From 5a40016b1512406bfda74c16b3ddb11347050981 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 12 May 2020 00:50:59 -0400 Subject: Added new code for measurements, and utility for investigating the Kauzmann paradox. --- kauzmann.cpp | 151 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ quantity.hpp | 110 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 261 insertions(+) create mode 100644 kauzmann.cpp create mode 100644 quantity.hpp 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 +#include +#include + +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::max(),'\n'); + inFile.ignore(std::numeric_limits::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 τ = 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 +#include +#include +#include +#include +#include + +class Quantity { +private: + unsigned n; + std::list hist; + double total; + double total2; + std::vector 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 ρ() const { + double C0 = C.front() / n; + double avg2 = pow(total / n, 2); + + std::vector ρtmp; + + for (double Ct : C) { + ρtmp.push_back((Ct / n - avg2) / (C0 - avg2)); + } + + return ρtmp; + } + + std::array τ() const { + double τtmp = 0.5; + unsigned M = 1; + double c = 8.0; + + std::vector ρ_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; } +}; + -- cgit v1.2.3-70-g09d2