cleanup, and change to initalization with walsh matrices

3 files changed, 129 insertions, 92 deletions

diff --git a/hadamard.cpp b/hadamard.cpp
index 6b9a64a..15f73b1 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -1,57 +1,54 @@
 
 #include "hadamard_pt.hpp"
 
-#include <iostream>
 #include <fstream>
+#include <iostream>
 
 class MeasureEnergy : public Measurement {
-  public:
-    unsigned N;
-    double totalE;
-    MeasureEnergy() {
-      N = 0;
-      totalE = 0;
-    }
+public:
+  unsigned N;
+  double totalE;
 
-    void after_sweep(double, double E, const Orthogonal&) override {
-      totalE += E;
-      N++;
-    }
+  MeasureEnergy() {
+    N = 0;
+    totalE = 0;
+  }
 
-    double energy() const {
-      return totalE / N;
-    }
+  void after_sweep(double, double E, const Orthogonal&) override {
+    totalE += E;
+    N++;
+  }
+
+  double energy() const { return totalE / N; }
 };
 
 class MeasureTransitionRates : public ParallelMeasurement {
-  public:
-    std::vector<std::vector<unsigned>> nAccepted;
-    unsigned total_steps;
-
-    MeasureTransitionRates(unsigned n) : nAccepted(n - 1) {
-      total_steps = 0;
-      for (unsigned i = 0; i < n - 1; i++) {
-        nAccepted[i].resize(i + 1);
-      }
+public:
+  std::vector<std::vector<unsigned>> nAccepted;
+  unsigned total_steps;
+
+  MeasureTransitionRates(unsigned n) : nAccepted(n - 1) {
+    total_steps = 0;
+    for (unsigned i = 0; i < n - 1; i++) {
+      nAccepted[i].resize(i + 1);
     }
+  }
 
-    void after_step(bool accepted, unsigned i, unsigned j, double, double, const MCMC&, const MCMC&) override {
-      if (accepted) {
-        nAccepted[j - 1][i]++;
-      }
-    }
+  void after_step(bool accepted, unsigned i, unsigned j, double, double, const MCMC&,
+                  const MCMC&) override {
+    if (accepted)
+      nAccepted[j - 1][i]++;
+  }
 
-    void after_sweep(const std::vector<MCMC>&) override {
-      total_steps++;
-    }
+  void after_sweep(const std::vector<MCMC>&) override { total_steps++; }
 };
 
 int main(int argc, char* argv[]) {
-  unsigned trials = 1e6;
+  unsigned n_tuning = 1e2;
   double β0 = 1.00;
   double βf = 40.00;
   unsigned num = 40;
-  unsigned n = 20;
+  unsigned k = 2;
   double ε = 0.01;
 
   unsigned M = 10;
@@ -59,37 +56,39 @@ int main(int argc, char* argv[]) {
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "d:b:c:n:t:N:M:e:")) != -1) {
+  while ((opt = getopt(argc, argv, "k:b:c:n:t:N:M:e:")) != -1) {
     switch (opt) {
-      case 'd':
-        n = atoi(optarg);
-        break;
-      case 'b':
-        β0 = atof(optarg);
-        break;
-      case 'c':
-        βf = atof(optarg);
-        break;
-      case 'e':
-        ε = atof(optarg);
-        break;
-      case 'n':
-        num = atof(optarg);
-        break;
-      case 't':
-        trials = (unsigned)atof(optarg);
-        break;
-      case 'N':
-        N = (unsigned)atof(optarg);
-        break;
-      case 'M':
-        M = (unsigned)atof(optarg);
-        break;
-      default:
-        exit(1);
+    case 'k':
+      k = atoi(optarg);
+      break;
+    case 'b':
+      β0 = atof(optarg);
+      break;
+    case 'c':
+      βf = atof(optarg);
+      break;
+    case 'e':
+      ε = atof(optarg);
+      break;
+    case 'n':
+      num = atof(optarg);
+      break;
+    case 't':
+      n_tuning = (unsigned)atof(optarg);
+      break;
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'M':
+      M = (unsigned)atof(optarg);
+      break;
+    default:
+      exit(1);
     }
   }
 
+  unsigned n = pow(2, k);
+
   std::vector<Measurement*> As(num);
   for (Measurement*& A : As) {
     A = new MeasureEnergy();
@@ -98,13 +97,19 @@ int main(int argc, char* argv[]) {
 
   PT p(β0, βf, num, n, B, As);
 
-  std::cout << "Beginning " << trials << " tuning steps for " << n << ".\n";
-  p.tune(trials, ε);
-  std::cout << "Finished tuning, beginning " << N << " tempering updates of " << M << " sweeps each.\n";
+  for (MCMC& sim : p.Ms) {
+    sim.M = walsh(k);
+  }
+
+  std::cout << "Beginning " << n_tuning << " tuning steps for " << n << ".\n";
+  p.tune(n_tuning, ε);
+  std::cout << "Finished tuning, beginning " << N << " tempering updates of " << M
+            << " sweeps each.\n";
   p.run(N, M);
   std::cout << "Finished " << n << ".\n";
 
-  std::string filename = "probs_" + std::to_string(n) + "_" + std::to_string(β0) + "_" + std::to_string(βf) + "_" + std::to_string(num) +  ".dat";
+  std::string filename = "probs_" + std::to_string(n) + "_" + std::to_string(β0) + "_" +
+                         std::to_string(βf) + "_" + std::to_string(num) + ".dat";
   std::ifstream file(filename);
 
   unsigned N_old = 0;
@@ -141,7 +146,8 @@ int main(int argc, char* argv[]) {
 
   file_out.close();
 
-  std::string efilename = "energies_" + std::to_string(n) + "_" + std::to_string(β0) + "_" + std::to_string(βf) + "_" + std::to_string(num) +  ".dat";
+  std::string efilename = "energies_" + std::to_string(n) + "_" + std::to_string(β0) + "_" +
+                          std::to_string(βf) + "_" + std::to_string(num) + ".dat";
   std::ifstream efile(efilename);
 
   unsigned Ne_old = 0;
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 + ε;
       }
     }
   }
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index 4006e0e..404129c 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -9,8 +9,8 @@ void swap(MCMC& s1, MCMC& s2) {
 
 class ParallelMeasurement {
 public:
-  virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&) {};
-  virtual void after_sweep(const std::vector<MCMC>&) {};
+  virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&){};
+  virtual void after_sweep(const std::vector<MCMC>&){};
 };
 
 class PT {
@@ -22,7 +22,8 @@ public:
   ParallelMeasurement& B;
   std::vector<Measurement*>& As;
 
-  PT(double β0, double β1, unsigned N, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As)
+  PT(double β0, double β1, unsigned N, unsigned n, ParallelMeasurement& B,
+     std::vector<Measurement*>& As)
       : B(B), As(As) {
     Ms.reserve(N);
     for (unsigned i = 0; i < N; i++) {
@@ -44,12 +45,10 @@ public:
 
     bool accepted = Δβ * ΔE > 0 || exp(Δβ * ΔE) > rng.uniform((double)0.0, 1.0);
 
-    if (accepted) {
+    if (accepted)
       swap(Ms[i], Ms[j]);
-    }
 
     B.after_step(accepted, i, j, Δβ, ΔE, Ms[i], Ms[j]);
-
     return accepted;
   }
 
@@ -68,6 +67,7 @@ public:
         Ms[j].run(m);
       }
       this->sweep();
+      B.after_sweep(this->Ms);
     }
   }
 };