diff options
| -rw-r--r-- | hadamard.cpp | 10 | ||||
| -rw-r--r-- | hadamard_mcmc.hpp | 12 | ||||
| -rw-r--r-- | hadamard_pt.hpp | 47 |
3 files changed, 47 insertions, 22 deletions
diff --git a/hadamard.cpp b/hadamard.cpp index 1d99a2e..1afb270 100644 --- a/hadamard.cpp +++ b/hadamard.cpp @@ -92,9 +92,6 @@ int main(int argc, char* argv[]) { case 'M': M = (unsigned)atof(optarg); break; - case 'm': - m = (unsigned)atof(optarg); - break; default: exit(1); } @@ -135,11 +132,10 @@ int main(int argc, char* argv[]) { sim.E = sim.M.energy(); } - std::cout << "Beginning " << n_tuning << " tuning steps for " << n << ".\n"; - p.tune(n_tuning, ε); - std::cout << "Finished tuning, beginning " << N << " dry tempering updates of " << M + std::cout << "Beginning simulation of " << n << ".\n"; + std::cout << "Beginning " << n_tuning << " tuning tempering updates of " << M << " sweeps each.\n"; - p.run(m, M, true); + p.tune(n_tuning, M, ε); std::cout << "Finished tuning, beginning " << N << " measurement tempering updates of " << M << " sweeps each.\n"; p.run(N, M); diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp index 9f0d8c0..c9fa2da 100644 --- a/hadamard_mcmc.hpp +++ b/hadamard_mcmc.hpp @@ -140,6 +140,12 @@ public: virtual void after_sweep(double, double, const Orthogonal&){}; }; +typedef enum { + none, + up, + down +} color; + class MCMC { private: randutils::mt19937_rng rng; @@ -150,10 +156,14 @@ public: const double β; double E; Orthogonal M; + unsigned tag; + color c; - MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0) { + MCMC(unsigned n, double β0, Measurement& A, unsigned tag = 0) : 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 a3d4657..fba97b5 100644 --- a/hadamard_pt.hpp +++ b/hadamard_pt.hpp @@ -6,6 +6,8 @@ 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 { @@ -29,23 +31,41 @@ public: ParallelMeasurement& B; std::vector<Measurement*>& As; - PT(std::list<range> ranges, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As) + PT(double β₀, double β₁, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As) : B(B), As(As) { - unsigned count = 0; - for (range r : ranges) { - for (unsigned i = 1; i <= r.N; i++) { - double β = r.β0 + i * (r.β1 - r.β0) / r.N; - Ms.push_back(MCMC(n, β, *As[count])); - count++; - } + 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, double ε) { + 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 i = 0; i < Ms.size(); i++) { - Ms[i].tune(N, ε); + 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]++; + } + } } + + for ( } bool step(unsigned i, unsigned j, bool dry = false) { @@ -64,10 +84,9 @@ public: } void sweep(bool dry = false) { + for (unsigned i = 0; i < Ms.size() - 1; i++) { - for (unsigned j = i + 1; j < Ms.size(); j++) { - this->step(i, j, dry); - } + this->step(i, i + 1, dry); } } |