mvpa.misc.transformers

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 """Simply functors that transform something.""" 10 11 _DEV_DOC = """ 12 Follow the convetion that functions start with lower case, and classes 13 with uppercase letter. 14 """ 15 16 __docformat__ = 'restructuredtext' 17 18 19 import numpy as N 20 21 from mvpa.base import externals, warning 22 from mvpa.misc.state import StateVariable, ClassWithCollections 23 24 if __debug__: 25 from mvpa.base import debug 26 27

28 -def Absolute(x):

29 """ 30 Returns the elementwise absolute of any argument. 31 32 :Parameter: 33 x: scalar | sequence 34 35 """ 36 return N.absolute(x)

37 38

39 -def OneMinus(x):

40 """Returns elementwise '1 - x' of any argument.""" 41 return 1 - x

42 43

44 -def Identity(x):

45 """Return whatever it was called with.""" 46 return x

47 48

49 -def FirstAxisMean(x):

50 """Mean computed along the first axis.""" 51 return N.mean(x, axis=0)

52

53 -def FirstAxisSumNotZero(x):

54 """Sum computed over first axis of whether the values are not 55 equal to zero.""" 56 return (N.asarray(x)!=0).sum(axis=0)

57 58

59 -def SecondAxisMean(x):

60 """Mean across 2nd axis 61 62 Use cases: 63 - to combine multiple sensitivities to get sense about 64 mean sensitivity across splits 65 """ 66 return N.mean(x, axis=1)

67 68

69 -def SecondAxisSumOfAbs(x):

70 """Sum of absolute values along the 2nd axis 71 72 Use cases: 73 - to combine multiple sensitivities to get sense about 74 what features are most influential 75 """ 76 return N.abs(x).sum(axis=1)

77 78

79 -def SecondAxisMaxOfAbs(x):

80 """Max of absolute values along the 2nd axis 81 """ 82 return N.abs(x).max(axis=1)

83 84

85 -def GrandMean(x):

86 """Just what the name suggests.""" 87 return N.mean(x)

88 89

90 -def L2Normed(x, norm=1.0, reverse=False):

91 """Norm the values so that regular vector norm becomes `norm` 92 93 More verbose: Norm that the sum of the squared elements of the 94 returned vector becomes `norm`. 95 """ 96 xnorm = N.linalg.norm(x) 97 return x * (norm/xnorm)

98

99 -def L1Normed(x, norm=1.0, reverse=False):

100 """Norm the values so that L_1 norm (sum|x|) becomes `norm`""" 101 xnorm = N.sum(N.abs(x)) 102 return x * (norm/xnorm)

103 104

105 -def RankOrder(x, reverse=False):

106 """Rank-order by value. Highest gets 0""" 107 108 # XXX was Yarik on drugs? please simplify this beast 109 nelements = len(x) 110 indexes = N.arange(nelements) 111 t_indexes = indexes 112 if not reverse: 113 t_indexes = indexes[::-1] 114 tosort = zip(x, indexes) 115 tosort.sort() 116 ztosort = zip(tosort, t_indexes) 117 rankorder = N.empty(nelements, dtype=int) 118 rankorder[ [x[0][1] for x in ztosort] ] = \ 119 [x[1] for x in ztosort] 120 return rankorder

121 122

123 -def ReverseRankOrder(x):

124 """Convinience functor""" 125 return RankOrder(x, reverse=True)

126 127

128 -class OverAxis(object):

129 """Helper to apply transformer over specific axis 130 """ 131

132 - def __init__(self, transformer, axis=None):

133 """Initialize transformer wrapper with an axis. 134 135 :Parameters: 136 transformer 137 A callable to be used 138 axis : None or int 139 If None -- apply transformer across all the data. If some 140 int -- over that axis 141 """ 142 self.transformer = transformer 143 # sanity check 144 if not (axis is None or isinstance(axis, int)): 145 raise ValueError, "axis must be specified by integer" 146 self.axis = axis

147 148

149 - def __call__(self, x, *args, **kwargs):

