src/math/math_GlobOptMin.cxx

   1 // Created on: 2014-01-20
   2 // Created by: Alexaner Malyshev
   3 // Copyright (c) 2014-2015 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement
  15
  16 #include <math_GlobOptMin.hxx>
  17
  18 #include <math_BFGS.hxx>
  19 #include <math_Matrix.hxx>
  20 #include <math_MultipleVarFunctionWithGradient.hxx>
  21 #include <math_MultipleVarFunctionWithHessian.hxx>
  22 #include <math_NewtonMinimum.hxx>
  23 #include <math_Powell.hxx>
  24 #include <Standard_Integer.hxx>
  25 #include <Standard_Real.hxx>
  26 #include <Precision.hxx>
  27
  28
  29 //=======================================================================
  30 //function : math_GlobOptMin
  31 //purpose  : Constructor
  32 //=======================================================================
  33 math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
  34                                  const math_Vector& theA,
  35                                  const math_Vector& theB,
  36                                  const Standard_Real theC,
  37                                  const Standard_Real theDiscretizationTol,
  38                                  const Standard_Real theSameTol)
  39 : myN(theFunc->NbVariables()),
  40   myA(1, myN),
  41   myB(1, myN),
  42   myGlobA(1, myN),
  43   myGlobB(1, myN),
  44   myIsConstLocked(Standard_False),
  45   myX(1, myN),
  46   myTmp(1, myN),
  47   myV(1, myN),
  48   myMaxV(1, myN),
  49   myExpandCoeff(1, myN),
  50   myCellSize(0, myN - 1),
  51   myFilter(theFunc->NbVariables()),
  52   myCont(2)
  53 {
  54   Standard_Integer i;
  55
  56   myFunc = theFunc;
  57   myC = theC;
  58   myInitC = theC;
  59   myIsFindSingleSolution = Standard_False;
  60   myFunctionalMinimalValue = -Precision::Infinite();
  61   myZ = -1;
  62   mySolCount = 0;
  63
  64   for(i = 1; i <= myN; i++)
  65   {
  66     myGlobA(i) = theA(i);
  67     myGlobB(i) = theB(i);
  68
  69     myA(i) = theA(i);
  70     myB(i) = theB(i);
  71   }
  72
  73   for(i = 1; i <= myN; i++)
  74   {
  75     myMaxV(i) = (myB(i) - myA(i)) / 3.0;
  76   }
  77
  78   myExpandCoeff(1) = 1.0;
  79   for(i = 2; i <= myN; i++)
  80   {
  81     myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
  82   }
  83
  84   myTol = theDiscretizationTol;
  85   mySameTol = theSameTol;
  86
  87   const Standard_Integer aMaxSquareSearchSol = 200;
  88   Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
  89   myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
  90   initCellSize();
  91   ComputeInitSol();
  92
  93   myDone = Standard_False;
  94 }
  95
  96 //=======================================================================
  97 //function : SetGlobalParams
  98 //purpose  : Set parameters without memory allocation.
  99 //=======================================================================
 100 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
 101                                       const math_Vector& theA,
 102                                       const math_Vector& theB,
 103                                       const Standard_Real theC,
 104                                       const Standard_Real theDiscretizationTol,
 105                                       const Standard_Real theSameTol)
 106 {
 107   Standard_Integer i;
 108
 109   myFunc = theFunc;
 110   myC = theC;
 111   myInitC = theC;
 112   myZ = -1;
 113   mySolCount = 0;
 114
 115   for(i = 1; i <= myN; i++)
 116   {
 117     myGlobA(i) = theA(i);
 118     myGlobB(i) = theB(i);
 119
 120     myA(i) = theA(i);
 121     myB(i) = theB(i);
 122   }
 123
 124   for(i = 1; i <= myN; i++)
 125   {
 126     myMaxV(i) = (myB(i) - myA(i)) / 3.0;
 127   }
 128
 129   myExpandCoeff(1) = 1.0;
 130   for(i = 2; i <= myN; i++)
 131   {
 132     myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
 133   }
 134
 135   myTol = theDiscretizationTol;
 136   mySameTol = theSameTol;
 137
 138   initCellSize();
 139   ComputeInitSol();
 140
 141   myDone = Standard_False;
 142 }
 143
 144 //=======================================================================
 145 //function : SetLocalParams
 146 //purpose  : Set parameters without memory allocation.
 147 //=======================================================================
 148 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
 149                                      const math_Vector& theLocalB)
 150 {
 151   Standard_Integer i;
 152
 153   myZ = -1;
 154   for(i = 1; i <= myN; i++)
 155   {
 156     myA(i) = theLocalA(i);
 157     myB(i) = theLocalB(i);
 158   }
 159
 160   for(i = 1; i <= myN; i++)
 161   {
 162     myMaxV(i) = (myB(i) - myA(i)) / 3.0;
 163   }
 164
 165   myExpandCoeff(1) = 1.0;
 166   for(i = 2; i <= myN; i++)
 167   {
 168     myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
 169   }
 170
 171   myDone = Standard_False;
 172 }
 173
 174 //=======================================================================
 175 //function : SetTol
 176 //purpose  : Set algorithm tolerances.
 177 //=======================================================================
 178 void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
 179                              const Standard_Real theSameTol)
 180 {
 181   myTol = theDiscretizationTol;
 182   mySameTol = theSameTol;
 183 }
 184
 185 //=======================================================================
 186 //function : GetTol
 187 //purpose  : Get algorithm tolerances.
 188 //=======================================================================
 189 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
 190                              Standard_Real& theSameTol)
 191 {
 192   theDiscretizationTol = myTol;
 193   theSameTol = mySameTol;
 194 }
 195
 196 //=======================================================================
 197 //function : ~math_GlobOptMin
 198 //purpose  :
 199 //=======================================================================
 200 math_GlobOptMin::~math_GlobOptMin()
 201 {
 202 }
 203
 204 //=======================================================================
 205 //function : Perform
 206 //purpose  : Compute Global extremum point
 207 //=======================================================================
 208 // In this algo indexes started from 1, not from 0.
 209 void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
 210 {
 211   Standard_Integer i;
 212
 213   // Compute parameters range
 214   Standard_Real minLength = RealLast();
 215   Standard_Real maxLength = RealFirst();
 216   for(i = 1; i <= myN; i++)
 217   {
 218     Standard_Real currentLength = myB(i) - myA(i);
 219     if (currentLength < minLength)
 220       minLength = currentLength;
 221     if (currentLength > maxLength)
 222       maxLength = currentLength;
 223   }
 224
 225   if (minLength < Precision::PConfusion())
 226   {
 227     #ifdef OCCT_DEBUG
 228     cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << endl;
 229     #endif
 230
 231     return;
 232   }
 233
 234   if (!myIsConstLocked)
 235   {
 236     // Compute initial value for myC.
 237     computeInitialValues();
 238   }
 239
 240   myE1 = minLength * myTol;
 241   myE2 = maxLength * myTol;
 242
 243   myIsFindSingleSolution = isFindSingleSolution;
 244   if (isFindSingleSolution)
 245   {
 246     // Run local optimization if current value better than optimal.
 247     myE3 = 0.0;
 248   }
 249   else
 250   {
 251     if (myC > 1.0)
 252       myE3 = - maxLength * myTol / 4.0;
 253     else
 254       myE3 = - maxLength * myTol * myC / 4.0;
 255   }
 256
 257   // Search single solution and current solution in its neighborhood.
 258   if (CheckFunctionalStopCriteria())
 259   {
 260     myDone = Standard_True;
 261     return;
 262   }
 263
 264   myLastStep = 0.0;
 265   isFirstCellFilterInvoke = Standard_True;
 266   computeGlobalExtremum(myN);
 267
 268   myDone = Standard_True;
 269 }
 270
 271 //=======================================================================
 272 //function : computeLocalExtremum
 273 //purpose  :
 274 //=======================================================================
 275 Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
 276                                                        Standard_Real& theVal,
 277                                                        math_Vector& theOutPnt)
 278 {
 279   Standard_Integer i;
 280
 281   //Newton method
 282   if (myCont >= 2 &&
 283       dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
 284   {
 285     math_MultipleVarFunctionWithHessian* aTmp =
 286       dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
 287     math_NewtonMinimum newtonMinimum(*aTmp);
 288     newtonMinimum.SetBoundary(myGlobA, myGlobB);
 289     newtonMinimum.Perform(*aTmp, thePnt);
 290
 291     if (newtonMinimum.IsDone())
 292     {
 293       newtonMinimum.Location(theOutPnt);
 294       theVal = newtonMinimum.Minimum();
 295     }
 296     else return Standard_False;
 297   } else
 298
 299   // BFGS method used.
 300   if (myCont >= 1 &&
 301       dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
 302   {
 303     math_MultipleVarFunctionWithGradient* aTmp =
 304       dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
 305     math_BFGS bfgs(aTmp->NbVariables());
 306     bfgs.Perform(*aTmp, thePnt);
 307     if (bfgs.IsDone())
 308     {
 309       bfgs.Location(theOutPnt);
 310       theVal = bfgs.Minimum();
 311     }
 312     else return Standard_False;
 313   } else
 314
 315   // Powell method used.
 316   if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
 317   {
 318     math_Matrix m(1, myN, 1, myN, 0.0);
 319     for(i = 1; i <= myN; i++)
 320       m(1, 1) = 1.0;
 321
 322     math_Powell powell(*myFunc, 1e-10);
 323     powell.Perform(*myFunc, thePnt, m);
 324
 325     if (powell.IsDone())
 326     {
 327       powell.Location(theOutPnt);
 328       theVal = powell.Minimum();
 329     }
 330     else return Standard_False;
 331   }
 332
 333   if (isInside(theOutPnt))
 334     return Standard_True;
 335   else
 336     return Standard_False;
 337 }
 338
 339 //=======================================================================
 340 //function : computeInitialValues
 341 //purpose  :
 342 //=======================================================================
 343 void math_GlobOptMin::computeInitialValues()
 344 {
 345   Standard_Integer i;
 346   math_Vector aCurrPnt(1, myN);
 347   math_Vector aBestPnt(1, myN);
 348   math_Vector aParamStep(1, myN);
 349   Standard_Real aCurrVal = RealLast();
 350
 351   // Lipchitz const approximation.
 352   Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
 353   Standard_Integer aPntNb = 13;
 354   myFunc->Value(myA, aPrevValDiag);
 355   aPrevValProj = aPrevValDiag;
 356   Standard_Real aStep = (myB - myA).Norm() / aPntNb;
 357   aParamStep = (myB - myA) / aPntNb;
 358   for(i = 1; i <= aPntNb; i++)
 359   {
 360     aCurrPnt = myA + aParamStep * i;
 361
 362     // Walk over diagonal.
 363     myFunc->Value(aCurrPnt, aCurrVal);
 364     aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
 365     aPrevValDiag = aCurrVal;
 366
 367     // Walk over diag in projected space aPnt(1) = myA(1) = const.
 368     aCurrPnt(1) = myA(1);
 369     myFunc->Value(aCurrPnt, aCurrVal);
 370     aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
 371     aPrevValProj = aCurrVal;
 372   }
 373
 374   myC = myInitC;
 375   aLipConst *= Sqrt(myN) / aStep;
 376   if (aLipConst < myC * 0.1)
 377     myC = Max(aLipConst * 0.1, 0.01);
 378   else if (aLipConst > myC * 5)
 379     myC = Min(myC * 5, 50.0);
 380
 381   // Clear all solutions except one.
 382   if (myY.Size() != myN)
 383   {
 384     for(i = 1; i <= myN; i++)
 385       aBestPnt(i) = myY(i);
 386     myY.Clear();
 387     for(i = 1; i <= myN; i++)
 388       myY.Append(aBestPnt(i));
 389   }
 390   mySolCount = 1;
 391 }
 392
 393 //=======================================================================
 394 //function : ComputeGlobalExtremum
 395 //purpose  :
 396 //=======================================================================
 397 void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
 398 {
 399   Standard_Integer i;
 400   Standard_Real d; // Functional in moved point.
 401   Standard_Real val = RealLast(); // Local extrema computed in moved point.
 402   Standard_Real aStepBestValue = RealLast();
 403   math_Vector aStepBestPoint(1, myN);
 404   Standard_Boolean isInside = Standard_False;
 405   Standard_Real r;
 406   Standard_Boolean isReached = Standard_False;
 407
 408
 409   for(myX(j) = myA(j) + myE1;
 410      (myX(j) < myB(j) + myE1) && (!isReached);
 411       myX(j) += myV(j))
 412   {
 413     if (myX(j) > myB(j))
 414     {
 415       myX(j) = myB(j);
 416       isReached = Standard_True;
 417     }
 418
 419     if (CheckFunctionalStopCriteria())
 420       return; // Best possible value is obtained.
 421
 422     if (j == 1)
 423     {
 424       isInside = Standard_False;
 425       myFunc->Value(myX, d);
 426       r = (d + myZ * myC * myLastStep - myF) * myZ;
 427       if(r > myE3)
 428       {
 429         isInside = computeLocalExtremum(myX, val, myTmp);
 430       }
 431       aStepBestValue = (isInside && (val < d))? val : d;
 432       aStepBestPoint = (isInside && (val < d))? myTmp : myX;
 433
 434       // Solutions are close to each other
 435       // and it is allowed to have more than one solution.
 436       if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
 437           !myIsFindSingleSolution)
 438       {
 439         if (!isStored(aStepBestPoint))
 440         {
 441           if ((aStepBestValue - myF) * myZ > 0.0)
 442             myF = aStepBestValue;
 443           for(i = 1; i <= myN; i++)
 444             myY.Append(aStepBestPoint(i));
 445           mySolCount++;
 446         }
 447       }
 448
 449       // New best solution:
 450       // new point is out of (mySameTol * 0.01) surrounding or
 451       // new point is better than old + single point search.
 452       Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
 453       if (aFunctionalDelta > mySameTol * 0.01 ||
 454          (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
 455       {
 456         mySolCount = 0;
 457         myF = aStepBestValue;
 458         myY.Clear();
 459         for(i = 1; i <= myN; i++)
 460           myY.Append(aStepBestPoint(i));
 461         mySolCount++;
 462
 463         isFirstCellFilterInvoke = Standard_True;
 464       }
 465
 466       if (CheckFunctionalStopCriteria())
 467         return; // Best possible value is obtained.
 468
 469       myV(1) = Min(myE2 + Abs(myF - d) / myC, myMaxV(1));
 470       myLastStep = myV(1);
 471     }
 472     else
 473     {
 474       myV(j) = RealLast() / 2.0;
 475       computeGlobalExtremum(j - 1);
 476
 477       // Nullify steps on lower dimensions.
 478       for(i = 1; i < j; i++)
 479         myV(i) = 0.0;
 480     }
 481     // Compute step in (j + 1) dimension according to scale.
 482     if (j < myN)
 483     {
 484       Standard_Real aUpperDimStep =  myV(j) * myExpandCoeff(j + 1);
 485       if (myV(j + 1) > aUpperDimStep)
 486       {
 487         if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
 488           myV(j + 1) = myMaxV(j + 1);
 489         else
 490           myV(j + 1) = aUpperDimStep;
 491       }
 492     }
 493   }
 494 }
 495
 496 //=======================================================================
 497 //function : IsInside
 498 //purpose  :
 499 //=======================================================================
 500 Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
 501 {
 502   Standard_Integer i;
 503
 504  for(i = 1; i <= myN; i++)
 505  {
 506    if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
 507      return Standard_False;
 508  }
 509
 510  return Standard_True;
 511 }
 512 //=======================================================================
 513 //function : IsStored
 514 //purpose  :
 515 //=======================================================================
 516 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
 517 {
 518   Standard_Integer i,j;
 519   Standard_Boolean isSame = Standard_True;
 520   math_Vector aTol(1, myN);
 521   aTol = (myB -  myA) * mySameTol;
 522
 523   // C1 * n^2 = C2 * 3^dim * n
 524   if (mySolCount < myMinCellFilterSol)
 525   {
 526     for(i = 0; i < mySolCount; i++)
 527     {
 528       isSame = Standard_True;
 529       for(j = 1; j <= myN; j++)
 530       {
 531         if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
 532         {
 533           isSame = Standard_False;
 534           break;
 535         }
 536       }
 537       if (isSame == Standard_True)
 538         return Standard_True;
 539     }
 540   }
 541   else
 542   {
 543     NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion());
 544     if (isFirstCellFilterInvoke)
 545     {
 546       myFilter.Reset(myCellSize);
 547
 548       // Copy initial data into cell filter.
 549       for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
 550       {
 551         math_Vector aVec(1, myN);
 552         for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
 553           aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
 554
 555         myFilter.Add(aVec, aVec);
 556       }
 557     }
 558
 559     isFirstCellFilterInvoke = Standard_False;
 560
 561     math_Vector aLow(1, myN), anUp(1, myN);
 562     anInspector.Shift(thePnt, myCellSize, aLow, anUp);
 563
 564     anInspector.ClearFind();
 565     anInspector.SetCurrent(thePnt);
 566     myFilter.Inspect(aLow, anUp, anInspector);
 567     if (!anInspector.isFind())
 568     {
 569       // Point is out of close cells, add new one.
 570       myFilter.Add(thePnt, thePnt);
 571     }
 572   }
 573   return Standard_False;
 574 }
 575
 576 //=======================================================================
 577 //function : NbExtrema
 578 //purpose  :
 579 //=======================================================================
 580 Standard_Integer math_GlobOptMin::NbExtrema()
 581 {
 582   return mySolCount;
 583 }
 584
 585 //=======================================================================
 586 //function : GetF
 587 //purpose  :
 588 //=======================================================================
 589 Standard_Real math_GlobOptMin::GetF()
 590 {
 591   return myF;
 592 }
 593
 594 //=======================================================================
 595 //function : SetFunctionalMinimalValue
 596 //purpose  :
 597 //=======================================================================
 598 void math_GlobOptMin::SetFunctionalMinimalValue(const Standard_Real theMinimalValue)
 599 {
 600   myFunctionalMinimalValue = theMinimalValue;
 601 }
 602
 603 //=======================================================================
 604 //function : GetFunctionalMinimalValue
 605 //purpose  :
 606 //=======================================================================
 607 Standard_Real math_GlobOptMin::GetFunctionalMinimalValue()
 608 {
 609   return myFunctionalMinimalValue;
 610 }
 611
 612 //=======================================================================
 613 //function : IsDone
 614 //purpose  :
 615 //=======================================================================
 616 Standard_Boolean math_GlobOptMin::isDone()
 617 {
 618   return myDone;
 619 }
 620
 621 //=======================================================================
 622 //function : Points
 623 //purpose  :
 624 //=======================================================================
 625 void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
 626 {
 627   Standard_Integer j;
 628
 629   for(j = 1; j <= myN; j++)
 630     theSol(j) = myY((theIndex - 1) * myN + j);
 631 }
 632
 633 //=======================================================================
 634 //function : initCellSize
 635 //purpose  :
 636 //=======================================================================
 637 void math_GlobOptMin::initCellSize()
 638 {
 639   for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
 640   {
 641     myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
 642       * Precision::PConfusion() / (2.0 * Sqrt(2.0));
 643   }
 644 }
 645
 646 //=======================================================================
 647 //function : CheckFunctionalStopCriteria
 648 //purpose  :
 649 //=======================================================================
 650 Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
 651 {
 652   // Search single solution and current solution in its neighborhood.
 653   if (myIsFindSingleSolution &&
 654       Abs (myF - myFunctionalMinimalValue) < mySameTol * 0.01)
 655     return Standard_True;
 656
 657   return Standard_False;
 658 }
 659
 660 //=======================================================================
 661 //function : ComputeInitSol
 662 //purpose  :
 663 //=======================================================================
 664 void math_GlobOptMin::ComputeInitSol()
 665 {
 666   Standard_Real aCurrVal, aBestVal;
 667   math_Vector aCurrPnt(1, myN);
 668   math_Vector aBestPnt(1, myN);
 669   math_Vector aParamStep(1, myN);
 670   // Check functional value in midpoint, lower and upper border points and
 671   // in each point try to perform local optimization.
 672   aBestPnt = (myGlobA + myGlobB) * 0.5;
 673   myFunc->Value(aBestPnt, aBestVal);
 674
 675   Standard_Integer i;
 676   for(i = 1; i <= 3; i++)
 677   {
 678     aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
 679
 680     if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
 681     {
 682       // Local search tries to find better solution than current point.
 683       if (aCurrVal < aBestVal)
 684       {
 685         aBestVal = aCurrVal;
 686         aBestPnt = aCurrPnt;
 687       }
 688     }
 689   }
 690
 691   myF = aBestVal;
 692   myY.Clear();
 693   for(i = 1; i <= myN; i++)
 694     myY.Append(aBestPnt(i));
 695   mySolCount = 1;
 696
 697   myDone = Standard_False;
 698 }
 699
 700 //=======================================================================
 701 //function : SetLipConstState
 702 //purpose  :
 703 //=======================================================================
 704 void math_GlobOptMin::SetLipConstState(const Standard_Boolean theFlag)
 705 {
 706   myIsConstLocked = theFlag;
 707 }