OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-03 17:56:21 +03:00
Code Packages Releases Activity

0028175: Bad result of curve-curve extrema

Browse Source 
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
This commit is contained in:
aml 2016-12-15 16:32:52 +03:00 committed by apn
parent 640d5fe219
commit f79b19a17e
17 changed files with 634 additions and 174 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/BRepTest/BRepTest_ExtremaCommands.cxx
View File

@ -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);

277
src/math/math_BFGS.cxx
View File

@ -27,6 +27,7 @@
#include <math_MultipleVarFunctionWithGradient.hxx>
#include <Standard_DimensionError.hxx>
#include <StdFail_NotDone.hxx>
#include <Precision.hxx>
#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;
}

17
src/math/math_BFGS.hxx
View File

@ -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:

151
src/math/math_BracketMinimum.cxx
View File

@ -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);
}

57
src/math/math_BracketMinimum.hxx
View File

@ -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;

41
src/math/math_BracketMinimum.lxx
View File

@ -14,6 +14,21 @@
// 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,
@ -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;
}

6
src/math/math_FRPR.cxx
View File

@ -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,

23
src/math/math_GlobOptMin.cxx
View File

@ -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<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))
@ -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;
}
//=======================================================================

33
src/math/math_NewtonMinimum.cxx
View File

@ -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

7
src/math/math_NewtonMinimum.hxx
View File

@ -21,6 +21,7 @@
#include <Standard_DefineAlloc.hxx>
#include <Standard_Handle.hxx>
#include <Precision.hxx>
#include <Standard_Boolean.hxx>
#include <math_Status.hxx>
#include <math_Vector.hxx>
@ -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);

5
src/math/math_Powell.cxx
View File

@ -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);
}

27
tests/bugs/fclasses/bug25635_1
View File

@ -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}

8
tests/bugs/modalg_5/bug23706_10
View File

@ -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"
}

48
tests/bugs/moddata_3/bug28175 Normal file
View File

@ -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"
}

26
tests/bugs/moddata_3/bug28175_1 Normal file
View File

@ -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"
}

39
tests/bugs/moddata_3/bug28182 Normal file
View File

@ -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"
}

41
tests/bugs/moddata_3/bug28183 Normal file
View File

@ -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"
}