Reworked linear_prediction and revcorr to handle the new slicestim

pyret/filtertools.py

@@ -198,7 +198,7 @@ def lowranksta(f_orig, k=10):
     Notes
     -----
 
-    This method requires that the linear filter be 3D. To decompose a 
+    This method requires that the linear filter be 3D. To decompose a
     linear filter into a temporal and 1-dimensional spatial filter, simply
     promote the filter to 3D before calling this method.
 
@@ -445,8 +445,8 @@ def get_contours(spatial_filter, threshold=10.0): # pragma: no cover
     Gets iso-value contours of a 2D spatial filter.
 
     This returns a list of arrays of shape (n, 2). Each array in the list
-    gives the indices into the spatial filter of one contour, and each 
-    column of the contour gives the indices along the two dimesions of 
+    gives the indices into the spatial filter of one contour, and each
+    column of the contour gives the indices along the two dimesions of
     the filter.
 
     Usage
@@ -473,9 +473,9 @@ def get_contours(spatial_filter, threshold=10.0): # pragma: no cover
 
 def get_regionprops(spatial_filter, threshold=10.0): # pragma: no cover
     """
-    Gets region properties of a 2D spatial filter. 
-    
-    This returns various attributes of the non-zero area of the given 
+    Gets region properties of a 2D spatial filter.
+
+    This returns various attributes of the non-zero area of the given
     spatial filter, such as its area, centroid, eccentricity, etc.
 
     Usage
@@ -543,7 +543,7 @@ def get_ellipse(spatial_filter, sigma=2.):
                                           p0=pinit)
 
     # return ellipse parameters, scaled by the appropriate scale factor
-    
+
     return _popt_to_ellipse(*popt, sigma)
 
 
@@ -592,46 +592,38 @@ def linear_prediction(filt, stim):
 
     Parameters
     ----------
-
     filt : array_like
         The linear filter whose response is to be computed. The array should
-        have shape ``(t, ...)``, where ``t`` is the number of time points in the 
+        have shape ``(t, ...)``, where ``t`` is the number of time points in the
         filter and the ellipsis indicates any remaining spatial dimenions.
         The number of dimensions and the sizes of the spatial dimensions
         must match that of ``stim``.
 
     stim : array_like
         The stimulus to which the predicted response is computed. The array
-        should have shape (T,...), where ``T`` is the number of time points 
+        should have shape (T,...), where ``T`` is the number of time points
         in the stimulus and the ellipsis indicates any remaining spatial
         dimensions. The number of dimensions and the sizes of the spatial
         dimenions must match that of ``filt``.
 
     Returns
     -------
-
     pred : array_like
-        The predicted linear response. The shape is (T,) where T is the 
+        The predicted linear response. The shape is (T,) where T is the
         number of time points in the input stimulus array.
 
     Raises
     ------
-
     ValueError : If the number of dimensions of ``stim`` and ``filt`` do not
         match, or if the spatial dimensions differ.
-
     """
-
     if (filt.ndim != stim.ndim) or (filt.shape[1:] != stim.shape[1:]):
         raise ValueError("The filter and stimulus must have the same " +
-                "number of dimensions and match in size along spatial dimensions")
+                         "number of dimensions and match in size along spatial dimensions")
 
     slices = slicestim(stim, filt.shape[0])
-    dim_start = ord('i')
-    indices = ''.join(map(chr, range(dim_start, dim_start + slices.ndim)))
-    subscripts = '{0},{1}{2}->{3}'.format(indices, indices[0], 
-            indices[2:], indices[1])
-    return np.einsum(subscripts, slices, filt)
+    T = slices.shape[0]
+    return np.einsum('tx,x->t', slices.reshape(T, -1), filt.ravel())
 
 
 def revcorr(response, stimulus, filter_length):
@@ -645,12 +637,11 @@ def revcorr(response, stimulus, filter_length):
 
     Parameters
     ----------
-
     response : array_like
         A continuous output response correlated with the stimulus. Must
         be one-dimensional.
 
-    stimulus : array_like 
+    stimulus : array_like
         A input stimulus correlated with the ``response``. Must be of shape
         (t, ...), where t is the time and ... indicates any spatial dimensions.
 
@@ -660,39 +651,33 @@ def revcorr(response, stimulus, filter_length):
 
     Returns
     -------
-
     filt : array_like
         An array of shape ``(filter_length, ...)`` containing the best-fitting
         linear filter which predicts the response from the stimulus. The ellipses
         indicates spatial dimensions of the filter.
 
     Raises
     ------
-
     ValueError : If the ``stimulus`` and ``response`` arrays are of different
     shapes.
 
     Notes
     -----
