[Numpy-discussion] single precision patch for arrayfnsmodule.c :interp()

Joe Van Andel vanandel at atd.ucar.edu
Tue Mar 28 12:00:12 CST 2000


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 at ucar.edu
-------------- next part --------------
*** arrayfnsmodule.c	1999/04/14 22:58:32	1.1
--- arrayfnsmodule.c	1999/08/12 23:27:28
***************
*** 6,11 ****
--- 6,13 ----
  #include <stdio.h>
  #include <stdlib.h>
  
+ #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__},


More information about the Numpy-discussion mailing list