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