//**************************************************************************
// RATING SIMULATIONS         BY ROYAL C. JONES            (C)1997
// simulation of ratings by three methods:  elo, linear, and ratio
// with statistics for deviation from predicted scores
// and correlation of calulated ranking with predefined ranking
//**************************************************************************

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#define PLAYERS 8000            /* number of players in population */
#define INTERVAL 10             /* number of games in sequential sample per
                                 * player */
#define GAMES 240000            /* total number of pairings in simulation */
#define GAMEGROUP 70            /* number of sequential samples per player */
#define RATINGGROUP 30          /* number of groups in rating distributions */
#define ELOSAMPLE 50            /* cutoff for linear and ratio performance
                                 * ratings */
#define LINEARSAMPLE 50
#define RATIOSAMPLE 50


double          totSqDevByGam_Lin[GAMEGROUP];   // total squared deviations
double          totSqDevByGam_Elo[GAMEGROUP];   // arranged by games played
 // per player
double          totSqDevByGam_Rat[GAMEGROUP];

double          stdByGam_Elo[GAMEGROUP];        // standard deviation
double          stdByGam_Rat[GAMEGROUP];
double          stdByGam_Lin[GAMEGROUP];

//grand total of squared deviations
double          totAllSqDev_Lin,
                totAllSqDev_Elo,
                totAllSqDev_Rat;

//grand total of standard deviations
double          totAllSTD_Lin,
                totAllSTD_Elo,
                totAllSTD_Rat;

int             cntByGam_Lin[GAMEGROUP];        // count by games played per
 // player
int             cntByGam_Elo[GAMEGROUP];
int             cntByGam_Rat[GAMEGROUP];

long            totCount;       // total count

int             distr_Elo[RATINGGROUP]; // distributions
int             distr_Lin[RATINGGROUP];
int             distr_Rat[RATINGGROUP];

//*************************************************************************
// player class:        stores individual ratings, score and game totals, and
//                opposition rating totals; for each game played, computes
//                expected scores, deviations, and ratings.
//*************************************************************************

class player {

 public:

    // ratings
    double          ideal;
    double          elo;
    double          linear;
    double          ratio;

    double          oldElo;     // pre-event rating
    double          oppRat;

    int             score,
                    rank;

    int             w;          // total points scored
    int             n;          // total games played

    // total of opposition ratings
    double          totOppElo;
    double          totOppLinear;
    double          totOppRatio;

    // predicted scores
    double          pctExpElo;
    double          pctExpLin;
    double          pctExpRat;

    double          deriv;      // derivative with respect to percentage
    // expectancy

                    player();
    void            FindSqDev_Elo();    // calculate squared deviation for
    // game
    void            FindSqDev_Lin();
    void            FindSqDev_Rat();
    void            RateElo();  // calculate rating
    void            RateLinear();
    void            RateRatio();

}               p[PLAYERS + 10];// array of player class

player::player()
{

    elo = 1400.0;
    linear = ratio = 0.5;
    n = w = 0;
    totOppElo = totOppLinear = totOppRatio = 0.0;


    ideal = 0;
    oldElo = 0;
    oppRat = 0;
    score = 0;
    rank = 0;
    pctExpElo = 0;
    pctExpLin = 0;
    pctExpRat = 0;
    deriv = 0;

}

void            player::FindSqDev_Elo()
{
    // calculate and record squared deviation for elo rating
    int             j;
    double          sqDev = (score - pctExpElo) * (score - pctExpElo);

//find slot (number of games) for squared deviation
    j = (n - 1) / INTERVAL;
    assert(j < GAMEGROUP);

    totSqDevByGam_Elo[j] += sqDev;      // add into global array
    cntByGam_Elo[j]++;          // bump count for this slot

    totAllSqDev_Elo += sqDev;   // add into grand total
}

void            player::FindSqDev_Lin()
{
    // calculate and record squared deviation for linear rating
    // see comments for FindSqDev_Elo()
    int             j;
    double          sqDev = (score - pctExpLin) * (score - pctExpLin);

    j = (n - 1) / INTERVAL;
    assert(j < GAMEGROUP);
    totSqDevByGam_Lin[j] += sqDev;
    cntByGam_Lin[j]++;

    totAllSqDev_Lin += sqDev;
}

