src/math/math_GlobOptMin.cxx

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