diff --git a/c/splines/linear/linear_splines.c b/c/splines/linear/linear_splines.c index 7a26223..2294cbe 100644 --- a/c/splines/linear/linear_splines.c +++ b/c/splines/linear/linear_splines.c @@ -4,80 +4,48 @@ double linear_spline_eval_double( double x - , interpol_splines_linear_InterPolator * interpolator) + , interpol_splines_linear_Interpolator * self) { unsigned long int i = 0; - double * X = PyArray_DATA(interpolator->x); - double * Y = PyArray_DATA(interpolator->y); + 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; + } - // Somehow this has to be handled explicitly. - if(x == X[0]) - { - i = 1; - } - else - { - while(x < X[i] && i < interpolator->grid_points) - { - i++; - } - } + 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]); - } -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 +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); + 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 +interpol_splines_linear_Interpolator_init + ( interpol_splines_linear_Interpolator * self , PyObject * args) { if(!PyArg_ParseTuple(args, "O!O!" @@ -104,17 +72,15 @@ interpol_splines_linear_InterPolator_init PyErr_SetString(PyExc_TypeError, "x and y have to have the same data type"); return -1; } - if(PyArray_TYPE(self->x) == NPY_FLOAT) + + + if(PyArray_TYPE(self->x) == NPY_DOUBLE) { 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"); + PyErr_SetString(PyExc_TypeError, "x and y have to be double"); return -1; } @@ -142,8 +108,8 @@ interpol_splines_linear_InterPolator_init } static PyObject * -interpol_splines_linear_InterPolator_eval_pyfloat - ( interpol_splines_linear_InterPolator * self +interpol_splines_linear_Interpolator_eval_pyfloat + ( interpol_splines_linear_Interpolator * self , PyObject * args ) { @@ -156,110 +122,14 @@ interpol_splines_linear_InterPolator_eval_pyfloat result = (double) (self->eval)(x, self); return Py_BuildValue("d", result); } -PyDoc_STRVAR(interpol_splines_linear_InterPolator_eval_pyfloat_doc +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 +interpol_splines_linear_Interpolator_equals_pyfloat_search_double + ( interpol_splines_linear_Interpolator * self , i_deque_t ** deque , double y) { @@ -270,7 +140,8 @@ interpol_splines_linear_InterPolator_equals_pyfloat_search_double long unsigned int i; for(i = 1; i < self->grid_points; i++) { - if(y_data[i - 1] <= y && y_data[i] >= y) + 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]) @@ -296,8 +167,8 @@ interpol_splines_linear_InterPolator_equals_pyfloat_search_double return 0; } static PyObject * -interpol_splines_linear_InterPolator_equals_pyfloat - ( interpol_splines_linear_InterPolator * self +interpol_splines_linear_Interpolator_equals_pyfloat + ( interpol_splines_linear_Interpolator * self , PyObject * args) { double y; @@ -307,31 +178,24 @@ interpol_splines_linear_InterPolator_equals_pyfloat } 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); - } + 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 +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[] = +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} + {"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 = +static PyTypeObject interpol_splines_linear_InterpolatorType = { PyVarObject_HEAD_INIT(NULL, 0) "interpol.splines.linear.do.Interpolator", - sizeof(interpol_splines_linear_InterPolator), + sizeof(interpol_splines_linear_Interpolator), 0, /* tp_itemsize */ 0, /* tp_dealloc */ 0, /* tp_print */ @@ -357,7 +221,7 @@ static PyTypeObject interpol_splines_linear_InterPolatorType = 0, 0, 0, - interpol_splines_linear_InterPolator_Methods, // methods + interpol_splines_linear_Interpolator_Methods, // methods 0, // members 0, 0, @@ -365,9 +229,9 @@ static PyTypeObject interpol_splines_linear_InterPolatorType = 0, 0, 0, - (initproc) interpol_splines_linear_InterPolator_init, + (initproc) interpol_splines_linear_Interpolator_init, 0, - interpol_splines_linear_InterPolator_new, + interpol_splines_linear_Interpolator_new, }; static PyModuleDef interpol_splines_linear_do_module = @@ -382,7 +246,7 @@ static PyModuleDef interpol_splines_linear_do_module = PyMODINIT_FUNC PyInit_do(void) { PyObject * module; - if(PyType_Ready(&interpol_splines_linear_InterPolatorType) < 0) + if(PyType_Ready(&interpol_splines_linear_InterpolatorType) < 0) { return NULL; } @@ -391,8 +255,8 @@ PyMODINIT_FUNC PyInit_do(void) { return NULL; } - Py_INCREF(&interpol_splines_linear_InterPolatorType); - PyModule_AddObject(module, "InterPolator", (PyObject *) &interpol_splines_linear_InterPolatorType); + 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 index d943783..ad10f34 100644 --- a/c/splines/linear/linear_splines.h +++ b/c/splines/linear/linear_splines.h @@ -4,7 +4,7 @@ #include #include -typedef struct _interpol_splines_linear_InterPolator_s +typedef struct _interpol_splines_linear_Interpolator_s { PyObject_HEAD long unsigned int grid_points; @@ -14,17 +14,17 @@ typedef struct _interpol_splines_linear_InterPolator_s // 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 (* eval)(double, struct _interpol_splines_linear_Interpolator_s *); +} interpol_splines_linear_Interpolator; double linear_spline_eval_double ( double x - , interpol_splines_linear_InterPolator * interpolator + , interpol_splines_linear_Interpolator * interpolator ); float linear_spline_eval_float ( float x - , interpol_splines_linear_InterPolator * interpolator + , interpol_splines_linear_Interpolator * interpolator ); #ifndef interpol_splines_linear_c diff --git a/test.py b/test.py index 9caaedf..6a89886 100644 --- a/test.py +++ b/test.py @@ -1,13 +1,17 @@ -import interpol.spline.linear.do as do +# coding: utf-8 +from interpol.spline.linear.do import Interpolator import numpy as np - import matplotlib.pyplot as plt +X = np.arange(0, 30, 0.3) +Y = np.cos(X) +i = Interpolator(X, Y) +X_bar = np.arange(0, 30, 0.01) +Y_bar = np.cos(X_bar) +interpolated = np.array([i.eval_float((f,)) for f in X_bar]) -x = np.arange(0, 4, 0.3) -y = np.cos(x) -ip = do.InterPolator(x, y) - -x_bar = np.arange(0, 4, 0.01) -plt.plot(x_bar, np.cos(x_bar)) -plt.plot(x_bar, [ip.eval_float((i,)) for i in x_bar]) +results_half = i.equals_float((0.5,)) +y_half = [0.5 for x in results_half] +plt.plot(X_bar, interpolated, "r-") +plt.plot(X_bar, Y_bar, "g-") +plt.plot(results_half, y_half, "bo") plt.show()