1
2 import numpy
3 import math
4 pylab = None
5
6
8 """Add arrays a and b in the overlap region.
9 Return (data, start).
10 If a, b are time series, they are assumed to have the
11 same sampling rate.
12 Astart and Bstart apply to the zeroth index.
13 All other indices are assumed to match start and length.
14 """
15 start = max(astart, bstart)
16 end = min(astart+a.shape[0], bstart+b.shape[0])
17
18 out = a[start-astart:end-astart] + b[start-bstart:end-bstart]
19
20 return (out, start)
21
22
23
24
25
26
27
28 asinh = numpy.arcsinh
29
30
31
32
33 acosh = numpy.arccosh
34
35
37 if padfactor == 0:
38 return numpy.array(d, copy=True)
39 elif padfactor<=0:
40 raise ValueError, 'pad factor <= 0'
41 assert len(d.shape) == 1
42 n = d.shape[0]
43 npad = int(round(n*padfactor))
44 return numpy.concatenate( (d, numpy.zeros((npad,), numpy.double)) )
45
46
48 """Return a Poisson random integer, whose distribution has
49 a mean = nbar.
50 """
51 import random
52 assert nbar >= 0.0
53 if nbar < 20.0:
54 L = math.exp(-nbar)
55 k = -1
56 p = 1.0
57 while p >= L:
58 k += 1
59 p *= random.random()
60 else:
61 lp = 0.0
62 lL = -nbar
63 k = 0
64 chunk = min( int(round (1 + nbar + 3*math.sqrt(nbar) )), 10000)
65 while lp >= lL:
66 ptmp = numpy.log(numpy.random.uniform(0.0, 1.0, chunk) )
67 lpp = numpy.add.accumulate(ptmp) + lp
68 fsm = numpy.nonzero(numpy.less(lpp, lL))[0]
69 if fsm.shape[0] != 0:
70 k += fsm[0]
71 break
72 k += ptmp.shape[0]
73 lp = lpp[-1]
74 return k
75
76
77
79 Nbar = 0.01
80 N = int(round(1000/Nbar))
81 s = 0
82 for i in range(N):
83 s += Poisson(Nbar)
84 assert abs(s-N*Nbar) < 5*math.sqrt(N*Nbar)
85
86
87 Nbar = 25.0
88 N = 1000
89 s = 0
90 for i in range(N):
91 s += Poisson(Nbar)
92 assert abs(s-N*Nbar) < 5*math.sqrt(N*Nbar)
93
94
95
96 -def bevel_concat(a, b, bevel=0, bevel_overlap=1.0, delay=0,
97 ta=None, tb=None):
98 """Concatenate two time series. Bevel the edges,
99 and overlap them slightly.
100
101 Bevel_overlap controls the fractional overlap of the two bevels,
102 and delay specifies an extra delay for b.
103
104 If ta and/or tb are specified, return a tuple of
105 (concatenated_time_series, tma, tmb) where tma and tmb are the locations
106 corresponding to ta and tb in the corresponding input arrays.
107 """
108 assert bevel >= 0
109 if bevel > 0:
110 bev = (0.5 + numpy.arange(bevel))/float(bevel)
111 else:
112 bev = 1
113 bc = numpy.array(b, copy=True)
114 bev = (0.5 + numpy.arange(bevel))/float(bevel)
115 ans = a.shape[0]-1
116 bns = bc.shape[0]-1
117 bos = a.shape[0] - int(round(bevel_overlap*bevel)) + delay
118 boe = bos + bc.shape[0]
119 o = numpy.zeros((boe,), numpy.double)
120 o[:a.shape[0]] = a
121 if bevel > 0:
122 numpy.multiply(o[:bevel], bev, o[:bevel])
123 numpy.multiply(o[ans:ans-bevel:-1], bev, o[ans:ans-bevel:-1])
124 numpy.multiply(bc[:bevel], bev, bc[:bevel])
125 numpy.multiply(bc[bns:bns-bevel:-1], bev, bc[bns:bns-bevel:-1])
126 numpy.add(o[bos:boe], bc, o[bos:boe])
127 if ta is not None or tb is not None:
128 if tb is not None:
129 tb += bos
130 return (o, ta, tb)
131 else:
132 return o
133
134
136 i = numpy.argmax(a, axis=None)
137 o = []
138 for n in reversed(a.shape):
139 o.append( i % n )
140 i = i//n
141 o.reverse()
142 return tuple(o)
143
144
146 x = numpy.array([1,2,3,4,3])
147 assert argmax(x) == (3,)
148 x = numpy.array([[1,2,3,4,3], [0,1,0,1,0]])
149 assert x[argmax(x)] == 4
150 x = numpy.array([[1,2,3,4,3], [0,1,0,5,0]])
151 assert x[argmax(x)] == 5
152 x = numpy.array([[1,2], [3,4], [5,6], [0,1],[0,5],[0,7]])
153 assert x[argmax(x)] == 7
154 x = numpy.zeros((5,4,3,3,1,2), numpy.double)
155 x[2,1,0,1,0,1] = 100.0
156 assert argmax(x) == (2,1,0,1,0,1)
157 x[0,0,0,0,0,0] = 200.0
158 assert argmax(x) == (0,0,0,0,0,0)
159 x[4,3,2,2,0,1] = 300.0
160 assert argmax(x) == (4,3,2,2,0,1)
161 x[4,3,2,2,0,0] = 400.0
162 assert argmax(x) == (4,3,2,2,0,0)
163 x[4,3,2,0,0,0] = 500.0
164 assert argmax(x) == (4,3,2,0,0,0)
165
166
168 assert len(a.shape) == 1
169 return a[numpy.argmax(a)].item()
170
171
173 assert len(a.shape) == 1
174 return a[numpy.argmin(a)].item()
175
176
177
178
180 assert 0 <= fr <= 1.0
181 tmp = numpy.sort(a)
182 n = tmp.shape[0]
183 assert n > 0, "Zero-sized array: cannot compute rank."
184 return tmp[int(round((n-1)*fr))].item()
185
186
188 """Mean absolute deviation. For a multi-dimensional array,
189 it takes the MAD along the first axis, so
190 N_mean_ad(x)[0]==N_mean_ad(x[:,0]).
191 """
192 ctr = numpy.median(a)
193 diff = numpy.subtract(a, ctr)
194 return numpy.sum(numpy.absolute(diff), axis=0)/(diff.shape[0]-1)
195
196
198 x = numpy.zeros((2,1), numpy.double)
199 x[0,0] = 1
200 x[1,0] = 0
201 assert numpy.allclose(N_mean_ad(x), [1.0])
202 x = numpy.zeros((5,7), numpy.double)
203 x[0,0] = 1
204 y = N_mean_ad(x)
205 assert y.shape == (7,)
206 assert numpy.allclose(y, [(1.0-0.0)/4.0, 0, 0, 0, 0, 0, 0])
207 assert abs(N_mean_ad(x)[0]-N_mean_ad(x[:,0])) < 0.001
208
209
222
223
224
239
242
243
244 variance = numpy.var
245 stdev = numpy.std
246 make_diag = numpy.diag
247
248
250 """Set the diagonal of a matrix x to be the vector a.
251 If a is shorter than the diagonal of x, just set the beginning."""
252
253 assert len(a.shape) == 1
254 assert len(x.shape) == 2
255 n = a.shape[0]
256 assert x.shape[0] >= n
257 assert x.shape[1] >= n
258 for i in range(n):
259 x[i,i] = a[i]
260
261
262
263
273
274
276 return numpy.clip(x, low, high)
277
278
279
281 import gpkavg
282 """Returns an element-by-element trimmed_mean of a list of Numeric vectors.
283 For instance, the first component of the output vector is the trimmed mean
284 of the first component of all of the input vectors.
285 """
286 import gpkavg
287 assert len(list_of_vec[0].shape) == 1
288 om = numpy.zeros(list_of_vec[0].shape, numpy.double)
289 osig = numpy.zeros(list_of_vec[0].shape, numpy.double)
290 tmp = numpy.zeros((len(list_of_vec),), numpy.double)
291 for t in list_of_vec:
292 assert t.shape == om.shape
293 for i in range(om.shape[0]):
294 for (j, v) in enumerate(list_of_vec):
295 tmp[j] = v[i]
296 om[i], osig[i] = gpkavg.avg(tmp, weights, clip)
297 return (om, osig)
298
301
304
305
307 """Take a component-by-component variance of a list of vectors."""
308 n = len(x)
309 if n < 2:
310 raise ValueError, "Cannot take variance unless len(x)>1"
311 x0 = x[0]
312 sh = x0.shape
313 s = numpy.zeros(sh, numpy.double)
314 ss = numpy.zeros(sh, numpy.double)
315 for xi in x:
316 assert xi.shape == sh
317 tmp = xi - x0
318 numpy.add(s, tmp, s)
319 numpy.add(ss, numpy.square(tmp), ss)
320 return numpy.maximum(ss-numpy.square(s)/n, 0.0)/(n-1)
321
322
323
324
340
341
342
344 d1 = numpy.asarray(d1)
345 d2 = numpy.asarray(d2)
346 assert len(d1.shape) == 1, 'd1.shape=%s' % str(d1.shape)
347 assert len(d2.shape) == 1, 'd2.shape=%s' % str(d2.shape)
348 if w1 is None:
349 w1 = numpy.ones(d1.shape, numpy.double)
350 if w2 is None:
351 w2 = numpy.ones(d2.shape, numpy.double)
352 ws1 = numpy.sum(w1)
353 ws2 = numpy.sum(w2)
354 i1 = numpy.argsort(d1)
355 i2 = numpy.argsort(d2)
356 c1 = 0.0
357 c2 = 0.0
358 j1 = 0
359 j2 = 0
360 maxdiff = 0.0
361 while True:
362 if abs(c1-c2) > maxdiff:
363 maxdiff = abs(c1 - c2)
364 if j1 < i1.shape[0]-1 and j2<i2.shape[0]-1:
365 if d1[i1[j1]] < d2[i2[j2]]:
366 j1 += 1
367 c1 += w1[i1[j1]]/ws1
368 elif d1[i1[j1]] > d2[i2[j2]]:
369 j2 += 1
370 c2 += w2[i2[j2]]/ws2
371 else:
372 j1 += 1
373 c1 += w1[i1[j1]]/ws1
374 j2 += 1
375 c2 += w2[i2[j2]]/ws2
376 elif j1==i1.shape[0]-1 and j2==i2.shape[0]-1:
377 break
378 elif j1 < i1.shape[0]-1:
379 j1 += 1
380 c1 += w1[i1[j1]]/ws1
381 else:
382 j2 += 1
383 c2 += w2[i2[j2]]/ws2
384 return maxdiff
385
386
388 assert KolmogorovSmirnov([1, 2, 3.01, 3, 4, 5], [1, 2, 3.01, 3, 4, 5]) < 0.001
389 assert abs( KolmogorovSmirnov([1, 2, 3.01, 3, 4, 5], [1, 2, 3, 4, 5, 6]) - 0.16667) < 0.001
390 print KolmogorovSmirnov([1, 2, 3.01, 3, 4, 5], [1, 2, 3, 4, 5])
391
392
393
394
396 """Interpolate array a to floating point indices, t,
397 via nearest-neighbor interpolation.
398 Returns a Numeric array.
399 """
400 ii = numpy.around(t)
401 n = a.shape[0]
402 if not numpy.alltrue((ii>=0) * (ii<=n-1)):
403 raise IndexError, "Index out of range."
404 iii = ii.astype(numpy.int)
405 return numpy.take(a, iii, axis=0)
406
407
408
410 """Interpolate to a specified time axis.
411 This does a linear interpolation.
412 A is a Numpy array, and t is an array of times.
413 @return: interpolated values
414 @rtype: numpy array.
415 """
416 n = a.shape[0]
417 nt = t.shape[0]
418 las = len(a.shape)
419
420 if las == 1:
421
422 a = numpy.transpose([a])
423 m = a.shape[1]
424 ii = numpy.around(t)
425 assert ii.shape == t.shape
426
427
428
429 if not numpy.alltrue((ii>=0) * (ii<=n-1)):
430 idxmin = numpy.argmin(t)
431 idxmax = numpy.argmax(t)
432 raise IndexError, "Index out of range: min=%g[index=%d,int=%d] max=%g[index=%d,int=%d] vs. %d" % (t[idxmin], idxmin, ii[idxmin], t[idxmax], idxmax, ii[idxmax], n)
433
434 nearestT = ii.astype(numpy.int)
435 assert nearestT.shape == t.shape
436
437 nearestA = numpy.take(a, nearestT, axis=0)
438 assert nearestA.shape == (nt,m)
439 deltaT = t - nearestT
440 assert deltaT.shape == nearestT.shape
441 isupport = numpy.where((deltaT>=0)*(nearestT<n-1)+(nearestT<=0), 1, -1) + nearestT
442 assert isupport.shape == (nearestA.shape[0],)
443 assert isupport.shape == nearestT.shape
444 assert isupport.shape == deltaT.shape
445
446
447 support = numpy.take(a, isupport, axis=0)
448 assert support.shape == (nt,m)
449 if len(a.shape) > 1:
450 assert nearestA.shape == (nt,m)
451 deltaA = (deltaT/(isupport-nearestT))[:,numpy.newaxis] * (support-nearestA)
452 assert deltaA.shape[1] == a.shape[1]
453 else:
454 deltaA = (deltaT/(isupport-nearestT)) * (support-nearestA)
455 assert deltaA.shape == nearestA.shape
456 rv = nearestA + deltaA
457 if las == 1:
458 return rv[:,0]
459 return rv
460
461
463 a = numpy.array([[1.0], [2.0]])
464
465 t = numpy.array([0.5])
466 q = interp(a, t)
467 assert q.shape == (1,1)
468 assert numpy.sum(numpy.absolute(interp(a, t) - [1.5])) < 1e-3
469
470
472 """This reports when the signal is above the threshold.
473 @param x: a signal
474 @type x: L{numpy.ndarray}, one-dimensional.
475 @param threshold: a threshold.
476 @type threshold: float
477 @returns: [(start, stop), ...] for each region ("clump") where C{x>threshold}.
478 """
479 assert len(x.shape) == 1
480 gt = numpy.greater(x, threshold)
481 chg = gt[1:] - gt[:-1]
482 nz = numpy.nonzero(chg)[0]
483 rv = []
484 if gt[0]:
485 last = 0
486 else:
487 last = None
488 for i in nz:
489 if last is None:
490 last = i+1
491 else:
492 if minsize > 0:
493 rv.append( (last, i+1) )
494 last = None
495 if last is not None and gt.shape[0]>=last+minsize:
496 rv.append( (last, gt.shape[0]) )
497 return rv
498
499
501 tmp = split_into_clumps(numpy.array([0,0,0,1,0,0]), 0.5)
502 assert tmp == [(3,4)]
503 tmp = split_into_clumps(numpy.array([1,0,0]), 0.5)
504 assert tmp == [(0,1)]
505 tmp = split_into_clumps(numpy.array([0,0,0,1]), 0.5)
506 assert tmp == [(3,4)]
507 tmp = split_into_clumps(numpy.array([1]), 0.5)
508 assert tmp == [(0,1)]
509 tmp = split_into_clumps(numpy.array([0]), 0.5)
510 assert tmp == []
511 tmp = split_into_clumps(numpy.array([1]), 0.5, minsize=2)
512 assert tmp == []
513 tmp = split_into_clumps(numpy.array([0,0,0,1,1,0]), 0.5, minsize=2)
514 assert tmp == [(3,5)]
515
516
517 BLOCK = 8192
518
520 """This is just a alternative implementation of
521 the standard deviation of each channel, but it is designed
522 in a block-wise fashion so the total memory usage is
523 not large.
524 """
525 m = (x.shape[0]+BLOCK-1)//BLOCK
526 avg = numpy.zeros((m, x.shape[1],))
527 sum2 = numpy.zeros((x.shape[1],))
528 nn = numpy.zeros((m,))
529 i = 0
530 j = 0
531 while i < x.shape[0]:
532 n = min(BLOCK, x.shape[0]-i)
533 tmp = numpy.sum(x[i:i+BLOCK], axis=0)/n
534 avg[j,:] = tmp
535 numpy.add(sum2, numpy.sum(numpy.square(x[i:i+BLOCK]-tmp), axis=0), sum2)
536 nn[j] = n
537 i += BLOCK
538 j += 1
539 assert j == m
540 avgavg = numpy.sum(nn[:,numpy.newaxis]*avg, axis=0)/float(x.shape[0])
541 for i in range(m):
542 numpy.add(sum2, nn[i]*numpy.square(avg[i,:]-avgavg), sum2)
543 return numpy.sqrt(sum2/float(x.shape[0]-1))
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
562 """This is basically like numpy.convolve(x, kernel, 1) except that it
563 properly handles the case where the kernel is longer than the data.
564 """
565 assert len(kernel.shape) == 1
566 assert len(x.shape) == 1, "Have not tried multidimensional data yet."
567 if x.shape[0] > kernel.shape[0]:
568 return numpy.convolve(x, kernel, 1)
569 n = kernel.shape[0] - x.shape[0] + 1
570 z = numpy.zeros((n,))
571 xx = numpy.concatenate((z, x, z), axis=0)
572 return numpy.convolve(xx, kernel, 1)[n:n+x.shape[0]]
573
574
578
579
580 -def edge_window(n, eleft, eright, typeleft="linear", typeright="linear", norm=None):
581 """Creates a window which is basically flat, but tapers off on the left and
582 right edges. The widths of the tapers can be controlled, as can the shapes.
583 @param n: the width of the window
584 @type n: C{int}
585 @param eleft: the width of the left taper (i.e. at zero index, in samples).
586 @type eleft: C{int}
587 @param eright: the width of the right taper (i.e. at index near C{n}, in samples).
588 @type eright: C{int}
589 @param typeleft: what kind of taper on the left? (Defaults to C{"linear"}).
590 @type typeleft: C{'linear'} or C{'cos'}
591 @param typeright: what kind of taper on the right? (Defaults to C{"linear"}).
592 @type typeright: C{'linear'} or C{'cos'}
593 @param norm: How to normalize? The default is L{None}, which means no normalization.
594 Providing a number C{x} will normalize the window so that the average of the
595 C{window**x}==1.
596 @type norm: L{None} or C{float != 0}.
597 """
598 assert int(n) >= 0
599 if not (eleft <= n and eright <= n):
600 raise EdgesTooWide, "n=%d, eleft=%d eright=%d" % (n, eleft, eright)
601 o = numpy.ones((n,))
602 if eleft > 0:
603 frac = (numpy.arange(eleft)+0.5)/eleft
604 if typeleft == "cos":
605 frac = (1.0 - numpy.cos(frac*math.pi))/2.0
606 o[:eleft] = frac
607 del frac
608 if eright > 0:
609 frac = ((eright-numpy.arange(eright))-0.5)/eright
610 if typeright == "cos":
611 frac = (1.0 - numpy.cos(frac*math.pi))/2.0
612 numpy.multiply(o[n-eright:], frac, o[n-eright:])
613 if norm is not None:
614 numpy.divide(o, numpy.average(o**norm)**(1.0/norm), o)
615 return o
616
617
618
619 -def edge_window_t(t, eleft, eright, typeleft="linear", typeright="linear", norm=None):
620 """Computes a window which is basically flat, but tapers off on the left and
621 right edges. The widths of the tapers can be controlled, as can the shapes.
622 @param t: an array of time values. They are required to be monotonically increasing,
623 and assumed to be linearly spaced.
624 @type t: C{numpy.ndarray}
625 @param eleft: the width of the left taper (i.e. at zero index, in time units).
626 @type eleft: C{float}
627 @param eright: the width of the right taper (i.e. at index near C{n}, in time units).
628 @type eright: C{float}
629 @param typeleft: what kind of taper on the left? (Defaults to C{"linear"}).
630 @type typeleft: C{'linear'} or C{'cos'}
631 @param typeright: what kind of taper on the right? (Defaults to C{"linear"}).
632 @type typeright: C{'linear'} or C{'cos'}
633 @param norm: How to normalize? The default is L{None}, which means no normalization.
634 Providing a number C{x} will normalize the window so that the average of the
635 C{window**x}==1.
636 @type norm: L{None} or C{float != 0}.
637 """
638 if not len(t)>1:
639 raise ValueError, "Time axis too short: len(t)=%d" % len(t)
640 t0 = float(t[0])
641 te = float(t[-1])
642 assert te > t0, "t=%s" % str(t)
643 w = te - t0
644 assert w > 0.0
645 if not ( eleft <= w and eright <= w):
646 return EdgesTooWide, "eleft=%s eright=%s w=%s" % (eleft, eright, w)
647 o = numpy.ones(t.shape)
648 if eleft > 0:
649 nedge = int(math.ceil((eleft/w)*t.shape[0]))
650 wedge = (nedge+1)*eleft/float(nedge)
651 tx = t0 - 0.5*wedge/(nedge+1)
652 if typeleft == 'cos':
653 f = (1.0 - numpy.cos((math.pi/wedge)*(t[:nedge]-tx)))/2.0
654 else:
655 f = (t[:nedge]-tx)/wedge
656 o[:nedge] = f
657 del f
658 if eright > 0:
659 nedge = int(math.ceil((eright/w)*t.shape[0]))
660 wedge = (nedge+1)*eright/float(nedge)
661 tx = te + 0.5*wedge/(nedge+1)
662 y = o.shape[0] - nedge
663 if typeright == 'cos':
664 f = (1.0 - numpy.cos((math.pi/wedge)*(tx-t[y:])))/2.0
665 else:
666 f = (tx-t[y:])/wedge
667 numpy.multiply(o[y:], f, o[y:])
668 if norm is not None:
669 numpy.divide(o, numpy.average(o**norm)**(1.0/norm), o)
670 return o
671
672
680
681
682 if __name__ == '__main__':
683 _test_split_into_clumps()
684 _test_argmax()
685 _test_interp1()
686 _test_Poisson()
687 _test_N_mean_ad()
688 _test_N_median_across()
689
690 _test_split_into_clumps()
691 import pylab
692 _test_ew()
693