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(F, ax, xx, bx, 1.0e-10, 100);
  92        if(Sol.IsDone()) {
  93          Standard_Real Scale = Sol.Location();
  94          Result = Sol.Minimum();
  95          Dir.Multiply(Scale);
  96          P.Add(Dir);
  97          return Standard_True;
  98        }
  99      }
 100      return Standard_False;
 101   }
 102
 103
 104 void  math_FRPR::Perform(math_MultipleVarFunctionWithGradient& F,
 105                          const math_Vector& StartingPoint) {
 106
 107        Standard_Boolean Good;
 108        Standard_Integer n = TheLocation.Length();
 109        Standard_Integer j, its;
 110        Standard_Real gg, gam, dgg;
 111
 112        math_Vector g(1, n), h(1, n);
 113
 114        math_Vector Temp1(1, n);
 115        math_Vector Temp2(1, n);
 116        math_Vector Temp3(1, n);
 117        DirFunctionTer F_Dir(Temp1, Temp2, Temp3, F);
 118
 119        TheLocation = StartingPoint;
 120        Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
 121        if(!Good) {
 122          Done = Standard_False;
 123          TheStatus = math_FunctionError;
 124          return;
 125        }
 126
 127        g = -TheGradient;
 128        h = g;
 129        TheGradient = g;
 130
 131        for(its = 1; its <= Itermax; its++) {
 132          Iter = its;
 133
 134          Standard_Boolean IsGood = MinimizeDirection(TheLocation,
 135                                             TheGradient, TheMinimum, F_Dir);
 136          if(!IsGood) {
 137            Done = Standard_False;
 138            TheStatus = math_DirectionSearchError;
 139            return;
 140          }
 141          if(IsSolutionReached(F)) {
 142            Done = Standard_True;
 143            State = F.GetStateNumber();
 144            TheStatus = math_OK;
 145            return;
 146          }
 147          Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
 148          if(!Good) {
 149            Done = Standard_False;
 150            TheStatus = math_FunctionError;
 151            return;
 152          }
 153
 154          dgg =0.0;
 155          gg = 0.0;
 156
 157          for(j = 1; j<= n; j++) {
 158            gg += g(j)*g(j);
 159 //         dgg += TheGradient(j)*TheGradient(j);  //for Fletcher-Reeves
 160            dgg += (TheGradient(j)+g(j)) * TheGradient(j);  //for Polak-Ribiere
 161          }
 162
 163          if (gg == 0.0) {
 164            //Unlikely. If gradient is exactly 0 then we are already done.
 165            Done = Standard_False;
 166            TheStatus = math_FunctionError;
 167            return;
 168          }
 169
 170          gam = dgg/gg;
 171          g = -TheGradient;
 172          TheGradient = g + gam*h;
 173          h = TheGradient;
 174        }
 175        Done = Standard_False;
 176        TheStatus = math_TooManyIterations;
 177        return;
 178
 179      }
 180
 181
 182
 183     Standard_Boolean math_FRPR::IsSolutionReached(
 184 //                           math_MultipleVarFunctionWithGradient& F) {
 185                            math_MultipleVarFunctionWithGradient& ) {
 186
 187        return (2.0 * fabs(TheMinimum - PreviousMinimum)) <=
 188               XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ);
 189     }
 190
 191     math_FRPR::math_FRPR(math_MultipleVarFunctionWithGradient& F,
 192                          const math_Vector& StartingPoint,
 193                          const Standard_Real        Tolerance,
 194                          const Standard_Integer     NbIterations,
 195                          const Standard_Real        ZEPS)
 196                          : TheLocation(1, StartingPoint.Length()),
 197                            TheGradient(1, StartingPoint.Length()) {
 198
 199        XTol = Tolerance;
 200        EPSZ = ZEPS;
 201        Itermax = NbIterations;
 202        Perform(F, StartingPoint);
 203     }
 204
 205
 206     math_FRPR::math_FRPR(math_MultipleVarFunctionWithGradient& F,
 207                          const Standard_Real        Tolerance,
 208                          const Standard_Integer     NbIterations,
 209                          const Standard_Real        ZEPS)
 210                          : TheLocation(1, F.NbVariables()),
 211                            TheGradient(1, F.NbVariables()) {
 212
 213        XTol = Tolerance;
 214        EPSZ = ZEPS;
 215        Itermax = NbIterations;
 216     }
 217
 218
 219     math_FRPR::~math_FRPR()
 220     {
 221     }
 222
 223     void math_FRPR::Dump(Standard_OStream& o) const {
 224
 225        o << "math_FRPR ";
 226        if(Done) {
 227          o << " Status = Done \n";
 228          o << " Location Vector = "<< TheLocation << "\n";
 229          o << " Minimum value = " << TheMinimum <<"\n";
 230          o << " Number of iterations = " << Iter <<"\n";
 231        }
 232        else {
 233          o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
 234        }
 235     }
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250