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