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 | |
42cf5bc1 |
15 | |
16 | #include <math_BrentMinimum.hxx> |
7fd59977 |
17 | #include <math_Function.hxx> |
42cf5bc1 |
18 | #include <StdFail_NotDone.hxx> |
7fd59977 |
19 | |
f542b7bb |
20 | static const Standard_Real CGOLD = 0.3819660; //0.5*(3 - sqrt(5)); |
21 | |
22 | |
23 | //======================================================================= |
24 | //function : SHFT |
25 | //purpose : Shifts arguments |
26 | //======================================================================= |
27 | inline void SHFT(Standard_Real &theA, Standard_Real &theB, |
28 | Standard_Real &theC, Standard_Real &theD) |
29 | { |
30 | theA = theB; |
31 | theB = theC; |
32 | theC = theD; |
33 | } |
7fd59977 |
34 | |
859a47c3 |
35 | //======================================================================= |
36 | //function : math_BrentMinimum |
37 | //purpose : Constructor |
38 | //======================================================================= |
39 | math_BrentMinimum::math_BrentMinimum(const Standard_Real theTolX, |
40 | const Standard_Integer theNbIterations, |
41 | const Standard_Real theZEPS) |
42 | : a (0.0), |
43 | b (0.0), |
44 | x (0.0), |
45 | fx (0.0), |
46 | fv (0.0), |
47 | fw (0.0), |
48 | XTol(theTolX), |
49 | EPSZ(theZEPS), |
50 | Done (Standard_False), |
51 | iter (0), |
52 | Itermax(theNbIterations), |
53 | myF (Standard_False) |
54 | { |
55 | } |
56 | |
57 | //======================================================================= |
58 | //function : math_BrentMinimum |
59 | //purpose : Constructor |
60 | //======================================================================= |
61 | math_BrentMinimum::math_BrentMinimum(const Standard_Real theTolX, |
62 | const Standard_Real theFbx, |
63 | const Standard_Integer theNbIterations, |
64 | const Standard_Real theZEPS) |
65 | : a (0.0), |
66 | b (0.0), |
67 | x (0.0), |
68 | fx (theFbx), |
69 | fv (0.0), |
70 | fw (0.0), |
71 | XTol(theTolX), |
72 | EPSZ(theZEPS), |
73 | Done (Standard_False), |
74 | iter (0), |
75 | Itermax(theNbIterations), |
76 | myF (Standard_True) |
77 | { |
78 | } |
79 | |
80 | //======================================================================= |
81 | //function : ~math_BrentMinimum |
82 | //purpose : Destructor |
83 | //======================================================================= |
6da30ff1 |
84 | math_BrentMinimum::~math_BrentMinimum() |
85 | { |
86 | } |
87 | |
859a47c3 |
88 | //======================================================================= |
89 | //function : Perform |
90 | //purpose : |
91 | //======================================================================= |
7fd59977 |
92 | void math_BrentMinimum::Perform(math_Function& F, |
859a47c3 |
93 | const Standard_Real ax, |
94 | const Standard_Real bx, |
95 | const Standard_Real cx) |
96 | { |
7fd59977 |
97 | Standard_Boolean OK; |
98 | Standard_Real etemp, fu, p, q, r; |
99 | Standard_Real tol1, tol2, u, v, w, xm; |
100 | Standard_Real e = 0.0; |
101 | Standard_Real d = RealLast(); |
f542b7bb |
102 | |
7fd59977 |
103 | a = ((ax < cx) ? ax : cx); |
104 | b = ((ax > cx) ? ax : cx); |
105 | x = w = v = bx; |
106 | if (!myF) { |
107 | OK = F.Value(x, fx); |
f542b7bb |
108 | if (!OK) return; |
7fd59977 |
109 | } |
110 | fw = fv = fx; |
f542b7bb |
111 | for (iter = 1; iter <= Itermax; iter++) { |
7fd59977 |
112 | xm = 0.5 * (a + b); |
113 | tol1 = XTol * fabs(x) + EPSZ; |
114 | tol2 = 2.0 * tol1; |
f542b7bb |
115 | if (IsSolutionReached(F)) { |
7fd59977 |
116 | Done = Standard_True; |
117 | return; |
118 | } |
f542b7bb |
119 | if (fabs(e) > tol1) { |
7fd59977 |
120 | r = (x - w) * (fx - fv); |
121 | q = (x - v) * (fx - fw); |
122 | p = (x - v) * q - (x - w) * r; |
123 | q = 2.0 * (q - r); |
f542b7bb |
124 | if (q > 0.0) p = -p; |
7fd59977 |
125 | q = fabs(q); |
126 | etemp = e; |
127 | e = d; |
f542b7bb |
128 | if (fabs(p) >= fabs(0.5 * q * etemp) |
129 | || p <= q * (a - x) || p >= q * (b - x)) { |
130 | e = (x >= xm ? a - x : b - x); |
131 | d = CGOLD * e; |
7fd59977 |
132 | } |
133 | else { |
f542b7bb |
134 | d = p / q; |
135 | u = x + d; |
136 | if (u - a < tol2 || b - u < tol2) d = Sign(tol1, xm - x); |
7fd59977 |
137 | } |
138 | } |
139 | else { |
140 | e = (x >= xm ? a - x : b - x); |
141 | d = CGOLD * e; |
142 | } |
f542b7bb |
143 | u = (fabs(d) >= tol1 ? x + d : x + Sign(tol1, d)); |
7fd59977 |
144 | OK = F.Value(u, fu); |
f542b7bb |
145 | if (!OK) return; |
146 | if (fu <= fx) { |
147 | if (u >= x) a = x; else b = x; |
7fd59977 |
148 | SHFT(v, w, x, u); |
149 | SHFT(fv, fw, fx, fu); |
150 | } |
151 | else { |
f542b7bb |
152 | if (u < x) a = u; else b = u; |
153 | if (fu <= fw || w == x) { |
154 | v = w; |
155 | w = u; |
156 | fv = fw; |
157 | fw = fu; |
7fd59977 |
158 | } |
f542b7bb |
159 | else if (fu <= fv || v == x || v == w) { |
160 | v = u; |
161 | fv = fu; |
7fd59977 |
162 | } |
163 | } |
164 | } |
165 | Done = Standard_False; |
166 | return; |
167 | } |
168 | |
859a47c3 |
169 | //======================================================================= |
170 | //function : Dump |
171 | //purpose : |
172 | //======================================================================= |
173 | void math_BrentMinimum::Dump(Standard_OStream& o) const |
174 | { |
175 | o << "math_BrentMinimum "; |
176 | if(Done) { |
177 | o << " Status = Done \n"; |
178 | o << " Location value = " << x <<"\n"; |
179 | o << " Minimum value = " << fx << "\n"; |
180 | o << " Number of iterations = " << iter <<"\n"; |
7fd59977 |
181 | } |
859a47c3 |
182 | else { |
183 | o << " Status = not Done \n"; |
184 | } |
185 | } |