summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--hadamard.cpp10
-rw-r--r--hadamard_mcmc.hpp12
-rw-r--r--hadamard_pt.hpp47
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);
}
}