From a8a8695211c3a5eb8873b0818c30a86f3284b136 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 9 Sep 2019 15:18:27 -0400
Subject: initial commit

---
 .gitmodules    |   6 ++
 metastable.cpp | 172 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 randutils      |   1 +
 wolff          |   1 +
 4 files changed, 180 insertions(+)
 create mode 100644 .gitmodules
 create mode 100644 metastable.cpp
 create mode 160000 randutils
 create mode 160000 wolff

diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..a509f25
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,6 @@
+[submodule "randutils"]
+	path = randutils
+	url = https://gist.github.com/imneme/540829265469e673d045
+[submodule "wolff"]
+	path = wolff
+	url = git:research/wolff/code/c++
diff --git a/metastable.cpp b/metastable.cpp
new file mode 100644
index 0000000..cb47ffc
--- /dev/null
+++ b/metastable.cpp
@@ -0,0 +1,172 @@
+
+#define WOLFF_USE_FINITE_STATES
+
+#include <iostream>
+#include <cstring>
+#include <fstream>
+
+
+#include "randutils/randutils.hpp"
+#include "wolff/lib/wolff_models/ising.hpp"
+
+using namespace wolff;
+
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, std::string model_string) {
+  std::string filename = model_string + id + ".dat";
+  std::ifstream file(filename);
+
+  std::vector<uint64_t> data_old(data.size(), 0);
+
+  if (file.is_open()) {
+    for (unsigned i = 0; i < data.size(); i++) {
+      uint64_t num;
+      file >> num;
+      data_old[i] = num;
+    }
+
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    file_out <<std::fixed<< data_old[i] + data[i] << " ";
+  }
+
+  file_out.close();
+}
+
+void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data, std::string model_string) {
+  std::string filename = model_string + id + ".dat";
+  std::ifstream file(filename);
+
+  std::vector<std::vector<uint64_t>> data_old(data.size());
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    data_old[i].resize(data[0].size());
+  }
+
+  if (file.is_open()) {
+    for (unsigned i = 0; i < data.size(); i++) {
+      for (unsigned j = 0; j < data[0].size(); j++) {
+        uint64_t num;
+        file >> num;
+        data_old[i][j] = num;
+      }
+    }
+
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    for (unsigned j = 0; j < data[0].size(); j++) {
+      file_out << std::fixed << data_old[i][j] + data[i][j] << " ";
+    }
+    file_out << "\n";
+  }
+
+  file_out.close();
+}
+
+class meas : public measurement<ising_t, ising_t, graph<>> {
+private:
+  std::vector<uint64_t> mag_dist;
+  std::vector<std::vector<uint64_t>> energy_mag_dist;
+
+  signed E;
+  signed M;
+
+  std::string model_string;
+
+public:
+  meas(const system<ising_t, ising_t, graph<>>& S, double H) : mag_dist(S.nv + 1, 0), energy_mag_dist(S.nv + 1) {
+    M = S.nv * S.s[0];
+    E = S.ne;
+    for (std::vector<uint64_t>& d : energy_mag_dist) {
+      d.resize(S.ne + 1, 0);
+    }
+    model_string = "metastable_" + std::to_string(S.G.D) + "_" + std::to_string(S.G.L) + "_" + std::to_string(S.T) + "_" + std::to_string(H) + "_";
+  }
+
+  ~meas() {
+    update_distribution_file("m", mag_dist, model_string);
+    update_distribution_file("e", energy_mag_dist, model_string);
+  }
+
+  void plain_bond_visited(const system<ising_t, ising_t, graph<>>&, const typename graph<>::halfedge&, const ising_t&,
+                          double dE) override {
+    if (dE > 0) {
+      E -= 2;
+    } else {
+      E += 2;
+    }
+  }
+
+  void ghost_bond_visited(const system<ising_t, ising_t, graph<>>&, const typename graph<>::vertex&,
+                          const ising_t& s_old, const ising_t& s_new, double dE) override {
+    M += s_new - s_old;
+  }
+
+  void post_cluster(unsigned, unsigned, const system<ising_t, ising_t, graph<>>& S) override {
+    mag_dist[(S.nv + M) / 2]++;
+    energy_mag_dist[(S.nv + M) / 2][(E + S.ne) / 2]++;
+  }
+};
+
+int main(int argc, char* argv[]) {
+  unsigned N = (unsigned)1e4;
+  unsigned D = 2;
+  unsigned L = 128;
+  double T = 2.26918531421;
+  double H = 0.0;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
+    switch (opt) {
+    case 'N': // number of steps
+      N = (unsigned)atof(optarg);
+      break;
+    case 'D': // dimension
+      D = atoi(optarg);
+      break;
+    case 'L': // linear size
+      L = atoi(optarg);
+      break;
+    case 'T': // temperature
+      T = atof(optarg);
+      break;
+    case 'H': // external field
+      H = atof(optarg);
+      break;
+    default:
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  // define the spin-spin coupling
+  std::function<double(const ising_t&, const ising_t&)> Z =
+      [](const ising_t& s1, const ising_t& s2) -> double { return (double)(s1 * s2); };
+
+  // define the spin-field coupling
+  std::function<double(const ising_t&)> B = [=](const ising_t& s) -> double { return H * s; };
+
+  // initialize the lattice
+  graph<> G(D, L);
+
+  // initialize the system
+  system<ising_t, ising_t, graph<>> S(G, T, Z, B);
+
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng(seeds);
+
+  meas A(S, H);
+
+  // run wolff N times
+  S.run_wolff(N, gen_ising<graph<>>, A, rng);
+
+  return 0;
+}
+
diff --git a/randutils b/randutils
new file mode 160000
index 0000000..8486a61
--- /dev/null
+++ b/randutils
@@ -0,0 +1 @@
+Subproject commit 8486a610a954a8248c12485fb4cfc390a5f5f854
diff --git a/wolff b/wolff
new file mode 160000
index 0000000..4991418
--- /dev/null
+++ b/wolff
@@ -0,0 +1 @@
+Subproject commit 4991418a416b75680737ec8f1e47b355b219728d
-- 
cgit v1.2.3-54-g00ecf