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()),
48 myExpandCoeff(1, myN),
49 myCellSize(0, myN - 1),
50 myFilter(theFunc->NbVariables())
56 myIsFindSingleSolution = Standard_False;
60 for(i = 1; i <= myN; i++)
69 for(i = 1; i <= myN; i++)
71 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
74 myExpandCoeff(1) = 1.0;
75 for(i = 2; i <= myN; i++)
77 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
80 myTol = theDiscretizationTol;
81 mySameTol = theSameTol;
83 const Standard_Integer aMaxSquareSearchSol = 200;
84 Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
85 myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
88 myDone = Standard_False;
91 //=======================================================================
92 //function : SetGlobalParams
93 //purpose : Set params without memory allocation.
94 //=======================================================================
95 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
96 const math_Vector& theA,
97 const math_Vector& theB,
98 const Standard_Real theC,
99 const Standard_Real theDiscretizationTol,
100 const Standard_Real theSameTol)
109 for(i = 1; i <= myN; i++)
111 myGlobA(i) = theA(i);
112 myGlobB(i) = theB(i);
118 for(i = 1; i <= myN; i++)
120 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
123 myExpandCoeff(1) = 1.0;
124 for(i = 2; i <= myN; i++)
126 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
129 myTol = theDiscretizationTol;
130 mySameTol = theSameTol;
134 myDone = Standard_False;
137 //=======================================================================
138 //function : SetLocalParams
139 //purpose : Set params without memory allocation.
140 //=======================================================================
141 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
142 const math_Vector& theLocalB)
149 for(i = 1; i <= myN; i++)
151 myA(i) = theLocalA(i);
152 myB(i) = theLocalB(i);
155 for(i = 1; i <= myN; i++)
157 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
160 myExpandCoeff(1) = 1.0;
161 for(i = 2; i <= myN; i++)
163 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
166 myDone = Standard_False;
169 //=======================================================================
171 //purpose : Set algorithm tolerances.
172 //=======================================================================
173 void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
174 const Standard_Real theSameTol)
176 myTol = theDiscretizationTol;
177 mySameTol = theSameTol;
180 //=======================================================================
182 //purpose : Get algorithm tolerances.
183 //=======================================================================
184 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
185 Standard_Real& theSameTol)
187 theDiscretizationTol = myTol;
188 theSameTol = mySameTol;
191 //=======================================================================
192 //function : ~math_GlobOptMin
194 //=======================================================================
195 math_GlobOptMin::~math_GlobOptMin()
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)
208 // Compute parameters range
209 Standard_Real minLength = RealLast();
210 Standard_Real maxLength = RealFirst();
211 for(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;
220 if (minLength < Precision::PConfusion())
223 cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << endl;
229 // Compute initial values for myF, myY, myC.
230 computeInitialValues();
232 myE1 = minLength * myTol;
233 myE2 = maxLength * myTol;
235 myIsFindSingleSolution = isFindSingleSolution;
236 if (isFindSingleSolution)
238 // Run local optimization
239 // if current value better than optimal.
245 myE3 = - maxLength * myTol / 4.0;
247 myE3 = - maxLength * myTol * myC / 4.0;
250 isFirstCellFilterInvoke = Standard_True;
251 computeGlobalExtremum(myN);
253 myDone = Standard_True;
256 //=======================================================================
257 //function : computeLocalExtremum
259 //=======================================================================
260 Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
261 Standard_Real& theVal,
262 math_Vector& theOutPnt)
267 if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
269 math_MultipleVarFunctionWithHessian* myTmp =
270 dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
271 math_NewtonMinimum newtonMinimum(*myTmp);
272 newtonMinimum.SetBoundary(myGlobA, myGlobB);
273 newtonMinimum.Perform(*myTmp, thePnt);
275 if (newtonMinimum.IsDone())
277 newtonMinimum.Location(theOutPnt);
278 theVal = newtonMinimum.Minimum();
280 else return Standard_False;
284 if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
286 math_MultipleVarFunctionWithGradient* myTmp =
287 dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
288 math_BFGS bfgs(myTmp->NbVariables());
289 bfgs.Perform(*myTmp, thePnt);
292 bfgs.Location(theOutPnt);
293 theVal = bfgs.Minimum();
295 else return Standard_False;
298 // Powell method used.
299 if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
301 math_Matrix m(1, myN, 1, myN, 0.0);
302 for(i = 1; i <= myN; i++)
305 math_Powell powell(*myFunc, 1e-10);
306 powell.Perform(*myFunc, thePnt, m);
310 powell.Location(theOutPnt);
311 theVal = powell.Minimum();
313 else return Standard_False;
316 if (isInside(theOutPnt))
317 return Standard_True;
319 return Standard_False;
322 //=======================================================================
323 //function : computeInitialValues
325 //=======================================================================
326 void math_GlobOptMin::computeInitialValues()
329 math_Vector aCurrPnt(1, myN);
330 math_Vector aBestPnt(1, myN);
331 math_Vector aParamStep(1, myN);
332 Standard_Real aCurrVal = RealLast();
333 Standard_Real aBestVal = RealLast();
335 // Check functional value in midpoint, low and upp point border and
336 // in each point try to perform local optimization.
337 aBestPnt = (myA + myB) * 0.5;
338 myFunc->Value(aBestPnt, aBestVal);
340 for(i = 1; i <= 3; i++)
342 aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
344 if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
346 // Local Extremum finds better solution than current point.
347 if (aCurrVal < aBestVal)
357 for(i = 1; i <= myN; i++)
358 myY.Append(aBestPnt(i));
361 // Lipschitz const approximation
362 Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
363 Standard_Integer aPntNb = 13;
364 myFunc->Value(myA, aPrevValDiag);
365 aPrevValProj = aPrevValDiag;
366 Standard_Real aStep = (myB - myA).Norm() / aPntNb;
367 aParamStep = (myB - myA) / aPntNb;
368 for(i = 1; i <= aPntNb; i++)
370 aCurrPnt = myA + aParamStep * i;
372 // Walk over diagonal.
373 myFunc->Value(aCurrPnt, aCurrVal);
374 aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
375 aPrevValDiag = aCurrVal;
377 // Walk over diag in projected space aPnt(1) = myA(1) = const.
378 aCurrPnt(1) = myA(1);
379 myFunc->Value(aCurrPnt, aCurrVal);
380 aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
381 aPrevValProj = aCurrVal;
384 aLipConst *= Sqrt(myN) / aStep;
386 if (aLipConst < myC * 0.1)
388 myC = Max(aLipConst * 0.1, 0.01);
390 else if (aLipConst > myC * 10)
392 myC = Min(myC * 2, 30.0);
396 //=======================================================================
397 //function : ComputeGlobalExtremum
399 //=======================================================================
400 void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
403 Standard_Real d; // Functional in moved point.
404 Standard_Real val = RealLast(); // Local extrema computed in moved point.
405 Standard_Real aStepBestValue = RealLast();
406 Standard_Real aRealStep = 0.0;
407 math_Vector aStepBestPoint(1, myN);
408 Standard_Boolean isInside = Standard_False;
411 for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j))
418 isInside = Standard_False;
419 myFunc->Value(myX, d);
420 r = (d + myZ * myC * aRealStep - myF) * myZ;
423 isInside = computeLocalExtremum(myX, val, myTmp);
425 aStepBestValue = (isInside && (val < d))? val : d;
426 aStepBestPoint = (isInside && (val < d))? myTmp : myX;
428 // Solutions are close to each other
429 // and it is allowed to have more than one solution.
430 if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
431 !myIsFindSingleSolution)
433 if (!isStored(aStepBestPoint))
435 if ((aStepBestValue - myF) * myZ > 0.0)
436 myF = aStepBestValue;
437 for(i = 1; i <= myN; i++)
438 myY.Append(aStepBestPoint(i));
443 // New best solution:
444 // new point is out of (mySameTol * 0.01) surrounding or
445 // new point is better than old + single point search.
446 Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
447 if (aFunctionalDelta > mySameTol * 0.01 ||
448 (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
451 myF = aStepBestValue;
453 for(i = 1; i <= myN; i++)
454 myY.Append(aStepBestPoint(i));
457 isFirstCellFilterInvoke = Standard_True;
460 aRealStep = myE2 + Abs(myF - d) / myC;
461 myV(1) = Min(aRealStep, myMaxV(1));
465 myV(j) = RealLast() / 2.0;
466 computeGlobalExtremum(j - 1);
468 // Nullify steps on lower dimensions.
469 for(i = 1; i < j; i++)
472 // Compute step in (j + 1) dimension according to scale.
475 Standard_Real aUpperDimStep = myV(j) * myExpandCoeff(j + 1);
476 if (myV(j + 1) > aUpperDimStep)
478 if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
479 myV(j + 1) = myMaxV(j + 1);
481 myV(j + 1) = aUpperDimStep;
487 //=======================================================================
488 //function : IsInside
490 //=======================================================================
491 Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
495 for(i = 1; i <= myN; i++)
497 if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
498 return Standard_False;
501 return Standard_True;
503 //=======================================================================
504 //function : IsStored
506 //=======================================================================
507 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
509 Standard_Integer i,j;
510 Standard_Boolean isSame = Standard_True;
511 math_Vector aTol(1, myN);
512 aTol = (myB - myA) * mySameTol;
514 // C1 * n^2 = C2 * 3^dim * n
515 if (mySolCount < myMinCellFilterSol)
517 for(i = 0; i < mySolCount; i++)
519 isSame = Standard_True;
520 for(j = 1; j <= myN; j++)
522 if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
524 isSame = Standard_False;
528 if (isSame == Standard_True)
529 return Standard_True;
534 NCollection_CellFilter_NDimInspector anInspector(myN, Precision::PConfusion());
535 if (isFirstCellFilterInvoke)
537 myFilter.Reset(myCellSize);
539 // Copy initial data into cell filter.
540 for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
542 math_Vector aVec(1, myN);
543 for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
544 aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
546 myFilter.Add(aVec, aVec);
550 isFirstCellFilterInvoke = Standard_False;
552 math_Vector aLow(1, myN), anUp(1, myN);
553 anInspector.Shift(thePnt, myCellSize, aLow, anUp);
555 anInspector.ClearFind();
556 anInspector.SetCurrent(thePnt);
557 myFilter.Inspect(aLow, anUp, anInspector);
558 if (!anInspector.isFind())
560 // Point is out of close cells, add new one.
561 myFilter.Add(thePnt, thePnt);
564 return Standard_False;
567 //=======================================================================
568 //function : NbExtrema
570 //=======================================================================
571 Standard_Integer math_GlobOptMin::NbExtrema()
576 //=======================================================================
579 //=======================================================================
580 Standard_Real math_GlobOptMin::GetF()
585 //=======================================================================
588 //=======================================================================
589 Standard_Boolean math_GlobOptMin::isDone()
594 //=======================================================================
597 //=======================================================================
598 void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
602 for(j = 1; j <= myN; j++)
603 theSol(j) = myY((theIndex - 1) * myN + j);
606 //=======================================================================
607 //function : initCellSize
609 //=======================================================================
610 void math_GlobOptMin::initCellSize()
612 for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
614 myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
615 * Precision::PConfusion() / (2.0 * Sqrt(2.0));