[SciPy-user] step function

Robert Kern robert.kern at gmail.com
Mon Mar 13 16:41:58 CST 2006


Jonathan Taylor wrote:
> believe it or not, i just discovered "searchsorted"... my earlier hack 
> was brutal, using "bisect"....

Ouch! Yeah, that may be why it was slow.  :-)

> this is my current version and it seems to work fine.

It looks good. I'll just make a few comments inline.

> class StepFunction:
>     '''A basic step function: values at the ends are handled in the 
> simplest way possible: everything to the left of x[0] is set to ival; 
> everything to the right of x[-1] is set to y[-1].
>    
>     >>>
>     >>> from numpy import *
>     >>>
>     >>> x = arange(20)
>     >>> y = arange(20)
>     >>>
>     >>> f = StepFunction(x, y)
>     >>>
>     >>> print f(3.2)
>     3
>     >>> print f([[3.2,4.5],[24,-3.1]])
>     [[ 3  4]
>      [19  0]]
>     >>>
> 
>     '''
> 
>     def __init__(self, x, y, ival=0., sorted=False):
>            
>         _x = N.array(x)
>         _y = N.array(y)

N.asarray() will be a tad more efficient.

>         if _x.shape != _y.shape:
>             raise ValueError, 'in StepFunction: x and y do not have the 
> same shape'
>         if len(_x.shape) != 1:
>             raise ValueError, 'in StepFunction: x and y must be 
> 1-dimensional'
> 
>         self.x = N.array([-N.inf] + list(x))
>         self.y = N.array([ival] + list(y))

You can do this more efficiently by using numpy.concatenate() (or numpy.hstack()
which happens to be equivalent in this particular case).

self.x = N.hstack(([-N.inf], _x))

>         if not sorted:
>             asort = N.argsort(self.x)
>             self.x = N.take(self.x, asort)
>             self.y = N.take(self.y, asort)
>         self.n = self.x.shape[0]
> 
>     def __call__(self, time):
> 
>         tind = scipy.searchsorted(self.x, time) - 1
>         _shape = tind.shape
>         if tind.shape is ():
>             return self.y[tind]
>         else:
>             tmp = N.take(self.y, tind)
>             tmp.shape = _shape
>             return tmp

self.y[tind] will give the correct answer regardless. In numpy, you can index
into an array using another array. The shape of the returned array will have the
shape of the index array.

-- 
Robert Kern
robert.kern at gmail.com

"I have come to believe that the whole world is an enigma, a harmless enigma
 that is made terrible by our own mad attempt to interpret it as though it had
 an underlying truth."
  -- Umberto Eco



More information about the SciPy-user mailing list