diff options
Diffstat (limited to 'hadamard_pt.hpp')
-rw-r--r-- | hadamard_pt.hpp | 74 |
1 files changed, 74 insertions, 0 deletions
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp new file mode 100644 index 0000000..4006e0e --- /dev/null +++ b/hadamard_pt.hpp @@ -0,0 +1,74 @@ + +#pragma once +#include "hadamard_mcmc.hpp" + +void swap(MCMC& s1, MCMC& s2) { + std::swap(s1.M, s2.M); + std::swap(s1.E, s2.E); +} + +class ParallelMeasurement { +public: + virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&) {}; + virtual void after_sweep(const std::vector<MCMC>&) {}; +}; + +class PT { +private: + randutils::mt19937_rng rng; + +public: + std::vector<MCMC> Ms; + ParallelMeasurement& B; + std::vector<Measurement*>& As; + + PT(double β0, double β1, unsigned N, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As) + : B(B), As(As) { + Ms.reserve(N); + for (unsigned i = 0; i < N; i++) { + double β = β0 + i * (β1 - β0) / (N - 1); + Ms.push_back(MCMC(n, β, *As[i])); + } + } + + void tune(unsigned N, double ε) { +#pragma omp parallel for + for (unsigned i = 0; i < Ms.size(); i++) { + Ms[i].tune(N, ε); + } + } + + bool step(unsigned i, unsigned j) { + double Δβ = Ms[i].β - Ms[j].β; + double ΔE = Ms[i].E - Ms[j].E; + + bool accepted = Δβ * ΔE > 0 || exp(Δβ * ΔE) > rng.uniform((double)0.0, 1.0); + + if (accepted) { + swap(Ms[i], Ms[j]); + } + + B.after_step(accepted, i, j, Δβ, ΔE, Ms[i], Ms[j]); + + return accepted; + } + + void sweep() { + for (unsigned i = 0; i < Ms.size() - 1; i++) { + for (unsigned j = i + 1; j < Ms.size(); j++) { + this->step(i, j); + } + } + } + + void run(unsigned n, unsigned m) { + for (unsigned i = 0; i < n; i++) { +#pragma omp parallel for + for (unsigned j = 0; j < Ms.size(); j++) { + Ms[j].run(m); + } + this->sweep(); + } + } +}; + |