mvpa.misc.support

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 """Support function -- little helpers in everyday life""" 10 11 __docformat__ = 'restructuredtext' 12 13 import numpy as N 14 import re, os 15 16 from mvpa.base import warning 17 from mvpa.support.copy import copy, deepcopy 18 from operator import isSequenceType 19 20 if __debug__: 21 from mvpa.base import debug 22 23

24 -def reuseAbsolutePath(file1, file2, force=False):

25 """Use path to file1 as the path to file2 is no absolute 26 path is given for file2 27 28 :Parameters: 29 force : bool 30 if True, force it even if the file2 starts with / 31 """ 32 if not file2.startswith('/') or force: 33 # lets reuse path to file1 34 return os.path.join(os.path.dirname(file1), file2.lstrip('/')) 35 else: 36 return file2

37 38

39 -def transformWithBoxcar(data, startpoints, boxlength, offset=0, fx=N.mean):

40 """This function extracts boxcar windows from an array. Such a boxcar is 41 defined by a starting point and the size of the window along the first axis 42 of the array (`boxlength`). Afterwards a customizable function is applied 43 to each boxcar individually (Default: averaging). 44 45 :param data: An array with an arbitrary number of dimensions. 46 :type data: array 47 :param startpoints: Boxcar startpoints as index along the first array axis 48 :type startpoints: sequence 49 :param boxlength: Length of the boxcar window in #array elements 50 :type boxlength: int 51 :param offset: Optional offset between the configured starting point and the 52 actual begining of the boxcar window. 53 :type offset: int 54 :rtype: array (len(startpoints) x data.shape[1:]) 55 """ 56 if boxlength < 1: 57 raise ValueError, "Boxlength lower than 1 makes no sense." 58 59 # check for illegal boxes 60 for sp in startpoints: 61 if ( sp + offset + boxlength - 1 > len(data)-1 ) \ 62 or ( sp + offset < 0 ): 63 raise ValueError, \ 64 'Illegal box: start: %i, offset: %i, length: %i' \ 65 % (sp, offset, boxlength) 66 67 # build a list of list where each sublist contains the indexes of to be 68 # averaged data elements 69 selector = [ range( i + offset, i + offset + boxlength ) \ 70 for i in startpoints ] 71 72 # average each box 73 selected = [ fx( data[ N.array(box) ], axis=0 ) for box in selector ] 74 75 return N.array( selected )

76 77

78 -def _getUniqueLengthNCombinations_lt3(data, n):

79 """Generates a list of lists containing all combinations of 80 elements of data of length 'n' without repetitions. 81 82 data: list 83 n: integer 84 85 This function is adapted from a Java version posted in some forum on 86 the web as an answer to the question 'How can I generate all possible 87 combinations of length n?'. Unfortunately I cannot remember which 88 forum it was. 89 90 NOTE: implementation is broken for n>2 91 """ 92 93 if n > 2: 94 raise ValueError, "_getUniqueLengthNCombinations_lt3 " \ 95 "is broken for n>2, hence should not be used directly." 96 97 # to be returned 98 combos = [] 99 100 # local function that will be called recursively to collect the 101 # combination elements 102 def take(data, occupied, depth, taken): 103 for i, d in enumerate(data): 104 # only do something if this element hasn't been touch yet 105 if occupied[i] == False: 106 # see whether will reached the desired length 107 if depth < n-1: 108 # flag the current element as touched 109 occupied[i] = True 110 # next level 111 take(data, occupied, depth+1, taken + [d]) 112 # if the current element would be set 'free', it would 113 # results in ALL combinations of elements (obeying order 114 # of elements) and not just in the unique sets of 115 # combinations (without order) 116 #occupied[i] = False 117 else: 118 # store the final combination 119 combos.append(taken + [d])

120 # some kind of bitset that stores the status of each element 121 # (contained in combination or not) 122 occupied = [False] * len(data) 123 # get the combinations 124 take(data, occupied, 0, []) 125 126 # return the result 127 return combos 128

