diff options
Diffstat (limited to 'src/wolff_vector.c')
-rw-r--r-- | src/wolff_vector.c | 87 |
1 files changed, 82 insertions, 5 deletions
diff --git a/src/wolff_vector.c b/src/wolff_vector.c index bbfc4da..b1584e6 100644 --- a/src/wolff_vector.c +++ b/src/wolff_vector.c @@ -47,11 +47,14 @@ int main(int argc, char *argv[]) { 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:q:T:H:m:e:s")) != -1) { + while ((opt = getopt(argc, argv, "N:n:D:L:q:T:H:m:e:saS:W:")) != -1) { switch (opt) { case 'N': N = (count_t)atof(optarg); @@ -84,6 +87,15 @@ int main(int argc, char *argv[]) { 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); } @@ -120,10 +132,9 @@ int main(int argc, char *argv[]) { s->M[0] = 1.0 * (double)h->nv; s->E = - ((double)h->ne) * s->J(1.0) - s->H(s->n, s->H_info, s->M); - printf("%g %g\n", s->E, s->M[0]); - double diff = 1e31; count_t n_runs = 0; + count_t n_steps = 0; meas_t *E, *clust, **M; @@ -135,6 +146,13 @@ int main(int argc, char *argv[]) { 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 @@ -151,7 +169,16 @@ int main(int argc, char *argv[]) { free(step); n_flips += tmp_flips; - update_meas(clust, tmp_flips); + if (n_runs > 0) { + n_steps++; + update_meas(clust, tmp_flips); + } + + if (record_autocorrelation && n_runs > 0) { + if (n_steps % ac_skip == 0) { + update_autocorr(autocorr, s->E); + } + } } for (q_t i = 0; i < q; i++) { @@ -170,6 +197,56 @@ int main(int argc, char *argv[]) { ", 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); + 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 ",q->%" PRIq ",T->%.15f", D, L, q, T); @@ -201,7 +278,7 @@ int main(int argc, char *argv[]) { fprintf(outfile, ","); } } - fprintf(outfile, "}|>\n"); + 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, clust->dx / h->nv, clust->c / h->nv, clust->dc / h->nv,tau); fclose(outfile); |