From 49ac78a6c04e215950bc9c0f97368605e63da15b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 14 Jan 2019 15:47:59 -0500 Subject: Large refactoring around changes in the graph class. - Graphs now use lists of references instead of vectors of indicies. - Vertices and edges have associated classes that can be given arbitrary properties via template specification. - All essential library headers have been combined into one, wolff.hpp. --- doc/compile.rst | 5 +- doc/finite_states.rst | 2 +- doc/graph.rst | 63 ++++-- doc/index.rst | 1 - doc/measurement.rst | 51 +++-- doc/models.rst | 2 +- doc/system.rst | 30 +-- doc/types.rst | 36 ---- examples/On.cpp | 22 +- examples/clock.cpp | 30 ++- examples/continuous_gaussian.cpp | 33 ++- examples/discrete_gaussian.cpp | 41 ++-- examples/ising.cpp | 18 +- examples/ising_animation.cpp | 84 ++++---- examples/ising_no_field.cpp | 20 +- examples/ising_random_field.cpp | 31 ++- examples/ising_standalone.cpp | 25 ++- examples/potts.cpp | 36 ++-- examples/simple_measurement.hpp | 36 ++-- lib/CMakeLists.txt | 16 +- lib/include/wolff.hpp | 32 --- lib/include/wolff/cluster.hpp | 120 ----------- lib/include/wolff/graph.hpp | 30 --- lib/include/wolff/measurement.hpp | 28 --- lib/include/wolff/models/dihedral.hpp | 56 ------ lib/include/wolff/models/dihedral_inf.hpp | 57 ------ lib/include/wolff/models/height.hpp | 57 ------ lib/include/wolff/models/ising.hpp | 93 --------- lib/include/wolff/models/orthogonal.hpp | 208 ------------------- lib/include/wolff/models/potts.hpp | 68 ------- lib/include/wolff/models/symmetric.hpp | 59 ------ lib/include/wolff/models/vector.hpp | 115 ----------- lib/include/wolff/system.hpp | 105 ---------- lib/include/wolff/types.h | 32 --- lib/src/graph.cpp | 55 ----- lib/wolff.hpp | 320 ++++++++++++++++++++++++++++++ lib/wolff_models/dihedral.hpp | 54 +++++ lib/wolff_models/dihedral_inf.hpp | 55 +++++ lib/wolff_models/height.hpp | 56 ++++++ lib/wolff_models/ising.hpp | 92 +++++++++ lib/wolff_models/orthogonal.hpp | 206 +++++++++++++++++++ lib/wolff_models/potts.hpp | 66 ++++++ lib/wolff_models/symmetric.hpp | 57 ++++++ lib/wolff_models/vector.hpp | 113 +++++++++++ 44 files changed, 1297 insertions(+), 1419 deletions(-) delete mode 100644 doc/types.rst delete mode 100644 lib/include/wolff.hpp delete mode 100644 lib/include/wolff/cluster.hpp delete mode 100644 lib/include/wolff/graph.hpp delete mode 100644 lib/include/wolff/measurement.hpp delete mode 100644 lib/include/wolff/models/dihedral.hpp delete mode 100644 lib/include/wolff/models/dihedral_inf.hpp delete mode 100644 lib/include/wolff/models/height.hpp delete mode 100644 lib/include/wolff/models/ising.hpp delete mode 100644 lib/include/wolff/models/orthogonal.hpp delete mode 100644 lib/include/wolff/models/potts.hpp delete mode 100644 lib/include/wolff/models/symmetric.hpp delete mode 100644 lib/include/wolff/models/vector.hpp delete mode 100644 lib/include/wolff/system.hpp delete mode 100644 lib/include/wolff/types.h delete mode 100644 lib/src/graph.cpp create mode 100644 lib/wolff.hpp create mode 100644 lib/wolff_models/dihedral.hpp create mode 100644 lib/wolff_models/dihedral_inf.hpp create mode 100644 lib/wolff_models/height.hpp create mode 100644 lib/wolff_models/ising.hpp create mode 100644 lib/wolff_models/orthogonal.hpp create mode 100644 lib/wolff_models/potts.hpp create mode 100644 lib/wolff_models/symmetric.hpp create mode 100644 lib/wolff_models/vector.hpp diff --git a/doc/compile.rst b/doc/compile.rst index e4960a3..d335dff 100644 --- a/doc/compile.rst +++ b/doc/compile.rst @@ -23,7 +23,7 @@ Building With 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:any:`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. + When defined, the bond coupling must be a function of the form :cpp:any:`double Z(const G_t::halfedge&, const X_t&, const X_t&)`, where the first argument is the edge being considered. A function of this type is passed to :cpp:class:`system` in place of the original bond coupling. Building With Site Dependence @@ -33,5 +33,6 @@ Building With 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:any:`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. + When defined, the field coupling must be a function of the form :cpp:any:`double B(const G_t::vertex&, const X_t&)`, where the first argument is the vertex the spin is on. A function of this type is passed to :cpp:class:`system` in place of the original field coupling. + An example of a system of this type can be found in :file:`examples/ising_random_field.cpp`, which uses a non-trivial vertex property to communicate vertex dependence to the field coupling function. diff --git a/doc/finite_states.rst b/doc/finite_states.rst index c36881e..66e34d2 100644 --- a/doc/finite_states.rst +++ b/doc/finite_states.rst @@ -5,7 +5,7 @@ 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 pair of spins. Once the appropriate things have been defined for your model, the compile definition :c:macro:`WOLFF_USE_FINITE_STATES` can be set to automate this process. The provided model headers :file:`wolff/models/ising.hpp` and :file:`wolff/models/potts.hpp` demonstrate the expected usage. +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 pair of spins. Once the appropriate things have been defined for your model, the compile definition :c:macro:`WOLFF_USE_FINITE_STATES` can be set 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 ==================== diff --git a/doc/graph.rst b/doc/graph.rst index 8a96e50..dfe684c 100644 --- a/doc/graph.rst +++ b/doc/graph.rst @@ -1,50 +1,83 @@ .. default-domain:: cpp +.. _graphs: ****** Graphs ****** -This class is defined in the header :file:`wolff/graph.hpp`. +This class is defined in the header :file:`lib/wolff.hpp`. -.. class:: graph +.. class:: \template , class edge_prop = std::tuple\<>> graph - Lattices are described by objects of class :class:`graph`, a minimal implementation of graphs. + Lattices are described by objects of class :class:`graph`, a minimal implementation of graphs. Can be called with `graph<>` if no properties need to be associated with vertices or edges. Otherwise, those properties can be supplied as classes. - .. member:: D_t graph::D + .. member:: unsigned graph::D The dimension of the graph. This property is unused by the core library. - .. member:: L_t graph::L + .. member:: unsigned graph::L The linear size of the graph. This property is unused by the core library. - .. member:: v_t graph::ne + .. member:: unsigned graph::ne The number of edges in the graph. This property is unused by the core library. - .. member:: v_t graph::nv + .. member:: unsigned graph::nv The number of vertices in the graph. - .. member:: std::vector> graph::adjacency + .. class:: graph::vertex - The adjacency list for the graph. The :math:`i\text{th}` element of :member:`adjacency` is a standard library vector containing the indices of all vertices adjacent to vertex :math:`i`. + This class describes the vertices on the graph. - .. member:: std::vector> graph::coordinate + .. member:: unsigned ind - The coordinates of the graph vertices. The :math:`i\text{th}` element of :var:`coordinate` is a standard library vector of length :var:`D` containing the spatial coordinates of vertex :math:`i`. This property is unused by the core library. + The index of the vertex, which is also its position in the list :var:`vertices`. + + .. member:: std::list edges + + The list of edges emanating from this vertex. + + .. member:: vertex_prop prop + + Template-defined class which stores optional properties of the vertex. + + .. class:: graph::halfedge + + This class describes the halfedges on the graph. + + .. member:: vertex& self + + A reference to the vertex this halfedge is emanating from. + + .. member:: vertex& neighbor + + A reference to the vertex this halfedge is going to. + + .. member:: edge_prop prop + + Template-defined class which stores optional properties of the edge. + + .. function:: halfedge(vertex &self, vertex &neighbor) + + Constructor which sets self and neighbor from supplied vertices. + + .. member:: std::vector graph::vertices + + A list of all vertices in the graph. .. function:: graph::graph() - The default constructor. Initializes an empty graph, i.e., :var:`D`, :var:`L`, :var:`ne`, and :var:`nv` are all zero and :var:`adjacency` and :var:`coordinate` are uninitialized. + The default constructor. Initializes an empty graph, i.e., :var:`D`, :var:`L`, :var:`ne`, and :var:`nv` are all zero and :var:`vertices` is uninitialized. - .. function:: graph::graph(D_t D, L_t L) + .. function:: graph::graph(unsigned D, unsigned L) Initializes a graph of a :var:`D`-dimensional hypercubic lattice with :var:`L` vertices per side. This is the only nontrivial graph constructor supplied by the core library. The library will work with arbitrary graphs, and if a different lattice is needed consider calling the default constructor and populating the member objects youself before handing the graph to the :class:`system` constructor. - :param D_t D: The dimension of space. - :param L_t L: The number of vertices per edge of the hypercube. + :param unsigned D: The dimension of space. + :param unsigned L: The number of vertices per edge of the hypercube. .. function:: void graph::add_ghost() diff --git a/doc/index.rst b/doc/index.rst index 7b156d3..c379dd7 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -11,7 +11,6 @@ Wolff Library system measurement graph - types finite_states compile diff --git a/doc/measurement.rst b/doc/measurement.rst index 1b4381c..1ca41cc 100644 --- a/doc/measurement.rst +++ b/doc/measurement.rst @@ -5,64 +5,63 @@ 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`. +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:`lib/wolff.hpp`. -.. class:: template measurement +.. 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. - .. function:: void measurement::pre_cluster(N_t n, N_t N, const system& S, v_t i0, const R_t& r) + .. function:: void measurement::pre_cluster(unsigned n, unsigned N, const system& S, const G_t::vertex& v, const R_t& r) This hook is run prior to cluster formation. - :param N_t n: The number of clusters already flipped. - :param N_t N: The total number of clusters to flip. - :param const system& S: The current system state. - :param v_t i0: The index of the seed vertex for the cluster. + :param unsigned n: The number of clusters already flipped. + :param unsigned N: The total number of clusters to flip. + :param const system& S: The current system state. + :param const G_t\:\:vertex& v: The seed vertex for the cluster. :param const R_t& r: The transformation by which the cluster will be flipped. - .. function:: void measurement::plain_bond_visited(const system& S, v_t i1, const X_t& s1_new, v_t i2, double dE) + .. function:: void measurement::plain_bond_visited(const system& S, const G_t::halfedge& e, const X_t& s1_new, double dE) This hook is run when an ordinary edge is visited during cluster formation, whether it is activated or not. - :param const system& S: The current system state. - :param v_t i1: The index of the vertex soon to be flipped. - :param const X_t& s1_new: The state vertex :code:`i1` will be flipped to. - :param v_t i2: The other vertex sharing the edge. - :param double dE: The change in energy that will result when the spin at site :code:`i1` is flipped. + :param const system& S: The current system state. + :param const G_t\:\:halfedge& e: The edge being considered, with e.self referencing the vertex soon to be flipped. + :param const X_t& s1_new: The state the vertex will be flipped to. + :param double dE: The change in energy that will result when the spin is flipped. - .. cpp:function:: void measurement::plain_site_transformed(const system& S, v_t i, const X_t& si_new) + .. cpp:function:: void measurement::plain_site_transformed(const system& S, const G_t::vertex& v, const X_t& si_new) This hook is run when an ordinary site is about to be flipped. - :param const system& S: The current system state. - :param v_t i: The index of the vertex soon to be flipped. - :param const X_t& si_new: The state vertex :code:`i` will be flipped to. + :param const system& S: The current system state. + :param const G_t\:\:vertex& v: The vertex soon to be flipped. + :param const X_t& si_new: The state that vertex 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) + .. cpp:function:: void measurement::ghost_bond_visited(const system& S, const G_t::vertex& v, 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. - :param const system& S: The current system state. - :param v_t i: The index of the non-ghost site in this edge. + :param const system& S: The current system state. + :param const G_t\:\:vertex& v: The non-ghost site in this edge. :param 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. :param 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. :param 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) + .. 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. - :param const system& S: The current system state. + :param const system& S: The current system state. :param 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) + .. cpp:function:: void measurement::post_cluster(unsigned n, unsigned N, const system& S) This hook is run after a cluster has been flipped. - :param N_t n: The number of clusters already flipped. - :param N_t N: The total number of clusters to flip. - :param const system& S: The current system state. + :param unsigned n: The number of clusters already flipped. + :param unsigned N: The total number of clusters to flip. + :param const system& S: The current system state. diff --git a/doc/models.rst b/doc/models.rst index e376e2e..8ba7e51 100644 --- a/doc/models.rst +++ b/doc/models.rst @@ -3,7 +3,7 @@ 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. Many examples of models can be found in the headers contained in :file:`wolff/models/`. +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. Many examples of models can be found in the headers contained in :file:`lib/wolff_models/`. Spin States =========== diff --git a/doc/system.rst b/doc/system.rst index 928d85b..9cd4a90 100644 --- a/doc/system.rst +++ b/doc/system.rst @@ -5,25 +5,25 @@ Constructing a System & Running Wolff ************************************* -This class and associated member functions are defined in the header files :file:`wolff.hpp`, :file:`wolff/system.hpp`, and :file:`wolff/cluster.hpp`. +This class and associated member functions are defined in the header file :file:`lib/wolff.hpp`. -.. class:: template\ system +.. class:: template\> system The core of the library lies in the :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`. - .. member:: v_t system::nv + .. member:: unsigned system::nv Stores the number of ordinary sites in the model. - .. member:: v_t system::ne + .. member:: unsigned system::ne Stores the number of ordinary bonds in the model. - .. member:: graph system::G + .. member:: G_t system::G - Stores the graph describing the lattice, including the ghost site. + Stores the graph describing the lattice, including the ghost site. See :ref:`graphs` for details on this class. The template parameter has a default option corresponding to a graph type with no special vertex or edge properties. .. member:: double system::T @@ -45,33 +45,33 @@ This class and associated member functions are defined in the header files :file The function that returns the coupling to the field. - .. function:: system::system(graph G, double T, std::function Z, std::function B) + .. function:: system::system(G_t G, double T, std::function Z, std::function B) The constructor for systems. - :param graph G: A lattice graph *without* the ghost spin added. + :param G_t G: A lattice graph *without* the ghost spin added. :param double T: The temperature. :param std\:\:function Z: The bond coupling. :param std\:\:function B: The field coupling. The states of the spins and ghost site are initialized using the default constructors for :type:`X_t` and :type:`R_t`, respectively. :any:`nv` and :any:`ne` are taken directly from :any:`G`, after which the ghost site is added to :any:`G`. - .. function:: system::flip_cluster(v_t i0, const R_t& r, std::mt19937& rng, measurement& A) + .. function:: system::flip_cluster(const G_t::vertex& v0, const R_t& r, std::mt19937& rng, measurement& A) Performs one Wolff cluster flip to the system. - :param v_t i0: The index of the seed site. + :param const G_t\:\:vertex& v0: The vertex of the seed site. :param const R_t& r: The transformation by which the cluster is flipped. :param std\:\:mt19937& rng: A random number generator. - :param measurement& A: Object whose class inherits :class:`measurement` and provides relevant measurement hooks. + :param measurement& A: Object whose class inherits :class:`measurement` and provides relevant measurement hooks. - .. function:: system::run_wolff(N_t N, std::function &, v_t)> r_gen, measurement& A, std::mt19937& rng) + .. function:: system::run_wolff(unsigned N, std::function &, const G_t::vertex&)> r_gen, measurement& A, std::mt19937& rng) Performs :any:`N` Wolff cluster flips to the system. - :param N_t N: Number of clusters to flip. - :param std\:\:function &, v_t>)> r_gen: Generator of transformations for the cluster flips. - :param measurement& A: Object whose class inherits :class:`measurement` and provides relevant measurement hooks. + :param unsigned N: Number of clusters to flip. + :param std\:\:function &, const G_t\:\:vertex&>)> r_gen: Generator of transformations for the cluster flips. + :param measurement& A: Object whose class inherits :class:`measurement` and provides relevant measurement hooks. :param std\:\:mt19937& rng: A random number generator. diff --git a/doc/types.rst b/doc/types.rst deleted file mode 100644 index b729ac3..0000000 --- a/doc/types.rst +++ /dev/null @@ -1,36 +0,0 @@ - -************ -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/On.cpp b/examples/On.cpp index 6c3f0fd..6885d2e 100644 --- a/examples/On.cpp +++ b/examples/On.cpp @@ -3,23 +3,21 @@ #include #include -#include "simple_measurement.hpp" - -#include -#include +#include +#include -#include +#include "simple_measurement.hpp" int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 0.8; vector_t H; H.fill(0.0); - q_t Hi = 0; + unsigned Hi = 0; int opt; @@ -27,7 +25,7 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) { switch (opt) { case 'N': // number of steps - N = (N_t)atof(optarg); + N = (unsigned)atof(optarg); break; case 'D': // dimension D = atoi(optarg); @@ -58,12 +56,12 @@ int main(int argc, char *argv[]) { }; // initialize the lattice - graph G(D, L); + graph<> G(D, L); // initialize the system system, vector_t> S(G, T, Z, B); - std::function (std::mt19937&, const system, vector_t>&, v_t)> gen_R = generate_rotation_uniform; + std::function (std::mt19937&, const system, vector_t, graph<>>&, const graph<>::vertex)> gen_R = generate_rotation_uniform>; // initailze the measurement object simple_measurement A(S); diff --git a/examples/clock.cpp b/examples/clock.cpp index 0cba8f6..3403f23 100644 --- a/examples/clock.cpp +++ b/examples/clock.cpp @@ -5,26 +5,24 @@ #define WOLFF_USE_FINITE_STATES #define WOLFF_FINITE_STATES_N WOLFF_POTTSQ - -#include -#include +#include +#include +#include #include "simple_measurement.hpp" -#include - using namespace wolff; int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 2.26918531421; vector_t<2, double> H; H.fill(0.0); - q_t Hi = 0; + unsigned Hi = 0; int opt; @@ -32,7 +30,7 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) { switch (opt) { case 'N': // number of steps - N = (N_t)atof(optarg); + N = (unsigned)atof(optarg); break; case 'D': // dimension D = atoi(optarg); @@ -63,22 +61,22 @@ int main(int argc, char *argv[]) { }; // initialize the lattice - graph G(D, L); + graph<> G(D, L); // initialize the system - system, potts_t> S(G, T, Z, B); + system, potts_t, graph<>> 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 (std::mt19937&, const system, potts_t>&, v_t)> gen_r = [] (std::mt19937& r, const system, potts_t>& S, v_t i0) -> dihedral_t { + std::function (std::mt19937&, const system, potts_t, graph<>>&, const graph<>::vertex&)> gen_r = [] (std::mt19937& r, const system, potts_t, graph<>>& S, const graph<>::vertex& v) -> dihedral_t { dihedral_t rot; rot.is_reflection = true; - std::uniform_int_distribution dist(0, WOLFF_POTTSQ - 2); - q_t x = dist(r); - rot.x = (2 * S.s[i0].x + x + 1) % WOLFF_POTTSQ; + std::uniform_int_distribution dist(0, WOLFF_POTTSQ - 2); + unsigned x = dist(r); + rot.x = (2 * S.s[v.ind].x + x + 1) % WOLFF_POTTSQ; return rot; }; diff --git a/examples/continuous_gaussian.cpp b/examples/continuous_gaussian.cpp index ef08247..eda22cb 100644 --- a/examples/continuous_gaussian.cpp +++ b/examples/continuous_gaussian.cpp @@ -3,19 +3,18 @@ #include #include -#include "simple_measurement.hpp" - -#include -#include - +#include +#include #include +#include "simple_measurement.hpp" + int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 0.8; double H = 0.0; @@ -25,7 +24,7 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) { switch (opt) { case 'N': // number of steps - N = (N_t)atof(optarg); + N = (unsigned)atof(optarg); break; case 'D': // dimension D = atoi(optarg); @@ -55,25 +54,25 @@ int main(int argc, char *argv[]) { }; // initialize the lattice - graph G(D, L); + graph<> G(D, L); // initialize the system - system, height_t> S(G, T, Z, B); + system, height_t, graph<>> S(G, T, Z, B); bool odd_run = false; - std::function (std::mt19937&, const system, height_t>&, v_t)> gen_R_IH = [&](std::mt19937& r, const system, height_t>& S, v_t i0) -> dihedral_inf_t { + std::function (std::mt19937&, const system, height_t, graph<>>&, const graph<>::vertex&)> gen_R_IH = [&](std::mt19937& r, const system, height_t, graph<>>& S, const graph<>::vertex& v) -> dihedral_inf_t { dihedral_inf_t rot; rot.is_reflection = true; if (odd_run) { - std::uniform_int_distribution dist(0, S.nv - 2); - v_t j = i0; + std::uniform_int_distribution dist(0, S.nv - 2); + unsigned j = v.ind; //while (S.s[j].x == S.s[i0].x) { - v_t tmp = dist(r); + unsigned tmp = dist(r); - if (tmp < i0) { + if (tmp < v.ind) { j = tmp; } else { j = tmp + 1; @@ -83,7 +82,7 @@ int main(int argc, char *argv[]) { rot.x = 2 * S.s[j].x; } else { std::normal_distribution dist(0.0,1.0); - rot.x = 2 * S.s[i0].x + dist(r); + rot.x = 2 * S.s[v.ind].x + dist(r); } odd_run = !odd_run; diff --git a/examples/discrete_gaussian.cpp b/examples/discrete_gaussian.cpp index 75ea0f9..a6d6ceb 100644 --- a/examples/discrete_gaussian.cpp +++ b/examples/discrete_gaussian.cpp @@ -3,19 +3,18 @@ #include #include -#include "simple_measurement.hpp" - -#include -#include - +#include +#include #include +#include "simple_measurement.hpp" + int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 0.8; double H = 0.0; @@ -25,7 +24,7 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) { switch (opt) { case 'N': // number of steps - N = (N_t)atof(optarg); + N = (unsigned)atof(optarg); break; case 'D': // dimension D = atoi(optarg); @@ -45,35 +44,35 @@ int main(int argc, char *argv[]) { } // define the spin-spin coupling - std::function &, const height_t&)> Z = [] (const height_t& s1, const height_t& s2) -> double { + std::function &, const height_t&)> Z = [] (const height_t& s1, const height_t& s2) -> double { return - pow(s1.x - s2.x, 2); }; // define the spin-field coupling - std::function &)> B = [&] (const height_t& s) -> double { + std::function &)> B = [&] (const height_t& s) -> double { return - H * pow(s.x, 2); }; // initialize the lattice - graph G(D, L); + graph<> G(D, L); // initialize the system - system, height_t> S(G, T, Z, B); + system, height_t, graph<>> S(G, T, Z, B); bool odd_run = false; - std::function (std::mt19937&, const system, height_t>&, v_t)> gen_R_IH = [&](std::mt19937& r, const system, height_t>& S, v_t i0) -> dihedral_inf_t { - dihedral_inf_t rot; + std::function (std::mt19937&, const system, height_t, graph<>>&, const graph<>::vertex&)> gen_R_IH = [&](std::mt19937& r, const system, height_t, graph<>>& S, const graph<>::vertex &v) -> dihedral_inf_t { + dihedral_inf_t rot; rot.is_reflection = true; if (odd_run) { - std::uniform_int_distribution dist(0, S.nv - 2); - v_t j = i0; + std::uniform_int_distribution dist(0, S.nv - 2); + unsigned j = v.ind; //while (S.s[j].x == S.s[i0].x) { - v_t tmp = dist(r); + unsigned tmp = dist(r); - if (tmp < i0) { + if (tmp < v.ind) { j = tmp; } else { j = tmp + 1; @@ -85,9 +84,9 @@ int main(int argc, char *argv[]) { std::uniform_int_distribution dist(0, 1); int j = dist(r); if (j) { - rot.x = 2 * S.s[i0].x + 1; + rot.x = 2 * S.s[v.ind].x + 1; } else { - rot.x = 2 * S.s[i0].x - 1; + rot.x = 2 * S.s[v.ind].x - 1; } } diff --git a/examples/ising.cpp b/examples/ising.cpp index c6e75fb..ecb296b 100644 --- a/examples/ising.cpp +++ b/examples/ising.cpp @@ -5,9 +5,7 @@ #define WOLFF_USE_FINITE_STATES -#include - -#include +#include #include "simple_measurement.hpp" @@ -16,9 +14,9 @@ using namespace wolff; int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 2.26918531421; double H = 0.0; @@ -28,7 +26,7 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) { switch (opt) { case 'N': // number of steps - N = (N_t)atof(optarg); + N = (unsigned)atof(optarg); break; case 'D': // dimension D = atoi(optarg); @@ -58,17 +56,17 @@ int main(int argc, char *argv[]) { }; // initialize the lattice - graph G(D, L); + graph<> G(D, L); // initialize the system - system S(G, 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 &, v_t)> gen_r = gen_ising; + std::function >&, const graph<>::vertex&)> gen_r = gen_ising>; // initailze the measurement object simple_measurement A(S); diff --git a/examples/ising_animation.cpp b/examples/ising_animation.cpp index 7bafcf8..bcaa589 100644 --- a/examples/ising_animation.cpp +++ b/examples/ising_animation.cpp @@ -6,19 +6,18 @@ #include #define WOLFF_USE_FINITE_STATES - -#include - -#include +#include using namespace wolff; -class draw_ising : public measurement { +typedef system> sys; + +class draw_ising : public measurement> { private: unsigned int frame_skip; - v_t C; + unsigned C; public: - draw_ising(const system& S, unsigned int window_size, unsigned int frame_skip, int argc, char *argv[]) : frame_skip(frame_skip){ + draw_ising(const sys& 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); @@ -29,9 +28,9 @@ class draw_ising : public measurement { 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&) { + void pre_cluster(unsigned, unsigned, const sys& S, const graph<>::vertex&, const ising_t&) override { glClear(GL_COLOR_BUFFER_BIT); - for (v_t i = 0; i < pow(S.G.L, 2); i++) { + for (unsigned 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 { @@ -43,9 +42,9 @@ class draw_ising : public measurement { C = 0; } - void plain_site_transformed(const system& S, v_t i, const ising_t&) { + void plain_site_transformed(const sys& S, const graph<>::vertex& v, const ising_t&) override { 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); + glRecti(v.ind / S.G.L, v.ind % S.G.L, (v.ind / S.G.L) + 1, (v.ind % S.G.L) + 1); C++; if (C % frame_skip == 0) { glFlush(); @@ -56,9 +55,9 @@ class draw_ising : public measurement { int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 2.26918531421; double H = 0.0; unsigned int window_size = 512; @@ -69,29 +68,29 @@ int main(int argc, char *argv[]) { // 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); + case 'N': // number of steps + N = (unsigned)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); } } @@ -106,15 +105,10 @@ int main(int argc, char *argv[]) { }; // initialize the lattice - graph G(D, L); + 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); - }; + sys S(G, T, Z, B); // initailze the measurement object draw_ising A(S, window_size, frame_skip, argc, argv); @@ -124,7 +118,7 @@ int main(int argc, char *argv[]) { std::mt19937 rng{seed}; // run wolff N times - S.run_wolff(N, gen_R, A, rng); + S.run_wolff(N, gen_ising>, A, rng); // exit return 0; diff --git a/examples/ising_no_field.cpp b/examples/ising_no_field.cpp index 28d7c70..0a5b722 100644 --- a/examples/ising_no_field.cpp +++ b/examples/ising_no_field.cpp @@ -5,22 +5,18 @@ #define WOLFF_NO_FIELD #define WOLFF_USE_FINITE_STATES - -#include - -#include +#include #include "simple_measurement.hpp" - using namespace wolff; int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 2.26918531421; int opt; @@ -29,7 +25,7 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:D:L:T:")) != -1) { switch (opt) { case 'N': // number of steps - N = (N_t)atof(optarg); + N = (unsigned)atof(optarg); break; case 'D': // dimension D = atoi(optarg); @@ -51,17 +47,17 @@ int main(int argc, char *argv[]) { }; // initialize the lattice - graph G(D, L); + graph<> G(D, L); // initialize the system - system S(G, T, Z); + system> S(G, T, Z); // 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 &, v_t)> gen_r = gen_ising; + std::function >&, const graph<>::vertex&)> gen_r = gen_ising>; // initailze the measurement object simple_measurement A(S); diff --git a/examples/ising_random_field.cpp b/examples/ising_random_field.cpp index b37b0ce..9284797 100644 --- a/examples/ising_random_field.cpp +++ b/examples/ising_random_field.cpp @@ -4,21 +4,18 @@ #include #define WOLFF_SITE_DEPENDENCE +#include #include "simple_measurement.hpp" -#include - -#include - using namespace wolff; int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 2.26918531421; double H = 0.0; @@ -28,7 +25,7 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) { switch (opt) { case 'N': // number of steps - N = (N_t)atof(optarg); + N = (unsigned)atof(optarg); break; case 'D': // dimension D = atoi(optarg); @@ -52,29 +49,29 @@ int main(int argc, char *argv[]) { return (double)(s1 * s2); }; - // initialize the lattice - graph G(D, L); + // initialize the lattice. vertex_prop is set to double and will contain the + // value of the random field at that site + graph G(D, L); // initialize the random number generator auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count(); std::mt19937 rng{seed}; // define the spin-field coupling - std::vector H_vals(G.nv); std::normal_distribution distribution(0.0, H); - for (v_t i = 0; i < G.nv; i++) { - H_vals[i] = distribution(rng); + for (auto &v : G.vertices) { + v.prop = distribution(rng); } - std::function B = [&] (v_t i, const ising_t& s) -> double { - return H_vals[i] * s; + std::function ::vertex&, const ising_t&)> B = [] (const graph::vertex& v, const ising_t& s) -> double { + return v.prop * s; }; // initialize the system - system S(G, T, Z, B); + system> S(G, T, Z, B); // define function that generates self-inverse rotations - std::function &, v_t)> gen_r = gen_ising; + std::function >&, const graph::vertex&)> gen_r = gen_ising>; // initailze the measurement object simple_measurement A(S); diff --git a/examples/ising_standalone.cpp b/examples/ising_standalone.cpp index 40572fd..6863ba5 100644 --- a/examples/ising_standalone.cpp +++ b/examples/ising_standalone.cpp @@ -21,27 +21,30 @@ class ising_t { } }; +typedef graph<> G_t; +typedef system sys; + class measure_clusters : public measurement { private: - v_t C; + unsigned C; public: double Ctotal; measure_clusters() { Ctotal = 0; } - void pre_cluster(N_t, N_t, const system&, v_t, const ising_t&) { C = 0; } + void pre_cluster(unsigned, unsigned, const sys&, const G_t::vertex&, const ising_t&) override { C = 0; } - void plain_site_transformed(const system&, v_t, const ising_t&) { C++; } + void plain_site_transformed(const sys&, const G_t::vertex&, const ising_t&) override { C++; } - void post_cluster(N_t, N_t, const system&) { Ctotal += C; } + void post_cluster(unsigned, unsigned, const sys&) override { Ctotal += C; } }; int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e3; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e3; + unsigned D = 2; + unsigned L = 128; double T = 2.26918531421; double H = 0.01; @@ -58,14 +61,14 @@ int main(int argc, char *argv[]) { }; // initialize the lattice - graph G(D, L); + G_t G(D, L); // initialize the system - system S(G, T, Z, B); + sys 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 { + std::function gen_R = + [] (std::mt19937&, const sys&, const G_t::vertex&) -> ising_t { return ising_t(-1); }; diff --git a/examples/potts.cpp b/examples/potts.cpp index 24be2b3..c15de8d 100644 --- a/examples/potts.cpp +++ b/examples/potts.cpp @@ -5,26 +5,24 @@ #define WOLFF_USE_FINITE_STATES #define WOLFF_FINITE_STATES_N WOLFF_POTTSQ - -#include -#include +#include +#include +#include #include "simple_measurement.hpp" -#include - using namespace wolff; int main(int argc, char *argv[]) { // set defaults - N_t N = (N_t)1e4; - D_t D = 2; - L_t L = 128; + unsigned N = (unsigned)1e4; + unsigned D = 2; + unsigned L = 128; double T = 2.26918531421; vector_t H; H.fill(0.0); - q_t Hi = 0; + unsigned Hi = 0; int opt; @@ -32,7 +30,7 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) { switch (opt) { case 'N': // number of steps - N = (N_t)atof(optarg); + N = (unsigned)atof(optarg); break; case 'D': // dimension D = atoi(optarg); @@ -67,30 +65,30 @@ int main(int argc, char *argv[]) { }; // initialize the lattice - graph G(D, L); + graph<> G(D, L); // initialize the system - system, potts_t> S(G, T, Z, B); + system, potts_t, graph<>> 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 (std::mt19937&, const system, potts_t>&, v_t)> gen_r = [] (std::mt19937& r, const system, potts_t>& S, v_t i0) -> symmetric_t { + std::function (std::mt19937&, const system, potts_t, graph<>>&, const graph<>::vertex&)> gen_r = [] (std::mt19937& r, const system, potts_t, graph<>>& S, const graph<>::vertex& v) -> symmetric_t { symmetric_t rot; - std::uniform_int_distribution dist(0, WOLFF_POTTSQ - 2); - q_t j = dist(r); - q_t swap_v; - if (j < S.s[i0].x) { + std::uniform_int_distribution dist(0, WOLFF_POTTSQ - 2); + unsigned j = dist(r); + unsigned swap_v; + if (j < S.s[v.ind].x) { swap_v = j; } else { swap_v = j + 1; } - rot[S.s[i0].x] = swap_v; - rot[swap_v] = S.s[i0].x; + rot[S.s[v.ind].x] = swap_v; + rot[swap_v] = S.s[v.ind].x; return rot; }; diff --git a/examples/simple_measurement.hpp b/examples/simple_measurement.hpp index 217f4f0..140da3b 100644 --- a/examples/simple_measurement.hpp +++ b/examples/simple_measurement.hpp @@ -1,31 +1,31 @@ -#include +#include using namespace wolff; -template -class simple_measurement : public measurement { +template +class simple_measurement : public measurement { private: - N_t n; + unsigned n; double E; typename X_t::M_t M; - v_t C; + unsigned C; double totalE; typename X_t::F_t totalM; double totalC; public: - simple_measurement(const system& S) { + simple_measurement(const system& S) { n = 0; M = S.nv * S.s[0]; E = 0; #ifdef WOLFF_BOND_DEPENDENCE - for (v_t i = 0; i < S.nv; i++) { - for (v_t j : S.G.adjacency[i]) { - E -= 0.5 * S.Z(i, S.s[i], j, S.s[j]); + for (unsigned i = 0; i < S.nv; i++) { + for (const typename G_t::halfedge& e : S.G.vertices[i].edges) { + E -= 0.5 * S.Z(e, S.s[i], S.s[j]); } } #else @@ -34,8 +34,8 @@ class simple_measurement : public measurement { #ifndef WOLFF_NO_FIELD #ifdef WOLFF_SITE_DEPENDENCE - for (v_t i = 0; i < S.nv; i++) { - E -= S.B(i, S.s[i]); + for (unsigned i = 0; i < S.nv; i++) { + E -= S.B(S.G.vertices[i], S.s[i]); } #else E -= S.nv * S.B(S.s[0]); @@ -47,28 +47,30 @@ class simple_measurement : public measurement { totalC = 0; } - void pre_cluster(N_t, N_t, const system&, v_t, const R_t&) { + void pre_cluster(unsigned, unsigned, const system&, const typename G_t::vertex&, const R_t&) override { C = 0; } - void plain_bond_visited(const system&, v_t, const X_t&, v_t, double dE) { + void plain_bond_visited(const system&, const typename G_t::halfedge&, const X_t&, double dE) override { E += dE; } - void ghost_bond_visited(const system&, v_t, const X_t& s_old, const X_t& s_new, double dE) { +#ifndef WOLFF_NO_FIELD + void ghost_bond_visited(const system&, const typename G_t::vertex&, const X_t& s_old, const X_t& s_new, double dE) override { E += dE; M += s_new - s_old; } +#endif - void plain_site_transformed(const system& S, v_t i, const X_t& si_new) { + void plain_site_transformed(const system& S, const typename G_t::vertex& v, const X_t& si_new) override { C++; #ifdef WOLFF_NO_FIELD - M += si_new - S.s[i]; + M += si_new - S.s[v.ind]; #endif } - void post_cluster(N_t, N_t, const system&) { + void post_cluster(unsigned, unsigned, const system&) override { totalE += E; totalM += M; totalC += C; diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 1df9807..445e68f 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -1,14 +1,10 @@ -project(libwolff LANGUAGES C CXX) +project(libwolff LANGUAGES CXX) -add_library(wolff SHARED src/graph.cpp) +add_library(wolff INTERFACE) +target_include_directories(wolff INTERFACE .) +target_compile_features(wolff INTERFACE cxx_std_17) -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(FILES "wolff.hpp" DESTINATION include) +install(DIRECTORY "wolff_models" DESTINATION include) diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp deleted file mode 100644 index 7e909c8..0000000 --- a/lib/include/wolff.hpp +++ /dev/null @@ -1,32 +0,0 @@ - -#ifndef WOLFF_H -#define WOLFF_H - -#include "wolff/cluster.hpp" -#include "wolff/measurement.hpp" - -namespace wolff{ - -template -void system::run_wolff(N_t N, - std::function &, v_t)> r_gen, - measurement& A, std::mt19937& rng) { - - 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, *this, i0); - - A.pre_cluster(n, N, *this, i0, r); - - this->flip_cluster(i0, r, rng, A); - - A.post_cluster(n, N, *this); - } -} - -} - -#endif - diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp deleted file mode 100644 index 3912446..0000000 --- a/lib/include/wolff/cluster.hpp +++ /dev/null @@ -1,120 +0,0 @@ - -#ifndef WOLFF_CLUSTER_H -#define WOLFF_CLUSTER_H - -#include -#include -#include -#include - -#include "system.hpp" - -namespace wolff { - -template -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(G.nv, false); - - while (!queue.empty()) { - v_t i = queue.front(); - queue.pop(); - - if (!marks[i]) { // don't reprocess anyone we've already visited! - marks[i] = true; - - X_t si_new; -#ifndef WOLFF_NO_FIELD - R_t s0_new; - - bool we_are_ghost = (i == nv); - - if (we_are_ghost) { - s0_new = r.act(s0); - } else -#endif - { - si_new = r.act(s[i]); - } - - for (const v_t &j : G.adjacency[i]) { - double dE, p; - -#ifndef WOLFF_NO_FIELD - bool neighbor_is_ghost = (j == nv); - - if (we_are_ghost || neighbor_is_ghost) { - X_t s0s_old, s0s_new; - v_t non_ghost; - - if (neighbor_is_ghost) { - non_ghost = i; - s0s_old = s0.act_inverse(s[i]); - s0s_new = s0.act_inverse(si_new); - } else { - non_ghost = j; - s0s_old = s0.act_inverse(s[j]); - s0s_new = s0_new.act_inverse(s[j]); - } - -#ifdef WOLFF_SITE_DEPENDENCE - dE = B(non_ghost, s0s_old) - B(non_ghost, s0s_new); -#else - dE = B(s0s_old) - B(s0s_new); -#endif - -#ifdef WOLFF_USE_FINITE_STATES - p = Bp[s0s_old.enumerate()][s0s_new.enumerate()]; -#endif - - // run measurement hooks for encountering a ghost bond - 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 = Z(i, s[i], j, s[j]) - Z(i, si_new, j, s[j]); -#else - dE = Z(s[i], s[j]) - Z(si_new, s[j]); -#endif - -#ifdef WOLFF_USE_FINITE_STATES - p = Zp[s[i].enumerate()][si_new.enumerate()][s[j].enumerate()]; -#endif - - // run measurement hooks for encountering a plain bond - A.plain_bond_visited(*this, i, si_new, j, dE); - } - -#ifndef WOLFF_USE_FINITE_STATES - p = 1.0 - exp(-dE / T); -#endif - - if (dist(rng) < p) { - queue.push(j); // push the neighboring vertex to the queue - } - } - -#ifndef WOLFF_NO_FIELD - if (we_are_ghost) { - A.ghost_site_transformed(*this, s0_new); - s0 = s0_new; - } else -#endif - { - A.plain_site_transformed(*this, i, si_new); - s[i] = si_new; - } - } - } -} - -} - -#endif - diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp deleted file mode 100644 index 65b3941..0000000 --- a/lib/include/wolff/graph.hpp +++ /dev/null @@ -1,30 +0,0 @@ - -#ifndef WOLFF_GRAPH_H -#define WOLFF_GRAPH_H - -#include -#include - -namespace wolff { - -#include "types.h" - -class graph { - public: - D_t D; - L_t L; - v_t ne; - v_t nv; - std::vector> adjacency; - std::vector> coordinate; - - 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 deleted file mode 100644 index 6e352a5..0000000 --- a/lib/include/wolff/measurement.hpp +++ /dev/null @@ -1,28 +0,0 @@ - -#ifndef WOLFF_MEASUREMENTS_H -#define WOLFF_MEASUREMENTS_H - -#include "system.hpp" - -namespace wolff { - -template -class measurement { - public: - virtual void pre_cluster(N_t, N_t, const system&, v_t, const R_t&) {}; - - virtual void plain_bond_visited(const system&, v_t, const X_t&, v_t, double) {}; - virtual void plain_site_transformed(const system&, v_t, const X_t&) {}; - -#ifndef WOLFF_NO_FIELD - virtual void ghost_bond_visited(const system&, v_t, const X_t&, const X_t&, double) {}; - virtual void ghost_site_transformed(const system&, const R_t&) {}; -#endif - - virtual void post_cluster(N_t, N_t, const system&) {}; -}; - -} - -#endif - diff --git a/lib/include/wolff/models/dihedral.hpp b/lib/include/wolff/models/dihedral.hpp deleted file mode 100644 index cf4b067..0000000 --- a/lib/include/wolff/models/dihedral.hpp +++ /dev/null @@ -1,56 +0,0 @@ - -#ifndef WOLFF_MODELS_DIHEDRAL_H -#define WOLFF_MODELS_DIHEDRAL_H - -#include "potts.hpp" - -namespace wolff { - -#include "../types.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); - } - } -}; - -} - -#endif - diff --git a/lib/include/wolff/models/dihedral_inf.hpp b/lib/include/wolff/models/dihedral_inf.hpp deleted file mode 100644 index 33d5cec..0000000 --- a/lib/include/wolff/models/dihedral_inf.hpp +++ /dev/null @@ -1,57 +0,0 @@ - -#ifndef WOLFF_MODELS_DIHEDRAL_INF_H -#define WOLFF_MODELS_DIHEDRAL_INF_H - -#include -#include "height.hpp" - -namespace wolff { - -#include "../types.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); - } - } -}; - -} - -#endif - diff --git a/lib/include/wolff/models/height.hpp b/lib/include/wolff/models/height.hpp deleted file mode 100644 index 9d7e87c..0000000 --- a/lib/include/wolff/models/height.hpp +++ /dev/null @@ -1,57 +0,0 @@ - -#pragma once - -#include -#include - -namespace wolff { - -template -struct height_t { - T x; - - height_t() : x(0) {} - height_t(T x) : x(x) {} - - typedef T M_t; - typedef double F_t; - - 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 typename height_t::M_t operator*(v_t a, height_t h) { - return h * a; -} - -template -inline typename height_t::F_t operator*(double a, height_t h) { - return h * a; -} - -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; -} - -} - diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp deleted file mode 100644 index d2f6f91..0000000 --- a/lib/include/wolff/models/ising.hpp +++ /dev/null @@ -1,93 +0,0 @@ - -#ifndef WOLFF_MODELS_ISING_H -#define WOLFF_MODELS_ISING_H - -#define WOLFF_FINITE_STATES_N 2 - -#include "../system.hpp" - -namespace wolff { - -#include "../types.h" - -class ising_t { - public: - bool x; - - ising_t() : x(false) {} - - ising_t(bool x) : x(x) {} - ising_t(q_t x) : x((bool)x) {} - - ising_t act(const ising_t& s) const { - if (x) { - return ising_t(!s.x); - } else { - return ising_t(s.x); - } - } - - ising_t act_inverse(const ising_t& s) const { - return this->act(s); - } - - typedef int M_t; - typedef double F_t; - - inline M_t operator*(v_t a) const { - if (x) { - return -(M_t)a; - } else { - return (M_t)a; - } - } - - inline F_t operator*(double a) const { - if (x) { - return -a; - } else { - return a; - } - } - - inline M_t operator-(const ising_t &s) const { - if (x == s.x) { - return 0; - } else { - if (x) { - return -2; - } else { - return 2; - } - } - } - - inline int operator*(const ising_t& s) const { - if (x == s.x) { - return 1; - } else { - return -1; - } - } - - q_t enumerate() const { - return (q_t)x; - } -}; - -inline ising_t::M_t operator*(v_t a, const ising_t& s) { - return s * a; -} - -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); -}; - -} - -#endif - diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp deleted file mode 100644 index 514c88a..0000000 --- a/lib/include/wolff/models/orthogonal.hpp +++ /dev/null @@ -1,208 +0,0 @@ - -#ifndef WOLFF_MODELS_ORTHOGONAL_H -#define WOLFF_MODELS_ORTHOGONAL_H - -#include -#include - -#include "../system.hpp" -#include "vector.hpp" - -namespace wolff { - -#include "../types.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 (std::mt19937& r, const system, vector_t>&, v_t) { - 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 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; - - 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] = 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] = S.s[i0][i]; - } - } else { - v = S.s[i0]; - } - - 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; -} - -} - -#endif - diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp deleted file mode 100644 index a0cb368..0000000 --- a/lib/include/wolff/models/potts.hpp +++ /dev/null @@ -1,68 +0,0 @@ - -#ifndef WOLFF_MODELS_POTTS_H -#define WOLFF_MODELS_POTTS_H - -#include - -#include "vector.hpp" - -namespace wolff { - -#include "../types.h" - -template -class potts_t { - public: - q_t x; - - potts_t() : x(0) {} - potts_t(q_t x) : x(x) {} - - typedef vector_t M_t; - typedef vector_t F_t; - - 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; - } - - q_t enumerate() const { - return x; - } -}; - -template -inline typename potts_t::M_t operator*(v_t a, const potts_t& s) { - return s * a; -} - -template -inline typename potts_t::F_t operator*(double a, const potts_t& s) { - return s * a; -} - -} - -#endif - diff --git a/lib/include/wolff/models/symmetric.hpp b/lib/include/wolff/models/symmetric.hpp deleted file mode 100644 index d15a7ba..0000000 --- a/lib/include/wolff/models/symmetric.hpp +++ /dev/null @@ -1,59 +0,0 @@ - -#ifndef WOLFF_MODELS_SYMMETRIC_H -#define WOLFF_MODELS_SYMMETRIC_H - -#include - -#include "potts.hpp" - -namespace wolff { - -#include "../types.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; - } -}; - -} - -#endif - diff --git a/lib/include/wolff/models/vector.hpp b/lib/include/wolff/models/vector.hpp deleted file mode 100644 index e4c4a1c..0000000 --- a/lib/include/wolff/models/vector.hpp +++ /dev/null @@ -1,115 +0,0 @@ - -#ifndef WOLFF_MODELS_VECTOR_H -#define WOLFF_MODELS_VECTOR_H - -#include -#include -#include - -namespace wolff { - -#include - -template -class vector_t : public std::array { - public: - - vector_t() { - this->fill((T)0); - (*this)[0] = (T)1; - } - - vector_t(const T *x) { - for (q_t i = 0; i < q; i++) { - (*this)[i] = x[i]; - } - } - - typedef vector_t M_t; - typedef vector_t F_t; - - template - inline vector_t& operator+=(const vector_t &v) { - for (q_t i = 0; i < q; i++) { - (*this)[i] += (T)v[i]; - } - return *this; - } - - template - inline vector_t& operator-=(const vector_t &v) { - for (q_t i = 0; i < q; i++) { - (*this)[i] -= (T)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; - } - - inline T operator*(const vector_t& v) const { - double prod = 0; - - for (q_t i = 0; i < q; i++) { - prod += v[i] * (*this)[i]; - } - - return prod; - } - - template - inline vector_t operator/(U a) const { - vector_t result; - for (q_t i = 0; i < q; i++) { - result[i] = (*this)[i] / a; - } - - return result; - } -}; - -template -inline vector_t operator*(v_t a, const vector_t&v) { - return v * a; -} - -template -inline vector_t operator*(double a, const vector_t&v) { - return v * a; -} - -template -std::ostream& operator<<(std::ostream& os, const vector_t&v) { - os << "( "; - for (T vi : v) { - os << vi << " "; - } - os << ")"; - return os; -} - -} - -#endif - diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp deleted file mode 100644 index fa9a1b0..0000000 --- a/lib/include/wolff/system.hpp +++ /dev/null @@ -1,105 +0,0 @@ - -#ifndef WOLFF_SYSTEM_H -#define WOLFF_SYSTEM_H - -#include -#include -#include - -#include "graph.hpp" - -namespace wolff { - -#include "types.h" - -template -class measurement; - -template -class system { - public: - v_t nv; // number of vertices - v_t ne; // number of edges - graph G; // the graph defining the lattice with ghost - double T; // the temperature - std::vector s; // the state of the ordinary spins - -#ifdef WOLFF_BOND_DEPENDENCE - std::function Z; // coupling between sites -#else - std::function Z; // coupling between sites -#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 - std::function B; // coupling with the external field -#endif -#endif - -#ifdef WOLFF_USE_FINITE_STATES - std::array, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp; -#ifndef WOLFF_NO_FIELD - std::array, WOLFF_FINITE_STATES_N> Bp; -#endif -#endif - - system(graph g, double T, -#ifdef WOLFF_BOND_DEPENDENCE - std::function Z -#else - std::function Z -#endif -#ifndef WOLFF_NO_FIELD -#ifdef WOLFF_SITE_DEPENDENCE - , std::function B -#else - , std::function B -#endif -#endif - ) : G(g), T(T), Z(Z) -#ifndef WOLFF_NO_FIELD - , s0(), B(B) -#endif - { - nv = G.nv; - ne = G.ne; - s.resize(nv); -#ifndef WOLFF_NO_FIELD - G.add_ghost(); -#endif -#ifdef WOLFF_USE_FINITE_STATES - this->finite_states_init(); -#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); - -#ifdef WOLFF_USE_FINITE_STATES - void finite_states_init() { -#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++) { - Bp[i][j] = 1.0 - exp(-(B(X_t(i)) - B(X_t(j))) / T); - } - } -#endif - for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) { - for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) { - for (q_t k = 0; k < WOLFF_FINITE_STATES_N; k++) { - Zp[i][j][k] = 1.0 - exp(-(Z(X_t(i), X_t(k)) - Z(X_t(j), X_t(k))) / T); - } - } - } - } -#endif - -}; - -} - -#endif - diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h deleted file mode 100644 index 5bbaa6d..0000000 --- a/lib/include/wolff/types.h +++ /dev/null @@ -1,32 +0,0 @@ - -#ifndef WOLFF_TYPES_H -#define WOLFF_TYPES_H - -#include - -typedef uint_fast32_t v_t; // vertex and edge indices -typedef uint_fast8_t q_t; // enumerated states -typedef uint_fast8_t D_t; // dimension -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_D UINT_FAST8_MAX -#define MAX_L UINT_FAST16_MAX -#define MAX_N UINT_FAST64_MAX - -#define PRIv PRIuFAST32 -#define PRIq PRIuFAST8 -#define PRID PRIuFAST8 -#define PRIL PRIuFAST16 -#define PRIN PRIuFAST64 - -#define SCNv SCNuFAST32 -#define SCNq SCNuFAST8 -#define SCND SCNuFAST8 -#define SCNL SCNuFAST16 -#define SCNN SCNuFAST64 - -#endif - diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp deleted file mode 100644 index 4f91be4..0000000 --- a/lib/src/graph.cpp +++ /dev/null @@ -1,55 +0,0 @@ - -#include - -namespace wolff { - -graph::graph() { - D = 0; - L = 0; - nv = 0; - ne = 0; -} - -graph::graph(D_t D, L_t L) : D(D), L(L) { - nv = pow(L, D); - ne = D * nv; - - adjacency.resize(nv); - coordinate.resize(nv); - - for (std::vector adj_i : adjacency) { - 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; - - 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::add_ghost() { - for (std::vector& adj_i : adjacency) { - adj_i.push_back(nv); - } - - adjacency.resize(nv + 1); - coordinate.resize(nv + 1); - adjacency[nv].reserve(nv); - - for (v_t i = 0; i < nv; i++) { - adjacency[nv].push_back(i); - } - - coordinate[nv].resize(coordinate[0].size()); - - ne += nv; - nv++; -} - -} - diff --git a/lib/wolff.hpp b/lib/wolff.hpp new file mode 100644 index 0000000..d21d788 --- /dev/null +++ b/lib/wolff.hpp @@ -0,0 +1,320 @@ + +#ifndef WOLFF_H +#define WOLFF_H + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace wolff{ + + template , class edge_prop = std::tuple<>> + class graph { + /* the graph class describes a lattice on which wolff is run. for most + * purposes, using the default square lattice constructor is sufficient, + * but arbitrary lattices can be constructed + */ + public: + unsigned D; // dimension of space + unsigned L; // linear size + unsigned ne; // number of edges + unsigned nv; // number of vertices + + struct _vertex; + + typedef struct _halfedge { + struct _vertex &self; // reference to the vertex this edge comes from + struct _vertex &neighbor; // reference to the vertex this edge goes to + edge_prop prop; // optional properties + + _halfedge(struct _vertex &v1, struct _vertex &v2) : self(v1), neighbor(v2) {}; + } halfedge; + + typedef struct _vertex { + unsigned ind; // index of the vertex + std::list edges; // list of edges incident on this vertex + vertex_prop prop; // optional properties + } vertex; + + std::vector vertices; /* vector of vertices, length nv, with + * vertices[i].ind = i for all 0 <= i < nv + */ + graph() { + // default constructor for empty graph. use it to build your own! + D = 0; + L = 0; + nv = 0; + ne = 0; + }; + + graph(unsigned D, unsigned L) : D(D), L(L) { + // default constructor for square lattice graph + nv = pow(L, D); + ne = D * nv; + + vertices.resize(nv); + + for (unsigned i = 0; i < nv; i++) { + vertices[i].ind = i; + for (unsigned j = 0; j < D; j++) { + unsigned n1 = pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) + + fmod(i + pow(L, j), pow(L, j + 1)); + unsigned n2 = pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) + + fmod(pow(L, j + 1) + i - pow(L, j), pow(L, j + 1)); + + halfedge f(vertices[i], vertices[n1]); + halfedge b(vertices[i], vertices[n2]); + + vertices[i].edges.push_back(f); + vertices[i].edges.push_back(b); + } + } + }; + + void add_ghost() { + // adds a ghost site to any graph + vertices.resize(nv + 1); + + for (auto it = vertices.begin(); it != std::prev(vertices.end()); it++) { + halfedge e1(*it, vertices[nv]); + it->edges.push_back(e1); + + halfedge e2(vertices[nv], *it); + vertices[nv].edges.push_back(e2); + } + + vertices[nv].ind = nv; + + ne += nv; + nv++; + }; + }; + + template > + class measurement; + + template > + class system { + public: + unsigned nv; // number of vertices + unsigned ne; // number of edges + G_t G; // the graph defining the lattice with ghost + double T; // the temperature + std::vector s; // the state of the ordinary spins + +#ifdef WOLFF_BOND_DEPENDENCE + std::function Z; // coupling between sites +#else + std::function Z; // coupling between sites +#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 + std::function B; // coupling with the external field +#endif +#endif + +#ifdef WOLFF_USE_FINITE_STATES + std::array, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp; +#ifndef WOLFF_NO_FIELD + std::array, WOLFF_FINITE_STATES_N> Bp; +#endif +#endif + + system(G_t g, double T, +#ifdef WOLFF_BOND_DEPENDENCE + std::function Z +#else + std::function Z +#endif +#ifndef WOLFF_NO_FIELD +#ifdef WOLFF_SITE_DEPENDENCE + , std::function B +#else + , std::function B +#endif +#endif + ) : G(g), T(T), Z(Z) +#ifndef WOLFF_NO_FIELD + , s0(), B(B) +#endif + { + nv = G.nv; + ne = G.ne; + s.resize(nv); +#ifndef WOLFF_NO_FIELD + G.add_ghost(); +#endif +#ifdef WOLFF_USE_FINITE_STATES + this->finite_states_init(); +#endif + } + + void flip_cluster(typename G_t::vertex& v0, const R_t& r, std::mt19937& rng, + measurement& A) { + std::uniform_real_distribution dist(0.0, 1.0); + + std::queue queue; + queue.push(v0.ind); + + std::vector visited(G.nv, false); + + while (!queue.empty()) { + unsigned i = queue.front(); + queue.pop(); + + if (!visited[i]) { // don't reprocess anyone we've already visited! + visited[i] = true; + + X_t si_new; +#ifndef WOLFF_NO_FIELD + R_t s0_new; + + bool we_are_ghost = (i == nv); + + if (we_are_ghost) { + s0_new = r.act(s0); + } else +#endif + { + si_new = r.act(s[i]); + } + + for (const typename G_t::halfedge &e : G.vertices[i].edges) { + double dE, p; + unsigned j = e.neighbor.ind; + +#ifndef WOLFF_NO_FIELD + bool neighbor_is_ghost = (j == nv); + + if (we_are_ghost || neighbor_is_ghost) { + X_t s0s_old, s0s_new; + unsigned non_ghost; + + if (neighbor_is_ghost) { + non_ghost = i; + s0s_old = s0.act_inverse(s[i]); + s0s_new = s0.act_inverse(si_new); + } else { + non_ghost = j; + s0s_old = s0.act_inverse(s[j]); + s0s_new = s0_new.act_inverse(s[j]); + } + +#ifdef WOLFF_SITE_DEPENDENCE + dE = B(G.vertices[non_ghost], s0s_old) - B(G.vertices[non_ghost], s0s_new); +#else + dE = B(s0s_old) - B(s0s_new); +#endif + +#ifdef WOLFF_USE_FINITE_STATES + p = Bp[s0s_old.enumerate()][s0s_new.enumerate()]; +#endif + + // run measurement hooks for encountering a ghost bond + A.ghost_bond_visited(*this, G.vertices[non_ghost], s0s_old, s0s_new, dE); + } else // this is a perfectly normal bond! +#endif + { +#ifdef WOLFF_BOND_DEPENDENCE + dE = Z(e, s[i], s[j]) - Z(e, si_new, s[j]); +#else + dE = Z(s[i], s[j]) - Z(si_new, s[j]); +#endif + +#ifdef WOLFF_USE_FINITE_STATES + p = Zp[s[i].enumerate()][si_new.enumerate()][s[j].enumerate()]; +#endif + + // run measurement hooks for encountering a plain bond + A.plain_bond_visited(*this, e, si_new, dE); + } + +#ifndef WOLFF_USE_FINITE_STATES + p = 1.0 - exp(-dE / T); +#endif + + if (dist(rng) < p) { + queue.push(j); // push the neighboring vertex to the queue + } + } + +#ifndef WOLFF_NO_FIELD + if (we_are_ghost) { + A.ghost_site_transformed(*this, s0_new); + s0 = s0_new; + } else +#endif + { + A.plain_site_transformed(*this, G.vertices[i], si_new); + s[i] = si_new; + } + } + } + } + + void run_wolff(unsigned N, + std::function &, const typename G_t::vertex&)> r_gen, measurement& A, std::mt19937& rng) { + std::uniform_int_distribution dist(0, nv - 1); + + for (unsigned n = 0; n < N; n++) { + unsigned i0 = dist(rng); + R_t r = r_gen(rng, *this, G.vertices[i0]); + + A.pre_cluster(n, N, *this, G.vertices[i0], r); + + this->flip_cluster(G.vertices[i0], r, rng, A); + + A.post_cluster(n, N, *this); + } + } + +#ifdef WOLFF_USE_FINITE_STATES + void finite_states_init() { +#ifndef WOLFF_NO_FIELD + for (unsigned i = 0; i < WOLFF_FINITE_STATES_N; i++) { + for (unsigned j = 0; j < WOLFF_FINITE_STATES_N; j++) { + Bp[i][j] = 1.0 - exp(-(B(X_t(i)) - B(X_t(j))) / T); + } + } +#endif + for (unsigned i = 0; i < WOLFF_FINITE_STATES_N; i++) { + for (unsigned j = 0; j < WOLFF_FINITE_STATES_N; j++) { + for (unsigned k = 0; k < WOLFF_FINITE_STATES_N; k++) { + Zp[i][j][k] = 1.0 - exp(-(Z(X_t(i), X_t(k)) - Z(X_t(j), X_t(k))) / T); + } + } + } + } +#endif + + }; + + template + class measurement { + public: + virtual void pre_cluster(unsigned, unsigned, const system&, const typename G_t::vertex& v, const R_t&) {}; + + virtual void plain_bond_visited(const system&, const typename G_t::halfedge& e, const X_t&, double) {}; + virtual void plain_site_transformed(const system&, const typename G_t::vertex& v, const X_t&) {}; + +#ifndef WOLFF_NO_FIELD + virtual void ghost_bond_visited(const system&, const typename G_t::vertex& v, const X_t&, const X_t&, double) {}; + virtual void ghost_site_transformed(const system&, const R_t&) {}; +#endif + + virtual void post_cluster(unsigned, unsigned, const system&) {}; + }; + +} + +#endif + diff --git a/lib/wolff_models/dihedral.hpp b/lib/wolff_models/dihedral.hpp new file mode 100644 index 0000000..aef3be6 --- /dev/null +++ b/lib/wolff_models/dihedral.hpp @@ -0,0 +1,54 @@ + +#ifndef WOLFF_MODELS_DIHEDRAL_H +#define WOLFF_MODELS_DIHEDRAL_H + +#include "potts.hpp" + +namespace wolff { + + template + class dihedral_t { + public: + bool is_reflection; + unsigned x; + + dihedral_t() : is_reflection(false), x(0) {} + dihedral_t(bool x, unsigned 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); + } + } + }; + +} + +#endif + diff --git a/lib/wolff_models/dihedral_inf.hpp b/lib/wolff_models/dihedral_inf.hpp new file mode 100644 index 0000000..0b09c59 --- /dev/null +++ b/lib/wolff_models/dihedral_inf.hpp @@ -0,0 +1,55 @@ + +#ifndef WOLFF_MODELS_DIHEDRAL_INF_H +#define WOLFF_MODELS_DIHEDRAL_INF_H + +#include +#include "height.hpp" + +namespace wolff { + + 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); + } + } + }; + +} + +#endif + diff --git a/lib/wolff_models/height.hpp b/lib/wolff_models/height.hpp new file mode 100644 index 0000000..ae3ad82 --- /dev/null +++ b/lib/wolff_models/height.hpp @@ -0,0 +1,56 @@ + +#pragma once + +#include + +namespace wolff { + + template + struct height_t { + T x; + + height_t() : x(0) {} + height_t(T x) : x(x) {} + + typedef T M_t; + typedef double F_t; + + inline T operator*(unsigned 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 typename height_t::M_t operator*(unsigned a, height_t h) { + return h * a; + } + + template + inline typename height_t::F_t operator*(double a, height_t h) { + return h * a; + } + + 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; + } + +} + diff --git a/lib/wolff_models/ising.hpp b/lib/wolff_models/ising.hpp new file mode 100644 index 0000000..345eeff --- /dev/null +++ b/lib/wolff_models/ising.hpp @@ -0,0 +1,92 @@ + +#ifndef WOLFF_MODELS_ISING_H +#define WOLFF_MODELS_ISING_H + +#define WOLFF_FINITE_STATES_N 2 + +#include "wolff.hpp" + +namespace wolff { + + class ising_t { + public: + bool x; + + ising_t() : x(false) {} + + ising_t(bool x) : x(x) {} + ising_t(unsigned x) : x((bool)x) {} + + ising_t act(const ising_t& s) const { + if (x) { + return ising_t(!s.x); + } else { + return ising_t(s.x); + } + } + + ising_t act_inverse(const ising_t& s) const { + return this->act(s); + } + + typedef int M_t; + typedef double F_t; + + inline M_t operator*(unsigned a) const { + if (x) { + return -(M_t)a; + } else { + return (M_t)a; + } + } + + inline F_t operator*(double a) const { + if (x) { + return -a; + } else { + return a; + } + } + + inline M_t operator-(const ising_t &s) const { + if (x == s.x) { + return 0; + } else { + if (x) { + return -2; + } else { + return 2; + } + } + } + + inline int operator*(const ising_t& s) const { + if (x == s.x) { + return 1; + } else { + return -1; + } + } + + unsigned enumerate() const { + return (unsigned)x; + } + }; + + inline ising_t::M_t operator*(unsigned a, const ising_t& s) { + return s * a; + } + + inline ising_t::F_t operator*(double a, const ising_t& s) { + return s * a; + } + + template + ising_t gen_ising(std::mt19937&, const system&, const typename G_t::vertex&) { + return ising_t(true); + }; + +} + +#endif + diff --git a/lib/wolff_models/orthogonal.hpp b/lib/wolff_models/orthogonal.hpp new file mode 100644 index 0000000..ea9af14 --- /dev/null +++ b/lib/wolff_models/orthogonal.hpp @@ -0,0 +1,206 @@ + +#ifndef WOLFF_MODELS_ORTHOGONAL_H +#define WOLFF_MODELS_ORTHOGONAL_H + +#include +#include + +#include +#include "vector.hpp" + +namespace wolff { + + template + class orthogonal_t : public std::array, q> { + public : + bool is_reflection; + + orthogonal_t() : is_reflection(false) { + for (unsigned 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 (unsigned i = 0; i < q; i++) { + prod += v[i] * (*this)[0][i]; + } + for (unsigned i = 0; i < q; i++) { + v_rot[i] = v[i] - 2 * prod * (*this)[0][i]; + } + } else { + for (unsigned i = 0; i < q; i++) { + for (unsigned 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 (unsigned i = 0; i < q; i++) { + double akOki = 0; + + for (unsigned k = 0; k < q; k++) { + akOki += (*this)[0][k] * m[k][i]; + } + + for (unsigned j = 0; j < q; j++) { + m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][j]; + } + } + } else { + for (unsigned i = 0; i < q; i++) { + m_rot[i].fill(0); + for (unsigned j = 0; j < q; j++) { + for (unsigned 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 (unsigned i = 0; i < q; i++) { + for (unsigned 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 (unsigned i = 0; i < q; i++) { + m_rot[i].fill(0); + for (unsigned j = 0; j < q; j++) { + for (unsigned 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 system, vector_t, G_t>&, const typename G_t::vertex&) { + std::normal_distribution dist(0.0,1.0); + orthogonal_t ptr; + ptr.is_reflection = true; + + double v2 = 0; + + for (unsigned i = 0; i < q; i++) { + ptr[0][i] = dist(r); + v2 += ptr[0][i] * ptr[0][i]; + } + + double mag_v = sqrt(v2); + + for (unsigned i = 0; i < q; i++) { + ptr[0][i] /= mag_v; + } + + return ptr; + } + + template + orthogonal_t generate_rotation_perturbation (std::mt19937& r, const system, vector_t, G_t>& S, const typename G_t::vertex& 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] = S.s[v0.ind][0] * cosr - S.s[v0.ind][1] * sinr; + v[1] = S.s[v0.ind][1] * cosr + S.s[v0.ind][0] * sinr; + + for (unsigned i = 2; i < q; i++) { + v[i] = S.s[v0.ind][i]; + } + } else { + v = S.s[v0.ind]; + } + + double m_dot_v = 0; + + for (unsigned 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 (unsigned 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 (unsigned i = 0; i < q; i++) { + m[0][i] /= mag_v; // normalize + } + + v2 = 0; + + double factor = epsilon * dist(r); + + for (unsigned 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 (unsigned i = 0; i < q; i++) { + m[0][i] /= mag_v; // normalize + } + + return m; + } + +} + +#endif + diff --git a/lib/wolff_models/potts.hpp b/lib/wolff_models/potts.hpp new file mode 100644 index 0000000..dd56134 --- /dev/null +++ b/lib/wolff_models/potts.hpp @@ -0,0 +1,66 @@ + +#ifndef WOLFF_MODELS_POTTS_H +#define WOLFF_MODELS_POTTS_H + +#include + +#include "vector.hpp" + +namespace wolff { + + template + class potts_t { + public: + unsigned x; + + potts_t() : x(0) {} + potts_t(unsigned x) : x(x) {} + + typedef vector_t M_t; + typedef vector_t F_t; + + inline vector_t operator*(unsigned 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; + } + + unsigned enumerate() const { + return x; + } + }; + + template + inline typename potts_t::M_t operator*(unsigned a, const potts_t& s) { + return s * a; + } + + template + inline typename potts_t::F_t operator*(double a, const potts_t& s) { + return s * a; + } + +} + +#endif + diff --git a/lib/wolff_models/symmetric.hpp b/lib/wolff_models/symmetric.hpp new file mode 100644 index 0000000..98e53e3 --- /dev/null +++ b/lib/wolff_models/symmetric.hpp @@ -0,0 +1,57 @@ + +#ifndef WOLFF_MODELS_SYMMETRIC_H +#define WOLFF_MODELS_SYMMETRIC_H + +#include + +#include "potts.hpp" + +namespace wolff { + + template + class symmetric_t : public std::array { + public: + + symmetric_t() { + for (unsigned 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 (unsigned i = 0; i < q; i++) { + r_rot[i] = (*this)[r[i]]; + } + + return r_rot; + } + + potts_t act_inverse(const potts_t& s) const { + for (unsigned 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 (unsigned i = 0; i < q; i++) { + r_rot[(*this)[i]] = r[i]; + } + + return r_rot; + } + }; + +} + +#endif + diff --git a/lib/wolff_models/vector.hpp b/lib/wolff_models/vector.hpp new file mode 100644 index 0000000..4106e62 --- /dev/null +++ b/lib/wolff_models/vector.hpp @@ -0,0 +1,113 @@ + +#ifndef WOLFF_MODELS_VECTOR_H +#define WOLFF_MODELS_VECTOR_H + +#include +#include +#include + +namespace wolff { + + template + class vector_t : public std::array { + public: + + vector_t() { + this->fill((T)0); + (*this)[0] = (T)1; + } + + vector_t(const T *x) { + for (unsigned i = 0; i < q; i++) { + (*this)[i] = x[i]; + } + } + + typedef vector_t M_t; + typedef vector_t F_t; + + template + inline vector_t& operator+=(const vector_t &v) { + for (unsigned i = 0; i < q; i++) { + (*this)[i] += (T)v[i]; + } + return *this; + } + + template + inline vector_t& operator-=(const vector_t &v) { + for (unsigned i = 0; i < q; i++) { + (*this)[i] -= (T)v[i]; + } + return *this; + } + + inline vector_t operator*(unsigned x) const { + vector_t result; + for (unsigned i = 0; i < q; i++) { + result[i] = x * (*this)[i]; + } + + return result; + } + + inline vector_t operator*(double x) const { + vector_t result; + for (unsigned 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; + } + + inline T operator*(const vector_t& v) const { + double prod = 0; + + for (unsigned i = 0; i < q; i++) { + prod += v[i] * (*this)[i]; + } + + return prod; + } + + template + inline vector_t operator/(U a) const { + vector_t result; + for (unsigned i = 0; i < q; i++) { + result[i] = (*this)[i] / a; + } + + return result; + } + }; + + template + inline vector_t operator*(unsigned a, const vector_t&v) { + return v * a; + } + + template + inline vector_t operator*(double a, const vector_t&v) { + return v * a; + } + + template + std::ostream& operator<<(std::ostream& os, const vector_t&v) { + os << "( "; + for (T vi : v) { + os << vi << " "; + } + os << ")"; + return os; + } + +} + +#endif + -- cgit v1.2.3-70-g09d2