summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--hadamard.cpp158
-rw-r--r--hadamard_pt.hpp21
-rw-r--r--hadamard_test.cpp6
3 files changed, 44 insertions, 141 deletions
diff --git a/hadamard.cpp b/hadamard.cpp
index 1afb270..8eba184 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -3,6 +3,7 @@
#include <fstream>
#include <iostream>
+#include <chrono>
class MeasureEnergy : public Measurement {
public:
@@ -33,20 +34,17 @@ public:
class MeasureTransitionRates : public ParallelMeasurement {
public:
- std::vector<std::vector<unsigned>> nAccepted;
+ std::vector<unsigned> nAccepted;
unsigned total_steps;
- MeasureTransitionRates(unsigned n) : nAccepted(n - 1) {
+ MeasureTransitionRates(unsigned n) : nAccepted(n - 1, 0) {
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&,
+ void after_step(bool accepted, unsigned i, double, double, const MCMC&,
const MCMC&) override {
if (accepted)
- nAccepted[j - 1][i]++;
+ nAccepted[i]++;
}
void after_sweep(const std::vector<MCMC>&) override { total_steps++; }
@@ -54,14 +52,13 @@ public:
int main(int argc, char* argv[]) {
unsigned n_tuning = 1e2;
- std::list<double> β0s;
- std::list<double> β1s;
- std::list<unsigned> Ns;
+ double β₀ = 0.1;
+ double β₁ = 10;
+ unsigned N = 16;
unsigned k = 2;
double ε = 0.01;
unsigned M = 10;
- unsigned N = 1e4;
unsigned m = 1e4;
int opt;
@@ -72,16 +69,16 @@ int main(int argc, char* argv[]) {
k = atoi(optarg);
break;
case 'b':
- β0s.push_back(atof(optarg));
+ β₀ = atof(optarg);
break;
case 'c':
- β1s.push_back(atof(optarg));
+ β₁ = atof(optarg);
break;
case 'e':
ε = atof(optarg);
break;
case 'n':
- Ns.push_back((unsigned)atof(optarg));
+ m = (unsigned)atof(optarg);
break;
case 't':
n_tuning = (unsigned)atof(optarg);
@@ -99,33 +96,13 @@ int main(int argc, char* argv[]) {
unsigned n = pow(2, k);
- std::list<range> rs;
- unsigned num = 0;
-
- if ((β0s.size() != β1s.size()) || (β0s.size() != Ns.size())) {
- std::cout << "You need the same number of ranges!\n";
- exit(0);
- } else {
- auto it0 = β0s.begin();
- auto it1 = β1s.begin();
- auto itN = Ns.begin();
-
- while (it0 != β0s.end()) {
- rs.push_back({*it0, *it1, *itN});
- num += *itN;
- it0++;
- it1++;
- itN++;
- }
- }
-
- std::vector<Measurement*> As(num);
+ std::vector<Measurement*> As(N);
for (Measurement*& A : As) {
A = new MeasureEnergy();
}
- MeasureTransitionRates B(num);
+ MeasureTransitionRates B(N);
- PT p(rs, n, B, As);
+ PT p(β₀, β₁, N, n, B, As);
for (MCMC& sim : p.Ms) {
sim.M = walsh(k);
@@ -135,115 +112,42 @@ int main(int argc, char* argv[]) {
std::cout << "Beginning simulation of " << n << ".\n";
std::cout << "Beginning " << n_tuning << " tuning tempering updates of " << M
<< " sweeps each.\n";
- p.tune(n_tuning, M, ε);
- std::cout << "Finished tuning, beginning " << N << " measurement tempering updates of " << M
+ std::vector<double> f = p.tune(n_tuning, M, ε, 0.05);
+ std::cout << "Finished tuning, beginning " << m << " measurement tempering updates of " << M
<< " sweeps each.\n";
- p.run(N, M);
+ p.run(m, M);
std::cout << "Finished " << n << ".\n";
- std::string rs_string = "";
-
- for (range r : rs) {
- rs_string +=
- "_" + std::to_string(r.β0) + "_" + std::to_string(r.β1) + "_" + std::to_string(r.N);
- }
+ auto tag = std::chrono::high_resolution_clock::now();
- std::string filename = "probs_" + std::to_string(n) + rs_string + ".dat";
- std::ifstream file(filename);
+ std::string filename = "hmm_" + std::to_string(n) + "_" + std::to_string(β₀) + "_" + std::to_string(β₁) + "_" + std::to_string(N) + "_" + std::to_string(tag.time_since_epoch().count()) + ".dat";
- unsigned N_old = 0;
- std::vector<std::vector<unsigned long>> data_old(B.nAccepted.size());
+ std::ofstream file(filename);
- for (unsigned i = 0; i < B.nAccepted.size(); i++) {
- data_old[i].resize(B.nAccepted[i].size());
+ for (const MCMC& M : p.Ms) {
+ file << M.β << " ";
}
-
- if (file.is_open()) {
- file >> N_old;
-
- for (unsigned i = 0; i < B.nAccepted.size(); i++) {
- for (unsigned j = 0; j < B.nAccepted[i].size(); j++) {
- double num;
- file >> num;
- data_old[i][j] = num;
- }
- }
-
- file.close();
- }
-
- std::ofstream file_out(filename);
-
- file_out << N_old + B.total_steps << "\n";
+ file << std::endl;
for (unsigned i = 0; i < B.nAccepted.size(); i++) {
- for (unsigned j = 0; j < B.nAccepted[i].size(); j++) {
- file_out << std::fixed << data_old[i][j] + B.nAccepted[i][j] << " ";
- }
- file_out << "\n";
+ file << std::fixed << B.nAccepted[i] / (double)B.total_steps << " ";
}
-
- file_out.close();
-
- std::string efilename = "energies_" + std::to_string(n) + rs_string + ".dat";
- std::ifstream efile(efilename);
-
- unsigned Ne_old = 0;
- std::vector<double> edata_old(p.As.size());
-
- if (efile.is_open()) {
- efile >> Ne_old;
-
- for (unsigned i = 0; i < p.As.size(); i++) {
- double num;
- efile >> num;
- edata_old[i] = num;
- }
-
- efile.close();
- }
-
- std::ofstream efile_out(efilename);
-
- efile_out << Ne_old + ((MeasureEnergy*)As[0])->N << "\n";
-
- for (unsigned i = 0; i < As.size(); i++) {
- efile_out << std::fixed << edata_old[i] + ((MeasureEnergy*)As[i])->totalE << " ";
- }
-
- efile_out.close();
-
- std::string ρfilename = "rhos_" + std::to_string(n) + rs_string + ".dat";
- std::ifstream ρfile(ρfilename);
-
- std::vector<std::vector<unsigned>> ρdata_old(As.size());
+ file << std::endl;
for (unsigned i = 0; i < As.size(); i++) {
- ρdata_old[i].resize(((MeasureEnergy*)As[0])->ρ_dist.size());
- }
-
- if (ρfile.is_open()) {
- for (unsigned i = 0; i < As.size(); i++) {
- for (unsigned j = 0; j < ((MeasureEnergy*)As[0])->ρ_dist.size(); j++) {
- unsigned num;
- ρfile >> num;
- ρdata_old[i][j] = num;
- }
- }
-
- ρfile.close();
+ file << std::fixed << ((MeasureEnergy*)As[i])->totalE / ((MeasureEnergy*)As[i])->N << " ";
}
- std::ofstream ρfile_out(ρfilename);
+ file << std::endl;
for (unsigned i = 0; i < As.size(); i++) {
- for (unsigned j = 0; j < ((MeasureEnergy*)As[0])->ρ_dist.size(); j++) {
- ρfile_out << std::fixed << ρdata_old[i][j] + ((MeasureEnergy*)As[i])->ρ_dist[j] << " ";
+ for (unsigned j = 0; j < ((MeasureEnergy*)As[i])->ρ_dist.size(); j++) {
+ file << std::fixed << ((MeasureEnergy*)As[i])->ρ_dist[j] << " ";
}
- ρfile_out << "\n";
+ file << std::endl;
}
- ρfile_out.close();
+ file.close();
return 0;
}
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index f1cdf17..5500f57 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -2,8 +2,6 @@
#pragma once
#include "hadamard_mcmc.hpp"
#include <list>
-#include <fstream>
-#include <iostream>
void swap(MCMC& s1, MCMC& s2) {
std::swap(s1.M, s2.M);
@@ -12,7 +10,7 @@ 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_step(bool, unsigned, double, double, const MCMC&, const MCMC&){};
virtual void after_sweep(const std::vector<MCMC>&){};
};
@@ -54,7 +52,7 @@ public:
swap(Ms[i], Ms[j]);
if (!dry) {
- B.after_step(accepted, i, j, Δβ, ΔE, Ms[i], Ms[j]);
+ B.after_step(accepted, i, Δβ, ΔE, Ms[i], Ms[j]);
}
return accepted;
}
@@ -66,7 +64,7 @@ public:
}
}
- void tune(unsigned n0, unsigned m, double ε, double ε2) {
+ std::vector<double> tune(unsigned n0, unsigned m, double ε, double ε2) {
unsigned n = n0;
@@ -112,18 +110,12 @@ public:
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];
}
}
- 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]].β);
@@ -134,7 +126,7 @@ public:
double C = 0;
for (unsigned i = 0; i < Ms.size() - 1; i++) {
- η[η.size() - i - 1] = sqrt((f[i + 1] - f[i])) / (T(i) - T(i + 1));
+ η[η.size() - i - 1] = sqrt((f[i + 1] - f[i])) / (1 / Ms[i].β - 1 / Ms[i+1].β);
}
std::vector<double> T1(Ms.size() - 2);
@@ -156,9 +148,14 @@ public:
}
for (unsigned i = 0; i < T1.size(); i++) {
+ err += fabs(T1[i] - 1/ Ms[Ms.size() - i - 2].β) / T1[i];
Ms[Ms.size() - i - 2].β = 1 / T1[i];
}
+ if (err / T1.size() < ε2) {
+ return f;
+ }
+
n *= 2;
}
}
diff --git a/hadamard_test.cpp b/hadamard_test.cpp
index e505d79..a27b117 100644
--- a/hadamard_test.cpp
+++ b/hadamard_test.cpp
@@ -1,9 +1,10 @@
#include "hadamard_pt.hpp"
+#include <iostream>
int main() {
double β0 = 1.0;
double β1 = 10.0;
- unsigned N = 4;
+ unsigned N = 16;
unsigned n = 4;
ParallelMeasurement m1;
@@ -18,6 +19,7 @@ int main() {
a.tune(1000, 10, 0.1, 0.1);
for (MCMC& M : a.Ms) {
- std::cout << M.β;
+ std::cout << M.β << " ";
}
+ std::cout << "\n";
}