split up some files

4 files changed, 180 insertions, 284 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index 29dd0cb..a20912e 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -20,6 +20,7 @@
 #include "graph.h"
 #include "tree.h"
 #include "measurement.h"
+#include "vector.h"
 #include "orthogonal.h"
 #include "dihedral.h"
 #include "dihinf.h"
@@ -95,182 +96,6 @@ class state_t {
     }
 };
 
-template <q_t q, class T>
-struct vector_t { T *x; };
-
-template <q_t q, class T>
-void init(vector_t <q, T> *ptr) {
-  ptr->x = (T *)calloc(q, sizeof(T));
-
-  ptr->x[0] = (T)1;
-}
-
-template <q_t q, class T>
-vector_t <q, T> copy (vector_t <q, T> v) {
-  vector_t <q, T> v_copy;
- 
- v_copy.x = (T *)calloc(q, sizeof(T));
-
-  for (q_t i = 0; i < q; i++) {
-    v_copy.x[i] = v.x[i];
-  }
-
-  return v_copy;
-}
-
-template <q_t q, class T>
-void add (vector_t <q, T> v1, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1.x[i] += v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-void subtract (vector_t <q, T> v1, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1.x[i] -= v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
-  vector_t <q, T> multiple;
-  multiple.x = (T *)malloc(q * sizeof(T));
-  for (q_t i = 0; i < q; i++) {
-    multiple.x[i] = a * v.x[i];
-  }
-
-  return multiple;
-}
-
-template <q_t q, class T>
-T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
-  T prod = 0;
-
-  for (q_t i = 0; i < q; i++) {
-    prod += v1.x[i] * v2.x[i];
-  }
-
-  return prod;
-}
-
-template <q_t q, class T>
-void free_spin (vector_t <q, T> v) {
-  free(v.x);
-}
-
-template <q_t q, class T>
-struct orthogonal_t { T *x; };
-
-template <q_t q, class T>
-void init(orthogonal_t <q, T> *ptr) {
-  ptr->x = (T *)calloc(q * q, sizeof(T));
-
-  for (q_t i = 0; i < q; i++) {
-    ptr->x[q * i + i] = (T)1;
-  }
-}
-
-template <q_t q, class T>
-orthogonal_t <q, T> copy (orthogonal_t <q, T> m) {
-  orthogonal_t <q, T> m_copy;
-  m_copy.x = (T *)calloc(q * q, sizeof(T));
-
-  for (q_t i = 0; i < q * q; i++) {
-    m_copy.x[i] = m.x[i];
-  }
-
-  return m_copy;
-}
-
-template <q_t q, class T>
-void free_spin (orthogonal_t <q, T> m) {
-  free(m.x);
-}
-
-template <q_t q, class T>
-vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) {
-  vector_t <q, T> v_rot;
-  v_rot.x = (T *)calloc(q, sizeof(T));
-
-  for (q_t i = 0; i < q; i++) {
-    for (q_t j = 0; j < q; j++) {
-      v_rot.x[i] += m.x[q * i + j] * v.x[j];
-    }
-  }
-
-  return v_rot;
-}
-
-
-template <q_t q, class T>
-orthogonal_t <q, T> act (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
-  orthogonal_t <q, T> m2_rot;
-  m2_rot.x = (T *)calloc(q * q, sizeof(T));
-
-  for (q_t i = 0; i < q; i++) {
-    for (q_t j = 0; j < q; j++) {
-      for (q_t k = 0; k < q; k++) {
-        m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k];
-      }
-    }
-  }
-
-  return m2_rot;
-}
-
-template <q_t q, class T>
-vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) {
-  vector_t <q, T> v_rot;
-  v_rot.x = (T *)calloc(q, sizeof(T));
-
-  for (q_t i = 0; i < q; i++) {
-    for (q_t j = 0; j < q; j++) {
-      v_rot.x[i] += m.x[q * j + i] * v.x[j];
-    }
-  }
-
-  return v_rot;
-}
-
-template <q_t q, class T>
-orthogonal_t <q, T> act_inverse (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
-  orthogonal_t <q, T> m2_rot;
-  m2_rot.x = (T *)calloc(q * q, sizeof(T));
-
-  for (q_t i = 0; i < q; i++) {
-    for (q_t j = 0; j < q; j++) {
-      for (q_t k = 0; k < q; k++) {
-        m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k];
-      }
-    }
-  }
-
-  return m2_rot;
-}
-
-template <q_t q>
-void generate_rotation (gsl_rng *r, orthogonal_t <q, double> *ptr) {
-  double *v = (double *)malloc(q * sizeof(double));
-  double v2 = 0;
-
-  for (q_t i = 0; i < q; i++) {
-    v[i] = gsl_ran_ugaussian(r);
-    v2 += v[i] * v[i];
-  }
-
-  ptr->x = (double *)calloc(q * q, sizeof(double));
-  
-  for (q_t i = 0; i < q; i++) {
-    ptr->x[q * i + i] = 1.0;
-    for (q_t j = 0; j < q; j++) {
-      ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2;
-    }
-  }
-
-  free(v);
-}
-
 template <class R_t, class X_t>
 v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
   v_t nv = 0;
