Extrema between curves has been made producing correct result for the cases of solution located near bounds.
- Class math_GlobOptMin has been improved to use lower order methods of local optimization when high-order methods are failed.
- Add support of conditional optimization (in bounds) in the classes math_BFGS and math_BracketMinimum.
- Turn on conditional optimization in the case of usage of math_BFGS in the class math_GlobOptMin.
- Correct mistake in distmini command, which caused incorrect reading of deflection parameter.
- To avoid possible FPE signals, ensure initialization of fields in the class math/math_BracketMinimum.
- In the algorithms math_BFGS, math_Powell and math_FRPR, take into account that the function math_MultipleVarFunction can return failure status (e.g. when computing D0 out of bounds).
New test cases have been added.
Tests cases are updated.
// correct test case
Standard_Real aDeflection = Precision::Confusion();
if (n == 5)
- aDeflection = Draw::Atoi(a[4]);
+ aDeflection = Draw::Atof(a[4]);
BRepExtrema_DistShapeShape dst(S1 ,S2, aDeflection);
#include <math_MultipleVarFunctionWithGradient.hxx>
#include <Standard_DimensionError.hxx>
#include <StdFail_NotDone.hxx>
+#include <Precision.hxx>
#define R 0.61803399
#define C (1.0-R)
{
*P = *Dir;
P->Multiply(x);
- P->Add(*P0);
- F->Value(*P, fval);
- return Standard_True;
+ P->Add(*P0);
+ fval = 0.;
+ return F->Value(*P, fval);
}
Standard_Boolean DirFunction::Values(const Standard_Real x,
{
*P = *Dir;
P->Multiply(x);
- P->Add(*P0);
- F->Values(*P, fval, *G);
- D = (*G).Multiplied(*Dir);
- return Standard_True;
+ P->Add(*P0);
+ fval = D = 0.;
+ if (F->Values(*P, fval, *G))
+ {
+ D = (*G).Multiplied(*Dir);
+ return Standard_True;
+ }
+ return Standard_False;
}
Standard_Boolean DirFunction::Derivative(const Standard_Real x,
Standard_Real& D)
P->Multiply(x);
P->Add(*P0);
Standard_Real fval;
- F->Values(*P, fval, *G);
- D = (*G).Multiplied(*Dir);
- return Standard_True;
+ D = 0.;
+ if (F->Values(*P, fval, *G))
+ {
+ D = (*G).Multiplied(*Dir);
+ return Standard_True;
+ }
+ return Standard_False;
}
+//=======================================================================
+//function : ComputeInitScale
+//purpose : Compute the appropriate initial value of scale factor to apply
+// to the direction to approach to the minimum of the function
+//=======================================================================
+static Standard_Boolean ComputeInitScale(const Standard_Real theF0,
+ const math_Vector& theDir,
+ const math_Vector& theGr,
+ Standard_Real& theScale)
+{
+ Standard_Real dy1 = theGr * theDir;
+ if (Abs(dy1) < RealSmall())
+ return Standard_False;
+ Standard_Real aHnr1 = theDir.Norm2();
+ Standard_Real alfa = 0.7*(-theF0) / dy1;
+ theScale = 0.015 / Sqrt(aHnr1);
+ if (theScale > alfa)
+ theScale = alfa;
+ return Standard_True;
+}
+
+//=======================================================================
+//function : ComputeMinMaxScale
+//purpose : For a given point and direction, and bounding box,
+// find min and max scale factors with which the point reaches borders
+// if we apply translation Point+Dir*Scale.
+//return : True if found, False if point is out of bounds.
+//=======================================================================
+static Standard_Boolean ComputeMinMaxScale(const math_Vector& thePoint,
+ const math_Vector& theDir,
+ const math_Vector& theLeft,
+ const math_Vector& theRight,
+ Standard_Real& theMinScale,
+ Standard_Real& theMaxScale)
+{
+ Standard_Integer anIdx;
+ for (anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
+ {
+ Standard_Real aLeft = theLeft(anIdx) - thePoint(anIdx);
+ Standard_Real aRight = theRight(anIdx) - thePoint(anIdx);
+ if (Abs(theDir(anIdx)) > RealSmall())
+ {
+ // use PConfusion to get off a little from the bounds to prevent
+ // possible refuse in Value function.
+ Standard_Real aLScale = (aLeft + Precision::PConfusion()) / theDir(anIdx);
+ Standard_Real aRScale = (aRight - Precision::PConfusion()) / theDir(anIdx);
+ if (Abs(aLeft) < Precision::PConfusion())
+ {
+ // point is on the left border
+ theMaxScale = Min(theMaxScale, Max(0., aRScale));
+ theMinScale = Max(theMinScale, Min(0., aRScale));
+ }
+ else if (Abs(aRight) < Precision::PConfusion())
+ {
+ // point is on the right border
+ theMaxScale = Min(theMaxScale, Max(0., aLScale));
+ theMinScale = Max(theMinScale, Min(0., aLScale));
+ }
+ else if (aLeft * aRight < 0)
+ {
+ // point is inside allowed range
+ theMaxScale = Min(theMaxScale, Max(aLScale, aRScale));
+ theMinScale = Max(theMinScale, Min(aLScale, aRScale));
+ }
+ else
+ // point is out of bounds
+ return Standard_False;
+ }
+ else
+ {
+ // Direction is parallel to the border.
+ // Check that the point is not out of bounds
+ if (aLeft > Precision::PConfusion() ||
+ aRight < -Precision::PConfusion())
+ {
+ return Standard_False;
+ }
+ }
+ }
+ return Standard_True;
+}
static Standard_Boolean MinimizeDirection(math_Vector& P,
- Standard_Real F0,
- math_Vector& Gr,
- math_Vector& Dir,
- Standard_Real& Result,
- DirFunction& F) {
-
-
-
- Standard_Real ax, xx, bx, Fax, Fxx, Fbx, F1;
- F.Initialize(P, Dir);
-
- Standard_Real dy1, Hnr1, lambda, alfa=0;
- dy1 = Gr*Dir;
- if (dy1 != 0) {
- Hnr1 = Dir.Norm2();
- alfa = 0.7*(-F0)/dy1;
- lambda = 0.015/Sqrt(Hnr1);
- }
- else lambda = 1.0;
- if (lambda > alfa) lambda = alfa;
- F.Value(lambda, F1);
- math_BracketMinimum Bracket(F, 0.0, lambda, F0, F1);
- if(Bracket.IsDone()) {
- Bracket.Values(ax, xx, bx);
- Bracket.FunctionValues(Fax, Fxx, Fbx);
-
- Standard_Integer niter = 100;
- Standard_Real tol = 1.e-03;
- math_BrentMinimum Sol(tol, Fxx, niter, 1.e-08);
- Sol.Perform(F, ax, xx, bx);
- if(Sol.IsDone()) {
- Standard_Real Scale = Sol.Location();
- Result = Sol.Minimum();
- Dir.Multiply(Scale);
- P.Add(Dir);
- return Standard_True;
- }
- }
- return Standard_False;
- }
+ Standard_Real F0,
+ math_Vector& Gr,
+ math_Vector& Dir,
+ Standard_Real& Result,
+ DirFunction& F,
+ Standard_Boolean isBounds,
+ const math_Vector& theLeft,
+ const math_Vector& theRight)
+{
+ Standard_Real lambda;
+ if (!ComputeInitScale(F0, Dir, Gr, lambda))
+ return Standard_False;
+
+ // by default the scaling range is unlimited
+ Standard_Real aMinLambda = -Precision::Infinite();
+ Standard_Real aMaxLambda = Precision::Infinite();
+ if (isBounds)
+ {
+ // limit the scaling range taking into account the bounds
+ if (!ComputeMinMaxScale(P, Dir, theLeft, theRight, aMinLambda, aMaxLambda))
+ return Standard_False;
+
+ if (aMinLambda > -Precision::PConfusion() && aMaxLambda < Precision::PConfusion())
+ {
+ // Point is on the border and the direction shows outside.
+ // Make direction to go along the border
+ for (Standard_Integer anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
+ {
+ if (Abs(P(anIdx) - theRight(anIdx)) < Precision::PConfusion() ||
+ Abs(P(anIdx) - theLeft(anIdx)) < Precision::PConfusion())
+ {
+ Dir(anIdx) = 0.0;
+ }
+ }
+
+ // re-compute scale values with new direction
+ if (!ComputeInitScale(F0, Dir, Gr, lambda))
+ return Standard_False;
+ if (!ComputeMinMaxScale(P, Dir, theLeft, theRight, aMinLambda, aMaxLambda))
+ return Standard_False;
+ }
+ lambda = Min(lambda, aMaxLambda);
+ lambda = Max(lambda, aMinLambda);
+ }
+
+ F.Initialize(P, Dir);
+ Standard_Real F1;
+ if (!F.Value(lambda, F1))
+ return Standard_False;
+ math_BracketMinimum Bracket(0.0, lambda);
+ if (isBounds)
+ Bracket.SetLimits(aMinLambda, aMaxLambda);
+ Bracket.SetFA(F0);
+ Bracket.SetFB(F1);
+ Bracket.Perform(F);
+ if (Bracket.IsDone()) {
+ // find minimum inside the bracket
+ Standard_Real ax, xx, bx, Fax, Fxx, Fbx;
+ Bracket.Values(ax, xx, bx);
+ Bracket.FunctionValues(Fax, Fxx, Fbx);
+
+ Standard_Integer niter = 100;
+ Standard_Real tol = 1.e-03;
+ math_BrentMinimum Sol(tol, Fxx, niter, 1.e-08);
+ Sol.Perform(F, ax, xx, bx);
+ if (Sol.IsDone()) {
+ Standard_Real Scale = Sol.Location();
+ Result = Sol.Minimum();
+ Dir.Multiply(Scale);
+ P.Add(Dir);
+ return Standard_True;
+ }
+ }
+ else if (isBounds)
+ {
+ // Bracket definition is failure. If the bounds are defined then
+ // set current point to intersection with bounds
+ Standard_Real aFMin, aFMax;
+ if (!F.Value(aMinLambda, aFMin))
+ return Standard_False;
+ if (!F.Value(aMaxLambda, aFMax))
+ return Standard_False;
+ Standard_Real aBestLambda;
+ if (aFMin < aFMax)
+ {
+ aBestLambda = aMinLambda;
+ Result = aFMin;
+ }
+ else
+ {
+ aBestLambda = aMaxLambda;
+ Result = aFMax;
+ }
+ Dir.Multiply(aBestLambda);
+ P.Add(Dir);
+ return Standard_True;
+ }
+ return Standard_False;
+}
void math_BFGS::Perform(math_MultipleVarFunctionWithGradient& F,
for(nbiter = 1; nbiter <= Itermax; nbiter++) {
TheMinimum = PreviousMinimum;
Standard_Boolean IsGood = MinimizeDirection(TheLocation, TheMinimum,
- TheGradient,
- xi, TheMinimum, F_Dir);
+ TheGradient,
+ xi, TheMinimum, F_Dir,
+ myIsBoundsDefined, myLeft, myRight);
if(!IsGood) {
Done = Standard_False;
TheStatus = math_DirectionSearchError;
const Standard_Real Tolerance,
const Standard_Integer NbIterations,
const Standard_Real ZEPS)
- : TheLocation(1, NbVariables),
- TheGradient(1, NbVariables) {
-
- XTol = Tolerance;
- EPSZ = ZEPS;
- Itermax = NbIterations;
+ : TheStatus(math_OK),
+ TheLocation(1, NbVariables),
+ TheGradient(1, NbVariables),
+ PreviousMinimum(0.),
+ TheMinimum(0.),
+ XTol(Tolerance),
+ EPSZ(ZEPS),
+ nbiter(0),
+ myIsBoundsDefined(Standard_False),
+ myLeft(1, NbVariables, 0.0),
+ myRight(1, NbVariables, 0.0),
+ Done(Standard_False),
+ Itermax(NbIterations)
+ {
}
}
}
-
-
+//=======================================================================
+//function : SetBoundary
+//purpose : Set boundaries for conditional optimization
+//=======================================================================
+void math_BFGS::SetBoundary(const math_Vector& theLeftBorder,
+ const math_Vector& theRightBorder)
+{
+ myLeft = theLeftBorder;
+ myRight = theRightBorder;
+ myIsBoundsDefined = Standard_True;
+}
//! This class implements the Broyden-Fletcher-Goldfarb-Shanno variant of
//! Davidson-Fletcher-Powell minimization algorithm of a function of
//! multiple variables.Knowledge of the function's gradient is required.
+//!
+//! It is possible to solve conditional optimization problem on hyperparallelepiped.
+//! Method SetBoundary is used to define hyperparallelepiped borders. With boundaries
+//! defined, the algorithm will not make evaluations of the function outside of the
+//! borders.
class math_BFGS
{
public:
//! initialization to effectively compute the minimum of the
//! function F.
Standard_EXPORT math_BFGS(const Standard_Integer NbVariables, const Standard_Real Tolerance = 1.0e-8, const Standard_Integer NbIterations = 200, const Standard_Real ZEPS = 1.0e-12);
-Standard_EXPORT virtual ~math_BFGS();
-
+
+ Standard_EXPORT virtual ~math_BFGS();
+
+ //! Set boundaries for conditional optimization.
+ //! The expected indices range of vectors is [1, NbVariables].
+ Standard_EXPORT void SetBoundary(const math_Vector& theLeftBorder, const math_Vector& theRightBorder);
+
//! Given the starting point StartingPoint,
//! minimization is done on the function F.
//! The solution F = Fi is found when :
Standard_Real XTol;
Standard_Real EPSZ;
Standard_Integer nbiter;
+ Standard_Boolean myIsBoundsDefined;
+ math_Vector myLeft;
+ math_Vector myRight;
private:
#define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d)
+Standard_Boolean math_BracketMinimum::LimitAndMayBeSwap
+ (math_Function& F,
+ const Standard_Real theA,
+ Standard_Real& theB,
+ Standard_Real& theFB,
+ Standard_Real& theC,
+ Standard_Real& theFC) const
+{
+ theC = Limited(theC);
+ if (Abs(theB - theC) < Precision::PConfusion())
+ return Standard_False;
+ Standard_Boolean OK = F.Value(theC, theFC);
+ if (!OK)
+ return Standard_False;
+ // check that B is between A and C
+ if ((theA - theB) * (theB - theC) < 0)
+ {
+ // swap B and C
+ Standard_Real dum;
+ SHFT(dum, theB, theC, dum);
+ SHFT(dum, theFB, theFC, dum);
+ }
+ return Standard_True;
+}
- void math_BracketMinimum::Perform(math_Function& F,
- const Standard_Real A,
- const Standard_Real B) {
+ void math_BracketMinimum::Perform(math_Function& F)
+ {
Standard_Boolean OK;
- Standard_Real ulim, u, r, q, f, fu, dum;
+ Standard_Real ulim, u, r, q, fu, dum;
Done = Standard_False;
- Ax = A;
- Bx = B;
Standard_Real Lambda = GOLD;
if (!myFA) {
OK = F.Value(Ax, FAx);
SHFT(dum, Ax, Bx, dum);
SHFT(dum, FBx, FAx, dum);
}
+
+ // get next prob after (A, B)
Cx = Bx + Lambda * (Bx - Ax);
- OK = F.Value(Cx, FCx);
- if(!OK) return;
+ if (myIsLimited)
+ {
+ OK = LimitAndMayBeSwap(F, Ax, Bx, FBx, Cx, FCx);
+ if (!OK)
+ return;
+ }
+ else
+ {
+ OK = F.Value(Cx, FCx);
+ if (!OK)
+ return;
+ }
+
while(FBx > FCx) {
r = (Bx - Ax) * (FBx -FCx);
q = (Bx - Cx) * (FBx -FAx);
u = Bx - ((Bx - Cx) * q - (Bx - Ax) * r) /
(2.0 * SIGN(MAX(fabs(q - r), TINY), q - r));
ulim = Bx + GLIMIT * (Cx - Bx);
- if((Bx - u) * (u - Cx) > 0.0) {
+ if (myIsLimited)
+ ulim = Limited(ulim);
+ if ((Bx - u) * (u - Cx) > 0.0) {
+ // u is between B and C
OK = F.Value(u, fu);
if(!OK) return;
if(fu < FCx) {
+ // solution is found (B, u, c)
Ax = Bx;
Bx = u;
FAx = FBx;
return;
}
else if(fu > FBx) {
+ // solution is found (A, B, u)
Cx = u;
FCx = fu;
Done = Standard_True;
return;
}
+ // get next prob after (B, C)
u = Cx + Lambda * (Cx - Bx);
- OK = F.Value(u, fu);
- if(!OK) return;
+ if (myIsLimited)
+ {
+ OK = LimitAndMayBeSwap(F, Bx, Cx, FCx, u, fu);
+ if (!OK)
+ return;
+ }
+ else
+ {
+ OK = F.Value(u, fu);
+ if (!OK)
+ return;
+ }
}
else if((Cx - u) * (u - ulim) > 0.0) {
+ // u is beyond C but between C and limit
OK = F.Value(u, fu);
if(!OK) return;
- if(fu < FCx) {
- SHFT(Bx, Cx, u, Cx + GOLD * (Cx - Bx));
- OK = F.Value(u, f);
- if(!OK) return;
- SHFT(FBx, FCx, fu, f);
- }
}
else if ((u - ulim) * (ulim - Cx) >= 0.0) {
+ // u is beyond limit
u = ulim;
OK = F.Value(u, fu);
if(!OK) return;
}
else {
+ // u tends to approach to the side of A,
+ // so reset it to the next prob after (B, C)
u = Cx + GOLD * (Cx - Bx);
- OK = F.Value(u, fu);
- if(!OK) return;
+ if (myIsLimited)
+ {
+ OK = LimitAndMayBeSwap(F, Bx, Cx, FCx, u, fu);
+ if (!OK)
+ return;
+ }
+ else
+ {
+ OK = F.Value(u, fu);
+ if (!OK)
+ return;
+ }
}
SHFT(Ax, Bx, Cx, u);
SHFT(FAx, FBx, FCx, fu);
math_BracketMinimum::math_BracketMinimum(math_Function& F,
const Standard_Real A,
- const Standard_Real B) {
-
- myFA = Standard_False;
- myFB = Standard_False;
- Perform(F, A, B);
+ const Standard_Real B)
+ : Done(Standard_False),
+ Ax(A), Bx(B), Cx(0.),
+ FAx(0.), FBx(0.), FCx(0.),
+ myLeft(-Precision::Infinite()),
+ myRight(Precision::Infinite()),
+ myIsLimited(Standard_False),
+ myFA(Standard_False),
+ myFB (Standard_False)
+ {
+ Perform(F);
}
math_BracketMinimum::math_BracketMinimum(math_Function& F,
const Standard_Real A,
const Standard_Real B,
- const Standard_Real FA) {
- FAx = FA;
- myFA = Standard_True;
- myFB = Standard_False;
- Perform(F, A, B);
+ const Standard_Real FA)
+ : Done(Standard_False),
+ Ax(A), Bx(B), Cx(0.),
+ FAx(FA), FBx(0.), FCx(0.),
+ myLeft(-Precision::Infinite()),
+ myRight(Precision::Infinite()),
+ myIsLimited(Standard_False),
+ myFA(Standard_True),
+ myFB (Standard_False)
+ {
+ Perform(F);
}
math_BracketMinimum::math_BracketMinimum(math_Function& F,
const Standard_Real A,
const Standard_Real B,
const Standard_Real FA,
- const Standard_Real FB) {
- FAx = FA;
- FBx = FB;
- myFA = Standard_True;
- myFB = Standard_True;
- Perform(F, A, B);
+ const Standard_Real FB)
+ : Done(Standard_False),
+ Ax(A), Bx(B), Cx(0.),
+ FAx(FA), FBx(FB), FCx(0.),
+ myLeft(-Precision::Infinite()),
+ myRight(Precision::Infinite()),
+ myIsLimited(Standard_False),
+ myFA(Standard_True),
+ myFB(Standard_True)
+ {
+ Perform(F);
}
//! Given two distinct initial points, BracketMinimum
//! implements the computation of three points (a, b, c) which
//! bracket the minimum of the function and verify A less than
-//! B, B less than C and F(A) less than F(B), F(B) less than (C).
-class math_BracketMinimum
+//! B, B less than C and F(B) less than F(A), F(B) less than F(C).
+//!
+//! The algorithm supports conditional optimization. By default no limits are
+//! applied to the parameter change. The method SetLimits defines the allowed range.
+//! If no minimum is found in limits then IsDone() will return false. The user
+//! is in charge of providing A and B to be in limits.
+class math_BracketMinimum
{
public:
DEFINE_STANDARD_ALLOC
-
+ //! Constructor preparing A and B parameters only. It does not perform the job.
+ Standard_EXPORT math_BracketMinimum(const Standard_Real A, const Standard_Real B);
//! Given two initial values this class computes a
//! bracketing triplet of abscissae Ax, Bx, Cx
//! on the function F.
//! This constructor has to be used if F(A) and F(B) are known.
Standard_EXPORT math_BracketMinimum(math_Function& F, const Standard_Real A, const Standard_Real B, const Standard_Real FA, const Standard_Real FB);
-
+
+ //! Set limits of the parameter. By default no limits are applied to the parameter change.
+ //! If no minimum is found in limits then IsDone() will return false. The user
+ //! is in charge of providing A and B to be in limits.
+ void SetLimits(const Standard_Real theLeft, const Standard_Real theRight);
+
+ //! Set function value at A
+ void SetFA(const Standard_Real theValue);
+
+ //! Set function value at B
+ void SetFB(const Standard_Real theValue);
+
+ //! The method performing the job. It is called automatically by constructors with the function.
+ Standard_EXPORT void Perform(math_Function& F);
+
//! Returns true if the computations are successful, otherwise returns false.
- Standard_Boolean IsDone() const;
+ Standard_Boolean IsDone() const;
//! Returns the bracketed triplet of abscissae.
//! Exceptions
//! Is used to redefine the operator <<.
Standard_EXPORT void Dump (Standard_OStream& o) const;
+private:
+ //! Limit the given value to become within the range [myLeft, myRight].
+ Standard_Real Limited(const Standard_Real theValue) const;
-
-protected:
-
-
- //! Is used internally by the constructors.
- Standard_EXPORT void Perform (math_Function& F, const Standard_Real A, const Standard_Real B);
-
-
-
+ //! Limit the value of C (see Limited) and compute the function in it.
+ //! If C occurs to be between A and B then swap parameters and function
+ //! values of B and C.
+ //! Return false in the case of C becomes equal to B or function calculation
+ //! failure.
+ Standard_Boolean LimitAndMayBeSwap(math_Function& F, const Standard_Real theA,
+ Standard_Real& theB, Standard_Real& theFB,
+ Standard_Real& theC, Standard_Real& theFC) const;
private:
-
-
Standard_Boolean Done;
Standard_Real Ax;
Standard_Real Bx;
Standard_Real FAx;
Standard_Real FBx;
Standard_Real FCx;
+ Standard_Real myLeft;
+ Standard_Real myRight;
+ Standard_Boolean myIsLimited;
Standard_Boolean myFA;
Standard_Boolean myFB;
// math_BracketMinimum.lxx
+#include <Precision.hxx>
+
+inline math_BracketMinimum::math_BracketMinimum(const Standard_Real A,
+ const Standard_Real B)
+: Done(Standard_False),
+ Ax(A), Bx(B), Cx(0.),
+ FAx(0.), FBx(0.), FCx(0.),
+ myLeft(-Precision::Infinite()),
+ myRight(Precision::Infinite()),
+ myIsLimited(Standard_False),
+ myFA(Standard_False),
+ myFB(Standard_False)
+{
+}
+
inline Standard_Boolean math_BracketMinimum::IsDone() const { return Done; }
inline Standard_OStream& operator<<(Standard_OStream& o,
Br.Dump(o);
return o;
}
+
+inline void math_BracketMinimum::SetLimits(const Standard_Real theLeft,
+ const Standard_Real theRight)
+{
+ myLeft = theLeft;
+ myRight = theRight;
+ myIsLimited = Standard_True;
+}
+
+inline void math_BracketMinimum::SetFA(const Standard_Real theValue)
+{
+ FAx = theValue;
+ myFA = Standard_True;
+}
+
+inline void math_BracketMinimum::SetFB(const Standard_Real theValue)
+{
+ FBx = theValue;
+ myFB = Standard_True;
+}
+
+inline Standard_Real math_BracketMinimum::Limited(const Standard_Real theValue) const
+{
+ return theValue < myLeft ? myLeft :
+ theValue > myRight ? myRight : theValue;
+}
*P = *Dir;
P->Multiply(x);
- P->Add(*P0);
- F->Value(*P, fval);
- return Standard_True;
+ P->Add(*P0);
+ fval = 0.;
+ return F->Value(*P, fval);
}
static Standard_Boolean MinimizeDirection(math_Vector& P,
{
newtonMinimum.Location(theOutPnt);
theVal = newtonMinimum.Minimum();
+
+ if (isInside(theOutPnt))
+ return Standard_True;
}
- else return Standard_False;
- } else
+ }
// BFGS method used.
if (myCont >= 1 &&
math_MultipleVarFunctionWithGradient* aTmp =
dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
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))
{
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;
}
//=======================================================================
return;
}
-
-
MinEigenValue = CalculVP.Values() ( CalculVP.Values().Min());
- if ( MinEigenValue < CTol) {
- Convex = Standard_False;
- if (NoConvexTreatement) {
- Standard_Real Delta = CTol+0.1*Abs(MinEigenValue) -MinEigenValue ;
- for (ii=1; ii<=TheGradient.Length(); ii++) {
- TheHessian(ii, ii) += Delta;
- }
- }
- else {
- Done = Standard_False;
- TheStatus = math_FunctionError;
- return;
- }
- }
+ if ( MinEigenValue < CTol)
+ {
+ Convex = Standard_False;
+ if (NoConvexTreatement && // Treatment is allowed.
+ Abs (MinEigenValue) > CTol) // Treatment will have effect.
+ {
+ Standard_Real Delta = CTol + 0.1 * Abs(MinEigenValue) - MinEigenValue;
+ for (ii=1; ii<=TheGradient.Length(); ii++)
+ TheHessian(ii, ii) += Delta;
+ }
+ else
+ {
+ Done = Standard_False;
+ TheStatus = math_FunctionError;
+ return;
+ }
+ }
// Schemas de Newton
#include <Standard_DefineAlloc.hxx>
#include <Standard_Handle.hxx>
+#include <Precision.hxx>
#include <Standard_Boolean.hxx>
#include <math_Status.hxx>
#include <math_Vector.hxx>
//! positive (if the smaller eigenvalue of H < Convexity)
//! or IsConverged() returns True for 2 successives Iterations.
//! Warning: This constructor does not perform computation.
- Standard_EXPORT math_NewtonMinimum(const math_MultipleVarFunctionWithHessian& theFunction, const Standard_Real theTolerance = 1.0e-7, const Standard_Integer theNbIterations = 40, const Standard_Real theConvexity = 1.0e-6, const Standard_Boolean theWithSingularity = Standard_True);
+ Standard_EXPORT math_NewtonMinimum(const math_MultipleVarFunctionWithHessian& theFunction,
+ const Standard_Real theTolerance = Precision::Confusion(),
+ const Standard_Integer theNbIterations = 40,
+ const Standard_Real theConvexity = 1.0e-6,
+ const Standard_Boolean theWithSingularity = Standard_True);
//! Search the solution.
Standard_EXPORT void Perform (math_MultipleVarFunctionWithHessian& theFunction, const math_Vector& theStartingPoint);
*P = *Dir;
P->Multiply(x);
P->Add(*P0);
- F->Value(*P, fval);
- return Standard_True;
+
+ fval = 0.0;
+ return F->Value(*P, fval);
}
set tol_abs 7.e-5
set tol_rel 0.01
-#-1
+# Check result distance.
regexp "dist 1: +(\[-0-9.+eE\]+)" ${info} full dist_1
set expected_dist_1 0.
checkreal "Distance" ${dist_1} ${expected_dist_1} ${tol_abs} ${tol_rel}
-
-#-2
-set dump_list [dump ext_1]
-
-regexp { *Parameters *: *([-0-9.+eE]+) *([-0-9.+eE]+)} ${dump_list} full Parameter1 Parameter2
-regexp { *Origin *:([-0-9.+eE]+), *([-0-9.+eE]+) } ${dump_list} full OriginX OriginY
-regexp { *Axis *:([-0-9.+eE]+), *([-0-9.+eE]+) } ${dump_list} full AxisX AxisY
-
-set expected_Parameter1 0.
-checkreal "Parameter1" ${Parameter1} ${expected_Parameter1} ${tol_abs} ${tol_rel}
-
-set expected_Parameter2 0.
-checkreal "Parameter2" ${Parameter2} ${expected_Parameter2} ${tol_abs} ${tol_rel}
-
-set expected_OriginX 2.
-checkreal "OriginX" ${OriginX} ${expected_OriginX} ${tol_abs} ${tol_rel}
-
-set expected_OriginY 0.0
-checkreal "OriginY" ${OriginY} ${expected_OriginY} ${tol_abs} ${tol_rel}
-
-set expected_AxisX 1.
-checkreal "AxisX" ${AxisX} ${expected_AxisX} ${tol_abs} ${tol_rel}
-
-set expected_AxisY 0.
-checkreal "AxisY" ${AxisY} ${expected_AxisY} ${tol_abs} ${tol_rel}
set info [extrema r3 r4]
-regexp {Extrema 1 is point : } $info full pp
+regexp {Infinite number of extremas, distance = +([-0-9.+eE]+)} $info full dist
-if { $pp == "4 8 3"} {
- puts "Error : Point of extrema is wrong"
+if { $dist > 4.0e-13 } {
+ puts "Error : Extrema distance is too big"
} else {
- puts "OK: Point of extrema is valid"
+ puts "OK: Extrema distance is good"
}
--- /dev/null
+puts "============"
+puts "CR28175"
+puts "==========="
+puts ""
+###############################################################################
+# Bad result of curve-curve extrema
+###############################################################################
+
+pload MODELING
+
+# Prepare curves.
+restore [locate_data_file bug28175_border.brep] b
+restore [locate_data_file bug28175_segment.brep] s
+explode b
+explode s
+mkcurve cb b_6
+mkcurve cs s_1
+
+# Check result of forward arguments order.
+extrema cb cs 1
+set info [length ext_1]
+regexp {The length ext_1 is +([-0-9.+eE]+)} $info full len1
+if { $len1 > 1e-7 } {
+ puts "Error: extrema finds distance $len1 (parameters [dval prm_1_1] and [dval prm_1_3]"
+ puts "while should be 4.5985092243989615e-008 (parameters 23.548772710035486
+ and 585.69374860054825"
+} else {
+ puts "OK: Correct extrema point is found for forward arguments order"
+}
+
+# Check result of reversed arguments order.
+extrema cs cb 1
+set info [length ext_1]
+regexp {The length ext_1 is +([-0-9.+eE]+)} $info full len2
+if { $len2 > 1e-7 } {
+ puts "Error: extrema finds distance $len2 (parameters [dval prm_1_1] and [dval prm_1_3]"
+ puts "while should be 4.5985092243989615e-008 (parameters 23.548772710035486
+ and 585.69374860054825"
+} else {
+ puts "OK: Correct extrema point is found for reversed arguments order"
+}
+
+# Check that order not influence minimum value.
+if { abs ($len1 - $len2) > 1e-4 * ($len1 + $len2) } {
+ puts "Error: distance between cb to cs ($len1) is not equal to distance between cs and cb ($len2)"
+} else {
+ puts "OK: Extrema values are equal for forward and reversed arguments order"
+}
--- /dev/null
+puts "============"
+puts "CR28175"
+puts "==========="
+puts ""
+###############################################################################
+# Bad result of curve-curve extrema
+###############################################################################
+
+# Set signals on.
+pload MODELING
+dsetsignal 1
+
+# Prepare input data.
+restore [locate_data_file bug28175_2.brep] c
+explode c
+
+# Compute minimal distance
+distmini d c_1 c_2
+set dist [dval d_val]
+
+# Check extrema distance
+if { $dist > 1e-7 } {
+ puts "ERROR: Extrema distance is too big"
+} else {
+ puts "OK: Correct extrema distance"
+}
\ No newline at end of file
--- /dev/null
+puts "============"
+puts "CR28182"
+puts "==========="
+puts ""
+###############################################################################
+# BRepExtrema_DistShapeShape returns bad result of non-default deflection is used
+###############################################################################
+
+pload MODELING
+
+restore [locate_data_file bug28175_borders2.brep] b
+restore [locate_data_file bug28175_segments2_diff.brep] s
+explode s
+donly s_198 b
+
+set ref_nbsol 4
+set defl 0.0001
+
+set res [distmini r s_198 b $defl]
+
+set redges [lrange [lindex [split $res \n] 1] 1 end]
+set nbsol [llength $redges]
+set dist [dval r_val]
+
+don b s_198
+foreach q $redges { display $q; foreach v [explode $q] { display $v } }
+fit
+
+if { $dist > $defl } {
+ puts "Error: too big distance is reported: $dist"
+} else {
+ puts "OK: reported distance $dist is below $defl"
+}
+
+if {$nbsol != $ref_nbsol} {
+ puts "Error: $ref_nbsol solutions expected, but $nbsol found"
+} else {
+ puts "OK: $ref_nbsol solutions are found"
+}
--- /dev/null
+puts "============"
+puts "CR28183"
+puts "==========="
+puts ""
+###############################################################################
+# BRepExtrema_DistShapeShape does not find all solutions
+###############################################################################
+
+puts "TODO #28183 ALL: Error: .* solutions expected"
+
+pload MODELING
+
+restore [locate_data_file bug28175_borders2.brep] b
+restore [locate_data_file bug28175_segments2.brep] s
+explode b
+explode s
+
+set ref_nbsol 4
+set ref_dist 1e-7
+
+# find extremum points
+set res [distmini r s_511 b_2]
+set redges [lrange [lindex [split $res \n] 1] 1 end]
+set nbsol [llength $redges]
+set dist [dval r_val]
+
+# display curves and points
+don b_2 s_511
+foreach q $redges { display $q; foreach v [explode $q] { display $v } }
+
+if { $dist > $ref_dist } {
+ puts "Error: too big distance is reported: $dist"
+} else {
+ puts "OK: reported distance $dist is below $ref_dist"
+}
+
+if {$nbsol != $ref_nbsol} {
+ puts "Error: $ref_nbsol solutions expected, but $nbsol found"
+} else {
+ puts "OK: $ref_nbsol solutions are found"
+}
projects / occt.git / commitdiff