From 3eb67e3bca774eb0441db60158e1968ad901273b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 28 Jun 2018 14:20:25 -0400 Subject: cleaned up the operation of the finite-group wolff code --- lib/cluster_finite.h | 4 + lib/graph.h | 2 + lib/initial_finite.c | 211 +++++++++++++++++++++++++++++++++++++++++++++++++++ lib/initial_finite.h | 24 ++++++ lib/measurement.h | 5 -- lib/wolff_finite.c | 70 ----------------- 6 files changed, 241 insertions(+), 75 deletions(-) create mode 100644 lib/initial_finite.c create mode 100644 lib/initial_finite.h delete mode 100644 lib/wolff_finite.c (limited to 'lib') diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h index abdc8fc..ad45ed3 100644 --- a/lib/cluster_finite.h +++ b/lib/cluster_finite.h @@ -23,6 +23,10 @@ #include "yule_walker.h" typedef struct { + D_t D; + L_t L; + v_t nv; + v_t ne; graph_t *g; q_t q; R_t n_transformations; diff --git a/lib/graph.h b/lib/graph.h index 9c80dd6..cb47faa 100644 --- a/lib/graph.h +++ b/lib/graph.h @@ -1,4 +1,6 @@ +#pragma once + #include #include #include diff --git a/lib/initial_finite.c b/lib/initial_finite.c new file mode 100644 index 0000000..f286dcc --- /dev/null +++ b/lib/initial_finite.c @@ -0,0 +1,211 @@ + +#include "initial_finite.h" + +double *Jprobs_from_J(q_t q, double T, double *J) { + double *J_probs = (double *)calloc(pow(q, 2), sizeof(double)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + J_probs[q * i + j] = 1.0 - exp((J[i] - J[j]) / T); + } + } + + return J_probs; +} + +q_t *initialize_R(q_t q) { + q_t *R = (q_t *)malloc(q * sizeof(q_t)); + + for (q_t i = 0; i < q; i++) { + R[i] = i; + } + + return R; +} + +state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H) { + state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t)); + + s->D = D; + s->L = L; + + { + graph_t *g = graph_create_square(D, L); + s->nv = g->nv; + s->ne = g->ne; + s->g = graph_add_ext(g); + graph_free(g); + } + + s->q = 2; + s->n_transformations = 1; + + s->transformations = (q_t *)malloc(2 * sizeof(q_t)); + s->transformations[0] = 1; + s->transformations[1] = 0; + + s->T = T; + s->J = (double *)malloc(2 * sizeof(double)); + s->J[0] = 1.0; + s->J[1] = -1.0; + s->H = (double *)malloc(2 * sizeof(double)); + s->H[0] = H[0]; + s->H[1] = -H[0]; + + s->J_probs = Jprobs_from_J(2, T, s->J); + s->H_probs = Jprobs_from_J(2, T, s->H); + + s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); + s->R = initialize_R(2); + + s->E = - ((double)s->ne) * s->J[0] - ((double)s->nv) * s->H[0]; + s->M = (v_t *)calloc(2, sizeof(v_t)); + s->M[0] = s->nv; // everyone starts in state 0, remember? + + return s; +} + +state_finite_t *initial_finite_prepare_potts(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)); + + s->D = D; + s->L = L; + + { + graph_t *g = graph_create_square(D, L); + s->nv = g->nv; + s->ne = g->ne; + s->g = graph_add_ext(g); + graph_free(g); + } + + s->q = q; + s->n_transformations = q; + s->transformations = dihedral_gen_transformations(q); + + s->T = T; + s->J = (double *)calloc(q, sizeof(double)); + s->J[0] = 1.0; + + s->H = (double *)malloc(q * sizeof(double)); + for (q_t i = 0; i < q; i++) { + s->H[i] = H[i]; + } + + s->J_probs = Jprobs_from_J(q, T, s->J); + s->H_probs = Jprobs_from_J(q, T, s->H); + + s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); + s->R = initialize_R(q); + + s->E = - ((double)s->ne) * s->J[0] - ((double)s->nv) * s->H[0]; + s->M = (v_t *)calloc(q, sizeof(v_t)); + s->M[0] = s->nv; // everyone starts in state 0, remember? + + return s; +} + +state_finite_t *initial_finite_prepare_clock(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)); + + s->D = D; + s->L = L; + + { + graph_t *g = graph_create_square(D, L); + s->nv = g->nv; + s->ne = g->ne; + s->g = graph_add_ext(g); + graph_free(g); + } + + s->q = q; + s->n_transformations = q; + s->transformations = dihedral_gen_transformations(q); + + s->T = T; + s->J = (double *)malloc(q * sizeof(double)); + + for (q_t i = 0; i < q; i++) { + s->J[i] = cos(2 * M_PI * i / ((double)q)); + } + + + s->H = (double *)malloc(q * sizeof(double)); + for (q_t i = 0; i < q; i++) { + s->H[i] = H[i]; + } + + s->J_probs = Jprobs_from_J(q, T, s->J); + s->H_probs = Jprobs_from_J(q, T, s->H); + + s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); + s->R = initialize_R(q); + + s->E = - ((double)s->ne) * s->J[0] - ((double)s->nv) * s->H[0]; + s->M = (v_t *)calloc(q, sizeof(v_t)); + s->M[0] = s->nv; // everyone starts in state 0, remember? + + 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)); + + s->D = D; + s->L = L; + + { + graph_t *g = graph_create_square(D, L); + s->nv = g->nv; + s->ne = g->ne; + s->g = graph_add_ext(g); + graph_free(g); + } + + s->q = q; + s->n_transformations = q; + s->transformations = dihedral_gen_transformations(q); + + s->T = T; + s->J = (double *)malloc(q * sizeof(double)); + + for (q_t i = 0; i < q / 2 + 1; i++) { + s->J[i] = -pow(i, 2); + } + for (q_t i = 1; i < (q + 1) / 2; i++) { + s->J[q - i] = -pow(i, 2); + } + + s->H = (double *)malloc(q * sizeof(double)); + for (q_t i = 0; i < q; i++) { + s->H[i] = H[i]; + } + + s->J_probs = Jprobs_from_J(q, T, s->J); + s->H_probs = Jprobs_from_J(q, T, s->H); + + s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); + s->R = initialize_R(q); + + s->E = - ((double)s->ne) * s->J[0] - ((double)s->nv) * s->H[0]; + s->M = (v_t *)calloc(q, sizeof(v_t)); + s->M[0] = s->nv; // everyone starts in state 0, remember? + + return s; +} + +void state_finite_free(state_finite_t *s) { + graph_free(s->g); + free(s->J); + free(s->H); + free(s->J_probs); + free(s->H_probs); + free(s->spins); + free(s->R); + free(s->M); + free(s->transformations); + free(s); +} + diff --git a/lib/initial_finite.h b/lib/initial_finite.h new file mode 100644 index 0000000..65414cd --- /dev/null +++ b/lib/initial_finite.h @@ -0,0 +1,24 @@ + +#pragma once + +#include + +#include "types.h" +#include "dihedral.h" +#include "cluster_finite.h" + +typedef enum { + ISING, + POTTS, + CLOCK, + DGM +} finite_model_t; + +state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H); +state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, double *H); +state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, double *H); +state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double *H); + +void state_finite_free(state_finite_t *s); + + diff --git a/lib/measurement.h b/lib/measurement.h index 46c034f..eaa260b 100644 --- a/lib/measurement.h +++ b/lib/measurement.h @@ -24,11 +24,6 @@ typedef struct { double O2; } autocorr_t; -typedef struct { - void (*f)(state_finite_t *, void *); - void *data; -} measurement_t; - void meas_update(meas_t *m, double x); double meas_dx(const meas_t *m); diff --git a/lib/wolff_finite.c b/lib/wolff_finite.c deleted file mode 100644 index 64de9ba..0000000 --- a/lib/wolff_finite.c +++ /dev/null @@ -1,70 +0,0 @@ - -#include "cluster_finite.h" - -void wolff_finite(state_finite_t *s, count_t sweeps, count_t sweeps_per_measurement, count_t n_measurements, measurement_t *measurements) { - for (count_t i = 0; i < sweeps; i++) { - - count_t n_flips = 0; - - while (n_flips / h->nv < sweeps_per_measurement) { - v_t v0 = gsl_rng_uniform_int(r, h->nv); - R_t step; - - bool changed = false; - while (!changed) { - step = gsl_rng_uniform_int(r, s->n_transformations); - if v(symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) { - changed = true; - } - } - - v_t tmp_flips = flip_cluster_finite(s, v0, step, r); - n_flips += tmp_flips; - - if (n_runs > 0) { - n_steps++; - meas_update(clust, tmp_flips); - - if (record_autocorrelation && n_steps % ac_skip == 0) { - update_autocorr(autocorr, s->E); - } - - } - - } - - for (q_t i = 0; i < q; i++) { - meas_update(M[i], s->M[i]); - } - meas_update(E, s->E); - - q_t n_at_max = 0; - q_t max_M_i = 0; - v_t max_M = 0; - - for (q_t i = 0; i < q; i++) { - if (s->M[i] > max_M) { - n_at_max = 1; - max_M_i = i; - max_M = s->M[i]; - } else if (s->M[i] == max_M) { - n_at_max++; - } - } - - if (record_distribution) { - mag_dist[s->M[0]]++; - } - - if (n_at_max == 1) { - for (q_t i = 0; i < q; i++) { - meas_update(sM[max_M_i][i], s->M[i]); - } - meas_update(sE[max_M_i], s->E); - freqs[max_M_i]++; - } - - diff = fabs(meas_dx(clust) / clust->x); - } -} - -- cgit v1.2.3-54-g00ecf From 78d8de381f0b1e99ad98364709cbf876689628b2 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 29 Jun 2018 15:43:10 -0400 Subject: completely removed measurement during the simulation, opting to just save binary data points to files throughout --- lib/cluster_finite.c | 4 +- lib/cluster_finite.h | 2 +- lib/initial_finite.c | 22 ++- lib/initial_finite.h | 3 + lib/measurement.c | 145 +++++++++++++++++-- lib/measurement.h | 18 ++- src/wolff_finite.c | 398 +++++++++------------------------------------------ 7 files changed, 242 insertions(+), 350 deletions(-) (limited to 'lib') diff --git a/lib/cluster_finite.c b/lib/cluster_finite.c index f11a3ea..71396e0 100644 --- a/lib/cluster_finite.c +++ b/lib/cluster_finite.c @@ -62,14 +62,14 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) { s->M[rot_s_old]--; s->M[rot_s_new]++; - s->E += - s->H[rot_s_new] + s->H[rot_s_old]; } else { q_t diff_old = (s_old + s->q - sn) % s->q; q_t diff_new = (s_new + s->q - sn) % s->q; prob = s->J_probs[diff_new * s->q + diff_old]; - s->E += - s->J[diff_new] + s->J[diff_old]; + s->B[diff_old]--; + s->B[diff_new]++; } if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h index ad45ed3..701c197 100644 --- a/lib/cluster_finite.h +++ b/lib/cluster_finite.h @@ -38,7 +38,7 @@ typedef struct { double *H_probs; q_t *spins; q_t *R; - double E; + v_t *B; v_t *M; } state_finite_t; diff --git a/lib/initial_finite.c b/lib/initial_finite.c index f286dcc..fb120f0 100644 --- a/lib/initial_finite.c +++ b/lib/initial_finite.c @@ -58,9 +58,10 @@ state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H) s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); s->R = initialize_R(2); - s->E = - ((double)s->ne) * s->J[0] - ((double)s->nv) * s->H[0]; s->M = (v_t *)calloc(2, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? + s->B = (v_t *)calloc(2, sizeof(v_t)); + s->B[0] = s->ne; return s; } @@ -98,9 +99,10 @@ state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, doub s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); s->R = initialize_R(q); - s->E = - ((double)s->ne) * s->J[0] - ((double)s->nv) * s->H[0]; s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? + s->B = (v_t *)calloc(q, sizeof(v_t)); + s->B[0] = s->ne; // everyone starts in state 0, remember? return s; } @@ -142,9 +144,10 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); s->R = initialize_R(q); - s->E = - ((double)s->ne) * s->J[0] - ((double)s->nv) * s->H[0]; s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? + s->B = (v_t *)calloc(q, sizeof(v_t)); + s->B[0] = s->ne; // everyone starts in state 0, remember? return s; } @@ -189,13 +192,23 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); s->R = initialize_R(q); - s->E = - ((double)s->ne) * s->J[0] - ((double)s->nv) * s->H[0]; s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? return s; } +double state_finite_energy(state_finite_t *s) { + double E = 0; + + for (q_t i = 0; i < s->q; i++) { + E += s->J[i] * s->B[i]; + E += s->H[i] * s->M[i]; + } + + return -E; +} + void state_finite_free(state_finite_t *s) { graph_free(s->g); free(s->J); @@ -205,6 +218,7 @@ void state_finite_free(state_finite_t *s) { free(s->spins); free(s->R); free(s->M); + free(s->B); free(s->transformations); free(s); } diff --git a/lib/initial_finite.h b/lib/initial_finite.h index 65414cd..542f923 100644 --- a/lib/initial_finite.h +++ b/lib/initial_finite.h @@ -7,6 +7,8 @@ #include "dihedral.h" #include "cluster_finite.h" +static char *finite_model_t_strings[] = {"ISING", "POTTS", "CLOCK", "DGM"}; + typedef enum { ISING, POTTS, @@ -21,4 +23,5 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double void state_finite_free(state_finite_t *s); +double state_finite_energy(state_finite_t *s); diff --git a/lib/measurement.c b/lib/measurement.c index ad824f6..b30cf6b 100644 --- a/lib/measurement.c +++ b/lib/measurement.c @@ -1,6 +1,15 @@ +#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); } @@ -10,24 +19,42 @@ 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); - /* - 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); + 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 (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 sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2))); + return 2 * get_tau(m) * Cxx(m, 0) / m->n; +// return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2))); } double meas_c(const meas_t *m) { @@ -74,3 +101,105 @@ 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 eaa260b..d9bd52e 100644 --- a/lib/measurement.h +++ b/lib/measurement.h @@ -3,16 +3,21 @@ #include #include +#include #include "types.h" #include "stack.h" typedef struct { - uint64_t n; + count_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 { @@ -36,3 +41,14 @@ 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); + diff --git a/src/wolff_finite.c b/src/wolff_finite.c index 47fcc88..e41c326 100644 --- a/src/wolff_finite.c +++ b/src/wolff_finite.c @@ -1,93 +1,60 @@ +#include #include #include int main(int argc, char *argv[]) { - L_t L = 128; count_t N = (count_t)1e7; - count_t min_runs = 10; - count_t n = 3; + + finite_model_t model = ISING; + q_t q = 2; D_t D = 2; + L_t L = 128; double T = 2.26918531421; double *J = (double *)calloc(MAX_Q, sizeof(double)); J[0] = 1.0; double *H = (double *)calloc(MAX_Q, sizeof(double)); - double eps = 0; - bool silent = false; - bool snapshots = false; - bool snapshot = false; - bool record_autocorrelation = false; - bool record_distribution = false; - count_t W = 10; - count_t ac_skip = 1; - finite_model_t model = ISING; + bool silent = false; int opt; q_t J_ind = 0; q_t H_ind = 0; - while ((opt = getopt(argc, argv, "N:n:D:L:q:T:J:H:m:e:IpsSPak:W:drt:")) != -1) { + while ((opt = getopt(argc, argv, "N:t:q:D:L:T:J:H:s")) != -1) { switch (opt) { - case 'N': + case 'N': // number of steps N = (count_t)atof(optarg); break; - case 'n': - n = (count_t)atof(optarg); + case 't': // type of simulation + model = (finite_model_t)atoi(optarg); + break; + case 'q': // number of states, if relevant + q = atoi(optarg); break; - case 'D': + case 'D': // dimension D = atoi(optarg); break; - case 'L': + case 'L': // linear size L = atoi(optarg); break; - case 'q': - q = atoi(optarg); - break; - case 'T': + case 'T': // temperature T = atof(optarg); break; - case 'J': + case 'J': // couplings, if relevant. nth call couples states i and i + n J[J_ind] = atof(optarg); J_ind++; break; - case 'H': + case 'H': // external field. nth call couples to state n H[H_ind] = atof(optarg); H_ind++; break; - case 'm': - min_runs = atoi(optarg); - break; - case 'e': - eps = atof(optarg); - break; - case 's': + case 's': // don't print anything during simulation. speeds up slightly silent = true; break; - case 'S': - snapshots = true; - break; - case 'P': - snapshot = true; - break; - case 'a': - record_autocorrelation = true; - break; - case 'k': - ac_skip = (count_t)atof(optarg); - break; - case 'W': - W = (count_t)atof(optarg); - break; - case 'd': - record_distribution = true; - break; - case 't': - model = (finite_model_t)atoi(optarg); - break; default: exit(EXIT_FAILURE); } @@ -95,9 +62,6 @@ int main(int argc, char *argv[]) { state_finite_t *s; - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - switch (model) { case ISING: s = initial_finite_prepare_ising(D, L, T, H); @@ -113,318 +77,84 @@ int main(int argc, char *argv[]) { break; default: printf("Not a valid model!\n"); - return 1; + free(J); + free(H); + exit(EXIT_FAILURE); } free(J); free(H); + // initialize random number generator + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(r, rand_seed()); - double diff = 1e31; - count_t n_runs = 0; - count_t n_steps = 0; + unsigned long timestamp = (unsigned long)time(NULL); - meas_t *E, *clust, **M, **sE, ***sM; + char *filename_M = (char *)malloc(256 * sizeof(char)); + char *filename_B = (char *)malloc(256 * sizeof(char)); + char *filename_S = (char *)malloc(256 * sizeof(char)); - M = (meas_t **)malloc(q * sizeof(meas_t *)); + sprintf(filename_M, "wolff_%s_%" PRIq "_%" PRID "_%" PRIL "_%.8f", finite_model_t_strings[model], q, D, L, T); for (q_t i = 0; i < q; i++) { - M[i] = (meas_t *)calloc(1, sizeof(meas_t)); + sprintf(filename_M + strlen(filename_M), "_%.8f", s->H[i]); } - E = calloc(1, sizeof(meas_t)); - clust = calloc(1, sizeof(meas_t)); + strcpy(filename_B, filename_M); + strcpy(filename_S, filename_M); - sE = (meas_t **)malloc(q * sizeof(meas_t *)); - sM = (meas_t ***)malloc(q * sizeof(meas_t **)); + sprintf(filename_M + strlen(filename_M), "_%lu_M.dat", timestamp); + sprintf(filename_B + strlen(filename_B), "_%lu_B.dat", timestamp); + sprintf(filename_S + strlen(filename_S), "_%lu_S.dat", timestamp); - for (q_t i = 0; i < q; i++) { - sE[i] = (meas_t *)calloc(1, sizeof(meas_t)); - sM[i] = (meas_t **)malloc(q * sizeof(meas_t *)); - for (q_t j = 0; j < q; j++) { - sM[i][j] = (meas_t *)calloc(1, sizeof(meas_t)); - } - } + FILE *outfile_M = fopen(filename_M, "wb"); + FILE *outfile_B = fopen(filename_B, "wb"); + FILE *outfile_S = fopen(filename_S, "wb"); - count_t *freqs = (count_t *)calloc(q, sizeof(count_t)); - q_t cur_M = 0; + free(filename_M); + free(filename_B); + free(filename_S); - autocorr_t *autocorr; - if (record_autocorrelation) { - autocorr = (autocorr_t *)calloc(1, sizeof(autocorr_t)); - autocorr->W = 2 * W + 1; - autocorr->OO = (double *)calloc(2 * W + 1, sizeof(double)); - } - - count_t *mag_dist; - if (record_distribution) { - mag_dist = (count_t *)calloc(s->nv + 1, sizeof(count_t)); - } + v_t cluster_size = 0; if (!silent) printf("\n"); - while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) { - if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 - ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", - n_runs, fabs(meas_dx(E) / E->x), meas_dx(M[0]) / M[0]->x, meas_dc(E) / meas_c(E), meas_dc(M[0]) / meas_c(M[0]), s->nv / clust->x); - - count_t n_flips = 0; + for (count_t steps = 0; steps < N; steps++) { + if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, B_0 = %" PRIv ", M_0 = %" PRIv ", S = %" PRIv "\n", steps, N, state_finite_energy(s), s->B[0], s->M[0], cluster_size); - while (n_flips / s->nv < n) { - v_t v0 = gsl_rng_uniform_int(r, s->nv); - R_t step; + v_t v0 = gsl_rng_uniform_int(r, s->nv); + R_t step; - bool changed = false; - while (!changed) { - step = gsl_rng_uniform_int(r, s->n_transformations); - if (symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) { - changed = true; - } + bool changed = false; + while (!changed) { + step = gsl_rng_uniform_int(r, s->n_transformations); + if (symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) { + changed = true; } - - v_t tmp_flips = flip_cluster_finite(s, v0, step, r); - n_flips += tmp_flips; - - if (n_runs > 0) { - n_steps++; - meas_update(clust, tmp_flips); - - if (record_autocorrelation && n_steps % ac_skip == 0) { - update_autocorr(autocorr, s->E); - } - - } - } - for (q_t i = 0; i < q; i++) { - meas_update(M[i], s->M[i]); - } - meas_update(E, s->E); - - q_t n_at_max = 0; - q_t max_M_i = 0; - v_t max_M = 0; - - for (q_t i = 0; i < q; i++) { - if (s->M[i] > max_M) { - n_at_max = 1; - max_M_i = i; - max_M = s->M[i]; - } else if (s->M[i] == max_M) { - n_at_max++; - } - } + cluster_size = flip_cluster_finite(s, v0, step, r); - if (record_distribution) { - mag_dist[s->M[0]]++; - } + // v_t is never going to be bigger than 32 bits, but since it's specified + // as a fast time many machines will actually have it be 64 bits. we cast + // it down here to halve space. - if (n_at_max == 1) { - for (q_t i = 0; i < q; i++) { - meas_update(sM[max_M_i][i], s->M[i]); - } - meas_update(sE[max_M_i], s->E); - freqs[max_M_i]++; + for (q_t i = 0; i < q - 1; i++) { + fwrite(&(s->M[i]), sizeof(uint32_t), 1, outfile_M); // if we know the occupation of the first q - 1 states, we know the occupation of the last + fwrite(&(s->B[i]), sizeof(uint32_t), 1, outfile_B); // if we know the occupation of the first q - 1 states, we know the occupation of the last } - diff = fabs(meas_dx(clust) / clust->x); + fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S); - n_runs++; } if (!silent) { printf("\033[F\033[J"); } - printf("WOLFF: sweep %" PRIu64 - ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", - n_runs, fabs(meas_dx(E) / E->x), meas_dx(M[0]) / M[0]->x, meas_dc(E) / meas_c(E), meas_dc(M[0]) / meas_c(M[0]), s->nv / clust->x); - - if (snapshots) { - FILE *snapfile = fopen("snapshots.m", "a"); - fprintf(snapfile, "\n"); - } - - if (snapshot) { - q_t *R_inv = symmetric_invert(q, s->R); - FILE *snapfile = fopen("snapshot.m", "a"); - fprintf(snapfile, "{{"); - for (L_t i = 0; i < L; i++) { - fprintf(snapfile, "{"); - for (L_t j = 0; j < L; j++) { - fprintf(snapfile, "%" PRIq, symmetric_act(R_inv, s->spins[L * i + j])); - if (j != L - 1) { - fprintf(snapfile, ","); - } - } - fprintf(snapfile, "}"); - if (i != L - 1) { - fprintf(snapfile, ","); - } - } - fprintf(snapfile, "}}\n"); - fclose(snapfile); - } + printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, B_0 = %" PRIv ", M_0 = %" PRIv ", S = %" PRIv "\n", N, N, state_finite_energy(s), s->B[0], s->M[0], cluster_size); - double tau = 0; - int tau_failed = 0; + fclose(outfile_M); + fclose(outfile_B); + fclose(outfile_S); - if (record_autocorrelation) { - double *Gammas = (double *)malloc((W + 1) * sizeof(double)); - - Gammas[0] = 1 + rho(autocorr, 0); - for (uint64_t i = 0; i < W; i++) { - Gammas[1 + i] = rho(autocorr, 2 * i + 1) + rho(autocorr, 2 * i + 2); - } - - uint64_t n; - for (n = 0; n < W + 1; n++) { - if (Gammas[n] <= 0) { - break; - } - } - - if (n == W + 1) { - printf("WARNING: correlation function never hit the noise floor.\n"); - tau_failed = 1; - } - - if (n < 2) { - printf("WARNING: correlation function only has one nonnegative term.\n"); - tau_failed = 2; - } - - double *conv_Gamma = get_convex_minorant(n, Gammas); - - double ttau = - 0.5; - - for (uint64_t i = 0; i < n + 1; i++) { - ttau += conv_Gamma[i]; - } - - tau = ttau * ac_skip * clust->x / s->nv; - - free(Gammas); - free(autocorr->OO); - while (autocorr->Op != NULL) { - stack_pop_d(&(autocorr->Op)); - } - free(autocorr); - } - - if (tau_failed) { - //tau = 0; - } - - { - FILE *outfile = fopen("out.m", "a"); - - fprintf(outfile, "<|N->%" PRIcount ",n->%" PRIcount ",D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,J->{", N, n, D, L, q, T); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", s->J[i]); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},H->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", s->H[i]); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->{", E->x / s->nv, meas_dx(E) / s->nv, meas_c(E) / s->nv, meas_dc(E) / s->nv); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", M[i]->x / s->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},\\[Delta]M->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", meas_dx(M[i]) / s->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},\\[Chi]->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", meas_c(M[i]) / s->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},\\[Delta]\\[Chi]->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", meas_dc(M[i]) / s->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "},Subscript[E,%" PRIq "]->%.15f,Subscript[\\[Delta]E,%" PRIq "]->%.15f,Subscript[C,%" PRIq "]->%.15f,Subscript[\\[Delta]C,%" PRIq "]->%.15f,Subscript[M,%" PRIq "]->{", i, sE[i]->x / s->nv, i, meas_dx(sE[i]) / s->nv, i, meas_c(sE[i]) / s->nv, i, meas_dc(sE[i]) / s->nv, i); - for (q_t j = 0; j < q; j++) { - fprintf(outfile, "%.15f", sM[i][j]->x / s->nv); - if (j != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},Subscript[\\[Delta]M,%" PRIq "]->{", i); - for (q_t j = 0; j < q; j++) { - fprintf(outfile, "%.15f", meas_dx(sM[i][j]) / s->nv); - if (j != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},Subscript[\\[Chi],%" PRIq "]->{", i); - for (q_t j = 0; j < q; j++) { - fprintf(outfile, "%.15f", meas_c(sM[i][j]) / s->nv); - if (j != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},Subscript[\\[Delta]\\[Chi],%" PRIq "]->{", i); - for (q_t j = 0; j < q; j++) { - fprintf(outfile, "%.15f", meas_dc(sM[i][j]) / s->nv); - if (j != q-1) { - fprintf(outfile, ","); - } - } - } - fprintf(outfile,"}"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, ",Subscript[f,%" PRIq "]->%.15f,Subscript[\\[Delta]f,%" PRIq "]->%.15f", i, (double)freqs[i] / (double)n_runs, i, sqrt(freqs[i]) / (double)n_runs); - } - fprintf(outfile, ",Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f,\\[Tau]s->%d", clust->x / s->nv, meas_dx(clust) / s->nv, meas_c(clust) / s->nv, meas_dc(clust) / s->nv,tau,tau_failed); - if (record_distribution) { - fprintf(outfile, ",S->{"); - for (v_t i = 0; i < s->nv + 1; i++) { - fprintf(outfile, "%" PRIcount, mag_dist[i]); - if (i != s->nv) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "}"); - free(mag_dist); - } - fprintf(outfile, "|>\n"); - - fclose(outfile); - } - - free(E); - free(clust); - for (q_t i = 0; i < q; i++) { - free(M[i]); - for (q_t j = 0; j < q; j++) { - free(sM[i][j]); - } - free(sM[i]); - } - free(M); - free(sM); - for (q_t i = 0; i < q; i++) { - free(sE[i]); - } - free(freqs); - free(sE); state_finite_free(s); gsl_rng_free(r); -- cgit v1.2.3-54-g00ecf From 90ae915b5a7961a36e6a33509b16229244c6615a Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 1 Jul 2018 00:52:31 -0400 Subject: fixed both the system for determining bond energy and how global state is tracked --- lib/cluster.c | 91 -------------------------------- lib/cluster_finite.c | 10 ++-- lib/cluster_finite.h | 7 ++- lib/dihedral.c | 23 ++++++-- lib/dihedral.h | 4 +- lib/initial_finite.c | 145 +++++++++++++++++++++++++++++++++++++++++++-------- lib/symmetric.c | 43 +++++++++++++++ lib/symmetric.h | 2 + src/wolff_finite.c | 12 ++--- 9 files changed, 207 insertions(+), 130 deletions(-) (limited to 'lib') diff --git a/lib/cluster.c b/lib/cluster.c index 7274eb9..96225a2 100644 --- a/lib/cluster.c +++ b/lib/cluster.c @@ -1,97 +1,6 @@ #include "cluster.h" -v_t flip_cluster(ising_state_t *s, v_t v0, q_t rot, gsl_rng *r) { - v_t nv = 0; - - ll_t *stack = NULL; // create a new stack - stack_push(&stack, v0); // push the initial vertex to the stack - - bool *marks = (bool *)calloc(s->g->nv, sizeof(bool)); - - while (stack != NULL) { - v_t v = stack_pop(&stack); - - if (!marks[v]) { - q_t s_old, s_new; - dihedral_t *R_new; - bool external_flipped; - - marks[v] = true; - - if (v == s->g->nv - 1) { - R_new = dihedral_compose(s->q, rot, s->R); - external_flipped = true; - } else { - s_old = s->spins[v]; - s_new = dihedral_act(s->q, rot, s_old); - external_flipped = false; - } - - 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; - double prob; - bool external_neighbor = false; - - v_t vn = s->g->v_adj[s->g->v_i[v] + i]; - - if (vn == s->g->nv - 1) { - external_neighbor = true; - } else { - sn = s->spins[vn]; - } - - if (external_flipped || external_neighbor) { - q_t rot_s_old, rot_s_new; - - if (external_neighbor) { - rot_s_old = dihedral_inverse_act(s->q, s->R, s_old); - rot_s_new = dihedral_inverse_act(s->q, s->R, s_new); - } else { - rot_s_old = dihedral_inverse_act(s->q, s->R, sn); - rot_s_new = dihedral_inverse_act(s->q, R_new, sn); - } - - prob = s->H_probs[rot_s_new * s->q + rot_s_old]; - - s->M[rot_s_old]--; - s->M[rot_s_new]++; - - s->E += - s->H[rot_s_new] + s->H[rot_s_old]; - } else { - q_t diff_old = (s_old + s->q - sn) % s->q; - q_t diff_new = (s_new + s->q - sn) % s->q; - - prob = s->J_probs[diff_new * s->q + diff_old]; - - s->E += - s->J[diff_new] + s->J[diff_old]; - } - - if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... - stack_push(&stack, vn); // push the neighboring vertex to the stack - } - } - - if (external_flipped) { - free(s->R); - s->R = R_new; - } else { - s->spins[v] = s_new; - } - - if (v != s->g->nv - 1) { // count the number of non-external sites that flip - nv++; - } - } - } - - free(marks); - - return nv; -} - v_t flip_cluster_dgm(dgm_state_t *s, v_t v0, h_t rot, gsl_rng *r) { v_t nv = 0; diff --git a/lib/cluster_finite.c b/lib/cluster_finite.c index 71396e0..9392cf8 100644 --- a/lib/cluster_finite.c +++ b/lib/cluster_finite.c @@ -1,8 +1,8 @@ #include "cluster_finite.h" -v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) { - q_t *rot = s->transformations + s->q * rot_ind; +v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) { + q_t *rot = s->transformations + s->q * s->involutions[rot_ind]; q_t *R_inv = symmetric_invert(s->q, s->R); v_t nv = 0; @@ -63,10 +63,10 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) { s->M[rot_s_new]++; } else { - q_t diff_old = (s_old + s->q - sn) % s->q; - q_t diff_new = (s_new + s->q - sn) % s->q; + 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]]; - prob = s->J_probs[diff_new * s->q + diff_old]; + prob = s->J_probs[diff_new * s->n_bond_types + diff_old]; s->B[diff_old]--; s->B[diff_new]++; diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h index 701c197..b2d764e 100644 --- a/lib/cluster_finite.h +++ b/lib/cluster_finite.h @@ -31,6 +31,11 @@ typedef struct { q_t q; R_t n_transformations; q_t *transformations; + R_t n_involutions; + R_t *involutions; + R_t *transform_site_to_zero; + q_t n_bond_types; + q_t *bond_with_zero_type; double T; double *J; double *H; @@ -42,5 +47,5 @@ typedef struct { v_t *M; } state_finite_t; -v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot, gsl_rng *r); +v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot, gsl_rng *r); diff --git a/lib/dihedral.c b/lib/dihedral.c index ac74a23..8158b43 100644 --- a/lib/dihedral.c +++ b/lib/dihedral.c @@ -11,10 +11,14 @@ dihedral_t *dihedral_compose(q_t q, q_t g1i, const dihedral_t *g2) { return g3; } -q_t dihedral_act(q_t q, q_t gi, q_t s) { +q_t dihedral_act(q_t q, q_t gi, bool r, q_t s) { // we only need to consider the action of reflections - return (gi + q - s) % q; + if (r) { + return (gi + q - s) % q; + } else { + return (gi + s) % q; + } } q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) { @@ -26,15 +30,26 @@ q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) { } q_t *dihedral_gen_transformations(q_t q) { - q_t *transformations = (q_t *)malloc(q * q * sizeof(q_t)); + q_t *transformations = (q_t *)malloc(2 * q * q * sizeof(q_t)); for (q_t i = 0; i < q; i++) { for (q_t j = 0; j < q; j++) { - transformations[q * i + j] = dihedral_act(q, i, j); + transformations[q * i + j] = dihedral_act(q, i, false, j); + transformations[q * q + q * i + j] = dihedral_act(q, i, true, j); } } return transformations; } +R_t *dihedral_gen_involutions(q_t q) { + R_t *transformations = (R_t *)malloc(q * sizeof(R_t)); + + for (q_t i = 0; i < q; i++) { + transformations[i] = q + i; + } + + return transformations; +} + diff --git a/lib/dihedral.h b/lib/dihedral.h index e5e4cbd..c95b23a 100644 --- a/lib/dihedral.h +++ b/lib/dihedral.h @@ -11,9 +11,11 @@ typedef struct { dihedral_t *dihedral_compose(q_t q, q_t gti, const dihedral_t *g2); -q_t dihedral_act(q_t q, q_t gi, q_t s); +q_t dihedral_act(q_t q, q_t gi, bool r, q_t s); q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s); q_t *dihedral_gen_transformations(q_t q); +R_t *dihedral_gen_involutions(q_t q); +R_t factorial(q_t); diff --git a/lib/initial_finite.c b/lib/initial_finite.c index fb120f0..6ea76ef 100644 --- a/lib/initial_finite.c +++ b/lib/initial_finite.c @@ -23,6 +23,40 @@ q_t *initialize_R(q_t q) { return R; } +R_t *transformation_bringing_to_zero(q_t q, R_t n_transformations, q_t *transformations) { + R_t *destination = (R_t *)malloc(q * sizeof(R_t)); + + for (q_t i = 0; i < q; i++) { + for (R_t j = 0; j < n_transformations; j++) { + if (transformations[q * j + i] == 0) { + destination[i] = j; + } + } + } + + return destination; +} + +R_t find_involutions(R_t *destination, q_t q, R_t n_transformations, q_t *transformations) { + R_t n_involutions = 0; + + for (R_t i = 1; i < n_transformations; i++) { + bool is_involution = true; + for (q_t j = 0; j < q; j++) { + if (j != transformations[q * i + transformations[q * i + j]]) { + is_involution = false; + break; + } + } + if (is_involution) { + destination[n_involutions] = i; + n_involutions++; + } + } + + return n_involutions; +} + state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H) { state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t)); @@ -38,11 +72,26 @@ state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H) } s->q = 2; - s->n_transformations = 1; - s->transformations = (q_t *)malloc(2 * sizeof(q_t)); - s->transformations[0] = 1; - s->transformations[1] = 0; + s->n_transformations = 2; + s->transformations = (q_t *)malloc(2 * 2 * sizeof(q_t)); + s->transformations[0] = 0; + s->transformations[1] = 1; + s->transformations[2] = 1; + s->transformations[3] = 0; + + s->n_involutions = 1; + s->involutions = (R_t *)malloc(1 * sizeof(R_t)); + s->involutions[0] = 1; + + s->transform_site_to_zero = (R_t *)malloc(2 * sizeof(R_t)); + s->transform_site_to_zero[0] = 0; + s->transform_site_to_zero[1] = 1; + + s->n_bond_types = 2; + s->bond_with_zero_type = (q_t *)malloc(2 * sizeof(q_t)); + s->bond_with_zero_type[0] = 0; + s->bond_with_zero_type[1] = 1; s->T = T; s->J = (double *)malloc(2 * sizeof(double)); @@ -81,19 +130,34 @@ state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, doub } s->q = q; - s->n_transformations = q; - s->transformations = dihedral_gen_transformations(q); + s->n_transformations = factorial(q); + s->transformations = symmetric_gen_transformations(q); + s->involutions = (R_t *)malloc(s->n_transformations * sizeof(R_t)); + s->n_involutions = find_involutions(s->involutions, q, s->n_transformations, s->transformations); + + s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations); + + s->n_bond_types = 2; + + s->bond_with_zero_type = (q_t *)malloc(q * sizeof(q_t)); + + s->bond_with_zero_type[0] = 0; + + for (q_t i = 1; i < q; i++) { + s->bond_with_zero_type[i] = 1; + } s->T = T; - s->J = (double *)calloc(q, sizeof(double)); + s->J = (double *)calloc(2, sizeof(double)); s->J[0] = 1.0; + s->J[1] = 0.0; s->H = (double *)malloc(q * sizeof(double)); for (q_t i = 0; i < q; i++) { s->H[i] = H[i]; } - s->J_probs = Jprobs_from_J(q, T, s->J); + s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J); s->H_probs = Jprobs_from_J(q, T, s->H); s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); @@ -101,7 +165,7 @@ state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, doub s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? - s->B = (v_t *)calloc(q, sizeof(v_t)); + s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); s->B[0] = s->ne; // everyone starts in state 0, remember? return s; @@ -122,13 +186,30 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub } s->q = q; - s->n_transformations = q; + + s->n_transformations = 2 * q; s->transformations = dihedral_gen_transformations(q); + s->n_involutions = q; + s->involutions = dihedral_gen_involutions(q); + + s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations); + s->bond_with_zero_type = malloc(q * sizeof(q_t)); + + s->n_bond_types = q / 2 + 1; + + for (q_t i = 0; i < q / 2 + 1; i++) { + s->bond_with_zero_type[i] = i; + } + + for (q_t i = 1; i < (q + 1) / 2; i++) { + s->bond_with_zero_type[q - i] = i; + } + s->T = T; - s->J = (double *)malloc(q * sizeof(double)); + s->J = (double *)malloc(s->n_bond_types * sizeof(double)); - for (q_t i = 0; i < q; i++) { + for (q_t i = 0; i < s->n_bond_types; i++) { s->J[i] = cos(2 * M_PI * i / ((double)q)); } @@ -138,7 +219,7 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub s->H[i] = H[i]; } - s->J_probs = Jprobs_from_J(q, T, s->J); + s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J); s->H_probs = Jprobs_from_J(q, T, s->H); s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); @@ -146,7 +227,7 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? - s->B = (v_t *)calloc(q, sizeof(v_t)); + s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); s->B[0] = s->ne; // everyone starts in state 0, remember? return s; @@ -168,17 +249,30 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double } s->q = q; - s->n_transformations = q; + + s->n_transformations = 2 * q; s->transformations = dihedral_gen_transformations(q); + s->n_involutions = q; + s->involutions = dihedral_gen_involutions(q); - s->T = T; - s->J = (double *)malloc(q * sizeof(double)); + s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations); + s->bond_with_zero_type = malloc(q * sizeof(q_t)); + + s->n_bond_types = q / 2 + 1; for (q_t i = 0; i < q / 2 + 1; i++) { - s->J[i] = -pow(i, 2); + s->bond_with_zero_type[i] = i; } - for (q_t i = 1; i < (q + 1) / 2; i++) { - s->J[q - i] = -pow(i, 2); + + for (q_t i = 1; i < (q + 1) / 2; i++) { + s->bond_with_zero_type[(int)q - (int)i] = i; + } + + s->T = T; + s->J = (double *)malloc(s->n_bond_types * sizeof(double)); + + for (q_t i = 0; i < s->n_bond_types; i++) { + s->J[i] = -pow(i, 2); } s->H = (double *)malloc(q * sizeof(double)); @@ -186,7 +280,7 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double s->H[i] = H[i]; } - s->J_probs = Jprobs_from_J(q, T, s->J); + s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J); s->H_probs = Jprobs_from_J(q, T, s->H); s->spins = (q_t *)calloc(s->nv, sizeof(q_t)); @@ -194,6 +288,8 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double s->M = (v_t *)calloc(q, sizeof(v_t)); s->M[0] = s->nv; // everyone starts in state 0, remember? + s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t)); + s->B[0] = s->nv; // everyone starts in state 0, remember? return s; } @@ -201,8 +297,10 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double double state_finite_energy(state_finite_t *s) { double E = 0; - for (q_t i = 0; i < s->q; i++) { + for (q_t i = 0; i < s->n_bond_types; i++) { E += s->J[i] * s->B[i]; + } + for (q_t i = 0; i < s->q; i++) { E += s->H[i] * s->M[i]; } @@ -220,6 +318,9 @@ void state_finite_free(state_finite_t *s) { free(s->M); free(s->B); free(s->transformations); + free(s->involutions); + free(s->transform_site_to_zero); + free(s->bond_with_zero_type); free(s); } diff --git a/lib/symmetric.c b/lib/symmetric.c index 729b38c..4487538 100644 --- a/lib/symmetric.c +++ b/lib/symmetric.c @@ -25,3 +25,46 @@ q_t *symmetric_invert(q_t q, const q_t *g) { return g_inv; } +void swap(q_t *q1, q_t *q2) { + q_t temp = *q1; + *q1 = *q2; + *q2 = temp; +} + +R_t factorial(q_t q) { + if (q == 0) { + return 1; + } else { + return q * factorial(q - 1); + } +} + +void permute(q_t *a, q_t l, q_t r, R_t pos, q_t *transformations) { + if (l == r - 1) { + for (q_t i = 0; i < r; i++) { + transformations[r * pos + i] = a[i]; + } + } else { + for (q_t i = l; i < r; i++) { + swap((a+l), (a+i)); + permute(a, l+1, r, pos + (i - l) * factorial(r - l - 1), transformations); + swap((a+l), (a+i)); + } + } +} + +q_t *symmetric_gen_transformations(q_t q) { + q_t *transformations = (q_t *)malloc(q * factorial(q) * sizeof(q_t)); + q_t *tmp = (q_t *)malloc(q * sizeof(q_t)); + + for (q_t i = 0; i < q; i++) { + tmp[i] = i; + } + + permute(tmp, 0, q, 0, transformations); + + free(tmp); + + return transformations; +} + diff --git a/lib/symmetric.h b/lib/symmetric.h index 6e00f52..c71521d 100644 --- a/lib/symmetric.h +++ b/lib/symmetric.h @@ -11,3 +11,5 @@ q_t symmetric_act(const q_t *g, q_t s); q_t *symmetric_invert(q_t q, const q_t *g); +q_t *symmetric_gen_transformations(q_t q); + diff --git a/src/wolff_finite.c b/src/wolff_finite.c index e41c326..4bf96b9 100644 --- a/src/wolff_finite.c +++ b/src/wolff_finite.c @@ -96,7 +96,7 @@ int main(int argc, char *argv[]) { char *filename_S = (char *)malloc(256 * sizeof(char)); sprintf(filename_M, "wolff_%s_%" PRIq "_%" PRID "_%" PRIL "_%.8f", finite_model_t_strings[model], q, D, L, T); - for (q_t i = 0; i < q; i++) { + for (q_t i = 0; i < s->q; i++) { sprintf(filename_M + strlen(filename_M), "_%.8f", s->H[i]); } @@ -126,8 +126,8 @@ int main(int argc, char *argv[]) { bool changed = false; while (!changed) { - step = gsl_rng_uniform_int(r, s->n_transformations); - if (symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) { + step = gsl_rng_uniform_int(r, s->n_involutions); + if (symmetric_act(s->transformations + s->q * s->involutions[step], s->spins[v0]) != s->spins[v0]) { changed = true; } } @@ -138,9 +138,9 @@ int main(int argc, char *argv[]) { // as a fast time many machines will actually have it be 64 bits. we cast // it down here to halve space. - for (q_t i = 0; i < q - 1; i++) { - fwrite(&(s->M[i]), sizeof(uint32_t), 1, outfile_M); // if we know the occupation of the first q - 1 states, we know the occupation of the last - fwrite(&(s->B[i]), sizeof(uint32_t), 1, outfile_B); // if we know the occupation of the first q - 1 states, we know the occupation of the last + for (q_t i = 0; i < s->q - 1; i++) { // if we know the occupation of the first q - 1 states, we know the occupation of the last + fwrite(&(s->M[i]), sizeof(uint32_t), 1, outfile_M); + fwrite(&(s->B[i]), sizeof(uint32_t), 1, outfile_B); } fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S); -- cgit v1.2.3-54-g00ecf From 2d8fcebf2f56efd1c3913ba49eaff6520ffdb33d Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 6 Jul 2018 14:42:44 -0400 Subject: rewrote wolff in c++ with templates so that any system can be run with it --- CMakeLists.txt | 15 +- lib/cluster.c | 277 ------------------------------------- lib/cluster.h | 391 +++++++++++++++++++++++++++++++++++++++++++++-------- lib/graph.h | 10 +- lib/rand.h | 9 ++ lib/stack.h | 9 ++ src/wolff.cpp | 137 +++++++++++++++++++ src/wolff_dgm.c | 247 --------------------------------- src/wolff_vector.c | 377 --------------------------------------------------- 9 files changed, 507 insertions(+), 965 deletions(-) delete mode 100644 lib/cluster.c create mode 100644 src/wolff.cpp delete mode 100644 src/wolff_dgm.c delete mode 100644 src/wolff_vector.c (limited to 'lib') diff --git a/CMakeLists.txt b/CMakeLists.txt index 4b46bcd..7c1bc32 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,15 +4,15 @@ project(wolff) set(CMAKE_CXX_FLAGS_DEBUG "-g") set(CMAKE_CXX_FLAGS_RELEASE "-O3") -set (CMAKE_CXX_STANDARD 11) +set (CMAKE_CXX_STANDARD 17) include_directories(lib ~/.local/include) link_directories(~/.local/lib) -file(GLOB SOURCES lib/*.c) -add_executable(wolff_finite src/wolff_finite.c ${SOURCES}) -add_executable(wolff_vector src/wolff_vector.c ${SOURCES}) -add_executable(wolff_dgm src/wolff_dgm.c ${SOURCES}) +file(GLOB CSOURCES lib/*.c) +file(GLOB CPPSOURCES lib/*.cpp) +add_executable(wolff_finite src/wolff_finite.c ${CSOURCES}) +add_executable(wolff src/wolff.cpp ${CPPSOURCES} ${CSOURCES}) find_package(OpenMP) if (OPENMP_FOUND) @@ -21,8 +21,7 @@ if (OPENMP_FOUND) endif() target_link_libraries(wolff_finite gsl m cblas fftw3) -target_link_libraries(wolff_vector gsl m cblas fftw3) -target_link_libraries(wolff_dgm gsl m cblas fftw3) +target_link_libraries(wolff gsl m cblas fftw3) -install(TARGETS wolff_finite wolff_vector wolff_dgm DESTINATION bin) +install(TARGETS wolff_finite wolff DESTINATION bin) diff --git a/lib/cluster.c b/lib/cluster.c deleted file mode 100644 index 96225a2..0000000 --- a/lib/cluster.c +++ /dev/null @@ -1,277 +0,0 @@ - -#include "cluster.h" - -v_t flip_cluster_dgm(dgm_state_t *s, v_t v0, h_t rot, gsl_rng *r) { - v_t nv = 0; - - ll_t *stack = NULL; // create a new stack - stack_push(&stack, v0); // push the initial vertex to the stack - - bool *marks = (bool *)calloc(s->g->nv, sizeof(bool)); - - while (stack != NULL) { - v_t v = stack_pop(&stack); - - if (!marks[v]) { - h_t s_old, s_new; - dihinf_t *R_new; - bool external_flipped; - - marks[v] = true; - - if (v == s->g->nv - 1) { - R_new = dihinf_compose(rot, s->R); - external_flipped = true; - } else { - s_old = s->spins[v]; - s_new = dihinf_act(rot, s_old); - external_flipped = false; - } - - v_t nn = s->g->v_i[v + 1] - s->g->v_i[v]; - - for (v_t i = 0; i < nn; i++) { - h_t sn; - double prob; - bool external_neighbor = false; - - v_t vn = s->g->v_adj[s->g->v_i[v] + i]; - - if (vn == s->g->nv - 1) { - external_neighbor = true; - } else { - sn = s->spins[vn]; - } - - if (external_flipped || external_neighbor) { - h_t rot_s_old, rot_s_new; - - if (external_neighbor) { - rot_s_old = dihinf_inverse_act(s->R, s_old); - rot_s_new = dihinf_inverse_act(s->R, s_new); - } else { - rot_s_old = dihinf_inverse_act(s->R, sn); - rot_s_new = dihinf_inverse_act(R_new, sn); - } - - double dE = s->H(s->H_info, rot_s_old) - s->H(s->H_info, rot_s_new); - prob = 1.0 - exp(-dE / s->T); - - s->M += rot_s_new - rot_s_old; - s->E += dE; - } else { - double dE = (s->J)(s_old - sn) - (s->J)(s_new - sn); - prob = 1.0 - exp(-dE / s->T); - s->E += dE; - } - - if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... - stack_push(&stack, vn); // push the neighboring vertex to the stack - } - } - - if (external_flipped) { - free(s->R); - s->R = R_new; - } else { - s->spins[v] = s_new; - } - - if (v != s->g->nv - 1) { // count the number of non-external sites that flip - nv++; - } - } - } - - free(marks); - - return nv; -} - -v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) { - v_t nv = 0; - - ll_t *stack = NULL; // create a new stack - stack_push(&stack, v0); // push the initial vertex to the stack - - //node_t *T = NULL; - bool *marks = (bool *)calloc(s->g->nv, sizeof(bool)); - - while (stack != NULL) { - v_t v = stack_pop(&stack); - -// if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped - if (!marks[v]) { - bool v_is_external = false; - double *s_old, *s_new, *R_tmp; - - if (v == s->g->nv - 1) { - v_is_external = true; - } - - //tree_insert(&T, v); - marks[v] = true; - - if (v == s->g->nv - 1) { - R_tmp = orthogonal_rotate(s->n, rot, s->R); - } else { - s_old = &(s->spins[s->n * v]); // don't free me! I'm a pointer within array s->spins - s_new = vector_rotate(s->n, rot, s_old); // free me! I'm a new vector - } - - v_t nn = s->g->v_i[v + 1] - s->g->v_i[v]; - - for (v_t i = 0; i < nn; i++) { - v_t vn = s->g->v_adj[s->g->v_i[v] + i]; - - bool vn_is_external = false; - - if (vn == s->g->nv - 1) { - vn_is_external = true; - } - - double *sn; - - if (!vn_is_external) { - sn = &(s->spins[s->n * vn]); - } - - double prob; - - if (v_is_external || vn_is_external) { - double *rs_old, *rs_new; - if (vn_is_external) { - rs_old = vector_rotate_inverse(s->n, s->R, s_old); - rs_new = vector_rotate_inverse(s->n, s->R, s_new); - } else { - rs_old = vector_rotate_inverse(s->n, s->R, sn); - rs_new = vector_rotate_inverse(s->n, R_tmp, sn); - } - double dE = s->H(s->n, s->H_info, rs_old) - s->H(s->n, s->H_info, rs_new); - prob = 1.0 - exp(-dE / s->T); - vector_subtract(s->n, s->M, rs_old); - vector_add(s->n, s->M, rs_new); - s->E += dE; - - free(rs_old); - free(rs_new); - } else { - double dE = (s->J)(vector_dot(s->n, sn, s_old)) - (s->J)(vector_dot(s->n, sn, s_new)); - prob = 1.0 - exp(-dE / s->T); - s->E += dE; - } - - if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... - stack_push(&stack, vn); // push the neighboring vertex to the stack - } - } - - if (v == s->g->nv - 1) { - free(s->R); - s->R = R_tmp; - } else { - vector_replace(s->n, s_old, s_new); - free(s_new); - } - - if (v != s->g->nv - 1) { // count the number of non-external sites that flip - nv++; - } - } - } - - //tree_freeNode(T); - free(marks); - - return nv; -} - -/*G -template -v_t flip_cluster(state_t *state, v_t v0, R_t *r, gsl_rng *rand) { - v_t nv = 0; - - 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)); - - while (stack != NULL) { - v_t v = stack_pop(&stack); - - if (!marks[v]) { - X_t *si_old, *si_new; - R_t *s0_old, *s0_new; - - si_old = state->s[v]; - s0_old = state->s0; - - marks[v] = true; - - if (v == state->g->nv - 1) { - s0_new = act (r, s0_old); - } else { - si_new = act (r, si_old); - } - - v_t nn = state->g->v_i[v + 1] - state->g->v_i[v]; - - for (v_t i = 0; i < nn; i++) { - v_t vn = state->g->v_adj[state->g->v_i[v] + i]; - - X_t *sj; - - if (vn != state->g->nv - 1) { - sj = state->s[vn]; - } - - double prob; - - bool is_ext = (v == state->g->nv - 1 || vn == state->g->nv - 1); - - if (is_ext) { - X_t *rs_old, *rs_new; - if (vn == state->g->nv - 1) { - rs_old = inverse_act (s0_old, si_old); - rs_new = inverse_act (s0_old, si_new); - } else { - rs_old = inverse_act (s0_old, sj); - rs_new = inverse_act (s0_new, sj); - } - double dE = state->B(rs_old) - state->B(rs_new); - prob = 1.0 - exp(-dE / state->T); - update_magnetization (state->M, rs_old, rs_new); - state->E += dE; - - free_X (rs_old); - free_X (rs_new); - } else { - double dE = state->Z(si_old, sj) - state->Z(si_new, sj); - prob = 1.0 - exp(-dE / state->T); - state->E += dE; - } - - if (gsl_rng_uniform(rand) < prob) { // and with probability... - stack_push(&stack, vn); // push the neighboring vertex to the stack - } - } - - if (v == state->g->nv - 1) { - free_R (state->s0); - state->s0 = s0_new; - } else { - free_X (state->s[v]); - state->s[v] = si_new; - } - - if (v != state->g->nv - 1) { // count the number of non-external sites that flip - nv++; - } - } - } - - free(marks); - - return nv; -} -*/ diff --git a/lib/cluster.h b/lib/cluster.h index d118735..29dd0cb 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -1,13 +1,14 @@ #pragma once +#include #include #include #include #include #include #include -#include +#include #include #include #include @@ -24,57 +25,339 @@ #include "dihinf.h" #include "yule_walker.h" -typedef struct { - graph_t *g; - q_t *spins; - double T; - double *J; - double *H; - double *J_probs; - double *H_probs; - dihedral_t *R; - double E; - v_t *M; - q_t q; -} ising_state_t; - -typedef struct { - graph_t *g; - h_t *spins; - double T; - double (*J)(h_t); - double (*H)(double *, h_t); - double *H_info; - dihinf_t *R; - double E; - h_t M; -} dgm_state_t; - -typedef struct { - graph_t *g; - double *spins; - double T; - double (*J)(double); - double (*H)(q_t, double *, double *); - double *H_info; - double *R; - double E; - double *M; - q_t n; -} vector_state_t; - -typedef enum { - VECTOR, - MODULATED, - CUBIC, - QUADRATIC -} vector_field_t; - -v_t flip_cluster(ising_state_t *s, v_t v0, q_t s1, gsl_rng *r); - -v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r); - -v_t flip_cluster_dgm(dgm_state_t *s, v_t v0, h_t rot, gsl_rng *r); - -graph_t *graph_add_ext(const graph_t *g); +template +void init(T*); + +template +T scalar_multiple(v_t a, T b); + +template +X_t act(R_t a, X_t b); + +template +X_t act_inverse(R_t a, X_t b); + +template +T copy(T a); + +template +void free_spin(T a); + +template +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 +struct vector_t { T *x; }; + +template +void init(vector_t *ptr) { + ptr->x = (T *)calloc(q, sizeof(T)); + + ptr->x[0] = (T)1; +} + +template +vector_t copy (vector_t v) { + vector_t v_copy; + + v_copy.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + v_copy.x[i] = v.x[i]; + } + + 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 subtract (vector_t v1, vector_t v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] -= v2.x[i]; + } +} + +template +vector_t scalar_multiple(v_t a, vector_t v) { + vector_t multiple; + multiple.x = (T *)malloc(q * sizeof(T)); + for (q_t i = 0; i < q; i++) { + multiple.x[i] = a * v.x[i]; + } + + return multiple; +} + +template +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 +void free_spin (vector_t v) { + free(v.x); +} + +template +struct orthogonal_t { T *x; }; + +template +void init(orthogonal_t *ptr) { + ptr->x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + ptr->x[q * i + i] = (T)1; + } +} + +template +orthogonal_t copy (orthogonal_t m) { + orthogonal_t m_copy; + m_copy.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q * q; i++) { + m_copy.x[i] = m.x[i]; + } + + return m_copy; +} + +template +void free_spin (orthogonal_t m) { + free(m.x); +} + +template +vector_t act (orthogonal_t m, vector_t v) { + vector_t v_rot; + v_rot.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot.x[i] += m.x[q * i + j] * v.x[j]; + } + } + + return v_rot; +} + + +template +orthogonal_t act (orthogonal_t m1, orthogonal_t m2) { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k]; + } + } + } + + return m2_rot; +} + +template +vector_t act_inverse (orthogonal_t m, vector_t v) { + vector_t v_rot; + v_rot.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot.x[i] += m.x[q * j + i] * v.x[j]; + } + } + + return v_rot; +} + +template +orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k]; + } + } + } + + return m2_rot; +} + +template +void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { + double *v = (double *)malloc(q * sizeof(double)); + double v2 = 0; + + for (q_t i = 0; i < q; i++) { + v[i] = gsl_ran_ugaussian(r); + v2 += v[i] * v[i]; + } + + ptr->x = (double *)calloc(q * q, sizeof(double)); + + for (q_t i = 0; i < q; i++) { + ptr->x[q * i + i] = 1.0; + for (q_t j = 0; j < q; j++) { + ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2; + } + } + + free(v); +} + +template +v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { + v_t nv = 0; + + 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)); + + while (stack != NULL) { + v_t v = stack_pop(&stack); + + if (!marks[v]) { + 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) { + R_new = act (r, R_old); + } else { + si_new = act (r, si_old); + } + + v_t nn = state->g->v_i[v + 1] - state->g->v_i[v]; + + for (v_t i = 0; i < nn; i++) { + v_t vn = state->g->v_adj[state->g->v_i[v] + i]; + + X_t sj; + + if (vn != state->g->nv - 1) { + sj = state->spins[vn]; + } + + double prob; + + bool is_ext = (v == state->g->nv - 1 || vn == state->g->nv - 1); + + if (is_ext) { + X_t rs_old, rs_new; + if (vn == state->g->nv - 1) { + rs_old = act_inverse (R_old, si_old); + rs_new = act_inverse (R_old, si_new); + } else { + rs_old = act_inverse (R_old, sj); + rs_new = act_inverse (R_new, sj); + } + 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); + state->E += dE; + + free_spin (rs_old); + free_spin (rs_new); + } else { + double dE = state->J(si_old, sj) - state->J(si_new, sj); + prob = 1.0 - exp(-dE / state->T); + state->E += dE; + } + + if (gsl_rng_uniform(rand) < prob) { // and with probability... + stack_push(&stack, vn); // push the neighboring vertex to the stack + } + } + + if (v == state->g->nv - 1) { + free_spin(state->R); + state->R = R_new; + } else { + free_spin(state->spins[v]); + state->spins[v] = si_new; + } + + if (v != state->g->nv - 1) { // count the number of non-external sites that flip + nv++; + } + } + } + + free(marks); + + return nv; +} diff --git a/lib/graph.h b/lib/graph.h index cb47faa..beb7f4c 100644 --- a/lib/graph.h +++ b/lib/graph.h @@ -7,6 +7,10 @@ #include "types.h" +#ifdef __cplusplus +extern "C" { +#endif + typedef struct { v_t ne; v_t nv; @@ -15,8 +19,10 @@ typedef struct { } graph_t; graph_t *graph_create_square(D_t D, L_t L); - graph_t *graph_add_ext(const graph_t *G); - void graph_free(graph_t *h); +#ifdef __cplusplus +} +#endif + diff --git a/lib/rand.h b/lib/rand.h index 2354f6a..7bb5354 100644 --- a/lib/rand.h +++ b/lib/rand.h @@ -4,4 +4,13 @@ #include #include +#ifdef __cplusplus +extern "C" { +#endif + unsigned long int rand_seed(); + +#ifdef __cplusplus +} +#endif + diff --git a/lib/stack.h b/lib/stack.h index a354ab5..8d25aff 100644 --- a/lib/stack.h +++ b/lib/stack.h @@ -8,6 +8,11 @@ #include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + typedef struct ll_tag { v_t x; struct ll_tag *next; @@ -24,3 +29,7 @@ void stack_push_d(dll_t **q, double x); v_t stack_pop(ll_t **q); double stack_pop_d(dll_t **q); +#ifdef __cplusplus +} +#endif + diff --git a/src/wolff.cpp b/src/wolff.cpp new file mode 100644 index 0000000..85df357 --- /dev/null +++ b/src/wolff.cpp @@ -0,0 +1,137 @@ + +#include +#include + +#include + +double H_vector(vector_t <2, double> v1, double *H) { + vector_t <2, double> H_vec; + H_vec.x = H; + return dot <2, double> (v1, H_vec); +} + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e7; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double *H = (double *)calloc(MAX_Q, sizeof(double)); + + bool silent = false; + + int opt; + q_t J_ind = 0; + q_t H_ind = 0; + + while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:s")) != -1) { + switch (opt) { + case 'N': // number of steps + N = (count_t)atof(optarg); + break; + case 'D': // dimension + D = atoi(optarg); + break; + case 'L': // linear size + L = atoi(optarg); + break; + case 'T': // temperature + T = atof(optarg); + break; + case 'H': // external field. nth call couples to state n + H[H_ind] = atof(optarg); + H_ind++; + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + state_t , vector_t <2, double>> s(D, L, T, dot <2, double>, std::bind(H_vector, std::placeholders::_1, H)); + + // initialize random number generator + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(r, rand_seed()); + + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); + + fprintf(outfile_info, "<| \"ID\" -> %lu, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, D, L, s.nv, s.ne, T); + + for (q_t i = 0; i < 2; i++) { + fprintf(outfile_info, "%.15f", H[i]); + if (i < 2 - 1) { + fprintf(outfile_info, ", "); + } + } + + fprintf(outfile_info, "} |>\n"); + + fclose(outfile_info); + + char *filename_M = (char *)malloc(255 * sizeof(char)); + char *filename_E = (char *)malloc(255 * sizeof(char)); + char *filename_S = (char *)malloc(255 * sizeof(char)); + + sprintf(filename_M, "wolff_%lu_M.dat", timestamp); + sprintf(filename_E, "wolff_%lu_E.dat", timestamp); + sprintf(filename_S, "wolff_%lu_S.dat", timestamp); + + FILE *outfile_M = fopen(filename_M, "wb"); + FILE *outfile_E = fopen(filename_E, "wb"); + FILE *outfile_S = fopen(filename_S, "wb"); + + free(filename_M); + free(filename_E); + free(filename_S); + + v_t cluster_size = 0; + + if (!silent) printf("\n"); + for (count_t steps = 0; steps < N; steps++) { + if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, M_0 = %.2f, S = %" PRIv "\n", steps, N, s.E, s.M.x[0], cluster_size); + + v_t v0 = gsl_rng_uniform_int(r, s.nv); + + orthogonal_t <2, double> step; + generate_rotation<2>(r, &step); + + printf("(%g %g) . (%g %g) = %g or %g, H = %g\n\n", s.spins[0].x[0], s.spins[0].x[1], s.spins[1].x[0], s.spins[1].x[1], dot(s.spins[0], s.spins[1]), s.J(s.spins[0],s.spins[1]), s.H(s.spins[0])); + + getchar(); + cluster_size = flip_cluster , vector_t <2, double>> (&s, v0, step, r); + + free_spin(step); + + fwrite(&(s.E), sizeof(double), 1, outfile_E); + fwrite(s.M.x, sizeof(double), 2, outfile_M); + fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S); + + } + if (!silent) { + printf("\033[F\033[J"); + } + printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, M_0 = %.2f, S = %" PRIv "\n", N, N, s.E, s.M.x[0], cluster_size); + + fclose(outfile_M); + fclose(outfile_E); + fclose(outfile_S); + + gsl_rng_free(r); + + free(H); + + return 0; +} + diff --git a/src/wolff_dgm.c b/src/wolff_dgm.c deleted file mode 100644 index f11b296..0000000 --- a/src/wolff_dgm.c +++ /dev/null @@ -1,247 +0,0 @@ - -#include - -#include - -double identity(h_t x) { - return -pow(x, 2); -} - -double basic_H(double *H, h_t x) { - return -H[0] * pow(x, 2); -} - -int main(int argc, char *argv[]) { - - L_t L = 128; - count_t N = (count_t)1e7; - count_t min_runs = 10; - count_t n = 3; - D_t D = 2; - double T = 2.26918531421; - double *H = (double *)calloc(MAX_Q, sizeof(double)); - double eps = 0; - bool silent = false; - bool record_autocorrelation = false; - count_t ac_skip = 1; - count_t W = 10; - - int opt; - q_t H_ind = 0; - - while ((opt = getopt(argc, argv, "N:n:D:L:T:H:m:e:saS:W:")) != -1) { - switch (opt) { - case 'N': - N = (count_t)atof(optarg); - break; - case 'n': - n = (count_t)atof(optarg); - break; - case 'D': - D = atoi(optarg); - break; - case 'L': - L = atoi(optarg); - break; - case 'T': - T = atof(optarg); - break; - case 'H': - H[H_ind] = atof(optarg); - H_ind++; - break; - case 'm': - min_runs = atoi(optarg); - break; - case 'e': - eps = atof(optarg); - break; - case 's': - silent = true; - break; - case 'a': - record_autocorrelation = true; - break; - case 'S': - ac_skip = (count_t)atof(optarg); - break; - case 'W': - W = (count_t)atof(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - - dgm_state_t *s = (dgm_state_t *)calloc(1, sizeof(dgm_state_t)); - - graph_t *h = graph_create_square(D, L); - s->g = graph_add_ext(h); - - s->spins = (h_t *)calloc(h->nv, sizeof(h_t)); - - s->H_info = H; - s->T = T; - s->H = basic_H; - s->J = identity; - - s->R = (dihinf_t *)calloc(1, sizeof(dihinf_t)); - - s->M = 0; - s->E = 0; - - double diff = 1e31; - count_t n_runs = 0; - count_t n_steps = 0; - - meas_t *E, *clust, *M, *dM; - - M = (meas_t *)calloc(1, sizeof(meas_t )); - dM = (meas_t *)calloc(1, sizeof(meas_t )); - - E = calloc(1, sizeof(meas_t)); - clust = calloc(1, sizeof(meas_t)); - - autocorr_t *autocorr; - if (record_autocorrelation) { - autocorr = (autocorr_t *)calloc(1, sizeof(autocorr_t)); - autocorr->W = 2 * W + 1; - autocorr->OO = (double *)calloc(2 * W + 1, sizeof(double)); - } - - if (!silent) printf("\n"); - while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) { - if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 - ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", - n_runs, fabs(meas_dx(E) / E->x), meas_dx(M) / M->x, meas_dc(E) / meas_c(E), meas_dc(M) / meas_c(M), h->nv / clust->x); - - count_t n_flips = 0; - - while (n_flips / h->nv < n) { - v_t v0 = gsl_rng_uniform_int(r, h->nv); - h_t step = round((((double)s->M) / h->nv) + gsl_ran_gaussian(r, 5)); - - v_t tmp_flips = flip_cluster_dgm(s, v0, step, r); - n_flips += tmp_flips; - - if (n_runs > 0) { - n_steps++; - meas_update(clust, tmp_flips); - } - - if (record_autocorrelation && n_runs > 0) { - if (n_steps % ac_skip == 0) { - update_autocorr(autocorr, s->E); - } - } - } - - meas_update(M, s->M); - h_t min_h, max_h; - min_h = MAX_H; - max_h = MIN_H; - for (v_t i = 0; i < h->nv; i++) { - if (s->spins[i] < min_h) { - min_h = s->spins[i]; - } else if (s->spins[i] > max_h) { - max_h = s->spins[i]; - } - } - meas_update(dM, max_h - min_h); - meas_update(E, s->E); - - diff = fabs(meas_dc(E) / meas_c(E)); - - n_runs++; - } - if (!silent) { - printf("\033[F\033[J"); - } - printf("WOLFF: sweep %" PRIu64 - ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", - n_runs, fabs(meas_dx(E) / E->x), meas_dx(M) / M->x, meas_dc(E) / meas_c(E), meas_dc(M) / meas_c(M), h->nv / clust->x); - - double tau = 0; - bool tau_failed = false; - - if (record_autocorrelation) { - double *Gammas = (double *)malloc((W + 1) * sizeof(double)); - - Gammas[0] = 1 + rho(autocorr, 0); - for (uint64_t i = 0; i < W; i++) { - Gammas[1 + i] = rho(autocorr, 2 * i + 1) + rho(autocorr, 2 * i + 2); - } - - uint64_t n; - for (n = 0; n < W + 1; n++) { - if (Gammas[n] <= 0) { - break; - } - } - - if (n == W + 1) { - printf("WARNING: correlation function never hit the noise floor.\n"); - tau_failed = true; - } - - if (n < 2) { - printf("WARNING: correlation function only has one nonnegative term.\n"); - tau_failed = true; - } - - double *conv_Gamma = get_convex_minorant(n, Gammas); - - double ttau = - 0.5; - - for (uint64_t i = 0; i < n + 1; i++) { - ttau += conv_Gamma[i]; - } - - free(Gammas); - free(autocorr->OO); - while (autocorr->Op != NULL) { - stack_pop_d(&(autocorr->Op)); - } - free(autocorr); - - tau = ttau * ac_skip * clust->x / h->nv; - } - - if (tau_failed) { - tau = 0; - } - - FILE *outfile = fopen("out.m", "a"); - - fprintf(outfile, "<|D->%" PRID ",L->%" PRIL ",T->%.15f", D, L, T); - fprintf(outfile, ",E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->%.15f", E->x / h->nv, meas_dx(E) / h->nv, meas_c(E) / h->nv, meas_dc(E) / h->nv, M->x / h->nv); - fprintf(outfile, ",\\[Delta]M->%.15f", meas_dx(M) / h->nv); - fprintf(outfile, ",\\[Chi]->%.15f", meas_c(M) / h->nv); - fprintf(outfile, ",\\[Delta]\\[Chi]->%.15f", meas_dc(M) / h->nv); - fprintf(outfile, ",w->%.15f,\\[Delta]w->%.15f,wc->%.15f,\\[Delta]wc->%.15f,Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", dM->x, meas_dx(dM), meas_c(dM), meas_dc(dM), clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau); - - fclose(outfile); - - FILE *image = fopen("out.dat", "a"); - for (v_t i = 0; i < h->nv; i++) { - fprintf(image, "%" PRIh " ", s->spins[i]); - } - fprintf(image, "\n"); - fclose(image); - - free(E); - free(clust); - free(H); - free(s->R); - free(s->spins); - graph_free(s->g); - free(s); - graph_free(h); - gsl_rng_free(r); - - return 0; -} - diff --git a/src/wolff_vector.c b/src/wolff_vector.c deleted file mode 100644 index c5ebcb5..0000000 --- a/src/wolff_vector.c +++ /dev/null @@ -1,377 +0,0 @@ - -#include - -#include - -double identity(double x) { - return x; -} - -double zero(q_t n, double *H, double *x) { - return 0.0; -} - -double dot(q_t n, double *H, double *x) { - double total = 0; - for (q_t i = 0; i < n; i++) { - total += H[i] * x[i]; - } - return total; -} - -double theta(double x, double y) { - double val = atan(y / x); - - if (x < 0.0 && y > 0.0) { - return M_PI + val; - } else if ( x < 0.0 && y < 0.0 ) { - return - M_PI + val; - } else { - return val; - } -} - -double modulated(q_t n, double *H_info, double *x) { - return H_info[0] * cos(H_info[1] * theta(x[0], x[1])); -} - -double cubic(q_t n, double *H_info, double *x) { - double v_sum = 0; - - for (q_t i = 0; i < n; i++) { - v_sum += pow(x[i], 4); - } - - return - H_info[0] * v_sum; -} - -double quadratic(q_t n, double *H_info, double *x) { - double tmp = 0; - - tmp += pow(x[0], 2); - - for (q_t i = 1; i < n; i++) { - tmp += - 1.0 / (n - 1.0) * pow(x[i], 2); - } - - return - 0.5 * H_info[0] * tmp; -} - -int main(int argc, char *argv[]) { - - L_t L = 128; - count_t N = (count_t)1e7; - count_t min_runs = 10; - count_t n = 3; - q_t q = 2; - D_t D = 2; - double T = 2.26918531421; - double *H = (double *)calloc(MAX_Q, sizeof(double)); - double eps = 0; - bool silent = false; - bool record_autocorrelation = false; - vector_field_t H_type = VECTOR; - count_t ac_skip = 1; - count_t W = 10; - - int opt; - q_t H_ind = 0; - - while ((opt = getopt(argc, argv, "N:n:D:L:q:T:H:m:e:saS:W:f:")) != -1) { - switch (opt) { - case 'N': - N = (count_t)atof(optarg); - break; - case 'n': - n = (count_t)atof(optarg); - break; - case 'D': - D = atoi(optarg); - break; - case 'L': - L = atoi(optarg); - break; - case 'q': - q = atoi(optarg); - break; - case 'T': - T = atof(optarg); - break; - case 'H': - H[H_ind] = atof(optarg); - H_ind++; - break; - case 'm': - min_runs = atoi(optarg); - break; - case 'e': - eps = atof(optarg); - break; - case 's': - silent = true; - break; - case 'a': - record_autocorrelation = true; - break; - case 'S': - ac_skip = (count_t)atof(optarg); - break; - case 'W': - W = (count_t)atof(optarg); - break; - case 'f': - switch (atoi(optarg)) { - case 0: - H_type = VECTOR; - break; - case 1: - H_type = MODULATED; - break; - case 2: - H_type = CUBIC; - break; - case 3: - H_type = QUADRATIC; - break; - } - break; - default: - exit(EXIT_FAILURE); - } - } - - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - - vector_state_t *s = (vector_state_t *)calloc(1, sizeof(vector_state_t)); - - graph_t *h = graph_create_square(D, L); - s->g = graph_add_ext(h); - - s->n = q; - - s->spins = (double *)calloc(n * h->nv, sizeof(double)); - - for (v_t i = 0; i < h->nv; i++) { - s->spins[q * i] = 1.0; - } - - s->H_info = H; - s->T = T; - switch (H_type) { - case VECTOR: - s->H = dot; - break; - case MODULATED: - s->H = modulated; - break; - case CUBIC: - s->H = cubic; - break; - case QUADRATIC: - s->H = quadratic; - break; - } - s->J = identity; - - s->R = (double *)calloc(q * q, sizeof(double)); - - for (q_t i = 0; i < q; i++) { - s->R[q * i + i] = 1.0; - } - - s->M = (double *)calloc(q, sizeof(double)); - s->M[0] = 1.0; - s->E = - ((double)h->ne) * s->J(1.0) - (double)h->nv * s->H(s->n, s->H_info, s->M); - s->M[0] *= (double)h->nv; - - double diff = 1e31; - count_t n_runs = 0; - count_t n_steps = 0; - - meas_t *E, *clust, **M, *aM; - - M = (meas_t **)malloc(q * sizeof(meas_t *)); - aM = (meas_t *)calloc(q, sizeof(meas_t )); - for (q_t i = 0; i < q; i++) { - M[i] = (meas_t *)calloc(1, sizeof(meas_t)); - } - - E = calloc(1, sizeof(meas_t)); - clust = calloc(1, sizeof(meas_t)); - - autocorr_t *autocorr; - if (record_autocorrelation) { - autocorr = (autocorr_t *)calloc(1, sizeof(autocorr_t)); - autocorr->W = 2 * W + 1; - autocorr->OO = (double *)calloc(2 * W + 1, sizeof(double)); - } - - if (!silent) printf("\n"); - while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) { - if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 - ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", - n_runs, fabs(meas_dx(E) / E->x), meas_dx(aM) / aM->x, meas_dc(E) / meas_c(E), meas_dc(aM) / meas_c(aM), h->nv / clust->x); - - count_t n_flips = 0; - - while (n_flips / h->nv < n) { - v_t v0 = gsl_rng_uniform_int(r, h->nv); - double *step = gen_rot(r, q); - - v_t tmp_flips = flip_cluster_vector(s, v0, step, r); - free(step); - n_flips += tmp_flips; - - if (n_runs > 0) { - n_steps++; - meas_update(clust, tmp_flips); - - if (record_autocorrelation && n_steps % ac_skip == 0) { - update_autocorr(autocorr, s->E); - } - } - } - - double aM_val = 0; - - for (q_t i = 0; i < q; i++) { - meas_update(M[i], s->M[i]); - aM_val += s->M[i] * s->M[i]; - } - - meas_update(aM, sqrt(aM_val)); - meas_update(E, s->E); - - diff = fabs(meas_dx(clust) / clust->x); - - n_runs++; - } - - if (!silent) { - printf("\033[F\033[J"); - } - printf("WOLFF: sweep %" PRIu64 - ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", - n_runs, fabs(meas_dx(E) / E->x), meas_dx(M[0]) / M[0]->x, meas_dc(E) / meas_c(E), meas_dc(M[0]) / meas_c(M[0]), h->nv / clust->x); - - double tau = 0; - bool tau_failed = false; - - if (record_autocorrelation) { - double *Gammas = (double *)malloc((W + 1) * sizeof(double)); - - Gammas[0] = 1 + rho(autocorr, 0); - for (uint64_t i = 0; i < W; i++) { - Gammas[1 + i] = rho(autocorr, 2 * i + 1) + rho(autocorr, 2 * i + 2); - } - - uint64_t n; - for (n = 0; n < W + 1; n++) { - if (Gammas[n] <= 0) { - break; - } - } - - if (n == W + 1) { - printf("WARNING: correlation function never hit the noise floor.\n"); - tau_failed = true; - } - - if (n < 2) { - printf("WARNING: correlation function only has one nonnegative term.\n"); - tau_failed = true; - } - - double *conv_Gamma = get_convex_minorant(n, Gammas); - - double ttau = - 0.5; - - for (uint64_t i = 0; i < n + 1; i++) { - ttau += conv_Gamma[i]; - } - - FILE *autocorr_file = fopen("autocorr.dat", "a"); - - printf("%g %g\n", Gammas[0], conv_Gamma[0]); - - for (count_t i = 0; i < n+1; i++) { - fprintf(autocorr_file, "%g ", conv_Gamma[i]); - } - fprintf(autocorr_file, "\n"); - - fclose(autocorr_file); - - free(Gammas); - free(autocorr->OO); - while (autocorr->Op != NULL) { - stack_pop_d(&(autocorr->Op)); - } - free(autocorr); - - tau = ttau * ac_skip * clust->x / h->nv; - } - - if (tau_failed) { - tau = 0; - } - - FILE *outfile = fopen("out.m", "a"); - - fprintf(outfile, "<|N->%" PRIcount ",n->%" PRIcount ",D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,H->{", N, n, D, L, q, T); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", H[i]); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->{", E->x / h->nv, meas_dx(E) / h->nv, meas_c(E) / h->nv, meas_dc(E) / h->nv); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", M[i]->x / h->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},\\[Delta]M->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", meas_dx(M[i]) / h->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},\\[Chi]->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", meas_c(M[i]) / h->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},\\[Delta]\\[Chi]->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", meas_dc(M[i]) / h->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},aM->%.15f,\\[Delta]aM->%.15f,a\\[Chi]->%.15f,\\[Delta]a\\[Chi]->%.15f,Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", aM->x / h->nv, meas_dx(aM) / h->nv, meas_c(aM) / h->nv, meas_dc(aM) / h->nv, clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau); - - fclose(outfile); - - free(E); - free(clust); - for (q_t i = 0; i < q; i++) { - free(M[i]); - } - free(M); - free(H); - free(s->M); - free(s->R); - free(s->spins); - graph_free(s->g); - free(s); - graph_free(h); - gsl_rng_free(r); - - return 0; -} - -- cgit v1.2.3-54-g00ecf From 31f4244352b5e68eed770090419541d469f7f999 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 6 Jul 2018 14:51:29 -0400 Subject: split up some files --- lib/cluster.h | 177 +------------------------------------------------------ lib/orthogonal.c | 99 ------------------------------- lib/orthogonal.h | 116 +++++++++++++++++++++++++++++++++--- lib/vector.h | 72 ++++++++++++++++++++++ 4 files changed, 180 insertions(+), 284 deletions(-) delete mode 100644 lib/orthogonal.c create mode 100644 lib/vector.h (limited to 'lib') diff --git a/lib/cluster.h b/lib/cluster.h index 29dd0cb..a20912e 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -20,6 +20,7 @@ #include "graph.h" #include "tree.h" #include "measurement.h" +#include "vector.h" #include "orthogonal.h" #include "dihedral.h" #include "dihinf.h" @@ -95,182 +96,6 @@ class state_t { } }; -template -struct vector_t { T *x; }; - -template -void init(vector_t *ptr) { - ptr->x = (T *)calloc(q, sizeof(T)); - - ptr->x[0] = (T)1; -} - -template -vector_t copy (vector_t v) { - vector_t v_copy; - - v_copy.x = (T *)calloc(q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - v_copy.x[i] = v.x[i]; - } - - 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 subtract (vector_t v1, vector_t v2) { - for (q_t i = 0; i < q; i++) { - v1.x[i] -= v2.x[i]; - } -} - -template -vector_t scalar_multiple(v_t a, vector_t v) { - vector_t multiple; - multiple.x = (T *)malloc(q * sizeof(T)); - for (q_t i = 0; i < q; i++) { - multiple.x[i] = a * v.x[i]; - } - - return multiple; -} - -template -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 -void free_spin (vector_t v) { - free(v.x); -} - -template -struct orthogonal_t { T *x; }; - -template -void init(orthogonal_t *ptr) { - ptr->x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - ptr->x[q * i + i] = (T)1; - } -} - -template -orthogonal_t copy (orthogonal_t m) { - orthogonal_t m_copy; - m_copy.x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q * q; i++) { - m_copy.x[i] = m.x[i]; - } - - return m_copy; -} - -template -void free_spin (orthogonal_t m) { - free(m.x); -} - -template -vector_t act (orthogonal_t m, vector_t v) { - vector_t v_rot; - v_rot.x = (T *)calloc(q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - v_rot.x[i] += m.x[q * i + j] * v.x[j]; - } - } - - return v_rot; -} - - -template -orthogonal_t act (orthogonal_t m1, orthogonal_t m2) { - orthogonal_t m2_rot; - m2_rot.x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k]; - } - } - } - - return m2_rot; -} - -template -vector_t act_inverse (orthogonal_t m, vector_t v) { - vector_t v_rot; - v_rot.x = (T *)calloc(q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - v_rot.x[i] += m.x[q * j + i] * v.x[j]; - } - } - - return v_rot; -} - -template -orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) { - orthogonal_t m2_rot; - m2_rot.x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k]; - } - } - } - - return m2_rot; -} - -template -void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { - double *v = (double *)malloc(q * sizeof(double)); - double v2 = 0; - - for (q_t i = 0; i < q; i++) { - v[i] = gsl_ran_ugaussian(r); - v2 += v[i] * v[i]; - } - - ptr->x = (double *)calloc(q * q, sizeof(double)); - - for (q_t i = 0; i < q; i++) { - ptr->x[q * i + i] = 1.0; - for (q_t j = 0; j < q; j++) { - ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2; - } - } - - free(v); -} - template v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { v_t nv = 0; diff --git a/lib/orthogonal.c b/lib/orthogonal.c deleted file mode 100644 index 87569ae..0000000 --- a/lib/orthogonal.c +++ /dev/null @@ -1,99 +0,0 @@ - -#include "orthogonal.h" - -void vector_replace(q_t n, double *v1, const double *v2) { - // writes vector v2 of length n to memory located at v1 - for (q_t i = 0; i < n; i++) { - v1[i] = v2[i]; - } -} - -void vector_add(q_t n, double *v1, const double *v2) { - // adds vector v2 of length n to vector v1 - for (q_t i = 0; i < n; i++) { - v1[i] += v2[i]; - } -} - -void vector_subtract(q_t n, double *v1, const double *v2) { - // subtracts vector v2 of length n from vector v1 - for (q_t i = 0; i < n; i++) { - v1[i] -= v2[i]; - } -} - -double *vector_rotate(q_t n, const double *rot, const double *vec) { - // multiplies n by n rotation matrix rot to vector vec - double *rot_vec = (double *)malloc(n * sizeof(double)); - - double prod = 0.0; - for (q_t i = 0; i < n; i++) { - prod += rot[i] * vec[i]; - } - - for (q_t i = 0; i < n; i++) { - rot_vec[i] = vec[i] - 2 * prod * rot[i]; - } - - return rot_vec; -} - -double *vector_rotate_inverse(q_t n, const double *rot, const double *vec) { - double *rot_vec = (double *)calloc(n, sizeof(double)); - - for (q_t i = 0; i < n; i++) { - for (q_t j = 0; j < n; j++) { - rot_vec[i] += rot[n * j + i] * vec[j]; - } - } - - return rot_vec; -} - -double vector_dot(q_t n, const double *v1, const double *v2) { - double dot = 0; - - for (q_t i = 0; i < n; i++) { - dot += v1[i] * v2[i]; - } - - return dot; -} - -double *orthogonal_rotate(q_t n, const double *r, const double *m) { - double *mul = (double *)calloc(n * n, sizeof(double)); - - for (q_t i = 0; i < n; i++) { - double akOki = 0; - - for (q_t k = 0; k < n; k++) { - akOki += r[k] * m[n * k + i]; - } - - for (q_t j = 0; j < n; j++) { - mul[n * j + i] = m[n * j + i] - 2 * akOki * r[j]; - } - } - - return mul; -} - -double *gen_rot(gsl_rng *r, q_t n) { - double *v = (double *)malloc(n * sizeof(double)); - - double v2 = 0; - - for (q_t i = 0; i < n; i++) { - v[i] = gsl_ran_ugaussian(r); - v2 += v[i] * v[i]; - } - - double magv = sqrt(v2); - - for (q_t i = 0; i < n; i++) { - v[i] /= magv; - } - - return v; -} - diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 60d5f49..0b2fdd5 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -1,24 +1,122 @@ +#pragma once + #include #include #include -#include +#include #include "types.h" -void vector_replace(q_t n, double *v1, const double *v2); +template +struct orthogonal_t { T *x; }; + +template +void init(orthogonal_t *ptr) { + ptr->x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + ptr->x[q * i + i] = (T)1; + } +} + +template +orthogonal_t copy (orthogonal_t m) { + orthogonal_t m_copy; + m_copy.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q * q; i++) { + m_copy.x[i] = m.x[i]; + } + + return m_copy; +} + +template +void free_spin (orthogonal_t m) { + free(m.x); +} + +template +vector_t act (orthogonal_t m, vector_t v) { + vector_t v_rot; + v_rot.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot.x[i] += m.x[q * i + j] * v.x[j]; + } + } + + return v_rot; +} + +template +orthogonal_t act (orthogonal_t m1, orthogonal_t m2) { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k]; + } + } + } + + return m2_rot; +} + +template +vector_t act_inverse (orthogonal_t m, vector_t v) { + vector_t v_rot; + v_rot.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot.x[i] += m.x[q * j + i] * v.x[j]; + } + } + + return v_rot; +} + +template +orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); -void vector_add(q_t n, double *v1, const double *v2); + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k]; + } + } + } -void vector_subtract(q_t n, double *v1, const double *v2); + return m2_rot; +} -double *vector_rotate(q_t n, const double *rot, const double *vec); +template +void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { + double *v = (double *)malloc(q * sizeof(double)); + double v2 = 0; -double *vector_rotate_inverse(q_t n, const double *rot, const double *vec); + for (q_t i = 0; i < q; i++) { + v[i] = gsl_ran_ugaussian(r); + v2 += v[i] * v[i]; + } -double vector_dot(q_t n, const double *v1, const double *v2); + ptr->x = (double *)calloc(q * q, sizeof(double)); + + for (q_t i = 0; i < q; i++) { + ptr->x[q * i + i] = 1.0; + for (q_t j = 0; j < q; j++) { + ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2; + } + } -double *orthogonal_rotate(q_t n, const double *m1, const double *m2); + free(v); +} -double *gen_rot(gsl_rng *r, q_t n); diff --git a/lib/vector.h b/lib/vector.h new file mode 100644 index 0000000..c7f459c --- /dev/null +++ b/lib/vector.h @@ -0,0 +1,72 @@ + +#pragma once + +#include +#include + +#include "types.h" + +template +struct vector_t { T *x; }; + +template +void init(vector_t *ptr) { + ptr->x = (T *)calloc(q, sizeof(T)); + + ptr->x[0] = (T)1; +} + +template +vector_t copy (vector_t v) { + vector_t v_copy; + + v_copy.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + v_copy.x[i] = v.x[i]; + } + + 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 subtract (vector_t v1, vector_t v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] -= v2.x[i]; + } +} + +template +vector_t scalar_multiple(v_t a, vector_t v) { + vector_t multiple; + multiple.x = (T *)malloc(q * sizeof(T)); + for (q_t i = 0; i < q; i++) { + multiple.x[i] = a * v.x[i]; + } + + return multiple; +} + +template +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 +void free_spin (vector_t v) { + free(v.x); +} + -- cgit v1.2.3-54-g00ecf From 466812e61e2ccec7750c791835111b402938411c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 9 Jul 2018 14:19:16 -0400 Subject: wolff run from own function, called with types to run --- CMakeLists.txt | 8 +-- lib/orthogonal.h | 116 +++++++++++++++++++++++++++------------- lib/wolff.h | 70 +++++++++++++++++++++++++ src/wolff.cpp | 134 ----------------------------------------------- src/wolff_heisenberg.cpp | 77 +++++++++++++++++++++++++++ src/wolff_planar.cpp | 77 +++++++++++++++++++++++++++ 6 files changed, 308 insertions(+), 174 deletions(-) create mode 100644 lib/wolff.h delete mode 100644 src/wolff.cpp create mode 100644 src/wolff_heisenberg.cpp create mode 100644 src/wolff_planar.cpp (limited to 'lib') diff --git a/CMakeLists.txt b/CMakeLists.txt index 7c1bc32..ffbde47 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,7 +12,8 @@ link_directories(~/.local/lib) file(GLOB CSOURCES lib/*.c) file(GLOB CPPSOURCES lib/*.cpp) add_executable(wolff_finite src/wolff_finite.c ${CSOURCES}) -add_executable(wolff src/wolff.cpp ${CPPSOURCES} ${CSOURCES}) +add_executable(wolff_heisenberg src/wolff_heisenberg.cpp ${CPPSOURCES} ${CSOURCES}) +add_executable(wolff_planar src/wolff_planar.cpp ${CPPSOURCES} ${CSOURCES}) find_package(OpenMP) if (OPENMP_FOUND) @@ -21,7 +22,8 @@ if (OPENMP_FOUND) endif() target_link_libraries(wolff_finite gsl m cblas fftw3) -target_link_libraries(wolff gsl m cblas fftw3) +target_link_libraries(wolff_heisenberg gsl m cblas fftw3) +target_link_libraries(wolff_planar gsl m cblas fftw3) -install(TARGETS wolff_finite wolff DESTINATION bin) +install(TARGETS wolff_finite wolff_heisenberg wolff_planar DESTINATION bin) diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 0b2fdd5..0a2b5c7 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -9,10 +9,11 @@ #include "types.h" template -struct orthogonal_t { T *x; }; +struct orthogonal_t { bool is_reflection; T *x; }; template void init(orthogonal_t *ptr) { + ptr->is_reflection = false; ptr->x = (T *)calloc(q * q, sizeof(T)); for (q_t i = 0; i < q; i++) { @@ -23,9 +24,19 @@ void init(orthogonal_t *ptr) { template orthogonal_t copy (orthogonal_t m) { orthogonal_t m_copy; - m_copy.x = (T *)calloc(q * q, sizeof(T)); + m_copy.is_reflection = m.is_reflection; - for (q_t i = 0; i < q * q; i++) { + q_t size; + + if (m.is_reflection) { + size = q; + } else { + size = q * q; + } + + m_copy.x = (T *)calloc(size, sizeof(T)); + + for (q_t i = 0; i < size; i++) { m_copy.x[i] = m.x[i]; } @@ -42,9 +53,19 @@ vector_t act (orthogonal_t m, vector_t v) { vector_t v_rot; v_rot.x = (T *)calloc(q, sizeof(T)); - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - v_rot.x[i] += m.x[q * i + j] * v.x[j]; + if (m.is_reflection) { + double prod = 0; + for (q_t i = 0; i < q; i++) { + prod += v.x[i] * m.x[i]; + } + for (q_t i = 0; i < q; i++) { + v_rot.x[i] = v.x[i] - 2 * prod * m.x[i]; + } + } else { + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot.x[i] += m.x[q * i + j] * v.x[j]; + } } } @@ -54,12 +75,28 @@ vector_t act (orthogonal_t m, vector_t v) { template orthogonal_t act (orthogonal_t m1, orthogonal_t m2) { orthogonal_t m2_rot; + + m2_rot.is_reflection = false; m2_rot.x = (T *)calloc(q * q, sizeof(T)); - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { + if (m1.is_reflection) { + for (q_t i = 0; i < q; i++) { + double akOki = 0; + for (q_t k = 0; k < q; k++) { - m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k]; + akOki += m1.x[k] * m2.x[q * k + i]; + } + + for (q_t j = 0; j < q; j++) { + m2_rot.x[q * j + i] = m2.x[q * j + i] - 2 * akOki * m1.x[j]; + } + } + } else { + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k]; + } } } } @@ -69,54 +106,59 @@ orthogonal_t act (orthogonal_t m1, orthogonal_t m2) { template vector_t act_inverse (orthogonal_t m, vector_t v) { - vector_t v_rot; - v_rot.x = (T *)calloc(q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - v_rot.x[i] += m.x[q * j + i] * v.x[j]; + if (m.is_reflection) { + return act(m, v); // reflections are their own inverse + } else { + vector_t v_rot; + v_rot.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + v_rot.x[i] += m.x[q * j + i] * v.x[j]; + } } - } - return v_rot; + return v_rot; + } } template orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) { - orthogonal_t m2_rot; - m2_rot.x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k]; + if (m1.is_reflection) { + return act(m1, m2); // reflections are their own inverse + } else { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k]; + } } } - } - return m2_rot; + return m2_rot; + } } template void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { - double *v = (double *)malloc(q * sizeof(double)); + ptr->is_reflection = true; + ptr->x = (double *)calloc(q, sizeof(double)); + double v2 = 0; for (q_t i = 0; i < q; i++) { - v[i] = gsl_ran_ugaussian(r); - v2 += v[i] * v[i]; + ptr->x[i] = gsl_ran_ugaussian(r); + v2 += ptr->x[i] * ptr->x[i]; } - ptr->x = (double *)calloc(q * q, sizeof(double)); - + double mag_v = sqrt(v2); + for (q_t i = 0; i < q; i++) { - ptr->x[q * i + i] = 1.0; - for (q_t j = 0; j < q; j++) { - ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2; - } + ptr->x[i] /= mag_v; } - - free(v); } diff --git a/lib/wolff.h b/lib/wolff.h new file mode 100644 index 0000000..81830ee --- /dev/null +++ b/lib/wolff.h @@ -0,0 +1,70 @@ + +#include +#include + +#include + +template +double H_vector(vector_t v1, T *H) { + vector_t H_vec; + H_vec.x = H; + return (double)(dot (v1, H_vec)); +} + +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) { + + 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()); + + char *filename_M = (char *)malloc(255 * sizeof(char)); + char *filename_E = (char *)malloc(255 * sizeof(char)); + char *filename_S = (char *)malloc(255 * sizeof(char)); + + sprintf(filename_M, "wolff_%lu_M.dat", timestamp); + sprintf(filename_E, "wolff_%lu_E.dat", timestamp); + sprintf(filename_S, "wolff_%lu_S.dat", timestamp); + + FILE *outfile_M = fopen(filename_M, "wb"); + FILE *outfile_E = fopen(filename_E, "wb"); + FILE *outfile_S = fopen(filename_S, "wb"); + + free(filename_M); + free(filename_E); + free(filename_S); + + v_t cluster_size = 0; + + if (!silent) printf("\n"); + for (count_t steps = 0; steps < N; steps++) { + if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, M_0 = %.2f, S = %" PRIv "\n", steps, N, s.E, s.M.x[0], cluster_size); + + v_t v0 = gsl_rng_uniform_int(r, s.nv); + + R_t step; + generate_rotation(r, &step); + + cluster_size = flip_cluster (&s, v0, step, r); + + free_spin(step); + + fwrite(&(s.E), sizeof(double), 1, outfile_E); + fwrite(s.M.x, sizeof(double), 2, outfile_M); + fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S); + + } + if (!silent) { + printf("\033[F\033[J"); + } + printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, M_0 = %.2f, S = %" PRIv "\n", N, N, s.E, s.M.x[0], cluster_size); + + fclose(outfile_M); + fclose(outfile_E); + fclose(outfile_S); + + gsl_rng_free(r); +} + diff --git a/src/wolff.cpp b/src/wolff.cpp deleted file mode 100644 index f685129..0000000 --- a/src/wolff.cpp +++ /dev/null @@ -1,134 +0,0 @@ - -#include -#include - -#include - -double H_vector(vector_t <2, double> v1, double *H) { - vector_t <2, double> H_vec; - H_vec.x = H; - return dot <2, double> (v1, H_vec); -} - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e7; - - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double *H = (double *)calloc(MAX_Q, sizeof(double)); - - bool silent = false; - - int opt; - q_t J_ind = 0; - q_t H_ind = 0; - - while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:s")) != -1) { - switch (opt) { - case 'N': // number of steps - N = (count_t)atof(optarg); - break; - case 'D': // dimension - D = atoi(optarg); - break; - case 'L': // linear size - L = atoi(optarg); - break; - case 'T': // temperature - T = atof(optarg); - break; - case 'H': // external field. nth call couples to state n - H[H_ind] = atof(optarg); - H_ind++; - break; - case 's': // don't print anything during simulation. speeds up slightly - silent = true; - break; - default: - exit(EXIT_FAILURE); - } - } - - state_t , vector_t <2, double>> s(D, L, T, dot <2, double>, std::bind(H_vector, std::placeholders::_1, H)); - - // initialize random number generator - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - - unsigned long timestamp; - - { - struct timespec spec; - clock_gettime(CLOCK_REALTIME, &spec); - timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; - } - - FILE *outfile_info = fopen("wolff_metadata.txt", "a"); - - fprintf(outfile_info, "<| \"ID\" -> %lu, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, D, L, s.nv, s.ne, T); - - for (q_t i = 0; i < 2; i++) { - fprintf(outfile_info, "%.15f", H[i]); - if (i < 2 - 1) { - fprintf(outfile_info, ", "); - } - } - - fprintf(outfile_info, "} |>\n"); - - fclose(outfile_info); - - char *filename_M = (char *)malloc(255 * sizeof(char)); - char *filename_E = (char *)malloc(255 * sizeof(char)); - char *filename_S = (char *)malloc(255 * sizeof(char)); - - sprintf(filename_M, "wolff_%lu_M.dat", timestamp); - sprintf(filename_E, "wolff_%lu_E.dat", timestamp); - sprintf(filename_S, "wolff_%lu_S.dat", timestamp); - - FILE *outfile_M = fopen(filename_M, "wb"); - FILE *outfile_E = fopen(filename_E, "wb"); - FILE *outfile_S = fopen(filename_S, "wb"); - - free(filename_M); - free(filename_E); - free(filename_S); - - v_t cluster_size = 0; - - if (!silent) printf("\n"); - for (count_t steps = 0; steps < N; steps++) { - if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, M_0 = %.2f, S = %" PRIv "\n", steps, N, s.E, s.M.x[0], cluster_size); - - v_t v0 = gsl_rng_uniform_int(r, s.nv); - - orthogonal_t <2, double> step; - generate_rotation<2>(r, &step); - - cluster_size = flip_cluster , vector_t <2, double>> (&s, v0, step, r); - - free_spin(step); - - fwrite(&(s.E), sizeof(double), 1, outfile_E); - fwrite(s.M.x, sizeof(double), 2, outfile_M); - fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S); - - } - if (!silent) { - printf("\033[F\033[J"); - } - printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, M_0 = %.2f, S = %" PRIv "\n", N, N, s.E, s.M.x[0], cluster_size); - - fclose(outfile_M); - fclose(outfile_E); - fclose(outfile_S); - - gsl_rng_free(r); - - free(H); - - return 0; -} - diff --git a/src/wolff_heisenberg.cpp b/src/wolff_heisenberg.cpp new file mode 100644 index 0000000..d1ebd48 --- /dev/null +++ b/src/wolff_heisenberg.cpp @@ -0,0 +1,77 @@ + +#include + +#include + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e7; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double *H = (double *)calloc(MAX_Q, sizeof(double)); + + bool silent = false; + + int opt; + q_t J_ind = 0; + q_t H_ind = 0; + + while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:s")) != -1) { + switch (opt) { + case 'N': // number of steps + N = (count_t)atof(optarg); + break; + case 'D': // dimension + D = atoi(optarg); + break; + case 'L': // linear size + L = atoi(optarg); + break; + case 'T': // temperature + T = atof(optarg); + break; + case 'H': // external field. nth call couples to state n + H[H_ind] = atof(optarg); + H_ind++; + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); + + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"HEISENBERG\", q -> \"3\", \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, D, L, L * L, D * L * L, T); + + for (q_t i = 0; i < 2; i++) { + fprintf(outfile_info, "%.15f", H[i]); + if (i < 2 - 1) { + fprintf(outfile_info, ", "); + } + } + + fprintf(outfile_info, "} |>\n"); + + fclose(outfile_info); + + + wolff , vector_t <3, double>> (N, D, L, T, dot <3, double>, std::bind(H_vector <3, double>, std::placeholders::_1, H), timestamp, silent); + + free(H); + + return 0; +} + diff --git a/src/wolff_planar.cpp b/src/wolff_planar.cpp new file mode 100644 index 0000000..02ededc --- /dev/null +++ b/src/wolff_planar.cpp @@ -0,0 +1,77 @@ + +#include + +#include + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e7; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double *H = (double *)calloc(MAX_Q, sizeof(double)); + + bool silent = false; + + int opt; + q_t J_ind = 0; + q_t H_ind = 0; + + while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:s")) != -1) { + switch (opt) { + case 'N': // number of steps + N = (count_t)atof(optarg); + break; + case 'D': // dimension + D = atoi(optarg); + break; + case 'L': // linear size + L = atoi(optarg); + break; + case 'T': // temperature + T = atof(optarg); + break; + case 'H': // external field. nth call couples to state n + H[H_ind] = atof(optarg); + H_ind++; + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); + + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"PLANAR\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, D, L, L * L, D * L * L, T); + + for (q_t i = 0; i < 2; i++) { + fprintf(outfile_info, "%.15f", H[i]); + if (i < 2 - 1) { + fprintf(outfile_info, ", "); + } + } + + fprintf(outfile_info, "} |>\n"); + + fclose(outfile_info); + + + wolff , vector_t <2, double>> (N, D, L, T, dot <2, double>, std::bind(H_vector <2, double>, std::placeholders::_1, H), timestamp, silent); + + free(H); + + return 0; +} + -- cgit v1.2.3-54-g00ecf From 88086fc222e1802d2a68edf5f1097cf9685ec42c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 9 Jul 2018 14:29:16 -0400 Subject: removed some references to the specific form of the magnetization for vector models from wolff.h --- lib/orthogonal.h | 12 ++++++++++++ lib/wolff.h | 11 +++++++---- 2 files changed, 19 insertions(+), 4 deletions(-) (limited to 'lib') diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 0a2b5c7..2d0e1a1 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -161,4 +161,16 @@ void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { } } +template +void write_magnetization(vector_t M, FILE *outfile) { + fwrite(M.x, sizeof(double), q, outfile); +} + +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++) { + float M_tmp = (float)M.x[i]; + fwrite(&M_tmp, sizeof(float), 1, outfile); + } +} diff --git a/lib/wolff.h b/lib/wolff.h index 81830ee..caf413b 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -40,7 +40,7 @@ void wolff(count_t N, D_t D, L_t L, double T, std::function J if (!silent) printf("\n"); for (count_t steps = 0; steps < N; steps++) { - if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, M_0 = %.2f, S = %" PRIv "\n", steps, N, s.E, s.M.x[0], cluster_size); + if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, cluster_size); v_t v0 = gsl_rng_uniform_int(r, s.nv); @@ -51,15 +51,18 @@ void wolff(count_t N, D_t D, L_t L, double T, std::function J free_spin(step); - fwrite(&(s.E), sizeof(double), 1, outfile_E); - fwrite(s.M.x, sizeof(double), 2, outfile_M); + { + float smaller_E = (float)s.E; + fwrite(&smaller_E, sizeof(float), 1, outfile_E); + } + write_magnetization(s.M, outfile_M); fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S); } if (!silent) { printf("\033[F\033[J"); } - printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, M_0 = %.2f, S = %" PRIv "\n", N, N, s.E, s.M.x[0], cluster_size); + printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, cluster_size); fclose(outfile_M); fclose(outfile_E); -- cgit v1.2.3-54-g00ecf From 1810103bc9ac4c9a8d432d113f5ca6eae6560fb4 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 10 Jul 2018 11:17:03 -0400 Subject: minor reorganizing --- lib/orthogonal.h | 13 ------------- lib/vector.h | 13 +++++++++++++ 2 files changed, 13 insertions(+), 13 deletions(-) (limited to 'lib') diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 2d0e1a1..85f7a20 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -161,16 +161,3 @@ void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { } } -template -void write_magnetization(vector_t M, FILE *outfile) { - fwrite(M.x, sizeof(double), q, outfile); -} - -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++) { - float M_tmp = (float)M.x[i]; - fwrite(&M_tmp, sizeof(float), 1, outfile); - } -} - diff --git a/lib/vector.h b/lib/vector.h index c7f459c..62ce59e 100644 --- a/lib/vector.h +++ b/lib/vector.h @@ -70,3 +70,16 @@ void free_spin (vector_t v) { free(v.x); } +template +void write_magnetization(vector_t M, FILE *outfile) { + fwrite(M.x, sizeof(double), q, outfile); +} + +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++) { + float M_tmp = (float)M.x[i]; + fwrite(&M_tmp, sizeof(float), 1, outfile); + } +} + -- cgit v1.2.3-54-g00ecf