daknuett/interpol
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

some bugfixes

Browse Source
This commit is contained in:
Daniel Knüttel 2019-12-05 08:38:14 +01:00
parent 91cd4c24c3
commit edf6b032c7
3 changed files with 64 additions and 196 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

214
c/splines/linear/linear_splines.c
View File

@ -4,80 +4,48 @@
double linear_spline_eval_double( double linear_spline_eval_double(
double x double x
, interpol_splines_linear_InterPolator * interpolator) , interpol_splines_linear_Interpolator * self)
{ {
unsigned long int i = 0; unsigned long int i = 0;
double * X = PyArray_DATA(interpolator->x); double * X = PyArray_DATA(self->x);
double * Y = PyArray_DATA(interpolator->y); double * Y = PyArray_DATA(self->y);
// Outside the known data. // Outside the known data.
if(x < X[0]) if(x < X[0])
{ {
return 0; return 0;
} }
if(x > X[self->grid_points - 1])
// 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; return 0;
} }
// Somehow this has to be handled explicitly. i = 0;
if(x == X[0]) while(i < self->grid_points - 1)
{ {
i = 1; if((X[i] <= x) && (x <= X[i + 1]))
{
break;
} }
else
{
while(x < X[i] && i < interpolator->grid_points)
{
i++; i++;
} }
}
return Y[i - 1] + (x - X[i - 1]) / (X[i] - X[i - 1]) * (Y[i] - Y[i - 1]); return Y[i - 1] + (x - X[i - 1]) / (X[i] - X[i - 1]) * (Y[i] - Y[i - 1]);
} }
static PyObject * static PyObject *
interpol_splines_linear_InterPolator_new interpol_splines_linear_Interpolator_new
( PyTypeObject * type ( PyTypeObject * type
, PyObject * args , PyObject * args
, PyObject * kwds) , 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; return (PyObject *) self;
} }
static int static int
interpol_splines_linear_InterPolator_init interpol_splines_linear_Interpolator_init
( interpol_splines_linear_InterPolator * self ( interpol_splines_linear_Interpolator * self
, PyObject * args) , PyObject * args)
{ {
if(!PyArg_ParseTuple(args, "O!O!" 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"); PyErr_SetString(PyExc_TypeError, "x and y have to have the same data type");
return -1; return -1;
} }
if(PyArray_TYPE(self->x) == NPY_FLOAT)
if(PyArray_TYPE(self->x) == NPY_DOUBLE)
{ {
self->eval = linear_spline_eval_double; self->eval = linear_spline_eval_double;
} }
else if(PyArray_TYPE(self->x) == NPY_DOUBLE)
{
self->eval = linear_spline_eval_float;
}
else 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; return -1;
} }
@ -142,8 +108,8 @@ interpol_splines_linear_InterPolator_init
} }
static PyObject * static PyObject *
interpol_splines_linear_InterPolator_eval_pyfloat interpol_splines_linear_Interpolator_eval_pyfloat
( interpol_splines_linear_InterPolator * self ( interpol_splines_linear_Interpolator * self
, PyObject * args , PyObject * args
) )
{ {
@ -156,110 +122,14 @@ interpol_splines_linear_InterPolator_eval_pyfloat
result = (double) (self->eval)(x, self); result = (double) (self->eval)(x, self);
return Py_BuildValue("d", result); 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"); , "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 int
interpol_splines_linear_InterPolator_equals_pyfloat_search_float interpol_splines_linear_Interpolator_equals_pyfloat_search_double
( interpol_splines_linear_InterPolator * self ( 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 , i_deque_t ** deque
, double y) , double y)
{ {
@ -270,7 +140,8 @@ interpol_splines_linear_InterPolator_equals_pyfloat_search_double
long unsigned int i; long unsigned int i;
for(i = 1; i < self->grid_points; 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 // just avoid zero division
if(y_data[i - 1] == y_data[i]) if(y_data[i - 1] == y_data[i])
@ -296,8 +167,8 @@ interpol_splines_linear_InterPolator_equals_pyfloat_search_double
return 0; return 0;
} }
static PyObject * static PyObject *
interpol_splines_linear_InterPolator_equals_pyfloat interpol_splines_linear_Interpolator_equals_pyfloat
( interpol_splines_linear_InterPolator * self ( interpol_splines_linear_Interpolator * self
, PyObject * args) , PyObject * args)
{ {
double y; double y;
@ -307,31 +178,24 @@ interpol_splines_linear_InterPolator_equals_pyfloat
} }
i_deque_t * deque = i_deque_t_start(); i_deque_t * deque = i_deque_t_start();
if(PyArray_TYPE(self->y) == NPY_FLOAT) interpol_splines_linear_Interpolator_equals_pyfloat_search_double(self, &deque, y);
{
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); 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"); , "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} {"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_O, interpol_splines_linear_InterPolator_equals_pyfloat_doc} , {"equals_float", (PyCFunction) interpol_splines_linear_Interpolator_equals_pyfloat, METH_VARARGS, interpol_splines_linear_Interpolator_equals_pyfloat_doc}
, {NULL, NULL, 0, NULL} , {NULL, NULL, 0, NULL}
}; };
static PyTypeObject interpol_splines_linear_InterPolatorType = static PyTypeObject interpol_splines_linear_InterpolatorType =
{ {
PyVarObject_HEAD_INIT(NULL, 0) PyVarObject_HEAD_INIT(NULL, 0)
"interpol.splines.linear.do.Interpolator", "interpol.splines.linear.do.Interpolator",
sizeof(interpol_splines_linear_InterPolator), sizeof(interpol_splines_linear_Interpolator),
0, /* tp_itemsize */ 0, /* tp_itemsize */
0, /* tp_dealloc */ 0, /* tp_dealloc */
0, /* tp_print */ 0, /* tp_print */
@ -357,7 +221,7 @@ static PyTypeObject interpol_splines_linear_InterPolatorType =
0, 0,
0, 0,
0, 0,
interpol_splines_linear_InterPolator_Methods, // methods interpol_splines_linear_Interpolator_Methods, // methods
0, // members 0, // members
0, 0,
0, 0,
@ -365,9 +229,9 @@ static PyTypeObject interpol_splines_linear_InterPolatorType =
0, 0,
0, 0,
0, 0,
(initproc) interpol_splines_linear_InterPolator_init, (initproc) interpol_splines_linear_Interpolator_init,
0, 0,
interpol_splines_linear_InterPolator_new, interpol_splines_linear_Interpolator_new,
}; };
static PyModuleDef interpol_splines_linear_do_module = static PyModuleDef interpol_splines_linear_do_module =
@ -382,7 +246,7 @@ static PyModuleDef interpol_splines_linear_do_module =
PyMODINIT_FUNC PyInit_do(void) PyMODINIT_FUNC PyInit_do(void)
{ {
PyObject * module; PyObject * module;
if(PyType_Ready(&interpol_splines_linear_InterPolatorType) < 0) if(PyType_Ready(&interpol_splines_linear_InterpolatorType) < 0)
{ {
return NULL; return NULL;
} }
@ -391,8 +255,8 @@ PyMODINIT_FUNC PyInit_do(void)
{ {
return NULL; return NULL;
} }
Py_INCREF(&interpol_splines_linear_InterPolatorType); Py_INCREF(&interpol_splines_linear_InterpolatorType);
PyModule_AddObject(module, "InterPolator", (PyObject *) &interpol_splines_linear_InterPolatorType); PyModule_AddObject(module, "Interpolator", (PyObject *) &interpol_splines_linear_InterpolatorType);
import_array(); import_array();
return module; return module;

10
c/splines/linear/linear_splines.h
View File

@ -4,7 +4,7 @@
#include <Python.h> #include <Python.h>
#include <numpy/arrayobject.h> #include <numpy/arrayobject.h>
typedef struct _interpol_splines_linear_InterPolator_s typedef struct _interpol_splines_linear_Interpolator_s
{ {
PyObject_HEAD PyObject_HEAD
long unsigned int grid_points; long unsigned int grid_points;
@ -14,17 +14,17 @@ typedef struct _interpol_splines_linear_InterPolator_s
// This is either // This is either
// linear_spline_eval_double or linear_spline_eval_float and // linear_spline_eval_double or linear_spline_eval_float and
// is used to evaluate the spline at a given x efficiently. // is used to evaluate the spline at a given x efficiently.
double (* eval)(double, struct _interpol_splines_linear_InterPolator_s *); double (* eval)(double, struct _interpol_splines_linear_Interpolator_s *);
} interpol_splines_linear_InterPolator; } interpol_splines_linear_Interpolator;
double linear_spline_eval_double double linear_spline_eval_double
( double x ( double x
, interpol_splines_linear_InterPolator * interpolator , interpol_splines_linear_Interpolator * interpolator
); );
float linear_spline_eval_float float linear_spline_eval_float
( float x ( float x
, interpol_splines_linear_InterPolator * interpolator , interpol_splines_linear_Interpolator * interpolator
); );
#ifndef interpol_splines_linear_c #ifndef interpol_splines_linear_c

22
test.py
View File

@ -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 numpy as np
import matplotlib.pyplot as plt 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) results_half = i.equals_float((0.5,))
y = np.cos(x) y_half = [0.5 for x in results_half]
ip = do.InterPolator(x, y) plt.plot(X_bar, interpolated, "r-")
plt.plot(X_bar, Y_bar, "g-")
x_bar = np.arange(0, 4, 0.01) plt.plot(results_half, y_half, "bo")
plt.plot(x_bar, np.cos(x_bar))
plt.plot(x_bar, [ip.eval_float((i,)) for i in x_bar])
plt.show() plt.show()