X-Git-Url: https://git.sesse.net/?p=wloh;a=blobdiff_plain;f=bayeswf.cpp;h=eee6882e4e18f59dcc9c2b8f37f950e22b3d6e13;hp=1a4893050f4d3520479785f6ccec9cfa0819bad0;hb=fde909c294de9806dd6337f5acb0ed87c41557c6;hpb=e6ad3fff4d7f562c9c68948b340777c90d4867c0

diff --git a/bayeswf.cpp b/bayeswf.cpp
index 1a48930..eee6882 100644
--- a/bayeswf.cpp
+++ b/bayeswf.cpp
@@ -2,6 +2,8 @@
 #include <math.h>
 #include <string.h>
 #include <stdlib.h>
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
 
 #include <map>
 #include <vector>
@@ -9,11 +11,38 @@
 #include <algorithm>
 
 using namespace std;
+using namespace Eigen;
 
+#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
 
@@ -29,43 +58,51 @@ float sigma[MAX_PLAYERS];
  */
 
 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 1
+#if USE_DB
 	for (int i = 0; i < num_players; ++i) {
-		fprintf(stderr, "%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);
 	}
-	fprintf(stderr, "\n");
 #else
 	for (int i = 0; i < num_players; ++i) {
-		fprintf(stderr, "%5.1f ", mu[i]);
+		printf("%f %f %d\n", mu[i], mu_stddev[i], players[i]);
 	}
-	fprintf(stderr, "\n");
 #endif
 }
 
 /*
- * diff(logL, mu1) = -(mu1 - mu2 - x) / sigma_c^2
- * maximizer for mu1 is given by: sum_i[ (1/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
- *                                sum_i[ (1/sigma_c_i)^2 mu1 ] = sum_i [ (1/sigma_c_i)^2 ( mu2_i + x_i ) ]
- *                                mu1 = sum_i [ (1/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (1/sigma_c_i)^2 ]
+ * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2
+ * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
+ *                                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;
 	}
 
 	float nom = 0.0f, denom = 0.0f;
+
+	// Prior.
+	{
+		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 = 1.0f / (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;
@@ -74,97 +111,189 @@ void update_mu(float *mu, float *sigma, int player_num, const vector<match> &mat
 	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);
 }
 
-void renormalize(float *mu, float *sigma, int num_players)
+/*
+ * 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_prior_sigma(float *mu, int num_players)
 {
-	float avg = 0.0f;
+	float nom = 0.0f, denom = 0.0f;
 	for (int i = 0; i < num_players; ++i) {
-		avg += mu[i];
+		float mu1 = mu[i];
+
+		nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
+		denom += 1.0f;
 	}
-	float corr = 1500.0f - avg / num_players;
+
+	prior_sigma = sqrt(nom / denom);
+	if (!(prior_sigma > 40.0f)) {
+		prior_sigma = 40.0f;
+	}
+}
+
+float compute_logl(float z)
+{
+	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) {
-		mu[i] += corr;
+		total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
 	}
+
+	// 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;
 }
 
 /*
- * Compute Fisher information matrix of the log-likelihood, evaluated at the MLE,
-c
- * ie. M_ij = E[ (D_i logL) (D_j logL) ] = - sum( ( x - (mu_1 - mu_2) )² / sigma_c⁴ )  for i != j
- *                                       = - sum( 1 / sigma_c² )                     for i == j
+ * Compute Hessian matrix of the negative log-likelihood, ie. for each term in logL:
  *
- * The Hessian matrix is generally zero and thus not as interesting.
+ * M_ij = D_i D_j (- logL) = -w / sigma²                                for i != j
+ *                            w / sigma²                                for i == j
+ *
+ * Note that this does not depend on mu or the margin at all.
  */
-void construct_fim(const float *mu, const float *sigma, int num_players)
+Matrix<float, Dynamic, Dynamic> hessian;
+void construct_hessian(const float *mu, int num_players)
 {
-	float fim[MAX_PLAYERS][MAX_PLAYERS];
-	memset(fim, 0, sizeof(fim));
+	hessian = Matrix<float, Dynamic, Dynamic>(num_players, num_players);
+	hessian.fill(0.0f);
 
 	for (int i = 0; i < num_players; ++i) {
-		float mu1 = mu[i];
-		float 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;
-			float mu2 = mu[j];
-			float sigma2 = sigma[j];
+		int p1 = m.player;
+		int p2 = m.other_player;
 
-			float x = matches_for_player[i][k].margin;
-			float sigma_sq = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
+		double sigma_sq = global_sigma * global_sigma;
+		float w = m.weight;
 
-			fprintf(stderr, "exp_diff_sq=%f  sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * sigma_sq);
+		hessian(p1, p2) -= w / sigma_sq;
+		hessian(p2, p1) -= w / sigma_sq;
 
-#if 1
-			fim[i][i] += (x - (mu1 - mu2)) * (x - (mu1 - mu2)) / (sigma_sq * sigma_sq);
-			fim[i][j] -= (x - (mu1 - mu2)) * (x - (mu1 - mu2)) / (sigma_sq * sigma_sq);
-#else
-			fim[i][i] -= 1.0f / sigma_sq;
-			fim[i][j] += 1.0f / sigma_sq;
-#endif
-		}
+		hessian(p1, p1) += w / sigma_sq;
+		hessian(p2, p2) += w / sigma_sq;
+	}
+}
+
+// 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));
+	}
 
-		for (int j = 0; j < num_players; ++j) {
-			printf("%f ", fim[i][j]);
+#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);
 		}
-		printf("\n");
 	}
+#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));
+		}
+	}
+#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);
 	}
@@ -174,90 +303,190 @@ int main(int argc, char **argv)
 		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", pl1, pl2, &score1, &score2) != 4) {
-			fprintf(stderr, "Read %d matches.\n", num_matches);
+		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);
-
-	dump_scores(players, mu, sigma, num_players);
 
-	for (int j = 0; j < 100; ++j) {
+	int num_iterations = -1;
+	for (int j = 0; j < 1000; ++j) {
 		float old_mu[MAX_PLAYERS];
-		float old_sigma[MAX_PLAYERS];
-		memcpy(old_mu, mu, sizeof(float) * MAX_PLAYERS);
-		memcpy(old_sigma, sigma, sizeof(float) * MAX_PLAYERS);
+		float old_global_sigma = global_sigma;
+		float old_prior_sigma = prior_sigma;
+		memcpy(old_mu, mu, sizeof(mu));
 		for (int i = 0; i < num_players; ++i) {
-			update_mu(mu, sigma, i, matches_for_player[i]);
-			renormalize(mu, sigma, num_players);
-			dump_scores(players, mu, sigma, num_players);
+			update_mu(mu, i, matches_for_player[i]);
 		}
+		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);
 		} */
-		bool any_difference = false;
+
+		float sumdiff = 0.0f;
 		for (int i = 0; i < num_players; ++i) {
-			if (fabs(mu[i] - old_mu[i]) > EPSILON ||
-			    fabs(sigma[i] - old_sigma[i]) > EPSILON) {
-				any_difference = true;
-				break;
-			}
+			sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
 		}
-		if (!any_difference) {
-			fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
+		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);
+			num_iterations = j + 1;
 			break;
 		}
 	}
 
-//	construct_fim(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));
+	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
+
+	txn.commit();
 }