From 2af9351db3aa97da9b0d3f23d53a593bc96c8a8e Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 2 Feb 2018 18:33:22 -0500 Subject: does potts now, no external libraries --- lib/convex.c | 8 +- lib/convex.h | 6 +- lib/graph.c | 67 ++++++++++++ lib/graph.h | 18 ++++ lib/queue.c | 50 --------- lib/queue.h | 25 ----- lib/rand.c | 20 ++++ lib/rand.h | 7 ++ lib/stack.c | 41 ++++++++ lib/stack.h | 26 +++++ lib/tree.c | 81 +++++++++++++++ lib/tree.h | 22 ++++ lib/types.h | 16 +++ lib/wolff.h | 44 +++----- lib/wolff_tools.c | 300 +++++++++++------------------------------------------- 15 files changed, 379 insertions(+), 352 deletions(-) create mode 100644 lib/graph.c create mode 100644 lib/graph.h delete mode 100644 lib/queue.c delete mode 100644 lib/queue.h create mode 100644 lib/rand.c create mode 100644 lib/rand.h create mode 100644 lib/stack.c create mode 100644 lib/stack.h create mode 100644 lib/tree.c create mode 100644 lib/tree.h create mode 100644 lib/types.h (limited to 'lib') diff --git a/lib/convex.c b/lib/convex.c index 096504d..df506f9 100644 --- a/lib/convex.c +++ b/lib/convex.c @@ -5,7 +5,7 @@ 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) { +double *get_convex_minorant(count_t n, double *Gammas) { if (n < 2) { return Gammas; } @@ -17,7 +17,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) { list_t *pos = L; - for (uint64_t i = 1; i < n; i++) { + for (count_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; @@ -69,7 +69,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) { double *g = (double *)calloc(n + 1, sizeof(double)); double rho = 0; - for (uint64_t i = 0; i < n + 1; i++) { + for (count_t i = 0; i < n + 1; i++) { if (i > pos->next->p->x) { pos = pos->next; } @@ -87,7 +87,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) { } } - for (uint64_t i = 0; i < n + 1; i++) { + for (count_t i = 0; i < n + 1; i++) { g[i] += rho / 2; } diff --git a/lib/convex.h b/lib/convex.h index 29f4ae0..5a405d4 100644 --- a/lib/convex.h +++ b/lib/convex.h @@ -6,8 +6,10 @@ #include #include +#include "types.h" + typedef struct { - uint64_t x; + count_t x; double y; } point_t; @@ -17,5 +19,5 @@ typedef struct list_tag { point_t *p; } list_t; -double *get_convex_minorant(uint64_t n, double *Gammas); +double *get_convex_minorant(count_t n, double *Gammas); diff --git a/lib/graph.c b/lib/graph.c new file mode 100644 index 0000000..a1edf13 --- /dev/null +++ b/lib/graph.c @@ -0,0 +1,67 @@ + +#include "graph.h" + +graph_t *graph_create_square(D_t D, L_t L) { + v_t nv = pow(L, D); + v_t ne = D * nv; + + v_t *v_i = (v_t *)malloc((nv + 1) * sizeof(v_t)); + + for (v_t i = 0; i < nv + 1; i++) { + v_i[i] = 2 * D * i; + } + + v_t *v_adj = (v_t *)malloc(2 * D * nv * sizeof(v_t)); + + for (v_t i = 0; i < nv; i++) { + for (D_t j = 0; j < D; j++) { + v_adj[v_i[i] + 2 * j] = pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)); + v_adj[v_i[i] + 2 * j + 1] = pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)); + } + } + + graph_t *g = (graph_t *)malloc(sizeof(graph_t)); + + g->ne = ne; + g->nv = nv; + g->v_i = v_i; + g->v_adj = v_adj; + + return g; +} + +graph_t *graph_add_ext(const graph_t *G) { + graph_t *tG = (graph_t *)calloc(1, sizeof(graph_t)); + + tG->nv = G->nv + 1; + tG->ne = G->ne + G->nv; + + tG->v_i = (v_t *)malloc((tG->nv + 1) * sizeof(v_t)); + tG->v_adj = (v_t *)malloc(2 * tG->ne * sizeof(v_t)); + + for (v_t i = 0; i < G->nv + 1; i++) { + tG->v_i[i] = G->v_i[i] + i; + } + + tG->v_i[tG->nv] = 2 * tG->ne; + + for (v_t i = 0; i < G->nv; i++) { + v_t nn = G->v_i[i + 1] - G->v_i[i]; + + for (v_t j = 0; j < nn; j++) { + tG->v_adj[tG->v_i[i] + j] = G->v_adj[G->v_i[i] + j]; + } + + tG->v_adj[tG->v_i[i] + nn] = G->nv; + tG->v_adj[tG->v_i[G->nv] + i] = i; + } + + return tG; +} + +void graph_free(graph_t *g) { + free(g->v_i); + free(g->v_adj); + free(g); +} + diff --git a/lib/graph.h b/lib/graph.h new file mode 100644 index 0000000..eda7630 --- /dev/null +++ b/lib/graph.h @@ -0,0 +1,18 @@ + +#include +#include +#include + +#include "types.h" + +typedef struct { + v_t ne; + v_t nv; + v_t *v_i; + v_t *v_adj; +} graph_t; + +graph_t *graph_create_square(D_t D, L_t L); + +graph_t *graph_add_ext(const graph_t *G); + diff --git a/lib/queue.c b/lib/queue.c deleted file mode 100644 index 0823518..0000000 --- a/lib/queue.c +++ /dev/null @@ -1,50 +0,0 @@ - -#include "queue.h" - -void stack_push(ll_t **q, uint32_t x) { - ll_t *nq = malloc(sizeof(ll_t)); - nq->x = x; - nq->next = *q; - - *q = nq; -} - -void stack_push_d(dll_t **q, double x) { - dll_t *nq = malloc(sizeof(dll_t)); - nq->x = x; - nq->next = *q; - - *q = nq; -} - -uint32_t stack_pop(ll_t **q) { - ll_t *old_q = *q; - - *q = old_q->next; - uint32_t x = old_q->x; - - free(old_q); - - return x; -} - -double stack_pop_d(dll_t **q) { - dll_t *old_q = *q; - - *q = old_q->next; - double x = old_q->x; - - free(old_q); - - return x; -} - -bool stack_contains(const ll_t *q, uint32_t x) { - if (q == NULL) { - return false; - } else if (q->x == x) { - return true; - } else { - return stack_contains(q->next, x); - } -} diff --git a/lib/queue.h b/lib/queue.h deleted file mode 100644 index 23834ba..0000000 --- a/lib/queue.h +++ /dev/null @@ -1,25 +0,0 @@ - -#pragma once - -#include -#include -#include -#include - -typedef struct ll_tag { - uint32_t x; - struct ll_tag *next; -} ll_t; - -typedef struct dll_tag { - double x; - struct dll_tag *next; -} dll_t; - -void stack_push(ll_t **q, uint32_t x); -void stack_push_d(dll_t **q, double x); - -uint32_t stack_pop(ll_t **q); -double stack_pop_d(dll_t **q); - -bool stack_contains(const ll_t *q, uint32_t x); diff --git a/lib/rand.c b/lib/rand.c new file mode 100644 index 0000000..76f537d --- /dev/null +++ b/lib/rand.c @@ -0,0 +1,20 @@ + +#include +#include + +unsigned long int rand_seed() { + /* + returns a random unsigned long integer read from the standard unix + pseudorandom device /dev/urandom + */ + + FILE *f = fopen("/dev/urandom", "r"); + assert(f != NULL); + + unsigned long int seed; + fread(&seed, sizeof(unsigned long int), 1, f); + + fclose(f); + + return seed; +} diff --git a/lib/rand.h b/lib/rand.h new file mode 100644 index 0000000..2354f6a --- /dev/null +++ b/lib/rand.h @@ -0,0 +1,7 @@ + +#pragma once + +#include +#include + +unsigned long int rand_seed(); diff --git a/lib/stack.c b/lib/stack.c new file mode 100644 index 0000000..9fde1e6 --- /dev/null +++ b/lib/stack.c @@ -0,0 +1,41 @@ + +#include "stack.h" + +void stack_push(ll_t **q, v_t x) { + ll_t *nq = malloc(sizeof(ll_t)); + nq->x = x; + nq->next = *q; + + *q = nq; +} + +void stack_push_d(dll_t **q, double x) { + dll_t *nq = malloc(sizeof(dll_t)); + nq->x = x; + nq->next = *q; + + *q = nq; +} + +v_t stack_pop(ll_t **q) { + ll_t *old_q = *q; + + *q = old_q->next; + v_t x = old_q->x; + + free(old_q); + + return x; +} + +double stack_pop_d(dll_t **q) { + dll_t *old_q = *q; + + *q = old_q->next; + double x = old_q->x; + + free(old_q); + + return x; +} + diff --git a/lib/stack.h b/lib/stack.h new file mode 100644 index 0000000..a354ab5 --- /dev/null +++ b/lib/stack.h @@ -0,0 +1,26 @@ + +#pragma once + +#include +#include +#include +#include + +#include "types.h" + +typedef struct ll_tag { + v_t x; + struct ll_tag *next; +} ll_t; + +typedef struct dll_tag { + double x; + struct dll_tag *next; +} dll_t; + +void stack_push(ll_t **q, v_t x); +void stack_push_d(dll_t **q, double x); + +v_t stack_pop(ll_t **q); +double stack_pop_d(dll_t **q); + diff --git a/lib/tree.c b/lib/tree.c new file mode 100644 index 0000000..081f1d1 --- /dev/null +++ b/lib/tree.c @@ -0,0 +1,81 @@ + +#include "tree.h" + +node_t *tree_createNode(v_t value, v_t level, node_t *left, node_t *right) { + node_t *node = (node_t *)malloc(sizeof(node_t)); + + node->value = value; + node->level = level; + node->left = left; + node->right = right; + + return node; +} + +void tree_freeNode(node_t *node) { + if (node != NULL) { + tree_freeNode(node->left); + tree_freeNode(node->right); + free(node); + } +} + +void tree_skew(node_t **T) { + if (*T != NULL) { + if ((*T)->left != NULL) { + if ((*T)->left->level == (*T)->level) { + node_t *L = (*T)->left; + (*T)->left = L->right; + L->right = (*T); + + *T = L; + } + } + } +} + +void tree_split(node_t **T) { + if (*T != NULL) { + if ((*T)->right != NULL) { + if ((*T)->right->right != NULL) { + if ((*T)->level == (*T)->right->right->level) { + node_t *R = (*T)->right; + (*T)->right = R->left; + R->left = *T; + R->level = R->level + 1; + *T = R; + } + } + } + } +} + +void tree_insert(node_t **T, v_t x) { + if (*T == NULL) { + *T = tree_createNode(x, 1, NULL, NULL); + } else if (x < (*T)->value) { + node_t *L = (*T)->left; + tree_insert(&L, x); + (*T)->left = L; + } else if (x > (*T)->value) { + node_t *R = (*T)->right; + tree_insert(&R, x); + (*T)->right = R; + } + + tree_skew(T); + tree_split(T); +} + +bool tree_contains(node_t *T, v_t x) { + if (T == NULL) { + return false; + } else if (x < T->value) { + return tree_contains(T->left, x); + } else if (x > T->value) { + return tree_contains(T->right, x); + } else { + return true; + } +} + diff --git a/lib/tree.h b/lib/tree.h new file mode 100644 index 0000000..dc22c2d --- /dev/null +++ b/lib/tree.h @@ -0,0 +1,22 @@ + +#pragma once + +#include +#include +#include + +#include "types.h" + +typedef struct node_t { + v_t value; + v_t level; + struct node_t *left; + struct node_t *right; +} node_t; + +void tree_insert(node_t **T, v_t x); + +bool tree_contains(node_t *T, v_t x); + +void tree_freeNode(node_t *T); + diff --git a/lib/types.h b/lib/types.h new file mode 100644 index 0000000..982db47 --- /dev/null +++ b/lib/types.h @@ -0,0 +1,16 @@ + +#include + +typedef uint_fast32_t v_t; +typedef uint_fast8_t q_t; +typedef uint_fast8_t D_t; +typedef uint_fast16_t L_t; +typedef uint_fast64_t count_t; + +#define MAX_Q 256 + +#define PRID PRIu8 +#define PRIL PRIu16 +#define PRIcount PRIu64 +#define PRIq PRIu8 + diff --git a/lib/wolff.h b/lib/wolff.h index dfe7b1c..dae907a 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -13,33 +13,25 @@ #include #include -#include -#include - -#include "queue.h" +#include "types.h" +#include "rand.h" +#include "stack.h" #include "convex.h" - -typedef enum { - WOLFF, - WOLFF_GHOST, - METROPOLIS -} sim_t; +#include "graph.h" +#include "tree.h" typedef struct { graph_t *g; - bool *spins; - int32_t M; - double H; + q_t *spins; + double T; + double *H; + double T_prob; + double *H_probs; + double E; + v_t *M; + q_t q; } ising_state_t; -typedef struct { - uint32_t nv; - int32_t dJb; - int32_t dHb; - bool hit_ghost; - ll_t *spins; -} cluster_t; - typedef struct { uint64_t n; double x; @@ -60,21 +52,13 @@ typedef struct { double O2; } autocorr_t; -int8_t sign(double x); - -cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost, - gsl_rng *r); +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); -uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r, - double *ps); - void update_meas(meas_t *m, double x); void update_autocorr(autocorr_t *OO, double O); -double add_to_avg(double mx, double x, uint64_t n); - double rho(autocorr_t *o, uint64_t i); diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c index a417d4e..416a498 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -1,257 +1,94 @@ -#include "queue.h" #include "wolff.h" -int8_t spin_to_sign(bool spin) { - /* takes a spin (represented by a bool) and converts it to an integer (plus - * or minus one). our convention takes true to be negative spin and false to - * be positive spin - */ - if (spin) { - return -1; - } else { - return 1; - } -} - -int8_t sign(double x) { - // the sign function. zero returns positive sign. - return x >= 0 ? 1 : -1; -} - -graph_t *graph_add_ext(const graph_t *G) { - /* takes a graph object G and returns tG, the same graph with an extra vertex - * that is connected to every other vertex. - */ - graph_t *tG = (graph_t *)calloc(1, sizeof(graph_t)); - - tG->nv = G->nv + 1; - tG->ne = G->ne + G->nv; - - tG->ev = (uint32_t *)malloc(2 * tG->ne * sizeof(uint32_t)); - tG->vei = (uint32_t *)malloc((tG->nv + 1) * sizeof(uint32_t)); - tG->ve = (uint32_t *)malloc(2 * tG->ne * sizeof(uint32_t)); - tG->vx = (double *)malloc(2 * tG->nv * sizeof(double)); - tG->bq = (bool *)malloc(tG->nv * sizeof(bool)); - - memcpy(tG->ev, G->ev, 2 * G->ne * sizeof(uint32_t)); - memcpy(tG->vx, G->vx, 2 * G->nv * sizeof(double)); - memcpy(tG->bq, G->bq, G->nv * sizeof(bool)); - - tG->vx[2 * G->nv] = -1; - tG->vx[2 * G->nv + 1] = -0.5; - tG->bq[G->nv] = false; - - for (uint32_t i = 0; i < G->nv; i++) { - tG->ev[2 * G->ne + 2 * i] = i; - tG->ev[2 * G->ne + 2 * i + 1] = G->nv; - } - - for (uint32_t i = 0; i < G->nv; i++) { - tG->vei[i] = G->vei[i] + i; - - for (uint32_t j = 0; j < G->vei[i + 1] - G->vei[i]; j++) { - tG->ve[tG->vei[i] + j] = G->ve[G->vei[i] + j]; - } - - tG->ve[tG->vei[i] + G->vei[i + 1] - G->vei[i]] = G->ne + i; - } - - tG->vei[G->nv] = G->vei[G->nv] + G->nv; - tG->vei[G->nv + 1] = tG->vei[G->nv] + G->nv; - - for (uint32_t i = 0; i < G->nv; i++) { - tG->ve[tG->vei[G->nv] + i] = G->ne + i; - } - - return tG; -} - -uint32_t get_neighbor(const graph_t *g, uint32_t v, uint32_t i) { - // returns the index of the ith neighbor of vertex v on graph g. - assert(i < g->vei[v + 1] - g->vei[v]); // don't request a neighbor over the total number of neighbors! - - uint32_t e, v1, v2; - - e = g->ve[g->vei[v] + i]; // select the ith bond connected to site - v1 = g->ev[2 * e]; - v2 = g->ev[2 * e + 1]; - - return v == v1 ? v2 : v1; // distinguish neighboring site from site itself +q_t delta(q_t s0, q_t s1) { + return s0 == s1 ? 1 : 0; } -cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *s, bool stop_on_ghost, - gsl_rng *r) { - /* flips a wolff cluster on the graph g with initial state s. s is always - * changed by the routine. the probability of adding a normal bond to the - * cluster is passed by ps[0], and the probability of adding a bond with the - * external field spin to the cluster is passed by ps[1]. if stop_on_ghost is - * true, adding the external field spin to the cluster stops execution of the - * routine. flip_cluster returns an object of type cluster_t, which encodes - * information about the flipped cluster. - */ - uint32_t v0; - int32_t n_h_bonds, n_bonds; - bool s0; - cluster_t *c; - - v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex - s0 = s[v0]; // record its orientation +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 stack_push(&stack, v0); // push the initial vertex to the stack - // initiate the data structure for returning flip information - c = (cluster_t *)calloc(1, sizeof(cluster_t)); + node_t *T = NULL; while (stack != NULL) { - uint32_t v; - uint32_t nn; - - v = stack_pop(&stack); - nn = g->vei[v + 1] - g->vei[v]; + v_t v = stack_pop(&stack); - if (s[v] == s0 && !(c->hit_ghost && stop_on_ghost)) { // if the vertex hasn't already been flipped - s[v] = !s[v]; // flip the vertex - if (stop_on_ghost) { - stack_push(&(c->spins), v); - } - - for (uint32_t i = 0; i < nn; i++) { - bool is_ext; - uint32_t vn; - int32_t *bond_counter; - double prob; + if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped + q_t s_old = s->spins[v]; + q_t s_new = (s->spins[v] + step) % s->q; - vn = get_neighbor(g, v, i); + s->spins[v] = s_new; // flip the vertex + tree_insert(&T, v); - is_ext = (v == g->nv - 1 || vn == g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph + v_t nn = s->g->v_i[v + 1] - s->g->v_i[v]; - bond_counter = is_ext ? &(c->dHb) : &(c->dJb); - prob = is_ext ? ps[1] : ps[0]; + 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]; + double prob; - if (s[vn] == - s0) { // if the neighboring site matches the flipping cluster... - (*bond_counter)++; + bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1); - if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... - if (is_ext && stop_on_ghost) { - c->hit_ghost = true; - } - stack_push(&stack, vn); // push the neighboring vertex to the stack + if (is_ext) { + q_t M_ind_0; + q_t M_ind_1; + if (vn == s->g->nv - 1) { + M_ind_0 = (s_old + s->q - s->spins[s->g->nv - 1]) % s->q; + M_ind_1 = (s_new + s->q - s->spins[s->g->nv - 1]) % s->q; + } else { + M_ind_0 = (sn + s->q - s_old) % s->q; + M_ind_1 = (sn + s->q - s_new) % s->q; } + 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]; } else { - (*bond_counter)--; + prob = s->T_prob * delta(s0, sn); + s->E += - ((double)delta(s->spins[v], sn)) + ((double)delta(s0, sn)); } - } - - if (v != g->nv - 1) { // count the number of non-external sites that flip - c->nv++; - } - } - } - - return c; -} - -uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r, - double *ps) { - uint32_t n_flips; - bool no_r, no_ps; - no_r = false; - no_ps = false; - - if (r == NULL) { - no_r = true; - r = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(r, jst_rand_seed()); - } - - if (ps == NULL) { // computing exponentials is relatively expensive, so calling functions are given the option to supply values calculated once in the entire runtime - no_ps = true; - ps = (double *)malloc(2 * sizeof(double)); - ps[0] = 1 - exp(-2 / T); - ps[1] = 1 - exp(-2 * fabs(H) / T); - } - - switch (sim) { - case METROPOLIS: { - uint32_t v0 = gsl_rng_uniform_int(r, s->g->nv); // pick a random vertex - uint32_t nn = s->g->vei[v0 + 1] - s->g->vei[v0]; - - double dE = 0; - for (uint32_t i = 0; i < nn; i++) { - bool is_ext; - uint32_t vn; - int32_t *bond_counter; - double prob; - vn = get_neighbor(s->g, v0, i); - - is_ext = (v0 == s->g->nv - 1 || vn == s->g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph - - if (is_ext) { - dE += 2 * H * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]); - } else { - dE += 2 * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]); + if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... + stack_push(&stack, vn); // push the neighboring vertex to the stack } } - double p = exp(-dE / T); - if (gsl_rng_uniform(r) < p) { - s->M += - sign(H) * 2 * spin_to_sign(s->spins[v0]) * spin_to_sign(s->spins[s->g->nv - 1]); - s->H += dE; - s->spins[v0] = !s->spins[v0]; + if (v != s->g->nv - 1) { // count the number of non-external sites that flip + nv++; } + } + } - n_flips = 1; - } break; - - case WOLFF: { - cluster_t *c = flip_cluster(s->g, ps, s->spins, false, r); - s->M += - sign(H) * 2 * c->dHb; - s->H += 2 * (c->dJb + sign (H) * H * c->dHb); - n_flips = c->nv; + tree_freeNode(T); - free(c); - } break; - case WOLFF_GHOST: { - cluster_t *c = flip_cluster(s->g, ps, s->spins, true, r); + return nv; +} - if (c->hit_ghost) { - while (c->spins != NULL) { - uint32_t v = stack_pop(&(c->spins)); - s->spins[v] = !s->spins[v]; - // if we hit the external spin, undo the cluster flip - } - } else { - while (c->spins != NULL) { - stack_pop(&(c->spins)); - // we have to clear the memory on the stack anyway... - } - s->M += - sign(H) * 2 * c->dHb; - s->H += 2 * (c->dJb + sign (H) * H * c->dHb); - } +double add_to_avg(double mx, double x, count_t n) { + return mx * (n / (n + 1.0)) + x * 1.0 / (n + 1.0); +} - n_flips = c->nv; +void update_meas(meas_t *m, double x) { + count_t n = m->n; - free(c); - } break; - } + m->x = add_to_avg(m->x, x, n); + m->x2 = add_to_avg(m->x2, pow(x, 2), n); - if (no_ps) { - free(ps); - } + 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 (no_r) { - gsl_rng_free(r); + 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); } - return n_flips; -} - -double add_to_avg(double mx, double x, uint64_t n) { - return mx * (n / (n + 1.)) + x * 1. / (n + 1.); + (m->n)++; } void update_autocorr(autocorr_t *OO, double O) { @@ -260,7 +97,7 @@ void update_autocorr(autocorr_t *OO, double O) { dll_t *Otmp = OO->Op; dll_t *Osave; - uint64_t t = 0; + count_t t = 0; while (Otmp != NULL) { OO->OO[t] = add_to_avg(OO->OO[t], O * (Otmp->x), OO->n - t - 1); @@ -286,25 +123,6 @@ void update_autocorr(autocorr_t *OO, double O) { OO->n++; } -double rho(autocorr_t *o, uint64_t i) { +double rho(autocorr_t *o, count_t i) { return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2)); } - -void update_meas(meas_t *m, double x) { - uint64_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)++; -} -- cgit v1.2.3-70-g09d2