From 2d8fcebf2f56efd1c3913ba49eaff6520ffdb33d Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 6 Jul 2018 14:42:44 -0400 Subject: rewrote wolff in c++ with templates so that any system can be run with it --- lib/cluster.c | 277 ---------------------------------------------------------- 1 file changed, 277 deletions(-) delete mode 100644 lib/cluster.c (limited to 'lib/cluster.c') diff --git a/lib/cluster.c b/lib/cluster.c deleted file mode 100644 index 96225a2..0000000 --- a/lib/cluster.c +++ /dev/null @@ -1,277 +0,0 @@ - -#include "cluster.h" - -v_t flip_cluster_dgm(dgm_state_t *s, v_t v0, h_t 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 - - bool *marks = (bool *)calloc(s->g->nv, sizeof(bool)); - - while (stack != NULL) { - v_t v = stack_pop(&stack); - - if (!marks[v]) { - h_t s_old, s_new; - dihinf_t *R_new; - bool external_flipped; - - marks[v] = true; - - if (v == s->g->nv - 1) { - R_new = dihinf_compose(rot, s->R); - external_flipped = true; - } else { - s_old = s->spins[v]; - s_new = dihinf_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++) { - h_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) { - h_t rot_s_old, rot_s_new; - - if (external_neighbor) { - rot_s_old = dihinf_inverse_act(s->R, s_old); - rot_s_new = dihinf_inverse_act(s->R, s_new); - } else { - rot_s_old = dihinf_inverse_act(s->R, sn); - rot_s_new = dihinf_inverse_act(R_new, sn); - } - - double dE = s->H(s->H_info, rot_s_old) - s->H(s->H_info, rot_s_new); - prob = 1.0 - exp(-dE / s->T); - - s->M += rot_s_new - rot_s_old; - s->E += dE; - } else { - double dE = (s->J)(s_old - sn) - (s->J)(s_new - sn); - prob = 1.0 - exp(-dE / s->T); - 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 (external_flipped) { - free(s->R); - s->R = R_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); - - 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]) { - bool v_is_external = false; - double *s_old, *s_new, *R_tmp; - - if (v == s->g->nv - 1) { - v_is_external = true; - } - - //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]; - - bool vn_is_external = false; - - if (vn == s->g->nv - 1) { - vn_is_external = true; - } - - double *sn; - - if (!vn_is_external) { - sn = &(s->spins[s->n * vn]); - } - - double prob; - - if (v_is_external || vn_is_external) { - double *rs_old, *rs_new; - if (vn_is_external) { - 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(s->n, s->H_info, rs_old) - s->H(s->n, s->H_info, 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); - 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++; - } - } - } - - //tree_freeNode(T); - free(marks); - - return nv; -} - -/*G -template -v_t flip_cluster(state_t *state, v_t v0, R_t *r, gsl_rng *rand) { - 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(state->g->nv, sizeof(bool)); - - while (stack != NULL) { - v_t v = stack_pop(&stack); - - if (!marks[v]) { - X_t *si_old, *si_new; - R_t *s0_old, *s0_new; - - si_old = state->s[v]; - s0_old = state->s0; - - marks[v] = true; - - if (v == state->g->nv - 1) { - s0_new = act (r, s0_old); - } else { - si_new = act (r, si_old); - } - - v_t nn = state->g->v_i[v + 1] - state->g->v_i[v]; - - for (v_t i = 0; i < nn; i++) { - v_t vn = state->g->v_adj[state->g->v_i[v] + i]; - - X_t *sj; - - if (vn != state->g->nv - 1) { - sj = state->s[vn]; - } - - double prob; - - bool is_ext = (v == state->g->nv - 1 || vn == state->g->nv - 1); - - if (is_ext) { - X_t *rs_old, *rs_new; - if (vn == state->g->nv - 1) { - rs_old = inverse_act (s0_old, si_old); - rs_new = inverse_act (s0_old, si_new); - } else { - rs_old = inverse_act (s0_old, sj); - rs_new = inverse_act (s0_new, sj); - } - double dE = state->B(rs_old) - state->B(rs_new); - prob = 1.0 - exp(-dE / state->T); - update_magnetization (state->M, rs_old, rs_new); - state->E += dE; - - free_X (rs_old); - free_X (rs_new); - } else { - double dE = state->Z(si_old, sj) - state->Z(si_new, sj); - prob = 1.0 - exp(-dE / state->T); - state->E += dE; - } - - if (gsl_rng_uniform(rand) < prob) { // and with probability... - stack_push(&stack, vn); // push the neighboring vertex to the stack - } - } - - if (v == state->g->nv - 1) { - free_R (state->s0); - state->s0 = s0_new; - } else { - free_X (state->s[v]); - state->s[v] = si_new; - } - - if (v != state->g->nv - 1) { // count the number of non-external sites that flip - nv++; - } - } - } - - free(marks); - - return nv; -} -*/ -- cgit v1.2.3-70-g09d2