400 lines
8.9 KiB
C
400 lines
8.9 KiB
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;
|
||
|
|
}
|