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>
28 //=======================================================================
29 //function : DistanceToBorder
31 //=======================================================================
32 static Standard_Real DistanceToBorder(const math_Vector & theX,
33 const math_Vector & theMin,
34 const math_Vector & theMax)
36 Standard_Real aDist = RealLast();
38 for (Standard_Integer anIdx = theMin.Lower(); anIdx <= theMin.Upper(); ++anIdx)
40 const Standard_Real aDist1 = Abs (theX(anIdx) - theMin(anIdx));
41 const Standard_Real aDist2 = Abs (theX(anIdx) - theMax(anIdx));
43 aDist = Min (aDist, Min (aDist1, aDist2));
50 //=======================================================================
51 //function : math_GlobOptMin
52 //purpose : Constructor
53 //=======================================================================
54 math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
55 const math_Vector& theA,
56 const math_Vector& theB,
57 const Standard_Real theC,
58 const Standard_Real theDiscretizationTol,
59 const Standard_Real theSameTol)
60 : myN(theFunc->NbVariables()),
65 myIsConstLocked(Standard_False),
70 myCellSize(0, myN - 1),
71 myFilter(theFunc->NbVariables()),
79 myIsFindSingleSolution = Standard_False;
80 myFunctionalMinimalValue = -Precision::Infinite();
84 for(i = 1; i <= myN; i++)
93 for(i = 1; i <= myN; i++)
95 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
98 myTol = theDiscretizationTol;
99 mySameTol = theSameTol;
101 const Standard_Integer aMaxSquareSearchSol = 200;
102 Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
103 myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
107 myDone = Standard_False;
110 //=======================================================================
111 //function : SetGlobalParams
112 //purpose : Set parameters without memory allocation.
113 //=======================================================================
114 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
115 const math_Vector& theA,
116 const math_Vector& theB,
117 const Standard_Real theC,
118 const Standard_Real theDiscretizationTol,
119 const Standard_Real theSameTol)
129 for(i = 1; i <= myN; i++)
131 myGlobA(i) = theA(i);
132 myGlobB(i) = theB(i);
138 for(i = 1; i <= myN; i++)
140 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
143 myTol = theDiscretizationTol;
144 mySameTol = theSameTol;
149 myDone = Standard_False;
152 //=======================================================================
153 //function : SetLocalParams
154 //purpose : Set parameters without memory allocation.
155 //=======================================================================
156 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
157 const math_Vector& theLocalB)
162 for(i = 1; i <= myN; i++)
164 myA(i) = theLocalA(i);
165 myB(i) = theLocalB(i);
168 for(i = 1; i <= myN; i++)
170 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
173 myDone = Standard_False;
176 //=======================================================================
178 //purpose : Set algorithm tolerances.
179 //=======================================================================
180 void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
181 const Standard_Real theSameTol)
183 myTol = theDiscretizationTol;
184 mySameTol = theSameTol;
187 //=======================================================================
189 //purpose : Get algorithm tolerances.
190 //=======================================================================
191 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
192 Standard_Real& theSameTol)
194 theDiscretizationTol = myTol;
195 theSameTol = mySameTol;
198 //=======================================================================
200 //purpose : Compute Global extremum point
201 //=======================================================================
202 // In this algo indexes started from 1, not from 0.
203 void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
205 myDone = Standard_False;
207 // Compute parameters range
208 Standard_Real minLength = RealLast();
209 Standard_Real maxLength = RealFirst();
210 for(Standard_Integer i = 1; i <= myN; i++)
212 Standard_Real currentLength = myB(i) - myA(i);
213 if (currentLength < minLength)
214 minLength = currentLength;
215 if (currentLength > maxLength)
216 maxLength = currentLength;
221 if (minLength < Precision::PConfusion())
224 cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << endl;
230 if (!myIsConstLocked)
232 // Compute initial value for myC.
233 computeInitialValues();
236 myE1 = minLength * myTol;
237 myE2 = maxLength * myTol;
239 myIsFindSingleSolution = isFindSingleSolution;
240 if (isFindSingleSolution)
242 // Run local optimization if current value better than optimal.
248 myE3 = - maxLength * myTol / 4.0;
250 myE3 = - maxLength * myTol * myC / 4.0;
253 // Search single solution and current solution in its neighborhood.
254 if (CheckFunctionalStopCriteria())
256 myDone = Standard_True;
261 isFirstCellFilterInvoke = Standard_True;
262 computeGlobalExtremum(myN);
264 myDone = Standard_True;
267 //=======================================================================
268 //function : computeLocalExtremum
270 //=======================================================================
271 Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
272 Standard_Real& theVal,
273 math_Vector& theOutPnt)
279 dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
281 math_MultipleVarFunctionWithHessian* aTmp =
282 dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
283 math_NewtonMinimum newtonMinimum(*aTmp);
284 newtonMinimum.SetBoundary(myGlobA, myGlobB);
285 newtonMinimum.Perform(*aTmp, thePnt);
287 if (newtonMinimum.IsDone())
289 newtonMinimum.Location(theOutPnt);
290 theVal = newtonMinimum.Minimum();
292 else return Standard_False;
297 dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
299 math_MultipleVarFunctionWithGradient* aTmp =
300 dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
301 math_BFGS bfgs(aTmp->NbVariables());
302 bfgs.Perform(*aTmp, thePnt);
305 bfgs.Location(theOutPnt);
306 theVal = bfgs.Minimum();
308 else return Standard_False;
311 // Powell method used.
312 if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
314 math_Matrix m(1, myN, 1, myN, 0.0);
315 for(i = 1; i <= myN; i++)
318 math_Powell powell(*myFunc, 1e-10);
319 powell.Perform(*myFunc, thePnt, m);
323 powell.Location(theOutPnt);
324 theVal = powell.Minimum();
326 else return Standard_False;
329 if (isInside(theOutPnt))
330 return Standard_True;
332 return Standard_False;
335 //=======================================================================
336 //function : computeInitialValues
338 //=======================================================================
339 void math_GlobOptMin::computeInitialValues()
342 math_Vector aCurrPnt(1, myN);
343 math_Vector aBestPnt(1, myN);
344 math_Vector aParamStep(1, myN);
345 Standard_Real aCurrVal = RealLast();
347 // Lipchitz const approximation.
348 Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
349 Standard_Integer aPntNb = 13;
350 myFunc->Value(myA, aPrevValDiag);
351 aPrevValProj = aPrevValDiag;
352 Standard_Real aStep = (myB - myA).Norm() / aPntNb;
353 aParamStep = (myB - myA) / aPntNb;
354 for(i = 1; i <= aPntNb; i++)
356 aCurrPnt = myA + aParamStep * i;
358 // Walk over diagonal.
359 myFunc->Value(aCurrPnt, aCurrVal);
360 aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
361 aPrevValDiag = aCurrVal;
363 // Walk over diag in projected space aPnt(1) = myA(1) = const.
364 aCurrPnt(1) = myA(1);
365 myFunc->Value(aCurrPnt, aCurrVal);
366 aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
367 aPrevValProj = aCurrVal;
371 aLipConst *= Sqrt(myN) / aStep;
372 if (aLipConst < myC * 0.1)
373 myC = Max(aLipConst * 0.1, 0.01);
374 else if (aLipConst > myC * 5)
375 myC = Min(myC * 5, 50.0);
377 // Clear all solutions except one.
378 if (myY.Size() != myN)
380 for(i = 1; i <= myN; i++)
381 aBestPnt(i) = myY(i);
383 for(i = 1; i <= myN; i++)
384 myY.Append(aBestPnt(i));
389 //=======================================================================
390 //function : ComputeGlobalExtremum
392 //=======================================================================
393 void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
396 Standard_Real d = RealLast(), aPrevVal; // Functional in original and moved points.
397 Standard_Real val = RealLast(); // Local extrema computed in moved point.
398 Standard_Real aStepBestValue = RealLast();
399 math_Vector aStepBestPoint(1, myN);
400 Standard_Boolean isInside = Standard_False,
401 isReached = Standard_False;
403 Standard_Real r1, r2, r;
405 for(myX(j) = myA(j) + myE1; !isReached; myX(j) += myV(j))
410 isReached = Standard_True;
413 if (CheckFunctionalStopCriteria())
414 return; // Best possible value is obtained.
418 isInside = Standard_False;
420 myFunc->Value(myX, d);
421 r1 = (d + myZ * myC * myLastStep - myF) * myZ; // Evtushenko estimation.
422 r2 = ((d + aPrevVal - myC * myLastStep) * 0.5 - myF) * myZ; // Shubert / Piyavsky estimation.
426 Standard_Real aSaveParam = myX(1);
428 // Piyavsky midpoint estimation.
429 Standard_Real aParam = (2 * myX(1) - myV(1) ) * 0.5 + (aPrevVal - d) * 0.5 / myC;
430 if (Precision::IsInfinite(aPrevVal))
431 aParam = myX(1) - myV(1) * 0.5; // Protection from upper dimension step.
434 Standard_Real aVal = 0;
435 myFunc->Value(myX, aVal);
438 if ( (aVal < d && aVal < aPrevVal) ||
439 DistanceToBorder(myX, myA, myB) < myE1 ) // Condition optimization case near the border.
441 isInside = computeLocalExtremum(myX, val, myTmp);
444 aStepBestValue = (isInside && (val < d))? val : d;
445 aStepBestPoint = (isInside && (val < d))? myTmp : myX;
447 // Solutions are close to each other
448 // and it is allowed to have more than one solution.
449 if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
450 !myIsFindSingleSolution)
452 if (!isStored(aStepBestPoint))
454 if ((aStepBestValue - myF) * myZ > 0.0)
455 myF = aStepBestValue;
456 for(i = 1; i <= myN; i++)
457 myY.Append(aStepBestPoint(i));
462 // New best solution:
463 // new point is out of (mySameTol * 0.01) surrounding or
464 // new point is better than old + single point search.
465 Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
466 if (aFunctionalDelta > mySameTol * 0.01 ||
467 (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
470 myF = aStepBestValue;
472 for(i = 1; i <= myN; i++)
473 myY.Append(aStepBestPoint(i));
476 isFirstCellFilterInvoke = Standard_True;
479 if (CheckFunctionalStopCriteria())
480 return; // Best possible value is obtained.
482 myV(1) = Min(myE2 + Abs(myF - d) / myC, myMaxV(1));
487 myV(j) = RealLast() / 2.0;
488 computeGlobalExtremum(j - 1);
490 // Nullify steps on lower dimensions.
491 for(i = 1; i < j; i++)
496 Standard_Real aUpperDimStep = Max(myV(j), myE2);
497 if (myV(j + 1) > aUpperDimStep)
499 if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
500 myV(j + 1) = myMaxV(j + 1);
502 myV(j + 1) = aUpperDimStep;
508 //=======================================================================
509 //function : IsInside
511 //=======================================================================
512 Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
516 for(i = 1; i <= myN; i++)
518 if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
519 return Standard_False;
522 return Standard_True;
524 //=======================================================================
525 //function : IsStored
527 //=======================================================================
528 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
530 Standard_Integer i,j;
531 Standard_Boolean isSame = Standard_True;
532 math_Vector aTol(1, myN);
533 aTol = (myB - myA) * mySameTol;
535 // C1 * n^2 = C2 * 3^dim * n
536 if (mySolCount < myMinCellFilterSol)
538 for(i = 0; i < mySolCount; i++)
540 isSame = Standard_True;
541 for(j = 1; j <= myN; j++)
543 if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
545 isSame = Standard_False;
549 if (isSame == Standard_True)
550 return Standard_True;
555 NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion());
556 if (isFirstCellFilterInvoke)
558 myFilter.Reset(myCellSize);
560 // Copy initial data into cell filter.
561 for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
563 math_Vector aVec(1, myN);
564 for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
565 aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
567 myFilter.Add(aVec, aVec);
571 isFirstCellFilterInvoke = Standard_False;
573 math_Vector aLow(1, myN), anUp(1, myN);
574 anInspector.Shift(thePnt, myCellSize, aLow, anUp);
576 anInspector.ClearFind();
577 anInspector.SetCurrent(thePnt);
578 myFilter.Inspect(aLow, anUp, anInspector);
579 if (!anInspector.isFind())
581 // Point is out of close cells, add new one.
582 myFilter.Add(thePnt, thePnt);
585 return Standard_False;
588 //=======================================================================
591 //=======================================================================
592 void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
596 for(j = 1; j <= myN; j++)
597 theSol(j) = myY((theIndex - 1) * myN + j);
600 //=======================================================================
601 //function : initCellSize
603 //=======================================================================
604 void math_GlobOptMin::initCellSize()
606 for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
608 myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
609 * Precision::PConfusion() / (2.0 * Sqrt(2.0));
613 //=======================================================================
614 //function : CheckFunctionalStopCriteria
616 //=======================================================================
617 Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
619 // Search single solution and current solution in its neighborhood.
620 if (myIsFindSingleSolution &&
621 Abs (myF - myFunctionalMinimalValue) < mySameTol * 0.01)
622 return Standard_True;
624 return Standard_False;
627 //=======================================================================
628 //function : ComputeInitSol
630 //=======================================================================
631 void math_GlobOptMin::ComputeInitSol()
633 Standard_Real aCurrVal, aBestVal;
634 math_Vector aCurrPnt(1, myN);
635 math_Vector aBestPnt(1, myN);
636 math_Vector aParamStep(1, myN);
637 // Check functional value in midpoint, lower and upper border points and
638 // in each point try to perform local optimization.
639 aBestPnt = (myGlobA + myGlobB) * 0.5;
640 myFunc->Value(aBestPnt, aBestVal);
643 for(i = 1; i <= 3; i++)
645 aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
647 if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
649 // Local search tries to find better solution than current point.
650 if (aCurrVal < aBestVal)
660 for(i = 1; i <= myN; i++)
661 myY.Append(aBestPnt(i));