From c960d23a9e147c2043df852460427b95bf752583 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 1 Mar 2018 01:43:33 -0500
Subject: more cosmetic changes

---
 lib/cluster.c      | 182 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 lib/cluster.h      |  57 +++++++++++++++++
 lib/wolff.h        |  58 -----------------
 lib/wolff_tools.c  | 174 --------------------------------------------------
 src/wolff_potts.c  |   4 +-
 src/wolff_vector.c |   4 +-
 6 files changed, 245 insertions(+), 234 deletions(-)
 create mode 100644 lib/cluster.c
 create mode 100644 lib/cluster.h
 delete mode 100644 lib/wolff.h
 delete mode 100644 lib/wolff_tools.c

diff --git a/lib/cluster.c b/lib/cluster.c
new file mode 100644
index 0000000..970620e
--- /dev/null
+++ b/lib/cluster.c
@@ -0,0 +1,182 @@
+
+#include "cluster.h"
+
+v_t flip_cluster(ising_state_t *s, v_t v0, q_t rot, gsl_rng *r) {
+  v_t nv = 0;
+
+  ll_t *stack = NULL;     // create a new stack
+  stack_push(&stack, v0); // push the initial vertex to the stack
+
+  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
+
+  while (stack != NULL) {
+    v_t v = stack_pop(&stack);
+
+    if (!marks[v]) {
+      q_t s_old, s_new;
+      dihedral_t *R_new; 
+      bool external_flipped;
+
+      marks[v] = true;
+
+      if (v == s->g->nv - 1) {
+        R_new = dihedral_compose(s->q, rot, s->R);
+        external_flipped = true;
+      } else {
+        s_old = s->spins[v];
+        s_new = dihedral_act(s->q, rot, s_old);
+        external_flipped = false;
+      }
+
+      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
+
+      for (v_t i = 0; i < nn; i++) {
+        q_t sn;
+        double prob;
+        bool external_neighbor = false;
+
+        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
+
+        if (vn == s->g->nv - 1) {
+          external_neighbor = true;
+        } else {
+          sn = s->spins[vn];
+        }
+
+        if (external_flipped || external_neighbor) {
+          q_t rot_s_old, rot_s_new;
+
+          if (external_neighbor) {
+            rot_s_old = dihedral_inverse_act(s->q, s->R, s_old);
+            rot_s_new = dihedral_inverse_act(s->q, s->R, s_new);
+          } else {
+            rot_s_old = dihedral_inverse_act(s->q, s->R, sn);
+            rot_s_new = dihedral_inverse_act(s->q, R_new, sn);
+          }
+
+          prob = s->H_probs[rot_s_new * s->q + rot_s_old];
+
+          s->M[rot_s_old]--;
+          s->M[rot_s_new]++;
+
+          s->E += - s->H[rot_s_new] + s->H[rot_s_old];
+        } else {
+          q_t diff_old = (s_old + s->q - sn) % s->q;
+          q_t diff_new = (s_new + s->q - sn) % s->q;
+
+          prob = s->J_probs[diff_new * s->q + diff_old];
+
+          s->E += - s->J[diff_new] + s->J[diff_old];
+        }
+
+        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
+          stack_push(&stack, vn); // push the neighboring vertex to the stack
+        }
+      }
+
+      if (external_flipped) {
+        free(s->R);
+        s->R = R_new;
+      } else {
+        s->spins[v] = s_new;
+      }
+
+      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
+        nv++;
+      }
+    }
+  }
+
+  free(marks);
+
+  return nv;
+}
+
+v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) {
+  v_t nv = 0;
+
+  ll_t *stack = NULL;     // create a new stack
+  stack_push(&stack, v0); // push the initial vertex to the stack
+
+  //node_t *T = NULL;
+  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
+
+  while (stack != NULL) {
+    v_t v = stack_pop(&stack);
+
+//    if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
+    if (!marks[v]) {
+      double *s_old, *s_new, *R_tmp; 
+
+      //tree_insert(&T, v);
+      marks[v] = true;
+
+      if (v == s->g->nv - 1) {
+        R_tmp = orthogonal_rotate(s->n, rot, s->R);
+      } else {
+        s_old = &(s->spins[s->n * v]); // don't free me! I'm a pointer within array s->spins
+        s_new = vector_rotate(s->n, rot, s_old); // free me! I'm a new vector
+      }
+
+      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
+
+      for (v_t i = 0; i < nn; i++) {
+        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
+        double *sn;
+        if (vn != s->g->nv - 1) {
+          sn = &(s->spins[s->n * vn]);
+        }
+        double prob;
+
+        bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
+
+        if (is_ext) {
+          double *rs_old, *rs_new;
+          if (vn == s->g->nv - 1) {
+            rs_old = vector_rotate_inverse(s->n, s->R, s_old);
+            rs_new = vector_rotate_inverse(s->n, s->R, s_new);
+          } else {
+            rs_old = vector_rotate_inverse(s->n, s->R, sn);
+            rs_new = vector_rotate_inverse(s->n, R_tmp, sn);
+          }
+          double dE = s->H(s->n, s->H_info, rs_old) - s->H(s->n, s->H_info, rs_new);
+          prob = 1.0 - exp(-dE / s->T);
+          vector_subtract(s->n, s->M, rs_old);
+          vector_add(s->n, s->M, rs_new);
+          s->E += dE;
+
+          free(rs_old);
+          free(rs_new);
+        } else {
+          double dE = (s->J)(vector_dot(s->n, sn, s_old)) - (s->J)(vector_dot(s->n, sn, s_new));
+          prob = 1.0 - exp(-dE / s->T);
+          //printf("(%g %g) (%g %g) (%g %g) %g\n", s_old[0], s_old[1], s_new[0], s_new[1], sn[0], sn[1], dE);
+          //getchar();
+          s->E += dE;
+        }
+
+        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
+          stack_push(&stack, vn); // push the neighboring vertex to the stack
+        }
+      }
+
+      if (v == s->g->nv - 1) {
+        free(s->R);
+        s->R = R_tmp;
+      } else {
+        vector_replace(s->n, s_old, s_new);
+        free(s_new);
+      }
+
+      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
+        nv++;
+      }
+    }
+  }
+
+  //tree_freeNode(T);
+  free(marks);
+
+  return nv;
+}
+
diff --git a/lib/cluster.h b/lib/cluster.h
new file mode 100644
index 0000000..095ed61
--- /dev/null
+++ b/lib/cluster.h
@@ -0,0 +1,57 @@
+
+#pragma once
+
+#include <assert.h>
+#include <fftw3.h>
+#include <float.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_rng.h>
+#include <inttypes.h>
+#include <math.h>
+#include <stdbool.h>
+#include <string.h>
+#include <sys/types.h>
+
+#include "types.h"
+#include "rand.h"
+#include "stack.h"
+#include "convex.h"
+#include "graph.h"
+#include "tree.h"
+#include "measurement.h"
+#include "orthogonal.h"
+#include "dihedral.h"
+
+typedef struct {
+  graph_t *g;
+  q_t *spins;
+  double T;
+  double *J;
+  double *H;
+  double *J_probs;
+  double *H_probs;
+  dihedral_t *R;
+  double E;
+  v_t *M;
+  q_t q;
+} ising_state_t;
+
+typedef struct {
+  graph_t *g;
+  double *spins;
+  double T;
+  double (*J)(double);
+  double (*H)(q_t, double *, double *);
+  double *H_info;
+  double *R;
+  double E;
+  double *M;
+  q_t n;
+} vector_state_t;
+
+v_t flip_cluster(ising_state_t *s, v_t v0, q_t s1, gsl_rng *r);
+
+v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r);
+
+graph_t *graph_add_ext(const graph_t *g);
+
diff --git a/lib/wolff.h b/lib/wolff.h
deleted file mode 100644
index 18ae6dd..0000000
--- a/lib/wolff.h
+++ /dev/null
@@ -1,58 +0,0 @@
-
-#pragma once
-
-#include <assert.h>
-#include <fftw3.h>
-#include <float.h>
-#include <getopt.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_rng.h>
-#include <inttypes.h>
-#include <math.h>
-#include <stdbool.h>
-#include <string.h>
-#include <sys/types.h>
-
-#include "types.h"
-#include "rand.h"
-#include "stack.h"
-#include "convex.h"
-#include "graph.h"
-#include "tree.h"
-#include "measurement.h"
-#include "orthogonal.h"
-#include "dihedral.h"
-
-typedef struct {
-  graph_t *g;
-  q_t *spins;
-  double T;
-  double *J;
-  double *H;
-  double *J_probs;
-  double *H_probs;
-  dihedral_t *R;
-  double E;
-  v_t *M;
-  q_t q;
-} ising_state_t;
-
-typedef struct {
-  graph_t *g;
-  double *spins;
-  double T;
-  double (*J)(double);
-  double (*H)(q_t, double *, double *);
-  double *H_info;
-  double *R;
-  double E;
-  double *M;
-  q_t n;
-} vector_state_t;
-
-v_t flip_cluster(ising_state_t *s, v_t v0, q_t s1, gsl_rng *r);
-
-v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r);
-
-graph_t *graph_add_ext(const graph_t *g);
-
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
deleted file mode 100644
index 7a52546..0000000
--- a/lib/wolff_tools.c
+++ /dev/null
@@ -1,174 +0,0 @@
-
-#include "wolff.h"
-
-v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
-  v_t nv = 0;
-
-  ll_t *stack = NULL;     // create a new stack
-  stack_push(&stack, v0); // push the initial vertex to the stack
-
-  //node_t *T = NULL;
-  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
-
-  while (stack != NULL) {
-    v_t v = stack_pop(&stack);
-
-//    if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
-    if (!marks[v]) {
-      q_t s_old, s_new;
-      dihedral_t *R_tmp; 
-
-      //tree_insert(&T, v);
-      marks[v] = true;
-
-      if (v == s->g->nv - 1) {
-        R_tmp = dihedral_compose(s->q, step, s->R);
-      } else {
-        s_old = s->spins[v];
-        s_new = dihedral_act(s->q, step, s_old); // free me! I'm a new vector
-      }
-
-      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
-
-      for (v_t i = 0; i < nn; i++) {
-        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
-        q_t sn;
-        if (vn != s->g->nv - 1) {
-          sn = s->spins[vn];
-        }
-        double prob;
-
-        bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
-
-        if (is_ext) {
-          q_t rs_old, rs_new;
-          if (vn == s->g->nv - 1) {
-            rs_old = dihedral_inverse_act(s->q, s->R, s_old);
-            rs_new = dihedral_inverse_act(s->q, s->R, s_new);
-          } else {
-            rs_old = dihedral_inverse_act(s->q, s->R, sn);
-            rs_new = dihedral_inverse_act(s->q, R_tmp, sn);
-          }
-          prob = s->H_probs[rs_new * s->q + rs_old];
-          s->M[rs_old]--;
-          s->M[rs_new]++;
-          s->E += - s->H[rs_new] + s->H[rs_old];
-        } else {
-          q_t M_ind_0 = (s_old + s->q - sn) % s->q;
-          q_t M_ind_1 = (s_new + s->q - sn) % s->q;
-          prob = s->J_probs[M_ind_1 * s->q + M_ind_0];
-          s->E += - s->J[M_ind_1] + s->J[M_ind_0];
-        }
-
-        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
-          stack_push(&stack, vn); // push the neighboring vertex to the stack
-        }
-      }
-
-      if (v == s->g->nv - 1) {
-        free(s->R);
-        s->R = R_tmp;
-      } else {
-        s->spins[v] = s_new;
-      }
-
-      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
-        nv++;
-      }
-    }
-  }
-
-  //tree_freeNode(T);
-  free(marks);
-
-  return nv;
-}
-
-v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) {
-  v_t nv = 0;
-
-  ll_t *stack = NULL;     // create a new stack
-  stack_push(&stack, v0); // push the initial vertex to the stack
-
-  //node_t *T = NULL;
-  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
-
-  while (stack != NULL) {
-    v_t v = stack_pop(&stack);
-
-//    if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
-    if (!marks[v]) {
-      double *s_old, *s_new, *R_tmp; 
-
-      //tree_insert(&T, v);
-      marks[v] = true;
-
-      if (v == s->g->nv - 1) {
-        R_tmp = orthogonal_rotate(s->n, rot, s->R);
-      } else {
-        s_old = &(s->spins[s->n * v]); // don't free me! I'm a pointer within array s->spins
-        s_new = vector_rotate(s->n, rot, s_old); // free me! I'm a new vector
-      }
-
-      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
-
-      for (v_t i = 0; i < nn; i++) {
-        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
-        double *sn;
-        if (vn != s->g->nv - 1) {
-          sn = &(s->spins[s->n * vn]);
-        }
-        double prob;
-
-        bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
-
-        if (is_ext) {
-          double *rs_old, *rs_new;
-          if (vn == s->g->nv - 1) {
-            rs_old = vector_rotate_inverse(s->n, s->R, s_old);
-            rs_new = vector_rotate_inverse(s->n, s->R, s_new);
-          } else {
-            rs_old = vector_rotate_inverse(s->n, s->R, sn);
-            rs_new = vector_rotate_inverse(s->n, R_tmp, sn);
-          }
-          double dE = s->H(s->n, s->H_info, rs_old) - s->H(s->n, s->H_info, rs_new);
-          prob = 1.0 - exp(-dE / s->T);
-          vector_subtract(s->n, s->M, rs_old);
-          vector_add(s->n, s->M, rs_new);
-          s->E += dE;
-
-          free(rs_old);
-          free(rs_new);
-        } else {
-          double dE = (s->J)(vector_dot(s->n, sn, s_old)) - (s->J)(vector_dot(s->n, sn, s_new));
-          prob = 1.0 - exp(-dE / s->T);
-          //printf("(%g %g) (%g %g) (%g %g) %g\n", s_old[0], s_old[1], s_new[0], s_new[1], sn[0], sn[1], dE);
-          //getchar();
-          s->E += dE;
-        }
-
-        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
-          stack_push(&stack, vn); // push the neighboring vertex to the stack
-        }
-      }
-
-      if (v == s->g->nv - 1) {
-        free(s->R);
-        s->R = R_tmp;
-      } else {
-        vector_replace(s->n, s_old, s_new);
-        free(s_new);
-      }
-
-      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
-        nv++;
-      }
-    }
-  }
-
-  //tree_freeNode(T);
-  free(marks);
-
-  return nv;
-}
-
diff --git a/src/wolff_potts.c b/src/wolff_potts.c
index 7d9beb3..8dc4471 100644
--- a/src/wolff_potts.c
+++ b/src/wolff_potts.c
@@ -1,5 +1,7 @@
 
-#include <wolff.h>
+#include <getopt.h>
+
+#include <cluster.h>
 
 int main(int argc, char *argv[]) {
 
diff --git a/src/wolff_vector.c b/src/wolff_vector.c
index 904fcfd..bbfc4da 100644
--- a/src/wolff_vector.c
+++ b/src/wolff_vector.c
@@ -1,5 +1,7 @@
 
-#include <wolff.h>
+#include <getopt.h>
+
+#include <cluster.h>
 
 double identity(double x) {
   return x;
-- 
cgit v1.2.3-54-g00ecf