summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--examples/CMakeLists.txt3
-rw-r--r--examples/include/correlation.hpp10
-rw-r--r--examples/include/measure.hpp164
-rw-r--r--examples/simple_ising.cpp192
-rw-r--r--examples/src/models/On/wolff_On.cpp25
-rw-r--r--examples/src/models/ising/ising.hpp33
-rw-r--r--examples/src/models/ising/wolff_ising.cpp25
-rw-r--r--examples/src/models/potts/symmetric.hpp1
-rw-r--r--examples/src/models/potts/wolff_clock.cpp22
-rw-r--r--examples/src/models/potts/wolff_potts.cpp25
-rw-r--r--lib/include/wolff.hpp15
-rw-r--r--lib/include/wolff/cluster.hpp20
-rw-r--r--lib/include/wolff/graph.hpp9
-rw-r--r--lib/include/wolff/meas.h19
-rw-r--r--lib/include/wolff/state.hpp81
-rw-r--r--lib/src/graph.cpp60
16 files changed, 491 insertions, 213 deletions
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index baf0a10..a435b29 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -3,5 +3,8 @@ add_library(wolff_examples INTERFACE)
target_include_directories(wolff_examples INTERFACE include)
+add_executable(simple_ising simple_ising.cpp)
+target_link_libraries(simple_ising wolff)
+
add_subdirectory(src)
diff --git a/examples/include/correlation.hpp b/examples/include/correlation.hpp
index 042cff3..da8ab7f 100644
--- a/examples/include/correlation.hpp
+++ b/examples/include/correlation.hpp
@@ -6,18 +6,18 @@
#include <fftw3.h>
-template <class R_t, class X_t>
-double correlation_length(const state_t <R_t, X_t>& s) {
+template <class X_t>
+double correlation_length(const std::vector<typename X_t::F_t>& ReF, const std::vector<typename X_t::F_t>& ImF, D_t D) {
double total = 0;
#ifdef DIMENSION
for (D_t j = 0; j < DIMENSION; j++) {
#else
- for (D_t j = 0; j < s.D; j++) {
+ for (D_t j = 0; j < D; j++) {
#endif
- total += norm_squared(s.ReF[j]) + norm_squared(s.ImF[j]);
+ total += norm_squared(ReF[j]) + norm_squared(ImF[j]);
}
- return total / s.D;
+ return total / D;
}
diff --git a/examples/include/measure.hpp b/examples/include/measure.hpp
index e20353c..bf8baab 100644
--- a/examples/include/measure.hpp
+++ b/examples/include/measure.hpp
@@ -2,6 +2,7 @@
#pragma once
#include <wolff/state.hpp>
+#include <wolff/meas.h>
#include "correlation.hpp"
#include <functional>
@@ -13,51 +14,146 @@ const unsigned char measurement_fourierZero = 1 << 3;
char const *measurement_labels[] = {"E", "S", "M", "F"};
-FILE **measure_setup_files(unsigned char flags, unsigned long timestamp) {
- FILE **files = (FILE **)calloc(POSSIBLE_MEASUREMENTS, sizeof(FILE *));
+template <class R_t, class X_t>
+class wolff_research_measurements : public wolff_measurement<R_t, X_t> {
+ private:
+ std::vector<std::vector<double>> precomputed_sin;
+ std::vector<std::vector<double>> precomputed_cos;
+ FILE **files;
+ D_t D;
+ unsigned char flags;
+ std::function <void(const state_t <R_t, X_t>&, const wolff_research_measurements<R_t, X_t>&)> other_f;
+ bool silent;
+ public:
+ v_t last_cluster_size;
+ double E;
+ typename X_t::M_t M;
+ std::vector<typename X_t::F_t> ReF;
+ std::vector<typename X_t::F_t> ImF;
+
+ wolff_research_measurements(unsigned char flags, unsigned long timestamp, std::function <void(const state_t <R_t, X_t>&, const wolff_research_measurements<R_t, X_t>&)> other_f, state_t<R_t, X_t>& s, bool silent) : flags(flags), D(s.D), other_f(other_f), silent(silent) {
+ FILE **files = (FILE **)calloc(POSSIBLE_MEASUREMENTS, sizeof(FILE *));
+
+ for (uint8_t i = 0; i < POSSIBLE_MEASUREMENTS; i++) {
+ if (flags & (1 << i)) {
+ char *filename = (char *)malloc(255 * sizeof(char));
+ sprintf(filename, "wolff_%lu_%s.dat", timestamp, measurement_labels[i]);
+ files[i] = fopen(filename, "wb");
+ free(filename);
+ }
+ }
+
+#ifdef BOND_DEPENDENCE
+ E = 0;
+ for (v_t v = 0; v < s.nv; v++) {
+ for (const v_t &vn : s.g.v_adj[v]) {
+ if (v < vn) {
+ E -= s.J(v, s.spins[v], vn, s.spins[vn]);
+ }
+ }
+ }
+#else
+ E = - (double)s.ne * s.J(s.spins[0], s.spins[0]);
+#endif
+
+#ifndef NOFIELD
+#ifdef SITE_DEPENDENCE
+ for (v_t i = 0; i < s.nv; i++) {
+ E -= s.H(i, s.spins[i]);
+ }
+#else
+ E -= (double)s.nv * s.H(s.spins[0]);
+#endif
+#endif
+
+ M = s.spins[0] * s.nv;
+
+ ReF.resize(s.D);
+ ImF.resize(s.D);
+ for (D_t i = 0; i < s.D; i++) {
+ ReF[i] = s.spins[0] * 0.0;
+ ImF[i] = s.spins[0] * 0.0;
+ }
+ precomputed_cos.resize(s.nv);
+ precomputed_sin.resize(s.nv);
+ for (v_t i = 0; i < s.nv; i++) {
+ precomputed_cos[i].resize(s.D);
+ precomputed_sin[i].resize(s.D);
+ for (D_t j = 0; j < s.D; j++) {
+ precomputed_cos[i][j] = cos(2 * M_PI * s.g.coordinate[i][j] / (double)s.L);
+ precomputed_sin[i][j] = sin(2 * M_PI * s.g.coordinate[i][j] / (double)s.L);
+ }
+ }
+
+ if (!silent) printf("\n");
- for (uint8_t i = 0; i < POSSIBLE_MEASUREMENTS; i++) {
- if (flags & (1 << i)) {
- char *filename = (char *)malloc(255 * sizeof(char));
- sprintf(filename, "wolff_%lu_%s.dat", timestamp, measurement_labels[i]);
- files[i] = fopen(filename, "wb");
- free(filename);
}
- }
- return files;
-}
+ void pre_cluster(const state_t<R_t, X_t>& s, count_t step, count_t N, v_t v0, const R_t& R) {
+ last_cluster_size = 0;
+ if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", step, N, E, last_cluster_size);
-template <class R_t, class X_t>
-std::function <void(const state_t <R_t, X_t>&)> measure_function_write_files(unsigned char flags, FILE **files, std::function <void(const state_t <R_t, X_t>&)> other_f) {
- return [=] (const state_t <R_t, X_t>& s) {
- if (flags & measurement_energy) {
- float smaller_E = (float)s.E;
- fwrite(&smaller_E, sizeof(float), 1, files[0]);
}
- if (flags & measurement_clusterSize) {
- fwrite(&(s.last_cluster_size), sizeof(uint32_t), 1, files[1]);
+
+ void plain_bond_added(v_t vi, const X_t& si_old, const X_t& si_new, v_t vn, const X_t& sn, double dE) {
+ E += dE;
}
- if (flags & measurement_magnetization) {
- write_magnetization(s.M, files[2]);
+
+ void ghost_bond_added(v_t v, const X_t& rs_old, const X_t& rs_new, double dE) {
+ E += dE;
+ M += rs_new - rs_old;
+
+#ifdef DIMENSION
+ for (D_t i = 0; i < DIMENSION; i++)
+#else
+ for (D_t i = 0; i < D; i++)
+#endif
+ {
+ ReF[i] += (rs_new - rs_old) * precomputed_cos[v][i];
+ ImF[i] += (rs_new - rs_old) * precomputed_sin[v][i];
+ }
}
- if (flags & measurement_fourierZero) {
- float smaller_X = (float)correlation_length(s);
- fwrite(&smaller_X, sizeof(float), 1, files[3]);
+
+ void plain_site_transformed(v_t v, const X_t& s_old, const X_t& s_new) {
+ last_cluster_size++;
}
- other_f(s);
- };
-}
+ void ghost_site_transformed(const R_t& r_old, const R_t& r_new) {
+ }
+
+ void post_cluster(const state_t<R_t, X_t>& s, count_t step, count_t N) {
+ if (flags & measurement_energy) {
+ float smaller_E = (float)E;
+ fwrite(&smaller_E, sizeof(float), 1, files[0]);
+ }
+ if (flags & measurement_clusterSize) {
+ fwrite(&(last_cluster_size), sizeof(uint32_t), 1, files[1]);
+ }
+ if (flags & measurement_magnetization) {
+ write_magnetization(M, files[2]);
+ }
+ if (flags & measurement_fourierZero) {
+ float smaller_X = (float)correlation_length<X_t>(ReF, ImF, D);
+ fwrite(&smaller_X, sizeof(float), 1, files[3]);
+ }
+
+ other_f(s, *this);
-void measure_free_files(unsigned char flags, FILE **files) {
- for (uint8_t i = 0; i < POSSIBLE_MEASUREMENTS; i++) {
- if (flags & (1 << i)) {
- fclose(files[i]);
}
- }
- free(files);
-}
+ ~wolff_research_measurements() {
+ for (uint8_t i = 0; i < POSSIBLE_MEASUREMENTS; i++) {
+ if (flags & (1 << i)) {
+ fclose(files[i]);
+ }
+ }
+ if (!silent) {
+ printf("\033[F\033[J");
+ }
+ printf("WOLFF COMPLETE: E = %.2f, S = %" PRIv "\n", E, last_cluster_size);
+
+ free(files);
+ }
+};
diff --git a/examples/simple_ising.cpp b/examples/simple_ising.cpp
new file mode 100644
index 0000000..24e4ae5
--- /dev/null
+++ b/examples/simple_ising.cpp
@@ -0,0 +1,192 @@
+
+#include <getopt.h>
+#include <iostream>
+
+#include "include/randutils/randutils.hpp"
+
+#include <wolff.hpp>
+
+// define your R_t and X_t. here both are the same, ising_t
+class ising_t {
+ public:
+ bool x;
+
+ // both X_t and R_t need a default constructor (and destructor, if relevant)
+ ising_t() : x(false) {}
+ ising_t(bool x) : x(x) {}
+
+ // R_t needs the member functions
+ // X_t act(const X_t& s) const {}
+ // R_t act(const R_t& s) const {}
+ // to define action on both spins and other transformations
+ ising_t act(const ising_t& s) const {
+ if (x) {
+ return ising_t(!(s.x));
+ } else {
+ return ising_t(s.x);
+ }
+ }
+
+ // R_t needs the member functions
+ // X_t act_inverse(const X_t& s) const {}
+ // R_t act_inverse(const R_t& s) const {}
+ // to define action of its inverse on both spins and other transformations
+ ising_t act_inverse(const ising_t& s) const {
+ return this->act(s);
+ }
+};
+
+// define how measurements should be taken by importing the interface wolff_measurement<R_t, X_t>
+class ising_measurements : public wolff_measurement<ising_t, ising_t> {
+ private:
+ count_t n;
+
+ double E;
+ int M;
+ v_t S;
+
+ double totalE;
+ double totalM;
+ double totalS;
+
+ public:
+ ising_measurements(D_t D, L_t L, double H) {
+ n = 0;
+ M = (int)pow(L, D);
+ E = -D * pow(L, D) - H * pow(L, D);
+
+ totalE = 0;
+ totalM = 0;
+ totalS = 0;
+ }
+
+ void pre_cluster(const state_t<ising_t, ising_t>& s, count_t step, count_t N, v_t v0, const ising_t& R) {
+ S = 0;
+ }
+
+ void plain_bond_added(v_t v, const ising_t& s_old, const ising_t& s_new, v_t vn, const ising_t& sn, double dE) {
+ E += dE;
+ }
+
+ void ghost_bond_added(v_t v, const ising_t& s_old, const ising_t& s_new, double dE) {
+ E += dE;
+
+ if (s_old.x) {
+ M++;
+ } else {
+ M--;
+ }
+
+ if (s_new.x) {
+ M--;
+ } else {
+ M++;
+ }
+ }
+
+ void plain_site_transformed(v_t v, const ising_t& s_old, const ising_t& s_new) {
+ S++;
+ }
+
+ void ghost_site_transformed(const ising_t& R_old, const ising_t& R_new) {
+ }
+
+ void post_cluster(const state_t<ising_t, ising_t>& s, count_t step, count_t N) {
+ totalE += E;
+ totalM += M;
+ totalS += S;
+ n++;
+ }
+
+ double avgE() {
+ return totalE / n;
+ }
+
+ double avgM() {
+ return totalM / n;
+ }
+
+ double avgS() {
+ return totalS / n;
+ }
+};
+
+int main(int argc, char *argv[]) {
+
+ // set defaults
+ count_t N = (count_t)1e4;
+ D_t D = 2;
+ L_t L = 128;
+ double T = 2.26918531421;
+ double H = 0.0;
+
+ int opt;
+
+ // take command line arguments
+ while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
+ switch (opt) {
+ case 'N': // number of steps
+ N = (count_t)atof(optarg);
+ break;
+ case 'D': // dimension
+ D = atoi(optarg);
+ break;
+ case 'L': // linear size
+ L = atoi(optarg);
+ break;
+ case 'T': // temperature
+ T = atof(optarg);
+ break;
+ case 'H': // external field
+ H = atof(optarg);
+ break;
+ default:
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ // define the spin-spin coupling
+ std::function <double(const ising_t&, const ising_t&)> Z = [] (const ising_t& s1, const ising_t& s2) -> double {
+ if (s1.x == s2.x) {
+ return 1.0;
+ } else {
+ return -1.0;
+ }
+ };
+
+ // define the spin-field coupling
+ std::function <double(const ising_t&)> B = [=] (const ising_t& s) -> double {
+ if (s.x) {
+ return -H;
+ } else {
+ return H;
+ }
+ };
+
+ // initialize the system
+ state_t<ising_t, ising_t> s(D, L, T, Z, B);
+
+ // initialize the random number generator
+ randutils::auto_seed_128 seeds;
+ std::mt19937 rng{seeds};
+
+ // define function that generates self-inverse rotations
+ std::function <ising_t(std::mt19937&, const ising_t&)> gen_R = [] (std::mt19937&, const ising_t& s) -> ising_t {
+ return ising_t(true);
+ };
+
+ // initailze the measurement object
+ ising_measurements m(D, L, H);
+
+ // run wolff N times
+ wolff<ising_t, ising_t>(N, s, gen_R, m, rng);
+
+ // print the result of our measurements
+ std::cout << "Wolff complete!\nThe average energy per site was " << m.avgE() / s.nv
+ << ".\nThe average magnetization per site was " << m.avgM() / s.nv
+ << ".\nThe average cluster size per site was " << m.avgS() / s.nv << ".\n";
+
+ // exit
+ return 0;
+}
+
diff --git a/examples/src/models/On/wolff_On.cpp b/examples/src/models/On/wolff_On.cpp
index ad2ac77..67f28a5 100644
--- a/examples/src/models/On/wolff_On.cpp
+++ b/examples/src/models/On/wolff_On.cpp
@@ -188,14 +188,12 @@ int main(int argc, char *argv[]) {
fclose(outfile_info);
- FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
-
- std::function <void(const On_t&)> other_f;
+ std::function <void(const On_t&, const wolff_research_measurements<orthogonal_R_t, vector_R_t>&)> other_f;
uint64_t sum_of_clusterSize = 0;
if (N_is_sweeps) {
- other_f = [&] (const On_t& s) {
- sum_of_clusterSize += s.last_cluster_size;
+ other_f = [&] (const On_t& s, const wolff_research_measurements<orthogonal_R_t, vector_R_t>& m) {
+ sum_of_clusterSize += m.last_cluster_size;
};
} else if (draw) {
#ifdef HAVE_GLUT
@@ -209,7 +207,7 @@ int main(int argc, char *argv[]) {
glLoadIdentity();
gluOrtho2D(0.0, L, 0.0, L);
- other_f = [&] (const On_t& s) {
+ other_f = [&] (const On_t& s, const wolff_research_measurements<orthogonal_R_t, vector_R_t>& m) {
glClear(GL_COLOR_BUFFER_BIT);
for (v_t i = 0; i < pow(L, 2); i++) {
#ifdef NOFIELD
@@ -228,11 +226,9 @@ int main(int argc, char *argv[]) {
};
#endif
} else {
- other_f = [] (const On_t& s) {};
+ other_f = [] (const On_t& s, const wolff_research_measurements<orthogonal_R_t, vector_R_t>& m) {};
}
- std::function <void(const On_t&)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f);
-
std::function <double(const vector_R_t&)> H;
if (modulated_field) {
@@ -251,20 +247,21 @@ int main(int argc, char *argv[]) {
state_t <orthogonal_R_t, vector_R_t> s(D, L, T, dot <N_COMP, double>);
#endif
+ wolff_research_measurements<orthogonal_R_t, vector_R_t> m(measurement_flags, timestamp, other_f, s, silent);
+
if (N_is_sweeps) {
count_t N_rounds = 0;
printf("\n");
while (sum_of_clusterSize < N * s.nv) {
- printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
- wolff <orthogonal_R_t, vector_R_t> (N, s, gen_R, measurements, rng, silent);
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, m.E, m.last_cluster_size);
+ wolff <orthogonal_R_t, vector_R_t> (N, s, gen_R, m, rng);
N_rounds++;
}
- printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, m.E, m.last_cluster_size);
} else {
- wolff <orthogonal_R_t, vector_R_t> (N, s, gen_R, measurements, rng, silent);
+ wolff <orthogonal_R_t, vector_R_t> (N, s, gen_R, m, rng);
}
- measure_free_files(measurement_flags, outfiles);
free(H_vec);
return 0;
diff --git a/examples/src/models/ising/ising.hpp b/examples/src/models/ising/ising.hpp
index ae20840..73b06ed 100644
--- a/examples/src/models/ising/ising.hpp
+++ b/examples/src/models/ising/ising.hpp
@@ -5,17 +5,31 @@
#include <wolff/types.h>
+// all that is required to use wolff.hpp is a default constructor
class ising_t {
public:
bool x;
- typedef int M_t;
- typedef double F_t;
-
ising_t() : x(false) {}
- ising_t(bool x) : x(x) {}
+
+ // optional constructors for syntactic sugar
+ ising_t(bool x) : x(x) {}
ising_t(int x) : x((bool)x) {}
+ /* below this comment is code required only for using measure.hpp in the
+ * examples folder, which provides an interface for measuring several
+ * generic features of models. these require
+ *
+ * - an M_t, representing the magnetization or sum of all spins
+ * - an F_t, representing a double-weighted version of the magnetization
+ * - the overloaded operator *, which takes a v_t (unsigned int) and returns an M_t
+ * - the overloaded operator *, which takes a double and returns an F_t
+ * - the overloaded operator -, which takes another X_t and returns an M_t
+ */
+
+ typedef int M_t;
+ typedef double F_t;
+
inline int operator*(v_t a) const {
if (x) {
return -(int)a;
@@ -45,6 +59,12 @@ class ising_t {
}
};
+/* using measure.hpp additionally requires a norm_squared function which takes
+ * an F_t to a double, and a write_magnetization function, which takes an M_t
+ * and a FILE pointer and appropriately records the contents of the former to
+ * the latter.
+ */
+
double norm_squared(double s) {
return pow(s, 2);
}
@@ -53,6 +73,11 @@ void write_magnetization(int M, FILE *outfile) {
fwrite(&M, sizeof(int), 1, outfile);
}
+/* these definitions allow wolff/finite_states.hpp to be invoked and provide
+ * much faster performance for models whose number of possible spin
+ * configurations is finite.
+ */
+
#define N_STATES 2
const ising_t states[2] = {ising_t(0), ising_t(1)};
q_t state_to_ind(ising_t state) { return (q_t)state.x; }
diff --git a/examples/src/models/ising/wolff_ising.cpp b/examples/src/models/ising/wolff_ising.cpp
index f7a3ada..de04f32 100644
--- a/examples/src/models/ising/wolff_ising.cpp
+++ b/examples/src/models/ising/wolff_ising.cpp
@@ -127,14 +127,12 @@ int main(int argc, char *argv[]) {
return z2_t(true);
};
- FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
-
- std::function <void(const state_t<z2_t, ising_t>&)> other_f;
+ std::function <void(const state_t<z2_t, ising_t>&, const wolff_research_measurements<z2_t, ising_t>&)> other_f;
uint64_t sum_of_clusterSize = 0;
if (N_is_sweeps) {
- other_f = [&] (const state_t<z2_t, ising_t>& s) {
- sum_of_clusterSize += s.last_cluster_size;
+ other_f = [&] (const state_t<z2_t, ising_t>& s, const wolff_research_measurements<z2_t, ising_t>& meas) {
+ sum_of_clusterSize += meas.last_cluster_size;
};
} else if (draw) {
#ifdef HAVE_GLUT
@@ -148,7 +146,7 @@ int main(int argc, char *argv[]) {
glLoadIdentity();
gluOrtho2D(0.0, L, 0.0, L);
- other_f = [] (const state_t <z2_t, ising_t>& s) {
+ other_f = [] (const state_t <z2_t, ising_t>& s, const wolff_research_measurements<z2_t, ising_t>& meas) {
glClear(GL_COLOR_BUFFER_BIT);
for (v_t i = 0; i < pow(s.L, 2); i++) {
#ifdef NOFIELD
@@ -166,10 +164,10 @@ int main(int argc, char *argv[]) {
};
#endif
} else {
- other_f = [] (const state_t<z2_t, ising_t>& s) {};
+ other_f = [] (const state_t<z2_t, ising_t>& s, const wolff_research_measurements<z2_t, ising_t>& meas) {};
}
- std::function <void(const state_t<z2_t, ising_t>&)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f);
+ wolff_research_measurements<z2_t, ising_t> m(measurement_flags, timestamp, other_f, s, silent);
// add line to metadata file with run info
{
@@ -185,18 +183,15 @@ int main(int argc, char *argv[]) {
count_t N_rounds = 0;
printf("\n");
while (sum_of_clusterSize < N * s.nv) {
- printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
- wolff(N, s, gen_R, measurements, rng, silent);
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, m.E, m.last_cluster_size);
+ wolff(N, s, gen_R, m, rng);
N_rounds++;
}
- printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, m.E, m.last_cluster_size);
} else {
- wolff(N, s, gen_R, measurements, rng, silent);
+ wolff(N, s, gen_R, m, rng);
}
- measure_free_files(measurement_flags, outfiles);
-
return 0;
-
}
diff --git a/examples/src/models/potts/symmetric.hpp b/examples/src/models/potts/symmetric.hpp
index 8636f15..bc8673f 100644
--- a/examples/src/models/potts/symmetric.hpp
+++ b/examples/src/models/potts/symmetric.hpp
@@ -36,6 +36,7 @@ class symmetric_t : public std::array<q_t, q> {
}
}
+ printf("Your spin wasn't a valid state!", s.x);
exit(EXIT_FAILURE);
}
diff --git a/examples/src/models/potts/wolff_clock.cpp b/examples/src/models/potts/wolff_clock.cpp
index 020415d..0706cc5 100644
--- a/examples/src/models/potts/wolff_clock.cpp
+++ b/examples/src/models/potts/wolff_clock.cpp
@@ -9,6 +9,7 @@
#include "dihedral.hpp"
#include "potts.hpp"
#include <colors.h>
+#include <measure.hpp>
// hack to speed things up considerably
#define N_STATES POTTSQ
@@ -95,7 +96,7 @@ int main(int argc, char *argv[]) {
std::function <dihedral_t<POTTSQ>(std::mt19937&, potts_t<POTTSQ>)> gen_R = [] (std::mt19937& r, potts_t<POTTSQ> v) -> dihedral_t<POTTSQ> {
dihedral_t<POTTSQ> rot;
rot.is_reflection = true;
- std::uniform_int_distribution<q_t> dist(0, POTTSQ - 1);
+ std::uniform_int_distribution<q_t> dist(0, POTTSQ - 2);
q_t x = dist(r);
rot.x = (2 * v.x + x + 1) % POTTSQ;
@@ -103,13 +104,11 @@ int main(int argc, char *argv[]) {
};
// define function that updates any number of measurements
- std::function <void(const sim_t&)> measurement;
+ std::function <void(const sim_t&, const wolff_research_measurements<dihedral_t<POTTSQ>, potts_t<POTTSQ>>&)> measurement;
- double average_M = 0;
if (!draw) {
// a very simple example: measure the average magnetization
- measurement = [&] (const sim_t& s) {
- average_M += (double)s.M[0] / (double)N / (double)s.nv;
+ measurement = [&] (const sim_t& s, const wolff_research_measurements<dihedral_t<POTTSQ>, potts_t<POTTSQ>>&) {
};
} else {
// a more complex example: measure the average magnetization, and draw the spin configuration to the screen
@@ -125,8 +124,7 @@ int main(int argc, char *argv[]) {
glLoadIdentity();
gluOrtho2D(0.0, L, 0.0, L);
- measurement = [&] (const sim_t& s) {
- average_M += (double)s.M[0] / (double)N / (double)s.nv;
+ measurement = [&] (const sim_t& s, const wolff_research_measurements<dihedral_t<POTTSQ>, potts_t<POTTSQ>>&) {
glClear(GL_COLOR_BUFFER_BIT);
for (v_t i = 0; i < pow(L, 2); i++) {
potts_t<POTTSQ> tmp_s = s.R.act_inverse(s.spins[i]);
@@ -138,17 +136,13 @@ int main(int argc, char *argv[]) {
#endif
}
- // run wolff for N cluster flips
- wolff(N, s, gen_R, measurement, rng, silent);
+ wolff_research_measurements<dihedral_t<POTTSQ>, potts_t<POTTSQ>> m(0, 0, measurement, s, silent);
- // tell us what we found!
- printf("%" PRIcount " %d-Potts runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, <M> = %g\n", N, POTTSQ, D, L, T, H_vec[0], average_M);
+ // run wolff for N cluster flips
+ wolff(N, s, gen_R, m, rng);
// free the random number generator
- if (draw) {
- }
-
return 0;
}
diff --git a/examples/src/models/potts/wolff_potts.cpp b/examples/src/models/potts/wolff_potts.cpp
index a1e9284..7b92ac1 100644
--- a/examples/src/models/potts/wolff_potts.cpp
+++ b/examples/src/models/potts/wolff_potts.cpp
@@ -119,7 +119,7 @@ int main(int argc, char *argv[]) {
std::function <symmetric_t<POTTSQ>(std::mt19937&, potts_t<POTTSQ>)> gen_R = [] (std::mt19937& r, potts_t<POTTSQ> v) -> symmetric_t<POTTSQ> {
symmetric_t<POTTSQ> rot;
- std::uniform_int_distribution<q_t> dist(0, POTTSQ - 1);
+ std::uniform_int_distribution<q_t> dist(0, POTTSQ - 2);
q_t j = dist(r);
q_t swap_v;
if (j < v.x) {
@@ -134,14 +134,12 @@ int main(int argc, char *argv[]) {
return rot;
};
- FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
-
- std::function <void(const sim_t&)> other_f;
+ std::function <void(const sim_t&, const wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>>&)> other_f;
uint64_t sum_of_clusterSize = 0;
if (N_is_sweeps) {
- other_f = [&] (const sim_t& s) {
- sum_of_clusterSize += s.last_cluster_size;
+ other_f = [&] (const sim_t& s, const wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>>& m) {
+ sum_of_clusterSize += m.last_cluster_size;
};
} else if (draw) {
#ifdef HAVE_GLUT
@@ -155,7 +153,7 @@ int main(int argc, char *argv[]) {
glLoadIdentity();
gluOrtho2D(0.0, L, 0.0, L);
- other_f = [] (const sim_t& s) {
+ other_f = [] (const sim_t& s, const wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>>& m) {
glClear(GL_COLOR_BUFFER_BIT);
for (v_t i = 0; i < pow(s.L, 2); i++) {
potts_t<POTTSQ> tmp_s = s.R.act_inverse(s.spins[i]);
@@ -166,10 +164,10 @@ int main(int argc, char *argv[]) {
};
#endif
} else {
- other_f = [] (const sim_t& s) {};
+ other_f = [] (const sim_t& s, const wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>>& m) {};
}
- std::function <void(const sim_t&)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f);
+ wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>> m(measurement_flags, timestamp, other_f, s, silent);
// add line to metadata file with run info
{
@@ -194,18 +192,17 @@ int main(int argc, char *argv[]) {
count_t N_rounds = 0;
printf("\n");
while (sum_of_clusterSize < N * s.nv) {
- printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
- wolff(N, s, gen_R, measurements, rng, silent);
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, m.E, m.last_cluster_size);
+ wolff(N, s, gen_R, m, rng);
N_rounds++;
}
- printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, m.E, m.last_cluster_size);
} else {
- wolff(N, s, gen_R, measurements, rng, silent);
+ wolff(N, s, gen_R, m, rng);
}
// free the random number generator
free(H_vec);
- measure_free_files(measurement_flags, outfiles);
return 0;
diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp
index c10a211..b730c8d 100644
--- a/lib/include/wolff.hpp
+++ b/lib/include/wolff.hpp
@@ -3,7 +3,7 @@
#include "wolff/state.hpp"
template <class R_t, class X_t>
-void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(std::mt19937&, X_t)> gen_R, std::function <void(const state_t <R_t, X_t>&)> measurements, std::mt19937& r, bool silent) {
+void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(std::mt19937&, X_t)> gen_R, wolff_measurement<R_t, X_t>& m, std::mt19937& r) {
#ifdef FINITE_STATES
#ifdef NOFIELD
@@ -15,21 +15,16 @@ void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(std::mt19937&, X
std::uniform_int_distribution<v_t> dist(0, s.nv);
- if (!silent) printf("\n");
for (count_t steps = 0; steps < N; steps++) {
- if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, s.last_cluster_size);
-
v_t v0 = dist(r);
R_t step = gen_R(r, s.spins[v0]);
- flip_cluster <R_t, X_t> (s, v0, step, r);
- measurements(s);
- }
+ m.pre_cluster(s, steps, N, v0, step);
+
+ flip_cluster<R_t, X_t>(s, v0, step, r, m);
- if (!silent) {
- printf("\033[F\033[J");
+ m.post_cluster(s, steps, N);
}
- printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size);
}
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
index e9dff7b..805e2c3 100644
--- a/lib/include/wolff/cluster.hpp
+++ b/lib/include/wolff/cluster.hpp
@@ -9,11 +9,11 @@
#include "types.h"
#include "state.hpp"
#include "graph.hpp"
+#include "meas.h"
template <class R_t, class X_t>
-void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand) {
- std::uniform_real_distribution<double> dist(0.0,1.0);
- v_t nv = 0;
+void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand, wolff_measurement<R_t, X_t>& m) {
+ std::uniform_real_distribution<double> dist(0.0, 1.0);
std::stack<v_t> stack;
stack.push(v0);
@@ -75,8 +75,8 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand
prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
#endif
- s.update_magnetization(rs_old, rs_new);
- s.update_fourierZero(non_ghost, rs_old, rs_new);
+ // run measurement hooks for encountering a ghost bond
+ m.ghost_bond_added(non_ghost, rs_old, rs_new, dE);
} else // this is a perfectly normal bond!
#endif
{
@@ -90,9 +90,10 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand
#ifdef FINITE_STATES
prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])];
#endif
- }
- s.update_energy(dE);
+ // run measurement hooks for encountering a plain bond
+ m.plain_bond_added(v, s.spins[v], si_new, vn, s.spins[vn], dE);
+ }
#ifndef FINITE_STATES
prob = 1.0 - exp(-dE / s.T);
@@ -105,16 +106,15 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand
#ifndef NOFIELD
if (v_is_ghost) {
+ m.ghost_site_transformed(s.R, R_new);
s.R = R_new;
} else
#endif
{
+ m.plain_site_transformed(v, s.spins[v], si_new);
s.spins[v] = si_new;
- nv++;
}
}
}
-
- s.last_cluster_size = nv;
}
diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp
index 165544a..0aeb6af 100644
--- a/lib/include/wolff/graph.hpp
+++ b/lib/include/wolff/graph.hpp
@@ -1,13 +1,16 @@
#pragma once
-#include <inttypes.h>
#include <cmath>
-#include <stdlib.h>
#include <vector>
#include "types.h"
+typedef enum lattice_t {
+ SQUARE_LATTICE,
+ DIAGONAL_LATTICE
+} lattice_t;
+
class graph_t {
public:
v_t ne;
@@ -15,7 +18,7 @@ class graph_t {
std::vector<std::vector<v_t>> v_adj;
std::vector<std::vector<double>> coordinate;
- graph_t(D_t D, L_t L);
+ graph_t(D_t D, L_t L, lattice_t lat = SQUARE_LATTICE);
void add_ext();
};
diff --git a/lib/include/wolff/meas.h b/lib/include/wolff/meas.h
new file mode 100644
index 0000000..4895eac
--- /dev/null
+++ b/lib/include/wolff/meas.h
@@ -0,0 +1,19 @@
+
+#pragma once
+
+#include "types.h"
+
+template <class R_t, class X_t>
+class wolff_measurement {
+ public:
+ virtual void pre_cluster(const state_t<R_t, X_t>&, count_t, count_t, v_t, const R_t&) = 0;
+
+ virtual void plain_bond_added(v_t, const X_t&, const X_t&, v_t, const X_t&, double) = 0;
+ virtual void ghost_bond_added(v_t, const X_t&, const X_t&, double) = 0;
+
+ virtual void plain_site_transformed(v_t, const X_t&, const X_t&) = 0;
+ virtual void ghost_site_transformed(const R_t&, const R_t&) = 0;
+
+ virtual void post_cluster(const state_t<R_t, X_t>&, count_t, count_t) = 0;
+};
+
diff --git a/lib/include/wolff/state.hpp b/lib/include/wolff/state.hpp
index e7c0ac3..4909881 100644
--- a/lib/include/wolff/state.hpp
+++ b/lib/include/wolff/state.hpp
@@ -9,26 +9,17 @@
template <class R_t, class X_t>
class state_t {
- private:
- // updating fourier terms F requires many cos and sin calls, faster to do it beforehand.
- std::vector<std::vector<double>> precomputed_cos;
- std::vector<std::vector<double>> precomputed_sin;
public:
- D_t D;
- L_t L;
- v_t nv; // the number of vertices in the lattice
- v_t ne; // the number of edges in the lattice
- graph_t g; // the graph defining the lattice without ghost
+ D_t D; // the dimension of the system
+ L_t L; // the linear size of the lattice
+ v_t nv; // the number of vertices in the original lattice
+ v_t ne; // the number of edges in the original lattice
+ graph_t g; // the graph defining the lattice with ghost
double T; // the temperature
std::vector<X_t> spins; // the state of the ordinary spins
#ifndef NOFIELD
R_t R; // the current state of the ghost site
#endif
- double E; // the system's total energy
- typename X_t::M_t M; // the "sum" of the spins, like the total magnetization
- v_t last_cluster_size; // the size of the last cluster
- std::vector<typename X_t::F_t> ReF;
- std::vector<typename X_t::F_t> ImF;
#ifdef BOND_DEPENDENCE
std::function <double(v_t, const X_t&, v_t, const X_t&)> J; // coupling between sites
@@ -57,7 +48,7 @@ class state_t {
, std::function <double(const X_t&)> H
#endif
#endif
- ) : D(D), L(L), g(D, L), T(T),
+ , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), g(D, L, lat), T(T),
#ifndef NOFIELD
R(),
#endif
@@ -69,69 +60,9 @@ class state_t {
nv = g.nv;
ne = g.ne;
spins.resize(nv);
-#ifdef BOND_DEPENDENCE
- E = 0;
- for (v_t v = 0; v < nv; v++) {
- for (const v_t &vn : g.v_adj[v]) {
- if (v < vn) {
- E -= J(v, spins[v], vn, spins[vn]);
- }
- }
- }
-#else
- E = - (double)ne * J(spins[0], spins[0]);
-#endif
-
#ifndef NOFIELD
g.add_ext();
-#ifdef SITE_DEPENDENCE
- for (v_t i = 0; i < nv; i++) {
- E -= H(i, spins[i]);
- }
-#else
- E -= (double)nv * H(spins[0]);
-#endif
-#endif
-
- M = spins[0] * nv;
- last_cluster_size = 0;
- ReF.resize(D);
- ImF.resize(D);
- for (D_t i = 0; i < D; i++) {
- ReF[i] = spins[0] * 0.0;
- ImF[i] = spins[0] * 0.0;
- }
- precomputed_cos.resize(nv);
- precomputed_sin.resize(nv);
- for (v_t i = 0; i < nv; i++) {
- precomputed_cos[i].resize(D);
- precomputed_sin[i].resize(D);
- for (D_t j = 0; j < D; j++) {
- precomputed_cos[i][j] = cos(2 * M_PI * g.coordinate[i][j] / (double)L);
- precomputed_sin[i][j] = sin(2 * M_PI * g.coordinate[i][j] / (double)L);
- }
- }
- }
-
- void update_magnetization(const X_t& s_old, const X_t& s_new) {
- M += s_new - s_old;
- }
-
- void update_energy(const double& dE) {
- E += dE;
- }
-
- void update_fourierZero(v_t v, const X_t& s_old, const X_t& s_new) {
-#ifdef DIMENSION
- for (D_t i = 0; i < DIMENSION; i++) {
-#else
- for (D_t i = 0; i < D; i++) {
#endif
- ReF[i] += (s_new - s_old) * precomputed_cos[v][i];
- ImF[i] += (s_new - s_old) * precomputed_sin[v][i];
- }
}
};
-
-
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 4043413..c6f0ba6 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -1,25 +1,55 @@
#include <wolff/graph.hpp>
-graph_t::graph_t(D_t D, L_t L) {
- nv = pow(L, D);
- ne = D * nv;
+graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
+ switch (lat) {
+ case SQUARE_LATTICE: {
+ nv = pow(L, D);
+ ne = D * nv;
- v_adj.resize(nv);
- coordinate.resize(nv);
+ v_adj.resize(nv);
+ coordinate.resize(nv);
- for (std::vector<v_t> v_adj_i : v_adj) {
- v_adj_i.reserve(2 * D);
- }
+ for (std::vector<v_t> v_adj_i : v_adj) {
+ v_adj_i.reserve(2 * D);
+ }
- for (v_t i = 0; i < nv; i++) {
- coordinate[i].resize(D);
- for (D_t j = 0; j < D; j++) {
- coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
+ for (v_t i = 0; i < nv; i++) {
+ coordinate[i].resize(D);
+ for (D_t j = 0; j < D; j++) {
+ coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
+
+ v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
+ v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
+ }
+ }
+ break;
+ }
+ case DIAGONAL_LATTICE: {
+ nv = D * pow(L, D);
+ ne = D * nv;
+
+ v_adj.resize(nv);
+ coordinate.resize(nv);
+
+ for (std::vector<v_t> v_adj_i : v_adj) {
+ v_adj_i.reserve(4 * (D - 1));
+ }
+
+ for (D_t i = 0; i < D; i++) {
+ v_t sb = i * pow(L, D);
+
+ for (v_t j = 0; j < pow(L, D); j++) {
+ v_t vc = sb + j;
- v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
- v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
- }
+ v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + j);
+ v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L);
+ v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
+ v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
+ }
+ }
+ break;
+ }
}
}