diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/cluster.c | 24 | ||||
| -rw-r--r-- | lib/measurement.c | 20 | ||||
| -rw-r--r-- | lib/measurement.h | 13 | ||||
| -rw-r--r-- | lib/orthogonal.c | 10 | ||||
| -rw-r--r-- | lib/orthogonal.h | 6 |
5 files changed, 52 insertions, 21 deletions
diff --git a/lib/cluster.c b/lib/cluster.c index fa417fa..7274eb9 100644 --- a/lib/cluster.c +++ b/lib/cluster.c @@ -193,8 +193,13 @@ v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) { // if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped if (!marks[v]) { + bool v_is_external = false; double *s_old, *s_new, *R_tmp; + if (v == s->g->nv - 1) { + v_is_external = true; + } + //tree_insert(&T, v); marks[v] = true; @@ -209,17 +214,24 @@ v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) { for (v_t i = 0; i < nn; i++) { v_t vn = s->g->v_adj[s->g->v_i[v] + i]; + + bool vn_is_external = false; + + if (vn == s->g->nv - 1) { + vn_is_external = true; + } + double *sn; - if (vn != s->g->nv - 1) { + + if (!vn_is_external) { sn = &(s->spins[s->n * vn]); } - double prob; - bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1); + double prob; - if (is_ext) { + if (v_is_external || vn_is_external) { double *rs_old, *rs_new; - if (vn == s->g->nv - 1) { + if (vn_is_external) { rs_old = vector_rotate_inverse(s->n, s->R, s_old); rs_new = vector_rotate_inverse(s->n, s->R, s_new); } else { @@ -237,8 +249,6 @@ v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) { } else { double dE = (s->J)(vector_dot(s->n, sn, s_old)) - (s->J)(vector_dot(s->n, sn, s_new)); prob = 1.0 - exp(-dE / s->T); - //printf("(%g %g) (%g %g) (%g %g) %g\n", s_old[0], s_old[1], s_new[0], s_new[1], sn[0], sn[1], dE); - //getchar(); s->E += dE; } diff --git a/lib/measurement.c b/lib/measurement.c index 435c749..ad824f6 100644 --- a/lib/measurement.c +++ b/lib/measurement.c @@ -2,10 +2,10 @@ #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); + return mx * (n / (n + 1.0)) + x / (n + 1.0); } -void update_meas(meas_t *m, double x) { +void meas_update(meas_t *m, double x) { count_t n = m->n; m->x = add_to_avg(m->x, x, n); @@ -14,16 +14,30 @@ void update_meas(meas_t *m, double x) { 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)++; } +double meas_dx(const meas_t *m) { + return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2))); +} + +double meas_c(const meas_t *m) { + return m->n / (m->n - 1.) * (m->x2 - pow(m->x, 2)); +} + +double meas_dc(const meas_t *m) { + return sqrt((m->m4 - (m->n - 3.)/(m->n - 1.) * pow(m->m2, 2)) / 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); @@ -56,7 +70,7 @@ void update_autocorr(autocorr_t *OO, double O) { OO->n++; } -double rho(autocorr_t *o, count_t i) { +double rho(const 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 index 4ced4dc..a24df9f 100644 --- a/lib/measurement.h +++ b/lib/measurement.h @@ -8,12 +8,9 @@ typedef struct { uint64_t n; double x; - double dx; double x2; double m2; double m4; - double c; - double dc; } meas_t; typedef struct { @@ -25,9 +22,15 @@ typedef struct { double O2; } autocorr_t; -void update_meas(meas_t *m, double x); +void meas_update(meas_t *m, double x); + +double meas_dx(const meas_t *m); + +double meas_c(const meas_t *m); + +double meas_dc(const meas_t *m); void update_autocorr(autocorr_t *OO, double O); -double rho(autocorr_t *o, uint64_t i); +double rho(const autocorr_t *o, uint64_t i); diff --git a/lib/orthogonal.c b/lib/orthogonal.c index 7e0a71b..87569ae 100644 --- a/lib/orthogonal.c +++ b/lib/orthogonal.c @@ -2,24 +2,28 @@ #include "orthogonal.h" void vector_replace(q_t n, double *v1, const double *v2) { + // writes vector v2 of length n to memory located at v1 for (q_t i = 0; i < n; i++) { v1[i] = v2[i]; } } void vector_add(q_t n, double *v1, const double *v2) { + // adds vector v2 of length n to vector v1 for (q_t i = 0; i < n; i++) { v1[i] += v2[i]; } } void vector_subtract(q_t n, double *v1, const double *v2) { + // subtracts vector v2 of length n from vector v1 for (q_t i = 0; i < n; i++) { v1[i] -= v2[i]; } } -double *vector_rotate(q_t n, double *rot, double *vec) { +double *vector_rotate(q_t n, const double *rot, const double *vec) { + // multiplies n by n rotation matrix rot to vector vec double *rot_vec = (double *)malloc(n * sizeof(double)); double prod = 0.0; @@ -46,7 +50,7 @@ double *vector_rotate_inverse(q_t n, const double *rot, const double *vec) { return rot_vec; } -double vector_dot(q_t n, double *v1, double *v2) { +double vector_dot(q_t n, const double *v1, const double *v2) { double dot = 0; for (q_t i = 0; i < n; i++) { @@ -56,7 +60,7 @@ double vector_dot(q_t n, double *v1, double *v2) { return dot; } -double *orthogonal_rotate(q_t n, double *r, double *m) { +double *orthogonal_rotate(q_t n, const double *r, const double *m) { double *mul = (double *)calloc(n * n, sizeof(double)); for (q_t i = 0; i < n; i++) { diff --git a/lib/orthogonal.h b/lib/orthogonal.h index a763b08..60d5f49 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -12,13 +12,13 @@ void vector_add(q_t n, double *v1, const double *v2); void vector_subtract(q_t n, double *v1, const double *v2); -double *vector_rotate(q_t n, double *rot, double *vec); +double *vector_rotate(q_t n, const double *rot, const double *vec); double *vector_rotate_inverse(q_t n, const double *rot, const double *vec); -double vector_dot(q_t n, double *v1, double *v2); +double vector_dot(q_t n, const double *v1, const double *v2); -double *orthogonal_rotate(q_t n, double *m1, double *m2); +double *orthogonal_rotate(q_t n, const double *m1, const double *m2); double *gen_rot(gsl_rng *r, q_t n); |