[frames] | no frames]

# Source Code for Module gmisclib.tsops

```  1  """Do mathematical operations on time series, where
2  the two operands don't necessarily have the same
3  sampling times.     It finds a common sampling time
4  and sampling interval, then interpolates as necessary
5  to bring the data onto the common time axis.
6  """
7
8  import math
9  from gmisclib import Num
10  from gmisclib import Numeric_gpk as NG
11  import gpkimgclass
12
13
14 -class axis:
15          """This class represents a time axis for a time series.
16          Indices of the underlying array are assumed to be zero-based.
17          """
18
19 -        def __init__(self, start=None, dt=None, n=None, end=None, crpix=0):
20                  if start is not None and end is not None and n is not None:
21                          self.crval = start
22                          self.n = int(round(n))
23                          self.cdelt = (end-start)/float(n-1)
24                  elif start is not None and n is not None and dt is not None:
25                          self.crval = start
26                          self.n = n
27                          self.cdelt = dt
28                  elif start is not None and dt is not None and end is not None:
29                          self.crval = start
30                          self.n = 1 + abs(int(math.floor((end-start)/dt)))
31                          # print '# self.n = 1+round(', (end-start)/dt, ')=', self.n
32                          if end > start:
33                                  self.cdelt = dt
34                          else:
35                                  self.cdelt = -dt
36                  else:
37                          raise ValueError, "Either silly values or not implemented."
38                  assert self.n >= 0, "Silly number of data"
39                  self.crpix = crpix
40                  # print 'info=', self.crval, 'at', self.crpix, 'dt=', self.cdelt, 'n=', self.n
41                  assert round(self.index(self.start())) == 0, 'start fail: %s' % self
42                  assert round(self.index(self.end())) == self.n-1, 'end fail: %s' % self
43
44 -        def N(self):
45                  return self.n
46
47 -        def coord(self, index):
48                  """What is the i^th coordinate if the i^th index==index?
49                  More plainly, this function gives you the coordinate that
50                  corresponds to the index."""
51                  assert self.n >= 0
52                  return self.crval + (index-self.crpix) * self.cdelt
53
54
55 -        def coords(self):
56                  """Generate an array of all the time values."""
57                  assert self.n >= 0
58                  index = Num.arrayrange(self.n)
59                  return index*self.cdelt + (self.crval-self.crpix*self.cdelt)
60
61
62 -        def index(self, t, limit=False, error=True):
63                  assert self.n > 0
64                  assert self.cdelt != 0
65                  tmp = self.crpix + (t-self.crval)/self.cdelt
66                  if limit:
67                          rtmp = round(tmp)
68                          if not ( rtmp >= 0 ):
69                                  tmp = 0.0
70                          elif not (rtmp < self.n):
71                                  tmp = self.n - 1.0
72                  elif error and not 0 <= round(tmp) < self.n:
73                          raise IndexError, "coordinate out of range: %g->%.1f" % (t, tmp)
74                  return tmp
75
76
77 -        def indices(self, t, limit=False, error=True):
78                  assert self.n > 0
79                  assert self.cdelt != 0
80                  tmp = self.crpix + (t-self.crval)/self.cdelt
81                  if limit:
82                          Num.maximum(tmp, 0.0, tmp)
83                          Num.minimum(tmp, self.n-1, tmp)
84                  elif error:
85                          nbad = Num.sum(Num.less(tmp, -0.5)) + Num.sum(Num.greater(tmp,self.n-0.5))
87                                  raise IndexError, "%d coordinates out of range" % nbad
88                  return tmp
89
90
91 -        def dt(self):
92                  return self.cdelt
93
94
95 -        def start(self):
96                  return self.coord(0)
97
98 -        def end(self):
99                  return self.coord(self.n-1)
100
101 -        def __str__(self):
102                  return '<axis start=%g end=%g dt=%g n=%d>' % (self.start(), self.end(),
103                                                                  self.dt(), self.N()
104                                                                  )
105          __repr__ = __str__
106
107
108 -def time(datasets, start=None, end=None):
109          dts = 0.0
110          dtn = 0
111          for x in datasets:
112                  tdt = x.dt()
113                  dts += math.log(tdt)
114                  dtn += 1
115          dt = math.exp(dts/dtn)
116
117          _s = [x.start() for x in datasets ]
118          if start is not None:
119                  _s.append(start)
120          _start = max(_s)
121          _e = [x.end() for x in datasets]
122          if end is not None:
123                  _e.append(end)
124          _end = min(_e)
125
126          n = int(math.floor( ( _end - _start ) / dt ) )
127          if n < 0:
128                  n = 0
129          return axis( start=_start, dt=dt, n=n )
130
131
132 -def time2(a, b):
133          dt = 0.5*( a.dt() + b.dt() )
134          start = max(a.start(), b.start())
135          end0 = min(a.end(), b.end())
136          n = int(math.floor((end0-start)/dt))
137          return axis(start=start, dt=dt, n=n)
138
139
140 -def _fill_interp_guts(a, t, fill, interpolator):
141          """Returns a Numeric array."""
142          assert len(t.shape) == 1
143          if a.n[2] == 0:
144                  # This case catches empty datasets.
145                  return Num.zeros((t.shape[0], a.n[1]), Num.Float) + fill
146          noneedfill = Num.greater_equal(t, a.start()) * Num.less_equal(t, a.end())
147          assert len(noneedfill.shape) == 1
148          # print 'nnf', noneedfill
149          nofillidx = Num.nonzero(noneedfill)[0]
150          # print 'nfi', nofillidx
151          if nofillidx.shape != t.shape:
152                  # print 't.shape=', t.shape
153                  ainofill = interpolator(a, Num.take(t, nofillidx, axis=0))
154                  # print 'ainf', ainofill
155                  out = Num.zeros((t.shape[0], a.n[1]), Num.Float) + fill
156                  # print 'out1=', out
157                  Num.put(out, nofillidx, ainofill)
158                  # print 'out2=', out
159                  # print 'outshape=', out.shape
160          else:
161                  out = interpolator(a, t)
162                  # print 'Outshape=', out.shape
163          return out
164
165
166 -def test_fig():
167          print "TEST_FIG"
168          a = axis(start=2.0, dt=1.0, n=10)
169          a.n = (0, 1, 10)
170          t = Num.arrayrange(20)
171          fill = -1
172          ifcn = lambda a, b: 100
173          q = _fill_interp_guts(a, t, fill, ifcn)
174          print q
175          print q.shape
176
177
178 -def interp_fill(a, t, fill):
179          return _fill_interp_guts(a, t, fill, interp)
180
181 -def interpN_fill(a, t, fill):
182          return _fill_interp_guts(a, t, fill, interpN)
183
184
185
186 -def interp(a, t):
187          """Interpolate to a specified time axis.
188          This does a linear interpolation.
189          @param a: data to be interpolated (a time series)
190          @type a: gpkimgclass.gpk_img
191          @type t: an array of times.
192          @rtype: numpy array.
193          @return: data interpolated onto the specified time values.
194          """
195          idx = a.t_index(t)
196          return NG.interp(a.d, idx)
197
198
199
200 -def test_interp1():
201          x = Num.array([0, 1], Num.Float)
202          xx = gpkimgclass.gpk_img({'CDELT2':1.0, 'CRPIX2':1, 'CRVAL2':0}, x)
203          t = Num.array([0.0, 0.7, 0.99, 1.0], Num.Float)
204          q = interp(xx, t)
205          print q.shape, t.shape, x.shape, xx.d.shape
206          print q
207          assert q.shape[0] == t.shape[0]
208          assert len(q.shape) == len(xx.d.shape)
209          assert Num.sum(Num.absolute(q-[[0.0], [0.7], [0.99], [1.0]]))<1e-6
210          qq = gpkimgclass.gpk_img({'CDELT2':1.0, 'CRPIX2':1, 'CRVAL2':0}, q)
211          print 'qq.n', qq.n, xx.n
212          assert qq.n[1] == xx.n[1] and qq.n[2] == t.shape[0]
213          qn = interpN(xx, t)
214          print 'qn=', qn
215          assert Num.sum(Num.absolute(qn-[[0.0], [1.0], [1.0], [1.0]]))<1e-6
216
217 -def test_interp2():
218          x = Num.array([0, 1, 1.5], Num.Float)
219          xx = gpkimgclass.gpk_img({'CDELT2':1.0, 'CRPIX2':1, 'CRVAL2':0}, x)
220          t = Num.array([0.0, 0.7, 0.99, 1.0, 1.1, 2.0], Num.Float)
221          q = interp(xx, t)
222          assert Num.sum(Num.absolute(q-[[0.0], [0.7], [0.99], [1.0], [1.05], [1.5]]))<1e-6
223          qn = interp(xx, t)
224          assert Num.sum(Num.absolute(qn-[[0.0], [0.7], [0.99], [1.0], [1.05], [1.5]]))<1e-6
225
226
227 -def test_interp3():
228          x = Num.array([[0, 100], [1, 101], [1.5, 101.5]], Num.Float)
229          xx = gpkimgclass.gpk_img({'CDELT2':1.0, 'CRPIX2':1, 'CRVAL2':0}, x)
230          t = Num.array([0.0, 0.7, 0.99, 1.0, 1.1, 2.0], Num.Float)
231          q = interp(xx, t)
232          print q
233          qq = gpkimgclass.gpk_img({'CDELT2':1.0, 'CRPIX2':1, 'CRVAL2':0}, q)
234          print 'qq.n', qq.n, xx.n
235          assert qq.n[1] == xx.n[1] and qq.n[2] == t.shape[0]
236          assert Num.sum(Num.absolute(q-[[0.0, 100.0], [0.7, 100.7], [0.99, 100.99],
237                                          [1.0, 101.0], [1.05, 101.05], [1.5, 101.5]]))<1e-6
238          qn = interp(xx, t)
239          assert Num.sum(Num.absolute(qn-[[0.0, 100.0], [0.7, 100.7], [0.99, 100.99],
240                                          [1.0, 101.0], [1.05, 101.05], [1.5, 101.5]]))<1e-6
241
242
243
244 -def interpN(a, t):
245          """Interpolate to a specified time axis via
246          nearest-neighbor interpolation.
247          A is a gpkimgclass, and t is an array of times.
248          Returns a Numeric array, not a gpkimgclass.
249          """
250          # print "interp.a=", a
251          # print "interp.t=", t
252          idx = a.t_index(t)
253          return NG.interpN(a.d, idx)
254
255
256
257 -def common(data_sets, start=None, end=None):
258          """Computes a common time axis for several datasets.
259          Linearly interpolate as needed.
260          The data sets need not be the same width, and need not have the
261          same sampling interval or a common starting time.
262          @param data_sets: this is a list of the data to be put on a common time axis.
263          @type data_sets: L{list}(L{gpkimgclass.gpk_img})
264          @param start: this allows you to restrict the output data to a smaller region.
265          @param end: this allows you to restrict the output data to a smaller region.
266          @type start: L{float} or L{None}
267          @type end: L{float} or L{None}
268          @rtype: list(numpy.ndarray) where the first ndarray is 1-D and the rest are two dimensional.
269          @return: the time_axis (as a 1-D numpy array of time values),
270                  followed by a 2-D numpy array for each of the input data sets.
271          """
272          tt = time(data_sets, start=start, end=end)
273          t = tt.coords()
274
275          return [tt] + [ interp(x, t) for x in data_sets ]
276
277
278
279 -def commonN(data_sets, start=None, end=None):
280          """Put several data sets on a common time axis.
281          Interpolate by choosing nearest neighbor.
282          """
283          tt = time(data_sets, start=start, end=end)
284          t = tt.coords()
285
286          return [t] + [ interpN(x, t) for x in data_sets ]
287
288
289 -def mul(a, b, hdr_op=None):
290          tt = time((a, b))
291          t = tt.coords()
292          # print "t=", t
293          ai = interp(a, t)
294          bi = interp(b, t)
295          # print "ai=", ai
296          # print "bi=", bi
297          c = ai * bi
298          # print "c=", c
299
300          if hdr_op is None:
301                  h = {}
302          else:
303                  h = hdr_op(a.hdr, b.hdr)
304
305          h['CDELT2'] = tt.dt()
306          h['CRVAL2'] = tt.start()
307          h['CRPIX2'] = 1
308
309          return gpkimgclass.gpk_img(h, c)
310
311
312
313 -def copy_interval(a, t0, t1, hdr_op=None, mode="rr"):
314          """This copies the part of the time-series in a
315          where t0 < t < t1.
316          'a' is a gpk_img object.
317          """
318
319          if hdr_op is None:
320                  h = a.hdr.copy()
321          else:
322                  h = hdr_op(a.hdr, t0, t1)
323          i0 = a.t_index(t0)
324          i1 = a.t_index(t1)
325          if mode[0] == "r":
326                  i0 = int(round(i0))
327          elif mode[0] == "w":
328                  i0 = int(math.floor(i0))
329          elif mode[0] == "n":
330                  i0 = int(math.ceil(i0))
331          if mode[1] == "r":
332                  i1 = int(round(i1))
333          elif mode[1] == "w":
334                  i1 = int(math.floor(i1))
335          elif mode[1] == "n":
336                  i1 = int(math.ceil(i1))
337
338          d = Num.array(a.d[i0:i1,:], Num.Float, copy=True)
339          h['CRVAL2'] = a.time(i0)
340          h['CRPIX2'] = 1
341          return gpkimgclass.gpk_img(h, d)
342
343
344 -def apply(fcn, a, hdrfcn=lambda x:x):
345          """Apply a function, point-by-point to the data in a."""
346          return gpkimgclass.gpk_img(hdrfcn(a.hdr), fcn(a.d))
347
348
349 -def resample(a, dt):
350          import samplerate
351          ratio = a.dt()/dt
352          d = samplerate.resample(a.d, ratio)
353          assert len(d.shape)==len(a.d.shape)
354          assert d.shape[1]==a.d.shape[1], "d.shape[0]: %d -> %d" % (a.d.shape[0], d.shape[0])
355          assert abs(float(d.shape[0])/float(a.d.shape[0]) - ratio) < 0.1
356          h = a.hdr.copy()
357          h['CDELT2'] = dt
358          return gpkimgclass.gpk_img(h, d)
359
360
361 -def test():
362          a = gpkimgclass.gpk_img({'CDELT2':1.0, 'CRPIX2':1, 'CRVAL2':0.0},
363                                  Num.arrayrange(10))
364          b = gpkimgclass.gpk_img({'CDELT2':0.5, 'CRPIX2':1, 'CRVAL2':0.4},
365                                  Num.arrayrange(10)*0.5+0.4)
366          # print 'a=', a.d
367          # print 'b=', b.d
368          # print 'b.time()', b.time()
369          ab = mul(a, b)
370          # print 'ab=', ab.d
371          # print 'ab.time()=', ab.time()
372          err = ab.d - Num.transpose([ab.time()])**2
373          # print 'err=', err
374          assert Num.sum(err**2) < 1e-5
375
376
377
378  if __name__ == '__main__':
379          test_interp1()
380          test_interp2()
381          test_interp3()
382          test()
383
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:12 2011 http://epydoc.sourceforge.net