Permalink
Browse files

pep8...pep8...pep8...

  • Loading branch information...
1 parent 746e477 commit e53dc701c97f258bafd0dcae5abc229b6d007242 @bthirion committed with matthew-brett May 22, 2012
@@ -20,18 +20,18 @@
# --------- Get the data -----------------------------------
#-----------------------------------------------------------
-fmri_files = [example_data.get_filename('fiac', 'fiac0', run)
+fmri_files = [example_data.get_filename('fiac', 'fiac0', run)
for run in ['run1.nii.gz', 'run2.nii.gz']]
-design_files = [example_data.get_filename('fiac', 'fiac0', run)
+design_files = [example_data.get_filename('fiac', 'fiac0', run)
for run in ['run1_design.npz', 'run2_design.npz']]
-mask_file = example_data.get_filename('fiac', 'fiac0', 'mask.nii.gz')
+mask_file = example_data.get_filename('fiac', 'fiac0', 'mask.nii.gz')
affine = load(mask_file).get_affine()
# Get design matrix as numpy array
print('Loading design matrices...')
X = [np.load(f)['X'] for f in design_files]
-# Get multi-session fMRI data
+# Get multi-session fMRI data
print('Loading fmri data...')
Y = [load(f) for f in fmri_files]
@@ -40,28 +40,31 @@
mask = load(mask_file)
mask_array = mask.get_data() > 0
-# GLM fitting
+# GLM fitting
print('Starting fit...')
results = []
for x, y in zip(X, Y):
data = y.get_data()[mask_array].T
mean = data.mean(0)
- data = 100 * (data / mean - 1)
+ data = 100 * (data / mean - 1)
results.append(glm_fit(x, data))
# make a mean volume for display
wmean = mask_array.astype(np.int16)
wmean[mask_array] = mean
+
def make_fiac_contrasts():
"""Specify some constrasts for the FIAC experiment"""
con = {}
# the design matrices of bothe sessions comprise 13 columns
- # the first 5 columns of the design matrices correpond to the following
+ # the first 5 columns of the design matrices correpond to the following
# conditions: ["SSt-SSp", "SSt-DSp", "DSt-SSp", "DSt-DSp", "FirstSt"]
p = 13
+
def length_p_vector(con, p):
return np.hstack((con, np.zeros(p - len(con))))
+
con["SStSSp_minus_DStDSp"] = length_p_vector([1, 0, 0, - 1], p)
con["DStDSp_minus_SStSSp"] = length_p_vector([- 1, 0, 0, 1], p)
con["DSt_minus_SSt"] = length_p_vector([- 1, - 1, 1, 1], p)
@@ -82,7 +85,7 @@ def length_p_vector(con, p):
index + 1, len(contrasts), contrast_id)
contrast_path = op.join(write_dir, '%s_z_map.nii' % contrast_id)
write_array = mask_array.astype(np.float)
- ffx_z_map = (results[0].contrast(contrast_val) +
+ ffx_z_map = (results[0].contrast(contrast_val) +
results[1].contrast(contrast_val)).z_score()
write_array[mask_array] = ffx_z_map
contrast_image = Nifti1Image(write_array, affine)
@@ -95,7 +98,7 @@ def length_p_vector(con, p):
vmin=- vmax,
vmax=vmax,
figure=1,
- threshold=2.5,
+ threshold=2.5,
black_bg=True)
pylab.savefig(op.join(write_dir, '%s_z_map.png' % contrast_id))
pylab.clf()
@@ -175,7 +175,7 @@ def vcov(self, matrix=None, column=None, dispersion=None, other=None):
value(s) for the dispersion parameters
other:array of shape (dim, self.theta.shape[0]), optional
alternative contrast specification (?)
-
+
Returns
=======
cov: array of shape(dim, dim) or (n_voxels, dim, dim),
@@ -257,7 +257,6 @@ def Tcontrast(self, matrix, store=('t', 'effect', 'sd'), dispersion=None):
return TContrastResults(effect=st_effect, t=st_t, sd=st_sd,
df_den=self.df_resid)
-
def Fcontrast(self, matrix, dispersion=None, invcov=None):
"""
Compute an Fcontrast for a contrast matrix.
@@ -288,7 +287,7 @@ def Fcontrast(self, matrix, dispersion=None, invcov=None):
-------
f_res : ``FContrastResults`` instance
with attributes F, df_den, df_num
-
+
Note
----
For F contrasts, we now specify an effect and covariance
@@ -314,7 +313,7 @@ def Fcontrast(self, matrix, dispersion=None, invcov=None):
F = np.squeeze(F)
return FContrastResults(
effect=ctheta, covariance=self.vcov(
- matrix=matrix, dispersion=dispersion[np.newaxis]),
+ matrix=matrix, dispersion=dispersion[np.newaxis]),
F=F, df_den=self.df_resid, df_num=invcov.shape[0])
def conf_int(self, alpha=.05, cols=None, dispersion=None):
@@ -418,4 +417,4 @@ def __array__(self):
def __str__(self):
return '<F contrast: F=%s, df_den=%d, df_num=%d>' % \
- (`self.F`, self.df_den, self.df_num)
+ (repr(self.F), self.df_den, self.df_num)
@@ -6,8 +6,8 @@
It contains the GLM and contrast classes that are meant to be the main objects
of fMRI data analyses.
-It is important to note that the GLM is meant as a one-session
-General Linear Model. But inference can be performed on multiple sessions
+It is important to note that the GLM is meant as a one-session
+General Linear Model. But inference can be performed on multiple sessions
by computing fixed effects on contrasts
@@ -32,21 +32,21 @@
def glm_fit(X, Y, model='ar1', steps=100):
- """ GLM fitting of a dataset using 'ols' regression or the two-pass
+ """ GLM fitting of a dataset using 'ols' regression or the two-pass
Parameters
----------
- X: array of shape(n_time_points, n_regressors),
+ X: array of shape(n_time_points, n_regressors),
the design matrix
Y: array of shape(n_time_points, n_samples), the fMRI data
- model: string, to be chosen in ['ar1', 'ols'], optional,
- the temporal variance model. Defaults to 'ar1'
+ model: string, to be chosen in ['ar1', 'ols'], optional,
+ the temporal variance model. Defaults to 'ar1'
steps: int, optional,
Maximum number of discrete steps for the AR(1) coef histogram
Returns
-------
- result: a GLMResults instance yielding the result of the GLM
+ result: a GLMResults instance yielding the result of the GLM
"""
if model not in ['ar1', 'ols']:
raise ValueError('Unknown model')
@@ -56,7 +56,7 @@ def glm_fit(X, Y, model='ar1', steps=100):
if Y.shape[0] != X.shape[0]:
raise ValueError('Response and predictors are inconsistent')
-
+
# fit the OLS model
ols_result = OLSModel(X).fit(Y)
@@ -93,74 +93,85 @@ def __init__(self, glm_results, labels):
----------
glm_results: dictionary of nipy.fixes.scipy.stats.models.regression,
describing results of a GLM fit
- labels: array of shape(n_voxels),
+ labels: array of shape(n_voxels),
labels that associate each voxel with a results key
"""
if glm_results.keys().sort() != labels.copy().sort():
raise ValueError('results and labels are inconsistant')
+
+ n_voxels = sum([grv.theta.shape[1] for grv in glm_results.values()])
+ if labels.size != n_voxels:
+ raise ValueError('The results do not match the labels size')
+
+ dim = glm_results.values()[0].theta.shape[0]
+ if np.array([grv.theta.shape[0] != dim
+ for grv in glm_results.values()]).any():
+ raise ValueError('The models do not have consistant dimension')
+
self.results = glm_results
self.labels = labels
# fixme: check that the number of voxels corresponds
# to the number of items in glm_results
- def contrast(self, c, contrast_type=None):
+ def contrast(self, con_val, contrast_type=None):
""" Specify and estimate a linear contrast
-
+
Parameters
----------
- c: numpy.ndarray of shape (p) or (q, p),
- where q = number of contrast vectors and p = number of regressors
- contrast_type: string, optional, either 't' or 'F',
+ con_val: numpy.ndarray of shape (p) or (q, p),
+ where q = number of contrast vectors
+ and p = number of regressors
+ contrast_type: string, optional, either 't' or 'F',
type of the contrast
-
+
Returns
-------
con: Contrast instance
"""
- c = np.asarray(c)
- if c.ndim == 1:
+ con_val = np.asarray(con_val)
+ if con_val.ndim == 1:
dim = 1
else:
- dim = c.shape[0]
+ dim = con_val.shape[0]
if contrast_type is None:
if dim == 1:
contrast_type = 't'
else:
contrast_type = 'F'
if contrast_type not in ['t', 'F']:
raise ValueError('Unknown contrast type: %s' % contrast_type)
-
+
effect_ = np.zeros((dim, self.labels.size), dtype=np.float)
var_ = np.zeros((dim, dim, self.labels.size), dtype=np.float)
#stat_ = np.zeros(self.labels.size, dtype=np.float)
if contrast_type == 't':
for l in self.results.keys():
- resl = self.results[l].Tcontrast(c)
+ resl = self.results[l].Tcontrast(con_val)
effect_[:, self.labels == l] = resl.effect.T
var_[:, :, self.labels == l] = (resl.sd ** 2).T
#stat_[self.labels == l] = resl.t
else:
for l in self.results.keys():
- resl = self.results[l].Fcontrast(c)
+ resl = self.results[l].Fcontrast(con_val)
effect_[:, self.labels == l] = resl.effect
var_[:, :, self.labels == l] = resl.covariance
dof_ = self.results[l].df_resid
- return Contrast(effect=effect_, variance=var_, dof=dof_,
+ return Contrast(effect=effect_, variance=var_, dof=dof_,
contrast_type=contrast_type)
class Contrast(object):
""" The contrast class handles the estimation of statistical contrasts
After application of the GLM.
- The important feature is that it supports addition,
+ The important feature is that it supports addition,
thus opening the possibility of fixed-effects models.
-
- The current implementation is meant to be simple,
- and could be enhanced in the future on the computational side
+
+ The current implementation is meant to be simple,
+ and could be enhanced in the future on the computational side
(high-dimensional F constrasts may lead to memory breakage)
"""
- def __init__(self, effect, variance, dof=DEF_DOFMAX, contrast_type='t',
+ def __init__(self, effect, variance, dof=DEF_DOFMAX, contrast_type='t',
tiny=DEF_TINY, dofmax=DEF_DOFMAX):
"""
Parameters
@@ -178,7 +189,7 @@ def __init__(self, effect, variance, dof=DEF_DOFMAX, contrast_type='t',
raise ValueError('Variance array should have 2 dimensions')
if variance.shape[0] != variance.shape[1]:
raise ValueError('Inconsistant shape for the variance estimate')
- if ((variance.shape[1] != effect.shape[0]) or
+ if ((variance.shape[1] != effect.shape[0]) or
(variance.shape[2] != effect.shape[1])):
raise ValueError('Effect and variance have inconsistant shape')
self.effect = effect
@@ -198,7 +209,7 @@ def __init__(self, effect, variance, dof=DEF_DOFMAX, contrast_type='t',
def stat(self, baseline=0.0):
""" Return the decision statistic associated with the test of the
null hypothesis: (H0) 'contrast equals baseline'
-
+
Parameters
==========
baseline: float, optional,
@@ -220,7 +231,8 @@ def stat(self, baseline=0.0):
self.effect = self.effect[np.newaxis]
if self.variance.ndim == 1:
self.variance = self.variance[np.newaxis, np.newaxis]
- stat = mahalanobis(self.effect - baseline, self.variance) / self.dim
+ stat = (mahalanobis(self.effect - baseline, self.variance)
+ / self.dim)
# Case: tmin (conjunctions)
elif self.contrast_type == 'tmin':
vdiag = self.variance.reshape([self.dim ** 2] + list(
@@ -260,7 +272,7 @@ def p_value(self, baseline=0.0):
def z_score(self, baseline=0.0):
"""Return a parametric estimation of the z-score associated
with the null hypothesis: (H0) 'contrast equals baseline'
-
+
Parameters
==========
baseline: float, optional,
@@ -273,6 +285,7 @@ def z_score(self, baseline=0.0):
return zscore(self.p_value_)
def __add__(self, other):
+ """Addition of selfwith others, Yields an new Contrast instance"""
if self.contrast_type != other.contrast_type:
raise ValueError(
'The two contrasts do not have consistant type dimensions')
@@ -282,17 +295,16 @@ def __add__(self, other):
effect_ = self.effect + other.effect
variance_ = self.variance + other.variance
dof_ = self.dof + other.dof
- return Contrast(effect=effect_, variance=variance_, dof=dof_,
+ return Contrast(effect=effect_, variance=variance_, dof=dof_,
contrast_type=self.contrast_type)
-
def __rmul__(self, scalar):
"""Multiplication of the contrast by a scalar"""
scalar = float(scalar)
effect_ = self.effect * scalar
- variance_ = self.variance * scalar ** 2
+ variance_ = self.variance * scalar ** 2
dof_ = self.dof
- return Contrast(effect=effect_, variance=variance_, dof=dof_,
+ return Contrast(effect=effect_, variance=variance_, dof=dof_,
contrast_type=self.contrast_type)
__mul__ = __rmul__
Oops, something went wrong.

0 comments on commit e53dc70

Please sign in to comment.