diff --git a/lib/orthogonal.c b/lib/orthogonal.c
deleted file mode 100644
index 87569ae..0000000
--- a/lib/orthogonal.c
+++ /dev/null
@@ -1,99 +0,0 @@
-
-#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, 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;
-  for (q_t i = 0; i < n; i++) {
-    prod += rot[i] * vec[i];
-  }
-
-  for (q_t i = 0; i < n; i++) {
-    rot_vec[i] = vec[i] - 2 * prod * rot[i];
-  }
-
-  return rot_vec;
-}
-
-double *vector_rotate_inverse(q_t n, const double *rot, const double *vec) {
-  double *rot_vec = (double *)calloc(n, sizeof(double));
-
-  for (q_t i = 0; i < n; i++) {
-    for (q_t j = 0; j < n; j++) {
-      rot_vec[i] += rot[n * j + i] * vec[j];
-    }
-  }
-
-  return rot_vec;
-}
-
-double vector_dot(q_t n, const double *v1, const double *v2) {
-  double dot = 0;
-
-  for (q_t i = 0; i < n; i++) {
-    dot += v1[i] * v2[i];
-  }
-
-  return dot;
-}
-
-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++) {
-    double akOki = 0;
-
-    for (q_t k = 0; k < n; k++) {
-      akOki += r[k] * m[n * k + i];
-    }
-
-    for (q_t j = 0; j < n; j++) {
-      mul[n * j + i] = m[n * j + i] - 2 * akOki * r[j];
-    }
-  }
-
-  return mul;
-}
-
-double *gen_rot(gsl_rng *r, q_t n) {
-  double *v = (double *)malloc(n * sizeof(double));
-
-  double v2 = 0;
-
-  for (q_t i = 0; i < n; i++) {
-    v[i] = gsl_ran_ugaussian(r);
-    v2 += v[i] * v[i];
-  }
-
-  double magv = sqrt(v2);
-
-  for (q_t i = 0; i < n; i++) {
-    v[i] /= magv;
-  }
-
-  return v;
-}
-
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 60d5f49..0b2fdd5 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -1,24 +1,122 @@
 
+#pragma once
+
 #include <stdlib.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
-#include <math.h>
+#include <cmath>
 
 #include "types.h"
 
-void vector_replace(q_t n, double *v1, const double *v2);
+template <q_t q, class T>
+struct orthogonal_t { T *x; };
+
+template <q_t q, class T>
+void init(orthogonal_t <q, T> *ptr) {
+  ptr->x = (T *)calloc(q * q, sizeof(T));
+
+  for (q_t i = 0; i < q; i++) {
+    ptr->x[q * i + i] = (T)1;
+  }
+}
+
+template <q_t q, class T>
+orthogonal_t <q, T> copy (orthogonal_t <q, T> m) {
+  orthogonal_t <q, T> m_copy;
+  m_copy.x = (T *)calloc(q * q, sizeof(T));
+
+  for (q_t i = 0; i < q * q; i++) {
+    m_copy.x[i] = m.x[i];
+  }
+
+  return m_copy;
+}
+
+template <q_t q, class T>
+void free_spin (orthogonal_t <q, T> m) {
+  free(m.x);
+}
+
+template <q_t q, class T>
+vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) {
+  vector_t <q, T> v_rot;
+  v_rot.x = (T *)calloc(q, sizeof(T));
+
+  for (q_t i = 0; i < q; i++) {
+    for (q_t j = 0; j < q; j++) {
+      v_rot.x[i] += m.x[q * i + j] * v.x[j];
+    }
+  }
+
+  return v_rot;
+}
+
+template <q_t q, class T>
+orthogonal_t <q, T> act (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
+  orthogonal_t <q, T> m2_rot;
+  m2_rot.x = (T *)calloc(q * q, sizeof(T));
+
+  for (q_t i = 0; i < q; i++) {
+    for (q_t j = 0; j < q; j++) {
+      for (q_t k = 0; k < q; k++) {
+        m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k];
+      }
+    }
+  }
+
+  return m2_rot;
+}
+
+template <q_t q, class T>
+vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) {
+  vector_t <q, T> v_rot;
+  v_rot.x = (T *)calloc(q, sizeof(T));
+
+  for (q_t i = 0; i < q; i++) {
+    for (q_t j = 0; j < q; j++) {
+      v_rot.x[i] += m.x[q * j + i] * v.x[j];
+    }
+  }
+
+  return v_rot;
+}
+
+template <q_t q, class T>
+orthogonal_t <q, T> act_inverse (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
+  orthogonal_t <q, T> m2_rot;
+  m2_rot.x = (T *)calloc(q * q, sizeof(T));
 
-void vector_add(q_t n, double *v1, const double *v2);
+  for (q_t i = 0; i < q; i++) {
+    for (q_t j = 0; j < q; j++) {
+      for (q_t k = 0; k < q; k++) {
+        m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k];
+      }
+    }
+  }
 