129 -def xuniqueCombinations(L, n):

130 """Generator of unique combinations form a list L of objects in 131 groups of size n. 132 133 # XXX EO: I guess they are already sorted. 134 # XXX EO: It seems to work well for n>20 :) 135 136 :Parameters: 137 L : list 138 list of unique ids 139 n : int 140 grouping size 141 142 Adopted from Li Daobing 143 http://code.activestate.com/recipes/190465/ 144 (MIT license, according to activestate.com's policy) 145 """ 146 if n==0: yield [] 147 else: 148 for i in xrange(len(L)-n+1): 149 for cc in xuniqueCombinations(L[i+1:],n-1): 150 yield [L[i]]+cc

151

152 -def _getUniqueLengthNCombinations_binary(L, n=None, sort=True):

153 """Find all subsets of data 154 155 :Parameters: 156 L : list 157 list of unique ids 158 n : None or int 159 If None, all possible subsets are returned. If n is specified (int), 160 then only the ones of the length n are returned 161 sort : bool 162 Either to sort the resultant sequence 163 164 Adopted from Alex Martelli: 165 http://mail.python.org/pipermail/python-list/2001-January/067815.html 166 """ 167 N = len(L) 168 if N > 20 and n == 1: 169 warning("getUniqueLengthNCombinations_binary should not be used for " 170 "large N") 171 result = [] 172 for X in range(2**N): 173 x = [ L[i] for i in range(N) if X & (1L<<i) ] 174 if n is None or len(x) == n: 175 # yield x # if we wanted to use it as a generator 176 result.append(x) 177 result.sort() 178 # if __debug__ and n is not None: 179 # # verify the result 180 # # would need scipy... screw it 181 # assert(len(result) == ...) 182 return result

183 184

185 -def getUniqueLengthNCombinations(L, n=None, sort=True):

186 """Find all subsets of data 187 188 :Parameters: 189 L : list 190 list of unique ids 191 n : None or int 192 If None, all possible subsets are returned. If n is specified (int), 193 then only the ones of the length n are returned 194 195 TODO: work out single stable implementation -- probably just by fixing 196 _getUniqueLengthNCombinations_lt3 197 """ 198 if n == 1: 199 return _getUniqueLengthNCombinations_lt3(L, n) 200 elif n==None: 201 return _getUniqueLengthNCombinations_binary(L, n, sort=True) 202 else: 203 # XXX EO: Is it safe to remove "list" and use generator 204 # directly to save memory for large number of 205 # combinations? It seems OK according to tests but 206 # I'd like a second opinion. 207 # XXX EO: what about inserting a warning if number of 208 # combinations > 10000? No one would run it... 209 return list(xuniqueCombinations(L, n))

210 211

212 -def _getUniqueLengthNCombinations_binary(L, n=None, sort=True):

213 """Find all subsets of data 214 215 :Parameters: 216 L : list 217 list of unique ids 218 n : None or int 219 If None, all possible subsets are returned. If n is specified (int), 220 then only the ones of the length n are returned 221 sort : bool 222 Either to sort the resultant sequence 223 224 Adopted from Alex Martelli: 225 http://mail.python.org/pipermail/python-list/2001-January/067815.html 226 """ 227 N = len(L) 228 if N > 20 and n == 1: 229 warning("getUniqueLengthNCombinations_binary should not be used for " 230 "large N") 231 result = [] 232 for X in range(2**N): 233 x = [ L[i] for i in range(N) if X & (1L<<i) ] 234 if n is None or len(x) == n: 235 # yield x # if we wanted to use it as a generator 236 result.append(x) 237 result.sort() 238 # if __debug__ and n is not None: 239 # # verify the result 240 # # would need scipy... screw it 241 # assert(len(result) == ...) 242 return result

243 244

245 -def getUniqueLengthNCombinations(L, n=None, sort=True):

