[frames] | no frames]

# Source Code for Module gmisclib.Numeric_gpk

```  1
2  import numpy
3  import math
4  pylab = None
5
6
7 -def add_overlap(a, astart, b, bstart):
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  # def asinh(x):
26          # sign = numpy.sign(x)
27          # return sign*numpy.log(sign*x + numpy.sqrt(x*x + 1))
28  asinh = numpy.arcsinh
29
30  # def acosh(x):
31          # assert numpy.all(numpy.less_equal(x, 1))
32          # return numpy.log(x - numpy.sqrt(x*x-1))
33  acosh = numpy.arccosh
34
35
38                  return numpy.array(d, copy=True)
40                  raise ValueError, 'pad factor <= 0'
41          assert len(d.shape) == 1
42          n = d.shape[0]
44          return numpy.concatenate( (d, numpy.zeros((npad,),  numpy.double)) )
45
46
47 -def Poisson(nbar):
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
78 -def _test_Poisson():
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])
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
135 -def argmax(a):
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
145 -def _test_argmax():
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
167 -def N_maximum(a):
168          assert len(a.shape) == 1
169          return a[numpy.argmax(a)].item()
170
171
172 -def N_minimum(a):
173          assert len(a.shape) == 1
174          return a[numpy.argmin(a)].item()
175
176
177
178
179 -def N_frac_rank(a, fr):
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
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
202          x = numpy.zeros((5,7), numpy.double)
203          x[0,0] = 1
205          assert y.shape == (7,)
206          assert numpy.allclose(y, [(1.0-0.0)/4.0, 0, 0, 0, 0, 0, 0])
208
209
210 -def median_across(list_of_vec):
211          """Returns an element-by-element median of a list of Numeric vectors."""
212          assert len(list_of_vec[0].shape) == 1
213          o = numpy.zeros(list_of_vec[0].shape, numpy.double)
214          tmp = numpy.zeros((len(list_of_vec),), numpy.double)
215          for t in list_of_vec:
216                  assert t.shape == o.shape
217          for i in range(o.shape[0]):
218                  for (j, v) in enumerate(list_of_vec):
219                          tmp[j] = v[i]
220                  o[i] = numpy.median(tmp)
221          return o
222
223
224  # Median along the first axis.
225 -def N_median(a, axis=0):
226          """Returns an element-by-element median of a list of Numeric vectors."""
227          if len(a.shape) ==  1:
228                  return numpy.median(a)
229          assert len(a.shape) == 2
230          if axis == 1:
231                  o = numpy.zeros((a.shape[0],), numpy.double)
232                  for i in range(o.shape[0]):
233                          o[i] = numpy.median(a[i])
234          o = numpy.zeros((a.shape[1],), numpy.double)
235          a = numpy.array(a)
236          for i in range(o.shape[0]):
237                  o[i] = numpy.median(a[:,i])
238          return o
239
240 -def N_median_across(a):
241          return N_median(a, axis=0)
242
243
244  variance = numpy.var
245  stdev = numpy.std
246  make_diag = numpy.diag
247
248
249 -def set_diag(x, a):
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  # Take the median along the first coordinate of a 2-D matrix,
263  # so if m.shape=(2,4) then the output will have shape (4,)
264 -def _test_N_median_across():
265          x = numpy.zeros((3,2), numpy.double)
266          x[0,0] = 1
267          x[1,0] = 2
268          x[2,0] = 3
269          y = N_median(x, axis=0)
270          # y = N_median_across(x)
271          print 'median_across=', y
272          assert numpy.allclose(y, [2.0, 0.0])
273
274
275 -def limit(low, x, high):
276          return numpy.clip(x, low, high)
277
278
279
280 -def trimmed_mean_sigma_across(list_of_vec, weights, clip):
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
299 -def trimmed_mean_across(list_of_vec, weights, clip):
300          return trimmed_mean_sigma_across(list_of_vec, weights, clip)[0]
301
302 -def trimmed_stdev_across(list_of_vec, weights, clip):
303          return trimmed_mean_sigma_across(list_of_vec, weights, clip)[1]
304
305
306 -def vec_variance(x):
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
320          return numpy.maximum(ss-numpy.square(s)/n, 0.0)/(n-1)
321
322
323
324
325 -def qform(vec, mat):
327          or vecs*mat*transpose(vecs)"""
328          if len(vec.shape) != 1 or len(mat.shape) != 2:
329                  raise ValueError, ': '.join(["Can't handle input",
330                                               "requires vector*matrix(vector)",
331                                               "shapes are %s and %s" % (
332                                                          vec.shape, mat.shape)
333                                                  ])
334          if mat.shape != (vec.shape[0], vec.shape[0]):
335                  raise ValueError, ': '.join([
336                                          "Matrix must be square and match the length of vector",
337                                          "shapes are %s and %s" % (vec.shape, mat.shape)
338                                          ])
339          return numpy.dot(vec, numpy.dot(mat, vec))
340
341
342
343 -def KolmogorovSmirnov(d1, d2, w1=None, w2=None):
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
387 -def _testKS():
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
395 -def interpN(a, t):
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
409 -def interp(a, t):
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          # print 'a=', a.shape, a
420          if las == 1:
421                  # print 'a.shape=', a.shape
422                  a = numpy.transpose([a])
423          m = a.shape[1]
424          ii = numpy.around(t)
425          assert ii.shape == t.shape
426
427          # print 't=', t
428          # print 'n=', n, 'ii=', ii
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          # print 'a=',a
446          # print 'isupport=', isupport
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
462 -def _test_interp1():
463          a = numpy.array([[1.0], [2.0]])
464          # a = numpy.array([1.0, 2.0], numpy.double)
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
471 -def split_into_clumps(x, threshold, minsize=1):
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
500 -def _test_split_into_clumps():
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
519 -def block_stdev(x):
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
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):
543          return numpy.sqrt(sum2/float(x.shape[0]-1))
544
545
546  # def _test_block_stdev():
547          # n = BLOCK//2
548          # while n < 4*BLOCK:
549                  # x = numpy.zeros((n,2))
550                  # x[0,0] = 1.0
551                  # x[0,1] = 2.0
552                  # x[1,0] = -1.0
553                  # avg2 = 2.0/float(n)
554                  # a2 = math.sqrt(((2-avg2)**2+(n-1)*avg2**2)/float(n-1))
555                  # answer = [math.sqrt(2/float(n-1)), a2]
556                  # calc = avg_dev(x)
557                  # assert numpy.sum(numpy.absolute(numpy.log(calc/answer))) < 0.001
558                  # n = (n*3)//2 + 17
559
560
561 -def convolve(x, kernel):
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
575 -class EdgesTooWide(Exception):
576 -        def __init__(self, *s):
577                  Exception.__init__(self, *s)
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
673 -def _test_ew():
674          t = 1.0 + 0.1*numpy.arange(100)
675          for i in range(0,60,4):
676                  e = 0.0 + 0.1*i
677                  w = edge_window_t(t, e, 1.432*e, typeright="cos")
678                  pylab.plot(t, w)
679          pylab.show()
680
681
682  if __name__ == '__main__':
683          _test_split_into_clumps()
684          _test_argmax()
685          _test_interp1()
686          _test_Poisson()