[frames] | no frames]

# Source Code for Module gmisclib.ortho_poly

```  1  import math
2  import numpy
3  from gmisclib import die
4
5
6  """Implements a set of orthonormal polynomials, based on a recurrence relation
7  on the order of the polynomial.    Each particular type of polynomial just
8  needs to specify that recurrence relation.
9  """
10
11
12 -def _orthogonalize(x, others, wt):
13          """Makes x orthonormal to all others.
14          Others are assumed to be orthonormal to eachother."""
15          against = others[:]
16          against.reverse()
17          ipass = 0
18          x = numpy.array(x, numpy.float, copy=True)
19          n = x.shape[0]
20          while ipass <= 2:
21                  ipass += 1
22                  do_all_again = False
23                  do_again = []
24                  renorm = True
25                  for t in against:
26                          if renorm > 0:
27                                  numpy.multiply(x, 1.0/math.sqrt(numpy.average( wt * x**2)), x)
28                                  renorm = False
29                          d = numpy.average( x * wt * t)
30                          # print 'i=', i, 'd=', d
31                          if abs(d) >= 1e-6:
32                                  do_again.append(t)
33                                  renorm = True
34
35                          if abs(d) > 0.1*n:
36                                  do_all_again = True
37                                  if ipass > 1 :
38                                          die.note('x', x)
39                                          die.note('t', t)
40                                          die.note('wt', wt)
41                                          die.note('d', d)
42                                          die.warn('Poor orthogality')
43                                          raise RuntimeError, "Poor orthogonality: %d %g %d/%d" % (len(wt), d, ipass, len(others))
44                          assert x.shape == t.shape
45                          # assert type(d) == type(0.1)
46                          numpy.subtract(x, t*d, x)
47
48                  if do_all_again:
49                          against = others
50                  elif do_again:
51                          against = do_again
52                  else:
53                          break
54          if ipass > 2:
55                  raise RuntimeError, "Can't orthogonalize"
56          return x * (1.0/math.sqrt(numpy.average(wt*x**2)))
57
58
59 -class ortho:
60          registry = {}
61
62 -        def __init__(self, n=None, x=None):
63                  """Argument n is how many points you want the polynomials evaluated at.
64                  The ordinates of the points are placed in self.x.
65                  You may alternatively specify the points as argument x, in which case
66                  the x_() function will never be called."""
67                  if x is not None:
68                          lx = len(x)
69                          if not lx > 0:
70                                  raise ValueError, "Need len(x)>0"
71                          if n is not None:
72                                  assert int(n) == len(x), "X and N mismatch: %d vs %s" % (len(x), n)
73                          self.x = numpy.array(x, numpy.float, copy=True)
74                  elif n is not None:
75                          n = int(n)
76                          if not n>0:
77                                  raise ValueError, "Need n>0"
78                          self.x = self.x_(n)
79                  else:
80                          raise RuntimeError, "Need to specify either x or n."
81
82                  self.o = []
83                  self._wt = None
84                  self._numwt = None
85
86
87 -        def wt(self):
88                  """Weighting function to get orthonormality.
89                  This sets self._numwt."""
90                  if self._wt is None:
91                          self._wt = self.wt_()
92                          self._numwt = numpy.sum( self._wt > 0.0 )
93                          if not numpy.alltrue( self._wt >= 0.0 ):
94                                  raise ValueError, 'Weights must be non-negative.'
95                          if not self._numwt > 0:
96                                  raise ValueError, 'Need some positive weights.'
97                  # Again, it might make sense to return a copy, as _wt gets used
98                  # internally after being passed out...
99                  return self._wt
100
101
102 -        def x_(self, m):
103                  """Calculates the points at which the function is evaluated,
104                          if you want m evenly spaced points.
105                          By default, the function is assumed to
106                          range over the open interval (-1, 1).
107                          You can override this function in a subclass
108                          to get a different range."""
109                  return (numpy.arange(m)-0.5*(m-1))*(2.0/m)
110
111
112 -        def P(self, i):
113                  """Self.P(i) is the ith orthogonal polynomial.
114                  The result is normalized so that
115                  numpy.sum(self.wt() * self.P(i)**2) == 1.
116                  It returns an un-shared array, so the array can be
117                  modified without affecting future calls to P(i).
118                  """
119                  assert i>=0
120                  if i<len(self.o):
121                          return numpy.array(self.o[i], copy=True)
122                  wt = self.wt()  # Sets self._numwt
123                  if i >= self._numwt:
124                          raise ValueError, 'Request for higher order polynomial than the number of positive weight points.'
125                  while i >= len(self.o):
126                          self.o.append(_orthogonalize(self.compute(len(self.o)), self.o, wt))
127                  # We return a copy, because it is possible that the user
128                  # might modify the data.  If so, it would screw up the orthogalization
129                  # step for higher polynomials.  Copy-on-write would be nice.
130                  return numpy.array(self.o[-1], numpy.float, copy=True)
131
132
133 -        def expand(self, c):
134                  o = numpy.zeros(self.x.shape, numpy.float)
135                  for (i, ci) in enumerate(c):
136                          tmp = self.P(i)
137                          numpy.multiply(tmp, ci, tmp)
139                  return o
140
141 -        def compute(self, n):
142                  raise RuntimeError, 'Virtual Function'
143
144
145 -        def wt_(self):
146                  raise RuntimeError, "Virtual Function"
147
148
149 -class ortho_poly(ortho):
150          __doc__ = """Virtual base class.
151                  This is a superclass of all orthorgonal polynomials.
152                  This does the recurrence [Abramowitz and Stegun p.782],
153                  and ensures that the generated polynomials are all
154                  orthonormal to each other.
155                  The derived classed need to specify two functions:
156                  recurse() and wt_().
157                  Recurse() is the recursion relation from P(i) and P(i-1)
158                  to P(i+1), and
159                  wt_() is the weighting function for the polynomial.
160                  Wt_() is needed to check orthogonality, and really
161                  defines the polynomial.
162                  These polynomials are guarenteed to be
163                  orthogonal to eachother when summed over the
164                  supplied set of x points.   Thus, if you change
165                  x_() or wt_(), you get a different set of functions."""
166
167
168 -        def __init__(self, n=None, x=None):
169                  ortho.__init__(self, n, x)
170
171
172 -        def P(self, i):
173                  """Self.P(i) is the ith orthogonal polynomial.
174                  The result is normalized so that
175                  numpy.sum(self.P(i)**2) == 1.  """
176
177                  assert i>=0
178                  if i<len(self.o):
179                          return numpy.array(self.o[i], copy=True)
180                  wt = self.wt()  # Sets self._numwt
181                  if i >= self._numwt:
182                          raise ValueError, 'Request for higher order polynomial than the number of positive weight points.'
183                  while i>=len(self.o):
184                          if len(self.o)>1:
185                                  m2 = self.o[-2]
186                          else:
187                                  m2 = None
188                          if len(self.o)>0:
189                                  m1 = self.o[-1]
190                          else:
191                                  m1 = None
192
193                          tmp = self.recurse(len(self.o), m1, m2)
194                          self.o.append(_orthogonalize(tmp, self.o, wt))
195                  # We return a copy, because it is possible that the user
196                  # might modify the data.  If so, it would screw up the orthogalization
197                  # step for higher polynomials.  Copy-on-write would be nice.
198                  return numpy.array(self.o[-1], numpy.float, copy=True)
199
200
201 -        def recurse(self, n, fn, fnm1):
202                  """Does the recurrence relation, evaluating f_(n+1) as a function
203                  of n, f_n, and f_(n-1).   For n=0 and n=1, fn and fnm1 may be None."""
204                  raise RuntimeError, "Virtual function"
205
206 -        def compute(self, n):
207                  raise RuntimeError, "Compute is not needed and should not be used for subclasses of ortho_poly."
208
209
210 -class Legendre(ortho_poly):
211          __doc__ = """Legendre polynomials,
212                  orthonormal over (-1, 1) with weight 1.
213                  Called Pn(x) in Abramowitz and Stegun.
214                  """
215
216          name = 'Legendre'
217
218 -        def __init__(self, n=None, x=None):
219                  ortho_poly.__init__(self, n, x)
220
221
222 -        def recurse(self, n, fn, fnm1):
223                  if fn is None:
224                          assert n == 0
225                          return numpy.ones(self.x.shape, numpy.float)
226                  if fnm1 is None:
227                          assert n == 1
228                          return numpy.array(self.x, numpy.float, copy=True)
229                  # tmp1 = (2*n-1)/float(n)
230                  # tmp2 = (n-1)/float(n)
231                  # return tmp1*self.x*fn - tmp2*fnm1
232                  return self.x*fn - 0.08*fnm1
233
234
235
236 -        def wt_(self):
237                  return numpy.ones(self.x.shape, numpy.float)
238
239  ortho.registry[Legendre.name] = Legendre
240
241
242 -class Chebyshev(ortho_poly):
243          __doc__ = """Chebyshev polynomials of the first kind,
244                  orthonormal over (-1, 1) with weight (1-x^2)^(-1/2).
245                  These are the equi-ripple polynomials.
246                  Called Tn(x) in Abramowitz and Stegun.
247                  """
248
249          name = 'Chebyshev'
250
251 -        def __init__(self, n=None, x=None):
252                  ortho_poly.__init__(self, n, x)
253
254
255 -        def recurse(self, n, fn, fnm1):
256                  if fn is None:
257                          assert n == 0
258                          return numpy.ones(self.x.shape, numpy.float)
259                  if fnm1 is None:
260                          assert n == 1
261                          return numpy.array(self.x, numpy.float, copy=True)
262                  # return 2*self.x*fn - fnm1
263                  return self.x*fn - 0.052*fnm1
264
265
266 -        def wt_(self):
267                  return numpy.sqrt(1/(1.0-self.x**2))
268
269  ortho.registry[Chebyshev.name] = Chebyshev
270
271
272 -class Chebyshev2(ortho_poly):
273          __doc__ = """Chebyshev polynomials of the second kind,
274                  orthonormal over (-1, 1) with weight (1-x^2)^(1/2).
275                  Called Un(x) in Abramowitz and Stegun.
276                  """
277
278          name = 'Chebyshev2'
279
280 -        def __init__(self, n=None, x=None):
281                  ortho_poly.__init__(self, n, x)
282
283
284 -        def recurse(self, n, fn, fnm1):
285                  if fn is None:
286                          assert n == 0
287                          return numpy.ones(self.x.shape, numpy.float)
288                  if fnm1 is None:
289                          assert n == 1
290                          return 2*self.x
291                  # return 2*self.x*fn - fnm1
292                  return self.x*fn - 0.445*fnm1
293
294
295 -        def wt_(self):
296                  return numpy.sqrt(1.0-self.x**2)
297
298  ortho.registry[Chebyshev2.name] = Chebyshev2
299
300
301 -class SinCos(ortho):
302          __doc__ = """1, sin(pi*x), cos(pi*x), sin(2*pi*x), cos(2*pi*x)...
303                  orthonormal over (-1, 1) with weight 1."""
304
305          name = 'SinCos'
306
307 -        def __init__(self, n=None, x=None):
308                  ortho.__init__(self, n, x)
309
310
311 -        def compute(self, n):
312                  # print 'SinCos.recurse', n
313                  if n == 0:
314                          return numpy.ones(self.x.shape, numpy.float)
315                  elif n%2 == 1:
316                          return numpy.sin((math.pi*(n+1)/2) * self.x)
317                  return numpy.cos((math.pi*(n/2))*self.x)
318
319
320 -        def wt_(self):
321                  return numpy.ones(self.x.shape, numpy.float)
322
323  ortho.registry[SinCos.name] = SinCos
324
325
326
327
328
329 -class SLTB(ortho):
330          __doc__ = """Smooth Local Trigonometric Basis.
331                  From Bj\:orn Jawerth, Yi Liu, Wim Sweldens,
332                  "Signal Compression with Smooth Local Basis Functions".
333                  """
334
335          name = 'SLTB'
336
337 -        def __init__(self, n=None, x=None, eps=0.5):
338                  assert eps <= 0.5 and eps>0.0
339                  self.eps = eps
340                  ortho.__init__(self, n, x)
341
342
343 -        def compute(self, n):
344                  omega = (2*n+1)*math.pi/2.0
345                  s = numpy.sin(omega * self.x)
346                  return s * self.b(self.x)
347
348
349 -        def r(self, x):
350                  """This must have l**2+r**2==1"""
351                  x = numpy.clip(x, -self.eps, self.eps)
352                  s = numpy.sin(x * math.pi/(2.0*self.eps))
353                  r = 1 + s
354                  l = 1 - s
355                  return r/numpy.hypot(r, l)
356
357
358 -        def b(self, x):
359                  # print 'x=', x
360                  # print 'b=', self.l(x)*self.l(1.0-x)
361                  return self.r(x) * self.r(1.0-x)
362
363
364 -        def wt_(self):
365                  return numpy.ones(self.x.shape, numpy.float)
366
367
368 -        def x_(self, m):
369                  """Calculates the points at which the function is evaluated,
370                          if you want m evenly spaced points.
371                          The function is assumed to
372                          range over the open interval (-eps, 1+eps).
373                          You can override this function in a subclass
374                          to get a different range."""
375                  return (0.5+numpy.arange(m))*((1+2*self.eps)/m) - self.eps
376
377
378 -        def test():
379                  M = 20
380                  N = 3*M
381                  x = (numpy.arange(N)+0.5)/float(N) * 3.0 - 0.5
382                  q = SLTB(x=x)
383                  for i in range(M//3):
384                          for j in range(M//3):
385                                  a = q.P(i)
386                                  b = q.P(j)
387                                  # print numpy.dot(a[:-M], b[M:])
388                                  assert abs(numpy.dot(a[:-M], b[M:])) < 0.001
389
390          test = staticmethod(test)
391
392  ortho.registry[SLTB.name] = SLTB
393
394
395
396 -def test_a_name(name):
397          N = 96
398          q = F(name, N)
399          for i in range(N//2):
400                  for j in range(10):
401                          dot = numpy.sum(q.P(i)*q.P(j)*q.wt())
402                          # print 'p(i=', i, ')', q.P(i)
403                          # print 'p(j=', j, ')', q.P(j)
404                          if i != j:
405                                  assert abs(dot) < 0.001
406                          else:
407                                  assert abs(dot - N) < 0.001*N
408
409          if hasattr(q, 'test'):
410                  q.test()
411
412
413
414
415
416
417
418
419
420
421
422 -def F(name, n=None, x=None):
423          """This is a factory function.  Get any kind of ortho poly
424          you want, as selected by the first argument."""
425
426          try:
427                  return ortho.registry[name](n, x)
428          except KeyError:
429                  raise RuntimeError, "Unimplemented orthogonal function name: %s" % name
430
431
432 -def test():
433          for name in ortho.registry.keys():
434                  test_a_name(name)
435
436
437  if __name__ == '__main__':
438          test()
439          # q = F('SLTB', n=100)
440          # for i in range(100):
441                  # print i, q.P(0)[i]
442
<!--
expandto(location.href);
// -->

```

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