246 """Find all subsets of data 247 248 :Parameters: 249 L : list 250 list of unique ids 251 n : None or int 252 If None, all possible subsets are returned. If n is specified (int), 253 then only the ones of the length n are returned 254 255 TODO: work out single stable implementation -- probably just by fixing 256 _getUniqueLengthNCombinations_lt3 257 """ 258 if n == 1: 259 return _getUniqueLengthNCombinations_lt3(L, n) 260 else: 261 return _getUniqueLengthNCombinations_binary(L, n, sort=True)

262 263

264 -def indentDoc(v):

265 """Given a `value` returns a string where each line is indented 266 267 Needed for a cleaner __repr__ output 268 `v` - arbitrary 269 """ 270 return re.sub('\n', '\n ', str(v))

271 272

273 -def idhash(val):

274 """Craft unique id+hash for an object 275 """ 276 res = "%s" % id(val) 277 if isinstance(val, list): 278 val = tuple(val) 279 try: 280 res += ":%s" % hash(buffer(val)) 281 except: 282 try: 283 res += ":%s" % hash(val) 284 except: 285 pass 286 pass 287 return res

288

289 -def isSorted(items):

290 """Check if listed items are in sorted order. 291 292 :Parameters: 293 `items`: iterable container 294 295 :return: `True` if were sorted. Otherwise `False` + Warning 296 """ 297 items_sorted = deepcopy(items) 298 items_sorted.sort() 299 equality = items_sorted == items 300 # XXX yarik forgotten analog to isiterable 301 if hasattr(equality, '__iter__'): 302 equality = N.all(equality) 303 return equality

304 305

306 -def isInVolume(coord, shape):

307 """For given coord check if it is within a specified volume size. 308 309 Returns True/False. Assumes that volume coordinates start at 0. 310 No more generalization (arbitrary minimal coord) is done to save 311 on performance 312 """ 313 for i in xrange(len(coord)): 314 if coord[i] < 0 or coord[i] >= shape[i]: 315 return False 316 return True

317 318

319 -def getBreakPoints(items, contiguous=True):

320 """Return a list of break points. 321 322 :Parameters: 323 items : iterable 324 list of items, such as chunks 325 contiguous : bool 326 if `True` (default) then raise Value Error if items are not 327 contiguous, i.e. a label occur in multiple contiguous sets 328 329 :raises: ValueError 330 331 :return: list of indexes for every new set of items 332 """ 333 prev = None # pylint happiness event! 334 known = [] 335 """List of items which was already seen""" 336 result = [] 337 """Resultant list""" 338 for index in xrange(len(items)): 339 item = items[index] 340 if item in known: 341 if index > 0: 342 if prev != item: # breakpoint 343 if contiguous: 344 raise ValueError, \ 345 "Item %s was already seen before" % str(item) 346 else: 347 result.append(index) 348 else: 349 known.append(item) 350 result.append(index) 351 prev = item 352 return result

353 354

355 -def RFEHistory2maps(history):

356 """Convert history generated by RFE into the array of binary maps 357 358 Example: 359 history2maps(N.array( [ 3,2,1,0 ] )) 360 results in 361 array([[ 1., 1., 1., 1.], 362 [ 1., 1., 1., 0.], 363 [ 1., 1., 0., 0.], 364 [ 1., 0., 0., 0.]]) 365 """ 366 367 # assure that it is an array 368 history = N.array(history) 369 nfeatures, steps = len(history), max(history) - min(history) + 1 370 history_maps = N.zeros((steps, nfeatures)) 371 372 for step in xrange(steps): 373 history_maps[step, history >= step] = 1 374 375 return history_maps

376 377

378 -class MapOverlap(object):

379 """Compute some overlap stats from a sequence of binary maps. 380 381 When called with a sequence of binary maps (e.g. lists or arrays) the 382 fraction of mask elements that are non-zero in a customizable proportion 383 of the maps is returned. By default this threshold is set to 1.0, i.e. 384 such an element has to be non-zero in *all* maps. 385 386 Three additional maps (same size as original) are computed: 387 388 * overlap_map: binary map which is non-zero for each overlapping element. 389 * spread_map: binary map which is non-zero for each element that is 390 non-zero in any map, but does not exceed the overlap 391 threshold. 392 * ovstats_map: map of float with the raw elementwise fraction of overlap. 393 394 All maps are available via class members. 395 """

