1 // Created on: 2014-01-20
2 // Created by: Alexaner Malyshev
3 // Copyright (c) 2014-2015 OPEN CASCADE SAS
5 // This file is part of Open CASCADE Technology software library.
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.
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement
16 #include <math_GlobOptMin.hxx>
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>
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()),
44 myIsConstLocked(Standard_False),
49 myExpandCoeff(1, myN),
50 myCellSize(0, myN - 1),
51 myFilter(theFunc->NbVariables())
58 myIsFindSingleSolution = Standard_False;
59 myFunctionalMinimalValue = -Precision::Infinite();
63 for(i = 1; i <= myN; i++)
72 for(i = 1; i <= myN; i++)
74 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
77 myExpandCoeff(1) = 1.0;
78 for(i = 2; i <= myN; i++)
80 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
83 myTol = theDiscretizationTol;
84 mySameTol = theSameTol;
86 const Standard_Integer aMaxSquareSearchSol = 200;
87 Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
88 myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
92 myDone = Standard_False;
95 //=======================================================================
96 //function : SetGlobalParams
97 //purpose : Set parameters without memory allocation.
98 //=======================================================================
99 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
100 const math_Vector& theA,
101 const math_Vector& theB,
102 const Standard_Real theC,
103 const Standard_Real theDiscretizationTol,
104 const Standard_Real theSameTol)
114 for(i = 1; i <= myN; i++)
116 myGlobA(i) = theA(i);
117 myGlobB(i) = theB(i);
123 for(i = 1; i <= myN; i++)
125 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
128 myExpandCoeff(1) = 1.0;
129 for(i = 2; i <= myN; i++)
131 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
134 myTol = theDiscretizationTol;
135 mySameTol = theSameTol;
140 myDone = Standard_False;
143 //=======================================================================
144 //function : SetLocalParams
145 //purpose : Set parameters without memory allocation.
146 //=======================================================================
147 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
148 const math_Vector& theLocalB)
153 for(i = 1; i <= myN; i++)
155 myA(i) = theLocalA(i);
156 myB(i) = theLocalB(i);
159 for(i = 1; i <= myN; i++)
161 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
164 myExpandCoeff(1) = 1.0;
165 for(i = 2; i <= myN; i++)
167 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
170 myDone = Standard_False;
173 //=======================================================================
175 //purpose : Set algorithm tolerances.
176 //=======================================================================
177 void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
178 const Standard_Real theSameTol)
180 myTol = theDiscretizationTol;
181 mySameTol = theSameTol;
184 //=======================================================================
186 //purpose : Get algorithm tolerances.
187 //=======================================================================
188 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
189 Standard_Real& theSameTol)
191 theDiscretizationTol = myTol;
192 theSameTol = mySameTol;
195 //=======================================================================
196 //function : ~math_GlobOptMin
198 //=======================================================================
199 math_GlobOptMin::~math_GlobOptMin()
203 //=======================================================================
205 //purpose : Compute Global extremum point
206 //=======================================================================
207 // In this algo indexes started from 1, not from 0.
208 void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
212 // Compute parameters range
213 Standard_Real minLength = RealLast();
214 Standard_Real maxLength = RealFirst();
215 for(i = 1; i <= myN; i++)
217 Standard_Real currentLength = myB(i) - myA(i);
218 if (currentLength < minLength)
219 minLength = currentLength;
220 if (currentLength > maxLength)
221 maxLength = currentLength;
224 if (minLength < Precision::PConfusion())
227 cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << endl;
233 if (!myIsConstLocked)
235 // Compute initial value for myC.
236 computeInitialValues();
239 myE1 = minLength * myTol;
240 myE2 = maxLength * myTol;
242 myIsFindSingleSolution = isFindSingleSolution;
243 if (isFindSingleSolution)
245 // Run local optimization if current value better than optimal.
251 myE3 = - maxLength * myTol / 4.0;
253 myE3 = - maxLength * myTol * myC / 4.0;
256 // Search single solution and current solution in its neighborhood.
257 if (CheckFunctionalStopCriteria())
259 myDone = Standard_True;
264 isFirstCellFilterInvoke = Standard_True;
265 computeGlobalExtremum(myN);
267 myDone = Standard_True;
270 //=======================================================================
271 //function : computeLocalExtremum
273 //=======================================================================
274 Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
275 Standard_Real& theVal,
276 math_Vector& theOutPnt)
281 if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
283 math_MultipleVarFunctionWithHessian* aTmp =
284 dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
285 math_NewtonMinimum newtonMinimum(*aTmp);
286 newtonMinimum.SetBoundary(myGlobA, myGlobB);
287 newtonMinimum.Perform(*aTmp, thePnt);
289 if (newtonMinimum.IsDone())
291 newtonMinimum.Location(theOutPnt);
292 theVal = newtonMinimum.Minimum();
294 else return Standard_False;
298 if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
300 math_MultipleVarFunctionWithGradient* aTmp =
301 dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
302 math_BFGS bfgs(aTmp->NbVariables());
303 bfgs.Perform(*aTmp, thePnt);
306 bfgs.Location(theOutPnt);
307 theVal = bfgs.Minimum();
309 else return Standard_False;
312 // Powell method used.
313 if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
315 math_Matrix m(1, myN, 1, myN, 0.0);
316 for(i = 1; i <= myN; i++)
319 math_Powell powell(*myFunc, 1e-10);
320 powell.Perform(*myFunc, thePnt, m);
324 powell.Location(theOutPnt);
325 theVal = powell.Minimum();
327 else return Standard_False;
330 if (isInside(theOutPnt))
331 return Standard_True;
333 return Standard_False;
336 //=======================================================================
337 //function : computeInitialValues
339 //=======================================================================
340 void math_GlobOptMin::computeInitialValues()
343 math_Vector aCurrPnt(1, myN);
344 math_Vector aBestPnt(1, myN);
345 math_Vector aParamStep(1, myN);
346 Standard_Real aCurrVal = RealLast();
348 // Lipchitz const approximation.
349 Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
350 Standard_Integer aPntNb = 13;
351 myFunc->Value(myA, aPrevValDiag);
352 aPrevValProj = aPrevValDiag;
353 Standard_Real aStep = (myB - myA).Norm() / aPntNb;
354 aParamStep = (myB - myA) / aPntNb;
355 for(i = 1; i <= aPntNb; i++)
357 aCurrPnt = myA + aParamStep * i;
359 // Walk over diagonal.
360 myFunc->Value(aCurrPnt, aCurrVal);
361 aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
362 aPrevValDiag = aCurrVal;
364 // Walk over diag in projected space aPnt(1) = myA(1) = const.
365 aCurrPnt(1) = myA(1);
366 myFunc->Value(aCurrPnt, aCurrVal);
367 aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
368 aPrevValProj = aCurrVal;
372 aLipConst *= Sqrt(myN) / aStep;
373 if (aLipConst < myC * 0.1)
374 myC = Max(aLipConst * 0.1, 0.01);
375 else if (aLipConst > myC * 5)
376 myC = Min(myC * 5, 50.0);
378 // Clear all solutions except one.
379 if (myY.Size() != myN)
381 for(i = 1; i <= myN; i++)
382 aBestPnt(i) = myY(i);
384 for(i = 1; i <= myN; i++)
385 myY.Append(aBestPnt(i));
390 //=======================================================================
391 //function : ComputeGlobalExtremum
393 //=======================================================================
394 void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
397 Standard_Real d; // Functional in moved point.
398 Standard_Real val = RealLast(); // Local extrema computed in moved point.
399 Standard_Real aStepBestValue = RealLast();
400 math_Vector aStepBestPoint(1, myN);
401 Standard_Boolean isInside = Standard_False;
403 Standard_Boolean isReached = Standard_False;
406 for(myX(j) = myA(j) + myE1;
407 (myX(j) < myB(j) + myE1) && (!isReached);
413 isReached = Standard_True;
416 if (CheckFunctionalStopCriteria())
417 return; // Best possible value is obtained.
421 isInside = Standard_False;
422 myFunc->Value(myX, d);
423 r = (d + myZ * myC * myLastStep - myF) * myZ;
426 isInside = computeLocalExtremum(myX, val, myTmp);
428 aStepBestValue = (isInside && (val < d))? val : d;
429 aStepBestPoint = (isInside && (val < d))? myTmp : myX;
431 // Solutions are close to each other
432 // and it is allowed to have more than one solution.
433 if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
434 !myIsFindSingleSolution)
436 if (!isStored(aStepBestPoint))
438 if ((aStepBestValue - myF) * myZ > 0.0)
439 myF = aStepBestValue;
440 for(i = 1; i <= myN; i++)
441 myY.Append(aStepBestPoint(i));
446 // New best solution:
447 // new point is out of (mySameTol * 0.01) surrounding or
448 // new point is better than old + single point search.
449 Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
450 if (aFunctionalDelta > mySameTol * 0.01 ||
451 (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
454 myF = aStepBestValue;
456 for(i = 1; i <= myN; i++)
457 myY.Append(aStepBestPoint(i));
460 isFirstCellFilterInvoke = Standard_True;
463 if (CheckFunctionalStopCriteria())
464 return; // Best possible value is obtained.
466 myV(1) = Min(myE2 + Abs(myF - d) / myC, myMaxV(1));
471 myV(j) = RealLast() / 2.0;
472 computeGlobalExtremum(j - 1);
474 // Nullify steps on lower dimensions.
475 for(i = 1; i < j; i++)
478 // Compute step in (j + 1) dimension according to scale.
481 Standard_Real aUpperDimStep = myV(j) * myExpandCoeff(j + 1);
482 if (myV(j + 1) > aUpperDimStep)
484 if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
485 myV(j + 1) = myMaxV(j + 1);
487 myV(j + 1) = aUpperDimStep;
493 //=======================================================================
494 //function : IsInside
496 //=======================================================================
497 Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
501 for(i = 1; i <= myN; i++)
503 if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
504 return Standard_False;
507 return Standard_True;
509 //=======================================================================
510 //function : IsStored
512 //=======================================================================
513 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
515 Standard_Integer i,j;
516 Standard_Boolean isSame = Standard_True;
517 math_Vector aTol(1, myN);
518 aTol = (myB - myA) * mySameTol;
520 // C1 * n^2 = C2 * 3^dim * n
521 if (mySolCount < myMinCellFilterSol)
523 for(i = 0; i < mySolCount; i++)
525 isSame = Standard_True;
526 for(j = 1; j <= myN; j++)
528 if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
530 isSame = Standard_False;
534 if (isSame == Standard_True)
535 return Standard_True;
540 NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion());
541 if (isFirstCellFilterInvoke)
543 myFilter.Reset(myCellSize);
545 // Copy initial data into cell filter.
546 for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
548 math_Vector aVec(1, myN);
549 for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
550 aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
552 myFilter.Add(aVec, aVec);
556 isFirstCellFilterInvoke = Standard_False;
558 math_Vector aLow(1, myN), anUp(1, myN);
559 anInspector.Shift(thePnt, myCellSize, aLow, anUp);
561 anInspector.ClearFind();
562 anInspector.SetCurrent(thePnt);
563 myFilter.Inspect(aLow, anUp, anInspector);
564 if (!anInspector.isFind())
566 // Point is out of close cells, add new one.
567 myFilter.Add(thePnt, thePnt);
570 return Standard_False;
573 //=======================================================================
574 //function : NbExtrema
576 //=======================================================================
577 Standard_Integer math_GlobOptMin::NbExtrema()
582 //=======================================================================
585 //=======================================================================
586 Standard_Real math_GlobOptMin::GetF()
591 //=======================================================================
592 //function : SetFunctionalMinimalValue
594 //=======================================================================
595 void math_GlobOptMin::SetFunctionalMinimalValue(const Standard_Real theMinimalValue)
597 myFunctionalMinimalValue = theMinimalValue;
600 //=======================================================================
601 //function : GetFunctionalMinimalValue
603 //=======================================================================
604 Standard_Real math_GlobOptMin::GetFunctionalMinimalValue()
606 return myFunctionalMinimalValue;
609 //=======================================================================
612 //=======================================================================
613 Standard_Boolean math_GlobOptMin::isDone()
618 //=======================================================================
621 //=======================================================================
622 void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
626 for(j = 1; j <= myN; j++)
627 theSol(j) = myY((theIndex - 1) * myN + j);
630 //=======================================================================
631 //function : initCellSize
633 //=======================================================================
634 void math_GlobOptMin::initCellSize()
636 for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
638 myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
639 * Precision::PConfusion() / (2.0 * Sqrt(2.0));
643 //=======================================================================
644 //function : CheckFunctionalStopCriteria
646 //=======================================================================
647 Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
649 // Search single solution and current solution in its neighborhood.
650 if (myIsFindSingleSolution &&
651 Abs (myF - myFunctionalMinimalValue) < mySameTol * 0.01)
652 return Standard_True;
654 return Standard_False;
657 //=======================================================================
658 //function : ComputeInitSol
660 //=======================================================================
661 void math_GlobOptMin::ComputeInitSol()
663 Standard_Real aCurrVal, aBestVal;
664 math_Vector aCurrPnt(1, myN);
665 math_Vector aBestPnt(1, myN);
666 math_Vector aParamStep(1, myN);
667 // Check functional value in midpoint, lower and upper border points and
668 // in each point try to perform local optimization.
669 aBestPnt = (myGlobA + myGlobB) * 0.5;
670 myFunc->Value(aBestPnt, aBestVal);
673 for(i = 1; i <= 3; i++)
675 aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
677 if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
679 // Local search tries to find better solution than current point.
680 if (aCurrVal < aBestVal)
690 for(i = 1; i <= myN; i++)
691 myY.Append(aBestPnt(i));
694 myDone = Standard_False;
697 //=======================================================================
698 //function : SetLipConstState
700 //=======================================================================
701 void math_GlobOptMin::SetLipConstState(const Standard_Boolean theFlag)
703 myIsConstLocked = theFlag;