strf

Classes and functions for fitting Gaussian population encoding models

SpatioTemporalFit(model, data, grids, ...[, ...]) A spatiotemporal population receptive field fit class
SpatioTemporalModel(stimulus) A Gaussian population receptive field model class
compute_model_ts(x, y, spatial_sigma, ...) The objective function for GaussianFi class.
parallel_fit(args) This is a convenience function for parallelizing the fitting procedure.
recast_estimation_results(output, grid_parent) Recasts the output of the prf estimation into two nifti_gz volumes.

PopulationFit

class popeye.strf.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.strf.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()

SpatioTemporalFit

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

Bases: popeye.base.PopulationFit

A spatiotemporal population receptive field fit class

__init__(model, data, grids, bounds, Ns, tr_length, voxel_index=(1, 2, 3), auto_fit=True, verbose=0)
data : ndarray
An array containing the measured BOLD signal.
model : GaussianModel class instance containing the representation
of the visual stimulus.
search_bounds : 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 fit_bounds=((-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.

fit_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.
tr_length : float
The length of the repetition time in seconds.
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 : bool
A flag for printing some summary information about the model estiamte after the fitting procedures have completed.
OLS()
baseline()
baseline0()
beta()
beta0()
coefficient()
generate_prediction(x, y, spatial_sigma, temporal_sigma, weight, beta, baseline)
hemodynamic_response()
msg()
prediction()
response_center()
rho()
rsquared()
rss()
spatial_rf()
spatial_rf_norm()
spatial_s0()
spatial_sigma()
stderr()
sustained_response()
sustained_response_norm()
temporal_s0()
temporal_sigma()
theta()
tr_timescale()
transient_response()
transient_response_norm()
weight()
weight0()
x()
x0()
y()
y0()

SpatioTemporalModel

class popeye.strf.SpatioTemporalModel(stimulus)

Bases: popeye.base.PopulationModel

A Gaussian population receptive field model class

__init__(stimulus)

A spatiotemporal population receptive field model.

stimulus : VisualStimulus class object
A class instantiation of the VisualStimulus class containing a representation of the visual stimulus.

auto_attr

popeye.strf.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

compute_model_ts

popeye.strf.compute_model_ts(x, y, spatial_sigma, temporal_sigma, weight, beta, baseline, deg_x, deg_y, stim_arr, tr_length, projector_hz)

The objective function for GaussianFi class.

x : float
The model estimate along the horizontal dimensions of the display.
y : float
The model estimate along the vertical dimensions of the display.
sigma : float
The model estimate of the dispersion across the the display.
hrf_delay : float
The model estimate of the relative delay of the HRF. The canonical HRF is assumed to be 5 s post-stimulus [1]_.
beta : float
The model estimate of the amplitude of the BOLD signal.
tr_length : float
The length of the repetition time in seconds.

model : ndarray The model prediction time-series.

[1]Glover, GH. (1999). Deconvolution of impulse response in

event-related BOLD fMRI. NeuroImage 9: 416-429.

fftconvolve

popeye.strf.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.

generate_og_receptive_field

popeye.strf.generate_og_receptive_field()

Generate a Gaussian.

x : float
The x coordinate of the center of the Gaussian (degrees)
y : float
The y coordinate of the center of the Gaussian (degrees)
s : float
The dispersion of the Gaussian (degrees)
beta : float
The amplitude of the Gaussian
deg_x : 2D array
The coordinate matrix along the horizontal dimension of the display (degrees)
deg_y : 2D array
The coordinate matrix along the vertical dimension of the display (degrees)

Returns

stim : ndarray
The 1D array containing the stimulus energies given the Gaussian coordinates

generate_rf_timeseries

popeye.strf.generate_rf_timeseries()

linregress

popeye.strf.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

parallel_fit

popeye.strf.parallel_fit(args)

This is a convenience function for parallelizing the fitting procedure. Each call is handed a tuple or list containing all the necessary inputs for instantiaing a GaussianFit class object and estimating the model parameters.

args : list/tuple
A list or tuple containing all the necessary inputs for fitting the Gaussian pRF model.
fit : GaussianFit class object
A fit object that contains all the inputs and outputs of the Gaussian pRF model estimation for a single voxel.

recast_estimation_results

popeye.strf.recast_estimation_results(output, grid_parent, polar=False)

Recasts the output of the prf estimation into two nifti_gz volumes.

Takes output, a list of multiprocessing.Queue objects containing the output of the prf estimation for each voxel. The prf estimates are expressed in both polar and Cartesian coordinates. If the default value for the write parameter is set to False, then the function returns the arrays without writing the nifti files to disk. Otherwise, if write is True, then the two nifti files are written to disk.

Each voxel contains the following metrics:

0 x / polar angle 1 y / eccentricity 2 sigma 3 HRF delay 4 RSS error of the model fit 5 correlation of the model fit
output : list
A list of PopulationFit objects.
grid_parent : nibabel object
A nibabel object to use as the geometric basis for the statmap. The grid_parent (x,y,z) dim and pixdim will be used.
cartes_filename : string
The absolute path of the recasted prf estimation output in Cartesian coordinates.
plar_filename : string
The absolute path of the recasted prf estimation output in polar coordinates.

simps

popeye.strf.simps(y, x=None, dx=1, axis=-1, even='avg')

Integrate y(x) using samples along the given axis and the composite Simpson’s rule. If x is None, spacing of dx is assumed.

If there are an even number of samples, N, then there are an odd number of intervals (N-1), but Simpson’s rule requires an even number of intervals. The parameter ‘even’ controls how this is handled.

y : array_like
Array to be integrated.
x : array_like, optional
If given, the points at which y is sampled.
dx : int, optional
Spacing of integration points along axis of y. Only used when x is None. Default is 1.
axis : int, optional
Axis along which to integrate. Default is the last axis.
even : {‘avg’, ‘first’, ‘str’}, optional
‘avg’ : Average two results:1) use the first N-2 intervals with
a trapezoidal rule on the last interval and 2) use the last N-2 intervals with a trapezoidal rule on the first interval.
‘first’ : Use Simpson’s rule for the first N-2 intervals with
a trapezoidal rule on the last interval.
‘last’ : Use Simpson’s rule for the last N-2 intervals with a
trapezoidal rule on the first interval.

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

For an odd number of samples that are equally spaced the result is exact if the function is a polynomial of order 3 or less. If the samples are not equally spaced, then the result is exact only if the function is a polynomial of order 2 or less.

trapz

popeye.strf.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.])