erange 4.0a dev release with integrated cistematic

[erange.git] / cistematic / cisstat / statextension.c

diff --git a/cistematic/cisstat/statextension.c b/cistematic/cisstat/statextension.c
new file mode 100644 (file)
index 0000000..a21703a
--- /dev/null
+++ b/cistematic/cisstat/statextension.c
@@ -0,0 +1,112 @@
+/*##########################################################################*/
+/*#                                                                         #*/
+/*# C O P Y R I G H T   N O T I C E                                         #*/
+/*#  Copyright (c) 2003-10 by:                                              #*/
+/*#    * California Institute of Technology                                 #*/
+/*#                                                                         #*/
+/*#    All Rights Reserved.                                                 #*/
+/*#                                                                         #*/
+/*# Permission is hereby granted, free of charge, to any person             #*/
+/*# obtaining a copy of this software and associated documentation files    #*/
+/*# (the "Software"), to deal in the Software without restriction,          #*/
+/*# including without limitation the rights to use, copy, modify, merge,    #*/
+/*# publish, distribute, sublicense, and/or sell copies of the Software,    #*/
+/*# and to permit persons to whom the Software is furnished to do so,       #*/
+/*# subject to the following conditions:                                    #*/
+/*#                                                                         #*/
+/*# The above copyright notice and this permission notice shall be          #*/
+/*# included in all copies or substantial portions of the Software.         #*/
+/*#                                                                         #*/
+/*# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #*/
+/*# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #*/
+/*# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #*/
+/*# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #*/
+/*# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #*/
+/*# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #*/
+/*# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #*/
+/*# SOFTWARE.                                                               #*/
+/*###########################################################################*/
+/*# C extension for commonly used stat methods.                             */
+
+#include <Python.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+static char module_doc[] = 
+"This module only implements the pearsonCorrelation in C for now.";
+
+static PyObject*
+the_func(PyObject *self, PyObject *args)
+{
+    PyObject *a, *b;
+    double c, numerator;
+    double *pwmA, meanA, denominatorA;
+    double *pwmB, meanB, denominatorB;
+    long listLen, index;
+    
+    if (!PyArg_UnpackTuple(args, "func", 2, 2, &a, &b)) {
+        return NULL;
+    }
+    
+    listLen = PyList_Size(a);
+    pwmA = malloc(listLen * sizeof(double));
+    for (index = 0; index < listLen; index++) {
+        pwmA[index] = PyFloat_AsDouble(PyList_GetItem(a, index));
+    }
+        
+    pwmB = malloc(PyList_Size(b) * sizeof(double));
+    for (index = 0; index < listLen; index++) {
+        pwmB[index] = PyFloat_AsDouble(PyList_GetItem(b, index));
+    }
+    
+    if (listLen > PyList_Size(b)) {
+        listLen = PyList_Size(b);
+    }
+    
+    meanA = 0.0;
+    meanB = 0.0;
+    
+    for (index = 0; index < listLen; index++) {
+        meanA += pwmA[index];
+        meanB += pwmB[index];
+    }
+    
+    meanA /= listLen;
+    meanB /= listLen;
+    
+    denominatorA = 0.0;
+    denominatorB = 0.0;
+    numerator = 0.0;
+    
+    for (index = 0; index < listLen; index++) {
+        numerator += (pwmA[index] - meanA) * (pwmB[index] - meanB);
+        denominatorA += (pwmA[index] - meanA) * (pwmA[index] - meanA);
+        denominatorB += (pwmB[index] - meanB) * (pwmB[index] - meanB);
+    }
+    
+    free(pwmA);
+    free(pwmB);
+    
+    if (denominatorA == 0.0 || denominatorB == 0.0) {
+        c = 0.0;
+    } else {
+        c = numerator / sqrt(denominatorA * denominatorB);
+    }
+    
+    return PyFloat_FromDouble(c);
+}
+
+static char the_func_doc[] =
+"pearsonCorrelation(a,b)\n\nReturns the pearsonCorrelation of a and b.";
+
+static PyMethodDef module_methods[] = {
+       {"pearsonCorrelation", the_func, METH_VARARGS, the_func_doc},
+       {NULL, NULL}
+};
+
+PyMODINIT_FUNC
+init_stat(void)
+{
+       Py_InitModule3("_stat", module_methods, module_doc);
+}