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 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 }
  34
  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 //=======================================================================
  84 math_BrentMinimum::~math_BrentMinimum()
  85 {
  86 }
  87
  88 //=======================================================================
  89 //function : Perform
  90 //purpose  :
  91 //=======================================================================
  92 void math_BrentMinimum::Perform(math_Function& F,
  93                                 const Standard_Real    ax,
  94                                 const Standard_Real    bx,
  95                                 const Standard_Real    cx)
  96 {
  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();
 102
 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);
 108     if (!OK) return;
 109   }
 110   fw = fv = fx;
 111   for (iter = 1; iter <= Itermax; iter++) {
 112     xm = 0.5 * (a + b);
 113     tol1 = XTol * fabs(x) + EPSZ;
 114     tol2 = 2.0 * tol1;
 115     if (IsSolutionReached(F)) {
 116       Done = Standard_True;
 117       return;
 118     }
 119     if (fabs(e) > tol1) {
 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);
 124       if (q > 0.0) p = -p;
 125       q = fabs(q);
 126       etemp = e;
 127       e = d;
 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;
 132       }
 133       else {
 134         d = p / q;
 135         u = x + d;
 136         if (u - a < tol2 || b - u < tol2) d = Sign(tol1, xm - x);
 137       }
 138     }
 139     else {
 140       e = (x >= xm ? a - x : b - x);
 141       d = CGOLD * e;
 142     }
 143     u = (fabs(d) >= tol1 ? x + d : x + Sign(tol1, d));
 144     OK = F.Value(u, fu);
 145     if (!OK) return;
 146     if (fu <= fx) {
 147       if (u >= x) a = x; else b = x;
 148       SHFT(v, w, x, u);
 149       SHFT(fv, fw, fx, fu);
 150     }
 151     else {
 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;
 158       }
 159       else if (fu <= fv || v == x || v == w) {
 160         v = u;
 161         fv = fu;
 162       }
 163     }
 164   }
 165   Done = Standard_False;
 166   return;
 167 }
 168
 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";
 181   }
 182   else {
 183     o << " Status = not Done \n";
 184   }
 185 }