From 1fbcb4dd4e52daeb53becba33827f8e48c5606b2 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 28 Feb 2018 20:33:41 -0500 Subject: fixed major mistake in process, also got n-component version fully working --- lib/dihedral.c | 28 +++++++++++ lib/dihedral.h | 17 +++++++ lib/orthogonal.c | 95 +++++++++++++++++++++++++++++++++++++ lib/orthogonal.h | 24 ++++++++++ lib/wolff.h | 17 +++++++ lib/wolff_tools.c | 137 +++++++++++++++++++++++++++++++++++++++++++++++------- 6 files changed, 300 insertions(+), 18 deletions(-) create mode 100644 lib/dihedral.c create mode 100644 lib/dihedral.h create mode 100644 lib/orthogonal.c create mode 100644 lib/orthogonal.h (limited to 'lib') 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 +#include + +#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 +#include +#include +#include + +#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++; } -- cgit v1.2.3-70-g09d2