From fbbc4d9655835c0d6fdf58f231e59d7007a99407 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 23 Apr 2018 18:32:49 -0400 Subject: lots of changes --- lib/cluster.h | 8 ++++++ lib/convex.c | 1 - lib/measurement.h | 2 ++ lib/yule_walker.c | 41 +++++++++++++++++++++++++++++++ lib/yule_walker.h | 12 +++++++++ src/wolff_potts.c | 37 ++++++++++++++++++++++++---- src/wolff_vector.c | 71 +++++++++++++++++++++++++++++++++++++++++++++++++++--- 7 files changed, 162 insertions(+), 10 deletions(-) create mode 100644 lib/yule_walker.c create mode 100644 lib/yule_walker.h diff --git a/lib/cluster.h b/lib/cluster.h index 2de17e5..d118735 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -22,6 +22,7 @@ #include "orthogonal.h" #include "dihedral.h" #include "dihinf.h" +#include "yule_walker.h" typedef struct { graph_t *g; @@ -62,6 +63,13 @@ typedef struct { 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); diff --git a/lib/convex.c b/lib/convex.c index df506f9..7255ad4 100644 --- a/lib/convex.c +++ b/lib/convex.c @@ -93,7 +93,6 @@ double *get_convex_minorant(count_t n, double *Gammas) { while (L != NULL) { free(L->p); - list_t *L_save = L; L = L->next; free(L); } diff --git a/lib/measurement.h b/lib/measurement.h index a24df9f..eaa260b 100644 --- a/lib/measurement.h +++ b/lib/measurement.h @@ -1,4 +1,6 @@ +#pragma once + #include #include diff --git a/lib/yule_walker.c b/lib/yule_walker.c new file mode 100644 index 0000000..77b3e96 --- /dev/null +++ b/lib/yule_walker.c @@ -0,0 +1,41 @@ + +#include "yule_walker.h" + +double yule_walker(const autocorr_t *a) { + gsl_vector *rhos = gsl_vector_alloc(a->W); + gsl_matrix *mat = gsl_matrix_alloc(a->W, a->W); + gsl_vector *phis = gsl_vector_alloc(a->W); + + for (count_t i = 0; i < a->W; i++) { + gsl_vector_set(rhos, i, rho(a, i)); + } + + for (count_t i = 0; i < a->W; i++) { + gsl_matrix_set(mat, i, i, 1.0); + + for (count_t j = 0; j < i; j++) { + gsl_matrix_set(mat, i, j, gsl_vector_get(rhos, i - 1 - j)); + } + + for (count_t j = 0; j < a->W - 1 - i; j++) { + gsl_matrix_set(mat, i, i + 1 + j, gsl_vector_get(rhos, j)); + } + } + + int out = gsl_linalg_cholesky_solve(mat, rhos, phis); + + double rhophi = 0; + double onephi = 0; + + for (count_t i = 0; i < a->W; i++) { + rhophi += gsl_vector_get(rhos, i) * gsl_vector_get(phis, i); + onephi += gsl_vector_get(phis, i); + } + + gsl_vector_free(rhos); + gsl_matrix_free(mat); + gsl_vector_free(phis); + + return (1 - rhophi) / pow(1 - onephi, 2); +} + diff --git a/lib/yule_walker.h b/lib/yule_walker.h new file mode 100644 index 0000000..883ea03 --- /dev/null +++ b/lib/yule_walker.h @@ -0,0 +1,12 @@ + +#pragma once + +#include +#include +#include + +#include "types.h" +#include "measurement.h" + +double yule_walker(const autocorr_t *a); + diff --git a/src/wolff_potts.c b/src/wolff_potts.c index e209043..d200b35 100644 --- a/src/wolff_potts.c +++ b/src/wolff_potts.c @@ -22,6 +22,7 @@ int main(int argc, char *argv[]) { bool snapshots = false; bool snapshot = false; bool record_autocorrelation = false; + bool record_distribution = false; count_t W = 10; count_t ac_skip = 1; @@ -29,7 +30,7 @@ int main(int argc, char *argv[]) { 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:")) != -1) { + while ((opt = getopt(argc, argv, "N:n:D:L:q:T:J:H:m:e:IpsSPak:W:d")) != -1) { switch (opt) { case 'N': N = (count_t)atof(optarg); @@ -87,6 +88,9 @@ int main(int argc, char *argv[]) { case 'W': W = (count_t)atof(optarg); break; + case 'd': + record_distribution = true; + break; default: exit(EXIT_FAILURE); } @@ -173,6 +177,11 @@ int main(int argc, char *argv[]) { autocorr->OO = (double *)calloc(2 * W + 1, sizeof(double)); } + count_t *cluster_dist; + if (record_distribution) { + cluster_dist = (count_t *)calloc(h->nv, sizeof(count_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 @@ -201,6 +210,10 @@ int main(int argc, char *argv[]) { if (record_autocorrelation && n_steps % ac_skip == 0) { update_autocorr(autocorr, s->E); } + + if (record_distribution) { + cluster_dist[tmp_flips - 1]++; + } } } @@ -269,6 +282,7 @@ int main(int argc, char *argv[]) { } double tau = 0; + double tauyw; bool tau_failed = false; if (record_autocorrelation) { @@ -304,14 +318,15 @@ int main(int argc, char *argv[]) { ttau += conv_Gamma[i]; } + tau = ttau * ac_skip * clust->x / h->nv; + tauyw = yule_walker(autocorr) * ac_skip * clust->x / h->nv; + 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) { @@ -320,7 +335,7 @@ int main(int argc, char *argv[]) { FILE *outfile = fopen("out.m", "a"); - fprintf(outfile, "<|N->%" PRIcount ",D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,J->{", N, D, L, q, T); + 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", J[i]); if (i != q-1) { @@ -396,7 +411,19 @@ int main(int argc, char *argv[]) { 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|>\n", clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau); + fprintf(outfile, ",Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f,\\[Tau]yw->%.15f", clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau, tauyw); + if (record_distribution) { + fprintf(outfile, ",S->{"); + for (v_t i = 0; i < h->nv; i++) { + fprintf(outfile, "%" PRIcount, cluster_dist[i]); + if (i != h->nv - 1) { + fprintf(outfile, ","); + } + } + fprintf(outfile, "}"); + free(cluster_dist); + } + fprintf(outfile, "|>\n"); fclose(outfile); diff --git a/src/wolff_vector.c b/src/wolff_vector.c index c87e60d..da93a15 100644 --- a/src/wolff_vector.c +++ b/src/wolff_vector.c @@ -35,6 +35,28 @@ 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 tmp = 0; + + for (q_t i = 0; i < n; i++) { + tmp += pow(x[i], 4); + } + + return - H_info[0] * tmp; +} + +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; @@ -48,13 +70,14 @@ int main(int argc, char *argv[]) { 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:")) != -1) { + 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); @@ -96,6 +119,22 @@ int main(int argc, char *argv[]) { 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); } @@ -119,7 +158,20 @@ int main(int argc, char *argv[]) { s->H_info = H; s->T = T; - s->H = dot; + 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)); @@ -190,7 +242,7 @@ int main(int argc, char *argv[]) { meas_update(aM, sqrt(aM_val)); meas_update(E, s->E); - diff = fabs(meas_dc(aM) / meas_c(aM)); + diff = fabs(meas_dx(clust) / clust->x); n_runs++; } @@ -237,6 +289,17 @@ int main(int argc, char *argv[]) { 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); @@ -254,7 +317,7 @@ int main(int argc, char *argv[]) { FILE *outfile = fopen("out.m", "a"); - fprintf(outfile, "<|N->%" PRIcount ",D->%" PRID ",L->%" PRIL ",q->%" PRIq ",T->%.15f,H->{", N, D, L, q, T); + 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) { -- cgit v1.2.3-54-g00ecf