From 5ffaf0a1bb0f0b47d57d0f24ee1134659775dacb Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 20 Jul 2018 22:57:39 -0400 Subject: added ising example to cpp collection --- lib/cluster.h | 14 +++--- lib/ising.h | 80 ++++++++++++++++++++++++++++++ lib/measure.h | 11 ++++- lib/measurement.c | 145 +++--------------------------------------------------- lib/measurement.h | 25 +++------- lib/state.h | 13 ++--- lib/vector.h | 98 ++++++++++++++++++------------------ lib/wolff.h | 21 +++----- lib/z2.h | 64 ++++++++++++++++++++++++ 9 files changed, 243 insertions(+), 228 deletions(-) create mode 100644 lib/ising.h create mode 100644 lib/z2.h (limited to 'lib') diff --git a/lib/cluster.h b/lib/cluster.h index 2225e2f..4b3c2f6 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -23,6 +23,8 @@ #include "measurement.h" #include "vector.h" #include "orthogonal.h" +#include "ising.h" +#include "z2.h" #include "dihedral.h" #include "dihinf.h" #include "yule_walker.h" @@ -107,17 +109,17 @@ void flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { double dE = state->H(rs_old) - state->H(rs_new); prob = 1.0 - exp(-dE / state->T); - subtract (&(state->M), rs_old); - add (&(state->M), rs_new); + add(&(state->M), -1, rs_old); + add(&(state->M), 1, rs_new); state->E += dE; for (D_t i = 0; i < state->D; i++) { double x = (double)((non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L) / (double)state->L; - scalar_subtract(&(state->ReF[i]), cos(2 * M_PI * x), rs_old); - scalar_add(&(state->ReF[i]), cos(2 * M_PI * x), rs_new); + add(&(state->ReF[i]), -cos(2 * M_PI * x), rs_old); + add(&(state->ReF[i]), cos(2 * M_PI * x), rs_new); - scalar_subtract(&(state->ImF[i]), sin(2 * M_PI * x), rs_old); - scalar_add(&(state->ImF[i]), sin(2 * M_PI * x), rs_new); + add(&(state->ImF[i]), -sin(2 * M_PI * x), rs_old); + add(&(state->ImF[i]), sin(2 * M_PI * x), rs_new); } free_spin (rs_old); diff --git a/lib/ising.h b/lib/ising.h new file mode 100644 index 0000000..4e5164b --- /dev/null +++ b/lib/ising.h @@ -0,0 +1,80 @@ +#pragma once + +#include +#include + +#include "types.h" + +class ising_t { + public: + bool x; + + typedef int M_t; + typedef double F_t; +}; + +void init(ising_t *p) { + p->x = false; +} + +void free_spin(ising_t s) { + // do nothing! +} + +void free_spin(int s) { + // do nothing +} + +void free_spin(double s) { + // do nothing +} + +ising_t copy(ising_t s) { + return s; +} + +template +void add(T *s1, T a, ising_t s2) { + if (s2.x) { + *s1 -= a; + } else { + *s1 += a; + } +} + +int scalar_multiple(int factor, ising_t s) { + if (s.x) { + return -factor; + } else { + return factor; + } +} + + +double norm_squared(double s) { + return pow(s, 2); +} + +template +void write_magnetization(T M, FILE *outfile) { + fwrite(&M, sizeof(T), 1, outfile); +} + +// below this line is unnecessary, but convenient + +double ising_dot(ising_t s1, ising_t s2) { + if (s1.x == s2.x) { + return 1.0; + } else { + return -1.0; + } +} + +double scalar_field(ising_t s, double H) { + if (s.x) { + return -H; + } else { + return H; + } +} + diff --git a/lib/measure.h b/lib/measure.h index d20c081..52e43af 100644 --- a/lib/measure.h +++ b/lib/measure.h @@ -1,7 +1,9 @@ #pragma once -#define POSSIBLE_MEASUREMENTS 4 +#include "measurement.h" + +#define POSSIBLE_MEASUREMENTS 5 const unsigned char measurement_energy = 1 << 0; const unsigned char measurement_clusterSize = 1 << 1; const unsigned char measurement_magnetization = 1 << 2; @@ -43,6 +45,13 @@ std::function *)> measurement_fourier_file(FILE * }; } +template +std::function *)> measurement_average_cluster(meas_t *x) { + return [=](const state_t *s) { + meas_update(x, s->last_cluster_size); + }; +} + #endif diff --git a/lib/measurement.c b/lib/measurement.c index b30cf6b..ad824f6 100644 --- a/lib/measurement.c +++ b/lib/measurement.c @@ -1,15 +1,6 @@ -#include "convex.h" #include "measurement.h" -meas_t *meas_initialize(count_t W) { - meas_t *m = (meas_t *)calloc(1, sizeof(meas_t)); - m->W = W; - m->xx = (double *)calloc(2 * W + 1, sizeof(double)); - - return m; -} - double add_to_avg(double mx, double x, count_t n) { return mx * (n / (n + 1.0)) + x / (n + 1.0); } @@ -19,42 +10,24 @@ void meas_update(meas_t *m, double x) { m->x = add_to_avg(m->x, x, n); m->x2 = add_to_avg(m->x2, pow(x, 2), n); - m->x4 = add_to_avg(m->x4, pow(x, 4), n); m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n); m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n); - dll_t *tmp_window = m->x_window; - dll_t *pos_save; - count_t t = 0; - - while (tmp_window != NULL) { - m->xx[t] = add_to_avg(m->xx[t], x * (tmp_window->x), m->n - t - 1); - t++; - if (t == 2 * (m->W)) { - pos_save = tmp_window; - } - tmp_window = tmp_window->next; + /* + if (n > 1) { + double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2)); + m->dx = sqrt(s2 / n); + m->c = s2; + m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n); } - - if (t == 2 * (m->W) + 1) { - if (2 * (m->W) + 1 == 1) { - free(m->x_window); - m->x_window = NULL; - } else { - free(pos_save->next); - pos_save->next = NULL; - } - } - - stack_push_d(&(m->x_window), x); + */ (m->n)++; } double meas_dx(const meas_t *m) { - return 2 * get_tau(m) * Cxx(m, 0) / m->n; -// return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2))); + return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2))); } double meas_c(const meas_t *m) { @@ -101,105 +74,3 @@ double rho(const autocorr_t *o, count_t i) { return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2)); } -double Cxx(const meas_t *m, count_t t) { - return m->xx[t] - pow(m->x, 2); -} - -double rho_m(const meas_t *m, count_t t) { - return Cxx(m, t) / Cxx(m, 0); -} - -double get_tau(const meas_t *m) { - double *Gammas = (double *)malloc((m->W + 1) * sizeof(double)); - - Gammas[0] = 1 + rho_m(m, 0); - for (uint64_t i = 0; i < m->W; i++) { - Gammas[1 + i] = rho_m(m, 2 * i + 1) + rho_m(m, 2 * i + 2); - } - - uint64_t n; - for (n = 0; n < m->W + 1; n++) { - if (Gammas[n] <= 0) { - break; - } - } - - double *conv_Gamma = get_convex_minorant(n, Gammas); - - double tau = - 0.5; - - for (uint64_t i = 0; i < n + 1; i++) { - tau += conv_Gamma[i]; - } - - free(Gammas); - - return tau; -} - -void print_meas(const meas_t *m, const char *sym, FILE *outfile) { - fprintf(outfile, "%s-><|n->%" PRIcount ",x->%.15f,x^2->%.15f,x^4->%.15f,xx->{", sym, m->n, m->x, m->x2, m->x4); - for (count_t i = 0; i < 2 * (m->W) + 1; i++) { - fprintf(outfile, "%.15f", m->xx[i]); - if (i < 2 * (m->W)) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "}|>"); -} - -void print_vec_meas(q_t q, const meas_t **m, const char *sym, FILE *outfile) { - fprintf(outfile, "%s-><|n->{", sym); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%" PRIcount, m[i]->n); - if (i < q - 1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},x->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", m[i]->x); - if (i < q - 1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},x^2->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", m[i]->x2); - if (i < q - 1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},x^4->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", m[i]->x4); - if (i < q - 1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},xx->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "{"); - for (count_t j = 0; j < 2 * (m[i]->W) + 1; j++) { - fprintf(outfile, "%.15f", m[i]->xx[j]); - if (j < 2 * (m[i]->W)) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "}"); - if (i < q - 1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "}|>"); -} - -void free_meas(meas_t *m) { - free(m->xx); - while (m->x_window != NULL) { - stack_pop_d(&(m->x_window)); - } - free(m); -} - - diff --git a/lib/measurement.h b/lib/measurement.h index d9bd52e..78fa51b 100644 --- a/lib/measurement.h +++ b/lib/measurement.h @@ -3,21 +3,20 @@ #include #include -#include #include "types.h" #include "stack.h" +#ifdef __cplusplus +extern "C" { +#endif + typedef struct { - count_t n; + uint64_t n; double x; double x2; - double x4; double m2; double m4; - count_t W; - double *xx; - dll_t *x_window; } meas_t; typedef struct { @@ -41,14 +40,6 @@ void update_autocorr(autocorr_t *OO, double O); double rho(const autocorr_t *o, uint64_t i); -void print_meas(const meas_t *m, const char *sym, FILE *outfile); -void print_vec_meas(q_t q, const meas_t **m, const char *sym, FILE *outfile); - -void free_meas(meas_t *m); - -meas_t *meas_initialize(count_t W); - -double get_tau(const meas_t *m); - -double Cxx(const meas_t *m, count_t t); - +#ifdef __cplusplus +} +#endif diff --git a/lib/state.h b/lib/state.h index c70459f..b40b85c 100644 --- a/lib/state.h +++ b/lib/state.h @@ -5,6 +5,7 @@ #include "types.h" #include "graph.h" +#include "ising.h" template class state_t { @@ -18,10 +19,10 @@ class state_t { X_t *spins; R_t R; double E; - X_t M; // the "sum" of the spins, like the total magnetization + typename X_t::M_t M; // the "sum" of the spins, like the total magnetization v_t last_cluster_size; - X_t *ReF; - X_t *ImF; + typename X_t::F_t *ReF; + typename X_t::F_t *ImF; std::function J; std::function H; @@ -38,10 +39,10 @@ class state_t { } init (&R); E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]); - M = scalar_multiple(nv, spins[0]); + M = scalar_multiple((int)nv, spins[0]); last_cluster_size = 0; - ReF = (X_t *)malloc(D * sizeof(X_t)); - ImF = (X_t *)malloc(D * sizeof(X_t)); + ReF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t)); + ImF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t)); for (D_t i = 0; i < D; i++) { ReF[i] = scalar_multiple(0, spins[0]); ImF[i] = scalar_multiple(0, spins[0]); diff --git a/lib/vector.h b/lib/vector.h index a1084e2..c478618 100644 --- a/lib/vector.h +++ b/lib/vector.h @@ -6,16 +6,48 @@ #include "types.h" +/* The following is the minimum definition of a spin class. + * + * The class must contain an M_t and an F_t for holding the sum of an + * integer number of spins and a double-weighted number of spins, + * respectively. + * + * void init(X_t *p); + * void free_spin(X_t p); + * X_t copy(X_t x); + * void add(M_t *x1, int factor, X_t x2); + * void add(F_t *x1, double factor, X_t x2); + * M_t scalar_multiple(int factor, X_t x); + * double norm_squared(F_t x); + * void write_magnetization(M_t M, FILE *outfile); + * + */ + template -struct vector_t { T *x; }; +class vector_t { + public: + T *x; + + // M_t needs to hold the sum of nv spins + typedef vector_t M_t; + + // F_t needs to hold the double-weighted sum of spins + typedef vector_t F_t; +}; template void init(vector_t *ptr) { ptr->x = (T *)calloc(q, sizeof(T)); + // initialize a unit vector ptr->x[0] = (T)1; } +template +void free_spin (vector_t v) { + free(v.x); +} + template vector_t copy (vector_t v) { vector_t v_copy; @@ -29,46 +61,15 @@ vector_t copy (vector_t v) { return v_copy; } -template -void add (vector_t *v1, vector_t v2) { - for (q_t i = 0; i < q; i++) { - v1->x[i] += v2.x[i]; - } -} - -template -void scalar_add (vector_t *v1, T a, vector_t v2) { - for (q_t i = 0; i < q; i++) { - v1->x[i] += a * v2.x[i]; - } -} - -template -void scalar_subtract (vector_t *v1, T a, vector_t v2) { - for (q_t i = 0; i < q; i++) { - v1->x[i] -= a * v2.x[i]; - } -} - -template -void subtract (vector_t *v1, vector_t v2) { - for (q_t i = 0; i < q; i++) { - v1->x[i] -= v2.x[i]; - } -} - -template -T norm_squared (vector_t v) { - T tmp = 0; +template +void add(vector_t *v1, V a, vector_t v2) { for (q_t i = 0; i < q; i++) { - tmp += v.x[i] * v.x[i]; + v1->x[i] += (U)(a * v2.x[i]); } - - return tmp; } template -vector_t scalar_multiple(v_t a, vector_t v) { +vector_t scalar_multiple(int a, vector_t v) { vector_t multiple; multiple.x = (T *)malloc(q * sizeof(T)); for (q_t i = 0; i < q; i++) { @@ -79,19 +80,14 @@ vector_t scalar_multiple(v_t a, vector_t v) { } template -T dot(vector_t v1, vector_t v2) { - T prod = 0; +double norm_squared (vector_t v) { + double tmp = 0; for (q_t i = 0; i < q; i++) { - prod += v1.x[i] * v2.x[i]; + tmp += pow(v.x[i], 2); } - return prod; -} - -template -void free_spin (vector_t v) { - free(v.x); + return tmp; } template @@ -99,6 +95,8 @@ void write_magnetization(vector_t M, FILE *outfile) { fwrite(M.x, sizeof(T), q, outfile); } +// below functions and definitions are unnecessary for wolff.h but useful. + template // save some space and don't write whole doubles void write_magnetization(vector_t M, FILE *outfile) { for (q_t i = 0; i < q; i++) { @@ -108,8 +106,14 @@ void write_magnetization(vector_t M, FILE *outfile) { } template -double correlation_component(vector_t v) { - return (double)v.x[0]; +T dot(vector_t v1, vector_t v2) { + T prod = 0; + + for (q_t i = 0; i < q; i++) { + prod += v1.x[i] * v2.x[i]; + } + + return prod; } template diff --git a/lib/wolff.h b/lib/wolff.h index eb357c7..21e88d9 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -3,33 +3,26 @@ #include "state.h" template -void wolff(count_t N, D_t D, L_t L, double T, std::function J, std::function H, std::function *)> gen_R, unsigned int n_measurements, std::function *)> *measurements, bool silent) { - - state_t s(D, L, T, J, H); - - // initialize random number generator - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); +void wolff(count_t N, state_t *s, std::function *)> gen_R, unsigned int n_measurements, std::function *)> *measurements, gsl_rng *r, bool silent) { if (!silent) printf("\n"); for (count_t steps = 0; steps < N; steps++) { - if (!silent) printf("\033[F\033[JWOLFF: sweep %" 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); - flip_cluster (&s, v0, step, r); + v_t v0 = gsl_rng_uniform_int(r, s->nv); + R_t step = gen_R(r, s); + flip_cluster (s, v0, step, r); free_spin(step); for (unsigned int i = 0; i < n_measurements; i++) { - measurements[i](&s); + measurements[i](s); } } if (!silent) { printf("\033[F\033[J"); } - printf("WOLFF: sweep %" 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); - gsl_rng_free(r); } diff --git a/lib/z2.h b/lib/z2.h new file mode 100644 index 0000000..599a6a5 --- /dev/null +++ b/lib/z2.h @@ -0,0 +1,64 @@ + +#pragma once + +#include +#include "types.h" +#include "ising.h" +#include "state.h" + +struct z2_t { bool x; }; + +void init(z2_t *p) { + p->x = false; +} + +void free_spin(z2_t p) { + // do nothing! +} + +z2_t copy(z2_t x) { + return x; +} + +ising_t act(z2_t r, ising_t s) { + ising_t rs; + + if (r.x) { + rs.x = !s.x; + return rs; + } else { + rs.x = s.x; + return rs; + } +} + +z2_t act(z2_t r1, z2_t r2) { + z2_t r3; + + if (r1.x) { + r3.x = !r2.x; + return r3; + } else { + r3.x = r2.x; + return r3; + } +} + +ising_t act_inverse(z2_t r, ising_t s) { + return act(r, s); +} + +z2_t act_inverse(z2_t r1, z2_t r2) { + return act(r1, r2); +} + +// these are all functions necessary for wolff.h + +z2_t generate_ising_rotation(gsl_rng *r, const state_t *s) { + z2_t rot; + rot.x = true; + return rot; +} + + + -- cgit v1.2.3-54-g00ecf