summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/cluster_finite.c100
-rw-r--r--lib/cluster_finite.h42
-rw-r--r--lib/dihedral.c12
-rw-r--r--lib/dihedral.h2
-rw-r--r--lib/measurement.h5
-rw-r--r--lib/symmetric.c27
-rw-r--r--lib/symmetric.h13
-rw-r--r--lib/types.h3
-rw-r--r--lib/wolff_finite.c70
9 files changed, 274 insertions, 0 deletions
diff --git a/lib/cluster_finite.c b/lib/cluster_finite.c
new file mode 100644
index 0000000..f11a3ea
--- /dev/null
+++ b/lib/cluster_finite.c
@@ -0,0 +1,100 @@
+
+#include "cluster_finite.h"
+
+v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) {
+ q_t *rot = s->transformations + s->q * rot_ind;
+ q_t *R_inv = symmetric_invert(s->q, s->R);
+ v_t nv = 0;
+
+ ll_t *stack = NULL; // create a new stack
+ stack_push(&stack, v0); // push the initial vertex to the stack
+
+ bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
+
+ while (stack != NULL) {
+ v_t v = stack_pop(&stack);
+
+ if (!marks[v]) {
+ q_t s_old, s_new;
+ q_t *R_new, *R_inv_new;
+ bool external_flipped;
+
+ marks[v] = true;
+
+ if (v == s->g->nv - 1) {
+ R_new = symmetric_compose(s->q, rot, s->R);
+ R_inv_new = symmetric_invert(s->q, R_new);
+ external_flipped = true;
+ } else {
+ s_old = s->spins[v];
+ s_new = symmetric_act(rot, s_old);
+ external_flipped = false;
+ }
+
+ v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
+
+ for (v_t i = 0; i < nn; i++) {
+ q_t sn;
+ double prob;
+ bool external_neighbor = false;
+
+ v_t vn = s->g->v_adj[s->g->v_i[v] + i];
+
+ if (vn == s->g->nv - 1) {
+ external_neighbor = true;
+ } else {
+ sn = s->spins[vn];
+ }
+
+ if (external_flipped || external_neighbor) {
+ q_t rot_s_old, rot_s_new;
+
+ if (external_neighbor) {
+ rot_s_old = symmetric_act(R_inv, s_old);
+ rot_s_new = symmetric_act(R_inv, s_new);
+ } else {
+ rot_s_old = symmetric_act(R_inv, sn);
+ rot_s_new = symmetric_act(R_inv_new, sn);
+ }
+
+ prob = s->H_probs[rot_s_new * s->q + rot_s_old];
+
+ s->M[rot_s_old]--;
+ s->M[rot_s_new]++;
+
+ s->E += - s->H[rot_s_new] + s->H[rot_s_old];
+ } else {
+ q_t diff_old = (s_old + s->q - sn) % s->q;
+ q_t diff_new = (s_new + s->q - sn) % s->q;
+
+ prob = s->J_probs[diff_new * s->q + diff_old];
+
+ s->E += - s->J[diff_new] + s->J[diff_old];
+ }
+
+ if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
+ stack_push(&stack, vn); // push the neighboring vertex to the stack
+ }
+ }
+
+ if (external_flipped) {
+ free(s->R);
+ free(R_inv);
+ s->R = R_new;
+ R_inv = R_inv_new;
+ } else {
+ s->spins[v] = s_new;
+ }
+
+ if (v != s->g->nv - 1) { // count the number of non-external sites that flip
+ nv++;
+ }
+ }
+ }
+
+ free(marks);
+ free(R_inv);
+
+ return nv;
+}
+
diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h
new file mode 100644
index 0000000..abdc8fc
--- /dev/null
+++ b/lib/cluster_finite.h
@@ -0,0 +1,42 @@
+
+#pragma once
+
+#include <assert.h>
+#include <fftw3.h>
+#include <float.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_rng.h>
+#include <inttypes.h>
+#include <math.h>
+#include <stdbool.h>
+#include <string.h>
+#include <sys/types.h>
+
+#include "types.h"
+#include "rand.h"
+#include "stack.h"
+#include "convex.h"
+#include "graph.h"
+#include "tree.h"
+#include "measurement.h"
+#include "symmetric.h"
+#include "yule_walker.h"
+
+typedef struct {
+ graph_t *g;
+ q_t q;
+ R_t n_transformations;
+ q_t *transformations;
+ double T;
+ double *J;
+ double *H;
+ double *J_probs;
+ double *H_probs;
+ q_t *spins;
+ q_t *R;
+ double E;
+ v_t *M;
+} state_finite_t;
+
+v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot, gsl_rng *r);
+
diff --git a/lib/dihedral.c b/lib/dihedral.c
index 623041a..ac74a23 100644
--- a/lib/dihedral.c
+++ b/lib/dihedral.c
@@ -25,4 +25,16 @@ q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) {
}
}
+q_t *dihedral_gen_transformations(q_t q) {
+ q_t *transformations = (q_t *)malloc(q * q * sizeof(q_t));
+
+ for (q_t i = 0; i < q; i++) {
+ for (q_t j = 0; j < q; j++) {
+ transformations[q * i + j] = dihedral_act(q, i, j);
+ }
+ }
+
+ return transformations;
+}
+
diff --git a/lib/dihedral.h b/lib/dihedral.h
index 813e796..e5e4cbd 100644
--- a/lib/dihedral.h
+++ b/lib/dihedral.h
@@ -15,3 +15,5 @@ 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);
+q_t *dihedral_gen_transformations(q_t q);
+
diff --git a/lib/measurement.h b/lib/measurement.h
index eaa260b..46c034f 100644
--- a/lib/measurement.h
+++ b/lib/measurement.h
@@ -24,6 +24,11 @@ typedef struct {
double O2;
} autocorr_t;
+typedef struct {
+ void (*f)(state_finite_t *, void *);
+ void *data;
+} measurement_t;
+
void meas_update(meas_t *m, double x);
double meas_dx(const meas_t *m);
diff --git a/lib/symmetric.c b/lib/symmetric.c
new file mode 100644
index 0000000..729b38c
--- /dev/null
+++ b/lib/symmetric.c
@@ -0,0 +1,27 @@
+
+#include "symmetric.h"
+
+q_t *symmetric_compose(q_t q, const q_t *g1, const q_t *g2) {
+ q_t *g3 = (q_t *)malloc(q * sizeof(q_t));
+
+ for (q_t i = 0; i < q; i++) {
+ g3[i] = g1[g2[i]];
+ }
+
+ return g3;
+}
+
+q_t symmetric_act(const q_t *g, q_t s) {
+ return g[s];
+}
+
+q_t *symmetric_invert(q_t q, const q_t *g) {
+ q_t *g_inv = (q_t *)malloc(q * sizeof(q_t));
+
+ for (q_t i = 0; i < q; i++) {
+ g_inv[g[i]] = i;
+ }
+
+ return g_inv;
+}
+
diff --git a/lib/symmetric.h b/lib/symmetric.h
new file mode 100644
index 0000000..6e00f52
--- /dev/null
+++ b/lib/symmetric.h
@@ -0,0 +1,13 @@
+
+#pragma once
+
+#include <stdlib.h>
+
+#include "types.h"
+
+q_t *symmetric_compose(q_t q, const q_t *g1, const q_t *g2);
+
+q_t symmetric_act(const q_t *g, q_t s);
+
+q_t *symmetric_invert(q_t q, const q_t *g);
+
diff --git a/lib/types.h b/lib/types.h
index fcc2ce7..572ec64 100644
--- a/lib/types.h
+++ b/lib/types.h
@@ -3,6 +3,7 @@
typedef uint_fast32_t v_t;
typedef uint_fast8_t q_t;
+typedef uint_fast16_t R_t;
typedef uint_fast8_t D_t;
typedef uint_fast16_t L_t;
typedef uint_fast64_t count_t;
@@ -10,6 +11,7 @@ typedef int_fast64_t h_t;
#define MAX_v UINT_FAST32_MAX
#define MAX_Q UINT_FAST8_MAX
+#define MAX_R UINT_FAST16_MAX
#define MAX_D UINT_FAST8_MAX
#define MAX_L UINT_FAST16_MAX
#define MAX_COUNT UINT_FAST64_MAX
@@ -18,6 +20,7 @@ typedef int_fast64_t h_t;
#define PRIv PRIuFAST32
#define PRIq PRIuFAST8
+#define PRIR PRIuFAST16
#define PRID PRIuFAST8
#define PRIL PRIuFAST16
#define PRIcount PRIuFAST64
diff --git a/lib/wolff_finite.c b/lib/wolff_finite.c
new file mode 100644
index 0000000..64de9ba
--- /dev/null
+++ b/lib/wolff_finite.c
@@ -0,0 +1,70 @@
+
+#include "cluster_finite.h"
+
+void wolff_finite(state_finite_t *s, count_t sweeps, count_t sweeps_per_measurement, count_t n_measurements, measurement_t *measurements) {
+ for (count_t i = 0; i < sweeps; i++) {
+
+ count_t n_flips = 0;
+
+ while (n_flips / h->nv < sweeps_per_measurement) {
+ v_t v0 = gsl_rng_uniform_int(r, h->nv);
+ R_t step;
+
+ bool changed = false;
+ while (!changed) {
+ step = gsl_rng_uniform_int(r, s->n_transformations);
+ if v(symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) {
+ changed = true;
+ }
+ }
+
+ v_t tmp_flips = flip_cluster_finite(s, v0, step, r);
+ n_flips += tmp_flips;
+
+ if (n_runs > 0) {
+ n_steps++;
+ meas_update(clust, tmp_flips);
+
+ if (record_autocorrelation && n_steps % ac_skip == 0) {
+ update_autocorr(autocorr, s->E);
+ }
+
+ }
+
+ }
+
+ for (q_t i = 0; i < q; i++) {
+ meas_update(M[i], s->M[i]);
+ }
+ meas_update(E, s->E);
+
+ q_t n_at_max = 0;
+ q_t max_M_i = 0;
+ v_t max_M = 0;
+
+ for (q_t i = 0; i < q; i++) {
+ if (s->M[i] > max_M) {
+ n_at_max = 1;
+ max_M_i = i;
+ max_M = s->M[i];
+ } else if (s->M[i] == max_M) {
+ n_at_max++;
+ }
+ }
+
+ if (record_distribution) {
+ mag_dist[s->M[0]]++;
+ }
+
+ if (n_at_max == 1) {
+ for (q_t i = 0; i < q; i++) {
+ meas_update(sM[max_M_i][i], s->M[i]);
+ }
+ meas_update(sE[max_M_i], s->E);
+ freqs[max_M_i]++;
+ }
+
+ diff = fabs(meas_dx(clust) / clust->x);
+ }
+}
+