From 8a6109a78a6f7b9e9585cd89e8a10d1f626b3af1 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 17 Jul 2018 22:37:57 -0400 Subject: made many changes to implement alternate reflection generators, this is broken --- lib/cluster.h | 60 +++++++------------------------------------------------- lib/orthogonal.h | 54 ++++++++++++++++++++++++++++++++++++++++++++------ lib/state.h | 53 +++++++++++++++++++++++++++++++++++++++++++++++++ lib/wolff.h | 8 ++++---- 4 files changed, 112 insertions(+), 63 deletions(-) create mode 100644 lib/state.h (limited to 'lib') diff --git a/lib/cluster.h b/lib/cluster.h index a20912e..ccd8725 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -13,6 +13,7 @@ #include #include +#include "state.h" #include "types.h" #include "rand.h" #include "stack.h" @@ -50,52 +51,6 @@ 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 v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { v_t nv = 0; @@ -103,7 +58,7 @@ v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { 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)); + bool *marks = (bool *)calloc(state->nv, sizeof(bool)); while (stack != NULL) { v_t v = stack_pop(&stack); @@ -112,14 +67,13 @@ v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { 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) { + if (v == state->nv) { + R_old = state->R; R_new = act (r, R_old); } else { + si_old = state->spins[v]; si_new = act (r, si_old); } @@ -130,13 +84,13 @@ v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { X_t sj; - if (vn != state->g->nv - 1) { + if (vn != state->nv) { sj = state->spins[vn]; } double prob; - bool is_ext = (v == state->g->nv - 1 || vn == state->g->nv - 1); + bool is_ext = (v == state->nv || vn == state->nv); if (is_ext) { X_t rs_old, rs_new; diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 85f7a20..6166223 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -6,6 +6,7 @@ #include #include +#include "state.h" #include "types.h" template @@ -143,21 +144,62 @@ orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) } template -void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { - ptr->is_reflection = true; - ptr->x = (double *)calloc(q, sizeof(double)); +orthogonal_t generate_rotation_uniform (gsl_rng *r, const state_t , vector_t > *s) { + orthogonal_t ptr; + ptr.is_reflection = true; + ptr.x = (double *)calloc(q, sizeof(double)); double v2 = 0; for (q_t i = 0; i < q; i++) { - ptr->x[i] = gsl_ran_ugaussian(r); - v2 += ptr->x[i] * ptr->x[i]; + ptr.x[i] = gsl_ran_ugaussian(r); + v2 += ptr.x[i] * ptr.x[i]; } double mag_v = sqrt(v2); for (q_t i = 0; i < q; i++) { - ptr->x[i] /= mag_v; + ptr.x[i] /= mag_v; } + + return ptr; +} + +template +orthogonal_t generate_rotation_perturbation (gsl_rng *r, const state_t , vector_t > *s, double epsilon) { + orthogonal_t ptr; + vector_t tmp_v; + ptr.is_reflection = true; + + tmp_v.x = (double *)malloc(q * sizeof(double)); + + double tmp2 = 0; + double M2 = 0; + double tmpM = 0; + + for (q_t i = 0; i < q; i++) { + tmp_v.x[i] = gsl_ran_ugaussian(r); + M2 += pow(s->M.x[i], 2); + tmpM += tmp_v.x[i] * s->M.x[i]; + } + + double v2 = 0; + + for (q_t i = 0; i < q; i++) { + tmp_v.x[i] = (tmp_v.x[i] - tmpM * s->M.x[i] / M2) + epsilon * gsl_ran_ugaussian(r); + v2 += tmp_v.x[i] * tmp_v.x[i]; + } + + double mag_v = sqrt(v2); + + for (q_t i = 0; i < q; i++) { + tmp_v.x[i] /= mag_v; + } + + vector_t v = act(s->R, tmp_v); + free(tmp_v.x); + ptr.x = v.x; + + return ptr; } diff --git a/lib/state.h b/lib/state.h new file mode 100644 index 0000000..8d4265b --- /dev/null +++ b/lib/state.h @@ -0,0 +1,53 @@ + +#pragma once + +#include + +#include "types.h" +#include "graph.h" + +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); + } +}; + + + diff --git a/lib/wolff.h b/lib/wolff.h index caf413b..90470d5 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -2,7 +2,8 @@ #include #include -#include +#include "cluster.h" +#include "state.h" template double H_vector(vector_t v1, T *H) { @@ -12,7 +13,7 @@ double H_vector(vector_t v1, T *H) { } template -void wolff(count_t N, D_t D, L_t L, double T, std::function J, std::function H, unsigned long timestamp, bool silent) { +void wolff(count_t N, D_t D, L_t L, double T, std::function J, std::function H, std::function *)> gen_R, unsigned long timestamp, bool silent) { state_t s(D, L, T, J, H); @@ -44,8 +45,7 @@ void wolff(count_t N, D_t D, L_t L, double T, std::function J v_t v0 = gsl_rng_uniform_int(r, s.nv); - R_t step; - generate_rotation(r, &step); + R_t step = gen_R(r, &s); cluster_size = flip_cluster (&s, v0, step, r); -- cgit v1.2.3-70-g09d2