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()),
73 myF(Precision::Infinite())
80 myIsFindSingleSolution = Standard_False;
81 myFunctionalMinimalValue = -Precision::Infinite();
85 for(i = 1; i <= myN; i++)
94 for(i = 1; i <= myN; i++)
96 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
99 myTol = theDiscretizationTol;
100 mySameTol = theSameTol;
102 const Standard_Integer aMaxSquareSearchSol = 200;
103 Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
104 myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
108 myDone = Standard_False;
111 //=======================================================================
112 //function : SetGlobalParams
113 //purpose : Set parameters without memory allocation.
114 //=======================================================================
115 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
116 const math_Vector& theA,
117 const math_Vector& theB,
118 const Standard_Real theC,
119 const Standard_Real theDiscretizationTol,
120 const Standard_Real theSameTol)
130 for(i = 1; i <= myN; i++)
132 myGlobA(i) = theA(i);
133 myGlobB(i) = theB(i);
139 for(i = 1; i <= myN; i++)
141 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
144 myTol = theDiscretizationTol;
145 mySameTol = theSameTol;
150 myDone = Standard_False;
153 //=======================================================================
154 //function : SetLocalParams
155 //purpose : Set parameters without memory allocation.
156 //=======================================================================
157 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
158 const math_Vector& theLocalB)
163 for(i = 1; i <= myN; i++)
165 myA(i) = theLocalA(i);
166 myB(i) = theLocalB(i);
169 for(i = 1; i <= myN; i++)
171 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
174 myDone = Standard_False;
177 //=======================================================================
179 //purpose : Set algorithm tolerances.
180 //=======================================================================
181 void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
182 const Standard_Real theSameTol)
184 myTol = theDiscretizationTol;
185 mySameTol = theSameTol;
188 //=======================================================================
190 //purpose : Get algorithm tolerances.
191 //=======================================================================
192 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
193 Standard_Real& theSameTol)
195 theDiscretizationTol = myTol;
196 theSameTol = mySameTol;
199 //=======================================================================
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)
206 myDone = Standard_False;
208 // Compute parameters range
209 Standard_Real minLength = RealLast();
210 Standard_Real maxLength = RealFirst();
211 for(Standard_Integer i = 1; i <= myN; i++)
213 Standard_Real currentLength = myB(i) - myA(i);
214 if (currentLength < minLength)
215 minLength = currentLength;
216 if (currentLength > maxLength)
217 maxLength = currentLength;
222 if (minLength < Precision::PConfusion())
225 std::cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << std::endl;
231 if (!myIsConstLocked)
233 // Compute initial value for myC.
234 computeInitialValues();
237 myE1 = minLength * myTol;
238 myE2 = maxLength * myTol;
240 myIsFindSingleSolution = isFindSingleSolution;
241 if (isFindSingleSolution)
243 // Run local optimization if current value better than optimal.
249 myE3 = - maxLength * myTol / 4.0;
251 myE3 = - maxLength * myTol * myC / 4.0;
254 // Search single solution and current solution in its neighborhood.
255 if (CheckFunctionalStopCriteria())
257 myDone = Standard_True;
262 isFirstCellFilterInvoke = Standard_True;
263 computeGlobalExtremum(myN);
265 myDone = Standard_True;
268 //=======================================================================
269 //function : computeLocalExtremum
271 //=======================================================================
272 Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
273 Standard_Real& theVal,
274 math_Vector& theOutPnt)
280 dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
282 math_MultipleVarFunctionWithHessian* aTmp =
283 dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
284 math_NewtonMinimum newtonMinimum(*aTmp);
285 newtonMinimum.SetBoundary(myGlobA, myGlobB);
286 newtonMinimum.Perform(*aTmp, thePnt);
288 if (newtonMinimum.IsDone())
290 newtonMinimum.Location(theOutPnt);
291 theVal = newtonMinimum.Minimum();
293 if (isInside(theOutPnt))
294 return Standard_True;
300 dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
302 math_MultipleVarFunctionWithGradient* aTmp =
303 dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
304 math_BFGS bfgs(aTmp->NbVariables());
305 bfgs.SetBoundary(myGlobA, myGlobB);
306 bfgs.Perform(*aTmp, thePnt);
310 bfgs.Location(theOutPnt);
311 theVal = bfgs.Minimum();
313 if (isInside(theOutPnt))
314 return Standard_True;
318 // Powell method used.
319 if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
321 math_Matrix m(1, myN, 1, myN, 0.0);
322 for(i = 1; i <= myN; i++)
325 math_Powell powell(*myFunc, 1e-10);
326 powell.Perform(*myFunc, thePnt, m);
330 powell.Location(theOutPnt);
331 theVal = powell.Minimum();
333 if (isInside(theOutPnt))
334 return Standard_True;
338 return Standard_False;
341 //=======================================================================
342 //function : computeInitialValues
344 //=======================================================================
345 void math_GlobOptMin::computeInitialValues()
347 const Standard_Real aMinLC = 0.01;
348 const Standard_Real aMaxLC = 1000.;
349 const Standard_Real aMinEps = 0.1;
350 const Standard_Real aMaxEps = 100.;
352 math_Vector aCurrPnt(1, myN);
353 math_Vector aBestPnt(1, myN);
354 math_Vector aParamStep(1, myN);
355 Standard_Real aCurrVal = RealLast();
357 // Lipchitz const approximation.
358 Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
359 Standard_Integer aPntNb = 13;
360 myFunc->Value(myA, aPrevValDiag);
361 aPrevValProj = aPrevValDiag;
362 Standard_Real aStep = (myB - myA).Norm() / aPntNb;
363 aParamStep = (myB - myA) / aPntNb;
364 for(i = 1; i <= aPntNb; i++)
366 aCurrPnt = myA + aParamStep * i;
368 // Walk over diagonal.
369 myFunc->Value(aCurrPnt, aCurrVal);
370 aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
371 aPrevValDiag = aCurrVal;
373 // Walk over diag in projected space aPnt(1) = myA(1) = const.
374 aCurrPnt(1) = myA(1);
375 myFunc->Value(aCurrPnt, aCurrVal);
376 aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
377 aPrevValProj = aCurrVal;
381 aLipConst *= Sqrt(myN) / aStep;
382 if (aLipConst < myC * aMinEps)
383 myC = Max(aLipConst * aMinEps, aMinLC);
384 else if (aLipConst > myC * aMaxEps)
385 myC = Min(myC * aMaxEps, aMaxLC);
388 //=======================================================================
389 //function : ComputeGlobalExtremum
391 //=======================================================================
392 void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
395 Standard_Real d = RealLast(), aPrevVal; // Functional in original and moved points.
396 Standard_Real val = RealLast(); // Local extrema computed in moved point.
397 Standard_Real aStepBestValue = RealLast();
398 math_Vector aStepBestPoint(1, myN);
399 Standard_Boolean isInside = Standard_False,
400 isReached = Standard_False;
402 Standard_Real r1, r2, r;
404 for(myX(j) = myA(j) + myE1; !isReached; myX(j) += myV(j))
409 isReached = Standard_True;
412 if (CheckFunctionalStopCriteria())
413 return; // Best possible value is obtained.
417 isInside = Standard_False;
419 myFunc->Value(myX, d);
420 r1 = (d + myZ * myC * myLastStep - myF) * myZ; // Evtushenko estimation.
421 r2 = ((d + aPrevVal - myC * myLastStep) * 0.5 - myF) * myZ; // Shubert / Piyavsky estimation.
425 Standard_Real aSaveParam = myX(1);
427 // Piyavsky midpoint estimation.
428 Standard_Real aParam = (2 * myX(1) - myV(1) ) * 0.5 + (aPrevVal - d) * 0.5 / myC;
429 if (Precision::IsInfinite(aPrevVal))
430 aParam = myX(1) - myV(1) * 0.5; // Protection from upper dimension step.
433 Standard_Real aVal = 0;
434 myFunc->Value(myX, aVal);
437 if ( (aVal < d && aVal < aPrevVal) ||
438 DistanceToBorder(myX, myA, myB) < myE1 ) // Condition optimization case near the border.
440 isInside = computeLocalExtremum(myX, val, myTmp);
443 aStepBestValue = (isInside && (val < d))? val : d;
444 aStepBestPoint = (isInside && (val < d))? myTmp : myX;
446 // Check point and value on the current step to be optimal.
447 checkAddCandidate(aStepBestPoint, aStepBestValue);
449 if (CheckFunctionalStopCriteria())
450 return; // Best possible value is obtained.
452 myV(1) = Min(myE2 + Abs(myF - d) / myC, myMaxV(1));
457 myV(j) = RealLast() / 2.0;
458 computeGlobalExtremum(j - 1);
460 // Nullify steps on lower dimensions.
461 for(i = 1; i < j; i++)
466 Standard_Real aUpperDimStep = Max(myV(j), myE2);
467 if (myV(j + 1) > aUpperDimStep)
469 if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
470 myV(j + 1) = myMaxV(j + 1);
472 myV(j + 1) = aUpperDimStep;
478 //=======================================================================
479 //function : IsInside
481 //=======================================================================
482 Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
486 for(i = 1; i <= myN; i++)
488 if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
489 return Standard_False;
492 return Standard_True;
494 //=======================================================================
495 //function : IsStored
497 //=======================================================================
498 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
500 Standard_Integer i,j;
501 Standard_Boolean isSame = Standard_True;
502 math_Vector aTol(1, myN);
503 aTol = (myB - myA) * mySameTol;
505 // C1 * n^2 = C2 * 3^dim * n
506 if (mySolCount < myMinCellFilterSol)
508 for(i = 0; i < mySolCount; i++)
510 isSame = Standard_True;
511 for(j = 1; j <= myN; j++)
513 if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
515 isSame = Standard_False;
519 if (isSame == Standard_True)
520 return Standard_True;
525 NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion());
526 if (isFirstCellFilterInvoke)
528 myFilter.Reset(myCellSize);
530 // Copy initial data into cell filter.
531 for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
533 math_Vector aVec(1, myN);
534 for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
535 aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
537 myFilter.Add(aVec, aVec);
541 isFirstCellFilterInvoke = Standard_False;
543 math_Vector aLow(1, myN), anUp(1, myN);
544 anInspector.Shift(thePnt, myCellSize, aLow, anUp);
546 anInspector.ClearFind();
547 anInspector.SetCurrent(thePnt);
548 myFilter.Inspect(aLow, anUp, anInspector);
549 if (!anInspector.isFind())
551 // Point is out of close cells, add new one.
552 myFilter.Add(thePnt, thePnt);
555 return Standard_False;
558 //=======================================================================
561 //=======================================================================
562 void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
566 for(j = 1; j <= myN; j++)
567 theSol(j) = myY((theIndex - 1) * myN + j);
570 //=======================================================================
571 //function : initCellSize
573 //=======================================================================
574 void math_GlobOptMin::initCellSize()
576 for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
578 myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
579 * Precision::PConfusion() / (2.0 * Sqrt(2.0));
583 //=======================================================================
584 //function : CheckFunctionalStopCriteria
586 //=======================================================================
587 Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
589 // Search single solution and current solution in its neighborhood.
590 if (myIsFindSingleSolution &&
591 Abs (myF - myFunctionalMinimalValue) < mySameTol * 0.01)
592 return Standard_True;
594 return Standard_False;
597 //=======================================================================
598 //function : ComputeInitSol
600 //=======================================================================
601 void math_GlobOptMin::ComputeInitSol()
604 math_Vector aPnt(1, myN);
606 // Check functional value in midpoint. It is necessary since local optimization
607 // algorithm may fail and return nothing. This is a protection from uninitialized
609 aPnt = (myGlobA + myGlobB) * 0.5;
610 myFunc->Value(aPnt, aVal);
611 checkAddCandidate(aPnt, aVal);
613 // Run local optimization from lower corner, midpoint, and upper corner.
614 for(Standard_Integer i = 1; i <= 3; i++)
616 aPnt = myA + (myB - myA) * (i - 1) / 2.0;
618 if(computeLocalExtremum(aPnt, aVal, aPnt))
619 checkAddCandidate(aPnt, aVal);
623 //=======================================================================
624 //function : checkAddCandidate
626 //=======================================================================
627 void math_GlobOptMin::checkAddCandidate(const math_Vector& thePnt,
628 const Standard_Real theValue)
630 if (Abs(theValue - myF) < mySameTol * 0.01 && // Value in point is close to optimal value.
631 !myIsFindSingleSolution) // Several optimal solutions are allowed.
633 if (!isStored(thePnt))
635 if ((theValue - myF) * myZ > 0.0)
637 for (Standard_Integer j = 1; j <= myN; j++)
638 myY.Append(thePnt(j));
643 // New best solution:
644 // new point is out of (mySameTol * 0.01) surrounding or
645 // new point is better than old and single point search.
646 Standard_Real aDelta = (theValue - myF) * myZ;
647 if (aDelta > mySameTol * 0.01 ||
648 (aDelta > 0.0 && myIsFindSingleSolution))
652 for (Standard_Integer j = 1; j <= myN; j++)
653 myY.Append(thePnt(j));
656 isFirstCellFilterInvoke = Standard_True;