src/math/math_BFGS.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_BFGS.hxx>
  23 #include <math_BracketMinimum.hxx>
  24 #include <math_BrentMinimum.hxx>
  25 #include <math_FunctionWithDerivative.hxx>
  26 #include <math_Matrix.hxx>
  27 #include <math_MultipleVarFunctionWithGradient.hxx>
  28 #include <Standard_DimensionError.hxx>
  29 #include <StdFail_NotDone.hxx>
  30
  31 #define R 0.61803399
  32 #define C (1.0-R)
  33 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
  34 #define SIGN(a, b) ((b) > 0.0 ? fabs(a) : -fabs(a))
  35 #define MOV3(a,b,c, d, e, f) (a)=(d); (b)= (e); (c)=(f);
  36
  37
  38 // l'utilisation de math_BrentMinumim pur trouver un minimum dans une direction
  39 // donnee n'est pas du tout optimale. voir peut etre interpolation cubique
  40 // classique et aussi essayer "recherche unidimensionnelle economique"
  41 // PROGRAMMATION MATHEMATIQUE (theorie et algorithmes) tome1 page 82.
  42
  43 class DirFunction : public math_FunctionWithDerivative {
  44
  45      math_Vector *P0;
  46      math_Vector *Dir;
  47      math_Vector *P;
  48      math_Vector *G;
  49      math_MultipleVarFunctionWithGradient *F;
  50
  51 public :
  52
  53      DirFunction(math_Vector& V1,
  54                  math_Vector& V2,
  55                  math_Vector& V3,
  56                  math_Vector& V4,
  57                  math_MultipleVarFunctionWithGradient& f) ;
  58
  59      void Initialize(const math_Vector& p0, const math_Vector& dir) const;
  60      void TheGradient(math_Vector& Grad);
  61      virtual Standard_Boolean Value(const Standard_Real x, Standard_Real& fval) ;
  62      virtual Standard_Boolean Values(const Standard_Real x, Standard_Real& fval, Standard_Real& D) ;
  63      virtual Standard_Boolean Derivative(const Standard_Real x, Standard_Real& D) ;
  64
  65
  66 };
  67
  68      DirFunction::DirFunction(math_Vector& V1,
  69                               math_Vector& V2,
  70                               math_Vector& V3,
  71                               math_Vector& V4,
  72                               math_MultipleVarFunctionWithGradient& f) {
  73
  74         P0  = &V1;
  75         Dir = &V2;
  76         P   = &V3;
  77         F   = &f;
  78         G =   &V4;
  79      }
  80
  81      void DirFunction::Initialize(const math_Vector& p0,
  82                                   const math_Vector& dir)  const{
  83
  84         *P0 = p0;
  85         *Dir = dir;
  86      }
  87
  88      void DirFunction::TheGradient(math_Vector& Grad) {
  89        Grad = *G;
  90      }
  91
  92
  93      Standard_Boolean DirFunction::Value(const Standard_Real x,
  94                                          Standard_Real& fval)
  95      {
  96         *P = *Dir;
  97         P->Multiply(x);
  98         P->Add(*P0);
  99         F->Value(*P, fval);
 100         return Standard_True;
 101      }
 102
 103      Standard_Boolean DirFunction::Values(const Standard_Real x,
 104                                           Standard_Real& fval,
 105                                           Standard_Real& D)
 106      {
 107         *P = *Dir;
 108         P->Multiply(x);
 109         P->Add(*P0);
 110         F->Values(*P, fval, *G);
 111         D = (*G).Multiplied(*Dir);
 112         return Standard_True;
 113      }
 114      Standard_Boolean DirFunction::Derivative(const Standard_Real x,
 115                                               Standard_Real& D)
 116      {
 117         *P = *Dir;
 118         P->Multiply(x);
 119         P->Add(*P0);
 120         Standard_Real fval;
 121         F->Values(*P, fval, *G);
 122         D = (*G).Multiplied(*Dir);
 123         return Standard_True;
 124      }
 125
 126
 127 static Standard_Boolean MinimizeDirection(math_Vector&   P,
 128                                           Standard_Real  F0,
 129                                           math_Vector&   Gr,
 130                                           math_Vector&   Dir,
 131                                           Standard_Real& Result,
 132                                           DirFunction&   F) {
 133
 134
 135
 136      Standard_Real ax, xx, bx, Fax, Fxx, Fbx, F1;
 137      F.Initialize(P, Dir);
 138
 139      Standard_Real dy1, Hnr1, lambda, alfa=0;
 140      dy1 = Gr*Dir;
 141      if (dy1 != 0) {
 142        Hnr1 = Dir.Norm2();
 143        alfa = 0.7*(-F0)/dy1;
 144        lambda = 0.015/Sqrt(Hnr1);
 145      }
 146      else lambda = 1.0;
 147      if (lambda > alfa) lambda = alfa;
 148      F.Value(lambda, F1);
 149      math_BracketMinimum Bracket(F, 0.0, lambda, F0, F1);
 150      if(Bracket.IsDone()) {
 151        Bracket.Values(ax, xx, bx);
 152        Bracket.FunctionValues(Fax, Fxx, Fbx);
 153
 154        Standard_Integer niter = 100;
 155        Standard_Real tol = 1.e-03;
 156        math_BrentMinimum Sol(tol, Fxx, niter, 1.e-08);
 157        Sol.Perform(F, ax, xx, bx);
 158        if(Sol.IsDone()) {
 159          Standard_Real Scale = Sol.Location();
 160          Result = Sol.Minimum();
 161          Dir.Multiply(Scale);
 162          P.Add(Dir);
 163          return Standard_True;
 164        }
 165      }
 166      return Standard_False;
 167    }
 168
 169
 170 void  math_BFGS::Perform(math_MultipleVarFunctionWithGradient& F,
 171                          const math_Vector& StartingPoint) {
 172
 173        Standard_Boolean Good;
 174        Standard_Integer n = TheLocation.Length();
 175        Standard_Integer j, i;
 176        Standard_Real fae, fad, fac;
 177
 178        math_Vector xi(1, n), dg(1, n), hdg(1, n);
 179        math_Matrix hessin(1, n, 1, n);
 180        hessin.Init(0.0);
 181
 182        math_Vector Temp1(1, n);
 183        math_Vector Temp2(1, n);
 184        math_Vector Temp3(1, n);
 185        math_Vector Temp4(1, n);
 186        DirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
 187
 188        TheLocation = StartingPoint;
 189        Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
 190        if(!Good) {
 191          Done = Standard_False;
 192          TheStatus = math_FunctionError;
 193          return;
 194        }
 195        for(i = 1; i <= n; i++) {
 196          hessin(i, i) = 1.0;
 197          xi(i) = -TheGradient(i);
 198        }
 199
 200
 201        for(nbiter = 1; nbiter <= Itermax; nbiter++) {
 202          TheMinimum = PreviousMinimum;
 203          Standard_Boolean IsGood = MinimizeDirection(TheLocation, TheMinimum,
 204                                                      TheGradient,
 205                                                      xi, TheMinimum, F_Dir);
 206          if(!IsGood) {
 207            Done = Standard_False;
 208            TheStatus = math_DirectionSearchError;
 209            return;
 210          }
 211          if(IsSolutionReached(F)) {
 212            Done = Standard_True;
 213            TheStatus = math_OK;
 214            return;
 215          }
 216          if (nbiter == Itermax) {
 217            Done = Standard_False;
 218            TheStatus = math_TooManyIterations;
 219            return;
 220          }
 221          PreviousMinimum = TheMinimum;
 222
 223          dg = TheGradient;
 224
 225          Good = F.Values(TheLocation, TheMinimum, TheGradient);
 226          if(!Good) {
 227            Done = Standard_False;
 228            TheStatus = math_FunctionError;
 229            return;
 230          }
 231
 232          for(i = 1; i <= n; i++) {
 233            dg(i) = TheGradient(i) - dg(i);
 234          }
 235          for(i = 1; i <= n; i++) {
 236            hdg(i) = 0.0;
 237            for (j = 1; j <= n; j++)
 238              hdg(i) += hessin(i, j) * dg(j);
 239          }
 240          fac = fae = 0.0;
 241          for(i = 1; i <= n; i++) {
 242            fac += dg(i) * xi(i);
 243            fae += dg(i) * hdg(i);
 244          }
 245          fac = 1.0 / fac;
 246          fad = 1.0 / fae;
 247
 248          for(i = 1; i <= n; i++)
 249            dg(i) = fac * xi(i) - fad * hdg(i);
 250
 251          for(i = 1; i <= n; i++)
 252            for(j = 1; j <= n; j++)
 253              hessin(i, j) += fac * xi(i) * xi(j)
 254                - fad * hdg(i) * hdg(j) + fae * dg(i) * dg(j);
 255
 256          for(i = 1; i <= n; i++) {
 257            xi(i) = 0.0;
 258            for (j = 1; j <= n; j++) xi(i) -= hessin(i, j) * TheGradient(j);
 259          }
 260        }
 261        Done = Standard_False;
 262        TheStatus = math_TooManyIterations;
 263        return;
 264    }
 265
 266     Standard_Boolean math_BFGS::IsSolutionReached(
 267                            math_MultipleVarFunctionWithGradient&) const {
 268
 269        return 2.0 * fabs(TheMinimum - PreviousMinimum) <=
 270               XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ);
 271     }
 272
 273     math_BFGS::math_BFGS(const Standard_Integer     NbVariables,
 274                          const Standard_Real        Tolerance,
 275                          const Standard_Integer     NbIterations,
 276                          const Standard_Real        ZEPS)
 277                          : TheLocation(1, NbVariables),
 278                            TheGradient(1, NbVariables) {
 279
 280        XTol = Tolerance;
 281        EPSZ = ZEPS;
 282        Itermax = NbIterations;
 283     }
 284
 285
 286     math_BFGS::~math_BFGS()
 287     {
 288     }
 289
 290     void math_BFGS::Dump(Standard_OStream& o) const {
 291
 292        o<< "math_BFGS resolution: ";
 293        if(Done) {
 294          o << " Status = Done \n";
 295          o <<" Location Vector = " << Location() << "\n";
 296          o <<" Minimum value = "<< Minimum()<<"\n";
 297          o <<" Number of iterations = "<<NbIterations() <<"\n";;
 298        }
 299        else {
 300          o<< " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
 301        }
 302     }
 303
 304
 305