summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--hadamard.cpp138
-rw-r--r--hadamard_mcmc.hpp71
-rw-r--r--hadamard_pt.hpp12
3 files changed, 129 insertions, 92 deletions
diff --git a/hadamard.cpp b/hadamard.cpp
index 6b9a64a..15f73b1 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -1,57 +1,54 @@
#include "hadamard_pt.hpp"
-#include <iostream>
#include <fstream>
+#include <iostream>
class MeasureEnergy : public Measurement {
- public:
- unsigned N;
- double totalE;
- MeasureEnergy() {
- N = 0;
- totalE = 0;
- }
+public:
+ unsigned N;
+ double totalE;
- void after_sweep(double, double E, const Orthogonal&) override {
- totalE += E;
- N++;
- }
+ MeasureEnergy() {
+ N = 0;
+ totalE = 0;
+ }
- double energy() const {
- return totalE / N;
- }
+ void after_sweep(double, double E, const Orthogonal&) override {
+ totalE += E;
+ N++;
+ }
+
+ double energy() const { return totalE / N; }
};
class MeasureTransitionRates : public ParallelMeasurement {
- public:
- std::vector<std::vector<unsigned>> nAccepted;
- unsigned total_steps;
-
- MeasureTransitionRates(unsigned n) : nAccepted(n - 1) {
- total_steps = 0;
- for (unsigned i = 0; i < n - 1; i++) {
- nAccepted[i].resize(i + 1);
- }
+public:
+ std::vector<std::vector<unsigned>> nAccepted;
+ unsigned total_steps;
+
+ MeasureTransitionRates(unsigned n) : nAccepted(n - 1) {
+ total_steps = 0;
+ for (unsigned i = 0; i < n - 1; i++) {
+ nAccepted[i].resize(i + 1);
}
+ }
- void after_step(bool accepted, unsigned i, unsigned j, double, double, const MCMC&, const MCMC&) override {
- if (accepted) {
- nAccepted[j - 1][i]++;
- }
- }
+ void after_step(bool accepted, unsigned i, unsigned j, double, double, const MCMC&,
+ const MCMC&) override {
+ if (accepted)
+ nAccepted[j - 1][i]++;
+ }
- void after_sweep(const std::vector<MCMC>&) override {
- total_steps++;
- }
+ void after_sweep(const std::vector<MCMC>&) override { total_steps++; }
};
int main(int argc, char* argv[]) {
- unsigned trials = 1e6;
+ unsigned n_tuning = 1e2;
double β0 = 1.00;
double βf = 40.00;
unsigned num = 40;
- unsigned n = 20;
+ unsigned k = 2;
double ε = 0.01;
unsigned M = 10;
@@ -59,37 +56,39 @@ int main(int argc, char* argv[]) {
int opt;
- while ((opt = getopt(argc, argv, "d:b:c:n:t:N:M:e:")) != -1) {
+ while ((opt = getopt(argc, argv, "k:b:c:n:t:N:M:e:")) != -1) {
switch (opt) {
- case 'd':
- n = atoi(optarg);
- break;
- case 'b':
- β0 = atof(optarg);
- break;
- case 'c':
- βf = atof(optarg);
- break;
- case 'e':
- ε = atof(optarg);
- break;
- case 'n':
- num = atof(optarg);
- break;
- case 't':
- trials = (unsigned)atof(optarg);
- break;
- case 'N':
- N = (unsigned)atof(optarg);
- break;
- case 'M':
- M = (unsigned)atof(optarg);
- break;
- default:
- exit(1);
+ case 'k':
+ k = atoi(optarg);
+ break;
+ case 'b':
+ β0 = atof(optarg);
+ break;
+ case 'c':
+ βf = atof(optarg);
+ break;
+ case 'e':
+ ε = atof(optarg);
+ break;
+ case 'n':
+ num = atof(optarg);
+ break;
+ case 't':
+ n_tuning = (unsigned)atof(optarg);
+ break;
+ case 'N':
+ N = (unsigned)atof(optarg);
+ break;
+ case 'M':
+ M = (unsigned)atof(optarg);
+ break;
+ default:
+ exit(1);
}
}
+ unsigned n = pow(2, k);
+
std::vector<Measurement*> As(num);
for (Measurement*& A : As) {
A = new MeasureEnergy();
@@ -98,13 +97,19 @@ int main(int argc, char* argv[]) {
PT p(β0, βf, num, n, B, As);
- std::cout << "Beginning " << trials << " tuning steps for " << n << ".\n";
- p.tune(trials, ε);
- std::cout << "Finished tuning, beginning " << N << " tempering updates of " << M << " sweeps each.\n";
+ for (MCMC& sim : p.Ms) {
+ sim.M = walsh(k);
+ }
+
+ std::cout << "Beginning " << n_tuning << " tuning steps for " << n << ".\n";
+ p.tune(n_tuning, ε);
+ std::cout << "Finished tuning, beginning " << N << " tempering updates of " << M
+ << " sweeps each.\n";
p.run(N, M);
std::cout << "Finished " << n << ".\n";
- std::string filename = "probs_" + std::to_string(n) + "_" + std::to_string(β0) + "_" + std::to_string(βf) + "_" + std::to_string(num) + ".dat";
+ std::string filename = "probs_" + std::to_string(n) + "_" + std::to_string(β0) + "_" +
+ std::to_string(βf) + "_" + std::to_string(num) + ".dat";
std::ifstream file(filename);
unsigned N_old = 0;
@@ -141,7 +146,8 @@ int main(int argc, char* argv[]) {
file_out.close();
- std::string efilename = "energies_" + std::to_string(n) + "_" + std::to_string(β0) + "_" + std::to_string(βf) + "_" + std::to_string(num) + ".dat";
+ std::string efilename = "energies_" + std::to_string(n) + "_" + std::to_string(β0) + "_" +
+ std::to_string(βf) + "_" + std::to_string(num) + ".dat";
std::ifstream efile(efilename);
unsigned Ne_old = 0;
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index 4647396..5ef062f 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -37,6 +37,28 @@ public:
}
};
+Orthogonal walsh(unsigned k) {
+ if (k == 0) {
+ return Orthogonal(1);
+ } else {
+ Orthogonal s = walsh(k - 1);
+ Orthogonal t = Orthogonal(2 * s.size());
+
+ for (unsigned i = 0; i < s.size(); i++) {
+ for (unsigned j = 0; j < s.size(); j++) {
+ double sij = s(i, j);
+
+ t(i, j) = sij;
+ t(s.size() + i, j) = sij;
+ t(i, s.size() + j) = sij;
+ t(s.size() + i, s.size() + j) = -sij;
+ }
+ }
+
+ return t;
+ }
+}
+
class Givens {
private:
Orthogonal& m;
@@ -137,36 +159,45 @@ public:
bool step(Givens& g) {
double ΔE = g.tryRotation();
- if (ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0)) {
+ bool accepted = ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0);
+
+ if (accepted) {
E += ΔE;
g.acceptRotation();
- A.after_step(true, g, θ0, E, ΔE, M);
- return true;
- } else {
- return false;
- A.after_step(false, g, θ0, E, ΔE, M);
}
- }
- void tune(unsigned N, double ε) {
- for (unsigned i = 0; i < N; i++) {
- Givens g(M, θ0, rng);
- bool stepAccepted = this->step(g);
- if (stepAccepted) {
- θ0 *= 1 + ε;
- } else {
- θ0 /= 1 + ε;
- }
- }
+ A.after_step(accepted, g, θ0, E, ΔE, M);
+ return accepted;
}
- void sweep() {
+ double sweep() {
+ unsigned total = 0;
+ unsigned accepted = 0;
+
for (unsigned i = 0; i < M.size() - 1; i++) {
for (unsigned j = i + 1; j < M.size(); j++) {
Givens g1(M, false, i, j, θ0, rng);
- this->step(g1);
Givens g2(M, true, i, j, θ0, rng);
- this->step(g2);
+
+ if (this->step(g1))
+ accepted++;
+ if (this->step(g2))
+ accepted++;
+
+ total += 2;
+ }
+ }
+
+ return (double)accepted / (double)total;
+ }
+
+ void tune(unsigned N, double ε) {
+ for (unsigned i = 0; i < N; i++) {
+ double ratio_accepted = this->sweep();
+ if (ratio_accepted > 0.5) {
+ θ0 *= 1 + ε;
+ } else {
+ θ0 /= 1 + ε;
}
}
}
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index 4006e0e..404129c 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -9,8 +9,8 @@ void swap(MCMC& s1, MCMC& s2) {
class ParallelMeasurement {
public:
- virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&) {};
- virtual void after_sweep(const std::vector<MCMC>&) {};
+ virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&){};
+ virtual void after_sweep(const std::vector<MCMC>&){};
};
class PT {
@@ -22,7 +22,8 @@ public:
ParallelMeasurement& B;
std::vector<Measurement*>& As;
- PT(double β0, double β1, unsigned N, unsigned n, 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++) {
@@ -44,12 +45,10 @@ public:
bool accepted = Δβ * ΔE > 0 || exp(Δβ * ΔE) > rng.uniform((double)0.0, 1.0);
- if (accepted) {
+ if (accepted)
swap(Ms[i], Ms[j]);
- }
B.after_step(accepted, i, j, Δβ, ΔE, Ms[i], Ms[j]);
-
return accepted;
}
@@ -68,6 +67,7 @@ public:
Ms[j].run(m);
}
this->sweep();
+ B.after_sweep(this->Ms);
}
}
};