#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_WEIGHT 1.0
*
* Note that this does not depend on mu or the margin at all.
*/
-double hessian[MAX_PLAYERS][MAX_PLAYERS];
+double hessian[MAX_PLAYERS * MAX_PLAYERS];
void construct_hessian(const float *mu, int num_players)
{
memset(hessian, 0, sizeof(hessian));
for (int i = 0; i < num_players; ++i) {
- hessian[i][i] += 1.0f / (prior_sigma * prior_sigma);
+ hessian[i * num_players + i] += 1.0f / (prior_sigma * prior_sigma);
}
for (unsigned i = 0; i < all_matches.size(); ++i) {
const match &m = all_matches[i];
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 * num_players + p2] -= w / sigma_sq;
+ hessian[p2 * num_players + p1] -= w / sigma_sq;
- hessian[p1][p1] += w / sigma_sq;
- hessian[p2][p2] += w / sigma_sq;
+ hessian[p1 * num_players + p1] += w / sigma_sq;
+ hessian[p2 * num_players + p2] += w / sigma_sq;
}
}
-// Compute uncertainty (stddev) of mu estimates, which is 1/sqrt(H_ii),
+// 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, int num_players)
{
- memset(mu_stddev, 0, sizeof(mu_stddev));
-
- // Temporarily use mu_stddev to store the diagonal of the Hessian.
-
- // Prior.
+ Matrix<float, Dynamic, Dynamic> h(num_players, num_players);
for (int i = 0; i < num_players; ++i) {
- mu_stddev[i] += 1.0f / (prior_sigma * prior_sigma);
- }
-
- // Matches.
- for (unsigned i = 0; i < all_matches.size(); ++i) {
- const match &m = all_matches[i];
-
- int p1 = m.player;
- int p2 = m.other_player;
-
- double sigma_sq = global_sigma * global_sigma;
- float w = m.weight;
-
- mu_stddev[p1] += w / sigma_sq;
- mu_stddev[p2] += w / sigma_sq;
+ for (int j = 0; j < num_players; ++j) {
+ h(i, j) = hessian[i * num_players + j];
+ }
}
- // Now convert to standard deviation.
+ // FIXME: Use pseudoinverse if applicable.
+ Matrix<float, Dynamic, Dynamic> ih = h.inverse();
for (int i = 0; i < num_players; ++i) {
- mu_stddev[i] = 1.0f / sqrt(mu_stddev[i]);
+ mu_stddev[i] = sqrt(ih(i, i));
}
}
#if DUMP_RAW
dump_raw(mu, num_players);
#else
+ construct_hessian(mu, num_players);
compute_mu_uncertainty(mu, num_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));
float total_logl = compute_total_logl(mu, num_players);
printf("%f 0 -4\n", total_logl);
-
-// construct_hessian(mu, sigma, num_players);
#endif
}