From c48fd16fe1554c88c79a1f0d50e81c803da8f61f Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 24 Jul 2018 13:12:35 -0400 Subject: implemented updating the first fourier moment in wolff_finite, but also make wolff_finite obselete by adding a hacky preprocessor method for making wolff.h as efficient --- lib/cluster.h | 8 ++++++++ lib/cluster_finite.c | 14 +++++++++++++- lib/cluster_finite.h | 4 ++++ lib/finite_states.h | 33 +++++++++++++++++++++++++++++++++ lib/initial_finite.c | 48 +++++++++++++++++++++++++++++++++++++++++++++++- lib/wolff.h | 4 ++++ src/wolff_ising.cpp | 5 +++++ src/wolff_potts.cpp | 5 +++++ 8 files changed, 119 insertions(+), 2 deletions(-) create mode 100644 lib/finite_states.h diff --git a/lib/cluster.h b/lib/cluster.h index 5124061..b8c98e5 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -64,7 +64,11 @@ void flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { non_ghost = vn; } double dE = state->H(rs_old) - state->H(rs_new); +#ifdef FINITE_STATES + prob = H_probs[N_STATES * state_to_ind(rs_old) + state_to_ind(rs_new)]; +#else prob = 1.0 - exp(-dE / state->T); +#endif add(&(state->M), -1, rs_old); add(&(state->M), 1, rs_new); @@ -84,7 +88,11 @@ void flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { free_spin (rs_new); } else { double dE = state->J(si_old, sj) - state->J(si_new, sj); +#ifdef FINITE_STATES + prob = J_probs[N_STATES * N_STATES * state_to_ind(si_old) + N_STATES * state_to_ind(si_new) + state_to_ind(sj)]; +#else prob = 1.0 - exp(-dE / state->T); +#endif state->E += dE; } diff --git a/lib/cluster_finite.c b/lib/cluster_finite.c index 9392cf8..2041df0 100644 --- a/lib/cluster_finite.c +++ b/lib/cluster_finite.c @@ -34,7 +34,7 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) { v_t nn = s->g->v_i[v + 1] - s->g->v_i[v]; for (v_t i = 0; i < nn; i++) { - q_t sn; + q_t sn, non_ghost; double prob; bool external_neighbor = false; @@ -42,8 +42,10 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) { if (vn == s->g->nv - 1) { external_neighbor = true; + non_ghost = v; } else { sn = s->spins[vn]; + non_ghost = vn; } if (external_flipped || external_neighbor) { @@ -62,6 +64,16 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) { s->M[rot_s_old]--; s->M[rot_s_new]++; + for (D_t i = 0; i < s->D; i++) { + L_t x = (non_ghost / (v_t)pow(s->L, s->D - i - 1)) % s->L; + + s->ReF[s->D * i + rot_s_old] -= s->precomputed_cos[i]; + s->ReF[s->D * i + rot_s_new] += s->precomputed_cos[i]; + + s->ImF[s->D * i + rot_s_old] -= s->precomputed_sin[i]; + s->ImF[s->D * i + rot_s_new] += s->precomputed_sin[i]; + } + } else { q_t diff_old = s->bond_with_zero_type[s->transformations[s->q * s->transform_site_to_zero[sn] + s_old]]; q_t diff_new = s->bond_with_zero_type[s->transformations[s->q * s->transform_site_to_zero[sn] + s_new]]; diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h index b2d764e..117087e 100644 --- a/lib/cluster_finite.h +++ b/lib/cluster_finite.h @@ -45,6 +45,10 @@ typedef struct { q_t *R; v_t *B; v_t *M; + double *ReF; + double *ImF; + double *precomputed_cos; + double *precomputed_sin; } state_finite_t; v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot, gsl_rng *r); diff --git a/lib/finite_states.h b/lib/finite_states.h new file mode 100644 index 0000000..1e8800c --- /dev/null +++ b/lib/finite_states.h @@ -0,0 +1,33 @@ + +#pragma once + +#include +#include + +#define FINITE_STATES + +// must have N_STATES, states[N_STATES], and state_to_ind defined before +// invoking header + +double J_probs[N_STATES * N_STATES * N_STATES]; +double H_probs[N_STATES * N_STATES]; + +template +void initialize_probs(std::function J, std::function H, double T) { + for (q_t i = 0; i < N_STATES; i++) { + for (q_t j = 0; j < N_STATES; j++) { + for (q_t k = 0; k < N_STATES; k++) { + J_probs[i * N_STATES * N_STATES + j * N_STATES +k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T); + } + } + } + for (q_t i = 0; i < N_STATES; i++) { + for (q_t j = 0; j < N_STATES; j++) { + H_probs[i * N_STATES + j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T); + } + } +} + + + + diff --git a/lib/initial_finite.c b/lib/initial_finite.c index 6ea76ef..59e7ec4 100644 --- a/lib/initial_finite.c +++ b/lib/initial_finite.c @@ -1,6 +1,24 @@ #include "initial_finite.h" +double *precompute_cos(L_t L) { + double *x = (double *)malloc(L * sizeof(double)); + for (L_t i = 0; i < L; i++) { + x[i] = cos(2 * M_PI * (double)i / (double)L); + } + + return x; +} + +double *precompute_sin(L_t L) { + double *x = (double *)malloc(L * sizeof(double)); + for (L_t i = 0; i < L; i++) { + x[i] = sin(2 * M_PI * (double)i / (double)L); + } + + return x; +} + double *Jprobs_from_J(q_t q, double T, double *J) { double *J_probs = (double *)calloc(pow(q, 2), sizeof(double)); @@ -112,6 +130,13 @@ state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H) s->B = (v_t *)calloc(2, sizeof(v_t)); s->B[0] = s->ne; + s->ReF = (double *)calloc(D * 2, sizeof(double)); + s->ImF = (double *)calloc(D * 2, sizeof(double)); + + + s->precomputed_cos = precompute_cos(L); + s->precomputed_sin = precompute_sin(L); + return s; } @@ -168,6 +193,12 @@ state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, doub s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); s->B[0] = s->ne; // everyone starts in state 0, remember? + s->ReF = (double *)calloc(D * q, sizeof(double)); + s->ImF = (double *)calloc(D * q, sizeof(double)); + + s->precomputed_cos = precompute_cos(L); + s->precomputed_sin = precompute_sin(L); + return s; } @@ -230,10 +261,15 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); s->B[0] = s->ne; // everyone starts in state 0, remember? + s->ReF = (double *)calloc(D * q, sizeof(double)); + s->ImF = (double *)calloc(D * q, sizeof(double)); + + s->precomputed_cos = precompute_cos(L); + s->precomputed_sin = precompute_sin(L); + return s; } - state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double *H) { state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t)); @@ -291,6 +327,12 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); s->B[0] = s->nv; // everyone starts in state 0, remember? + s->ReF = (double *)calloc(D * q, sizeof(double)); + s->ImF = (double *)calloc(D * q, sizeof(double)); + + s->precomputed_cos = precompute_cos(L); + s->precomputed_sin = precompute_sin(L); + return s; } @@ -321,6 +363,10 @@ void state_finite_free(state_finite_t *s) { free(s->involutions); free(s->transform_site_to_zero); free(s->bond_with_zero_type); + free(s->ReF); + free(s->ImF); + free(s->precomputed_cos); + free(s->precomputed_sin); free(s); } diff --git a/lib/wolff.h b/lib/wolff.h index 8286024..da4b7b6 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -5,6 +5,10 @@ template void wolff(count_t N, state_t *s, std::function gen_R, std::function *)> measurements, gsl_rng *r, bool silent) { +#ifdef FINITE_STATES + initialize_probs(s->J, s->H, s->T); +#endif + if (!silent) printf("\n"); for (count_t steps = 0; steps < N; steps++) { if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s->E, s->last_cluster_size); diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp index 9866248..f614707 100644 --- a/src/wolff_ising.cpp +++ b/src/wolff_ising.cpp @@ -8,6 +8,11 @@ #include #include +#define N_STATES 2 +const ising_t states[N_STATES] = {false, true}; +q_t state_to_ind(ising_t state) { return (q_t)state.x; } +#include + // include wolff.h #include #include diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp index e3259e4..6b6f602 100644 --- a/src/wolff_potts.cpp +++ b/src/wolff_potts.cpp @@ -10,6 +10,11 @@ #include #include +#define N_STATES POTTSQ +const potts_t states[8] = {0, 1, 2, 3, 4, 5, 6, 7}; +q_t state_to_ind(potts_t state) { return (q_t)state.x; } +#include + // include wolff.h #include #include -- cgit v1.2.3-70-g09d2