[Numpy-svn] r3645 - in trunk/numpy: . core/src lib lib/src
numpy-svn at scipy.org
numpy-svn at scipy.org
Mon Apr 2 16:22:11 EDT 2007
Author: oliphant
Date: 2007-04-02 15:21:57 -0500 (Mon, 02 Apr 2007)
New Revision: 3645
Modified:
trunk/numpy/add_newdocs.py
trunk/numpy/core/src/arraymethods.c
trunk/numpy/lib/function_base.py
trunk/numpy/lib/src/_compiled_base.c
Log:
Add interp to numpy.lib adapted from arrayfns. Add an unfinished arrayfns compatibility module to old_numeric.
Modified: trunk/numpy/add_newdocs.py
===================================================================
--- trunk/numpy/add_newdocs.py 2007-04-02 18:25:50 UTC (rev 3644)
+++ trunk/numpy/add_newdocs.py 2007-04-02 20:21:57 UTC (rev 3645)
@@ -1235,3 +1235,4 @@
Type can be either a new sub-type object or a data-descriptor object
"""))
+
Modified: trunk/numpy/core/src/arraymethods.c
===================================================================
--- trunk/numpy/core/src/arraymethods.c 2007-04-02 18:25:50 UTC (rev 3644)
+++ trunk/numpy/core/src/arraymethods.c 2007-04-02 20:21:57 UTC (rev 3645)
@@ -479,7 +479,9 @@
PyErr_SetString(PyExc_ValueError, "invalid integer");
return NULL;
}
- if (value >= PyArray_SIZE(self)) {
+ factor = PyArray_SIZE(self);
+ if (value < 0) value += factor;
+ if ((value >= factor) || (value < 0)) {
PyErr_SetString(PyExc_ValueError,
"index out of bounds");
return NULL;
Modified: trunk/numpy/lib/function_base.py
===================================================================
--- trunk/numpy/lib/function_base.py 2007-04-02 18:25:50 UTC (rev 3644)
+++ trunk/numpy/lib/function_base.py 2007-04-02 20:21:57 UTC (rev 3645)
@@ -8,7 +8,8 @@
'histogram', 'histogramdd', 'bincount', 'digitize', 'cov',
'corrcoef', 'msort', 'median', 'sinc', 'hamming', 'hanning',
'bartlett', 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc',
- 'add_docstring', 'meshgrid', 'delete', 'insert', 'append'
+ 'add_docstring', 'meshgrid', 'delete', 'insert', 'append',
+ 'interp'
]
import types
@@ -24,7 +25,7 @@
from numpy.lib.shape_base import atleast_1d, atleast_2d
from numpy.lib.twodim_base import diag
from _compiled_base import _insert, add_docstring
-from _compiled_base import digitize, bincount
+from _compiled_base import digitize, bincount, interp
from arraysetops import setdiff1d
#end Fernando's utilities
@@ -592,6 +593,25 @@
except RuntimeError:
pass
+try:
+ add_docstring(interp,
+r"""interp(x, xp, fp, left=None, right=None)
+
+Return the value of a piecewise-linear function at each value in x.
+
+The piecewise-linear function, f, is defined by the known data-points fp=f(xp).
+The xp points must be sorted in increasing order but this is not checked.
+
+For values of x < xp[0] return the value given by left. If left is None, then
+return fp[0].
+For values of x > xp[-1] return the value given by right. If right is None, then
+return fp[-1].
+"""
+ )
+except RuntimeError:
+ pass
+
+
def angle(z, deg=0):
"""Return the angle of the complex argument z.
"""
Modified: trunk/numpy/lib/src/_compiled_base.c
===================================================================
--- trunk/numpy/lib/src/_compiled_base.c 2007-04-02 18:25:50 UTC (rev 3644)
+++ trunk/numpy/lib/src/_compiled_base.c 2007-04-02 20:21:57 UTC (rev 3645)
@@ -337,7 +337,129 @@
return NULL;
}
+static npy_intp
+binary_search(double dval, double dlist [], npy_intp len)
+{
+ /* binary_search accepts three arguments: a numeric value and
+ * a numeric array and its length. It assumes that the array is sorted in
+ * increasing order. It returns the index of the array's
+ * largest element which is <= the value. It will return -1 if
+ * the value is less than the least element of the array. */
+ /* self is not used */
+ npy_intp bottom , top , middle, result;
+ if (dval < dlist [0])
+ result = -1 ;
+ else {
+ bottom = 0;
+ top = len - 1;
+ while (bottom < top) {
+ middle = (top + bottom) / 2 ;
+ if (dlist [middle] < dval)
+ bottom = middle + 1 ;
+ else if (dlist [middle] > dval)
+ top = middle - 1 ;
+ else
+ return middle ;
+ }
+ if (dlist [bottom] > dval)
+ result = bottom - 1 ;
+ else
+ result = bottom ;
+ }
+
+ return result ;
+}
+
+static PyObject *
+arr_interp(PyObject *self, PyObject *args, PyObject *kwdict)
+{
+
+ PyObject *fp, *xp, *x;
+ PyObject *left=NULL, *right=NULL;
+ PyArrayObject *afp=NULL, *axp=NULL, *ax=NULL, *af=NULL;
+ npy_intp i, lenx, lenxp, indx;
+ double lval, rval;
+ double *dy, *dx, *dz, *dres, *slopes;
+
+ static char *kwlist[] = {"x", "xp", "fp", "left", "right", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwdict, "OOO|OO", kwlist,
+ &x, &xp, &fp, &left, &right))
+ return NULL;
+
+
+ afp = (NPY_AO*)PyArray_ContiguousFromAny(fp, NPY_DOUBLE, 1, 1);
+ if (afp == NULL) return NULL;
+ axp = (NPY_AO*)PyArray_ContiguousFromAny(xp, NPY_DOUBLE, 1, 1);
+ if (axp == NULL) goto fail;
+ ax = (NPY_AO*)PyArray_ContiguousFromAny(x, NPY_DOUBLE, 1, 0);
+ if (ax == NULL) goto fail;
+
+ lenxp = axp->dimensions[0];
+ if (afp->dimensions[0] != lenxp) {
+ PyErr_SetString(PyExc_ValueError, "interp: fp and xp are not the same length.");
+ goto fail;
+ }
+
+ af = (NPY_AO*)PyArray_SimpleNew(ax->nd, ax->dimensions, NPY_DOUBLE);
+ if (af == NULL) goto fail;
+
+ lenx = PyArray_SIZE(ax);
+
+ dy = (double *)PyArray_DATA(afp);
+ dx = (double *)PyArray_DATA(axp);
+ dz = (double *)PyArray_DATA(ax);
+ dres = (double *)PyArray_DATA(af);
+
+ /* Get left and right fill values. */
+ if ((left == NULL) || (left == Py_None)) {
+ lval = dy[0];
+ }
+ else {
+ lval = PyFloat_AsDouble(left);
+ if ((lval==-1) && PyErr_Occurred())
+ goto fail;
+ }
+ if ((right == NULL) || (right == Py_None)) {
+ rval = dy[lenxp-1];
+ }
+ else {
+ rval = PyFloat_AsDouble(right);
+ if ((rval==-1) && PyErr_Occurred())
+ goto fail;
+ }
+
+ slopes = (double *) PyDataMem_NEW((lenxp-1)*sizeof(double));
+ for (i=0; i < lenxp-1; i++) {
+ slopes[i] = (dy[i+1] - dy[i])/(dx[i+1]-dx[i]);
+ }
+ for (i=0; i<lenx; i++) {
+ indx = binary_search(dz[i], dx, lenxp);
+ if (indx < 0)
+ dres[i] = lval;
+ else if (indx >= lenxp - 1)
+ dres[i] = rval;
+ else
+ dres[i] = slopes[indx]*(dz[i]-dx[indx]) + dy[indx];
+ }
+
+ PyDataMem_FREE(slopes);
+ Py_DECREF(afp);
+ Py_DECREF(axp);
+ Py_DECREF(ax);
+ return (PyObject *)af;
+
+ fail:
+ Py_XDECREF(afp);
+ Py_XDECREF(axp);
+ Py_XDECREF(ax);
+ Py_XDECREF(af);
+ return NULL;
+}
+
+
+
static PyTypeObject *PyMemberDescr_TypePtr=NULL;
static PyTypeObject *PyGetSetDescr_TypePtr=NULL;
static PyTypeObject *PyMethodDescr_TypePtr=NULL;
@@ -408,8 +530,10 @@
arr_insert__doc__},
{"bincount", (PyCFunction)arr_bincount,
METH_VARARGS | METH_KEYWORDS, NULL},
- {"digitize", (PyCFunction)arr_digitize, METH_VARARGS | METH_KEYWORDS,
+ {"digitize", (PyCFunction)arr_digitize, METH_VARARGS | METH_KEYWORDS,
NULL},
+ {"interp", (PyCFunction)arr_interp, METH_VARARGS | METH_KEYWORDS,
+ NULL},
{"add_docstring", (PyCFunction)arr_add_docstring, METH_VARARGS,
NULL},
{NULL, NULL} /* sentinel */
@@ -435,7 +559,7 @@
return;
}
-/* Initialization function for the module (*must* be called initArray) */
+/* Initialization function for the module (*must* be called init<name>) */
PyMODINIT_FUNC init_compiled_base(void) {
PyObject *m, *d, *s;
More information about the Numpy-svn
mailing list