changed the way that the finite states system works

11 files changed, 72 insertions, 74 deletions

diff --git a/doc/finite_states.rst b/doc/finite_states.rst
index 26178cd..c36881e 100644
--- a/doc/finite_states.rst
+++ b/doc/finite_states.rst
@@ -5,20 +5,24 @@
 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 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.
+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
 ====================
 
+.. c:macro:: WOLFF_USE_FINITE_STATES
+
+   This macro must defined before :file:`wolff.hpp` or any of the other header files are invoked.
+
 .. 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.
+   This macro must be defined and given a value equal to the number of states your model can take. It must be defined before :file:`wolff.hpp` or any of the other header files are invoked.
 
-.. cpp:var:: const X_t finite_states_possible[WOLFF_FINITE_STATES_N]
+.. cpp:function:: X_t::X_t(q_t)
 
-   You must supply a constant array which lists each of the possible states of your individual spins. 
+   Your spin class :cpp:class:`X_t` must have a constructor defined that takes a :cpp:type:`q_t` and returns a unique state for all arguments less than :c:macro:`WOLFF_FINITE_STATES_N`.
 
-.. cpp:function:: q_t finite_states_enum(const X_t&)
+.. cpp:function:: q_t X_t::enumerate()
 
-   You must define a function which takes each state and returns its index in :cpp:var:`finite_states_possible`.
+   Your spin class :cpp:class:`X_t` must have a function defined that returns the index associated with a given state. This must be the inverse function of the constructor above.
 
diff --git a/examples/clock.cpp b/examples/clock.cpp
index 4b0ffb8..0cba8f6 100644
--- a/examples/clock.cpp
+++ b/examples/clock.cpp
@@ -3,11 +3,13 @@
 #include <iostream>
 #include <chrono>
 
-#include "simple_measurement.hpp"
+#define WOLFF_USE_FINITE_STATES
+#define WOLFF_FINITE_STATES_N WOLFF_POTTSQ
 
 #include <wolff/models/potts.hpp>
 #include <wolff/models/dihedral.hpp>
-#include <wolff/finite_states.hpp>
+
+#include "simple_measurement.hpp"
 
 #include <wolff.hpp>
 
diff --git a/examples/ising.cpp b/examples/ising.cpp
index fe21279..f187081 100644
--- a/examples/ising.cpp
+++ b/examples/ising.cpp
@@ -3,13 +3,14 @@
 #include <iostream>
 #include <chrono>
 
-#include "simple_measurement.hpp"
+#define WOLFF_USE_FINITE_STATES
 
 #include <wolff/models/ising.hpp>
-#include <wolff/finite_states.hpp>
 
 #include <wolff.hpp>
 
