1 // Created on: 2014-01-20
2 // Created by: Alexaner Malyshev
3 // Copyright (c) 2014-2014 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()),
54 myIsFindSingleSolution = Standard_False;
58 for(i = 1; i <= myN; i++)
67 for(i = 1; i <= myN; i++)
69 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
72 myExpandCoeff(1) = 1.0;
73 for(i = 2; i <= myN; i++)
75 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
78 myTol = theDiscretizationTol;
79 mySameTol = theSameTol;
81 myDone = Standard_False;
84 //=======================================================================
85 //function : SetGlobalParams
86 //purpose : Set params without memory allocation.
87 //=======================================================================
88 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
89 const math_Vector& theA,
90 const math_Vector& theB,
91 const Standard_Real theC,
92 const Standard_Real theDiscretizationTol,
93 const Standard_Real theSameTol)
102 for(i = 1; i <= myN; i++)
104 myGlobA(i) = theA(i);
105 myGlobB(i) = theB(i);
111 for(i = 1; i <= myN; i++)
113 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
116 myExpandCoeff(1) = 1.0;
117 for(i = 2; i <= myN; i++)
119 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
122 myTol = theDiscretizationTol;
123 mySameTol = theSameTol;
125 myDone = Standard_False;
128 //=======================================================================
129 //function : SetLocalParams
130 //purpose : Set params without memory allocation.
131 //=======================================================================
132 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
133 const math_Vector& theLocalB)
140 for(i = 1; i <= myN; i++)
142 myA(i) = theLocalA(i);
143 myB(i) = theLocalB(i);
146 for(i = 1; i <= myN; i++)
148 myMaxV(i) = (myB(i) - myA(i)) / 3.0;
151 myExpandCoeff(1) = 1.0;
152 for(i = 2; i <= myN; i++)
154 myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
157 myDone = Standard_False;
160 //=======================================================================
162 //purpose : Set algorithm tolerances.
163 //=======================================================================
164 void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
165 const Standard_Real theSameTol)
167 myTol = theDiscretizationTol;
168 mySameTol = theSameTol;
171 //=======================================================================
173 //purpose : Get algorithm tolerances.
174 //=======================================================================
175 void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
176 Standard_Real& theSameTol)
178 theDiscretizationTol = myTol;
179 theSameTol = mySameTol;
182 //=======================================================================
183 //function : ~math_GlobOptMin
185 //=======================================================================
186 math_GlobOptMin::~math_GlobOptMin()
190 //=======================================================================
192 //purpose : Compute Global extremum point
193 //=======================================================================
194 // In this algo indexes started from 1, not from 0.
195 void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
199 // Compute parameters range
200 Standard_Real minLength = RealLast();
201 Standard_Real maxLength = RealFirst();
202 for(i = 1; i <= myN; i++)
204 Standard_Real currentLength = myB(i) - myA(i);
205 if (currentLength < minLength)
206 minLength = currentLength;
207 if (currentLength > maxLength)
208 maxLength = currentLength;
211 if (minLength < Precision::PConfusion())
214 cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << endl;
220 // Compute initial values for myF, myY, myC.
221 computeInitialValues();
223 myE1 = minLength * myTol;
224 myE2 = maxLength * myTol;
226 myIsFindSingleSolution = isFindSingleSolution;
227 if (isFindSingleSolution)
229 // Run local optimization
230 // if current value better than optimal.
236 myE3 = - maxLength * myTol / 4.0;
238 myE3 = - maxLength * myTol * myC / 4.0;
241 computeGlobalExtremum(myN);
243 myDone = Standard_True;
246 //=======================================================================
247 //function : computeLocalExtremum
249 //=======================================================================
250 Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
251 Standard_Real& theVal,
252 math_Vector& theOutPnt)
257 if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
259 math_MultipleVarFunctionWithHessian* myTmp =
260 dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
261 math_NewtonMinimum newtonMinimum(*myTmp);
262 newtonMinimum.SetBoundary(myGlobA, myGlobB);
263 newtonMinimum.Perform(*myTmp, thePnt);
265 if (newtonMinimum.IsDone())
267 newtonMinimum.Location(theOutPnt);
268 theVal = newtonMinimum.Minimum();
270 else return Standard_False;
274 if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
276 math_MultipleVarFunctionWithGradient* myTmp =
277 dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
278 math_BFGS bfgs(myTmp->NbVariables());
279 bfgs.Perform(*myTmp, thePnt);
282 bfgs.Location(theOutPnt);
283 theVal = bfgs.Minimum();
285 else return Standard_False;
288 // Powell method used.
289 if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
291 math_Matrix m(1, myN, 1, myN, 0.0);
292 for(i = 1; i <= myN; i++)
295 math_Powell powell(*myFunc, 1e-10);
296 powell.Perform(*myFunc, thePnt, m);
300 powell.Location(theOutPnt);
301 theVal = powell.Minimum();
303 else return Standard_False;
306 if (isInside(theOutPnt))
307 return Standard_True;
309 return Standard_False;
312 //=======================================================================
313 //function : computeInitialValues
315 //=======================================================================
316 void math_GlobOptMin::computeInitialValues()
319 math_Vector aCurrPnt(1, myN);
320 math_Vector aBestPnt(1, myN);
321 math_Vector aParamStep(1, myN);
322 Standard_Real aCurrVal = RealLast();
323 Standard_Real aBestVal = RealLast();
325 // Check functional value in midpoint, low and upp point border and
326 // in each point try to perform local optimization.
327 aBestPnt = (myA + myB) * 0.5;
328 myFunc->Value(aBestPnt, aBestVal);
330 for(i = 1; i <= 3; i++)
332 aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
334 if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
336 // Local Extremum finds better solution than current point.
337 if (aCurrVal < aBestVal)
347 for(i = 1; i <= myN; i++)
348 myY.Append(aBestPnt(i));
351 // Lipschitz const approximation
352 Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
353 Standard_Integer aPntNb = 13;
354 myFunc->Value(myA, aPrevValDiag);
355 aPrevValProj = aPrevValDiag;
356 Standard_Real aStep = (myB - myA).Norm() / aPntNb;
357 aParamStep = (myB - myA) / aPntNb;
358 for(i = 1; i <= aPntNb; i++)
360 aCurrPnt = myA + aParamStep * i;
362 // Walk over diagonal.
363 myFunc->Value(aCurrPnt, aCurrVal);
364 aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
365 aPrevValDiag = aCurrVal;
367 // Walk over diag in projected space aPnt(1) = myA(1) = const.
368 aCurrPnt(1) = myA(1);
369 myFunc->Value(aCurrPnt, aCurrVal);
370 aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
371 aPrevValProj = aCurrVal;
374 aLipConst *= Sqrt(myN) / aStep;
376 if (aLipConst < myC * 0.1)
378 myC = Max(aLipConst * 0.1, 0.01);
380 else if (aLipConst > myC * 10)
382 myC = Min(myC * 2, 30.0);
386 //=======================================================================
387 //function : ComputeGlobalExtremum
389 //=======================================================================
390 void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
393 Standard_Real d; // Functional in moved point.
394 Standard_Real val = RealLast(); // Local extrema computed in moved point.
395 Standard_Real aStepBestValue = RealLast();
396 Standard_Real aRealStep = 0.0;
397 math_Vector aStepBestPoint(1, myN);
398 Standard_Boolean isInside = Standard_False;
402 for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j))
409 isInside = Standard_False;
410 myFunc->Value(myX, d);
411 r = (d + myZ * myC * aRealStep - myF) * myZ;
414 isInside = computeLocalExtremum(myX, val, myTmp);
416 aStepBestValue = (isInside && (val < d))? val : d;
417 aStepBestPoint = (isInside && (val < d))? myTmp : myX;
419 // Solutions are close to each other
420 // and it is allowed to have more than one solution.
421 if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
422 !myIsFindSingleSolution)
424 if (!isStored(aStepBestPoint))
426 if ((aStepBestValue - myF) * myZ > 0.0)
427 myF = aStepBestValue;
428 for(i = 1; i <= myN; i++)
429 myY.Append(aStepBestPoint(i));
434 // New best solution:
435 // new point is out of (mySameTol * 0.01) surrounding or
436 // new point is better than old + single point search.
437 Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
438 if (aFunctionalDelta > mySameTol * 0.01 ||
439 (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
442 myF = aStepBestValue;
444 for(i = 1; i <= myN; i++)
445 myY.Append(aStepBestPoint(i));
449 aRealStep = myE2 + Abs(myF - d) / myC;
450 myV(1) = Min(aRealStep, myMaxV(1));
454 myV(j) = RealLast() / 2.0;
455 computeGlobalExtremum(j - 1);
457 // Nullify steps on lower dimensions.
458 for(i = 1; i < j; i++)
461 // Compute step in (j + 1) dimension according to scale.
464 Standard_Real aUpperDimStep = myV(j) * myExpandCoeff(j + 1);
465 if (myV(j + 1) > aUpperDimStep)
467 if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
468 myV(j + 1) = myMaxV(j + 1);
470 myV(j + 1) = aUpperDimStep;
476 //=======================================================================
477 //function : IsInside
479 //=======================================================================
480 Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
484 for(i = 1; i <= myN; i++)
486 if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
487 return Standard_False;
490 return Standard_True;
492 //=======================================================================
493 //function : IsStored
495 //=======================================================================
496 Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
498 Standard_Integer i,j;
499 Standard_Boolean isSame = Standard_True;
501 for(i = 0; i < mySolCount; i++)
503 isSame = Standard_True;
504 for(j = 1; j <= myN; j++)
506 if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) - myA(j)) * mySameTol)
508 isSame = Standard_False;
512 if (isSame == Standard_True)
513 return Standard_True;
516 return Standard_False;
519 //=======================================================================
520 //function : NbExtrema
522 //=======================================================================
523 Standard_Integer math_GlobOptMin::NbExtrema()
528 //=======================================================================
531 //=======================================================================
532 Standard_Real math_GlobOptMin::GetF()
537 //=======================================================================
540 //=======================================================================
541 Standard_Boolean math_GlobOptMin::isDone()
546 //=======================================================================
549 //=======================================================================
550 void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
554 for(j = 1; j <= myN; j++)
555 theSol(j) = myY((theIndex - 1) * myN + j);