made some simplifications to the computation of magnetization and energy

3 files changed, 11 insertions, 22 deletions

diff --git a/lib/wolff.h b/lib/wolff.h
index 6beb03e..ca7635a 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -29,6 +29,7 @@ typedef struct ll_tag {
 typedef struct {
   int32_t nv;
   double dH;
+  int32_t dM;
 } cluster_t;
 
 double get_hamiltonian(graph_t *g, double *coupling, bool *x);
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index adef40f..d361344 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -123,10 +123,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g
   }
 
   c->dH = n_bonds + H * n_h_bonds;
-
-  if (x0) {
-    c->nv = -c->nv;
-  }
+  c->dM = n_h_bonds;
 
   return c;
 }
@@ -137,7 +134,7 @@ double hh(double th) {
 
 double *get_bond_probs(double T, double H, ising_state_t *s) {
   double p = 1 - exp(-2 / T);
-  double q = 1 - exp(-2 * H / T);
+  double q = 1 - exp(-2 * fabs(H) / T);
 
   double *ps = (double *)malloc(s->g->ne * sizeof(double));
 
@@ -167,13 +164,11 @@ int32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps)
 
   cluster_t *c = flip_cluster(s->g, ps, H, s->spins, r);
 
-  s->M += 2 * c->nv;
+  s->M += -2 * c->dM;
   s->H += 2 * c->dH;
 
-  int32_t n_flipped = c->nv;
-
   free(c);
 
-  return n_flipped;
+  return c->nv;
 }
 
diff --git a/src/wolff.c b/src/wolff.c
index 22e6857..5f06665 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -56,7 +56,7 @@ int main(int argc, char *argv[]) {
 
   s->g = graph_add_ext(h);
   s->spins = (bool *)calloc(h->nv + 1, sizeof(bool));
-  s->M = -h->nv;
+  s->M = h->nv;
   s->H = -(1.0 * h->ne) - H * h->nv;
 
   double *bond_probs = get_bond_probs(T, H, s);
@@ -82,24 +82,17 @@ int main(int argc, char *argv[]) {
       n_clust++;
     }
 
-    int32_t MM;
-    if (s->spins[h->nv]) {
-      MM = s->M;
-    } else {
-      MM = -s->M;
-    }
-
     E1 = E1 * (n_runs / (n_runs + 1.)) + s->H * 1. / (n_runs + 1.);
-    M1 = M1 * (n_runs / (n_runs + 1.)) + MM * 1. / (n_runs + 1.);
+    M1 = M1 * (n_runs / (n_runs + 1.)) + s->M * 1. / (n_runs + 1.);
     E2 = E2 * (n_runs / (n_runs + 1.)) + pow(s->H, 2) * 1. / (n_runs + 1.);
-    M2 = M2 * (n_runs / (n_runs + 1.)) + pow(MM, 2) * 1. / (n_runs + 1.);
+    M2 = M2 * (n_runs / (n_runs + 1.)) + pow(s->M, 2) * 1. / (n_runs + 1.);
 
-    Mmu2 = Mmu2 * (n_runs / (n_runs + 1.)) + pow(MM - M1, 2) * 1. / (n_runs + 1.);
-    Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) + pow(MM - M1, 4) * 1. / (n_runs + 1.);
+    Mmu2 = Mmu2 * (n_runs / (n_runs + 1.)) + pow(s->M - M1, 2) * 1. / (n_runs + 1.);
+    Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) + pow(s->M - M1, 4) * 1. / (n_runs + 1.);
     Emu2 = Emu2 * (n_runs / (n_runs + 1.)) + pow(s->H - E1, 2) * 1. / (n_runs + 1.);
     Emu4 = Emu4 * (n_runs / (n_runs + 1.)) + pow(s->H - E1, 4) * 1. / (n_runs + 1.);
 
-    if (n_runs > 1){
+    if (n_runs > 1) {
       double Msigma2 = n_runs / (n_runs - 1) * (M2 - pow(M1, 2));
       X = Msigma2 / T;
       dX = sqrt((Mmu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Mmu2, 2)) / n_runs) / T;