irregularlysampledsignal.py 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494
  1. # -*- coding: utf-8 -*-
  2. '''
  3. This module implements :class:`IrregularlySampledSignal`, an array of analog
  4. signals with samples taken at arbitrary time points.
  5. :class:`IrregularlySampledSignal` inherits from :class:`basesignal.BaseSignal`
  6. and derives from :class:`BaseNeo`, from :module:`neo.core.baseneo`,
  7. and from :class:`quantities.Quantity`, which inherits from :class:`numpy.array`.
  8. Inheritance from :class:`numpy.array` is explained here:
  9. http://docs.scipy.org/doc/numpy/user/basics.subclassing.html
  10. In brief:
  11. * Initialization of a new object from constructor happens in :meth:`__new__`.
  12. This is where user-specified attributes are set.
  13. * :meth:`__array_finalize__` is called for all new objects, including those
  14. created by slicing. This is where attributes are copied over from
  15. the old object.
  16. '''
  17. # needed for Python 3 compatibility
  18. from __future__ import absolute_import, division, print_function
  19. import numpy as np
  20. import quantities as pq
  21. from neo.core.baseneo import BaseNeo, MergeError, merge_annotations
  22. from neo.core import basesignal
  23. def _new_IrregularlySampledSignal(cls, times, signal, units=None, time_units=None, dtype=None,
  24. copy=True, name=None, file_origin=None, description=None,
  25. annotations=None, segment=None, channel_index=None):
  26. '''
  27. A function to map IrregularlySampledSignal.__new__ to function that
  28. does not do the unit checking. This is needed for pickle to work.
  29. '''
  30. iss = cls(times=times, signal=signal, units=units, time_units=time_units,
  31. dtype=dtype, copy=copy, name=name, file_origin=file_origin,
  32. description=description, **annotations)
  33. iss.segment = segment
  34. iss.channel_index = channel_index
  35. return iss
  36. class IrregularlySampledSignal(basesignal.BaseSignal):
  37. '''
  38. An array of one or more analog signals with samples taken at arbitrary time points.
  39. A representation of one or more continuous, analog signals acquired at time
  40. :attr:`t_start` with a varying sampling interval. Each channel is sampled
  41. at the same time points.
  42. *Usage*::
  43. >>> from neo.core import IrregularlySampledSignal
  44. >>> from quantities import s, nA
  45. >>>
  46. >>> irsig0 = IrregularlySampledSignal([0.0, 1.23, 6.78], [1, 2, 3],
  47. ... units='mV', time_units='ms')
  48. >>> irsig1 = IrregularlySampledSignal([0.01, 0.03, 0.12]*s,
  49. ... [[4, 5], [5, 4], [6, 3]]*nA)
  50. *Required attributes/properties*:
  51. :times: (quantity array 1D, numpy array 1D, or list)
  52. The time of each data point. Must have the same size as :attr:`signal`.
  53. :signal: (quantity array 2D, numpy array 2D, or list (data, channel))
  54. The data itself.
  55. :units: (quantity units)
  56. Required if the signal is a list or NumPy array, not if it is
  57. a :class:`Quantity`.
  58. :time_units: (quantity units) Required if :attr:`times` is a list or
  59. NumPy array, not if it is a :class:`Quantity`.
  60. *Recommended attributes/properties*:.
  61. :name: (str) A label for the dataset
  62. :description: (str) Text description.
  63. :file_origin: (str) Filesystem path or URL of the original data file.
  64. *Optional attributes/properties*:
  65. :dtype: (numpy dtype or str) Override the dtype of the signal array.
  66. (times are always floats).
  67. :copy: (bool) True by default.
  68. Note: Any other additional arguments are assumed to be user-specific
  69. metadata and stored in :attr:`annotations`.
  70. *Properties available on this object*:
  71. :sampling_intervals: (quantity array 1D) Interval between each adjacent
  72. pair of samples.
  73. (``times[1:] - times[:-1]``)
  74. :duration: (quantity scalar) Signal duration, read-only.
  75. (``times[-1] - times[0]``)
  76. :t_start: (quantity scalar) Time when signal begins, read-only.
  77. (``times[0]``)
  78. :t_stop: (quantity scalar) Time when signal ends, read-only.
  79. (``times[-1]``)
  80. *Slicing*:
  81. :class:`IrregularlySampledSignal` objects can be sliced. When this
  82. occurs, a new :class:`IrregularlySampledSignal` (actually a view) is
  83. returned, with the same metadata, except that :attr:`times` is also
  84. sliced in the same way.
  85. *Operations available on this object*:
  86. == != + * /
  87. '''
  88. _single_parent_objects = ('Segment', 'ChannelIndex')
  89. _quantity_attr = 'signal'
  90. _necessary_attrs = (('times', pq.Quantity, 1),
  91. ('signal', pq.Quantity, 2))
  92. def __new__(cls, times, signal, units=None, time_units=None, dtype=None,
  93. copy=True, name=None, file_origin=None,
  94. description=None,
  95. **annotations):
  96. '''
  97. Construct a new :class:`IrregularlySampledSignal` instance.
  98. This is called whenever a new :class:`IrregularlySampledSignal` is
  99. created from the constructor, but not when slicing.
  100. '''
  101. if units is None:
  102. if hasattr(signal, "units"):
  103. units = signal.units
  104. else:
  105. raise ValueError("Units must be specified")
  106. elif isinstance(signal, pq.Quantity):
  107. # could improve this test, what if units is a string?
  108. if units != signal.units:
  109. signal = signal.rescale(units)
  110. if time_units is None:
  111. if hasattr(times, "units"):
  112. time_units = times.units
  113. else:
  114. raise ValueError("Time units must be specified")
  115. elif isinstance(times, pq.Quantity):
  116. # could improve this test, what if units is a string?
  117. if time_units != times.units:
  118. times = times.rescale(time_units)
  119. # should check time units have correct dimensions
  120. obj = pq.Quantity.__new__(cls, signal, units=units,
  121. dtype=dtype, copy=copy)
  122. if obj.ndim == 1:
  123. obj = obj.reshape(-1, 1)
  124. if len(times) != obj.shape[0]:
  125. raise ValueError("times array and signal array must "
  126. "have same length")
  127. obj.times = pq.Quantity(times, units=time_units,
  128. dtype=float, copy=copy)
  129. obj.segment = None
  130. obj.channel_index = None
  131. return obj
  132. def __init__(self, times, signal, units=None, time_units=None, dtype=None,
  133. copy=True, name=None, file_origin=None, description=None,
  134. **annotations):
  135. '''
  136. Initializes a newly constructed :class:`IrregularlySampledSignal`
  137. instance.
  138. '''
  139. BaseNeo.__init__(self, name=name, file_origin=file_origin,
  140. description=description, **annotations)
  141. def __reduce__(self):
  142. '''
  143. Map the __new__ function onto _new_IrregularlySampledSignal, so that pickle
  144. works
  145. '''
  146. return _new_IrregularlySampledSignal, (self.__class__,
  147. self.times,
  148. np.array(self),
  149. self.units,
  150. self.times.units,
  151. self.dtype,
  152. True,
  153. self.name,
  154. self.file_origin,
  155. self.description,
  156. self.annotations,
  157. self.segment,
  158. self.channel_index)
  159. def __array_finalize__(self, obj):
  160. '''
  161. This is called every time a new :class:`IrregularlySampledSignal` is
  162. created.
  163. It is the appropriate place to set default values for attributes
  164. for :class:`IrregularlySampledSignal` constructed by slicing or
  165. viewing.
  166. User-specified values are only relevant for construction from
  167. constructor, and these are set in __new__. Then they are just
  168. copied over here.
  169. '''
  170. super(IrregularlySampledSignal, self).__array_finalize__(obj)
  171. self.times = getattr(obj, 'times', None)
  172. # The additional arguments
  173. self.annotations = getattr(obj, 'annotations', None)
  174. # Globally recommended attributes
  175. self.name = getattr(obj, 'name', None)
  176. self.file_origin = getattr(obj, 'file_origin', None)
  177. self.description = getattr(obj, 'description', None)
  178. # Parent objects
  179. self.segment = getattr(obj, 'segment', None)
  180. self.channel_index = getattr(obj, 'channel_index', None)
  181. def __repr__(self):
  182. '''
  183. Returns a string representing the :class:`IrregularlySampledSignal`.
  184. '''
  185. return '<%s(%s at times %s)>' % (self.__class__.__name__,
  186. super(IrregularlySampledSignal,
  187. self).__repr__(), self.times)
  188. def __getslice__(self, i, j):
  189. '''
  190. Get a slice from :attr:`i` to :attr:`j`.
  191. Doesn't get called in Python 3, :meth:`__getitem__` is called instead
  192. '''
  193. return self.__getitem__(slice(i, j))
  194. def __getitem__(self, i):
  195. '''
  196. Get the item or slice :attr:`i`.
  197. '''
  198. obj = super(IrregularlySampledSignal, self).__getitem__(i)
  199. if isinstance(i, (int, np.integer)): # a single point in time across all channels
  200. obj = pq.Quantity(obj.magnitude, units=obj.units)
  201. elif isinstance(i, tuple):
  202. j, k = i
  203. if isinstance(j, (int, np.integer)): # a single point in time across some channels
  204. obj = pq.Quantity(obj.magnitude, units=obj.units)
  205. else:
  206. if isinstance(j, slice):
  207. obj.times = self.times.__getitem__(j)
  208. elif isinstance(j, np.ndarray):
  209. raise NotImplementedError("Arrays not yet supported")
  210. else:
  211. raise TypeError("%s not supported" % type(j))
  212. if isinstance(k, (int, np.integer)):
  213. obj = obj.reshape(-1, 1)
  214. elif isinstance(i, slice):
  215. obj.times = self.times.__getitem__(i)
  216. else:
  217. raise IndexError("index should be an integer, tuple or slice")
  218. return obj
  219. @property
  220. def duration(self):
  221. '''
  222. Signal duration.
  223. (:attr:`times`[-1] - :attr:`times`[0])
  224. '''
  225. return self.times[-1] - self.times[0]
  226. @property
  227. def t_start(self):
  228. '''
  229. Time when signal begins.
  230. (:attr:`times`[0])
  231. '''
  232. return self.times[0]
  233. @property
  234. def t_stop(self):
  235. '''
  236. Time when signal ends.
  237. (:attr:`times`[-1])
  238. '''
  239. return self.times[-1]
  240. def __eq__(self, other):
  241. '''
  242. Equality test (==)
  243. '''
  244. return (super(IrregularlySampledSignal, self).__eq__(other).all() and
  245. (self.times == other.times).all())
  246. def _check_consistency(self, other):
  247. '''
  248. Check if the attributes of another :class:`IrregularlySampledSignal`
  249. are compatible with this one.
  250. '''
  251. # if not an array, then allow the calculation
  252. if not hasattr(other, 'ndim'):
  253. return
  254. # if a scalar array, then allow the calculation
  255. if not other.ndim:
  256. return
  257. # dimensionality should match
  258. if self.ndim != other.ndim:
  259. raise ValueError('Dimensionality does not match: %s vs %s' %
  260. (self.ndim, other.ndim))
  261. # if if the other array does not have a times property,
  262. # then it should be okay to add it directly
  263. if not hasattr(other, 'times'):
  264. return
  265. # if there is a times property, the times need to be the same
  266. if not (self.times == other.times).all():
  267. raise ValueError('Times do not match: %s vs %s' %
  268. (self.times, other.times))
  269. def _copy_data_complement(self, other):
  270. '''
  271. Copy the metadata from another :class:`IrregularlySampledSignal`.
  272. '''
  273. for attr in ("times", "name", "file_origin",
  274. "description", "annotations"):
  275. setattr(self, attr, getattr(other, attr, None))
  276. def __rsub__(self, other, *args):
  277. '''
  278. Backwards subtraction (other-self)
  279. '''
  280. return self.__mul__(-1) + other
  281. def _repr_pretty_(self, pp, cycle):
  282. '''
  283. Handle pretty-printing the :class:`IrregularlySampledSignal`.
  284. '''
  285. pp.text("{cls} with {channels} channels of length {length}; "
  286. "units {units}; datatype {dtype} ".format(
  287. cls=self.__class__.__name__,
  288. channels=self.shape[1],
  289. length=self.shape[0],
  290. units=self.units.dimensionality.string,
  291. dtype=self.dtype))
  292. if self._has_repr_pretty_attrs_():
  293. pp.breakable()
  294. self._repr_pretty_attrs_(pp, cycle)
  295. def _pp(line):
  296. pp.breakable()
  297. with pp.group(indent=1):
  298. pp.text(line)
  299. for line in ["sample times: {0}".format(self.times)]:
  300. _pp(line)
  301. @property
  302. def sampling_intervals(self):
  303. '''
  304. Interval between each adjacent pair of samples.
  305. (:attr:`times[1:]` - :attr:`times`[:-1])
  306. '''
  307. return self.times[1:] - self.times[:-1]
  308. def mean(self, interpolation=None):
  309. '''
  310. Calculates the mean, optionally using interpolation between sampling
  311. times.
  312. If :attr:`interpolation` is None, we assume that values change
  313. stepwise at sampling times.
  314. '''
  315. if interpolation is None:
  316. return (self[:-1]*self.sampling_intervals.reshape(-1, 1)).sum()/self.duration
  317. else:
  318. raise NotImplementedError
  319. def resample(self, at=None, interpolation=None):
  320. '''
  321. Resample the signal, returning either an :class:`AnalogSignal` object
  322. or another :class:`IrregularlySampledSignal` object.
  323. Arguments:
  324. :at: either a :class:`Quantity` array containing the times at
  325. which samples should be created (times must be within the
  326. signal duration, there is no extrapolation), a sampling rate
  327. with dimensions (1/Time) or a sampling interval
  328. with dimensions (Time).
  329. :interpolation: one of: None, 'linear'
  330. '''
  331. # further interpolation methods could be added
  332. raise NotImplementedError
  333. def rescale(self, units):
  334. '''
  335. Return a copy of the :class:`IrregularlySampledSignal` converted to the
  336. specified units
  337. '''
  338. to_dims = pq.quantity.validate_dimensionality(units)
  339. if self.dimensionality == to_dims:
  340. to_u = self.units
  341. signal = np.array(self)
  342. else:
  343. to_u = pq.Quantity(1.0, to_dims)
  344. from_u = pq.Quantity(1.0, self.dimensionality)
  345. try:
  346. cf = pq.quantity.get_conversion_factor(from_u, to_u)
  347. except AssertionError:
  348. raise ValueError('Unable to convert between units of "%s" \
  349. and "%s"' % (from_u._dimensionality,
  350. to_u._dimensionality))
  351. signal = cf * self.magnitude
  352. new = self.__class__(times=self.times, signal=signal, units=to_u)
  353. new._copy_data_complement(self)
  354. new.channel_index = self.channel_index
  355. new.segment = self.segment
  356. new.annotations.update(self.annotations)
  357. return new
  358. def merge(self, other):
  359. '''
  360. Merge another :class:`IrregularlySampledSignal` with this one, and return the
  361. merged signal.
  362. The :class:`IrregularlySampledSignal` objects are concatenated horizontally
  363. (column-wise, :func:`np.hstack`).
  364. If the attributes of the two :class:`IrregularlySampledSignal` are not
  365. compatible, a :class:`MergeError` is raised.
  366. '''
  367. if not np.array_equal(self.times, other.times):
  368. raise MergeError("Cannot merge these two signals as the sample times differ.")
  369. if self.segment != other.segment:
  370. raise MergeError("Cannot merge these two signals as they belong to different segments.")
  371. if hasattr(self, "lazy_shape"):
  372. if hasattr(other, "lazy_shape"):
  373. if self.lazy_shape[0] != other.lazy_shape[0]:
  374. raise MergeError("Cannot merge signals of different length.")
  375. merged_lazy_shape = (self.lazy_shape[0], self.lazy_shape[1] + other.lazy_shape[1])
  376. else:
  377. raise MergeError("Cannot merge a lazy object with a real object.")
  378. if other.units != self.units:
  379. other = other.rescale(self.units)
  380. stack = np.hstack(map(np.array, (self, other)))
  381. kwargs = {}
  382. for name in ("name", "description", "file_origin"):
  383. attr_self = getattr(self, name)
  384. attr_other = getattr(other, name)
  385. if attr_self == attr_other:
  386. kwargs[name] = attr_self
  387. else:
  388. kwargs[name] = "merge(%s, %s)" % (attr_self, attr_other)
  389. merged_annotations = merge_annotations(self.annotations,
  390. other.annotations)
  391. kwargs.update(merged_annotations)
  392. signal = IrregularlySampledSignal(self.times, stack, units=self.units,
  393. dtype=self.dtype, copy=False,
  394. **kwargs)
  395. signal.segment = self.segment
  396. if hasattr(self, "lazy_shape"):
  397. signal.lazy_shape = merged_lazy_shape
  398. return signal
  399. def time_slice (self, t_start, t_stop):
  400. '''
  401. Creates a new :class:`IrregularlySampledSignal` corresponding to the time slice of
  402. the original :class:`IrregularlySampledSignal` between times
  403. `t_start` and `t_stop`. Either parameter can also be None
  404. to use infinite endpoints for the time interval.
  405. '''
  406. _t_start = t_start
  407. _t_stop = t_stop
  408. if t_start is None:
  409. _t_start = -np.inf
  410. if t_stop is None:
  411. _t_stop = np.inf
  412. indices = (self.times >= _t_start) & (self.times <= _t_stop)
  413. count = 0
  414. id_start = None
  415. id_stop = None
  416. for i in indices :
  417. if id_start == None :
  418. if i == True :
  419. id_start = count
  420. else :
  421. if i == False :
  422. id_stop = count
  423. break
  424. count += 1
  425. new_st = self[id_start:id_stop]
  426. return new_st