partially class-ified, ising and On work but potts and height do not

6 files changed, 110 insertions, 149 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index 8aae02a..427c3c8 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -68,22 +68,20 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
           prob = 1.0 - exp(-dE / state->T);
 #endif
 
-          add(&(state->M), -1, rs_old);
-          add(&(state->M),  1, rs_new);
+          state->M -= rs_old;
+          state->M += rs_new;
+
           state->E += dE;
 
           for (D_t i = 0; i < state->D; i++) {
             L_t x = (non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L;
 
-            add(&(state->ReF[i]), -state->precomputed_cos[x], rs_old);
-            add(&(state->ReF[i]),  state->precomputed_cos[x], rs_new);
+            state->ReF[i] -= rs_old * state->precomputed_cos[x];
+            state->ReF[i] += rs_new * state->precomputed_cos[x];
 
-            add(&(state->ImF[i]), -state->precomputed_sin[x], rs_old);
-            add(&(state->ImF[i]),  state->precomputed_sin[x], rs_new);
+            state->ImF[i] -= rs_old * state->precomputed_sin[x];
+            state->ImF[i] += rs_new * state->precomputed_sin[x];
           }
-
-          free_spin (rs_old);
-          free_spin (rs_new);
         } else {
           double dE = state->J(si_old, sj) - state->J(si_new, sj);
 #ifdef FINITE_STATES
@@ -103,7 +101,6 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
         free_spin(state->R);
         state->R = R_new;
       } else {
-        free_spin(state->spins[v]);
         state->spins[v] = si_new;
       }
 
diff --git a/lib/correlation.h b/lib/correlation.h
index 26c3a99..49c6ff2 100644
--- a/lib/correlation.h
+++ b/lib/correlation.h
@@ -11,7 +11,7 @@ double correlation_length(const state_t <R_t, X_t> *s) {
   double total = 0;
 
   for (D_t j = 0; j < s->D; j++) {
-    total += norm_squared(s->ReF[j]) + norm_squared(s->ImF[j]);
+    total += norm_squared(s->ReF[j]) + norm_squared(s->ReF[j]);
   }
 
   return total / s->D;
diff --git a/lib/ising.h b/lib/ising.h
index f08ae58..5932ecf 100644
--- a/lib/ising.h
+++ b/lib/ising.h
@@ -31,58 +31,49 @@ class ising_t {
 
     typedef int M_t;
     typedef double F_t;
-};
 
-void init(ising_t *p) {
-  p->x = false;
-}
+    ising_t() {
+      x = false;
+    }
 
-void free_spin(ising_t s) {
-  // do nothing!
-}
+    ising_t(bool x) : x(x) {}
 
-void free_spin(int s) {
-  // do nothing
-}
+    inline int operator*(v_t a) const {
+      if (x) {
+        return -(int)a;
+      } else {
+        return (int)a;
+      }
+    }
 
-void free_spin(double s) {
-  // do nothing
-}
+    inline double operator*(double a) const {
+      if (x) {
+        return -a;
+      } else {
+        return a;
+      }
+    }
 
-ising_t copy(ising_t s) {
-  return s;
-}
+};
 
