interpol/c/splines/linear/linear_splines.c

#define interpol_splines_linear_c
#include "linear_splines.h"
#include "../../internal_deque.h"

double linear_spline_eval_double(
		double x
		, interpol_splines_linear_Interpolator * self)
{
	unsigned long int i = 0;
	double * X = PyArray_DATA(self->x);
	double * Y = PyArray_DATA(self->y);

	// Outside the known data.
	if(x < X[0])
	{
		return 0;
	}
    if(x > X[self->grid_points - 1])
    {
        return 0;
    }

    i = 0;
    while(i < self->grid_points - 1)
    {
        if((X[i] <= x) && (x <= X[i + 1]))
        {
            break;
        }
        i++;
    }

	return Y[i - 1] + (x - X[i - 1]) / (X[i] - X[i - 1]) * (Y[i] - Y[i - 1]);
}

static PyObject * 
interpol_splines_linear_Interpolator_new
	( PyTypeObject * type
	, PyObject * args
	, PyObject * kwds)
{
	interpol_splines_linear_Interpolator * self = (interpol_splines_linear_Interpolator *) type->tp_alloc(type, 0);
	return (PyObject *) self;
}

static int 
interpol_splines_linear_Interpolator_init
	( interpol_splines_linear_Interpolator * self
	, PyObject * args)
{
	if(!PyArg_ParseTuple(args, "O!O!"
				, &PyArray_Type, &(self->x)
				, &PyArray_Type, &(self->y)))
	{
		return -1;
	}

	// Check for proper data types & dimensions first.
	if(self->x->nd != 1)
	{
		PyErr_SetString(PyExc_ValueError, "x has to have 1 dimension");
		return -1;
	}
	if(self->y->nd != 1)
	{
		PyErr_SetString(PyExc_ValueError, "y has to have 1 dimension");
		return -1;
	}

	if(PyArray_TYPE(self->x) != PyArray_TYPE(self->y))
	{
		PyErr_SetString(PyExc_TypeError, "x and y have to have the same data type");
		return -1;
	}


	if(PyArray_TYPE(self->x) == NPY_DOUBLE)
	{
		self->eval = linear_spline_eval_double;
	}
	else
	{
		PyErr_SetString(PyExc_TypeError, "x and y have to be double");
		return -1;
	}

	if(self->x->dimensions[0] != self->y->dimensions[0])
	{
		PyErr_SetString(PyExc_ValueError, "x and y must have the same length");
		return -1;
	}
	if((PyArray_FLAGS(self->x) & NPY_CARRAY_RO) != NPY_CARRAY_RO)
	{
		PyErr_SetString(PyExc_ValueError, "x must be a readable C Array");
		return -1;
	}
	if((PyArray_FLAGS(self->y) & NPY_CARRAY_RO) != NPY_CARRAY_RO)
	{
		PyErr_SetString(PyExc_ValueError, "y must be a readable C Array");
		return -1;
	}

	self->grid_points = self->x->dimensions[0];
	Py_INCREF(self->x);
	Py_INCREF(self->y);

	return 0;
}

static PyObject *
interpol_splines_linear_Interpolator_eval_pyfloat
	( interpol_splines_linear_Interpolator * self
	, PyObject * args
	)
{
	double x;
	double result;
	if(!PyArg_ParseTuple(args, "d", &x))
	{
		return NULL;
	}
	result = (double) (self->eval)(x, self);
	return Py_BuildValue("d", result);
}
PyDoc_STRVAR(interpol_splines_linear_Interpolator_eval_pyfloat_doc
	, "evaluate the spline at the given point");



