commit 5292197cbadff2b7103bc2e3dfd72d9c3d8b7208 Author: Daniel Knüttel Date: Wed Jul 24 19:19:41 2019 +0200 initial diff --git a/c/internal_deque.c b/c/internal_deque.c new file mode 100644 index 0000000..d0e4cdd --- /dev/null +++ b/c/internal_deque.c @@ -0,0 +1,70 @@ +#include "internal_deque.h" +#include + +i_deque_t * +i_deque_t_new + ( PyObject * content + , i_deque_t * next) +{ + i_deque_t * deque = malloc(sizeof(i_deque_t)); + if(!deque) + { + return NULL; + } + + deque->content = content; + deque->next = next; + + if(next) + { + deque->following = next->following + 1; + } + else + { + deque->following = 0; + } + return deque; +} + +void +i_deque_t_del + (i_deque_t * deque) +{ + if(deque) + { + free(deque); + } +} + +PyObject * +i_deque_t_to_List + (i_deque_t * deque) +{ + if(!deque) + { + return PyList_New(0); + } + PyObject * list = PyList_New(deque->following + 1); + if(!list) + { + return NULL; + } + i_deque_t * this_node; + PyObject * this_content; + + size_t i; + size_t content_length = deque->following + 1; + for(i = 0; i < content_length; i++) + { + this_node = deque; + deque = deque->next; + this_content = this_node->content; + i_deque_t_del(this_node); + Py_INCREF(this_content); + + PyList_SetItem(list, content_length - 1 - i, this_content); + } + return list; + +} + diff --git a/c/internal_deque.h b/c/internal_deque.h new file mode 100644 index 0000000..488ebed --- /dev/null +++ b/c/internal_deque.h @@ -0,0 +1,38 @@ +#ifndef interpol_internal_deque_h +#define interpol_internal_deque_h +#include + +/* + * + * This module contains a linked list implementation that is used internally + * as a deque. The internal deques are then converted to python lists. + * + * */ + +typedef struct i_deque_s +{ + struct i_deque_s * next; + PyObject * content; + size_t following; +} i_deque_t; + + +i_deque_t * +i_deque_t_new + ( PyObject * content + , i_deque_t * next); + +void +i_deque_t_del + (i_deque_t * deque); + +PyObject * +i_deque_t_to_List + (i_deque_t * deque); + + +#define i_deque_t_start() NULL +#define i_deque_t_insert(deque, content) i_deque_t_new(content, deque) + + +#endif diff --git a/c/splines/linear/linear_splines.c b/c/splines/linear/linear_splines.c new file mode 100644 index 0000000..7a26223 --- /dev/null +++ b/c/splines/linear/linear_splines.c @@ -0,0 +1,399 @@ +#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; +} diff --git a/c/splines/linear/linear_splines.h b/c/splines/linear/linear_splines.h new file mode 100644 index 0000000..d943783 --- /dev/null +++ b/c/splines/linear/linear_splines.h @@ -0,0 +1,35 @@ +#ifndef interpol_splines_linear_h +#define interpol_splines_linear_h + +#include +#include + +typedef struct _interpol_splines_linear_InterPolator_s +{ + PyObject_HEAD + long unsigned int grid_points; + PyArrayObject * x; + PyArrayObject * y; + + // This is either + // linear_spline_eval_double or linear_spline_eval_float and + // is used to evaluate the spline at a given x efficiently. + double (* eval)(double, struct _interpol_splines_linear_InterPolator_s *); +} interpol_splines_linear_InterPolator; + + +double linear_spline_eval_double + ( double x + , interpol_splines_linear_InterPolator * interpolator + ); +float linear_spline_eval_float + ( float x + , interpol_splines_linear_InterPolator * interpolator + ); + +#ifndef interpol_splines_linear_c + + +#endif + +#endif diff --git a/py/interpol/__init__.py b/py/interpol/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/py/interpol/__pycache__/__init__.cpython-35.pyc b/py/interpol/__pycache__/__init__.cpython-35.pyc new file mode 100644 index 0000000..b93baa2 Binary files /dev/null and b/py/interpol/__pycache__/__init__.cpython-35.pyc differ diff --git a/py/interpol/__pycache__/linear.cpython-35.pyc b/py/interpol/__pycache__/linear.cpython-35.pyc new file mode 100644 index 0000000..72941b2 Binary files /dev/null and b/py/interpol/__pycache__/linear.cpython-35.pyc differ diff --git a/py/interpol/linear.py b/py/interpol/linear.py new file mode 100644 index 0000000..a31b5d3 --- /dev/null +++ b/py/interpol/linear.py @@ -0,0 +1,39 @@ +import numpy as np +from collections import deque + +class HatLinearSplineInterpolator(object): + def __init__(self, x, y): + if(len(x) != len(y)): + raise ValueError("x and y have to have the same length") + + self._x = x + self._y = y + + def eval_at_point(self, x): + if(x < self._x[0]): + return 0 + if(x > self._x[-1]): + return 0 + + for i, x_test in enumerate(self._x[:-1]): + if(x_test <= x <= self._x[i + 1]): + break + + return self._y[i - 1] + (x - self._x[i - 1]) / (self._x[i] - self._x[i - 1]) * (self._y[i] - self._y[i - 1]) + + def equals(self, y): + result = deque() + + for i, y_test in enumerate(self._y): + if(not i): + continue + if((self._y[i - 1] <= y <= y_test) + or (self._y[i - 1] >= y >= y_test)): + if(self._y[i - 1] == y_test): + result.append(self._x[i - 1]) + continue + result.append(self._x[i - 1] + (y - self._y[i - 1]) / (self._y[i] - self._y[i - 1]) * (self._x[i] - self._x[i - 1])) + + return np.array(result) + + diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..dd325cd --- /dev/null +++ b/setup.py @@ -0,0 +1,17 @@ +from distutils.core import setup, Extension + + +interpol_spline_linear_do = Extension("interpol.spline.linear.do", + sources = ["c/splines/linear/linear_splines.c", "c/internal_deque.c"]) + +setup(name = "interpol", + version = "0.0.1", + description = "a numerical library for working with interpolation", + ext_modules = [interpol_spline_linear_do], + packages = [ + "interpol" + ], + package_dir = {"interpol": "py/interpol"}, + url="https://github.com/daknuett/python3-interpol", + author = "Daniel Knüttel", + author_email = "daniel.knuettel@daknuett.eu") diff --git a/test.py b/test.py new file mode 100644 index 0000000..ad95b62 --- /dev/null +++ b/test.py @@ -0,0 +1,5 @@ +import interpol.spline.linear.do as do +import numpy as np +x = np.arange(0, 4, 0.3) +y = np.cos(x) +ip = do.InterPolator(x, y)