git.sesse.net Git - foosball/blob - foosrank.cpp

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <assert.h>
   4
   5 #include <vector>
   6 #include <algorithm>
   7
   8 #include <complex>
   9 #include <fftw3.h>
  10
  11 // step sizes
  12 static const double int_step_size = 75.0;
  13
  14 // rating constant (see below)
  15 static const double rating_constant = 455.0;
  16
  17 using namespace std;
  18
  19 double prob_score(int k, double a, double rd);
  20 double prob_score_real(int k, double a, double binomial, double rd_norm);
  21 double prodai(int k, double a);
  22 double fac(int x);
  23
  24
  25 // probability of match ending k-a (k>a) when winnerR - loserR = RD
  26 //
  27 //   +inf
  28 //     /
  29 //    |
  30 //    | Poisson[lambda1, t](a) * Erlang[lambda2, k](t) dt
  31 //    |
  32 //   /
  33 // -inf
  34 //
  35 // where lambda1 = 1.0, lambda2 = 2^(rd/455)
  36 //
  37 // The constant of 455 is chosen carefully so to match with the
  38 // Glicko/Bradley-Terry assumption that a player rated 400 points over
  39 // his/her opponent will win with a probability of 10/11 =~ 0.90909.
  40 //
  41 double prob_score(int k, double a, double rd)
  42 {
  43         return prob_score_real(k, a, prodai(k, a) / fac(k-1), rd/rating_constant);
  44 }
  45
  46 // Same, but takes in binomial(a+k-1, k-1) as an argument in
  47 // addition to a. Faster if you already have that precomputed, and assumes rd
  48 // is already divided by 455.
  49 double prob_score_real(int k, double a, double binomial, double rd_norm)
  50 {
  51         double nom = binomial * pow(2.0, rd_norm * a);
  52         double denom = pow(1.0 + pow(2.0, rd_norm), k+a);
  53         return nom/denom;
  54 }
  55
  56 // Calculates Product(a+i, i=1..k-1) (see above).
  57 double prodai(int k, double a)
  58 {
  59         double prod = 1.0;
  60         for (int i = 1; i < k; ++i)
  61                 prod *= (a+i);
  62         return prod;
  63 }
  64
  65 double fac(int x)
  66 {
  67         double prod = 1.0;
  68         for (int i = 2; i <= x; ++i)
  69                 prod *= i;
  70         return prod;
  71 }
  72
  73 //
  74 // Computes the integral
  75 //
  76 //   +inf
  77 //    /
  78 //    |
  79 //    | ProbScore[a] (r1-r2) Gaussian[mu2, sigma2] (r2) dr2
  80 //    |
  81 //   /
  82 // -inf
  83 //
  84 // For practical reasons, -inf and +inf are replaced by 0 and 3000, which
  85 // is reasonable in the this context.
  86 //
  87 // The Gaussian is not normalized.
  88 //
  89 // Set the last parameter to 1.0 if player 1 won, or -1.0 if player 2 won.
  90 // In the latter case, ProbScore will be given (r2-r1) instead of (r1-r2).
  91 //
  92 class ProbScoreEvaluator {
  93 private:
  94         int k;
  95         double a;
  96         double binomial_precompute, r1, mu2, sigma2, winfac;
  97
  98 public:
  99         ProbScoreEvaluator(int k, double a, double binomial_precompute, double r1, double mu2, double sigma2, double winfac)
 100                 : k(k), a(a), binomial_precompute(binomial_precompute), r1(r1), mu2(mu2), sigma2(sigma2), winfac(winfac) {}
 101         inline double operator() (double x) const
 102         {
 103                 double probscore = prob_score_real(k, a, binomial_precompute, (r1 - x)*winfac);
 104                 double z = (x - mu2)/sigma2;
 105                 double gaussian = exp(-(z*z/2.0));
 106                 return probscore * gaussian;
 107         }
 108 };
 109
 110 void convolve(int size)
 111 {
 112 }
 113
 114 void compute_opponent_rating_pdf(int k, double a, double mu2, double sigma2, double winfac, vector<pair<double, double> > &result)
 115 {
 116         double binomial_precompute = prodai(k, a) / fac(k-1);
 117         winfac /= rating_constant;
 118
 119         int sz = (6000.0 - 0.0) / int_step_size;
 120         double h = (6000.0 - 0.0) / sz;
 121
 122         fftw_plan f1, f2, b;
 123         complex<double> *func1, *func2, *res;
 124
 125         func1 = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
 126         func2 = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
 127         res = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
 128         f1 = fftw_plan_dft_1d(sz*2,
 129                 reinterpret_cast<fftw_complex*>(func1),
 130                 reinterpret_cast<fftw_complex*>(func1),
 131                 FFTW_FORWARD,
 132                 FFTW_MEASURE);
 133         f2 = fftw_plan_dft_1d(sz*2,
 134                 reinterpret_cast<fftw_complex*>(func2),
 135                 reinterpret_cast<fftw_complex*>(func2),
 136                 FFTW_FORWARD,
 137                 FFTW_MEASURE);
 138         b = fftw_plan_dft_1d(sz*2,
 139                 reinterpret_cast<fftw_complex*>(res),
 140                 reinterpret_cast<fftw_complex*>(res),
 141                 FFTW_BACKWARD,
 142                 FFTW_MEASURE);
 143
 144         // start off by zero
 145         for (int i = 0; i < sz*2; ++i) {
 146                 func1[i].real() = func1[i].imag() = func2[i].real() = func2[i].imag() = 0.0;
 147         }
 148
 149         for (int i = 0; i < sz; ++i) {
 150                 double x1 = 0.0 + h*i;
 151                 double z = (x1 - mu2)/sigma2;
 152                 func1[i].real() = exp(-(z*z/2.0));
 153
 154                 double x2 = -3000.0 + h*i;
 155                 func2[(i - sz/2 + sz*2)%(sz*2)].real() = prob_score_real(k, a, binomial_precompute, x2*winfac);
 156         }
 157
 158         result.reserve(sz*2);
 159
 160         // convolve
 161         fftw_execute(f1);
 162         fftw_execute(f2);
 163         for (int i = 0; i < sz*2; ++i) {
 164                 res[i] = func1[i] * func2[i];
 165         }
 166         fftw_execute(b);
 167         for (int i = 0; i < sz; ++i) {
 168                 double r1 = i*h;
 169                 result.push_back(make_pair(r1, abs(res[i])));
 170         }
 171 }
 172
 173 // normalize the curve so we know that A ~= 1
 174 void normalize(vector<pair<double, double> > &curve)
 175 {
 176         double peak = 0.0;
 177         for (vector<pair<double, double> >::const_iterator i = curve.begin(); i != curve.end(); ++i) {
 178                 peak = max(peak, i->second);
 179         }
 180
 181         double invpeak = 1.0 / peak;
 182         for (vector<pair<double, double> >::iterator i = curve.begin(); i != curve.end(); ++i) {
 183                 i->second *= invpeak;
 184         }
 185 }
 186
 187 // computes matA * matB
 188 void mat_mul(double *matA, unsigned ah, unsigned aw,
 189              double *matB, unsigned bh, unsigned bw,
 190              double *result)
 191 {
 192         assert(aw == bh);
 193         for (unsigned y = 0; y < bw; ++y) {
 194                 for (unsigned x = 0; x < ah; ++x) {
 195                         double sum = 0.0;
 196                         for (unsigned c = 0; c < aw; ++c) {
 197                                 sum += matA[c*ah + x] * matB[y*bh + c];
 198                         }
 199                         result[y*bw + x] = sum;
 200                 }
 201         }
 202 }
 203
 204 // computes matA^T * matB
 205 void mat_mul_trans(double *matA, unsigned ah, unsigned aw,
 206                    double *matB, unsigned bh, unsigned bw,
 207                    double *result)
 208 {
 209         assert(ah == bh);
 210         for (unsigned y = 0; y < bw; ++y) {
 211                 for (unsigned x = 0; x < aw; ++x) {
 212                         double sum = 0.0;
 213                         for (unsigned c = 0; c < ah; ++c) {
 214                                 sum += matA[x*ah + c] * matB[y*bh + c];
 215                         }
 216                         result[y*bw + x] = sum;
 217                 }
 218         }
 219 }
 220
 221 void print3x3(double *M)
 222 {
 223         printf("%f %f %f\n", M[0], M[3], M[6]);
 224         printf("%f %f %f\n", M[1], M[4], M[7]);
 225         printf("%f %f %f\n", M[2], M[5], M[8]);
 226 }
 227
 228 void print3x1(double *M)
 229 {
 230         printf("%f\n", M[0]);
 231         printf("%f\n", M[1]);
 232         printf("%f\n", M[2]);
 233 }
 234
 235 // solves Ax = B by Gauss-Jordan elimination, where A is a 3x3 matrix,
 236 // x is a column vector of length 3 and B is a row vector of length 3.
 237 // Destroys its input in the process.
 238 void solve3x3(double *A, double *x, double *B)
 239 {
 240         // row 1 -= row 0 * (a1/a0)
 241         {
 242                 double f = A[1] / A[0];
 243                 A[1] = 0.0;
 244                 A[4] -= A[3] * f;
 245                 A[7] -= A[6] * f;
 246
 247                 B[1] -= B[0] * f;
 248         }
 249
 250         // row 2 -= row 0 * (a2/a0)
 251         {
 252                 double f = A[2] / A[0];
 253                 A[2] = 0.0;
 254                 A[5] -= A[3] * f;
 255                 A[8] -= A[6] * f;
 256
 257                 B[2] -= B[0] * f;
 258         }
 259
 260         // row 2 -= row 1 * (a5/a4)
 261         {
 262                 double f = A[5] / A[4];
 263                 A[5] = 0.0;
 264                 A[8] -= A[7] * f;
 265
 266                 B[2] -= B[1] * f;
 267         }
 268
 269         // back substitute:
 270
 271         // row 1 -= row 2 * (a7/a8)
 272         {
 273                 double f = A[7] / A[8];
 274                 A[7] = 0.0;
 275
 276                 B[1] -= B[2] * f;
 277         }
 278
 279         // row 0 -= row 2 * (a6/a8)
 280         {
 281                 double f = A[6] / A[8];
 282                 A[6] = 0.0;
 283
 284                 B[0] -= B[2] * f;
 285         }
 286
 287         // row 0 -= row 1 * (a3/a4)
 288         {
 289                 double f = A[3] / A[4];
 290                 A[3] = 0.0;
 291
 292                 B[0] -= B[1] * f;
 293         }
 294
 295         // normalize
 296         x[0] = B[0] / A[0];
 297         x[1] = B[1] / A[4];
 298         x[2] = B[2] / A[8];
 299 }
 300
 301 // Give an OK starting estimate for the least squares, by numerical integration
 302 // of statistical moments.
 303 void estimate_musigma(vector<pair<double, double> > &curve, double &mu_result, double &sigma_result)
 304 {
 305         double h = (curve.back().first - curve.front().first) / (curve.size() - 1);
 306
 307         double area = curve.front().second;
 308         double ex = curve.front().first * curve.front().second;
 309         double ex2 = curve.front().first * curve.front().first * curve.front().second;
 310
 311         for (unsigned i = 1; i < curve.size() - 1; i += 2) {
 312                 double x = curve[i].first;
 313                 double y = curve[i].second;
 314                 area += 4.0 * y;
 315                 ex += 4.0 * x * y;
 316                 ex2 += 4.0 * x * x * y;
 317         }
 318         for (unsigned i = 2; i < curve.size() - 1; i += 2) {
 319                 double x = curve[i].first;
 320                 double y = curve[i].second;
 321                 area += 2.0 * y;
 322                 ex += 2.0 * x * y;
 323                 ex2 += 2.0 * x * x * y;
 324         }
 325
 326         area += curve.back().second;
 327         ex += curve.back().first * curve.back().second;
 328         ex2 += curve.back().first * curve.back().first * curve.back().second;
 329
 330         area = (h/3.0) * area;
 331         ex = (h/3.0) * ex / area;
 332         ex2 = (h/3.0) * ex2 / area;
 333
 334         mu_result = ex;
 335         sigma_result = sqrt(ex2 - ex * ex);
 336 }
 337
 338 // Find best fit of the data in curves to a Gaussian pdf, based on the
 339 // given initial estimates. Works by nonlinear least squares, iterating
 340 // until we're below a certain threshold.
 341 //
 342 // Note that the algorithm blows up quite hard if the initial estimate is
 343 // not good enough. Use estimate_musigma to get a reasonable starting
 344 // estimate.
 345 void least_squares(vector<pair<double, double> > &curve, double mu1, double sigma1, double &mu_result, double &sigma_result)
 346 {
 347         double A = 1.0;
 348         double mu = mu1;
 349         double sigma = sigma1;
 350
 351         // column-major
 352         double matA[curve.size() * 3];  // N x 3
 353         double dbeta[curve.size()];     // N x 1
 354
 355         // A^T * A: 3xN * Nx3 = 3x3
 356         double matATA[3*3];
 357
 358         // A^T * dβ: 3xN * Nx1 = 3x1
 359         double matATdb[3];
 360
 361         double dlambda[3];
 362
 363         for ( ;; ) {
 364                 //printf("A=%f mu=%f sigma=%f\n", A, mu, sigma);
 365
 366                 // fill in A (depends only on x_i, A, mu, sigma -- not y_i)
 367                 for (unsigned i = 0; i < curve.size(); ++i) {
 368                         double x = curve[i].first;
 369
 370                         // df/dA(x_i)
 371                         matA[i + 0 * curve.size()] =
 372                                 exp(-(x-mu)*(x-mu)/(2.0*sigma*sigma));
 373
 374                         // df/dµ(x_i)
 375                         matA[i + 1 * curve.size()] =
 376                                 A * (x-mu)/(sigma*sigma) * matA[i + 0 * curve.size()];
 377
 378                         // df/dσ(x_i)
 379                         matA[i + 2 * curve.size()] =
 380                                 matA[i + 1 * curve.size()] * (x-mu)/sigma;
 381                 }
 382
 383                 // find dβ
 384                 for (unsigned i = 0; i < curve.size(); ++i) {
 385                         double x = curve[i].first;
 386                         double y = curve[i].second;
 387
 388                         dbeta[i] = y - A * exp(- (x-mu)*(x-mu)/(2.0*sigma*sigma));
 389                 }
 390
 391                 // compute a and b
 392                 mat_mul_trans(matA, curve.size(), 3, matA, curve.size(), 3, matATA);
 393                 mat_mul_trans(matA, curve.size(), 3, dbeta, curve.size(), 1, matATdb);
 394
 395                 // solve
 396                 solve3x3(matATA, dlambda, matATdb);
 397
 398                 A += dlambda[0];
 399                 mu += dlambda[1];
 400                 sigma += dlambda[2];
 401
 402                 // terminate when we're down to three digits
 403                 if (fabs(dlambda[0]) <= 1e-3 && fabs(dlambda[1]) <= 1e-3 && fabs(dlambda[2]) <= 1e-3)
 404                         break;
 405         }
 406
 407         mu_result = mu;
 408         sigma_result = sigma;
 409 }
 410
 411 void compute_new_rating(double mu1, double sigma1, double mu2, double sigma2, int score1, int score2, double &mu, double &sigma)
 412 {
 413         vector<pair<double, double> > curve;
 414
 415         if (score1 > score2) {
 416                 compute_opponent_rating_pdf(score1, score2, mu2, sigma2, -1.0, curve);
 417         } else {
 418                 compute_opponent_rating_pdf(score2, score1, mu2, sigma2, 1.0, curve);
 419         }
 420
 421         // multiply in the gaussian
 422         for (unsigned i = 0; i < curve.size(); ++i) {
 423                 double r1 = curve[i].first;
 424                 double z = (r1 - mu1) / sigma1;
 425                 double gaussian = exp(-(z*z/2.0));
 426                 curve[i].second *= gaussian;
 427         }
 428
 429         double mu_est, sigma_est;
 430         normalize(curve);
 431         estimate_musigma(curve, mu_est, sigma_est);
 432         least_squares(curve, mu_est, sigma_est, mu, sigma);
 433 }
 434
 435 void compute_new_double_rating(double mu1, double sigma1, double mu2, double sigma2, double mu3, double sigma3, double mu4, double sigma4, int score1, int score2, double &mu, double &sigma)
 436 {
 437         vector<pair<double, double> > curve, newcurve;
 438         double mu_t = mu3 + mu4;
 439         double sigma_t = sqrt(sigma3*sigma3 + sigma4*sigma4);
 440
 441         if (score1 > score2) {
 442                 compute_opponent_rating_pdf(score1, score2, mu_t, sigma_t, -1.0, curve);
 443         } else {
 444                 compute_opponent_rating_pdf(score2, score1, mu_t, sigma_t, 1.0, curve);
 445         }
 446
 447         // iterate over r1
 448         double h = 3000.0 / curve.size();
 449         for (unsigned i = 0; i < curve.size(); ++i) {
 450                 double sum = 0.0;
 451
 452                 // could be anything, but this is a nice start
 453                 //double r1 = curve[i].first;
 454                 double r1 = i * h;
 455
 456                 // iterate over r2
 457                 for (unsigned j = 0; j < curve.size(); ++j) {
 458                         double r1plusr2 = curve[j].first;
 459                         double r2 = r1plusr2 - r1;
 460
 461                         double z = (r2 - mu2) / sigma2;
 462                         double gaussian = exp(-(z*z/2.0));
 463                         sum += curve[j].second * gaussian;
 464                 }
 465
 466                 double z = (r1 - mu1) / sigma1;
 467                 double gaussian = exp(-(z*z/2.0));
 468                 newcurve.push_back(make_pair(r1, gaussian * sum));
 469         }
 470
 471
 472         double mu_est, sigma_est;
 473         normalize(newcurve);
 474         estimate_musigma(newcurve, mu_est, sigma_est);
 475         least_squares(newcurve, mu_est, sigma_est, mu, sigma);
 476 }
 477
 478 int main(int argc, char **argv)
 479 {
 480         FILE *fp = fopen("fftw-wisdom", "rb");
 481         if (fp != NULL) {
 482                 fftw_import_wisdom_from_file(fp);
 483                 fclose(fp);
 484         }
 485
 486         double mu1 = atof(argv[1]);
 487         double sigma1 = atof(argv[2]);
 488         double mu2 = atof(argv[3]);
 489         double sigma2 = atof(argv[4]);
 490
 491         if (argc > 10) {
 492                 double mu3 = atof(argv[5]);
 493                 double sigma3 = atof(argv[6]);
 494                 double mu4 = atof(argv[7]);
 495                 double sigma4 = atof(argv[8]);
 496                 int score1 = atoi(argv[9]);
 497                 int score2 = atoi(argv[10]);
 498                 double mu, sigma;
 499                 compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, score1, score2, mu, sigma);
 500                 printf("%f %f\n", mu, sigma);
 501         } else if (argc > 8) {
 502                 double mu3 = atof(argv[5]);
 503                 double sigma3 = atof(argv[6]);
 504                 double mu4 = atof(argv[7]);
 505                 double sigma4 = atof(argv[8]);
 506                 int k = atoi(argv[9]);
 507
 508                 // assess all possible scores
 509                 for (int i = 0; i < k; ++i) {
 510                         double newmu1_1, newmu1_2, newmu2_1, newmu2_2;
 511                         double newsigma1_1, newsigma1_2, newsigma2_1, newsigma2_2;
 512                         compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, k, i, newmu1_1, newsigma1_1);
 513                         compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, k, i, newmu1_2, newsigma1_2);
 514                         compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, i, k, newmu2_1, newsigma2_1);
 515                         compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, i, k, newmu2_2, newsigma2_2);
 516                         printf("%u-%u,%f,%+f,%+f,%+f,%+f\n",
 517                                 k, i, prob_score(k, i, mu3+mu4-(mu1+mu2)), newmu1_1-mu1, newmu1_2-mu2,
 518                                 newmu2_1-mu3, newmu2_2-mu4);
 519                 }
 520                 for (int i = k; i --> 0; ) {
 521                         double newmu1_1, newmu1_2, newmu2_1, newmu2_2;
 522                         double newsigma1_1, newsigma1_2, newsigma2_1, newsigma2_2;
 523                         compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, i, k, newmu1_1, newsigma1_1);
 524                         compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, i, k, newmu1_2, newsigma1_2);
 525                         compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, k, i, newmu2_1, newsigma2_1);
 526                         compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, k, i, newmu2_2, newsigma2_2);
 527                         printf("%u-%u,%f,%+f,%+f,%+f,%+f\n",
 528                                 i, k, prob_score(k, i, mu1+mu2-(mu3+mu4)), newmu1_1-mu1, newmu1_2-mu2,
 529                                 newmu2_1-mu3, newmu2_2-mu4);
 530                 }
 531         } else if (argc > 6) {
 532                 int score1 = atoi(argv[5]);
 533                 int score2 = atoi(argv[6]);
 534                 double mu, sigma;
 535                 compute_new_rating(mu1, sigma1, mu2, sigma2, score1, score2, mu, sigma);
 536                 printf("%f %f\n", mu, sigma);
 537         } else {
 538                 int k = atoi(argv[5]);
 539
 540                 // assess all possible scores
 541                 for (int i = 0; i < k; ++i) {
 542                         double newmu1, newmu2, newsigma1, newsigma2;
 543                         compute_new_rating(mu1, sigma1, mu2, sigma2, k, i, newmu1, newsigma1);
 544                         compute_new_rating(mu2, sigma2, mu1, sigma1, i, k, newmu2, newsigma2);
 545                         printf("%u-%u,%f,%+f,%+f\n",
 546                                 k, i, prob_score(k, i, mu2-mu1), newmu1-mu1, newmu2-mu2);
 547                 }
 548                 for (int i = k; i --> 0; ) {
 549                         double newmu1, newmu2, newsigma1, newsigma2;
 550                         compute_new_rating(mu1, sigma1, mu2, sigma2, i, k, newmu1, newsigma1);
 551                         compute_new_rating(mu2, sigma2, mu1, sigma1, k, i, newmu2, newsigma2);
 552                         printf("%u-%u,%f,%+f,%+f\n",
 553                                 i, k, prob_score(k, i, mu1-mu2), newmu1-mu1, newmu2-mu2);
 554                 }
 555         }
 556
 557         fp = fopen("fftw-wisdom", "wb");
 558         if (fp != NULL) {
 559                 fftw_export_wisdom_to_file(fp);
 560                 fclose(fp);
 561         }
 562 }
 563