396 - def __init__(self, overlap_threshold=1.0):

397 """Nothing to be seen here. 398 """ 399 self.__overlap_threshold = overlap_threshold 400 401 # pylint happiness block 402 self.overlap_map = None 403 self.spread_map = None 404 self.ovstats_map = None

405 406

407 - def __call__(self, maps):

408 """Returns fraction of overlapping elements. 409 """ 410 ovstats = N.mean(maps, axis=0) 411 412 self.overlap_map = (ovstats >= self.__overlap_threshold ) 413 self.spread_map = N.logical_and(ovstats > 0.0, 414 ovstats < self.__overlap_threshold) 415 self.ovstats_map = ovstats 416 417 return N.mean(ovstats >= self.__overlap_threshold)

418 419

420 -class Event(dict):

421 """Simple class to define properties of an event. 422 423 The class is basically a dictionary. Any properties can 424 be passed as keyword arguments to the constructor, e.g.: 425 426 >>> ev = Event(onset=12, duration=2.45) 427 428 Conventions for keys: 429 430 `onset` 431 The onset of the event in some unit. 432 `duration` 433 The duration of the event in the same unit as `onset`. 434 `label` 435 E.g. the condition this event is part of. 436 `chunk` 437 Group this event is part of (if any), e.g. experimental run. 438 `features` 439 Any amount of additional features of the event. This might include 440 things like physiological measures, stimulus intensity. Must be a mutable 441 sequence (e.g. list), if present. 442 """ 443 _MUSTHAVE = ['onset'] 444

445 - def __init__(self, **kwargs):

446 # store everything 447 dict.__init__(self, **kwargs) 448 449 # basic checks 450 for k in Event._MUSTHAVE: 451 if not self.has_key(k): 452 raise ValueError, "Event must have '%s' defined." % k

453 454

455 - def asDescreteTime(self, dt, storeoffset=False):

456 """Convert `onset` and `duration` information into descrete timepoints. 457 458 :Parameters: 459 dt: float 460 Temporal distance between two timepoints in the same unit as `onset` 461 and `duration`. 462 storeoffset: bool 463 If True, the temporal offset between original `onset` and 464 descretized `onset` is stored as an additional item in `features`. 465 466 :Return: 467 A copy of the original `Event` with `onset` and optionally `duration` 468 replaced by their corresponding descrete timepoint. The new onset will 469 correspond to the timepoint just before or exactly at the original 470 onset. The new duration will be the number of timepoints covering the 471 event from the computed onset timepoint till the timepoint exactly at 472 the end, or just after the event. 473 474 Note again, that the new values are expressed as #timepoint and not 475 in their original unit! 476 """ 477 dt = float(dt) 478 onset = self['onset'] 479 out = deepcopy(self) 480 481 # get the timepoint just prior the onset 482 out['onset'] = int(N.floor(onset / dt)) 483 484 if storeoffset: 485 # compute offset 486 offset = onset - (out['onset'] * dt) 487 488 if out.has_key('features'): 489 out['features'].append(offset) 490 else: 491 out['features'] = [offset] 492 493 if out.has_key('duration'): 494 # how many timepoint cover the event (from computed onset 495 # to the one timepoint just after the end of the event 496 out['duration'] = int(N.ceil((onset + out['duration']) / dt) \ 497 - out['onset']) 498 499 return out

500 501 502

503 -class HarvesterCall(object):

504 - def __init__(self, call, attribs=None, argfilter=None, expand_args=True, 505 copy_attribs=True):

