// Created on: 2014-01-20
// Created by: Alexaner Malyshev
-// Copyright (c) 2014-2014 OPEN CASCADE SAS
+// Copyright (c) 2014-2015 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
#include <math_Powell.hxx>
#include <Standard_Integer.hxx>
#include <Standard_Real.hxx>
+#include <Precision.hxx>
-const Handle(Standard_Type)& STANDARD_TYPE(math_GlobOptMin)
+//=======================================================================
+//function : DistanceToBorder
+//purpose :
+//=======================================================================
+static Standard_Real DistanceToBorder(const math_Vector & theX,
+ const math_Vector & theMin,
+ const math_Vector & theMax)
{
- static Handle(Standard_Type) _atype = new Standard_Type ("math_GlobOptMin", sizeof (math_GlobOptMin));
- return _atype;
+ Standard_Real aDist = RealLast();
+
+ for (Standard_Integer anIdx = theMin.Lower(); anIdx <= theMin.Upper(); ++anIdx)
+ {
+ const Standard_Real aDist1 = Abs (theX(anIdx) - theMin(anIdx));
+ const Standard_Real aDist2 = Abs (theX(anIdx) - theMax(anIdx));
+
+ aDist = Min (aDist, Min (aDist1, aDist2));
+ }
+
+ return aDist;
}
+
//=======================================================================
//function : math_GlobOptMin
//purpose : Constructor
myB(1, myN),
myGlobA(1, myN),
myGlobB(1, myN),
+ myIsConstLocked(Standard_False),
myX(1, myN),
myTmp(1, myN),
myV(1, myN),
- myMaxV(1, myN)
+ myMaxV(1, myN),
+ myCellSize(0, myN - 1),
+ myFilter(theFunc->NbVariables()),
+ myCont(2),
+ myF(Precision::Infinite())
{
Standard_Integer i;
myFunc = theFunc;
myC = theC;
+ myInitC = theC;
+ myIsFindSingleSolution = Standard_False;
+ myFunctionalMinimalValue = -Precision::Infinite();
myZ = -1;
mySolCount = 0;
myTol = theDiscretizationTol;
mySameTol = theSameTol;
+ const Standard_Integer aMaxSquareSearchSol = 200;
+ Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
+ myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
+ initCellSize();
+ ComputeInitSol();
+
myDone = Standard_False;
}
//=======================================================================
//function : SetGlobalParams
-//purpose : Set params without memory allocation.
+//purpose : Set parameters without memory allocation.
//=======================================================================
void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
const math_Vector& theA,
myFunc = theFunc;
myC = theC;
+ myInitC = theC;
myZ = -1;
mySolCount = 0;
myB(i) = theB(i);
}
+ for(i = 1; i <= myN; i++)
+ {
+ myMaxV(i) = (myB(i) - myA(i)) / 3.0;
+ }
+
myTol = theDiscretizationTol;
mySameTol = theSameTol;
+ initCellSize();
+ ComputeInitSol();
+
myDone = Standard_False;
}
//=======================================================================
//function : SetLocalParams
-//purpose : Set params without memory allocation.
+//purpose : Set parameters without memory allocation.
//=======================================================================
void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
const math_Vector& theLocalB)
Standard_Integer i;
myZ = -1;
- mySolCount = 0;
-
for(i = 1; i <= myN; i++)
{
myA(i) = theLocalA(i);
theSameTol = mySameTol;
}
-//=======================================================================
-//function : ~math_GlobOptMin
-//purpose :
-//=======================================================================
-math_GlobOptMin::~math_GlobOptMin()
-{
-}
-
//=======================================================================
//function : Perform
//purpose : Compute Global extremum point
//=======================================================================
// In this algo indexes started from 1, not from 0.
-void math_GlobOptMin::Perform()
+void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
{
- Standard_Integer i;
+ myDone = Standard_False;
// Compute parameters range
Standard_Real minLength = RealLast();
Standard_Real maxLength = RealFirst();
- for(i = 1; i <= myN; i++)
+ for(Standard_Integer i = 1; i <= myN; i++)
{
Standard_Real currentLength = myB(i) - myA(i);
if (currentLength < minLength)
minLength = currentLength;
if (currentLength > maxLength)
maxLength = currentLength;
+
+ myV(i) = 0.0;
}
- myE1 = minLength * myTol / myC;
- myE2 = maxLength * myTol * 2.0 / myC;
- myE3 = - maxLength * myTol / 4.0;
+ if (minLength < Precision::PConfusion())
+ {
+ #ifdef OCCT_DEBUG
+ std::cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << std::endl;
+ #endif
- // Compute start point.
- math_Vector aPnt(1,myN);
- for(i = 1; i <= myN; i++)
+ return;
+ }
+
+ if (!myIsConstLocked)
{
- Standard_Real currCentral = (myA(i) + myB(i)) / 2.0;
- aPnt(i) = currCentral;
+ // Compute initial value for myC.
+ computeInitialValues();
}
- myFunc->Value(aPnt, myF);
+ myE1 = minLength * myTol;
+ myE2 = maxLength * myTol;
- math_Vector aExtremumPoint(1,myN);
- Standard_Real aExtremumValue = RealLast();
- if (computeLocalExtremum(aPnt, aExtremumValue, aExtremumPoint))
+ myIsFindSingleSolution = isFindSingleSolution;
+ if (isFindSingleSolution)
{
- // Local Extremum finds better solution than midpoint.
- if (aExtremumValue < myF)
- {
- myF = aExtremumValue;
- aPnt = aExtremumPoint;
- }
+ // Run local optimization if current value better than optimal.
+ myE3 = 0.0;
+ }
+ else
+ {
+ if (myC > 1.0)
+ myE3 = - maxLength * myTol / 4.0;
+ else
+ myE3 = - maxLength * myTol * myC / 4.0;
}
- myY.Clear();
- for(i = 1; i <= myN; i++)
- myY.Append(aPnt(i));
- mySolCount++;
+ // Search single solution and current solution in its neighborhood.
+ if (CheckFunctionalStopCriteria())
+ {
+ myDone = Standard_True;
+ return;
+ }
+ myLastStep = 0.0;
+ isFirstCellFilterInvoke = Standard_True;
computeGlobalExtremum(myN);
myDone = Standard_True;
Standard_Integer i;
//Newton method
- if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
+ if (myCont >= 2 &&
+ dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
{
- math_MultipleVarFunctionWithHessian* myTmp =
+ math_MultipleVarFunctionWithHessian* aTmp =
dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
-
- math_NewtonMinimum newtonMinimum(*myTmp, thePnt);
+ math_NewtonMinimum newtonMinimum(*aTmp);
+ newtonMinimum.SetBoundary(myGlobA, myGlobB);
+ newtonMinimum.Perform(*aTmp, thePnt);
+
if (newtonMinimum.IsDone())
{
newtonMinimum.Location(theOutPnt);
theVal = newtonMinimum.Minimum();
+
+ if (isInside(theOutPnt))
+ return Standard_True;
}
- else return Standard_False;
- } else
+ }
// BFGS method used.
- if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
+ if (myCont >= 1 &&
+ dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
{
- math_MultipleVarFunctionWithGradient* myTmp =
+ math_MultipleVarFunctionWithGradient* aTmp =
dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
- math_BFGS bfgs(*myTmp, thePnt);
+ math_BFGS bfgs(aTmp->NbVariables());
+ bfgs.SetBoundary(myGlobA, myGlobB);
+ bfgs.Perform(*aTmp, thePnt);
+
if (bfgs.IsDone())
{
bfgs.Location(theOutPnt);
theVal = bfgs.Minimum();
+
+ if (isInside(theOutPnt))
+ return Standard_True;
}
- else return Standard_False;
- } else
+ }
// Powell method used.
if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
{
math_Matrix m(1, myN, 1, myN, 0.0);
for(i = 1; i <= myN; i++)
- m(1, 1) = 1.0;
+ m(i, i) = 1.0;
- math_Powell powell(*myFunc, thePnt, m, 1e-10);
+ math_Powell powell(*myFunc, 1e-10);
+ powell.Perform(*myFunc, thePnt, m);
if (powell.IsDone())
{
powell.Location(theOutPnt);
theVal = powell.Minimum();
+
+ if (isInside(theOutPnt))
+ return Standard_True;
}
- else return Standard_False;
}
- if (isInside(theOutPnt))
- return Standard_True;
- else
- return Standard_False;
+ return Standard_False;
+}
+
+//=======================================================================
+//function : computeInitialValues
+//purpose :
+//=======================================================================
+void math_GlobOptMin::computeInitialValues()
+{
+ const Standard_Real aMinLC = 0.01;
+ const Standard_Real aMaxLC = 1000.;
+ const Standard_Real aMinEps = 0.1;
+ const Standard_Real aMaxEps = 100.;
+ Standard_Integer i;
+ math_Vector aCurrPnt(1, myN);
+ math_Vector aBestPnt(1, myN);
+ math_Vector aParamStep(1, myN);
+ Standard_Real aCurrVal = RealLast();
+
+ // Lipchitz const approximation.
+ Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
+ Standard_Integer aPntNb = 13;
+ myFunc->Value(myA, aPrevValDiag);
+ aPrevValProj = aPrevValDiag;
+ Standard_Real aStep = (myB - myA).Norm() / aPntNb;
+ aParamStep = (myB - myA) / aPntNb;
+ for(i = 1; i <= aPntNb; i++)
+ {
+ aCurrPnt = myA + aParamStep * i;
+
+ // Walk over diagonal.
+ myFunc->Value(aCurrPnt, aCurrVal);
+ aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
+ aPrevValDiag = aCurrVal;
+
+ // Walk over diag in projected space aPnt(1) = myA(1) = const.
+ aCurrPnt(1) = myA(1);
+ myFunc->Value(aCurrPnt, aCurrVal);
+ aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
+ aPrevValProj = aCurrVal;
+ }
+
+ myC = myInitC;
+ aLipConst *= Sqrt(myN) / aStep;
+ if (aLipConst < myC * aMinEps)
+ myC = Max(aLipConst * aMinEps, aMinLC);
+ else if (aLipConst > myC * aMaxEps)
+ myC = Min(myC * aMaxEps, aMaxLC);
}
//=======================================================================
void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
{
Standard_Integer i;
- Standard_Real d; // Functional in moved point.
+ Standard_Real d = RealLast(), aPrevVal; // Functional in original and moved points.
Standard_Real val = RealLast(); // Local extrema computed in moved point.
- Standard_Real stepBestValue = RealLast();
- Standard_Real realStep = RealLast();
- math_Vector stepBestPoint(1, myN);
- Standard_Boolean isInside = Standard_False;
- Standard_Real r;
+ Standard_Real aStepBestValue = RealLast();
+ math_Vector aStepBestPoint(1, myN);
+ Standard_Boolean isInside = Standard_False,
+ isReached = Standard_False;
+ Standard_Real r1, r2, r;
- for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j))
+ for(myX(j) = myA(j) + myE1; !isReached; myX(j) += myV(j))
{
if (myX(j) > myB(j))
+ {
myX(j) = myB(j);
+ isReached = Standard_True;
+ }
+
+ if (CheckFunctionalStopCriteria())
+ return; // Best possible value is obtained.
if (j == 1)
{
isInside = Standard_False;
+ aPrevVal = d;
myFunc->Value(myX, d);
- r = (d - myF) * myZ;
+ r1 = (d + myZ * myC * myLastStep - myF) * myZ; // Evtushenko estimation.
+ r2 = ((d + aPrevVal - myC * myLastStep) * 0.5 - myF) * myZ; // Shubert / Piyavsky estimation.
+ r = Min(r1, r2);
if(r > myE3)
{
- isInside = computeLocalExtremum(myX, val, myTmp);
- }
- stepBestValue = (isInside && (val < d))? val : d;
- stepBestPoint = (isInside && (val < d))? myTmp : myX;
+ Standard_Real aSaveParam = myX(1);
- // Solutions are close to each other.
- if (Abs(stepBestValue - myF) < mySameTol * 0.01)
- {
- if (!isStored(stepBestPoint))
+ // Piyavsky midpoint estimation.
+ Standard_Real aParam = (2 * myX(1) - myV(1) ) * 0.5 + (aPrevVal - d) * 0.5 / myC;
+ if (Precision::IsInfinite(aPrevVal))
+ aParam = myX(1) - myV(1) * 0.5; // Protection from upper dimension step.
+
+ myX(1) = aParam;
+ Standard_Real aVal = 0;
+ myFunc->Value(myX, aVal);
+ myX(1) = aSaveParam;
+
+ if ( (aVal < d && aVal < aPrevVal) ||
+ DistanceToBorder(myX, myA, myB) < myE1 ) // Condition optimization case near the border.
{
- if ((stepBestValue - myF) * myZ > 0.0)
- myF = stepBestValue;
- for(i = 1; i <= myN; i++)
- myY.Append(stepBestPoint(i));
- mySolCount++;
+ isInside = computeLocalExtremum(myX, val, myTmp);
}
}
+ aStepBestValue = (isInside && (val < d))? val : d;
+ aStepBestPoint = (isInside && (val < d))? myTmp : myX;
- // New best solution.
- if ((stepBestValue - myF) * myZ > mySameTol * 0.01)
- {
- mySolCount = 0;
- myF = stepBestValue;
- myY.Clear();
- for(i = 1; i <= myN; i++)
- myY.Append(stepBestPoint(i));
- mySolCount++;
- }
+ // Check point and value on the current step to be optimal.
+ checkAddCandidate(aStepBestPoint, aStepBestValue);
+
+ if (CheckFunctionalStopCriteria())
+ return; // Best possible value is obtained.
- realStep = myE2 + Abs(myF - d) / myC;
- myV(1) = Min(realStep, myMaxV(1));
+ myV(1) = Min(myE2 + Abs(myF - d) / myC, myMaxV(1));
+ myLastStep = myV(1);
}
else
{
myV(j) = RealLast() / 2.0;
computeGlobalExtremum(j - 1);
+
+ // Nullify steps on lower dimensions.
+ for(i = 1; i < j; i++)
+ myV(i) = 0.0;
}
- if ((j < myN) && (myV(j + 1) > realStep))
+ if (j < myN)
{
- if (realStep > myMaxV(j + 1)) // Case of too big step.
- myV(j + 1) = myMaxV(j + 1);
- else
- myV(j + 1) = realStep;
+ Standard_Real aUpperDimStep = Max(myV(j), myE2);
+ if (myV(j + 1) > aUpperDimStep)
+ {
+ if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
+ myV(j + 1) = myMaxV(j + 1);
+ else
+ myV(j + 1) = aUpperDimStep;
+ }
}
}
}
{
Standard_Integer i,j;
Standard_Boolean isSame = Standard_True;
+ math_Vector aTol(1, myN);
+ aTol = (myB - myA) * mySameTol;
- for(i = 0; i < mySolCount; i++)
+ // C1 * n^2 = C2 * 3^dim * n
+ if (mySolCount < myMinCellFilterSol)
+ {
+ for(i = 0; i < mySolCount; i++)
+ {
+ isSame = Standard_True;
+ for(j = 1; j <= myN; j++)
+ {
+ if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
+ {
+ isSame = Standard_False;
+ break;
+ }
+ }
+ if (isSame == Standard_True)
+ return Standard_True;
+ }
+ }
+ else
{
- isSame = Standard_True;
- for(j = 1; j <= myN; j++)
+ NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion());
+ if (isFirstCellFilterInvoke)
{
- if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) - myA(j)) * mySameTol)
+ myFilter.Reset(myCellSize);
+
+ // Copy initial data into cell filter.
+ for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
{
- isSame = Standard_False;
- break;
+ math_Vector aVec(1, myN);
+ for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
+ aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
+
+ myFilter.Add(aVec, aVec);
}
}
- if (isSame == Standard_True)
- return Standard_True;
+ isFirstCellFilterInvoke = Standard_False;
+
+ math_Vector aLow(1, myN), anUp(1, myN);
+ anInspector.Shift(thePnt, myCellSize, aLow, anUp);
+
+ anInspector.ClearFind();
+ anInspector.SetCurrent(thePnt);
+ myFilter.Inspect(aLow, anUp, anInspector);
+ if (!anInspector.isFind())
+ {
+ // Point is out of close cells, add new one.
+ myFilter.Add(thePnt, thePnt);
+ }
}
return Standard_False;
}
//=======================================================================
-//function : NbExtrema
+//function : Points
+//purpose :
+//=======================================================================
+void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
+{
+ Standard_Integer j;
+
+ for(j = 1; j <= myN; j++)
+ theSol(j) = myY((theIndex - 1) * myN + j);
+}
+
+//=======================================================================
+//function : initCellSize
//purpose :
//=======================================================================
-Standard_Integer math_GlobOptMin::NbExtrema()
+void math_GlobOptMin::initCellSize()
{
- return mySolCount;
+ for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
+ {
+ myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
+ * Precision::PConfusion() / (2.0 * Sqrt(2.0));
+ }
}
//=======================================================================
-//function : GetF
+//function : CheckFunctionalStopCriteria
//purpose :
//=======================================================================
-Standard_Real math_GlobOptMin::GetF()
+Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
{
- return myF;
+ // Search single solution and current solution in its neighborhood.
+ if (myIsFindSingleSolution &&
+ Abs (myF - myFunctionalMinimalValue) < mySameTol * 0.01)
+ return Standard_True;
+
+ return Standard_False;
}
//=======================================================================
-//function : IsDone
+//function : ComputeInitSol
//purpose :
//=======================================================================
-Standard_Boolean math_GlobOptMin::isDone()
+void math_GlobOptMin::ComputeInitSol()
{
- return myDone;
+ Standard_Real aVal;
+ math_Vector aPnt(1, myN);
+
+ // Check functional value in midpoint. It is necessary since local optimization
+ // algorithm may fail and return nothing. This is a protection from uninitialized
+ // variables.
+ aPnt = (myGlobA + myGlobB) * 0.5;
+ myFunc->Value(aPnt, aVal);
+ checkAddCandidate(aPnt, aVal);
+
+ // Run local optimization from lower corner, midpoint, and upper corner.
+ for(Standard_Integer i = 1; i <= 3; i++)
+ {
+ aPnt = myA + (myB - myA) * (i - 1) / 2.0;
+
+ if(computeLocalExtremum(aPnt, aVal, aPnt))
+ checkAddCandidate(aPnt, aVal);
+ }
}
//=======================================================================
-//function : Points
+//function : checkAddCandidate
//purpose :
//=======================================================================
-void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
+void math_GlobOptMin::checkAddCandidate(const math_Vector& thePnt,
+ const Standard_Real theValue)
{
- Standard_Integer j;
+ if (Abs(theValue - myF) < mySameTol * 0.01 && // Value in point is close to optimal value.
+ !myIsFindSingleSolution) // Several optimal solutions are allowed.
+ {
+ if (!isStored(thePnt))
+ {
+ if ((theValue - myF) * myZ > 0.0)
+ myF = theValue;
+ for (Standard_Integer j = 1; j <= myN; j++)
+ myY.Append(thePnt(j));
+ mySolCount++;
+ }
+ }
- for(j = 1; j <= myN; j++)
- theSol(j) = myY((theIndex - 1) * myN + j);
+ // New best solution:
+ // new point is out of (mySameTol * 0.01) surrounding or
+ // new point is better than old and single point search.
+ Standard_Real aDelta = (theValue - myF) * myZ;
+ if (aDelta > mySameTol * 0.01 ||
+ (aDelta > 0.0 && myIsFindSingleSolution))
+ {
+ myF = theValue;
+ myY.Clear();
+ for (Standard_Integer j = 1; j <= myN; j++)
+ myY.Append(thePnt(j));
+ mySolCount = 1;
+
+ isFirstCellFilterInvoke = Standard_True;
+ }
}