pylal.metaarray

1 """ 2 MetaArray - a subclass of ndarray that holds metadata and preserves it across 3 array operations. 4 Metadata - a class for metadata stored in MetaArray 5 MetaArrayList - a subclass of list that ollows element-wise MetaArray 6 operations 7 8 Spectrum - a subclass of MetaArray as demonstration 9 SpectrumMetadata - a subclass of MetaData that strictly checks for metadata 10 compatibility between two Spectra 11 SpectrumList - subclass of MetaArrayList that has some nice features specific 12 to SpectrumMetadata 13 """ 14 from __future__ import division 15 16 __author__ = "Nickolas Fotopoulos <nvf@gravity.phys.uwm.edu>" 17 18 import operator 19 import os 20 import sys 21 import copy 22 23 import numpy 24 from glue import segments 25 26 ############################################################################## 27 # Method wrappers 28 ############################################################################## 29

30 -class _arraymethod(object):

31 """ 32 Used to attach methods to MetaArrays. Wrap ndarray methods that return 33 a single array. Merge metadata of all input Spectra. 34 35 __init__ gets called when we define the MetaArray (sub-)class and attach 36 methods. 37 __get__ gets called on the object at method call time. 38 __call__ is called immediately after __get__. 39 """

40 - def __init__(self, methname):

41 self._name = methname 42 self.obj = None 43 self.__doc__ = getattr(numpy.ndarray, self._name, None).__doc__

44

45 - def __get__(self, obj, objtype=None):

46 self.obj = obj 47 return self

48

49 - def __call__(self, *args, **params):

50 data = self.obj.A 51 cls = type(self.obj) 52 result = getattr(data, self._name)(*args, **params).view(cls) 53 result.metadata = self.obj.metadata.copy() 54 for arg in args: 55 result.metadata |= getattr(arg, 'metadata', None) 56 return result

57

58 -class _elementwise_method(object):

59 """ 60 Used to attach methods to MetaArrayLists. Wrap ndarray methods that 61 operate upon arrays and apply them to each element of the list. 62 63 __init__ gets called when we define the MetaArrayList class and attach 64 methods. 65 __get__ gets called on the list object at method call time. 66 __call__ is called immediately after __get__. 67 """

68 - def __init__(self, methname):

69 self._name = methname 70 self.list = None 71 self.listtype = None 72 self.itemtype = None

73 # I don't think I can get the docstring, unfortunately 74 # self.__doc__ = getattr(itemtype, self._name, None).__doc__ 75

76 - def __get__(self, list, listtype=None):

77 self.list = list 78 self.listtype = type(list) 79 self.itemtype = list._itemtype 80 return self

81

82 - def __call__(self, *others, **params):

83 result = self.listtype([]) 84 cls = self.itemtype 85 # if we have two alike lists, just take pairwise action 86 if len(others) == 1 and isinstance(others[0], self.listtype): 87 result.extend([getattr(a, self._name)(b, **params).view(cls) \ 88 for a, b in zip(self.list, others[0])]) 89 # otherwise, just try passing the args along 90 else: 91 result.extend([getattr(a, self._name)(*others, **params).view(cls) for a in self.list]) 92 return result

93 94 ############################################################################## 95 # Generic classes 96 ############################################################################## 97

98 -class Metadata(object):

99 """ 100 Abstract class to hold metadata 101 """ 102 # specify all allowed metadata here; type is required 103 typemap = {} 104 __slots__ = typemap.keys() + ['typemap'] 105

106 - def __new__(cls, metadata=None):

107 if isinstance(metadata, Metadata): 108 return metadata 109 return super(Metadata, cls).__new__(cls)

110

111 - def __init__(self, metadata):

112 # Recycle if possible 113 if isinstance(metadata, Metadata): 114 return 115 # user forgot to provide metadata 116 elif metadata is None: 117 return None 118 # create Metadata from dictionary 119 elif isinstance(metadata, dict): 120 for key, val in metadata.items(): 121 try: 122 setattr(self, key, self.typemap[key](val)) 123 except KeyError, e: 124 raise KeyError, \ 125 "This key is not in the typemap of %s: %s" \ 126 % (type(self), str(e)) 127 128 # all fields must be filled 129 for slot in self.__slots__: 130 if not hasattr(self, slot): 131 raise ValueError, \ 132 "Not enough metadata supplied; missing %s" % slot 133 else: 134 raise NotImplementedError

135

136 - def __str__(self):

137 return str(dict([(slot, getattr(self, slot)) for slot \ 138 in self.__slots__ if slot != "typemap"]))

