1 files changed, 94 insertions, 0 deletions
diff --git a/hmm_correlation.cpp b/hmm_correlation.cpp
new file mode 100644
index 0000000..80a6927
--- /dev/null
+++ b/hmm_correlation.cpp
@@ -0,0 +1,94 @@
+
+#include "hadamard_mcmc.hpp"
+#include "quantity.hpp"
+#include "matrices.hpp"
+
+#include <chrono>
+#include <fstream>
+#include <iostream>
+
+std::string getFilename(unsigned n, double β, double θ₀, unsigned lag, unsigned skip) {
+  return "correlation_" + std::to_string(n) + "_" + std::to_string(β) + "_" +
+    std::to_string(θ₀) + "_" + std::to_string(lag) + "_" + std::to_string(skip)
+    + ".dat";
+}
+
+class MeasureCorrelation : public Measurement {
+public:
+  Quantity<Orthogonal> Q;
+
+  MeasureCorrelation(unsigned lag, unsigned skip) : Q(lag, skip) {}
+
+  void after_sweep(double, double E, const Orthogonal& M) override {
+    Q.add(M);
+  }
+};
+
+int main(int argc, char* argv[]) {
+  // model parameters
+  unsigned n = 2; // matrix size over four
+
+  // simulation settings
+  double β = 1; // temperature
+  unsigned m = 1e4; // number of relaxation sweeps
+  unsigned N = 1e4; // number of measurement sweeps
+  double θ₀ = 0.05; // standard deviation of step size
+
+  // measurement settings
+  unsigned l = 1e3; // lag in autocorrelation function
+  unsigned s = 1; // skip in autocorrelation function
+
+  bool loadDataFromFile = false;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "n:b:m:N:l:t:i:Ls:")) != -1) {
+    switch (opt) {
+    case 'n':
+      n = atoi(optarg);
+      break;
+    case 'b':
+      β = atof(optarg);
+      break;
+    case 'm':
+      m = (unsigned)atof(optarg);
+      break;
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'l':
+      l = (unsigned)atof(optarg);
+      break;
+    case 's':
+      s = (unsigned)atof(optarg);
+      break;
+    case 't':
+      θ₀ = atof(optarg);
+      break;
+    case 'L':
+      loadDataFromFile = true;
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  MeasureCorrelation Q(l, s);
+  MCMC simulation(n, β, Q, θ₀);
+
+  if (loadDataFromFile) {
+    Q.Q.read(getFilename(n, β, θ₀, l, s));
+    std::cout << "Imported data from file, number of steps " << Q.Q.num_added() << "." << std::endl;
+  }
+
+  simulation.run(m, true);
+  std::cout << "Finished initial relaxation at " << β << "." << std::endl;
+  simulation.run(N);
+
+  std::cout << "Finished simulation of " << β << ", writing output file." << std::endl;
+  std::string filename = getFilename(n, β, θ₀, l, s);
+  Q.Q.write(filename);
+
+  return 0;
+}
+