Added new code for measurements, and utility for investigating the Kauzmann paradox.

2 files changed, 261 insertions, 0 deletions

diff --git a/kauzmann.cpp b/kauzmann.cpp
new file mode 100644
index 0000000..fe0c242
--- /dev/null
+++ b/kauzmann.cpp
@@ -0,0 +1,151 @@
+
+#include "hadamard_mcmc.hpp"
+#include "quantity.hpp"
+#include "matrices.hpp"
+
+#include <chrono>
+#include <fstream>
+#include <iostream>
+
+std::string getFilename(unsigned n, double β, double θ₀) {
+  return "kauzmann_" + std::to_string(n) + "_" + std::to_string(β) + "_" + std::to_string(θ₀) + ".dat";
+}
+
+class MeasureEnergy : public Measurement {
+public:
+  Quantity Eavg;
+
+  MeasureEnergy(unsigned lag) : Eavg(lag) {}
+
+  void after_sweep(double, double E, const Orthogonal& M) override {
+    Eavg.add(E);
+  }
+};
+
+int main(int argc, char* argv[]) {
+  // model parameters
+  unsigned n = 2; // matrix size over four
+
+  // temperatures to simulate
+  double β₀ = 1; // starting temperature
+  double Δβ = 0.5; // step size
+  double β₁ = 50; // ending temperature
+
+  // simulation settings
+  unsigned m = 1e4; // number of tuning sweeps at each temperature
+  unsigned N = 1e3; // number of autocorrelation times to simulate
+  double ε = 0.1; // uncertainty on autocorrelation time
+  double θ₀ = 0.05; // standard deviation of step size
+
+  // measurement settings
+  unsigned l = 1e3; // lag in autocorrelation function
+
+  bool loadInitialFromFile = false;
+  std::string inFilename;
+
+  bool loadDataFromFile = false;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "n:b:d:c:m:N:e:l:t:i:L")) != -1) {
+    switch (opt) {
+    case 'n':
+      n = atoi(optarg);
+      break;
+    case 'b':
+      β₀ = atof(optarg);
+      break;
+    case 'd':
+      Δβ = atof(optarg);
+      break;
+    case 'c':
+      β₁ = atof(optarg);
+      break;
+    case 'm':
+      m = (unsigned)atof(optarg);
+      break;
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'e':
+      ε = atof(optarg);
+      break;
+    case 'l':
+      l = (unsigned)atof(optarg);
+      break;
+    case 't':
+      θ₀ = atof(optarg);
+      break;
+    case 'i':
+      loadInitialFromFile = true;
+      inFilename = optarg;
+      break;
+    case 'L':
+      loadDataFromFile = true;
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  double β = β₀;
+
+  MeasureEnergy energyMeasurement(l);
+  MCMC simulation(n, β₀, energyMeasurement, θ₀);
+
+  if (loadInitialFromFile) {
+    std::ifstream inFile(inFilename);
+    inFile.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
+    inFile.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
+    for (unsigned i = 0; i < n; i++) {
+      for (unsigned j = 0; j < n; j++) {
+        inFile >> simulation.M(i, j);
+      }
+    }
+
+    simulation.E = simulation.M.energy();
+
+    std::cout << "Kauzmann: Imported matrix from file, energy " << simulation.E << "." << std::endl;
+  }
+
+  if (loadDataFromFile) {
+    energyMeasurement.Eavg.read(getFilename(n, β, θ₀));
+    std::cout << "Kauzmann: Imported data from file, number of steps " << energyMeasurement.Eavg.num_added() << "." << std::endl;
+  }
+
+  while (β <= β₁) {
+    std::cout << "Kauzmann: Starting simulation of " << β << "." << std::endl;
+    double rat = 0;
+    unsigned num = 0;
+    simulation.run(m, true);
+    simulation.run(l);
+    while (true) {
+      Quantity EavgBak = energyMeasurement.Eavg;
+      Orthogonal Mbak = simulation.M;
+      rat += simulation.run(m);
+      num++;
+      std::array<double, 2> τ = energyMeasurement.Eavg.τ();
+      std::cout << "Kauzmann: After " << energyMeasurement.Eavg.num_added() << " steps the autocorrelation time is " << τ[0] << " ± " << τ[1] << "." << std::endl;
+      if (τ[1] / τ[0] < ε && τ[0] * N <= energyMeasurement.Eavg.num_added()) {
+        break;
+      }
+    }
+    std::cout << "Kauzmann: Finished simulation of " << β << ", " << rat / num << " acceptance ratio, writing output file." << std::endl;
+    std::string filename = getFilename(n, β, θ₀);
+    energyMeasurement.Eavg.write(filename);
+    std::ofstream outFile(filename, std::ios::app);
+    for (unsigned i = 0; i < n; i++) {
+      for (unsigned j = 0; j < n; j++) {
+        outFile << simulation.M(i, j) << " ";
+      }
+    }
+    outFile << std::endl;
+    outFile.close();
+    energyMeasurement.Eavg.reset();
+    β += Δβ;
+    simulation.β = β;
+  }
+
+  return 0;
+}
+
diff --git a/quantity.hpp b/quantity.hpp
new file mode 100644
index 0000000..296a6e8
--- /dev/null
+++ b/quantity.hpp
@@ -0,0 +1,110 @@
+
+#pragma once
+
+#include <array>
+#include <cmath>
+#include <list>
+#include <vector>
+#include <fstream>
+#include <iomanip>
+
+class Quantity {
+private:
+  unsigned n;
+  std::list<double> hist;
+  double total;
+  double total2;
+  std::vector<double> C;
+
+public:
+  Quantity(unsigned lag) : C(lag) {
+    n = 0;
+    total = 0;
+    total2 = 0;
+  }
+
+  void read(std::string filename) {
+    std::ifstream inFile(filename);
+    inFile >> n;
+    inFile >> total;
+    inFile >> total2;
+
+    for (double& Ci : C) {
+      inFile >> Ci;
+    }
+  }
+
+  void write(std::string filename) const {
+    std::ofstream outFile(filename);
+    outFile << std::setprecision(15) << n << " " << total << " " << total2 << std::endl;
+    for (double Ci : C) {
+      outFile << Ci << " ";
+    }
+    outFile << std::endl;
+    outFile.close();
+  }
+
+  void reset() {
+    total = 0;
+    total2 = 0;
+    std::fill(C.begin(), C.end(), 0);
+    n = 0;
+    hist = {};
+  }
+
+  void add(double x) {
+    hist.push_front(x);
+    if (hist.size() > C.size()) {
+      hist.pop_back();
+      unsigned t = 0;
+      for (double a : hist) {
+        C[t] += a * x;
+        t++;
+      }
+      total += x;
+      total2 += pow(x, 2);
+      n++;
+    }
+  }
+
+  double avg() const { return total / n; }
+
+  double avg2() const { return total2 / n; }
+
+  std::vector<double> ρ() const {
+    double C0 = C.front() / n;
+    double avg2 = pow(total / n, 2);
+
+    std::vector<double> ρtmp;
+
+    for (double Ct : C) {
+      ρtmp.push_back((Ct / n - avg2) / (C0 - avg2));
+    }
+
+    return ρtmp;
+  }
+
+  std::array<double, 2> τ() const {
+    double τtmp = 0.5;
+    unsigned M = 1;
+    double c = 8.0;
+
+    std::vector<double> ρ_tmp = this->ρ();
+
+    while (c * τtmp > M && M < C.size()) {
+      τtmp += ρ_tmp[M];
+      M++;
+    }
+
+    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
+  }
+
+  double σ() const {
+    return 2.0 / n * this->τ()[0] * (C[0] / n - pow(this->avg(), 2));
+  }
+
+  double serr() const { return sqrt(this->σ()); }
+
+  unsigned num_added() const { return n; }
+};
+