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
   6 #include <map>
   7 #include <vector>
   8 #include <string>
   9 #include <algorithm>
  10
  11 using namespace std;
  12
  13 #define PRIOR_MU 1500
  14 #define PRIOR_WEIGHT 1.0
  15 #define MAX_PLAYERS 4096
  16 #define DUMP_RAW 0
  17
  18 float mu[MAX_PLAYERS];
  19 float global_sigma = 70.0f;
  20 float prior_sigma = 70.0f;
  21
  22 #define EPSILON 1e-3
  23
  24 /*
  25  * L(mu_vec, sigma_vec, matches) = product[ L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ]
  26  * log-likelihood = sum[ log( L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ) ]
  27  *
  28  * L(mu1, sigma1, mu2, sigma2, score2 - score1) = sigmoid(mu2 - mu1, sqrt(sigma1² + sigma2²), (score2 - score1))
  29  *
  30  * pdf := 1/(sigma * sqrt(2*Pi)) * exp(-(x - mu)^2 / (2 * sigma^2));
  31  * pdfs := subs({ mu = mu1 - mu2, sigma = sqrt(sigma1^2 + sigma2^2) }, pdf);
  32  * diff(log(pdfs), mu1);
  33  */
  34
  35 struct match {
  36         int player, other_player;
  37         int margin;
  38         float weight;
  39 };
  40 map<int, vector<match> > matches_for_player;
  41 vector<match> all_matches;
  42
  43 void dump_scores(const vector<string> &players, const float *mu, int num_players)
  44 {
  45 #if 0
  46         for (int i = 0; i < num_players; ++i) {
  47                 printf("%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]);
  48         }
  49         printf("\n");
  50 #elif 0
  51         for (int i = 0; i < num_players; ++i) {
  52                 printf("%5.1f ", mu[i]);
  53         }
  54         printf("\n");
  55 #else
  56         for (int i = 0; i < num_players; ++i) {
  57                 printf("%f %s\n", mu[i], players[i].c_str());
  58         }
  59 #endif
  60 }
  61
  62 /*
  63  * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2
  64  * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
  65  *                                sum_i[ (w_i/sigma_c_i)^2 mu1 ] = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ]
  66  *                                mu1 = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (w_i/sigma_c_i)^2 ]
  67  */
  68 void update_mu(float *mu, int player_num, const vector<match> &matches)
  69 {
  70         if (matches.empty()) {
  71                 return;
  72         }
  73
  74         float nom = 0.0f, denom = 0.0f;
  75
  76         // Prior.
  77         {
  78                 float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
  79                 nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
  80                 denom += PRIOR_WEIGHT * inv_sigma2;
  81         }
  82
  83         // All matches.
  84         for (unsigned i = 0; i < matches.size(); ++i) {
  85                 float inv_sigma_c2 = matches[i].weight / (global_sigma * global_sigma);
  86                 float x = matches[i].margin; // / 70.0f;
  87
  88                 nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
  89                 denom += inv_sigma_c2;
  90         }
  91         mu[player_num] = nom / denom;
  92 }
  93
  94 void dump_raw(const float *mu, int num_players)
  95 {
  96         for (unsigned i = 0; i < all_matches.size(); ++i) {
  97                 const match& m = all_matches[i];
  98
  99                 float mu1 = mu[m.player];
 100                 float mu2 = mu[m.other_player];
 101                 float sigma = global_sigma;
 102                 float mu = mu1 - mu2;
 103                 float x = m.margin;
 104                 float w = m.weight;
 105
 106                 printf("%f %f\n", (x - mu) / sigma, w);
 107         }
 108 }
 109
 110 /*
 111  * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
 112  * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
 113  *                                   sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0                            |: sigma != 0
 114  *                                   sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
 115  *                                   sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
 116  */
 117 void update_global_sigma(float *mu, int num_players)
 118 {
 119         float nom = 0.0f, denom = 0.0f;
 120         for (unsigned i = 0; i < all_matches.size(); ++i) {
 121                 const match& m = all_matches[i];
 122
 123                 float mu1 = mu[m.player];
 124                 float mu2 = mu[m.other_player];
 125                 float mu = mu1 - mu2;
 126                 float x = m.margin;
 127                 float w = m.weight;
 128
 129                 nom += w * ((x - mu) * (x - mu));
 130                 denom += w;
 131         }
 132
 133         global_sigma = sqrt(nom / denom);
 134 }
 135
 136 /*
 137  * diff(priorlogL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
 138  * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
 139  *                                   sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0                            |: sigma != 0
 140  *                                   sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
 141  *                                   sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
 142  */
 143 void update_prior_sigma(float *mu, int num_players)
 144 {
 145         float nom = 0.0f, denom = 0.0f;
 146         for (int i = 0; i < num_players; ++i) {
 147                 float mu1 = mu[i];
 148
 149                 nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
 150                 denom += 1.0f;
 151         }
 152
 153         prior_sigma = sqrt(nom / denom);
 154         if (!(prior_sigma > 40.0f)) {
 155                 prior_sigma = 40.0f;
 156         }
 157 }
 158
 159 float compute_logl(float z)
 160 {
 161         return -0.5 * (log(2.0f / M_PI) + z * z);
 162 }
 163
 164 float compute_total_logl(float *mu, int num_players)
 165 {
 166         float total_logl = 0.0f;
 167
 168         // Prior.
 169         for (int i = 0; i < num_players; ++i) {
 170                 total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
 171         }
 172
 173         // Matches.
 174         for (unsigned i = 0; i < all_matches.size(); ++i) {
 175                 const match& m = all_matches[i];
 176
 177                 float mu1 = mu[m.player];
 178                 float mu2 = mu[m.other_player];
 179                 float sigma = global_sigma;
 180                 float mu = mu1 - mu2;
 181                 float x = m.margin;
 182                 float w = m.weight;
 183
 184                 total_logl += w * compute_logl((x - mu) / sigma);
 185         }
 186
 187         return total_logl;
 188 }
 189
 190 /*
 191  * Compute Hessian matrix of the negative log-likelihood, ie. for each term in logL:
 192  *
 193  * M_ij = D_i D_j (- logL) = -w / sigma²                                for i != j
 194  *                            w / sigma²                                for i == j
 195  *
 196  * Note that this does not depend on mu or the margin at all.
 197  */
 198 double hessian[MAX_PLAYERS][MAX_PLAYERS];
 199 void construct_hessian(const float *mu, const float *sigma, int num_players)
 200 {
 201         memset(hessian, 0, sizeof(hessian));
 202
 203         for (unsigned i = 0; i < all_matches.size(); ++i) {
 204                 const match &m = all_matches[i];
 205
 206                 int p1 = m.player;
 207                 int p2 = m.other_player;
 208
 209                 double sigma_sq = global_sigma * global_sigma;
 210                 float w = m.weight;
 211
 212                 hessian[p1][p2] -= w / sigma_sq;
 213                 hessian[p2][p1] -= w / sigma_sq;
 214
 215                 hessian[p1][p1] += w / sigma_sq;
 216                 hessian[p2][p2] += w / sigma_sq;
 217         }
 218
 219         for (int i = 0; i < num_players; ++i) {
 220                 for (int j = 0; j < num_players; ++j) {
 221                         printf("%.12f ", hessian[i][j]);
 222                 }
 223                 printf("\n");
 224         }
 225 }
 226
 227 int main(int argc, char **argv)
 228 {
 229         int num_players;
 230         if (scanf("%d", &num_players) != 1) {
 231                 fprintf(stderr, "Could't read number of players\n");
 232                 exit(1);
 233         }
 234
 235         if (num_players > MAX_PLAYERS) {
 236                 fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS);
 237                 exit(1);
 238         }
 239
 240         vector<string> players;
 241         map<string, int> player_map;
 242
 243         for (int i = 0; i < num_players; ++i) {
 244                 char buf[256];
 245                 if (scanf("%s", buf) != 1) {
 246                         fprintf(stderr, "Couldn't read player %d\n", i);
 247                         exit(1);
 248                 }
 249
 250                 players.push_back(buf);
 251                 player_map[buf] = i;
 252         }
 253
 254         int num_matches = 0;
 255         for ( ;; ) {
 256                 char pl1[256], pl2[256];
 257                 int score1, score2;
 258                 float weight;
 259
 260                 if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) {
 261                         //fprintf(stderr, "Read %d matches.\n", num_matches);
 262                         break;
 263                 }
 264
 265                 ++num_matches;
 266
 267                 if (player_map.count(pl1) == 0) {
 268                         fprintf(stderr, "Unknown player '%s'\n", pl1);
 269                         exit(1);
 270                 }
 271                 if (player_map.count(pl2) == 0) {
 272                         fprintf(stderr, "Unknown player '%s'\n", pl2);
 273                         exit(1);
 274                 }
 275
 276                 match m1;
 277                 m1.player = player_map[pl1];
 278                 m1.other_player = player_map[pl2];
 279                 m1.margin = score1 - score2;
 280                 m1.weight = weight;
 281                 matches_for_player[player_map[pl1]].push_back(m1);
 282
 283                 match m2;
 284                 m2.player = player_map[pl2];
 285                 m2.other_player = player_map[pl1];
 286                 m2.margin = score2 - score1;
 287                 m2.weight = weight;
 288                 matches_for_player[player_map[pl2]].push_back(m2);
 289
 290                 all_matches.push_back(m1);
 291         }
 292
 293         float mu[MAX_PLAYERS];
 294
 295         for (int i = 0; i < num_players; ++i) {
 296                 mu[i] = PRIOR_MU;
 297         }
 298
 299         for (int j = 0; j < 1000; ++j) {
 300                 float old_mu[MAX_PLAYERS];
 301                 float old_global_sigma = global_sigma;
 302                 float old_prior_sigma = prior_sigma;
 303                 memcpy(old_mu, mu, sizeof(mu));
 304                 for (int i = 0; i < num_players; ++i) {
 305                         update_mu(mu, i, matches_for_player[i]);
 306                 }
 307                 update_global_sigma(mu, num_players);
 308                 update_prior_sigma(mu, num_players);
 309                 /* for (int i = 0; i < num_players; ++i) {
 310                         update_sigma(mu, sigma, i, matches_for_player[i]);
 311                         dump_scores(players, mu, sigma, num_players);
 312                 } */
 313
 314                 float sumdiff = 0.0f;
 315                 for (int i = 0; i < num_players; ++i) {
 316                         sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
 317                 }
 318                 sumdiff += (prior_sigma - old_prior_sigma) * (prior_sigma - old_prior_sigma);
 319                 sumdiff += (global_sigma - old_global_sigma) * (global_sigma - old_global_sigma);
 320                 if (sumdiff < EPSILON) {
 321                         //fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
 322                         printf("%d -1\n", j + 1);
 323                         break;
 324                 }
 325         }
 326
 327 #if DUMP_RAW
 328         dump_raw(mu, num_players);
 329 #else
 330         dump_scores(players, mu, num_players);
 331         //fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
 332         printf("%f -2\n", global_sigma / sqrt(2.0f));
 333         printf("%f -3\n", prior_sigma);
 334
 335         float total_logl = compute_total_logl(mu, num_players);
 336         printf("%f -4\n", total_logl);
 337
 338 //      construct_hessian(mu, sigma, num_players);
 339 #endif
 340 }