summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-24 13:12:35 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-24 13:12:35 -0400
commitc48fd16fe1554c88c79a1f0d50e81c803da8f61f (patch)
tree2b8a07de4e0fc795854fc4c1ac0ff095783218fa /lib
parent8d96c4d30214a2c27561740b7b3f7e1e3b0bbfe4 (diff)
downloadc++-c48fd16fe1554c88c79a1f0d50e81c803da8f61f.tar.gz
c++-c48fd16fe1554c88c79a1f0d50e81c803da8f61f.tar.bz2
c++-c48fd16fe1554c88c79a1f0d50e81c803da8f61f.zip
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
Diffstat (limited to 'lib')
-rw-r--r--lib/cluster.h8
-rw-r--r--lib/cluster_finite.c14
-rw-r--r--lib/cluster_finite.h4
-rw-r--r--lib/finite_states.h33
-rw-r--r--lib/initial_finite.c48
-rw-r--r--lib/wolff.h4
6 files changed, 109 insertions, 2 deletions
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 <R_t, X_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 <R_t, X_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 <cmath>
+#include <functional>
+
+#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 <class X_t>
+void initialize_probs(std::function <double(X_t, X_t)> J, std::function <double(X_t)> 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 <class R_t, class X_t>
void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s0)> gen_R, std::function <void(const state_t <R_t, X_t> *)> 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);