139

140 - def __repr__(self):

141 return repr(dict([(slot, getattr(self, slot)) for slot \ 142 in self.__slots__ if slot != "typemap"]))

143

144 - def __ior__(self, other):

145 """ 146 Merge metadata; this must be subclassed. 147 """ 148 raise NotImplementedError

149

150 - def __or__(self, other):

151 return self.copy().__ior__(other)

152

153 - def __ror__(self, other):

154 return self.copy().__ior__(other)

155

156 - def __eq__(self, other):

157 for slot in self.__slots__: 158 if getattr(self, slot) != getattr(other, slot): 159 return False 160 return True

161

162 - def copy(self):

163 return type(self)(dict([(slot, getattr(self, slot)) for slot \ 164 in self.__slots__ if slot != "typemap"]))

165

166 - def todict(self):

167 return dict([(key, getattr(self, key)) for key in self.__slots__ \ 168 if key != "typemap"])

169

170 -class MetaArray(numpy.ndarray):

171 """ 172 An array containing a data and metadata. Intended to be 173 subclassed. 174 175 On b = MetaArray(a) where a is a MetaArray, metadata is copied, but data 176 is not. 177 """ 178 __array_priority__ = 10.1 # ufuncs mixing array types return MetaArray 179 _metadata_type = Metadata 180

181 - def __new__(subtype, data=None, metadata=None, dtype=None, copy=False, 182 subok=True):

183 # Case 1: data is an MetaArray, it has metadata already. 184 if isinstance(data, MetaArray): 185 if dtype is None: 186 dtype = data.dtype 187 else: 188 dtype = numpy.dtype(dtype) 189 if not copy and dtype==data.dtype and metadata is None: 190 return data 191 elif metadata is not None: 192 # share memory, but not metadata 193 new = numpy.array(data, dtype=dtype, copy=copy, subok=True) 194 new.metadata = subtype._metadata_type(metadata) 195 return new 196 else: 197 # copy, update metadata, then return 198 new = data.astype(dtype) # always copies 199 new.metadata = data.metadata 200 new._baseclass = data._baseclass 201 return new 202 203 # All other cases, create a new array and attach metadata 204 # Unless you specified otherwise, we'll reuse memory from existing 205 # arrays. 206 new = numpy.array(data, dtype=dtype, copy=copy, subok=True) 207 _baseclass = type(new) 208 new = new.view(subtype) 209 new.metadata = subtype._metadata_type(metadata) 210 new._baseclass = _baseclass 211 return new

212

213 - def __array_finalize__(self, obj):

214 """ 215 Called anytime a MetaArray is returned; make sure that metadata 216 is set to something. 217 """ 218 self.metadata = getattr(obj, "metadata", None) 219 self._baseclass = getattr(obj, "_baseclass", type(obj))

220

221 - def __array_wrap__(self, obj):

222 """ 223 Called anytime a ufunc operates on a MetaArray and another object. 224 The result of the ufunc is obj. The MetaArray operand is self. 225 """ 226 result = obj.view(type(self)) 227 result.metadata = self.metadata.copy() 228 return result

229

230 - def __repr__(self):

231 return "%s, %s)" % (repr(self.view(numpy.ndarray))[:-1], repr(self.metadata))

232

233 - def __str__(self):

234 return "%s %s" % (str(self.view(numpy.ndarray)), str(self.metadata))

