// 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 <Standard_Real.hxx>
#include <Precision.hxx>
+//=======================================================================
+//function : DistanceToBorder
+//purpose :
+//=======================================================================
+static Standard_Real DistanceToBorder(const math_Vector & theX,
+ const math_Vector & theMin,
+ const math_Vector & theMax)
+{
+ 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
myB(1, myN),
myGlobA(1, myN),
myGlobB(1, myN),
+ myIsConstLocked(Standard_False),
myX(1, myN),
myTmp(1, myN),
myV(1, myN),
myMaxV(1, myN),
- myExpandCoeff(1, myN)
+ myCellSize(0, myN - 1),
+ myFilter(theFunc->NbVariables()),
+ myCont(2)
{
Standard_Integer i;
myFunc = theFunc;
myC = theC;
+ myInitC = theC;
+ myIsFindSingleSolution = Standard_False;
+ myFunctionalMinimalValue = -Precision::Infinite();
myZ = -1;
mySolCount = 0;
myMaxV(i) = (myB(i) - myA(i)) / 3.0;
}
- myExpandCoeff(1) = 1.0;
- for(i = 2; i <= myN; i++)
- {
- myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
- }
-
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;
myMaxV(i) = (myB(i) - myA(i)) / 3.0;
}
- myExpandCoeff(1) = 1.0;
- for(i = 2; i <= myN; i++)
- {
- myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
- }
-
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);
myMaxV(i) = (myB(i) - myA(i)) / 3.0;
}
- myExpandCoeff(1) = 1.0;
- for(i = 2; i <= myN; i++)
- {
- myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
- }
-
myDone = Standard_False;
}
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;
}
if (minLength < Precision::PConfusion())
return;
}
- // Compute initial values for myF, myY, myC.
- computeInitialValues();
+ if (!myIsConstLocked)
+ {
+ // Compute initial value for myC.
+ computeInitialValues();
+ }
myE1 = minLength * myTol;
myE2 = maxLength * myTol;
- if (myC > 1.0)
- myE3 = - maxLength * myTol / 4.0;
+
+ myIsFindSingleSolution = isFindSingleSolution;
+ if (isFindSingleSolution)
+ {
+ // Run local optimization if current value better than optimal.
+ myE3 = 0.0;
+ }
else
- myE3 = - maxLength * myTol * myC / 4.0;
+ {
+ if (myC > 1.0)
+ myE3 = - maxLength * myTol / 4.0;
+ else
+ myE3 = - maxLength * myTol * myC / 4.0;
+ }
+
+ // 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);
- newtonMinimum.Perform(*myTmp, thePnt);
+ math_NewtonMinimum newtonMinimum(*aTmp);
+ newtonMinimum.SetBoundary(myGlobA, myGlobB);
+ newtonMinimum.Perform(*aTmp, thePnt);
if (newtonMinimum.IsDone())
{
} 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->NbVariables());
- bfgs.Perform(*myTmp, thePnt);
+ math_BFGS bfgs(aTmp->NbVariables());
+ bfgs.Perform(*aTmp, thePnt);
if (bfgs.IsDone())
{
bfgs.Location(theOutPnt);
math_Vector aBestPnt(1, myN);
math_Vector aParamStep(1, myN);
Standard_Real aCurrVal = RealLast();
- Standard_Real aBestVal = RealLast();
-
- // Check functional value in midpoint, low and upp point border and
- // in each point try to perform local optimization.
- aBestPnt = (myA + myB) * 0.5;
- myFunc->Value(aBestPnt, aBestVal);
-
- for(i = 1; i <= 3; i++)
- {
- aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
- if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
- {
- // Local Extremum finds better solution than current point.
- if (aCurrVal < aBestVal)
- {
- aBestVal = aCurrVal;
- aBestPnt = aCurrPnt;
- }
- }
- }
-
- myF = aBestVal;
- myY.Clear();
- for(i = 1; i <= myN; i++)
- myY.Append(aBestPnt(i));
- mySolCount++;
-
- // Lipschitz const approximation
+ // Lipchitz const approximation.
Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
Standard_Integer aPntNb = 13;
myFunc->Value(myA, aPrevValDiag);
aPrevValProj = aCurrVal;
}
+ myC = myInitC;
aLipConst *= Sqrt(myN) / aStep;
-
if (aLipConst < myC * 0.1)
- {
myC = Max(aLipConst * 0.1, 0.01);
- }
- else if (aLipConst > myC * 10)
+ else if (aLipConst > myC * 5)
+ myC = Min(myC * 5, 50.0);
+
+ // Clear all solutions except one.
+ if (myY.Size() != myN)
{
- myC = Min(myC * 2, 30.0);
+ for(i = 1; i <= myN; i++)
+ aBestPnt(i) = myY(i);
+ myY.Clear();
+ for(i = 1; i <= myN; i++)
+ myY.Append(aBestPnt(i));
}
+ mySolCount = 1;
}
//=======================================================================
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 aStepBestValue = RealLast();
- Standard_Real aRealStep = 0.0;
math_Vector aStepBestPoint(1, myN);
- Standard_Boolean isInside = Standard_False;
- Standard_Real r;
+ 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 + myZ * myC * aRealStep - 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);
+ Standard_Real aSaveParam = myX(1);
+
+ // 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.
+ {
+ isInside = computeLocalExtremum(myX, val, myTmp);
+ }
}
aStepBestValue = (isInside && (val < d))? val : d;
aStepBestPoint = (isInside && (val < d))? myTmp : myX;
- // Solutions are close to each other.
- if (Abs(aStepBestValue - myF) < mySameTol * 0.01)
+ // Solutions are close to each other
+ // and it is allowed to have more than one solution.
+ if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
+ !myIsFindSingleSolution)
{
if (!isStored(aStepBestPoint))
{
}
}
- // New best solution.
- if ((aStepBestValue - myF) * myZ > mySameTol * 0.01)
+ // New best solution:
+ // new point is out of (mySameTol * 0.01) surrounding or
+ // new point is better than old + single point search.
+ Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
+ if (aFunctionalDelta > mySameTol * 0.01 ||
+ (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
{
mySolCount = 0;
myF = aStepBestValue;
for(i = 1; i <= myN; i++)
myY.Append(aStepBestPoint(i));
mySolCount++;
+
+ isFirstCellFilterInvoke = Standard_True;
}
- aRealStep = myE2 + Abs(myF - d) / myC;
- myV(1) = Min(aRealStep, myMaxV(1));
+ if (CheckFunctionalStopCriteria())
+ return; // Best possible value is obtained.
+
+ myV(1) = Min(myE2 + Abs(myF - d) / myC, myMaxV(1));
+ myLastStep = myV(1);
}
else
{
for(i = 1; i < j; i++)
myV(i) = 0.0;
}
- // Compute step in (j + 1) dimension according to scale.
if (j < myN)
{
- Standard_Real aUpperDimStep = myV(j) * myExpandCoeff(j + 1);
+ Standard_Real aUpperDimStep = Max(myV(j), myE2);
if (myV(j + 1) > aUpperDimStep)
{
if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
{
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)
{
- isSame = Standard_True;
- for(j = 1; j <= myN; j++)
+ for(i = 0; i < mySolCount; i++)
{
- if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) - myA(j)) * mySameTol)
+ isSame = Standard_True;
+ for(j = 1; j <= myN; j++)
{
- isSame = Standard_False;
- break;
+ if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
+ {
+ isSame = Standard_False;
+ break;
+ }
}
+ if (isSame == Standard_True)
+ return Standard_True;
}
- if (isSame == Standard_True)
- return Standard_True;
+ }
+ else
+ {
+ NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion());
+ if (isFirstCellFilterInvoke)
+ {
+ myFilter.Reset(myCellSize);
+
+ // Copy initial data into cell filter.
+ for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
+ {
+ math_Vector aVec(1, myN);
+ for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
+ aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
+
+ myFilter.Add(aVec, aVec);
+ }
+ }
+
+ 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 :
//=======================================================================
-Standard_Integer math_GlobOptMin::NbExtrema()
+void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
{
- return mySolCount;
+ Standard_Integer j;
+
+ for(j = 1; j <= myN; j++)
+ theSol(j) = myY((theIndex - 1) * myN + j);
}
//=======================================================================
-//function : GetF
+//function : initCellSize
//purpose :
//=======================================================================
-Standard_Real math_GlobOptMin::GetF()
+void math_GlobOptMin::initCellSize()
{
- return myF;
+ for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
+ {
+ myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
+ * Precision::PConfusion() / (2.0 * Sqrt(2.0));
+ }
}
//=======================================================================
-//function : IsDone
+//function : CheckFunctionalStopCriteria
//purpose :
//=======================================================================
-Standard_Boolean math_GlobOptMin::isDone()
+Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
{
- return myDone;
+ // Search single solution and current solution in its neighborhood.
+ if (myIsFindSingleSolution &&
+ Abs (myF - myFunctionalMinimalValue) < mySameTol * 0.01)
+ return Standard_True;
+
+ return Standard_False;
}
//=======================================================================
-//function : Points
+//function : ComputeInitSol
//purpose :
//=======================================================================
-void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
+void math_GlobOptMin::ComputeInitSol()
{
- Standard_Integer j;
+ Standard_Real aCurrVal, aBestVal;
+ math_Vector aCurrPnt(1, myN);
+ math_Vector aBestPnt(1, myN);
+ math_Vector aParamStep(1, myN);
+ // Check functional value in midpoint, lower and upper border points and
+ // in each point try to perform local optimization.
+ aBestPnt = (myGlobA + myGlobB) * 0.5;
+ myFunc->Value(aBestPnt, aBestVal);
- for(j = 1; j <= myN; j++)
- theSol(j) = myY((theIndex - 1) * myN + j);
+ Standard_Integer i;
+ for(i = 1; i <= 3; i++)
+ {
+ aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
+
+ if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
+ {
+ // Local search tries to find better solution than current point.
+ if (aCurrVal < aBestVal)
+ {
+ aBestVal = aCurrVal;
+ aBestPnt = aCurrPnt;
+ }
+ }
+ }
+
+ myF = aBestVal;
+ myY.Clear();
+ for(i = 1; i <= myN; i++)
+ myY.Append(aBestPnt(i));
+ mySolCount = 1;
}
projects / occt.git / blobdiff