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
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
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
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
60 registry = {}
61
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
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
98
99 return self._wt
100
101
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
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()
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
128
129
130 return numpy.array(self.o[-1], numpy.float, copy=True)
131
132
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)
138 numpy.add(o, tmp, o)
139 return o
140
142 raise RuntimeError, 'Virtual Function'
143
144
146 raise RuntimeError, "Virtual Function"
147
148
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
170
171
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()
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
196
197
198 return numpy.array(self.o[-1], numpy.float, copy=True)
199
200
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
207 raise RuntimeError, "Compute is not needed and should not be used for subclasses of ortho_poly."
208
209
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
220
221
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
230
231
232 return self.x*fn - 0.08*fnm1
233
234
235
237 return numpy.ones(self.x.shape, numpy.float)
238
239 ortho.registry[Legendre.name] = Legendre
240
241
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
253
254
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
263 return self.x*fn - 0.052*fnm1
264
265
267 return numpy.sqrt(1/(1.0-self.x**2))
268
269 ortho.registry[Chebyshev.name] = Chebyshev
270
271
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
282
283
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
292 return self.x*fn - 0.445*fnm1
293
294
296 return numpy.sqrt(1.0-self.x**2)
297
298 ortho.registry[Chebyshev2.name] = Chebyshev2
299
300
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
309
310
312
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
321 return numpy.ones(self.x.shape, numpy.float)
322
323 ortho.registry[SinCos.name] = SinCos
324
325
326
327
328
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
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
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
359
360
361 return self.r(x) * self.r(1.0-x)
362
363
365 return numpy.ones(self.x.shape, numpy.float)
366
367
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
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
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
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
403
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
435
436
437 if __name__ == '__main__':
438 test()
439
440
441
442