From 643baba78eb15a685d959aae718ee3eeade2f806 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 19 Oct 2018 01:26:14 -0400 Subject: big library overhual, fixed all examples, added documentation with sphinx --- .gitignore | 4 +- doc/Makefile | 19 ++++ doc/compile.rst | 51 +++++++++ doc/conf.py | 177 ++++++++++++++++++++++++++++++++ doc/finite_states.rst | 24 +++++ doc/graph.rst | 47 +++++++++ doc/index.rst | 15 +++ doc/measurement.rst | 84 +++++++++++++++ doc/models.rst | 61 +++++++++++ doc/system.rst | 61 +++++++++++ doc/types.rst | 36 +++++++ examples/CMakeLists.txt | 4 +- examples/On.cpp | 10 +- examples/animate_ising.cpp | 133 ------------------------ examples/ising.cpp | 14 ++- examples/ising_animation.cpp | 139 +++++++++++++++++++++++++ examples/simple_measurement.hpp | 18 ++-- lib/include/wolff.hpp | 21 ++-- lib/include/wolff/cluster.hpp | 57 +++++----- lib/include/wolff/finite_states.hpp | 6 +- lib/include/wolff/graph.hpp | 21 ++-- lib/include/wolff/measurement.hpp | 18 ++-- lib/include/wolff/models/ising.hpp | 11 +- lib/include/wolff/models/orthogonal.hpp | 15 +-- lib/include/wolff/models/potts.hpp | 2 +- lib/include/wolff/system.hpp | 30 ++++-- lib/include/wolff/types.h | 2 +- lib/src/graph.cpp | 79 ++++++-------- 28 files changed, 882 insertions(+), 277 deletions(-) create mode 100644 doc/Makefile create mode 100644 doc/compile.rst create mode 100644 doc/conf.py create mode 100644 doc/finite_states.rst create mode 100644 doc/graph.rst create mode 100644 doc/index.rst create mode 100644 doc/measurement.rst create mode 100644 doc/models.rst create mode 100644 doc/system.rst create mode 100644 doc/types.rst delete mode 100644 examples/animate_ising.cpp create mode 100644 examples/ising_animation.cpp diff --git a/.gitignore b/.gitignore index bf5bfbe..e7bdd38 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,2 @@ build/* -wolfram_link/*.o -wolfram_link/*.c -wolfram_link/convexminorant +doc/_build/* diff --git a/doc/Makefile b/doc/Makefile new file mode 100644 index 0000000..298ea9e --- /dev/null +++ b/doc/Makefile @@ -0,0 +1,19 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) \ No newline at end of file diff --git a/doc/compile.rst b/doc/compile.rst new file mode 100644 index 0000000..2a3b4a2 --- /dev/null +++ b/doc/compile.rst @@ -0,0 +1,51 @@ + +.. _compile: + +*************** +Compile Options +*************** + +There are several features that are not active by default, but can be used by defining a compiler flag. These will modify the constructor of :cpp:class:`system`. These can be used in concert when applicable. Note that at the moment, there is no :ref:`finite_states` support for systems with bond or site dependence. + +Building Without A Field +======================== + +.. c:macro:: WOLFF_NO_FIELD + +When this flag is defined Wolff will not use a field. When a field isn't necessary, this will improve performance: no ghost site will be initialize and no time will be wasted checking the energy change with respect to a uncoupled ghost site. + +When defined, the constructor for :cpp:class:`system` becomes + +.. cpp:namespace:: cflags + +.. cpp:function:: system(graph, double, std::function ) + +The resulting :cpp:class:`system` object does not have member objects :cpp:member:`B` or :cpp:member:`s0`, and :cpp:member:`system::G` does not have a ghost site initialized. + +Building With Bond Dependence +============================= + +.. c:macro:: WOLFF_BOND_DEPENDENCE + +When this flag is defined Wolff will ask for the indices of the spins on each side a bond, allowing the implementation of random bonds or anisotropic interactions. + +When defined, the bond coupling must be a function of the form + +.. cpp:function:: double Z(v_t, const X_t&, v_t, const X_t&) + +where the first argument is the index of the first spin and the third argument is the index of the second spin. A function of this type is passed to :cpp:class:`system` in place of the original bond coupling. + + +Building With Site Dependence +============================= + +.. c:macro:: WOLFF_SITE_DEPENDENCE + +When this flag is defined Wolff will ask for the indices of the spin when measuring the external field, allowing the implementation of random fields or to emulate boundaries. + +When defined, the field coupling must be a function of the form + +.. cpp:function:: double B(v_t, const X_t&) + +where the first argument is the index of the spin. A function of this type is passed to :cpp:class:`system` in place of the original field coupling. + diff --git a/doc/conf.py b/doc/conf.py new file mode 100644 index 0000000..112d4b7 --- /dev/null +++ b/doc/conf.py @@ -0,0 +1,177 @@ +# -*- coding: utf-8 -*- +# +# Configuration file for the Sphinx documentation builder. +# +# This file does only contain a selection of the most common options. For a +# full list see the documentation: +# http://www.sphinx-doc.org/en/master/config + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +# import os +# import sys +# sys.path.insert(0, os.path.abspath('.')) + + +# -- Project information ----------------------------------------------------- + +project = 'wolff' +copyright = '2018, Jaron Kent-Dobias' +author = 'Jaron Kent-Dobias' + +# The short X.Y version +version = '' +# The full version, including alpha/beta/rc tags +release = '' + + +# -- General configuration --------------------------------------------------- + +# If your documentation needs a minimal Sphinx version, state it here. +# +# needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.mathjax', +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# +# source_suffix = ['.rst', '.md'] +source_suffix = '.rst' + +# The master toctree document. +master_doc = 'index' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = None + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = None + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'alabaster' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +# html_theme_options = {} + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +# Custom sidebar templates, must be a dictionary that maps document names +# to template names. +# +# The default sidebars (for documents that don't match any pattern) are +# defined by theme itself. Builtin themes are using these templates by +# default: ``['localtoc.html', 'relations.html', 'sourcelink.html', +# 'searchbox.html']``. +# +# html_sidebars = {} + + +# -- Options for HTMLHelp output --------------------------------------------- + +# Output file base name for HTML help builder. +htmlhelp_basename = 'wolffdoc' + + +# -- Options for LaTeX output ------------------------------------------------ + +latex_elements = { + # The paper size ('letterpaper' or 'a4paper'). + # + # 'papersize': 'letterpaper', + + # The font size ('10pt', '11pt' or '12pt'). + # + # 'pointsize': '10pt', + + # Additional stuff for the LaTeX preamble. + # + # 'preamble': '', + + # Latex figure (float) alignment + # + # 'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'wolff.tex', 'wolff Documentation', + 'Jaron Kent-Dobias', 'manual'), +] + + +# -- Options for manual page output ------------------------------------------ + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [ + (master_doc, 'wolff', 'wolff Documentation', + [author], 1) +] + + +# -- Options for Texinfo output ---------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + (master_doc, 'wolff', 'wolff Documentation', + author, 'wolff', 'One line description of project.', + 'Miscellaneous'), +] + + +# -- Options for Epub output ------------------------------------------------- + +# Bibliographic Dublin Core info. +epub_title = project + +# The unique identifier of the text. This can be a ISBN number +# or the project homepage. +# +# epub_identifier = '' + +# A unique identification for the text. +# +# epub_uid = '' + +# A list of files that should not be packed into the epub file. +epub_exclude_files = ['search.html'] + + +# -- Extension configuration ------------------------------------------------- \ No newline at end of file diff --git a/doc/finite_states.rst b/doc/finite_states.rst new file mode 100644 index 0000000..5d1fce8 --- /dev/null +++ b/doc/finite_states.rst @@ -0,0 +1,24 @@ + +.. _finite_states: + +************* +Finite States +************* + +One of the slower steps in running the algorithm is taking the exponent involved in computing the bond activation probabilities for each prospective bond. When the spins in your system have a finite number of possible states, the algorithm can be sped up considerably by precomputing the bond activation probabilities for every possible combination pair spins. Once the appropriate things have been defined for your model, the header :file:`wolff/finite_states.hpp` can be invoked to automate this process. The provided model headers :file:`wolff/models/ising.hpp` and :file:`wolff/models/potts.hpp` demonstrate the expected usage. + +Required Definitions +==================== + +.. c:macro:: WOLFF_FINITE_STATES_N + +This macro must be defined and given a value equal to the number of states your model can take. + +.. cpp:var:: const X_t finite_states_possible[WOLFF_FINITE_STATES_N] + +You must supply a constant array which lists each of the possible states of your individual spins. + +.. cpp:function:: q_t finite_states_enum(const X_t&) + +You must define a function which takes each state and returns its index in :cpp:var:`finite_states_possible`. + diff --git a/doc/graph.rst b/doc/graph.rst new file mode 100644 index 0000000..e44a14f --- /dev/null +++ b/doc/graph.rst @@ -0,0 +1,47 @@ + +****** +Graphs +****** + +This class is defined in the header :file:`wolff/graph.hpp`. + +.. cpp:class:: graph + + Lattices are described by objects of class :cpp:class:`graph`, a minimal implementation of graphs. + +.. cpp:member:: D_t graph::D + + The dimension of the graph. This property is unused by the core library. + +.. cpp:member:: L_t graph::L + + The linear size of the graph. This property is unused by the core library. + +.. cpp:member:: v_t graph::ne + + The number of edges in the graph. This property is unused by the core library. + +.. cpp:member:: v_t graph::nv + + The number of vertices in the graph. + +.. cpp:member:: std::vector> graph::adjacency + + The adjacency list for the graph. The :math:`i\text{th}` element of :cpp:member:`adjacency` is a standard library vector containing the indices of all vertices adjacent to vertex :math:`i`. + +.. cpp:member:: std::vector> graph::coordinate + + The coordinates of the graph vertices. The :math:`i\text{th}` element of :cpp:var:`coordinate` is a standard library vector of length :cpp:var:`D` containing the spatial coordinates of vertex :math:`i`. This property is unused by the core library. + +.. cpp:function:: graph::graph() + + The default constructor. Initializes an empty graph, i.e., :cpp:var:`D`, :cpp:var:`L`, :cpp:var:`ne`, and :cpp:var:`nv` are all zero and :cpp:var:`adjacency` and :cpp:var:`coordinate` are uninitialized. + +.. cpp:function:: graph::graph(D_t D, L_t L) + + Calling the constructor with :cpp:var:`D` and :cpp:var:`L` will initialize the graph of a :cpp:var:`D`-dimensional hypercubic lattice with :cpp:var:`L` vertices per side. This is the only nontrivial graph constructor supplied by the core library. + +.. cpp:function:: void graph::add_ghost() + + Calling this function on a graph object will add a ghost site to the graph. Explicitly, a new vertex is added that is adjacent to every other vertex in the graph. This vertex will have the last index, which is equal to number of vertices in the original graph. + diff --git a/doc/index.rst b/doc/index.rst new file mode 100644 index 0000000..4e58e46 --- /dev/null +++ b/doc/index.rst @@ -0,0 +1,15 @@ + +Wolff Libraries ++++++++++++++++ + +.. toctree :: + :maxdepth: 2 + + models + system + measurement + graph + types + finite_states + compile + diff --git a/doc/measurement.rst b/doc/measurement.rst new file mode 100644 index 0000000..e71010d --- /dev/null +++ b/doc/measurement.rst @@ -0,0 +1,84 @@ + +************ +Measurements +************ + +One of the most complicated parts of using this library is setting up measurements to take during a run, but once understood they provide a powerful way to quickly measure arbitrary properties. The class is defined in the header :file:`wolff/measurement.hpp`. + +.. cpp:class:: template\ measurement + + This class defines several virtual member functions. To use the library, you must supply your own measurement class that inherits this one and defines those functions, which may be trivial. + + Each member function defines a hook that is run at a specific time during execution. These hooks can be used to modify member objects of your inheritor measurement class, and thereby extract information from the simulation. + +.. cpp:function:: template& S, v_t i0, const R_t& r) + + This hook is run prior to cluster formation. The arguments supplied are: + + =========================== === + :cpp:var:`n` the number of clusters already flipped + :cpp:var:`N` the total number of clusters to flip + :cpp:var:`S` the current system state + :cpp:var:`i0` the index of the seed vertex for the cluster + :cpp:var:`r` the transformation by which the cluster will be flipped + =========================== === + +.. cpp:function:: void measurement::plain_bond_visited(const system& S, v_t i1, const X_t& s1_new, v_t i2, double dE) + + This hook is run when an ordinary edge is visited during cluster formation, whether it is activated or not. The arguments supplied are: + + ================================= === + :code:`const system& S` the current system state + :code:`v_t i1` the index of the vertex soon to be flipped + :code:`const X_t& s1_new` the state vertex :cpp:var:`i1` will be flipped to + :code:`v_t i2` the other vertex sharing the edge + :code:`double dE` the change in energy that will result when the spin at site :cpp:var:`i1` is flipped + ================================= === + +.. cpp:function:: void measurement::plain_site_transformed(const system& S, v_t i, const X_t& si_new) + + This hook is run when an ordinary site is about to be flipped. The arguments supplied are: + + ================================= === + :code:`const system& S` the current system state + :code:`v_t i` the index of the vertex soon to be flipped + :code:`const X_t& si_new` the state vertex :cpp:var:`i` will be flipped to + ================================= === + +| + +.. cpp:function:: void measurement::ghost_bond_visited(const system& S, v_t i, const X_t& s0si_old, const X_t& s0si_new, double dE) + + This hook is run when an edge containing the ghost site is visited during cluster formation, whether activated or not. The arguments supplied are: + + ================================= === + :code:`const system& S` the current system state + :code:`v_t i` the index of the non-ghost site in this edge + :code:`const X_t& s0si_old` the state that the spin on the non-ghost site has before the transformation is applied, rotated by the inverse action of the ghost site + :code:`const X_t& s0si_new` the state that the spin on the non-ghost site will have after the transformation is applied, rotated by the inverse action of the ghost site + :code:`double dE` the change in energy that will result when one site is transformed + ================================= === + +| + +.. cpp:function:: void measurement::ghost_site_transformed(const system& S, const R_t& R_new) + + This hook is run when the ghost site is about to be flipped. The arguments supplied are: + + ================================= === + :code:`const system& S` the current system state + :code:`const R_t& R_new` the state the ghost site will be flipped to + ================================= === + +| + +.. cpp:function:: void measurement::post_cluster(N_t n, N_t N, const system& S) + + This hook is run after a cluster has been flipped. The arguments supplied are: + + ================================= === + :code:`N_t n` the number of clusters already flipped + :code:`N_t N` the total number of clusters to flip + :code:`const system& S` the current system state + ================================= === + diff --git a/doc/models.rst b/doc/models.rst new file mode 100644 index 0000000..aa80ee6 --- /dev/null +++ b/doc/models.rst @@ -0,0 +1,61 @@ + +**************** +Defining a Model +**************** + +To define your own model, you need a class for spin states, a class for spin transformations, a bond and site coupling, and a generator of transformations. + +Spin States +=========== + +.. cpp:class:: X_t + + The spin type, which we have been denoting with :cpp:class:`X_t`, only needs to have a default constructor. If your spins can take only finitely many values, consider following the instructions in :ref:`finite_states` to significantly speed the algorithm. + +Transformations +=============== + +.. cpp:class:: R_t + + The transformation type must have a default constructor, and the following member functions: + + .. cpp:function:: X_t act(const X_t&) + + The action of the transformation on a spin, returning the transformed spin. + + .. cpp:function:: R_t act(const R_t&) + + The action of the transformation on another transformation, returning the transformed transformation. + + .. cpp:function:: X_t act_inverse(const X_t&) + + The action of the inverse transformation on a spin, returning the transformed spin. + + .. cpp:function:: R_t act_inverse(const R_t&) + + The action of the inverse transformation on another transformation, returning the transformed transformation. + +Couplings +========= + +When a :cpp:class:`system` object is initialized it needs to be given a bond and field coupling, to resemble the Hamiltonian + +.. math:: + \mathcal H = -\sum_{\langle ij\rangle}Z(s_i,s_j)-\sum_iB(s_i) + +Note that building with certain compile flags will change the form that these coupling functions must take, as outlined in :ref:`compile`. + +Bond Coupling +------------- + +.. cpp:function:: double Z(const X_t&, const X_t&) + + A function that takes two spins and outputs the coupling between them. For traditional spin models, this is typically something like a dot product. + +Field Coupling +-------------- + +.. cpp:function:: double B(const X_t&) + + A function that takes one spin and outputs the coupling between it and an external field. + diff --git a/doc/system.rst b/doc/system.rst new file mode 100644 index 0000000..1751e2b --- /dev/null +++ b/doc/system.rst @@ -0,0 +1,61 @@ + +************************************* +Constructing a System & Running Wolff +************************************* + +The core of the library lies in the :cpp:class:`system` class and its member functions. Here, the state of your model is stored and cluster flip Monte Carlo can be carried out in various ways. + +Note that the member objects and functions described here will change when compiled with certain compiler flags active, as described in :ref:`compile`. + +.. cpp:class:: template\ system + +Member Objects +============== + +.. cpp:member:: v_t system::nv + + Stores the number of ordinary sites in the model. + +.. cpp:member:: v_t system::ne + + Stores the number of ordinary bonds in the model. + +.. cpp:member:: graph system::G + + Stores the graph describing the lattice, including the ghost site. + +.. cpp:member:: double system::T + + Stores the temperature of the model. + +.. cpp:member:: std::vector system::s + + The :math:`i\text{th}` component stores the spin state of site :math:`i`. + +.. cpp:member:: R_t system::s0 + + Stores the state of the ghost site. + +.. cpp:member:: std::function system::Z + + The function that returns the coupling between neighboring sites. + +.. cpp:member:: std::function system::B + + The function that returns the coupling to the field. + +Member Functions +================ + +.. cpp:function:: system::system(graph G, double T, std::function Z, std::function B) + + The constructor for systems. The arguments given are a graph :cpp:any:`G` *without* the ghost spin added, the temperature :cpp:any:`T`, and the coupling functions :cpp:any:`Z` and :cpp:any:`B`. The states of the spins and ghost site are initialized using the default constructors for :cpp:type:`X_t` and :cpp:type:`R_t`, respectively. :cpp:any:`nv` and :cpp:any:`ne` are taken directly from :cpp:any:`G`, after which the ghost site is added to :cpp:any:`G`. + +.. cpp:function:: system::flip_cluster(v_t i0, const R_t& r, std::mt19937& rng, measurement& A) + + Performs one Wolff cluster flip to the system. The cluster is seeded at vertex :cpp:any:`i0` and the spins added are transformed by :cpp:any:`r`. A random number generator :cpp:any:`rng` provides required random numbers during, and the relevant measurement hooks defined in the inherited class of :cpp:any:`A` are run during the cluster formation. + +.. cpp:function:: system::run_wolff(N_t N, std::function &, v_t)> r_gen, measurement& A, std::mt19937& rng) + + Performs :cpp:any:`N` Wolff cluster flips to the system. One must provide a function :cpp:func:`r_gen` that takes a random number generator, the system state, and the index of the seed site and returns a transformation for each flip. Also required are an object from a class which inherits :cpp:class:`measurement` to run measurements, and a random number generator :cpp:any:`rng`. + diff --git a/doc/types.rst b/doc/types.rst new file mode 100644 index 0000000..b729ac3 --- /dev/null +++ b/doc/types.rst @@ -0,0 +1,36 @@ + +************ +Custom Types +************ + +The Wolff library uses several custom types, defined to promote memory and cache effiency, promote shorter lines, and associate natural scope to certain parameters. They are all aliases for `standard C99 fixed width integer types`_. These are defined in the header file :file:`wolff/types.h`. Here is a list of the types and what they are used to hold: + +.. cpp:type:: uint_fast32_t v_t + + Holds indicies for lattice vertices and edges. + +.. cpp:type:: uint_fast8_t q_t + + Holds indicies for enumerating states, as in :math:`q`-state Potts. + +.. cpp:type:: uint_fast8_t D_t + + Holds the dimension of space. + +.. cpp:type:: uint_fast16_t L_t + + Holds the linear extent of the lattice. + +.. cpp:type:: uint_fast64_t N_t + + Holds a count of Monte Carlo steps. + +Other Definitions +================= + +All types are unsigned integers of a minimum size appropriate for most realistic purposes, with the :c:type:`fast` directives meaning that your compiler may choose a larger type for efficiency's sake. + +For each type :c:type:`x_t`, the macro :c:macro:`MAX_x` supplies the maximum value the type may hold, the macro :c:macro:`PRIx` supplies the format constants for use with the :c:func:`fprintf` family of functions, and the macro :c:macro:`SCNx` supplies the format constants for use with the :c:func:`fscanf` family of functions. + +.. _standard C99 fixed width integer types: https://en.cppreference.com/w/c/types/integer + diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 73922c6..a4a4a5c 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,11 +1,11 @@ add_executable(ising ising.cpp) -add_executable(animate_ising animate_ising.cpp) +add_executable(ising_animation ising_animation.cpp) add_executable(xy On.cpp) add_compile_definitions(xy WOLFF_N=2) target_link_libraries(ising wolff) -target_link_libraries(animate_ising wolff GL GLU glut) +target_link_libraries(ising_animation wolff GL GLU glut) target_link_libraries(xy wolff) diff --git a/examples/On.cpp b/examples/On.cpp index e045f52..6c3f0fd 100644 --- a/examples/On.cpp +++ b/examples/On.cpp @@ -7,6 +7,7 @@ #include #include + #include int main(int argc, char *argv[]) { @@ -56,10 +57,13 @@ int main(int argc, char *argv[]) { return H * s; }; + // initialize the lattice + graph G(D, L); + // initialize the system - wolff_system, vector_t> S(D, L, T, Z, B); + system, vector_t> S(G, T, Z, B); - std::function (std::mt19937&, const vector_t&)> gen_R = generate_rotation_uniform; + std::function (std::mt19937&, const system, vector_t>&, v_t)> gen_R = generate_rotation_uniform; // initailze the measurement object simple_measurement A(S); @@ -69,7 +73,7 @@ int main(int argc, char *argv[]) { std::mt19937 rng{seed}; // run wolff N times - wolff, vector_t>(N, S, gen_R, A, rng); + S.run_wolff(N, gen_R, A, rng); // print the result of our measurements std::cout << "Wolff complete!\nThe average energy per site was " << A.avgE() / S.nv diff --git a/examples/animate_ising.cpp b/examples/animate_ising.cpp deleted file mode 100644 index e33736e..0000000 --- a/examples/animate_ising.cpp +++ /dev/null @@ -1,133 +0,0 @@ - -#include -#include -#include - -#include - -#include -#include -#include - -class draw_ising : public wolff_measurement { - private: - unsigned int frame_skip; - v_t C; - public: - draw_ising(const wolff_system& S, unsigned int window_size, unsigned int frame_skip, int argc, char *argv[]) : frame_skip(frame_skip){ - 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, S.L, 0.0, S.L); - } - - void pre_cluster(N_t, N_t, const wolff_system& S, v_t, const ising_t&) { - glClear(GL_COLOR_BUFFER_BIT); - for (v_t i = 0; i < pow(S.L, 2); i++) { - if (S.s[i].x == S.s0.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(); - C = 0; - } - - void plain_bond_visited(const wolff_system&, v_t, const ising_t&, v_t, double dE) {} - - void ghost_bond_visited(const wolff_system&, v_t, const ising_t& s_old, const ising_t& s_new, double dE) {} - - void plain_site_transformed(const wolff_system& S, v_t i, const ising_t&) { - glColor3f(1.0, 0.0, 0.0); - glRecti(i / S.L, i % S.L, (i / S.L) + 1, (i % S.L) + 1); - C++; - if (C % frame_skip == 0) { - glFlush(); - } - } - - void ghost_site_transformed(const wolff_system&, const ising_t&) {} - - void post_cluster(N_t, N_t, const wolff_system&) {} -}; - -int main(int argc, char *argv[]) { - - // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double H = 0.0; - unsigned int window_size = 512; - unsigned int frame_skip = 1; - - int opt; - - // take command line arguments - while ((opt = getopt(argc, argv, "N:D:L:T:H:w:f:")) != -1) { - switch (opt) { - case 'N': // number of steps - N = (N_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 'w': - window_size = atoi(optarg); - break; - case 'f': - frame_skip = atoi(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - // define the spin-spin coupling - std::function Z = [] (const ising_t& s1, const ising_t& s2) -> double { - return (double)(s1 * s2); - }; - - // define the spin-field coupling - std::function B = [=] (const ising_t& s) -> double { - return H * s; - }; - - // initialize the system - wolff_system S(D, L, T, Z, B); - - // define function that generates self-inverse rotations - std::function gen_R = [] (std::mt19937&, const ising_t& s) -> ising_t { - return ising_t(true); - }; - - // initailze the measurement object - draw_ising A(S, window_size, frame_skip, argc, argv); - - // initialize the random number generator - auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count(); - std::mt19937 rng{seed}; - - // run wolff N times - wolff(N, S, gen_R, A, rng); - - // exit - return 0; -} - diff --git a/examples/ising.cpp b/examples/ising.cpp index ffcb7e4..fe21279 100644 --- a/examples/ising.cpp +++ b/examples/ising.cpp @@ -7,8 +7,11 @@ #include #include + #include +using namespace wolff; + int main(int argc, char *argv[]) { // set defaults @@ -53,23 +56,24 @@ int main(int argc, char *argv[]) { return H * s; }; + // initialize the lattice + graph G(D, L); + // initialize the system - wolff_system S(D, L, T, Z, B); + system S(G, T, Z, B); // initialize the random number generator auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count(); std::mt19937 rng{seed}; // define function that generates self-inverse rotations - std::function gen_R = [] (std::mt19937&, const ising_t& s) -> ising_t { - return ising_t(true); - }; + std::function &, v_t)> gen_r = gen_ising; // initailze the measurement object simple_measurement A(S); // run wolff N times - wolff(N, S, gen_R, A, rng); + S.run_wolff(N, gen_r, A, rng); // print the result of our measurements std::cout << "Wolff complete!\nThe average energy per site was " << A.avgE() / S.nv diff --git a/examples/ising_animation.cpp b/examples/ising_animation.cpp new file mode 100644 index 0000000..2104e91 --- /dev/null +++ b/examples/ising_animation.cpp @@ -0,0 +1,139 @@ + +#include +#include +#include + +#include + +#include +#include + +#include + +using namespace wolff; + +class draw_ising : public measurement { + private: + unsigned int frame_skip; + v_t C; + public: + draw_ising(const system& S, unsigned int window_size, unsigned int frame_skip, int argc, char *argv[]) : frame_skip(frame_skip){ + 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, S.G.L, 0.0, S.G.L); + } + + void pre_cluster(N_t, N_t, const system& S, v_t, const ising_t&) { + glClear(GL_COLOR_BUFFER_BIT); + for (v_t i = 0; i < pow(S.G.L, 2); i++) { + if (S.s[i].x == S.s0.x) { + glColor3f(0.0, 0.0, 0.0); + } else { + glColor3f(1.0, 1.0, 1.0); + } + glRecti(i / S.G.L, i % S.G.L, (i / S.G.L) + 1, (i % S.G.L) + 1); + } + glFlush(); + C = 0; + } + + void plain_bond_visited(const system&, v_t, const ising_t&, v_t, double dE) {} + + void ghost_bond_visited(const system&, v_t, const ising_t& s_old, const ising_t& s_new, double dE) {} + + void plain_site_transformed(const system& S, v_t i, const ising_t&) { + glColor3f(1.0, 0.0, 0.0); + glRecti(i / S.G.L, i % S.G.L, (i / S.G.L) + 1, (i % S.G.L) + 1); + C++; + if (C % frame_skip == 0) { + glFlush(); + } + } + + void ghost_site_transformed(const system&, const ising_t&) {} + + void post_cluster(N_t, N_t, const system&) {} +}; + +int main(int argc, char *argv[]) { + + // set defaults + N_t N = (N_t)1e4; + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double H = 0.0; + unsigned int window_size = 512; + unsigned int frame_skip = 1; + + int opt; + + // take command line arguments + while ((opt = getopt(argc, argv, "N:D:L:T:H:w:f:")) != -1) { + switch (opt) { + case 'N': // number of steps + N = (N_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 'w': + window_size = atoi(optarg); + break; + case 'f': + frame_skip = atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + // define the spin-spin coupling + std::function Z = [] (const ising_t& s1, const ising_t& s2) -> double { + return (double)(s1 * s2); + }; + + // define the spin-field coupling + std::function B = [=] (const ising_t& s) -> double { + return H * s; + }; + + // initialize the lattice + graph G(D, L); + + // initialize the system + system S(G, T, Z, B); + + // define function that generates self-inverse rotations + std::function &, v_t)> gen_R = [] (std::mt19937&, const system&, v_t) -> ising_t { + return ising_t(true); + }; + + // initailze the measurement object + draw_ising A(S, window_size, frame_skip, argc, argv); + + // initialize the random number generator + auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count(); + std::mt19937 rng{seed}; + + // run wolff N times + S.run_wolff(N, gen_R, A, rng); + + // exit + return 0; +} + diff --git a/examples/simple_measurement.hpp b/examples/simple_measurement.hpp index 2287c58..518631c 100644 --- a/examples/simple_measurement.hpp +++ b/examples/simple_measurement.hpp @@ -1,8 +1,10 @@ #include +using namespace wolff; + template -class simple_measurement : public wolff_measurement { +class simple_measurement : public measurement { private: N_t n; @@ -15,7 +17,7 @@ class simple_measurement : public wolff_measurement { double totalC; public: - simple_measurement(const wolff_system& S) { + simple_measurement(const system& S) { n = 0; M = S.nv * S.s[0]; E = - (S.ne * S.Z(S.s[0], S.s[0]) + S.nv * S.B(S.s[0])); @@ -25,26 +27,26 @@ class simple_measurement : public wolff_measurement { totalC = 0; } - void pre_cluster(N_t, N_t, const wolff_system&, v_t, const R_t&) { + void pre_cluster(N_t, N_t, const system&, v_t, const R_t&) { C = 0; } - void plain_bond_visited(const wolff_system&, v_t, const X_t&, v_t, double dE) { + void plain_bond_visited(const system&, v_t, const X_t&, v_t, double dE) { E += dE; } - void ghost_bond_visited(const wolff_system&, v_t, const X_t& s_old, const X_t& s_new, double dE) { + void ghost_bond_visited(const system&, v_t, const X_t& s_old, const X_t& s_new, double dE) { E += dE; M += s_new - s_old; } - void plain_site_transformed(const wolff_system&, v_t, const X_t&) { + void plain_site_transformed(const system&, v_t, const X_t&) { C++; } - void ghost_site_transformed(const wolff_system&, const R_t&) {} + void ghost_site_transformed(const system&, const R_t&) {} - void post_cluster(N_t, N_t, const wolff_system&) { + void post_cluster(N_t, N_t, const system&) { totalE += E; totalM += M; totalC += C; diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp index b78cd4b..7e909c8 100644 --- a/lib/include/wolff.hpp +++ b/lib/include/wolff.hpp @@ -3,25 +3,30 @@ #define WOLFF_H #include "wolff/cluster.hpp" +#include "wolff/measurement.hpp" + +namespace wolff{ template -void wolff(N_t N, wolff_system& S, - std::function r_gen, - wolff_measurement& A, std::mt19937& rng) { +void system::run_wolff(N_t N, + std::function &, v_t)> r_gen, + measurement& A, std::mt19937& rng) { - std::uniform_int_distribution dist(0, S.nv - 1); + std::uniform_int_distribution dist(0, nv - 1); for (N_t n = 0; n < N; n++) { v_t i0 = dist(rng); - R_t r = r_gen(rng, S.s[i0]); + R_t r = r_gen(rng, *this, i0); - A.pre_cluster(n, N, S, i0, r); + A.pre_cluster(n, N, *this, i0, r); - wolff_cluster_flip(S, i0, r, rng, A); + this->flip_cluster(i0, r, rng, A); - A.post_cluster(n, N, S); + A.post_cluster(n, N, *this); } } +} + #endif diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp index 055cdf3..b66f367 100644 --- a/lib/include/wolff/cluster.hpp +++ b/lib/include/wolff/cluster.hpp @@ -7,20 +7,19 @@ #include #include -#include "types.h" #include "system.hpp" -#include "graph.hpp" -#include "measurement.hpp" + +namespace wolff { template -void wolff_cluster_flip(wolff_system& S, v_t i0, const R_t& r, - std::mt19937& rng, wolff_measurement& A) { +void system::flip_cluster(v_t i0, const R_t& r, + std::mt19937& rng, measurement& A) { std::uniform_real_distribution dist(0.0, 1.0); std::queue queue; queue.push(i0); - std::vector marks(S.G.nv, false); + std::vector marks(G.nv, false); while (!queue.empty()) { v_t i = queue.front(); @@ -33,21 +32,21 @@ void wolff_cluster_flip(wolff_system& S, v_t i0, const R_t& r, #ifndef WOLFF_NO_FIELD R_t s0_new; - bool we_are_ghost = (i == S.nv); + bool we_are_ghost = (i == nv); if (we_are_ghost) { - s0_new = r.act(S.s0); + s0_new = r.act(s0); } else #endif { - si_new = r.act(S.s[i]); + si_new = r.act(s[i]); } - for (const v_t &j : S.G.adj[i]) { + for (const v_t &j : G.adjacency[i]) { double dE, p; #ifndef WOLFF_NO_FIELD - bool neighbor_is_ghost = (j == S.nv); + bool neighbor_is_ghost = (j == nv); if (we_are_ghost || neighbor_is_ghost) { X_t s0s_old, s0s_new; @@ -55,18 +54,18 @@ void wolff_cluster_flip(wolff_system& S, v_t i0, const R_t& r, if (neighbor_is_ghost) { non_ghost = i; - s0s_old = S.s0.act_inverse(S.s[i]); - s0s_new = S.s0.act_inverse(si_new); + s0s_old = s0.act_inverse(s[i]); + s0s_new = s0.act_inverse(si_new); } else { non_ghost = j; - s0s_old = S.s0.act_inverse(S.s[j]); - s0s_new = s0_new.act_inverse(S.s[j]); + s0s_old = s0.act_inverse(s[j]); + s0s_new = s0_new.act_inverse(s[j]); } #ifdef WOLFF_SITE_DEPENDENCE - dE = S.B(non_ghost, s0s_old) - S.B(non_ghost, s0s_new); + dE = B(non_ghost, s0s_old) - B(non_ghost, s0s_new); #else - dE = S.B(s0s_old) - S.B(s0s_new); + dE = B(s0s_old) - B(s0s_new); #endif #ifdef WOLFF_FINITE_STATES @@ -75,28 +74,28 @@ void wolff_cluster_flip(wolff_system& S, v_t i0, const R_t& r, #endif // run measurement hooks for encountering a ghost bond - A.ghost_bond_visited(S, non_ghost, s0s_old, s0s_new, dE); + A.ghost_bond_visited(*this, non_ghost, s0s_old, s0s_new, dE); } else // this is a perfectly normal bond! #endif { #ifdef WOLFF_BOND_DEPENDENCE - dE = S.Z(i, S.s[i], j, S.s[j]) - S.Z(i, si_new, j, S.s[j]); + dE = Z(i, s[i], j, s[j]) - Z(i, si_new, j, s[j]); #else - dE = S.Z(S.s[i], S.s[j]) - S.Z(si_new, S.s[j]); + dE = Z(s[i], s[j]) - Z(si_new, s[j]); #endif #ifdef WOLFF_FINITE_STATES - p = finite_states_Zp[finite_states_enum(S.s[i])] + p = finite_states_Zp[finite_states_enum(s[i])] [finite_states_enum(si_new)] - [finite_states_enum(S.s[j])]; + [finite_states_enum(s[j])]; #endif // run measurement hooks for encountering a plain bond - A.plain_bond_visited(S, i, si_new, j, dE); + A.plain_bond_visited(*this, i, si_new, j, dE); } #ifndef FINITE_STATES - p = 1.0 - exp(-dE / S.T); + p = 1.0 - exp(-dE / T); #endif if (dist(rng) < p) { @@ -106,17 +105,19 @@ void wolff_cluster_flip(wolff_system& S, v_t i0, const R_t& r, #ifndef WOLFF_NO_FIELD if (we_are_ghost) { - A.ghost_site_transformed(S, s0_new); - S.s0 = s0_new; + A.ghost_site_transformed(*this, s0_new); + s0 = s0_new; } else #endif { - A.plain_site_transformed(S, i, si_new); - S.s[i] = si_new; + A.plain_site_transformed(*this, i, si_new); + s[i] = si_new; } } } } +} + #endif diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp index 4cca69e..7e54eff 100644 --- a/lib/include/wolff/finite_states.hpp +++ b/lib/include/wolff/finite_states.hpp @@ -8,13 +8,15 @@ #include "system.hpp" +namespace wolff { + std::array, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Zp; #ifndef WOLFF_NO_FIELD std::array, WOLFF_FINITE_STATES_N> finite_states_Bp; #endif template -void finite_states_init(const wolff_system& S) { +void finite_states_init(const system& S) { #ifndef WOLFF_NO_FIELD for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) { for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) { @@ -31,5 +33,7 @@ void finite_states_init(const wolff_system& S) { } } +} + #endif diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp index fcd73fa..65b3941 100644 --- a/lib/include/wolff/graph.hpp +++ b/lib/include/wolff/graph.hpp @@ -5,23 +5,26 @@ #include #include -#include "types.h" +namespace wolff { -typedef enum lattice_t { - SQUARE_LATTICE, - DIAGONAL_LATTICE -} lattice_t; +#include "types.h" -class graph_t { +class graph { public: + D_t D; + L_t L; v_t ne; v_t nv; - std::vector> adj; + std::vector> adjacency; std::vector> coordinate; - graph_t(D_t D, L_t L, lattice_t lat = SQUARE_LATTICE); - void add_ext(); + graph(); + graph(D_t D, L_t L); + + void add_ghost(); }; +} + #endif diff --git a/lib/include/wolff/measurement.hpp b/lib/include/wolff/measurement.hpp index a6a5f79..604eaf5 100644 --- a/lib/include/wolff/measurement.hpp +++ b/lib/include/wolff/measurement.hpp @@ -4,21 +4,25 @@ #include "system.hpp" +namespace wolff { + template -class wolff_measurement { +class measurement { public: - virtual void pre_cluster(N_t, N_t, const wolff_system&, v_t, const R_t&) = 0; + virtual void pre_cluster(N_t, N_t, const system&, v_t, const R_t&) = 0; - virtual void plain_bond_visited(const wolff_system&, v_t, const X_t&, v_t, double) = 0; - virtual void plain_site_transformed(const wolff_system&, v_t, const X_t&) = 0; + virtual void plain_bond_visited(const system&, v_t, const X_t&, v_t, double) = 0; + virtual void plain_site_transformed(const system&, v_t, const X_t&) = 0; #ifndef WOLFF_NO_FIELD - virtual void ghost_bond_visited(const wolff_system&, v_t, const X_t&, const X_t&, double) = 0; - virtual void ghost_site_transformed(const wolff_system&, const R_t&) = 0; + virtual void ghost_bond_visited(const system&, v_t, const X_t&, const X_t&, double) = 0; + virtual void ghost_site_transformed(const system&, const R_t&) = 0; #endif - virtual void post_cluster(N_t, N_t, const wolff_system&) = 0; + virtual void post_cluster(N_t, N_t, const system&) = 0; }; +} + #endif diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp index c6a3a47..824c063 100644 --- a/lib/include/wolff/models/ising.hpp +++ b/lib/include/wolff/models/ising.hpp @@ -3,6 +3,9 @@ #define WOLFF_MODELS_ISING #include "../types.h" +#include "../system.hpp" + +namespace wolff { class ising_t { public: @@ -73,9 +76,15 @@ inline ising_t::F_t operator*(double a, const ising_t& s) { return s * a; } +ising_t gen_ising(std::mt19937&, const system&, v_t) { + return ising_t(true); +}; + #define WOLFF_FINITE_STATES_N 2 const ising_t finite_states_possible[2] = {ising_t(0), ising_t(1)}; -q_t finite_states_enum(ising_t state) { return (q_t)state.x; } +q_t finite_states_enum(const ising_t& state) { return (q_t)state.x; } + +} #endif diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp index 9dd5ddd..58203b7 100644 --- a/lib/include/wolff/models/orthogonal.hpp +++ b/lib/include/wolff/models/orthogonal.hpp @@ -4,7 +4,8 @@ #include #include -#include +#include "../types.h" +#include "../system.hpp" #include "vector.hpp" template @@ -113,7 +114,7 @@ class orthogonal_t : public std::array, q> { }; template -orthogonal_t generate_rotation_uniform (std::mt19937& r, const vector_t & v) { +orthogonal_t generate_rotation_uniform (std::mt19937& r, const system, vector_t>&, v_t) { std::normal_distribution dist(0.0,1.0); orthogonal_t ptr; ptr.is_reflection = true; @@ -135,7 +136,7 @@ orthogonal_t generate_rotation_uniform (std::mt19937& r, const vecto } template -orthogonal_t generate_rotation_perturbation (std::mt19937& r, const vector_t & v0, double epsilon, unsigned int n) { +orthogonal_t generate_rotation_perturbation (std::mt19937& r, const system, vector_t>& S, v_t i0, double epsilon, unsigned int n) { std::normal_distribution dist(0.0,1.0); orthogonal_t m; m.is_reflection = true; @@ -149,14 +150,14 @@ orthogonal_t generate_rotation_perturbation (std::mt19937& r, const 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; + v[0] = S.s[i0][0] * cosr - S.s[i0][1] * sinr; + v[1] = S.s[i0][1] * cosr + S.s[i0][0] * sinr; for (q_t i = 2; i < q; i++) { - v[i] = v0[i]; + v[i] = S.s[i0][i]; } } else { - v = v0; + v = S.s[i0]; } double m_dot_v = 0; diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp index 903f25f..10b440e 100644 --- a/lib/include/wolff/models/potts.hpp +++ b/lib/include/wolff/models/potts.hpp @@ -46,5 +46,5 @@ class potts_t { 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 finite_states_enum(potts_t state) { return (q_t)state.x; } +q_t finite_states_enum(const potts_t& state) { return (q_t)state.x; } diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp index ac82f44..32ad38e 100644 --- a/lib/include/wolff/system.hpp +++ b/lib/include/wolff/system.hpp @@ -4,23 +4,25 @@ #include #include +#include -#include "types.h" #include "graph.hpp" +namespace wolff { + +#include "types.h" + +template +class measurement; + template -class wolff_system { +class system { public: - D_t D; // dimension - L_t L; // linear size v_t nv; // number of vertices v_t ne; // number of edges - graph_t G; // the graph defining the lattice with ghost + graph G; // the graph defining the lattice with ghost double T; // the temperature std::vector s; // the state of the ordinary spins -#ifndef WOLFF_NO_FIELD - R_t s0; // the current state of the ghost site -#endif #ifdef WOLFF_BOND_DEPENDENCE std::function Z; // coupling between sites @@ -29,6 +31,7 @@ class wolff_system { #endif #ifndef WOLFF_NO_FIELD + R_t s0; // the current state of the ghost site #ifdef WOLFF_SITE_DEPENDENCE std::function B; // coupling with the external field #else @@ -36,7 +39,7 @@ class wolff_system { #endif #endif - wolff_system(D_t D, L_t L, double T, + system(graph g, double T, #ifdef WOLFF_BOND_DEPENDENCE std::function Z #else @@ -49,7 +52,7 @@ class wolff_system { , std::function B #endif #endif - , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), G(D, L, lat), T(T), Z(Z) + ) : G(g), T(T), Z(Z) #ifndef WOLFF_NO_FIELD , s0(), B(B) #endif @@ -58,13 +61,18 @@ class wolff_system { ne = G.ne; s.resize(nv); #ifndef WOLFF_NO_FIELD - G.add_ext(); + G.add_ghost(); #endif #ifdef WOLFF_FINITE_STATES finite_states_init(*this); #endif } + + void flip_cluster(v_t, const R_t&, std::mt19937&, measurement&); + void run_wolff(N_t, std::function &, v_t)> r_gen, measurement& A, std::mt19937& rng); }; +} + #endif diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h index 68bd09e..1bc8c4d 100644 --- a/lib/include/wolff/types.h +++ b/lib/include/wolff/types.h @@ -8,7 +8,7 @@ typedef uint_fast16_t L_t; // linear size typedef uint_fast64_t N_t; // cycle iterator #define MAX_v UINT_FAST32_MAX -#define MAX_Q UINT_FAST8_MAX +#define MAX_q UINT_FAST8_MAX #define MAX_D UINT_FAST8_MAX #define MAX_L UINT_FAST16_MAX #define MAX_N UINT_FAST64_MAX diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index b42aa78..4f91be4 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -1,74 +1,55 @@ #include -graph_t::graph_t(D_t D, L_t L, lattice_t lat) { - switch (lat) { - case SQUARE_LATTICE: { - nv = pow(L, D); - ne = D * nv; +namespace wolff { - adj.resize(nv); - coordinate.resize(nv); - - for (std::vector adj_i : adj) { - 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; - - adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1))); - adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1))); - } - } - break; - } - case DIAGONAL_LATTICE: { - nv = D * pow(L, D); - ne = 2 * nv; +graph::graph() { + D = 0; + L = 0; + nv = 0; + ne = 0; +} - adj.resize(nv); - coordinate.resize(nv); +graph::graph(D_t D, L_t L) : D(D), L(L) { + nv = pow(L, D); + ne = D * nv; - for (std::vector adj_i : adj) { - adj_i.reserve(4 * (D - 1)); - } + adjacency.resize(nv); + coordinate.resize(nv); - for (D_t i = 0; i < D; i++) { - v_t sb = i * pow(L, D); + for (std::vector adj_i : adjacency) { + adj_i.reserve(2 * D); + } - for (v_t j = 0; j < pow(L, D); j++) { - v_t vc = sb + j; + 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; - adj[vc].push_back(((i + 1) % D) * pow(L, D) + j); - adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L); - adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L); - adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L); - } - } - break; - } + adjacency[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1))); + adjacency[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& adj_i : adj) { +void graph::add_ghost() { + for (std::vector& adj_i : adjacency) { adj_i.push_back(nv); } - adj.resize(nv + 1); + adjacency.resize(nv + 1); coordinate.resize(nv + 1); - adj[nv].reserve(nv); + adjacency[nv].reserve(nv); for (v_t i = 0; i < nv; i++) { - adj[nv].push_back(i); + adjacency[nv].push_back(i); } coordinate[nv].resize(coordinate[0].size()); ne += nv; - nv += 1; + nv++; +} + } -- cgit v1.2.3-54-g00ecf