diff options
Diffstat (limited to 'lib')
-rw-r--r-- | lib/measurement.c | 62 | ||||
-rw-r--r-- | lib/measurement.h | 33 | ||||
-rw-r--r-- | lib/wolff.h | 27 | ||||
-rw-r--r-- | lib/wolff_tools.c | 58 |
4 files changed, 96 insertions, 84 deletions
diff --git a/lib/measurement.c b/lib/measurement.c new file mode 100644 index 0000000..435c749 --- /dev/null +++ b/lib/measurement.c @@ -0,0 +1,62 @@ + +#include "measurement.h" + +double add_to_avg(double mx, double x, count_t n) { + return mx * (n / (n + 1.0)) + x * 1.0 / (n + 1.0); +} + +void update_meas(meas_t *m, double x) { + count_t n = m->n; + + m->x = add_to_avg(m->x, x, n); + m->x2 = add_to_avg(m->x2, pow(x, 2), n); + + m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n); + m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n); + + if (n > 1) { + double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2)); + m->dx = sqrt(s2 / n); + m->c = s2; + m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n); + } + + (m->n)++; +} + +void update_autocorr(autocorr_t *OO, double O) { + OO->O = add_to_avg(OO->O, O, OO->n); + OO->O2 = add_to_avg(OO->O2, pow(O, 2), OO->n); + + dll_t *Otmp = OO->Op; + dll_t *Osave; + count_t t = 0; + + while (Otmp != NULL) { + OO->OO[t] = add_to_avg(OO->OO[t], O * (Otmp->x), OO->n - t - 1); + t++; + if (t == OO->W - 1) { + Osave = Otmp; + } + Otmp = Otmp->next; + } + + if (t == OO->W) { + if (OO->W == 1) { + free(OO->Op); + OO->Op = NULL; + } else { + free(Osave->next); + Osave->next = NULL; + } + } + + stack_push_d(&(OO->Op), O); + + OO->n++; +} + +double rho(autocorr_t *o, count_t i) { + return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2)); +} + diff --git a/lib/measurement.h b/lib/measurement.h new file mode 100644 index 0000000..4ced4dc --- /dev/null +++ b/lib/measurement.h @@ -0,0 +1,33 @@ + +#include <math.h> +#include <stdlib.h> + +#include "types.h" +#include "stack.h" + +typedef struct { + uint64_t n; + double x; + double dx; + double x2; + double m2; + double m4; + double c; + double dc; +} meas_t; + +typedef struct { + uint64_t n; + uint64_t W; + double *OO; + dll_t *Op; + double O; + double O2; +} autocorr_t; + +void update_meas(meas_t *m, double x); + +void update_autocorr(autocorr_t *OO, double O); + +double rho(autocorr_t *o, uint64_t i); + diff --git a/lib/wolff.h b/lib/wolff.h index 44e5926..ce7040a 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -19,6 +19,7 @@ #include "convex.h" #include "graph.h" #include "tree.h" +#include "measurement.h" typedef struct { graph_t *g; @@ -33,33 +34,7 @@ typedef struct { q_t q; } ising_state_t; -typedef struct { - uint64_t n; - double x; - double dx; - double x2; - double m2; - double m4; - double c; - double dc; -} meas_t; - -typedef struct { - uint64_t n; - uint64_t W; - double *OO; - dll_t *Op; - double O; - double O2; -} autocorr_t; - v_t flip_cluster(ising_state_t *s, v_t v0, q_t s1, gsl_rng *r); graph_t *graph_add_ext(const graph_t *g); -void update_meas(meas_t *m, double x); - -void update_autocorr(autocorr_t *OO, double O); - -double rho(autocorr_t *o, uint64_t i); - diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c index 20bc529..698d58c 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -71,61 +71,3 @@ v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) { return nv; } -double add_to_avg(double mx, double x, count_t n) { - return mx * (n / (n + 1.0)) + x * 1.0 / (n + 1.0); -} - -void update_meas(meas_t *m, double x) { - count_t n = m->n; - - m->x = add_to_avg(m->x, x, n); - m->x2 = add_to_avg(m->x2, pow(x, 2), n); - - m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n); - m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n); - - if (n > 1) { - double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2)); - m->dx = sqrt(s2 / n); - m->c = s2; - m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n); - } - - (m->n)++; -} - -void update_autocorr(autocorr_t *OO, double O) { - OO->O = add_to_avg(OO->O, O, OO->n); - OO->O2 = add_to_avg(OO->O2, pow(O, 2), OO->n); - - dll_t *Otmp = OO->Op; - dll_t *Osave; - count_t t = 0; - - while (Otmp != NULL) { - OO->OO[t] = add_to_avg(OO->OO[t], O * (Otmp->x), OO->n - t - 1); - t++; - if (t == OO->W - 1) { - Osave = Otmp; - } - Otmp = Otmp->next; - } - - if (t == OO->W) { - if (OO->W == 1) { - free(OO->Op); - OO->Op = NULL; - } else { - free(Osave->next); - Osave->next = NULL; - } - } - - stack_push_d(&(OO->Op), O); - - OO->n++; -} - -double rho(autocorr_t *o, count_t i) { - return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2)); -} |