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 <Precision.hxx>
  29
  30 // l'utilisation de math_BrentMinumim pur trouver un minimum dans une direction
  31 // donnee n'est pas du tout optimale. voir peut etre interpolation cubique
  32 // classique et aussi essayer "recherche unidimensionnelle economique"
  33 // PROGRAMMATION MATHEMATIQUE (theorie et algorithmes) tome1 page 82.
  34
  35 // Target function for 1D problem, point and direction are known.
  36 class DirFunction : public math_FunctionWithDerivative
  37 {
  38
  39   math_Vector *P0;
  40   math_Vector *Dir;
  41   math_Vector *P;
  42   math_Vector *G;
  43   math_MultipleVarFunctionWithGradient *F;
  44
  45 public:
  46
  47   //! Ctor.
  48   DirFunction(math_Vector& V1,
  49               math_Vector& V2,
  50               math_Vector& V3,
  51               math_Vector& V4,
  52               math_MultipleVarFunctionWithGradient& f)
  53   : P0(&V1),
  54     Dir(&V2),
  55     P(&V3),
  56     G(&V4),
  57     F(&f)
  58   {}
  59
  60   //! Sets point and direction.
  61   void Initialize(const math_Vector& p0,
  62                   const math_Vector& dir) const
  63   {
  64     *P0 = p0;
  65     *Dir = dir;
  66   }
  67
  68   void TheGradient(math_Vector& Grad)
  69   {
  70     Grad = *G;
  71   }
  72
  73   virtual Standard_Boolean Value(const Standard_Real x,
  74                                  Standard_Real&      fval)
  75   {
  76     *P = *Dir;
  77     P->Multiply(x);
  78     P->Add(*P0);
  79     fval = 0.;
  80     return F->Value(*P, fval);
  81   }
  82
  83   virtual Standard_Boolean Values(const Standard_Real x,
  84                                   Standard_Real&      fval,
  85                                   Standard_Real&      D)
  86   {
  87     *P = *Dir;
  88     P->Multiply(x);
  89     P->Add(*P0);
  90     fval = D = 0.;
  91     if (F->Values(*P, fval, *G))
  92     {
  93       D = (*G).Multiplied(*Dir);
  94       return Standard_True;
  95     }
  96
  97     return Standard_False;
  98   }
  99   virtual Standard_Boolean Derivative(const Standard_Real x,
 100                                       Standard_Real&      D)
 101   {
 102     *P = *Dir;
 103     P->Multiply(x);
 104     P->Add(*P0);
 105     Standard_Real fval;
 106     D = 0.;
 107     if (F->Values(*P, fval, *G))
 108     {
 109       D = (*G).Multiplied(*Dir);
 110       return Standard_True;
 111     }
 112
 113     return Standard_False;
 114   }
 115
 116
 117 };
 118
 119 //=============================================================================
 120 //function : ComputeInitScale
 121 //purpose  : Compute the appropriate initial value of scale factor to apply
 122 //           to the direction to approach to the minimum of the function
 123 //=============================================================================
 124 static Standard_Boolean ComputeInitScale(const Standard_Real theF0,
 125                                          const math_Vector&  theDir,
 126                                          const math_Vector&  theGr,
 127                                          Standard_Real&      theScale)
 128 {
 129   const Standard_Real dy1 = theGr * theDir;
 130   if (Abs(dy1) < RealSmall())
 131     return Standard_False;
 132
 133   const Standard_Real aHnr1 = theDir.Norm2();
 134   const Standard_Real alfa = 0.7*(-theF0) / dy1;
 135   theScale = 0.015 / Sqrt(aHnr1);
 136   if (theScale > alfa)
 137     theScale = alfa;
 138
 139   return Standard_True;
 140 }
 141
 142 //=============================================================================
 143 //function : ComputeMinMaxScale
 144 //purpose  : For a given point and direction, and bounding box,
 145 //           find min and max scale factors with which the point reaches borders
 146 //           if we apply translation Point+Dir*Scale.
 147 //return   : True if found, False if point is out of bounds.
 148 //=============================================================================
 149 static Standard_Boolean ComputeMinMaxScale(const math_Vector& thePoint,
 150                                            const math_Vector& theDir,
 151                                            const math_Vector& theLeft,
 152                                            const math_Vector& theRight,
 153                                            Standard_Real&     theMinScale,
 154                                            Standard_Real&     theMaxScale)
 155 {
 156   for (Standard_Integer anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
 157   {
 158     const Standard_Real aLeft = theLeft(anIdx) - thePoint(anIdx);
 159     const Standard_Real aRight = theRight(anIdx) - thePoint(anIdx);
 160     if (Abs(theDir(anIdx)) > RealSmall())
 161     {
 162       // Use PConfusion to get off a little from the bounds to prevent
 163       // possible refuse in Value function.
 164       const Standard_Real aLScale = (aLeft + Precision::PConfusion()) / theDir(anIdx);
 165       const Standard_Real aRScale = (aRight - Precision::PConfusion()) / theDir(anIdx);
 166       if (Abs(aLeft) < Precision::PConfusion())
 167       {
 168         // Point is on the left border.
 169         theMaxScale = Min(theMaxScale, Max(0., aRScale));
 170         theMinScale = Max(theMinScale, Min(0., aRScale));
 171       }
 172       else if (Abs(aRight) < Precision::PConfusion())
 173       {
 174         // Point is on the right border.
 175         theMaxScale = Min(theMaxScale, Max(0., aLScale));
 176         theMinScale = Max(theMinScale, Min(0., aLScale));
 177       }
 178       else if (aLeft * aRight < 0)
 179       {
 180         // Point is inside allowed range.
 181         theMaxScale = Min(theMaxScale, Max(aLScale, aRScale));
 182         theMinScale = Max(theMinScale, Min(aLScale, aRScale));
 183       }
 184       else
 185         // point is out of bounds
 186         return Standard_False;
 187     }
 188     else
 189     {
 190       // Direction is parallel to the border.
 191       // Check that the point is not out of bounds
 192       if (aLeft  >  Precision::PConfusion() ||
 193           aRight < -Precision::PConfusion())
 194       {
 195         return Standard_False;
 196       }
 197     }
 198   }
 199   return Standard_True;
 200 }
 201
 202 //=============================================================================
 203 //function : MinimizeDirection
 204 //purpose  : Solves 1D minimization problem when point and directions
 205 //           are known.
 206 //=============================================================================
 207 static Standard_Boolean MinimizeDirection(math_Vector&       P,
 208                                           Standard_Real      F0,
 209                                           math_Vector&       Gr,
 210                                           math_Vector&       Dir,
 211                                           Standard_Real&     Result,
 212                                           DirFunction&       F,
 213                                           Standard_Boolean   isBounds,
 214                                           const math_Vector& theLeft,
 215                                           const math_Vector& theRight)
 216 {
 217   Standard_Real lambda;
 218   if (!ComputeInitScale(F0, Dir, Gr, lambda))
 219     return Standard_False;
 220
 221   // by default the scaling range is unlimited
 222   Standard_Real aMinLambda = -Precision::Infinite();
 223   Standard_Real aMaxLambda = Precision::Infinite();
 224   if (isBounds)
 225   {
 226     // limit the scaling range taking into account the bounds
 227     if (!ComputeMinMaxScale(P, Dir, theLeft, theRight, aMinLambda, aMaxLambda))
 228       return Standard_False;
 229
 230     if (aMinLambda > -Precision::PConfusion() && aMaxLambda < Precision::PConfusion())
 231     {
 232       // Point is on the border and the direction shows outside.
 233       // Make direction to go along the border
 234       for (Standard_Integer anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
 235       {
 236         if ((Abs(P(anIdx) - theRight(anIdx)) < Precision::PConfusion() && Dir(anIdx) > 0.0) ||
 237             (Abs(P(anIdx) - theLeft(anIdx))  < Precision::PConfusion() && Dir(anIdx) < 0.0))
 238         {
 239           Dir(anIdx) = 0.0;
 240         }
 241       }
 242
 243       // re-compute scale values with new direction
 244       if (!ComputeInitScale(F0, Dir, Gr, lambda))
 245         return Standard_False;
 246       if (!ComputeMinMaxScale(P, Dir, theLeft, theRight, aMinLambda, aMaxLambda))
 247         return Standard_False;
 248     }
 249     lambda = Min(lambda, aMaxLambda);
 250     lambda = Max(lambda, aMinLambda);
 251   }
 252
 253   F.Initialize(P, Dir);
 254   Standard_Real F1;
 255   if (!F.Value(lambda, F1))
 256     return Standard_False;
 257
 258   math_BracketMinimum Bracket(0.0, lambda);
 259   if (isBounds)
 260     Bracket.SetLimits(aMinLambda, aMaxLambda);
 261   Bracket.SetFA(F0);
 262   Bracket.SetFB(F1);
 263   Bracket.Perform(F);
 264   if (Bracket.IsDone())
 265   {
 266     // find minimum inside the bracket
 267     Standard_Real ax, xx, bx, Fax, Fxx, Fbx;
 268     Bracket.Values(ax, xx, bx);
 269     Bracket.FunctionValues(Fax, Fxx, Fbx);
 270
 271     Standard_Integer niter = 100;
 272     Standard_Real tol = 1.e-03;
 273     math_BrentMinimum Sol(tol, Fxx, niter, 1.e-08);
 274     Sol.Perform(F, ax, xx, bx);
 275     if (Sol.IsDone())
 276     {
 277       Standard_Real Scale = Sol.Location();
 278       Result = Sol.Minimum();
 279       Dir.Multiply(Scale);
 280       P.Add(Dir);
 281       return Standard_True;
 282     }
 283   }
 284   else if (isBounds)
 285   {
 286     // Bracket definition is failure. If the bounds are defined then
 287     // set current point to intersection with bounds
 288     Standard_Real aFMin, aFMax;
 289     if (!F.Value(aMinLambda, aFMin))
 290       return Standard_False;
 291     if (!F.Value(aMaxLambda, aFMax))
 292       return Standard_False;
 293     Standard_Real aBestLambda;
 294     if (aFMin < aFMax)
 295     {
 296       aBestLambda = aMinLambda;
 297       Result = aFMin;
 298     }
 299     else
 300     {
 301       aBestLambda = aMaxLambda;
 302       Result = aFMax;
 303     }
 304     Dir.Multiply(aBestLambda);
 305     P.Add(Dir);
 306     return Standard_True;
 307   }
 308   return Standard_False;
 309 }
 310
 311 //=============================================================================
 312 //function : Perform
 313 //purpose  : Performs minimization problem using BFGS method.
 314 //=============================================================================
 315 void  math_BFGS::Perform(math_MultipleVarFunctionWithGradient& F,
 316                          const math_Vector&                    StartingPoint)
 317 {
 318   const Standard_Integer n = TheLocation.Length();
 319   Standard_Boolean Good = Standard_True;
 320   Standard_Integer j, i;
 321   Standard_Real fae, fad, fac;
 322
 323   math_Vector xi(1, n), dg(1, n), hdg(1, n);
 324   math_Matrix hessin(1, n, 1, n);
 325   hessin.Init(0.0);
 326
 327   math_Vector Temp1(1, n);
 328   math_Vector Temp2(1, n);
 329   math_Vector Temp3(1, n);
 330   math_Vector Temp4(1, n);
 331   DirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
 332
 333   TheLocation = StartingPoint;
 334   Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
 335   if (!Good)
 336   {
 337     Done = Standard_False;
 338     TheStatus = math_FunctionError;
 339     return;
 340   }
 341   for (i = 1; i <= n; i++)
 342   {
 343     hessin(i, i) = 1.0;
 344     xi(i) = -TheGradient(i);
 345   }
 346
 347
 348   for (nbiter = 1; nbiter <= Itermax; nbiter++)
 349   {
 350     TheMinimum = PreviousMinimum;
 351     const Standard_Boolean IsGood = MinimizeDirection(TheLocation, TheMinimum, TheGradient,
 352                                                       xi, TheMinimum, F_Dir, myIsBoundsDefined,
 353                                                       myLeft, myRight);
 354
 355     if (IsSolutionReached(F))
 356     {
 357       Done = Standard_True;
 358       TheStatus = math_OK;
 359       return;
 360     }
 361
 362     if (!IsGood)
 363     {
 364       Done = Standard_False;
 365       TheStatus = math_DirectionSearchError;
 366       return;
 367     }
 368     PreviousMinimum = TheMinimum;
 369
 370     dg = TheGradient;
 371
 372     Good = F.Values(TheLocation, TheMinimum, TheGradient);
 373     if (!Good)
 374     {
 375       Done = Standard_False;
 376       TheStatus = math_FunctionError;
 377       return;
 378     }
 379
 380     for (i = 1; i <= n; i++)
 381       dg(i) = TheGradient(i) - dg(i);
 382
 383     for (i = 1; i <= n; i++)
 384     {
 385       hdg(i) = 0.0;
 386       for (j = 1; j <= n; j++)
 387         hdg(i) += hessin(i, j) * dg(j);
 388     }
 389
 390     fac = fae = 0.0;
 391     for (i = 1; i <= n; i++)
 392     {
 393       fac += dg(i) * xi(i);
 394       fae += dg(i) * hdg(i);
 395     }
 396     fac = 1.0 / fac;
 397     fad = 1.0 / fae;
 398
 399     for (i = 1; i <= n; i++)
 400       dg(i) = fac * xi(i) - fad * hdg(i);
 401
 402     for (i = 1; i <= n; i++)
 403     {
 404       for (j = 1; j <= n; j++)
 405         hessin(i, j) += fac * xi(i) * xi(j) - fad * hdg(i) * hdg(j) + fae * dg(i) * dg(j);
 406     }
 407
 408     for (i = 1; i <= n; i++)
 409     {
 410       xi(i) = 0.0;
 411       for (j = 1; j <= n; j++)
 412         xi(i) -= hessin(i, j) * TheGradient(j);
 413     }
 414   }
 415   Done = Standard_False;
 416   TheStatus = math_TooManyIterations;
 417   return;
 418 }
 419
 420 //=============================================================================
 421 //function : IsSolutionReached
 422 //purpose  : Checks whether solution reached or not.
 423 //=============================================================================
 424 Standard_Boolean math_BFGS::IsSolutionReached(math_MultipleVarFunctionWithGradient&) const
 425 {
 426
 427   return 2.0 * fabs(TheMinimum - PreviousMinimum) <=
 428     XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ);
 429 }
 430
 431 //=============================================================================
 432 //function : math_BFGS
 433 //purpose  : Constructor.
 434 //=============================================================================
 435 math_BFGS::math_BFGS(const Standard_Integer     NbVariables,
 436                      const Standard_Real        Tolerance,
 437                      const Standard_Integer     NbIterations,
 438                      const Standard_Real        ZEPS)
 439 : TheStatus(math_OK),
 440   TheLocation(1, NbVariables),
 441   TheGradient(1, NbVariables),
 442   PreviousMinimum(0.),
 443   TheMinimum(0.),
 444   XTol(Tolerance),
 445   EPSZ(ZEPS),
 446   nbiter(0),
 447   myIsBoundsDefined(Standard_False),
 448   myLeft(1, NbVariables, 0.0),
 449   myRight(1, NbVariables, 0.0),
 450   Done(Standard_False),
 451   Itermax(NbIterations)
 452 {
 453 }
 454
 455 //=============================================================================
 456 //function : ~math_BFGS
 457 //purpose  : Destructor.
 458 //=============================================================================
 459 math_BFGS::~math_BFGS()
 460 {
 461 }
 462
 463 //=============================================================================
 464 //function : Dump
 465 //purpose  : Prints dump.
 466 //=============================================================================
 467 void math_BFGS::Dump(Standard_OStream& o) const
 468 {
 469
 470   o << "math_BFGS resolution: ";
 471   if (Done)
 472   {
 473     o << " Status = Done \n";
 474     o << " Location Vector = " << Location() << "\n";
 475     o << " Minimum value = " << Minimum() << "\n";
 476     o << " Number of iterations = " << NbIterations() << "\n";
 477   }
 478   else
 479     o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
 480 }
 481
 482 //=============================================================================
 483 //function : SetBoundary
 484 //purpose  : Set boundaries for conditional optimization
 485 //=============================================================================
 486 void math_BFGS::SetBoundary(const math_Vector& theLeftBorder,
 487                             const math_Vector& theRightBorder)
 488 {
 489   myLeft = theLeftBorder;
 490   myRight = theRightBorder;
 491   myIsBoundsDefined = Standard_True;
 492 }