Commit 7ba5756f authored by Paul McCarthy's avatar Paul McCarthy
Browse files

Unit tests for featdesign module. Still need to test voxelwise EVs

parent 2b5d433f
#!/usr/bin/env python
#
# test_featdesign.py -
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""Test data sets (in testdata/test_feat) were generated with FSL 5.0.9, and
then 'cleaned' to remove unnecessary files and reduce size, with the following
commands:
First level analyses:
find . -name "*png" -delete
find . -name "*ppm" -delete
find . -name "*html" -delete
rm -r logs .files tsplot .ramp.gif
for f in `find . -name "*nii.gz"`; do
echo $f > $f
done
Second level analyses:
find . -name "*png" -delete
find . -name "*ppm" -delete
find . -name "*html" -delete
rm -r logs .files .ramp.gif
for f in `find . -name "*nii.gz"`; do
echo $f > $f
done
Second level cope?.feats:
rm -r logs .files tsplot .ramp.gif
`1stlevel_1.feat`
- 45 time points
- 10 EVs in total:
- 2 stimulus EVs
- 2 temporal derivative EVs
- 6 Standard motion parameters
- 2 contrasts - one on each stimulus EV
`1stlevel_2.feat`
- 45 time points
- 11 EVs in total:
- 2 stimulus EVs
- 2 temporal derivative EVs
- 1 voxelwise EV
- 6 Standard motion parameters
- 2 contrasts - one on each stimulus EV
`1stlevel_3.feat`
- 45 time points
- 32 EVs in total:
- 2 stimulus EVs
- 1 Temporal derivative EV on first
- 2 gamma basis functions on second
- 24 Standard+extended motion parameters
- 2 Confound EVs
- 1 Voxelwise confound EV
- 2 contrasts - one on each stimulus EV
`2ndlevel_1.feat`
- Three inputs
- Two copes
- One main EV - group average
`2ndlevel_1.feat`
- Three inputs
- Two copes
- One main EV - group average
- One voxelwise EV
"""
import os.path as op
import pytest
import fsl.data.featdesign as featdesign
import fsl.data.featanalysis as featanalysis
datadir = op.join(op.dirname(__file__), 'testdata', 'test_feat')
def test_FeatFSFDesign():
featdirs = ['1stlevel_1.feat', '1stlevel_2.feat', '1stlevel_3.feat',
'2ndlevel_1.gfeat', '2ndlevel_2.gfeat']
nevs = [10, 11, 32, 1, 2]
shapes = [(45, 10), (45, 11), (45, 32), (3, 1), (3, 2)]
for featdir, nev, shape in zip(featdirs, nevs, shapes):
featdir = op.join(datadir, featdir)
settings = featanalysis.loadSettings(featdir)
# We can't load the voxelwise EVs
# here, because all of the .nii.gz
# files in the test directory are
# stubs.
#
# See the test_firstLevelVoxelwiseEV
# function
des = featdesign.FEATFSFDesign(featdir,
loadVoxelwiseEVs=False)
# Can also specify the design.fsf settings
featdesign.FEATFSFDesign(featdir,
settings=settings,
loadVoxelwiseEVs=False)
assert len(des.getEVs()) == nev
assert des.getDesign().shape == shape
assert des.getDesign((10, 10, 3)).shape == shape
def test_getFirstLevelEVs_1():
featdir = op.join(datadir, '1stlevel_1.feat')
settings = featanalysis.loadSettings(featdir)
matrix = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))
expected = [(featdesign.NormalEV, {'index' : 0, 'origIndex' : 0}),
(featdesign.TemporalDerivativeEV, {'index' : 1}),
(featdesign.NormalEV, {'index' : 2, 'origIndex' : 1}),
(featdesign.TemporalDerivativeEV, {'index' : 3}),
(featdesign.MotionParameterEV, {'index' : 4, 'motionIndex' : 0}),
(featdesign.MotionParameterEV, {'index' : 5, 'motionIndex' : 1}),
(featdesign.MotionParameterEV, {'index' : 6, 'motionIndex' : 2}),
(featdesign.MotionParameterEV, {'index' : 7, 'motionIndex' : 3}),
(featdesign.MotionParameterEV, {'index' : 8, 'motionIndex' : 4}),
(featdesign.MotionParameterEV, {'index' : 9, 'motionIndex' : 5})]
evs = featdesign.getFirstLevelEVs(featdir, settings, matrix)
assert len(evs) == 10
for i, (evtype, atts) in enumerate(expected):
assert isinstance(evs[i], evtype)
for k, v in atts.items():
assert getattr(evs[i], k) == v
def test_getFirstLevelEVs_2():
featdir = op.join(datadir, '1stlevel_2.feat')
settings = featanalysis.loadSettings(featdir)
matrix = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))
expected = [(featdesign.NormalEV, {'index' : 0, 'origIndex' : 0}),
(featdesign.TemporalDerivativeEV, {'index' : 1}),
(featdesign.NormalEV, {'index' : 2, 'origIndex' : 1}),
(featdesign.TemporalDerivativeEV, {'index' : 3}),
(featdesign.VoxelwiseEV, {'index' : 4, 'origIndex' : 2}),
(featdesign.MotionParameterEV, {'index' : 5, 'motionIndex' : 0}),
(featdesign.MotionParameterEV, {'index' : 6, 'motionIndex' : 1}),
(featdesign.MotionParameterEV, {'index' : 7, 'motionIndex' : 2}),
(featdesign.MotionParameterEV, {'index' : 8, 'motionIndex' : 3}),
(featdesign.MotionParameterEV, {'index' : 9, 'motionIndex' : 4}),
(featdesign.MotionParameterEV, {'index' : 10, 'motionIndex' : 5})]
evs = featdesign.getFirstLevelEVs(featdir, settings, matrix)
assert len(evs) == 11
for i, (evtype, atts) in enumerate(expected):
assert isinstance(evs[i], evtype)
for k, v in atts.items():
assert getattr(evs[i], k) == v
def test_getFirstLevelEVs_3():
featdir = op.join(datadir, '1stlevel_3.feat')
settings = featanalysis.loadSettings(featdir)
matrix = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))
expected = [(featdesign.NormalEV, {'index' : 0, 'origIndex' : 0}),
(featdesign.TemporalDerivativeEV, {'index' : 1}),
(featdesign.NormalEV, {'index' : 2, 'origIndex' : 1}),
(featdesign.BasisFunctionEV, {'index' : 3}),
(featdesign.BasisFunctionEV, {'index' : 4}),
(featdesign.VoxelwiseConfoundEV, {'index' : 5, 'voxIndex' : 0}),
(featdesign.MotionParameterEV, {'index' : 6, 'motionIndex' : 0}),
(featdesign.MotionParameterEV, {'index' : 7, 'motionIndex' : 1}),
(featdesign.MotionParameterEV, {'index' : 8, 'motionIndex' : 2}),
(featdesign.MotionParameterEV, {'index' : 9, 'motionIndex' : 3}),
(featdesign.MotionParameterEV, {'index' : 10, 'motionIndex' : 4}),
(featdesign.MotionParameterEV, {'index' : 11, 'motionIndex' : 5}),
(featdesign.MotionParameterEV, {'index' : 12, 'motionIndex' : 6}),
(featdesign.MotionParameterEV, {'index' : 13, 'motionIndex' : 7}),
(featdesign.MotionParameterEV, {'index' : 14, 'motionIndex' : 8}),
(featdesign.MotionParameterEV, {'index' : 15, 'motionIndex' : 9}),
(featdesign.MotionParameterEV, {'index' : 16, 'motionIndex' : 10}),
(featdesign.MotionParameterEV, {'index' : 17, 'motionIndex' : 11}),
(featdesign.MotionParameterEV, {'index' : 18, 'motionIndex' : 12}),
(featdesign.MotionParameterEV, {'index' : 19, 'motionIndex' : 13}),
(featdesign.MotionParameterEV, {'index' : 20, 'motionIndex' : 14}),
(featdesign.MotionParameterEV, {'index' : 21, 'motionIndex' : 15}),
(featdesign.MotionParameterEV, {'index' : 22, 'motionIndex' : 16}),
(featdesign.MotionParameterEV, {'index' : 23, 'motionIndex' : 17}),
(featdesign.MotionParameterEV, {'index' : 24, 'motionIndex' : 18}),
(featdesign.MotionParameterEV, {'index' : 25, 'motionIndex' : 19}),
(featdesign.MotionParameterEV, {'index' : 26, 'motionIndex' : 20}),
(featdesign.MotionParameterEV, {'index' : 27, 'motionIndex' : 21}),
(featdesign.MotionParameterEV, {'index' : 28, 'motionIndex' : 22}),
(featdesign.MotionParameterEV, {'index' : 29, 'motionIndex' : 23}),
(featdesign.ConfoundEV, {'index' : 30, 'confIndex' : 0}),
(featdesign.ConfoundEV, {'index' : 31, 'confIndex' : 1})]
evs = featdesign.getFirstLevelEVs(featdir, settings, matrix)
assert len(evs) == 32
for i, (evtype, atts) in enumerate(expected):
print(i, evs[i])
assert isinstance(evs[i], evtype)
for k, v in atts.items():
assert getattr(evs[i], k) == v
def test_getHigherLevelEVs_1():
featdir = op.join(datadir, '2ndlevel_1.gfeat')
settings = featanalysis.loadSettings(featdir)
matrix = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))
evs = featdesign.getHigherLevelEVs(featdir, settings, matrix)
assert len(evs) == 1
assert isinstance(evs[0], featdesign.NormalEV)
assert evs[0].index == 0
assert evs[0].origIndex == 0
def test_getHigherLevelEVs_2():
featdir = op.join(datadir, '2ndlevel_2.gfeat')
settings = featanalysis.loadSettings(featdir)
matrix = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))
evs = featdesign.getHigherLevelEVs(featdir, settings, matrix)
assert len(evs) == 2
assert isinstance(evs[0], featdesign.NormalEV)
assert evs[0].index == 0
assert evs[0].origIndex == 0
assert isinstance(evs[1], featdesign.VoxelwiseEV)
def test_loadDesignMat():
analyses = ['1stlevel_1.feat', '1stlevel_2.feat', '1stlevel_3.feat',
'2ndlevel_1.gfeat', '2ndlevel_2.gfeat']
shapes = [(45, 10), (45, 11), (45, 32), (3, 1), (3, 2)]
for analysis, shape in zip(analyses, shapes):
featdir = op.join(datadir, analysis)
fname = op.join(featdir, 'design.mat')
mat = featdesign.loadDesignMat(fname)
assert mat.shape == shape
nonfile = op.join(datadir, 'non-existent-file')
badfile = op.join(datadir, '1stlevel_1.feat', 'design.fsf')
with pytest.raises(Exception):
featdesign.loadDesignMat(nonfile)
with pytest.raises(Exception):
featdesign.loadDesignMat(badfile)
Cluster Index Voxels P -log10(P) Z-MAX Z-MAX X (vox) Z-MAX Y (vox) Z-MAX Z (vox) Z-COG X (vox) Z-COG Y (vox) Z-COG Z (vox) COPE-MAX COPE-MAX X (vox) COPE-MAX Y (vox) COPE-MAX Z (vox) COPE-MEAN
1 296 1.79e-27 26.7 6.03 34 10 1 31.4 12.3 1.72 612 34 10 1 143
Cluster Index Voxels P -log10(P) Z-MAX Z-MAX X (mm) Z-MAX Y (mm) Z-MAX Z (mm) Z-COG X (mm) Z-COG Y (mm) Z-COG Z (mm) COPE-MAX COPE-MAX X (mm) COPE-MAX Y (mm) COPE-MAX Z (mm) COPE-MEAN
1 296 1.79e-27 26.7 6.03 -8.41 -89.9 -2.42 3.01 -81.2 2.15 612 -8.41 -89.9 -2.42 143
Cluster Index Voxels P -log10(P) Z-MAX Z-MAX X (vox) Z-MAX Y (vox) Z-MAX Z (vox) Z-COG X (vox) Z-COG Y (vox) Z-COG Z (vox) COPE-MAX COPE-MAX X (vox) COPE-MAX Y (vox) COPE-MAX Z (vox) COPE-MEAN
5 271 7.77e-26 25.1 6.27 47 27 2 44.3 29.7 1.72 468 46 25 3 113
4 169 1.53e-18 17.8 5.75 20 32 2 23.9 28.4 1.45 304 19 31 1 100
3 24 0.000341 3.47 4.58 17 39 3 18.8 38.3 2.23 157 22 39 0 62
2 19 0.00202 2.69 4.88 25 44 0 26.3 42.1 0 197 26 43 0 118
1 13 0.0213 1.67 3.73 26 19 2 26.6 21.3 1.85 97.1 26 19 2 69.2
Cluster Index Voxels P -log10(P) Z-MAX Z-MAX X (mm) Z-MAX Y (mm) Z-MAX Z (mm) Z-COG X (mm) Z-COG Y (mm) Z-COG Z (mm) COPE-MAX COPE-MAX X (mm) COPE-MAX Y (mm) COPE-MAX Z (mm) COPE-MEAN
5 271 7.77e-26 25.1 6.27 -63.6 -21.8 4 -51.3 -11.2 1.84 468 -59.6 -29.9 10.8 113
4 169 1.53e-18 17.8 5.75 55.4 -4.27 2.45 37.9 -18.1 -0.842 304 59.8 -8.19 -4.25 100
3 24 0.000341 3.47 4.58 69.3 23 8.65 61.3 20.6 3.62 157 47.6 23.8 -11.3 62
2 19 0.00202 2.69 4.88 35 43.8 -11.4 29.3 36.5 -11.2 197 30.5 39.9 -11.3 118
1 13 0.0213 1.67 3.73 27.6 -55.1 3.43 25.3 -46.1 2.31 97.1 27.6 -55.1 3.43 69.2
0 -0.000264589 0.000183168 0.0217952 -0.0424268 -0.0254739
0.00070214 -0.000526306 0.000397315 0.0337431 -0.0193044 -0.0283846
0.000498598 -0.000809397 0 0.0395197 -0.0131337 -0.0452989
0.000579988 -0.000221715 8.21135e-05 0.0258738 0.0348006 -0.0129759
0.000381893 -0.000506425 0 0.0167397 0.0709324 -0.0453051
0.000388283 -0.000809397 0 0.0215977 -0.00161064 -0.0651091
0.00028076 -0.000604851 0 0.0241783 -0.0191748 -0.0536596
0 -0.000155116 0 0.0128543 -0.0191748 -0.0452136
0.000384396 -7.84812e-05 0 0.00734347 0.021921 -0.0415896
0.000321737 -0.000393791 0 0.00375303 -0.0230888 -0.0415753
0 0.000111998 0 -0.000902616 0.0727951 -0.0415698
0.000215285 -0.000136386 0 -0.00282623 -0.0354858 -0.0415777
0.000261125 -9.06374e-05 0 0.000629449 -0.0731166 -0.0415577
0.000486193 -0.000133878 0 0.00314586 0.000714271 -0.0452225
5.40494e-05 -2.19574e-05 0 0.00489143 -0.0056629 0.0214406
0.000459199 0.0004553 -0.000374243 -0.031751 -0.075842 0.0214637
0.000347496 0.000667638 -0.000316943 -0.0241288 -0.0550503 0.052287
0.000153477 -2.19574e-05 -0.000196968 -0.0113862 -0.0806421 0.0129587
0.000459199 0.000326759 -0.000257516 -0.018521 -0.0673269 0.0310445
-0.00010422 -0.000709562 -0.000551959 -0.0343409 -0.0673334 0.0190272
-4.78084e-05 -0.000709562 -0.000167275 -0.00978313 -0.0746766 -0.00060184
0.000459199 -0.000265855 -8.30117e-05 -0.00650845 -0.0269721 -0.00502952
-7.98919e-06 -0.000265855 9.6357e-06 -0.00184664 0.000679341 0.000189414
8.2695e-05 7.07945e-05 9.6357e-06 0.00353303 -0.0269518 -0.00502371
0.000151626 -0.000265855 -9.28369e-05 -0.00184551 -0.0269531 -0.012929
0.000381967 -0.000265855 -4.45114e-07 0.00220544 -0.0269575 -0.00503093
0.000198744 0.000212547 -4.45114e-07 -0.00184187 -0.0334269 0.0113516
-7.04892e-05 -0.000265855 -0.000111352 -0.00183592 -0.0434752 -0.0156616
0.000118964 -0.000265855 -4.45114e-07 -0.00184005 -0.0223293 -0.00503604
0.000458349 -0.000265855 -4.45114e-07 0.0217036 -0.00490347 -0.00509171
0.000571088 0.000180565 -0.000332005 0.0027614 -0.0223381 -0.00508851
0.000362825 0.000104016 -4.45114e-07 -0.0050014 -0.0505434 -0.00502399
0.00036589 -7.31897e-05 3.36877e-05 0.00123403 -0.0879296 -0.0285882
0.000133338 -0.000211688 -0.000144184 -0.00941317 -0.0856607 -0.00688414
0.000177832 0.000215669 -0.000185755 -0.0128922 -0.0856624 0.0114737
-0.000130942 -8.65635e-05 -1.11437e-05 -0.0021262 -0.0156379 -0.0286166
4.3281e-05 0.000387841 -0.00013645 -0.00962933 -0.0459189 0.0222458
-2.70728e-05 0.000729944 -0.000153679 -0.0179646 -0.0459249 0.0222224
4.3281e-05 0.000247161 -0.000452448 -0.0306068 -0.0621157 0.00586292
0.000128363 0.000237582 -0.00027413 -0.0180748 -0.065616 -0.00337989
4.3281e-05 0.000239086 -0.000195342 -0.018503 -0.0773948 -0.0033793
0.00038521 0.000105191 -0.000219467 -0.0197604 -0.0525178 -0.00335665
0.000186606 -0.000411083 -0.000184065 -0.0228653 -0.059242 -0.00333508
0.000446345 -0.000305936 -0.000125382 -0.0135317 -0.0429802 0.00123055
0.000197933 -0.000176617 -0.000408 -0.0193382 -0.0681025 1.65356e-05
/ContrastName1 auditory
/ContrastName2 visual
/NumWaves 10
/NumContrasts 2
/PPheights 1.470856e+00 1.389712e+00
/RequiredEffect 3.840 3.963
/Matrix
1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
# FEAT version number
set fmri(version) 6.00
# Are we in MELODIC?
set fmri(inmelodic) 0
# Analysis level
# 1 : First-level analysis
# 2 : Higher-level analysis
set fmri(level) 1
# Which stages to run
# 0 : No first-level analysis (registration and/or group stats only)
# 7 : Full first-level analysis
# 1 : Pre-processing
# 2 : Statistics
set fmri(analysis) 7
# Use relative filenames
set fmri(relative_yn) 0
# Balloon help
set fmri(help_yn) 1
# Run Featwatcher
set fmri(featwatcher_yn) 1
# Cleanup first-level standard-space images
set fmri(sscleanup_yn) 0
# Output directory
set fmri(outputdir) "/Users/paulmc/Projects/fslpy/tests/testdata/test_feat/fmri"
# TR(s)
set fmri(tr) 3.000000
# Total volumes
set fmri(npts) 45
# Delete volumes
set fmri(ndelete) 0
# Perfusion tag/control order
set fmri(tagfirst) 1
# Number of first-level analyses
set fmri(multiple) 3
# Higher-level input type
# 1 : Inputs are lower-level FEAT directories
# 2 : Inputs are cope images from FEAT directories
set fmri(inputtype) 2
# Carry out pre-stats processing?
set fmri(filtering_yn) 1
# Brain/background threshold, %
set fmri(brain_thresh) 10
# Critical z for design efficiency calculation
set fmri(critical_z) 5.3
# Noise level
set fmri(noise) 0.66
# Noise AR(1)
set fmri(noisear) 0.34
# Motion correction
# 0 : None
# 1 : MCFLIRT
set fmri(mc) 1
# Spin-history (currently obsolete)
set fmri(sh_yn) 0
# B0 fieldmap unwarping?
set fmri(regunwarp_yn) 0
# EPI dwell time (ms)
set fmri(dwell) 0.7
# EPI TE (ms)
set fmri(te) 35
# % Signal loss threshold
set fmri(signallossthresh) 10
# Unwarp direction
set fmri(unwarp_dir) y-
# Slice timing correction
# 0 : None
# 1 : Regular up (0, 1, 2, 3, ...)
# 2 : Regular down
# 3 : Use slice order file
# 4 : Use slice timings file
# 5 : Interleaved (0, 2, 4 ... 1, 3, 5 ... )
set fmri(st) 0
# Slice timings file
set fmri(st_file) ""
# BET brain extraction
set fmri(bet_yn) 1
# Spatial smoothing FWHM (mm)
set fmri(smooth) 5
# Intensity normalization
set fmri(norm_yn) 0
# Perfusion subtraction
set fmri(perfsub_yn) 0
# Highpass temporal filtering
set fmri(temphp_yn) 1
# Lowpass temporal filtering
set fmri(templp_yn) 0
# MELODIC ICA data exploration
set fmri(melodic_yn) 0
# Carry out main stats?
set fmri(stats_yn) 1
# Carry out prewhitening?
set fmri(prewhiten_yn) 1
# Add motion parameters to model
# 0 : No
# 1 : Yes
set fmri(motionevs) 1
set fmri(motionevsbeta) ""
set fmri(scriptevsbeta) ""
# Robust outlier detection in FLAME?
set fmri(robust_yn) 0
# Higher-level modelling
# 3 : Fixed effects
# 0 : Mixed Effects: Simple OLS
# 2 : Mixed Effects: FLAME 1
# 1 : Mixed Effects: FLAME 1+2
set fmri(mixed_yn) 2
# Number of EVs
set fmri(evs_orig) 2
set fmri(evs_real) 4
set fmri(evs_vox) 0
# Number of contrasts
set fmri(ncon_orig) 2
set fmri(ncon_real) 2
# Number of F-tests
set fmri(nftests_orig) 0
set fmri(nftests_real) 0
# Add constant column to design matrix? (obsolete)
set fmri(constcol) 0
# Carry out post-stats steps?
set fmri(poststats_yn) 1
# Pre-threshold masking?
set fmri(threshmask) ""
# Thresholding
# 0 : None
# 1 : Uncorrected
# 2 : Voxel
# 3 : Cluster
set fmri(thresh) 3
# P threshold
set fmri(prob_thresh) 0.05
# Z threshold
set fmri(z_thresh) 2.3
# Z min/max for colour rendering
# 0 : Use actual Z min/max
# 1 : Use preset Z min/max
set fmri(zdisplay) 0
# Z min in colour rendering
set fmri(zmin) 2
# Z max in colour rendering
set fmri(zmax) 8
# Colour rendering type
# 0 : Solid blobs
# 1 : Transparent blobs
set fmri(rendertype) 1
# Background image for higher-level stats overlays
# 1 : Mean highres
# 2 : First highres
# 3 : Mean functional
# 4 : First functional
# 5 : Standard space template
set fmri(bgimage) 1
# Create time series plots
set fmri(tsplot_yn) 1
# Registration to initial structural
set fmri(reginitial_highres_yn) 1
# Search space for registration to initial structural
# 0 : No search
# 90 : Normal search
# 180 : Full search
set fmri(reginitial_highres_search) 90
# Degrees of Freedom for registration to initial structural
set fmri(reginitial_highres_dof) 3
# Registration to main structural
set fmri(reghighres_yn) 1
# Search space for registration to main structural
# 0 : No search
# 90 : Normal search
# 180 : Full search
set fmri(reghighres_search) 90
# Degrees of Freedom for registration to main structural
set fmri(reghighres_dof) BBR