From 4669226698d1b37ba05de8332751a5c1710fef95 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 19 Dec 2019 23:01:48 -0500 Subject: working pt with tuning --- hadamard.cpp | 158 +++++++++++------------------------------------------- hadamard_pt.hpp | 21 ++++---- hadamard_test.cpp | 6 ++- 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 #include +#include class MeasureEnergy : public Measurement { public: @@ -33,20 +34,17 @@ public: class MeasureTransitionRates : public ParallelMeasurement { public: - std::vector> nAccepted; + std::vector 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&) override { total_steps++; } @@ -54,14 +52,13 @@ public: int main(int argc, char* argv[]) { unsigned n_tuning = 1e2; - std::list β0s; - std::list β1s; - std::list 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 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 As(num); + std::vector 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 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> 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 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> ρ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 -#include -#include 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&){}; }; @@ -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 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 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 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"; } -- cgit v1.2.3-54-g00ecf