+#include "simple_measurement.hpp"
+
 using namespace wolff;
 
 int main(int argc, char *argv[]) {
diff --git a/examples/ising_animation.cpp b/examples/ising_animation.cpp
index 2104e91..ab10585 100644
--- a/examples/ising_animation.cpp
+++ b/examples/ising_animation.cpp
@@ -5,8 +5,9 @@
 
 #include <GL/glut.h>
 
+#define WOLFF_USE_FINITE_STATES
+
 #include <wolff/models/ising.hpp>
-#include <wolff/finite_states.hpp>
 
 #include <wolff.hpp>
 
diff --git a/examples/ising_no_field.cpp b/examples/ising_no_field.cpp
index eebd672..28d7c70 100644
--- a/examples/ising_no_field.cpp
+++ b/examples/ising_no_field.cpp
@@ -4,13 +4,15 @@
 #include <chrono>
 
 #define WOLFF_NO_FIELD
-#include "simple_measurement.hpp"
+#define WOLFF_USE_FINITE_STATES
 
 #include <wolff/models/ising.hpp>
-#include <wolff/finite_states.hpp>
 
 #include <wolff.hpp>
 
+#include "simple_measurement.hpp"
+
+
 using namespace wolff;
 
 int main(int argc, char *argv[]) {
diff --git a/examples/potts.cpp b/examples/potts.cpp
index 84494e2..24be2b3 100644
--- a/examples/potts.cpp
+++ b/examples/potts.cpp
@@ -3,11 +3,13 @@
 #include <iostream>
 #include <chrono>
 
-#include "simple_measurement.hpp"
+#define WOLFF_USE_FINITE_STATES
+#define WOLFF_FINITE_STATES_N WOLFF_POTTSQ
 
 #include <wolff/models/potts.hpp>
 #include <wolff/models/symmetric.hpp>
-#include <wolff/finite_states.hpp>
+
+#include "simple_measurement.hpp"
 
 #include <wolff.hpp>
 
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
index b66f367..3912446 100644
--- a/lib/include/wolff/cluster.hpp
+++ b/lib/include/wolff/cluster.hpp
@@ -68,9 +68,8 @@ void system<R_t, X_t>::flip_cluster(v_t i0, const R_t& r,
           dE = B(s0s_old) - B(s0s_new);
 #endif
 
-#ifdef WOLFF_FINITE_STATES
-          p = finite_states_Bp[finite_states_enum(s0s_old)]
-                              [finite_states_enum(s0s_new)];
+#ifdef WOLFF_USE_FINITE_STATES
+          p = Bp[s0s_old.enumerate()][s0s_new.enumerate()];
 #endif
 
           // run measurement hooks for encountering a ghost bond
@@ -84,17 +83,15 @@ void system<R_t, X_t>::flip_cluster(v_t i0, const R_t& r,
           dE = Z(s[i], s[j]) - Z(si_new, s[j]);
 #endif
 
-#ifdef WOLFF_FINITE_STATES
-          p = finite_states_Zp[finite_states_enum(s[i])]
-                              [finite_states_enum(si_new)]
-                              [finite_states_enum(s[j])];
+#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 FINITE_STATES
+#ifndef WOLFF_USE_FINITE_STATES
         p = 1.0 - exp(-dE / T);
 #endif
 
diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp
deleted file mode 100644
index 7e54eff..0000000
--- a/lib/include/wolff/finite_states.hpp
+++ /dev/null
@@ -1,39 +0,0 @@
-
-
-#ifndef WOLFF_FINITE_STATES
-#define WOLFF_FINITE_STATES
-
-#include <cmath>
-#include <array>
-
-#include "system.hpp"
-
-namespace wolff {
-
-std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Zp;
-#ifndef WOLFF_NO_FIELD
-std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Bp;
-#endif
-
-template <class R_t, class X_t>
-void finite_states_init(const system<R_t, X_t>& 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++) {
-      finite_states_Bp[i][j] = 1.0 - exp(-(S.B(finite_states_possible[i]) - S.B(finite_states_possible[j])) / S.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++) {
-        finite_states_Zp[i][j][k] = 1.0 - exp(-(S.Z(finite_states_possible[i], finite_states_possible[k]) - S.Z(finite_states_possible[j], finite_states_possible[k])) / S.T);
-      }
-    }
-  }
-}
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp
index 824c063..39d0ad6 100644
--- a/lib/include/wolff/models/ising.hpp
+++ b/lib/include/wolff/models/ising.hpp
@@ -2,6 +2,8 @@
 #ifndef WOLFF_MODELS_ISING
 #define WOLFF_MODELS_ISING
 
+#define WOLFF_FINITE_STATES_N 2
+
 #include "../types.h"
 #include "../system.hpp"
 
@@ -14,7 +16,7 @@ class ising_t {
     ising_t() : x(false) {}
 
     ising_t(bool x) : x(x) {} 
-    ising_t(int x) : x((bool)x) {}
+    ising_t(q_t x) : x((bool)x) {}
 
     ising_t act(const ising_t& s) const {
       if (x) {
@@ -66,6 +68,10 @@ class ising_t {
         return -1;
       }
     }
+
+    q_t enumerate() const {
+      return (q_t)x;
+    }
 };
 
 inline ising_t::M_t operator*(v_t a, const ising_t& s) {
@@ -80,10 +86,6 @@ ising_t gen_ising(std::mt19937&, const system<ising_t, ising_t>&, 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(const ising_t& state) { return (q_t)state.x; }
-
 }
 
 #endif
diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp
index ce419ce..3d0a0e9 100644
--- a/lib/include/wolff/models/potts.hpp
+++ b/lib/include/wolff/models/potts.hpp
@@ -44,6 +44,10 @@ class potts_t {
 
       return result;
     }
+
+    q_t enumerate() const {
+      return x;
+    }
 };
 
 template<q_t q>
@@ -56,10 +60,5 @@ inline typename potts_t<q>::F_t operator*(double a, const potts_t<q>& s) {
   return s * a;
 }
 
-#define WOLFF_FINITE_STATES_N WOLFF_POTTSQ
-const potts_t<WOLFF_POTTSQ> finite_states_possible[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 q>
-q_t finite_states_enum(const potts_t<q>& state) { return (q_t)state.x; }
-
 }
 
diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp
index 32ad38e..4aa8ee4 100644
--- a/lib/include/wolff/system.hpp
+++ b/lib/include/wolff/system.hpp
@@ -39,6 +39,13 @@ class system {
 #endif
 #endif
 
+#ifdef WOLFF_USE_FINITE_STATES
+    std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp;
+#ifndef WOLFF_NO_FIELD
+    std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Bp;
+#endif
+#endif
+
     system(graph g, double T,
 #ifdef WOLFF_BOND_DEPENDENCE
         std::function <double(v_t, const X_t&, v_t, const X_t&)> Z
@@ -63,13 +70,33 @@ class system {
 #ifndef WOLFF_NO_FIELD
       G.add_ghost();
 #endif
-#ifdef WOLFF_FINITE_STATES
-      finite_states_init(*this);
+#ifdef WOLFF_USE_FINITE_STATES
+      this->finite_states_init();
 #endif
     }
 
     void flip_cluster(v_t, const R_t&, std::mt19937&, measurement<R_t, X_t>&);
     void run_wolff(N_t, std::function <R_t(std::mt19937&, const system<R_t, X_t>&, v_t)> r_gen, measurement<R_t, X_t>& 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
+
 };
 
 }