int
interpol_splines_linear_Interpolator_equals_pyfloat_search_double
	( interpol_splines_linear_Interpolator * self
	, i_deque_t ** deque
	, double y)
{
	double * x_data = (double *) PyArray_DATA(self->x),
	      * y_data = (double *) PyArray_DATA(self->y);
	PyObject * this_result;
	
	long unsigned int i;
	for(i = 1; i < self->grid_points; i++)
	{
		if(((y_data[i - 1]) <= y && (y <= y_data[i]))
                || ((y_data[i - 1] >= y) && (y >= y_data[i])))
		{
			// just avoid zero division
			if(y_data[i - 1] == y_data[i])
			{
				this_result = PyFloat_FromDouble(x_data[i - 1]);
				if(!this_result)
				{
					return -1;
				}
				*deque = i_deque_t_insert(*deque, this_result);
				continue;
			}
			this_result = PyFloat_FromDouble(
					(y - y_data[i - 1]) / (y_data[i] - y_data[i - 1]) * (x_data[i] - x_data[i - 1]) + x_data[i - 1]);
			if(!this_result)
			{
				return -1;
			}
			*deque = i_deque_t_insert(*deque, this_result);

		}
	}
	return 0;	
}
static PyObject *
interpol_splines_linear_Interpolator_equals_pyfloat
	( interpol_splines_linear_Interpolator * self
	, PyObject * args)
{
	double y;
	if(!PyArg_ParseTuple(args, "d", &y))
	{
		return NULL;
	}
	i_deque_t * deque = i_deque_t_start();

    interpol_splines_linear_Interpolator_equals_pyfloat_search_double(self, &deque, y);
	return i_deque_t_to_List(deque);
}
PyDoc_STRVAR(interpol_splines_linear_Interpolator_equals_pyfloat_doc
	, "return a deque to all the points where s(x) = y");

static PyMethodDef interpol_splines_linear_Interpolator_Methods[] = 
{
	{"eval_float", (PyCFunction) interpol_splines_linear_Interpolator_eval_pyfloat, METH_VARARGS, interpol_splines_linear_Interpolator_eval_pyfloat_doc}
	, {"equals_float", (PyCFunction) interpol_splines_linear_Interpolator_equals_pyfloat, METH_VARARGS, interpol_splines_linear_Interpolator_equals_pyfloat_doc}
	, {NULL, NULL, 0, NULL}
};

static PyTypeObject interpol_splines_linear_InterpolatorType = 
{
	PyVarObject_HEAD_INIT(NULL, 0)
	"interpol.splines.linear.do.Interpolator",
	sizeof(interpol_splines_linear_Interpolator),
	0,                 /* tp_itemsize */
	0,                 /* tp_dealloc */
	0,                 /* tp_print */
	0,                 /* tp_getattr */
	0,                 /* tp_setattr */
	0,                 /* tp_reserved */
	0,                 /* tp_repr */
	0,                 /* tp_as_number */
	0,                 /* tp_as_sequence */
	0,                 /* tp_as_mapping */
	0,                 /* tp_hash  */
	0,                 /* tp_call */
	0,                 /* tp_str */
	0,                 /* tp_getattro */
	0,                 /* tp_setattro */
	0,                 /* tp_as_buffer */
	Py_TPFLAGS_DEFAULT ,/* tp_flags */
	"The backend for linear splines.",

	0,
	0,
	0,
	0,
	0,
	0,
	interpol_splines_linear_Interpolator_Methods, // methods
	0, // members
	0,
	0,
	0,
	0,
	0,
	0,
	(initproc) interpol_splines_linear_Interpolator_init,
	0,
	interpol_splines_linear_Interpolator_new,
};

static PyModuleDef interpol_splines_linear_do_module = 
{
	PyModuleDef_HEAD_INIT,
	"interpol.splines.linear.do",
	"module containing the backend to linear splines",
	-1,
	NULL,NULL,NULL,NULL,NULL
};

PyMODINIT_FUNC PyInit_do(void)
{
	PyObject * module;
	if(PyType_Ready(&interpol_splines_linear_InterpolatorType) < 0)
	{
		return NULL;
	}
	module = PyModule_Create(&interpol_splines_linear_do_module);
	if(!module)
	{
		return NULL;
	}
	Py_INCREF(&interpol_splines_linear_InterpolatorType);
	PyModule_AddObject(module, "Interpolator", (PyObject *) &interpol_splines_linear_InterpolatorType);

	import_array();
	return module;
}