[Scipy-svn] r4689 - in trunk/scipy/interpolate: . src
scipy-svn@scip...
scipy-svn@scip...
Fri Sep 5 07:50:13 CDT 2008
Author: oliphant
Date: 2008-09-05 07:50:12 -0500 (Fri, 05 Sep 2008)
New Revision: 4689
Added:
trunk/scipy/interpolate/src/_interpolate.cpp
trunk/scipy/interpolate/src/interpolate.h
Modified:
trunk/scipy/interpolate/setup.py
Log:
Add _interpolate functions.
Modified: trunk/scipy/interpolate/setup.py
===================================================================
--- trunk/scipy/interpolate/setup.py 2008-09-04 20:22:25 UTC (rev 4688)
+++ trunk/scipy/interpolate/setup.py 2008-09-05 12:50:12 UTC (rev 4689)
@@ -21,6 +21,11 @@
libraries=['fitpack'],
)
+ config.add_extension('_interpolate',
+ sources=['src/_interpolate.cpp'],
+ include_dirs = ['src'],
+ depends = ['src/interpolate.h'])
+
config.add_data_dir('tests')
return config
Added: trunk/scipy/interpolate/src/_interpolate.cpp
===================================================================
--- trunk/scipy/interpolate/src/_interpolate.cpp 2008-09-04 20:22:25 UTC (rev 4688)
+++ trunk/scipy/interpolate/src/_interpolate.cpp 2008-09-05 12:50:12 UTC (rev 4689)
@@ -0,0 +1,236 @@
+#include "Python.h"
+#include <stdlib.h>
+
+#include "interpolate.h"
+#include "numpy/arrayobject.h"
+
+using namespace std;
+
+extern "C" {
+
+static PyObject* linear_method(PyObject*self, PyObject* args, PyObject* kywds)
+{
+ static char *kwlist[] = {"x","y","new_x","new_y", NULL};
+ PyObject *py_x, *py_y, *py_new_x, *py_new_y;
+ py_x = py_y = py_new_x = py_new_y = NULL;
+ PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
+ arr_x = arr_y = arr_new_x = arr_new_y = NULL;
+
+ if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:linear_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
+ return NULL;
+ arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_x) {
+ PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_y) {
+ PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_new_x) {
+ PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
+ if (!arr_new_y) {
+ PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
+ goto fail;
+ }
+
+ linear((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
+ PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
+ (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
+
+ Py_DECREF(arr_x);
+ Py_DECREF(arr_y);
+ Py_DECREF(arr_new_x);
+ Py_DECREF(arr_new_y);
+
+ Py_RETURN_NONE;
+
+fail:
+ Py_XDECREF(arr_x);
+ Py_XDECREF(arr_y);
+ Py_XDECREF(arr_new_x);
+ Py_XDECREF(arr_new_y);
+ return NULL;
+}
+
+static PyObject* loginterp_method(PyObject*self, PyObject* args, PyObject* kywds)
+{
+ static char *kwlist[] = {"x","y","new_x","new_y", NULL};
+ PyObject *py_x, *py_y, *py_new_x, *py_new_y;
+ py_x = py_y = py_new_x = py_new_y = NULL;
+ PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
+ arr_x = arr_y = arr_new_x = arr_new_y = NULL;
+
+ if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
+ return NULL;
+ arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_x) {
+ PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_y) {
+ PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_new_x) {
+ PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
+ if (!arr_new_y) {
+ PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
+ goto fail;
+ }
+
+ loginterp((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
+ PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
+ (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
+
+ Py_DECREF(arr_x);
+ Py_DECREF(arr_y);
+ Py_DECREF(arr_new_x);
+ Py_DECREF(arr_new_y);
+
+ Py_RETURN_NONE;
+
+fail:
+ Py_XDECREF(arr_x);
+ Py_XDECREF(arr_y);
+ Py_XDECREF(arr_new_x);
+ Py_XDECREF(arr_new_y);
+ return NULL;
+}
+
+static PyObject* window_average_method(PyObject*self, PyObject* args, PyObject* kywds)
+{
+ static char *kwlist[] = {"x","y","new_x","new_y", NULL};
+ PyObject *py_x, *py_y, *py_new_x, *py_new_y;
+ py_x = py_y = py_new_x = py_new_y = NULL;
+ PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
+ arr_x = arr_y = arr_new_x = arr_new_y = NULL;
+ double width;
+
+ if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOOd:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y, &width))
+ return NULL;
+ arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_x) {
+ PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_y) {
+ PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_new_x) {
+ PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
+ if (!arr_new_y) {
+ PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
+ goto fail;
+ }
+
+ window_average((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
+ PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
+ (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0), width);
+
+ Py_DECREF(arr_x);
+ Py_DECREF(arr_y);
+ Py_DECREF(arr_new_x);
+ Py_DECREF(arr_new_y);
+
+ Py_RETURN_NONE;
+
+fail:
+ Py_XDECREF(arr_x);
+ Py_XDECREF(arr_y);
+ Py_XDECREF(arr_new_x);
+ Py_XDECREF(arr_new_y);
+ return NULL;
+}
+
+static PyObject* block_average_above_method(PyObject*self, PyObject* args, PyObject* kywds)
+{
+ static char *kwlist[] = {"x","y","new_x","new_y", NULL};
+ PyObject *py_x, *py_y, *py_new_x, *py_new_y;
+ py_x = py_y = py_new_x = py_new_y = NULL;
+ PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
+ arr_x = arr_y = arr_new_x = arr_new_y = NULL;
+
+ if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
+ return NULL;
+ arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_x) {
+ PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_y) {
+ PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_new_x) {
+ PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
+ if (!arr_new_y) {
+ PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
+ goto fail;
+ }
+
+ block_average_above((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
+ PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
+ (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
+
+ Py_DECREF(arr_x);
+ Py_DECREF(arr_y);
+ Py_DECREF(arr_new_x);
+ Py_DECREF(arr_new_y);
+
+ Py_RETURN_NONE;
+
+fail:
+ Py_XDECREF(arr_x);
+ Py_XDECREF(arr_y);
+ Py_XDECREF(arr_new_x);
+ Py_XDECREF(arr_new_y);
+ return NULL;
+}
+
+static PyMethodDef interpolate_methods[] = {
+ {"linear_dddd", (PyCFunction)linear_method, METH_VARARGS|METH_KEYWORDS,
+ ""},
+ {"loginterp_dddd", (PyCFunction)loginterp_method, METH_VARARGS|METH_KEYWORDS,
+ ""},
+ {"window_average_ddddd", (PyCFunction)window_average_method, METH_VARARGS|METH_KEYWORDS,
+ ""},
+ {"block_average_above_dddd", (PyCFunction)block_average_above_method, METH_VARARGS|METH_KEYWORDS,
+ ""},
+ {NULL, NULL, 0, NULL}
+};
+
+
+PyMODINIT_FUNC init_interpolate(void)
+{
+ PyObject* m;
+ m = Py_InitModule3("_interpolate", interpolate_methods,
+ "A few interpolation routines.\n"
+ );
+ if (m == NULL)
+ return;
+ import_array();
+}
+
+} // extern "C"
Added: trunk/scipy/interpolate/src/interpolate.h
===================================================================
--- trunk/scipy/interpolate/src/interpolate.h 2008-09-04 20:22:25 UTC (rev 4688)
+++ trunk/scipy/interpolate/src/interpolate.h 2008-09-05 12:50:12 UTC (rev 4689)
@@ -0,0 +1,205 @@
+#include <time.h>
+#include <math.h>
+#include <iostream>
+#include <algorithm>
+
+template <class T>
+void linear(T* x_vec, T* y_vec, int len,
+ T* new_x_vec, T* new_y_vec, int new_len)
+{
+ for (int i=0;i<new_len;i++)
+ {
+ T new_x = new_x_vec[i];
+ int index;
+ if (new_x <= x_vec[0])
+ index = 0;
+ else if (new_x >=x_vec[len-1])
+ index = len-2;
+ else
+ {
+ T* which = std::lower_bound(x_vec, x_vec+len, new_x);
+ index = which - x_vec-1;
+ }
+
+ if(new_x == x_vec[index])
+ {
+ // exact value
+ new_y_vec[i] = y_vec[index];
+ }
+ else
+ {
+ //interpolate
+ double x_lo = x_vec[index];
+ double x_hi = x_vec[index+1];
+ double y_lo = y_vec[index];
+ double y_hi = y_vec[index+1];
+ double slope = (y_hi-y_lo)/(x_hi-x_lo);
+ new_y_vec[i] = slope * (new_x-x_lo) + y_lo;
+ }
+ }
+}
+
+template <class T>
+void loginterp(T* x_vec, T* y_vec, int len,
+ T* new_x_vec, T* new_y_vec, int new_len)
+{
+ for (int i=0;i<new_len;i++)
+ {
+ T new_x = new_x_vec[i];
+ int index;
+ if (new_x <= x_vec[0])
+ index = 0;
+ else if (new_x >=x_vec[len-1])
+ index = len-2;
+ else
+ {
+ T* which = std::lower_bound(x_vec, x_vec+len, new_x);
+ index = which - x_vec-1;
+ }
+
+ if(new_x == x_vec[index])
+ {
+ // exact value
+ new_y_vec[i] = y_vec[index];
+ }
+ else
+ {
+ //interpolate
+ double x_lo = x_vec[index];
+ double x_hi = x_vec[index+1];
+ double y_lo = log10(y_vec[index]);
+ double y_hi = log10(y_vec[index+1]);
+ double slope = (y_hi-y_lo)/(x_hi-x_lo);
+ new_y_vec[i] = pow(10.0, (slope * (new_x-x_lo) + y_lo));
+ }
+ }
+}
+
+template <class T>
+int block_average_above(T* x_vec, T* y_vec, int len,
+ T* new_x_vec, T* new_y_vec, int new_len)
+{
+ int bad_index = -1;
+ int start_index = 0;
+ T last_y = 0.0;
+ T thickness = 0.0;
+
+ for(int i=0;i<new_len;i++)
+ {
+ T new_x = new_x_vec[i];
+ if ((new_x < x_vec[0]) || (new_x > x_vec[len-1]))
+ {
+ bad_index = i;
+ break;
+ }
+ else if (new_x == x_vec[0])
+ {
+ // for the first sample, just return the cooresponding y value
+ new_y_vec[i] = y_vec[0];
+ }
+ else
+ {
+ T* which = std::lower_bound(x_vec, x_vec+len, new_x);
+ int index = which - x_vec-1;
+
+ // calculate weighted average
+
+ // Start off with "residue" from last interval in case last x
+ // was between to samples.
+ T weighted_y_sum = last_y * thickness;
+ T thickness_sum = thickness;
+ for(int j=start_index; j<=index; j++)
+ {
+ if (x_vec[j+1] < new_x)
+ thickness = x_vec[j+1] - x_vec[j];
+ else
+ thickness = new_x -x_vec[j];
+ weighted_y_sum += y_vec[j] * thickness;
+ thickness_sum += thickness;
+ }
+ new_y_vec[i] = weighted_y_sum/thickness_sum;
+
+ // Store the thickness between the x value and the next sample
+ // to add to the next weighted average.
+ last_y = y_vec[index];
+ thickness = x_vec[index+1] - new_x;
+
+ // start next weighted average at next sample
+ start_index =index+1;
+ }
+ }
+ return bad_index;
+}
+
+template <class T>
+int window_average(T* x_vec, T* y_vec, int len,
+ T* new_x_vec, T* new_y_vec, int new_len,
+ T width)
+{
+ for(int i=0;i<new_len;i++)
+ {
+ T new_x = new_x_vec[i];
+ T bottom = new_x - width/2;
+ T top = new_x + width/2;
+
+ T* which = std::lower_bound(x_vec, x_vec+len, bottom);
+ int bottom_index = which - x_vec;
+ if (bottom_index < 0)
+ {
+ //bottom = x_vec[0];
+ bottom_index = 0;
+ }
+
+ which = std::lower_bound(x_vec, x_vec+len, top);
+ int top_index = which - x_vec;
+ if (top_index >= len)
+ {
+ //top = x_vec[len-1];
+ top_index = len-1;
+ }
+ //std::cout << std::endl;
+ //std::cout << bottom_index << " " << top_index << std::endl;
+ //std::cout << bottom << " " << top << std::endl;
+ // calculate weighted average
+ T thickness =0.0;
+ T thickness_sum =0.0;
+ T weighted_y_sum =0.0;
+ for(int j=bottom_index; j < top_index; j++)
+ {
+ thickness = x_vec[j+1] - bottom;
+ weighted_y_sum += y_vec[j] * thickness;
+ thickness_sum += thickness;
+ bottom = x_vec[j+1];
+ /*
+ std::cout << "iter: " << j - bottom_index << " " <<
+ "index: " << j << " " <<
+ "bottom: " << bottom << " " <<
+ "x+1: " << x_vec[j+1] << " " <<
+ "x: " << x_vec[j] << " " <<
+ "y: " << y_vec[j] << " " <<
+ "weighted_sum: " << weighted_y_sum <<
+ "thickness: " << thickness << " " <<
+ "thickness_sum: " << thickness_sum << std::endl;
+ */
+ //std::cout << x_vec[j] << " ";
+ //std::cout << thickness << " ";
+ }
+
+ // last element
+ thickness = top - bottom;
+ weighted_y_sum += y_vec[top_index] * thickness;
+ thickness_sum += thickness;
+ /*
+ std::cout << "iter: last" << " " <<
+ "index: " << top_index << " " <<
+ "x: " << x_vec[top_index] << " " <<
+ "y: " << y_vec[top_index] << " " <<
+ "weighted_sum: " << weighted_y_sum <<
+ "thickness: " << thickness << " " <<
+ "thickness_sum: " << thickness_sum << std::endl;
+ */
+ //std::cout << x_vec[top_index] << " " << thickness_sum << std::endl;
+ new_y_vec[i] = weighted_y_sum/thickness_sum;
+ }
+ return -1;
+}
More information about the Scipy-svn
mailing list