/*
 *  untitled.c
 *  SOM
 *
 *  Created by Ali Mortazavi on 12/28/09.
 *  Copyright 2009 __MyCompanyName__. All rights reserved.
 *
 */

#include <Python.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static char module_doc[] = 
"This module only implements only some of the Self-Organizing Map code in C for now.";

static PyObject*
cProp(PyObject *self, PyObject *args)
{
    PyObject *SOM, *SOMrow, *SOMunit, *pyVec, *pyRows, *pyCols, *pyDim;
    long nrows, ncols, dim, winRow, winCol;
    double *inputVec, currentVal;
    int row, col, index; 
    float bestScore, currentDist;   
    
    if (!PyArg_UnpackTuple(args, "func", 5, 5, &SOM, &pyVec, &pyRows, &pyCols, &pyDim)) {
        return NULL;
    }
    
    nrows = PyInt_AsLong(pyRows);
    ncols = PyInt_AsLong(pyCols);
    dim = PyInt_AsLong(pyDim);
    bestScore = 16777216.0;
    winRow = 0;
    winCol = 0;
    
    inputVec = malloc(dim * sizeof(double));
    for (index = 0; index < dim; index++) {
        inputVec[index] = PyFloat_AsDouble(PyList_GetItem(pyVec, index));
    }
    
    for (row = 0; row < nrows; row++) {
        SOMrow = PyList_GetItem(SOM, row);
        for (col = 0; col < ncols; col++) {
            currentDist = 0.0;
            SOMunit = PyList_GetItem(SOMrow, col);
            for (index = 0; index < dim; index++) {
                /*currentVal = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(PyList_GetItem(SOM, row), col),index));*/
                currentVal = PyFloat_AsDouble(PyList_GetItem(SOMunit, index));
                currentDist = currentDist + (currentVal - inputVec[index]) * (currentVal - inputVec[index]);
            }
            /* fprintf(stdout,"row %d col %d current %f best %f\n", row, col, currentDist, bestScore); */
            if (currentDist < bestScore) {
                winRow = row;
                winCol = col;
                bestScore = currentDist;
            }
        }
    }
    
    free(inputVec);
    
    return Py_BuildValue("(ii)", winRow, winCol);
}

static PyObject*
cEuclidean(PyObject *self, PyObject *args)
{
    PyObject *vecA, *vecB;
    double distance, a, b;
    int dim, index; 
    
    if (!PyArg_UnpackTuple(args, "func", 2, 2, &vecA, &vecB)) {
        return NULL;
    }
    
    dim = PyList_Size(vecA);
    distance = 0.0;
    
    for (index = 0; index < dim; index++) {
        a = PyFloat_AsDouble(PyList_GetItem(vecA, index));
        b = PyFloat_AsDouble(PyList_GetItem(vecB, index));
        distance = distance + (a - b) * (a - b);
    }
    
    distance = sqrt(distance);
    
    return PyFloat_FromDouble(distance);
}

static PyObject*
cUpdateUnit(PyObject *self, PyObject *args)
{
    PyObject *oldVec, *inVec, *pyRate;
    PyObject *outVec;
    double old_i, in_i, rate;
    int dim, index; 
    
    if (!PyArg_UnpackTuple(args, "func", 3, 3, &oldVec, &inVec, &pyRate)) {
        return NULL;
    }
    
    dim = PyList_Size(oldVec);
    rate = PyFloat_AsDouble(pyRate);
    outVec = PyList_New(dim);
    
    for (index = 0; index < dim; index++) {
        old_i = PyFloat_AsDouble(PyList_GetItem(oldVec, index));
        in_i = PyFloat_AsDouble(PyList_GetItem(inVec, index));
        PyList_SET_ITEM(outVec, index, PyFloat_FromDouble(old_i + rate * (in_i - old_i)));
    }
    
    return outVec;
}


static char cProp_doc[] =
"cPropagate(SOM,inputVec, nrows, ncols, dim)\n\nReturns the winning unit on the SOM of nrows x ncols for inputVec of 1 x dim.";

static char cEuclidean_doc[] =
"cEuclidean(vecA, vecB)\n\nReturns the eucledian distance between vecA and vecB.";

static char cUpdateUnit_doc[] =
"cUpdateUnit(origVec, inputVec, rate)\n\nReturns origVec + rate * (inputVec - origVec).";


static PyMethodDef module_methods[] = {
	{"cPropagate", cProp, METH_VARARGS, cProp_doc},
        {"cEuclidean", cEuclidean, METH_VARARGS, cEuclidean_doc},
        {"cUpdateUnit", cUpdateUnit, METH_VARARGS, cUpdateUnit_doc},
	{NULL, NULL}
};

PyMODINIT_FUNC
init_csom(void)
{
	Py_InitModule3("_csom", module_methods, module_doc);
}