everything compiles now, and fixed a major bug in the graph code

7 files changed, 57 insertions, 85 deletions

diff --git a/lib/graph.cpp b/lib/graph.cpp
index 021873d..8c97274 100644
--- a/lib/graph.cpp
+++ b/lib/graph.cpp
@@ -20,7 +20,7 @@ graph_t::graph_t(D_t D, L_t L) {
 }
 
 void graph_t::add_ext() {
-  for (std::vector<v_t> v_adj_i : v_adj) {
+  for (std::vector<v_t>& v_adj_i : v_adj) {
     v_adj_i.push_back(nv);
   }
 
diff --git a/lib/height.h b/lib/height.h
index 81d3a2d..a4cbbed 100644
--- a/lib/height.h
+++ b/lib/height.h
@@ -32,50 +32,32 @@ struct height_t {
 
   typedef T M_t;
   typedef double F_t;
-};
-
-template <class T>
-void init(height_t<T> *ptr) {
-  ptr->x = (T)0;
-}
 
-template <class T>
-void free_spin(height_t<T> h) {
-  // do nothing!
-}
+  height_t() : x(0) {}
 
-template <class T>
-void free_spin(T h) {
-  // do nothing!
-}
+  height_t(T x) : x(x) {}
 
-void free_spin(double h) {
-  // do nothing!
-}
+  inline T operator*(v_t a) const {
+    return x * a;
+  }
 
-template <class T>
-height_t<T> copy(height_t<T> h) {
-  return h;
-}
+  inline double operator*(double a) const {
+    return x * a;
+  }
+};
 
 template <class T>
-void add(T *h1, int a, height_t<T> h2) {
-  (*h1) += a * h2.x;
-}
+inline T& operator+=(T& M, const height_t<T> &h) {
+  M += h.x;
 
-template <class T>
-void add(double *h1, double a, height_t<T> h2) {
-  (*h1) += a * h2.x;
+  return M;
 }
 
 template <class T>
-T scalar_multiple(int factor, height_t<T> h) {
-  return factor * h.x;
-}
+inline T& operator-=(T& M, const height_t<T> &h) {
+  M -= h.x;
 
-template <class T>
-double scalar_multiple(double factor, height_t<T> h) {
-  return factor * h.x;
+  return M;
 }
 
 double norm_squared(double h) {
diff --git a/lib/ising.h b/lib/ising.h
index 5932ecf..473a18a 100644
--- a/lib/ising.h
+++ b/lib/ising.h
@@ -53,7 +53,6 @@ class ising_t {
         return a;
       }
     }
-
 };
 
 inline int& operator+=(int& M, const ising_t &s) {
diff --git a/lib/potts.h b/lib/potts.h
index 3ab08bd..e6ea636 100644
--- a/lib/potts.h
+++ b/lib/potts.h
@@ -33,47 +33,40 @@ class potts_t {
 
     typedef vector_t<q, int> M_t;
     typedef vector_t<q, double> F_t;
-};
 
-template <q_t q>
-void init(potts_t <q> *p) {
-  p->x = 0;
-}
+    potts_t() : x(0) {}
 
-template <q_t q>
-void free_spin(potts_t <q> s) {
-  // do nothing!
-}
+    potts_t(q_t x) : x(x) {}
 
-template <q_t q>
-potts_t <q> copy(potts_t <q> s) {
-  return s;
-}
+    inline vector_t<q, int> operator*(v_t a) const {
+      vector_t<q, int> result;
+      result.fill(0);
+      result[x] = (int)a;
 
-template <q_t q>
-void add(vector_t<q, int> *s1, int a, potts_t <q> s2) {
-  (s1->x)[s2.x] += a;
-}
+      return result;
+    }
 
-template <q_t q>
-void add(vector_t<q, double> *s1, double a, potts_t <q> s2) {
-  (s1->x)[s2.x] += a;
-}
+    inline vector_t<q, double> operator*(double a) const {
+      vector_t<q, double> result;
+      result.fill(0.0);
+      result[x] = a;
+
+      return result;
+    }
+};
 
 template <q_t q>
-vector_t<q, int> scalar_multiple(int factor, potts_t <q> s) {
-  vector_t<q, int> M;
-  M.x = (int *)calloc(q, sizeof(int));
-  M.x[s.x] += factor;
+inline vector_t<q, int>& operator+=(vector_t<q, int>& M, const potts_t<q> &s) {
+  M[s.x]++;
+
   return M;
 }
 
 template <q_t q>
-vector_t<q, double> scalar_multiple(double factor, potts_t <q> s) {
-  vector_t<q, double> F;
-  F.x = (double *)calloc(q, sizeof(double));
-  F.x[s.x] += factor;
-  return F;
+inline vector_t<q, int>& operator-=(vector_t<q, int>& M, const potts_t<q> &s) {
+  M[s.x]--;
+
+  return M;
 }
 
 // we could inherit norm_squared from vector.h, but convention dictates that
@@ -81,8 +74,8 @@ vector_t<q, double> scalar_multiple(double factor, potts_t <q> s) {
 template <q_t q>
 double norm_squared(vector_t<q, double> s) {
   double total = 0;
-  for (q_t i = 0; i < q; i++) {
-    total += pow(s.x[i], 2);
+  for (double& x : s) {
+    total += pow(x, 2);
   }
 
   return total * (double)q / ((double)q - 1.0);
@@ -92,7 +85,9 @@ double norm_squared(vector_t<q, double> s) {
 // to nv we don't need to write the last element
 template <q_t q>
 void write_magnetization(vector_t<q, int> M, FILE *outfile) {
-  fwrite(M.x, sizeof(int), q - 1, outfile);
+  for (int& x : M) {
+    fwrite(&x, sizeof(int), q - 1, outfile);
+  }
 }
 
 // knock yourself out
diff --git a/lib/state.h b/lib/state.h
index 360e8f8..b57a9e2 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -21,11 +21,11 @@ class state_t {
     double E;
     typename X_t::M_t M; // the "sum" of the spins, like the total magnetization
     v_t last_cluster_size;
-    typename X_t::F_t *ReF;
-    typename X_t::F_t *ImF;
+    std::vector<typename X_t::F_t> ReF;
+    std::vector<typename X_t::F_t> ImF;
     // updating fourier terms F requires many cos and sin calls, faster to do it beforehand.
-    double *precomputed_cos;
-    double *precomputed_sin;
+    std::vector<double> precomputed_cos;
+    std::vector<double> precomputed_sin;
 
     std::function <double(X_t, X_t)> J;
     std::function <double(X_t)> H;
@@ -39,14 +39,14 @@ class state_t {
       E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(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));
+      ReF.resize(D);
+      ImF.resize(D);
       for (D_t i = 0; i < D; i++) {
         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));
+      precomputed_cos.resize(L);
+      precomputed_sin.resize(L);
       for (L_t i = 0; i < L; i++) {
         precomputed_cos[i] = cos(2 * M_PI * (double)i / (double)L);
         precomputed_sin[i] = sin(2 * M_PI * (double)i / (double)L);
@@ -55,10 +55,6 @@ class state_t {
 
     ~state_t() {
       free_spin(R);
-      free(ReF);
-      free(ImF);
-      free(precomputed_sin);
-      free(precomputed_cos);
     }
 };
 
diff --git a/lib/vector.h b/lib/vector.h
index d5ced03..beee1a7 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -81,10 +81,10 @@ class vector_t : public std::array<T, q> {
 };
 
 
-template<q_t q>
-double norm_squared(vector_t<q, double> v) {
+template<q_t q, class T>
+double norm_squared(vector_t<q, T> v) {
   double tmp = 0;
-  for (double &x : v) {
+  for (T &x : v) {
     tmp += pow(x, 2);
   }
 
diff --git a/src/wolff_clock.cpp b/src/wolff_clock.cpp
index 86badfe..e186c44 100644
--- a/src/wolff_clock.cpp
+++ b/src/wolff_clock.cpp
@@ -107,7 +107,7 @@ int main(int argc, char *argv[]) {
   if (!draw) {
     // a very simple example: measure the average magnetization
     measurement = [&] (const sim_t *s) {
-      average_M += (double)s->M.x[0] / (double)N / (double)s->nv;
+      average_M += (double)s->M[0] / (double)N / (double)s->nv;
     };
   } else {
     // a more complex example: measure the average magnetization, and draw the spin configuration to the screen
@@ -124,7 +124,7 @@ int main(int argc, char *argv[]) {
     gluOrtho2D(0.0, L, 0.0, L);
 
     measurement = [&] (const sim_t *s) {
-      average_M += (double)s->M.x[0] / (double)N / (double)s->nv;
+      average_M += (double)s->M[0] / (double)N / (double)s->nv;
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(L, 2); i++) {
         potts_t<POTTSQ> tmp_s = act_inverse(s->R, s->spins[i]);