summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/cluster.h65
-rw-r--r--lib/correlation.h8
-rw-r--r--lib/measure.h10
-rw-r--r--lib/wolff.h12
-rw-r--r--lib/z2.h8
5 files changed, 53 insertions, 50 deletions
diff --git a/lib/cluster.h b/lib/cluster.h
index f948586..f3271b0 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -12,45 +12,45 @@
#include "graph.h"
template <class R_t, class X_t>
-void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
+void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, gsl_rng *rand) {
v_t nv = 0;
std::stack<v_t> stack;
stack.push(v0);
- std::vector<bool> marks(state->g.nv, false);
+ std::vector<bool> marks(s.g.nv, false);
while (!stack.empty()) {
v_t v = stack.top();
stack.pop();
- if (!marks[v]) {
+ if (!marks[v]) { // don't reprocess anyone we've already visited!
X_t si_old, si_new;
R_t R_old, R_new;
- R_old = state->R;
+ R_old = s.R;
marks[v] = true;
- bool v_is_ghost = (v == state->nv);
+ bool v_is_ghost = (v == s.nv); // ghost site has the last index
if (v_is_ghost) {
- R_new = r.act(R_old);
+ R_new = r.act(R_old); // if we are, then we're moving the transformation
} else {
- si_old = state->spins[v];
- si_new = r.act(si_old);
+ si_old = s.spins[v];
+ si_new = r.act(si_old); // otherwise, we're moving the spin at our site
}
- for (const v_t &vn : state->g.v_adj[v]) {
+ for (const v_t &vn : s.g.v_adj[v]) {
X_t sj;
- bool vn_is_ghost = (vn == state->nv);
+ bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
if (!vn_is_ghost) {
- sj = state->spins[vn];
+ sj = s.spins[vn];
}
double prob;
- if (v_is_ghost || vn_is_ghost) {
+ if (v_is_ghost || vn_is_ghost) { // if this is a ghost-involved bond...
X_t rs_old, rs_new;
v_t non_ghost;
if (vn_is_ghost) {
@@ -63,51 +63,54 @@ 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);
+ double dE = s.H(rs_old) - s.H(rs_new);
+
#ifdef FINITE_STATES
prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
#else
- prob = 1.0 - exp(-dE / state->T);
+ prob = 1.0 - exp(-dE / s.T);
#endif
- state->M -= rs_old;
- state->M += rs_new;
+ s.M -= rs_old;
+ s.M += rs_new;
- state->E += dE;
+ s.E += dE;
- for (D_t i = 0; i < state->D; i++) {
- L_t x = (non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L;
+ 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;
- state->ReF[i] -= rs_old * state->precomputed_cos[x];
- state->ReF[i] += rs_new * state->precomputed_cos[x];
+ s.ReF[i] -= rs_old * s.precomputed_cos[x];
+ s.ReF[i] += rs_new * s.precomputed_cos[x];
- state->ImF[i] -= rs_old * state->precomputed_sin[x];
- state->ImF[i] += rs_new * state->precomputed_sin[x];
+ s.ImF[i] -= rs_old * s.precomputed_sin[x];
+ s.ImF[i] += rs_new * s.precomputed_sin[x];
}
- } else {
- double dE = state->J(si_old, sj) - state->J(si_new, sj);
+ } else { // otherwise, we're at a perfectly normal bond!
+ double dE = s.J(si_old, sj) - s.J(si_new, sj);
+
#ifdef FINITE_STATES
prob = J_probs[state_to_ind(si_old)][state_to_ind(si_new)][state_to_ind(sj)];
#else
- prob = 1.0 - exp(-dE / state->T);
+ prob = 1.0 - exp(-dE / s.T);
#endif
- state->E += dE;
+
+ s.E += dE;
}
- if (gsl_rng_uniform(rand) < prob) { // and with probability...
+ if (gsl_rng_uniform(rand) < prob) {
stack.push(vn); // push the neighboring vertex to the stack
}
}
if (v_is_ghost) {
- state->R = R_new;
+ s.R = R_new;
} else {
- state->spins[v] = si_new;
+ s.spins[v] = si_new;
nv++;
}
}
}
- state->last_cluster_size = nv;
+ s.last_cluster_size = nv;
}
diff --git a/lib/correlation.h b/lib/correlation.h
index 49c6ff2..23b1bd9 100644
--- a/lib/correlation.h
+++ b/lib/correlation.h
@@ -7,13 +7,13 @@
#include <fftw3.h>
template <class R_t, class X_t>
-double correlation_length(const state_t <R_t, X_t> *s) {
+double correlation_length(const state_t <R_t, X_t>& s) {
double total = 0;
- for (D_t j = 0; j < s->D; j++) {
- total += norm_squared(s->ReF[j]) + norm_squared(s->ReF[j]);
+ for (D_t j = 0; j < s.D; j++) {
+ total += norm_squared(s.ReF[j]) + norm_squared(s.ReF[j]);
}
- return total / s->D;
+ return total / s.D;
}
diff --git a/lib/measure.h b/lib/measure.h
index 8082dd2..2c5ffb7 100644
--- a/lib/measure.h
+++ b/lib/measure.h
@@ -29,17 +29,17 @@ FILE **measure_setup_files(unsigned char flags, unsigned long timestamp) {
}
template <class R_t, class X_t>
-std::function <void(const state_t <R_t, X_t> *)> measure_function_write_files(unsigned char flags, FILE **files, std::function <void(const state_t <R_t, X_t> *)> other_f) {
- return [=] (const state_t <R_t, X_t> *s) {
+std::function <void(const state_t <R_t, X_t>&)> measure_function_write_files(unsigned char flags, FILE **files, std::function <void(const state_t <R_t, X_t>&)> other_f) {
+ return [=] (const state_t <R_t, X_t>& s) {
if (flags & measurement_energy) {
- float smaller_E = (float)s->E;
+ float smaller_E = (float)s.E;
fwrite(&smaller_E, sizeof(float), 1, files[0]);
}
if (flags & measurement_clusterSize) {
- fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, files[1]);
+ fwrite(&(s.last_cluster_size), sizeof(uint32_t), 1, files[1]);
}
if (flags & measurement_magnetization) {
- write_magnetization(s->M, files[2]);
+ write_magnetization(s.M, files[2]);
}
if (flags & measurement_fourierZero) {
float smaller_X = (float)correlation_length(s);
diff --git a/lib/wolff.h b/lib/wolff.h
index a4a663c..498f7f3 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -3,18 +3,18 @@
#include "state.h"
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) {
+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);
+ 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);
+ if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, s.last_cluster_size);
- v_t v0 = gsl_rng_uniform_int(r, s->nv);
- R_t step = gen_R(r, s->spins[v0]);
+ v_t v0 = gsl_rng_uniform_int(r, s.nv);
+ R_t step = gen_R(r, s.spins[v0]);
flip_cluster <R_t, X_t> (s, v0, step, r);
measurements(s);
@@ -23,7 +23,7 @@ void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s
if (!silent) {
printf("\033[F\033[J");
}
- printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s->E, s->last_cluster_size);
+ printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size);
}
diff --git a/lib/z2.h b/lib/z2.h
index fd4722e..a18d740 100644
--- a/lib/z2.h
+++ b/lib/z2.h
@@ -25,7 +25,7 @@ class z2_t {
z2_t(bool x) : x(x) {}
- ising_t act(const ising_t& s) {
+ ising_t act(const ising_t& s) const {
if (x) {
return ising_t(!s.x);
} else {
@@ -33,7 +33,7 @@ class z2_t {
}
}
- z2_t act(const z2_t& r) {
+ z2_t act(const z2_t& r) const {
if (x) {
return z2_t(!r.x);
} else {
@@ -41,11 +41,11 @@ class z2_t {
}
}
- ising_t act_inverse(const ising_t& s) {
+ ising_t act_inverse(const ising_t& s) const {
return this->act(s);
}
- z2_t act_inverse(const z2_t& r) {
+ z2_t act_inverse(const z2_t& r) const {
return this->act(r);
}
};