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