Convert index.pl to templates.

[wloh] / bayeswf.cpp

diff --git a/bayeswf.cpp b/bayeswf.cpp
index 6214d281fd24ff83571cf0e272203af5c328456b..cc591832a556f023ea6aa364a44ce5422a0b0740 100644 (file)
--- 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,13 +11,15 @@
 #include <algorithm>
 
 using namespace std;
+using namespace Eigen;
 
-#define PRIOR_MU 1500
+#define PRIOR_MU 500
 #define PRIOR_WEIGHT 1.0
 #define MAX_PLAYERS 4096
 #define DUMP_RAW 0
 
 float mu[MAX_PLAYERS];
+float mu_stddev[MAX_PLAYERS];
 float global_sigma = 70.0f;
 float prior_sigma = 70.0f;
 
@@ -40,7 +44,7 @@ struct match {
 map<int, vector<match> > matches_for_player;
 vector<match> all_matches;
 
-void dump_scores(const vector<string> &players, const float *mu, int num_players)
+void dump_scores(const vector<string> &players, const float *mu, const float *mu_stddev, int num_players)
 {
 #if 0
        for (int i = 0; i < num_players; ++i) {
@@ -54,7 +58,7 @@ void dump_scores(const vector<string> &players, const float *mu, int num_players
        printf("\n");
 #else
        for (int i = 0; i < num_players; ++i) {
-               printf("%f %s\n", mu[i], players[i].c_str());
+               printf("%f %f %s\n", mu[i], mu_stddev[i], players[i].c_str());
        }
 #endif
 }
@@ -195,11 +199,15 @@ float compute_total_logl(float *mu, int num_players)
  *
  * 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) {
+               hessian(i, i) += 1.0f / (prior_sigma * prior_sigma);
+       }
        for (unsigned i = 0; i < all_matches.size(); ++i) {
                const match &m = all_matches[i];
 
@@ -209,18 +217,31 @@ void construct_hessian(const float *mu, const float *sigma, int num_players)
                double sigma_sq = global_sigma * global_sigma;
                float w = m.weight;
 
-               hessian[p1][p2] -= w / sigma_sq;
-               hessian[p2][p1] -= w / sigma_sq;
+               hessian(p1, p2) -= w / sigma_sq;
+               hessian(p2, p1) -= w / sigma_sq;
 
-               hessian[p1][p1] += w / sigma_sq;
-               hessian[p2][p2] += 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<string> &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 (unsigned i = 0; i < players.size(); ++i) {
+               for (unsigned j = 0; j < players.size(); ++j) {
+                       printf("covariance %s %s %f\n",
+                              players[i].c_str(),
+                              players[j].c_str(),
+                              ih(i, j));
                }
-               printf("\n");
        }
 }
 
@@ -319,7 +340,7 @@ int main(int argc, char **argv)
                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);
+                       printf("aux_param num_iterations %d\n", j + 1);
                        break;
                }
        }
@@ -327,14 +348,14 @@ int main(int argc, char **argv)
 #if DUMP_RAW
        dump_raw(mu, num_players);
 #else
-       dump_scores(players, mu, num_players);
+       construct_hessian(mu, num_players);
+       compute_mu_uncertainty(mu, players);
+       dump_scores(players, mu, mu_stddev, num_players);
        //fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
-       printf("%f -2\n", global_sigma / sqrt(2.0f));
-       printf("%f -3\n", prior_sigma);
+       printf("aux_param score_stddev %f\n", global_sigma / sqrt(2.0f));
+       printf("aux_param rating_prior_stddev %f\n", prior_sigma);
 
        float total_logl = compute_total_logl(mu, num_players);
-       printf("%f -4\n", total_logl);
-
-//     construct_hessian(mu, sigma, num_players);
+       printf("aux_param total_log_likelihood %f\n", total_logl);
 #endif
 }