diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/cluster.c | 277 | ||||
| -rw-r--r-- | lib/cluster.h | 391 | ||||
| -rw-r--r-- | lib/graph.h | 10 | ||||
| -rw-r--r-- | lib/rand.h | 9 | ||||
| -rw-r--r-- | lib/stack.h | 9 |
5 files changed, 363 insertions, 333 deletions
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 <class R_t, class X_t> -v_t flip_cluster(state_t <R_t, X_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_t, R_t> (r, s0_old); - } else { - si_new = act <R_t, X_t> (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 <class R_t, class X_t> (s0_old, si_old); - rs_new = inverse_act <class R_t, class X_t> (s0_old, si_new); - } else { - rs_old = inverse_act <class R_t, class X_t> (s0_old, sj); - rs_new = inverse_act <class R_t, class X_t> (s0_new, sj); - } - double dE = state->B(rs_old) - state->B(rs_new); - prob = 1.0 - exp(-dE / state->T); - update_magnetization <X_t> (state->M, rs_old, rs_new); - state->E += dE; - - free_X <X_t> (rs_old); - free_X <X_t> (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 <R_t> (state->s0); - state->s0 = s0_new; - } else { - free_X <X_t> (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; -} -*/ 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 <functional> #include <assert.h> #include <fftw3.h> #include <float.h> #include <gsl/gsl_randist.h> #include <gsl/gsl_rng.h> #include <inttypes.h> -#include <math.h> +#include <cmath> #include <stdbool.h> #include <string.h> #include <sys/types.h> @@ -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 <class T> +void init(T*); + +template <class T> +T scalar_multiple(v_t a, T b); + +template <class R_t, class X_t> +X_t act(R_t a, X_t b); + +template <class R_t, class X_t> +X_t act_inverse(R_t a, X_t b); + +template <class T> +T copy(T a); + +template <class T> +void free_spin(T a); + +template <class T> +T add(T, T); + +template <class T> +T subtract(T, T); + +template <class T> +T gen_rot(gsl_rng *r); + +template <class R_t, class X_t> +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 <double(X_t, X_t)> J; + std::function <double(X_t)> H; + + state_t(D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> 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 <q_t q, class T> +struct vector_t { T *x; }; + +template <q_t q, class T> +void init(vector_t <q, T> *ptr) { + ptr->x = (T *)calloc(q, sizeof(T)); + + ptr->x[0] = (T)1; +} + +template <q_t q, class T> +vector_t <q, T> copy (vector_t <q, T> v) { + vector_t <q, 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 <q_t q, class T> +void add (vector_t <q, T> v1, vector_t <q, T> v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] += v2.x[i]; + } +} + +template <q_t q, class T> +void subtract (vector_t <q, T> v1, vector_t <q, T> v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] -= v2.x[i]; + } +} + +template <q_t q, class T> +vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) { + vector_t <q, 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 <q_t q, class T> +T dot(vector_t <q, T> v1, vector_t <q, T> v2) { + T prod = 0; + + for (q_t i = 0; i < q; i++) { + prod += v1.x[i] * v2.x[i]; + } + + return prod; +} + +template <q_t q, class T> +void free_spin (vector_t <q, T> v) { + free(v.x); +} + +template <q_t q, class T> +struct orthogonal_t { T *x; }; + +template <q_t q, class T> +void init(orthogonal_t <q, 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 <q_t q, class T> +orthogonal_t <q, T> copy (orthogonal_t <q, T> m) { + orthogonal_t <q, 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 <q_t q, class T> +void free_spin (orthogonal_t <q, T> m) { + free(m.x); +} + +template <q_t q, class T> +vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) { + vector_t <q, 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 <q_t q, class T> +orthogonal_t <q, T> act (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) { + orthogonal_t <q, 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 <q_t q, class T> +vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) { + vector_t <q, 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 <q_t q, class T> +orthogonal_t <q, T> act_inverse (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) { + orthogonal_t <q, 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 <q_t q> +void generate_rotation (gsl_rng *r, orthogonal_t <q, double> *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 <class R_t, class X_t> +v_t flip_cluster(state_t <R_t, X_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; +} diff --git a/lib/graph.h b/lib/graph.h index cb47faa..beb7f4c 100644 --- a/lib/graph.h +++ b/lib/graph.h @@ -7,6 +7,10 @@ #include "types.h" +#ifdef __cplusplus +extern "C" { +#endif + typedef struct { v_t ne; v_t nv; @@ -15,8 +19,10 @@ typedef struct { } graph_t; graph_t *graph_create_square(D_t D, L_t L); - graph_t *graph_add_ext(const graph_t *G); - void graph_free(graph_t *h); +#ifdef __cplusplus +} +#endif + @@ -4,4 +4,13 @@ #include <assert.h> #include <stdio.h> +#ifdef __cplusplus +extern "C" { +#endif + unsigned long int rand_seed(); + +#ifdef __cplusplus +} +#endif + diff --git a/lib/stack.h b/lib/stack.h index a354ab5..8d25aff 100644 --- a/lib/stack.h +++ b/lib/stack.h @@ -8,6 +8,11 @@ #include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + typedef struct ll_tag { v_t x; struct ll_tag *next; @@ -24,3 +29,7 @@ void stack_push_d(dll_t **q, double x); v_t stack_pop(ll_t **q); double stack_pop_d(dll_t **q); +#ifdef __cplusplus +} +#endif + |