mvpa.datasets.miscfx

1 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 2 # vi: set ft=python sts=4 ts=4 sw=4 et: 3 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 4 # 5 # See COPYING file distributed along with the PyMVPA package for the 6 # copyright and license terms. 7 # 8 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 9 """Misc function performing operations on datasets which are based on scipy 10 """ 11 12 __docformat__ = 'restructuredtext' 13 14 from mvpa.base import externals 15 16 import numpy as N 17 np = N # just for easy merging of changes into >=0.5 18 19 from operator import isSequenceType 20 21 from mvpa.datasets.base import Dataset, datasetmethod 22 from mvpa.misc.support import getBreakPoints 23 24 if externals.exists('scipy', raiseException=True): 25 from scipy import signal 26 from scipy.linalg import lstsq 27 from scipy.special import legendre

28 29 - def legendre_(n, x):

30 # Scipy 0.8.0 (and possibly later) has regression of reporting 31 # 'inf's for negative boundary. Lets guard against it for now 32 leg = legendre(n) 33 r = leg(x) 34 infs = np.isinf(r) 35 if np.any(infs): 36 r[infs] = leg(x[infs] + 1e-10) # offset to try to overcome problems 37 return r

38

39 @datasetmethod 40 -def detrend(dataset, perchunk=False, model='linear', 41 polyord=None, opt_reg=None):

42 """ 43 Given a dataset, detrend the data inplace either entirely 44 or per each chunk 45 46 :Parameters: 47 dataset : Dataset 48 dataset to operate on 49 perchunk : bool 50 either to operate on whole dataset at once or on each chunk 51 separately 52 model 53 Type of detrending model to run. If 'linear' or 'constant', 54 scipy.signal.detrend is used to perform a linear or demeaning 55 detrend. Polynomial detrending is activated when 'regress' is 56 used or when polyord or opt_reg are specified. 57 polyord : int or list 58 Order of the Legendre polynomial to remove from the data. This 59 will remove every polynomial up to and including the provided 60 value. For example, 3 will remove 0th, 1st, 2nd, and 3rd order 61 polynomials from the data. N.B.: The 0th polynomial is the 62 baseline shift, the 1st is the linear trend. 63 If you specify a single int and perchunk is True, then this value 64 is used for each chunk. You can also specify a different polyord 65 value for each chunk by providing a list or ndarray of polyord 66 values the length of the number of chunks. 67 opt_reg : ndarray 68 Optional ndarray of additional information to regress out from the 69 dataset. One example would be to regress out motion parameters. 70 As with the data, time is on the first axis. 71 72 """ 73 if polyord is not None or opt_reg is not None: 74 model='regress' 75 76 if model in ['linear', 'constant']: 77 # perform scipy detrend 78 bp = 0 # no break points by default 79 80 if perchunk: 81 try: 82 bp = getBreakPoints(dataset.chunks) 83 except ValueError, e: 84 raise ValueError, \ 85 "Failed to assess break points between chunks. Often " \ 86 "that is due to discontinuities within a chunk, which " \ 87 "ruins idea of per-chunk detrending. Original " \ 88 "exception was: %s" % str(e) 89 90 dataset.samples[:] = signal.detrend(dataset.samples, axis=0, 91 type=model, bp=bp) 92 elif model in ['regress']: 93 # perform regression-based detrend 94 return __detrend_regress(dataset, perchunk=perchunk, 95 polyord=polyord, opt_reg=opt_reg) 96 else: 97 # raise exception because not found 98 raise ValueError('Specified model type (%s) is unknown.' 99 % (model))

100

101 102 103 -def __detrend_regress(dataset, perchunk=True, polyord=None, opt_reg=None):

104 """ 105 Given a dataset, perform a detrend inplace, regressing out polynomial 106 terms as well as optional regressors, such as motion parameters. 107 108 :Parameters: 109 dataset : Dataset 110 Dataset to operate on 111 perchunk : bool 112 Either to operate on whole dataset at once or on each chunk 113 separately. If perchunk is True, all the samples within a 114 chunk should be contiguous and the chunks should be sorted in 115 order from low to high. 116 polyord : int 117 Order of the Legendre polynomial to remove from the data. This 118 will remove every polynomial up to and including the provided 119 value. For example, 3 will remove 0th, 1st, 2nd, and 3rd order 120 polynomials from the data. N.B.: The 0th polynomial is the 121 baseline shift, the 1st is the linear trend. 122 If you specify a single int and perchunk is True, then this value 123 is used for each chunk. You can also specify a different polyord 124 value for each chunk by providing a list or ndarray of polyord 125 values the length of the number of chunks. 126 opt_reg : ndarray 127 Optional ndarray of additional information to regress out from the 128 dataset. One example would be to regress out motion parameters. 129 As with the data, time is on the first axis. 130 """ 131 132 # Process the polyord to be a list with length of the number of chunks 133 if not polyord is None: 134 if not isSequenceType(polyord): 135 # repeat to be proper length 136 polyord = [polyord]*len(dataset.uniquechunks) 137 elif perchunk and len(polyord) != len(dataset.uniquechunks): 138 raise ValueError("If you specify a sequence of polyord values " + \ 139 "they sequence length must match the " + \ 140 "number of unique chunks in the dataset.") 141 142 # loop over chunks if necessary 143 if perchunk: 144 # get the unique chunks 145 uchunks = dataset.uniquechunks 146 147 # loop over each chunk 148 reg = [] 149 for n, chunk in enumerate(uchunks): 150 # get the indices for that chunk 151 cinds = dataset.chunks == chunk 152 153 # see if add in polyord values 154 if not polyord is None: 155 # create the timespan 156 x = N.linspace(-1, 1, cinds.sum()) 157 # create each polyord with the value for that chunk 158 for n in range(polyord[n] + 1): 159 newreg = N.zeros((dataset.nsamples, 1)) 160 newreg[cinds, 0] = legendre_(n, x) 161 reg.append(newreg) 162 else: 163 # take out mean over entire dataset 164 reg = [] 165 # see if add in polyord values 166 if not polyord is None: 167 # create the timespan 168 x = N.linspace(-1, 1, dataset.nsamples) 169 for n in range(polyord[0] + 1): 170 reg.append(legendre_(n, x)[:, N.newaxis]) 171 172 # see if add in optional regs 173 if not opt_reg is None: 174 # add in the optional regressors, too 175 reg.append(opt_reg) 176 177 # combine the regs 178 if len(reg) > 0: 179 if len(reg) > 1: 180 regs = N.hstack(reg) 181 else: 182 regs = reg[0] 183 else: 184 # no regs to remove 185 raise ValueError("You must specify at least one " + \ 186 "regressor to regress out.") 187 188 # perform the regression 189 res = lstsq(regs, dataset.samples) 190 191 # remove all but the residuals 192 yhat = N.dot(regs, res[0]) 193 dataset.samples -= yhat 194 195 # return the results 196 return res, regs

197

Source Code for Module mvpa.datasets.miscfx_sp