1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
#pragma once
#include "hadamard_mcmc.hpp"
#include <list>
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 {
public:
virtual void after_step(bool, unsigned, unsigned, double, double, const MCMC&, const MCMC&){};
virtual void after_sweep(const std::vector<MCMC>&){};
};
typedef struct range {
double β0;
double β1;
unsigned N;
} range;
class PT {
private:
randutils::mt19937_rng rng;
public:
std::vector<MCMC> Ms;
ParallelMeasurement& B;
std::vector<Measurement*>& As;
PT(double β₀, double β₁, unsigned n, ParallelMeasurement& B, std::vector<Measurement*>& As)
: B(B), As(As) {
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, 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 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) {
double Δβ = Ms[i].β - Ms[j].β;
double ΔE = Ms[i].E - Ms[j].E;
bool accepted = Δβ * ΔE > 0 || exp(Δβ * ΔE) > rng.uniform((double)0.0, 1.0);
if (accepted)
swap(Ms[i], Ms[j]);
if (!dry) {
B.after_step(accepted, i, j, Δβ, ΔE, Ms[i], Ms[j]);
}
return accepted;
}
void sweep(bool dry = false) {
for (unsigned i = 0; i < Ms.size() - 1; i++) {
this->step(i, i + 1, dry);
}
}
void run(unsigned n, unsigned m, bool dry = false) {
for (unsigned i = 0; i < n; i++) {
#pragma omp parallel for
for (unsigned j = 0; j < Ms.size(); j++) {
Ms[j].run(m, dry);
}
this->sweep(dry);
if (!dry) {
B.after_sweep(this->Ms);
}
}
}
};
|
