b311480e |
1 | // Copyright (c) 1997-1999 Matra Datavision |
973c2be1 |
2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
3 | // |
973c2be1 |
4 | // This file is part of Open CASCADE Technology software library. |
b311480e |
5 | // |
d5f74e42 |
6 | // This library is free software; you can redistribute it and/or modify it under |
7 | // the terms of the GNU Lesser General Public License version 2.1 as published |
973c2be1 |
8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
10 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
11 | // |
973c2be1 |
12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. |
b311480e |
14 | |
7fd59977 |
15 | |
42cf5bc1 |
16 | #include <math_BracketMinimum.hxx> |
7fd59977 |
17 | #include <math_Function.hxx> |
42cf5bc1 |
18 | #include <StdFail_NotDone.hxx> |
7fd59977 |
19 | |
42cf5bc1 |
20 | // waiting for NotDone Exception |
7fd59977 |
21 | #define GOLD 1.618034 |
22 | #define CGOLD 0.3819660 |
23 | #define GLIMIT 100.0 |
24 | #define TINY 1.0e-20 |
25 | #ifdef MAX |
26 | #undef MAX |
27 | #endif |
28 | #define MAX(a,b) ((a) > (b) ? (a) : (b)) |
29 | #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) |
30 | #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d) |
31 | |
f79b19a1 |
32 | Standard_Boolean math_BracketMinimum::LimitAndMayBeSwap |
33 | (math_Function& F, |
34 | const Standard_Real theA, |
35 | Standard_Real& theB, |
36 | Standard_Real& theFB, |
37 | Standard_Real& theC, |
38 | Standard_Real& theFC) const |
39 | { |
40 | theC = Limited(theC); |
41 | if (Abs(theB - theC) < Precision::PConfusion()) |
42 | return Standard_False; |
43 | Standard_Boolean OK = F.Value(theC, theFC); |
44 | if (!OK) |
45 | return Standard_False; |
46 | // check that B is between A and C |
47 | if ((theA - theB) * (theB - theC) < 0) |
48 | { |
49 | // swap B and C |
50 | Standard_Real dum; |
51 | SHFT(dum, theB, theC, dum); |
52 | SHFT(dum, theFB, theFC, dum); |
53 | } |
54 | return Standard_True; |
55 | } |
7fd59977 |
56 | |
f79b19a1 |
57 | void math_BracketMinimum::Perform(math_Function& F) |
58 | { |
7fd59977 |
59 | |
60 | Standard_Boolean OK; |
f79b19a1 |
61 | Standard_Real ulim, u, r, q, fu, dum; |
7fd59977 |
62 | |
63 | Done = Standard_False; |
7fd59977 |
64 | Standard_Real Lambda = GOLD; |
65 | if (!myFA) { |
66 | OK = F.Value(Ax, FAx); |
67 | if(!OK) return; |
68 | } |
69 | if (!myFB) { |
70 | OK = F.Value(Bx, FBx); |
71 | if(!OK) return; |
72 | } |
73 | if(FBx > FAx) { |
74 | SHFT(dum, Ax, Bx, dum); |
75 | SHFT(dum, FBx, FAx, dum); |
76 | } |
f79b19a1 |
77 | |
78 | // get next prob after (A, B) |
7fd59977 |
79 | Cx = Bx + Lambda * (Bx - Ax); |
f79b19a1 |
80 | if (myIsLimited) |
81 | { |
82 | OK = LimitAndMayBeSwap(F, Ax, Bx, FBx, Cx, FCx); |
83 | if (!OK) |
84 | return; |
85 | } |
86 | else |
87 | { |
88 | OK = F.Value(Cx, FCx); |
89 | if (!OK) |
90 | return; |
91 | } |
92 | |
7fd59977 |
93 | while(FBx > FCx) { |
94 | r = (Bx - Ax) * (FBx -FCx); |
95 | q = (Bx - Cx) * (FBx -FAx); |
96 | u = Bx - ((Bx - Cx) * q - (Bx - Ax) * r) / |
97 | (2.0 * SIGN(MAX(fabs(q - r), TINY), q - r)); |
98 | ulim = Bx + GLIMIT * (Cx - Bx); |
f79b19a1 |
99 | if (myIsLimited) |
100 | ulim = Limited(ulim); |
101 | if ((Bx - u) * (u - Cx) > 0.0) { |
102 | // u is between B and C |
7fd59977 |
103 | OK = F.Value(u, fu); |
104 | if(!OK) return; |
105 | if(fu < FCx) { |
f79b19a1 |
106 | // solution is found (B, u, c) |
7fd59977 |
107 | Ax = Bx; |
108 | Bx = u; |
109 | FAx = FBx; |
110 | FBx = fu; |
111 | Done = Standard_True; |
112 | return; |
113 | } |
114 | else if(fu > FBx) { |
f79b19a1 |
115 | // solution is found (A, B, u) |
7fd59977 |
116 | Cx = u; |
117 | FCx = fu; |
118 | Done = Standard_True; |
119 | return; |
120 | } |
f79b19a1 |
121 | // get next prob after (B, C) |
7fd59977 |
122 | u = Cx + Lambda * (Cx - Bx); |
f79b19a1 |
123 | if (myIsLimited) |
124 | { |
125 | OK = LimitAndMayBeSwap(F, Bx, Cx, FCx, u, fu); |
126 | if (!OK) |
127 | return; |
128 | } |
129 | else |
130 | { |
131 | OK = F.Value(u, fu); |
132 | if (!OK) |
133 | return; |
134 | } |
7fd59977 |
135 | } |
136 | else if((Cx - u) * (u - ulim) > 0.0) { |
f79b19a1 |
137 | // u is beyond C but between C and limit |
7fd59977 |
138 | OK = F.Value(u, fu); |
139 | if(!OK) return; |
7fd59977 |
140 | } |
141 | else if ((u - ulim) * (ulim - Cx) >= 0.0) { |
f79b19a1 |
142 | // u is beyond limit |
7fd59977 |
143 | u = ulim; |
144 | OK = F.Value(u, fu); |
145 | if(!OK) return; |
146 | } |
147 | else { |
f79b19a1 |
148 | // u tends to approach to the side of A, |
149 | // so reset it to the next prob after (B, C) |
7fd59977 |
150 | u = Cx + GOLD * (Cx - Bx); |
f79b19a1 |
151 | if (myIsLimited) |
152 | { |
153 | OK = LimitAndMayBeSwap(F, Bx, Cx, FCx, u, fu); |
154 | if (!OK) |
155 | return; |
156 | } |
157 | else |
158 | { |
159 | OK = F.Value(u, fu); |
160 | if (!OK) |
161 | return; |
162 | } |
7fd59977 |
163 | } |
164 | SHFT(Ax, Bx, Cx, u); |
165 | SHFT(FAx, FBx, FCx, fu); |
166 | } |
167 | Done = Standard_True; |
168 | } |
169 | |
170 | |
171 | |
172 | |
173 | math_BracketMinimum::math_BracketMinimum(math_Function& F, |
174 | const Standard_Real A, |
f79b19a1 |
175 | const Standard_Real B) |
176 | : Done(Standard_False), |
177 | Ax(A), Bx(B), Cx(0.), |
178 | FAx(0.), FBx(0.), FCx(0.), |
179 | myLeft(-Precision::Infinite()), |
180 | myRight(Precision::Infinite()), |
181 | myIsLimited(Standard_False), |
182 | myFA(Standard_False), |
183 | myFB (Standard_False) |
184 | { |
185 | Perform(F); |
7fd59977 |
186 | } |
187 | |
188 | math_BracketMinimum::math_BracketMinimum(math_Function& F, |
189 | const Standard_Real A, |
190 | const Standard_Real B, |
f79b19a1 |
191 | const Standard_Real FA) |
192 | : Done(Standard_False), |
193 | Ax(A), Bx(B), Cx(0.), |
194 | FAx(FA), FBx(0.), FCx(0.), |
195 | myLeft(-Precision::Infinite()), |
196 | myRight(Precision::Infinite()), |
197 | myIsLimited(Standard_False), |
198 | myFA(Standard_True), |
199 | myFB (Standard_False) |
200 | { |
201 | Perform(F); |
7fd59977 |
202 | } |
203 | |
204 | math_BracketMinimum::math_BracketMinimum(math_Function& F, |
205 | const Standard_Real A, |
206 | const Standard_Real B, |
207 | const Standard_Real FA, |
f79b19a1 |
208 | const Standard_Real FB) |
209 | : Done(Standard_False), |
210 | Ax(A), Bx(B), Cx(0.), |
211 | FAx(FA), FBx(FB), FCx(0.), |
212 | myLeft(-Precision::Infinite()), |
213 | myRight(Precision::Infinite()), |
214 | myIsLimited(Standard_False), |
215 | myFA(Standard_True), |
216 | myFB(Standard_True) |
217 | { |
218 | Perform(F); |
7fd59977 |
219 | } |
220 | |
221 | |
222 | void math_BracketMinimum::Values(Standard_Real& A, Standard_Real& B, Standard_Real& C) const{ |
223 | |
224 | StdFail_NotDone_Raise_if(!Done, " "); |
225 | A = Ax; |
226 | B = Bx; |
227 | C = Cx; |
228 | } |
229 | |
230 | void math_BracketMinimum::FunctionValues(Standard_Real& FA, Standard_Real& FB, Standard_Real& FC) const{ |
231 | |
232 | StdFail_NotDone_Raise_if(!Done, " "); |
233 | FA = FAx; |
234 | FB = FBx; |
235 | FC = FCx; |
236 | } |
237 | |
238 | void math_BracketMinimum::Dump(Standard_OStream& o) const { |
239 | |
240 | o << "math_BracketMinimum "; |
241 | if(Done) { |
242 | o << " Status = Done \n"; |
04232180 |
243 | o << " The bracketed triplet is: " << std::endl; |
244 | o << Ax << ", " << Bx << ", " << Cx << std::endl; |
245 | o << " The corresponding function values are: "<< std::endl; |
246 | o << FAx << ", " << FBx << ", " << FCx << std::endl; |
7fd59977 |
247 | } |
248 | else { |
249 | o << " Status = not Done \n"; |
250 | } |
251 | } |
252 | |