-
     The ``response`` and ``stimulus`` arrays must share the same sampling
     rate. As the stimulus often has a lower sampling rate, one can use
     ``stimulustools.upsamplestim`` to upsample it.
-
     """
-
     if response.ndim > 1:
         raise ValueError("The `response` must be 1-dimensional")
-    if response.size != (stimulus.shape[0] - filter_length):
-        raise ValueError(("`stimulus` must have {:#d} time points " + 
-                "(`response.size` + `filter_length`").format(response.size + filter_length))
-    
+    if response.size != (stimulus.shape[0] - filter_length + 1):
+        msg = "`stimulus` must have {:#d} time points (`response.size` + `filter_length`)"
+        raise ValueError(msg.format(response.size + filter_length + 1))
+
     slices = slicestim(stimulus, filter_length)
-    dim_start = ord('i')
-    indices = ''.join(map(chr, range(dim_start, dim_start + slices.ndim)))
-    subscripts = '{0},{1}->{2}{3}'.format(indices, indices[1], indices[0], indices[2:])
-    recovered = np.einsum(subscripts, slices, response)
-    return recovered
+    T = slices.shape[0]
+    recovered = np.einsum('tx,t->x', slices.reshape(T, -1), response)
+    return recovered.reshape(slices.shape[1:])
+
 
 def _gaussian_function(data, x0, y0, a, b, c):
     """
@@ -794,7 +779,7 @@ def _initial_gaussian_params(xm, ym, z, width=5):
     ----------
     xm : array_like
         The x-points for the filter.
-        
+
     ym : array_like
         The y-points for the filter.
 

pyret/stimulustools.py

@@ -120,8 +120,8 @@ def slicestim(stimulus, history):
 
     Examples
     --------
-    >>> x=np.arange(10).reshape((2,5))
-    >>> rolling_window(x, 3)
+    >>> x=np.arange(10).reshape((5, 2))
+    >>> slicestim(x, 3)
     array([[[0, 1, 2], [1, 2, 3], [2, 3, 4]],
            [[5, 6, 7], [6, 7, 8], [7, 8, 9]]])
 
@@ -131,18 +131,19 @@ def slicestim(stimulus, history):
     array([[ 1.,  2.,  3.],
            [ 6.,  7.,  8.]])
     """
-
     if not (1 <= history <= stimulus.shape[0]):
-        raise ValueError(
-                '`history` must be between 1 and {0:#d}'.format(
-                stimulus.shape[0]))
-    history = int(history)
+        msg = '`history` must be between 1 and {0:#d}'.format(stimulus.shape[0])
+        raise ValueError(msg)
+    elif not isinstance(history, int):
+        raise ValueError("`history` must be an integer")
 
     # Use strides to create view onto array
-    shape = (history, stimulus.shape[0] - history) + stimulus.shape[1:]
-    strides = (stimulus.strides[0],) + stimulus.strides
-    return np.lib.stride_tricks.as_strided(stimulus, 
-            shape=shape, strides=strides)
+    shape = (history, stimulus.shape[0] - history + 1) + stimulus.shape[1:]
+    stride = (stimulus.strides[0],) + stimulus.strides
+
+    # return the newly strided array
+    arr = np.lib.stride_tricks.as_strided(stimulus, shape=shape, strides=stride)
+    return np.rollaxis(arr, 1)
 
 
 def cov(stimulus, history, nsamples=None, verbose=False):
@@ -160,60 +161,13 @@ def cov(stimulus, history, nsamples=None, verbose=False):
     history : int
         Integer number of time points to keep in each slice.
 
-    nsamples : int_like, optional
-
-    verbose : boolean, optional
-        If True, print out progress of the computation. (defaults to False)
-
     Returns
     ------
     stim_cov : array_like
-        (n*n*t by n*n*t) Covariance matrix
-
+        Covariance matrix
     """