void            player::FindSqDev_Rat()
{
    // calculate and record squared deviation for rectangular rating
    // see comments for FindSqDev_Elo()
    int             j;
    double          sqDev = (score - pctExpRat) * (score - pctExpRat);

    j = (n - 1) / INTERVAL;
    assert(j < GAMEGROUP);
    totSqDevByGam_Rat[j] += sqDev;
    cntByGam_Rat[j]++;

    totAllSqDev_Rat += sqDev;
}

void            player::RateElo()
{
    // computes elo rating by official description of USCF system

    // performance
    assert(n > 0);
    if (n < 20)
        elo = (totOppElo + 400.0 * (w - (n - w))) / n;
    else
        elo += 16 * (score - pctExpElo);

// adjustment for crossing K factor
//              if(oldElo < 2100 && elo >= 2100) elo = 2100 + (elo - 2100)*.75;
//              else if(oldElo >= 2100 && elo < 2100) elo = 2100 + (elo - 2100)*1.33;
//              else if(oldElo < 2400 && elo >= 2400) elo = 2400 + (elo - 2400)*.66;
//              else if(oldElo >= 2400 && elo < 2400) elo = 2400 + (elo - 2400)*1.5;
}

void            player::RateLinear()
{
    // compute linear rating

    if (n < LINEARSAMPLE) {
        // performance
        assert(n > 0);
        linear = (totOppLinear + 2.0 * w) / n - 1;
        // R = Rc + 2P -1
    }
    // established
    else
        linear += .5 * (score - pctExpLin) / LINEARSAMPLE;
    // change(R) = (W - We)/25
}

void            player::RateRatio()
{
    // computes ratio rating

    assert(n > 0);
    // average of opposition ratings:
    double          aveOpp = totOppRatio / n;

    // percentage score:
    double          pct = (double) w / n;

    if (aveOpp == 0 || pct == 0)
        ratio = 0;
    else if (aveOpp == 1 || pct == 1)
        ratio = 1;
    else if (n < RATIOSAMPLE || ratio == 1 || ratio == 0) {     // performance
        if (aveOpp <= 1 - pct)
            ratio = aveOpp * ((pct / (1 - pct)));
        // R = Rc(P)/(1-P)
        else
            ratio = (1 - (1 - aveOpp) * sqrt(((1 - pct) / pct)));
        // R = 1 - (1-Rc)(1-P)/(P)
    } else {                    // established rating
        if (ratio < .5)
            ratio += 2 * ratio * (score - pctExpRat) / RATIOSAMPLE;
        // change(R) = 2(R)(W-We)/50
        else
            ratio += 2 * (1 - ratio) * (score - pctExpRat) / RATIOSAMPLE;
        // change(R) = 2(1-R)(W-We)/50
        assert(ratio >= 0);
        assert(ratio <= 1);
    }
}

//**************************************************************************
// pairing class:    selects pairing randomly, updates opposition totals
//                   for each player, calculates expected scores, generates
//                   a random result, and calls player functions for
//                   squared deviation and rating
//**************************************************************************

class pairing {

 public:

    player * p1, *p2;

    void            Select();   // select random pairing
    void            UpdateOpp();// update opposition totals for each player

    // calculate expected score
    void            EloOdds();  // Elo
    void            LinearOdds();       // linear
    void            RatioOdds();// ratio

    void            RandomResult();     // result according to predefined
    // expectations

    void            Rate();

}               pair;

void            pairing::EloOdds()
{
    // calculates predicted scores from Elo formula (Elo, p.141)

    p1->pctExpElo = 1 / (1 + (pow(10.0, (p2->elo - p1->elo) / 400.0)));
    p2->pctExpElo = 1 - p1->pctExpElo;
}

void            pairing::LinearOdds()
{
    // calculates predicted scores from linear formula

    p1->pctExpLin = .5 * (p1->linear - p2->linear + 1.0);
    // Pe = .5(R - Rc) + .5
    if (p1->pctExpLin > 1.0)
        p1->pctExpLin = 1.0;
    else if (p1->pctExpLin < 0.0)
        p1->pctExpLin = 0.0;
    // 0 <= Pe <= 1
    p2->pctExpLin = 1 - p1->pctExpLin;
    // Pe(R) + Pe(Rc) = 1
}

