src/math/math_FRPR.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 //#ifndef OCCT_DEBUG
  16 #define No_Standard_RangeError
  17 #define No_Standard_OutOfRange
  18 #define No_Standard_DimensionError
  19
  20 //#endif
  21
  22 #include <math_BracketMinimum.hxx>
  23 #include <math_BrentMinimum.hxx>
  24 #include <math_FRPR.hxx>
  25 #include <math_Function.hxx>
  26 #include <math_MultipleVarFunction.hxx>
  27 #include <math_MultipleVarFunctionWithGradient.hxx>
  28 #include <Standard_DimensionError.hxx>
  29 #include <StdFail_NotDone.hxx>
  30
  31 // l'utilisation de math_BrentMinumim pur trouver un minimum dans une direction
  32 // donnee n'est pas du tout optimale. voir peut etre interpolation cubique
  33 // classique et aussi essayer "recherche unidimensionnelle economique"
  34 // PROGRAMMATION MATHEMATIQUE (theorie et algorithmes) tome1 page 82.
  35 class DirFunctionTer : public math_Function {
  36
  37      math_Vector *P0;
  38      math_Vector *Dir;
  39      math_Vector *P;
  40      math_MultipleVarFunction *F;
  41
  42 public :
  43
  44      DirFunctionTer(math_Vector& V1,
  45                  math_Vector& V2,
  46                  math_Vector& V3,
  47                  math_MultipleVarFunction& f);
  48
  49      void Initialize(const math_Vector& p0, const math_Vector& dir);
  50
  51      virtual Standard_Boolean Value(const Standard_Real x, Standard_Real& fval);
  52 };
  53
  54      DirFunctionTer::DirFunctionTer(math_Vector& V1,
  55                               math_Vector& V2,
  56                               math_Vector& V3,
  57                               math_MultipleVarFunction& f) {
  58
  59         P0  = &V1;
  60         Dir = &V2;
  61         P   = &V3;
  62         F   = &f;
  63      }
  64
  65      void DirFunctionTer::Initialize(const math_Vector& p0,
  66                                   const math_Vector& dir) {
  67
  68         *P0 = p0;
  69         *Dir = dir;
  70      }
  71
  72      Standard_Boolean DirFunctionTer::Value(const Standard_Real x, Standard_Real& fval) {
  73
  74         *P = *Dir;
  75         P->Multiply(x);
  76         P->Add(*P0);
  77         fval = 0.;
  78         return F->Value(*P, fval);
  79      }
  80
  81 static Standard_Boolean MinimizeDirection(math_Vector& P,
  82                                  math_Vector& Dir,
  83                                  Standard_Real& Result,
  84                                  DirFunctionTer& F) {
  85
  86      Standard_Real ax, xx, bx;
  87
  88      F.Initialize(P, Dir);
  89      math_BracketMinimum Bracket(F, 0.0, 1.0);
  90      if(Bracket.IsDone()) {
  91        Bracket.Values(ax, xx, bx);
  92        math_BrentMinimum Sol(1.e-10);
  93        Sol.Perform(F, ax, xx, bx);
  94        if (Sol.IsDone()) {
  95          Standard_Real Scale = Sol.Location();
  96          Result = Sol.Minimum();
  97          Dir.Multiply(Scale);
  98          P.Add(Dir);
  99          return Standard_True;
 100        }
 101      }
 102      return Standard_False;
 103   }
 104
 105 //=======================================================================
 106 //function : math_FRPR
 107 //purpose  : Constructor
 108 //=======================================================================
 109 math_FRPR::math_FRPR(const math_MultipleVarFunctionWithGradient& theFunction,
 110                      const Standard_Real                         theTolerance,
 111                      const Standard_Integer                      theNbIterations,
 112                      const Standard_Real                         theZEPS)
 113
 114 : TheLocation(1, theFunction.NbVariables()),
 115   TheGradient(1, theFunction.NbVariables()),
 116   TheMinimum     (0.0),
 117   PreviousMinimum(0.0),
 118   XTol           (theTolerance),
 119   EPSZ           (theZEPS),
 120   Done           (Standard_False),
 121   Iter           (0),
 122   State          (0),
 123   TheStatus      (math_NotBracketed),
 124   Itermax        (theNbIterations)
 125 {
 126 }
 127
 128 //=======================================================================
 129 //function : ~math_FRPR
 130 //purpose  : Destructor
 131 //=======================================================================
 132 math_FRPR::~math_FRPR()
 133 {
 134 }
 135
 136
 137 //=======================================================================
 138 //function : Perform
 139 //purpose  :
 140 //=======================================================================
 141 void  math_FRPR::Perform(math_MultipleVarFunctionWithGradient& F,
 142                          const math_Vector& StartingPoint)
 143 {
 144        Standard_Boolean Good;
 145        Standard_Integer n = TheLocation.Length();
 146        Standard_Integer j, its;
 147        Standard_Real gg, gam, dgg;
 148
 149        math_Vector g(1, n), h(1, n);
 150
 151        math_Vector Temp1(1, n);
 152        math_Vector Temp2(1, n);
 153        math_Vector Temp3(1, n);
 154        DirFunctionTer F_Dir(Temp1, Temp2, Temp3, F);
 155
 156        TheLocation = StartingPoint;
 157        Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
 158        if(!Good) {
 159          Done = Standard_False;
 160          TheStatus = math_FunctionError;
 161          return;
 162        }
 163
 164        g = -TheGradient;
 165        h = g;
 166        TheGradient = g;
 167
 168        for(its = 1; its <= Itermax; its++) {
 169          Iter = its;
 170
 171          Standard_Boolean IsGood = MinimizeDirection(TheLocation,
 172                                             TheGradient, TheMinimum, F_Dir);
 173          if(!IsGood) {
 174            Done = Standard_False;
 175            TheStatus = math_DirectionSearchError;
 176            return;
 177          }
 178          if(IsSolutionReached(F)) {
 179            Done = Standard_True;
 180            State = F.GetStateNumber();
 181            TheStatus = math_OK;
 182            return;
 183          }
 184          Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
 185          if(!Good) {
 186            Done = Standard_False;
 187            TheStatus = math_FunctionError;
 188            return;
 189          }
 190
 191          dgg =0.0;
 192          gg = 0.0;
 193
 194          for(j = 1; j<= n; j++) {
 195            gg += g(j)*g(j);
 196 //         dgg += TheGradient(j)*TheGradient(j);  //for Fletcher-Reeves
 197            dgg += (TheGradient(j)+g(j)) * TheGradient(j);  //for Polak-Ribiere
 198          }
 199
 200          if (gg == 0.0) {
 201            //Unlikely. If gradient is exactly 0 then we are already done.
 202            Done = Standard_False;
 203            TheStatus = math_FunctionError;
 204            return;
 205          }
 206
 207          gam = dgg/gg;
 208          g = -TheGradient;
 209          TheGradient = g + gam*h;
 210          h = TheGradient;
 211        }
 212        Done = Standard_False;
 213        TheStatus = math_TooManyIterations;
 214        return;
 215 }
 216
 217 //=======================================================================
 218 //function : Dump
 219 //purpose  :
 220 //=======================================================================
 221 void math_FRPR::Dump(Standard_OStream& o) const
 222 {
 223   o << "math_FRPR ";
 224   if(Done) {
 225     o << " Status = Done \n";
 226     o << " Location Vector = "<< TheLocation << "\n";
 227     o << " Minimum value = " << TheMinimum <<"\n";
 228     o << " Number of iterations = " << Iter <<"\n";
 229   }
 230   else {
 231     o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
 232   }
 233 }
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248