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__},