235 236 # methods that return an array, wrapped to return a MetaArray 237 __abs__ = _arraymethod('__abs__') 238 __add__ = _arraymethod('__add__') 239 __and__ = _arraymethod('__and__') 240 __copy__ = _arraymethod('__copy__') 241 __deepcopy__ = _arraymethod('__deepcopy__') 242 __div__ = _arraymethod('__div__') 243 __divmod__ = _arraymethod('__divmod__') 244 __floordiv__ = _arraymethod('__floordiv__') 245 __hex__ = _arraymethod('__hex__') 246 __iadd__ = _arraymethod('__iadd__') 247 __iand__ = _arraymethod('__iand__') 248 __idiv__ = _arraymethod('__idiv__') 249 __ifloordiv__ = _arraymethod('__ifloordiv__') 250 __ilshift__ = _arraymethod('__ilshift__') 251 __imod__ = _arraymethod('__imod__') 252 __imul__ = _arraymethod('__imul__') 253 __invert__ = _arraymethod('__invert__') 254 __ior__ = _arraymethod('__ior__') 255 __ipow__ = _arraymethod('__ipow__') 256 __irshift__ = _arraymethod('__irshift__') 257 __isub__ = _arraymethod('__isub__') 258 __itruediv__ = _arraymethod('__itruediv__') 259 __ixor__ = _arraymethod('__ixor__') 260 __lshift__ = _arraymethod('__lshift__') 261 __mul__ = _arraymethod('__mul__') 262 __rmod__ = _arraymethod('__rmod__') 263 __rmul__ = _arraymethod('__rmul__') 264 __ror__ = _arraymethod('__ror__') 265 __rpow__ = _arraymethod('__rpow__') 266 __rrshift__ = _arraymethod('__rrshift__') 267 __rshift__ = _arraymethod('__rshift__') 268 __rsub__ = _arraymethod('__rsub__') 269 __rtruediv__ = _arraymethod('__rtruediv__') 270 __rxor__ = _arraymethod('__rxor__') 271 __sub__ = _arraymethod('__sub__') 272 __truediv__ = _arraymethod('__truediv__') 273 __xor__ = _arraymethod('__xor__') 274 astype = _arraymethod('astype') 275 byteswap = _arraymethod('byteswap') 276 choose = _arraymethod('choose') 277 clip = _arraymethod('clip') 278 compress = _arraymethod('compress') 279 conj = _arraymethod('conj') 280 conjugate = _arraymethod('conjugate') 281 copy = _arraymethod('copy') 282 cumprod = _arraymethod('cumprod') 283 cumsum = _arraymethod('cumsum') 284 diagonal = _arraymethod('diagonal') 285 fill = _arraymethod('fill') 286 flat = _arraymethod('flat') 287 flatten = _arraymethod('flatten') 288 repeat = _arraymethod('repeat') 289 squeeze = _arraymethod('squeeze') 290 transpose = _arraymethod('transpose') 291 292 T = property(fget = lambda self: self.transpose()) 293 H = property(fget = lambda self: self.T.conj()) # Hermitian transpose 294

295 - def _get_data(self):

296 return self.view(self._baseclass)

297 A = property(fget=_get_data) # get at the underlying Array data 298 299 # Pickling

300 - def __getstate__(self):

301 """ 302 Returns the internal state of the object, for pickling purposes. 303 """ 304 state = (1, 305 self.shape, 306 self.dtype, 307 self.flags.fnc, 308 self.A.tostring(), 309 self.metadata.todict(), 310 ) 311 return state

312

313 - def __setstate__(self, state):

314 """ 315 Restores the internal state of the masked array, for unpickling 316 purposes. `state` is typically the output of the ``__getstate__`` 317 output, and is a 5-tuple: 318 319 - class name 320 - a tuple giving the shape of the data 321 - a typecode for the data 322 - a binary string for the data 323 - a binary string for the mask. 324 """ 325 (ver, shp, typ, isf, raw, meta) = state 326 if ver != 1: raise NotImplementedError 327 numpy.ndarray.__setstate__(self, (shp, typ, isf, raw)) 328 self.metadata = self._metadata_type(meta)

329

330 - def __reduce__(self):

331 """ 332 Returns a 3-tuple for pickling a MetaArray 333 - reconstruction function 334 - tuple to pass reconstruction function 335 - state, which will be passed to __setstate__ 336 """ 337 return (_mareconstruct, 338 (self.__class__, self._baseclass, (0,), 'b', ), 339 self.__getstate__())

340

341 -def _mareconstruct(subtype, baseclass, baseshape, basetype):

342 """ 343 Internal function that builds a new MaskedArray from the information 344 stored in a pickle. 345 """ 346 return subtype.__new__(subtype, [], dtype=basetype)

347

348 -class MetaArrayList(list):

