From 075a5abe9b30b544f3ca51e2118b3fc84bdc3ad0 Mon Sep 17 00:00:00 2001 From: ifv Date: Wed, 12 Jul 2017 16:53:31 +0300 Subject: [PATCH] 0028856: Extrema between two curves gives wrong result, patch for 7.1.0 --- src/Extrema/Extrema_GenExtCC.gxx | 12 ++ src/math/math_BFGS.cxx | 287 ++++++++++++++++++++++++------- src/math/math_BFGS.hxx | 17 +- src/math/math_BracketMinimum.cxx | 145 ++++++++++++---- src/math/math_BracketMinimum.hxx | 55 ++++-- src/math/math_BracketMinimum.lxx | 41 +++++ src/math/math_GlobOptMin.cxx | 35 ++-- 7 files changed, 466 insertions(+), 126 deletions(-) diff --git a/src/Extrema/Extrema_GenExtCC.gxx b/src/Extrema/Extrema_GenExtCC.gxx index 63294ae1bd..250f610c15 100644 --- a/src/Extrema/Extrema_GenExtCC.gxx +++ b/src/Extrema/Extrema_GenExtCC.gxx @@ -211,6 +211,7 @@ void Extrema_GenExtCC::Perform() C2.Intervals(anIntervals2, aContinuity); // Lipchitz constant computation. + const Standard_Real aMaxLC = 10000.; Standard_Real aLC = 9.0; // Default value. const Standard_Real aMaxDer1 = 1.0 / C1.Resolution(1.0); const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0); @@ -220,6 +221,17 @@ void Extrema_GenExtCC::Perform() // Change constant value according to the concrete curve types. Standard_Boolean isConstLockedFlag = Standard_False; + //To prevent LipConst to became too small + const Standard_Real aCR = 0.001; + if (aMaxDer1 / aMaxDer < aCR || aMaxDer2 / aMaxDer < aCR) + { + isConstLockedFlag = Standard_True; + } + if (aMaxDer > aMaxLC) + { + aLC = aMaxLC; + isConstLockedFlag = Standard_True; + } if (C1.GetType() == GeomAbs_Line) { aMaxDer = 1.0 / C2.Resolution(1.0); 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 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, @@ -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; + 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); @@ -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,22 +103,22 @@ public: //! 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; @@ -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_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx index 8008ab33da..15107354c2 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; } //======================================================================= @@ -338,6 +343,10 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt //======================================================================= 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); @@ -369,10 +378,10 @@ void math_GlobOptMin::computeInitialValues() myC = myInitC; aLipConst *= Sqrt(myN) / aStep; - if (aLipConst < myC * 0.1) - myC = Max(aLipConst * 0.1, 0.01); - else if (aLipConst > myC * 5) - myC = Min(myC * 5, 50.0); + if (aLipConst < myC * aMinEps) + myC = Max(aLipConst * aMinEps, aMinLC); + else if (aLipConst > myC * aMaxEps) + myC = Min(myC * aMaxEps, aMaxLC); // Clear all solutions except one. if (myY.Size() != myN) -- 2.39.5