simplified bond probabilities

3 files changed, 28 insertions, 28 deletions

diff --git a/lib/wolff.h b/lib/wolff.h
index d3bc412..4270f33 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -36,7 +36,5 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
 
 graph_t *graph_add_ext(const graph_t *g);
 
-double *get_bond_probs(double T, double H, ising_state_t *s);
-
 uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
                     double *ps);
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index f45e1b5..74e5bd6 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -73,8 +73,10 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
       x[v] = !x[v];   // flip the vertex
 
       for (uint32_t i = 0; i < nn; i++) {
-        uint32_t e, v1, v2, vn;
+        bool is_ext;
+        uint32_t e, v1, v2, vn, ind;
         int32_t *bond_counter;
+        double prob;
 
         e = g->ve[g->vei[v] + i]; // select the ith bond connected to site
         v1 = g->ev[2 * e];
@@ -82,14 +84,16 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
 
         vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself
 
-        bond_counter =
-            (v1 == g->nv - 1 || v2 == g->nv - 1) ? &(c->dHb) : &(c->dJb);
+        is_ext = (v1 == g->nv - 1 || v2 == g->nv - 1);
+
+        bond_counter = is_ext ? &(c->dHb) : &(c->dJb);
+        prob = is_ext ? ps[1] : ps[0];
 
         if (x[vn] ==
             x0) { // if the neighboring site matches the flipping cluster...
           (*bond_counter)++;
 
-          if (gsl_rng_uniform(r) < ps[e]) { // and with probability ps[e]...
+          if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
             stack_push(&stack, vn); // push the neighboring vertex to the stack
           }
         } else {
@@ -106,35 +110,23 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
   return c;
 }
 
-double *get_bond_probs(double T, double H, ising_state_t *s) {
-  double p = 1 - exp(-2 / T);
-  double q = 1 - exp(-2 * fabs(H) / T);
-
-  double *ps = (double *)malloc(s->g->ne * sizeof(double));
-
-  for (uint32_t i = 0; i < s->g->ne; i++) {
-    uint32_t v1, v2;
-    v1 = s->g->ev[2 * i];
-    v2 = s->g->ev[2 * i + 1];
-    if (v1 == s->g->nv - 1 || v2 == s->g->nv - 1) {
-      ps[i] = q;
-    } else {
-      ps[i] = p;
-    }
-  }
-
-  return ps;
-}
-
 uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
                     double *ps) {
+  bool no_r, no_ps;
+  no_r = false;
+  no_ps = false;
+
   if (r == NULL) {
+    no_r = true;
     r = gsl_rng_alloc(gsl_rng_mt19937);
     gsl_rng_set(r, jst_rand_seed());
   }
 
   if (ps == NULL) {
-    ps = get_bond_probs(T, H, s);
+    no_ps = true;
+    ps = (double *)malloc(2 * sizeof(double));
+    ps[0] = 1 - exp(-2 / T);
+    ps[1] = 1 - exp(-2 * fabs(H) / T);
   }
 
   cluster_t *c = flip_cluster(s->g, ps, s->spins, r);
@@ -146,5 +138,13 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
 
   free(c);
 
+  if (no_ps) {
+    free(ps);
+  }
+
+  if (no_r) {
+    gsl_rng_free(r);
+  }
+
   return n_flips;
 }
diff --git a/src/wolff.c b/src/wolff.c
index 7675889..81906fa 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -87,7 +87,9 @@ int main(int argc, char *argv[]) {
   s->M = h->nv;
   s->H = -(1.0 * h->ne) - H * h->nv;
 
-  double *bond_probs = get_bond_probs(T, H, s);
+  double *bond_probs = (double *)malloc(2 * sizeof(double));
+  bond_probs[0] = 1 - exp(-2 / T);
+  bond_probs[1] = 1 - exp(-2 * fabs(H) / T);
 
   double diff = 1e31;
   uint64_t n_runs = 0;