349 _itemtype = MetaArray 350 351 # these methods will act upon each element of this list 352 __abs__ = _elementwise_method('__abs__') 353 __add__ = _elementwise_method('__add__') 354 __and__ = _elementwise_method('__and__') 355 __copy__ = _elementwise_method('__copy__') 356 __deepcopy__ = _elementwise_method('__deepcopy__') 357 __div__ = _elementwise_method('__div__') 358 __divmod__ = _elementwise_method('__divmod__') 359 __floordiv__ = _elementwise_method('__floordiv__') 360 __hex__ = _elementwise_method('__hex__') 361 __iadd__ = _elementwise_method('__iadd__') 362 __iand__ = _elementwise_method('__iand__') 363 __idiv__ = _elementwise_method('__idiv__') 364 __ifloordiv__ = _elementwise_method('__ifloordiv__') 365 __ilshift__ = _elementwise_method('__ilshift__') 366 __imod__ = _elementwise_method('__imod__') 367 __imul__ = _elementwise_method('__imul__') 368 __invert__ = _elementwise_method('__invert__') 369 __ior__ = _elementwise_method('__ior__') 370 __ipow__ = _elementwise_method('__ipow__') 371 __irshift__ = _elementwise_method('__irshift__') 372 __isub__ = _elementwise_method('__isub__') 373 __itruediv__ = _elementwise_method('__itruediv__') 374 __ixor__ = _elementwise_method('__ixor__') 375 __lshift__ = _elementwise_method('__lshift__') 376 __mul__ = _elementwise_method('__mul__') 377 __rmod__ = _elementwise_method('__rmod__') 378 __rmul__ = _elementwise_method('__rmul__') 379 __ror__ = _elementwise_method('__ror__') 380 __rpow__ = _elementwise_method('__rpow__') 381 __rrshift__ = _elementwise_method('__rrshift__') 382 __rshift__ = _elementwise_method('__rshift__') 383 __rsub__ = _elementwise_method('__rsub__') 384 __rtruediv__ = _elementwise_method('__rtruediv__') 385 __rxor__ = _elementwise_method('__rxor__') 386 __sub__ = _elementwise_method('__sub__') 387 __truediv__ = _elementwise_method('__truediv__') 388 __xor__ = _elementwise_method('__xor__') 389 astype = _elementwise_method('astype') 390 byteswap = _elementwise_method('byteswap') 391 choose = _elementwise_method('choose') 392 clip = _elementwise_method('clip') 393 compress = _elementwise_method('compress') 394 conj = _elementwise_method('conj') 395 conjugate = _elementwise_method('conjugate') 396 copy = _elementwise_method('copy') 397 cumprod = _elementwise_method('cumprod') 398 cumsum = _elementwise_method('cumsum') 399 diagonal = _elementwise_method('diagonal') 400 fill = _elementwise_method('fill') 401 flat = _elementwise_method('flat') 402 flatten = _elementwise_method('flatten') 403 repeat = _elementwise_method('repeat') 404 squeeze = _elementwise_method('squeeze') 405 sum = _elementwise_method('sum') 406 transpose = _elementwise_method('transpose') 407 408 T = property(fget = lambda self: self.transpose()) # Transpose 409 H = property(fget = lambda self: self.T.conj()) # Hermitian transpose 410 411 # special case a few attribute accessors

412 - def _get_real(self):

413 return type(self)([x.real for x in self])

414 real = property(fget=_get_real)

415 - def _get_imag(self):

416 return type(self)([x.imag for x in self])

417 imag = property(fget=_get_imag) 418

419 - def _get_data(self):

420 """ 421 Return a regular list of arrays. 422 """ 423 return [x.view(x._baseclass) for x in self]

424 A = property(fget=_get_data)

425 426 427 ############################################################################## 428 # Define TimeSeries and associated structures 429 ############################################################################## 430

431 -class TimeSeriesMetadata(Metadata):

432 """ 433 Hold the metadata associated with a spectrum, including frequency 434 resolution, a segmentlist indicating what times were involved in taking 435 the spectrum, a channel name, etc. 436 """ 437 # specify all allowed metadata here; type is required 438 typemap = {"name": str, 439 "dt": float, 440 "segments": segments.segmentlist, 441 "comments": list} 442 __slots__ = typemap.keys() + ['typemap'] 443

444 - def __ior__(self, other):

445 """ 446 Merge metadata. No repeats. Throw error on incompatible spectra. 447 Let None act as identity. 448 """ 449 if other is None: 450 return self 451 452 # check that metadata are compatible for merging 453 assert numpy.alltrue([getattr(self, attr) == getattr(other, attr) \ 454 for attr in self.__slots__ \ 455 if attr not in ('segments', 'comments', 'typemap')]) 456 457 # add, but do not join, segments 458 self.segments.extend([seg for seg in other.segments \ 459 if seg not in self.segments]) 460 self.segments.sort() 461 462 # add only new comments 463 self.comments.extend([comment for comment in other.comments \ 464 if comment not in self.comments]) 465 return self

466

467 -class TimeSeries(MetaArray):

468 """ 469 This is a MetaArray, but with the metadata typemap specified. 470 """ 471 _metadata_type = TimeSeriesMetadata 472

473 - def ordinates(self):

