From fbbc4d9655835c0d6fdf58f231e59d7007a99407 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 23 Apr 2018 18:32:49 -0400
Subject: lots of changes

---
 lib/cluster.h     |  8 ++++++++
 lib/convex.c      |  1 -
 lib/measurement.h |  2 ++
 lib/yule_walker.c | 41 +++++++++++++++++++++++++++++++++++++++++
 lib/yule_walker.h | 12 ++++++++++++
 5 files changed, 63 insertions(+), 1 deletion(-)
 create mode 100644 lib/yule_walker.c
 create mode 100644 lib/yule_walker.h

(limited to 'lib')

diff --git a/lib/cluster.h b/lib/cluster.h
index 2de17e5..d118735 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -22,6 +22,7 @@
 #include "orthogonal.h"
 #include "dihedral.h"
 #include "dihinf.h"
+#include "yule_walker.h"
 
 typedef struct {
   graph_t *g;
@@ -62,6 +63,13 @@ typedef struct {
   q_t n;
 } vector_state_t;
 
+typedef enum {
+  VECTOR,
+  MODULATED,
+  CUBIC,
+  QUADRATIC
+} vector_field_t;
+
 v_t flip_cluster(ising_state_t *s, v_t v0, q_t s1, gsl_rng *r);
 
 v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r);
diff --git a/lib/convex.c b/lib/convex.c
index df506f9..7255ad4 100644
--- a/lib/convex.c
+++ b/lib/convex.c
@@ -93,7 +93,6 @@ double *get_convex_minorant(count_t n, double *Gammas) {
 
   while (L != NULL) {
     free(L->p);
-    list_t *L_save = L;
     L = L->next;
     free(L);
   }
diff --git a/lib/measurement.h b/lib/measurement.h
index a24df9f..eaa260b 100644
--- a/lib/measurement.h
+++ b/lib/measurement.h
@@ -1,4 +1,6 @@
 
+#pragma once
+
 #include <math.h>
 #include <stdlib.h>
 
diff --git a/lib/yule_walker.c b/lib/yule_walker.c
new file mode 100644
index 0000000..77b3e96
--- /dev/null
+++ b/lib/yule_walker.c
@@ -0,0 +1,41 @@
+
+#include "yule_walker.h"
+
+double yule_walker(const autocorr_t *a) {
+  gsl_vector *rhos = gsl_vector_alloc(a->W);
+  gsl_matrix *mat = gsl_matrix_alloc(a->W, a->W);
+  gsl_vector *phis = gsl_vector_alloc(a->W);
+
+  for (count_t i = 0; i < a->W; i++) {
+    gsl_vector_set(rhos, i, rho(a, i));
+  }
+
+  for (count_t i = 0; i < a->W; i++) {
+    gsl_matrix_set(mat, i, i, 1.0);
+
+    for (count_t j = 0; j < i; j++) {
+      gsl_matrix_set(mat, i, j, gsl_vector_get(rhos, i - 1 - j));
+    }
+
+    for (count_t j = 0; j < a->W - 1 - i; j++) {
+      gsl_matrix_set(mat, i, i + 1 + j, gsl_vector_get(rhos, j));
+    }
+  }
+
+  int out = gsl_linalg_cholesky_solve(mat, rhos, phis);
+
+  double rhophi = 0;
+  double onephi = 0;
+
+  for (count_t i = 0; i < a->W; i++) {
+    rhophi += gsl_vector_get(rhos, i) * gsl_vector_get(phis, i);
+    onephi += gsl_vector_get(phis, i);
+  }
+
+  gsl_vector_free(rhos);
+  gsl_matrix_free(mat);
+  gsl_vector_free(phis);
+
+  return (1 - rhophi) / pow(1 - onephi, 2);
+}
+
diff --git a/lib/yule_walker.h b/lib/yule_walker.h
new file mode 100644
index 0000000..883ea03
--- /dev/null
+++ b/lib/yule_walker.h
@@ -0,0 +1,12 @@
+
+#pragma once
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+
+#include "types.h"
+#include "measurement.h"
+
+double yule_walker(const autocorr_t *a);
+
-- 
cgit v1.2.3-70-g09d2