lots of changes

7 files changed, 162 insertions, 10 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index 2de17e5..d118735 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -22,6 +22,7 @@
 #include "orthogonal.h"
 #include "dihedral.h"
 #include "dihinf.h"
+#include "yule_walker.h"
 
 typedef struct {
   graph_t *g;
@@ -62,6 +63,13 @@ typedef struct {
   q_t n;
 } vector_state_t;
 
+typedef enum {
+  VECTOR,
+  MODULATED,
+  CUBIC,
+  QUADRATIC
+} vector_field_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);
diff --git a/lib/convex.c b/lib/convex.c
index df506f9..7255ad4 100644
--- a/lib/convex.c
+++ b/lib/convex.c
@@ -93,7 +93,6 @@ double *get_convex_minorant(count_t n, double *Gammas) {
 
   while (L != NULL) {
     free(L->p);
-    list_t *L_save = L;
     L = L->next;
     free(L);
   }
diff --git a/lib/measurement.h b/lib/measurement.h
index a24df9f..eaa260b 100644
--- a/lib/measurement.h
+++ b/lib/measurement.h
@@ -1,4 +1,6 @@
 
+#pragma once
+
 #include <math.h>
 #include <stdlib.h>
 
diff --git a/lib/yule_walker.c b/lib/yule_walker.c
new file mode 100644
index 0000000..77b3e96
--- /dev/null
+++ b/lib/yule_walker.c
@@ -0,0 +1,41 @@
+
+#include "yule_walker.h"
+
+double yule_walker(const autocorr_t *a) {
+  gsl_vector *rhos = gsl_vector_alloc(a->W);
+  gsl_matrix *mat = gsl_matrix_alloc(a->W, a->W);
+  gsl_vector *phis = gsl_vector_alloc(a->W);
+
+  for (count_t i = 0; i < a->W; i++) {
+    gsl_vector_set(rhos, i, rho(a, i));
+  }
+
+  for (count_t i = 0; i < a->W; i++) {
+    gsl_matrix_set(mat, i, i, 1.0);
+
+    for (count_t j = 0; j < i; j++) {
+      gsl_matrix_set(mat, i, j, gsl_vector_get(rhos, i - 1 - j));
+    }
+
+    for (count_t j = 0; j < a->W - 1 - i; j++) {
+      gsl_matrix_set(mat, i, i + 1 + j, gsl_vector_get(rhos, j));
+    }
+  }
+
+  int out = gsl_linalg_cholesky_solve(mat, rhos, phis);
+
+  double rhophi = 0;
+  double onephi = 0;
+
+  for (count_t i = 0; i < a->W; i++) {
+    rhophi += gsl_vector_get(rhos, i) * gsl_vector_get(phis, i);
+    onephi += gsl_vector_get(phis, i);
+  }
+
+  gsl_vector_free(rhos);
+  gsl_matrix_free(mat);
+  gsl_vector_free(phis);
+
+  return (1 - rhophi) / pow(1 - onephi, 2);
+}
+
diff --git a/lib/yule_walker.h b/lib/yule_walker.h
new file mode 100644
index 0000000..883ea03
--- /dev/null
+++ b/lib/yule_walker.h
@@ -0,0 +1,12 @@
+
+#pragma once
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+
+#include "types.h"
+#include "measurement.h"
+
+double yule_walker(const autocorr_t *a);
+
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) {