150 transformer = self.transformer 151 axis = self.axis 152 if axis is None: 153 return transformer(x, *args, **kwargs) 154 155 x = N.asanyarray(x) 156 shape = x.shape 157 if axis >= len(shape): 158 raise ValueError, "Axis given in constructor %d is higher " \ 159 "than dimensionality of the data of shape %s" % (axis, shape) 160 161 # WRONG! ;-) 162 #for ind in xrange(shape[axis]): 163 # results.append(transformer(x.take([ind], axis=axis), 164 # *args, **kwargs)) 165 166 # TODO: more elegant/speedy solution? 167 shape_sweep = shape[:axis] + shape[axis+1:] 168 shrinker = None 169 """Either transformer reduces the dimensionality of the data""" 170 #results = N.empty(shape_out, dtype=x.dtype) 171 for index_sweep in N.ndindex(shape_sweep): 172 if axis > 0: 173 index = index_sweep[:axis] 174 else: 175 index = () 176 index = index + (slice(None),) + index_sweep[axis:] 177 x_sel = x[index] 178 x_t = transformer(x_sel, *args, **kwargs) 179 if shrinker is None: 180 if N.isscalar(x_t) or x_t.shape == shape_sweep: 181 results = N.empty(shape_sweep, dtype=x.dtype) 182 shrinker = True 183 elif x_t.shape == x_sel.shape: 184 results = N.empty(x.shape, dtype=x.dtype) 185 shrinker = False 186 else: 187 raise RuntimeError, 'Not handled by OverAxis kind of transformer' 188 189 if shrinker: 190 results[index_sweep] = x_t 191 else: 192 results[index] = x_t 193 194 return results

195 196

197 -class DistPValue(ClassWithCollections):

198 """Converts values into p-values under vague and non-scientific assumptions 199 """ 200 201 nulldist_number = StateVariable(enabled=True, 202 doc="Number of features within the estimated null-distribution") 203 204 positives_recovered = StateVariable(enabled=True, 205 doc="Number of features considered to be positives and which were recovered") 206 207

208 - def __init__(self, sd=0, distribution='rdist', fpp=None, nbins=400, **kwargs):

209 """L2-Norm the values, convert them to p-values of a given distribution. 210 211 :Parameters: 212 sd : int 213 Samples dimension (if len(x.shape)>1) on which to operate 214 distribution : string 215 Which distribution to use. Known are: 'rdist' (later normal should 216 be there as well) 217 fpp : float 218 At what p-value (both tails) if not None, to control for false 219 positives. It would iteratively prune the tails (tentative real positives) 220 until empirical p-value becomes less or equal to numerical. 221 nbins : int 222 Number of bins for the iterative pruning of positives 223 224 WARNING: Highly experimental/slow/etc: no theoretical grounds have been 225 presented in any paper, nor proven 226 """ 227 externals.exists('scipy', raiseException=True) 228 ClassWithCollections.__init__(self, **kwargs) 229 230 self.sd = sd 231 if not (distribution in ['rdist']): 232 raise ValueError, "Actually only rdist supported at the moment" \ 233 " got %s" % distribution 234 self.distribution = distribution 235 self.fpp = fpp 236 self.nbins = nbins

237 238

239 - def __call__(self, x):

