summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/dihedral.c28
-rw-r--r--lib/dihedral.h17
-rw-r--r--lib/orthogonal.c95
-rw-r--r--lib/orthogonal.h24
-rw-r--r--lib/wolff.h17
-rw-r--r--lib/wolff_tools.c137
6 files changed, 300 insertions, 18 deletions
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++;
}