/*##########################################################################*/
/*#                                                                         #*/
/*# 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);
}