summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--.gitignore1
-rw-r--r--.gitmodules3
-rw-r--r--hadamard.cpp174
-rw-r--r--hadamard_mcmc.hpp181
-rw-r--r--hadamard_pt.hpp74
-rw-r--r--hadamard_time.cpp34
m---------randutils0
7 files changed, 467 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..773a6df
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+*.dat
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..94345fd
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "randutils"]
+ path = randutils
+ url = https://gist.github.com/imneme/540829265469e673d045
diff --git a/hadamard.cpp b/hadamard.cpp
new file mode 100644
index 0000000..6b9a64a
--- /dev/null
+++ b/hadamard.cpp
@@ -0,0 +1,174 @@
+
+#include "hadamard_pt.hpp"
+
+#include <iostream>
+#include <fstream>
+
+class MeasureEnergy : public Measurement {
+ public:
+ unsigned N;
+ double totalE;
+ MeasureEnergy() {
+ N = 0;
+ totalE = 0;
+ }
+
+ void after_sweep(double, double E, const Orthogonal&) override {
+ totalE += E;
+ N++;
+ }
+
+ double energy() const {
+ return totalE / N;
+ }
+};
+
+class MeasureTransitionRates : public ParallelMeasurement {
+ public:
+ std::vector<std::vector<unsigned>> nAccepted;
+ unsigned total_steps;
+
+ MeasureTransitionRates(unsigned n) : nAccepted(n - 1) {
+ total_steps = 0;
+ for (unsigned i = 0; i < n - 1; i++) {
+ nAccepted[i].resize(i + 1);
+ }
+ }
+
+ void after_step(bool accepted, unsigned i, unsigned j, double, double, const MCMC&, const MCMC&) override {
+ if (accepted) {
+ nAccepted[j - 1][i]++;
+ }
+ }
+
+ void after_sweep(const std::vector<MCMC>&) override {
+ total_steps++;
+ }
+};
+
+int main(int argc, char* argv[]) {
+ unsigned trials = 1e6;
+ double β0 = 1.00;
+ double βf = 40.00;
+ unsigned num = 40;
+ unsigned n = 20;
+ double ε = 0.01;
+
+ unsigned M = 10;
+ unsigned N = 1e4;
+
+ int opt;
+
+ while ((opt = getopt(argc, argv, "d:b:c:n:t:N:M:e:")) != -1) {
+ switch (opt) {
+ case 'd':
+ n = atoi(optarg);
+ break;
+ case 'b':
+ β0 = atof(optarg);
+ break;
+ case 'c':
+ βf = atof(optarg);
+ break;
+ case 'e':
+ ε = atof(optarg);
+ break;
+ case 'n':
+ num = atof(optarg);
+ break;
+ case 't':
+ trials = (unsigned)atof(optarg);
+ break;
+ case 'N':
+ N = (unsigned)atof(optarg);
+ break;
+ case 'M':
+ M = (unsigned)atof(optarg);
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ std::vector<Measurement*> As(num);
+ for (Measurement*& A : As) {
+ A = new MeasureEnergy();
+ }
+ MeasureTransitionRates B(num);
+
+ PT p(β0, βf, num, n, B, As);
+
+ std::cout << "Beginning " << trials << " tuning steps for " << n << ".\n";
+ p.tune(trials, ε);
+ std::cout << "Finished tuning, beginning " << N << " tempering updates of " << M << " sweeps each.\n";
+ p.run(N, M);
+ std::cout << "Finished " << n << ".\n";
+
+ std::string filename = "probs_" + std::to_string(n) + "_" + std::to_string(β0) + "_" + std::to_string(βf) + "_" + std::to_string(num) + ".dat";
+ std::ifstream file(filename);
+
+ unsigned N_old = 0;
+ std::vector<std::vector<unsigned long>> data_old(B.nAccepted.size());
+
+ for (unsigned i = 0; i < B.nAccepted.size(); i++) {
+ data_old[i].resize(B.nAccepted[i].size());
+ }
+
+ if (file.is_open()) {
+ file >> N_old;
+
+ for (unsigned i = 0; i < B.nAccepted.size(); i++) {
+ for (unsigned j = 0; j < B.nAccepted[i].size(); j++) {
+ double num;
+ file >> num;
+ data_old[i][j] = num;
+ }
+ }
+
+ file.close();
+ }
+
+ std::ofstream file_out(filename);
+
+ file_out << N_old + B.total_steps << "\n";
+
+ for (unsigned i = 0; i < B.nAccepted.size(); i++) {
+ for (unsigned j = 0; j < B.nAccepted[i].size(); j++) {
+ file_out << std::fixed << data_old[i][j] + B.nAccepted[i][j] << " ";
+ }
+ file_out << "\n";
+ }
+
+ file_out.close();
+
+ std::string efilename = "energies_" + std::to_string(n) + "_" + std::to_string(β0) + "_" + std::to_string(βf) + "_" + std::to_string(num) + ".dat";
+ std::ifstream efile(efilename);
+
+ unsigned Ne_old = 0;
+ std::vector<double> edata_old(p.As.size());
+
+ if (efile.is_open()) {
+ efile >> N_old;
+
+ for (unsigned i = 0; i < p.As.size(); i++) {
+ double num;
+ efile >> num;
+ edata_old[i] = num;
+ }
+
+ efile.close();
+ }
+
+ std::ofstream efile_out(efilename);
+
+ efile_out << Ne_old + ((MeasureEnergy*)As[0])->N << "\n";
+
+ for (unsigned i = 0; i < As.size(); i++) {
+ efile_out << std::fixed << edata_old[i] + ((MeasureEnergy*)As[i])->totalE << " ";
+ }
+
+ efile_out.close();
+
+ return 0;
+}
+
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
new file mode 100644
index 0000000..4647396
--- /dev/null
+++ b/hadamard_mcmc.hpp
@@ -0,0 +1,181 @@
+
+#pragma once
+#include "randutils/randutils.hpp"
+#include <vector>
+
+inline double Ei(double x) { return -fabs(x); }
+
+class Orthogonal {
+private:
+ unsigned d;
+ std::vector<double> m;
+
+public:
+ Orthogonal(unsigned size) : m(size * size) {
+ d = size;
+ for (unsigned i = 0; i < size; i++) {
+ m[size * i + i] = sqrt(size);
+ }
+ }
+
+ unsigned size() const { return d; }
+
+ double& operator()(unsigned i, unsigned j) { return m[d * i + j]; }
+
+ const double& operator()(unsigned i, unsigned j) const { return m[d * i + j]; }
+
+ double energy() const {
+ double E = 0;
+
+ for (unsigned i = 0; i < this->size(); i++) {
+ for (unsigned j = 0; j < this->size(); j++) {
+ E += Ei(this->operator()(i, j));
+ }
+ }
+
+ return E;
+ }
+};
+
+class Givens {
+private:
+ Orthogonal& m;
+
+ bool transpose;
+ unsigned axis_1;
+ unsigned axis_2;
+ double Δθ;
+
+ std::vector<double> rows;
+
+public:
+ Givens(Orthogonal& m, bool t, unsigned a1, unsigned a2, double θ0, randutils::mt19937_rng& rng)
+ : m(m), rows(2 * m.size()) {
+ transpose = t;
+ axis_1 = a1;
+ axis_2 = a2;
+ Δθ = rng.uniform(-θ0, θ0);
+ }
+
+ Givens(Orthogonal& m, double θ0, randutils::mt19937_rng& rng) : m(m), rows(2 * m.size()) {
+ Δθ = rng.uniform(-θ0, θ0);
+ unsigned axis1axis2 = rng.uniform((unsigned)0, m.size() * (m.size() - 1) - 1);
+
+ axis_1 = axis1axis2 / (m.size() - 1);
+ axis_2 = axis1axis2 % (m.size() - 1);
+ transpose = axis_2 >= axis_1;
+ if (transpose) {
+ axis_2++;
+ }
+ }
+
+ double tryRotation() {
+ double ΔE = 0;
+ double c = cos(Δθ);
+ double s = sin(Δθ);
+
+ for (unsigned i = 0; i < m.size(); i++) {
+ double m1i, m2i, m1i_new, m2i_new;
+
+ if (transpose) {
+ m1i = m(i, axis_1);
+ m2i = m(i, axis_2);
+ } else {
+ m1i = m(axis_1, i);
+ m2i = m(axis_2, i);
+ }
+
+ ΔE -= Ei(m1i) + Ei(m2i);
+
+ m1i_new = c * m1i + s * m2i;
+ m2i_new = c * m2i - s * m1i;
+
+ ΔE += Ei(m1i_new) + Ei(m2i_new);
+
+ rows[i] = m1i_new;
+ rows[m.size() + i] = m2i_new;
+ }
+
+ return ΔE;
+ }
+
+ void acceptRotation() const {
+ for (unsigned i = 0; i < m.size(); i++) {
+ if (transpose) {
+ m(i, axis_1) = rows[i];
+ m(i, axis_2) = rows[m.size() + i];
+ } else {
+ m(axis_1, i) = rows[i];
+ m(axis_2, i) = rows[m.size() + i];
+ }
+ }
+ }
+};
+
+class Measurement {
+public:
+ virtual void after_step(bool, const Givens&, double, double, double, const Orthogonal&){};
+ virtual void after_sweep(double, double, const Orthogonal&){};
+};
+
+class MCMC {
+private:
+ randutils::mt19937_rng rng;
+ Measurement& A;
+ double θ0;
+
+public:
+ const double β;
+ double E;
+ Orthogonal M;
+
+ MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0) {
+ θ0 = M_PI;
+ E = M.energy();
+ }
+
+ bool step(Givens& g) {
+ double ΔE = g.tryRotation();
+
+ if (ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0)) {
+ E += ΔE;
+ g.acceptRotation();
+ A.after_step(true, g, θ0, E, ΔE, M);
+ return true;
+ } else {
+ return false;
+ A.after_step(false, g, θ0, E, ΔE, M);
+ }
+ }
+
+ void tune(unsigned N, double ε) {
+ for (unsigned i = 0; i < N; i++) {
+ Givens g(M, θ0, rng);
+ bool stepAccepted = this->step(g);
+ if (stepAccepted) {
+ θ0 *= 1 + ε;
+ } else {
+ θ0 /= 1 + ε;
+ }
+ }
+ }
+
+ void sweep() {
+ for (unsigned i = 0; i < M.size() - 1; i++) {
+ for (unsigned j = i + 1; j < M.size(); j++) {
+ Givens g1(M, false, i, j, θ0, rng);
+ this->step(g1);
+ Givens g2(M, true, i, j, θ0, rng);
+ this->step(g2);
+ }
+ }
+ }
+
+ void run(unsigned N) {
+ for (unsigned i = 0; i < N; i++) {
+ this->sweep();
+ A.after_sweep(θ0, E, M);
+ }
+ }
+};
+
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
new file mode 100644
index 0000000..4006e0e
--- /dev/null
+++ b/hadamard_pt.hpp
@@ -0,0 +1,74 @@
+
+#pragma once
+#include "hadamard_mcmc.hpp"
+
+void swap(MCMC& s1, MCMC& s2) {
+ std::swap(s1.M, s2.M);
+ std::swap(s1.E, s2.E);
+}
+
+class ParallelMeasurement {
+public:
+ virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&) {};
+ virtual void after_sweep(const std::vector<MCMC>&) {};
+};
+
+class PT {
+private:
+ randutils::mt19937_rng rng;
+
+public:
+ std::vector<MCMC> Ms;
+ ParallelMeasurement& B;
+ std::vector<Measurement*>& As;
+
+ PT(double β0, double β1, unsigned N, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As)
+ : B(B), As(As) {
+ Ms.reserve(N);
+ for (unsigned i = 0; i < N; i++) {
+ double β = β0 + i * (β1 - β0) / (N - 1);
+ Ms.push_back(MCMC(n, β, *As[i]));
+ }
+ }
+
+ void tune(unsigned N, double ε) {
+#pragma omp parallel for
+ for (unsigned i = 0; i < Ms.size(); i++) {
+ Ms[i].tune(N, ε);
+ }
+ }
+
+ bool step(unsigned i, unsigned j) {
+ double Δβ = Ms[i].β - Ms[j].β;
+ double ΔE = Ms[i].E - Ms[j].E;
+
+ bool accepted = Δβ * ΔE > 0 || exp(Δβ * ΔE) > rng.uniform((double)0.0, 1.0);
+
+ if (accepted) {
+ swap(Ms[i], Ms[j]);
+ }
+
+ B.after_step(accepted, i, j, Δβ, ΔE, Ms[i], Ms[j]);
+
+ return accepted;
+ }
+
+ void sweep() {
+ for (unsigned i = 0; i < Ms.size() - 1; i++) {
+ for (unsigned j = i + 1; j < Ms.size(); j++) {
+ this->step(i, j);
+ }
+ }
+ }
+
+ void run(unsigned n, unsigned m) {
+ for (unsigned i = 0; i < n; i++) {
+#pragma omp parallel for
+ for (unsigned j = 0; j < Ms.size(); j++) {
+ Ms[j].run(m);
+ }
+ this->sweep();
+ }
+ }
+};
+
diff --git a/hadamard_time.cpp b/hadamard_time.cpp
new file mode 100644
index 0000000..c9dd961
--- /dev/null
+++ b/hadamard_time.cpp
@@ -0,0 +1,34 @@
+#include "hadamard_mcmc.hpp"
+#include <iostream>
+
+class MeasureEnergy : public Measurement {
+ private:
+ unsigned N;
+ double totalE;
+
+ public:
+ MeasureEnergy() {
+ N = 0;
+ totalE = 0;
+ }
+
+ void after_sweep(double, double E, const Orthogonal&) override {
+ totalE += E;
+ N++;
+ }
+
+ double energy() const {
+ return totalE / N;
+ }
+};
+
+int main() {
+ MeasureEnergy A;
+ MCMC sim(20, 6.0, A);
+ sim.tune(1e4, 0.01);
+ sim.run(1e4);
+
+ std::cout << "The average energy was " << A.energy() << ".";
+
+ return 0;
+}
diff --git a/randutils b/randutils
new file mode 160000
+Subproject 8486a610a954a8248c12485fb4cfc390a5f5f85