src/math/math_Powell.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_Function.hxx>
  25 #include <math_Matrix.hxx>
  26 #include <math_MultipleVarFunction.hxx>
  27 #include <math_Powell.hxx>
  28 #include <Standard_DimensionError.hxx>
  29 #include <StdFail_NotDone.hxx>
  30
  31 namespace {
  32 static inline Standard_Real SQR (const Standard_Real a)
  33 {
  34     return a * a;
  35 }
  36 }
  37
  38 class DirFunctionBis : public math_Function {
  39
  40   math_Vector *P0;
  41   math_Vector *Dir;
  42   math_Vector *P;
  43   math_MultipleVarFunction *F;
  44
  45   public :
  46     DirFunctionBis(math_Vector& V1,
  47                 math_Vector& V2,
  48                 math_Vector& V3,
  49                 math_MultipleVarFunction& f);
  50
  51     void Initialize(const math_Vector& p0, const math_Vector& dir);
  52
  53     virtual Standard_Boolean Value(const Standard_Real x, Standard_Real& fval);
  54 };
  55
  56  DirFunctionBis::DirFunctionBis(math_Vector& V1,
  57                           math_Vector& V2,
  58                           math_Vector& V3,
  59                           math_MultipleVarFunction& f) {
  60    P0 = &V1;
  61    Dir = &V2;
  62    P = &V3;
  63    F = &f;
  64  }
  65
  66
  67 void DirFunctionBis::Initialize(const math_Vector& p0,
  68                              const math_Vector& dir) {
  69
  70   *P0 = p0;
  71   *Dir = dir;
  72 }
  73
  74 Standard_Boolean DirFunctionBis::Value(const Standard_Real x, Standard_Real& fval) {
  75   *P = *Dir;
  76   P->Multiply(x);
  77   P->Add(*P0);
  78   F->Value(*P, fval);
  79   return Standard_True;
  80 }
  81
  82
  83 static Standard_Boolean MinimizeDirection(math_Vector& P,
  84                                  math_Vector& Dir,
  85                                  Standard_Real& Result,
  86                                  DirFunctionBis& F) {
  87
  88   Standard_Real ax;
  89   Standard_Real xx;
  90   Standard_Real bx;
  91
  92
  93   F.Initialize(P, Dir);
  94
  95   math_BracketMinimum Bracket(F, 0.0, 1.0);
  96   if (Bracket.IsDone()) {
  97     Bracket.Values(ax, xx, bx);
  98     math_BrentMinimum Sol(1.0e-10);
  99     Sol.Perform(F, ax, xx, bx);
 100     if (Sol.IsDone()) {
 101       Standard_Real Scale = Sol.Location();
 102       Result = Sol.Minimum();
 103       Dir.Multiply(Scale);
 104       P.Add(Dir);
 105       return Standard_True;
 106     }
 107   }
 108   return Standard_False;
 109 }
 110
 111 //=======================================================================
 112 //function : math_Powell
 113 //purpose  : Constructor
 114 //=======================================================================
 115 math_Powell::math_Powell(const math_MultipleVarFunction& theFunction,
 116                          const Standard_Real             theTolerance,
 117                          const Standard_Integer          theNbIterations,
 118                          const Standard_Real             theZEPS)
 119 : TheLocation     (1, theFunction.NbVariables()),
 120   TheMinimum      (RealLast()),
 121   TheLocationError(RealLast()),
 122   PreviousMinimum (RealLast()),
 123   XTol            (theTolerance),
 124   EPSZ            (theZEPS),
 125   Done            (Standard_False),
 126   Iter            (0),
 127   TheStatus       (math_NotBracketed),
 128   TheDirections   (1, theFunction.NbVariables(), 1, theFunction.NbVariables()),
 129   State           (0),
 130   Itermax         (theNbIterations)
 131 {
 132 }
 133
 134 //=======================================================================
 135 //function : ~math_Powell
 136 //purpose  : Destructor
 137 //=======================================================================
 138 math_Powell::~math_Powell()
 139 {
 140 }
 141
 142 //=======================================================================
 143 //function : Perform
 144 //purpose  :
 145 //=======================================================================
 146 void math_Powell::Perform(math_MultipleVarFunction& F,
 147                           const math_Vector& StartingPoint,
 148                           const math_Matrix& StartingDirections)
 149 {
 150   Done = Standard_False;
 151   Standard_Integer i, ibig, j;
 152   Standard_Real t, fptt, del;
 153   Standard_Integer n = TheLocation.Length();
 154   math_Vector pt(1,n);
 155   math_Vector ptt(1,n);
 156   math_Vector xit(1,n);
 157   math_Vector Temp1(1, n);
 158   math_Vector Temp2(1, n);
 159   math_Vector Temp3(1, n);
 160   DirFunctionBis F_Dir(Temp1, Temp2, Temp3, F);
 161
 162   TheLocation = StartingPoint;
 163   TheDirections = StartingDirections;
 164   pt = TheLocation;  //sauvegarde du point initial
 165
 166
 167   for(Iter = 1; Iter<= Itermax; Iter++) {
 168     F.Value(TheLocation, PreviousMinimum);
 169     ibig = 0;
 170     del = 0.0;
 171     for (i = 1; i <= n; i++){
 172       for(j =1; j<= n; j++) xit(j) = TheDirections(j,i);
 173       F.Value(TheLocation, fptt);
 174       Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit,
 175                                          TheMinimum, F_Dir);
 176
 177       if (!IsGood) {
 178         Done = Standard_False;
 179         TheStatus = math_DirectionSearchError;
 180         return;
 181       }
 182
 183       if (fabs(fptt - TheMinimum)> del) {
 184         del = fabs(fptt- TheMinimum);
 185         ibig = i;
 186       }
 187     }
 188
 189     if (IsSolutionReached(F)) {
 190       //Termination criterion
 191       State = F.GetStateNumber();
 192       Done = Standard_True;
 193       TheStatus = math_OK;
 194       return;
 195     }
 196
 197     if (Iter == Itermax) {
 198       Done = Standard_False;
 199       TheStatus = math_TooManyIterations;
 200       return;
 201     }
 202
 203       ptt = 2.0 * TheLocation - pt;
 204       xit = TheLocation - pt;
 205       pt = TheLocation;
 206
 207     // Valeur de la fonction au point extrapole:
 208
 209     F.Value(ptt, fptt);
 210
 211     if (fptt < PreviousMinimum) {
 212       t = 2.0 *(PreviousMinimum -2.0*TheMinimum +fptt)*
 213         SQR(PreviousMinimum-TheMinimum -del)-del*
 214           SQR(PreviousMinimum-fptt);
 215       if (t <0.0) {
 216         //Minimisation along the direction
 217         Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit,
 218                                            TheMinimum, F_Dir);
 219         if(!IsGood) {
 220           Done = Standard_False;
 221           TheStatus = math_FunctionError;
 222           return;
 223         }
 224
 225         for(j =1; j <= n; j++) {
 226           TheDirections(j, ibig)=xit(j);
 227         }
 228       }
 229     }
 230   }
 231 }
 232
 233 //=======================================================================
 234 //function : Dump
 235 //purpose  :
 236 //=======================================================================
 237 void math_Powell::Dump(Standard_OStream& o) const
 238 {
 239   o << "math_Powell resolution:";
 240   if(Done) {
 241     o << " Status = Done \n";
 242     o << " Location Vector = "<< TheLocation << "\n";
 243     o << " Minimum value = " << TheMinimum <<"\n";
 244     o << " Number of iterations = " << Iter <<"\n";
 245   }
 246   else {
 247     o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
 248   }
 249 }