#include <math.h>
#include <string.h>
#include <stdlib.h>
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
#include <map>
#include <vector>
#include <algorithm>
using namespace std;
+using namespace Eigen;
-#define PRIOR_MU 1500
-#define PRIOR_SIGMA 100
+#define PRIOR_MU 500
+#define PRIOR_WEIGHT 1.0
#define MAX_PLAYERS 4096
+#define DUMP_RAW 0
+#define USE_DB 1
+
+#if USE_DB
+#include <pqxx/connection>
+#include <pqxx/tablewriter>
+#include <pqxx/transaction>
+#endif
float mu[MAX_PLAYERS];
-float sigma[MAX_PLAYERS];
+float mu_stddev[MAX_PLAYERS];
+float global_sigma;
+float prior_sigma;
+
+// Data waiting for insertion into the database.
+
+struct RatingDBTuple {
+ int player;
+ float mu, mu_stddev;
+};
+struct CovarianceDBTuple {
+ int player1, player2;
+ float covariance;
+};
+vector<RatingDBTuple> rating_db_tuples;
+vector<CovarianceDBTuple> covariance_db_tuples;
+map<pair<string, string>, float> aux_params;
#define EPSILON 1e-3
*/
struct match {
- int other_player;
+ int player, other_player;
int margin;
float weight;
};
map<int, vector<match> > matches_for_player;
+vector<match> all_matches;
-void dump_scores(const vector<string> &players, const float *mu, const float *sigma, int num_players)
+void dump_scores(const vector<int> &players, const float *mu, const float *mu_stddev, int num_players)
{
-#if 0
+#if USE_DB
for (int i = 0; i < num_players; ++i) {
- printf("%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]);
+ RatingDBTuple tuple = { players[i], mu[i], mu_stddev[i] };
+ rating_db_tuples.push_back(tuple);
}
- printf("\n");
-#elif 0
- for (int i = 0; i < num_players; ++i) {
- printf("%5.1f ", mu[i]);
- }
- printf("\n");
#else
for (int i = 0; i < num_players; ++i) {
- printf("%f %s\n", mu[i], players[i].c_str());
+ printf("%f %f %d\n", mu[i], mu_stddev[i], players[i]);
}
#endif
}
* sum_i[ (w_i/sigma_c_i)^2 mu1 ] = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ]
* mu1 = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (w_i/sigma_c_i)^2 ]
*/
-void update_mu(float *mu, float *sigma, int player_num, const vector<match> &matches)
+void update_mu(float *mu, int player_num, const vector<match> &matches)
{
if (matches.empty()) {
return;
// Prior.
{
- float inv_sigma2 = 1.0f / (PRIOR_SIGMA * PRIOR_SIGMA);
- nom += PRIOR_MU * inv_sigma2;
- denom += inv_sigma2;
+ float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
+ nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
+ denom += PRIOR_WEIGHT * inv_sigma2;
}
// All matches.
for (unsigned i = 0; i < matches.size(); ++i) {
- float sigma1 = sigma[player_num];
- float sigma2 = sigma[matches[i].other_player];
- float inv_sigma_c2 = matches[i].weight / (sigma1 * sigma1 + sigma2 * sigma2);
+ float inv_sigma_c2 = matches[i].weight / (global_sigma * global_sigma);
float x = matches[i].margin; // / 70.0f;
nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
mu[player_num] = nom / denom;
}
+void dump_raw(const float *mu, int num_players)
+{
+ for (unsigned i = 0; i < all_matches.size(); ++i) {
+ const match& m = all_matches[i];
+
+ float mu1 = mu[m.player];
+ float mu2 = mu[m.other_player];
+ float sigma = global_sigma;
+ float mu = mu1 - mu2;
+ float x = m.margin;
+ float w = m.weight;
+
+ printf("%f %f\n", (x - mu) / sigma, w);
+ }
+}
+
/*
- * diff(logL, sigma1) = sigma1 (-sigma1² - sigma2² + (x - mu)²) / sigma_c²
- * maximizer for sigma1 is given by: sum_i[ (1/sigma_c_i)² sigma1 ((x - mu)² - (sigma1² + sigma2²) ] = 0
- * sum_i[ (x - mu)² - sigma1² - sigma2² ] = 0 |: sigma1 != 0, sigma2 != 0
- * sum_i[ (x - mu)² - sigma2² ] = sum[ sigma1² ]
- * sigma1 = sqrt( sum_i[ (x - mu)² - sigma2² ] / N )
+ * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
+ * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
+ * sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0 |: sigma != 0
+ * sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
+ * sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
*/
-void update_sigma(float *mu, float *sigma, int player_num, const vector<match> &matches)
+void update_global_sigma(float *mu, int num_players)
{
- if (matches.size() < 2) {
- return;
- }
+ float nom = 0.0f, denom = 0.0f;
+ for (unsigned i = 0; i < all_matches.size(); ++i) {
+ const match& m = all_matches[i];
- float sum = 0.0f;
- for (unsigned i = 0; i < matches.size(); ++i) {
- float mu1 = mu[player_num];
- float mu2 = mu[matches[i].other_player];
+ float mu1 = mu[m.player];
+ float mu2 = mu[m.other_player];
float mu = mu1 - mu2;
- float sigma2 = sigma[matches[i].other_player];
- float x = matches[i].margin;
+ float x = m.margin;
+ float w = m.weight;
- //fprintf(stderr, "x=%f mu=%f sigma2=%f add %f-%f = %f\n", x, mu, sigma2, (x-mu)*(x-mu), sigma2*sigma2, (x - mu) * (x - mu) - sigma2 * sigma2);
- sum += (x - mu) * (x - mu) - sigma2 * sigma2;
+ nom += w * ((x - mu) * (x - mu));
+ denom += w;
}
- if (sum <= 0) {
- return;
- }
- //fprintf(stderr, "sum=%f\n", sum);
- sigma[player_num] = sqrt(sum / matches.size());
+ global_sigma = sqrt(nom / denom);
}
/*
- * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
+ * diff(priorlogL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
* maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
* sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0 |: sigma != 0
* sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
* sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
*/
-void update_global_sigma(float *mu, float *sigma, int num_players)
+void update_prior_sigma(float *mu, int num_players)
{
float nom = 0.0f, denom = 0.0f;
for (int i = 0; i < num_players; ++i) {
- for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
- const match& m = matches_for_player[i][j];
-
- // Only count each match once.
- if (m.other_player <= i) {
- continue;
- }
-
- float mu1 = mu[i];
- float mu2 = mu[m.other_player];
- float mu = mu1 - mu2;
- float x = m.margin;
- float w = m.weight;
-
- nom += w * ((x - mu) * (x - mu));
- denom += w;
- }
+ float mu1 = mu[i];
+
+ nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
+ denom += 1.0f;
}
- float best_sigma = sqrt(nom / denom) / sqrt(2.0f); // Divide evenly between the two players.
- for (int i = 0; i < num_players; ++i) {
- sigma[i] = best_sigma;
+ prior_sigma = sqrt(nom / denom);
+ if (!(prior_sigma > 40.0f)) {
+ prior_sigma = 40.0f;
}
}
-void renormalize(float *mu, float *sigma, int num_players)
+float compute_logl(float z)
{
- float avg = 0.0f;
+ return -0.5 * (log(2.0f / M_PI) + z * z);
+}
+
+float compute_total_logl(float *mu, int num_players)
+{
+ float total_logl = 0.0f;
+
+ // Prior.
for (int i = 0; i < num_players; ++i) {
- avg += mu[i];
+ total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
}
- float corr = 1500.0f - avg / num_players;
- for (int i = 0; i < num_players; ++i) {
- mu[i] += corr;
+
+ // Matches.
+ for (unsigned i = 0; i < all_matches.size(); ++i) {
+ const match& m = all_matches[i];
+
+ float mu1 = mu[m.player];
+ float mu2 = mu[m.other_player];
+ float sigma = global_sigma;
+ float mu = mu1 - mu2;
+ float x = m.margin;
+ float w = m.weight;
+
+ total_logl += w * compute_logl((x - mu) / sigma);
}
+
+ return total_logl;
}
/*
*
* Note that this does not depend on mu or the margin at all.
*/
-double hessian[MAX_PLAYERS][MAX_PLAYERS];
-void construct_hessian(const float *mu, const float *sigma, int num_players)
+Matrix<float, Dynamic, Dynamic> hessian;
+void construct_hessian(const float *mu, int num_players)
{
- memset(hessian, 0, sizeof(hessian));
+ hessian = Matrix<float, Dynamic, Dynamic>(num_players, num_players);
+ hessian.fill(0.0f);
for (int i = 0; i < num_players; ++i) {
- double sigma1 = sigma[i];
+ hessian(i, i) += 1.0f / (prior_sigma * prior_sigma);
+ }
+ for (unsigned i = 0; i < all_matches.size(); ++i) {
+ const match &m = all_matches[i];
- for (unsigned k = 0; k < matches_for_player[i].size(); ++k) {
- int j = matches_for_player[i][k].other_player;
+ int p1 = m.player;
+ int p2 = m.other_player;
- double sigma2 = sigma[j];
- double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2;
+ double sigma_sq = global_sigma * global_sigma;
+ float w = m.weight;
- float w = matches_for_player[i][k].weight;
+ hessian(p1, p2) -= w / sigma_sq;
+ hessian(p2, p1) -= w / sigma_sq;
- hessian[i][j] -= w / sigma_sq;
- hessian[i][i] += w / sigma_sq;
- }
+ hessian(p1, p1) += w / sigma_sq;
+ hessian(p2, p2) += w / sigma_sq;
}
+}
- for (int i = 0; i < num_players; ++i) {
- for (int j = 0; j < num_players; ++j) {
- printf("%.12f ", hessian[i][j]);
+// Compute uncertainty (stddev) of mu estimates, which is sqrt((H^-1)_ii),
+// where H is the Hessian (see construct_hessian()).
+void compute_mu_uncertainty(const float *mu, const vector<int> &players)
+{
+ // FIXME: Use pseudoinverse if applicable.
+ Matrix<float, Dynamic, Dynamic> ih = hessian.inverse();
+ for (unsigned i = 0; i < players.size(); ++i) {
+ mu_stddev[i] = sqrt(ih(i, i));
+ }
+
+#if USE_DB
+ for (unsigned i = 0; i < players.size(); ++i) {
+ for (unsigned j = 0; j < players.size(); ++j) {
+ CovarianceDBTuple tuple;
+ tuple.player1 = players[i];
+ tuple.player2 = players[j];
+ tuple.covariance = ih(i, j);
+ covariance_db_tuples.push_back(tuple);
+ }
+ }
+#else
+ for (unsigned i = 0; i < players.size(); ++i) {
+ for (unsigned j = 0; j < players.size(); ++j) {
+ printf("covariance %d %d %f\n",
+ players[i],
+ players[j],
+ ih(i, j));
}
- printf("\n");
}
+#endif
}
-int main(int argc, char **argv)
+void process_file(const char *filename)
{
+ global_sigma = 70.0f;
+ prior_sigma = 70.0f;
+ matches_for_player.clear();
+ all_matches.clear();
+
+ FILE *fp = fopen(filename, "r");
+ if (fp == NULL) {
+ perror(filename);
+ exit(1);
+ }
+
+ char locale[256];
+ if (fscanf(fp, "%s", locale) != 1) {
+ fprintf(stderr, "Could't read number of players\n");
+ exit(1);
+ }
+
int num_players;
- if (scanf("%d", &num_players) != 1) {
+ if (fscanf(fp,"%d", &num_players) != 1) {
fprintf(stderr, "Could't read number of players\n");
exit(1);
}
exit(1);
}
- vector<string> players;
- map<string, int> player_map;
+ vector<int> players;
+ map<int, int> player_map;
for (int i = 0; i < num_players; ++i) {
char buf[256];
- if (scanf("%s", buf) != 1) {
+ if (fscanf(fp, "%s", buf) != 1) {
fprintf(stderr, "Couldn't read player %d\n", i);
exit(1);
}
- players.push_back(buf);
- player_map[buf] = i;
+ players.push_back(atoi(buf));
+ player_map[atoi(buf)] = i;
}
int num_matches = 0;
for ( ;; ) {
- char pl1[256], pl2[256];
+ int pl1, pl2;
int score1, score2;
float weight;
- if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) {
+ if (fscanf(fp, "%d %d %d %d %f", &pl1, &pl2, &score1, &score2, &weight) != 5) {
//fprintf(stderr, "Read %d matches.\n", num_matches);
break;
}
++num_matches;
if (player_map.count(pl1) == 0) {
- fprintf(stderr, "Unknown player '%s'\n", pl1);
+ fprintf(stderr, "Unknown player '%d'\n", pl1);
exit(1);
}
if (player_map.count(pl2) == 0) {
- fprintf(stderr, "Unknown player '%s'\n", pl2);
+ fprintf(stderr, "Unknown player '%d'\n", pl2);
exit(1);
}
match m1;
+ m1.player = player_map[pl1];
m1.other_player = player_map[pl2];
m1.margin = score1 - score2;
m1.weight = weight;
matches_for_player[player_map[pl1]].push_back(m1);
match m2;
+ m2.player = player_map[pl2];
m2.other_player = player_map[pl1];
m2.margin = score2 - score1;
m2.weight = weight;
matches_for_player[player_map[pl2]].push_back(m2);
+
+ all_matches.push_back(m1);
}
+
+ fclose(fp);
float mu[MAX_PLAYERS];
- float sigma[MAX_PLAYERS];
for (int i = 0; i < num_players; ++i) {
- mu[i] = 1500.0f;
- sigma[i] = 70.0f / sqrt(2.0f);
+ mu[i] = PRIOR_MU;
}
- renormalize(mu, sigma, num_players);
+ int num_iterations = -1;
for (int j = 0; j < 1000; ++j) {
float old_mu[MAX_PLAYERS];
- float old_sigma[MAX_PLAYERS];
+ float old_global_sigma = global_sigma;
+ float old_prior_sigma = prior_sigma;
memcpy(old_mu, mu, sizeof(mu));
- memcpy(old_sigma, sigma, sizeof(sigma));
for (int i = 0; i < num_players; ++i) {
- update_mu(mu, sigma, i, matches_for_player[i]);
- renormalize(mu, sigma, num_players);
+ update_mu(mu, i, matches_for_player[i]);
}
- update_global_sigma(mu, sigma, num_players);
+ update_global_sigma(mu, num_players);
+ update_prior_sigma(mu, num_players);
/* for (int i = 0; i < num_players; ++i) {
update_sigma(mu, sigma, i, matches_for_player[i]);
dump_scores(players, mu, sigma, num_players);
float sumdiff = 0.0f;
for (int i = 0; i < num_players; ++i) {
sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
- sumdiff += (sigma[i] - old_sigma[i]) * (sigma[i] - old_sigma[i]);
}
+ sumdiff += (prior_sigma - old_prior_sigma) * (prior_sigma - old_prior_sigma);
+ sumdiff += (global_sigma - old_global_sigma) * (global_sigma - old_global_sigma);
if (sumdiff < EPSILON) {
//fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
- printf("%d -1\n", j + 1);
+ num_iterations = j + 1;
break;
}
}
- dump_scores(players, mu, sigma, num_players);
+
+ construct_hessian(mu, num_players);
+ aux_params[make_pair(locale, "num_iterations")] = num_iterations;
+ aux_params[make_pair(locale, "score_stddev")] = global_sigma / sqrt(2.0f);
+ aux_params[make_pair(locale, "rating_prior_stddev")] = prior_sigma;
+ aux_params[make_pair(locale, "total_log_likelihood")] = compute_total_logl(mu, num_players);
+
+ compute_mu_uncertainty(mu, players);
+ dump_scores(players, mu, mu_stddev, num_players);
+}
+
+int main(int argc, char **argv)
+{
+#if USE_DB
+ pqxx::connection conn("dbname=wloh host=127.0.0.1 user=wloh password=oto4iCh5");
+#endif
+
+ for (int i = 1; i < argc; ++i) {
+ process_file(argv[i]);
+ }
+
+#if DUMP_RAW
+ dump_raw(mu, num_players);
+ return 0;
+#endif
+
+#if USE_DB
+ pqxx::work txn(conn);
+ txn.exec("SET client_min_messages TO WARNING");
+
+ // Dump aux_params.
+ {
+ txn.exec("TRUNCATE aux_params");
+ pqxx::tablewriter writer(txn, "aux_params");
+ for (map<pair<string, string>, float>::const_iterator it = aux_params.begin(); it != aux_params.end(); ++it) {
+ char str[128];
+ snprintf(str, 128, "%f", it->second);
+
+ vector<string> tuple;
+ tuple.push_back(it->first.first); // locale
+ tuple.push_back(it->first.second); // parameter name
+ tuple.push_back(str);
+ writer.push_back(tuple);
+ }
+ writer.complete();
+ }
+
+ // Dump ratings.
+ {
+ txn.exec("TRUNCATE ratings");
+ pqxx::tablewriter writer(txn, "ratings");
+ for (unsigned i = 0; i < rating_db_tuples.size(); ++i) {
+ char player_str[128], mu_str[128], mu_stddev_str[128];
+ snprintf(player_str, 128, "%d", rating_db_tuples[i].player);
+ snprintf(mu_str, 128, "%f", rating_db_tuples[i].mu);
+ snprintf(mu_stddev_str, 128, "%f", rating_db_tuples[i].mu_stddev);
+
+ vector<string> tuple;
+ tuple.push_back(player_str);
+ tuple.push_back(mu_str);
+ tuple.push_back(mu_stddev_str);
+ writer.push_back(tuple);
+ }
+ writer.complete();
+ }
+
+ // Create a table new_covariance, and dump covariance into it.
+ {
+ txn.exec("CREATE TABLE new_covariance ( player1 smallint NOT NULL, player2 smallint NOT NULL, cov float NOT NULL )");
+ pqxx::tablewriter writer(txn, "new_covariance");
+ for (unsigned i = 0; i < covariance_db_tuples.size(); ++i) {
+ char player1_str[128], player2_str[128], cov_str[128];
+ snprintf(player1_str, 128, "%d", covariance_db_tuples[i].player1);
+ snprintf(player2_str, 128, "%d", covariance_db_tuples[i].player2);
+ snprintf(cov_str, 128, "%f", covariance_db_tuples[i].covariance);
+
+ vector<string> tuple;
+ tuple.push_back(player1_str);
+ tuple.push_back(player2_str);
+ tuple.push_back(cov_str);
+ writer.push_back(tuple);
+ }
+ writer.complete();
+ }
+
+ // Create index, and rename new_covariance on top of covariance.
+ txn.exec("ALTER TABLE new_covariance ADD PRIMARY KEY ( player1, player2 );");
+ txn.exec("DROP TABLE IF EXISTS covariance");
+ txn.exec("ALTER TABLE new_covariance RENAME TO covariance");
+#else
//fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
- printf("%f -2\n", sigma[0]);
+ for (map<pair<string, string>, float>::const_iterator it = aux_params.begin(); it != aux_params.end(); ++it) {
+ printf("%s: aux_param %s %f\n", it->first.first.c_str(), it->first.second.c_str(), it->second);
+ }
+#endif
-// construct_hessian(mu, sigma, num_players);
+ txn.commit();
}
projects / wloh / blobdiff