From f3d53bac095fd47f7c10f54007ffa11ab0f1ac1f Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 2 Nov 2017 23:00:48 -0400
Subject: added new method for computing autocorrelation

---
 lib/convex.c | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 93 insertions(+)
 create mode 100644 lib/convex.c

(limited to 'lib/convex.c')

diff --git a/lib/convex.c b/lib/convex.c
new file mode 100644
index 0000000..33af1b4
--- /dev/null
+++ b/lib/convex.c
@@ -0,0 +1,93 @@
+
+#include "convex.h"
+
+double slope(point_t *P, point_t *Q) {
+  return (Q->y - P->y) / ((double)(Q->x) - (double)(P->x));
+}
+
+double *get_convex_minorant(uint64_t n, double *Gammas) {
+  list_t *L = (list_t *)calloc(1, sizeof(list_t));
+  L->p = (point_t *)calloc(1, sizeof(point_t));
+  L->p->x = 0;
+  L->p->y = Gammas[0];
+
+  list_t *pos = L;
+
+  for (uint64_t i = 1; i < n; i++) {
+    pos->next = (list_t *)calloc(1, sizeof(list_t));
+    pos->next->p = (point_t *)calloc(1, sizeof(point_t));
+    pos->next->p->x = i;
+    pos->next->p->y = Gammas[i];
+    pos->next->prev = pos;
+    pos = pos->next;
+  }
+
+  pos->next = (list_t *)calloc(1, sizeof(list_t));
+  pos->next->p = (point_t *)calloc(1, sizeof(point_t));
+  pos->next->p->x = n;
+  pos->next->p->y = 0;
+  pos->next->prev = pos;
+
+  list_t *X = L;
+  list_t *Y = L->next;
+  list_t *Z = Y->next;
+
+  while (true) {
+    if (slope(X->p, Y->p) <= slope(Y->p, Z->p)) {
+      X = Y;
+      Y = Z;
+      if (Z->next == NULL) {
+        break;
+      } else {
+        Z = Z->next;
+      }
+    } else {
+      Y->prev->next = Y->next;
+      Y->next->prev = Y->prev;
+      free(Y->p);
+      free(Y);
+      if (X->prev != NULL) {
+        Y = X;
+        X = X->prev;
+      } else {
+        if (Z->next != NULL) {
+          Y = Z;
+          Z = Z->next;
+        } else {
+          break;
+        }
+      }
+    }
+  }
+
+  pos = L;
+
+  double *g = (double *)calloc(n, sizeof(double));
+  double rho = 0;
+
+  for (uint64_t i = 0; i < n; i++) {
+    if (i > pos->next->p->x) {
+      pos = pos->next;
+    }
+
+    g[i] = pos->p->y + ((double)i - (double)(pos->p->x)) *  (pos->next->p->y - pos->p->y) / ((double)(pos->next->p->x) - (double)(pos->p->x));
+
+    if (Gammas[i] - g[i] > rho) {
+      rho = Gammas[i] - g[i];
+    }
+  }
+
+  for (uint64_t i = 0; i < n; i++) {
+    g[i] += rho / 2;
+  }
+
+  while (L != NULL) {
+    free(L->p);
+    list_t *L_save = L;
+    L = L->next;
+    free(L);
+  }
+
+  return g;
+}
+
-- 
cgit v1.2.3-70-g09d2