From 2cc44d8981e42bfe663e8e7dee8dbb7914a9482f Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 20 Dec 2019 17:37:01 -0500
Subject: fixed it
---
hadamard.cpp | 7 ++++++-
hadamard_mcmc.hpp | 2 +-
hadamard_pt.hpp | 35 +++++++++++++++++++----------------
3 files changed, 26 insertions(+), 18 deletions(-)
diff --git a/hadamard.cpp b/hadamard.cpp
index 5154677..2f45dcc 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -116,7 +116,7 @@ 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, ε, ε2);
+ std::vector<double> f = p.tune(n_tuning, M, ε, ε2);
std::cout << "Finished tuning, beginning " << m << " measurement tempering updates of " << M
<< " sweeps each.\n";
p.run(m, M);
@@ -133,6 +133,11 @@ int main(int argc, char* argv[]) {
}
file << std::endl;
+ for (double ff : f) {
+ file << ff << " ";
+ }
+ file << std::endl;
+
for (unsigned i = 0; i < B.nAccepted.size(); i++) {
file << std::fixed << B.nAccepted[i] / (double)B.total_steps << " ";
}
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index 145454b..eab9916 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -148,11 +148,11 @@ typedef enum {
class MCMC {
private:
- randutils::mt19937_rng rng;
Measurement& A;
double θ0;
public:
+ randutils::mt19937_rng rng;
double β;
double E;
Orthogonal M;
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp
index 4cde211..de6590d 100644
--- a/hadamard_pt.hpp
+++ b/hadamard_pt.hpp
@@ -2,7 +2,6 @@
#pragma once
#include "hadamard_mcmc.hpp"
#include <list>
-#include <iostream>
void swap(MCMC& s1, MCMC& s2) {
std::swap(s1.M, s2.M);
@@ -65,15 +64,16 @@ public:
}
}
- void tune(unsigned n0, unsigned m, double ε, double ε2) {
+ std::vector<double> tune(unsigned n0, unsigned m, double ε, double ε2) {
unsigned n = n0;
while (true) {
std::vector<color> colors(Ms.size(), none);
+ colors.front() = down;
+ colors.back() = up;
+
std::vector<unsigned> nu(Ms.size(), 0);
std::vector<unsigned> nd(Ms.size(), 0);
- colors.front() = down;
- colors.back() = up;
for (unsigned i = 0; i < n; i++) {
#pragma omp parallel for
@@ -81,23 +81,25 @@ public:
Ms[j].tune(m, ε);
}
- for (unsigned j = 0; j < Ms.size() - 1; j++) {
+ for (unsigned k = 0; k < m * Ms.size() - 1; k++) {
+ unsigned j = Ms[0].rng.uniform((unsigned)0, (unsigned)(Ms.size() - 2));
+
if (this->step(j, j + 1, true)) {
std::swap(colors[j], colors[j + 1]);
colors.front() = down;
colors.back() = up;
}
- if (i > n / 2) {
- for (unsigned j = 0; j < Ms.size(); j++) {
- if (colors[j] == up) {
- nu[j]++;
- } else if (colors[j] == down) {
- nd[j]++;
- }
- }
+ }
+
+ for (unsigned j = 0; j < Ms.size(); j++) {
+ if (colors[j] == up) {
+ nu[j]++;
+ } else if (colors[j] == down) {
+ nd[j]++;
}
}
+
}
std::vector<double> f(Ms.size());
@@ -129,10 +131,11 @@ public:
C += η[i] * (T(i + 1) - T(i));
}
- double x= 0;
- unsigned j = 1;
std::vector<double> T1(Ms.size() - 2);
+ double x = 0;
+ unsigned j = 1;
+
for (unsigned i = 0; i < η.size(); i++) {
double xnew = x + η[i] * (T(i + 1) - T(i)) / C;
while (j < xnew * η.size() && j < η.size()) {
@@ -149,7 +152,7 @@ public:
}
if (err / T1.size() * Ms.size() / (1/Ms.front().β - 1/Ms.back().β) < ε2) {
- break;
+ return f;
}
n *= 2;
--
cgit v1.2.3-70-g09d2