working generator of rotations near the current state

2 files changed, 20 insertions, 8 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index ccd8725..5f770ad 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -58,7 +58,7 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
   ll_t *stack = NULL;     // create a new stack
   stack_push(&stack, v0); // push the initial vertex to the stack
 
-  bool *marks = (bool *)calloc(state->nv, sizeof(bool));
+  bool *marks = (bool *)calloc(state->nv + 1, sizeof(bool));
 
   while (stack != NULL) {
     v_t v = stack_pop(&stack);
@@ -67,10 +67,10 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
       X_t si_old, si_new;
       R_t R_old, R_new;
 
+      R_old = state->R;
       marks[v] = true;
 
       if (v == state->nv) {
-        R_old = state->R;
         R_new = act (r, R_old);
       } else {
         si_old = state->spins[v];
@@ -94,7 +94,7 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
 
         if (is_ext) {
           X_t rs_old, rs_new;
-          if (vn == state->g->nv - 1) {
+          if (vn == state->nv) {
             rs_old = act_inverse (R_old, si_old);
             rs_new = act_inverse (R_old, si_new);
           } else {
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 6166223..d8ad33d 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -167,9 +167,9 @@ orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, const state_t <o
 
 template <q_t q>
 orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const state_t <orthogonal_t <q, double>, vector_t <q, double>> *s, double epsilon) {
-  orthogonal_t <q, double> ptr;
+  orthogonal_t <q, double> m;
   vector_t <q, double> tmp_v;
-  ptr.is_reflection = true;
+  m.is_reflection = true;
 
   tmp_v.x = (double *)malloc(q * sizeof(double));
 
@@ -196,10 +196,22 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const state
     tmp_v.x[i] /= mag_v;
   }
 
-  vector_t <q, double> v = act(s->R, tmp_v);
+  vector_t <q, double> v = act_inverse(s->R, tmp_v);
   free(tmp_v.x);
-  ptr.x = v.x;
+  m.x = v.x;
 
-  return ptr;
+  v2 = 0;
+
+  for (q_t i = 0; i < q; i++) {
+    v2 += m.x[i] * m.x[i];
+  }
+
+  mag_v = sqrt(v2);
+
+  for (q_t i = 0; i < q; i++) {
+    m.x[i] /= mag_v;
+  }
+
+  return m;
 }