diff options
| -rw-r--r-- | hadamard_mcmc.hpp | 8 | ||||
| -rw-r--r-- | hadamard_pt.hpp | 144 | ||||
| -rw-r--r-- | hadamard_test.cpp | 21 |
3 files changed, 132 insertions, 41 deletions
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp index c9fa2da..145454b 100644 --- a/hadamard_mcmc.hpp +++ b/hadamard_mcmc.hpp @@ -153,17 +153,13 @@ private: double θ0; public: - const double β; + double β; double E; Orthogonal M; - unsigned tag; - color c; - MCMC(unsigned n, double β0, Measurement& A, unsigned tag = 0) : A(A), M(n), β(β0) { + MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0) { θ0 = M_PI; E = M.energy(); - tag = tag; - c = none; } bool step(Givens& g, bool dry = false) { diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp index fba97b5..33dff71 100644 --- a/hadamard_pt.hpp +++ b/hadamard_pt.hpp @@ -2,12 +2,12 @@ #pragma once #include "hadamard_mcmc.hpp" #include <list> +#include <fstream> +#include <iostream> void swap(MCMC& s1, MCMC& s2) { std::swap(s1.M, s2.M); std::swap(s1.E, s2.E); - std::swap(s1.tag, s2.tag); - std::swap(s1.c, s2.c); } class ParallelMeasurement { @@ -22,6 +22,7 @@ typedef struct range { unsigned N; } range; + class PT { private: randutils::mt19937_rng rng; @@ -31,41 +32,13 @@ public: ParallelMeasurement& B; std::vector<Measurement*>& As; - PT(double β₀, double β₁, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As) + PT(double β₀, double β₁, unsigned N, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As) : B(B), As(As) { - Ms.reserve(n); - for (unsigned i = 1; i <= n; i++) { - double β = β₀ + i * (β₁ - β₀) / n; - Ms.push_back(MCMC(n, β, *As[i - 1], i - 1)); - } - Ms[0].c = down; - Ms[n - 1].c = up; - } - - void tune(unsigned n, unsigned m, double ε) { - std::vector<unsigned> nu(Ms.size()); - std::vector<unsigned> nd(Ms.size()); - - for (unsigned i = 0; i < n; i++) { - Ms[0].c = down; - Ms[Ms.size() - 1].c = up; - -#pragma omp parallel for - for (unsigned j = 0; j < Ms.size(); j++) { - Ms[j].tune(m, ε); - } - this->sweep(true); - - for (unsigned j = 0; j < Ms.size(); j++) { - if (Ms[j].c == up) { - nu[j]++; - } else if (Ms[j].c == down) { - nd[j]++; - } - } + Ms.reserve(N); + for (unsigned i = 0; i < N; i++) { + double β = β₀ + i * (β₁ - β₀) / (N - 1); + Ms.push_back(MCMC(n, β, *As[i])); } - - for ( } bool step(unsigned i, unsigned j, bool dry = false) { @@ -90,6 +63,107 @@ public: } } + void tune(unsigned n, unsigned m, double ε) { + std::vector<color> colors(Ms.size(), none); + std::vector<unsigned> nu(Ms.size(), 0); + std::vector<unsigned> nd(Ms.size(), 0); + + colors.front() = down; + colors.back() = up; + + for (unsigned i = 0; i < n; i++) { +#pragma omp parallel for + for (unsigned j = 0; j < Ms.size(); j++) { + Ms[j].tune(m, ε); + } + + for (unsigned j = 0; j < Ms.size() - 1; j++) { + if (this->step(j, j + 1, true)) { + std::swap(colors[j], colors[j + 1]); + colors.front() = down; + colors.back() = up; + + if (colors[j] == up) { + nu[j]++; + } else if (colors[j] == down) { + nd[j]++; + } + if (colors[j + 1] == up) { + nu[j + 1]++; + } else if (colors[j + 1] == down) { + nd[j + 1]++; + } + } + } + } + + std::ofstream outfile("test.dat", std::ios::app); + + std::vector<double> f(Ms.size()); + + double δ = 0.5; + std::vector<unsigned> f_keep; + double f_last = 0; + for (unsigned i = 0; i < Ms.size(); i++) { + f[i] = nu[i] / (double)(nu[i] + nd[i]); + if (fabs(f[i] - (1 - (Ms.size() - i - 1) / (Ms.size() - 1.0))) < δ && f[i] >= f_last) { + f_keep.push_back(i); + f_last = f[i]; + } + } + + for (unsigned i = 0; i < f_keep.size() - 1; i++) { + for (unsigned j = f_keep[i]; j < f_keep[i + 1]; j++) { + f[j] = f[f_keep[i]] + (f[f_keep[i + 1]] - f[f_keep[i]]) / (Ms[f_keep[i+1]].β - Ms[f_keep[i]].β) * (Ms[j].β - Ms[f_keep[i]].β); + } + } + + + + std::vector<double> η(Ms.size() - 1); + double C = 0; + + for (unsigned i = 0; i < Ms.size() - 1; i++) { + η[η.size() - i - 1] = sqrt((f[i + 1] - f[i]) / pow((1 / Ms[i].β - 1 / Ms[i + 1].β), 2)); + } + + std::vector<double> T0(Ms.size()); + std::vector<double> T1(Ms.size()); + + for (unsigned i = 0; i < Ms.size(); i++) { + T0[i] = 1 / Ms[Ms.size() - i - 1].β; + } + T1.front() = 1 / Ms.back().β; + T1.back() = 1 / Ms.front().β; + + for (unsigned i = 0; i < η.size(); i++) { + C += η[i] * (T0[i + 1] - T0[i]); + } + + double x= 0; + unsigned j = 1; + + for (unsigned i = 0; i < η.size(); i++) { + double xnew = x + η[i] * (T0[i + 1] - T0[i]) / C; + while (j < xnew * η.size()) { + T1[j] = T0[i] + (j / (double)η.size() - x) / η[i] * C; + j++; + } + x = xnew; + } + + for (unsigned i = 0; i < T1.size(); i++) { + Ms[Ms.size() - i - 1].β = 1 / T1[i]; + } + + for (double ff : f) { + outfile << ff << " "; + } + outfile << "\n"; + outfile.close(); + + } + void run(unsigned n, unsigned m, bool dry = false) { for (unsigned i = 0; i < n; i++) { #pragma omp parallel for diff --git a/hadamard_test.cpp b/hadamard_test.cpp new file mode 100644 index 0000000..c1dd8cd --- /dev/null +++ b/hadamard_test.cpp @@ -0,0 +1,21 @@ +#include "hadamard_pt.hpp" + +int main() { + double β0 = 1.0; + double β1 = 10.0; + unsigned N = 8; + unsigned n = 4; + + ParallelMeasurement m1; + std::vector<Measurement> m2(N); + std::vector<Measurement*> m3; + for (Measurement& m : m2) { + m3.push_back(&m); + } + + PT a(β0, β1, N, n, m1, m3); + + while (true) { + a.tune(1000, 10, 0.1); + } +} |