474 """ 475 Return an single-precision ndarray containing the times at 476 which this TimeSeries is sampled. 477 """ 478 m = self.metadata 479 segs = m.segments 480 segs.coalesce() 481 482 # start with a uniformly spaced domain 483 result = numpy.arange(len(self), dtype=numpy.float32)*m.dt 484 485 # split domain by applying different offsets for each segment 486 offsets = [seg[0] for seg in segs] 487 seg_lengths = [int(abs(seg)/m.dt) for seg in segs] 488 assert sum(seg_lengths) == len(self) 489 boundaries = [0]+[sum(seg_lengths[:i]) for i in range(1, len(seg_lengths)+1)] 490 assert boundaries[-1] == len(self) 491 slices = [slice(a, b) for a,b in zip(boundaries[:-1], boundaries[1:])] 492 493 for sl, offset in zip(slices, offsets): 494 result[sl] += offset 495 return result

496

497 -class TimeSeriesList(MetaArrayList):

498 """ 499 A list of TimeSeries. 500 """ 501 _itemtype = TimeSeries 502

503 - def segments(self):

504 """ 505 Return the (uncoalesced) list of segments represented by the TimeSeries 506 in this TimeSeriesList. 507 """ 508 segs = segments.segmentlist() 509 for series in self: 510 segs.extend(series.metadata.segments) 511 return segs

512

513 - def merge_list(self):

514 """ 515 Concatenate the list into one single TimeSeries. 516 """ 517 meta = reduce(operator.or_, [ts.metadata for ts in self]) 518 return TimeSeries(numpy.concatenate(self.A), meta)

519 520 ############################################################################## 521 # Define Spectrum and associated structures 522 ############################################################################## 523

524 -class SpectrumMetadata(Metadata):

525 """ 526 Hold the metadata associated with a spectrum, including frequency 527 resolution, a segmentlist indicating what times were involved in taking 528 the spectrum, a channel name, etc. 529 """ 530 # specify all allowed metadata here; type is required 531 typemap = {"name": str, 532 "df": float, 533 "f_low": float, 534 "segments": segments.segmentlist, 535 "comments": list} 536 __slots__ = typemap.keys() + ['typemap'] 537

538 - def __ior__(self, other):

539 """ 540 Merge metadata. No repeats. Throw error on incompatible spectra. 541 Let None act as identity. 542 """ 543 if other is None: 544 return self 545 if self is None: 546 return other 547 548 # check that metadata are compatible for merging 549 assert numpy.alltrue([getattr(self, attr) == getattr(other, attr) \ 550 for attr in ("df", "f_low")]) 551 552 # add, but do not join segments 553 self.segments.extend([seg for seg in other.segments \ 554 if seg not in self.segments]) 555 self.segments.sort() 556 557 # add only new comments 558 self.comments.extend([comment for comment in other.comments \ 559 if comment not in self.comments]) 560 return self

561

562 -class Spectrum(MetaArray):

563 """ 564 This is a MetaArray, but with the metadata typemap specified. 565 """ 566 _metadata_type = SpectrumMetadata 567

568 - def ordinates(self):

569 """ 570 Return an single-precision ndarray containing the frequencies at 571 which this Spectrum is sampled. 572 """ 573 m = self.metadata 574 return numpy.arange(len(self), dtype=numpy.float32)*m.df + m.f_low

575

576 -class SpectrumList(MetaArrayList):

577 """ 578 A list of Spectra. 579 """ 580 _itemtype = Spectrum 581

582 - def segments(self):

583 """ 584 Return the (uncoalesced) list of segments represented by the Spectra 585 in this SpectrumList. 586 """ 587 segs = segments.segmentlist() 588 for spectrum in self: 589 segs.extend(spectrum.metadata.segments) 590 return segs

591

592 - def ordinates(self):

593 """ 594 Return an single-precision ndarray containing the frequencies at 595 which these Spectrums are sampled (they are required to match). 596 For an empty list, return a zero-length array. 597 """ 598 if len(self) == 0: 599 return numpy.array([], dtype=numpy.float32) 600 m = self[0].metadata 601 return numpy.arange(len(self[0]), dtype=numpy.float32)*m.df + m.f_low

602

603 - def sum_spectra(self):

604 """ 605 Sum across Spectrums, returning a single Spectrum 606 """ 607 if len(self) == 0: 608 raise ValueError 609 meta = reduce(operator.or_, [s.metadata for s in self]) 610 return Spectrum(numpy.sum(self.A, axis=0), meta)

611

612 -class SpectrumDict(dict):

613 """ 614 A dictionary allowing access to FFTs of different data streams and 615 providing convenient batch operations. 616 """ 617 pass

618

Source Code for Module pylal.metaarray