From 2928df0a4b4be30cde2b11bdcf2698875d5b9536 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 20 Mar 2018 19:06:50 -0400 Subject: new system --- lib/cluster.c | 87 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ lib/cluster.h | 15 +++++++++++ lib/dihinf.c | 28 +++++++++++++++++++ lib/dihinf.h | 17 ++++++++++++ lib/types.h | 14 +++++----- 5 files changed, 155 insertions(+), 6 deletions(-) create mode 100644 lib/dihinf.c create mode 100644 lib/dihinf.h (limited to 'lib') diff --git a/lib/cluster.c b/lib/cluster.c index c800a0f..fa417fa 100644 --- a/lib/cluster.c +++ b/lib/cluster.c @@ -92,6 +92,93 @@ v_t flip_cluster(ising_state_t *s, v_t v0, q_t rot, gsl_rng *r) { return nv; } +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; diff --git a/lib/cluster.h b/lib/cluster.h index 095ed61..2de17e5 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -21,6 +21,7 @@ #include "measurement.h" #include "orthogonal.h" #include "dihedral.h" +#include "dihinf.h" typedef struct { graph_t *g; @@ -36,6 +37,18 @@ typedef struct { 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; @@ -53,5 +66,7 @@ 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); diff --git a/lib/dihinf.c b/lib/dihinf.c new file mode 100644 index 0000000..4f88a7a --- /dev/null +++ b/lib/dihinf.c @@ -0,0 +1,28 @@ + +#include "dihinf.h" + +dihinf_t *dihinf_compose(h_t g1i, const dihinf_t *g2) { + // we only need to consider the action of reflections + dihinf_t *g3 = (dihinf_t *)malloc(1 * sizeof(dihinf_t)); + + g3->r = !g2->r; + g3->i = g1i - g2->i; + + return g3; +} + +h_t dihinf_act(h_t gi, h_t s) { + // we only need to consider the action of reflections + + return gi - s; +} + +h_t dihinf_inverse_act(const dihinf_t *g, h_t s) { + if (g->r) { + return g->i - s; + } else { + return s - g->i; + } +} + + diff --git a/lib/dihinf.h b/lib/dihinf.h new file mode 100644 index 0000000..2bc7dc2 --- /dev/null +++ b/lib/dihinf.h @@ -0,0 +1,17 @@ + +#include +#include + +#include "types.h" + +typedef struct { + h_t i; + bool r; +} dihinf_t; + +dihinf_t *dihinf_compose(h_t gti, const dihinf_t *g2); + +h_t dihinf_act(h_t gi, h_t s); + +h_t dihinf_inverse_act(const dihinf_t *g, h_t s); + diff --git a/lib/types.h b/lib/types.h index daac873..fcc2ce7 100644 --- a/lib/types.h +++ b/lib/types.h @@ -8,16 +8,18 @@ typedef uint_fast16_t L_t; typedef uint_fast64_t count_t; typedef int_fast64_t h_t; -#define MAX_v INT_FAST32_MAX -#define MAX_Q INT_FAST8_MAX -#define MAX_D INT_FAST8_MAX -#define MAX_L INT_FAST16_MAX -#define MAX_COUNT INT_FAST64_MAX +#define MAX_v UINT_FAST32_MAX +#define MAX_Q UINT_FAST8_MAX +#define MAX_D UINT_FAST8_MAX +#define MAX_L UINT_FAST16_MAX +#define MAX_COUNT UINT_FAST64_MAX +#define MAX_H INT_FAST64_MAX +#define MIN_H INT_FAST64_MIN #define PRIv PRIuFAST32 #define PRIq PRIuFAST8 #define PRID PRIuFAST8 #define PRIL PRIuFAST16 #define PRIcount PRIuFAST64 -#define PRIh PRIFAST64 +#define PRIh PRIdFAST64 -- cgit v1.2.3-70-g09d2