cleanup, and change to initalization with walsh matrices

1 files changed, 51 insertions, 20 deletions

diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index 4647396..5ef062f 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -37,6 +37,28 @@ public:
   }
 };
 
+Orthogonal walsh(unsigned k) {
+  if (k == 0) {
+    return Orthogonal(1);
+  } else {
+    Orthogonal s = walsh(k - 1);
+    Orthogonal t = Orthogonal(2 * s.size());
+
+    for (unsigned i = 0; i < s.size(); i++) {
+      for (unsigned j = 0; j < s.size(); j++) {
+        double sij = s(i, j);
+
+        t(i, j) = sij;
+        t(s.size() + i, j) = sij;
+        t(i, s.size() + j) = sij;
+        t(s.size() + i, s.size() + j) = -sij;
+      }
+    }
+
+    return t;
+  }
+}
+
 class Givens {
 private:
   Orthogonal& m;
@@ -137,36 +159,45 @@ public:
   bool step(Givens& g) {
     double ΔE = g.tryRotation();
 
-    if (ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0)) {
+    bool accepted = ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0);
+
+    if (accepted) {
       E += ΔE;
       g.acceptRotation();
-      A.after_step(true, g, θ0, E, ΔE, M);
-      return true;
-    } else {
-      return false;
-      A.after_step(false, g, θ0, E, ΔE, M);
     }
-  }
 
-  void tune(unsigned N, double ε) {
-    for (unsigned i = 0; i < N; i++) {
-      Givens g(M, θ0, rng);
-      bool stepAccepted = this->step(g);
-      if (stepAccepted) {
-        θ0 *= 1 + ε;
-      } else {
-        θ0 /= 1 + ε;
-      }
-    }
+    A.after_step(accepted, g, θ0, E, ΔE, M);
+    return accepted;
   }
 
-  void sweep() {
+  double sweep() {
+    unsigned total = 0;
+    unsigned accepted = 0;
+
     for (unsigned i = 0; i < M.size() - 1; i++) {
       for (unsigned j = i + 1; j < M.size(); j++) {
         Givens g1(M, false, i, j, θ0, rng);
-        this->step(g1);
         Givens g2(M, true, i, j, θ0, rng);
-        this->step(g2);
+
+        if (this->step(g1))
+          accepted++;
+        if (this->step(g2))
+          accepted++;
+
+        total += 2;
+      }
+    }
+
+    return (double)accepted / (double)total;
+  }
+
+  void tune(unsigned N, double ε) {
+    for (unsigned i = 0; i < N; i++) {
+      double ratio_accepted = this->sweep();
+      if (ratio_accepted > 0.5) {
+        θ0 *= 1 + ε;
+      } else {
+        θ0 /= 1 + ε;
       }
     }
   }