now using const passing of object reference instead of pointers

5 files changed, 53 insertions, 50 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index f948586..f3271b0 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -12,45 +12,45 @@
 #include "graph.h"
 
 template <class R_t, class X_t>
-void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
+void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, gsl_rng *rand) {
   v_t nv = 0;
 
   std::stack<v_t> stack;
   stack.push(v0);
 
-  std::vector<bool> marks(state->g.nv, false);
+  std::vector<bool> marks(s.g.nv, false);
 
   while (!stack.empty()) {
     v_t v = stack.top();
     stack.pop();
 
-    if (!marks[v]) {
+    if (!marks[v]) { // don't reprocess anyone we've already visited!
       X_t si_old, si_new;
       R_t R_old, R_new;
 
-      R_old = state->R;
+      R_old = s.R;
       marks[v] = true;
 
-      bool v_is_ghost = (v == state->nv);
+      bool v_is_ghost = (v == s.nv); // ghost site has the last index
 
       if (v_is_ghost) {
-        R_new = r.act(R_old);
+        R_new = r.act(R_old); // if we are, then we're moving the transformation
       } else {
-        si_old = state->spins[v];
-        si_new = r.act(si_old);
+        si_old = s.spins[v];
+        si_new = r.act(si_old); // otherwise, we're moving the spin at our site
       }
 
-      for (const v_t &vn : state->g.v_adj[v]) {
+      for (const v_t &vn : s.g.v_adj[v]) {
         X_t sj;
-        bool vn_is_ghost = (vn == state->nv);
+        bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
 
         if (!vn_is_ghost) {
-          sj = state->spins[vn];
+          sj = s.spins[vn];
         }
 
         double prob;
 
-        if (v_is_ghost || vn_is_ghost) {
+        if (v_is_ghost || vn_is_ghost) { // if this is a ghost-involved bond...
           X_t rs_old, rs_new;
           v_t non_ghost;
           if (vn_is_ghost) {
@@ -63,51 +63,54 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
             non_ghost = vn;
           }
 
-          double dE = state->H(rs_old) - state->H(rs_new);
+          double dE = s.H(rs_old) - s.H(rs_new);
+
 #ifdef FINITE_STATES
           prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
 #else
-          prob = 1.0 - exp(-dE / state->T);
+          prob = 1.0 - exp(-dE / s.T);
 #endif
 
-          state->M -= rs_old;
-          state->M += rs_new;
+          s.M -= rs_old;
+          s.M += rs_new;
 
-          state->E += dE;
+          s.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;
+          for (D_t i = 0; i < s.D; i++) {
+            L_t x = (non_ghost / (v_t)pow(s.L, s.D - i - 1)) % s.L;
 
-            state->ReF[i] -= rs_old * state->precomputed_cos[x];
-            state->ReF[i] += rs_new * state->precomputed_cos[x];
+            s.ReF[i] -= rs_old * s.precomputed_cos[x];
+            s.ReF[i] += rs_new * s.precomputed_cos[x];
 
-            state->ImF[i] -= rs_old * state->precomputed_sin[x];
-            state->ImF[i] += rs_new * state->precomputed_sin[x];
+            s.ImF[i] -= rs_old * s.precomputed_sin[x];
+            s.ImF[i] += rs_new * s.precomputed_sin[x];
           }
-        } else {
-          double dE = state->J(si_old, sj) - state->J(si_new, sj);
+        } else { // otherwise, we're at a perfectly normal bond!
+          double dE = s.J(si_old, sj) - s.J(si_new, sj);
+
 #ifdef FINITE_STATES
           prob = J_probs[state_to_ind(si_old)][state_to_ind(si_new)][state_to_ind(sj)];
 #else
-          prob = 1.0 - exp(-dE / state->T);
+          prob = 1.0 - exp(-dE / s.T);
 #endif
-          state->E += dE;
+
+          s.E += dE;
         }
 
-        if (gsl_rng_uniform(rand) < prob) { // and with probability...
+        if (gsl_rng_uniform(rand) < prob) {
           stack.push(vn); // push the neighboring vertex to the stack
         }
       }
 
       if (v_is_ghost) {
-        state->R = R_new;
+        s.R = R_new;
       } else {
-        state->spins[v] = si_new;
+        s.spins[v] = si_new;
         nv++;
       }
     }
   }
 
-  state->last_cluster_size = nv;
+  s.last_cluster_size = nv;
 }
 
diff --git a/lib/correlation.h b/lib/correlation.h
index 49c6ff2..23b1bd9 100644
--- a/lib/correlation.h
+++ b/lib/correlation.h
@@ -7,13 +7,13 @@
 #include <fftw3.h>
 
 template <class R_t, class X_t>
-double correlation_length(const state_t <R_t, X_t> *s) {
+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->ReF[j]);
+  for (D_t j = 0; j < s.D; j++) {
+    total += norm_squared(s.ReF[j]) + norm_squared(s.ReF[j]);
   }
 
-  return total / s->D;
+  return total / s.D;
 }
 
diff --git a/lib/measure.h b/lib/measure.h
index 8082dd2..2c5ffb7 100644
--- a/lib/measure.h
+++ b/lib/measure.h
@@ -29,17 +29,17 @@ FILE **measure_setup_files(unsigned char flags, unsigned long timestamp) {
 }
 
 template <class R_t, class X_t>
-std::function <void(const state_t <R_t, X_t> *)> measure_function_write_files(unsigned char flags, FILE **files, std::function <void(const state_t <R_t, X_t> *)> other_f) {
-  return [=] (const state_t <R_t, X_t> *s) {
+std::function <void(const state_t <R_t, X_t>&)> measure_function_write_files(unsigned char flags, FILE **files, std::function <void(const state_t <R_t, X_t>&)> other_f) {
+  return [=] (const state_t <R_t, X_t>& s) {
     if (flags & measurement_energy) {
-      float smaller_E = (float)s->E;
+      float smaller_E = (float)s.E;
       fwrite(&smaller_E, sizeof(float), 1, files[0]);
     }
     if (flags & measurement_clusterSize) {
-      fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, files[1]);
+      fwrite(&(s.last_cluster_size), sizeof(uint32_t), 1, files[1]);
     }
     if (flags & measurement_magnetization) {
-      write_magnetization(s->M, files[2]);
+      write_magnetization(s.M, files[2]);
     }
     if (flags & measurement_fourierZero) {
       float smaller_X = (float)correlation_length(s);
diff --git a/lib/wolff.h b/lib/wolff.h
index a4a663c..498f7f3 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -3,18 +3,18 @@
 #include "state.h"
 
 template <class R_t, class X_t>
-void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s0)> gen_R, std::function <void(const state_t <R_t, X_t> *)> measurements, gsl_rng *r, bool silent) {
+void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(gsl_rng *, X_t s0)> gen_R, std::function <void(const state_t <R_t, X_t>&)> measurements, gsl_rng *r, bool silent) {
 
 #ifdef FINITE_STATES
-  initialize_probs(s->J, s->H, s->T);
+  initialize_probs(s.J, s.H, s.T);
 #endif
 
   if (!silent) printf("\n");
   for (count_t steps = 0; steps < N; steps++) {
-    if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s->E, s->last_cluster_size);
+    if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, s.last_cluster_size);
 
-    v_t v0 = gsl_rng_uniform_int(r, s->nv);
-    R_t step = gen_R(r, s->spins[v0]);
+    v_t v0 = gsl_rng_uniform_int(r, s.nv);
+    R_t step = gen_R(r, s.spins[v0]);
     flip_cluster <R_t, X_t> (s, v0, step, r);
 
     measurements(s);
@@ -23,7 +23,7 @@ void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s
   if (!silent) {
     printf("\033[F\033[J");
   }
-  printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s->E, s->last_cluster_size);
+  printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size);
 
 }
 
diff --git a/lib/z2.h b/lib/z2.h
index fd4722e..a18d740 100644
--- a/lib/z2.h
+++ b/lib/z2.h
@@ -25,7 +25,7 @@ class z2_t {
 
   z2_t(bool x) : x(x) {}
 
-  ising_t act(const ising_t& s) {
+  ising_t act(const ising_t& s) const {
     if (x) {
       return ising_t(!s.x);
     } else {
@@ -33,7 +33,7 @@ class z2_t {
     }
   }
 
-  z2_t act(const z2_t& r) {
+  z2_t act(const z2_t& r) const {
     if (x) {
       return z2_t(!r.x);
     } else {
@@ -41,11 +41,11 @@ class z2_t {
     }
   }
 
-  ising_t act_inverse(const ising_t& s) {
+  ising_t act_inverse(const ising_t& s) const {
     return this->act(s);
   }
 
-  z2_t act_inverse(const z2_t& r) {
+  z2_t act_inverse(const z2_t& r) const {
     return this->act(r);
   }
 };