git.sesse.net Git - wloh/blob - bayeswf.cpp

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <string.h>
   4 #include <stdlib.h>
   5 #include <Eigen/Core>
   6 #include <Eigen/Eigenvalues>
   7
   8 #include <map>
   9 #include <vector>
  10 #include <string>
  11 #include <algorithm>
  12
  13 using namespace std;
  14 using namespace Eigen;
  15
  16 #define PRIOR_MU 500
  17 #define PRIOR_WEIGHT 1.0
  18 #define MAX_PLAYERS 8192
  19 #define DUMP_RAW 0
  20 #define USE_DB 1
  21
  22 #if USE_DB
  23 #include <pqxx/connection>
  24 #include <pqxx/tablewriter>
  25 #include <pqxx/transaction>
  26 #endif
  27
  28 float mu[MAX_PLAYERS];
  29 float mu_stddev[MAX_PLAYERS];
  30 float global_sigma;
  31 float prior_sigma;
  32
  33 // Data waiting for insertion into the database.
  34
  35 struct RatingDBTuple {
  36         int player;
  37         float mu, mu_stddev;
  38 };
  39 struct CovarianceDBTuple {
  40         int player1, player2;
  41         float covariance;
  42 };
  43 vector<RatingDBTuple> rating_db_tuples;
  44 vector<CovarianceDBTuple> covariance_db_tuples;
  45 map<pair<string, string>, float> aux_params;
  46
  47 #define EPSILON 1e-3
  48
  49 /*
  50  * L(mu_vec, sigma_vec, matches) = product[ L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ]
  51  * log-likelihood = sum[ log( L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ) ]
  52  *
  53  * L(mu1, sigma1, mu2, sigma2, score2 - score1) = sigmoid(mu2 - mu1, sqrt(sigma1² + sigma2²), (score2 - score1))
  54  *
  55  * pdf := 1/(sigma * sqrt(2*Pi)) * exp(-(x - mu)^2 / (2 * sigma^2));
  56  * pdfs := subs({ mu = mu1 - mu2, sigma = sqrt(sigma1^2 + sigma2^2) }, pdf);
  57  * diff(log(pdfs), mu1);
  58  */
  59
  60 struct match {
  61         int player, other_player;
  62         int margin;
  63         float weight;
  64 };
  65 map<int, vector<match> > matches_for_player;
  66 vector<match> all_matches;
  67
  68 void dump_scores(const vector<int> &players, const float *mu, const float *mu_stddev, int num_players)
  69 {
  70 #if USE_DB
  71         for (int i = 0; i < num_players; ++i) {
  72                 RatingDBTuple tuple = { players[i], mu[i], mu_stddev[i] };
  73                 rating_db_tuples.push_back(tuple);
  74         }
  75 #else
  76         for (int i = 0; i < num_players; ++i) {
  77                 printf("%f %f %d\n", mu[i], mu_stddev[i], players[i]);
  78         }
  79 #endif
  80 }
  81
  82 /*
  83  * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2
  84  * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
  85  *                                sum_i[ (w_i/sigma_c_i)^2 mu1 ] = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ]
  86  *                                mu1 = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (w_i/sigma_c_i)^2 ]
  87  */
  88 void update_mu(float *mu, int player_num, const vector<match> &matches)
  89 {
  90         if (matches.empty()) {
  91                 return;
  92         }
  93
  94         float nom = 0.0f, denom = 0.0f;
  95
  96         // Prior.
  97         {
  98                 float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
  99                 nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
 100                 denom += PRIOR_WEIGHT * inv_sigma2;
 101         }
 102
 103         // All matches.
 104         for (unsigned i = 0; i < matches.size(); ++i) {
 105                 float inv_sigma_c2 = matches[i].weight / (global_sigma * global_sigma);
 106                 float x = matches[i].margin; // / 70.0f;
 107
 108                 nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
 109                 denom += inv_sigma_c2;
 110         }
 111         mu[player_num] = nom / denom;
 112 }
 113
 114 void dump_raw(const float *mu, int num_players)
 115 {
 116         for (unsigned i = 0; i < all_matches.size(); ++i) {
 117                 const match& m = all_matches[i];
 118
 119                 float mu1 = mu[m.player];
 120                 float mu2 = mu[m.other_player];
 121                 float sigma = global_sigma;
 122                 float mu = mu1 - mu2;
 123                 float x = m.margin;
 124                 float w = m.weight;
 125
 126                 printf("%f %f\n", (x - mu) / sigma, w);
 127         }
 128 }
 129
 130 /*
 131  * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
 132  * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
 133  *                                   sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0                            |: sigma != 0
 134  *                                   sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
 135  *                                   sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
 136  */
 137 void update_global_sigma(float *mu, int num_players)
 138 {
 139         float nom = 0.0f, denom = 0.0f;
 140         for (unsigned i = 0; i < all_matches.size(); ++i) {
 141                 const match& m = all_matches[i];
 142
 143                 float mu1 = mu[m.player];
 144                 float mu2 = mu[m.other_player];
 145                 float mu = mu1 - mu2;
 146                 float x = m.margin;
 147                 float w = m.weight;
 148
 149                 nom += w * ((x - mu) * (x - mu));
 150                 denom += w;
 151         }
 152
 153         global_sigma = sqrt(nom / denom);
 154 }
 155
 156 /*
 157  * diff(priorlogL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
 158  * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
 159  *                                   sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0                            |: sigma != 0
 160  *                                   sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
 161  *                                   sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
 162  */
 163 void update_prior_sigma(float *mu, int num_players)
 164 {
 165         float nom = 0.0f, denom = 0.0f;
 166         for (int i = 0; i < num_players; ++i) {
 167                 float mu1 = mu[i];
 168
 169                 nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
 170                 denom += 1.0f;
 171         }
 172
 173         prior_sigma = sqrt(nom / denom);
 174         if (!(prior_sigma > 40.0f)) {
 175                 prior_sigma = 40.0f;
 176         }
 177 }
 178
 179 float compute_logl(float z)
 180 {
 181         return -0.5 * (log(2.0f / M_PI) + z * z);
 182 }
 183
 184 float compute_total_logl(float *mu, int num_players)
 185 {
 186         float total_logl = 0.0f;
 187
 188         // Prior.
 189         for (int i = 0; i < num_players; ++i) {
 190                 total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
 191         }
 192
 193         // Matches.
 194         for (unsigned i = 0; i < all_matches.size(); ++i) {
 195                 const match& m = all_matches[i];
 196
 197                 float mu1 = mu[m.player];
 198                 float mu2 = mu[m.other_player];
 199                 float sigma = global_sigma;
 200                 float mu = mu1 - mu2;
 201                 float x = m.margin;
 202                 float w = m.weight;
 203
 204                 total_logl += w * compute_logl((x - mu) / sigma);
 205         }
 206
 207         return total_logl;
 208 }
 209
 210 /*
 211  * Compute Hessian matrix of the negative log-likelihood, ie. for each term in logL:
 212  *
 213  * M_ij = D_i D_j (- logL) = -w / sigma²                                for i != j
 214  *                            w / sigma²                                for i == j
 215  *
 216  * Note that this does not depend on mu or the margin at all.
 217  */
 218 Matrix<float, Dynamic, Dynamic> hessian;
 219 void construct_hessian(const float *mu, int num_players)
 220 {
 221         hessian = Matrix<float, Dynamic, Dynamic>(num_players, num_players);
 222         hessian.fill(0.0f);
 223
 224         for (int i = 0; i < num_players; ++i) {
 225                 hessian(i, i) += 1.0f / (prior_sigma * prior_sigma);
 226         }
 227         for (unsigned i = 0; i < all_matches.size(); ++i) {
 228                 const match &m = all_matches[i];
 229
 230                 int p1 = m.player;
 231                 int p2 = m.other_player;
 232
 233                 double sigma_sq = global_sigma * global_sigma;
 234                 float w = m.weight;
 235
 236                 hessian(p1, p2) -= w / sigma_sq;
 237                 hessian(p2, p1) -= w / sigma_sq;
 238
 239                 hessian(p1, p1) += w / sigma_sq;
 240                 hessian(p2, p2) += w / sigma_sq;
 241         }
 242 }
 243
 244 // Compute uncertainty (stddev) of mu estimates, which is sqrt((H^-1)_ii),
 245 // where H is the Hessian (see construct_hessian()).
 246 void compute_mu_uncertainty(const float *mu, const vector<int> &players)
 247 {
 248         // FIXME: Use pseudoinverse if applicable.
 249         Matrix<float, Dynamic, Dynamic> ih = hessian.inverse();
 250         for (unsigned i = 0; i < players.size(); ++i) {
 251                 mu_stddev[i] = sqrt(ih(i, i));
 252         }
 253
 254 #if USE_DB
 255         for (unsigned i = 0; i < players.size(); ++i) {
 256                 for (unsigned j = 0; j < players.size(); ++j) {
 257                         CovarianceDBTuple tuple;
 258                         tuple.player1 = players[i];
 259                         tuple.player2 = players[j];
 260                         tuple.covariance = ih(i, j);
 261                         covariance_db_tuples.push_back(tuple);
 262                 }
 263         }
 264 #else
 265         for (unsigned i = 0; i < players.size(); ++i) {
 266                 for (unsigned j = 0; j < players.size(); ++j) {
 267                         printf("covariance %d %d %f\n",
 268                                players[i],
 269                                players[j],
 270                                ih(i, j));
 271                 }
 272         }
 273 #endif
 274 }
 275
 276 void process_file(const char *filename)
 277 {
 278         global_sigma = 70.0f;
 279         prior_sigma = 70.0f;
 280         matches_for_player.clear();
 281         all_matches.clear();
 282
 283         FILE *fp = fopen(filename, "r");
 284         if (fp == NULL) {
 285                 perror(filename);
 286                 exit(1);
 287         }
 288
 289         char locale[256];
 290         if (fscanf(fp, "%s", locale) != 1) {
 291                 fprintf(stderr, "Could't read number of players\n");
 292                 exit(1);
 293         }
 294
 295         int num_players;
 296         if (fscanf(fp,"%d", &num_players) != 1) {
 297                 fprintf(stderr, "Could't read number of players\n");
 298                 exit(1);
 299         }
 300
 301         if (num_players > MAX_PLAYERS) {
 302                 fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS);
 303                 exit(1);
 304         }
 305
 306         vector<int> players;
 307         map<int, int> player_map;
 308
 309         for (int i = 0; i < num_players; ++i) {
 310                 char buf[256];
 311                 if (fscanf(fp, "%s", buf) != 1) {
 312                         fprintf(stderr, "Couldn't read player %d\n", i);
 313                         exit(1);
 314                 }
 315
 316                 players.push_back(atoi(buf));
 317                 player_map[atoi(buf)] = i;
 318         }
 319
 320         int num_matches = 0;
 321         for ( ;; ) {
 322                 int pl1, pl2;
 323                 int score1, score2;
 324                 float weight;
 325
 326                 if (fscanf(fp, "%d %d %d %d %f", &pl1, &pl2, &score1, &score2, &weight) != 5) {
 327                         //fprintf(stderr, "Read %d matches.\n", num_matches);
 328                         break;
 329                 }
 330
 331                 ++num_matches;
 332
 333                 if (player_map.count(pl1) == 0) {
 334                         fprintf(stderr, "Unknown player '%d'\n", pl1);
 335                         exit(1);
 336                 }
 337                 if (player_map.count(pl2) == 0) {
 338                         fprintf(stderr, "Unknown player '%d'\n", pl2);
 339                         exit(1);
 340                 }
 341
 342                 match m1;
 343                 m1.player = player_map[pl1];
 344                 m1.other_player = player_map[pl2];
 345                 m1.margin = score1 - score2;
 346                 m1.weight = weight;
 347                 matches_for_player[player_map[pl1]].push_back(m1);
 348
 349                 match m2;
 350                 m2.player = player_map[pl2];
 351                 m2.other_player = player_map[pl1];
 352                 m2.margin = score2 - score1;
 353                 m2.weight = weight;
 354                 matches_for_player[player_map[pl2]].push_back(m2);
 355
 356                 all_matches.push_back(m1);
 357         }
 358
 359         fclose(fp);
 360
 361         float mu[MAX_PLAYERS];
 362
 363         for (int i = 0; i < num_players; ++i) {
 364                 mu[i] = PRIOR_MU;
 365         }
 366
 367         int num_iterations = -1;
 368         for (int j = 0; j < 1000; ++j) {
 369                 float old_mu[MAX_PLAYERS];
 370                 float old_global_sigma = global_sigma;
 371                 float old_prior_sigma = prior_sigma;
 372                 memcpy(old_mu, mu, sizeof(mu));
 373                 for (int i = 0; i < num_players; ++i) {
 374                         update_mu(mu, i, matches_for_player[i]);
 375                 }
 376                 update_global_sigma(mu, num_players);
 377                 update_prior_sigma(mu, num_players);
 378                 /* for (int i = 0; i < num_players; ++i) {
 379                         update_sigma(mu, sigma, i, matches_for_player[i]);
 380                         dump_scores(players, mu, sigma, num_players);
 381                 } */
 382
 383                 float sumdiff = 0.0f;
 384                 for (int i = 0; i < num_players; ++i) {
 385                         sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
 386                 }
 387                 sumdiff += (prior_sigma - old_prior_sigma) * (prior_sigma - old_prior_sigma);
 388                 sumdiff += (global_sigma - old_global_sigma) * (global_sigma - old_global_sigma);
 389                 if (sumdiff < EPSILON) {
 390                         //fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
 391                         num_iterations = j + 1;
 392                         break;
 393                 }
 394         }
 395
 396         construct_hessian(mu, num_players);
 397         aux_params[make_pair(locale, "num_iterations")] = num_iterations;
 398         aux_params[make_pair(locale, "score_stddev")] = global_sigma / sqrt(2.0f);
 399         aux_params[make_pair(locale, "rating_prior_stddev")] = prior_sigma;
 400         aux_params[make_pair(locale, "total_log_likelihood")] = compute_total_logl(mu, num_players);
 401
 402         compute_mu_uncertainty(mu, players);
 403         dump_scores(players, mu, mu_stddev, num_players);
 404 }
 405
 406 int main(int argc, char **argv)
 407 {
 408 #if USE_DB
 409         pqxx::connection conn("dbname=wloh host=127.0.0.1 user=wloh password=censored");
 410 #endif
 411
 412         for (int i = 1; i < argc; ++i) {
 413                 process_file(argv[i]);
 414         }
 415
 416 #if DUMP_RAW
 417         dump_raw(mu, num_players);
 418         return 0;
 419 #endif
 420
 421 #if USE_DB
 422         pqxx::work txn(conn);
 423         txn.exec("SET client_min_messages TO WARNING");
 424
 425         // Dump aux_params.
 426         {
 427                 txn.exec("TRUNCATE aux_params");
 428                 pqxx::tablewriter writer(txn, "aux_params");
 429                 for (map<pair<string, string>, float>::const_iterator it = aux_params.begin(); it != aux_params.end(); ++it) {
 430                         char str[128];
 431                         snprintf(str, 128, "%f", it->second);
 432
 433                         vector<string> tuple;
 434                         tuple.push_back(it->first.first);  // locale
 435                         tuple.push_back(it->first.second);  // parameter name
 436                         tuple.push_back(str);
 437                         writer.push_back(tuple);
 438                 }
 439                 writer.complete();
 440         }
 441
 442         // Dump ratings.
 443         {
 444                 txn.exec("TRUNCATE ratings");
 445                 pqxx::tablewriter writer(txn, "ratings");
 446                 for (unsigned i = 0; i < rating_db_tuples.size(); ++i) {
 447                         char player_str[128], mu_str[128], mu_stddev_str[128];
 448                         snprintf(player_str, 128, "%d", rating_db_tuples[i].player);
 449                         snprintf(mu_str, 128, "%f", rating_db_tuples[i].mu);
 450                         snprintf(mu_stddev_str, 128, "%f", rating_db_tuples[i].mu_stddev);
 451
 452                         vector<string> tuple;
 453                         tuple.push_back(player_str);
 454                         tuple.push_back(mu_str);
 455                         tuple.push_back(mu_stddev_str);
 456                         writer.push_back(tuple);
 457                 }
 458                 writer.complete();
 459         }
 460
 461         // Create a table new_covariance, and dump covariance into it.
 462         {
 463                 txn.exec("CREATE TABLE new_covariance ( player1 smallint NOT NULL, player2 smallint NOT NULL, cov float NOT NULL )");
 464                 pqxx::tablewriter writer(txn, "new_covariance");
 465                 for (unsigned i = 0; i < covariance_db_tuples.size(); ++i) {
 466                         char player1_str[128], player2_str[128], cov_str[128];
 467                         snprintf(player1_str, 128, "%d", covariance_db_tuples[i].player1);
 468                         snprintf(player2_str, 128, "%d", covariance_db_tuples[i].player2);
 469                         snprintf(cov_str, 128, "%f", covariance_db_tuples[i].covariance);
 470
 471                         vector<string> tuple;
 472                         tuple.push_back(player1_str);
 473                         tuple.push_back(player2_str);
 474                         tuple.push_back(cov_str);
 475                         writer.push_back(tuple);
 476                 }
 477                 writer.complete();
 478         }
 479
 480         // Create index, and rename new_covariance on top of covariance.
 481         txn.exec("ALTER TABLE new_covariance ADD PRIMARY KEY ( player1, player2 );");
 482         txn.exec("DROP TABLE IF EXISTS covariance");
 483         txn.exec("ALTER TABLE new_covariance RENAME TO covariance");
 484 #else
 485         //fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
 486         for (map<pair<string, string>, float>::const_iterator it = aux_params.begin(); it != aux_params.end(); ++it) {
 487                 printf("%s: aux_param %s %f\n", it->first.first.c_str(), it->first.second.c_str(), it->second);
 488         }
 489 #endif
 490
 491         txn.commit();
 492 }