summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--hadamard_pt.hpp155
-rw-r--r--hadamard_test.cpp8
2 files changed, 82 insertions, 81 deletions
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index 33dff71..f1cdf17 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -22,7 +22,6 @@ typedef struct range {
unsigned N;
} range;
-
class PT {
private:
randutils::mt19937_rng rng;
@@ -41,6 +40,10 @@ public:
}
}
+ double T(unsigned i) {
+ return 1 / Ms[Ms.size() - i - 1].β;
+ }
+
bool step(unsigned i, unsigned j, bool dry = false) {
double Δβ = Ms[i].β - Ms[j].β;
double ΔE = Ms[i].E - Ms[j].E;
@@ -63,105 +66,101 @@ 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);
+ void tune(unsigned n0, unsigned m, double ε, double ε2) {
- colors.front() = down;
- colors.back() = up;
+ unsigned n = n0;
- for (unsigned i = 0; i < n; i++) {
-#pragma omp parallel for
- for (unsigned j = 0; j < Ms.size(); j++) {
- Ms[j].tune(m, ε);
- }
+ while (true) {
+ std::vector<color> colors(Ms.size(), none);
+ std::vector<unsigned> nu(Ms.size(), 0);
+ std::vector<unsigned> nd(Ms.size(), 0);
- 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;
+ 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, ε);
+ }
- 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]++;
+ 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];
+ std::vector<double> f(Ms.size());
+
+ double err = 0;
+ 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]);
+ double δerr = fabs(f[i] - (1 - (Ms.size() - i - 1) / (Ms.size() - 1.0)));
+ err += δerr;
+ if (δerr < δ && 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::cout << err << "\n";
+ if (err < ε2) {
+ break;
}
- }
+ 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;
- 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() - 1; i++) {
+ η[η.size() - i - 1] = sqrt((f[i + 1] - f[i])) / (T(i) - T(i + 1));
+ }
- 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().β;
+ std::vector<double> T1(Ms.size() - 2);
- for (unsigned i = 0; i < η.size(); i++) {
- C += η[i] * (T0[i + 1] - T0[i]);
- }
+ for (unsigned i = 0; i < η.size(); i++) {
+ C += η[i] * (T(i + 1) - T(i));
+ }
- double x= 0;
- unsigned j = 1;
+ 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++;
+ for (unsigned i = 0; i < η.size(); i++) {
+ double xnew = x + η[i] * (T(i + 1) - T(i)) / C;
+ while (j < xnew * η.size() && j < η.size()) {
+ T1[j - 1] = T(i) + (j / (double)η.size() - x) / η[i] * C;
+ j++;
+ }
+ x = xnew;
}
- x = xnew;
- }
- for (unsigned i = 0; i < T1.size(); i++) {
- Ms[Ms.size() - i - 1].β = 1 / T1[i];
- }
+ for (unsigned i = 0; i < T1.size(); i++) {
+ Ms[Ms.size() - i - 2].β = 1 / T1[i];
+ }
- for (double ff : f) {
- outfile << ff << " ";
+ n *= 2;
}
- outfile << "\n";
- outfile.close();
-
}
void run(unsigned n, unsigned m, bool dry = false) {
diff --git a/hadamard_test.cpp b/hadamard_test.cpp
index c1dd8cd..e505d79 100644
--- a/hadamard_test.cpp
+++ b/hadamard_test.cpp
@@ -3,7 +3,7 @@
int main() {
double β0 = 1.0;
double β1 = 10.0;
- unsigned N = 8;
+ unsigned N = 4;
unsigned n = 4;
ParallelMeasurement m1;
@@ -15,7 +15,9 @@ int main() {
PT a(β0, β1, N, n, m1, m3);
- while (true) {
- a.tune(1000, 10, 0.1);
+ a.tune(1000, 10, 0.1, 0.1);
+
+ for (MCMC& M : a.Ms) {
+ std::cout << M.β;
}
}