cistematic/cisstat/statextension.c

   1 /*##########################################################################*/
   2 /*#                                                                         #*/
   3 /*# C O P Y R I G H T   N O T I C E                                         #*/
   4 /*#  Copyright (c) 2003-10 by:                                              #*/
   5 /*#    * California Institute of Technology                                 #*/
   6 /*#                                                                         #*/
   7 /*#    All Rights Reserved.                                                 #*/
   8 /*#                                                                         #*/
   9 /*# Permission is hereby granted, free of charge, to any person             #*/
  10 /*# obtaining a copy of this software and associated documentation files    #*/
  11 /*# (the "Software"), to deal in the Software without restriction,          #*/
  12 /*# including without limitation the rights to use, copy, modify, merge,    #*/
  13 /*# publish, distribute, sublicense, and/or sell copies of the Software,    #*/
  14 /*# and to permit persons to whom the Software is furnished to do so,       #*/
  15 /*# subject to the following conditions:                                    #*/
  16 /*#                                                                         #*/
  17 /*# The above copyright notice and this permission notice shall be          #*/
  18 /*# included in all copies or substantial portions of the Software.         #*/
  19 /*#                                                                         #*/
  20 /*# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #*/
  21 /*# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #*/
  22 /*# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #*/
  23 /*# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #*/
  24 /*# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #*/
  25 /*# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #*/
  26 /*# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #*/
  27 /*# SOFTWARE.                                                               #*/
  28 /*###########################################################################*/
  29 /*# C extension for commonly used stat methods.                             */
  30
  31 #include <Python.h>
  32 #include <stdio.h>
  33 #include <stdlib.h>
  34 #include <math.h>
  35
  36 static char module_doc[] =
  37 "This module only implements the pearsonCorrelation in C for now.";
  38
  39 static PyObject*
  40 the_func(PyObject *self, PyObject *args)
  41 {
  42     PyObject *a, *b;
  43     double c, numerator;
  44     double *pwmA, meanA, denominatorA;
  45     double *pwmB, meanB, denominatorB;
  46     long listLen, index;
  47
  48     if (!PyArg_UnpackTuple(args, "func", 2, 2, &a, &b)) {
  49         return NULL;
  50     }
  51
  52     listLen = PyList_Size(a);
  53     pwmA = malloc(listLen * sizeof(double));
  54     for (index = 0; index < listLen; index++) {
  55         pwmA[index] = PyFloat_AsDouble(PyList_GetItem(a, index));
  56     }
  57
  58     pwmB = malloc(PyList_Size(b) * sizeof(double));
  59     for (index = 0; index < listLen; index++) {
  60         pwmB[index] = PyFloat_AsDouble(PyList_GetItem(b, index));
  61     }
  62
  63     if (listLen > PyList_Size(b)) {
  64         listLen = PyList_Size(b);
  65     }
  66
  67     meanA = 0.0;
  68     meanB = 0.0;
  69
  70     for (index = 0; index < listLen; index++) {
  71         meanA += pwmA[index];
  72         meanB += pwmB[index];
  73     }
  74
  75     meanA /= listLen;
  76     meanB /= listLen;
  77
  78     denominatorA = 0.0;
  79     denominatorB = 0.0;
  80     numerator = 0.0;
  81
  82     for (index = 0; index < listLen; index++) {
  83         numerator += (pwmA[index] - meanA) * (pwmB[index] - meanB);
  84         denominatorA += (pwmA[index] - meanA) * (pwmA[index] - meanA);
  85         denominatorB += (pwmB[index] - meanB) * (pwmB[index] - meanB);
  86     }
  87
  88     free(pwmA);
  89     free(pwmB);
  90
  91     if (denominatorA == 0.0 || denominatorB == 0.0) {
  92         c = 0.0;
  93     } else {
  94         c = numerator / sqrt(denominatorA * denominatorB);
  95     }
  96
  97     return PyFloat_FromDouble(c);
  98 }
  99
 100 static char the_func_doc[] =
 101 "pearsonCorrelation(a,b)\n\nReturns the pearsonCorrelation of a and b.";
 102
 103 static PyMethodDef module_methods[] = {
 104         {"pearsonCorrelation", the_func, METH_VARARGS, the_func_doc},
 105         {NULL, NULL}
 106 };
 107
 108 PyMODINIT_FUNC
 109 init_stat(void)
 110 {
 111         Py_InitModule3("_stat", module_methods, module_doc);
 112 }