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 * interpolator)
{
	unsigned long int i = 0;
	double * X = PyArray_DATA(interpolator->x);
	double * Y = PyArray_DATA(interpolator->y);

	// Outside the known data.
	if(x < X[0])
	{
		return 0;
	}

	// Somehow this has to be handled explicitly.
	if(x == X[0])
	{
		i = 1;
	}
	else
	{
		while(x < X[i] && i < interpolator->grid_points)
		{
			i++;
		}
	}

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

}
float linear_spline_eval_float(
		float x
		, interpol_splines_linear_InterPolator * interpolator)
{
	unsigned long int i = 0;
	float * X = PyArray_DATA(interpolator->x);
	float * Y = PyArray_DATA(interpolator->y);

	// Outside the known data.
	if(x < X[0])
	{
		return 0;
	}

	// Somehow this has to be handled explicitly.
	if(x == X[0])
	{
		i = 1;
	}
	else
	{
		while(x < X[i] && i < interpolator->grid_points)
		{
			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_FLOAT)
	{
		self->eval = linear_spline_eval_double;
	}
	else if(PyArray_TYPE(self->x) == NPY_DOUBLE)
	{
		self->eval = linear_spline_eval_float;
	}
	else
	{
		PyErr_SetString(PyExc_TypeError, "x and y have to be either float or 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");

/*
static PyObject *
interpol_splines_linear_InterPolator_eval_ndarray
	( interpol_splines_linear_InterPolator * self
	, PyObject * args
	)
{
	PyArray * x;
	PyArray * y;
	if(!PyArg_ParseTuple(args, "O!O!", &PyArray_Type, &x, &PyArray_Type, &y))
	{
		return NULL;
	}
	if(PyArray_FLAGS(x) & NPY_CARRAY_RO != NPY_CARRAY_RO)
	{
		PyErr_SetString(PyExc_ValueError, "x must be a readable c array");
		return NULL;
	}
	if(PyArray_FLAGS(y) & NPY_CARRAY_RW != NPY_CARRAY_RW)
	{
		PyErr_SetString(PyExc_ValueError, "y must be a readable and writable c array");
		return NULL;
	}

	if(PyArray_TYPE(x) != NPY_FLOAT 
		&& PyArray_TYPE(x) != NPY_DOUBLE)
	{
		PyErr_SetString(PyExc_TypeError, "x must be either float or double");
		return NULL;
	}

	if(PyArray_TYPE(x) != PyArray_TYPE(y))
	{
		PyErr_SetString(PyExc_TypeError, "x and y must have same datatype");
		return NULL;
	}

	if(PyArray_Size(x) != PyArray_Size(y))
	{
		PyErr_SetString(PyExc_ValueError, "x and y must have same length");
	}

	npy_intp size = PyArray_NBYTES(x);
	npy_intp element_size = PyArray_ITEMSIZE(x);

	char * x_data = PyArray_BYTES(x);
	char * y_data = PyArray_BYTES(y);
	
	// FIXME:
	// I have to pay attention to datatypes here.
	while(x_data < size)
	{
		*y_data = self->eval(*x, self);
	}

}
*/


int
interpol_splines_linear_InterPolator_equals_pyfloat_search_float
	( interpol_splines_linear_InterPolator * self
	, i_deque_t ** deque
	, double y)
{
	float * x_data = (float *) PyArray_DATA(self->x),
	      * y_data = (float *) 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_data[i] >= y)
		{
			// 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;	
}

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_data[i] >= y)
		{
			// 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();

	if(PyArray_TYPE(self->y) == NPY_FLOAT)
	{
		interpol_splines_linear_InterPolator_equals_pyfloat_search_float(self, &deque, y);
	}
	else
	{
		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_O, interpol_splines_linear_InterPolator_eval_pyfloat_doc}
	, {"equals_float", (PyCFunction) interpol_splines_linear_InterPolator_equals_pyfloat, METH_O, 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;
}