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