void            pairing::RatioOdds()
{
    // calculates predicted scores from ratio formula

    if (p1->ratio == p2->ratio)
        p1->pctExpRat = .5;
    // if R = Rc, Pe(R) = Pe(Rc) = .5
    else if (p1->ratio + p2->ratio < 1.0)
        p1->pctExpRat = p1->ratio / (p1->ratio + p2->ratio);
    // if mean(R, Rc) < .5, Pe = (R)/(R - Rc)
    else
        p1->pctExpRat = (1 - p2->ratio) / (2 - p1->ratio - p2->ratio);
    // Pe = (1 - Rc)/((1 - R)+(1 - Rc))
    p2->pctExpRat = 1 - p1->pctExpRat;
    // Pe(R) + Pe(Rc) = 1
}

void            pairing::RandomResult()
{
    // random result according to ideal odds
    double          prob,
                    random;

    prob = 1 / (1 + pow(10.0, (p2->ideal - p1->ideal) / 400.0));        // (Elo, p.141)
    random = (double) rand() / RAND_MAX;        // in interval 0 < random <
    // 1;
    if (random < prob) {
        p1->score = 1;
        p2->score = 0;
        p1->w++;
    } else {                    // upset
        p1->score = 0;
        p2->score = 1;
        p2->w++;
    }
    p1->n++;                    // increment game totals
    p2->n++;
    totCount += 2;
}

void            pairing::UpdateOpp()
{
    // update opposition scores for performance ratings

    p1->totOppElo += p2->elo;
    p1->totOppLinear += p2->linear;
    p1->totOppRatio += p2->ratio;

    p2->totOppElo += p1->elo;
    p2->totOppLinear += p1->linear;
    p2->totOppRatio += p1->ratio;
}

void            pairing::Select()
{
    // select random pair; make sure players are distinct
    int             i,
                    j;

    i = rand() % PLAYERS;       // select player from full population
    j = rand() % (PLAYERS - 1); // select opponent from remaining population

    if (j >= i)
        j++;

    p1 = p + i;                 // convert index to pointer
    p2 = p + j;
}

void            pairing::Rate()
{

    p1->FindSqDev_Elo();        // find square deviations
    p2->FindSqDev_Elo();
    p1->FindSqDev_Lin();
    p2->FindSqDev_Lin();
    p1->FindSqDev_Rat();
    p2->FindSqDev_Rat();

    p1->RateElo();              // compute ratings
    p2->RateElo();
    p1->RateLinear();
    p2->RateLinear();
    p1->oppRat = p2->ratio;
    p2->oppRat = p1->ratio;
    p1->RateRatio();
    p2->RateRatio();
}

void            StdErr(void)
{
// compute standard deviation of expected from actual scores
// (this is not the standard error of linear correlation,
// but a useful statistic nonetheless)
    int             i;

    for (i = 0; i < GAMEGROUP; i++) {
        if (cntByGam_Elo[i] > 0)
            stdByGam_Elo[i] = sqrt(totSqDevByGam_Elo[i] / cntByGam_Elo[i]);
        if (cntByGam_Lin[i] > 0)
            stdByGam_Lin[i] = sqrt(totSqDevByGam_Lin[i] / cntByGam_Lin[i]);
        if (cntByGam_Rat[i] > 0)
            stdByGam_Rat[i] = sqrt(totSqDevByGam_Rat[i] / cntByGam_Rat[i]);
    }
    assert(totCount > 0);
    totAllSTD_Elo = sqrt(totAllSqDev_Elo / totCount);
    totAllSTD_Lin = sqrt(totAllSqDev_Lin / totCount);
    totAllSTD_Rat = sqrt(totAllSqDev_Rat / totCount);
}

