made the wolff hybrid more efficient

3 files changed, 11 insertions, 10 deletions

diff --git a/lib/wolff.h b/lib/wolff.h
index bd49bf0..84b3d22 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -36,6 +36,7 @@ typedef struct {
   int32_t dJb;
   int32_t dHb;
   bool hit_ghost;
+  ll_t *spins;
 } cluster_t;
 
 typedef struct {
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index 5e88a57..6b1e006 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -93,6 +93,9 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_o
 
     if (x[v] == x0 && !(c->hit_ghost && stop_on_ghost)) { // if the vertex hasn't already been flipped
       x[v] = !x[v];   // flip the vertex
+      if (stop_on_ghost) {
+        stack_push(&(c->spins), v);
+      }
 
       for (uint32_t i = 0; i < nn; i++) {
         bool is_ext;
@@ -195,19 +198,17 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
                 }
                 break;
     case WOLFF_GHOST: {
-      bool *spins_bak;
-      spins_bak = (bool *)malloc(s->g->nv * sizeof(bool));
-      memcpy(spins_bak, s->spins, s->g->nv * sizeof(bool));
-
       cluster_t *c = flip_cluster(s->g, ps, s->spins, true, r);
 
       if (c->hit_ghost) {
-        memcpy(s->spins, spins_bak, s->g->nv * sizeof(bool));
+        while (c->spins != NULL) {
+          uint32_t v = stack_pop(&(c->spins));
+          s->spins[v] = !s->spins[v];
+        }
       } else {
         s->M += - sign(H) * 2 * c->dHb;
         s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
       }
-      free(spins_bak);
 
       n_flips = c->nv;
 
diff --git a/src/wolff.c b/src/wolff.c
index 9082c81..77971e2 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -210,15 +210,14 @@ int main(int argc, char *argv[]) {
   FILE *outfile = fopen("out.dat", "a");
 
   double tau = batch_size * corrmE->c / corrE->c;
-  double dtau = tau * sqrt(pow(corrmE->dc / corrmE->c, 2) + pow(corrE->dc / corrE->c, 2));
 
   fprintf(outfile,
-          "%u %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L,
-          T, H, n_runs,
+          "%u %.15f %.15f %" PRIu64 " %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n",
+          L, T, H, n_runs, n_steps,
           E->x / h->nv, E->dx / h->nv, M->x / h->nv, M->dx / h->nv, E->c / h->nv, E->dc / h->nv, M->c / h->nv, M->dc / h->nv,
           eE->n, eE->x / h->nv, eE->dx / h->nv, eM->x / h->nv, eM->dx / h->nv, eE->c / h->nv, eE->dc / h->nv, eM->c / h->nv, eM->dc / h->nv,
           mE->n, mE->x / h->nv, mE->dx / h->nv, mM->x / h->nv, mM->dx / h->nv, mE->c / h->nv, mE->dc / h->nv, mM->c / h->nv, mM->dc / h->nv,
-          aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv, tau, dtau);
+          aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv, tau); 
   fclose(outfile);
 
   if (output_state) {