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

[Joe Van Andel]

I'm very sorry, I previously sent the wrong patch.  Here's the correct
patch.

(The other was a previous version of arrayfnsmodule.c, that added an
interpf() command.  This version adds an optional typecode to allow
specifying single precision results from interp)

-- 
Joe VanAndel  	          
National Center for Atmospheric Research
http://www.atd.ucar.edu/~vanandel/
Internet: vanandel at ucar.edu
-------------- next part --------------
Index: arrayfnsmodule.c
===================================================================
RCS file: /cvsroot/numpy/Numerical/Src/arrayfnsmodule.c,v
retrieving revision 1.1
diff -u -r1.1 arrayfnsmodule.c
--- arrayfnsmodule.c	2000/01/20 18:18:28	1.1
+++ arrayfnsmodule.c	2000/03/28 17:55:04
@@ -575,6 +575,96 @@
    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 ;
+}
+/* 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 ;
+
+   PyObject *tpo = Py_None;  /* unused argument, we've already parsed it*/
+
+   Py_Try(PyArg_ParseTuple(args, "OOO|O", &oy, &ox, &oz, &tpo));
+   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 char arr_interp__doc__[] =
 ""
 ;
@@ -592,8 +682,24 @@
    PyArrayObject * ay, * ax, * az , * _interp;
    double * dy, * dx, * dz , * dres, * slopes;
    int leny, lenz, i, left ;
+   PyObject *tpo = Py_None;
+   char type_char = 'd';
+   char *type = &type_char;
 
-   Py_Try(PyArg_ParseTuple(args, "OOO", &oy, &ox, &oz));
+   Py_Try(PyArg_ParseTuple(args, "OOO|O", &oy, &ox, &oz,&tpo));
+   if (tpo != Py_None) {
+       type = PyString_AsString(tpo);
+       if (!type)
+	   return NULL;
+       if(!*type)
+	   type = &type_char;
+   }
+   if (*type == 'f' ) {
+       return arr_interpf(self, args);
+   } else if (*type != 'd') {
+       SETERR ("interp: unimplemented typecode.");
+       return NULL;
+   }
    GET_ARR(ay,oy,PyArray_DOUBLE,1);
    GET_ARR(ax,ox,PyArray_DOUBLE,1);
    if ( (leny = A_SIZE (ay)) != A_SIZE (ax)) {