1 files changed, 125 insertions, 0 deletions
diff --git a/hmm_aging.cpp b/hmm_aging.cpp
new file mode 100644
index 0000000..522dac6
--- /dev/null
+++ b/hmm_aging.cpp
@@ -0,0 +1,125 @@
+
+#include "hadamard_mcmc.hpp"
+#include "matrices.hpp"
+
+#include <chrono>
+#include <fstream>
+#include <iostream>
+#include <list>
+
+std::string getFilename(unsigned n, double β1, double β2, double θ₀, unsigned m, unsigned s) {
+  auto tag = std::chrono::high_resolution_clock::now();
+  return "aging_" + std::to_string(n) + "_" + std::to_string(β1) + "_" + std::to_string(β2) + "_" +
+    std::to_string(θ₀) + "_" + std::to_string(m) + "_" + std::to_string(s) + "_" + std::to_string(tag.time_since_epoch().count()) + ".dat";
+}
+
+class Correlation {
+  public:
+    Orthogonal M0;
+    std::list<double> C;
+
+    Correlation(const Orthogonal& M) : M0(M) {}
+
+    void add(const Orthogonal& M) {
+      C.push_back(M0 * M);
+    }
+};
+
+class MeasureCorrelation : public Measurement {
+public:
+  unsigned s;
+  unsigned n;
+  unsigned N;
+  std::list<Correlation> Cs;
+
+  MeasureCorrelation(unsigned m, unsigned skip) : N(m), s(skip) {
+    n = 0;
+  }
+
+  void after_sweep(double, double, const Orthogonal& M) override {
+    if (n % N == 0) {
+      Cs.push_back(Correlation(M));
+    }
+
+    if (n % s == 0) {
+      for (Correlation& C : Cs) {
+        C.add(M);
+      }
+    }
+
+    n++;
+  }
+};
+
+int main(int argc, char* argv[]) {
+  // model parameters
+  unsigned n = 2; // matrix size over four
+
+  // simulation settings
+  double β1 = 5; // temperature
+  double β2 = 15; // temperature
+  unsigned M = 1e4; // number of relaxation sweeps
+  unsigned N = 1e4; // number of measurement sweeps
+  unsigned m = 1e3;
+  unsigned s = 10;
+  double θ₀ = 0.05; // standard deviation of step size
+
+  bool loadDataFromFile = false;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "n:b:c:M:N:m:t:s:")) != -1) {
+    switch (opt) {
+    case 'n':
+      n = atoi(optarg);
+      break;
+    case 'b':
+      β1 = atof(optarg);
+      break;
+    case 'c':
+      β2 = atof(optarg);
+      break;
+    case 'M':
+      M = (unsigned)atof(optarg);
+      break;
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'm':
+      m = (unsigned)atof(optarg);
+      break;
+    case 't':
+      θ₀ = atof(optarg);
+      break;
+    case 's':
+      s = (unsigned)atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  MeasureCorrelation Q(m, s);
+  MCMC simulation(n, β1, Q, θ₀);
+
+  simulation.run(M, true);
+  std::cout << "Finished initial relaxation at " << β1 << ", quenching to " << β2 << "." << std::endl;
+
+  simulation.β = β2;
+  simulation.run(N * m);
+
+  std::cout << "Finished simulation." << std::endl;
+
+  std::string filename = getFilename(n, β1, β2, θ₀, m, s);
+  std::ofstream file(filename);
+
+  for (const Correlation& C : Q.Cs) {
+    for (double x : C.C) {
+      file << x << " ";
+    }
+    file << std::endl;
+  }
+
+  return 0;
+}
+