b2c0f9ac413f883a29444f088b8ae9b1d79598e7
[occt.git] / src / math / math_BrentMinimum.cxx
1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
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
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.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15
16 #include <math_BrentMinimum.hxx>
17 #include <math_Function.hxx>
18 #include <StdFail_NotDone.hxx>
19
20 #define CGOLD         0.3819660
21 #ifdef MAX
22 #undef MAX
23 #endif
24 #define MAX(a,b)      ((a) > (b) ? (a) : (b))
25 #define SIGN(a,b)     ((b) > 0.0 ? fabs(a) : -fabs(a))
26 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d)
27
28 //=======================================================================
29 //function : math_BrentMinimum
30 //purpose  : Constructor
31 //=======================================================================
32 math_BrentMinimum::math_BrentMinimum(const Standard_Real    theTolX,
33                                      const Standard_Integer theNbIterations,
34                                      const Standard_Real    theZEPS)
35 : a   (0.0),
36   b   (0.0),
37   x   (0.0),
38   fx  (0.0),
39   fv  (0.0),
40   fw  (0.0),
41   XTol(theTolX),
42   EPSZ(theZEPS),
43   Done   (Standard_False),
44   iter   (0),
45   Itermax(theNbIterations),
46   myF    (Standard_False)
47 {
48 }
49
50 //=======================================================================
51 //function : math_BrentMinimum
52 //purpose  : Constructor
53 //=======================================================================
54 math_BrentMinimum::math_BrentMinimum(const Standard_Real    theTolX,
55                                      const Standard_Real    theFbx,
56                                      const Standard_Integer theNbIterations,
57                                      const Standard_Real    theZEPS)
58 : a   (0.0),
59   b   (0.0),
60   x   (0.0),
61   fx  (theFbx),
62   fv  (0.0),
63   fw  (0.0),
64   XTol(theTolX),
65   EPSZ(theZEPS),
66   Done   (Standard_False),
67   iter   (0),
68   Itermax(theNbIterations),
69   myF    (Standard_True)
70 {
71 }
72
73 //=======================================================================
74 //function : ~math_BrentMinimum
75 //purpose  : Destructor
76 //=======================================================================
77 math_BrentMinimum::~math_BrentMinimum()
78 {
79 }
80
81 //=======================================================================
82 //function : Perform
83 //purpose  : 
84 //=======================================================================
85 void math_BrentMinimum::Perform(math_Function& F,
86                                 const Standard_Real    ax,
87                                 const Standard_Real    bx,
88                                 const Standard_Real    cx)
89 {
90   Standard_Boolean OK;  
91   Standard_Real etemp, fu, p, q, r;  
92   Standard_Real tol1, tol2, u, v, w, xm;
93   Standard_Real e = 0.0;
94   Standard_Real d = RealLast();
95   
96   a = ((ax < cx) ? ax : cx);
97   b = ((ax > cx) ? ax : cx);
98   x = w = v = bx;
99   if (!myF) {
100     OK = F.Value(x, fx);
101     if(!OK) return;
102   }
103   fw = fv = fx;
104   for(iter = 1; iter <= Itermax; iter++) {
105     xm = 0.5 * (a + b);
106     tol1 = XTol * fabs(x) + EPSZ;
107     tol2 = 2.0 * tol1;
108     if(IsSolutionReached(F)) {
109       Done = Standard_True;
110       return;
111     }
112     if(fabs(e) > tol1) {
113       r = (x - w) * (fx - fv);
114       q = (x - v) * (fx - fw);
115       p = (x - v) * q - (x - w) * r;
116       q = 2.0 * (q - r);
117       if(q > 0.0) p = -p;
118       q = fabs(q);
119       etemp = e;
120       e = d;
121       if(fabs(p) >= fabs(0.5 * q * etemp) 
122          || p <= q * ( a - x) || p >= q * (b - x)) {
123         e = (x >= xm ? a - x : b - x);
124         d = CGOLD * e;
125       }
126       else {
127         d = p / q;
128         u = x + d;
129         if(u - a < tol2 || b - u < tol2) d = SIGN(tol1, xm - x);
130       }
131     }
132     else {
133       e = (x >= xm ? a - x : b - x);
134       d = CGOLD * e;
135     }
136     u = (fabs(d) >= tol1 ? x + d : x + SIGN(tol1, d));
137     OK = F.Value(u, fu);
138     if(!OK) return;
139     if(fu <= fx) {
140       if(u >= x) a = x; else b = x;
141       SHFT(v, w, x, u);
142       SHFT(fv, fw, fx, fu);
143     }
144     else {
145       if(u < x) a = u; else b = u;
146       if(fu <= fw || w == x) {
147         v = w;
148         w = u;
149         fv = fw;
150         fw = fu;
151       }
152       else if(fu <= fv || v == x || v == w) {
153         v = u;
154         fv = fu;
155       }
156     }
157   }
158   Done = Standard_False;
159   return;
160 }
161
162 //=======================================================================
163 //function : Dump
164 //purpose  : 
165 //=======================================================================
166 void math_BrentMinimum::Dump(Standard_OStream& o) const
167 {
168   o << "math_BrentMinimum ";
169   if(Done) {
170     o << " Status = Done \n";
171     o << " Location value = " << x <<"\n";
172     o << " Minimum value = " << fx << "\n";
173     o << " Number of iterations = " << iter <<"\n";
174   }
175   else {
176     o << " Status = not Done \n";
177   }
178 }