-void add(int *s1, int a, ising_t s2) {
-  if (s2.x) {
-    *s1 -= a;
+inline int& operator+=(int& M, const ising_t &s) {
+  if (s.x) {
+    M--;
   } else {
-    *s1 += a;
+    M++;
   }
-}
 
-void add(double *s1, double a, ising_t s2) {
-  if (s2.x) {
-    *s1 -= a;
-  } else {
-    *s1 += a;
-  }
+  return M;
 }
 
-int scalar_multiple(int factor, ising_t s) {
+inline int& operator-=(int& M, const ising_t &s) {
   if (s.x) {
-    return -factor;
+    M++;
   } else {
-    return factor;
+    M--;
   }
-}
 
-double scalar_multiple(double factor, ising_t s) {
-  if (s.x) {
-    return -factor;
-  } else {
-    return factor;
-  }
+  return M;
 }
 
 double norm_squared(double s) {
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 2a391ae..c0958f2 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -53,20 +53,19 @@ void free_spin (orthogonal_t <q, T> m) {
 template <q_t q, class T>
 vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) {
   vector_t <q, T> v_rot;
-  v_rot.x = (T *)calloc(q, sizeof(T));
 
   if (m.is_reflection) {
     double prod = 0;
     for (q_t i = 0; i < q; i++) {
-      prod += v.x[i] * m.x[i];
+      prod += v[i] * m.x[i];
     }
     for (q_t i = 0; i < q; i++) {
-      v_rot.x[i] = v.x[i] - 2 * prod * m.x[i];
+      v_rot[i] = v[i] - 2 * prod * m.x[i];
     }
   } else {
     for (q_t i = 0; i < q; i++) {
       for (q_t j = 0; j < q; j++) {
-        v_rot.x[i] += m.x[q * i + j] * v.x[j];
+        v_rot[i] += m.x[q * i + j] * v[j];
       }
     }
   }
@@ -112,11 +111,10 @@ vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) {
     return act(m, v); // reflections are their own inverse
   } else {
     vector_t <q, T> v_rot;
-    v_rot.x = (T *)calloc(q, sizeof(T));
 
     for (q_t i = 0; i < q; i++) {
       for (q_t j = 0; j < q; j++) {
-        v_rot.x[i] += m.x[q * j + i] * v.x[j];
+        v_rot[i] += m.x[q * j + i] * v[j];
       }
     }
 
@@ -176,33 +174,32 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q
 
   if (n > 1) {
     unsigned int rotation = gsl_rng_uniform_int(r, n);
-    v.x = (double *)malloc(q * sizeof(double));
 
     double cosr = cos(2 * M_PI * rotation / (double)n / 2.0);
     double sinr = sin(2 * M_PI * rotation / (double)n / 2.0);
 
-    v.x[0] = v0.x[0] * cosr - v0.x[1] * sinr;
-    v.x[1] = v0.x[1] * cosr + v0.x[0] * sinr;
+    v[0] = v0[0] * cosr - v0[1] * sinr;
+    v[1] = v0[1] * cosr + v0[0] * sinr;
 
     for (q_t i = 2; i < q; i++) {
-      v.x[i] = v0.x[i];
+      v[i] = v0[i];
     }
   } else {
-    v.x = v0.x;
+    v = v0;
   }
 
   double m_dot_v = 0;
 
   for (q_t i = 0; i < q; i++) {
     m.x[i] = gsl_ran_ugaussian(r);
-    m_dot_v += m.x[i] * v.x[i];
+    m_dot_v += m.x[i] * v[i];
   }
 
   double v2 = 0;
   double factor = epsilon * gsl_ran_ugaussian(r);
 
   for (q_t i = 0; i < q; i++) {
-    m.x[i] = m.x[i] - m_dot_v * v.x[i] + factor * v.x[i];
+    m.x[i] = m.x[i] - m_dot_v * v[i] + factor * v[i];
     v2 += pow(m.x[i], 2);
   }
 
@@ -212,10 +209,6 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q
     m.x[i] /= mag_v;
   }
 
-  if (n > 1) {
-    free(v.x);
-  }
-
   return m;
 }
 
diff --git a/lib/state.h b/lib/state.h
index 3cef157..360e8f8 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -2,6 +2,7 @@
 #pragma once
 
 #include <functional>
+#include <vector>
 
 #include "types.h"
 #include "graph.h"
@@ -15,7 +16,7 @@ class state_t {
     v_t ne;
     graph_t g;
     double T;
-    X_t *spins;
+    std::vector<X_t> spins;
     R_t R;
     double E;
     typename X_t::M_t M; // the "sum" of the spins, like the total magnetization
@@ -29,23 +30,20 @@ class state_t {
     std::function <double(X_t, X_t)> J;
     std::function <double(X_t)> H;
 
-    state_t(D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H) : D(D), L(L), T(T), J(J), H(H), g(D, L) {
+    state_t(D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H) : D(D), L(L), g(D, L), T(T), J(J), H(H) {
       nv = g.nv;
       ne = g.ne;
       g.add_ext();
-      spins = (X_t *)malloc(nv * sizeof(X_t));
-      for (v_t i = 0; i < nv; i++) {
-        init (&(spins[i]));
-      }
+      spins.resize(nv);
       init (&R);
       E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
-      M = scalar_multiple((int)nv, spins[0]);
+      M = spins[0] * nv;
       last_cluster_size = 0;
       ReF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t));
       ImF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t));
       for (D_t i = 0; i < D; i++) {
-        ReF[i] = scalar_multiple(0.0, spins[0]);
-        ImF[i] = scalar_multiple(0.0, spins[0]);
+        ReF[i] = spins[0] * 0.0;
+        ImF[i] = spins[0] * 0.0;
       }
       precomputed_cos = (double *)malloc(L * sizeof(double));
       precomputed_sin = (double *)malloc(L * sizeof(double));
@@ -56,16 +54,7 @@ class state_t {
     }
 
     ~state_t() {
-      for (v_t i = 0; i < nv; i++) {
-        free_spin(spins[i]);
-      }
-      free(spins);
       free_spin(R);
-      free_spin(M);
-      for (D_t i = 0; i < D; i++) {
-        free_spin(ReF[i]);
-        free_spin(ImF[i]);
-      }
       free(ReF);
       free(ImF);
       free(precomputed_sin);
diff --git a/lib/vector.h b/lib/vector.h
index 2fe6ab8..d5ced03 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -3,6 +3,7 @@
 
 #include <stdlib.h>
 #include <cmath>
+#include <array>
 
 #include "types.h"
 
@@ -24,78 +25,67 @@
  */
 
 template <q_t q, class T>
-class vector_t {
+class vector_t : public std::array<T, q> {
   public: 
-    T *x;
 
     // M_t needs to hold the sum of nv spins
     typedef vector_t <q, T> M_t;
 
     // F_t needs to hold the double-weighted sum of spins
     typedef vector_t <q, double> F_t;
-};
-
-template <q_t q, class T>
-void init(vector_t <q, T> *ptr) {
-  ptr->x = (T *)calloc(q, sizeof(T));
-
-  // initialize a unit vector
-  ptr->x[0] = (T)1;
-}
-
-template <q_t q, class T>
-void free_spin (vector_t <q, T> v) {
-  free(v.x);
-}
-
-template <q_t q, class T>
-vector_t <q, T> copy (vector_t <q, T> v) {
-  vector_t <q, T> v_copy;
- 
- v_copy.x = (T *)calloc(q, sizeof(T));
 
-  for (q_t i = 0; i < q; i++) {
-    v_copy.x[i] = v.x[i];
-  }
-
-  return v_copy;
-}
-
-template <q_t q, class T, class U, class V>
-void add(vector_t<q, U> *v1, V a, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] += (U)(a * v2.x[i]);
-  }
-}
-
-template <q_t q, class T>
-vector_t <q, T> scalar_multiple(int a, vector_t <q, T> v) {
-  vector_t <q, T> multiple;
-  multiple.x = (T *)malloc(q * sizeof(T));
-  for (q_t i = 0; i < q; i++) {
-    multiple.x[i] = a * v.x[i];
-  }
-
-  return multiple;
-}
-
-template <q_t q, class T>
-vector_t <q, T> scalar_multiple(double a, vector_t <q, T> v) {
-  vector_t <q, T> multiple;
-  multiple.x = (T *)malloc(q * sizeof(T));
-  for (q_t i = 0; i < q; i++) {
-    multiple.x[i] = a * v.x[i];
-  }
+    vector_t() {
+      this->fill((T)0);
+      (*this)[1] = (T)1;
+    }
+
+    vector_t(const T *x) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] = x[i];
+      }
+    }
+
+    template <class U>
+    inline vector_t<q, T>& operator+=(const vector_t<q, U> &v) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] += (U)v[i];
+      }
+      return *this;
+    }
+
+    template <class U>
+    inline vector_t<q, T>& operator-=(const vector_t<q, U> &v) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] -= (U)v[i];
+      }
+      return *this;
+    }
+
+    inline vector_t<q, T> operator*(v_t x) const {
+      vector_t<q, T> result;
+      for (q_t 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 (q_t i = 0; i < q; i++) {
+        result[i] = x * (*this)[i];
+      }
+
+      return result;
+    }
+};
 
