6 #include <Eigen/Eigenvalues>
14 using namespace Eigen;
17 #define PRIOR_WEIGHT 1.0
18 #define MAX_PLAYERS 4096
23 #include <pqxx/connection>
24 #include <pqxx/tablewriter>
25 #include <pqxx/transaction>
28 float mu[MAX_PLAYERS];
29 float mu_stddev[MAX_PLAYERS];
33 // Data waiting for insertion into the database.
35 struct RatingDBTuple {
39 struct CovarianceDBTuple {
43 vector<RatingDBTuple> rating_db_tuples;
44 vector<CovarianceDBTuple> covariance_db_tuples;
45 map<pair<string, string>, float> aux_params;
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) ) ]
53 * L(mu1, sigma1, mu2, sigma2, score2 - score1) = sigmoid(mu2 - mu1, sqrt(sigma1² + sigma2²), (score2 - score1))
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);
61 int player, other_player;
65 map<int, vector<match> > matches_for_player;
66 vector<match> all_matches;
68 void dump_scores(const vector<int> &players, const float *mu, const float *mu_stddev, int num_players)
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);
76 for (int i = 0; i < num_players; ++i) {
77 printf("%f %f %d\n", mu[i], mu_stddev[i], players[i]);
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 ]
88 void update_mu(float *mu, int player_num, const vector<match> &matches)
90 if (matches.empty()) {
94 float nom = 0.0f, denom = 0.0f;
98 float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
99 nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
100 denom += PRIOR_WEIGHT * inv_sigma2;
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;
108 nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
109 denom += inv_sigma_c2;
111 mu[player_num] = nom / denom;
114 void dump_raw(const float *mu, int num_players)
116 for (unsigned i = 0; i < all_matches.size(); ++i) {
117 const match& m = all_matches[i];
119 float mu1 = mu[m.player];
120 float mu2 = mu[m.other_player];
121 float sigma = global_sigma;
122 float mu = mu1 - mu2;
126 printf("%f %f\n", (x - mu) / sigma, w);
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] )
137 void update_global_sigma(float *mu, int num_players)
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];
143 float mu1 = mu[m.player];
144 float mu2 = mu[m.other_player];
145 float mu = mu1 - mu2;
149 nom += w * ((x - mu) * (x - mu));
153 global_sigma = sqrt(nom / denom);
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] )
163 void update_prior_sigma(float *mu, int num_players)
165 float nom = 0.0f, denom = 0.0f;
166 for (int i = 0; i < num_players; ++i) {
169 nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
173 prior_sigma = sqrt(nom / denom);
174 if (!(prior_sigma > 40.0f)) {
179 float compute_logl(float z)
181 return -0.5 * (log(2.0f / M_PI) + z * z);
184 float compute_total_logl(float *mu, int num_players)
186 float total_logl = 0.0f;
189 for (int i = 0; i < num_players; ++i) {
190 total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
194 for (unsigned i = 0; i < all_matches.size(); ++i) {
195 const match& m = all_matches[i];
197 float mu1 = mu[m.player];
198 float mu2 = mu[m.other_player];
199 float sigma = global_sigma;
200 float mu = mu1 - mu2;
204 total_logl += w * compute_logl((x - mu) / sigma);
211 * Compute Hessian matrix of the negative log-likelihood, ie. for each term in logL:
213 * M_ij = D_i D_j (- logL) = -w / sigma² for i != j
214 * w / sigma² for i == j
216 * Note that this does not depend on mu or the margin at all.
218 Matrix<float, Dynamic, Dynamic> hessian;
219 void construct_hessian(const float *mu, int num_players)
221 hessian = Matrix<float, Dynamic, Dynamic>(num_players, num_players);
224 for (int i = 0; i < num_players; ++i) {
225 hessian(i, i) += 1.0f / (prior_sigma * prior_sigma);
227 for (unsigned i = 0; i < all_matches.size(); ++i) {
228 const match &m = all_matches[i];
231 int p2 = m.other_player;
233 double sigma_sq = global_sigma * global_sigma;
236 hessian(p1, p2) -= w / sigma_sq;
237 hessian(p2, p1) -= w / sigma_sq;
239 hessian(p1, p1) += w / sigma_sq;
240 hessian(p2, p2) += w / sigma_sq;
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)
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));
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);
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",
276 void process_file(const char *filename)
278 global_sigma = 70.0f;
280 matches_for_player.clear();
283 FILE *fp = fopen(filename, "r");
290 if (fscanf(fp, "%s", locale) != 1) {
291 fprintf(stderr, "Could't read number of players\n");
296 if (fscanf(fp,"%d", &num_players) != 1) {
297 fprintf(stderr, "Could't read number of players\n");
301 if (num_players > MAX_PLAYERS) {
302 fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS);
307 map<int, int> player_map;
309 for (int i = 0; i < num_players; ++i) {
311 if (fscanf(fp, "%s", buf) != 1) {
312 fprintf(stderr, "Couldn't read player %d\n", i);
316 players.push_back(atoi(buf));
317 player_map[atoi(buf)] = i;
326 if (fscanf(fp, "%d %d %d %d %f", &pl1, &pl2, &score1, &score2, &weight) != 5) {
327 //fprintf(stderr, "Read %d matches.\n", num_matches);
333 if (player_map.count(pl1) == 0) {
334 fprintf(stderr, "Unknown player '%d'\n", pl1);
337 if (player_map.count(pl2) == 0) {
338 fprintf(stderr, "Unknown player '%d'\n", pl2);
343 m1.player = player_map[pl1];
344 m1.other_player = player_map[pl2];
345 m1.margin = score1 - score2;
347 matches_for_player[player_map[pl1]].push_back(m1);
350 m2.player = player_map[pl2];
351 m2.other_player = player_map[pl1];
352 m2.margin = score2 - score1;
354 matches_for_player[player_map[pl2]].push_back(m2);
356 all_matches.push_back(m1);
361 float mu[MAX_PLAYERS];
363 for (int i = 0; i < num_players; ++i) {
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]);
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);
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]);
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;
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);
402 compute_mu_uncertainty(mu, players);
403 dump_scores(players, mu, mu_stddev, num_players);
406 int main(int argc, char **argv)
409 pqxx::connection conn("dbname=wloh host=127.0.0.1 user=wloh password=oto4iCh5");
412 for (int i = 1; i < argc; ++i) {
413 process_file(argv[i]);
417 dump_raw(mu, num_players);
422 pqxx::work txn(conn);
423 txn.exec("SET client_min_messages TO WARNING");
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) {
431 snprintf(str, 128, "%f", it->second);
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);
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);
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);
461 // Create a table new_covariance, and dump covariance into it.
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);
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);
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");
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);