summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/cluster.c24
-rw-r--r--lib/measurement.c20
-rw-r--r--lib/measurement.h13
-rw-r--r--lib/orthogonal.c10
-rw-r--r--lib/orthogonal.h6
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);