diff options
| -rw-r--r-- | hadamard_pt.hpp | 155 | ||||
| -rw-r--r-- | hadamard_test.cpp | 8 |
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.β; } } |