[SciPy-user] Sample and hold

Robert Clewley rclewley at cam.cornell.edu
Sat Jun 3 19:21:08 CDT 2006

```Hi Martin,

On Sat, 3 Jun 2006, Martin Spacek wrote:
>
> I'd like the output to have regularly spaced data points. This means
> that the exact times of the original data points will be lost, but
> that's OK.

FWIW, you might be interested in a slightly more sophisticated solution:
abstracting the data and treating it like a true curve. You can have
arbitrary access to the curve for convenient resampling at your desired
interval length.

For this you could use the interp0d class that I adapted from scipy's
interp1d. You could also easily modify it to use the input points as the
left or right-hand side of the constant intervals (mine uses them as right
points, and averages the y values from the interval endpoints). Naturally,
this solution involves more overhead than the methods previously posted.

Below is the modified __call__ method of interp1d (from scipy 0.3.2) that
you need to create interp0d. Alternatively, the whole thing (including a
curve class over these interpolators) can be found in the PyDSTool
package. Note that the original data points are retained inside the class.
A demo of it working with independent time and data arrays in this package
also appears below.

-Rob

def __call__(self,x_new):
"""Find piecewise-constant interpolated y_new = <name>(x_new).

Inputs:
x_new -- New independent variables.

Outputs:
y_new -- Piecewise-constant interpolated values corresponding to x_new.
"""
# 1. Handle values in x_new that are outside of x.  Throw error,
#    or return a list of mask array indicating the outofbounds values.
#    The behavior is set by the bounds_error variable.
x_new_1d = atleast_1d(x_new)
out_of_bounds = self._check_bounds(x_new_1d)
# 2. Find where in the orignal data, the values to interpolate
#    would be inserted.
#    Note: If x_new[n] = x[m], then m is returned by searchsorted.
x_new_indices = searchsorted(self.x,x_new_1d)
# 3. Clip x_new_indices so that they are within the range of
#    self.x indices and at least 1.  Removes mis-interpolation
#    of x_new[n] = x[0]
x_new_indices = clip(x_new_indices,1,len(self.x)-1).astype(Int)
# 4. Calculate the region that each x_new value falls in.
lo = x_new_indices - 1; hi = x_new_indices

# !! take() should default to the last axis (IMHO) and remove
# !! the extra argument.
# 5. Calculate the actual value for each entry in x_new.
y_lo = take(self.y,lo,axis=self.interp_axis)
y_hi = take(self.y,hi,axis=self.interp_axis)
y_new = (y_lo+y_hi)/2.
# 6. Fill any values that were out of bounds with NaN
# !! Need to think about how to do this efficiently for
# !! mutli-dimensional Cases.
yshape = y_new.shape
y_new = y_new.flat
new_shape = list(yshape)
new_shape[self.interp_axis] = 1
sec_shape = [1]*len(new_shape)
sec_shape[self.interp_axis] = len(out_of_bounds)
out_of_bounds.shape = sec_shape
new_out = ones(new_shape)*out_of_bounds
y_new.shape = yshape
# Rotate the values of y_new back so that they correspond to the
# correct x_new values.
result = swapaxes(y_new,self.interp_axis,self.axis)
try:
len(x_new)
return result
except TypeError:
return result[0]
return result

Example of use in PyDSTool, using the InterpolateTable class:

>> timeData = linspace(0, 10, 30)
>> sindata = sin(timeData)
>> xData = {'sinx': sindata}
>> sin_pcwc_interp = InterpolateTable({'tdata': timeData,
'ics': xData,
'name': 'interp0d',
'method': 'constant',
'checklevel': 1,
'abseps': 1e-5
}).compute('sin_interp')
>> sin_pcwc_interp(2*pi)  # should be close to zero!
0.094554071417120258
>> t = linspace(0,10,300) # new time array, much finer resolution
>> y = sin_pcwc_interp(t)  # new data (as pointset, so use array(y) to get array)

```