7 files changed, 467 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..773a6df
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+*.dat
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..94345fd
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "randutils"]
+	path = randutils
+	url = https://gist.github.com/imneme/540829265469e673d045
diff --git a/hadamard.cpp b/hadamard.cpp
new file mode 100644
index 0000000..6b9a64a
--- /dev/null
+++ b/hadamard.cpp
@@ -0,0 +1,174 @@
+
+#include "hadamard_pt.hpp"
+
+#include <iostream>
+#include <fstream>
+
+class MeasureEnergy : public Measurement {
+  public:
+    unsigned N;
+    double totalE;
+    MeasureEnergy() {
+      N = 0;
+      totalE = 0;
+    }
+
+    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);
+      }
+    }
+
+    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++;
+    }
+};
+
+int main(int argc, char* argv[]) {
+  unsigned trials = 1e6;
+  double β0 = 1.00;
+  double βf = 40.00;
+  unsigned num = 40;
+  unsigned n = 20;
+  double ε = 0.01;
+
+  unsigned M = 10;
+  unsigned N = 1e4;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "d: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);
+    }
+  }
+
+  std::vector<Measurement*> As(num);
+  for (Measurement*& A : As) {
+    A = new MeasureEnergy();
+  }
+  MeasureTransitionRates B(num);
+
+  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";
+  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::ifstream file(filename);
+
+  unsigned N_old = 0;
+  std::vector<std::vector<unsigned long>> data_old(B.nAccepted.size());
+
+  for (unsigned i = 0; i < B.nAccepted.size(); i++) {
+    data_old[i].resize(B.nAccepted[i].size());
+  }
+
+  if (file.is_open()) {
+    file >> N_old;
+
+    for (unsigned i = 0; i < B.nAccepted.size(); i++) {
+      for (unsigned j = 0; j < B.nAccepted[i].size(); j++) {
+        double num;
+        file >> num;
+        data_old[i][j] = num;
+      }
+    }
+
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  file_out << N_old + B.total_steps << "\n";
+
+  for (unsigned i = 0; i < B.nAccepted.size(); i++) {
+    for (unsigned j = 0; j < B.nAccepted[i].size(); j++) {
+      file_out << std::fixed << data_old[i][j] + B.nAccepted[i][j] << " ";
+    }
+    file_out << "\n";
+  }
+
+  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::ifstream efile(efilename);
+
+  unsigned Ne_old = 0;
+  std::vector<double> edata_old(p.As.size());
+
+  if (efile.is_open()) {
+    efile >> N_old;
+
+    for (unsigned i = 0; i < p.As.size(); i++) {
+      double num;
+      efile >> num;
+      edata_old[i] = num;
+    }
+
+    efile.close();
+  }
+
+  std::ofstream efile_out(efilename);
+
+  efile_out << Ne_old + ((MeasureEnergy*)As[0])->N << "\n";
+
+  for (unsigned i = 0; i < As.size(); i++) {
+    efile_out << std::fixed << edata_old[i] + ((MeasureEnergy*)As[i])->totalE << " ";
+  }
+
+  efile_out.close();
+
+  return 0;
+}
+
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
new file mode 100644
index 0000000..4647396
--- /dev/null
+++ b/hadamard_mcmc.hpp
@@ -0,0 +1,181 @@
+
+#pragma once
+#include "randutils/randutils.hpp"
+#include <vector>
+
+inline double Ei(double x) { return -fabs(x); }
+
+class Orthogonal {
+private:
+  unsigned d;
+  std::vector<double> m;
+
+public:
+  Orthogonal(unsigned size) : m(size * size) {
+    d = size;
+    for (unsigned i = 0; i < size; i++) {
+      m[size * i + i] = sqrt(size);
+    }
+  }
+
+  unsigned size() const { return d; }
+
+  double& operator()(unsigned i, unsigned j) { return m[d * i + j]; }
+
+  const double& operator()(unsigned i, unsigned j) const { return m[d * i + j]; }
+
+  double energy() const {
+    double E = 0;
+
+    for (unsigned i = 0; i < this->size(); i++) {
+      for (unsigned j = 0; j < this->size(); j++) {
+        E += Ei(this->operator()(i, j));
+      }
+    }
+
+    return E;
+  }
+};
+
+class Givens {
+private:
+  Orthogonal& m;
+
+  bool transpose;
+  unsigned axis_1;
+  unsigned axis_2;
+  double Δθ;
+
+  std::vector<double> rows;
+
+public:
+  Givens(Orthogonal& m, bool t, unsigned a1, unsigned a2, double θ0, randutils::mt19937_rng& rng)
+      : m(m), rows(2 * m.size()) {
+    transpose = t;
+    axis_1 = a1;
+    axis_2 = a2;
+    Δθ = rng.uniform(-θ0, θ0);
+  }
+
+  Givens(Orthogonal& m, double θ0, randutils::mt19937_rng& rng) : m(m), rows(2 * m.size()) {
+    Δθ = rng.uniform(-θ0, θ0);
+    unsigned axis1axis2 = rng.uniform((unsigned)0, m.size() * (m.size() - 1) - 1);
+
+    axis_1 = axis1axis2 / (m.size() - 1);
+    axis_2 = axis1axis2 % (m.size() - 1);
+    transpose = axis_2 >= axis_1;
+    if (transpose) {
+      axis_2++;
+    }
+  }
+
+  double tryRotation() {
+    double ΔE = 0;
+    double c = cos(Δθ);
+    double s = sin(Δθ);
+
+    for (unsigned i = 0; i < m.size(); i++) {
+      double m1i, m2i, m1i_new, m2i_new;
+
+      if (transpose) {
+        m1i = m(i, axis_1);
+        m2i = m(i, axis_2);
+      } else {
+        m1i = m(axis_1, i);
+        m2i = m(axis_2, i);
+      }
+
+      ΔE -= Ei(m1i) + Ei(m2i);
+
+      m1i_new = c * m1i + s * m2i;
+      m2i_new = c * m2i - s * m1i;
+
+      ΔE += Ei(m1i_new) + Ei(m2i_new);
+
+      rows[i] = m1i_new;
+      rows[m.size() + i] = m2i_new;
+    }
+
+    return ΔE;
+  }
+
+  void acceptRotation() const {
+    for (unsigned i = 0; i < m.size(); i++) {
+      if (transpose) {
+        m(i, axis_1) = rows[i];
+        m(i, axis_2) = rows[m.size() + i];
+      } else {
+        m(axis_1, i) = rows[i];
+        m(axis_2, i) = rows[m.size() + i];
+      }
+    }
+  }
+};
+
+class Measurement {
+public:
+  virtual void after_step(bool, const Givens&, double, double, double, const Orthogonal&){};
+  virtual void after_sweep(double, double, const Orthogonal&){};
+};
+
+class MCMC {
+private:
+  randutils::mt19937_rng rng;
+  Measurement& A;
+  double θ0;
+
+public:
+  const double β;
+  double E;
+  Orthogonal M;
+
+  MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0) {
+    θ0 = M_PI;
+    E = M.energy();
+  }
+
+  bool step(Givens& g) {
+    double ΔE = g.tryRotation();
+
+    if (ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0)) {
+      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 + ε;
+      }
+    }
+  }
+
+  void sweep() {
+    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);
+      }
+    }
+  }
+
+  void run(unsigned N) {
+    for (unsigned i = 0; i < N; i++) {
+      this->sweep();
+      A.after_sweep(θ0, E, M);
+    }
+  }
+};
+
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();
+    }
+  }
+};
+
diff --git a/hadamard_time.cpp b/hadamard_time.cpp
new file mode 100644
index 0000000..c9dd961
--- /dev/null
+++ b/hadamard_time.cpp
@@ -0,0 +1,34 @@
+#include "hadamard_mcmc.hpp"
+#include <iostream>
+
+class MeasureEnergy : public Measurement {
+  private:
+    unsigned N;
+    double totalE;
+
+  public:
+    MeasureEnergy() {
+      N = 0;
+      totalE = 0;
+    }
+
+    void after_sweep(double, double E, const Orthogonal&) override {
+      totalE += E;
+      N++;
+    }
+
+    double energy() const {
+      return totalE / N;
+    }
+};
+
+int main() {
+  MeasureEnergy A;
+  MCMC sim(20, 6.0, A);
+  sim.tune(1e4, 0.01);
+  sim.run(1e4);
+
+  std::cout << "The average energy was " << A.energy() << ".";
+
+  return 0;
+}
diff --git a/randutils b/randutils
new file mode 160000
+Subproject 8486a610a954a8248c12485fb4cfc390a5f5f85