From a43ff1f98e9b9814f858bccb11c174b418458491 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 10 Oct 2018 21:45:32 -0400 Subject: big rearrangement of files to make libraries and example (research) files clearer, and changed to c++ std lib random numbers --- .gitmodules | 3 + CMakeLists.txt | 91 +---- examples/CMakeLists.txt | 7 + examples/include/colors.h | 34 ++ examples/include/correlation.hpp | 23 ++ examples/include/measure.hpp | 63 +++ examples/include/randutils | 1 + examples/src/CMakeLists.txt | 4 + examples/src/models/CMakeLists.txt | 6 + examples/src/models/On/CMakeLists.txt | 29 ++ examples/src/models/On/orthogonal.hpp | 202 ++++++++++ examples/src/models/On/vector.hpp | 118 ++++++ examples/src/models/On/wolff_On.cpp | 268 +++++++++++++ examples/src/models/ising/CMakeLists.txt | 28 ++ examples/src/models/ising/ising.hpp | 59 +++ examples/src/models/ising/wolff_ising.cpp | 198 ++++++++++ examples/src/models/ising/z2.hpp | 53 +++ examples/src/models/potts/CMakeLists.txt | 39 ++ examples/src/models/potts/dihedral.hpp | 48 +++ examples/src/models/potts/potts.hpp | 72 ++++ examples/src/models/potts/symmetric.hpp | 51 +++ examples/src/models/potts/wolff_clock.cpp | 155 ++++++++ examples/src/models/potts/wolff_potts.cpp | 213 +++++++++++ examples/src/models/roughening/CMakeLists.txt | 21 + examples/src/models/roughening/dihedral_inf.hpp | 47 +++ examples/src/models/roughening/height.hpp | 75 ++++ examples/src/models/roughening/wolff_cgm.cpp | 167 ++++++++ examples/src/models/roughening/wolff_dgm.cpp | 164 ++++++++ examples/src/tools/CMakeLists.txt | 7 + examples/src/tools/analyze_correlations.cpp | 486 ++++++++++++++++++++++++ lib/CMakeLists.txt | 17 + lib/angle.h | 48 --- lib/circle_group.h | 46 --- lib/cluster.h | 111 ------ lib/colors.h | 34 -- lib/convex.c | 102 ----- lib/convex.h | 23 -- lib/correlation.h | 23 -- lib/dihedral.h | 48 --- lib/dihedral_inf.h | 47 --- lib/finite_states.h | 40 -- lib/graph.cpp | 44 --- lib/graph.h | 22 -- lib/height.h | 75 ---- lib/include/wolff.hpp | 35 ++ lib/include/wolff/cluster.hpp | 111 ++++++ lib/include/wolff/finite_states.hpp | 40 ++ lib/include/wolff/graph.hpp | 21 + lib/include/wolff/state.hpp | 89 +++++ lib/include/wolff/types.h | 36 ++ lib/ising.h | 59 --- lib/measure.h | 63 --- lib/orthogonal.h | 200 ---------- lib/potts.h | 72 ---- lib/rand.c | 20 - lib/rand.h | 16 - lib/src/graph.cpp | 44 +++ lib/state.h | 89 ----- lib/symmetric.h | 51 --- lib/torus.h | 64 ---- lib/types.h | 36 -- lib/vector.h | 118 ------ lib/wolff.h | 33 -- lib/z2.h | 53 --- src/analyze_correlations.cpp | 486 ------------------------ src/wolff_On.cpp | 269 ------------- src/wolff_cgm.cpp | 166 -------- src/wolff_clock.cpp | 154 -------- src/wolff_dgm.cpp | 171 --------- src/wolff_ising.cpp | 201 ---------- src/wolff_potts.cpp | 213 ----------- wolfram_link/Makefile | 13 - wolfram_link/convexminorant.tm | 34 -- 73 files changed, 3039 insertions(+), 3330 deletions(-) create mode 100644 .gitmodules create mode 100644 examples/CMakeLists.txt create mode 100644 examples/include/colors.h create mode 100644 examples/include/correlation.hpp create mode 100644 examples/include/measure.hpp create mode 160000 examples/include/randutils create mode 100644 examples/src/CMakeLists.txt create mode 100644 examples/src/models/CMakeLists.txt create mode 100644 examples/src/models/On/CMakeLists.txt create mode 100644 examples/src/models/On/orthogonal.hpp create mode 100644 examples/src/models/On/vector.hpp create mode 100644 examples/src/models/On/wolff_On.cpp create mode 100644 examples/src/models/ising/CMakeLists.txt create mode 100644 examples/src/models/ising/ising.hpp create mode 100644 examples/src/models/ising/wolff_ising.cpp create mode 100644 examples/src/models/ising/z2.hpp create mode 100644 examples/src/models/potts/CMakeLists.txt create mode 100644 examples/src/models/potts/dihedral.hpp create mode 100644 examples/src/models/potts/potts.hpp create mode 100644 examples/src/models/potts/symmetric.hpp create mode 100644 examples/src/models/potts/wolff_clock.cpp create mode 100644 examples/src/models/potts/wolff_potts.cpp create mode 100644 examples/src/models/roughening/CMakeLists.txt create mode 100644 examples/src/models/roughening/dihedral_inf.hpp create mode 100644 examples/src/models/roughening/height.hpp create mode 100644 examples/src/models/roughening/wolff_cgm.cpp create mode 100644 examples/src/models/roughening/wolff_dgm.cpp create mode 100644 examples/src/tools/CMakeLists.txt create mode 100644 examples/src/tools/analyze_correlations.cpp create mode 100644 lib/CMakeLists.txt delete mode 100644 lib/angle.h delete mode 100644 lib/circle_group.h delete mode 100644 lib/cluster.h delete mode 100644 lib/colors.h delete mode 100644 lib/convex.c delete mode 100644 lib/convex.h delete mode 100644 lib/correlation.h delete mode 100644 lib/dihedral.h delete mode 100644 lib/dihedral_inf.h delete mode 100644 lib/finite_states.h delete mode 100644 lib/graph.cpp delete mode 100644 lib/graph.h delete mode 100644 lib/height.h create mode 100644 lib/include/wolff.hpp create mode 100644 lib/include/wolff/cluster.hpp create mode 100644 lib/include/wolff/finite_states.hpp create mode 100644 lib/include/wolff/graph.hpp create mode 100644 lib/include/wolff/state.hpp create mode 100644 lib/include/wolff/types.h delete mode 100644 lib/ising.h delete mode 100644 lib/measure.h delete mode 100644 lib/orthogonal.h delete mode 100644 lib/potts.h delete mode 100644 lib/rand.c delete mode 100644 lib/rand.h create mode 100644 lib/src/graph.cpp delete mode 100644 lib/state.h delete mode 100644 lib/symmetric.h delete mode 100644 lib/torus.h delete mode 100644 lib/types.h delete mode 100644 lib/vector.h delete mode 100644 lib/wolff.h delete mode 100644 lib/z2.h delete mode 100644 src/analyze_correlations.cpp delete mode 100644 src/wolff_On.cpp delete mode 100644 src/wolff_cgm.cpp delete mode 100644 src/wolff_clock.cpp delete mode 100644 src/wolff_dgm.cpp delete mode 100644 src/wolff_ising.cpp delete mode 100644 src/wolff_potts.cpp delete mode 100644 wolfram_link/Makefile delete mode 100644 wolfram_link/convexminorant.tm diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..66f0577 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "examples/include/randutils"] + path = examples/include/randutils + url = https://gist.github.com/imneme/540829265469e673d045 diff --git a/CMakeLists.txt b/CMakeLists.txt index 7bf604c..073ae14 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ -cmake_minimum_required(VERSION 3.0) -project(wolff) +cmake_minimum_required(VERSION 3.9) +project(wolff DESCRIPTION "a library for preforming the wolff algorithm") set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall") set(CMAKE_CXX_FLAGS_RELEASE "-O3") @@ -8,89 +8,8 @@ set(CMAKE_CXX_FLAGS_RELEASE "-O3") set (CMAKE_CXX_STANDARD 17) set (CMAKE_C_STANDARD 11) -include_directories(lib ~/.local/include) -link_directories(~/.local/lib) +include(GNUInstallDirs) -file(GLOB CSOURCES lib/*.c) -file(GLOB CPPSOURCES lib/*.cpp) - -add_executable(wolff_ising src/wolff_ising.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_ising_2D_0 src/wolff_ising.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_dgm src/wolff_dgm.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_cgm src/wolff_cgm.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_3potts src/wolff_potts.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_4potts src/wolff_potts.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_7potts src/wolff_potts.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_3clock src/wolff_clock.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_5clock src/wolff_clock.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_planar src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_planar_2D src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_planar_2D_0 src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_heisenberg src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(analyze_correlations src/analyze_correlations.cpp ${CPPSOURCES} ${CSOURCES}) - -SET_TARGET_PROPERTIES(wolff_ising_2D_0 PROPERTIES COMPILE_FLAGS "-DDIMENSION=2 -DNOFIELD") -SET_TARGET_PROPERTIES(wolff_3potts PROPERTIES COMPILE_FLAGS "-DPOTTSQ=3") -SET_TARGET_PROPERTIES(wolff_4potts PROPERTIES COMPILE_FLAGS "-DPOTTSQ=4") -SET_TARGET_PROPERTIES(wolff_7potts PROPERTIES COMPILE_FLAGS "-DPOTTSQ=7") -SET_TARGET_PROPERTIES(wolff_3clock PROPERTIES COMPILE_FLAGS "-DPOTTSQ=3") -SET_TARGET_PROPERTIES(wolff_5clock PROPERTIES COMPILE_FLAGS "-DPOTTSQ=5") -SET_TARGET_PROPERTIES(wolff_planar PROPERTIES COMPILE_FLAGS "-DN_COMP=2") -SET_TARGET_PROPERTIES(wolff_planar_2D PROPERTIES COMPILE_FLAGS "-DN_COMP=2 -DDIMENSION=2") -SET_TARGET_PROPERTIES(wolff_planar_2D_0 PROPERTIES COMPILE_FLAGS "-DN_COMP=2 -DDIMENSION=2 -DNOFIELD") -SET_TARGET_PROPERTIES(wolff_heisenberg PROPERTIES COMPILE_FLAGS "-DN_COMP=3") - -find_library(GSL REQUIRED NAMES gsl) -find_library(FFTW REQUIRED NAMES fftw3) -find_library(M REQUIRED NAMES m) -FIND_LIBRARY(GL NAMES GL) -FIND_LIBRARY(GLU NAMES GLU) -FIND_LIBRARY(GLUT NAMES glut) - -target_link_libraries(analyze_correlations cblas gsl fftw3 m) - -if (${GLUT} MATCHES "GLUT-NOTFOUND") - target_link_libraries(wolff_ising cblas gsl m) - target_link_libraries(wolff_ising_2D_0 cblas gsl m) - target_link_libraries(wolff_dgm cblas gsl m) - target_link_libraries(wolff_cgm cblas gsl m) - target_link_libraries(wolff_3potts cblas gsl m) - target_link_libraries(wolff_4potts cblas gsl m) - target_link_libraries(wolff_7potts cblas gsl m) - target_link_libraries(wolff_3clock cblas gsl m) - target_link_libraries(wolff_5clock cblas gsl m) - target_link_libraries(wolff_heisenberg cblas gsl m) - target_link_libraries(wolff_planar cblas gsl m) - target_link_libraries(wolff_planar_2D cblas gsl m) - target_link_libraries(wolff_planar_2D_0 cblas gsl m) -else() - target_link_libraries(wolff_ising cblas gsl m glut GL GLU) - target_link_libraries(wolff_ising_2D_0 cblas gsl m glut GL GLU) - target_link_libraries(wolff_dgm cblas gsl m glut GL GLU) - target_link_libraries(wolff_cgm cblas gsl m glut GL GLU) - target_link_libraries(wolff_3potts cblas gsl m glut GL GLU) - target_link_libraries(wolff_4potts cblas gsl m glut GL GLU) - target_link_libraries(wolff_7potts cblas gsl m glut GL GLU) - target_link_libraries(wolff_3clock cblas gsl m glut GL GLU) - target_link_libraries(wolff_5clock cblas gsl m glut GL GLU) - target_link_libraries(wolff_heisenberg cblas gsl m glut GL GLU) - target_link_libraries(wolff_planar cblas gsl m glut GL GLU) - target_link_libraries(wolff_planar_2D cblas gsl m glut GL GLU) - target_link_libraries(wolff_planar_2D_0 cblas gsl m glut GL GLU) - target_compile_definitions(wolff_ising PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_ising_2D_0 PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_dgm PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_cgm PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_3potts PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_4potts PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_7potts PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_3clock PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_5clock PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_planar PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_planar_2D PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_planar_2D_0 PUBLIC HAVE_GLUT) - target_compile_definitions(wolff_heisenberg PUBLIC HAVE_GLUT) -endif() - -install(TARGETS wolff_ising wolff_dgm wolff_cgm wolff_3potts wolff_4potts wolff_3clock wolff_heisenberg wolff_planar analyze_correlations DESTINATION bin) +add_subdirectory(lib) +add_subdirectory(examples) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt new file mode 100644 index 0000000..baf0a10 --- /dev/null +++ b/examples/CMakeLists.txt @@ -0,0 +1,7 @@ + +add_library(wolff_examples INTERFACE) + +target_include_directories(wolff_examples INTERFACE include) + +add_subdirectory(src) + diff --git a/examples/include/colors.h b/examples/include/colors.h new file mode 100644 index 0000000..abf137c --- /dev/null +++ b/examples/include/colors.h @@ -0,0 +1,34 @@ +#pragma once + +#include + +double hue_to_R(double theta) { + if (((M_PI / 3 <= theta) && (theta < 2 * M_PI / 3)) || ((4 * M_PI / 3 <= theta) && (theta < 5 * M_PI / 3))) { + return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); + } else if (((0 <= theta) && (theta < M_PI / 3)) || ((5 * M_PI / 3 <= theta) && (theta <= 2 * M_PI))) { + return 1.0; + } else { + return 0.0; + } +} + +double hue_to_G(double theta) { + if (((0 <= theta) && (theta < M_PI / 3)) || ((M_PI <= theta) && (theta < 4 * M_PI / 3))) { + return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); + } else if (((M_PI / 3 <= theta) && (theta < 2 * M_PI / 3)) || ((2 * M_PI / 3 <= theta) && (theta < M_PI))) { + return 1.0; + } else { + return 0.0; + } +} + +double hue_to_B(double theta) { + if (((2 * M_PI / 3 <= theta) && (theta < M_PI)) || ((5 * M_PI / 3 <= theta) && (theta <= 2 * M_PI))) { + return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); + } else if (((M_PI <= theta) && (theta < 4 * M_PI / 3)) || ((4 * M_PI / 3 <= theta) && (theta < 5 * M_PI / 3))) { + return 1.0; + } else { + return 0.0; + } +} + diff --git a/examples/include/correlation.hpp b/examples/include/correlation.hpp new file mode 100644 index 0000000..042cff3 --- /dev/null +++ b/examples/include/correlation.hpp @@ -0,0 +1,23 @@ + +#pragma once + +#include +#include + +#include + +template +double correlation_length(const state_t & s) { + double total = 0; + +#ifdef DIMENSION + for (D_t j = 0; j < DIMENSION; j++) { +#else + for (D_t j = 0; j < s.D; j++) { +#endif + total += norm_squared(s.ReF[j]) + norm_squared(s.ImF[j]); + } + + return total / s.D; +} + diff --git a/examples/include/measure.hpp b/examples/include/measure.hpp new file mode 100644 index 0000000..e20353c --- /dev/null +++ b/examples/include/measure.hpp @@ -0,0 +1,63 @@ + +#pragma once + +#include +#include "correlation.hpp" +#include + +#define POSSIBLE_MEASUREMENTS 4 +const unsigned char measurement_energy = 1 << 0; +const unsigned char measurement_clusterSize = 1 << 1; +const unsigned char measurement_magnetization = 1 << 2; +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 *)); + + 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; +} + +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]); + } + if (flags & measurement_magnetization) { + write_magnetization(s.M, files[2]); + } + if (flags & measurement_fourierZero) { + float smaller_X = (float)correlation_length(s); + fwrite(&smaller_X, sizeof(float), 1, files[3]); + } + + other_f(s); + }; +} + +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); +} + + diff --git a/examples/include/randutils b/examples/include/randutils new file mode 160000 index 0000000..8486a61 --- /dev/null +++ b/examples/include/randutils @@ -0,0 +1 @@ +Subproject commit 8486a610a954a8248c12485fb4cfc390a5f5f854 diff --git a/examples/src/CMakeLists.txt b/examples/src/CMakeLists.txt new file mode 100644 index 0000000..3397426 --- /dev/null +++ b/examples/src/CMakeLists.txt @@ -0,0 +1,4 @@ + +add_subdirectory(models) +add_subdirectory(tools) + diff --git a/examples/src/models/CMakeLists.txt b/examples/src/models/CMakeLists.txt new file mode 100644 index 0000000..0b0c111 --- /dev/null +++ b/examples/src/models/CMakeLists.txt @@ -0,0 +1,6 @@ + +add_subdirectory(ising) +add_subdirectory(On) +add_subdirectory(potts) +add_subdirectory(roughening) + diff --git a/examples/src/models/On/CMakeLists.txt b/examples/src/models/On/CMakeLists.txt new file mode 100644 index 0000000..26985b9 --- /dev/null +++ b/examples/src/models/On/CMakeLists.txt @@ -0,0 +1,29 @@ + +add_executable(wolff_planar wolff_On.cpp) +add_executable(wolff_planar_2d_no-field wolff_On.cpp) +add_executable(wolff_heisenberg wolff_On.cpp) + +set_target_properties(wolff_planar PROPERTIES COMPILE_FLAGS "-DN_COMP=2") +set_target_properties(wolff_planar_2d_no-field PROPERTIES COMPILE_FLAGS "-DN_COMP=2 -DDIMENSION=2 -DNOFIELD") +set_target_properties(wolff_heisenberg PROPERTIES COMPILE_FLAGS "-DN_COMP=3") + +find_library(GL NAMES GL) +find_library(GLU NAMES GLU) +find_library(GLUT NAMES glut) + +if (${GLUT} MATCHES "GLUT-NOTFOUND") + target_link_libraries(wolff_planar wolff wolff_examples) + target_link_libraries(wolff_planar_2d_no-field wolff wolff_examples) + target_link_libraries(wolff_heisenberg wolff wolff_examples) +else() + target_compile_definitions(wolff_planar PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_planar_2d_no-field PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_heisenberg PUBLIC HAVE_GLUT) + + target_link_libraries(wolff_planar wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_planar_2d_no-field wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_heisenberg wolff wolff_examples glut GL GLU) +endif() + +install(TARGETS wolff_planar wolff_planar_2d_no-field wolff_heisenberg DESTINATION ${CMAKE_INSTALL_BINDIR}) + diff --git a/examples/src/models/On/orthogonal.hpp b/examples/src/models/On/orthogonal.hpp new file mode 100644 index 0000000..f13357f --- /dev/null +++ b/examples/src/models/On/orthogonal.hpp @@ -0,0 +1,202 @@ + +#pragma once + +#include +#include +#include + +#include +#include +#include "vector.hpp" + +template +class orthogonal_t : public std::array, q> { + public : + bool is_reflection; + + orthogonal_t() : is_reflection(false) { + for (q_t i = 0; i < q; i++) { + (*this)[i].fill(0); + (*this)[i][i] = (T)1; + } + } + + vector_t act(const vector_t & v) const { + vector_t v_rot; + v_rot.fill(0); + + if (is_reflection) { + double prod = 0; + for (q_t i = 0; i < q; i++) { + prod += v[i] * (*this)[0][i]; + } + for (q_t i = 0; i < q; i++) { + v_rot[i] = v[i] - 2 * prod * (*this)[0][i]; + } + } else { + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot[i] += (*this)[i][j] * v[j]; + } + } + } + + return v_rot; + } + + orthogonal_t act(const orthogonal_t & m) const { + orthogonal_t m_rot; + + m_rot.is_reflection = false; + + if (is_reflection) { + for (q_t i = 0; i < q; i++) { + double akOki = 0; + + for (q_t k = 0; k < q; k++) { + akOki += (*this)[0][k] * m[k][i]; + } + + for (q_t j = 0; j < q; j++) { + m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][j]; + } + } + } else { + for (q_t i = 0; i < q; i++) { + m_rot[i].fill(0); + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m_rot[i][j] += (*this)[i][j] * m[j][k]; + } + } + } + } + + return m_rot; + } + + vector_t act_inverse(const vector_t & v) const { + if (is_reflection) { + return this->act(v); // reflections are their own inverse + } else { + vector_t v_rot; + v_rot.fill(0); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot[i] += (*this)[j][i] * v[j]; + } + } + + return v_rot; + } + } + + vector_t act_inverse(const orthogonal_t & m) const { + if (is_reflection) { + return this->act(m); // reflections are their own inverse + } else { + orthogonal_t m_rot; + m_rot.is_reflection = false; + + for (q_t i = 0; i < q; i++) { + m_rot[i].fill(0); + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m_rot[i][j] += (*this)[j][i] * m[j][k]; + } + } + } + + return m_rot; + } + } + +}; + + +template +orthogonal_t generate_rotation_uniform (std::mt19937& r, const vector_t & v) { + std::normal_distribution dist(0.0,1.0); + orthogonal_t ptr; + ptr.is_reflection = true; + + double v2 = 0; + + for (q_t i = 0; i < q; i++) { + ptr[0][i] = dist(r); + v2 += ptr[0][i] * ptr[0][i]; + } + + double mag_v = sqrt(v2); + + for (q_t i = 0; i < q; i++) { + ptr[0][i] /= mag_v; + } + + return ptr; +} + +template +orthogonal_t generate_rotation_perturbation (std::mt19937& r, const vector_t & v0, double epsilon, unsigned int n) { + std::normal_distribution dist(0.0,1.0); + orthogonal_t m; + m.is_reflection = true; + + vector_t v; + + if (n > 1) { + std::uniform_int_distribution udist(0, n); + unsigned int rotation = udist(r); + + double cosr = cos(2 * M_PI * rotation / (double)n / 2.0); + double sinr = sin(2 * M_PI * rotation / (double)n / 2.0); + + v[0] = v0[0] * cosr - v0[1] * sinr; + v[1] = v0[1] * cosr + v0[0] * sinr; + + for (q_t i = 2; i < q; i++) { + v[i] = v0[i]; + } + } else { + v = v0; + } + + double m_dot_v = 0; + + for (q_t i = 0; i < q; i++) { + m[0][i] = dist(r); // create a random vector + m_dot_v += m[0][i] * v[i]; + } + + double v2 = 0; + + for (q_t i = 0; i < q; i++) { + m[0][i] = m[0][i] - m_dot_v * v[i]; // find the component orthogonal to v + v2 += pow(m[0][i], 2); + } + + double mag_v = sqrt(v2); + + for (q_t i = 0; i < q; i++) { + m[0][i] /= mag_v; // normalize + } + + v2 = 0; + + double factor = epsilon * dist(r); + + for (q_t i = 0; i < q; i++) { + m[0][i] += factor * v[i]; // perturb orthogonal vector in original direction + v2 += pow(m[0][i], 2); + } + + mag_v = sqrt(v2); + + for (q_t i = 0; i < q; i++) { + m[0][i] /= mag_v; // normalize + } + + return m; +} + diff --git a/examples/src/models/On/vector.hpp b/examples/src/models/On/vector.hpp new file mode 100644 index 0000000..1cdb60a --- /dev/null +++ b/examples/src/models/On/vector.hpp @@ -0,0 +1,118 @@ + +#pragma once + +#include +#include +#include + +#include + +template +class vector_t : public std::array { + public: + + // M_t needs to hold the sum of nv spins + typedef vector_t M_t; + + // F_t needs to hold the double-weighted sum of spins + typedef vector_t F_t; + + vector_t() { + this->fill((T)0); + (*this)[1] = (T)1; + } + + vector_t(const T *x) { + for (q_t i = 0; i < q; i++) { + (*this)[i] = x[i]; + } + } + + template + inline vector_t& operator+=(const vector_t &v) { + for (q_t i = 0; i < q; i++) { + (*this)[i] += (U)v[i]; + } + return *this; + } + + template + inline vector_t& operator-=(const vector_t &v) { + for (q_t i = 0; i < q; i++) { + (*this)[i] -= (U)v[i]; + } + return *this; + } + + inline vector_t operator*(v_t x) const { + vector_t result; + for (q_t i = 0; i < q; i++) { + result[i] = x * (*this)[i]; + } + + return result; + } + + inline vector_t operator*(double x) const { + vector_t result; + for (q_t i = 0; i < q; i++) { + result[i] = x * (*this)[i]; + } + + return result; + } + + inline vector_t operator-(const vector_t& v) const { + vector_t diff = *this; + diff -= v; + return diff; + } +}; + + +template +double norm_squared(vector_t v) { + double tmp = 0; + for (T &x : v) { + tmp += pow(x, 2); + } + + return tmp; +} + +template +void write_magnetization(vector_t M, FILE *outfile) { + for (q_t i = 0; i < q; i++) { + fwrite(&(M[i]), sizeof(T), q, outfile); + } +} + +// below functions and definitions are unnecessary for wolff.h but useful. + +template // save some space and don't write whole doubles +void write_magnetization(vector_t M, FILE *outfile) { + for (q_t i = 0; i < q; i++) { + float M_tmp = (float)M[i]; + fwrite(&M_tmp, sizeof(float), 1, outfile); + } +} + +template +T dot(const vector_t & v1, const vector_t & v2) { + T prod = 0; + + for (q_t i = 0; i < q; i++) { + prod += v1[i] * v2[i]; + } + + return prod; +} + +template +double H_vector(const vector_t & v1, T *H) { + vector_t H_vec(H); + return (double)(dot (v1, H_vec)); +} + +char const *ON_strings[] = {"TRIVIAL", "ISING", "PLANAR", "HEISENBERG"}; + diff --git a/examples/src/models/On/wolff_On.cpp b/examples/src/models/On/wolff_On.cpp new file mode 100644 index 0000000..e3568c7 --- /dev/null +++ b/examples/src/models/On/wolff_On.cpp @@ -0,0 +1,268 @@ + +#include +#include + +#ifdef HAVE_GLUT +#include +#endif + +#include "orthogonal.hpp" +#include "vector.hpp" + +#include +#include +#include +#include + +typedef orthogonal_t orthogonal_R_t; +typedef vector_t vector_R_t; +typedef state_t On_t; + +// angle from the x-axis of a two-vector +double theta(vector_R_t v) { + double x = v[0]; + double y = v[1]; + + double val = atan(y / x); + + if (x < 0.0 && y > 0.0) { + return M_PI + val; + } else if ( x < 0.0 && y < 0.0 ) { + return - M_PI + val; + } else { + return val; + } +} + +double H_modulated(vector_R_t v, int order, double mag) { + return mag * cos(order * theta(v)); +} + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e7; + +#ifdef DIMENSION + D_t D = DIMENSION; +#else + D_t D = 2; +#endif + L_t L = 128; + double T = 2.26918531421; + double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); + + bool silent = false; + bool use_pert = false; + bool N_is_sweeps = false; + bool draw = false; + unsigned int window_size = 512; + + bool modulated_field = false; + unsigned int order = 1; + + int opt; + q_t H_ind = 0; + double epsilon = 1; + +// unsigned char measurement_flags = measurement_energy | measurement_clusterSize; + + unsigned char measurement_flags = 0; + + while ((opt = getopt(argc, argv, "N:D:L:T:H:spe:mo:M:Sdw:")) != -1) { + switch (opt) { + case 'N': // number of steps + N = (count_t)atof(optarg); + break; +#ifdef DIMENSION + case 'D': // dimension + printf("Dimension was specified at compile time, you can't change it now!\n"); + exit(EXIT_FAILURE); +#else + case 'D': // dimension + D = atoi(optarg); + break; +#endif + case 'L': // linear size + L = atoi(optarg); + break; + case 'T': // temperature + T = atof(optarg); + break; + case 'H': // external field. nth call couples to state n + H_vec[H_ind] = atof(optarg); + H_ind++; + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + case 'p': + use_pert = true; + break; + case 'e': + epsilon = atof(optarg); + break; + case 'm': + modulated_field = true; + break; + case 'M': + measurement_flags ^= 1 << atoi(optarg); + break; + case 'o': + order = atoi(optarg); + break; + case 'S': + N_is_sweeps = true; + break; + case 'd': +#ifdef HAVE_GLUT + draw = true; + break; +#else + printf("You didn't compile this with the glut library installed!\n"); + exit(EXIT_FAILURE); +#endif + case 'w': + window_size = atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + const char *pert_type; + + std::function gen_R; + + if (use_pert) { + double Hish; + if (modulated_field) { + Hish = fabs(H_vec[0]); + } else { + double H2 = 0; + for (q_t i = 0; i < N_COMP; i++) { + H2 += pow(H_vec[i], 2); + } + Hish = sqrt(H2); + } + + epsilon = sqrt((N_COMP - 1) * T / (D + Hish / 2)) / 2; + + gen_R = std::bind(generate_rotation_perturbation , std::placeholders::_1, std::placeholders::_2, epsilon, order); + pert_type = "PERTURB5"; + } else { + gen_R = generate_rotation_uniform ; + pert_type = "UNIFORM"; + } + + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); + + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"%s\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"FIELD_TYPE\" -> ", timestamp, ON_strings[N_COMP], N_COMP, D, L, (v_t)pow(L, D), D * (v_t)pow(L, D), T); + if (modulated_field) { + fprintf(outfile_info, "\"MODULATED\", \"ORDER\" -> %d, \"H\" -> %.15f, ", order, H_vec[0]); + } else { + fprintf(outfile_info, "\"VECTOR\", \"H\" -> {"); + for (q_t i = 0; i < N_COMP; i++) { + fprintf(outfile_info, "%.15f", H_vec[i]); + if (i < N_COMP - 1) { + fprintf(outfile_info, ", "); + } + } + fprintf(outfile_info, "}, "); + } + + fprintf(outfile_info, "\"GENERATOR\" -> \"%s\"", pert_type); + + if (use_pert) { + fprintf(outfile_info, ", \"EPS\" -> %g", epsilon); + } + + fprintf(outfile_info, " |>\n"); + + fclose(outfile_info); + + FILE **outfiles = measure_setup_files(measurement_flags, timestamp); + + 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; + }; + } else if (draw) { +#ifdef HAVE_GLUT + // initialize glut + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("wolff"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0.0, L, 0.0, L); + + other_f = [&] (const On_t& s) { + glClear(GL_COLOR_BUFFER_BIT); + for (v_t i = 0; i < pow(L, 2); i++) { + vector_R_t v_tmp = s.R.act_inverse(s.spins[i]); + double thetai = fmod(2 * M_PI + theta(v_tmp), 2 * M_PI); + double saturation = 0.7; + double value = 0.9; + double chroma = saturation * value; + glColor3f(chroma * hue_to_R(thetai) + (value - chroma), chroma * hue_to_G(thetai) + (value - chroma), chroma * hue_to_B(thetai) + (value - chroma)); + glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); + } + glFlush(); + }; +#endif + } else { + other_f = [] (const On_t& s) {}; + } + + std::function measurements = measure_function_write_files(measurement_flags, outfiles, other_f); + + std::function H; + + if (modulated_field) { + H = std::bind(H_modulated, std::placeholders::_1, order, H_vec[0]); + } else { + H = std::bind(H_vector , std::placeholders::_1, H_vec); + } + + // initialize random number generator + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + +#ifndef NOFIELD + state_t s(D, L, T, dot , H); +#else + state_t s(D, L, T, dot ); +#endif + + 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); + 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); + } else { + wolff (N, s, gen_R, measurements, rng, silent); + } + + measure_free_files(measurement_flags, outfiles); + free(H_vec); + + return 0; +} + diff --git a/examples/src/models/ising/CMakeLists.txt b/examples/src/models/ising/CMakeLists.txt new file mode 100644 index 0000000..e8fbc9a --- /dev/null +++ b/examples/src/models/ising/CMakeLists.txt @@ -0,0 +1,28 @@ + +add_executable(wolff_ising wolff_ising.cpp) +add_executable(wolff_ising_2d wolff_ising.cpp) +add_executable(wolff_ising_2d_no-field wolff_ising.cpp) + +set_target_properties(wolff_ising_2d PROPERTIES COMPILE_FLAGS "-DDIMENSION=2") +set_target_properties(wolff_ising_2d_no-field PROPERTIES COMPILE_FLAGS "-DDIMENSION=2 -DNOFIELD") + +find_library(GL NAMES GL) +find_library(GLU NAMES GLU) +find_library(GLUT NAMES glut) + +if (${GLUT} MATCHES "GLUT-NOTFOUND") + target_link_libraries(wolff_ising wolff wolff_examples) + target_link_libraries(wolff_ising_2d wolff wolff_examples) + target_link_libraries(wolff_ising_2d_no-field wolff wolff_examples) +else() + target_compile_definitions(wolff_ising PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_ising_2d PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_ising_2d_no-field PUBLIC HAVE_GLUT) + + target_link_libraries(wolff_ising wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_ising_2d wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_ising_2d_no-field wolff wolff_examples glut GL GLU) +endif() + +install(TARGETS wolff_ising wolff_ising_2d wolff_ising_2d_no-field DESTINATION ${CMAKE_INSTALL_BINDIR}) + diff --git a/examples/src/models/ising/ising.hpp b/examples/src/models/ising/ising.hpp new file mode 100644 index 0000000..ae20840 --- /dev/null +++ b/examples/src/models/ising/ising.hpp @@ -0,0 +1,59 @@ +#pragma once + +#include +#include + +#include + +class ising_t { + public: + bool x; + + typedef int M_t; + typedef double F_t; + + ising_t() : x(false) {} + ising_t(bool x) : x(x) {} + ising_t(int x) : x((bool)x) {} + + inline int operator*(v_t a) const { + if (x) { + return -(int)a; + } else { + return (int)a; + } + } + + inline double operator*(double a) const { + if (x) { + return -a; + } else { + return a; + } + } + + inline int operator-(const ising_t &s) const { + if (x == s.x) { + return 0; + } else { + if (x) { + return -2; + } else { + return 2; + } + } + } +}; + +double norm_squared(double s) { + return pow(s, 2); +} + +void write_magnetization(int M, FILE *outfile) { + fwrite(&M, sizeof(int), 1, outfile); +} + +#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 new file mode 100644 index 0000000..5bdaa82 --- /dev/null +++ b/examples/src/models/ising/wolff_ising.cpp @@ -0,0 +1,198 @@ + +#include +#include + +// if you have GLUT installed, you can see graphics! +#ifdef HAVE_GLUT +#include +#endif + +// include your group and spin space +#include "z2.hpp" +#include "ising.hpp" + +// finite_states.h can be included for spin types that have special variables +// defined, and it causes wolff execution to use precomputed bond probabilities +#include + +#include + +// measure.hpp contains useful functions for saving timeseries to files +#include + +// include wolff.hpp +#include + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e4; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double H = 0.0; + + bool silent = false; + bool draw = false; + bool N_is_sweeps = false; + unsigned int window_size = 512; + + // don't measure anything by default + unsigned char measurement_flags = 0; + + int opt; + + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -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; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + case 'S': + N_is_sweeps = true; + break; + case 'd': +#ifdef HAVE_GLUT + draw = true; + break; +#else + printf("You didn't compile this with the glut library installed!\n"); + exit(EXIT_FAILURE); +#endif + case 'w': + window_size = atoi(optarg); + break; + case 'M': + measurement_flags ^= 1 << atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + // get nanosecond timestamp for unique run id + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + // initialize random number generator + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + // define 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 spin-field coupling + std::function B = [=] (const ising_t& s) -> double { + if (s.x) { + return -H; + } else { + return H; + } + }; + + // initialize state object +#ifndef NOFIELD + state_t s(D, L, T, Z, B); +#else + state_t s(D, L, T, Z); +#endif + + // define function that generates self-inverse rotations + std::function gen_R = [] (std::mt19937&, const ising_t& s) -> z2_t { + return z2_t(true); + }; + + FILE **outfiles = measure_setup_files(measurement_flags, timestamp); + + std::function &)> 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; + }; + } else if (draw) { +#ifdef HAVE_GLUT + // initialize glut + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("wolff"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0.0, L, 0.0, L); + + other_f = [] (const state_t & s) { + glClear(GL_COLOR_BUFFER_BIT); + for (v_t i = 0; i < pow(s.L, 2); i++) { + if (s.spins[i].x == s.R.x) { + glColor3f(0.0, 0.0, 0.0); + } else { + glColor3f(1.0, 1.0, 1.0); + } + glRecti(i / s.L, i % s.L, (i / s.L) + 1, (i % s.L) + 1); + } + glFlush(); + }; +#endif + } else { + other_f = [] (const state_t& s) {}; + } + + std::function &)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f); + + // add line to metadata file with run info + { + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); + + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"ISING\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> %.15f |>\n", timestamp, s.D, s.L, s.nv, s.ne, T, H); + + fclose(outfile_info); + } + + // run wolff for N cluster flips + 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); + 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); + } else { + wolff(N, s, gen_R, measurements, rng, silent); + } + + measure_free_files(measurement_flags, outfiles); + + return 0; + +} + diff --git a/examples/src/models/ising/z2.hpp b/examples/src/models/ising/z2.hpp new file mode 100644 index 0000000..19b6c05 --- /dev/null +++ b/examples/src/models/ising/z2.hpp @@ -0,0 +1,53 @@ + +#pragma once + +#include +#include "ising.hpp" + +/* The minimum definition for a group type R_t to act on a spin type X_t is + * given by the following. + * + * void init(R_t *p); + * void free_spin(R_t r); + * R_t copy(R_t r); + * X_t act(R_t r, X_t x); + * R_t act(R_t r, R_t r); + * X_t act_inverse(R_t r, X_t x); + * R_t act_inverse(R_t r, R_t r); + * + */ + +class z2_t { + public: + bool x; + + z2_t() : x(false) {} + + z2_t(bool x) : x(x) {} + + ising_t act(const ising_t& s) const { + if (x) { + return ising_t(!s.x); + } else { + return ising_t(s.x); + } + } + + z2_t act(const z2_t& r) const { + if (x) { + return z2_t(!r.x); + } else { + return z2_t(r.x); + } + } + + ising_t act_inverse(const ising_t& s) const { + return this->act(s); + } + + z2_t act_inverse(const z2_t& r) const { + return this->act(r); + } +}; + + diff --git a/examples/src/models/potts/CMakeLists.txt b/examples/src/models/potts/CMakeLists.txt new file mode 100644 index 0000000..53133b9 --- /dev/null +++ b/examples/src/models/potts/CMakeLists.txt @@ -0,0 +1,39 @@ + +add_executable(wolff_3potts wolff_potts.cpp) +add_executable(wolff_4potts wolff_potts.cpp) +add_executable(wolff_7potts wolff_potts.cpp) +add_executable(wolff_3clock wolff_clock.cpp) +add_executable(wolff_5clock wolff_clock.cpp) + +set_target_properties(wolff_3potts PROPERTIES COMPILE_FLAGS "-DPOTTSQ=3") +set_target_properties(wolff_4potts PROPERTIES COMPILE_FLAGS "-DPOTTSQ=4") +set_target_properties(wolff_7potts PROPERTIES COMPILE_FLAGS "-DPOTTSQ=7") +set_target_properties(wolff_3clock PROPERTIES COMPILE_FLAGS "-DPOTTSQ=3") +set_target_properties(wolff_5clock PROPERTIES COMPILE_FLAGS "-DPOTTSQ=5") + +find_library(GL NAMES GL) +find_library(GLU NAMES GLU) +find_library(GLUT NAMES glut) + +if (${GLUT} MATCHES "GLUT-NOTFOUND") + target_link_libraries(wolff_3potts wolff wolff_examples) + target_link_libraries(wolff_4potts wolff wolff_examples) + target_link_libraries(wolff_7potts wolff wolff_examples) + target_link_libraries(wolff_3clock wolff wolff_examples) + target_link_libraries(wolff_5clock wolff wolff_examples) +else() + target_compile_definitions(wolff_3potts PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_4potts PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_7potts PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_3clock PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_5clock PUBLIC HAVE_GLUT) + + target_link_libraries(wolff_3potts wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_4potts wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_7potts wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_3clock wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_5clock wolff wolff_examples glut GL GLU) +endif() + +install(TARGETS wolff_3potts wolff_4potts wolff_7potts wolff_3clock wolff_5clock DESTINATION ${CMAKE_INSTALL_BINDIR}) + diff --git a/examples/src/models/potts/dihedral.hpp b/examples/src/models/potts/dihedral.hpp new file mode 100644 index 0000000..cbc5687 --- /dev/null +++ b/examples/src/models/potts/dihedral.hpp @@ -0,0 +1,48 @@ + +#pragma once + +#include +#include "potts.hpp" + +template +class dihedral_t { + public: + bool is_reflection; + q_t x; + + dihedral_t() : is_reflection(false), x(0) {} + dihedral_t(bool x, q_t y) : is_reflection(x), x(y) {} + + potts_t act(const potts_t& s) const { + if (this->is_reflection) { + return potts_t(((q + this->x) - s.x) % q); + } else { + return potts_t((this->x + s.x) % q); + } + } + + dihedral_t act(dihedral_t r) const { + if (this->is_reflection) { + return dihedral_t(!(r.is_reflection), ((q + this->x) - r.x) % q); + } else { + return dihedral_t(r.is_reflection, (this->x + r.x) % q); + } + } + + potts_t act_inverse(potts_t s) const { + if (this->is_reflection) { + return this->act(s); + } else { + return potts_t(((s.x + q) - this->x) % q); + } + } + + dihedral_t act_inverse(dihedral_t r) const { + if (this->is_reflection) { + return this->act(r); + } else { + return dihedral_t(r.is_reflection, ((r.x + q) - this->x) % q); + } + } +}; + diff --git a/examples/src/models/potts/potts.hpp b/examples/src/models/potts/potts.hpp new file mode 100644 index 0000000..f4765e2 --- /dev/null +++ b/examples/src/models/potts/potts.hpp @@ -0,0 +1,72 @@ +#pragma once + +#include +#include + +#include +#include "../On/vector.hpp" + +template +class potts_t { + public: + q_t x; + + typedef vector_t M_t; + typedef vector_t F_t; + + potts_t() : x(0) {} + potts_t(q_t x) : x(x) {} + + inline vector_t operator*(v_t a) const { + vector_t result; + result.fill(0); + result[x] = (int)a; + + return result; + } + + inline vector_t operator*(double a) const { + vector_t result; + result.fill(0.0); + result[x] = a; + + return result; + } + + inline vector_t operator-(const potts_t &s) const { + vector_t result; + result.fill(0); + + result[x]++; + result[s.x]--; + + return result; + } +}; + +// we could inherit norm_squared from vector.h, but convention dictates that +// potts norms be changed by a constant factor +template +double norm_squared(vector_t s) { + double total = 0; + for (double& x : s) { + total += pow(x, 2); + } + + return total * (double)q / ((double)q - 1.0); +} + +// we could inherit write_magnetization from vector.h, but since M.x must sum +// to nv we don't need to write the last element +template +void write_magnetization(vector_t M, FILE *outfile) { + for (int& x : M) { + fwrite(&x, sizeof(int), q - 1, outfile); + } +} + +// knock yourself out +const potts_t states[256] = {{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {21}, {22}, {23}, {24}, {25}, {26}, {27}, {28}, {29}, {30}, {31}, {32}, {33}, {34}, {35}, {36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46}, {47}, {48}, {49}, {50}, {51}, {52}, {53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}, {61}, {62}, {63}, {64}, {65}, {66}, {67}, {68}, {69}, {70}, {71}, {72}, {73}, {74}, {75}, {76}, {77}, {78}, {79}, {80}, {81}, {82}, {83}, {84}, {85}, {86}, {87}, {88}, {89}, {90}, {91}, {92}, {93}, {94}, {95}, {96}, {97}, {98}, {99}, {100}, {101}, {102}, {103}, {104}, {105}, {106}, {107}, {108}, {109}, {110}, {111}, {112}, {113}, {114}, {115}, {116}, {117}, {118}, {119}, {120}, {121}, {122}, {123}, {124}, {125}, {126}, {127}, {128}, {129}, {130}, {131}, {132}, {133}, {134}, {135}, {136}, {137}, {138}, {139}, {140}, {141}, {142}, {143}, {144}, {145}, {146}, {147}, {148}, {149}, {150}, {151}, {152}, {153}, {154}, {155}, {156}, {157}, {158}, {159}, {160}, {161}, {162}, {163}, {164}, {165}, {166}, {167}, {168}, {169}, {170}, {171}, {172}, {173}, {174}, {175}, {176}, {177}, {178}, {179}, {180}, {181}, {182}, {183}, {184}, {185}, {186}, {187}, {188}, {189}, {190}, {191}, {192}, {193}, {194}, {195}, {196}, {197}, {198}, {199}, {200}, {201}, {202}, {203}, {204}, {205}, {206}, {207}, {208}, {209}, {210}, {211}, {212}, {213}, {214}, {215}, {216}, {217}, {218}, {219}, {220}, {221}, {222}, {223}, {224}, {225}, {226}, {227}, {228}, {229}, {230}, {231}, {232}, {233}, {234}, {235}, {236}, {237}, {238}, {239}, {240}, {241}, {242}, {243}, {244}, {245}, {246}, {247}, {248}, {249}, {250}, {251}, {252}, {253}, {254}, {255}}; +template +q_t state_to_ind(potts_t state) { return (q_t)state.x; } + diff --git a/examples/src/models/potts/symmetric.hpp b/examples/src/models/potts/symmetric.hpp new file mode 100644 index 0000000..8636f15 --- /dev/null +++ b/examples/src/models/potts/symmetric.hpp @@ -0,0 +1,51 @@ + +#pragma once + +#include +#include +#include +#include "potts.hpp" + +template +class symmetric_t : public std::array { + public: + + symmetric_t() { + for (q_t i = 0; i < q; i++) { + (*this)[i] = i; + } + } + + potts_t act(const potts_t &s) const { + return potts_t((*this)[s.x]); + } + + symmetric_t act(const symmetric_t& r) const { + symmetric_t r_rot; + for (q_t i = 0; i < q; i++) { + r_rot[i] = (*this)[r[i]]; + } + + return r_rot; + } + + potts_t act_inverse(const potts_t& s) const { + for (q_t i = 0; i < q; i++) { + if ((*this)[i] == s.x) { + return potts_t(i); + } + } + + exit(EXIT_FAILURE); + } + + symmetric_t act_inverse(const symmetric_t& r) const { + symmetric_t r_rot; + for (q_t i = 0; i < q; i++) { + r_rot[(*this)[i]] = r[i]; + } + + return r_rot; + } +}; + diff --git a/examples/src/models/potts/wolff_clock.cpp b/examples/src/models/potts/wolff_clock.cpp new file mode 100644 index 0000000..020415d --- /dev/null +++ b/examples/src/models/potts/wolff_clock.cpp @@ -0,0 +1,155 @@ + +#include + +#ifdef HAVE_GLUT +#include +#endif + +// include your group and spin space +#include "dihedral.hpp" +#include "potts.hpp" +#include + +// hack to speed things up considerably +#define N_STATES POTTSQ +#include + +#include + +// include wolff.hpp +#include + +typedef state_t , potts_t> sim_t; + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e4; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); + + bool silent = false; + bool draw = false; + unsigned int window_size = 512; + + int opt; + q_t H_ind = 0; + + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:")) != -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. nth call couples to state n + H_vec[H_ind] = atof(optarg); + H_ind++; + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + case 'd': +#ifdef HAVE_GLUT + draw = true; + break; +#else + printf("You didn't compile this with the glut library installed!\n"); + exit(EXIT_FAILURE); +#endif + case 'w': + window_size = atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + // initialize random number generator + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + // define spin-spin coupling + std::function &, const potts_t&)> Z = [] (const potts_t& s1, const potts_t& s2) -> double { + return cos(2 * M_PI * (double)(s1.x + POTTSQ - s2.x) / (double)POTTSQ); + }; + + // define spin-field coupling + std::function &)> B = [=] (const potts_t& s) -> double { + return H_vec[s.x]; + }; + + // initialize state object + state_t , potts_t> s(D, L, T, Z, B); + + // define function that generates self-inverse rotations + 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); + q_t x = dist(r); + rot.x = (2 * v.x + x + 1) % POTTSQ; + + return rot; + }; + + // define function that updates any number of measurements + std::function 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; + }; + } else { + // a more complex example: measure the average magnetization, and draw the spin configuration to the screen + +#ifdef HAVE_GLUT + // initialize glut + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("wolff"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0.0, L, 0.0, L); + + measurement = [&] (const sim_t& s) { + average_M += (double)s.M[0] / (double)N / (double)s.nv; + 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]); + glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ)); + glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); + } + glFlush(); + }; +#endif + } + + // run wolff for N cluster flips + wolff(N, s, gen_R, measurement, rng, 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); + + // 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 new file mode 100644 index 0000000..a1e9284 --- /dev/null +++ b/examples/src/models/potts/wolff_potts.cpp @@ -0,0 +1,213 @@ + +#include +#include + +#ifdef HAVE_GLUT +#include +#endif + +// include your group and spin space +#include "symmetric.hpp" +#include "potts.hpp" + +// hack to speed things up considerably +#define N_STATES POTTSQ +#include + +// include wolff.h +#include +#include +#include +#include + +typedef state_t , potts_t> sim_t; + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e4; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); + + bool silent = false; + bool draw = false; + bool N_is_sweeps = false; + unsigned int window_size = 512; + + // don't measure anything by default + unsigned char measurement_flags = 0; + + int opt; + q_t H_ind = 0; + + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -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. nth call couples to state n + H_vec[H_ind] = atof(optarg); + H_ind++; + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + case 'S': + N_is_sweeps = true; + break; + case 'd': +#ifdef HAVE_GLUT + draw = true; + break; +#else + printf("You didn't compile this with the glut library installed!\n"); + exit(EXIT_FAILURE); +#endif + case 'w': + window_size = atoi(optarg); + break; + case 'M': + measurement_flags ^= 1 << atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + // get nanosecond timestamp for unique run id + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + // initialize random number generator + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + // define spin-spin coupling + std::function &, const potts_t&)> Z = [] (const potts_t& s1, const potts_t& s2) -> double { + if (s1.x == s2.x) { + return 1.0; + } else { + return 0.0; + } + }; + + // define spin-field coupling + std::function &)> B = [=] (const potts_t& s) -> double { + return H_vec[s.x]; + }; + + // initialize state object + state_t , potts_t> s(D, L, T, Z, B); + + // define function that generates self-inverse rotations + 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); + q_t j = dist(r); + q_t swap_v; + if (j < v.x) { + swap_v = j; + } else { + swap_v = j + 1; + } + + rot[v.x] = swap_v; + rot[swap_v] = v.x; + + return rot; + }; + + FILE **outfiles = measure_setup_files(measurement_flags, timestamp); + + std::function 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; + }; + } else if (draw) { +#ifdef HAVE_GLUT + // initialize glut + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("wolff"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0.0, L, 0.0, L); + + other_f = [] (const sim_t& s) { + 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]); + glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ)); + glRecti(i / s.L, i % s.L, (i / s.L) + 1, (i % s.L) + 1); + } + glFlush(); + }; +#endif + } else { + other_f = [] (const sim_t& s) {}; + } + + std::function measurements = measure_function_write_files(measurement_flags, outfiles, other_f); + + // add line to metadata file with run info + { + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); + + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"POTTS\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, POTTSQ, s.D, s.L, s.nv, s.ne, T); + + for (q_t i = 0; i < POTTSQ; i++) { + fprintf(outfile_info, "%.15f", H_vec[i]); + if (i < POTTSQ - 1) { + fprintf(outfile_info, ", "); + } + } + + fprintf(outfile_info, "} |>\n"); + + fclose(outfile_info); + } + + // run wolff for N cluster flips + 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); + 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); + } else { + wolff(N, s, gen_R, measurements, rng, silent); + } + + // free the random number generator + free(H_vec); + measure_free_files(measurement_flags, outfiles); + + return 0; + +} + diff --git a/examples/src/models/roughening/CMakeLists.txt b/examples/src/models/roughening/CMakeLists.txt new file mode 100644 index 0000000..51a8644 --- /dev/null +++ b/examples/src/models/roughening/CMakeLists.txt @@ -0,0 +1,21 @@ + +add_executable(wolff_dgm wolff_dgm.cpp) +add_executable(wolff_cgm wolff_cgm.cpp) + +find_library(GL NAMES GL) +find_library(GLU NAMES GLU) +find_library(GLUT NAMES glut) + +if (${GLUT} MATCHES "GLUT-NOTFOUND") + target_link_libraries(wolff_dgm wolff wolff_examples) + target_link_libraries(wolff_cgm wolff wolff_examples) +else() + target_compile_definitions(wolff_dgm PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_cgm PUBLIC HAVE_GLUT) + + target_link_libraries(wolff_dgm wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_cgm wolff wolff_examples glut GL GLU) +endif() + +install(TARGETS wolff_dgm wolff_cgm DESTINATION ${CMAKE_INSTALL_BINDIR}) + diff --git a/examples/src/models/roughening/dihedral_inf.hpp b/examples/src/models/roughening/dihedral_inf.hpp new file mode 100644 index 0000000..19fa195 --- /dev/null +++ b/examples/src/models/roughening/dihedral_inf.hpp @@ -0,0 +1,47 @@ + +#include +#include +#include "height.hpp" + +template +class dihedral_inf_t { + public: + bool is_reflection; + T x; + + dihedral_inf_t() : is_reflection(false), x(0) {} + dihedral_inf_t(bool x, T y) : is_reflection(x), x(y) {} + + height_t act(const height_t& h) const { + if (this->is_reflection) { + return height_t(this->x - h.x); + } else { + return height_t(this->x + h.x); + } + } + + dihedral_inf_t act(const dihedral_inf_t& r) const { + if (this->is_reflection) { + return dihedral_inf_t(!r.is_reflection, this->x - r.x); + } else { + return dihedral_inf_t(r.is_reflection, this->x + r.x); + } + } + + height_t act_inverse(const height_t& h) const { + if (this->is_reflection) { + return this->act(h); + } else { + return height_t(h.x - this->x); + } + } + + dihedral_inf_t act_inverse(const dihedral_inf_t& r) const { + if (this->is_reflection) { + return this->act(r); + } else { + return dihedral_inf_t(r.is_reflection, r.x - this->x); + } + } +}; + diff --git a/examples/src/models/roughening/height.hpp b/examples/src/models/roughening/height.hpp new file mode 100644 index 0000000..4023063 --- /dev/null +++ b/examples/src/models/roughening/height.hpp @@ -0,0 +1,75 @@ + +#pragma once + +#include +#include + +#include + +/* The following is the minimum definition of a spin class. + * + * The class must contain an M_t and an F_t for holding the sum of an + * integer number of spins and a double-weighted number of spins, + * respectively. + * + * void init(X_t *p); + * void free_spin(X_t p); + * void free_spin(M_t p); + * void free_spin(F_t p); + * X_t copy(X_t x); + * void add(M_t *x1, int factor, X_t x2); + * void add(F_t *x1, double factor, X_t x2); + * M_t scalar_multiple(int factor, X_t x); + * F_t scalar_multiple(double factor, X_t x); + * double norm_squared(F_t x); + * void write_magnetization(M_t M, FILE *outfile); + * + */ + +template +struct height_t { + T x; + + typedef T M_t; + typedef double F_t; + + height_t() : x(0) {} + + height_t(T x) : x(x) {} + + inline T operator*(v_t a) const { + return x * a; + } + + inline double operator*(double a) const { + return x * a; + } + + inline T operator-(const height_t& h) const { + return x - h.x; + } +}; + +template +inline T& operator+=(T& M, const height_t &h) { + M += h.x; + + return M; +} + +template +inline T& operator-=(T& M, const height_t &h) { + M -= h.x; + + return M; +} + +double norm_squared(double h) { + return pow(h, 2); +} + +template +void write_magnetization(T M, FILE *outfile) { + fwrite(&M, sizeof(T), 1, outfile); +} + diff --git a/examples/src/models/roughening/wolff_cgm.cpp b/examples/src/models/roughening/wolff_cgm.cpp new file mode 100644 index 0000000..65f8d66 --- /dev/null +++ b/examples/src/models/roughening/wolff_cgm.cpp @@ -0,0 +1,167 @@ + +#include + +#ifdef HAVE_GLUT +#include +#endif + +// include your group and spin space +#include "dihedral_inf.hpp" +#include "height.hpp" + +#include + +// include wolff.h +#include + +typedef state_t , height_t> sim_t; + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e4; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double H = 0; + + bool silent = false; + bool draw = false; + unsigned int window_size = 512; + double epsilon = 1; + + int opt; + + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:e:")) != -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. nth call couples to state n + H = atof(optarg); + break; + case 'e': // external field. nth call couples to state n + epsilon = atof(optarg); + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + case 'd': +#ifdef HAVE_GLUT + draw = true; + break; +#else + printf("You didn't compile this with the glut library installed!\n"); + exit(EXIT_FAILURE); +#endif + case 'w': + window_size = atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + // initialize random number generator + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + // define spin-spin coupling + std::function &, const height_t&)> Z = [] (const height_t& h1, const height_t& h2) -> double { + return -pow(h1.x - h2.x, 2); + }; + + // define spin-field coupling + std::function )> B = [=] (height_t h) -> double { + return -H * pow(h.x, 2);; + }; + + // initialize state object + sim_t s(D, L, T, Z, B); + + // define function that generates self-inverse rotations + std::function (std::mt19937&, height_t)> gen_R = [=] (std::mt19937& r, height_t h) -> dihedral_inf_t { + dihedral_inf_t rot; + rot.is_reflection = true; + std::normal_distribution dist(0.0,1.0); + + double amount = epsilon * dist(r); + + rot.x = 2 * h.x + amount; + + return rot; + }; + + // define function that updates any number of measurements + std::function 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 / (double)N / (double)s.nv; + }; + } else { + // a more complex example: measure the average magnetization, and draw the spin configuration to the screen + +#ifdef HAVE_GLUT + // initialize glut + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("wolff"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0.0, L, 0.0, L); + + measurement = [&] (const sim_t& s) { + average_M += (double)s.M / (double)N / (double)s.nv; + glClear(GL_COLOR_BUFFER_BIT); + double max_h = INT64_MIN; + double min_h = INT64_MAX; + for (v_t i = 0; i < pow(L, 2); i++) { + double cur_h = (s.R.act_inverse(s.spins[i])).x; + if (cur_h < min_h) { + min_h = cur_h; + } + if (cur_h > max_h) { + max_h = cur_h; + } + } + + for (v_t i = 0; i < pow(L, 2); i++) { + double cur_h = (s.R.act_inverse(s.spins[i])).x; + double mag = ((double)(cur_h - min_h)) / ((double)(max_h - min_h)); + glColor3f(mag, mag, mag); + glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); + } + glFlush(); + }; +#endif + } + + // run wolff for N cluster flips + wolff(N, s, gen_R, measurement, rng, silent); + + // tell us what we found! + printf("%" PRIcount " DGM runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, = %g\n", N, D, L, T, H, average_M); + + // free the random number generator + + if (draw) { + } + + return 0; + +} + diff --git a/examples/src/models/roughening/wolff_dgm.cpp b/examples/src/models/roughening/wolff_dgm.cpp new file mode 100644 index 0000000..8395382 --- /dev/null +++ b/examples/src/models/roughening/wolff_dgm.cpp @@ -0,0 +1,164 @@ + +#include + +#ifdef HAVE_GLUT +#include +#endif + +// include your group and spin space +#include "dihedral_inf.hpp" +#include "height.hpp" + +#include + +// include wolff.h +#include + +typedef state_t , height_t> sim_t; + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e4; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double H = 0; + + bool silent = false; + bool draw = false; + unsigned int window_size = 512; + uint64_t epsilon = 1; + + int opt; + + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:e:")) != -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. nth call couples to state n + H = atof(optarg); + break; + case 'e': // external field. nth call couples to state n + epsilon = atof(optarg); + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + case 'd': +#ifdef HAVE_GLUT + draw = true; + break; +#else + printf("You didn't compile this with the glut library installed!\n"); + exit(EXIT_FAILURE); +#endif + case 'w': + window_size = atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + // initialize random number generator + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + // define spin-spin coupling + std::function &, const height_t&)> Z = [] (const height_t& h1, const height_t& h2) -> double { + return -pow(h1.x - h2.x, 2); + }; + + // define spin-field coupling + std::function &)> B = [=] (const height_t& h) -> double { + return -H * pow(h.x, 2);; + }; + + // initialize state object + sim_t s(D, L, T, Z, B); + + // define function that generates self-inverse rotations + std::function (std::mt19937&, height_t)> gen_R = [=] (std::mt19937& r, height_t h) -> dihedral_inf_t { + dihedral_inf_t rot; + rot.is_reflection = true; + + std::uniform_int_distribution dist(-epsilon,epsilon); + + rot.x = 2 * h.x + dist(r); + + return rot; + }; + + // define function that updates any number of measurements + std::function 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 / (double)N / (double)s.nv; + }; + } else { + // a more complex example: measure the average magnetization, and draw the spin configuration to the screen + +#ifdef HAVE_GLUT + // initialize glut + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("wolff"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0.0, L, 0.0, L); + + measurement = [&] (const sim_t& s) { + average_M += (double)s.M / (double)N / (double)s.nv; + glClear(GL_COLOR_BUFFER_BIT); + int64_t max_h = INT64_MIN; + int64_t min_h = INT64_MAX; + for (v_t i = 0; i < pow(L, 2); i++) { + int64_t cur_h = (s.R.act_inverse(s.spins[i])).x; + if (cur_h < min_h) { + min_h = cur_h; + } + if (cur_h > max_h) { + max_h = cur_h; + } + } + + for (v_t i = 0; i < pow(L, 2); i++) { + int64_t cur_h = (s.R.act_inverse(s.spins[i])).x; + double mag = ((double)(cur_h - min_h)) / ((double)(max_h - min_h)); + glColor3f(mag, mag, mag); + glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); + } + glFlush(); + }; +#endif + } + + // run wolff for N cluster flips + wolff(N, s, gen_R, measurement, rng, silent); + + // tell us what we found! + printf("%" PRIcount " DGM runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, = %g\n", N, D, L, T, H, average_M); + + if (draw) { + } + + return 0; + +} + diff --git a/examples/src/tools/CMakeLists.txt b/examples/src/tools/CMakeLists.txt new file mode 100644 index 0000000..1c73c2d --- /dev/null +++ b/examples/src/tools/CMakeLists.txt @@ -0,0 +1,7 @@ + +find_library(fftw REQUIRED NAMES fftw3) + +add_executable(analyze_correlations analyze_correlations.cpp) + +target_link_libraries(analyze_correlations fftw3 wolff) + diff --git a/examples/src/tools/analyze_correlations.cpp b/examples/src/tools/analyze_correlations.cpp new file mode 100644 index 0000000..abeaff3 --- /dev/null +++ b/examples/src/tools/analyze_correlations.cpp @@ -0,0 +1,486 @@ + +#include +#include +#include +#include +#include +#include +#include + +template +double mean(int N, T *data) { + double total = 0; + for (int i = 0; i < N; i++) { + total += (double)data[i]; + } + + return total / N; +} + +double squared_mean(int N, double *data) { + double total = 0; + for (int i = 0; i < N; i++) { + total += pow(data[i], 2); + } + + return total / N; +} + +double central_moment(int N, double *data, double mean, int m) { + double total = 0; + for (int i = 0; i < N; i++) { + total += pow(data[i] - mean, m); + } + + return total / N; +} + +void compute_OO(int N, fftw_plan forward_plan, double *forward_data, fftw_plan reverse_plan, double *reverse_data) { + + fftw_execute(forward_plan); + + reverse_data[0] = forward_data[0] * forward_data[0]; + reverse_data[N / 2] = forward_data[N/2] * forward_data[N/2]; + + for (count_t i = 1; i < N / 2; i++) { + reverse_data[i] = pow(forward_data[i], 2) + pow(forward_data[N - i], 2); + reverse_data[N - i] = 0; + } + + fftw_execute(reverse_plan); + +} + +double finite_energy(q_t nb, double *J, q_t q, double *H, v_t nv, v_t ne, uint32_t *bo, uint32_t *so) { + double energy = 0; + + v_t tot = 0; + for (q_t i = 0; i < nb - 1; i++) { + energy -= J[i] * bo[i]; + tot += bo[i]; + } + + energy -= J[nb - 1] * (ne - tot); + + tot = 0; + for (q_t i = 0; i < q - 1; i++) { + energy -= H[i] * so[i]; + tot += so[i]; + } + + energy -= H[q - 1] * (nv - tot); + + return energy; +} + +int main (int argc, char *argv[]) { + count_t drop = (count_t)1e4; + count_t length = (count_t)1e4; + bool speedy_drop = false; + bool from_stdin = false; + bool oldstyle = false; + + int opt; + + while ((opt = getopt(argc, argv, "d:l:spo")) != -1) { + switch (opt) { + case 'd': + drop = (count_t)atof(optarg); + break; + case 'l': + length = (count_t)atof(optarg); + break; + case 's': + speedy_drop = true; + break; + case 'p': + from_stdin = true; + break; + case 'o': + oldstyle = true; + break; + default: + exit(EXIT_FAILURE); + } + } + FILE *metadata; + + fftw_set_timelimit(1); + + if (from_stdin) { + metadata = stdin; + } else { + metadata = fopen("wolff_metadata.txt", "r"); + } + + if (metadata == NULL) { + printf("Metadata file not found. Make sure you are in the correct directory!\n"); + exit(EXIT_FAILURE); + } + + unsigned long id; + char *model = (char *)malloc(32 * sizeof(char)); + + if (model == NULL) { + printf("Malloc failed.\n"); + exit(EXIT_FAILURE); + } + + q_t q; + D_t D; + L_t L; + v_t nv, ne; + + while (EOF != fscanf(metadata, "<| \"ID\" -> %lu, \"MODEL\" -> \"%[^\"]\", \"q\" -> %" SCNq ", \"D\" -> %" SCND ", \"L\" -> %" SCNL ", \"NV\" -> %" SCNv ", \"NE\" -> %" SCNv ", ", &id, model, &q, &D, &L, &nv, &ne)) { + + printf("%lu: Processing...\n", id); + +// bool is_finite = 0 == strcmp(model, "ISING") || 0 == strcmp(model, "POTTS") || 0 == strcmp(model, "CLOCK"); + + if (oldstyle) { + q_t nb; + double T; + fscanf(metadata, "\"NB\" -> %" SCNq ", \"T\" -> %lf, \"J\" -> {", &nb, &T); + double *J = (double *)malloc(nb * sizeof(double)); + double *H = (double *)malloc(q * sizeof(double)); + + if (J == NULL || H == NULL) { + printf("%lu: Malloc failed.\n", id); + break; + } + + for (q_t i = 0; i < nb - 1; i++) { + fscanf(metadata, "%lf, ", &(J[i])); + } + fscanf(metadata, "%lf}, \"H\" -> {", &(J[nb - 1])); + for (q_t i = 0; i < q - 1; i++) { + fscanf(metadata, "%lf, ", &(H[i])); + } + fscanf(metadata, "%lf} |>\n", &(H[q - 1])); + + char *filename_M = (char *)malloc(128 * sizeof(char)); + char *filename_B = (char *)malloc(128 * sizeof(char)); + char *filename_S = (char *)malloc(128 * sizeof(char)); + + if (filename_M == NULL || filename_B == NULL || filename_S == NULL) { + printf("%lu: Malloc failed.\n", id); + break; + } + + sprintf(filename_M, "wolff_%lu_M.dat", id); + sprintf(filename_B, "wolff_%lu_B.dat", id); + sprintf(filename_S, "wolff_%lu_S.dat", id); + + FILE *file_M = fopen(filename_M, "rb"); + FILE *file_B = fopen(filename_B, "rb"); + FILE *file_S = fopen(filename_S, "rb"); + + if (file_M == NULL || file_B == NULL || file_S == NULL) { + printf("%lu: Opening data file failed.\n", id); + break; + } + + fseek(file_S, 0, SEEK_END); + unsigned long N = ftell(file_S) / sizeof(uint32_t); + fseek(file_S, 0, SEEK_SET); + + if (speedy_drop) { + drop = N - pow(2, floor(log(N) / log(2))); + } else { + if (N % 2 == 1 && drop % 2 == 0) { + drop++; // make sure M is even + } + } + + if (N <= drop) { + printf("\033[F%lu: Number of steps %lu is less than %" PRIcount ", nothing done.\n", id, N, drop); + } else { + int M = N - drop; + + double M_f = (double)M; + + if (length > M) { + length = M; + } + + double *forward_data = (double *)fftw_malloc(M * sizeof(double)); + fftw_plan forward_plan = fftw_plan_r2r_1d(M, forward_data, forward_data, FFTW_R2HC, 0); + double *reverse_data = (double *)fftw_malloc(M * sizeof(double)); + fftw_plan reverse_plan = fftw_plan_r2r_1d(M, reverse_data, reverse_data, FFTW_HC2R, 0); + + + uint32_t *data_S = (uint32_t *)malloc(N * sizeof(uint32_t)); + fread(data_S, N, sizeof(uint32_t), file_S); + for (count_t i = 0; i < M; i++) { + forward_data[i] = (double)data_S[drop + i]; + } + free(data_S); + double mean_S = mean(M, forward_data); + double squaredMean_S = squared_mean(M, forward_data); + double moment2_S = central_moment(M, forward_data, mean_S, 2); + double moment4_S = central_moment(M, forward_data, mean_S, 4); + + compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); + + sprintf(filename_S, "wolff_%lu_S_OO.dat", id); + + FILE *file_S = fopen(filename_S, "wb"); + fwrite(&M_f, sizeof(double), 1, file_S); + fwrite(&mean_S, sizeof(double), 1, file_S); + fwrite(&squaredMean_S, sizeof(double), 1, file_S); + fwrite(&moment2_S, sizeof(double), 1, file_S); + fwrite(&moment4_S, sizeof(double), 1, file_S); + fwrite(reverse_data, sizeof(double), length, file_S); + fclose(file_S); + + uint32_t *data_B = (uint32_t *)malloc((nb - 1) * N * sizeof(uint32_t)); + uint32_t *data_M = (uint32_t *)malloc((q - 1) * N * sizeof(uint32_t)); + fread(data_B, N * (nb - 1), sizeof(uint32_t), file_B); + fread(data_M, N * (q - 1), sizeof(uint32_t), file_M); + + for (count_t i = 0; i < M; i++) { + forward_data[i] = finite_energy(nb, J, q, H, nv, ne, data_B + (nb - 1) * (drop + i), data_M + (q - 1) * (drop + i)); + } + + double mean_E = mean(M, forward_data); + double squaredMean_E = squared_mean(M, forward_data); + double moment2_E = central_moment(M, forward_data, mean_E, 2); + double moment4_E = central_moment(M, forward_data, mean_E, 4); + + free(data_B); + free(data_M); + + compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); + + sprintf(filename_B, "wolff_%lu_E_OO.dat", id); + + FILE *file_E = fopen(filename_B, "wb"); + fwrite(&M_f, sizeof(double), 1, file_E); + fwrite(&mean_E, sizeof(double), 1, file_E); + fwrite(&squaredMean_E, sizeof(double), 1, file_E); + fwrite(&moment2_E, sizeof(double), 1, file_E); + fwrite(&moment4_E, sizeof(double), 1, file_E); + fwrite(reverse_data, sizeof(double), length, file_E); + fclose(file_E); + + printf("\033[F%lu: Correlation functions for %d steps written.\n", id, M); + + fftw_destroy_plan(forward_plan); + fftw_destroy_plan(reverse_plan); + fftw_free(forward_data); + fftw_free(reverse_data); + + } + + fclose(file_M); + fclose(file_B); + fclose(file_S); + + free(J); + free(H); + + free(filename_S); + free(filename_B); + free(filename_M); + + } else { + char *junk = (char *)malloc(1024 * sizeof(char)); + fscanf(metadata, "%[^\n]\n", junk); // throw away the rest of the line, we don't need it + free(junk); + + char *filename_E = (char *)malloc(128 * sizeof(char)); + char *filename_F = (char *)malloc(128 * sizeof(char)); + char *filename_M = (char *)malloc(128 * sizeof(char)); + char *filename_S = (char *)malloc(128 * sizeof(char)); + + sprintf(filename_E, "wolff_%lu_E.dat", id); + sprintf(filename_F, "wolff_%lu_F.dat", id); + sprintf(filename_M, "wolff_%lu_M.dat", id); + sprintf(filename_S, "wolff_%lu_S.dat", id); + + FILE *file_E = fopen(filename_E, "rb"); + FILE *file_F = fopen(filename_F, "rb"); + FILE *file_M = fopen(filename_M, "rb"); + FILE *file_S = fopen(filename_S, "rb"); + + fseek(file_S, 0, SEEK_END); + unsigned long N = ftell(file_S) / sizeof(uint32_t); + fseek(file_S, 0, SEEK_SET); + + if (speedy_drop) { + drop = N - pow(2, floor(log(N) / log(2))); + } else { + if (N % 2 == 1 && drop % 2 == 0) { + drop++; // make sure M is even + } + } + + if (N <= drop) { + printf("\033[F%lu: Number of steps %lu is less than %" PRIcount ", nothing done.\n", id, N, drop); + } else { + int M = N - drop; + double M_f = (double)M; + + if (length > M) { + length = M; + } + + double *forward_data = (double *)fftw_malloc(M * sizeof(double)); + fftw_plan forward_plan = fftw_plan_r2r_1d(M, forward_data, forward_data, FFTW_R2HC, 0); + + double *reverse_data = (double *)fftw_malloc(M * sizeof(double)); + fftw_plan reverse_plan = fftw_plan_r2r_1d(M, reverse_data, reverse_data, FFTW_HC2R, 0); + + if (file_S != NULL) { + uint32_t *data_S = (uint32_t *)malloc(N * sizeof(uint32_t)); + + fread(data_S, sizeof(uint32_t), N, file_S); + fclose(file_S); + + for (int i = 0; i < M; i++) { + forward_data[i] = (double)data_S[drop + i]; + } + free(data_S); + + double mean_S = mean(M, forward_data); + double squaredMean_S = squared_mean(M, forward_data); + double moment2_S = central_moment(M, forward_data, mean_S, 2); + double moment4_S = central_moment(M, forward_data, mean_S, 4); + + compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); + + sprintf(filename_S, "wolff_%lu_S_OO.dat", id); + FILE *file_S_new = fopen(filename_S, "wb"); + fwrite(&M_f, sizeof(double), 1, file_S_new); + fwrite(&mean_S, sizeof(double), 1, file_S_new); + fwrite(&squaredMean_S, sizeof(double), 1, file_S_new); + fwrite(&moment2_S, sizeof(double), 1, file_S_new); + fwrite(&moment4_S, sizeof(double), 1, file_S_new); + fwrite(reverse_data, sizeof(double), length, file_S_new); + fclose(file_S_new); + } + if (file_F != NULL) { + float *data_F = (float *)malloc(N * sizeof(float)); + + fread(data_F, sizeof(float), N, file_F); + fclose(file_F); + + for (int i = 0; i < M; i++) { + forward_data[i] = (double)data_F[drop + i]; + } + free(data_F); + + double mean_F = mean(M, forward_data); + double squaredMean_F = squared_mean(M, forward_data); + double moment2_F = central_moment(M, forward_data, mean_F, 2); + double moment4_F = central_moment(M, forward_data, mean_F, 4); + + compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); + + sprintf(filename_F, "wolff_%lu_F_OO.dat", id); + FILE *file_F_new = fopen(filename_F, "wb"); + fwrite(&M_f, sizeof(double), 1, file_F_new); + fwrite(&mean_F, sizeof(double), 1, file_F_new); + fwrite(&squaredMean_F, sizeof(double), 1, file_F_new); + fwrite(&moment2_F, sizeof(double), 1, file_F_new); + fwrite(&moment4_F, sizeof(double), 1, file_F_new); + fwrite(reverse_data, sizeof(double), length, file_F_new); + fclose(file_F_new); + } + if (file_E != NULL) { + float *data_E = (float *)malloc(N * sizeof(float)); + + fread(data_E, sizeof(float), N, file_E); + fclose(file_E); + + for (int i = 0; i < M; i++) { + forward_data[i] = (double)data_E[drop + i]; + } + free(data_E); + + double mean_E = mean(M, forward_data); + double squaredMean_E = squared_mean(M, forward_data); + double moment2_E = central_moment(M, forward_data, mean_E, 2); + double moment4_E = central_moment(M, forward_data, mean_E, 4); + + compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); + + sprintf(filename_E, "wolff_%lu_E_OO.dat", id); + FILE *file_E_new = fopen(filename_E, "wb"); + fwrite(&M_f, sizeof(double), 1, file_E_new); + fwrite(&mean_E, sizeof(double), 1, file_E_new); + fwrite(&squaredMean_E, sizeof(double), 1, file_E_new); + fwrite(&moment2_E, sizeof(double), 1, file_E_new); + fwrite(&moment4_E, sizeof(double), 1, file_E_new); + fwrite(reverse_data, sizeof(double), length, file_E_new); + fclose(file_E_new); + } + if (file_M != NULL) { + if (0 == strcmp(model, "PLANAR")) { + float *data_M = (float *)malloc(2 * N * sizeof(float)); + fread(data_M, sizeof(float), 2 * N, file_M); + fclose(file_M); + for (int i = 0; i < M; i++) { + forward_data[i] = (double)sqrt(pow(data_M[2 * drop + 2 * i], 2) + pow(data_M[2 * drop + 2 * i + 1], 2)); + } + free(data_M); + } else if (0 == strcmp(model, "HEISENBERG")) { + float *data_M = (float *)malloc(3 * N * sizeof(float)); + fread(data_M, sizeof(float), 3 * N, file_M); + fclose(file_M); + for (int i = 0; i < M; i++) { + forward_data[i] = sqrt(pow(data_M[3 * drop + 3 * i], 2) + pow(data_M[3 * drop + 3 * i + 1], 2) + pow(data_M[3 * drop + 3 * i + 2], 2)); + } + free(data_M); + } else if (0 == strcmp(model, "ISING")) { + int *data_M = (int *)malloc(N * sizeof(float)); + fread(data_M, sizeof(int), N, file_M); + fclose(file_M); + for (int i = 0; i < M; i++) { + forward_data[i] = (double)data_M[i]; + } + free(data_M); + } else { + printf("UNKNOWN MODEL\n"); + exit(EXIT_FAILURE); + } + + double mean_M = mean(M, forward_data); + double squaredMean_M = squared_mean(M, forward_data); + double moment2_M = central_moment(M, forward_data, mean_M, 2); + double moment4_M = central_moment(M, forward_data, mean_M, 4); + + compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); + + sprintf(filename_M, "wolff_%lu_M_OO.dat", id); + FILE *file_M_new = fopen(filename_M, "wb"); + fwrite(&M_f, sizeof(double), 1, file_M_new); + fwrite(&mean_M, sizeof(double), 1, file_M_new); + fwrite(&squaredMean_M, sizeof(double), 1, file_M_new); + fwrite(&moment2_M, sizeof(double), 1, file_M_new); + fwrite(&moment4_M, sizeof(double), 1, file_M_new); + fwrite(reverse_data, sizeof(double), length, file_M_new); + fclose(file_M_new); + } + + printf("\033[F%lu: Correlation functions for %d steps written.\n", id, M); + fftw_destroy_plan(forward_plan); + fftw_destroy_plan(reverse_plan); + fftw_free(forward_data); + fftw_free(reverse_data); + + } + free(filename_E); + free(filename_S); + free(filename_F); + free(filename_M); + } + } + + free(model); + fclose(metadata); + fftw_cleanup(); + + return 0; +} + diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt new file mode 100644 index 0000000..a31500e --- /dev/null +++ b/lib/CMakeLists.txt @@ -0,0 +1,17 @@ + +project(libwolff LANGUAGES C CXX) + +add_library(wolff SHARED src/graph.cpp) + +target_include_directories(wolff PUBLIC $ + $) + +install(TARGETS wolff EXPORT wolffConfig + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) +install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) + +install(EXPORT wolffConfig DESTINATION share/wolff/cmake) + +export(TARGETS wolff FILE wolffConfig.cmake) diff --git a/lib/angle.h b/lib/angle.h deleted file mode 100644 index c3f128e..0000000 --- a/lib/angle.h +++ /dev/null @@ -1,48 +0,0 @@ -#pragma once - -#include "types.h" - -#include -#include "vector.h" - -class angle_t { - public: - double x; - - typedef vector_t<2, double> M_t; - typedef vector_t<2, double> F_t; - - angle_t() : x(0) {} - angle_t(double x) : x(x) {} - - inline vector_t<2, double> operator*(v_t a) const { - vector_t<2, double>M; - M[0] = a * cos(x); - M[1] = a * sin(x); - - return M; - } - - inline vector_t<2, double> operator*(double a) const { - vector_t<2, double>M; - M[0] = a * cos(x); - M[1] = a * sin(x); - - return M; - } -}; - -inline vector_t<2,double>& operator+=(vector_t<2,double>& M, const angle_t& theta) { - M[0] += cos(theta.x); - M[1] += sin(theta.x); - - return M; -} - -inline vector_t<2,double>& operator-=(vector_t<2,double>& M, const angle_t& theta) { - M[0] -= cos(theta.x); - M[1] -= sin(theta.x); - - return M; -} - diff --git a/lib/circle_group.h b/lib/circle_group.h deleted file mode 100644 index cb7cadd..0000000 --- a/lib/circle_group.h +++ /dev/null @@ -1,46 +0,0 @@ -#pragma once - -#include "angle.h" - -class circle_group_t { - public: - bool is_reflection; - double x; - - circle_group_t() : is_reflection(false), x(0) {} - circle_group_t(bool x, double y) : is_reflection(x), x(y) {} - - angle_t act(const angle_t& theta) const { - if (is_reflection) { - return angle_t(fmod(2 * M_PI + x - theta.x, 2 * M_PI)); - } else { - return angle_t(fmod(x + theta.x, 2 * M_PI)); - } - } - - circle_group_t act(const circle_group_t& r) const { - if (is_reflection) { - return circle_group_t(!r.is_reflection, fmod(2 * M_PI + x - r.x, 2 * M_PI)); - } else { - return circle_group_t(r.is_reflection, fmod(x + r.x, 2 * M_PI)); - } - } - - angle_t act_inverse(const angle_t& theta) const { - if (is_reflection) { - return act(theta); - } else { - return angle_t(fmod(2 * M_PI + theta.x - x, 2 * M_PI)); - } - } - - circle_group_t act_inverse(const circle_group_t& r) const { - if (is_reflection) { - return act(r); - } else { - return circle_group_t(r.is_reflection, fmod(2 * M_PI + r.x - x, 2 * M_PI)); - } - } -}; - - diff --git a/lib/cluster.h b/lib/cluster.h deleted file mode 100644 index f57bb68..0000000 --- a/lib/cluster.h +++ /dev/null @@ -1,111 +0,0 @@ - -#pragma once - -#include -#include -#include -#include -#include - -#include "types.h" -#include "state.h" -#include "graph.h" - -template -void flip_cluster(state_t& s, v_t v0, const R_t& r, gsl_rng *rand) { - v_t nv = 0; - - std::stack stack; - stack.push(v0); - - std::vector marks(s.g.nv, false); - - while (!stack.empty()) { - v_t v = stack.top(); - stack.pop(); - - if (!marks[v]) { // don't reprocess anyone we've already visited! - marks[v] = true; - - X_t si_new; -#ifndef NOFIELD - R_t R_new; - - bool v_is_ghost = (v == s.nv); // ghost site has the last index - - if (v_is_ghost) { - R_new = r.act(s.R); // if we are, then we're moving the transformation - } else -#endif - { - si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site - } - - for (const v_t &vn : s.g.v_adj[v]) { - double dE, prob; - -#ifndef NOFIELD - bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost - - if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond - X_t rs_old, rs_new; - v_t non_ghost; - - if (vn_is_ghost) { - // if our neighbor is the ghost, the current site is a normal - // spin - rotate it back! - rs_old = s.R.act_inverse(s.spins[v]); - rs_new = s.R.act_inverse(si_new); - non_ghost = v; - } else { - /* if we're the ghost, we need to rotate our neighbor back in - both the old and new ways */ - rs_old = s.R.act_inverse(s.spins[vn]); - rs_new = R_new.act_inverse(s.spins[vn]); - non_ghost = vn; - } - - dE = s.H(rs_old) - s.H(rs_new); - -#ifdef FINITE_STATES - 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); - } else // this is a perfectly normal bond! -#endif - { - dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]); - -#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); - -#ifndef FINITE_STATES - prob = 1.0 - exp(-dE / s.T); -#endif - - if (gsl_rng_uniform(rand) < prob) { - stack.push(vn); // push the neighboring vertex to the stack - } - } - -#ifndef NOFIELD - if (v_is_ghost) { - s.R = R_new; - } else -#endif - { - s.spins[v] = si_new; - nv++; - } - } - } - - s.last_cluster_size = nv; -} - diff --git a/lib/colors.h b/lib/colors.h deleted file mode 100644 index 04d39a8..0000000 --- a/lib/colors.h +++ /dev/null @@ -1,34 +0,0 @@ -#pragma once - -#include "types.h" - -double hue_to_R(double theta) { - if (((M_PI / 3 <= theta) && (theta < 2 * M_PI / 3)) || ((4 * M_PI / 3 <= theta) && (theta < 5 * M_PI / 3))) { - return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); - } else if (((0 <= theta) && (theta < M_PI / 3)) || ((5 * M_PI / 3 <= theta) && (theta <= 2 * M_PI))) { - return 1.0; - } else { - return 0.0; - } -} - -double hue_to_G(double theta) { - if (((0 <= theta) && (theta < M_PI / 3)) || ((M_PI <= theta) && (theta < 4 * M_PI / 3))) { - return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); - } else if (((M_PI / 3 <= theta) && (theta < 2 * M_PI / 3)) || ((2 * M_PI / 3 <= theta) && (theta < M_PI))) { - return 1.0; - } else { - return 0.0; - } -} - -double hue_to_B(double theta) { - if (((2 * M_PI / 3 <= theta) && (theta < M_PI)) || ((5 * M_PI / 3 <= theta) && (theta <= 2 * M_PI))) { - return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); - } else if (((M_PI <= theta) && (theta < 4 * M_PI / 3)) || ((4 * M_PI / 3 <= theta) && (theta < 5 * M_PI / 3))) { - return 1.0; - } else { - return 0.0; - } -} - diff --git a/lib/convex.c b/lib/convex.c deleted file mode 100644 index 7255ad4..0000000 --- a/lib/convex.c +++ /dev/null @@ -1,102 +0,0 @@ - -#include "convex.h" - -double slope(point_t *P, point_t *Q) { - return (Q->y - P->y) / ((double)(Q->x) - (double)(P->x)); -} - -double *get_convex_minorant(count_t n, double *Gammas) { - if (n < 2) { - return Gammas; - } - - list_t *L = (list_t *)calloc(1, sizeof(list_t)); - L->p = (point_t *)calloc(1, sizeof(point_t)); - L->p->x = 0; - L->p->y = Gammas[0]; - - list_t *pos = L; - - for (count_t i = 1; i < n; i++) { - pos->next = (list_t *)calloc(1, sizeof(list_t)); - pos->next->p = (point_t *)calloc(1, sizeof(point_t)); - pos->next->p->x = i; - pos->next->p->y = Gammas[i]; - pos->next->prev = pos; - pos = pos->next; - } - - pos->next = (list_t *)calloc(1, sizeof(list_t)); - pos->next->p = (point_t *)calloc(1, sizeof(point_t)); - pos->next->p->x = n; - pos->next->p->y = 0; - pos->next->prev = pos; - - list_t *X = L; - list_t *Y = L->next; - list_t *Z = Y->next; - - while (true) { - if (slope(X->p, Y->p) <= slope(Y->p, Z->p)) { - X = Y; - Y = Z; - if (Z->next == NULL) { - break; - } else { - Z = Z->next; - } - } else { - Y->prev->next = Y->next; - Y->next->prev = Y->prev; - free(Y->p); - free(Y); - if (X->prev != NULL) { - Y = X; - X = X->prev; - } else { - if (Z->next != NULL) { - Y = Z; - Z = Z->next; - } else { - break; - } - } - } - } - - pos = L; - - double *g = (double *)calloc(n + 1, sizeof(double)); - double rho = 0; - - for (count_t i = 0; i < n + 1; i++) { - if (i > pos->next->p->x) { - pos = pos->next; - } - - g[i] = pos->p->y + ((double)i - (double)(pos->p->x)) * (pos->next->p->y - pos->p->y) / ((double)(pos->next->p->x) - (double)(pos->p->x)); - - if (i rho) { - rho = Gammas[i] - g[i]; - } - } else { - if (0 - g[i] > rho) { - rho = 0 - g[i]; - } - } - } - - for (count_t i = 0; i < n + 1; i++) { - g[i] += rho / 2; - } - - while (L != NULL) { - free(L->p); - L = L->next; - free(L); - } - - return g; -} - diff --git a/lib/convex.h b/lib/convex.h deleted file mode 100644 index 5a405d4..0000000 --- a/lib/convex.h +++ /dev/null @@ -1,23 +0,0 @@ - -#pragma once - -#include -#include -#include -#include - -#include "types.h" - -typedef struct { - count_t x; - double y; -} point_t; - -typedef struct list_tag { - struct list_tag *prev; - struct list_tag *next; - point_t *p; -} list_t; - -double *get_convex_minorant(count_t n, double *Gammas); - diff --git a/lib/correlation.h b/lib/correlation.h deleted file mode 100644 index 29357a5..0000000 --- a/lib/correlation.h +++ /dev/null @@ -1,23 +0,0 @@ - -#pragma once - -#include "types.h" -#include "state.h" - -#include - -template -double correlation_length(const state_t & s) { - double total = 0; - -#ifdef DIMENSION - for (D_t j = 0; j < DIMENSION; j++) { -#else - for (D_t j = 0; j < s.D; j++) { -#endif - total += norm_squared(s.ReF[j]) + norm_squared(s.ImF[j]); - } - - return total / s.D; -} - diff --git a/lib/dihedral.h b/lib/dihedral.h deleted file mode 100644 index 8d0472b..0000000 --- a/lib/dihedral.h +++ /dev/null @@ -1,48 +0,0 @@ - -#pragma once - -#include "types.h" -#include "potts.h" - -template -class dihedral_t { - public: - bool is_reflection; - q_t x; - - dihedral_t() : is_reflection(false), x(0) {} - dihedral_t(bool x, q_t y) : is_reflection(x), x(y) {} - - potts_t act(const potts_t& s) const { - if (this->is_reflection) { - return potts_t(((q + this->x) - s.x) % q); - } else { - return potts_t((this->x + s.x) % q); - } - } - - dihedral_t act(dihedral_t r) const { - if (this->is_reflection) { - return dihedral_t(!(r.is_reflection), ((q + this->x) - r.x) % q); - } else { - return dihedral_t(r.is_reflection, (this->x + r.x) % q); - } - } - - potts_t act_inverse(potts_t s) const { - if (this->is_reflection) { - return this->act(s); - } else { - return potts_t(((s.x + q) - this->x) % q); - } - } - - dihedral_t act_inverse(dihedral_t r) const { - if (this->is_reflection) { - return this->act(r); - } else { - return dihedral_t(r.is_reflection, ((r.x + q) - this->x) % q); - } - } -}; - diff --git a/lib/dihedral_inf.h b/lib/dihedral_inf.h deleted file mode 100644 index a064b48..0000000 --- a/lib/dihedral_inf.h +++ /dev/null @@ -1,47 +0,0 @@ - -#include "types.h" -#include -#include "height.h" - -template -class dihedral_inf_t { - public: - bool is_reflection; - T x; - - dihedral_inf_t() : is_reflection(false), x(0) {} - dihedral_inf_t(bool x, T y) : is_reflection(x), x(y) {} - - height_t act(const height_t& h) const { - if (this->is_reflection) { - return height_t(this->x - h.x); - } else { - return height_t(this->x + h.x); - } - } - - dihedral_inf_t act(const dihedral_inf_t& r) const { - if (this->is_reflection) { - return dihedral_inf_t(!r.is_reflection, this->x - r.x); - } else { - return dihedral_inf_t(r.is_reflection, this->x + r.x); - } - } - - height_t act_inverse(const height_t& h) const { - if (this->is_reflection) { - return this->act(h); - } else { - return height_t(h.x - this->x); - } - } - - dihedral_inf_t act_inverse(const dihedral_inf_t& r) const { - if (this->is_reflection) { - return this->act(r); - } else { - return dihedral_inf_t(r.is_reflection, r.x - this->x); - } - } -}; - diff --git a/lib/finite_states.h b/lib/finite_states.h deleted file mode 100644 index 426edad..0000000 --- a/lib/finite_states.h +++ /dev/null @@ -1,40 +0,0 @@ - -#pragma once - -#include -#include -#include - -#define FINITE_STATES - -// must have N_STATES, states[N_STATES], and state_to_ind defined before -// invoking header - -std::array, N_STATES>, N_STATES> J_probs; -#ifndef NOFIELD -std::array, N_STATES> H_probs; -#endif - -template -#ifndef NOFIELD -void initialize_probs(std::function J, std::function H, double T) { - for (q_t i = 0; i < N_STATES; i++) { - for (q_t j = 0; j < N_STATES; j++) { - H_probs[i][j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T); - } - } -#else -void initialize_probs(std::function J, double T) { -#endif - for (q_t i = 0; i < N_STATES; i++) { - for (q_t j = 0; j < N_STATES; j++) { - for (q_t k = 0; k < N_STATES; k++) { - J_probs[i][j][k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T); - } - } - } -} - - - - diff --git a/lib/graph.cpp b/lib/graph.cpp deleted file mode 100644 index ca251f3..0000000 --- a/lib/graph.cpp +++ /dev/null @@ -1,44 +0,0 @@ - -#include "graph.h" - -graph_t::graph_t(D_t D, L_t L) { - nv = 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(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; - - 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))); - } - } -} - -void graph_t::add_ext() { - for (std::vector& v_adj_i : v_adj) { - v_adj_i.push_back(nv); - } - - v_adj.resize(nv + 1); - coordinate.resize(nv + 1); - v_adj[nv].reserve(nv); - - for (v_t i = 0; i < nv; i++) { - v_adj[nv].push_back(i); - } - - coordinate[nv].resize(coordinate[0].size()); - - ne += nv; - nv += 1; -} - diff --git a/lib/graph.h b/lib/graph.h deleted file mode 100644 index de06924..0000000 --- a/lib/graph.h +++ /dev/null @@ -1,22 +0,0 @@ - -#pragma once - -#include -#include -#include -#include - -#include "types.h" - -class graph_t { - public: - v_t ne; - v_t nv; - std::vector> v_adj; - std::vector> coordinate; - - graph_t(D_t D, L_t L); - void add_ext(); - -}; - diff --git a/lib/height.h b/lib/height.h deleted file mode 100644 index d2173fe..0000000 --- a/lib/height.h +++ /dev/null @@ -1,75 +0,0 @@ - -#pragma once - -#include -#include - -#include "types.h" - -/* The following is the minimum definition of a spin class. - * - * The class must contain an M_t and an F_t for holding the sum of an - * integer number of spins and a double-weighted number of spins, - * respectively. - * - * void init(X_t *p); - * void free_spin(X_t p); - * void free_spin(M_t p); - * void free_spin(F_t p); - * X_t copy(X_t x); - * void add(M_t *x1, int factor, X_t x2); - * void add(F_t *x1, double factor, X_t x2); - * M_t scalar_multiple(int factor, X_t x); - * F_t scalar_multiple(double factor, X_t x); - * double norm_squared(F_t x); - * void write_magnetization(M_t M, FILE *outfile); - * - */ - -template -struct height_t { - T x; - - typedef T M_t; - typedef double F_t; - - height_t() : x(0) {} - - height_t(T x) : x(x) {} - - inline T operator*(v_t a) const { - return x * a; - } - - inline double operator*(double a) const { - return x * a; - } - - inline T operator-(const height_t& h) const { - return x - h.x; - } -}; - -template -inline T& operator+=(T& M, const height_t &h) { - M += h.x; - - return M; -} - -template -inline T& operator-=(T& M, const height_t &h) { - M -= h.x; - - return M; -} - -double norm_squared(double h) { - return pow(h, 2); -} - -template -void write_magnetization(T M, FILE *outfile) { - fwrite(&M, sizeof(T), 1, outfile); -} - diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp new file mode 100644 index 0000000..c10a211 --- /dev/null +++ b/lib/include/wolff.hpp @@ -0,0 +1,35 @@ + +#include "wolff/cluster.hpp" +#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) { + +#ifdef FINITE_STATES +#ifdef NOFIELD + initialize_probs(s.J, s.T); +#else + initialize_probs(s.J, s.H, s.T); +#endif +#endif + + std::uniform_int_distribution 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); + } + + if (!silent) { + printf("\033[F\033[J"); + } + 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 new file mode 100644 index 0000000..104f3c2 --- /dev/null +++ b/lib/include/wolff/cluster.hpp @@ -0,0 +1,111 @@ + +#pragma once + +#include +#include +#include +#include + +#include "types.h" +#include "state.hpp" +#include "graph.hpp" + +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; + + std::stack stack; + stack.push(v0); + + std::vector marks(s.g.nv, false); + + while (!stack.empty()) { + v_t v = stack.top(); + stack.pop(); + + if (!marks[v]) { // don't reprocess anyone we've already visited! + marks[v] = true; + + X_t si_new; +#ifndef NOFIELD + R_t R_new; + + bool v_is_ghost = (v == s.nv); // ghost site has the last index + + if (v_is_ghost) { + R_new = r.act(s.R); // if we are, then we're moving the transformation + } else +#endif + { + si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site + } + + for (const v_t &vn : s.g.v_adj[v]) { + double dE, prob; + +#ifndef NOFIELD + bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost + + if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond + X_t rs_old, rs_new; + v_t non_ghost; + + if (vn_is_ghost) { + // if our neighbor is the ghost, the current site is a normal + // spin - rotate it back! + rs_old = s.R.act_inverse(s.spins[v]); + rs_new = s.R.act_inverse(si_new); + non_ghost = v; + } else { + /* if we're the ghost, we need to rotate our neighbor back in + both the old and new ways */ + rs_old = s.R.act_inverse(s.spins[vn]); + rs_new = R_new.act_inverse(s.spins[vn]); + non_ghost = vn; + } + + dE = s.H(rs_old) - s.H(rs_new); + +#ifdef FINITE_STATES + 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); + } else // this is a perfectly normal bond! +#endif + { + dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]); + +#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); + +#ifndef FINITE_STATES + prob = 1.0 - exp(-dE / s.T); +#endif + + if (dist(rand) < prob) { + stack.push(vn); // push the neighboring vertex to the stack + } + } + +#ifndef NOFIELD + if (v_is_ghost) { + s.R = R_new; + } else +#endif + { + s.spins[v] = si_new; + nv++; + } + } + } + + s.last_cluster_size = nv; +} + diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp new file mode 100644 index 0000000..426edad --- /dev/null +++ b/lib/include/wolff/finite_states.hpp @@ -0,0 +1,40 @@ + +#pragma once + +#include +#include +#include + +#define FINITE_STATES + +// must have N_STATES, states[N_STATES], and state_to_ind defined before +// invoking header + +std::array, N_STATES>, N_STATES> J_probs; +#ifndef NOFIELD +std::array, N_STATES> H_probs; +#endif + +template +#ifndef NOFIELD +void initialize_probs(std::function J, std::function H, double T) { + for (q_t i = 0; i < N_STATES; i++) { + for (q_t j = 0; j < N_STATES; j++) { + H_probs[i][j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T); + } + } +#else +void initialize_probs(std::function J, double T) { +#endif + for (q_t i = 0; i < N_STATES; i++) { + for (q_t j = 0; j < N_STATES; j++) { + for (q_t k = 0; k < N_STATES; k++) { + J_probs[i][j][k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T); + } + } + } +} + + + + diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp new file mode 100644 index 0000000..165544a --- /dev/null +++ b/lib/include/wolff/graph.hpp @@ -0,0 +1,21 @@ + +#pragma once + +#include +#include +#include +#include + +#include "types.h" + +class graph_t { + public: + v_t ne; + v_t nv; + std::vector> v_adj; + std::vector> coordinate; + + graph_t(D_t D, L_t L); + void add_ext(); +}; + diff --git a/lib/include/wolff/state.hpp b/lib/include/wolff/state.hpp new file mode 100644 index 0000000..1f5e359 --- /dev/null +++ b/lib/include/wolff/state.hpp @@ -0,0 +1,89 @@ + +#pragma once + +#include +#include + +#include "types.h" +#include "graph.hpp" + +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; + v_t ne; + graph_t g; + double T; + std::vector spins; + R_t R; + double E; + typename X_t::M_t M; // the "sum" of the spins, like the total magnetization + v_t last_cluster_size; + std::vector ReF; + std::vector ImF; + + std::function J; +#ifndef NOFIELD + std::function H; + + state_t(D_t D, L_t L, double T, std::function J, std::function H) : D(D), L(L), g(D, L), T(T), R(), J(J), H(H) { +#else + state_t(D_t D, L_t L, double T, std::function J) : D(D), L(L), g(D, L), T(T), R(), J(J) { +#endif + nv = g.nv; + ne = g.ne; + spins.resize(nv); +#ifndef NOFIELD + g.add_ext(); + E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]); +#else + E = - (double)ne * J(spins[0], spins[0]); +#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/include/wolff/types.h b/lib/include/wolff/types.h new file mode 100644 index 0000000..ec9efaf --- /dev/null +++ b/lib/include/wolff/types.h @@ -0,0 +1,36 @@ + +#include + +typedef uint_fast32_t v_t; +typedef uint_fast8_t q_t; +typedef uint_fast16_t R_t; +typedef uint_fast8_t D_t; +typedef uint_fast16_t L_t; +typedef uint_fast64_t count_t; +typedef int_fast64_t h_t; + +#define MAX_v UINT_FAST32_MAX +#define MAX_Q UINT_FAST8_MAX +#define MAX_R UINT_FAST16_MAX +#define MAX_D UINT_FAST8_MAX +#define MAX_L UINT_FAST16_MAX +#define MAX_COUNT UINT_FAST64_MAX +#define MAX_H INT_FAST64_MAX +#define MIN_H INT_FAST64_MIN + +#define PRIv PRIuFAST32 +#define PRIq PRIuFAST8 +#define PRIR PRIuFAST16 +#define PRID PRIuFAST8 +#define PRIL PRIuFAST16 +#define PRIcount PRIuFAST64 +#define PRIh PRIdFAST64 + +#define SCNv SCNuFAST32 +#define SCNq SCNuFAST8 +#define SCNR SCNuFAST16 +#define SCND SCNuFAST8 +#define SCNL SCNuFAST16 +#define SCNcount SCNuFAST64 +#define SCNh SCNdFAST64 + diff --git a/lib/ising.h b/lib/ising.h deleted file mode 100644 index 45058fb..0000000 --- a/lib/ising.h +++ /dev/null @@ -1,59 +0,0 @@ -#pragma once - -#include -#include - -#include "types.h" - -class ising_t { - public: - bool x; - - typedef int M_t; - typedef double F_t; - - ising_t() : x(false) {} - ising_t(bool x) : x(x) {} - ising_t(int x) : x((bool)x) {} - - inline int operator*(v_t a) const { - if (x) { - return -(int)a; - } else { - return (int)a; - } - } - - inline double operator*(double a) const { - if (x) { - return -a; - } else { - return a; - } - } - - inline int operator-(const ising_t &s) const { - if (x == s.x) { - return 0; - } else { - if (x) { - return -2; - } else { - return 2; - } - } - } -}; - -double norm_squared(double s) { - return pow(s, 2); -} - -void write_magnetization(int M, FILE *outfile) { - fwrite(&M, sizeof(int), 1, outfile); -} - -#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/lib/measure.h b/lib/measure.h deleted file mode 100644 index 2c5ffb7..0000000 --- a/lib/measure.h +++ /dev/null @@ -1,63 +0,0 @@ - -#pragma once - -#include "state.h" -#include "correlation.h" -#include - -#define POSSIBLE_MEASUREMENTS 4 -const unsigned char measurement_energy = 1 << 0; -const unsigned char measurement_clusterSize = 1 << 1; -const unsigned char measurement_magnetization = 1 << 2; -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 *)); - - 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; -} - -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]); - } - if (flags & measurement_magnetization) { - write_magnetization(s.M, files[2]); - } - if (flags & measurement_fourierZero) { - float smaller_X = (float)correlation_length(s); - fwrite(&smaller_X, sizeof(float), 1, files[3]); - } - - other_f(s); - }; -} - -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); -} - - diff --git a/lib/orthogonal.h b/lib/orthogonal.h deleted file mode 100644 index e1bf33a..0000000 --- a/lib/orthogonal.h +++ /dev/null @@ -1,200 +0,0 @@ - -#pragma once - -#include -#include -#include -#include - -#include "state.h" -#include "types.h" -#include "vector.h" - -template -class orthogonal_t : public std::array, q> { - public : - bool is_reflection; - - orthogonal_t() : is_reflection(false) { - for (q_t i = 0; i < q; i++) { - (*this)[i].fill(0); - (*this)[i][i] = (T)1; - } - } - - vector_t act(const vector_t & v) const { - vector_t v_rot; - v_rot.fill(0); - - if (is_reflection) { - double prod = 0; - for (q_t i = 0; i < q; i++) { - prod += v[i] * (*this)[0][i]; - } - for (q_t i = 0; i < q; i++) { - v_rot[i] = v[i] - 2 * prod * (*this)[0][i]; - } - } else { - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - v_rot[i] += (*this)[i][j] * v[j]; - } - } - } - - return v_rot; - } - - orthogonal_t act(const orthogonal_t & m) const { - orthogonal_t m_rot; - - m_rot.is_reflection = false; - - if (is_reflection) { - for (q_t i = 0; i < q; i++) { - double akOki = 0; - - for (q_t k = 0; k < q; k++) { - akOki += (*this)[0][k] * m[k][i]; - } - - for (q_t j = 0; j < q; j++) { - m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][j]; - } - } - } else { - for (q_t i = 0; i < q; i++) { - m_rot[i].fill(0); - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m_rot[i][j] += (*this)[i][j] * m[j][k]; - } - } - } - } - - return m_rot; - } - - vector_t act_inverse(const vector_t & v) const { - if (is_reflection) { - return this->act(v); // reflections are their own inverse - } else { - vector_t v_rot; - v_rot.fill(0); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - v_rot[i] += (*this)[j][i] * v[j]; - } - } - - return v_rot; - } - } - - vector_t act_inverse(const orthogonal_t & m) const { - if (is_reflection) { - return this->act(m); // reflections are their own inverse - } else { - orthogonal_t m_rot; - m_rot.is_reflection = false; - - for (q_t i = 0; i < q; i++) { - m_rot[i].fill(0); - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m_rot[i][j] += (*this)[j][i] * m[j][k]; - } - } - } - - return m_rot; - } - } - -}; - - -template -orthogonal_t generate_rotation_uniform (gsl_rng *r, const vector_t & v) { - orthogonal_t ptr; - ptr.is_reflection = true; - - double v2 = 0; - - for (q_t i = 0; i < q; i++) { - ptr[0][i] = gsl_ran_ugaussian(r); - v2 += ptr[0][i] * ptr[0][i]; - } - - double mag_v = sqrt(v2); - - for (q_t i = 0; i < q; i++) { - ptr[0][i] /= mag_v; - } - - return ptr; -} - -template -orthogonal_t generate_rotation_perturbation (gsl_rng *r, const vector_t & v0, double epsilon, unsigned int n) { - orthogonal_t m; - m.is_reflection = true; - - vector_t v; - - if (n > 1) { - unsigned int rotation = gsl_rng_uniform_int(r, n); - - double cosr = cos(2 * M_PI * rotation / (double)n / 2.0); - double sinr = sin(2 * M_PI * rotation / (double)n / 2.0); - - v[0] = v0[0] * cosr - v0[1] * sinr; - v[1] = v0[1] * cosr + v0[0] * sinr; - - for (q_t i = 2; i < q; i++) { - v[i] = v0[i]; - } - } else { - v = v0; - } - - double m_dot_v = 0; - - for (q_t i = 0; i < q; i++) { - m[0][i] = gsl_ran_ugaussian(r); // create a random vector - m_dot_v += m[0][i] * v[i]; - } - - double v2 = 0; - - for (q_t i = 0; i < q; i++) { - m[0][i] = m[0][i] - m_dot_v * v[i]; // find the component orthogonal to v - v2 += pow(m[0][i], 2); - } - - double mag_v = sqrt(v2); - - for (q_t i = 0; i < q; i++) { - m[0][i] /= mag_v; // normalize - } - - v2 = 0; - - double factor = gsl_ran_gaussian(r, epsilon); - - for (q_t i = 0; i < q; i++) { - m[0][i] += factor * v[i]; // perturb orthogonal vector in original direction - v2 += pow(m[0][i], 2); - } - - mag_v = sqrt(v2); - - for (q_t i = 0; i < q; i++) { - m[0][i] /= mag_v; // normalize - } - - return m; -} - diff --git a/lib/potts.h b/lib/potts.h deleted file mode 100644 index 1764e53..0000000 --- a/lib/potts.h +++ /dev/null @@ -1,72 +0,0 @@ -#pragma once - -#include -#include - -#include "types.h" -#include "vector.h" - -template -class potts_t { - public: - q_t x; - - typedef vector_t M_t; - typedef vector_t F_t; - - potts_t() : x(0) {} - potts_t(q_t x) : x(x) {} - - inline vector_t operator*(v_t a) const { - vector_t result; - result.fill(0); - result[x] = (int)a; - - return result; - } - - inline vector_t operator*(double a) const { - vector_t result; - result.fill(0.0); - result[x] = a; - - return result; - } - - inline vector_t operator-(const potts_t &s) const { - vector_t result; - result.fill(0); - - result[x]++; - result[s.x]--; - - return result; - } -}; - -// we could inherit norm_squared from vector.h, but convention dictates that -// potts norms be changed by a constant factor -template -double norm_squared(vector_t s) { - double total = 0; - for (double& x : s) { - total += pow(x, 2); - } - - return total * (double)q / ((double)q - 1.0); -} - -// we could inherit write_magnetization from vector.h, but since M.x must sum -// to nv we don't need to write the last element -template -void write_magnetization(vector_t M, FILE *outfile) { - for (int& x : M) { - fwrite(&x, sizeof(int), q - 1, outfile); - } -} - -// knock yourself out -const potts_t states[256] = {{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {21}, {22}, {23}, {24}, {25}, {26}, {27}, {28}, {29}, {30}, {31}, {32}, {33}, {34}, {35}, {36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46}, {47}, {48}, {49}, {50}, {51}, {52}, {53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}, {61}, {62}, {63}, {64}, {65}, {66}, {67}, {68}, {69}, {70}, {71}, {72}, {73}, {74}, {75}, {76}, {77}, {78}, {79}, {80}, {81}, {82}, {83}, {84}, {85}, {86}, {87}, {88}, {89}, {90}, {91}, {92}, {93}, {94}, {95}, {96}, {97}, {98}, {99}, {100}, {101}, {102}, {103}, {104}, {105}, {106}, {107}, {108}, {109}, {110}, {111}, {112}, {113}, {114}, {115}, {116}, {117}, {118}, {119}, {120}, {121}, {122}, {123}, {124}, {125}, {126}, {127}, {128}, {129}, {130}, {131}, {132}, {133}, {134}, {135}, {136}, {137}, {138}, {139}, {140}, {141}, {142}, {143}, {144}, {145}, {146}, {147}, {148}, {149}, {150}, {151}, {152}, {153}, {154}, {155}, {156}, {157}, {158}, {159}, {160}, {161}, {162}, {163}, {164}, {165}, {166}, {167}, {168}, {169}, {170}, {171}, {172}, {173}, {174}, {175}, {176}, {177}, {178}, {179}, {180}, {181}, {182}, {183}, {184}, {185}, {186}, {187}, {188}, {189}, {190}, {191}, {192}, {193}, {194}, {195}, {196}, {197}, {198}, {199}, {200}, {201}, {202}, {203}, {204}, {205}, {206}, {207}, {208}, {209}, {210}, {211}, {212}, {213}, {214}, {215}, {216}, {217}, {218}, {219}, {220}, {221}, {222}, {223}, {224}, {225}, {226}, {227}, {228}, {229}, {230}, {231}, {232}, {233}, {234}, {235}, {236}, {237}, {238}, {239}, {240}, {241}, {242}, {243}, {244}, {245}, {246}, {247}, {248}, {249}, {250}, {251}, {252}, {253}, {254}, {255}}; -template -q_t state_to_ind(potts_t state) { return (q_t)state.x; } - diff --git a/lib/rand.c b/lib/rand.c deleted file mode 100644 index 76f537d..0000000 --- a/lib/rand.c +++ /dev/null @@ -1,20 +0,0 @@ - -#include -#include - -unsigned long int rand_seed() { - /* - returns a random unsigned long integer read from the standard unix - pseudorandom device /dev/urandom - */ - - FILE *f = fopen("/dev/urandom", "r"); - assert(f != NULL); - - unsigned long int seed; - fread(&seed, sizeof(unsigned long int), 1, f); - - fclose(f); - - return seed; -} diff --git a/lib/rand.h b/lib/rand.h deleted file mode 100644 index 7bb5354..0000000 --- a/lib/rand.h +++ /dev/null @@ -1,16 +0,0 @@ - -#pragma once - -#include -#include - -#ifdef __cplusplus -extern "C" { -#endif - -unsigned long int rand_seed(); - -#ifdef __cplusplus -} -#endif - diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp new file mode 100644 index 0000000..4043413 --- /dev/null +++ b/lib/src/graph.cpp @@ -0,0 +1,44 @@ + +#include + +graph_t::graph_t(D_t D, L_t L) { + nv = 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(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; + + 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))); + } + } +} + +void graph_t::add_ext() { + for (std::vector& v_adj_i : v_adj) { + v_adj_i.push_back(nv); + } + + v_adj.resize(nv + 1); + coordinate.resize(nv + 1); + v_adj[nv].reserve(nv); + + for (v_t i = 0; i < nv; i++) { + v_adj[nv].push_back(i); + } + + coordinate[nv].resize(coordinate[0].size()); + + ne += nv; + nv += 1; +} + diff --git a/lib/state.h b/lib/state.h deleted file mode 100644 index cad453c..0000000 --- a/lib/state.h +++ /dev/null @@ -1,89 +0,0 @@ - -#pragma once - -#include -#include - -#include "types.h" -#include "graph.h" - -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; - v_t ne; - graph_t g; - double T; - std::vector spins; - R_t R; - double E; - typename X_t::M_t M; // the "sum" of the spins, like the total magnetization - v_t last_cluster_size; - std::vector ReF; - std::vector ImF; - - std::function J; -#ifndef NOFIELD - std::function H; - - state_t(D_t D, L_t L, double T, std::function J, std::function H) : D(D), L(L), g(D, L), T(T), R(), J(J), H(H) { -#else - state_t(D_t D, L_t L, double T, std::function J) : D(D), L(L), g(D, L), T(T), R(), J(J) { -#endif - nv = g.nv; - ne = g.ne; - spins.resize(nv); -#ifndef NOFIELD - g.add_ext(); - E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]); -#else - E = - (double)ne * J(spins[0], spins[0]); -#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/symmetric.h b/lib/symmetric.h deleted file mode 100644 index 9e9b6e4..0000000 --- a/lib/symmetric.h +++ /dev/null @@ -1,51 +0,0 @@ - -#pragma once - -#include -#include -#include "types.h" -#include "potts.h" - -template -class symmetric_t : public std::array { - public: - - symmetric_t() { - for (q_t i = 0; i < q; i++) { - (*this)[i] = i; - } - } - - potts_t act(const potts_t &s) const { - return potts_t((*this)[s.x]); - } - - symmetric_t act(const symmetric_t& r) const { - symmetric_t r_rot; - for (q_t i = 0; i < q; i++) { - r_rot[i] = (*this)[r[i]]; - } - - return r_rot; - } - - potts_t act_inverse(const potts_t& s) const { - for (q_t i = 0; i < q; i++) { - if ((*this)[i] == s.x) { - return potts_t(i); - } - } - - exit(EXIT_FAILURE); - } - - symmetric_t act_inverse(const symmetric_t& r) const { - symmetric_t r_rot; - for (q_t i = 0; i < q; i++) { - r_rot[(*this)[i]] = r[i]; - } - - return r_rot; - } -}; - diff --git a/lib/torus.h b/lib/torus.h deleted file mode 100644 index 2aead52..0000000 --- a/lib/torus.h +++ /dev/null @@ -1,64 +0,0 @@ - -#pragma once - -#include -#include -#include "types.h" - -template -class torus_t : public std::array { - public: - typedef std::array M_t; - typedef std::array F_t; - - torus_t() { - this->fill(0); - } - - inline torus_t operator*(v_t a) const { - torus_t x; - for (q_t i = 0; i < n; i++) { - x[i] = a * (*this)[i]; - } - - return x; - } - - inline torus_t operator*(double a) const { - torus_t x; - for (q_t i = 0; i < n; i++) { - x[i] = a * (*this)[i]; - } - - return x; - } - - inline torus_t& operator+=(const torus_t& x) { - for (q_t i = 0; i < n; i++) { - (*this)[i] += x[i]; - } - } - - inline torus_t& operator-=(const torus_t& x) { - for (q_t i = 0; i < n; i++) { - (*this)[i] -= x[i]; - } - } -}; - -template -double norm_squared(const torus_t& x) { - double tmp = 0; - for (const double& xi : x) { - tmp += pow(xi, 2); - } - return tmp; -} - -void write_magnetization(const torus_t& x, FILE *outfile) { - for (const double& xi : x) { - float tmp_xi = (float)xi; - fwrite(&tmp_xi, sizeof(float), 1, outfile); - } -} - diff --git a/lib/types.h b/lib/types.h deleted file mode 100644 index ec9efaf..0000000 --- a/lib/types.h +++ /dev/null @@ -1,36 +0,0 @@ - -#include - -typedef uint_fast32_t v_t; -typedef uint_fast8_t q_t; -typedef uint_fast16_t R_t; -typedef uint_fast8_t D_t; -typedef uint_fast16_t L_t; -typedef uint_fast64_t count_t; -typedef int_fast64_t h_t; - -#define MAX_v UINT_FAST32_MAX -#define MAX_Q UINT_FAST8_MAX -#define MAX_R UINT_FAST16_MAX -#define MAX_D UINT_FAST8_MAX -#define MAX_L UINT_FAST16_MAX -#define MAX_COUNT UINT_FAST64_MAX -#define MAX_H INT_FAST64_MAX -#define MIN_H INT_FAST64_MIN - -#define PRIv PRIuFAST32 -#define PRIq PRIuFAST8 -#define PRIR PRIuFAST16 -#define PRID PRIuFAST8 -#define PRIL PRIuFAST16 -#define PRIcount PRIuFAST64 -#define PRIh PRIdFAST64 - -#define SCNv SCNuFAST32 -#define SCNq SCNuFAST8 -#define SCNR SCNuFAST16 -#define SCND SCNuFAST8 -#define SCNL SCNuFAST16 -#define SCNcount SCNuFAST64 -#define SCNh SCNdFAST64 - diff --git a/lib/vector.h b/lib/vector.h deleted file mode 100644 index 7d0ee36..0000000 --- a/lib/vector.h +++ /dev/null @@ -1,118 +0,0 @@ - -#pragma once - -#include -#include -#include - -#include "types.h" - -template -class vector_t : public std::array { - public: - - // M_t needs to hold the sum of nv spins - typedef vector_t M_t; - - // F_t needs to hold the double-weighted sum of spins - typedef vector_t F_t; - - vector_t() { - this->fill((T)0); - (*this)[1] = (T)1; - } - - vector_t(const T *x) { - for (q_t i = 0; i < q; i++) { - (*this)[i] = x[i]; - } - } - - template - inline vector_t& operator+=(const vector_t &v) { - for (q_t i = 0; i < q; i++) { - (*this)[i] += (U)v[i]; - } - return *this; - } - - template - inline vector_t& operator-=(const vector_t &v) { - for (q_t i = 0; i < q; i++) { - (*this)[i] -= (U)v[i]; - } - return *this; - } - - inline vector_t operator*(v_t x) const { - vector_t result; - for (q_t i = 0; i < q; i++) { - result[i] = x * (*this)[i]; - } - - return result; - } - - inline vector_t operator*(double x) const { - vector_t result; - for (q_t i = 0; i < q; i++) { - result[i] = x * (*this)[i]; - } - - return result; - } - - inline vector_t operator-(const vector_t& v) const { - vector_t diff = *this; - diff -= v; - return diff; - } -}; - - -template -double norm_squared(vector_t v) { - double tmp = 0; - for (T &x : v) { - tmp += pow(x, 2); - } - - return tmp; -} - -template -void write_magnetization(vector_t M, FILE *outfile) { - for (q_t i = 0; i < q; i++) { - fwrite(&(M[i]), sizeof(T), q, outfile); - } -} - -// below functions and definitions are unnecessary for wolff.h but useful. - -template // save some space and don't write whole doubles -void write_magnetization(vector_t M, FILE *outfile) { - for (q_t i = 0; i < q; i++) { - float M_tmp = (float)M[i]; - fwrite(&M_tmp, sizeof(float), 1, outfile); - } -} - -template -T dot(const vector_t & v1, const vector_t & v2) { - T prod = 0; - - for (q_t i = 0; i < q; i++) { - prod += v1[i] * v2[i]; - } - - return prod; -} - -template -double H_vector(const vector_t & v1, T *H) { - vector_t H_vec(H); - return (double)(dot (v1, H_vec)); -} - -char const *ON_strings[] = {"TRIVIAL", "ISING", "PLANAR", "HEISENBERG"}; - diff --git a/lib/wolff.h b/lib/wolff.h deleted file mode 100644 index 141a5b2..0000000 --- a/lib/wolff.h +++ /dev/null @@ -1,33 +0,0 @@ - -#include "cluster.h" -#include "state.h" - -template -void wolff(count_t N, state_t & s, std::function gen_R, std::function &)> measurements, gsl_rng *r, bool silent) { - -#ifdef FINITE_STATES -#ifdef NOFIELD - initialize_probs(s.J, s.T); -#else - initialize_probs(s.J, s.H, s.T); -#endif -#endif - - 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 = gsl_rng_uniform_int(r, s.nv); - R_t step = gen_R(r, s.spins[v0]); - flip_cluster (s, v0, step, r); - - measurements(s); - } - - if (!silent) { - printf("\033[F\033[J"); - } - printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size); - -} - diff --git a/lib/z2.h b/lib/z2.h deleted file mode 100644 index a18d740..0000000 --- a/lib/z2.h +++ /dev/null @@ -1,53 +0,0 @@ - -#pragma once - -#include "types.h" -#include "ising.h" - -/* The minimum definition for a group type R_t to act on a spin type X_t is - * given by the following. - * - * void init(R_t *p); - * void free_spin(R_t r); - * R_t copy(R_t r); - * X_t act(R_t r, X_t x); - * R_t act(R_t r, R_t r); - * X_t act_inverse(R_t r, X_t x); - * R_t act_inverse(R_t r, R_t r); - * - */ - -class z2_t { - public: - bool x; - - z2_t() : x(false) {} - - z2_t(bool x) : x(x) {} - - ising_t act(const ising_t& s) const { - if (x) { - return ising_t(!s.x); - } else { - return ising_t(s.x); - } - } - - z2_t act(const z2_t& r) const { - if (x) { - return z2_t(!r.x); - } else { - return z2_t(r.x); - } - } - - ising_t act_inverse(const ising_t& s) const { - return this->act(s); - } - - z2_t act_inverse(const z2_t& r) const { - return this->act(r); - } -}; - - diff --git a/src/analyze_correlations.cpp b/src/analyze_correlations.cpp deleted file mode 100644 index 48ee426..0000000 --- a/src/analyze_correlations.cpp +++ /dev/null @@ -1,486 +0,0 @@ - -#include -#include -#include -#include -#include -#include -#include - -template -double mean(int N, T *data) { - double total = 0; - for (int i = 0; i < N; i++) { - total += (double)data[i]; - } - - return total / N; -} - -double squared_mean(int N, double *data) { - double total = 0; - for (int i = 0; i < N; i++) { - total += pow(data[i], 2); - } - - return total / N; -} - -double central_moment(int N, double *data, double mean, int m) { - double total = 0; - for (int i = 0; i < N; i++) { - total += pow(data[i] - mean, m); - } - - return total / N; -} - -void compute_OO(int N, fftw_plan forward_plan, double *forward_data, fftw_plan reverse_plan, double *reverse_data) { - - fftw_execute(forward_plan); - - reverse_data[0] = forward_data[0] * forward_data[0]; - reverse_data[N / 2] = forward_data[N/2] * forward_data[N/2]; - - for (count_t i = 1; i < N / 2; i++) { - reverse_data[i] = pow(forward_data[i], 2) + pow(forward_data[N - i], 2); - reverse_data[N - i] = 0; - } - - fftw_execute(reverse_plan); - -} - -double finite_energy(q_t nb, double *J, q_t q, double *H, v_t nv, v_t ne, uint32_t *bo, uint32_t *so) { - double energy = 0; - - v_t tot = 0; - for (q_t i = 0; i < nb - 1; i++) { - energy -= J[i] * bo[i]; - tot += bo[i]; - } - - energy -= J[nb - 1] * (ne - tot); - - tot = 0; - for (q_t i = 0; i < q - 1; i++) { - energy -= H[i] * so[i]; - tot += so[i]; - } - - energy -= H[q - 1] * (nv - tot); - - return energy; -} - -int main (int argc, char *argv[]) { - count_t drop = (count_t)1e4; - count_t length = (count_t)1e4; - bool speedy_drop = false; - bool from_stdin = false; - bool oldstyle = false; - - int opt; - - while ((opt = getopt(argc, argv, "d:l:spo")) != -1) { - switch (opt) { - case 'd': - drop = (count_t)atof(optarg); - break; - case 'l': - length = (count_t)atof(optarg); - break; - case 's': - speedy_drop = true; - break; - case 'p': - from_stdin = true; - break; - case 'o': - oldstyle = true; - break; - default: - exit(EXIT_FAILURE); - } - } - FILE *metadata; - - fftw_set_timelimit(1); - - if (from_stdin) { - metadata = stdin; - } else { - metadata = fopen("wolff_metadata.txt", "r"); - } - - if (metadata == NULL) { - printf("Metadata file not found. Make sure you are in the correct directory!\n"); - exit(EXIT_FAILURE); - } - - unsigned long id; - char *model = (char *)malloc(32 * sizeof(char)); - - if (model == NULL) { - printf("Malloc failed.\n"); - exit(EXIT_FAILURE); - } - - q_t q; - D_t D; - L_t L; - v_t nv, ne; - - while (EOF != fscanf(metadata, "<| \"ID\" -> %lu, \"MODEL\" -> \"%[^\"]\", \"q\" -> %" SCNq ", \"D\" -> %" SCND ", \"L\" -> %" SCNL ", \"NV\" -> %" SCNv ", \"NE\" -> %" SCNv ", ", &id, model, &q, &D, &L, &nv, &ne)) { - - printf("%lu: Processing...\n", id); - -// bool is_finite = 0 == strcmp(model, "ISING") || 0 == strcmp(model, "POTTS") || 0 == strcmp(model, "CLOCK"); - - if (oldstyle) { - q_t nb; - double T; - fscanf(metadata, "\"NB\" -> %" SCNq ", \"T\" -> %lf, \"J\" -> {", &nb, &T); - double *J = (double *)malloc(nb * sizeof(double)); - double *H = (double *)malloc(q * sizeof(double)); - - if (J == NULL || H == NULL) { - printf("%lu: Malloc failed.\n", id); - break; - } - - for (q_t i = 0; i < nb - 1; i++) { - fscanf(metadata, "%lf, ", &(J[i])); - } - fscanf(metadata, "%lf}, \"H\" -> {", &(J[nb - 1])); - for (q_t i = 0; i < q - 1; i++) { - fscanf(metadata, "%lf, ", &(H[i])); - } - fscanf(metadata, "%lf} |>\n", &(H[q - 1])); - - char *filename_M = (char *)malloc(128 * sizeof(char)); - char *filename_B = (char *)malloc(128 * sizeof(char)); - char *filename_S = (char *)malloc(128 * sizeof(char)); - - if (filename_M == NULL || filename_B == NULL || filename_S == NULL) { - printf("%lu: Malloc failed.\n", id); - break; - } - - sprintf(filename_M, "wolff_%lu_M.dat", id); - sprintf(filename_B, "wolff_%lu_B.dat", id); - sprintf(filename_S, "wolff_%lu_S.dat", id); - - FILE *file_M = fopen(filename_M, "rb"); - FILE *file_B = fopen(filename_B, "rb"); - FILE *file_S = fopen(filename_S, "rb"); - - if (file_M == NULL || file_B == NULL || file_S == NULL) { - printf("%lu: Opening data file failed.\n", id); - break; - } - - fseek(file_S, 0, SEEK_END); - unsigned long N = ftell(file_S) / sizeof(uint32_t); - fseek(file_S, 0, SEEK_SET); - - if (speedy_drop) { - drop = N - pow(2, floor(log(N) / log(2))); - } else { - if (N % 2 == 1 && drop % 2 == 0) { - drop++; // make sure M is even - } - } - - if (N <= drop) { - printf("\033[F%lu: Number of steps %lu is less than %" PRIcount ", nothing done.\n", id, N, drop); - } else { - int M = N - drop; - - double M_f = (double)M; - - if (length > M) { - length = M; - } - - double *forward_data = (double *)fftw_malloc(M * sizeof(double)); - fftw_plan forward_plan = fftw_plan_r2r_1d(M, forward_data, forward_data, FFTW_R2HC, 0); - double *reverse_data = (double *)fftw_malloc(M * sizeof(double)); - fftw_plan reverse_plan = fftw_plan_r2r_1d(M, reverse_data, reverse_data, FFTW_HC2R, 0); - - - uint32_t *data_S = (uint32_t *)malloc(N * sizeof(uint32_t)); - fread(data_S, N, sizeof(uint32_t), file_S); - for (count_t i = 0; i < M; i++) { - forward_data[i] = (double)data_S[drop + i]; - } - free(data_S); - double mean_S = mean(M, forward_data); - double squaredMean_S = squared_mean(M, forward_data); - double moment2_S = central_moment(M, forward_data, mean_S, 2); - double moment4_S = central_moment(M, forward_data, mean_S, 4); - - compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); - - sprintf(filename_S, "wolff_%lu_S_OO.dat", id); - - FILE *file_S = fopen(filename_S, "wb"); - fwrite(&M_f, sizeof(double), 1, file_S); - fwrite(&mean_S, sizeof(double), 1, file_S); - fwrite(&squaredMean_S, sizeof(double), 1, file_S); - fwrite(&moment2_S, sizeof(double), 1, file_S); - fwrite(&moment4_S, sizeof(double), 1, file_S); - fwrite(reverse_data, sizeof(double), length, file_S); - fclose(file_S); - - uint32_t *data_B = (uint32_t *)malloc((nb - 1) * N * sizeof(uint32_t)); - uint32_t *data_M = (uint32_t *)malloc((q - 1) * N * sizeof(uint32_t)); - fread(data_B, N * (nb - 1), sizeof(uint32_t), file_B); - fread(data_M, N * (q - 1), sizeof(uint32_t), file_M); - - for (count_t i = 0; i < M; i++) { - forward_data[i] = finite_energy(nb, J, q, H, nv, ne, data_B + (nb - 1) * (drop + i), data_M + (q - 1) * (drop + i)); - } - - double mean_E = mean(M, forward_data); - double squaredMean_E = squared_mean(M, forward_data); - double moment2_E = central_moment(M, forward_data, mean_E, 2); - double moment4_E = central_moment(M, forward_data, mean_E, 4); - - free(data_B); - free(data_M); - - compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); - - sprintf(filename_B, "wolff_%lu_E_OO.dat", id); - - FILE *file_E = fopen(filename_B, "wb"); - fwrite(&M_f, sizeof(double), 1, file_E); - fwrite(&mean_E, sizeof(double), 1, file_E); - fwrite(&squaredMean_E, sizeof(double), 1, file_E); - fwrite(&moment2_E, sizeof(double), 1, file_E); - fwrite(&moment4_E, sizeof(double), 1, file_E); - fwrite(reverse_data, sizeof(double), length, file_E); - fclose(file_E); - - printf("\033[F%lu: Correlation functions for %d steps written.\n", id, M); - - fftw_destroy_plan(forward_plan); - fftw_destroy_plan(reverse_plan); - fftw_free(forward_data); - fftw_free(reverse_data); - - } - - fclose(file_M); - fclose(file_B); - fclose(file_S); - - free(J); - free(H); - - free(filename_S); - free(filename_B); - free(filename_M); - - } else { - char *junk = (char *)malloc(1024 * sizeof(char)); - fscanf(metadata, "%[^\n]\n", junk); // throw away the rest of the line, we don't need it - free(junk); - - char *filename_E = (char *)malloc(128 * sizeof(char)); - char *filename_F = (char *)malloc(128 * sizeof(char)); - char *filename_M = (char *)malloc(128 * sizeof(char)); - char *filename_S = (char *)malloc(128 * sizeof(char)); - - sprintf(filename_E, "wolff_%lu_E.dat", id); - sprintf(filename_F, "wolff_%lu_F.dat", id); - sprintf(filename_M, "wolff_%lu_M.dat", id); - sprintf(filename_S, "wolff_%lu_S.dat", id); - - FILE *file_E = fopen(filename_E, "rb"); - FILE *file_F = fopen(filename_F, "rb"); - FILE *file_M = fopen(filename_M, "rb"); - FILE *file_S = fopen(filename_S, "rb"); - - fseek(file_S, 0, SEEK_END); - unsigned long N = ftell(file_S) / sizeof(uint32_t); - fseek(file_S, 0, SEEK_SET); - - if (speedy_drop) { - drop = N - pow(2, floor(log(N) / log(2))); - } else { - if (N % 2 == 1 && drop % 2 == 0) { - drop++; // make sure M is even - } - } - - if (N <= drop) { - printf("\033[F%lu: Number of steps %lu is less than %" PRIcount ", nothing done.\n", id, N, drop); - } else { - int M = N - drop; - double M_f = (double)M; - - if (length > M) { - length = M; - } - - double *forward_data = (double *)fftw_malloc(M * sizeof(double)); - fftw_plan forward_plan = fftw_plan_r2r_1d(M, forward_data, forward_data, FFTW_R2HC, 0); - - double *reverse_data = (double *)fftw_malloc(M * sizeof(double)); - fftw_plan reverse_plan = fftw_plan_r2r_1d(M, reverse_data, reverse_data, FFTW_HC2R, 0); - - if (file_S != NULL) { - uint32_t *data_S = (uint32_t *)malloc(N * sizeof(uint32_t)); - - fread(data_S, sizeof(uint32_t), N, file_S); - fclose(file_S); - - for (int i = 0; i < M; i++) { - forward_data[i] = (double)data_S[drop + i]; - } - free(data_S); - - double mean_S = mean(M, forward_data); - double squaredMean_S = squared_mean(M, forward_data); - double moment2_S = central_moment(M, forward_data, mean_S, 2); - double moment4_S = central_moment(M, forward_data, mean_S, 4); - - compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); - - sprintf(filename_S, "wolff_%lu_S_OO.dat", id); - FILE *file_S_new = fopen(filename_S, "wb"); - fwrite(&M_f, sizeof(double), 1, file_S_new); - fwrite(&mean_S, sizeof(double), 1, file_S_new); - fwrite(&squaredMean_S, sizeof(double), 1, file_S_new); - fwrite(&moment2_S, sizeof(double), 1, file_S_new); - fwrite(&moment4_S, sizeof(double), 1, file_S_new); - fwrite(reverse_data, sizeof(double), length, file_S_new); - fclose(file_S_new); - } - if (file_F != NULL) { - float *data_F = (float *)malloc(N * sizeof(float)); - - fread(data_F, sizeof(float), N, file_F); - fclose(file_F); - - for (int i = 0; i < M; i++) { - forward_data[i] = (double)data_F[drop + i]; - } - free(data_F); - - double mean_F = mean(M, forward_data); - double squaredMean_F = squared_mean(M, forward_data); - double moment2_F = central_moment(M, forward_data, mean_F, 2); - double moment4_F = central_moment(M, forward_data, mean_F, 4); - - compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); - - sprintf(filename_F, "wolff_%lu_F_OO.dat", id); - FILE *file_F_new = fopen(filename_F, "wb"); - fwrite(&M_f, sizeof(double), 1, file_F_new); - fwrite(&mean_F, sizeof(double), 1, file_F_new); - fwrite(&squaredMean_F, sizeof(double), 1, file_F_new); - fwrite(&moment2_F, sizeof(double), 1, file_F_new); - fwrite(&moment4_F, sizeof(double), 1, file_F_new); - fwrite(reverse_data, sizeof(double), length, file_F_new); - fclose(file_F_new); - } - if (file_E != NULL) { - float *data_E = (float *)malloc(N * sizeof(float)); - - fread(data_E, sizeof(float), N, file_E); - fclose(file_E); - - for (int i = 0; i < M; i++) { - forward_data[i] = (double)data_E[drop + i]; - } - free(data_E); - - double mean_E = mean(M, forward_data); - double squaredMean_E = squared_mean(M, forward_data); - double moment2_E = central_moment(M, forward_data, mean_E, 2); - double moment4_E = central_moment(M, forward_data, mean_E, 4); - - compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); - - sprintf(filename_E, "wolff_%lu_E_OO.dat", id); - FILE *file_E_new = fopen(filename_E, "wb"); - fwrite(&M_f, sizeof(double), 1, file_E_new); - fwrite(&mean_E, sizeof(double), 1, file_E_new); - fwrite(&squaredMean_E, sizeof(double), 1, file_E_new); - fwrite(&moment2_E, sizeof(double), 1, file_E_new); - fwrite(&moment4_E, sizeof(double), 1, file_E_new); - fwrite(reverse_data, sizeof(double), length, file_E_new); - fclose(file_E_new); - } - if (file_M != NULL) { - if (0 == strcmp(model, "PLANAR")) { - float *data_M = (float *)malloc(2 * N * sizeof(float)); - fread(data_M, sizeof(float), 2 * N, file_M); - fclose(file_M); - for (int i = 0; i < M; i++) { - forward_data[i] = (double)sqrt(pow(data_M[2 * drop + 2 * i], 2) + pow(data_M[2 * drop + 2 * i + 1], 2)); - } - free(data_M); - } else if (0 == strcmp(model, "HEISENBERG")) { - float *data_M = (float *)malloc(3 * N * sizeof(float)); - fread(data_M, sizeof(float), 3 * N, file_M); - fclose(file_M); - for (int i = 0; i < M; i++) { - forward_data[i] = sqrt(pow(data_M[3 * drop + 3 * i], 2) + pow(data_M[3 * drop + 3 * i + 1], 2) + pow(data_M[3 * drop + 3 * i + 2], 2)); - } - free(data_M); - } else if (0 == strcmp(model, "ISING")) { - int *data_M = (int *)malloc(N * sizeof(float)); - fread(data_M, sizeof(int), N, file_M); - fclose(file_M); - for (int i = 0; i < M; i++) { - forward_data[i] = (double)data_M[i]; - } - free(data_M); - } else { - printf("UNKNOWN MODEL\n"); - exit(EXIT_FAILURE); - } - - double mean_M = mean(M, forward_data); - double squaredMean_M = squared_mean(M, forward_data); - double moment2_M = central_moment(M, forward_data, mean_M, 2); - double moment4_M = central_moment(M, forward_data, mean_M, 4); - - compute_OO(M, forward_plan, forward_data, reverse_plan, reverse_data); - - sprintf(filename_M, "wolff_%lu_M_OO.dat", id); - FILE *file_M_new = fopen(filename_M, "wb"); - fwrite(&M_f, sizeof(double), 1, file_M_new); - fwrite(&mean_M, sizeof(double), 1, file_M_new); - fwrite(&squaredMean_M, sizeof(double), 1, file_M_new); - fwrite(&moment2_M, sizeof(double), 1, file_M_new); - fwrite(&moment4_M, sizeof(double), 1, file_M_new); - fwrite(reverse_data, sizeof(double), length, file_M_new); - fclose(file_M_new); - } - - printf("\033[F%lu: Correlation functions for %d steps written.\n", id, M); - fftw_destroy_plan(forward_plan); - fftw_destroy_plan(reverse_plan); - fftw_free(forward_data); - fftw_free(reverse_data); - - } - free(filename_E); - free(filename_S); - free(filename_F); - free(filename_M); - } - } - - free(model); - fclose(metadata); - fftw_cleanup(); - - return 0; -} - diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp deleted file mode 100644 index f6661af..0000000 --- a/src/wolff_On.cpp +++ /dev/null @@ -1,269 +0,0 @@ - -#include -#include - -#ifdef HAVE_GLUT -#include -#endif - -#include -#include - -#include -#include -#include -#include - -typedef orthogonal_t orthogonal_R_t; -typedef vector_t vector_R_t; -typedef state_t On_t; - -// angle from the x-axis of a two-vector -double theta(vector_R_t v) { - double x = v[0]; - double y = v[1]; - - double val = atan(y / x); - - if (x < 0.0 && y > 0.0) { - return M_PI + val; - } else if ( x < 0.0 && y < 0.0 ) { - return - M_PI + val; - } else { - return val; - } -} - -double H_modulated(vector_R_t v, int order, double mag) { - return mag * cos(order * theta(v)); -} - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e7; - -#ifdef DIMENSION - D_t D = DIMENSION; -#else - D_t D = 2; -#endif - L_t L = 128; - double T = 2.26918531421; - double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); - - bool silent = false; - bool use_pert = false; - bool N_is_sweeps = false; - bool draw = false; - unsigned int window_size = 512; - - bool modulated_field = false; - unsigned int order = 1; - - int opt; - q_t H_ind = 0; - double epsilon = 1; - -// unsigned char measurement_flags = measurement_energy | measurement_clusterSize; - - unsigned char measurement_flags = 0; - - while ((opt = getopt(argc, argv, "N:D:L:T:H:spe:mo:M:Sdw:")) != -1) { - switch (opt) { - case 'N': // number of steps - N = (count_t)atof(optarg); - break; -#ifdef DIMENSION - case 'D': // dimension - printf("Dimension was specified at compile time, you can't change it now!\n"); - exit(EXIT_FAILURE); -#else - case 'D': // dimension - D = atoi(optarg); - break; -#endif - case 'L': // linear size - L = atoi(optarg); - break; - case 'T': // temperature - T = atof(optarg); - break; - case 'H': // external field. nth call couples to state n - H_vec[H_ind] = atof(optarg); - H_ind++; - break; - case 's': // don't print anything during simulation. speeds up slightly - silent = true; - break; - case 'p': - use_pert = true; - break; - case 'e': - epsilon = atof(optarg); - break; - case 'm': - modulated_field = true; - break; - case 'M': - measurement_flags ^= 1 << atoi(optarg); - break; - case 'o': - order = atoi(optarg); - break; - case 'S': - N_is_sweeps = true; - break; - case 'd': -#ifdef HAVE_GLUT - draw = true; - break; -#else - printf("You didn't compile this with the glut library installed!\n"); - exit(EXIT_FAILURE); -#endif - case 'w': - window_size = atoi(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - unsigned long timestamp; - - { - struct timespec spec; - clock_gettime(CLOCK_REALTIME, &spec); - timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; - } - - const char *pert_type; - - std::function gen_R; - - if (use_pert) { - double Hish; - if (modulated_field) { - Hish = fabs(H_vec[0]); - } else { - double H2 = 0; - for (q_t i = 0; i < N_COMP; i++) { - H2 += pow(H_vec[i], 2); - } - Hish = sqrt(H2); - } - - epsilon = sqrt((N_COMP - 1) * T / (D + Hish / 2)) / 2; - - gen_R = std::bind(generate_rotation_perturbation , std::placeholders::_1, std::placeholders::_2, epsilon, order); - pert_type = "PERTURB5"; - } else { - gen_R = generate_rotation_uniform ; - pert_type = "UNIFORM"; - } - - FILE *outfile_info = fopen("wolff_metadata.txt", "a"); - - fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"%s\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"FIELD_TYPE\" -> ", timestamp, ON_strings[N_COMP], N_COMP, D, L, (v_t)pow(L, D), D * (v_t)pow(L, D), T); - if (modulated_field) { - fprintf(outfile_info, "\"MODULATED\", \"ORDER\" -> %d, \"H\" -> %.15f, ", order, H_vec[0]); - } else { - fprintf(outfile_info, "\"VECTOR\", \"H\" -> {"); - for (q_t i = 0; i < N_COMP; i++) { - fprintf(outfile_info, "%.15f", H_vec[i]); - if (i < N_COMP - 1) { - fprintf(outfile_info, ", "); - } - } - fprintf(outfile_info, "}, "); - } - - fprintf(outfile_info, "\"GENERATOR\" -> \"%s\"", pert_type); - - if (use_pert) { - fprintf(outfile_info, ", \"EPS\" -> %g", epsilon); - } - - fprintf(outfile_info, " |>\n"); - - fclose(outfile_info); - - FILE **outfiles = measure_setup_files(measurement_flags, timestamp); - - 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; - }; - } else if (draw) { -#ifdef HAVE_GLUT - // initialize glut - glutInit(&argc, argv); - glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(window_size, window_size); - glutCreateWindow("wolff"); - glClearColor(0.0,0.0,0.0,0.0); - glMatrixMode(GL_PROJECTION); - glLoadIdentity(); - gluOrtho2D(0.0, L, 0.0, L); - - other_f = [&] (const On_t& s) { - glClear(GL_COLOR_BUFFER_BIT); - for (v_t i = 0; i < pow(L, 2); i++) { - vector_R_t v_tmp = s.R.act_inverse(s.spins[i]); - double thetai = fmod(2 * M_PI + theta(v_tmp), 2 * M_PI); - double saturation = 0.7; - double value = 0.9; - double chroma = saturation * value; - glColor3f(chroma * hue_to_R(thetai) + (value - chroma), chroma * hue_to_G(thetai) + (value - chroma), chroma * hue_to_B(thetai) + (value - chroma)); - glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); - } - glFlush(); - }; -#endif - } else { - other_f = [] (const On_t& s) {}; - } - - std::function measurements = measure_function_write_files(measurement_flags, outfiles, other_f); - - std::function H; - - if (modulated_field) { - H = std::bind(H_modulated, std::placeholders::_1, order, H_vec[0]); - } else { - H = std::bind(H_vector , std::placeholders::_1, H_vec); - } - - // initialize random number generator - gsl_rng *r = gsl_rng_alloc(gsl_rng_taus2); - gsl_rng_set(r, rand_seed()); - -#ifndef NOFIELD - state_t s(D, L, T, dot , H); -#else - state_t s(D, L, T, dot ); -#endif - - 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, r, silent); - 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); - } else { - wolff (N, s, gen_R, measurements, r, silent); - } - - measure_free_files(measurement_flags, outfiles); - free(H_vec); - gsl_rng_free(r); - - return 0; -} - diff --git a/src/wolff_cgm.cpp b/src/wolff_cgm.cpp deleted file mode 100644 index ce91bf2..0000000 --- a/src/wolff_cgm.cpp +++ /dev/null @@ -1,166 +0,0 @@ - -#include - -#ifdef HAVE_GLUT -#include -#endif - -// include your group and spin space -#include -#include - -// include wolff.h -#include -#include - -typedef state_t , height_t> sim_t; - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e4; - - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double H = 0; - - bool silent = false; - bool draw = false; - unsigned int window_size = 512; - double epsilon = 1; - - int opt; - - while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:e:")) != -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. nth call couples to state n - H = atof(optarg); - break; - case 'e': // external field. nth call couples to state n - epsilon = atof(optarg); - break; - case 's': // don't print anything during simulation. speeds up slightly - silent = true; - break; - case 'd': -#ifdef HAVE_GLUT - draw = true; - break; -#else - printf("You didn't compile this with the glut library installed!\n"); - exit(EXIT_FAILURE); -#endif - case 'w': - window_size = atoi(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - // initialize random number generator - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - - // define spin-spin coupling - std::function &, const height_t&)> Z = [] (const height_t& h1, const height_t& h2) -> double { - return -pow(h1.x - h2.x, 2); - }; - - // define spin-field coupling - std::function )> B = [=] (height_t h) -> double { - return -H * pow(h.x, 2);; - }; - - // initialize state object - sim_t s(D, L, T, Z, B); - - // define function that generates self-inverse rotations - std::function (gsl_rng *, height_t)> gen_R = [=] (gsl_rng *r, height_t h) -> dihedral_inf_t { - dihedral_inf_t rot; - rot.is_reflection = true; - - double amount = epsilon * gsl_ran_ugaussian(r); - - rot.x = 2 * h.x + amount; - - return rot; - }; - - // define function that updates any number of measurements - std::function 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 / (double)N / (double)s.nv; - }; - } else { - // a more complex example: measure the average magnetization, and draw the spin configuration to the screen - -#ifdef HAVE_GLUT - // initialize glut - glutInit(&argc, argv); - glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(window_size, window_size); - glutCreateWindow("wolff"); - glClearColor(0.0,0.0,0.0,0.0); - glMatrixMode(GL_PROJECTION); - glLoadIdentity(); - gluOrtho2D(0.0, L, 0.0, L); - - measurement = [&] (const sim_t& s) { - average_M += (double)s.M / (double)N / (double)s.nv; - glClear(GL_COLOR_BUFFER_BIT); - double max_h = INT64_MIN; - double min_h = INT64_MAX; - for (v_t i = 0; i < pow(L, 2); i++) { - double cur_h = (s.R.act_inverse(s.spins[i])).x; - if (cur_h < min_h) { - min_h = cur_h; - } - if (cur_h > max_h) { - max_h = cur_h; - } - } - - for (v_t i = 0; i < pow(L, 2); i++) { - double cur_h = (s.R.act_inverse(s.spins[i])).x; - double mag = ((double)(cur_h - min_h)) / ((double)(max_h - min_h)); - glColor3f(mag, mag, mag); - glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); - } - glFlush(); - }; -#endif - } - - // run wolff for N cluster flips - wolff(N, s, gen_R, measurement, r, silent); - - // tell us what we found! - printf("%" PRIcount " DGM runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, = %g\n", N, D, L, T, H, average_M); - - // free the random number generator - gsl_rng_free(r); - - if (draw) { - } - - return 0; - -} - diff --git a/src/wolff_clock.cpp b/src/wolff_clock.cpp deleted file mode 100644 index 3dec284..0000000 --- a/src/wolff_clock.cpp +++ /dev/null @@ -1,154 +0,0 @@ - -#include - -#ifdef HAVE_GLUT -#include -#endif - -// include your group and spin space -#include -#include -#include - -// hack to speed things up considerably -#define N_STATES POTTSQ -#include - -// include wolff.h -#include -#include - -typedef state_t , potts_t> sim_t; - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e4; - - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); - - bool silent = false; - bool draw = false; - unsigned int window_size = 512; - - int opt; - q_t H_ind = 0; - - while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:")) != -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. nth call couples to state n - H_vec[H_ind] = atof(optarg); - H_ind++; - break; - case 's': // don't print anything during simulation. speeds up slightly - silent = true; - break; - case 'd': -#ifdef HAVE_GLUT - draw = true; - break; -#else - printf("You didn't compile this with the glut library installed!\n"); - exit(EXIT_FAILURE); -#endif - case 'w': - window_size = atoi(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - // initialize random number generator - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - - // define spin-spin coupling - std::function &, const potts_t&)> Z = [] (const potts_t& s1, const potts_t& s2) -> double { - return cos(2 * M_PI * (double)(s1.x + POTTSQ - s2.x) / (double)POTTSQ); - }; - - // define spin-field coupling - std::function &)> B = [=] (const potts_t& s) -> double { - return H_vec[s.x]; - }; - - // initialize state object - state_t , potts_t> s(D, L, T, Z, B); - - // define function that generates self-inverse rotations - std::function (gsl_rng *, potts_t)> gen_R = [] (gsl_rng *r, potts_t v) -> dihedral_t { - dihedral_t rot; - rot.is_reflection = true; - q_t x = gsl_rng_uniform_int(r, POTTSQ - 1); - rot.x = (2 * v.x + x + 1) % POTTSQ; - - return rot; - }; - - // define function that updates any number of measurements - std::function 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; - }; - } else { - // a more complex example: measure the average magnetization, and draw the spin configuration to the screen - -#ifdef HAVE_GLUT - // initialize glut - glutInit(&argc, argv); - glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(window_size, window_size); - glutCreateWindow("wolff"); - glClearColor(0.0,0.0,0.0,0.0); - glMatrixMode(GL_PROJECTION); - glLoadIdentity(); - gluOrtho2D(0.0, L, 0.0, L); - - measurement = [&] (const sim_t& s) { - average_M += (double)s.M[0] / (double)N / (double)s.nv; - 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]); - glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ)); - glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); - } - glFlush(); - }; -#endif - } - - // run wolff for N cluster flips - wolff(N, s, gen_R, measurement, r, 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); - - // free the random number generator - gsl_rng_free(r); - - if (draw) { - } - - return 0; - -} - diff --git a/src/wolff_dgm.cpp b/src/wolff_dgm.cpp deleted file mode 100644 index 8667fb5..0000000 --- a/src/wolff_dgm.cpp +++ /dev/null @@ -1,171 +0,0 @@ - -#include - -#ifdef HAVE_GLUT -#include -#endif - -// include your group and spin space -#include -#include - -// include wolff.h -#include -#include - -typedef state_t , height_t> sim_t; - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e4; - - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double H = 0; - - bool silent = false; - bool draw = false; - unsigned int window_size = 512; - uint64_t epsilon = 1; - - int opt; - - while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:e:")) != -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. nth call couples to state n - H = atof(optarg); - break; - case 'e': // external field. nth call couples to state n - epsilon = atof(optarg); - break; - case 's': // don't print anything during simulation. speeds up slightly - silent = true; - break; - case 'd': -#ifdef HAVE_GLUT - draw = true; - break; -#else - printf("You didn't compile this with the glut library installed!\n"); - exit(EXIT_FAILURE); -#endif - case 'w': - window_size = atoi(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - // initialize random number generator - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - - // define spin-spin coupling - std::function &, const height_t&)> Z = [] (const height_t& h1, const height_t& h2) -> double { - return -pow(h1.x - h2.x, 2); - }; - - // define spin-field coupling - std::function &)> B = [=] (const height_t& h) -> double { - return -H * pow(h.x, 2);; - }; - - // initialize state object - sim_t s(D, L, T, Z, B); - - // define function that generates self-inverse rotations - std::function (gsl_rng *, height_t)> gen_R = [=] (gsl_rng *r, height_t h) -> dihedral_inf_t { - dihedral_inf_t rot; - rot.is_reflection = true; - - int direction = gsl_rng_uniform_int(r, 2); - int64_t amount = gsl_rng_uniform_int(r, epsilon); - - if (direction == 0) { - rot.x = 2 * h.x + (1 + amount); - } else { - rot.x = 2 * h.x - (1 + amount); - } - - return rot; - }; - - // define function that updates any number of measurements - std::function 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 / (double)N / (double)s.nv; - }; - } else { - // a more complex example: measure the average magnetization, and draw the spin configuration to the screen - -#ifdef HAVE_GLUT - // initialize glut - glutInit(&argc, argv); - glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(window_size, window_size); - glutCreateWindow("wolff"); - glClearColor(0.0,0.0,0.0,0.0); - glMatrixMode(GL_PROJECTION); - glLoadIdentity(); - gluOrtho2D(0.0, L, 0.0, L); - - measurement = [&] (const sim_t& s) { - average_M += (double)s.M / (double)N / (double)s.nv; - glClear(GL_COLOR_BUFFER_BIT); - int64_t max_h = INT64_MIN; - int64_t min_h = INT64_MAX; - for (v_t i = 0; i < pow(L, 2); i++) { - int64_t cur_h = (s.R.act_inverse(s.spins[i])).x; - if (cur_h < min_h) { - min_h = cur_h; - } - if (cur_h > max_h) { - max_h = cur_h; - } - } - - for (v_t i = 0; i < pow(L, 2); i++) { - int64_t cur_h = (s.R.act_inverse(s.spins[i])).x; - double mag = ((double)(cur_h - min_h)) / ((double)(max_h - min_h)); - glColor3f(mag, mag, mag); - glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); - } - glFlush(); - }; -#endif - } - - // run wolff for N cluster flips - wolff(N, s, gen_R, measurement, r, silent); - - // tell us what we found! - printf("%" PRIcount " DGM runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, = %g\n", N, D, L, T, H, average_M); - - // free the random number generator - gsl_rng_free(r); - - if (draw) { - } - - return 0; - -} - diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp deleted file mode 100644 index a6f43b1..0000000 --- a/src/wolff_ising.cpp +++ /dev/null @@ -1,201 +0,0 @@ - -#include -#include - -// if you have GLUT installed, you can see graphics! -#ifdef HAVE_GLUT -#include -#endif - -// include your group and spin space -#include -#include - -// finite_states.h can be included for spin types that have special variables -// defined, and it causes wolff execution to use precomputed bond probabilities -#include - -// rand.h uses a unix-specific way to seed the random number generator -#include - -// measure.h contains useful functions for saving timeseries to files -#include - -// include wolff.h -#include - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e4; - - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double H = 0.0; - - bool silent = false; - bool draw = false; - bool N_is_sweeps = false; - unsigned int window_size = 512; - - // don't measure anything by default - unsigned char measurement_flags = 0; - - int opt; - - while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -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; - case 's': // don't print anything during simulation. speeds up slightly - silent = true; - break; - case 'S': - N_is_sweeps = true; - break; - case 'd': -#ifdef HAVE_GLUT - draw = true; - break; -#else - printf("You didn't compile this with the glut library installed!\n"); - exit(EXIT_FAILURE); -#endif - case 'w': - window_size = atoi(optarg); - break; - case 'M': - measurement_flags ^= 1 << atoi(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - // get nanosecond timestamp for unique run id - unsigned long timestamp; - - { - struct timespec spec; - clock_gettime(CLOCK_REALTIME, &spec); - timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; - } - - // initialize random number generator - gsl_rng *r = gsl_rng_alloc(gsl_rng_taus2); - gsl_rng_set(r, rand_seed()); - - // define 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 spin-field coupling - std::function B = [=] (const ising_t& s) -> double { - if (s.x) { - return -H; - } else { - return H; - } - }; - - // initialize state object -#ifndef NOFIELD - state_t s(D, L, T, Z, B); -#else - state_t s(D, L, T, Z); -#endif - - // define function that generates self-inverse rotations - std::function gen_R = [] (gsl_rng *, const ising_t& s) -> z2_t { - return z2_t(true); - }; - - FILE **outfiles = measure_setup_files(measurement_flags, timestamp); - - std::function &)> 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; - }; - } else if (draw) { -#ifdef HAVE_GLUT - // initialize glut - glutInit(&argc, argv); - glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(window_size, window_size); - glutCreateWindow("wolff"); - glClearColor(0.0,0.0,0.0,0.0); - glMatrixMode(GL_PROJECTION); - glLoadIdentity(); - gluOrtho2D(0.0, L, 0.0, L); - - other_f = [] (const state_t & s) { - glClear(GL_COLOR_BUFFER_BIT); - for (v_t i = 0; i < pow(s.L, 2); i++) { - if (s.spins[i].x == s.R.x) { - glColor3f(0.0, 0.0, 0.0); - } else { - glColor3f(1.0, 1.0, 1.0); - } - glRecti(i / s.L, i % s.L, (i / s.L) + 1, (i % s.L) + 1); - } - glFlush(); - }; -#endif - } else { - other_f = [] (const state_t& s) {}; - } - - std::function &)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f); - - // add line to metadata file with run info - { - FILE *outfile_info = fopen("wolff_metadata.txt", "a"); - - fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"ISING\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> %.15f |>\n", timestamp, s.D, s.L, s.nv, s.ne, T, H); - - fclose(outfile_info); - } - - // run wolff for N cluster flips - 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, r, silent); - 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); - } else { - wolff(N, s, gen_R, measurements, r, silent); - } - - // free the random number generator - gsl_rng_free(r); - measure_free_files(measurement_flags, outfiles); - - return 0; - -} - diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp deleted file mode 100644 index 9fe3ffe..0000000 --- a/src/wolff_potts.cpp +++ /dev/null @@ -1,213 +0,0 @@ - -#include -#include - -#ifdef HAVE_GLUT -#include -#endif - -// include your group and spin space -#include -#include - -// hack to speed things up considerably -#define N_STATES POTTSQ -#include - -// include wolff.h -#include -#include -#include -#include - -typedef state_t , potts_t> sim_t; - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e4; - - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); - - bool silent = false; - bool draw = false; - bool N_is_sweeps = false; - unsigned int window_size = 512; - - // don't measure anything by default - unsigned char measurement_flags = 0; - - int opt; - q_t H_ind = 0; - - while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -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. nth call couples to state n - H_vec[H_ind] = atof(optarg); - H_ind++; - break; - case 's': // don't print anything during simulation. speeds up slightly - silent = true; - break; - case 'S': - N_is_sweeps = true; - break; - case 'd': -#ifdef HAVE_GLUT - draw = true; - break; -#else - printf("You didn't compile this with the glut library installed!\n"); - exit(EXIT_FAILURE); -#endif - case 'w': - window_size = atoi(optarg); - break; - case 'M': - measurement_flags ^= 1 << atoi(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - // get nanosecond timestamp for unique run id - unsigned long timestamp; - - { - struct timespec spec; - clock_gettime(CLOCK_REALTIME, &spec); - timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; - } - - // initialize random number generator - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - - // define spin-spin coupling - std::function &, const potts_t&)> Z = [] (const potts_t& s1, const potts_t& s2) -> double { - if (s1.x == s2.x) { - return 1.0; - } else { - return 0.0; - } - }; - - // define spin-field coupling - std::function &)> B = [=] (const potts_t& s) -> double { - return H_vec[s.x]; - }; - - // initialize state object - state_t , potts_t> s(D, L, T, Z, B); - - // define function that generates self-inverse rotations - std::function (gsl_rng *, potts_t)> gen_R = [] (gsl_rng *r, potts_t v) -> symmetric_t { - symmetric_t rot; - - q_t j = gsl_rng_uniform_int(r, POTTSQ - 1); - q_t swap_v; - if (j < v.x) { - swap_v = j; - } else { - swap_v = j + 1; - } - - rot[v.x] = swap_v; - rot[swap_v] = v.x; - - return rot; - }; - - FILE **outfiles = measure_setup_files(measurement_flags, timestamp); - - std::function 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; - }; - } else if (draw) { -#ifdef HAVE_GLUT - // initialize glut - glutInit(&argc, argv); - glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(window_size, window_size); - glutCreateWindow("wolff"); - glClearColor(0.0,0.0,0.0,0.0); - glMatrixMode(GL_PROJECTION); - glLoadIdentity(); - gluOrtho2D(0.0, L, 0.0, L); - - other_f = [] (const sim_t& s) { - 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]); - glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ)); - glRecti(i / s.L, i % s.L, (i / s.L) + 1, (i % s.L) + 1); - } - glFlush(); - }; -#endif - } else { - other_f = [] (const sim_t& s) {}; - } - - std::function measurements = measure_function_write_files(measurement_flags, outfiles, other_f); - - // add line to metadata file with run info - { - FILE *outfile_info = fopen("wolff_metadata.txt", "a"); - - fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"POTTS\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, POTTSQ, s.D, s.L, s.nv, s.ne, T); - - for (q_t i = 0; i < POTTSQ; i++) { - fprintf(outfile_info, "%.15f", H_vec[i]); - if (i < POTTSQ - 1) { - fprintf(outfile_info, ", "); - } - } - - fprintf(outfile_info, "} |>\n"); - - fclose(outfile_info); - } - - // run wolff for N cluster flips - 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, r, silent); - 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); - } else { - wolff(N, s, gen_R, measurements, r, silent); - } - - // free the random number generator - gsl_rng_free(r); - free(H_vec); - measure_free_files(measurement_flags, outfiles); - - return 0; - -} - diff --git a/wolfram_link/Makefile b/wolfram_link/Makefile deleted file mode 100644 index 9d65623..0000000 --- a/wolfram_link/Makefile +++ /dev/null @@ -1,13 +0,0 @@ - -CC = clang -WSPREP = /opt/Mathematica/SystemFiles/Links/WSTP/DeveloperKit/Linux-x86-64/CompilerAdditions/wsprep -CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -flto -fopenmp=libiomp5 -I/usr/lib/gcc/x86_64-linux-gnu/4.8/include/ -march=native -I/opt/Mathematica/SystemFiles/Links/WSTP/DeveloperKit/Linux-x86-64/CompilerAdditions/ -I../lib/ -LDFLAGS = -L/opt/Mathematica/SystemFiles/Links/WSTP/DeveloperKit/Linux-x86-64/CompilerAdditions -lm -lpthread -lrt -lstdc++ -ldl -luuid -l WSTP64i4 - -convex.o: ../lib/convex.c - @${CC} -c -o $@ $< ${CFLAGS} - -convexminorant: convexminorant.tm convex.o - @${WSPREP} $< -o $@.c - @${CC} -o $@ convex.o $@.c ${CFLAGS} ${LDFLAGS} - diff --git a/wolfram_link/convexminorant.tm b/wolfram_link/convexminorant.tm deleted file mode 100644 index edca6f4..0000000 --- a/wolfram_link/convexminorant.tm +++ /dev/null @@ -1,34 +0,0 @@ - -#include -#include - -extern int WSMain(int, char **); - -void convexminorant(double *, int); - -:Begin: -:Function: convexminorant -:Pattern: GetConvexMinorant[ list:{___Real} ] -:Arguments: { list } -:ArgumentTypes: { Real64List } -:ReturnType: Manual -:End: - -:Evaluate: GetConvexMinorant[ sequence___Float]:= GetConvexMinorant[ {sequence} ] - -void convexminorant(double * Gammas, int len) { - int i; - for (i = 0; i < len; i++) { - if (Gammas[i] <= 0) { - break; - } - } - double *m = get_convex_minorant(i, Gammas); - WSPutReal64List(stdlink, m, i); - free(m); -} - -int main(int argc, char **argv) { - return WSMain(argc, argv); -} - -- cgit v1.2.3-54-g00ecf