diff options
-rw-r--r-- | hadamard_mcmc.hpp | 17 |
1 files changed, 8 insertions, 9 deletions
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp index 218f8d1..e64eacf 100644 --- a/hadamard_mcmc.hpp +++ b/hadamard_mcmc.hpp @@ -53,18 +53,16 @@ private: unsigned axis_2; double Δθ; - std::vector<double> rows; - public: Givens(Orthogonal& m, bool t, unsigned a1, unsigned a2, double θ0, Rng& rng) - : m(m), rows(2 * m.size()) { + : m(m) { transpose = t; axis_1 = a1; axis_2 = a2; Δθ = rng.uniform(-θ0, θ0); } - Givens(Orthogonal& m, double θ0, Rng& rng) : m(m), rows(2 * m.size()) { + Givens(Orthogonal& m, double θ0, Rng& rng) : m(m) { Δθ = rng.uniform(-θ0, θ0); unsigned axis1axis2 = rng.uniform((unsigned)0, m.size() * (m.size() - 1) - 1); @@ -76,7 +74,7 @@ public: } } - double tryRotation() { + double tryRotation(std::vector<double>& rows) const { double ΔE = 0; double c = cos(Δθ); double s = sin(Δθ); @@ -106,7 +104,7 @@ public: return ΔE; } - void acceptRotation() const { + void acceptRotation(std::vector<double>& rows) const { for (unsigned i = 0; i < m.size(); i++) { if (transpose) { m(i, axis_1) = rows[i]; @@ -133,6 +131,7 @@ typedef enum { class MCMC { private: + std::vector<double> row_storage; Measurement& A; double θ0; @@ -142,19 +141,19 @@ public: double E; Orthogonal M; - MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0) { + MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0), row_storage(2 * n) { θ0 = M_PI; E = M.energy(); } bool step(Givens& g, bool dry = false) { - double ΔE = g.tryRotation(); + double ΔE = g.tryRotation(row_storage); bool accepted = ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0); if (accepted) { E += ΔE; - g.acceptRotation(); + g.acceptRotation(row_storage); } if (!dry) { |