From 1343a3fe6bd17a2487f12a0d61be8dc83cd722a0 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 15 Oct 2018 22:57:17 -0400 Subject: many changes, including reworking the measurements system --- examples/CMakeLists.txt | 3 + examples/include/correlation.hpp | 10 +- examples/include/measure.hpp | 164 +++++++++++++++++++------ examples/simple_ising.cpp | 192 ++++++++++++++++++++++++++++++ examples/src/models/On/wolff_On.cpp | 25 ++-- examples/src/models/ising/ising.hpp | 33 ++++- examples/src/models/ising/wolff_ising.cpp | 25 ++-- examples/src/models/potts/symmetric.hpp | 1 + examples/src/models/potts/wolff_clock.cpp | 22 ++-- examples/src/models/potts/wolff_potts.cpp | 25 ++-- lib/include/wolff.hpp | 15 +-- lib/include/wolff/cluster.hpp | 20 ++-- lib/include/wolff/graph.hpp | 9 +- lib/include/wolff/meas.h | 19 +++ lib/include/wolff/state.hpp | 81 +------------ lib/src/graph.cpp | 60 +++++++--- 16 files changed, 491 insertions(+), 213 deletions(-) create mode 100644 examples/simple_ising.cpp create mode 100644 lib/include/wolff/meas.h 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 -template -double correlation_length(const state_t & s) { +template +double correlation_length(const std::vector& ReF, const std::vector& 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 +#include #include "correlation.hpp" #include @@ -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 wolff_research_measurements : public wolff_measurement { + private: + std::vector> precomputed_sin; + std::vector> precomputed_cos; + FILE **files; + D_t D; + unsigned char flags; + std::function &, const wolff_research_measurements&)> other_f; + bool silent; + public: + v_t last_cluster_size; + double E; + typename X_t::M_t M; + std::vector ReF; + std::vector ImF; + + wolff_research_measurements(unsigned char flags, unsigned long timestamp, std::function &, const wolff_research_measurements&)> other_f, state_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& 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 -std::function &)> measure_function_write_files(unsigned char flags, FILE **files, std::function &)> other_f) { - return [=] (const state_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& 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(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 +#include + +#include "include/randutils/randutils.hpp" + +#include + +// 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 +class ising_measurements : public wolff_measurement { + 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& 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& 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 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 B = [=] (const ising_t& s) -> double { + if (s.x) { + return -H; + } else { + return H; + } + }; + + // initialize the system + state_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 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(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 other_f; + std::function &)> 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& 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& 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& m) {}; } - std::function measurements = measure_function_write_files(measurement_flags, outfiles, other_f); - std::function H; if (modulated_field) { @@ -251,20 +247,21 @@ int main(int argc, char *argv[]) { state_t s(D, L, T, dot ); #endif + wolff_research_measurements 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 (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); 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 +// 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 &)> other_f; + std::function &, const wolff_research_measurements&)> other_f; uint64_t sum_of_clusterSize = 0; if (N_is_sweeps) { - other_f = [&] (const state_t& s) { - sum_of_clusterSize += s.last_cluster_size; + other_f = [&] (const state_t& s, const wolff_research_measurements& 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 & s) { + other_f = [] (const state_t & s, const wolff_research_measurements& 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& s) {}; + other_f = [] (const state_t& s, const wolff_research_measurements& meas) {}; } - std::function &)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f); + wolff_research_measurements 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 { } } + 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 +#include // hack to speed things up considerably #define N_STATES POTTSQ @@ -95,7 +96,7 @@ int main(int argc, char *argv[]) { std::function (std::mt19937&, potts_t)> gen_R = [] (std::mt19937& r, potts_t v) -> dihedral_t { dihedral_t rot; rot.is_reflection = true; - std::uniform_int_distribution dist(0, POTTSQ - 1); + std::uniform_int_distribution 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 measurement; + std::function , potts_t>&)> 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, potts_t>&) { }; } 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, potts_t>&) { glClear(GL_COLOR_BUFFER_BIT); for (v_t i = 0; i < pow(L, 2); i++) { potts_t 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, potts_t> 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, = %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 (std::mt19937&, potts_t)> gen_R = [] (std::mt19937& r, potts_t v) -> symmetric_t { symmetric_t rot; - std::uniform_int_distribution dist(0, POTTSQ - 1); + std::uniform_int_distribution 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 other_f; + std::function , potts_t>&)> 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, potts_t>& 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, potts_t>& m) { glClear(GL_COLOR_BUFFER_BIT); for (v_t i = 0; i < pow(s.L, 2); i++) { potts_t 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, potts_t>& m) {}; } - std::function measurements = measure_function_write_files(measurement_flags, outfiles, other_f); + wolff_research_measurements, potts_t> 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 -void wolff(count_t N, state_t & s, std::function gen_R, std::function &)> measurements, std::mt19937& r, bool silent) { +void wolff(count_t N, state_t & s, std::function gen_R, wolff_measurement& m, std::mt19937& r) { #ifdef FINITE_STATES #ifdef NOFIELD @@ -15,21 +15,16 @@ void wolff(count_t N, state_t & s, std::function 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 (s, v0, step, r); - measurements(s); - } + m.pre_cluster(s, steps, N, v0, step); + + flip_cluster(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 -void flip_cluster(state_t& s, v_t v0, const R_t& r, std::mt19937& rand) { - std::uniform_real_distribution dist(0.0,1.0); - v_t nv = 0; +void flip_cluster(state_t& s, v_t v0, const R_t& r, std::mt19937& rand, wolff_measurement& m) { + std::uniform_real_distribution dist(0.0, 1.0); std::stack stack; stack.push(v0); @@ -75,8 +75,8 @@ void flip_cluster(state_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& 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& 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 #include -#include #include #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> v_adj; std::vector> 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 wolff_measurement { + public: + virtual void pre_cluster(const state_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&, 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 state_t { - private: - // updating fourier terms F requires many cos and sin calls, faster to do it beforehand. - std::vector> precomputed_cos; - std::vector> 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 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 ReF; - std::vector ImF; #ifdef BOND_DEPENDENCE std::function J; // coupling between sites @@ -57,7 +48,7 @@ class state_t { , std::function 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 -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_adj_i : v_adj) { - v_adj_i.reserve(2 * D); - } + for (std::vector 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_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; + } } } -- cgit v1.2.3-54-g00ecf