lots of changes

2 files changed, 99 insertions, 9 deletions

diff --git a/src/wolff_potts.c b/src/wolff_potts.c
index e209043..d200b35 100644
--- a/src/wolff_potts.c
+++ b/src/wolff_potts.c
@@ -22,6 +22,7 @@ int main(int argc, char *argv[]) {
   bool snapshots = false;
   bool snapshot = false;
   bool record_autocorrelation = false;
+  bool record_distribution = false;
   count_t W = 10;
   count_t ac_skip = 1;
 
@@ -29,7 +30,7 @@ int main(int argc, char *argv[]) {
   q_t J_ind = 0;
   q_t H_ind = 0;
 
-  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:J:H:m:e:IpsSPak:W:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:J:H:m:e:IpsSPak:W:d")) != -1) {
     switch (opt) {
     case 'N':
       N = (count_t)atof(optarg);
@@ -87,6 +88,9 @@ int main(int argc, char *argv[]) {
     case 'W':
       W = (count_t)atof(optarg);
       break;
+    case 'd':
+      record_distribution = true;
+      break;
     default:
       exit(EXIT_FAILURE);
     }
@@ -173,6 +177,11 @@ int main(int argc, char *argv[]) {
     autocorr->OO = (double *)calloc(2 * W + 1, sizeof(double));
   }
 
+  count_t *cluster_dist;
+  if (record_distribution) {
+    cluster_dist = (count_t *)calloc(h->nv, sizeof(count_t));
+  }
+
   if (!silent) printf("\n");
   while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) {
     if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64
@@ -201,6 +210,10 @@ int main(int argc, char *argv[]) {
         if (record_autocorrelation && n_steps % ac_skip == 0) {
           update_autocorr(autocorr, s->E);
         }
+
+        if (record_distribution) {
+          cluster_dist[tmp_flips - 1]++;
+        }
       }
 
     }
@@ -269,6 +282,7 @@ int main(int argc, char *argv[]) {
   }
 
   double tau = 0;
+  double tauyw;
   bool tau_failed = false;
 
   if (record_autocorrelation) {
@@ -304,14 +318,15 @@ int main(int argc, char *argv[]) {
       ttau += conv_Gamma[i];
     }
     
+    tau = ttau * ac_skip * clust->x / h->nv;
+    tauyw = yule_walker(autocorr) * ac_skip * clust->x / h->nv;
+    
     free(Gammas);
     free(autocorr->OO);
     while (autocorr->Op != NULL) {
       stack_pop_d(&(autocorr->Op));
     }
     free(autocorr);
-    
-    tau = ttau * ac_skip * clust->x / h->nv;
   }
 
   if (tau_failed) {
@@ -320,7 +335,7 @@ int main(int argc, char *argv[]) {
 
   FILE *outfile = fopen("out.m", "a");
 
-  fprintf(outfile, "<|N->%" PRIcount ",D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,J->{", N, D, L, q, T);
+  fprintf(outfile, "<|N->%" PRIcount ",n->%" PRIcount ",D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,J->{", N, n, D, L, q, T);
   for (q_t i = 0; i < q; i++) {
     fprintf(outfile, "%.15f", J[i]);
     if (i != q-1) {
@@ -396,7 +411,19 @@ int main(int argc, char *argv[]) {
   for (q_t i = 0; i < q; i++) {
     fprintf(outfile, ",Subscript[f,%" PRIq "]->%.15f,Subscript[\\[Delta]f,%" PRIq "]->%.15f", i, (double)freqs[i] / (double)n_runs, i, sqrt(freqs[i]) / (double)n_runs);
   }
-  fprintf(outfile, ",Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau);
+  fprintf(outfile, ",Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f,\\[Tau]yw->%.15f", clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau, tauyw);
+  if (record_distribution) {
+    fprintf(outfile, ",S->{");
+    for (v_t i = 0; i < h->nv; i++) {
+      fprintf(outfile, "%" PRIcount, cluster_dist[i]);
+      if (i != h->nv - 1) {
+        fprintf(outfile, ",");
+      }
+    }
+    fprintf(outfile, "}");
+    free(cluster_dist);
+  }
+  fprintf(outfile, "|>\n");
 
   fclose(outfile);
 
diff --git a/src/wolff_vector.c b/src/wolff_vector.c
index c87e60d..da93a15 100644
--- a/src/wolff_vector.c
+++ b/src/wolff_vector.c
@@ -35,6 +35,28 @@ double modulated(q_t n, double *H_info, double *x) {
   return H_info[0] * cos(H_info[1] * theta(x[0], x[1]));
 }
 
+double cubic(q_t n, double *H_info, double *x) {
+  double tmp = 0;
+
+  for (q_t i = 0; i < n; i++) {
+    tmp += pow(x[i], 4);
+  }
+
+  return - H_info[0] * tmp;
+}
+
+double quadratic(q_t n, double *H_info, double *x) {
+  double tmp = 0;
+
+  tmp += pow(x[0], 2);
+
+  for (q_t i = 1; i < n; i++) {
+    tmp += - 1.0 / (n - 1.0) * pow(x[i], 2);
+  }
+
+  return - 0.5 *  H_info[0] * tmp;
+}
+
 int main(int argc, char *argv[]) {
 
   L_t L = 128;
@@ -48,13 +70,14 @@ int main(int argc, char *argv[]) {
   double eps = 0;
   bool silent = false;
   bool record_autocorrelation = false;
+  vector_field_t H_type = VECTOR;
   count_t ac_skip = 1;
   count_t W = 10;
 
   int opt;
   q_t H_ind = 0;
 
-  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:H:m:e:saS:W:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:H:m:e:saS:W:f:")) != -1) {
     switch (opt) {
     case 'N':
       N = (count_t)atof(optarg);
@@ -96,6 +119,22 @@ int main(int argc, char *argv[]) {
     case 'W':
       W = (count_t)atof(optarg);
       break;
+    case 'f':
+      switch (atoi(optarg)) {
+        case 0:
+          H_type = VECTOR;
+          break;
+        case 1:
+          H_type = MODULATED;
+          break;
+        case 2:
+          H_type = CUBIC;
+          break;
+        case 3:
+          H_type = QUADRATIC;
+          break;
+      }
+      break;
     default:
       exit(EXIT_FAILURE);
     }
@@ -119,7 +158,20 @@ int main(int argc, char *argv[]) {
 
   s->H_info = H;
   s->T = T;
-  s->H = dot;
+  switch (H_type) {
+    case VECTOR:
+      s->H = dot;
+      break;
+    case MODULATED:
+      s->H = modulated;
+      break;
+    case CUBIC:
+      s->H = cubic;
+      break;
+    case QUADRATIC:
+      s->H = quadratic;
+      break;
+  }
   s->J = identity;
 
   s->R = (double *)calloc(q * q, sizeof(double));
@@ -190,7 +242,7 @@ int main(int argc, char *argv[]) {
     meas_update(aM, sqrt(aM_val));
     meas_update(E, s->E);
 
-    diff = fabs(meas_dc(aM) / meas_c(aM));
+    diff = fabs(meas_dx(clust) / clust->x);
 
     n_runs++;
   }
@@ -237,6 +289,17 @@ int main(int argc, char *argv[]) {
     for (uint64_t i = 0; i < n + 1; i++) {
       ttau += conv_Gamma[i];
     }
+
+    FILE *autocorr_file = fopen("autocorr.dat", "a");
+
+    printf("%g %g\n", Gammas[0], conv_Gamma[0]);
+
+    for (count_t i = 0; i < n+1; i++) {
+      fprintf(autocorr_file, "%g ", conv_Gamma[i]);
+    }
+    fprintf(autocorr_file, "\n");
+
+    fclose(autocorr_file);
     
     free(Gammas);
     free(autocorr->OO);
@@ -254,7 +317,7 @@ int main(int argc, char *argv[]) {
 
   FILE *outfile = fopen("out.m", "a");
 
-  fprintf(outfile, "<|N->%" PRIcount ",D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,H->{", N, D, L, q, T);
+  fprintf(outfile, "<|N->%" PRIcount ",n->%" PRIcount ",D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,H->{", N, n, D, L, q, T);
   for (q_t i = 0; i < q; i++) {
     fprintf(outfile, "%.15f", H[i]);
     if (i != q-1) {