-  return multiple;
-}
 
-template <q_t q, class T>
-double norm_squared (vector_t <q, T> v) {
+template<q_t q>
+double norm_squared(vector_t<q, double> v) {
   double tmp = 0;
-
-  for (q_t i = 0; i < q; i++) {
-    tmp += pow(v.x[i], 2);
+  for (double &x : v) {
+    tmp += pow(x, 2);
   }
 
   return tmp;
@@ -103,7 +93,9 @@ double norm_squared (vector_t <q, T> v) {
 
 template <q_t q, class T>
 void write_magnetization(vector_t <q, T> M, FILE *outfile) {
-  fwrite(M.x, sizeof(T), q, outfile);
+  for (q_t i = 0; i < q; i++) {
+    fwrite(&(M[i]), sizeof(T), q, outfile);
+  }
 }
 
 // below functions and definitions are unnecessary for wolff.h but useful.
@@ -111,7 +103,7 @@ void write_magnetization(vector_t <q, T> M, FILE *outfile) {
 template <q_t q> // save some space and don't write whole doubles
 void write_magnetization(vector_t <q, double> M, FILE *outfile) {
   for (q_t i = 0; i < q; i++) {
-    float M_tmp = (float)M.x[i];
+    float M_tmp = (float)M[i];
     fwrite(&M_tmp, sizeof(float), 1, outfile);
   }
 }
@@ -121,7 +113,7 @@ T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
   T prod = 0;
 
   for (q_t i = 0; i < q; i++) {
-    prod += v1.x[i] * v2.x[i];
+    prod += v1[i] * v2[i];
   }
 
   return prod;
@@ -129,8 +121,7 @@ T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
 
 template <q_t q, class T>
 double H_vector(vector_t <q, T> v1, T *H) {
-  vector_t <q, T> H_vec;
-  H_vec.x = H;
+  vector_t <q, T> H_vec(H);
   return (double)(dot <q, T> (v1, H_vec));
 }