From a3a1a4029bf52b13376e5caeae2f6ad3966a6338 2022-07-12 22:32:12
From: Roberto Flores <rflores@igp.gob.pe>
Date: 2022-07-12 22:32:12
Subject: [PATCH] HS_outliers in C

---

diff --git a/schainc/_HS_algorithm.c b/schainc/_HS_algorithm.c
new file mode 100644
index 0000000..31bab54
--- /dev/null
+++ b/schainc/_HS_algorithm.c
@@ -0,0 +1,86 @@
+#include <Python.h>
+#include <numpy/arrayobject.h>
+#include <math.h>
+
+
+static PyObject *HS_algorithm(PyObject *self, PyObject *args) {
+  double navg;
+  PyObject *data_obj, *data_array;
+
+  if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
+      return NULL;
+  }
+
+  data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
+
+  if (data_array == NULL) {
+      Py_XDECREF(data_array);
+      Py_XDECREF(data_obj);
+      return NULL;
+  }
+  double *sortdata = (double*)PyArray_DATA(data_array);
+  int lenOfData = (int)PyArray_SIZE(data_array) ;
+  double nums_min = lenOfData*0.75;
+  if (nums_min <= 5) nums_min = 5;
+  double sump = 0;
+  double sumq = 0;
+  long j = 0;
+  int cont = 1;
+  double rtest = 0;
+  while ((cont == 1) && (j < lenOfData)) {
+    sump = sump + sortdata[j];
+    sumq = sumq + pow(sortdata[j], 2);
+    if (j > nums_min) {
+      rtest = (double)j/(j-1) + 1/navg;
+      //printf("%ld\n", j);
+      //printf("%f \n", rtest);
+      //printf("%f \n", sump);
+      if ((sumq*j) > (rtest*pow(sump, 2))) {
+        j = j - 1;
+        sump = sump - sortdata[j];
+        sumq = sumq - pow(sortdata[j],2);
+        cont = 0;
+      }
+    }
+    //printf("%ld\n", j);
+    j = j + 1;
+  }
+
+  double lnoise = sump / j;
+
+  Py_DECREF(data_array);
+
+  // return PyLong_FromLong(lnoise);
+  return PyFloat_FromDouble(j);
+}
+
+
+static PyMethodDef noiseMethods[] = {
+  { "HS_algorithm", HS_algorithm, METH_VARARGS, "Applies hildebrand_sekhon algorithm" },
+  { NULL, NULL, 0, NULL }
+};
+
+#if PY_MAJOR_VERSION >= 3
+
+static struct PyModuleDef noisemodule = {
+  PyModuleDef_HEAD_INIT,
+  "_HS_algorithm",
+  "Applies hildebrand_sekhon algorithm",
+  -1,
+  noiseMethods
+};
+
+#endif
+
+#if PY_MAJOR_VERSION >= 3
+  PyMODINIT_FUNC PyInit__HS_algorithm(void) {
+    Py_Initialize();
+    import_array();
+    return PyModule_Create(&noisemodule);
+  }
+#else
+  PyMODINIT_FUNC init_HS_algorithm() {
+    Py_InitModule("_HS_algorithm", noiseMethods);
+    import_array();
+  }
+#endif