From 8a6109a78a6f7b9e9585cd89e8a10d1f626b3af1 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 17 Jul 2018 22:37:57 -0400 Subject: made many changes to implement alternate reflection generators, this is broken --- lib/cluster.h | 60 ++++++-------------------------------------- lib/orthogonal.h | 54 ++++++++++++++++++++++++++++++++++----- lib/state.h | 53 ++++++++++++++++++++++++++++++++++++++ lib/wolff.h | 8 +++--- src/analyze_correlations.cpp | 6 ++++- src/wolff_heisenberg.cpp | 26 ++++++++++++++++--- src/wolff_planar.cpp | 27 +++++++++++++++++--- 7 files changed, 164 insertions(+), 70 deletions(-) create mode 100644 lib/state.h diff --git a/lib/cluster.h b/lib/cluster.h index a20912e..ccd8725 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -13,6 +13,7 @@ #include #include +#include "state.h" #include "types.h" #include "rand.h" #include "stack.h" @@ -50,52 +51,6 @@ 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 v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { v_t nv = 0; @@ -103,7 +58,7 @@ v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { 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)); + bool *marks = (bool *)calloc(state->nv, sizeof(bool)); while (stack != NULL) { v_t v = stack_pop(&stack); @@ -112,14 +67,13 @@ v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { 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) { + if (v == state->nv) { + R_old = state->R; R_new = act (r, R_old); } else { + si_old = state->spins[v]; si_new = act (r, si_old); } @@ -130,13 +84,13 @@ v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { X_t sj; - if (vn != state->g->nv - 1) { + if (vn != state->nv) { sj = state->spins[vn]; } double prob; - bool is_ext = (v == state->g->nv - 1 || vn == state->g->nv - 1); + bool is_ext = (v == state->nv || vn == state->nv); if (is_ext) { X_t rs_old, rs_new; diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 85f7a20..6166223 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -6,6 +6,7 @@ #include #include +#include "state.h" #include "types.h" template @@ -143,21 +144,62 @@ orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) } template -void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { - ptr->is_reflection = true; - ptr->x = (double *)calloc(q, sizeof(double)); +orthogonal_t generate_rotation_uniform (gsl_rng *r, const state_t , vector_t > *s) { + orthogonal_t ptr; + ptr.is_reflection = true; + ptr.x = (double *)calloc(q, sizeof(double)); double v2 = 0; for (q_t i = 0; i < q; i++) { - ptr->x[i] = gsl_ran_ugaussian(r); - v2 += ptr->x[i] * ptr->x[i]; + ptr.x[i] = gsl_ran_ugaussian(r); + v2 += ptr.x[i] * ptr.x[i]; } double mag_v = sqrt(v2); for (q_t i = 0; i < q; i++) { - ptr->x[i] /= mag_v; + ptr.x[i] /= mag_v; } + + return ptr; +} + +template +orthogonal_t generate_rotation_perturbation (gsl_rng *r, const state_t , vector_t > *s, double epsilon) { + orthogonal_t ptr; + vector_t tmp_v; + ptr.is_reflection = true; + + tmp_v.x = (double *)malloc(q * sizeof(double)); + + double tmp2 = 0; + double M2 = 0; + double tmpM = 0; + + for (q_t i = 0; i < q; i++) { + tmp_v.x[i] = gsl_ran_ugaussian(r); + M2 += pow(s->M.x[i], 2); + tmpM += tmp_v.x[i] * s->M.x[i]; + } + + double v2 = 0; + + for (q_t i = 0; i < q; i++) { + tmp_v.x[i] = (tmp_v.x[i] - tmpM * s->M.x[i] / M2) + epsilon * gsl_ran_ugaussian(r); + v2 += tmp_v.x[i] * tmp_v.x[i]; + } + + double mag_v = sqrt(v2); + + for (q_t i = 0; i < q; i++) { + tmp_v.x[i] /= mag_v; + } + + vector_t v = act(s->R, tmp_v); + free(tmp_v.x); + ptr.x = v.x; + + return ptr; } diff --git a/lib/state.h b/lib/state.h new file mode 100644 index 0000000..8d4265b --- /dev/null +++ b/lib/state.h @@ -0,0 +1,53 @@ + +#pragma once + +#include + +#include "types.h" +#include "graph.h" + +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); + } +}; + + + diff --git a/lib/wolff.h b/lib/wolff.h index caf413b..90470d5 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -2,7 +2,8 @@ #include #include -#include +#include "cluster.h" +#include "state.h" template double H_vector(vector_t v1, T *H) { @@ -12,7 +13,7 @@ double H_vector(vector_t v1, T *H) { } 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) { +void wolff(count_t N, D_t D, L_t L, double T, std::function J, std::function H, std::function *)> gen_R, unsigned long timestamp, bool silent) { state_t s(D, L, T, J, H); @@ -44,8 +45,7 @@ void wolff(count_t N, D_t D, L_t L, double T, std::function J v_t v0 = gsl_rng_uniform_int(r, s.nv); - R_t step; - generate_rotation(r, &step); + R_t step = gen_R(r, &s); cluster_size = flip_cluster (&s, v0, step, r); diff --git a/src/analyze_correlations.cpp b/src/analyze_correlations.cpp index 6b8d28d..cf8e740 100644 --- a/src/analyze_correlations.cpp +++ b/src/analyze_correlations.cpp @@ -239,7 +239,11 @@ int main (int argc, char *argv[]) { for (q_t i = 0; i < q - 1; i++) { fscanf(metadata, "%lf, ", &(H[i])); } - fscanf(metadata, "%lf} |>\n", &(H[q - 1])); + char *field = (char *)malloc(32 * sizeof(char)); + double epsilon; + fscanf(metadata, "%lf}, \"GENERATOR\" -> \"%[^\"]\", \"EPS\" -> %lf |>\n", &(H[q - 1]), field, &epsilon); + + free(field); free(H); char *filename_E = (char *)malloc(128 * sizeof(char)); diff --git a/src/wolff_heisenberg.cpp b/src/wolff_heisenberg.cpp index aed4156..83ea503 100644 --- a/src/wolff_heisenberg.cpp +++ b/src/wolff_heisenberg.cpp @@ -13,12 +13,14 @@ int main(int argc, char *argv[]) { double *H = (double *)calloc(MAX_Q, sizeof(double)); bool silent = false; + bool use_pert = false; int opt; q_t J_ind = 0; q_t H_ind = 0; + double epsilon = 1; - while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:s")) != -1) { + while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:")) != -1) { switch (opt) { case 'N': // number of steps N = (count_t)atof(optarg); @@ -39,6 +41,12 @@ int main(int argc, char *argv[]) { case 's': // don't print anything during simulation. speeds up slightly silent = true; break; + case 'p': + use_pert = true; + break; + case 'e': + epsilon = atof(optarg); + break; default: exit(EXIT_FAILURE); } @@ -52,6 +60,18 @@ int main(int argc, char *argv[]) { timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; } + std::function (gsl_rng *, const state_t , vector_t <3, double>> *)> gen_R; + + const char *pert_type; + + if (use_pert) { + gen_R = std::bind(generate_rotation_perturbation <3>, std::placeholders::_1, std::placeholders::_2, epsilon); + pert_type = "PERTURB"; + } else { + gen_R = generate_rotation_uniform <3>; + pert_type = "UNIFORM"; + } + 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); @@ -63,12 +83,12 @@ int main(int argc, char *argv[]) { } } - fprintf(outfile_info, "} |>\n"); + fprintf(outfile_info, "}, \"GENERATOR\" -> \"%s\", \"EPS\" -> %g |>\n", pert_type, epsilon); 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); + wolff , vector_t <3, double>> (N, D, L, T, dot <3, double>, std::bind(H_vector <3, double>, std::placeholders::_1, H), gen_R, timestamp, silent); free(H); diff --git a/src/wolff_planar.cpp b/src/wolff_planar.cpp index 02ededc..2a3d02b 100644 --- a/src/wolff_planar.cpp +++ b/src/wolff_planar.cpp @@ -13,12 +13,14 @@ int main(int argc, char *argv[]) { double *H = (double *)calloc(MAX_Q, sizeof(double)); bool silent = false; + bool use_pert = false; int opt; q_t J_ind = 0; q_t H_ind = 0; + double epsilon = 1; - while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:s")) != -1) { + while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:")) != -1) { switch (opt) { case 'N': // number of steps N = (count_t)atof(optarg); @@ -39,6 +41,12 @@ int main(int argc, char *argv[]) { case 's': // don't print anything during simulation. speeds up slightly silent = true; break; + case 'p': + use_pert = true; + break; + case 'e': + epsilon = atof(optarg); + break; default: exit(EXIT_FAILURE); } @@ -52,6 +60,19 @@ int main(int argc, char *argv[]) { timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; } + const char *pert_type; + + std::function (gsl_rng *, const state_t , vector_t <2, double>> *)> gen_R; + + if (use_pert) { + gen_R = std::bind(generate_rotation_perturbation <2>, std::placeholders::_1, std::placeholders::_2, epsilon); + pert_type = "PERTURB"; + } else { + gen_R = generate_rotation_uniform <2>; + pert_type = "UNIFORM"; + } + + 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); @@ -63,12 +84,12 @@ int main(int argc, char *argv[]) { } } - fprintf(outfile_info, "} |>\n"); + fprintf(outfile_info, "}, \"GENERATOR\" -> \"%s\", \"EPS\" -> %g |>\n", pert_type, epsilon); 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); + wolff , vector_t <2, double>> (N, D, L, T, dot <2, double>, std::bind(H_vector <2, double>, std::placeholders::_1, H), gen_R, timestamp, silent); free(H); -- cgit v1.2.3-70-g09d2