-void vector_subtract(q_t n, double *v1, const double *v2);
+  return m2_rot;
+}
 
-double *vector_rotate(q_t n, const double *rot, const double *vec);
+template <q_t q>
+void generate_rotation (gsl_rng *r, orthogonal_t <q, double> *ptr) {
+  double *v = (double *)malloc(q * sizeof(double));
+  double v2 = 0;
 
-double *vector_rotate_inverse(q_t n, const double *rot, const double *vec);
+  for (q_t i = 0; i < q; i++) {
+    v[i] = gsl_ran_ugaussian(r);
+    v2 += v[i] * v[i];
+  }
 
-double vector_dot(q_t n, const double *v1, const double *v2);
+  ptr->x = (double *)calloc(q * q, sizeof(double));
+  
+  for (q_t i = 0; i < q; i++) {
+    ptr->x[q * i + i] = 1.0;
+    for (q_t j = 0; j < q; j++) {
+      ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2;
+    }
+  }
 
-double *orthogonal_rotate(q_t n, const double *m1, const double *m2);
+  free(v);
+}
 
-double *gen_rot(gsl_rng *r, q_t n);
 
diff --git a/lib/vector.h b/lib/vector.h
new file mode 100644
index 0000000..c7f459c
--- /dev/null
+++ b/lib/vector.h
@@ -0,0 +1,72 @@
+
+#pragma once
+
+#include <stdlib.h>
+#include <cmath>
+
+#include "types.h"
+
+template <q_t q, class T>
+struct vector_t { T *x; };
+
+template <q_t q, class T>
+void init(vector_t <q, T> *ptr) {
+  ptr->x = (T *)calloc(q, sizeof(T));
+
+  ptr->x[0] = (T)1;
+}
+
+template <q_t q, class T>
+vector_t <q, T> copy (vector_t <q, T> v) {
+  vector_t <q, T> v_copy;
+ 
+ v_copy.x = (T *)calloc(q, sizeof(T));
+
+  for (q_t i = 0; i < q; i++) {
+    v_copy.x[i] = v.x[i];
+  }
+
+  return v_copy;
+}
+
+template <q_t q, class T>
+void add (vector_t <q, T> v1, vector_t <q, T> v2) {
+  for (q_t i = 0; i < q; i++) {
+    v1.x[i] += v2.x[i];
+  }
+}
+
+template <q_t q, class T>
+void subtract (vector_t <q, T> v1, vector_t <q, T> v2) {
+  for (q_t i = 0; i < q; i++) {
+    v1.x[i] -= v2.x[i];
+  }
+}
+
+template <q_t q, class T>
+vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
+  vector_t <q, T> multiple;
+  multiple.x = (T *)malloc(q * sizeof(T));
+  for (q_t i = 0; i < q; i++) {
+    multiple.x[i] = a * v.x[i];
+  }
+
+  return multiple;
+}
+
+template <q_t q, class T>
+T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
+  T prod = 0;
+
+  for (q_t i = 0; i < q; i++) {
+    prod += v1.x[i] * v2.x[i];
+  }
+
+  return prod;
+}
+
+template <q_t q, class T>
+void free_spin (vector_t <q, T> v) {
+  free(v.x);
+}
+