8 files changed, 699 insertions, 0 deletions
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 <unsigned q>
+  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<q> act(const potts_t<q>& s) const {
+        if (this->is_reflection) {
+          return potts_t<q>(((q + this->x) - s.x) % q);
+        } else {
+          return potts_t<q>((this->x + s.x) % q);
+        }
+      }
+
+      dihedral_t<q> act(dihedral_t<q> r) const {
+        if (this->is_reflection) {
+          return dihedral_t<q>(!(r.is_reflection), ((q + this->x) - r.x) % q); 
+        } else {
+          return dihedral_t<q>(r.is_reflection, (this->x + r.x) % q);
+        }
+      }
+
+      potts_t<q> act_inverse(potts_t<q> s) const {
+        if (this->is_reflection) {
+          return this->act(s);
+        } else {
+          return potts_t<q>(((s.x + q) - this->x) % q);
+        }
+      }
+
+      dihedral_t<q> act_inverse(dihedral_t<q> r) const {
+        if (this->is_reflection) {
+          return this->act(r);
+        } else {
+          return dihedral_t<q>(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 <cmath>
+#include "height.hpp"
+
+namespace wolff {
+
+  template <class T>
+  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<T> act(const height_t<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<T> act(const dihedral_inf_t<T>& r) const {
+      if (this->is_reflection) {
+        return dihedral_inf_t<T>(!r.is_reflection, this->x - r.x);
+      } else {
+        return dihedral_inf_t<T>(r.is_reflection, this->x + r.x);
+      }
+    }
+
+    height_t<T> act_inverse(const height_t<T>& h) const {
+      if (this->is_reflection) {
+        return this->act(h);
+      } else {
+        return height_t(h.x - this->x);
+      }
+    }
+
+    dihedral_inf_t<T> act_inverse(const dihedral_inf_t<T>& r) const {
+      if (this->is_reflection) {
+        return this->act(r);
+      } else {
+        return dihedral_inf_t<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 <cmath>
+
+namespace wolff {
+
+  template <class T>
+  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 <class T>
+  inline typename height_t<T>::M_t operator*(unsigned a, height_t<T> h) {
+    return h * a;
+  }
+
+  template <class T>
+  inline typename height_t<T>::F_t operator*(double a, height_t<T> h) {
+    return h * a;
+  }
+
+  template <class T>
+  inline T& operator+=(T& M, const height_t<T> &h) {
+    M += h.x;
+
+    return M;
+  }
+
+  template <class T>
+  inline T& operator-=(T& M, const height_t<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 <class G_t>
+  ising_t gen_ising(std::mt19937&, const system<ising_t, ising_t, G_t>&, 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 <random>
+#include <cmath>
+
+#include <wolff.hpp>
+#include "vector.hpp"
+
+namespace wolff {
+
+  template <unsigned q, class T>
+  class orthogonal_t : public std::array<std::array<T, q>, 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<q, T> act(const vector_t <q, T>& v) const {
+      vector_t <q, 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<q, T> act(const orthogonal_t <q, T>& m) const {
+      orthogonal_t <q, 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 <q, T> act_inverse(const vector_t <q, T>& v) const {
+      if (is_reflection) {
+        return this->act(v); // reflections are their own inverse
+      } else {
+        vector_t <q, 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 <q, T> act_inverse(const orthogonal_t <q, T>& m) const {
+      if (is_reflection) {
+        return this->act(m); // reflections are their own inverse
+      } else {
+        orthogonal_t <q, 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 <unsigned q, class G_t>
+  orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const system<orthogonal_t<q, double>, vector_t<q, double>, G_t>&, const typename G_t::vertex&) {
+    std::normal_distribution<double> dist(0.0,1.0);
+    orthogonal_t <q, double> 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 <unsigned q, class G_t>
+  orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const system<orthogonal_t<q, double>, vector_t<q, double>, G_t>& S, const typename G_t::vertex& v0, double epsilon, unsigned int n) {
+    std::normal_distribution<double> dist(0.0,1.0);
+    orthogonal_t <q, double> m;
+    m.is_reflection = true;
+
+    vector_t <q, double> v;
+
+    if (n > 1) {
+      std::uniform_int_distribution<unsigned int> 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 <cmath>
+
+#include "vector.hpp"
+
+namespace wolff {
+
+  template <unsigned q>
+  class potts_t {
+    public:
+      unsigned x;
+
+      potts_t() : x(0) {}
+      potts_t(unsigned x) : x(x) {}
+
+      typedef vector_t<q, int> M_t;
+      typedef vector_t<q, double> F_t;
+
+      inline vector_t<q, int> operator*(unsigned a) const {
+        vector_t<q, int> result;
+        result.fill(0);
+        result[x] = (int)a;
+
+        return result;
+      }
+
+      inline vector_t<q, double> operator*(double a) const {
+        vector_t<q, double> result;
+        result.fill(0.0);
+        result[x] = a;
+
+        return result;
+      }
+
+      inline vector_t<q, int> operator-(const potts_t<q> &s) const {
+        vector_t<q, int> result;
+        result.fill(0);
+
+        result[x]++;
+        result[s.x]--;
+
+        return result;
+      }
+
+      unsigned enumerate() const {
+        return x;
+      }
+  };
+
+  template<unsigned q>
+  inline typename potts_t<q>::M_t operator*(unsigned a, const potts_t<q>& s) {
+    return s * a;
+  }
+
+  template<unsigned q>
+  inline typename potts_t<q>::F_t operator*(double a, const potts_t<q>& 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 <array>
+
+#include "potts.hpp"
+
+namespace wolff {
+
+  template <unsigned q>
+  class symmetric_t : public std::array<unsigned, q> {
+    public:
+
+      symmetric_t() {
+        for (unsigned i = 0; i < q; i++) {
+          (*this)[i] = i;
+        }
+      }
+
+      potts_t<q> act(const potts_t<q> &s) const {
+        return potts_t<q>((*this)[s.x]);
+      }
+
+      symmetric_t<q> act(const symmetric_t<q>& r) const {
+        symmetric_t<q> r_rot;
+        for (unsigned i = 0; i < q; i++) {
+          r_rot[i] = (*this)[r[i]];
+        }
+
+        return r_rot;
+      }
+
+      potts_t<q> act_inverse(const potts_t<q>& s) const {
+        for (unsigned i = 0; i < q; i++) {
+          if ((*this)[i] == s.x) {
+            return potts_t<q>(i);
+          }
+        }
+
+        exit(EXIT_FAILURE);
+      }
+
+      symmetric_t<q> act_inverse(const symmetric_t<q>& r) const {
+        symmetric_t<q> 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 <cmath>
+#include <array>
+#include <iostream>
+
+namespace wolff {
+
+  template <unsigned q, class T>
+  class vector_t : public std::array<T, q> {
+    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 <q, T> M_t;
+      typedef vector_t <q, double> F_t;
+
+      template <class U>
+      inline vector_t<q, T>& operator+=(const vector_t<q, U> &v) {
+        for (unsigned i = 0; i < q; i++) {
+          (*this)[i] += (T)v[i];
+        }
+        return *this;
+      }
+
+      template <class U>
+      inline vector_t<q, T>& operator-=(const vector_t<q, U> &v) {
+        for (unsigned i = 0; i < q; i++) {
+          (*this)[i] -= (T)v[i];
+        }
+        return *this;
+      }
+
+      inline vector_t<q, T> operator*(unsigned x) const {
+        vector_t<q, T> result;
+        for (unsigned i = 0; i < q; i++) {
+          result[i] = x * (*this)[i];
+        }
+
+        return result;
+      }
+
+      inline vector_t<q, double> operator*(double x) const {
+        vector_t<q, double> result;
+        for (unsigned i = 0; i < q; i++) {
+          result[i] = x * (*this)[i];
+        }
+
+        return result;
+      }
+
+      inline vector_t<q, T> operator-(const vector_t<q, T>& v) const {
+        vector_t<q, T> diff = *this;
+        diff -= v;
+        return diff;
+      }
+
+      inline T operator*(const vector_t<q, T>& v) const {
+        double prod = 0;
+
+        for (unsigned i = 0; i < q; i++) {
+          prod += v[i] * (*this)[i];
+        }
+
+        return prod;
+      }
+
+      template <class U>
+      inline vector_t<q, T> operator/(U a) const {
+        vector_t<q, T> result;
+        for (unsigned i = 0; i < q; i++) {
+          result[i] = (*this)[i] / a;
+        }
+
+        return result;
+      }
+  };
+
+  template<unsigned q, class T>
+  inline vector_t<q, T> operator*(unsigned a, const vector_t<q, T>&v) {
+    return v * a;
+  }
+
+  template<unsigned q, class T>
+  inline vector_t<q, double> operator*(double a, const vector_t<q, T>&v) {
+    return v * a;
+  }
+
+  template<unsigned q, class T>
+  std::ostream& operator<<(std::ostream& os, const vector_t<q, T>&v) {
+    os << "( ";
+    for (T vi : v) {
+      os << vi << " ";
+    }
+    os << ")";
+    return os;
+  }
+
+}
+
+#endif
+