As I've mentioned earlier, I'm using interp() with very large datasets,
and I can't afford to use double precision arrays. Here's a patch that
lets interp() accept an optional typecode
argument. Passing 'f' calls the new single precision version, and
returns a single precision
array. Passing no argument or 'd' uses the previous, double precision
version. (No other array types are supported - an error is returned.)
I hope this can be added to the CVS version of Numeric.
Thanks!
--
Joe VanAndel
National Center for Atmospheric Research
http://www.atd.ucar.edu/~vanandel/
Internet: vanandel@ucar.edu
*** arrayfnsmodule.c 1999/04/14 22:58:32 1.1
--- arrayfnsmodule.c 1999/08/12 23:27:28
***************
*** 6,11 ****
--- 6,13 ----
#include
#include
+ #define MAX_INTERP_DIMS 6
+
static PyObject *ErrorObject;
/* Define 2 macros for error handling:
***************
*** 34,39 ****
--- 36,43 ----
#define A_DIM(a,i) (((PyArrayObject *)a)->dimensions[i])
#define GET_ARR(ap,op,type,dim) \
Py_Try(ap=(PyArrayObject *)PyArray_ContiguousFromObject(op,type,dim,dim))
+ #define GET_ARR2(ap,op,type,min,max) \
+ Py_Try(ap=(PyArrayObject *)PyArray_ContiguousFromObject(op,type,min,max))
#define ERRSS(s) ((PyObject *)(PyErr_SetString(ErrorObject,s),0))
#define SETERR(s) if(!PyErr_Occurred()) ERRSS(errstr ? errstr : s)
#define DECREF_AND_ZERO(p) do{Py_XDECREF(p);p=0;}while(0)
***************
*** 571,576 ****
--- 575,613 ----
return result ;
}
+ static int
+ binary_searchf(float dval, float dlist [], int 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 */
+ int 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 char arr_interp__doc__[] =
""
;
***************
*** 597,609 ****
Py_DECREF(ay);
Py_DECREF(ax);
return NULL ;}
! GET_ARR(az,oz,PyArray_DOUBLE,1);
lenz = A_SIZE (az);
dy = (double *) A_DATA (ay);
dx = (double *) A_DATA (ax);
dz = (double *) A_DATA (az);
! Py_Try (_interp = (PyArrayObject *) PyArray_FromDims (1, &lenz,
! PyArray_DOUBLE));
dres = (double *) A_DATA (_interp) ;
slopes = (double *) malloc ( (leny - 1) * sizeof (double)) ;
for (i = 0 ; i < leny - 1; i++) {
--- 634,647 ----
Py_DECREF(ay);
Py_DECREF(ax);
return NULL ;}
! GET_ARR2(az,oz,PyArray_DOUBLE,1,MAX_INTERP_DIMS);
lenz = A_SIZE (az);
dy = (double *) A_DATA (ay);
dx = (double *) A_DATA (ax);
dz = (double *) A_DATA (az);
! /* create output array with same size as 'Z' input array */
! Py_Try (_interp = (PyArrayObject *) PyArray_FromDims
! (A_NDIM(az), az->dimensions, PyArray_DOUBLE));
dres = (double *) A_DATA (_interp) ;
slopes = (double *) malloc ( (leny - 1) * sizeof (double)) ;
for (i = 0 ; i < leny - 1; i++) {
***************
*** 627,632 ****
--- 665,724 ----
return PyArray_Return (_interp);
}
+ /* return float, rather than double */
+
+ static PyObject *
+ arr_interpf(PyObject *self, PyObject *args)
+ {
+ /* interp (y, x, z) treats (x, y) as a piecewise linear function
+ * whose value is y [0] for x < x [0] and y [len (y) -1] for x >
+ * x [len (y) -1]. An array of floats the same length as z is
+ * returned, whose values are ordinates for the corresponding z
+ * abscissae interpolated into the piecewise linear function. */
+ /* self is not used */
+ PyObject * oy, * ox, * oz ;
+ PyArrayObject * ay, * ax, * az , * _interp;
+ float * dy, * dx, * dz , * dres, * slopes;
+ int leny, lenz, i, left ;
+
+ Py_Try(PyArg_ParseTuple(args, "OOO", &oy, &ox, &oz));
+ GET_ARR(ay,oy,PyArray_FLOAT,1);
+ GET_ARR(ax,ox,PyArray_FLOAT,1);
+ if ( (leny = A_SIZE (ay)) != A_SIZE (ax)) {
+ SETERR ("interp: x and y are not the same length.");
+ Py_DECREF(ay);
+ Py_DECREF(ax);
+ return NULL ;}
+ GET_ARR2(az,oz,PyArray_FLOAT,1,MAX_INTERP_DIMS);
+ lenz = A_SIZE (az);
+ dy = (float *) A_DATA (ay);
+ dx = (float *) A_DATA (ax);
+ dz = (float *) A_DATA (az);
+ /* create output array with same size as 'Z' input array */
+ Py_Try (_interp = (PyArrayObject *) PyArray_FromDims
+ (A_NDIM(az), az->dimensions, PyArray_FLOAT));
+ dres = (float *) A_DATA (_interp) ;
+ slopes = (float *) malloc ( (leny - 1) * sizeof (float)) ;
+ for (i = 0 ; i < leny - 1; i++) {
+ slopes [i] = (dy [i + 1] - dy [i]) / (dx [i + 1] - dx [i]) ;
+ }
+ for ( i = 0 ; i < lenz ; i ++ )
+ {
+ left = binary_searchf (dz [i], dx, leny) ;
+ if (left < 0)
+ dres [i] = dy [0] ;
+ else if (left >= leny - 1)
+ dres [i] = dy [leny - 1] ;
+ else
+ dres [i] = slopes [left] * (dz [i] - dx [left]) + dy [left];
+ }
+
+ free (slopes);
+ Py_DECREF(ay);
+ Py_DECREF(ax);
+ Py_DECREF(az);
+ return PyArray_Return (_interp);
+ }
static int incr_slot_ (float x, double *bins, int lbins)
{
int i ;
***************
*** 1295,1300 ****
--- 1387,1393 ----
{"array_set", arr_array_set, 1, arr_array_set__doc__},
{"index_sort", arr_index_sort, 1, arr_index_sort__doc__},
{"interp", arr_interp, 1, arr_interp__doc__},
+ {"interpf", arr_interpf, 1, arr_interp__doc__},
{"digitize", arr_digitize, 1, arr_digitize__doc__},
{"zmin_zmax", arr_zmin_zmax, 1, arr_zmin_zmax__doc__},
{"reverse", arr_reverse, 1, arr_reverse__doc__},
