summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--hadamard_mcmc.hpp11
-rw-r--r--hmm_correlation.cpp94
-rw-r--r--quantity.hpp39
3 files changed, 129 insertions, 15 deletions
diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index 107622a..634b822 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -68,6 +68,17 @@ public:
}
}
}
+
+ double operator*(const Orthogonal& o2) const {
+ const Orthogonal& o1 = *this;
+ double total = 0;
+ for (unsigned i = 0; i < this->size(); i++) {
+ for (unsigned j = 0; j < this->size(); j++) {
+ total += o1(i, j) * o2(i, j);
+ }
+ }
+ return total;
+ }
};
class Givens {
diff --git a/hmm_correlation.cpp b/hmm_correlation.cpp
new file mode 100644
index 0000000..80a6927
--- /dev/null
+++ b/hmm_correlation.cpp
@@ -0,0 +1,94 @@
+
+#include "hadamard_mcmc.hpp"
+#include "quantity.hpp"
+#include "matrices.hpp"
+
+#include <chrono>
+#include <fstream>
+#include <iostream>
+
+std::string getFilename(unsigned n, double β, double θ₀, unsigned lag, unsigned skip) {
+ return "correlation_" + std::to_string(n) + "_" + std::to_string(β) + "_" +
+ std::to_string(θ₀) + "_" + std::to_string(lag) + "_" + std::to_string(skip)
+ + ".dat";
+}
+
+class MeasureCorrelation : public Measurement {
+public:
+ Quantity<Orthogonal> Q;
+
+ MeasureCorrelation(unsigned lag, unsigned skip) : Q(lag, skip) {}
+
+ void after_sweep(double, double E, const Orthogonal& M) override {
+ Q.add(M);
+ }
+};
+
+int main(int argc, char* argv[]) {
+ // model parameters
+ unsigned n = 2; // matrix size over four
+
+ // simulation settings
+ double β = 1; // temperature
+ unsigned m = 1e4; // number of relaxation sweeps
+ unsigned N = 1e4; // number of measurement sweeps
+ double θ₀ = 0.05; // standard deviation of step size
+
+ // measurement settings
+ unsigned l = 1e3; // lag in autocorrelation function
+ unsigned s = 1; // skip in autocorrelation function
+
+ bool loadDataFromFile = false;
+
+ int opt;
+
+ while ((opt = getopt(argc, argv, "n:b:m:N:l:t:i:Ls:")) != -1) {
+ switch (opt) {
+ case 'n':
+ n = atoi(optarg);
+ break;
+ case 'b':
+ β = atof(optarg);
+ break;
+ case 'm':
+ m = (unsigned)atof(optarg);
+ break;
+ case 'N':
+ N = (unsigned)atof(optarg);
+ break;
+ case 'l':
+ l = (unsigned)atof(optarg);
+ break;
+ case 's':
+ s = (unsigned)atof(optarg);
+ break;
+ case 't':
+ θ₀ = atof(optarg);
+ break;
+ case 'L':
+ loadDataFromFile = true;
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ MeasureCorrelation Q(l, s);
+ MCMC simulation(n, β, Q, θ₀);
+
+ if (loadDataFromFile) {
+ Q.Q.read(getFilename(n, β, θ₀, l, s));
+ std::cout << "Imported data from file, number of steps " << Q.Q.num_added() << "." << std::endl;
+ }
+
+ simulation.run(m, true);
+ std::cout << "Finished initial relaxation at " << β << "." << std::endl;
+ simulation.run(N);
+
+ std::cout << "Finished simulation of " << β << ", writing output file." << std::endl;
+ std::string filename = getFilename(n, β, θ₀, l, s);
+ Q.Q.write(filename);
+
+ return 0;
+}
+
diff --git a/quantity.hpp b/quantity.hpp
index 296a6e8..400a0dd 100644
--- a/quantity.hpp
+++ b/quantity.hpp
@@ -8,16 +8,21 @@
#include <fstream>
#include <iomanip>
+template <class T>
class Quantity {
private:
- unsigned n;
- std::list<double> hist;
+ uint64_t N;
+ uint64_t n;
+ unsigned skip;
+ std::list<T> hist;
double total;
double total2;
std::vector<double> C;
public:
- Quantity(unsigned lag) : C(lag) {
+ Quantity(unsigned lag, unsigned s = 1) : C(lag) {
+ skip = s;
+ N = 0;
n = 0;
total = 0;
total2 = 0;
@@ -52,19 +57,23 @@ public:
hist = {};
}
- void add(double x) {
- hist.push_front(x);
- if (hist.size() > C.size()) {
- hist.pop_back();
- unsigned t = 0;
- for (double a : hist) {
- C[t] += a * x;
- t++;
+ void add(const T& x) {
+ if (N % skip == 0) {
+ hist.push_front(x);
+ if (hist.size() > C.size()) {
+ hist.pop_back();
+ unsigned t = 0;
+ for (T a : hist) {
+ C[t] += a * x;
+ t++;
+ }
+ double norm = x * x;
+ total += sqrt(norm);
+ total2 += norm;
+ n++;
}
- total += x;
- total2 += pow(x, 2);
- n++;
}
+ N++;
}
double avg() const { return total / n; }
@@ -96,7 +105,7 @@ public:
M++;
}
- return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
+ return {skip * τtmp, skip * 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
}
double σ() const {