initialized repo

1 files changed, 74 insertions, 0 deletions

diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
new file mode 100644
index 0000000..4006e0e
--- /dev/null
+++ b/hadamard_pt.hpp
@@ -0,0 +1,74 @@
+
+#pragma once
+#include "hadamard_mcmc.hpp"
+
+void swap(MCMC& s1, MCMC& s2) {
+  std::swap(s1.M, s2.M);
+  std::swap(s1.E, s2.E);
+}
+
+class ParallelMeasurement {
+public:
+  virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&) {};
+  virtual void after_sweep(const std::vector<MCMC>&) {};
+};
+
+class PT {
+private:
+  randutils::mt19937_rng rng;
+
+public:
+  std::vector<MCMC> Ms;
+  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++) {
+      double β = β0 + i * (β1 - β0) / (N - 1);
+      Ms.push_back(MCMC(n, β, *As[i]));
+    }
+  }
+
+  void tune(unsigned N, double ε) {
+#pragma omp parallel for
+    for (unsigned i = 0; i < Ms.size(); i++) {
+      Ms[i].tune(N, ε);
+    }
+  }
+
+  bool step(unsigned i, unsigned j) {
+    double Δβ = Ms[i].β - Ms[j].β;
+    double ΔE = Ms[i].E - Ms[j].E;
+
+    bool accepted = Δβ * ΔE > 0 || exp(Δβ * ΔE) > rng.uniform((double)0.0, 1.0);
+
+    if (accepted) {
+      swap(Ms[i], Ms[j]);
+    }
+
+    B.after_step(accepted, i, j, Δβ, ΔE, Ms[i], Ms[j]);
+
+    return accepted;
+  }
+
+  void sweep() {
+    for (unsigned i = 0; i < Ms.size() - 1; i++) {
+      for (unsigned j = i + 1; j < Ms.size(); j++) {
+        this->step(i, j);
+      }
+    }
+  }
+
+  void run(unsigned n, unsigned m) {
+    for (unsigned i = 0; i < n; i++) {
+#pragma omp parallel for
+      for (unsigned j = 0; j < Ms.size(); j++) {
+        Ms[j].run(m);
+      }
+      this->sweep();
+    }
+  }
+};
+