diff options
Diffstat (limited to 'hadamard_mcmc.hpp')
-rw-r--r-- | hadamard_mcmc.hpp | 71 |
1 files changed, 51 insertions, 20 deletions
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp index 4647396..5ef062f 100644 --- a/hadamard_mcmc.hpp +++ b/hadamard_mcmc.hpp @@ -37,6 +37,28 @@ public: } }; +Orthogonal walsh(unsigned k) { + if (k == 0) { + return Orthogonal(1); + } else { + Orthogonal s = walsh(k - 1); + Orthogonal t = Orthogonal(2 * s.size()); + + for (unsigned i = 0; i < s.size(); i++) { + for (unsigned j = 0; j < s.size(); j++) { + double sij = s(i, j); + + t(i, j) = sij; + t(s.size() + i, j) = sij; + t(i, s.size() + j) = sij; + t(s.size() + i, s.size() + j) = -sij; + } + } + + return t; + } +} + class Givens { private: Orthogonal& m; @@ -137,36 +159,45 @@ public: bool step(Givens& g) { double ΔE = g.tryRotation(); - if (ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0)) { + bool accepted = ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0); + + if (accepted) { E += ΔE; g.acceptRotation(); - A.after_step(true, g, θ0, E, ΔE, M); - return true; - } else { - return false; - A.after_step(false, g, θ0, E, ΔE, M); } - } - void tune(unsigned N, double ε) { - for (unsigned i = 0; i < N; i++) { - Givens g(M, θ0, rng); - bool stepAccepted = this->step(g); - if (stepAccepted) { - θ0 *= 1 + ε; - } else { - θ0 /= 1 + ε; - } - } + A.after_step(accepted, g, θ0, E, ΔE, M); + return accepted; } - void sweep() { + double sweep() { + unsigned total = 0; + unsigned accepted = 0; + for (unsigned i = 0; i < M.size() - 1; i++) { for (unsigned j = i + 1; j < M.size(); j++) { Givens g1(M, false, i, j, θ0, rng); - this->step(g1); Givens g2(M, true, i, j, θ0, rng); - this->step(g2); + + if (this->step(g1)) + accepted++; + if (this->step(g2)) + accepted++; + + total += 2; + } + } + + return (double)accepted / (double)total; + } + + void tune(unsigned N, double ε) { + for (unsigned i = 0; i < N; i++) { + double ratio_accepted = this->sweep(); + if (ratio_accepted > 0.5) { + θ0 *= 1 + ε; + } else { + θ0 /= 1 + ε; } } } |