From 1fbcb4dd4e52daeb53becba33827f8e48c5606b2 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 28 Feb 2018 20:33:41 -0500 Subject: fixed major mistake in process, also got n-component version fully working --- CMakeLists.txt | 11 +- lib/dihedral.c | 28 +++++ lib/dihedral.h | 17 +++ lib/orthogonal.c | 95 +++++++++++++++ lib/orthogonal.h | 24 ++++ lib/wolff.h | 17 +++ lib/wolff_tools.c | 137 ++++++++++++++++++--- src/wolff.c | 340 --------------------------------------------------- src/wolff_potts.c | 349 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/wolff_vector.c | 202 +++++++++++++++++++++++++++++++ 10 files changed, 858 insertions(+), 362 deletions(-) create mode 100644 lib/dihedral.c create mode 100644 lib/dihedral.h create mode 100644 lib/orthogonal.c create mode 100644 lib/orthogonal.h delete mode 100644 src/wolff.c create mode 100644 src/wolff_potts.c create mode 100644 src/wolff_vector.c diff --git a/CMakeLists.txt b/CMakeLists.txt index eb86ba6..3b0de87 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,11 +2,14 @@ cmake_minimum_required(VERSION 3.0) project(wolff) +set (CMAKE_CXX_STANDARD 11) + include_directories(lib ~/.local/include) link_directories(~/.local/lib) file(GLOB SOURCES lib/*.c) -add_executable(wolff src/wolff.c ${SOURCES}) +add_executable(wolff_potts src/wolff_potts.c ${SOURCES}) +add_executable(wolff_vector src/wolff_vector.c ${SOURCES}) find_package(OpenMP) if (OPENMP_FOUND) @@ -14,8 +17,8 @@ if (OPENMP_FOUND) set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") endif() -target_link_libraries(wolff gsl m cblas fftw3 tcmalloc profiler) -#target_link_libraries(pt_wolff gsl m jst cblas tcmalloc profiler) +target_link_libraries(wolff_potts gsl m cblas fftw3) +target_link_libraries(wolff_vector gsl m cblas fftw3) -install(TARGETS wolff DESTINATION bin) +install(TARGETS wolff_potts wolff_vector DESTINATION bin) diff --git a/lib/dihedral.c b/lib/dihedral.c new file mode 100644 index 0000000..623041a --- /dev/null +++ b/lib/dihedral.c @@ -0,0 +1,28 @@ + +#include "dihedral.h" + +dihedral_t *dihedral_compose(q_t q, q_t g1i, const dihedral_t *g2) { + // we only need to consider the action of reflections + dihedral_t *g3 = (dihedral_t *)malloc(1 * sizeof(dihedral_t)); + + g3->r = !g2->r; + g3->i = (g1i + q - g2->i) % q; + + return g3; +} + +q_t dihedral_act(q_t q, q_t gi, q_t s) { + // we only need to consider the action of reflections + + return (gi + q - s) % q; +} + +q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) { + if (g->r) { + return (q - ((q + s - g->i) % q)) % q; + } else { + return (q + s - g->i) % q; + } +} + + diff --git a/lib/dihedral.h b/lib/dihedral.h new file mode 100644 index 0000000..813e796 --- /dev/null +++ b/lib/dihedral.h @@ -0,0 +1,17 @@ + +#include +#include + +#include "types.h" + +typedef struct { + q_t i; + bool r; +} dihedral_t; + +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_inverse_act(q_t q, const dihedral_t *g, q_t s); + diff --git a/lib/orthogonal.c b/lib/orthogonal.c new file mode 100644 index 0000000..7e0a71b --- /dev/null +++ b/lib/orthogonal.c @@ -0,0 +1,95 @@ + +#include "orthogonal.h" + +void vector_replace(q_t n, double *v1, const double *v2) { + for (q_t i = 0; i < n; i++) { + v1[i] = v2[i]; + } +} + +void vector_add(q_t n, double *v1, const double *v2) { + for (q_t i = 0; i < n; i++) { + v1[i] += v2[i]; + } +} + +void vector_subtract(q_t n, double *v1, const double *v2) { + for (q_t i = 0; i < n; i++) { + v1[i] -= v2[i]; + } +} + +double *vector_rotate(q_t n, double *rot, double *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, double *v1, 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, double *r, 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 new file mode 100644 index 0000000..a763b08 --- /dev/null +++ b/lib/orthogonal.h @@ -0,0 +1,24 @@ + +#include +#include +#include +#include + +#include "types.h" + +void vector_replace(q_t n, double *v1, const double *v2); + +void vector_add(q_t n, double *v1, const double *v2); + +void vector_subtract(q_t n, double *v1, const double *v2); + +double *vector_rotate(q_t n, double *rot, double *vec); + +double *vector_rotate_inverse(q_t n, const double *rot, const double *vec); + +double vector_dot(q_t n, double *v1, double *v2); + +double *orthogonal_rotate(q_t n, double *m1, double *m2); + +double *gen_rot(gsl_rng *r, q_t n); + diff --git a/lib/wolff.h b/lib/wolff.h index ce7040a..29286ba 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -20,6 +20,8 @@ #include "graph.h" #include "tree.h" #include "measurement.h" +#include "orthogonal.h" +#include "dihedral.h" typedef struct { graph_t *g; @@ -29,12 +31,27 @@ typedef struct { 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; + double *spins; + double T; + double (*J)(double); + double (*H)(double *); + double *R; + double E; + double *M; + q_t n; +} vector_state_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); + graph_t *graph_add_ext(const graph_t *g); diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c index 698d58c..e63623c 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -2,7 +2,6 @@ #include "wolff.h" v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) { - q_t s0 = s->spins[v0]; v_t nv = 0; ll_t *stack = NULL; // create a new stack @@ -16,40 +15,47 @@ v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) { // if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped if (!marks[v]) { - q_t s_old = s->spins[v]; - q_t s_new = (s->spins[v] + step) % s->q; + q_t s_old, s_new; + dihedral_t *R_tmp; - s->spins[v] = s_new; // flip the vertex //tree_insert(&T, v); marks[v] = true; + if (v == s->g->nv - 1) { + R_tmp = dihedral_compose(s->q, step, s->R); + } else { + s_old = s->spins[v]; + s_new = dihedral_act(s->q, step, 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]; - q_t sn = s->spins[vn]; + q_t sn; + if (vn != s->g->nv - 1) { + sn = s->spins[vn]; + } double prob; bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1); - q_t M_ind_0; - q_t M_ind_1; - if (is_ext) { + q_t rs_old, rs_new; if (vn == s->g->nv - 1) { - M_ind_0 = (s_old + s->q - sn) % s->q; - M_ind_1 = (s_new + s->q - sn) % s->q; + rs_old = dihedral_inverse_act(s->q, s->R, s_old); + rs_new = dihedral_inverse_act(s->q, s->R, s_new); } else { - M_ind_0 = (sn + s->q - s_old) % s->q; - M_ind_1 = (sn + s->q - s_new) % s->q; + rs_old = dihedral_inverse_act(s->q, s->R, sn); + rs_new = dihedral_inverse_act(s->q, R_tmp, sn); } - prob = s->H_probs[M_ind_1 * s->q + M_ind_0]; - s->M[M_ind_0]--; - s->M[M_ind_1]++; - s->E += - s->H[M_ind_1] + s->H[M_ind_0]; + prob = s->H_probs[rs_new * s->q + rs_old]; + s->M[rs_old]--; + s->M[rs_new]++; + s->E += - s->H[rs_new] + s->H[rs_old]; } else { - M_ind_0 = (s_old + s->q - sn) % s->q; - M_ind_1 = (s_new + s->q - sn) % s->q; + q_t M_ind_0 = (s_old + s->q - sn) % s->q; + q_t M_ind_1 = (s_new + s->q - sn) % s->q; prob = s->J_probs[M_ind_1 * s->q + M_ind_0]; s->E += - s->J[M_ind_1] + s->J[M_ind_0]; } @@ -59,6 +65,101 @@ v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) { } } + if (v == s->g->nv - 1) { + free(s->R); + s->R = R_tmp; + } else { + s->spins[v] = 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; +} + +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]) { + double *s_old, *s_new, *R_tmp; + + //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]; + double *sn; + if (vn != s->g->nv - 1) { + sn = &(s->spins[s->n * vn]); + } + double prob; + + bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1); + + if (is_ext) { + double *rs_old, *rs_new; + if (vn == s->g->nv - 1) { + 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(rs_old) - s->H(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); + //printf("(%g %g) (%g %g) (%g %g) %g\n", s_old[0], s_old[1], s_new[0], s_new[1], sn[0], sn[1], dE); + //getchar(); + 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++; } diff --git a/src/wolff.c b/src/wolff.c deleted file mode 100644 index 7b6b56f..0000000 --- a/src/wolff.c +++ /dev/null @@ -1,340 +0,0 @@ - -#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; - q_t q = 2; - D_t D = 2; - 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 pretend_ising = false; - bool planar_potts = false; - 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:Ips")) != -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 'J': - J[J_ind] = atof(optarg); - J_ind++; - 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 'I': - pretend_ising = true; - break; - case 'p': - planar_potts = true; - break; - case 's': - silent = true; - break; - default: - exit(EXIT_FAILURE); - } - } - - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, rand_seed()); - - if (pretend_ising) { - q = 2; - H[1] = -H[0]; - J[1] = -J[0]; - } - - if (planar_potts) { - for (q_t i = 0; i < q; i++) { - J[i] = cos(2 * M_PI * i / ((double)q)); - } - } - - ising_state_t *s = (ising_state_t *)calloc(1, sizeof(ising_state_t)); - - graph_t *h = graph_create_square(D, L); - s->g = graph_add_ext(h); - - s->q = q; - - s->spins = (q_t *)calloc(h->nv + 1, sizeof(q_t)); - - s->T = T; - s->H = H; - s->J = J; - - s->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++) { - s->J_probs[q * i + j] = 1.0 - exp((s->J[i] - s->J[j]) / T); - } - } - s->H_probs = (double *)calloc(pow(q, 2), sizeof(double)); - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - s->H_probs[q * i + j] = 1.0 - exp((s->H[i] - s->H[j]) / T); - } - } - - s->M = (v_t *)calloc(q, sizeof(v_t)); - s->M[0] = h->nv; - s->E = - ((double)h->ne) * s->J[0] - ((double)h->nv) * s->H[0]; - - double diff = 1e31; - count_t n_runs = 0; - - meas_t *E, *clust, **M, **sE, ***sM, **lifetimes; - - M = (meas_t **)malloc(q * sizeof(meas_t *)); - lifetimes = (meas_t **)malloc(q * sizeof(meas_t *)); - for (q_t i = 0; i < q; i++) { - M[i] = (meas_t *)calloc(1, sizeof(meas_t)); - lifetimes[i] = (meas_t *)calloc(1, sizeof(meas_t)); - } - - E = calloc(1, sizeof(meas_t)); - clust = calloc(1, sizeof(meas_t)); - - sE = (meas_t **)malloc(q * sizeof(meas_t *)); - sM = (meas_t ***)malloc(q * sizeof(meas_t **)); - - 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)); - } - } - - count_t *freqs = (count_t *)calloc(q, sizeof(count_t)); - count_t lifetime_n = 0; - q_t cur_M = 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(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x); - - count_t n_flips = 0; - - q_t max_M_i; - v_t max_M; - q_t n_at_max; - - while (n_flips / h->nv < n) { - v_t v0 = gsl_rng_uniform_int(r, h->nv); - q_t step = 1 + gsl_rng_uniform_int(r, q - 1); - - v_t tmp_flips = flip_cluster(s, v0, step, r); - n_flips += tmp_flips; - - max_M_i = 0; - max_M = 0; - n_at_max = 0; - - for (q_t i = 0; i < q; i++) { - if (s->M[i] > max_M) { - max_M = s->M[i]; - max_M_i = i; - n_at_max = 1; - } else if (s->M[i] == max_M) { - n_at_max++; - } - } - - if (n_at_max == 1) { - if (max_M_i == cur_M) { - lifetime_n++; - } else { - if (cur_M != MAX_Q) { - update_meas(lifetimes[cur_M], lifetime_n); - } - lifetime_n = 0; - cur_M = max_M_i; - } - } else { - if (cur_M != MAX_Q) { - update_meas(lifetimes[cur_M], lifetime_n); - cur_M = MAX_Q; - } - } - - update_meas(clust, tmp_flips); - } - - for (q_t i = 0; i < q; i++) { - update_meas(M[i], s->M[i]); - } - update_meas(E, s->E); - - if (n_at_max == 1) { - for (q_t i = 0; i < q; i++) { - update_meas(sM[max_M_i][i], s->M[i]); - } - update_meas(sE[max_M_i], s->E); - freqs[max_M_i]++; - } - - diff = fabs(sM[0][0]->dc / sM[0][0]->c); - - 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(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x); - - FILE *outfile = fopen("out.m", "a"); - - fprintf(outfile, "<|D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,J->{", D, L, q, T); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", J[i]); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},H->{"); - 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, E->dx / h->nv, E->c / h->nv, E->dc / 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", M[i]->dx / h->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},\\[Chi]->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", M[i]->c / h->nv); - if (i != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},\\[Delta]\\[Chi]->{"); - for (q_t i = 0; i < q; i++) { - fprintf(outfile, "%.15f", M[i]->dc / h->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 / h->nv, i, sE[i]->dx / h->nv, i, sE[i]->c / h->nv, i, sE[i]->dc / h->nv, i); - for (q_t j = 0; j < q; j++) { - fprintf(outfile, "%.15f", sM[i][j]->x / h->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", sM[i][j]->dx / h->nv); - if (j != q-1) { - fprintf(outfile, ","); - } - } - fprintf(outfile, "},Subscript[\\[Chi],%" PRIq "]->{", i); - for (q_t j = 0; j < q; j++) { - fprintf(outfile, "%.15f", sM[i][j]->c / h->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", sM[i][j]->dc / h->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); - } - for (q_t i = 0; i < q; i++) { - fprintf(outfile, ",Subscript[t,%" PRIq "]->%.15f,Subscript[\\[Delta]t,%" PRIq "]->%.15f", i, lifetimes[i]->x, i, lifetimes[i]->dx); - } - fprintf(outfile, ",Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f|>\n", clust->x / h->nv, clust->dx / h->nv); - - 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(lifetimes[i]); - } - free(freqs); - free(sE); - free(s->H_probs); - free(s->J_probs); - free(s->M); - free(s->spins); - graph_free(s->g); - free(s); - free(H); - free(J); - graph_free(h); - gsl_rng_free(r); - - return 0; -} - diff --git a/src/wolff_potts.c b/src/wolff_potts.c new file mode 100644 index 0000000..7d9beb3 --- /dev/null +++ b/src/wolff_potts.c @@ -0,0 +1,349 @@ + +#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; + q_t q = 2; + D_t D = 2; + 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 pretend_ising = false; + bool planar_potts = false; + 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:Ips")) != -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 'J': + J[J_ind] = atof(optarg); + J_ind++; + 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 'I': + pretend_ising = true; + break; + case 'p': + planar_potts = true; + break; + case 's': + silent = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(r, rand_seed()); + + if (pretend_ising) { + q = 2; + H[1] = -H[0]; + J[1] = -J[0]; + } + + if (planar_potts) { + for (q_t i = 0; i < q; i++) { + J[i] = cos(2 * M_PI * i / ((double)q)); + } + } + + ising_state_t *s = (ising_state_t *)calloc(1, sizeof(ising_state_t)); + + graph_t *h = graph_create_square(D, L); + s->g = graph_add_ext(h); + + s->q = q; + + s->spins = (q_t *)calloc(h->nv, sizeof(q_t)); + + s->T = T; + s->H = H; + s->J = J; + s->R = (dihedral_t *)calloc(1, sizeof(dihedral_t)); + + s->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++) { + s->J_probs[q * i + j] = 1.0 - exp((s->J[i] - s->J[j]) / T); + } + } + s->H_probs = (double *)calloc(pow(q, 2), sizeof(double)); + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + s->H_probs[q * i + j] = 1.0 - exp((s->H[i] - s->H[j]) / T); + } + } + + s->M = (v_t *)calloc(q, sizeof(v_t)); + s->M[0] = h->nv; + s->E = - ((double)h->ne) * s->J[0] - ((double)h->nv) * s->H[0]; + + double diff = 1e31; + count_t n_runs = 0; + + meas_t *E, *clust, **M, **sE, ***sM, **lifetimes; + + M = (meas_t **)malloc(q * sizeof(meas_t *)); + lifetimes = (meas_t **)malloc(q * sizeof(meas_t *)); + for (q_t i = 0; i < q; i++) { + M[i] = (meas_t *)calloc(1, sizeof(meas_t)); + lifetimes[i] = (meas_t *)calloc(1, sizeof(meas_t)); + } + + E = calloc(1, sizeof(meas_t)); + clust = calloc(1, sizeof(meas_t)); + + sE = (meas_t **)malloc(q * sizeof(meas_t *)); + sM = (meas_t ***)malloc(q * sizeof(meas_t **)); + + 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)); + } + } + + count_t *freqs = (count_t *)calloc(q, sizeof(count_t)); + count_t lifetime_n = 0; + q_t cur_M = 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(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x); + + count_t n_flips = 0; + + q_t max_M_i; + v_t max_M; + q_t n_at_max; + + while (n_flips / h->nv < n) { + v_t v0 = gsl_rng_uniform_int(r, h->nv); + q_t step; + + if (q == 2) { + step = 1; + } else { + step = gsl_rng_uniform_int(r, q); + } + + v_t tmp_flips = flip_cluster(s, v0, step, r); + n_flips += tmp_flips; + + max_M_i = 0; + max_M = 0; + n_at_max = 0; + + for (q_t i = 0; i < q; i++) { + if (s->M[i] > max_M) { + max_M = s->M[i]; + max_M_i = i; + n_at_max = 1; + } else if (s->M[i] == max_M) { + n_at_max++; + } + } + + if (n_at_max == 1) { + if (max_M_i == cur_M) { + lifetime_n++; + } else { + if (cur_M != MAX_Q) { + update_meas(lifetimes[cur_M], lifetime_n); + } + lifetime_n = 0; + cur_M = max_M_i; + } + } else { + if (cur_M != MAX_Q) { + update_meas(lifetimes[cur_M], lifetime_n); + cur_M = MAX_Q; + } + } + + update_meas(clust, tmp_flips); + } + + for (q_t i = 0; i < q; i++) { + update_meas(M[i], s->M[i]); + } + update_meas(E, s->E); + + if (n_at_max == 1) { + for (q_t i = 0; i < q; i++) { + update_meas(sM[max_M_i][i], s->M[i]); + } + update_meas(sE[max_M_i], s->E); + freqs[max_M_i]++; + } + + diff = fabs(sM[0][0]->dc / sM[0][0]->c); + + 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(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x); + + FILE *outfile = fopen("out.m", "a"); + + fprintf(outfile, "<|D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,J->{", D, L, q, T); + for (q_t i = 0; i < q; i++) { + fprintf(outfile, "%.15f", J[i]); + if (i != q-1) { + fprintf(outfile, ","); + } + } + fprintf(outfile, "},H->{"); + 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, E->dx / h->nv, E->c / h->nv, E->dc / 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", M[i]->dx / h->nv); + if (i != q-1) { + fprintf(outfile, ","); + } + } + fprintf(outfile, "},\\[Chi]->{"); + for (q_t i = 0; i < q; i++) { + fprintf(outfile, "%.15f", M[i]->c / h->nv); + if (i != q-1) { + fprintf(outfile, ","); + } + } + fprintf(outfile, "},\\[Delta]\\[Chi]->{"); + for (q_t i = 0; i < q; i++) { + fprintf(outfile, "%.15f", M[i]->dc / h->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 / h->nv, i, sE[i]->dx / h->nv, i, sE[i]->c / h->nv, i, sE[i]->dc / h->nv, i); + for (q_t j = 0; j < q; j++) { + fprintf(outfile, "%.15f", sM[i][j]->x / h->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", sM[i][j]->dx / h->nv); + if (j != q-1) { + fprintf(outfile, ","); + } + } + fprintf(outfile, "},Subscript[\\[Chi],%" PRIq "]->{", i); + for (q_t j = 0; j < q; j++) { + fprintf(outfile, "%.15f", sM[i][j]->c / h->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", sM[i][j]->dc / h->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); + } + for (q_t i = 0; i < q; i++) { + fprintf(outfile, ",Subscript[t,%" PRIq "]->%.15f,Subscript[\\[Delta]t,%" PRIq "]->%.15f", i, lifetimes[i]->x, i, lifetimes[i]->dx); + } + fprintf(outfile, ",Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f|>\n", clust->x / h->nv, clust->dx / h->nv); + + 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(lifetimes[i]); + } + free(lifetimes); + free(freqs); + free(sE); + free(s->H_probs); + free(s->J_probs); + free(s->M); + free(s->spins); + free(s->R); + graph_free(s->g); + free(s); + free(H); + free(J); + graph_free(h); + gsl_rng_free(r); + + return 0; +} + diff --git a/src/wolff_vector.c b/src/wolff_vector.c new file mode 100644 index 0000000..c6ab695 --- /dev/null +++ b/src/wolff_vector.c @@ -0,0 +1,202 @@ + +#include + +double identity(double x) { + return x; +} + +double zero(double *x) { + return 0.0; +} + +double test(double *x) { + return -x[1]; +} + +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; + + int opt; + q_t H_ind = 0; + + while ((opt = getopt(argc, argv, "N:n:D:L:q:T:H:m:e:s")) != -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; + 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->T = T; + s->H = test; + 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 * (double)h->nv; + s->E = - ((double)h->ne) * s->J(1.0) - s->H(s->M); + + printf("%g %g\n", s->E, s->M[0]); + + double diff = 1e31; + count_t n_runs = 0; + + meas_t *E, *clust, **M; + + M = (meas_t **)malloc(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)); + + 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(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, 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; + + update_meas(clust, tmp_flips); + } + + for (q_t i = 0; i < q; i++) { + update_meas(M[i], s->M[i]); + } + update_meas(E, s->E); + + diff = fabs(E->dc / E->c); + + 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(E->dx / E->x), M[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x); + + FILE *outfile = fopen("out.m", "a"); + + fprintf(outfile, "<|D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f", D, L, q, T); + fprintf(outfile, "},E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->{", E->x / h->nv, E->dx / h->nv, E->c / h->nv, E->dc / 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", M[i]->dx / h->nv); + if (i != q-1) { + fprintf(outfile, ","); + } + } + fprintf(outfile, "},\\[Chi]->{"); + for (q_t i = 0; i < q; i++) { + fprintf(outfile, "%.15f", M[i]->c / h->nv); + if (i != q-1) { + fprintf(outfile, ","); + } + } + fprintf(outfile, "},\\[Delta]\\[Chi]->{"); + for (q_t i = 0; i < q; i++) { + fprintf(outfile, "%.15f", M[i]->dc / h->nv); + if (i != q-1) { + fprintf(outfile, ","); + } + } + fprintf(outfile, "}|>\n"); + + 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-70-g09d2