summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ciamarra.cpp1
-rw-r--r--distinguishable.cpp11
-rw-r--r--quantity.hpp119
3 files changed, 129 insertions, 2 deletions
diff --git a/ciamarra.cpp b/ciamarra.cpp
index d1087a4..d3acc38 100644
--- a/ciamarra.cpp
+++ b/ciamarra.cpp
@@ -118,7 +118,6 @@ class CiamarraSystem : public HardSystem<D, CiamarraState<D>> {
}
};
-
void print(const CiamarraSystem<2>& s) {
for (const Vertex<2, CiamarraState<2>>& v : s.vertices) {
if (v.state(0) == 1 && v.state(1) == 0) {
diff --git a/distinguishable.cpp b/distinguishable.cpp
index b739927..0aae0c1 100644
--- a/distinguishable.cpp
+++ b/distinguishable.cpp
@@ -1,6 +1,7 @@
#include <iostream>
#include "glass.hpp"
+#include "quantity.hpp"
class DistinguishableState {
public:
@@ -71,6 +72,7 @@ int main() {
}
}
+
unsigned n = 1;
for (double T = Tmax; T > Tmin; T -= δT) {
@@ -86,6 +88,7 @@ int main() {
std::cout << T << " ";
s.setInitialPosition();
auto start = std::chrono::high_resolution_clock::now();
+ Quantity<double> energy(1e3);
for (unsigned i = 0; i < 1e5; i++) {
for (unsigned j = 0; j < s.vertices.size(); j++) {
s.tryRandomSwap(T, r);
@@ -96,9 +99,10 @@ int main() {
s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
*/
if (i % 10 == 0) {
+ energy.add(s.energy());
auto stop = std::chrono::high_resolution_clock::now();
auto duration = duration_cast<std::chrono::microseconds>(stop - start);
- std::cout << duration.count() << " " << s.selfIntermediateScattering(ms) << " ";
+ std::cout /*<< duration.count() << " "*/ << s.selfIntermediateScattering(ms) << " ";
}
}
/*
@@ -110,6 +114,11 @@ int main() {
}
*/
std::cout << std::endl;
+ std::vector<double> rho = energy.ρ();
+ for (const double& x : rho) {
+ std::cout << x << " ";
+ }
+ std::cout << std::endl;
std::cerr << T << " " << s.energy() / N << std::endl;
// s.sweepLocal(r);
// s.sweepSwap(r);
diff --git a/quantity.hpp b/quantity.hpp
new file mode 100644
index 0000000..400a0dd
--- /dev/null
+++ b/quantity.hpp
@@ -0,0 +1,119 @@
+
+#pragma once
+
+#include <array>
+#include <cmath>
+#include <list>
+#include <vector>
+#include <fstream>
+#include <iomanip>
+
+template <class T>
+class Quantity {
+private:
+ uint64_t N;
+ uint64_t n;
+ unsigned skip;
+ std::list<T> hist;
+ double total;
+ double total2;
+ std::vector<double> C;
+
+public:
+ Quantity(unsigned lag, unsigned s = 1) : C(lag) {
+ skip = s;
+ N = 0;
+ n = 0;
+ total = 0;
+ total2 = 0;
+ }
+
+ void read(std::string filename) {
+ std::ifstream inFile(filename);
+ inFile >> n;
+ inFile >> total;
+ inFile >> total2;
+
+ for (double& Ci : C) {
+ inFile >> Ci;
+ }
+ }
+
+ void write(std::string filename) const {
+ std::ofstream outFile(filename);
+ outFile << std::setprecision(15) << n << " " << total << " " << total2 << std::endl;
+ for (double Ci : C) {
+ outFile << Ci << " ";
+ }
+ outFile << std::endl;
+ outFile.close();
+ }
+
+ void reset() {
+ total = 0;
+ total2 = 0;
+ std::fill(C.begin(), C.end(), 0);
+ n = 0;
+ hist = {};
+ }
+
+ 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++;
+ }
+ }
+ N++;
+ }
+
+ double avg() const { return total / n; }
+
+ double avg2() const { return total2 / n; }
+
+ std::vector<double> ρ() const {
+ double C0 = C.front() / n;
+ double avg2 = pow(total / n, 2);
+
+ std::vector<double> ρtmp;
+
+ for (double Ct : C) {
+ ρtmp.push_back((Ct / n - avg2) / (C0 - avg2));
+ }
+
+ return ρtmp;
+ }
+
+ std::array<double, 2> τ() const {
+ double τtmp = 0.5;
+ unsigned M = 1;
+ double c = 8.0;
+
+ std::vector<double> ρ_tmp = this->ρ();
+
+ while (c * τtmp > M && M < C.size()) {
+ τtmp += ρ_tmp[M];
+ M++;
+ }
+
+ return {skip * τtmp, skip * 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
+ }
+
+ double σ() const {
+ return 2.0 / n * this->τ()[0] * (C[0] / n - pow(this->avg(), 2));
+ }
+
+ double serr() const { return sqrt(this->σ()); }
+
+ unsigned num_added() const { return n; }
+};
+