spectrotemporal

Classes and functions for fitting population encoding models

SpectroTemporalFit(model, data, grids, ...)
SpectroTemporalModel(stimulus) Gaussian population receptive field model.
compute_model_ts(center_freq, sigma, ...)
gaussian_1D(freqs, center_freq, sigma)
gaussian_2D(X, Y, x0, y0, sigma_x, sigma_y, ...)
parallel_fit(args)
recast_estimation_results(output, grid_parent)

PopulationFit

class popeye.spectrotemporal.PopulationFit(model, data, grids, bounds, Ns, voxel_index, auto_fit, verbose)

Bases: object

Base class for all pRF model fits.

__init__(model, data, grids, bounds, Ns, voxel_index, auto_fit, verbose)

A class containing tools for fitting pRF models.

The PopulationFit class houses all the fitting tool that are associated with estimatinga pRF model. The PopulationFit takes a PoulationModel instance model and a time-series data. In addition, extent and sampling-rate of a brute-force grid-search is set with grids and Ns. Use bounds to set limits on the search space for each parameter.

model : AuditoryModel class instance
An object representing the 1D Gaussian model.
data : ndarray
An array containing the measured BOLD signal of a single voxel.
grids : tuple

A tuple indicating the search space for the brute-force grid-search. The tuple contains pairs of upper and lower bounds for exploring a given dimension. For example grids=((-10,10),(0,5),) will search the first dimension from -10 to 10 and the second from 0 to 5. These values cannot be None.

For more information, see scipy.optimize.brute.

bounds : tuple
A tuple containing the upper and lower bounds for each parameter in parameters. If a parameter is not bounded, simply use None. For example, fit_bounds=((0,None),(-10,10),) would bound the first parameter to be any positive number while the second parameter would be bounded between -10 and 10.
Ns : int

Number of samples per stimulus dimension to sample during the ballpark search.

For more information, see scipy.optimize.brute.

voxel_index : tuple
A tuple containing the index of the voxel being modeled. The fitting procedure does not require a voxel index, but collating the results across many voxels will does require voxel indices. With voxel indices, the brain volume can be reconstructed using the newly computed model estimates.
auto_fit : bool
A flag for automatically running the fitting procedures once the GaussianFit object is instantiated.
verbose : int
0 = silent 1 = print the final solution of an error-minimization 2 = print each error-minimization step
Jout()
allvecs()
ballpark()
brute_force()
direc()
estimate()
fopt()
funcalls()
fval()
gradient_descent()
grid()
iter()
msg()
overloaded_estimate()
prediction()
rsquared()
rss()

PopulationModel

class popeye.spectrotemporal.PopulationModel(stimulus, hrf_model, nuissance=None)

Bases: object

Base class for all pRF models.

__init__(stimulus, hrf_model, nuissance=None)

Base class for all pRF models.

