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   myX(1, myN),
  45   myTmp(1, myN),
  46   myV(1, myN),
  47   myMaxV(1, myN),
  48   myExpandCoeff(1, myN),
  49   myCellSize(0, myN - 1),
  50   myFilter(theFunc->NbVariables())
  51 {
  52   Standard_Integer i;
  53
  54   myFunc = theFunc;
  55   myC = theC;
  56   myIsFindSingleSolution = Standard_False;
  57   myZ = -1;
  58   mySolCount = 0;
  59
  60   for(i = 1; i <= myN; i++)
  61   {
  62     myGlobA(i) = theA(i);
  63     myGlobB(i) = theB(i);
  64
  65     myA(i) = theA(i);
  66     myB(i) = theB(i);
  67   }
  68
  69   for(i = 1; i <= myN; i++)
  70   {
  71     myMaxV(i) = (myB(i) - myA(i)) / 3.0;
  72   }
  73
  74   myExpandCoeff(1) = 1.0;
  75   for(i = 2; i <= myN; i++)
  76   {
  77     myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
  78   }
  79
  80   myTol = theDiscretizationTol;
  81   mySameTol = theSameTol;
  82
  83   const Standard_Integer aMaxSquareSearchSol = 200;
  84   Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
  85   myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
  86   initCellSize();
  87
  88   myDone = Standard_False;
  89 }
  90
  91 //=======================================================================
  92 //function : SetGlobalParams
  93 //purpose  : Set params without memory allocation.
  94 //=======================================================================
  95 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
  96                                       const math_Vector& theA,
  97                                       const math_Vector& theB,
  98                                       const Standard_Real theC,
  99                                       const Standard_Real theDiscretizationTol,
 100                                       const Standard_Real theSameTol)
 101 {
 102   Standard_Integer i;
 103
 104   myFunc = theFunc;
 105   myC = theC;
 106   myZ = -1;
 107   mySolCount = 0;
 108
 109   for(i = 1; i <= myN; i++)
 110   {
 111     myGlobA(i) = theA(i);
 112     myGlobB(i) = theB(i);
 113
 114     myA(i) = theA(i);
 115     myB(i) = theB(i);
 116   }
 117
 118   for(i = 1; i <= myN; i++)
 119   {
 120     myMaxV(i) = (myB(i) - myA(i)) / 3.0;
 121   }
 122
 123   myExpandCoeff(1) = 1.0;
 124   for(i = 2; i <= myN; i++)
 125   {
 126     myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
 127   }
 128
 129   myTol = theDiscretizationTol;
 130   mySameTol = theSameTol;
 131
 132   initCellSize();
 133
 134   myDone = Standard_False;
 135 }
 136
 137 //=======================================================================
 138 //function : SetLocalParams
 139 //purpose  : Set params without memory allocation.
 140 //=======================================================================
 141 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
 142                                      const math_Vector& theLocalB)
 143 {
 144   Standard_Integer i;
 145
 146   myZ = -1;
 147   mySolCount = 0;
 148
 149   for(i = 1; i <= myN; i++)
 150   {
 151     myA(i) = theLocalA(i);
 152     myB(i) = theLocalB(i);
 153   }
 154
 155   for(i = 1; i <= myN; i++)
 156   {
 157     myMaxV(i) = (myB(i) - myA(i)) / 3.0;
 158   }
 159
 160   myExpandCoeff(1) = 1.0;
 161   for(i = 2; i <= myN; i++)
 162   {
 163     myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
 164   }
 165
 166   myDone = Standard_False;
 167 }
 168
 169 //=======================================================================
 170 //function : SetTol
 171 //purpose  : Set algorithm tolerances.
 172 //=======================================================================
 173 void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
 174                              const Standard_Real theSameTol)
 175 {
 176   myTol = theDiscretizationTol;
 177   mySameTol = theSameTol;
 178 }
 179
 180 //=======================================================================
 181 //function : GetTol
 182 //purpose  : Get algorithm tolerances.
 183 //=======================================================================
 184 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
 185                              Standard_Real& theSameTol)
 186 {
 187   theDiscretizationTol = myTol;
 188   theSameTol = mySameTol;
 189 }
 190
 191 //=======================================================================
 192 //function : ~math_GlobOptMin
 193 //purpose  :
 194 //=======================================================================
 195 math_GlobOptMin::~math_GlobOptMin()
 196 {
 197 }
 198
 199 //=======================================================================
 200 //function : Perform
 201 //purpose  : Compute Global extremum point
 202 //=======================================================================
 203 // In this algo indexes started from 1, not from 0.
 204 void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
 205 {
 206   Standard_Integer i;
 207
 208   // Compute parameters range
 209   Standard_Real minLength = RealLast();
 210   Standard_Real maxLength = RealFirst();
 211   for(i = 1; i <= myN; i++)
 212   {
 213     Standard_Real currentLength = myB(i) - myA(i);
 214     if (currentLength < minLength)
 215       minLength = currentLength;
 216     if (currentLength > maxLength)
 217       maxLength = currentLength;
 218   }
 219
 220   if (minLength < Precision::PConfusion())
 221   {
 222     #ifdef OCCT_DEBUG
 223     cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << endl;
 224     #endif
 225
 226     return;
 227   }
 228
 229   // Compute initial values for myF, myY, myC.
 230   computeInitialValues();
 231
 232   myE1 = minLength * myTol;
 233   myE2 = maxLength * myTol;
 234
 235   myIsFindSingleSolution = isFindSingleSolution;
 236   if (isFindSingleSolution)
 237   {
 238     // Run local optimization
 239     // if current value better than optimal.
 240     myE3 = 0.0;
 241   }
 242   else
 243   {
 244     if (myC > 1.0)
 245       myE3 = - maxLength * myTol / 4.0;
 246     else
 247       myE3 = - maxLength * myTol * myC / 4.0;
 248   }
 249
 250   isFirstCellFilterInvoke = Standard_True;
 251   computeGlobalExtremum(myN);
 252
 253   myDone = Standard_True;
 254 }
 255
 256 //=======================================================================
 257 //function : computeLocalExtremum
 258 //purpose  :
 259 //=======================================================================
 260 Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
 261                                                        Standard_Real& theVal,
 262                                                        math_Vector& theOutPnt)
 263 {
 264   Standard_Integer i;
 265
 266   //Newton method
 267   if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
 268   {
 269     math_MultipleVarFunctionWithHessian* myTmp =
 270       dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
 271     math_NewtonMinimum newtonMinimum(*myTmp);
 272     newtonMinimum.SetBoundary(myGlobA, myGlobB);
 273     newtonMinimum.Perform(*myTmp, thePnt);
 274
 275     if (newtonMinimum.IsDone())
 276     {
 277       newtonMinimum.Location(theOutPnt);
 278       theVal = newtonMinimum.Minimum();
 279     }
 280     else return Standard_False;
 281   } else
 282
 283   // BFGS method used.
 284   if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
 285   {
 286     math_MultipleVarFunctionWithGradient* myTmp =
 287       dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
 288     math_BFGS bfgs(myTmp->NbVariables());
 289     bfgs.Perform(*myTmp, thePnt);
 290     if (bfgs.IsDone())
 291     {
 292       bfgs.Location(theOutPnt);
 293       theVal = bfgs.Minimum();
 294     }
 295     else return Standard_False;
 296   } else
 297
 298   // Powell method used.
 299   if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
 300   {
 301     math_Matrix m(1, myN, 1, myN, 0.0);
 302     for(i = 1; i <= myN; i++)
 303       m(1, 1) = 1.0;
 304
 305     math_Powell powell(*myFunc, 1e-10);
 306     powell.Perform(*myFunc, thePnt, m);
 307
 308     if (powell.IsDone())
 309     {
 310       powell.Location(theOutPnt);
 311       theVal = powell.Minimum();
 312     }
 313     else return Standard_False;
 314   }
 315
 316   if (isInside(theOutPnt))
 317     return Standard_True;
 318   else
 319     return Standard_False;
 320 }
 321
 322 //=======================================================================
 323 //function : computeInitialValues
 324 //purpose  :
 325 //=======================================================================
 326 void math_GlobOptMin::computeInitialValues()
 327 {
 328   Standard_Integer i;
 329   math_Vector aCurrPnt(1, myN);
 330   math_Vector aBestPnt(1, myN);
 331   math_Vector aParamStep(1, myN);
 332   Standard_Real aCurrVal = RealLast();
 333   Standard_Real aBestVal = RealLast();
 334
 335   // Check functional value in midpoint, low and upp point border and
 336   // in each point try to perform local optimization.
 337   aBestPnt = (myA + myB) * 0.5;
 338   myFunc->Value(aBestPnt, aBestVal);
 339
 340   for(i = 1; i <= 3; i++)
 341   {
 342     aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
 343
 344     if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
 345     {
 346       // Local Extremum finds better solution than current point.
 347       if (aCurrVal < aBestVal)
 348       {
 349         aBestVal = aCurrVal;
 350         aBestPnt = aCurrPnt;
 351       }
 352     }
 353   }
 354
 355   myF = aBestVal;
 356   myY.Clear();
 357   for(i = 1; i <= myN; i++)
 358     myY.Append(aBestPnt(i));
 359   mySolCount++;
 360
 361   // Lipschitz const approximation
 362   Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
 363   Standard_Integer aPntNb = 13;
 364   myFunc->Value(myA, aPrevValDiag);
 365   aPrevValProj = aPrevValDiag;
 366   Standard_Real aStep = (myB - myA).Norm() / aPntNb;
 367   aParamStep = (myB - myA) / aPntNb;
 368   for(i = 1; i <= aPntNb; i++)
 369   {
 370     aCurrPnt = myA + aParamStep * i;
 371
 372     // Walk over diagonal.
 373     myFunc->Value(aCurrPnt, aCurrVal);
 374     aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
 375     aPrevValDiag = aCurrVal;
 376
 377     // Walk over diag in projected space aPnt(1) = myA(1) = const.
 378     aCurrPnt(1) = myA(1);
 379     myFunc->Value(aCurrPnt, aCurrVal);
 380     aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
 381     aPrevValProj = aCurrVal;
 382   }
 383
 384   aLipConst *= Sqrt(myN) / aStep;
 385
 386   if (aLipConst < myC * 0.1)
 387   {
 388     myC = Max(aLipConst * 0.1, 0.01);
 389   }
 390   else if (aLipConst > myC * 10)
 391   {
 392     myC = Min(myC * 2, 30.0);
 393   }
 394 }
 395
 396 //=======================================================================
 397 //function : ComputeGlobalExtremum
 398 //purpose  :
 399 //=======================================================================
 400 void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
 401 {
 402   Standard_Integer i;
 403   Standard_Real d; // Functional in moved point.
 404   Standard_Real val = RealLast(); // Local extrema computed in moved point.
 405   Standard_Real aStepBestValue = RealLast();
 406   Standard_Real aRealStep = 0.0;
 407   math_Vector aStepBestPoint(1, myN);
 408   Standard_Boolean isInside = Standard_False;
 409   Standard_Real r;
 410
 411   for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j))
 412   {
 413     if (myX(j) > myB(j))
 414       myX(j) = myB(j);
 415
 416     if (j == 1)
 417     {
 418       isInside = Standard_False;
 419       myFunc->Value(myX, d);
 420       r = (d + myZ * myC * aRealStep - myF) * myZ;
 421       if(r > myE3)
 422       {
 423         isInside = computeLocalExtremum(myX, val, myTmp);
 424       }
 425       aStepBestValue = (isInside && (val < d))? val : d;
 426       aStepBestPoint = (isInside && (val < d))? myTmp : myX;
 427
 428       // Solutions are close to each other
 429       // and it is allowed to have more than one solution.
 430       if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
 431           !myIsFindSingleSolution)
 432       {
 433         if (!isStored(aStepBestPoint))
 434         {
 435           if ((aStepBestValue - myF) * myZ > 0.0)
 436             myF = aStepBestValue;
 437           for(i = 1; i <= myN; i++)
 438             myY.Append(aStepBestPoint(i));
 439           mySolCount++;
 440         }
 441       }
 442
 443       // New best solution:
 444       // new point is out of (mySameTol * 0.01) surrounding or
 445       // new point is better than old + single point search.
 446       Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
 447       if (aFunctionalDelta > mySameTol * 0.01 ||
 448          (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
 449       {
 450         mySolCount = 0;
 451         myF = aStepBestValue;
 452         myY.Clear();
 453         for(i = 1; i <= myN; i++)
 454           myY.Append(aStepBestPoint(i));
 455         mySolCount++;
 456
 457         isFirstCellFilterInvoke = Standard_True;
 458       }
 459
 460       aRealStep = myE2 + Abs(myF - d) / myC;
 461       myV(1) = Min(aRealStep, myMaxV(1));
 462     }
 463     else
 464     {
 465       myV(j) = RealLast() / 2.0;
 466       computeGlobalExtremum(j - 1);
 467
 468       // Nullify steps on lower dimensions.
 469       for(i = 1; i < j; i++)
 470         myV(i) = 0.0;
 471     }
 472     // Compute step in (j + 1) dimension according to scale.
 473     if (j < myN)
 474     {
 475       Standard_Real aUpperDimStep =  myV(j) * myExpandCoeff(j + 1);
 476       if (myV(j + 1) > aUpperDimStep)
 477       {
 478         if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
 479           myV(j + 1) = myMaxV(j + 1);
 480         else
 481           myV(j + 1) = aUpperDimStep;
 482       }
 483     }
 484   }
 485 }
 486
 487 //=======================================================================
 488 //function : IsInside
 489 //purpose  :
 490 //=======================================================================
 491 Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
 492 {
 493   Standard_Integer i;
 494
 495  for(i = 1; i <= myN; i++)
 496  {
 497    if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
 498      return Standard_False;
 499  }
 500
 501  return Standard_True;
 502 }
 503 //=======================================================================
 504 //function : IsStored
 505 //purpose  :
 506 //=======================================================================
 507 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
 508 {
 509   Standard_Integer i,j;
 510   Standard_Boolean isSame = Standard_True;
 511   math_Vector aTol(1, myN);
 512   aTol = (myB -  myA) * mySameTol;
 513
 514   // C1 * n^2 = C2 * 3^dim * n
 515   if (mySolCount < myMinCellFilterSol)
 516   {
 517     for(i = 0; i < mySolCount; i++)
 518     {
 519       isSame = Standard_True;
 520       for(j = 1; j <= myN; j++)
 521       {
 522         if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
 523         {
 524           isSame = Standard_False;
 525           break;
 526         }
 527       }
 528       if (isSame == Standard_True)
 529         return Standard_True;
 530     }
 531   }
 532   else
 533   {
 534     NCollection_CellFilter_NDimInspector anInspector(myN, Precision::PConfusion());
 535     if (isFirstCellFilterInvoke)
 536     {
 537       myFilter.Reset(myCellSize);
 538
 539       // Copy initial data into cell filter.
 540       for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
 541       {
 542         math_Vector aVec(1, myN);
 543         for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
 544           aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
 545
 546         myFilter.Add(aVec, aVec);
 547       }
 548     }
 549
 550     isFirstCellFilterInvoke = Standard_False;
 551
 552     math_Vector aLow(1, myN), anUp(1, myN);
 553     anInspector.Shift(thePnt, myCellSize, aLow, anUp);
 554
 555     anInspector.ClearFind();
 556     anInspector.SetCurrent(thePnt);
 557     myFilter.Inspect(aLow, anUp, anInspector);
 558     if (!anInspector.isFind())
 559     {
 560       // Point is out of close cells, add new one.
 561       myFilter.Add(thePnt, thePnt);
 562     }
 563   }
 564   return Standard_False;
 565 }
 566
 567 //=======================================================================
 568 //function : NbExtrema
 569 //purpose  :
 570 //=======================================================================
 571 Standard_Integer math_GlobOptMin::NbExtrema()
 572 {
 573   return mySolCount;
 574 }
 575
 576 //=======================================================================
 577 //function : GetF
 578 //purpose  :
 579 //=======================================================================
 580 Standard_Real math_GlobOptMin::GetF()
 581 {
 582   return myF;
 583 }
 584
 585 //=======================================================================
 586 //function : IsDone
 587 //purpose  :
 588 //=======================================================================
 589 Standard_Boolean math_GlobOptMin::isDone()
 590 {
 591   return myDone;
 592 }
 593
 594 //=======================================================================
 595 //function : Points
 596 //purpose  :
 597 //=======================================================================
 598 void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
 599 {
 600   Standard_Integer j;
 601
 602   for(j = 1; j <= myN; j++)
 603     theSol(j) = myY((theIndex - 1) * myN + j);
 604 }
 605
 606 //=======================================================================
 607 //function : initCellSize
 608 //purpose  :
 609 //=======================================================================
 610 void math_GlobOptMin::initCellSize()
 611 {
 612   for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
 613   {
 614     myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
 615       * Precision::PConfusion() / (2.0 * Sqrt(2.0));
 616   }
 617 }