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
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
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
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
46
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
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))
86 if nbad:
87 raise IndexError, "%d coordinates out of range" % nbad
88 return tmp
89
90
93
94
97
99 return self.coord(self.n-1)
100
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
138
139
141 """Returns a Numeric array."""
142 assert len(t.shape) == 1
143 if a.n[2] == 0:
144
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
149 nofillidx = Num.nonzero(noneedfill)[0]
150
151 if nofillidx.shape != t.shape:
152
153 ainofill = interpolator(a, Num.take(t, nofillidx, axis=0))
154
155 out = Num.zeros((t.shape[0], a.n[1]), Num.Float) + fill
156
157 Num.put(out, nofillidx, ainofill)
158
159
160 else:
161 out = interpolator(a, t)
162
163 return out
164
165
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
180
183
184
185
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
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
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
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
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
251
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
293 ai = interp(a, t)
294 bi = interp(b, t)
295
296
297 c = ai * bi
298
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
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
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
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
367
368
369 ab = mul(a, b)
370
371
372 err = ab.d - Num.transpose([ab.time()])**2
373
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