stimulus : StimulusModel class instance
An instance of the StimulusModel class containing the stim_arr and various stimulus parameters, and can represent various stimulus modalities, including visual, auditory, etc.
hrf_model : callable
A function that generates an HRF model given an HRF delay. For more information, see popeye.utilties.double_gamma_hrf_hrf and `popeye.utilities.spm_hrf
nuissance : ndarray
A nuissance regressor for removing effects of non-interest. You can regress out any nuissance effects from you data prior to fitting the model of interest. The nuissance model is a statsmodels.OLS compatible design matrix, and the user is expected to have already added any constants.
generate_ballpark_prediction()
generate_prediction()

SpectroTemporalFit

class popeye.spectrotemporal.SpectroTemporalFit(model, data, grids, bounds, Ns, tr_length, voxel_index=(1, 2, 3), auto_fit=True, verbose=0)

Bases: popeye.base.PopulationFit

__init__(model, data, grids, bounds, Ns, tr_length, voxel_index=(1, 2, 3), auto_fit=True, verbose=0)
OLS()
ballpark()
center_freq()
center_freq0()
coefficient()
errmsg()
estimate()
msg()
prediction()
receptive_field()
rsquared()
rss()
sigma()
sigma0()
stderr()

SpectroTemporalModel

class popeye.spectrotemporal.SpectroTemporalModel(stimulus)

Bases: popeye.base.PopulationModel

Gaussian population receptive field model.

__init__(stimulus)

interp1d

class popeye.spectrotemporal.interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=True, fill_value=nan)

Bases: scipy.interpolate.polyint._Interpolator1D

interp1d(x, y, kind=’linear’, axis=-1, copy=True, bounds_error=True,
fill_value=np.nan)

Interpolate a 1-D function.

x and y are arrays of values used to approximate some function f: y = f(x). This class returns a function whose call method uses interpolation to find the value of new points.

x : (N,) array_like
A 1-D array of monotonically increasing real values.
y : (...,N,...) array_like
A N-D array of real values. The length of y along the interpolation axis must be equal to the length of x.
kind : str or int, optional
Specifies the kind of interpolation as a string (‘linear’,’nearest’, ‘zero’, ‘slinear’, ‘quadratic, ‘cubic’) or as an integer specifying the order of the spline interpolator to use. Default is ‘linear’.
axis : int, optional
Specifies the axis of y along which to interpolate. Interpolation defaults to the last axis of y.
copy : bool, optional
If True, the class makes internal copies of x and y. If False, references to x and y are used. The default is to copy.
bounds_error : bool, optional
If True, a ValueError is raised any time interpolation is attempted on a value outside of the range of x (where extrapolation is necessary). If False, out of bounds values are assigned fill_value. By default, an error is raised.
fill_value : float, optional
If provided, then this value will be used to fill in for requested points outside of the data range. If not provided, then the default is NaN.

UnivariateSpline : A more recent wrapper of the FITPACK routines. splrep, splev

Spline interpolation based on FITPACK.

interp2d

>>> from scipy import interpolate
>>> x = np.arange(0, 10)
>>> y = np.exp(-x/3.0)
>>> f = interpolate.interp1d(x, y)

>>> xnew = np.arange(0,9, 0.1)
>>> ynew = f(xnew)   # use interpolation function returned by `interp1d`
>>> plt.plot(x, y, 'o', xnew, ynew, '-')
>>> plt.show()
__init__(x, y, kind='linear', axis=-1, copy=True, bounds_error=True, fill_value=nan)

Initialize a 1D linear interpolation class.

auto_attr

popeye.spectrotemporal.auto_attr(func)

Decorator to create OneTimeProperty attributes.

func : method
The method that will be called the first time to compute a value. Afterwards, the method’s name will be a standard attribute holding the value of this computation.
>>> class MagicProp(object):
...     @auto_attr
...     def a(self):
...         return 99
...
>>> x = MagicProp()
>>> 'a' in x.__dict__
False
>>> x.a
99
>>> 'a' in x.__dict__
True

brute

popeye.spectrotemporal.brute(func, ranges, args=(), Ns=20, full_output=0, finish=<function fmin>, disp=False)

Minimize a function over a given range by brute force.

func : callable f(x,*args)
Objective function to be minimized.
ranges : tuple
Each element is a tuple of parameters or a slice object to be handed to numpy.mgrid.
args : tuple
Extra arguments passed to function.
Ns : int
Default number of samples, if those are not provided.
full_output : bool
If True, return the evaluation grid.
finish : callable, optional
An optimization function that is called with the result of brute force minimization as initial guess. finish should take the initial guess as positional argument, and take take args, full_output and disp as keyword arguments. See Notes for more details.
disp : bool, optional
Set to True to print convergence messages.
x0 : ndarray
Value of arguments to func, giving minimum over the grid.
fval : int
Function value at minimum.
grid : tuple
Representation of the evaluation grid. It has the same length as x0.
Jout : ndarray
Function values over grid: Jout = func(*grid).

The range is respected by the brute force minimization, but if the finish keyword specifies another optimization function (including the default fmin), the returned value may still be (just) outside the range. In order to ensure the range is specified, use finish=None.

compute_model_ts

popeye.spectrotemporal.compute_model_ts(center_freq, sigma, spectrogram, freqs, target_times)

decimate

popeye.spectrotemporal.decimate(x, q, n=None, ftype='iir', axis=-1)

Downsample the signal by using a filter.

By default, an order 8 Chebyshev type I filter is used. A 30 point FIR filter with hamming window is used if ftype is ‘fir’.

x : ndarray
The signal to be downsampled, as an N-dimensional array.
q : int
The downsampling factor.
n : int, optional
The order of the filter (1 less than the length for ‘fir’).
ftype : str {‘iir’, ‘fir’}, optional
The type of the lowpass filter.
axis : int, optional
The axis along which to decimate.
y : ndarray
The down-sampled signal.

resample

fftconvolve

popeye.spectrotemporal.fftconvolve(in1, in2, mode='full')

Convolve two N-dimensional arrays using FFT.

Convolve in1 and in2 using the fast Fourier transform method, with the output size determined by the mode argument.

This is generally much faster than convolve for large arrays (n > ~500), but can be slower when only a few output values are needed, and can only output float arrays (int or object array inputs will be cast to float).

in1 : array_like
First input.
in2 : array_like
Second input. Should have the same number of dimensions as in1; if sizes of in1 and in2 are not equal then in1 has to be the larger array.
mode : str {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
out : array
An N-dimensional array containing a subset of the discrete linear convolution of in1 with in2.

fmin_powell

popeye.spectrotemporal.fmin_powell(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None, direc=None)

Minimize a function using modified Powell’s method. This method only uses function values, not derivatives.

func : callable f(x,*args)
Objective function to be minimized.
x0 : ndarray
Initial guess.
args : tuple, optional
Extra arguments passed to func.
callback : callable, optional
An optional user-supplied function, called after each iteration. Called as callback(xk), where xk is the current parameter vector.
direc : ndarray, optional
Initial direction set.
xtol : float, optional
Line-search error tolerance.
ftol : float, optional
Relative error in func(xopt) acceptable for convergence.
maxiter : int, optional
Maximum number of iterations to perform.
maxfun : int, optional
Maximum number of function evaluations to make.
full_output : bool, optional
If True, fopt, xi, direc, iter, funcalls, and warnflag are returned.
disp : bool, optional
If True, print convergence messages.
retall : bool, optional
If True, return a list of the solution at each iteration.
xopt : ndarray
Parameter which minimizes func.
fopt : number
Value of function at minimum: fopt = func(xopt).
direc : ndarray
Current direction set.
iter : int
Number of iterations.
funcalls : int
Number of function calls made.
warnflag : int
Integer warning flag:
1 : Maximum number of function evaluations. 2 : Maximum number of iterations.
allvecs : list
List of solutions at each iteration.
minimize: Interface to unconstrained minimization algorithms for
multivariate functions. See the ‘Powell’ method in particular.

Uses a modification of Powell’s method to find the minimum of a function of N variables. Powell’s method is a conjugate direction method.

The algorithm has two loops. The outer loop merely iterates over the inner loop. The inner loop minimizes over each current direction in the direction set. At the end of the inner loop, if certain conditions are met, the direction that gave the largest decrease is dropped and replaced with the difference between the current estiamted x and the estimated x from the beginning of the inner-loop.

The technical conditions for replacing the direction of greatest increase amount to checking that

  1. No further gain can be made along the direction of greatest increase from that iteration.
  2. The direction of greatest increase accounted for a large sufficient fraction of the decrease in the function value from that iteration of the inner loop.

Powell M.J.D. (1964) An efficient method for finding the minimum of a function of several variables without calculating derivatives, Computer Journal, 7 (2):155-162.

Press W., Teukolsky S.A., Vetterling W.T., and Flannery B.P.: Numerical Recipes (any edition), Cambridge University Press

gaussian_1D

popeye.spectrotemporal.gaussian_1D(freqs, center_freq, sigma)

gaussian_2D

popeye.spectrotemporal.gaussian_2D(X, Y, x0, y0, sigma_x, sigma_y, degrees, amplitude=1)

generate_rf_timeseries_1D

popeye.spectrotemporal.generate_rf_timeseries_1D()

linregress

popeye.spectrotemporal.linregress(x, y=None)

Calculate a regression line

This computes a least-squares regression for two sets of measurements.

x, y : array_like
two sets of measurements. Both arrays should have the same length. If only x is given (and y=None), then it must be a two-dimensional array where one dimension has length 2. The two sets of measurements are then found by splitting the array along the length-2 dimension.
slope : float
slope of the regression line
intercept : float
intercept of the regression line
r-value : float
correlation coefficient
p-value : float
two-sided p-value for a hypothesis test whose null hypothesis is that the slope is zero.
stderr : float
Standard error of the estimate
>>> from scipy import stats
>>> import numpy as np
>>> x = np.random.random(10)
>>> y = np.random.random(10)
>>> slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)

# To get coefficient of determination (r_squared)

>>> print "r-squared:", r_value**2
r-squared: 0.15286643777

median_filter

popeye.spectrotemporal.median_filter(input, size=None, footprint=None, output=None, mode='reflect', cval=0.0, origin=0)

Calculates a multidimensional median filter.

input : array_like
Input array to filter.
size : scalar or tuple, optional
See footprint, below
footprint : array, optional
Either size or footprint must be defined. size gives the shape that is taken from the input array, at every element position, to define the input to the filter function. footprint is a boolean array that specifies (implicitly) a shape, but also which of the elements within this shape will get passed to the filter function. Thus size=(n,m) is equivalent to footprint=np.ones((n,m)). We adjust size to the number of dimensions of the input array, so that, if the input array is shape (10,10,10), and size is 2, then the actual size used is (2,2,2).
output : array, optional
The output parameter passes an array in which to store the filter output.
mode : {‘reflect’, ‘constant’, ‘nearest’, ‘mirror’, ‘wrap’}, optional
The mode parameter determines how the array borders are handled, where cval is the value when mode is equal to ‘constant’. Default is ‘reflect’
cval : scalar, optional
Value to fill past edges of input if mode is ‘constant’. Default is 0.0
origin : scalar, optional
The origin parameter controls the placement of the filter. Default 0.0.
median_filter : ndarray
Return of same shape as input.

parallel_fit

popeye.spectrotemporal.parallel_fit(args)

recast_estimation_results

popeye.spectrotemporal.recast_estimation_results(output, grid_parent, write=True)

romb

popeye.spectrotemporal.romb(y, dx=1.0, axis=-1, show=False)

Romberg integration using samples of a function.

y : array_like
A vector of 2**k + 1 equally-spaced samples of a function.
dx : array_like, optional
The sample spacing. Default is 1.
axis : array_like?, optional
The axis along which to integrate. Default is -1 (last axis).
show : bool, optional
When y is a single 1-D array, then if this argument is True print the table showing Richardson extrapolation from the samples. Default is False.
ret : array_like?
The integrated result for each axis.

quad - adaptive quadrature using QUADPACK romberg - adaptive Romberg quadrature quadrature - adaptive Gaussian quadrature fixed_quad - fixed-order Gaussian quadrature dblquad, tplquad - double and triple integrals simps, trapz - integrators for sampled data cumtrapz - cumulative integration for sampled data ode, odeint - ODE integrators

trapz

popeye.spectrotemporal.trapz(y, x=None, dx=1.0, axis=-1)

Integrate along the given axis using the composite trapezoidal rule.

Integrate y (x) along given axis.

y : array_like
Input array to integrate.
x : array_like, optional
The sample points corresponding to the y values. If x is None, the sample points are assumed to be evenly spaced dx apart. The default is None.
dx : scalar, optional
The spacing between sample points when x is None. The default is 1.
axis : int, optional
The axis along which to integrate.
trapz : float
Definite integral as approximated by trapezoidal rule.

sum, cumsum

Image [2] illustrates trapezoidal rule – y-axis locations of points will be taken from y array, by default x-axis distances between points will be 1.0, alternatively they can be provided with x array or with dx scalar. Return value will be equal to combined area under the red lines.

[1]Wikipedia page: http://en.wikipedia.org/wiki/Trapezoidal_rule
[2]Illustration image: http://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
 
>>> np.trapz([1,2,3])
4.0
>>> np.trapz([1,2,3], x=[4,6,8])
8.0
>>> np.trapz([1,2,3], dx=2)
8.0
>>> a = np.arange(6).reshape(2, 3)
>>> a
array([[0, 1, 2],
       [3, 4, 5]])
>>> np.trapz(a, axis=0)
array([ 1.5,  2.5,  3.5])
>>> np.trapz(a, axis=1)
array([ 2.,  8.])