From f79b19a17e10f628ee0dbeb220ade1e98e27c28e Mon Sep 17 00:00:00 2001 From: aml Date: Thu, 15 Dec 2016 16:32:52 +0300 Subject: [PATCH] 0028175: Bad result of curve-curve extrema 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 --- src/BRepTest/BRepTest_ExtremaCommands.cxx | 2 +- src/math/math_BFGS.cxx | 277 +++++++++++++++++----- src/math/math_BFGS.hxx | 17 +- src/math/math_BracketMinimum.cxx | 151 +++++++++--- src/math/math_BracketMinimum.hxx | 57 +++-- src/math/math_BracketMinimum.lxx | 41 ++++ src/math/math_FRPR.cxx | 6 +- src/math/math_GlobOptMin.cxx | 23 +- src/math/math_NewtonMinimum.cxx | 33 +-- src/math/math_NewtonMinimum.hxx | 7 +- src/math/math_Powell.cxx | 5 +- tests/bugs/fclasses/bug25635_1 | 27 +-- tests/bugs/modalg_5/bug23706_10 | 8 +- tests/bugs/moddata_3/bug28175 | 48 ++++ tests/bugs/moddata_3/bug28175_1 | 26 ++ tests/bugs/moddata_3/bug28182 | 39 +++ tests/bugs/moddata_3/bug28183 | 41 ++++ 17 files changed, 634 insertions(+), 174 deletions(-) create mode 100644 tests/bugs/moddata_3/bug28175 create mode 100644 tests/bugs/moddata_3/bug28175_1 create mode 100644 tests/bugs/moddata_3/bug28182 create mode 100644 tests/bugs/moddata_3/bug28183 diff --git a/src/BRepTest/BRepTest_ExtremaCommands.cxx b/src/BRepTest/BRepTest_ExtremaCommands.cxx index 3fb1887f90..07c774014f 100644 --- a/src/BRepTest/BRepTest_ExtremaCommands.cxx +++ b/src/BRepTest/BRepTest_ExtremaCommands.cxx @@ -75,7 +75,7 @@ static Standard_Integer distmini(Draw_Interpretor& di, Standard_Integer n, const Standard_Real aDeflection = Precision::Confusion(); if (n == 5) - aDeflection = Draw::Atoi(a[4]); + aDeflection = Draw::Atof(a[4]); BRepExtrema_DistShapeShape dst(S1 ,S2, aDeflection); diff --git a/src/math/math_BFGS.cxx b/src/math/math_BFGS.cxx index 04076da132..6a5c450b34 100644 --- a/src/math/math_BFGS.cxx +++ b/src/math/math_BFGS.cxx @@ -27,6 +27,7 @@ #include #include #include +#include #define R 0.61803399 #define C (1.0-R) @@ -95,9 +96,9 @@ public : { *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, @@ -106,10 +107,14 @@ public : { *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) @@ -118,53 +123,197 @@ public : 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 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; + } + } - Standard_Real ax, xx, bx, Fax, Fxx, Fbx, F1; - F.Initialize(P, Dir); + // 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); + } - 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); + 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; - } - } - return Standard_False; - } + 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, @@ -201,8 +350,9 @@ 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; @@ -274,12 +424,20 @@ void math_BFGS::Perform(math_MultipleVarFunctionWithGradient& F, 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) + { } @@ -301,5 +459,14 @@ void math_BFGS::Perform(math_MultipleVarFunctionWithGradient& F, } } - - +//======================================================================= +//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; +} diff --git a/src/math/math_BFGS.hxx b/src/math/math_BFGS.hxx index 2a6f77c404..553e7c907f 100644 --- a/src/math/math_BFGS.hxx +++ b/src/math/math_BFGS.hxx @@ -36,6 +36,11 @@ class math_MultipleVarFunctionWithGradient; //! 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: @@ -51,8 +56,13 @@ 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 : @@ -120,6 +130,9 @@ protected: Standard_Real XTol; Standard_Real EPSZ; Standard_Integer nbiter; + Standard_Boolean myIsBoundsDefined; + math_Vector myLeft; + math_Vector myRight; private: diff --git a/src/math/math_BracketMinimum.cxx b/src/math/math_BracketMinimum.cxx index 3f1996c71e..4601b7ae0d 100644 --- a/src/math/math_BracketMinimum.cxx +++ b/src/math/math_BracketMinimum.cxx @@ -29,17 +29,38 @@ #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); @@ -53,19 +74,36 @@ 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; @@ -74,34 +112,54 @@ 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; - } - else if((Cx - u) * (u - ulim) > 0.0) { - 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); + 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; + } 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); @@ -114,33 +172,50 @@ 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); } diff --git a/src/math/math_BracketMinimum.hxx b/src/math/math_BracketMinimum.hxx index b92427bd1c..1d399aeb89 100644 --- a/src/math/math_BracketMinimum.hxx +++ b/src/math/math_BracketMinimum.hxx @@ -31,14 +31,20 @@ class math_Function; //! 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 @@ -64,9 +70,23 @@ public: //! 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 @@ -83,21 +103,21 @@ public: //! Is used to redefine the operator <<. Standard_EXPORT void Dump (Standard_OStream& o) const; - - - -protected: - - - //! Is used internally by the constructors. - Standard_EXPORT void Perform (math_Function& F, const Standard_Real A, const Standard_Real B); - - - - private: + //! Limit the given value to become within the range [myLeft, myRight]. + Standard_Real Limited(const Standard_Real theValue) const; + //! 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; @@ -106,6 +126,9 @@ private: Standard_Real FAx; Standard_Real FBx; Standard_Real FCx; + Standard_Real myLeft; + Standard_Real myRight; + Standard_Boolean myIsLimited; Standard_Boolean myFA; Standard_Boolean myFB; diff --git a/src/math/math_BracketMinimum.lxx b/src/math/math_BracketMinimum.lxx index 8edd5a722d..e7b656146f 100644 --- a/src/math/math_BracketMinimum.lxx +++ b/src/math/math_BracketMinimum.lxx @@ -14,6 +14,21 @@ // math_BracketMinimum.lxx +#include + +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, @@ -22,3 +37,29 @@ 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; +} diff --git a/src/math/math_FRPR.cxx b/src/math/math_FRPR.cxx index 80d459e860..26f0b5cb35 100644 --- a/src/math/math_FRPR.cxx +++ b/src/math/math_FRPR.cxx @@ -73,9 +73,9 @@ public : *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, diff --git a/src/math/math_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx index 8008ab33da..ace855134a 100644 --- a/src/math/math_GlobOptMin.cxx +++ b/src/math/math_GlobOptMin.cxx @@ -288,9 +288,11 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt { newtonMinimum.Location(theOutPnt); theVal = newtonMinimum.Minimum(); + + if (isInside(theOutPnt)) + return Standard_True; } - else return Standard_False; - } else + } // BFGS method used. if (myCont >= 1 && @@ -299,14 +301,18 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt math_MultipleVarFunctionWithGradient* aTmp = dynamic_cast (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(myFunc)) @@ -322,14 +328,13 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt { 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; } //======================================================================= diff --git a/src/math/math_NewtonMinimum.cxx b/src/math/math_NewtonMinimum.cxx index d8595833cd..dd326dc1a9 100644 --- a/src/math/math_NewtonMinimum.cxx +++ b/src/math/math_NewtonMinimum.cxx @@ -128,23 +128,24 @@ void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F, 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 diff --git a/src/math/math_NewtonMinimum.hxx b/src/math/math_NewtonMinimum.hxx index a344203a16..8ae683cfed 100644 --- a/src/math/math_NewtonMinimum.hxx +++ b/src/math/math_NewtonMinimum.hxx @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -47,7 +48,11 @@ public: //! 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); diff --git a/src/math/math_Powell.cxx b/src/math/math_Powell.cxx index e163070058..237c8521e1 100644 --- a/src/math/math_Powell.cxx +++ b/src/math/math_Powell.cxx @@ -75,8 +75,9 @@ Standard_Boolean DirFunctionBis::Value(const Standard_Real x, Standard_Real& fva *P = *Dir; P->Multiply(x); P->Add(*P0); - F->Value(*P, fval); - return Standard_True; + + fval = 0.0; + return F->Value(*P, fval); } diff --git a/tests/bugs/fclasses/bug25635_1 b/tests/bugs/fclasses/bug25635_1 index 0770219c53..995282f026 100755 --- a/tests/bugs/fclasses/bug25635_1 +++ b/tests/bugs/fclasses/bug25635_1 @@ -14,33 +14,8 @@ set info [2dextrema c1 c2] 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} diff --git a/tests/bugs/modalg_5/bug23706_10 b/tests/bugs/modalg_5/bug23706_10 index 78ce1cf761..2fc95e4d35 100755 --- a/tests/bugs/modalg_5/bug23706_10 +++ b/tests/bugs/modalg_5/bug23706_10 @@ -13,10 +13,10 @@ bsplinecurve r4 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 3 -1 2 3 1 1 11 3 1 3 9 3 1 3 9 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" } diff --git a/tests/bugs/moddata_3/bug28175 b/tests/bugs/moddata_3/bug28175 new file mode 100644 index 0000000000..468d3c4dae --- /dev/null +++ b/tests/bugs/moddata_3/bug28175 @@ -0,0 +1,48 @@ +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" +} diff --git a/tests/bugs/moddata_3/bug28175_1 b/tests/bugs/moddata_3/bug28175_1 new file mode 100644 index 0000000000..48ffa1b23d --- /dev/null +++ b/tests/bugs/moddata_3/bug28175_1 @@ -0,0 +1,26 @@ +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 diff --git a/tests/bugs/moddata_3/bug28182 b/tests/bugs/moddata_3/bug28182 new file mode 100644 index 0000000000..d4df323a49 --- /dev/null +++ b/tests/bugs/moddata_3/bug28182 @@ -0,0 +1,39 @@ +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" +} diff --git a/tests/bugs/moddata_3/bug28183 b/tests/bugs/moddata_3/bug28183 new file mode 100644 index 0000000000..5b52bc8579 --- /dev/null +++ b/tests/bugs/moddata_3/bug28183 @@ -0,0 +1,41 @@ +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" +}