From 802b63ddf121b520db7942fe330cce6004fbeb6d Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 25 Jul 2018 16:22:50 -0400 Subject: got everyone recording data, and fixed huge bug in the updating of ReF and ImF --- lib/cluster.h | 8 ++--- lib/orthogonal.h | 1 + lib/potts.h | 47 +++++++++++++-------------- src/wolff_On.cpp | 4 +-- src/wolff_clock.cpp | 4 +-- src/wolff_ising.cpp | 91 ++++++++++++++++++++++++++++++++++++++------------- src/wolff_potts.cpp | 93 ++++++++++++++++++++++++++++++++++++++++------------- 7 files changed, 169 insertions(+), 79 deletions(-) diff --git a/lib/cluster.h b/lib/cluster.h index b8c98e5..3261969 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -77,11 +77,11 @@ void flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { for (D_t i = 0; i < state->D; i++) { L_t x = (non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L; - add(&(state->ReF[i]), -state->precomputed_cos[i], rs_old); - add(&(state->ReF[i]), state->precomputed_cos[i], rs_new); + add(&(state->ReF[i]), -state->precomputed_cos[x], rs_old); + add(&(state->ReF[i]), state->precomputed_cos[x], rs_new); - add(&(state->ImF[i]), -state->precomputed_sin[i], rs_old); - add(&(state->ImF[i]), state->precomputed_sin[i], rs_new); + add(&(state->ImF[i]), -state->precomputed_sin[x], rs_old); + add(&(state->ImF[i]), state->precomputed_sin[x], rs_new); } free_spin (rs_old); diff --git a/lib/orthogonal.h b/lib/orthogonal.h index ce2d300..0ba5eee 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -8,6 +8,7 @@ #include "state.h" #include "types.h" +#include "vector.h" template struct orthogonal_t { bool is_reflection; T *x; }; diff --git a/lib/potts.h b/lib/potts.h index 732a76f..e61e4e1 100644 --- a/lib/potts.h +++ b/lib/potts.h @@ -4,6 +4,7 @@ #include #include "types.h" +#include "vector.h" /* The following is the minimum definition of a spin class. * @@ -30,8 +31,8 @@ class potts_t { public: q_t x; - typedef int *M_t; - typedef double *F_t; + typedef vector_t M_t; + typedef vector_t F_t; }; template @@ -44,56 +45,54 @@ void free_spin(potts_t s) { // do nothing! } -void free_spin(int *s) { - free(s); -} - -void free_spin(double *s) { - free(s); -} - template potts_t copy(potts_t s) { return s; } template -void add(typename potts_t::M_t *s1, int a, potts_t s2) { - (*s1)[s2.x] += a; +void add(vector_t *s1, int a, potts_t s2) { + (s1->x)[s2.x] += a; } template -void add(typename potts_t::F_t *s1, double a, potts_t s2) { - (*s1)[s2.x] += a; +void add(vector_t *s1, double a, potts_t s2) { + (s1->x)[s2.x] += a; } template -typename potts_t::M_t scalar_multiple(int factor, potts_t s) { - int *M = (int *)calloc(q, sizeof(int)); - M[s.x] += factor; +vector_t scalar_multiple(int factor, potts_t s) { + vector_t M; + M.x = (int *)calloc(q, sizeof(int)); + M.x[s.x] += factor; return M; } template -typename potts_t::F_t scalar_multiple(double factor, potts_t s) { - double *F = (double *)calloc(q, sizeof(double)); - F[s.x] += factor; +vector_t scalar_multiple(double factor, potts_t s) { + vector_t F; + F.x = (double *)calloc(q, sizeof(double)); + F.x[s.x] += factor; return F; } +// we could inherit norm_squared from vector.h, but convention dictates that +// potts norms be changed by a constant factor template -double norm_squared(typename potts_t::F_t s) { +double norm_squared(vector_t s) { double total = 0; for (q_t i = 0; i < q; i++) { - total += pow(s[i], 2); + total += pow(s.x[i], 2); } return total * (double)q / ((double)q - 1.0); } +// we could inherit write_magnetization from vector.h, but since M.x must sum +// to nv we don't need to write the last element template -void write_magnetization(typename potts_t::M_t M, FILE *outfile) { - fwrite(&M, sizeof(int), q, outfile); +void write_magnetization(vector_t M, FILE *outfile) { + fwrite(M.x, sizeof(int), q - 1, outfile); } // knock yourself out diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index 997ec09..a277e0f 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -6,11 +6,10 @@ #include #endif -#include #include +#include #include -#include #include #include #include @@ -235,7 +234,6 @@ int main(int argc, char *argv[]) { wolff (N, &s, gen_R, measurements, r, silent); } - measure_free_files(measurement_flags, outfiles); free(H_vec); gsl_rng_free(r); diff --git a/src/wolff_clock.cpp b/src/wolff_clock.cpp index e186c44..86badfe 100644 --- a/src/wolff_clock.cpp +++ b/src/wolff_clock.cpp @@ -107,7 +107,7 @@ int main(int argc, char *argv[]) { if (!draw) { // a very simple example: measure the average magnetization measurement = [&] (const sim_t *s) { - average_M += (double)s->M[0] / (double)N / (double)s->nv; + average_M += (double)s->M.x[0] / (double)N / (double)s->nv; }; } else { // a more complex example: measure the average magnetization, and draw the spin configuration to the screen @@ -124,7 +124,7 @@ int main(int argc, char *argv[]) { gluOrtho2D(0.0, L, 0.0, L); measurement = [&] (const sim_t *s) { - average_M += (double)s->M[0] / (double)N / (double)s->nv; + average_M += (double)s->M.x[0] / (double)N / (double)s->nv; glClear(GL_COLOR_BUFFER_BIT); for (v_t i = 0; i < pow(L, 2); i++) { potts_t tmp_s = act_inverse(s->R, s->spins[i]); diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp index 7492ebf..5e44cab 100644 --- a/src/wolff_ising.cpp +++ b/src/wolff_ising.cpp @@ -1,5 +1,8 @@ #include +#include + +// if you have GLUT installed, you can see graphics! #ifdef HAVE_GLUT #include #endif @@ -8,10 +11,17 @@ #include #include +// finite_states.h can be included for spin types that have special variables +// defined, and it causes wolff execution to use precomputed bond probabilities #include -// include wolff.h +// rand.h uses a unix-specific way to seed the random number generator #include + +// measure.h contains useful functions for saving timeseries to files +#include + +// include wolff.h #include int main(int argc, char *argv[]) { @@ -25,11 +35,15 @@ int main(int argc, char *argv[]) { bool silent = false; bool draw = false; + bool N_is_sweeps = false; unsigned int window_size = 512; + // don't measure anything by default + unsigned char measurement_flags = 0; + int opt; - while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:")) != -1) { + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -1) { switch (opt) { case 'N': // number of steps N = (count_t)atof(optarg); @@ -49,6 +63,9 @@ int main(int argc, char *argv[]) { case 's': // don't print anything during simulation. speeds up slightly silent = true; break; + case 'S': + N_is_sweeps = true; + break; case 'd': #ifdef HAVE_GLUT draw = true; @@ -60,11 +77,23 @@ int main(int argc, char *argv[]) { case 'w': window_size = atoi(optarg); break; + case 'M': + measurement_flags ^= 1 << atoi(optarg); + break; default: exit(EXIT_FAILURE); } } + // get nanosecond timestamp for unique run id + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + // initialize random number generator gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); gsl_rng_set(r, rand_seed()); @@ -97,18 +126,16 @@ int main(int argc, char *argv[]) { return rot; }; - // define function that updates any number of measurements - std::function *)> measurement; + FILE **outfiles = measure_setup_files(measurement_flags, timestamp); - double average_M = 0; - if (!draw) { - // a very simple example: measure the average magnetization - measurement = [&] (const state_t *s) { - average_M += (double)s->M / (double)N / (double)s->nv; - }; - } else { - // a more complex example: measure the average magnetization, and draw the spin configuration to the screen + std::function *)> other_f; + uint64_t sum_of_clusterSize = 0; + if (N_is_sweeps) { + other_f = [&] (const state_t *s) { + sum_of_clusterSize += s->last_cluster_size; + }; + } else if (draw) { #ifdef HAVE_GLUT // initialize glut glutInit(&argc, argv); @@ -120,34 +147,52 @@ int main(int argc, char *argv[]) { glLoadIdentity(); gluOrtho2D(0.0, L, 0.0, L); - measurement = [&] (const state_t *s) { - average_M += (double)s->M / (double)N / (double)s->nv; + other_f = [] (const state_t *s) { glClear(GL_COLOR_BUFFER_BIT); - for (v_t i = 0; i < pow(L, 2); i++) { + for (v_t i = 0; i < pow(s->L, 2); i++) { if (s->spins[i].x == s->R.x) { glColor3f(0.0, 0.0, 0.0); } else { glColor3f(1.0, 1.0, 1.0); } - glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); + glRecti(i / s->L, i % s->L, (i / s->L) + 1, (i % s->L) + 1); } glFlush(); }; #endif + } else { + other_f = [] (const state_t *s) {}; } - // run wolff for N cluster flips - wolff(N, &s, gen_R, measurement, r, silent); + std::function *)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f); - // tell us what we found! - printf("%" PRIcount " Ising runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, = %g\n", N, D, L, T, H, average_M); + // add line to metadata file with run info + { + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); - // free the random number generator - gsl_rng_free(r); + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"ISING\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> %.15f |>\n", timestamp, s.D, s.L, s.nv, s.ne, T, H); - if (draw) { + fclose(outfile_info); } + // run wolff for N cluster flips + if (N_is_sweeps) { + count_t N_rounds = 0; + printf("\n"); + while (sum_of_clusterSize < N * s.nv) { + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size); + wolff(N, &s, gen_R, measurements, r, silent); + N_rounds++; + } + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size); + } else { + wolff(N, &s, gen_R, measurements, r, silent); + } + + // free the random number generator + gsl_rng_free(r); + measure_free_files(measurement_flags, outfiles); + return 0; } diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp index f8f1523..cdc4c07 100644 --- a/src/wolff_potts.cpp +++ b/src/wolff_potts.cpp @@ -1,5 +1,6 @@ #include +#include #ifdef HAVE_GLUT #include @@ -8,13 +9,14 @@ // include your group and spin space #include #include -#include // hack to speed things up considerably #define N_STATES POTTSQ #include // include wolff.h +#include +#include #include #include @@ -31,12 +33,16 @@ int main(int argc, char *argv[]) { bool silent = false; bool draw = false; + bool N_is_sweeps = false; unsigned int window_size = 512; + // don't measure anything by default + unsigned char measurement_flags = 0; + int opt; q_t H_ind = 0; - while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:")) != -1) { + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -1) { switch (opt) { case 'N': // number of steps N = (count_t)atof(optarg); @@ -57,6 +63,9 @@ int main(int argc, char *argv[]) { case 's': // don't print anything during simulation. speeds up slightly silent = true; break; + case 'S': + N_is_sweeps = true; + break; case 'd': #ifdef HAVE_GLUT draw = true; @@ -68,11 +77,23 @@ int main(int argc, char *argv[]) { case 'w': window_size = atoi(optarg); break; + case 'M': + measurement_flags ^= 1 << atoi(optarg); + break; default: exit(EXIT_FAILURE); } } + // get nanosecond timestamp for unique run id + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + // initialize random number generator gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); gsl_rng_set(r, rand_seed()); @@ -113,18 +134,16 @@ int main(int argc, char *argv[]) { return rot; }; - // define function that updates any number of measurements - std::function measurement; + FILE **outfiles = measure_setup_files(measurement_flags, timestamp); - double average_M = 0; - if (!draw) { - // a very simple example: measure the average magnetization - measurement = [&] (const sim_t *s) { - average_M += (double)s->M[0] / (double)N / (double)s->nv; - }; - } else { - // a more complex example: measure the average magnetization, and draw the spin configuration to the screen + std::function other_f; + uint64_t sum_of_clusterSize = 0; + if (N_is_sweeps) { + other_f = [&] (const sim_t *s) { + sum_of_clusterSize += s->last_cluster_size; + }; + } else if (draw) { #ifdef HAVE_GLUT // initialize glut glutInit(&argc, argv); @@ -136,31 +155,59 @@ int main(int argc, char *argv[]) { glLoadIdentity(); gluOrtho2D(0.0, L, 0.0, L); - measurement = [&] (const sim_t *s) { - average_M += (double)s->M[0] / (double)N / (double)s->nv; + other_f = [] (const sim_t *s) { glClear(GL_COLOR_BUFFER_BIT); - for (v_t i = 0; i < pow(L, 2); i++) { + for (v_t i = 0; i < pow(s->L, 2); i++) { potts_t tmp_s = act_inverse(s->R, s->spins[i]); glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ)); - glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); + glRecti(i / s->L, i % s->L, (i / s->L) + 1, (i % s->L) + 1); } glFlush(); }; #endif + } else { + other_f = [] (const sim_t *s) {}; } - // run wolff for N cluster flips - wolff(N, &s, gen_R, measurement, r, silent); + std::function measurements = measure_function_write_files(measurement_flags, outfiles, other_f); - // tell us what we found! - printf("%" PRIcount " %d-Potts runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, = %g\n", N, POTTSQ, D, L, T, H_vec[0], average_M); + // add line to metadata file with run info + { + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); - // free the random number generator - gsl_rng_free(r); + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"POTTS\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, POTTSQ, s.D, s.L, s.nv, s.ne, T); + + for (q_t i = 0; i < POTTSQ; i++) { + fprintf(outfile_info, "%.15f", H_vec[i]); + if (i < POTTSQ - 1) { + fprintf(outfile_info, ", "); + } + } - if (draw) { + fprintf(outfile_info, "} |>\n"); + + fclose(outfile_info); } + // run wolff for N cluster flips + if (N_is_sweeps) { + count_t N_rounds = 0; + printf("\n"); + while (sum_of_clusterSize < N * s.nv) { + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size); + wolff(N, &s, gen_R, measurements, r, silent); + N_rounds++; + } + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size); + } else { + wolff(N, &s, gen_R, measurements, r, silent); + } + + // free the random number generator + gsl_rng_free(r); + free(H_vec); + measure_free_files(measurement_flags, outfiles); + return 0; } -- cgit v1.2.3-70-g09d2