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.h | 391 ++++++++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 337 insertions(+), 54 deletions(-) (limited to 'lib/cluster.h') diff --git a/lib/cluster.h b/lib/cluster.h index d118735..29dd0cb 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -1,13 +1,14 @@ #pragma once +#include #include #include #include #include #include #include -#include +#include #include #include #include @@ -24,57 +25,339 @@ #include "dihinf.h" #include "yule_walker.h" -typedef struct { - graph_t *g; - q_t *spins; - double T; - double *J; - 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; - h_t *spins; - double T; - double (*J)(h_t); - double (*H)(double *, h_t); - double *H_info; - dihinf_t *R; - double E; - h_t M; -} dgm_state_t; - -typedef struct { - graph_t *g; - double *spins; - double T; - double (*J)(double); - double (*H)(q_t, double *, double *); - double *H_info; - double *R; - double E; - double *M; - q_t n; -} vector_state_t; - -typedef enum { - VECTOR, - MODULATED, - CUBIC, - QUADRATIC -} vector_field_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); - -v_t flip_cluster_dgm(dgm_state_t *s, v_t v0, h_t rot, gsl_rng *r); - -graph_t *graph_add_ext(const graph_t *g); +template +void init(T*); + +template +T scalar_multiple(v_t a, T b); + +template +X_t act(R_t a, X_t b); + +template +X_t act_inverse(R_t a, X_t b); + +template +T copy(T a); + +template +void free_spin(T a); + +template +T add(T, T); + +template +T subtract(T, T); + +template +T gen_rot(gsl_rng *r); + +template +class state_t { + public: + D_t D; + L_t L; + v_t nv; + v_t ne; + graph_t *g; + double T; + X_t *spins; + R_t R; + double E; + X_t M; // the "sum" of the spins, like the total magnetization + + std::function J; + std::function H; + + state_t(D_t D, L_t L, double T, std::function J, std::function H) : D(D), L(L), T(T), J(J), H(H) { + graph_t *h = graph_create_square(D, L); + nv = h->nv; + ne = h->ne; + g = graph_add_ext(h); + graph_free(h); + spins = (X_t *)malloc(nv * sizeof(X_t)); + for (v_t i = 0; i < nv; i++) { + init (&(spins[i])); + } + init (&R); + E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]); + M = scalar_multiple (nv, spins[0]); + } + + ~state_t() { + graph_free(g); + for (v_t i = 0; i < nv; i++) { + free_spin(spins[i]); + } + free(spins); + free_spin(R); + free_spin(M); + } +}; + +template +struct vector_t { T *x; }; + +template +void init(vector_t *ptr) { + ptr->x = (T *)calloc(q, sizeof(T)); + + ptr->x[0] = (T)1; +} + +template +vector_t copy (vector_t v) { + vector_t v_copy; + + v_copy.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + v_copy.x[i] = v.x[i]; + } + + return v_copy; +} + +template +void add (vector_t v1, vector_t v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] += v2.x[i]; + } +} + +template +void subtract (vector_t v1, vector_t v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] -= v2.x[i]; + } +} + +template +vector_t scalar_multiple(v_t a, vector_t v) { + vector_t multiple; + multiple.x = (T *)malloc(q * sizeof(T)); + for (q_t i = 0; i < q; i++) { + multiple.x[i] = a * v.x[i]; + } + + return multiple; +} + +template +T dot(vector_t v1, vector_t v2) { + T prod = 0; + + for (q_t i = 0; i < q; i++) { + prod += v1.x[i] * v2.x[i]; + } + + return prod; +} + +template +void free_spin (vector_t v) { + free(v.x); +} + +template +struct orthogonal_t { T *x; }; + +template +void init(orthogonal_t *ptr) { + ptr->x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + ptr->x[q * i + i] = (T)1; + } +} + +template +orthogonal_t copy (orthogonal_t m) { + orthogonal_t m_copy; + m_copy.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q * q; i++) { + m_copy.x[i] = m.x[i]; + } + + return m_copy; +} + +template +void free_spin (orthogonal_t m) { + free(m.x); +} + +template +vector_t act (orthogonal_t m, vector_t v) { + vector_t v_rot; + v_rot.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot.x[i] += m.x[q * i + j] * v.x[j]; + } + } + + return v_rot; +} + + +template +orthogonal_t act (orthogonal_t m1, orthogonal_t m2) { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k]; + } + } + } + + return m2_rot; +} + +template +vector_t act_inverse (orthogonal_t m, vector_t v) { + vector_t v_rot; + v_rot.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot.x[i] += m.x[q * j + i] * v.x[j]; + } + } + + return v_rot; +} + +template +orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k]; + } + } + } + + return m2_rot; +} + +template +void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { + double *v = (double *)malloc(q * sizeof(double)); + double v2 = 0; + + for (q_t i = 0; i < q; i++) { + v[i] = gsl_ran_ugaussian(r); + v2 += v[i] * v[i]; + } + + ptr->x = (double *)calloc(q * q, sizeof(double)); + + for (q_t i = 0; i < q; i++) { + ptr->x[q * i + i] = 1.0; + for (q_t j = 0; j < q; j++) { + ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2; + } + } + + free(v); +} + +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 R_old, R_new; + + si_old = state->spins[v]; + R_old = state->R; + + marks[v] = true; + + if (v == state->g->nv - 1) { + R_new = act (r, R_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->spins[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 = act_inverse (R_old, si_old); + rs_new = act_inverse (R_old, si_new); + } else { + rs_old = act_inverse (R_old, sj); + rs_new = act_inverse (R_new, sj); + } + double dE = state->H(rs_old) - state->H(rs_new); + prob = 1.0 - exp(-dE / state->T); + + subtract (state->M, rs_old); + add (state->M, rs_new); + state->E += dE; + + free_spin (rs_old); + free_spin (rs_new); + } else { + double dE = state->J(si_old, sj) - state->J(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_spin(state->R); + state->R = R_new; + } else { + free_spin(state->spins[v]); + state->spins[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