Use the matrix class from Eigen directly instead of copying from an array. Seems...

diff --git a/bayeswf.cpp b/bayeswf.cpp
index ca205ef47e4eb9cf75f9b27370d280e60d431b46..a691b8fe915ac8f6bb8714e62bd2349c173201b1 100644 (file)
--- a/bayeswf.cpp
+++ b/bayeswf.cpp
@@ -199,13 +199,14 @@ 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];
+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 * num_players + i] += 1.0f / (prior_sigma * prior_sigma);
+               hessian(i, i) += 1.0f / (prior_sigma * prior_sigma);
        }
        for (unsigned i = 0; i < all_matches.size(); ++i) {
                const match &m = all_matches[i];
@@ -216,11 +217,11 @@ void construct_hessian(const float *mu, int num_players)
                double sigma_sq = global_sigma * global_sigma;
                float w = m.weight;
 
-               hessian[p1 * num_players + p2] -= w / sigma_sq;
-               hessian[p2 * num_players + p1] -= w / sigma_sq;
+               hessian(p1, p2) -= w / sigma_sq;
+               hessian(p2, p1) -= w / sigma_sq;
 
-               hessian[p1 * num_players + p1] += w / sigma_sq;
-               hessian[p2 * num_players + p2] += w / sigma_sq;
+               hessian(p1, p1) += w / sigma_sq;
+               hessian(p2, p2) += w / sigma_sq;
        }
 }
 
@@ -228,15 +229,8 @@ void construct_hessian(const float *mu, int num_players)
 // where H is the Hessian (see construct_hessian()).
 void compute_mu_uncertainty(const float *mu, int num_players)
 {
-       Matrix<float, Dynamic, Dynamic> h(num_players, num_players);
-       for (int i = 0; i < num_players; ++i) {
-               for (int j = 0; j < num_players; ++j) {
-                       h(i, j) = hessian[i * num_players + j];
-               }
-       }
-
        // FIXME: Use pseudoinverse if applicable.
-       Matrix<float, Dynamic, Dynamic> ih = h.inverse();
+       Matrix<float, Dynamic, Dynamic> ih = hessian.inverse();
        for (int i = 0; i < num_players; ++i) {
                mu_stddev[i] = sqrt(ih(i, i));
        }