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 From 31f4244352b5e68eed770090419541d469f7f999 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 6 Jul 2018 14:51:29 -0400 Subject: split up some files --- lib/cluster.h | 177 +------------------------------------------------------ lib/orthogonal.c | 99 ------------------------------- lib/orthogonal.h | 116 +++++++++++++++++++++++++++++++++--- lib/vector.h | 72 ++++++++++++++++++++++ 4 files changed, 180 insertions(+), 284 deletions(-) delete mode 100644 lib/orthogonal.c create mode 100644 lib/vector.h (limited to 'lib/cluster.h') diff --git a/lib/cluster.h b/lib/cluster.h index 29dd0cb..a20912e 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -20,6 +20,7 @@ #include "graph.h" #include "tree.h" #include "measurement.h" +#include "vector.h" #include "orthogonal.h" #include "dihedral.h" #include "dihinf.h" @@ -95,182 +96,6 @@ class state_t { } }; -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; diff --git a/lib/orthogonal.c b/lib/orthogonal.c deleted file mode 100644 index 87569ae..0000000 --- a/lib/orthogonal.c +++ /dev/null @@ -1,99 +0,0 @@ - -#include "orthogonal.h" - -void vector_replace(q_t n, double *v1, const double *v2) { - // writes vector v2 of length n to memory located at v1 - for (q_t i = 0; i < n; i++) { - v1[i] = v2[i]; - } -} - -void vector_add(q_t n, double *v1, const double *v2) { - // adds vector v2 of length n to vector v1 - for (q_t i = 0; i < n; i++) { - v1[i] += v2[i]; - } -} - -void vector_subtract(q_t n, double *v1, const double *v2) { - // subtracts vector v2 of length n from vector v1 - for (q_t i = 0; i < n; i++) { - v1[i] -= v2[i]; - } -} - -double *vector_rotate(q_t n, const double *rot, const double *vec) { - // multiplies n by n rotation matrix rot to vector 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, const double *v1, const 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, const double *r, const 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 index 60d5f49..0b2fdd5 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -1,24 +1,122 @@ +#pragma once + #include #include #include -#include +#include #include "types.h" -void vector_replace(q_t n, double *v1, const double *v2); +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)); -void vector_add(q_t n, double *v1, const double *v2); + 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]; + } + } + } -void vector_subtract(q_t n, double *v1, const double *v2); + return m2_rot; +} -double *vector_rotate(q_t n, const double *rot, const double *vec); +template +void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { + double *v = (double *)malloc(q * sizeof(double)); + double v2 = 0; -double *vector_rotate_inverse(q_t n, const double *rot, const double *vec); + for (q_t i = 0; i < q; i++) { + v[i] = gsl_ran_ugaussian(r); + v2 += v[i] * v[i]; + } -double vector_dot(q_t n, const double *v1, const double *v2); + 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; + } + } -double *orthogonal_rotate(q_t n, const double *m1, const double *m2); + free(v); +} -double *gen_rot(gsl_rng *r, q_t n); diff --git a/lib/vector.h b/lib/vector.h new file mode 100644 index 0000000..c7f459c --- /dev/null +++ b/lib/vector.h @@ -0,0 +1,72 @@ + +#pragma once + +#include +#include + +#include "types.h" + +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); +} + -- cgit v1.2.3-70-g09d2