1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
#!/usr/bin/env python
"""
1d model representation for polynomials
References::
[1] Storn, R. and Price, K. Differential Evolution - A Simple and Efficient
Heuristic for Global Optimization over Continuous Spaces. Journal of Global
Optimization 11: 341-359, 1997.
[2] Storn, R. and Price, K.
(Same title as above, but as a technical report.)
http://www.icsi.berkeley.edu/~storn/deshort1.ps
"""
from numpy import sum as numpysum
from numpy import polyval
def poly(x, c):
return polyval(c, x)
# coefficients for specific Chebyshev polynomials
chebyshev8coeffs = [128.0, 0.0, -256.0, 0.0, 160.0, 0.0, -32.0, 0.0, 1.0]
chebyshev16coeffs = [
32768.0,
0.0,
-131072.0,
0.0,
212992.0,
0.0,
-180224.0,
0.0,
84480.0,
0.0,
-21504.0,
0.0,
2688.0,
0.0,
-128.0,
0.0,
1,
]
class Chebyshev(Polynomial):
"""Chebyshev polynomial models and functions,
including specific methods for T8(z) & T16(z), Equation (27-33) of [2]
NOTE: default is T8(z)"""
def __init__(self, order=8, name="poly", metric=lambda x: numpysum(x * x), sigma=1.0):
Polynomial.__init__(self, name, metric, sigma)
if order == 8:
self.coeffs = chebyshev8coeffs
elif order == 16:
self.coeffs = chebyshev16coeffs
else:
raise NotImplementedError("provide self.coeffs 'by hand'")
return
def cost(self, trial, M=61):
"""The costfunction for order-n Chebyshev fitting.
M evaluation points between [-1, 1], and two end points""" # % (len(self.coeffs)-1)
# XXX: throw error when len(trial) != len(self.coeffs) ?
myCost = chebyshevcostfactory(self.coeffs)
return myCost(trial, M)
pass
# faster implementation
def chebyshevcostfactory(target):
def chebyshevcost(trial, M=61):
"""The costfunction for order-n Chebyshev fitting.
M evaluation points between [-1, 1], and two end points"""
result = 0.0
x = -1.0
dx = 2.0 / (M - 1)
for i in range(M):
px = polyeval(trial, x)
if px < -1 or px > 1:
result += (1 - px) * (1 - px)
x += dx
px = polyeval(trial, 1.2) - polyeval(target, 1.2)
if px < 0:
result += px * px
px = polyeval(trial, -1.2) - polyeval(target, -1.2)
if px < 0:
result += px * px
return result
return chebyshevcost
|
