From f3d53bac095fd47f7c10f54007ffa11ab0f1ac1f Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 2 Nov 2017 23:00:48 -0400 Subject: added new method for computing autocorrelation --- lib/convex.c | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ lib/convex.h | 21 ++++++++++++++ lib/wolff.h | 1 + src/wolff.c | 63 +++++++++++++++++++++++++++------------- 4 files changed, 159 insertions(+), 19 deletions(-) create mode 100644 lib/convex.c create mode 100644 lib/convex.h diff --git a/lib/convex.c b/lib/convex.c new file mode 100644 index 0000000..33af1b4 --- /dev/null +++ b/lib/convex.c @@ -0,0 +1,93 @@ + +#include "convex.h" + +double slope(point_t *P, point_t *Q) { + return (Q->y - P->y) / ((double)(Q->x) - (double)(P->x)); +} + +double *get_convex_minorant(uint64_t n, double *Gammas) { + list_t *L = (list_t *)calloc(1, sizeof(list_t)); + L->p = (point_t *)calloc(1, sizeof(point_t)); + L->p->x = 0; + L->p->y = Gammas[0]; + + list_t *pos = L; + + for (uint64_t i = 1; i < n; i++) { + pos->next = (list_t *)calloc(1, sizeof(list_t)); + pos->next->p = (point_t *)calloc(1, sizeof(point_t)); + pos->next->p->x = i; + pos->next->p->y = Gammas[i]; + pos->next->prev = pos; + pos = pos->next; + } + + pos->next = (list_t *)calloc(1, sizeof(list_t)); + pos->next->p = (point_t *)calloc(1, sizeof(point_t)); + pos->next->p->x = n; + pos->next->p->y = 0; + pos->next->prev = pos; + + list_t *X = L; + list_t *Y = L->next; + list_t *Z = Y->next; + + while (true) { + if (slope(X->p, Y->p) <= slope(Y->p, Z->p)) { + X = Y; + Y = Z; + if (Z->next == NULL) { + break; + } else { + Z = Z->next; + } + } else { + Y->prev->next = Y->next; + Y->next->prev = Y->prev; + free(Y->p); + free(Y); + if (X->prev != NULL) { + Y = X; + X = X->prev; + } else { + if (Z->next != NULL) { + Y = Z; + Z = Z->next; + } else { + break; + } + } + } + } + + pos = L; + + double *g = (double *)calloc(n, sizeof(double)); + double rho = 0; + + for (uint64_t i = 0; i < n; i++) { + if (i > pos->next->p->x) { + pos = pos->next; + } + + g[i] = pos->p->y + ((double)i - (double)(pos->p->x)) * (pos->next->p->y - pos->p->y) / ((double)(pos->next->p->x) - (double)(pos->p->x)); + + if (Gammas[i] - g[i] > rho) { + rho = Gammas[i] - g[i]; + } + } + + for (uint64_t i = 0; i < n; i++) { + g[i] += rho / 2; + } + + while (L != NULL) { + free(L->p); + list_t *L_save = L; + L = L->next; + free(L); + } + + return g; +} + diff --git a/lib/convex.h b/lib/convex.h new file mode 100644 index 0000000..29f4ae0 --- /dev/null +++ b/lib/convex.h @@ -0,0 +1,21 @@ + +#pragma once + +#include +#include +#include +#include + +typedef struct { + uint64_t x; + double y; +} point_t; + +typedef struct list_tag { + struct list_tag *prev; + struct list_tag *next; + point_t *p; +} list_t; + +double *get_convex_minorant(uint64_t n, double *Gammas); + diff --git a/lib/wolff.h b/lib/wolff.h index 66ce493..dfe7b1c 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -17,6 +17,7 @@ #include #include "queue.h" +#include "convex.h" typedef enum { WOLFF, diff --git a/src/wolff.c b/src/wolff.c index 1e62c8c..a806681 100644 --- a/src/wolff.c +++ b/src/wolff.c @@ -136,7 +136,7 @@ int main(int argc, char *argv[]) { uint64_t n_steps = 0; double clust_per_sweep = 0; - meas_t *M, *aM, *eM, *mM, *E, *eE, *mE, *clust; + meas_t *M, *aM, *eM, *mM, *E, *eE, *mE, *clust, *everyE, *blockE; M = calloc(1, sizeof(meas_t)); aM = calloc(1, sizeof(meas_t)); @@ -146,12 +146,18 @@ int main(int argc, char *argv[]) { eE = calloc(1, sizeof(meas_t)); mE = calloc(1, sizeof(meas_t)); clust = calloc(1, sizeof(meas_t)); + everyE = calloc(1, sizeof(meas_t)); + + blockE = calloc(1,sizeof(meas_t)); + + uint64_t blocksize = 1000; + double Etot = 0; autocorr_t *autocorr; if (record_autocorrelation) { autocorr = (autocorr_t *)calloc(1, sizeof(autocorr_t)); - autocorr->W = W; - autocorr->OO = (double *)calloc(W, sizeof(double)); + autocorr->W = 2 * W + 1; + autocorr->OO = (double *)calloc(2 * W + 1, sizeof(double)); } printf("\n"); @@ -167,9 +173,17 @@ int main(int argc, char *argv[]) { uint32_t tmp_flips = wolff_step(T, H, s, sim, r, bond_probs); n_flips += tmp_flips; n_clust++; - n_steps++; + if (n_runs > 0){ + n_steps++; + } if (record_autocorrelation && n_runs > 0) { + update_meas(everyE, s->H); + if (n_steps % blocksize == 0) { + update_meas(blockE, Etot / blocksize); + Etot = 0; + } + Etot += s->H; if (n_steps % ac_skip == 0) { update_autocorr(autocorr, s->H); } @@ -213,25 +227,36 @@ int main(int argc, char *argv[]) { double tau = 0; double dtau = 0; + if (record_autocorrelation) { - double kappa = rho(autocorr, W - 1) / rho(autocorr, W - 2); - double R = rho(autocorr, W - 1) / (1-kappa); - - double ttau = 0.5 + R; - double X = 0.5; - double Y = 0; - for (uint64_t i = 0; i < W - 1; i++) { - ttau += rho(autocorr, i); - X += pow(rho(autocorr, i), 2); - if (i > 0) { - Y += 0.5 * pow(rho(autocorr, i) - rho(autocorr, i-1), 2); - } else { - Y += 0.5 * pow(rho(autocorr, i) - 1, 2); + double *Gammas = (double *)malloc((W + 1) * sizeof(double)); + + Gammas[0] = 1 + rho(autocorr, 0); + printf("%g ", Gammas[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; } } - double dttau = sqrt(4.0 / autocorr->n * (pow(ttau, 2) * (W +(1+kappa)/(1-kappa)-0.5) + kappa * X / pow(1-kappa, 2) + pow(kappa, 2) * Y / pow(1-kappa, 4))); + if (n == W) { + printf("WARNING: correlation function never hit the noise floor.\n"); + } + + double *conv_Gamma = get_convex_minorant(n, Gammas); + + double ttau = - 0.5; + + for (uint64_t i = 0; i < n; i++) { + ttau += conv_Gamma[i]; + } + free(Gammas); free(autocorr->OO); while (autocorr->Op != NULL) { stack_pop_d(&(autocorr->Op)); @@ -239,7 +264,7 @@ int main(int argc, char *argv[]) { free(autocorr); tau = ttau * ac_skip * clust->x / h->nv; - dtau = tau * sqrt(pow(dttau / ttau, 2) + pow(clust->dx / clust->x, 2)); + dtau = 0; } fprintf(outfile, -- cgit v1.2.3-70-g09d2