-    if not (1 <= history < stimulus.shape[0]):
-        raise ValueError('`history` must be in [1, {0:#d})'.format(
-                stimulus.shape[0]))
-
-    # Collapse spatial dimensions
-    cstim = stimulus.reshape(stimulus.shape[0], -1)
-
-    # store mean + covariance matrix
-    mean = np.zeros(cstim.shape[-1] * history)
-    stim_cov = np.zeros((cstim.shape[-1] * history,) * 2)
-
-    # Update covariance matrix from random points in the stimulus
-    indices = np.arange(history, cstim.shape[0])
-    np.random.shuffle(indices)
-
-    # BLAS rank-1 covariance update function
-    blas_ger_fnc = get_blas_funcs(('ger',), (stim_cov,))[0]
-
-    if nsamples is None:
-        nsamples = indices.size
-
-    for j in range(nsamples):
-        idx = indices[j]
-        if verbose and np.mod(j, 100) == 0:
-            print('[%i of %i]' % (j, numpts))
-
-        stimslice = np.reshape(cstim[idx - history : idx, :], 
-                (-1, 1), order='F')
-
-        mean += np.squeeze(stimslice)
-
-        # Update covariance matrix
-        blas_ger_fnc(1, stimslice, stimslice, 
-                a=stim_cov.T, overwrite_a=True)
-
-    # Normalize and compute the mean outer product
-    mean /= nsamples
-    mean_op = mean.reshape(-1, 1).dot(mean.reshape(1, -1))
-
-    # Return normalized covariance matrix
-    stim_cov = (stim_cov / (nsamples - 1)) - mean_op
-
-    return stim_cov
+    stim = slicestim(stimulus, history)
+    return np.cov(stim.reshape(stim.shape[0], -1).T)
 
 
 def rolling_window(array, window, time_axis=0):

tests/test_filtertools.py

@@ -47,7 +47,7 @@ def test_empty_sta():
     spikes = np.array(())
     stimulus = np.random.randn(100,)
     filter_length = 5
-    
+
     sta, _ = flt.sta(time, stimulus, spikes, filter_length)
     assert np.all(np.isnan(sta))
 
@@ -154,13 +154,13 @@ def test_rfsize():
     filter_length = 100
     nx, ny = 10, 10
     true_filter = utils.create_spatiotemporal_filter(nx, ny, filter_length)[1]
-    
+
     xsize, ysize = flt.rfsize(true_filter, 1., 1.)
     assert np.allclose(xsize, 3., 0.1) # 1 SD is about 3 units
     assert np.allclose(ysize, 3., 0.1)
 
 def test_linear_prediction_1d():
-    """Test method for computing linear prediction from a 
+    """Test method for computing linear prediction from a
     filter to a one-dimensional stimulus.
     """
     np.random.seed(0)
@@ -169,10 +169,11 @@ def test_linear_prediction_1d():
     pred = flt.linear_prediction(filt, stim)
 
     sl = slicestim(stim, filt.shape[0])
-    assert np.allclose(filt.reshape(1, -1).dot(sl), pred)
+    assert np.allclose(sl.dot(filt), pred)
+
 
 def test_linear_prediction_nd():
-    """Test method for computing linear prediction from a 
+    """Test method for computing linear prediction from a
     filter to a multi-dimensional stimulus.
     """
     np.random.seed(0)
@@ -182,10 +183,9 @@ def test_linear_prediction_nd():
         pred = flt.linear_prediction(filt, stim)
 
         sl = slicestim(stim, filt.shape[0])
-        tmp = np.zeros(sl.shape[1])
-        filt_reshape = filt.reshape(1, -1)
+        tmp = np.zeros(sl.shape[0])
         for i in range(tmp.size):
-            tmp[i] = filt_reshape.dot(sl[:, i, :].reshape(-1, 1))
+            tmp[i] = np.inner(filt.ravel(), sl[i].ravel())
 
         assert np.allclose(tmp, pred)
 

tests/test_stimulustools.py

@@ -49,7 +49,7 @@ def test_slicestim_1d():
     sliced_stim = stimulustools.slicestim(stim, history)
 
     for i in range(stim_size - history):
-        assert np.all(sliced_stim[:, i] == stim[i : i + history]), 'slicing failed'
+        assert np.all(sliced_stim[i] == stim[i:i + history]), 'slicing failed'
 
 def test_slicestim_3d():
     """Test slicing a 3D stimulus into overlapping segments."""
@@ -60,19 +60,13 @@ def test_slicestim_3d():
 
     sliced_stim = stimulustools.slicestim(stim, history)
     assert sliced_stim.ndim == stim.ndim + 1
-    assert sliced_stim.shape[1] == stim.shape[0] - history
+    assert sliced_stim.shape[0] == stim.shape[0] - history + 1
 
     for i in range(stim_size[0] - history):
-        assert np.all(sliced_stim[:, i, ...] == stim[i : i + history, ...]), 'slicing failed'
+        assert np.all(sliced_stim[i] == stim[i:i + history, ...]), 'slicing failed'
 
-def test_getcov():
+def test_cov():
     """Test recovering a stimulus covariance matrix."""
     np.random.seed(0)
     stim = np.random.randn(10, 2)
     assert np.allclose(np.cov(stim.T), stimulustools.cov(stim, 1))
-
-def test_rolling_window_warning():
-    """Verify that calling rolling_window() raises a DeprecationWarning"""
-    with pytest.warns(DeprecationWarning):
-        stimulustools.rolling_window(np.random.randn(20, 3, 3), 5)
-