summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--hadamard.cpp21
-rw-r--r--hadamard_mcmc.hpp26
-rw-r--r--matrices.hpp139
3 files changed, 161 insertions, 25 deletions
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 <fstream>
#include <iostream>
@@ -9,6 +10,7 @@ class MeasureEnergy : public Measurement {
public:
unsigned N;
double totalE;
+ double totalE2;
unsigned n;
std::vector<unsigned> ρ_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<Measurement*> 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<double>& 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 <array>
+
+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<Orthogonal> hadamards = {m4, m8, m12, m16, m20, m24, m28, m32};
+
+