double          Correlate(int sys)
{
    // caclulate the Spearman rank-difference correlation statistic
    // for the indicated rating system against predefined ranking
    int             i,
                    j,
                    n,
                    ranki,
                    rankj,
                    j_outranks_i = 0;
    double          sumSq = 0;

    n = PLAYERS;
    for (i = 0; i < n; i++)
        p[i].rank = i;
    for (i = 0; i < n - 1; i++) {
        for (j = i + 1; j < n; j++) {   // insertion sort
            ranki = p[i].rank;
            rankj = p[j].rank;
            switch (sys) {
            case 0:
                j_outranks_i = p[rankj].elo < p[ranki].elo;
                break;
            case 1:
                j_outranks_i = p[rankj].linear < p[ranki].linear;
                break;
            case 2:
                j_outranks_i = p[rankj].ratio < p[ranki].ratio;
                break;
            default:
                assert(sys < 3);// no more than three systems
            }
            if (j_outranks_i) {
                p[i].rank = rankj;
                p[j].rank = ranki;
            }
        }
    }

    for (i = 0; i < n; i++)
        sumSq += (p[i].rank - i) * (p[i].rank - i);
    return 1 - 6.0 * sumSq / (n * (n * n - 1.0));
}

void            Distribution(void)
{
    // rough distribution of ratings
    // mainly for internal checking, not a test criterion
    int             i,
                    j;

    for (i = 0; i < PLAYERS; i++) {
        j = (int) (p[i].elo / 100);
        if (j >= 0 && j < RATINGGROUP)
            ++distr_Elo[j];
        j = (int) (p[i].linear * 10) + 10;
        if (j >= 0 && j < RATINGGROUP)
            ++distr_Lin[j];
        j = (int) (p[i].ratio * 10) + 10;
        if (j >= 0 && j < RATINGGROUP)
            ++distr_Rat[j];
    }
}

int             main(void)
{

    long            i;
    int             j,
                    n;

    for (i = 0; i < PLAYERS; i++)       // initialize ideal ratings
        p[i].ideal = 1000 + i;

    n = PLAYERS * INTERVAL;
    printf("\nSpearman Rank Difference Correlation\n");
    printf("\n          elo   linear   ratio\n");
    for (i = 0; i < GAMES; i += n) {
        for (j = 0; j < n; j++) {
            pair.Select();
            pair.RandomResult();
            pair.UpdateOpp();
            pair.EloOdds();
            pair.LinearOdds();
            pair.RatioOdds();
            pair.Rate();
        }
        // output correlations every n games
        printf("%6d %7.4f %7.4f %7.4f\n",
               i + n, Correlate(0), Correlate(1), Correlate(2));
    }
    printf("\nRating distributions\n");
    Distribution();

    for (i = 0; i < RATINGGROUP; i++)
        printf("%2d: %4d %4d %4d\n",
               i, distr_Elo[i], distr_Lin[i], distr_Rat[i]);

    StdErr();

    // output standard deviations at conclusion of input
    // representing sample size with respect to individuals

    printf("\nStandard deviations, actual from expected scores\n");
    for (i = 0; i < GAMEGROUP; i++)
        printf("%3i: %6.4f %6.4f %6.4f\n",
           10 * (i + 1), stdByGam_Elo[i], stdByGam_Lin[i], stdByGam_Rat[i]);

    printf("\nTotal Squared Deviations: Elo = %f", totAllSqDev_Elo);
    printf(" Linear = %f", totAllSqDev_Lin);
    printf(" Ratio = %f\n", totAllSqDev_Rat);

    printf("std Elo = %f", totAllSTD_Elo);
    printf(" std Linear = %f", totAllSTD_Lin);
    printf(" std Ratio = %f\n", totAllSTD_Rat);

//  exit(0);

    printf("\nElo Ratings\n");
    for (i = 0; i < PLAYERS; i += 10) {
        for (j = 0; j < 10; j++)
            printf("%4i ", (int) p[i + j].elo);
        printf("\n");
    }

    printf("Linear Ratings\n");
    for (i = 0; i < PLAYERS; i += 10) {
        for (j = 0; j < 10; j++)
            printf("%4i ", (int) (p[i + j].linear * 1000 + 500));
        printf("\n");
    }

    printf("Ratio Ratings\n");
    for (i = 0; i < PLAYERS; i += 10) {
        for (j = 0; j < 10; j++)
            printf("%4i ", (int) (p[i + j].ratio * 10000));
        printf("\n");
    }
    return 0;
}