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