#!/usr/bin/env python
#
# featdesign.py -
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module provides the :class:`FEATFSFDesign` class, which encapsulates
a FEAT design matrix.
The FEAT design matrix
----------------------
A FEAT design matrix may contain the following types of explanatory variables:
- *Normal* EVs. This is simply a column in the design matrix, defined by the
user.
- Temporal derivative of normal EVs. A column in the design matrix containing
the derivative of the previous normal EV. The presence of a temporal
derivative EV for a given normal EV can be determined by the ``deriv_yn``
flag in the ``design.fsf`` file.
- Basis function EV. One or more columns derived from a normal EV
(``basisfnumN``)
- Voxelwise EVs (``designVoxelwiseEV<N>.nii.gz``, EV number must be offset by
temporal derivative EVs)
- Confound EVs (... ?)
- Voxelwise confound EVs (``vef.dat`` and ``ven.dat``)
- Motion parameter EVs (``design.fsf:motionevs``, the last 6 or 24 columns of
the design matrix [i think])
For each voxelwise EV, the design matrix (in ``design.mat``) contains a
'dummy' column which contains the mean across all voxels.
For voxelwise EVs, the column number (1-indexed) is conatined in the file name
(``<N>`` in the above list entry). But this number does not take into account
the temporal derivative EVs of regular evs, so you need to offset this number
by the number of TD EVs in the design matrix, which come before the voxelwise
EV.
For voxelwise confound EVs, the column number mappings (1-indexed) are
contained in ``vef.dat`` and ``ven.dat``.
*Original* EV: An EV defined by the user
*Real* EV: A derived EV (either temporal derivative, or derived with
basis functions).
"""
import logging
import collections
import os.path as op
import numpy as np
from . import featanalysis
from . import image as fslimage
log = logging.getLogger(__name__)
class FSFError(Exception):
pass
class EV(object):
def __init__(self, index, title):
self.index = index
self.title = title
class NormalEV(EV):
def __init__(self, realIdx, origIdx, title):
EV.__init__(self, realIdx, title)
self.origIndex = origIdx
class TemporalDerivativeEV(NormalEV):
pass
class BasisFunctionEV(NormalEV):
pass
class VoxelwiseEV(NormalEV):
def __init__(self, realIdx, origIdx, title, filename):
NormalEV.__init__(self, realIdx, origIdx, title)
self.filename = filename
class ConfoundEV(EV):
def __init__(self, index, confIndex, title):
EV.__init__(self, index, title)
self.confIndex = confIndex
class MotionParameterEV(EV):
def __init__(self, index, motionIndex, title):
EV.__init__(self, index, title)
self.motionIndex = motionIndex
class VoxelwiseConfoundEV(EV):
def __init__(self, index, voxIndex, title, filename):
EV.__init__(self, index, title)
self.voxIndex = voxIndex
self.filename = filename
class FEATFSFDesign(object):
"""
"""
def __init__(self, featDir, settings, designMatrix):
# Get some information about the analysis
version = float(settings['version'])
level = int( settings['level'])
# Print a warning if we're
# using an old version of FEAT
if version < 6:
log.warning('Unsupported FEAT version: {}'.format(version))
# We need to parse the EVS a bit
# differently depending on whether
# this is a first level or higher
# level analysis.
if level == 1: getEVs = getFirstLevelEVs
else: getEVs = getHigherLevelEVs
self.__settings = collections.OrderedDict(settings.items())
self.__design = np.array(designMatrix)
self.__numEVs = self.__design.shape[1]
self.__evs = getEVs(featDir, self.__settings, self.__design)
for i, ev in enumerate(self.__evs):
print 'EV{}: {} [{}]'.format(
ev.index + 1,
ev.title,
type(ev).__name__)
if len(self.__evs) != self.__numEVs:
raise FSFError('Number of EVs does not match design.mat')
def getDesign(self, x, y, z):
"""Returns the design matrix for the specified voxel.
"""
# if no vox EVs, just
# return the design
pass
def getVoxelEVFile(self, idx):
return self.__evs[idx].filename
def getVoxelConfoundFile(self, idx):
return self.__evs[idx].filename
def getFirstLevelEVs(featDir, settings, designMat):
evs = []
origEVs = int(settings['evs_orig'])
# First, we loop through the EVs that
# are explicitly defined in design.fsf.
# This includes
# - normal EVs
# - temporal derivative EVs
# - basis function EVs
# - voxelwise EVs
for origIdx in range(origEVs):
title = settings[ 'evtitle{}' .format(origIdx + 1)]
shape = int(settings[ 'shape{}' .format(origIdx + 1)])
convolve = int(settings[ 'convolve{}' .format(origIdx + 1)])
deriv = int(settings[ 'deriv_yn{}' .format(origIdx + 1)])
basis = int(settings.get('basisfnum{}'.format(origIdx + 1), -1))
# Normal EV. This is just a column
# in the design matrix, defined by
# the user.
if shape != 9:
evs.append(NormalEV(len(evs), origIdx, title))
# Voxelwise EV. This is a 'normal' EV
# defined by the user, with different
# values for each voxel. The voxelwise
# values should be contained in the
# feat directory, in an image called
# designVoxelwiseEVN, where N is the
# original EV index.
else:
# The addExt function will
# raise an error if the
# file does not exist.
filename = op.join(
featDir, 'designVoxelwiseEV{}'.format(origIdx + 1))
filename = fslimage.addExt(filename, True)
evs.append(VoxelwiseEV(len(evs), origIdx, title, filename))
# This EV has been convolved with a set of basis
# functions. A set of N additional EVs have been
# added to the design matrix, immediately after
# the EV, where N is specified by the basisfnumN
# parameter in design.fsf.
if convolve in (4, 5, 6):
if basis == -1:
raise FSFError('Number of EVs is not specified '
'for basis function EV')
for i in range(basis - 1):
evs.append(BasisFunctionEV(len(evs), origIdx, title))
# A temporal derivative EV has been
# added for this EV - in the design
# matrix, it is the column immediately
# after this EV.
if deriv == 1:
evs.append(TemporalDerivativeEV(len(evs), origIdx, title))
# In the design matrix, after all EVs which
# have been explicilty defined, the rest of
# the EVs in the design matrix are confounds,
# in the following order:
#
# 1. Voxelwise confounds
# 2. Motion parameters
# 3. Other confounds
# Any voxelwise confounds are specified
# in two plain text files - vef.dat
# contains a comma-separated list of
# files, and ven.dat contains the column
# index of this confound in the design
# matrix (1-indexed). If these files
# don't exist, then it means that there
# are no voxelwise confounds.
#
# n.b. Even though the indices into the
# final design matrix are stored in ven.dat,
# I'm just assuming that the voxelwise
# confound columns are immediately after
# the 'real' EVs procesed above, in the
# order defined in vef.dat.
voxConfFiles = op.join(featDir, 'vef.dat')
voxConfLocs = op.join(featDir, 'ven.dat')
if op.exists(voxConfFiles) and op.exists(voxConfLocs):
with open(voxConfFiles, 'rt') as vcff:
voxConfFiles = vcff.read()
with open(voxConfLocs, 'rt') as vclf:
voxConfLocs = vclf.read()
voxConfFiles = voxConfFiles.strip().split(',')
voxConfLocs = voxConfLocs .strip().split(',')
if len(voxConfFiles) != len(voxConfLocs):
raise FSFError('vef.dat does not match ven.dat')
# An error will be raised if any of
# the files in vef.dat do not exist,
# or if any of the indices in
# ven.dat are not integers.
voxConfFiles = [op.join(featDir, f) for f in voxConfFiles]
voxConfFiles = [fslimage.addExt(f, True) for f in voxConfFiles]
voxConfLocs = [int(i) for i in voxConfLocs]
# Check to see if my assumption
# above, about the voxelwise
# confound EV locations, holds
startIdx = len(evs) + 1
if voxConfLocs != range(startIdx, startIdx + len(voxConfFiles)):
raise FSFError('Unsupported voxelwise confound ordering '
'({} -> {})'.format(startIdx, voxConfLocs))
# Create the voxelwise confound EVs
for i, (f, l) in enumerate(zip(voxConfFiles, voxConfLocs)):
evs.append(VoxelwiseConfoundEV(len(evs), i, 'voxconf', f))
# Have motion parameters been added
# as regressors to the design matrix?
motion = int(settings['motionevs'])
if motion == 0: numMotionEVs = 0
elif motion == 1: numMotionEVs = 6
elif motion == 2: numMotionEVs = 24
for i in range(numMotionEVs):
evs.append(MotionParameterEV(len(evs), i, 'motion'))
# Last step - any columns in the design
# matrix which have not yet been accounted
# for are other confounds, specified by
# the user with a text file.
numConfoundEVs = designMat.shape[1] - len(evs)
for i in range(numConfoundEVs):
evs.append(ConfoundEV(len(evs), i, 'confound'))
# Phew.
return evs
def getHigherLevelEVs(featDir, settings, designMat):
titleKeys = [s for s in settings.keys() if s.startswith('evtitle')]
evs = []
return evs