240 from mvpa.support.stats import scipy 241 import scipy.stats as stats 242 243 # some local bindings 244 distribution = self.distribution 245 sd = self.sd 246 fpp = self.fpp 247 nbins = self.nbins 248 249 x = N.asanyarray(x) 250 shape_orig = x.shape 251 ndims = len(shape_orig) 252 253 # (very) old numpy had different format of returned bins -- 254 # there were not edges but center points. We care here about 255 # center points, so we will transform edge points into center 256 # points for newer versions of numpy 257 numpy_center_points = externals.versions['numpy'] < (1, 1) 258 259 # XXX May be just utilize OverAxis transformer? 260 if ndims > 2: 261 raise NotImplementedError, \ 262 "TODO: add support for more than 2 dimensions" 263 elif ndims == 1: 264 x, sd = x[:, N.newaxis], 0 265 266 # lets transpose for convenience 267 if sd == 0: x = x.T 268 269 # Output p-values of x in null-distribution 270 pvalues = N.zeros(x.shape) 271 nulldist_number, positives_recovered = [], [] 272 273 # finally go through all data 274 nd = x.shape[1] 275 if __debug__: 276 if nd < x.shape[0]: 277 warning("Number of features in DistPValue lower than number of" 278 " items -- may be incorrect sd=%d was provided" % sd) 279 for i, xx in enumerate(x): 280 dist = stats.rdist(nd-1, 0, 1) 281 xx /= N.linalg.norm(xx) 282 283 if fpp is not None: 284 if __debug__: 285 debug('TRAN_', "starting adaptive adjustment i=%d" % i) 286 287 # Adaptive adjustment for false negatives: 288 Nxx, xxx, pN_emp_prev = len(xx), xx, 1.0 289 Nxxx = Nxx 290 indexes = N.arange(Nxx) 291 """What features belong to Null-distribution""" 292 while True: 293 hist, bins = N.histogram(xxx, bins=nbins, normed=False) 294 pdf = hist.astype(float)/Nxxx 295 if not numpy_center_points: 296 # if we obtain edge points for bins -- take centers 297 bins = 0.5 * (bins[0:-1] + bins[1:]) 298 bins_halfstep = (bins[1] - bins[2])/2.0 299 300 # theoretical CDF 301 # was really unstable -- now got better ;) 302 dist_cdf = dist.cdf(bins) 303 304 # otherwise just recompute manually 305 # dist_pdf = dist.pdf(bins) 306 # dist_pdf /= N.sum(dist_pdf) 307 308 # XXX can't recall the function... silly 309 # probably could use N.integrate 310 cdf = N.zeros(nbins, dtype=float) 311 #dist_cdf = cdf.copy() 312 dist_prevv = cdf_prevv = 0.0 313 for j in range(nbins): 314 cdf_prevv = cdf[j] = cdf_prevv + pdf[j] 315 #dist_prevv = dist_cdf[j] = dist_prevv + dist_pdf[j] 316 317 # what bins fall into theoretical 'positives' in both tails 318 p = (0.5 - N.abs(dist_cdf - 0.5) < fpp/2.0) 319 320 # amount in empirical tails -- if we match theoretical, we 321 # should have total >= p 322 323 pN_emp = N.sum(pdf[p]) # / (1.0 * nbins) 324 325 if __debug__: 326 debug('TRAN_', "empirical p=%.3f for theoretical p=%.3f" 327 % (pN_emp, fpp)) 328 329 if pN_emp <= fpp: 330 # we are done 331 break 332 333 if pN_emp > pN_emp_prev: 334 if __debug__: 335 debug('TRAN_', "Diverging -- thus keeping last result " 336 "with p=%.3f" % pN_emp_prev) 337 # we better restore previous result 338 indexes, xxx, dist = indexes_prev, xxx_prev, dist_prev 339 break 340 341 pN_emp_prev = pN_emp 342 # very silly way for now -- just proceed by 1 bin 343 keep = N.logical_and(xxx > bins[0], # + bins_halfstep, 344 xxx < bins[-1]) # - bins_halfstep) 345 if __debug__: 346 debug('TRAN_', "Keeping %d out of %d elements" % 347 (N.sum(keep), Nxxx)) 348 349 # Preserve them if we need to "roll back" 350 indexes_prev, xxx_prev, dist_prev = indexes, xxx, dist 351 352 # we should remove those which we assume to be positives and 353 # which should not belong to Null-dist 354 xxx, indexes = xxx[keep], indexes[keep] 355 # L2 renorm it 356 xxx = xxx / N.linalg.norm(xxx) 357 Nxxx = len(xxx) 358 dist = stats.rdist(Nxxx-1, 0, 1) 359 360 Nindexes = len(indexes) 361 Nrecovered = Nxx - Nindexes 362 363 nulldist_number += [Nindexes] 364 positives_recovered += [Nrecovered] 365 366 if __debug__: 367 if distribution == 'rdist': 368 assert(dist.args[0] == Nindexes-1) 369 debug('TRAN', "Positives recovery finished with %d out of %d " 370 "entries in Null-distribution, thus %d positives " 371 "were recovered" % (Nindexes, Nxx, Nrecovered)) 372 373 # And now we need to perform our duty -- assign p-values 374 #dist = stats.rdist(Nindexes-1, 0, 1) 375 pvalues[i, :] = dist.cdf(xx) 376 377 # XXX we might add an option to transform it to z-scores? 378 result = pvalues 379 380 # charge state variables 381 # XXX might want to populate them for non-adaptive handling as well 382 self.nulldist_number = nulldist_number 383 self.positives_recovered = positives_recovered 384 385 # transpose if needed 386 if sd == 0: 387 result = result.T 388 389 return result

390

Source Code for Module mvpa.misc.transformers