From 84edbcb1b0dd1c5ac5c1626e7dc8cf22bbf10139 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 14 Oct 2019 20:31:35 -0400 Subject: initialized repo --- .gitignore | 1 + .gitmodules | 3 + hadamard.cpp | 174 +++++++++++++++++++++++++++++++++++++++++++++++++++ hadamard_mcmc.hpp | 181 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ hadamard_pt.hpp | 74 ++++++++++++++++++++++ hadamard_time.cpp | 34 ++++++++++ randutils | 1 + 7 files changed, 468 insertions(+) create mode 100644 .gitignore create mode 100644 .gitmodules create mode 100644 hadamard.cpp create mode 100644 hadamard_mcmc.hpp create mode 100644 hadamard_pt.hpp create mode 100644 hadamard_time.cpp create mode 160000 randutils 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 +#include + +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> 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&) 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 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> 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 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 + +inline double Ei(double x) { return -fabs(x); } + +class Orthogonal { +private: + unsigned d; + std::vector 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 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&) {}; +}; + +class PT { +private: + randutils::mt19937_rng rng; + +public: + std::vector Ms; + ParallelMeasurement& B; + std::vector& As; + + PT(double β0, double β1, unsigned N, unsigned n, ParallelMeasurement& B, std::vector& 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 + +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 index 0000000..8486a61 --- /dev/null +++ b/randutils @@ -0,0 +1 @@ +Subproject commit 8486a610a954a8248c12485fb4cfc390a5f5f854 -- cgit v1.2.3-70-g09d2