summaryrefslogtreecommitdiff
path: root/hadamard_pt.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'hadamard_pt.hpp')
-rw-r--r--hadamard_pt.hpp144
1 files changed, 109 insertions, 35 deletions
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