1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """Quadratic and cubic interpolation method.
25
26 This file is part of the U{minfx optimisation library<https://sourceforge.net/projects/minfx>}.
27 """
28
29
30 from math import sqrt
31
32
34 """Cubic interpolation using f(a), f(b), g(a), and g(b).
35
36 Equations
37 =========
38
39 The equations used are::
40
41 f(a) = a'a**3 + b'a**2 + c'a + d'
42 f(b) = a'b**3 + b'b**2 + c'b + d'
43 g(a) = 3a'a**2 + 2b'a + c'
44 g(b) = 3a'b**2 + 2b'b + c'
45
46
47 Interpolation
48 =============
49
50 The extrema are the roots of the quadratic equation::
51
52 3a'*alpha**2 + 2b'*alpha + c' = 0
53
54 The cubic interpolant is given by the formula::
55
56 g(b) + beta2 - beta1
57 ac = b - (b - a) . ---------------------
58 g(b) - g(a) + 2*beta2
59
60 where::
61
62 f(a) - f(b)
63 beta1 = g(a) + g(b) - 3 . -----------
64 a - b
65
66 if a < b:
67
68 beta2 = sqrt(beta1**2 - g(a).g(b))
69
70 else:
71
72 beta2 = -sqrt(beta1**2 - g(a).g(b))
73 """
74
75
76 beta1 = ga + gb - 3.0*(fa - fb)/(a - b)
77 if a < b:
78 beta2 = sqrt(beta1**2 - ga*gb)
79 else:
80 beta2 = -sqrt(beta1**2 - ga*gb)
81
82 denom = gb - ga + 2.0*beta2
83 if denom == 0.0:
84 return -1e99
85 else:
86 return b - (b - a)*(gb + beta2 - beta1)/denom
87
88
89 -def cubic_ext(a, b, fa, fb, ga, gb, full_output=0):
90 """Cubic Extrapolation using f(a), f(b), g(a), and g(b).
91
92 Extrapolation
93 =============
94
95 The extrema are the roots of the quadratic equation::
96
97 3a'*alpha**2 + 2b'*alpha + c' = 0
98
99 The cubic extrapolant is given by the formula::
100
101 g(b) + beta2 - beta1
102 ac = b - (b - a) . ---------------------
103 g(b) - g(a) + 2*beta2
104
105 where::
106
107 f(a) - f(b)
108 beta1 = g(a) + g(b) - 3 . -----------
109 a - b
110
111 if a < b:
112
113 beta2 = sqrt(max(0.0, beta1**2 - g(a).g(b)))
114
115 else:
116
117 beta2 = -sqrt(max(0.0, beta1**2 - g(a).g(b)))
118 """
119
120
121 beta1 = ga + gb - 3.0*(fa - fb)/(a - b)
122 if a < b:
123 beta2 = sqrt(max(0.0, beta1**2 - ga*gb))
124 else:
125 beta2 = -sqrt(max(0.0, beta1**2 - ga*gb))
126
127 denom = gb - ga + 2.0*beta2
128 if denom == 0.0:
129 alpha = -1e99
130 else:
131 alpha = b - (b - a)*(gb + beta2 - beta1)/denom
132
133 if full_output:
134 return alpha, beta1, beta2
135 else:
136 return alpha
137
138
140 """Quadratic interpolation using f(a), f(b), and g(a).
141
142 The extremum of the quadratic is given by::
143
144 1 g(a)
145 aq = a + - . -------------------------
146 2 f(a) - f(b) - (a - b)g(a)
147 """
148
149 denom = fa - fb - (a - b)*ga
150 if denom == 0.0:
151 return 1e99
152 else:
153 return a + 0.5 * ga / denom
154
155
157 """Quadratic interpolation using g(a) and g(b).
158
159 The extremum of the quadratic is given by::
160
161 bg(a) - ag(b)
162 aq = -------------
163 g(a) - g(b)
164 """
165
166 denom = ga - gb
167 if denom == 0.0:
168 return 1e99
169 else:
170 return (b*ga - a*gb) / denom
171