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 | // |
973c2be1 |
6 | // This library is free software; you can redistribute it and / or modify it |
7 | // under the terms of the GNU Lesser General Public version 2.1 as published |
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 | #include <math_BrentMinimum.ixx> |
16 | #include <math_Function.hxx> |
17 | |
18 | #define CGOLD 0.3819660 |
19 | #ifdef MAX |
20 | #undef MAX |
21 | #endif |
22 | #define MAX(a,b) ((a) > (b) ? (a) : (b)) |
23 | #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) |
24 | #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d) |
25 | |
26 | void math_BrentMinimum::Perform(math_Function& F, |
27 | const Standard_Real ax, |
28 | const Standard_Real bx, |
29 | const Standard_Real cx) { |
30 | |
31 | Standard_Boolean OK; |
32 | Standard_Real etemp, fu, p, q, r; |
33 | Standard_Real tol1, tol2, u, v, w, xm; |
34 | Standard_Real e = 0.0; |
35 | Standard_Real d = RealLast(); |
36 | |
37 | a = ((ax < cx) ? ax : cx); |
38 | b = ((ax > cx) ? ax : cx); |
39 | x = w = v = bx; |
40 | if (!myF) { |
41 | OK = F.Value(x, fx); |
42 | if(!OK) return; |
43 | } |
44 | fw = fv = fx; |
45 | for(iter = 1; iter <= Itermax; iter++) { |
46 | xm = 0.5 * (a + b); |
47 | tol1 = XTol * fabs(x) + EPSZ; |
48 | tol2 = 2.0 * tol1; |
49 | if(IsSolutionReached(F)) { |
50 | Done = Standard_True; |
51 | return; |
52 | } |
53 | if(fabs(e) > tol1) { |
54 | r = (x - w) * (fx - fv); |
55 | q = (x - v) * (fx - fw); |
56 | p = (x - v) * q - (x - w) * r; |
57 | q = 2.0 * (q - r); |
58 | if(q > 0.0) p = -p; |
59 | q = fabs(q); |
60 | etemp = e; |
61 | e = d; |
62 | if(fabs(p) >= fabs(0.5 * q * etemp) |
63 | || p <= q * ( a - x) || p >= q * (b - x)) { |
64 | e = (x >= xm ? a - x : b - x); |
65 | d = CGOLD * e; |
66 | } |
67 | else { |
68 | d = p / q; |
69 | u = x + d; |
70 | if(u - a < tol2 || b - u < tol2) d = SIGN(tol1, xm - x); |
71 | } |
72 | } |
73 | else { |
74 | e = (x >= xm ? a - x : b - x); |
75 | d = CGOLD * e; |
76 | } |
77 | u = (fabs(d) >= tol1 ? x + d : x + SIGN(tol1, d)); |
78 | OK = F.Value(u, fu); |
79 | if(!OK) return; |
80 | if(fu <= fx) { |
81 | if(u >= x) a = x; else b = x; |
82 | SHFT(v, w, x, u); |
83 | SHFT(fv, fw, fx, fu); |
84 | } |
85 | else { |
86 | if(u < x) a = u; else b = u; |
87 | if(fu <= fw || w == x) { |
88 | v = w; |
89 | w = u; |
90 | fv = fw; |
91 | fw = fu; |
92 | } |
93 | else if(fu <= fv || v == x || v == w) { |
94 | v = u; |
95 | fv = fu; |
96 | } |
97 | } |
98 | } |
99 | Done = Standard_False; |
100 | return; |
101 | } |
102 | |
103 | |
104 | math_BrentMinimum::math_BrentMinimum(math_Function& F, |
105 | const Standard_Real Ax, |
106 | const Standard_Real Bx, |
107 | const Standard_Real Cx, |
108 | const Standard_Real TolX, |
109 | const Standard_Integer NbIterations, |
110 | const Standard_Real ZEPS) { |
111 | |
112 | XTol = TolX; |
113 | EPSZ = ZEPS; |
114 | Itermax = NbIterations; |
115 | myF = Standard_False; |
116 | Perform(F, Ax, Bx, Cx); |
117 | } |
118 | |
119 | |
120 | // Constructeur d'initialisation des champs. |
121 | |
122 | math_BrentMinimum::math_BrentMinimum(const Standard_Real TolX, |
123 | const Standard_Integer NbIterations, |
124 | const Standard_Real ZEPS) { |
125 | myF = Standard_False; |
126 | XTol = TolX; |
127 | EPSZ = ZEPS; |
128 | Itermax = NbIterations; |
129 | } |
130 | |
131 | math_BrentMinimum::math_BrentMinimum(const Standard_Real TolX, |
132 | const Standard_Real Fbx, |
133 | const Standard_Integer NbIterations, |
134 | const Standard_Real ZEPS) { |
135 | |
136 | fx = Fbx; |
137 | myF = Standard_True; |
138 | XTol = TolX; |
139 | EPSZ = ZEPS; |
140 | Itermax = NbIterations; |
141 | } |
142 | |
143 | |
144 | // Standard_Boolean math_BrentMinimum::IsSolutionReached(math_Function& F) { |
145 | Standard_Boolean math_BrentMinimum::IsSolutionReached(math_Function& ) { |
146 | |
6e6cd5d9 |
147 | // Standard_Real xm = 0.5 * (a + b); |
7fd59977 |
148 | // modified by NIZHNY-MKK Mon Oct 3 17:45:57 2005.BEGIN |
149 | // Standard_Real tol = XTol * fabs(x) + EPSZ; |
150 | // return fabs(x - xm) <= 2.0 * tol - 0.5 * (b - a); |
151 | Standard_Real TwoTol = 2.0 *(XTol * fabs(x) + EPSZ); |
152 | return ((x <= (TwoTol + a)) && (x >= (b - TwoTol))); |
153 | // modified by NIZHNY-MKK Mon Oct 3 17:46:00 2005.END |
154 | } |
155 | |
156 | |
157 | |
158 | void math_BrentMinimum::Dump(Standard_OStream& o) const { |
159 | o << "math_BrentMinimum "; |
160 | if(Done) { |
161 | o << " Status = Done \n"; |
162 | o << " Location value = " << x <<"\n"; |
163 | o << " Minimum value = " << fx << "\n"; |
164 | o << " Number of iterations = " << iter <<"\n"; |
165 | } |
166 | else { |
167 | o << " Status = not Done \n"; |
168 | } |
169 | } |