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