From dbee70fe5eb7abca31c01c046ec2174f2dde7665 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 30 Jan 2020 13:52:53 -0500 Subject: any orthogonal size that is a multiple of four is now available instead of powers of two --- hadamard.cpp | 21 +++++++-- hadamard_mcmc.hpp | 26 ++-------- matrices.hpp | 139 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 161 insertions(+), 25 deletions(-) create mode 100644 matrices.hpp diff --git a/hadamard.cpp b/hadamard.cpp index d9e5bfd..86c27d8 100644 --- a/hadamard.cpp +++ b/hadamard.cpp @@ -1,5 +1,6 @@ #include "hadamard_pt.hpp" +#include "matrices.hpp" #include #include @@ -9,6 +10,7 @@ class MeasureEnergy : public Measurement { public: unsigned N; double totalE; + double totalE2; unsigned n; std::vector ρ_dist; @@ -16,11 +18,13 @@ public: n = n_bins; N = 0; totalE = 0; + totalE2 = 0; } void after_sweep(double, double E, const Orthogonal& M) override { N++; totalE += E; + totalE2 += pow(E, 2); double max = sqrt(M.size()); for (unsigned i = 0; i < M.size(); i++) { for (unsigned j = 0; j < M.size(); j++) { @@ -30,6 +34,7 @@ public: } double energy() const { return totalE / N; } + double specific_heat() const { return totalE2 / N - pow(totalE / N, 2); } }; class MeasureTransitionRates : public ParallelMeasurement { @@ -65,6 +70,10 @@ int main(int argc, char* argv[]) { switch (opt) { case 'k': k = atoi(optarg); + if (k == 0 || k > 8) { + std::cout << "The size k must be an integer from 1 to 8!" << std::endl; + exit(1); + } break; case 'b': β₀ = atof(optarg); @@ -95,7 +104,7 @@ int main(int argc, char* argv[]) { } } - unsigned n = pow(2, k); + unsigned n = 4 * k; std::vector As(N); for (Measurement*& A : As) { @@ -106,7 +115,7 @@ int main(int argc, char* argv[]) { PT p(β₀, β₁, N, n, B, As); for (MCMC& sim : p.Ms) { - sim.M = walsh(k); + sim.M = hadamards[k - 1]; sim.E = sim.M.energy(); } @@ -141,7 +150,13 @@ int main(int argc, char* argv[]) { file << std::endl; for (unsigned i = 0; i < As.size(); i++) { - file << std::fixed << ((MeasureEnergy*)As[i])->totalE / ((MeasureEnergy*)As[i])->N << " "; + file << std::fixed << ((MeasureEnergy*)As[i])->energy() << " "; + } + + file << std::endl; + + for (unsigned i = 0; i < As.size(); i++) { + file << std::fixed << ((MeasureEnergy*)As[i])->specific_heat() << " "; } file << std::endl; diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp index eab9916..bca8ae2 100644 --- a/hadamard_mcmc.hpp +++ b/hadamard_mcmc.hpp @@ -18,6 +18,10 @@ public: } } + Orthogonal(const std::vector& x) : m(x) { + d = sqrt(x.size()); + } + unsigned size() const { return d; } double& operator()(unsigned i, unsigned j) { return m[d * i + j]; } @@ -37,28 +41,6 @@ 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; diff --git a/matrices.hpp b/matrices.hpp new file mode 100644 index 0000000..9dc24fd --- /dev/null +++ b/matrices.hpp @@ -0,0 +1,139 @@ +#pragma once + +#include "hadamard_mcmc.hpp" +#include + +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; + } +} + +Orthogonal m4 = walsh(2); +Orthogonal m8 = walsh(3); +Orthogonal m16 = walsh(4); +Orthogonal m32 = walsh(5); + +Orthogonal m12( + { + 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, + 1, 1,-1, 1,-1,-1,-1, 1, 1, 1,-1, 1, + 1, 1, 1,-1, 1,-1,-1,-1, 1, 1, 1,-1, + 1,-1, 1, 1,-1, 1,-1,-1,-1, 1, 1, 1, + 1, 1,-1, 1, 1,-1, 1,-1,-1,-1, 1, 1, + 1, 1, 1,-1, 1, 1,-1, 1,-1,-1,-1, 1, + 1, 1, 1, 1,-1, 1, 1,-1, 1,-1,-1,-1, + 1,-1, 1, 1, 1,-1, 1, 1,-1, 1,-1,-1, + 1,-1,-1, 1, 1, 1,-1, 1, 1,-1, 1,-1, + 1,-1,-1,-1, 1, 1, 1,-1, 1, 1,-1, 1, + 1, 1,-1,-1,-1, 1, 1, 1,-1, 1, 1,-1, + 1,-1, 1,-1,-1,-1, 1, 1, 1,-1, 1, 1 + } + ); + +Orthogonal m20( + { + 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, + 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, + 1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, + 1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1, + 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1, + 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, + 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, + 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, + 1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, + 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1, + 1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, + 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1, + 1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, + 1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1, + 1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1, + 1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1, + 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1, + 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, + 1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, + 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1 + } + ); + +Orthogonal m24( + { + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, + 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1, + 1, 1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1,-1,-1,-1, + 1, 1, 1, 1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, + 1, 1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1, + 1, 1,-1,-1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1, + 1, 1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1, + 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, + 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1, 1, 1, + 1, 1,-1,-1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1,-1,-1, 1, 1, + 1, 1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1, 1, 1, + 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, + 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1, + 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1, + 1,-1,-1, 1,-1, 1, 1,-1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1, 1,-1,-1, 1,-1, 1, + 1,-1, 1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, + 1,-1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1, + 1,-1,-1, 1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1,-1, 1, + 1,-1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1,-1, 1, 1,-1, + 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, + 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1, 1,-1, + 1,-1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1,-1, 1, 1,-1, + 1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1 + } + ); + +Orthogonal m28( + { + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1,-1,-1, 1, 1, 1,-1,-1,-1, 1, 1, 1,-1,-1,-1, 1, 1, 1,-1,-1,-1, 1, 1, 1, 1,-1,-1,-1, + 1,-1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1, 1,-1, 1, 1, 1, 1,-1,-1,-1, + 1,-1,-1,-1,-1,-1, 1, 1, 1,-1, 1, 1,-1,-1, 1,-1, 1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1, + 1,-1,-1, 1, 1,-1, 1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1,-1,-1,-1, 1, 1,-1, 1,-1, 1, 1,-1, + 1,-1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1, 1,-1,-1,-1, 1,-1, 1, 1,-1, 1,-1, 1, 1, 1,-1, + 1,-1,-1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1,-1,-1,-1, 1, 1, 1,-1,-1, 1, 1, 1, 1,-1, + -1, 1,-1, 1, 1,-1,-1,-1, 1,-1, 1,-1, 1,-1, 1,-1,-1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1, + -1, 1,-1,-1,-1, 1,-1, 1, 1,-1, 1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1, 1,-1, 1,-1, 1,-1,-1, + -1, 1,-1,-1,-1, 1, 1,-1, 1, 1,-1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1,-1, + -1, 1,-1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1, 1, 1,-1, 1, 1, 1,-1, 1,-1,-1, 1,-1, 1,-1, + -1, 1,-1,-1, 1,-1, 1, 1,-1, 1, 1,-1,-1, 1,-1,-1, 1, 1, 1,-1, 1,-1,-1, 1, 1,-1, 1,-1, + -1, 1,-1, 1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1, 1,-1,-1, 1, 1,-1, + -1,-1, 1,-1, 1, 1,-1, 1,-1,-1,-1, 1, 1, 1, 1,-1, 1,-1, 1,-1, 1, 1, 1,-1, 1,-1,-1,-1, + -1,-1, 1, 1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, + -1,-1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1, 1,-1, 1, 1,-1, 1,-1, 1,-1,-1,-1, 1, 1, 1,-1,-1, + -1,-1, 1, 1,-1,-1,-1, 1, 1, 1,-1, 1,-1, 1, 1, 1,-1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1, + -1,-1, 1, 1,-1,-1, 1,-1, 1, 1, 1, 1, 1,-1,-1,-1, 1,-1,-1, 1, 1, 1,-1,-1, 1,-1, 1,-1, + -1,-1, 1,-1, 1, 1,-1,-1, 1, 1, 1,-1, 1, 1,-1, 1,-1, 1, 1,-1,-1, 1,-1,-1,-1, 1, 1,-1, + 1, 1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1,-1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1,-1,-1, + 1, 1, 1, 1,-1,-1,-1, 1,-1, 1, 1,-1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1,-1, + 1, 1, 1, 1,-1,-1,-1,-1, 1,-1,-1, 1, 1, 1,-1,-1, 1, 1, 1,-1,-1,-1,-1, 1, 1, 1,-1,-1, + 1, 1, 1,-1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, + 1, 1, 1,-1,-1, 1, 1,-1,-1,-1, 1, 1, 1,-1,-1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1,-1, 1,-1, + 1, 1, 1,-1, 1,-1,-1,-1, 1,-1, 1,-1,-1, 1, 1, 1, 1,-1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1, + 1, 1,-1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, + 1,-1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1, + -1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1 + } + ); + +std::vector hadamards = {m4, m8, m12, m16, m20, m24, m28, m32}; + + -- cgit v1.2.3-70-g09d2