9 files changed, 520 insertions, 24 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index eb86ba6..3b0de87 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,11 +2,14 @@
 cmake_minimum_required(VERSION 3.0)
 project(wolff)
 
+set (CMAKE_CXX_STANDARD 11)
+
 include_directories(lib ~/.local/include)
 link_directories(~/.local/lib)
 
 file(GLOB SOURCES lib/*.c)
-add_executable(wolff src/wolff.c ${SOURCES})
+add_executable(wolff_potts src/wolff_potts.c ${SOURCES})
+add_executable(wolff_vector src/wolff_vector.c ${SOURCES})
 
 find_package(OpenMP)
 if (OPENMP_FOUND)
@@ -14,8 +17,8 @@ if (OPENMP_FOUND)
     set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
 endif()
 
-target_link_libraries(wolff gsl m cblas fftw3 tcmalloc profiler)
-#target_link_libraries(pt_wolff gsl m jst cblas tcmalloc profiler)
+target_link_libraries(wolff_potts gsl m cblas fftw3)
+target_link_libraries(wolff_vector gsl m cblas fftw3)
 
-install(TARGETS wolff DESTINATION bin)
+install(TARGETS wolff_potts wolff_vector DESTINATION bin)
 
diff --git a/lib/dihedral.c b/lib/dihedral.c
new file mode 100644
index 0000000..623041a
--- /dev/null
+++ b/lib/dihedral.c
@@ -0,0 +1,28 @@
+
+#include "dihedral.h"
+
+dihedral_t *dihedral_compose(q_t q, q_t g1i, const dihedral_t *g2) {
+  // we only need to consider the action of reflections
+  dihedral_t *g3 = (dihedral_t *)malloc(1 * sizeof(dihedral_t));
+
+  g3->r = !g2->r;
+  g3->i = (g1i + q - g2->i) % q;
+
+  return g3;
+}
+
+q_t dihedral_act(q_t q, q_t gi, q_t s) {
+  // we only need to consider the action of reflections
+
+  return (gi + q - s) % q;
+}
+
+q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) {
+  if (g->r) {
+    return (q - ((q + s - g->i) % q)) % q;
+  } else {
+    return (q + s - g->i) % q;
+  }
+}
+
+
diff --git a/lib/dihedral.h b/lib/dihedral.h
new file mode 100644
index 0000000..813e796
--- /dev/null
+++ b/lib/dihedral.h
@@ -0,0 +1,17 @@
+
+#include <stdbool.h>
+#include <stdlib.h>
+
+#include "types.h"
+
+typedef struct {
+  q_t i;
+  bool r;
+} dihedral_t;
+
+dihedral_t *dihedral_compose(q_t q, q_t gti, const dihedral_t *g2);
+
+q_t dihedral_act(q_t q, q_t gi, q_t s);
+
+q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s);
+
diff --git a/lib/orthogonal.c b/lib/orthogonal.c
new file mode 100644
index 0000000..7e0a71b
--- /dev/null
+++ b/lib/orthogonal.c
@@ -0,0 +1,95 @@
+
+#include "orthogonal.h"
+
+void vector_replace(q_t n, double *v1, const double *v2) {
+  for (q_t i = 0; i < n; i++) {
+    v1[i] = v2[i];
+  }
+}
+
+void vector_add(q_t n, double *v1, const double *v2) {
+  for (q_t i = 0; i < n; i++) {
+    v1[i] += v2[i];
+  }
+}
+
+void vector_subtract(q_t n, double *v1, const double *v2) {
+  for (q_t i = 0; i < n; i++) {
+    v1[i] -= v2[i];
+  }
+}
+
+double *vector_rotate(q_t n, double *rot, double *vec) {
+  double *rot_vec = (double *)malloc(n * sizeof(double));
+
+  double prod = 0.0;
+  for (q_t i = 0; i < n; i++) {
+    prod += rot[i] * vec[i];
+  }
+
+  for (q_t i = 0; i < n; i++) {
+    rot_vec[i] = vec[i] - 2 * prod * rot[i];
+  }
+
+  return rot_vec;
+}
+
+double *vector_rotate_inverse(q_t n, const double *rot, const double *vec) {
+  double *rot_vec = (double *)calloc(n, sizeof(double));
+
+  for (q_t i = 0; i < n; i++) {
+    for (q_t j = 0; j < n; j++) {
+      rot_vec[i] += rot[n * j + i] * vec[j];
+    }
+  }
+
+  return rot_vec;
+}
+
+double vector_dot(q_t n, double *v1, double *v2) {
+  double dot = 0;
+
+  for (q_t i = 0; i < n; i++) {
+    dot += v1[i] * v2[i];
+  }
+
+  return dot;
+}
+
+double *orthogonal_rotate(q_t n, double *r, double *m) {
+  double *mul = (double *)calloc(n * n, sizeof(double));
+
+  for (q_t i = 0; i < n; i++) {
+    double akOki = 0;
+
+    for (q_t k = 0; k < n; k++) {
+      akOki += r[k] * m[n * k + i];
+    }
+
+    for (q_t j = 0; j < n; j++) {
+      mul[n * j + i] = m[n * j + i] - 2 * akOki * r[j];
+    }
+  }
+
+  return mul;
+}
+
+double *gen_rot(gsl_rng *r, q_t n) {
+  double *v = (double *)malloc(n * sizeof(double));
+
+  double v2 = 0;
+
+  for (q_t i = 0; i < n; i++) {
+    v[i] = gsl_ran_ugaussian(r);
+    v2 += v[i] * v[i];
+  }
+
+  double magv = sqrt(v2);
+
+  for (q_t i = 0; i < n; i++) {
+    v[i] /= magv;
+  }
+
+  return v;
+}
+
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
new file mode 100644
index 0000000..a763b08
--- /dev/null
+++ b/lib/orthogonal.h
@@ -0,0 +1,24 @@
+
+#include <stdlib.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <math.h>
+
+#include "types.h"
+
+void vector_replace(q_t n, double *v1, const double *v2);
+
+void vector_add(q_t n, double *v1, const double *v2);
+
+void vector_subtract(q_t n, double *v1, const double *v2);
+
+double *vector_rotate(q_t n, double *rot, double *vec);
+
+double *vector_rotate_inverse(q_t n, const double *rot, const double *vec);
+
+double vector_dot(q_t n, double *v1, double *v2);
+
+double *orthogonal_rotate(q_t n, double *m1, double *m2);
+
+double *gen_rot(gsl_rng *r, q_t n);
+
diff --git a/lib/wolff.h b/lib/wolff.h
index ce7040a..29286ba 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -20,6 +20,8 @@
 #include "graph.h"
 #include "tree.h"
 #include "measurement.h"
+#include "orthogonal.h"
+#include "dihedral.h"
 
 typedef struct {
   graph_t *g;
@@ -29,12 +31,27 @@ typedef struct {
   double *H;
   double *J_probs;
   double *H_probs;
+  dihedral_t *R;
   double E;
   v_t *M;
   q_t q;
 } ising_state_t;
 
+typedef struct {
+  graph_t *g;
+  double *spins;
+  double T;
+  double (*J)(double);
+  double (*H)(double *);
+  double *R;
+  double E;
+  double *M;
+  q_t n;
+} vector_state_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);
+
 graph_t *graph_add_ext(const graph_t *g);
 
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index 698d58c..e63623c 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -2,7 +2,6 @@
 #include "wolff.h"
 
 v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
-  q_t s0 = s->spins[v0];
   v_t nv = 0;
 
   ll_t *stack = NULL;     // create a new stack
@@ -16,40 +15,47 @@ v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
 
 //    if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
     if (!marks[v]) {
-      q_t s_old = s->spins[v];
-      q_t s_new = (s->spins[v] + step) % s->q;
+      q_t s_old, s_new;
+      dihedral_t *R_tmp; 
 
-      s->spins[v] = s_new;   // flip the vertex
       //tree_insert(&T, v);
       marks[v] = true;
 
+      if (v == s->g->nv - 1) {
+        R_tmp = dihedral_compose(s->q, step, s->R);
+      } else {
+        s_old = s->spins[v];
+        s_new = dihedral_act(s->q, step, s_old); // free me! I'm a new vector
+      }
+
       v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
 
       for (v_t i = 0; i < nn; i++) {
         v_t vn = s->g->v_adj[s->g->v_i[v] + i];
-        q_t sn = s->spins[vn];
+        q_t sn;
+        if (vn != s->g->nv - 1) {
+          sn = s->spins[vn];
+        }
         double prob;
 
         bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
 
-        q_t M_ind_0;
-        q_t M_ind_1;
-
         if (is_ext) {
+          q_t rs_old, rs_new;
           if (vn == s->g->nv - 1) {
-            M_ind_0 = (s_old + s->q - sn) % s->q;
-            M_ind_1 = (s_new + s->q - sn) % s->q;
+            rs_old = dihedral_inverse_act(s->q, s->R, s_old);
+            rs_new = dihedral_inverse_act(s->q, s->R, s_new);
           } else {
-            M_ind_0 = (sn + s->q - s_old) % s->q;
-            M_ind_1 = (sn + s->q - s_new) % s->q;
+            rs_old = dihedral_inverse_act(s->q, s->R, sn);
+            rs_new = dihedral_inverse_act(s->q, R_tmp, sn);
           }
-          prob = s->H_probs[M_ind_1 * s->q + M_ind_0];
-          s->M[M_ind_0]--;
-          s->M[M_ind_1]++;
-          s->E += - s->H[M_ind_1] + s->H[M_ind_0];
+          prob = s->H_probs[rs_new * s->q + rs_old];
+          s->M[rs_old]--;
+          s->M[rs_new]++;
+          s->E += - s->H[rs_new] + s->H[rs_old];
         } else {
-          M_ind_0 = (s_old + s->q - sn) % s->q;
-          M_ind_1 = (s_new + s->q - sn) % s->q;
+          q_t M_ind_0 = (s_old + s->q - sn) % s->q;
+          q_t M_ind_1 = (s_new + s->q - sn) % s->q;
           prob = s->J_probs[M_ind_1 * s->q + M_ind_0];
           s->E += - s->J[M_ind_1] + s->J[M_ind_0];
         }
@@ -59,6 +65,101 @@ v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
         }
       }
 
+      if (v == s->g->nv - 1) {
+        free(s->R);
+        s->R = R_tmp;
+      } else {
+        s->spins[v] = s_new;
+      }
+
+      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
+        nv++;
+      }
+    }
+  }
+
+  //tree_freeNode(T);
+  free(marks);
+
+  return nv;
+}
+
+v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) {
+  v_t nv = 0;
+
+  ll_t *stack = NULL;     // create a new stack
+  stack_push(&stack, v0); // push the initial vertex to the stack
+
+  //node_t *T = NULL;
+  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
+
+  while (stack != NULL) {
+    v_t v = stack_pop(&stack);
+
+//    if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
+    if (!marks[v]) {
+      double *s_old, *s_new, *R_tmp; 
+
+      //tree_insert(&T, v);
+      marks[v] = true;
+
+      if (v == s->g->nv - 1) {
+        R_tmp = orthogonal_rotate(s->n, rot, s->R);
+      } else {
+        s_old = &(s->spins[s->n * v]); // don't free me! I'm a pointer within array s->spins
+        s_new = vector_rotate(s->n, rot, s_old); // free me! I'm a new vector
+      }
+
+      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
+
+      for (v_t i = 0; i < nn; i++) {
+        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
+        double *sn;
+        if (vn != s->g->nv - 1) {
+          sn = &(s->spins[s->n * vn]);
+        }
+        double prob;
+
+        bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
+
+        if (is_ext) {
+          double *rs_old, *rs_new;
+          if (vn == s->g->nv - 1) {
+            rs_old = vector_rotate_inverse(s->n, s->R, s_old);
+            rs_new = vector_rotate_inverse(s->n, s->R, s_new);
+          } else {
+            rs_old = vector_rotate_inverse(s->n, s->R, sn);
+            rs_new = vector_rotate_inverse(s->n, R_tmp, sn);
+          }
+          double dE = s->H(rs_old) - s->H(rs_new);
+          prob = 1.0 - exp(-dE / s->T);
+          vector_subtract(s->n, s->M, rs_old);
+          vector_add(s->n, s->M, rs_new);
+          s->E += dE;
+
+          free(rs_old);
+          free(rs_new);
+        } else {
+          double dE = (s->J)(vector_dot(s->n, sn, s_old)) - (s->J)(vector_dot(s->n, sn, s_new));
+          prob = 1.0 - exp(-dE / s->T);
+          //printf("(%g %g) (%g %g) (%g %g) %g\n", s_old[0], s_old[1], s_new[0], s_new[1], sn[0], sn[1], dE);
+          //getchar();
+          s->E += dE;
+        }
+
+        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
+          stack_push(&stack, vn); // push the neighboring vertex to the stack
+        }
+      }
+
+      if (v == s->g->nv - 1) {
+        free(s->R);
+        s->R = R_tmp;
+      } else {
+        vector_replace(s->n, s_old, s_new);
+        free(s_new);
+      }
+
       if (v != s->g->nv - 1) { // count the number of non-external sites that flip
         nv++;
       }
diff --git a/src/wolff.c b/src/wolff_potts.c
index 7b6b56f..7d9beb3 100644
--- a/src/wolff.c
+++ b/src/wolff_potts.c
@@ -92,11 +92,12 @@ int main(int argc, char *argv[]) {
 
   s->q = q;
 
-  s->spins = (q_t *)calloc(h->nv + 1, sizeof(q_t));
+  s->spins = (q_t *)calloc(h->nv, sizeof(q_t));
 
   s->T = T;
   s->H = H;
   s->J = J;
+  s->R = (dihedral_t *)calloc(1, sizeof(dihedral_t));
 
   s->J_probs = (double *)calloc(pow(q, 2), sizeof(double));
   for (q_t i = 0; i < q; i++) {
@@ -159,7 +160,13 @@ int main(int argc, char *argv[]) {
 
     while (n_flips / h->nv < n) {
       v_t v0 = gsl_rng_uniform_int(r, h->nv);
-      q_t step = 1 + gsl_rng_uniform_int(r, q - 1);
+      q_t step;
+     
+      if (q == 2) {
+        step = 1;
+      } else {
+        step = gsl_rng_uniform_int(r, q);
+      }
 
       v_t tmp_flips = flip_cluster(s, v0, step, r);
       n_flips += tmp_flips;
@@ -322,12 +329,14 @@ int main(int argc, char *argv[]) {
     free(sE[i]);
     free(lifetimes[i]);
   }
+  free(lifetimes);
   free(freqs);
   free(sE);
   free(s->H_probs);
   free(s->J_probs);
   free(s->M);
   free(s->spins);
+  free(s->R);
   graph_free(s->g);
   free(s);
   free(H);
diff --git a/src/wolff_vector.c b/src/wolff_vector.c
new file mode 100644
index 0000000..c6ab695
--- /dev/null
+++ b/src/wolff_vector.c
@@ -0,0 +1,202 @@
+
+#include <wolff.h>
+
+double identity(double x) {
+  return x;
+}
+
+double zero(double *x) {
+  return 0.0;
+}
+
+double test(double *x) {
+  return -x[1];
+}
+
+int main(int argc, char *argv[]) {
+
+  L_t L = 128;
+  count_t N = (count_t)1e7;
+  count_t min_runs = 10;
+  count_t n = 3;
+  q_t q = 2;
+  D_t D = 2;
+  double T = 2.26918531421;
+  double *H = (double *)calloc(MAX_Q, sizeof(double));
+  double eps = 0;
+  bool silent = false;
+
+  int opt;
+  q_t H_ind = 0;
+
+  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:H:m:e:s")) != -1) {
+    switch (opt) {
+    case 'N':
+      N = (count_t)atof(optarg);
+      break;
+    case 'n':
+      n = (count_t)atof(optarg);
+      break;
+    case 'D':
+      D = atoi(optarg);
+      break;
+    case 'L':
+      L = atoi(optarg);
+      break;
+    case 'q':
+      q = atoi(optarg);
+      break;
+    case 'T':
+      T = atof(optarg);
+      break;
+    case 'H':
+      H[H_ind] = atof(optarg);
+      H_ind++;
+      break;
+    case 'm':
+      min_runs = atoi(optarg);
+      break;
+    case 'e':
+      eps = atof(optarg);
+      break;
+    case 's':
+      silent = true;
+      break;
+    default:
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+  gsl_rng_set(r, rand_seed());
+
+  vector_state_t *s = (vector_state_t *)calloc(1, sizeof(vector_state_t));
+
+  graph_t *h = graph_create_square(D, L);
+  s->g = graph_add_ext(h);
+
+  s->n = q;
+
+  s->spins = (double *)calloc(n * h->nv, sizeof(double));
+
+  for (v_t i = 0; i < h->nv; i++) {
+    s->spins[q * i] = 1.0;
+  }
+
+  s->T = T;
+  s->H = test;
+  s->J = identity;
+
+  s->R = (double *)calloc(q * q, sizeof(double));
+
+  for (q_t i = 0; i < q; i++) {
+    s->R[q * i + i] = 1.0;
+  }
+
+  s->M = (double *)calloc(q, sizeof(double));
+  s->M[0] = 1.0 * (double)h->nv;
+  s->E = - ((double)h->ne) * s->J(1.0) - s->H(s->M);
+
+  printf("%g %g\n", s->E, s->M[0]);
+
+  double diff = 1e31;
+  count_t n_runs = 0;
+
+  meas_t *E, *clust, **M;
+
+  M = (meas_t **)malloc(q * sizeof(meas_t *));
+  for (q_t i = 0; i < q; i++) {
+    M[i] = (meas_t *)calloc(1, sizeof(meas_t));
+  }
+
+  E = calloc(1, sizeof(meas_t));
+  clust = calloc(1, sizeof(meas_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
+           ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n",
+           n_runs, fabs(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x);
+
+    count_t n_flips = 0;
+
+    while (n_flips / h->nv < n) {
+      v_t v0 = gsl_rng_uniform_int(r, h->nv);
+      double *step = gen_rot(r, q);
+
+      v_t tmp_flips = flip_cluster_vector(s, v0, step, r);
+      free(step);
+      n_flips += tmp_flips;
+
+      update_meas(clust, tmp_flips);
+    }
+
+    for (q_t i = 0; i < q; i++) {
+      update_meas(M[i], s->M[i]);
+    }
+    update_meas(E, s->E);
+
+    diff = fabs(E->dc / E->c);
+
+    n_runs++;
+  }
+  if (!silent) {
+    printf("\033[F\033[J");
+  }
+  printf("WOLFF: sweep %" PRIu64
+         ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n",
+         n_runs, fabs(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x);
+
+  FILE *outfile = fopen("out.m", "a");
+
+  fprintf(outfile, "<|D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f", D, L, q, T);
+  fprintf(outfile, "},E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->{", E->x / h->nv, E->dx / h->nv, E->c / h->nv, E->dc / h->nv);
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", M[i]->x / h->nv);
+    if (i != q-1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "},\\[Delta]M->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", M[i]->dx / h->nv);
+    if (i != q-1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "},\\[Chi]->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", M[i]->c / h->nv);
+    if (i != q-1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "},\\[Delta]\\[Chi]->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", M[i]->dc / h->nv);
+    if (i != q-1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "}|>\n");
+
+  fclose(outfile);
+
+  free(E);
+  free(clust);
+  for (q_t i = 0; i < q; i++) {
+    free(M[i]);
+  }
+  free(M);
+  free(H);
+  free(s->M);
+  free(s->R);
+  free(s->spins);
+  graph_free(s->g);
+  free(s);
+  graph_free(h);
+  gsl_rng_free(r);
+
+  return 0;
+}
+