506 """Initialize 507 508 :Parameters: 509 expand_args : bool 510 Either to expand the output of looper into a list of arguments for 511 call 512 attribs : list of basestr 513 What attributes of call to store and return later on? 514 copy_attribs : bool 515 Force copying values of attributes 516 """ 517 518 self.call = call 519 """Call which gets called in the harvester.""" 520 521 if attribs is None: 522 attribs = [] 523 if not isSequenceType(attribs): 524 raise ValueError, "'attribs' have to specified as a sequence." 525 526 if not (argfilter is None or isSequenceType(argfilter)): 527 raise ValueError, "'argfilter' have to be a sequence or None." 528 529 # now give it to me... 530 self.argfilter = argfilter 531 self.expand_args = expand_args 532 self.copy_attribs = copy_attribs 533 self.attribs = attribs

534 535 536

537 -class Harvester(object):

538 """World domination helper: do whatever it is asked and accumulate results 539 540 XXX Thinks about: 541 - Might we need to deepcopy attributes values? 542 - Might we need to specify what attribs to copy and which just to bind? 543 """ 544

545 - def __init__(self, source, calls, simplify_results=True):

546 """Initialize 547 548 :Parameters: 549 source 550 Generator which produce food for the calls. 551 calls : sequence of HarvesterCall instances 552 Calls which are processed in the loop. All calls are processed in 553 order of apperance in the sequence. 554 simplify_results: bool 555 Remove unecessary overhead in results if possible (nested lists 556 and dictionaries). 557 """ 558 if not isSequenceType(calls): 559 raise ValueError, "'calls' have to specified as a sequence." 560 561 self.__source = source 562 """Generator which feeds the harvester""" 563 564 self.__calls = calls 565 """Calls which gets called with each generated source""" 566 567 self.__simplify_results = simplify_results

568 569

570 - def __call__(self, *args, **kwargs):

571 """ 572 """ 573 # prepare complex result structure for all calls and their respective 574 # attributes: calls x dict(attributes x loop iterations) 575 results = [dict([('result', [])] + [(a, []) for a in c.attribs]) \ 576 for c in self.__calls] 577 578 # Lets do it! 579 for (i, X) in enumerate(self.__source(*args, **kwargs)): 580 for (c, call) in enumerate(self.__calls): 581 # sanity check 582 if i == 0 and call.expand_args and not isSequenceType(X): 583 raise RuntimeError, \ 584 "Cannot expand non-sequence result from %s" % \ 585 `self.__source` 586 587 # apply argument filter (and reorder) if requested 588 if call.argfilter: 589 filtered_args = [X[f] for f in call.argfilter] 590 else: 591 filtered_args = X 592 593 if call.expand_args: 594 result = call.call(*filtered_args) 595 else: 596 result = call.call(filtered_args) 597 598 # # XXX pylint doesn't like `` for some reason 599 # if __debug__: 600 # debug("LOOP", "Iteration %i on call %s. Got result %s" % 601 # (i, `self.__call`, `result`)) 602 603 604 results[c]['result'].append(result) 605 606 for attrib in call.attribs: 607 attrv = call.call.__getattribute__(attrib) 608 609 if call.copy_attribs: 610 attrv = copy(attrv) 611 612 results[c][attrib].append(attrv) 613 614 # reduce results structure 615 if self.__simplify_results: 616 # get rid of dictionary if just the results are requested 617 for (c, call) in enumerate(self.__calls): 618 if not len(call.attribs): 619 results[c] = results[c]['result'] 620 621 if len(self.__calls) == 1: 622 results = results[0] 623 624 return results

625 626 627 # XXX MH: this doesn't work in all cases, as you cannot have *args after a 628 # kwarg. 629 #def loop(looper, call, 630 # unroll=True, attribs=None, copy_attribs=True, *args, **kwargs): 631 # """XXX Loop twin brother 632 # 633 # Helper for those who just wants to do smth like 634 # loop(blah, bleh, grgr) 635 # instead of 636 # Loop(blah, bleh)(grgr) 637 # """ 638 # print looper, call, unroll, attribs, copy_attribs 639 # print args, kwargs 640 # return Loop(looper=looper, call=call, unroll=unroll, 641 # attribs=attribs, copy_attribs=copy_attribs)(*args, **kwargs) 642

Source Code for Module mvpa.misc.support