From: azn Date: Thu, 22 Jan 2015 12:19:05 +0000 (+0300) Subject: 0025720: Incorrect code of math classes can lead to unpredicted behavior of algorithms X-Git-Tag: V6_9_0_beta~101 X-Git-Url: http://git.dev.opencascade.org/gitweb/?p=occt.git;a=commitdiff_plain;h=859a47c3d1434d63911586a6e035c5b6097943e0 0025720: Incorrect code of math classes can lead to unpredicted behavior of algorithms The calling of virtual methods has been removed from constructors & destructors: math_BissecNewton math_BrentMinimum math_FRPR math_FunctionSetRoot math_NewtonFunctionSetRoot math_NewtonMinimum math_Powell --- diff --git a/src/Bisector/Bisector_BisecCC.cxx b/src/Bisector/Bisector_BisecCC.cxx index 165d4c5dfb..fec6bca926 100644 --- a/src/Bisector/Bisector_BisecCC.cxx +++ b/src/Bisector/Bisector_BisecCC.cxx @@ -647,17 +647,23 @@ gp_Pnt2d Bisector_BisecCC::ValueAndDist (const Standard_Real U, if (Abs(FInit) < EpsH) { U2 = VInit; } - else { - math_BissecNewton SolNew (H,VMin - EpsH100,VMax + EpsH100,EpsH,10); - if (SolNew.IsDone()) { - U2 = SolNew.Root(); + else + { + math_BissecNewton aNewSolution(EpsH); + aNewSolution.Perform(H, VMin - EpsH100, VMax + EpsH100, 10); + + if (aNewSolution.IsDone()) + { + U2 = aNewSolution.Root(); } - else { + else + { math_FunctionRoot SolRoot (H,VInit,EpsH,VMin - EpsH100,VMax + EpsH100); - if (SolRoot.IsDone()) { + + if (SolRoot.IsDone()) U2 = SolRoot.Root(); - } - else { Valid = Standard_False;} + else + Valid = Standard_False; } } } diff --git a/src/Bisector/Bisector_Inter.cxx b/src/Bisector/Bisector_Inter.cxx index ab19c88af5..c7b5bc9480 100644 --- a/src/Bisector/Bisector_Inter.cxx +++ b/src/Bisector/Bisector_Inter.cxx @@ -346,12 +346,14 @@ void Bisector_Inter::NeighbourPerform(const Handle(Bisector_BisecCC)& Bis1, if (UMin - Eps > UMax + Eps) {return;} // Solution F = 0 to find the common point. - Bisector_FunctionInter Fint (Guide,Bis1,BisTemp); - math_BissecNewton Sol (Fint,UMin,UMax,Tol,20); - if (Sol.IsDone()) { - USol = Sol.Root(); - } - else { return; } + Bisector_FunctionInter Fint(Guide,Bis1,BisTemp); + + math_BissecNewton aSolution(Tol); + aSolution.Perform(Fint, UMin, UMax, 20); + + if (aSolution.IsDone()) + USol = aSolution.Root(); + else return; PSol = BisTemp ->ValueAndDist(USol,U1,U2,Dist); diff --git a/src/Extrema/Extrema_ExtPExtS.cxx b/src/Extrema/Extrema_ExtPExtS.cxx index f178a4d784..f6f22976d2 100644 --- a/src/Extrema/Extrema_ExtPExtS.cxx +++ b/src/Extrema/Extrema_ExtPExtS.cxx @@ -338,7 +338,8 @@ void Extrema_ExtPExtS::Perform (const gp_Pnt& P) Pe = ProjectPnt (anOrtogSection, myDirection, E), V = gp_Vec(E,Pe) * gp_Vec(myDirection); UV(1) = U; UV(2) = V; - math_FunctionSetRoot aFSR (myF,UV,Tol,UVinf,UVsup); + math_FunctionSetRoot aFSR (myF, Tol); + aFSR.Perform(myF, UV, UVinf, UVsup); // for (Standard_Integer k=1 ; k <= myF.NbExt(); Standard_Integer k; for ( k=1 ; k <= myF.NbExt(); k++) { diff --git a/src/Extrema/Extrema_GenExtCS.cxx b/src/Extrema/Extrema_GenExtCS.cxx index c5e32544b7..07d309b51b 100644 --- a/src/Extrema/Extrema_GenExtCS.cxx +++ b/src/Extrema/Extrema_GenExtCS.cxx @@ -289,7 +289,8 @@ void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C, math_PSO aPSO(&aFunc, TUVinf, TUVsup, aStep); aPSO.Perform(aParticles, aNbParticles, aValue, TUV); - math_FunctionSetRoot anA (myF, TUV, Tol, TUVinf, TUVsup, 100, Standard_False); + math_FunctionSetRoot anA(myF, Tol); + anA.Perform(myF, TUV, TUVinf, TUVsup); myDone = Standard_True; } diff --git a/src/Extrema/Extrema_GenExtPS.cxx b/src/Extrema/Extrema_GenExtPS.cxx index dc281f3f3b..97cf856d3f 100644 --- a/src/Extrema/Extrema_GenExtPS.cxx +++ b/src/Extrema/Extrema_GenExtPS.cxx @@ -778,8 +778,8 @@ void Extrema_GenExtPS::FindSolution(const gp_Pnt& /*P*/, UVsup(1) = myusup; UVsup(2) = myvsup; - const Standard_Integer aNbMaxIter = 100; - math_FunctionSetRoot S (myF, UV, Tol, UVinf, UVsup, aNbMaxIter); + math_FunctionSetRoot S(myF, Tol); + S.Perform(myF, UV, UVinf, UVsup); myDone = Standard_True; } diff --git a/src/Extrema/Extrema_GenExtSS.cxx b/src/Extrema/Extrema_GenExtSS.cxx index 33090e96e8..7fa8236c19 100644 --- a/src/Extrema/Extrema_GenExtSS.cxx +++ b/src/Extrema/Extrema_GenExtSS.cxx @@ -266,14 +266,16 @@ b- Calcul des minima: UV(3) = U20 + (N2Umin - 1) * PasU2; UV(4) = V20 + (N2Vmin - 1) * PasV2; - math_FunctionSetRoot SR1 (myF,UV,Tol,UVinf,UVsup); + math_FunctionSetRoot SR1(myF, Tol); + SR1.Perform(myF, UV, UVinf, UVsup); UV(1) = U10 + (N1Umax - 1) * PasU1; UV(2) = V10 + (N1Vmax - 1) * PasV1; UV(3) = U20 + (N2Umax - 1) * PasU2; UV(4) = V20 + (N2Vmax - 1) * PasV2; - math_FunctionSetRoot SR2 (myF,UV,Tol,UVinf,UVsup); + math_FunctionSetRoot SR2(myF, Tol); + SR2.Perform(myF, UV, UVinf, UVsup); myDone = Standard_True; } diff --git a/src/Extrema/Extrema_GenLocateExtCC.gxx b/src/Extrema/Extrema_GenLocateExtCC.gxx index 77ac8229d2..5473856390 100644 --- a/src/Extrema/Extrema_GenLocateExtCC.gxx +++ b/src/Extrema/Extrema_GenLocateExtCC.gxx @@ -82,7 +82,9 @@ Methode: Uusup(1)=Usup; Uusup(2)=Vsup; - math_FunctionSetRoot S (F,Start,Tol,Uuinf,Uusup); + math_FunctionSetRoot S(F, Tol); + S.Perform(F, Start, Uuinf, Uusup); + if (S.IsDone() && F.NbExt() > 0) { mySqDist = F.SquareDistance(1); F.Points(1,myPoint1,myPoint2); diff --git a/src/Extrema/Extrema_GenLocateExtCS.cxx b/src/Extrema/Extrema_GenLocateExtCS.cxx index 4195f8fdfb..35719a403d 100644 --- a/src/Extrema/Extrema_GenLocateExtCS.cxx +++ b/src/Extrema/Extrema_GenLocateExtCS.cxx @@ -89,7 +89,8 @@ void Extrema_GenLocateExtCS::Perform(const Adaptor3d_Curve& C, BSup(2) = Usup; BSup(3) = Vsup; - math_FunctionSetRoot SR (F, Start,Tol, BInf, BSup); + math_FunctionSetRoot SR (F, Tol); + SR.Perform(F, Start, BInf, BSup); if (!SR.IsDone()) return; diff --git a/src/Extrema/Extrema_GenLocateExtPS.cxx b/src/Extrema/Extrema_GenLocateExtPS.cxx index 50afd94b68..4227ee9b87 100644 --- a/src/Extrema/Extrema_GenLocateExtPS.cxx +++ b/src/Extrema/Extrema_GenLocateExtPS.cxx @@ -75,7 +75,8 @@ Method: BSup(1) = Usup; BSup(2) = Vsup; - math_FunctionSetRoot SR (F, Start,Tol, BInf, BSup); + math_FunctionSetRoot SR (F, Tol); + SR.Perform(F, Start, BInf, BSup); if (!SR.IsDone()) return; diff --git a/src/Extrema/Extrema_GenLocateExtSS.cxx b/src/Extrema/Extrema_GenLocateExtSS.cxx index eb8cad3333..320f5eaa3b 100644 --- a/src/Extrema/Extrema_GenLocateExtSS.cxx +++ b/src/Extrema/Extrema_GenLocateExtSS.cxx @@ -97,7 +97,8 @@ void Extrema_GenLocateExtSS::Perform(const Adaptor3d_Surface& S1, BSup(3) = Usup2; BSup(4) = Vsup2; - math_FunctionSetRoot SR (F, Start,Tol, BInf, BSup); + math_FunctionSetRoot SR (F, Tol); + SR.Perform(F, Start, BInf, BSup); if (!SR.IsDone()) return; diff --git a/src/FairCurve/FairCurve_Newton.cdl b/src/FairCurve/FairCurve_Newton.cdl index b39e5eee60..1a62b0c46c 100644 --- a/src/FairCurve/FairCurve_Newton.cdl +++ b/src/FairCurve/FairCurve_Newton.cdl @@ -23,50 +23,30 @@ uses MultipleVarFunctionWithHessian from math is - Create(F: in out MultipleVarFunctionWithHessian; - StartingPoint: Vector; - SpatialTolerance : Real = 1.0e-7; - CriteriumTolerance: Real = 1.0e-2; - NbIterations: Integer=40; - Convexity: Real=1.0e-6; - WithSingularity : Boolean = Standard_True) - ---Purpose: -- Given the starting point StartingPoint, - -- The tolerance required on the solution is given by - -- Tolerance. - -- Iteration are stopped if - -- (!WithSingularity) and H(F(Xi)) is not definite - -- positive (if the smaller eigenvalue of H < Convexity) - -- or IsConverged() returns True for 2 successives Iterations. - -- Warning: Obsolete Constructor (because IsConverged can not be redefined - -- with this. ) - returns Newton; - - Create(F: in out MultipleVarFunctionWithHessian; - SpatialTolerance : Real = 1.0e-7; - Tolerance: Real=1.0e-7; - NbIterations: Integer=40; - Convexity: Real=1.0e-6; - WithSingularity : Boolean = Standard_True) + Create(theFunction: in MultipleVarFunctionWithHessian; + theSpatialTolerance: Real = 1.0e-7; + theCriteriumTolerance: Real=1.0e-7; + theNbIterations: Integer=40; + theConvexity: Real=1.0e-6; + theWithSingularity: Boolean = Standard_True) ---Purpose: - -- The tolerance required on the solution is given by - -- Tolerance. - -- Iteration are stopped if - -- (!WithSingularity) and H(F(Xi)) is not definite - -- positive (if the smaller eigenvalue of H < Convexity) - -- or IsConverged() returns True for 2 successives Iterations. - -- Warning: This constructor do not computation + -- The tolerance required on the solution is given by Tolerance. + -- Iteration are stopped if (!WithSingularity) and H(F(Xi)) is not definite + -- positive (if the smaller eigenvalue of H < Convexity) + -- or IsConverged() returns True for 2 successives Iterations. + -- Warning: This constructor do not computation returns Newton; IsConverged(me) - ---Purpose: This method is called at the end of each - -- iteration to check the convergence : - -- || Xi+1 - Xi || < SpatialTolerance/100 Or - -- || Xi+1 - Xi || < SpatialTolerance and - -- |F(Xi+1) - F(Xi)| < CriteriumTolerance * |F(xi)| - -- It can be redefined in a sub-class to implement a specific test. - - returns Boolean - is redefined; + ---Purpose: + -- This method is called at the end of each + -- iteration to check the convergence : + -- || Xi+1 - Xi || < SpatialTolerance/100 Or + -- || Xi+1 - Xi || < SpatialTolerance and + -- |F(Xi+1) - F(Xi)| < CriteriumTolerance * |F(xi)| + -- It can be redefined in a sub-class to implement a specific test. + + returns Boolean is redefined; fields mySpTol : Real; diff --git a/src/FairCurve/FairCurve_Newton.cxx b/src/FairCurve/FairCurve_Newton.cxx index 53100205e2..bee99ce134 100644 --- a/src/FairCurve/FairCurve_Newton.cxx +++ b/src/FairCurve/FairCurve_Newton.cxx @@ -16,44 +16,37 @@ #include -FairCurve_Newton::FairCurve_Newton(math_MultipleVarFunctionWithHessian& F, - const math_Vector& StartingPoint, - const Standard_Real SpatialTolerance, - const Standard_Real CriteriumTolerance, - const Standard_Integer NbIterations, - const Standard_Real Convexity, - const Standard_Boolean WithSingularity) - :math_NewtonMinimum(F, StartingPoint, CriteriumTolerance, - NbIterations, Convexity, WithSingularity), - mySpTol(SpatialTolerance) - +//======================================================================= +//function : FairCurve_Newton +//purpose : Constructor +//======================================================================= +FairCurve_Newton::FairCurve_Newton( + const math_MultipleVarFunctionWithHessian& theFunction, + const Standard_Real theSpatialTolerance, + const Standard_Real theCriteriumTolerance, + const Standard_Integer theNbIterations, + const Standard_Real theConvexity, + const Standard_Boolean theWithSingularity + ) +: math_NewtonMinimum(theFunction, + theCriteriumTolerance, + theNbIterations, + theConvexity, + theWithSingularity), + mySpTol(theSpatialTolerance) { -// Attention this writing is wrong as FairCurve_Newton::IsConverged() is not -// used in the constructor of NewtonMinimum !! -} - -FairCurve_Newton::FairCurve_Newton(math_MultipleVarFunctionWithHessian& F, - const Standard_Real SpatialTolerance, - const Standard_Real CriteriumTolerance, - const Standard_Integer NbIterations, - const Standard_Real Convexity, - const Standard_Boolean WithSingularity) - :math_NewtonMinimum(F, CriteriumTolerance, - NbIterations, Convexity, WithSingularity), - mySpTol(SpatialTolerance) - -{ -// It is much better } +//======================================================================= +//function : IsConverged +//purpose : Convert if the steps are too small or if the criterion +// progresses little with a reasonable step, this last +// requirement allows detecting infinite slidings +// (case when the criterion varies troo slowly). +//======================================================================= Standard_Boolean FairCurve_Newton::IsConverged() const -// Convert if the steps are too small -// or if the criterion progresses little with a reasonable step, this last requirement -// allows detecting infinite slidings, -// (case when the criterion varies troo slowly). { - Standard_Real N = TheStep.Norm(); - return ( (N <= mySpTol/100 ) || - ( Abs(TheMinimum-PreviousMinimum) <= XTol*Abs(PreviousMinimum) - && N<=mySpTol) ); + const Standard_Real N = TheStep.Norm(); + return ( N <= 0.01 * mySpTol ) || ( N <= mySpTol && + Abs(TheMinimum - PreviousMinimum) <= XTol * Abs(PreviousMinimum)); } diff --git a/src/Geom2dGcc/Geom2dGcc_Circ2d2TanOnIter.cxx b/src/Geom2dGcc/Geom2dGcc_Circ2d2TanOnIter.cxx index 66ba352ef0..626be8a949 100644 --- a/src/Geom2dGcc/Geom2dGcc_Circ2d2TanOnIter.cxx +++ b/src/Geom2dGcc/Geom2dGcc_Circ2d2TanOnIter.cxx @@ -91,7 +91,8 @@ Geom2dGcc_Circ2d2TanOnIter (const GccEnt_QualifiedLin& Qualified1 , gp_Pnt2d point3 = ElCLib::Value(Param3,OnLine); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(L1,Cu2,OnLine,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -192,7 +193,8 @@ Geom2dGcc_Circ2d2TanOnIter (const Geom2dGcc_QCurve& Qualified1 , gp_Pnt2d point3 = ElCLib::Value(Param3,OnLine); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(Cu1,Cu2,OnLine,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -288,7 +290,8 @@ Geom2dGcc_Circ2d2TanOnIter (const Geom2dGcc_QCurve& Qualified1 , gp_Pnt2d point3 = ElCLib::Value(Param2,OnLine); Ufirst(3) = (point3.Distance(Point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(Cu1,Point2,OnLine,Ufirst(3)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -381,7 +384,8 @@ Geom2dGcc_Circ2d2TanOnIter (const GccEnt_QualifiedCirc& Qualified1 , gp_Pnt2d point3 = ElCLib::Value(Param3,OnLine); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(C1,Cu2,OnLine,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -488,7 +492,8 @@ Geom2dGcc_Circ2d2TanOnIter (const GccEnt_QualifiedCirc& Qualified1 , gp_Pnt2d point3 = ElCLib::Value(Param3,OnCirc); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(C1,Cu2,OnCirc,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -592,7 +597,8 @@ Geom2dGcc_Circ2d2TanOnIter (const GccEnt_QualifiedLin& Qualified1 , gp_Pnt2d point3 = ElCLib::Value(Param3,OnCirc); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(L1,Cu2,OnCirc,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -696,7 +702,8 @@ Geom2dGcc_Circ2d2TanOnIter (const Geom2dGcc_QCurve& Qualified1 , gp_Pnt2d point3(OnCirc.Location().XY()+R1*gp_XY(Cos(Param3),Sin(Param3))); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(Cu1,Cu2,OnCirc,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -796,7 +803,8 @@ Geom2dGcc_Circ2d2TanOnIter (const Geom2dGcc_QCurve& Qualified1 , gp_Pnt2d point3 = ElCLib::Value(Param2,OnCirc); Ufirst(3) = (point3.Distance(Point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(Cu1,Point2,OnCirc,Ufirst(3)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -889,7 +897,8 @@ Geom2dGcc_Circ2d2TanOnIter (const Geom2dGcc_QCurve& Qualified1 , gp_Pnt2d point3 = Geom2dGcc_CurveTool::Value(OnCurv,Param3); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(Cu1,Cu2,OnCurv,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -994,7 +1003,8 @@ Geom2dGcc_Circ2d2TanOnIter (const GccEnt_QualifiedCirc& Qualified1 , gp_Pnt2d point3 = Geom2dGcc_CurveTool::Value(OnCurv,ParamOn); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(C1,Cu2,OnCurv,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -1095,7 +1105,8 @@ Geom2dGcc_Circ2d2TanOnIter (const GccEnt_QualifiedLin& Qualified1 , gp_Pnt2d point3 = Geom2dGcc_CurveTool::Value(OnCurv,ParamOn); Ufirst(4) = (point3.Distance(point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(L1,Cu2,OnCurv,Ufirst(4)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); @@ -1186,7 +1197,8 @@ Geom2dGcc_Circ2d2TanOnIter (const Geom2dGcc_QCurve& Qualified1 , gp_Pnt2d point3 = Geom2dGcc_CurveTool::Value(OnCurv,ParamOn); Ufirst(3) = (point3.Distance(Point2)+point3.Distance(point1))/2.; Geom2dGcc_FunctionTanCuCuOnCu Func(Cu1,Point2,OnCurv,Ufirst(3)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); Func.Value(Ufirst,Umin); if (Root.IsDone()) { Root.Root(Ufirst); diff --git a/src/Geom2dGcc/Geom2dGcc_Circ2d3TanIter.cxx b/src/Geom2dGcc/Geom2dGcc_Circ2d3TanIter.cxx index f50338005a..7a3a0433c9 100644 --- a/src/Geom2dGcc/Geom2dGcc_Circ2d3TanIter.cxx +++ b/src/Geom2dGcc/Geom2dGcc_Circ2d3TanIter.cxx @@ -89,7 +89,8 @@ Geom2dGcc_Circ2d3TanIter (const Geom2dGcc_QCurve& Qualified1 , tol(1) = Geom2dGcc_CurveTool::EpsX(Cu1,Abs(Tolerance)); tol(2) = Geom2dGcc_CurveTool::EpsX(Cu2,Abs(Tolerance)); tol(3) = Geom2dGcc_CurveTool::EpsX(Cu3,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); @@ -212,7 +213,8 @@ Geom2dGcc_Circ2d3TanIter (const GccEnt_QualifiedCirc& Qualified1 , tol(1) = 2.e-15*M_PI; tol(2) = Geom2dGcc_CurveTool::EpsX(Cu2,Abs(Tolerance)); tol(3) = Geom2dGcc_CurveTool::EpsX(Cu3,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); @@ -340,7 +342,8 @@ Geom2dGcc_Circ2d3TanIter (const GccEnt_QualifiedCirc& Qualified1 , tol(1) = 2.e-15*M_PI; tol(2) = 2.e-15*M_PI; tol(3) = Geom2dGcc_CurveTool::EpsX(Cu3,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); @@ -471,7 +474,8 @@ Geom2dGcc_Circ2d3TanIter (const GccEnt_QualifiedLin& Qualified1 , tol(1) = 1.e-15; tol(2) = Geom2dGcc_CurveTool::EpsX(Cu2,Abs(Tolerance)); tol(3) = Geom2dGcc_CurveTool::EpsX(Cu3,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); @@ -602,7 +606,8 @@ Geom2dGcc_Circ2d3TanIter (const GccEnt_QualifiedLin& Qualified1 , tol(1) = 1.e-15; tol(2) = 1.e-15; tol(3) = Geom2dGcc_CurveTool::EpsX(Cu3,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); @@ -726,7 +731,8 @@ Geom2dGcc_Circ2d3TanIter (const Geom2dGcc_QCurve& Qualified1 , tol(1) = 2.e-15*M_PI; tol(2) = Geom2dGcc_CurveTool::EpsX(Cu1,Abs(Tolerance)); tol(3) = Geom2dGcc_CurveTool::EpsX(Cu2,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); @@ -839,7 +845,8 @@ Geom2dGcc_Circ2d3TanIter (const Geom2dGcc_QCurve& Qualified1 , tol(1) = 2.e-15*M_PI; tol(2) = 2.e-15*M_PI; tol(3) = Geom2dGcc_CurveTool::EpsX(Cu1,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); @@ -949,7 +956,8 @@ Geom2dGcc_Circ2d3TanIter (const GccEnt_QualifiedLin& Qualified1 , tol(1) = 2.e-15; tol(2) = 1.e-15; tol(3) = Geom2dGcc_CurveTool::EpsX(Cu2,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); @@ -1068,7 +1076,8 @@ Geom2dGcc_Circ2d3TanIter (const GccEnt_QualifiedCirc& Qualified1 , tol(1) = 2.e-15*M_PI; tol(2) = 1.e-15; tol(3) = Geom2dGcc_CurveTool::EpsX(Cu3,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Func.Value(Ufirst,Umin); Root.Root(Ufirst); @@ -1195,7 +1204,8 @@ Geom2dGcc_Circ2d3TanIter (const GccEnt_QualifiedCirc& Qualified1 , tol(1) = 2.e-15*M_PI; tol(2) = 2.e-15*M_PI; tol(3) = Geom2dGcc_CurveTool::EpsX(Cu2,Abs(Tolerance)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); Func.Value(Ufirst,Umin); diff --git a/src/Geom2dGcc/Geom2dGcc_Lin2d2TanIter.cxx b/src/Geom2dGcc/Geom2dGcc_Lin2d2TanIter.cxx index b0e28a0ce4..cfc8efa378 100644 --- a/src/Geom2dGcc/Geom2dGcc_Lin2d2TanIter.cxx +++ b/src/Geom2dGcc/Geom2dGcc_Lin2d2TanIter.cxx @@ -144,7 +144,8 @@ Geom2dGcc_Lin2d2TanIter (const Geom2dGcc_QCurve& Qualified1 , Ufirst(2) = Param2; tol(1) = Geom2dGcc_CurveTool::EpsX(Cu1,Abs(Tolang)); tol(2) = Geom2dGcc_CurveTool::EpsX(Cu2,Abs(Tolang)); - math_FunctionSetRoot Root(Func,Ufirst,tol,Umin,Umax); + math_FunctionSetRoot Root(Func, tol); + Root.Perform(Func, Ufirst, Umin, Umax); if (Root.IsDone()) { Root.Root(Ufirst); // Modified by Sergey KHROMOV - Thu Apr 5 17:45:00 2001 Begin diff --git a/src/GeomFill/GeomFill_LocationDraft.cxx b/src/GeomFill/GeomFill_LocationDraft.cxx index 46d96af53b..e5b1b52a0d 100644 --- a/src/GeomFill/GeomFill_LocationDraft.cxx +++ b/src/GeomFill/GeomFill_LocationDraft.cxx @@ -310,30 +310,27 @@ GeomFill_LocationDraft::GeomFill_LocationDraft GeomFill_FunctionDraft E(mySurf , G); // resolution - math_NewtonFunctionSetRoot Result(E, X, XTol, FTol, Iter); - + math_NewtonFunctionSetRoot Result(E, XTol, FTol, Iter); + Result.Perform(E, X); + if (Result.IsDone()) - { - math_Vector R(1,3); - Result.Root(R); // solution - - gp_Pnt2d p (R(2), R(3)); // point sur la surface - gp_Pnt2d q (R(1), Param); // point de la courbe - Poles2d.SetValue(1,p); - Poles2d.SetValue(2,q); - - } + { + math_Vector R(1,3); + Result.Root(R); // solution + + gp_Pnt2d p (R(2), R(3)); // point sur la surface + gp_Pnt2d q (R(1), Param); // point de la courbe + Poles2d.SetValue(1,p); + Poles2d.SetValue(2,q); + } else { - return Standard_False; + return Standard_False; } - + }// if_Intersec - - }// if_Intersec - -// la generatrice n'intersecte pas la surface d'arret + // la generatrice n'intersecte pas la surface d'arret return Standard_True; -} + } //================================================================== //Function: D1 @@ -427,44 +424,44 @@ GeomFill_LocationDraft::GeomFill_LocationDraft // resolution - math_NewtonFunctionSetRoot Result (E, X, XTol, FTol, Iter); + math_NewtonFunctionSetRoot Result(E, XTol, FTol, Iter); + Result.Perform(E, X); if (Result.IsDone()) - { - math_Vector R(1,3); - Result.Root(R); // solution - - gp_Pnt2d p (R(2), R(3)); // point sur la surface - gp_Pnt2d q (R(1), Param); // point de la courbe - Poles2d.SetValue(1,p); - Poles2d.SetValue(2,q); - - // derivee de la fonction par rapport a Param - math_Vector DEDT(1,3,0); - E.DerivT(myTrimmed, Param, R(1), DN, myAngle, DEDT); // dE/dt => DEDT - - math_Vector DSDT (1,3,0); - math_Matrix DEDX (1,3,1,3,0); - E.Derivatives(R, DEDX); // dE/dx au point R => DEDX - - // resolution du syst. lin. : DEDX*DSDT = -DEDT - math_Gauss Ga(DEDX); - if (Ga.IsDone()) - { - Ga.Solve (DEDT.Opposite(), DSDT); // resolution du syst. lin. - gp_Vec2d dp (DSDT(2), DSDT(3)); // surface - gp_Vec2d dq (DSDT(1), 1); // courbe - DPoles2d.SetValue(1, dp); - DPoles2d.SetValue(2, dq); - }//if - - - }//if_Result + { + math_Vector R(1,3); + Result.Root(R); // solution + + gp_Pnt2d p (R(2), R(3)); // point sur la surface + gp_Pnt2d q (R(1), Param); // point de la courbe + Poles2d.SetValue(1,p); + Poles2d.SetValue(2,q); + + // derivee de la fonction par rapport a Param + math_Vector DEDT(1,3,0); + E.DerivT(myTrimmed, Param, R(1), DN, myAngle, DEDT); // dE/dt => DEDT + + math_Vector DSDT (1,3,0); + math_Matrix DEDX (1,3,1,3,0); + E.Derivatives(R, DEDX); // dE/dx au point R => DEDX + + // resolution du syst. lin. : DEDX*DSDT = -DEDT + math_Gauss Ga(DEDX); + if (Ga.IsDone()) + { + Ga.Solve (DEDT.Opposite(), DSDT); // resolution du syst. lin. + gp_Vec2d dp (DSDT(2), DSDT(3)); // surface + gp_Vec2d dq (DSDT(1), 1); // courbe + DPoles2d.SetValue(1, dp); + DPoles2d.SetValue(2, dq); + }//if + + + }//if_Result else {// la generatrice n'intersecte pas la surface d'arret - return Standard_False; + return Standard_False; } - }// if_Intersec - + }// if_Intersec return Standard_True; } @@ -564,7 +561,8 @@ GeomFill_LocationDraft::GeomFill_LocationDraft // resolution - math_NewtonFunctionSetRoot Result (E, X, XTol, FTol, Iter); + math_NewtonFunctionSetRoot Result (E, XTol, FTol, Iter); + Result.Perform(E, X); if (Result.IsDone()) { diff --git a/src/GeomFill/GeomFill_LocationGuide.cxx b/src/GeomFill/GeomFill_LocationGuide.cxx index 86241c8b7f..b26d958a1b 100644 --- a/src/GeomFill/GeomFill_LocationGuide.cxx +++ b/src/GeomFill/GeomFill_LocationGuide.cxx @@ -349,34 +349,34 @@ static void InGoodPeriod(const Standard_Real Prec, #endif Standard_Boolean SOS=Standard_False; if (ii>1) { - // Intersection de secour entre surf revol et guide - // equation - X(1) = myPoles2d->Value(1,ii-1).Y(); - X(2) = myPoles2d->Value(2,ii-1).X(); - X(3) = myPoles2d->Value(2,ii-1).Y(); - GeomFill_FunctionGuide E (mySec, myGuide, U); - E.SetParam(U, P, T.XYZ(), N.XYZ()); - // resolution => angle - math_FunctionSetRoot Result(E, X, TolRes, - Inf, Sup); - - if (Result.IsDone() && - (Result.FunctionSetErrors().Norm() < TolRes(1)*TolRes(1)) ) { + // Intersection de secour entre surf revol et guide + // equation + X(1) = myPoles2d->Value(1,ii-1).Y(); + X(2) = myPoles2d->Value(2,ii-1).X(); + X(3) = myPoles2d->Value(2,ii-1).Y(); + GeomFill_FunctionGuide E (mySec, myGuide, U); + E.SetParam(U, P, T.XYZ(), N.XYZ()); + // resolution => angle + math_FunctionSetRoot Result(E, TolRes); + Result.Perform(E, X, Inf, Sup); + + if (Result.IsDone() && + (Result.FunctionSetErrors().Norm() < TolRes(1)*TolRes(1)) ) { #ifdef OCCT_DEBUG - cout << "Ratrappage Reussi !" << endl; + cout << "Ratrappage Reussi !" << endl; #endif - SOS = Standard_True; - math_Vector RR(1,3); - Result.Root(RR); - PInt.SetValues(P, RR(2), RR(3), RR(1), IntCurveSurface_Out); - theU = PInt.U(); - theV = PInt.V(); - } - else { + SOS = Standard_True; + math_Vector RR(1,3); + Result.Root(RR); + PInt.SetValues(P, RR(2), RR(3), RR(1), IntCurveSurface_Out); + theU = PInt.U(); + theV = PInt.V(); + } + else { #ifdef OCCT_DEBUG - cout << "Echec du Ratrappage !" << endl; + cout << "Echec du Ratrappage !" << endl; #endif - } + } } if (!SOS) { myStatus = GeomFill_ImpossibleContact; @@ -614,8 +614,8 @@ static void InGoodPeriod(const Standard_Real Prec, GeomFill_FunctionGuide E (mySec, myGuide, U); E.SetParam(Param, P, t, n); // resolution => angle - math_FunctionSetRoot Result(E, X, TolRes, - Inf, Sup, Iter); + math_FunctionSetRoot Result(E, TolRes, Iter); + Result.Perform(E, X, Inf, Sup); if (Result.IsDone()) { // solution @@ -682,11 +682,11 @@ static void InGoodPeriod(const Standard_Real Prec, GeomFill_FunctionGuide E (mySec, myGuide, myFirstS + (Param-myCurve->FirstParameter())*ratio); E.SetParam(Param, P, t, n); - + // resolution - math_FunctionSetRoot Result(E, X, TolRes, - Inf, Sup, Iter); - + math_FunctionSetRoot Result(E, TolRes, Iter); + Result.Perform(E, X, Inf, Sup); + if (Result.IsDone()) { // solution Result.Root(R); diff --git a/src/IntCurve/IntCurve_ExactIntersectionPoint.gxx b/src/IntCurve/IntCurve_ExactIntersectionPoint.gxx index 5318961892..9159682d5d 100644 --- a/src/IntCurve/IntCurve_ExactIntersectionPoint.gxx +++ b/src/IntCurve/IntCurve_ExactIntersectionPoint.gxx @@ -198,12 +198,8 @@ void IntCurve_ExactIntersectionPoint::Roots(Standard_Real& U,Standard_Real& V) { void IntCurve_ExactIntersectionPoint::MathPerform(void) { - math_FunctionSetRoot Fct( FctDist - ,StartingPoint - ,ToleranceVector - ,BInfVector - ,BSupVector - ,60); + math_FunctionSetRoot Fct(FctDist, ToleranceVector, 60); + Fct.Perform(FctDist, StartingPoint, BInfVector, BSupVector); if(Fct.IsDone()) { Fct.Root(Root); nbroots = 1; diff --git a/src/ProjLib/ProjLib_PrjResolve.cxx b/src/ProjLib/ProjLib_PrjResolve.cxx index a8322021e6..8288888848 100644 --- a/src/ProjLib/ProjLib_PrjResolve.cxx +++ b/src/ProjLib/ProjLib_PrjResolve.cxx @@ -72,15 +72,17 @@ ProjLib_PrjResolve::ProjLib_PrjResolve(const Adaptor3d_Curve& C,const Adaptor3d_ // if (!S1.IsDone()) { return; } // } // else { - math_NewtonFunctionSetRoot SR (F, Tol, 1.e-10); - SR.Perform(F, Start, BInf, BSup); + math_NewtonFunctionSetRoot SR (F, Tol, 1.e-10); + SR.Perform(F, Start, BInf, BSup); // if (!SR.IsDone()) { return; } - if (!SR.IsDone()) { - math_FunctionSetRoot S1 (F, Start,Tol, BInf, BSup); - if (!S1.IsDone()) { return; } - } - + if (!SR.IsDone()) + { + math_FunctionSetRoot S1 (F, Tol); + S1.Perform(F, Start, BInf, BSup); + if (!S1.IsDone()) + return; + } mySolution.SetXY(F.Solution().XY()); diff --git a/src/math/math_BissecNewton.cdl b/src/math/math_BissecNewton.cdl index e32d1ff143..7e8f96b857 100644 --- a/src/math/math_BissecNewton.cdl +++ b/src/math/math_BissecNewton.cdl @@ -30,41 +30,34 @@ raises NotDone from StdFail is + Create(theXTolerance: in Real) returns BissecNewton; + ---Purpose: Constructor. + -- @param theXTolerance - algorithm tolerance. Perform(me: in out; F: out FunctionWithDerivative; - Bound1, Bound2: Real; - NbIterations: Integer) - - is static protected; - - - Create(F: in out FunctionWithDerivative; - Bound1, Bound2, TolX: Real; - NbIterations: Integer = 100) - ---Purpose: - -- A combination of Newton-Raphson and bissection methods is done to find - -- the root of the function F between the bounds Bound1 and Bound2. - -- on the function F. - -- The tolerance required on the root is given by TolX. - -- The solution is found when : - -- abs(Xi - Xi-1) <= TolX and F(Xi) * F(Xi-1) <= 0 - -- The maximum number of iterations allowed is given by NbIterations. - - returns BissecNewton; + Bound1, Bound2: Real; + NbIterations: Integer = 100) + ---Purpose: + -- A combination of Newton-Raphson and bissection methods is done to find + -- the root of the function F between the bounds Bound1 and Bound2 + -- on the function F. + -- The tolerance required on the root is given by TolX. + -- The solution is found when: + -- abs(Xi - Xi-1) <= TolX and F(Xi) * F(Xi-1) <= 0 + -- The maximum number of iterations allowed is given by NbIterations. + is static; - ---C++: alias " Standard_EXPORT virtual ~math_BissecNewton();" - IsSolutionReached(me: in out; F: out FunctionWithDerivative) - ---Purpose: - -- This method is called at the end of each iteration to check if the - -- solution has been found. - -- It can be redefined in a sub-class to implement a specific test to - -- stop the iterations. + IsSolutionReached(me: in out; theFunction: out FunctionWithDerivative) + ---Purpose: + -- This method is called at the end of each iteration to check if the + -- solution has been found. + -- It can be redefined in a sub-class to implement a specific test to + -- stop the iterations. + ---C++: inline + returns Boolean is virtual; - returns Boolean - is virtual; - IsDone(me) ---Purpose: Tests is the root has been successfully found. ---C++: inline @@ -99,7 +92,7 @@ is raises NotDone is static; - + Dump(me; o: in out OStream) ---Purpose: Prints on the stream o information on the current state -- of the object. @@ -108,6 +101,12 @@ is is static; + Delete(me) is static; + ---Purpose: Destructor alias. + ---C++: inline + ---C++: alias " Standard_EXPORT virtual ~math_BissecNewton();" + + fields Done: Boolean; diff --git a/src/math/math_BissecNewton.cxx b/src/math/math_BissecNewton.cxx index b2d0c71bf1..884deec0c4 100644 --- a/src/math/math_BissecNewton.cxx +++ b/src/math/math_BissecNewton.cxx @@ -15,14 +15,37 @@ #include #include +//======================================================================= +//function : math_BissecNewton +//purpose : Constructor +//======================================================================= +math_BissecNewton::math_BissecNewton(const Standard_Real theXTolerance) +: TheStatus(math_NotBracketed), + XTol (theXTolerance), + x (0.0), + dx (0.0), + f (0.0), + df (0.0), + Done (Standard_False) +{ +} + +//======================================================================= +//function : ~math_BissecNewton +//purpose : Destructor +//======================================================================= math_BissecNewton::~math_BissecNewton() { } +//======================================================================= +//function : Perform +//purpose : +//======================================================================= void math_BissecNewton::Perform(math_FunctionWithDerivative& F, - const Standard_Real Bound1, - const Standard_Real Bound2, - const Standard_Integer NbIterations) + const Standard_Real Bound1, + const Standard_Real Bound2, + const Standard_Integer NbIterations) { Standard_Boolean GOOD; Standard_Integer j; @@ -125,25 +148,10 @@ void math_BissecNewton::Perform(math_FunctionWithDerivative& F, return; } -Standard_Boolean math_BissecNewton::IsSolutionReached -//(math_FunctionWithDerivative& F) -(math_FunctionWithDerivative& ) -{ - return Abs(dx) <= XTol; -} - - -math_BissecNewton::math_BissecNewton(math_FunctionWithDerivative& F, - const Standard_Real Bound1, - const Standard_Real Bound2, - const Standard_Real TolX, - const Standard_Integer NbIterations) { - - XTol = TolX; - Perform(F, Bound1, Bound2, NbIterations); -} - - +//======================================================================= +//function : Dump +//purpose : +//======================================================================= void math_BissecNewton::Dump(Standard_OStream& o) const { o << "math_BissecNewton "; diff --git a/src/math/math_BissecNewton.lxx b/src/math/math_BissecNewton.lxx index aa154b10f3..06815b1cfc 100644 --- a/src/math/math_BissecNewton.lxx +++ b/src/math/math_BissecNewton.lxx @@ -14,6 +14,16 @@ #include +inline void math_BissecNewton::Delete() const +{ +} + +inline Standard_Boolean math_BissecNewton::IsSolutionReached(math_FunctionWithDerivative&) +{ + return Abs(dx) <= XTol; +} + + inline Standard_OStream& operator<<(Standard_OStream& o, const math_BissecNewton& Bi) { diff --git a/src/math/math_BrentMinimum.cdl b/src/math/math_BrentMinimum.cdl index 2501732a3a..81444a6ff8 100644 --- a/src/math/math_BrentMinimum.cdl +++ b/src/math/math_BrentMinimum.cdl @@ -32,63 +32,46 @@ is Create(TolX: Real; NbIterations: Integer = 100; ZEPS: Real = 1.0e-12) - ---Purpose: - -- This constructor should be used in a sub-class to initialize - -- correctly all the fields of this class. - + ---Purpose: + -- This constructor should be used in a sub-class to initialize + -- correctly all the fields of this class. returns BrentMinimum; Create(TolX: Real; Fbx: Real; NbIterations: Integer = 100; ZEPS: Real = 1.0e-12) - ---Purpose: - -- This constructor should be used in a sub-class to initialize - -- correctly all the fields of this class. - -- It has to be used if F(Bx) is known. - + ---Purpose: + -- This constructor should be used in a sub-class to initialize + -- correctly all the fields of this class. + -- It has to be used if F(Bx) is known. returns BrentMinimum; - - Create(F: in out Function; Ax, Bx, Cx, TolX: Real; - NbIterations: Integer = 100; ZEPS: Real=1.0e-12) - ---Purpose: - -- Given a bracketing triplet of abscissae Ax, Bx, Cx - -- (such as Bx is between Ax and Cx, F(Bx) is - -- less than both F(Bx) and F(Cx)) the Brent minimization is done - -- on the function F. - -- The tolerance required on F is given by Tolerance. - -- The solution is found when : - -- abs(Xi - Xi-1) <= TolX * abs(Xi) + ZEPS; - -- The maximum number of iterations allowed is given by NbIterations. - - returns BrentMinimum; - + Delete(me) is static; + ---Purpose: Destructor alias. + ---C++: inline ---C++: alias " Standard_EXPORT virtual ~math_BrentMinimum();" - Perform(me: in out; F: in out Function; - Ax, Bx, Cx: Real) - ---Purpose: - -- Brent minimization is performed on function F from a given - -- bracketing triplet of abscissas Ax, Bx, Cx (such that Bx is - -- between Ax and Cx, F(Bx) is less than both F(Bx) and F(Cx)) - -- Warning - -- The initialization constructors must have been called - -- before the call to the Perform method. + Perform(me: in out; F: in out Function; + Ax, Bx, Cx: Real) + ---Purpose: + -- Brent minimization is performed on function F from a given + -- bracketing triplet of abscissas Ax, Bx, Cx (such that Bx is + -- between Ax and Cx, F(Bx) is less than both F(Bx) and F(Cx)) + -- The solution is found when: abs(Xi - Xi-1) <= TolX * abs(Xi) + ZEPS; is static; - IsSolutionReached(me: in out; F: in out Function) - ---Purpose: - -- This method is called at the end of each iteration to check if the - -- solution is found. - -- It can be redefined in a sub-class to implement a specific test to - -- stop the iterations. - - returns Boolean - is virtual; + IsSolutionReached(me: in out; theFunction: in out Function) + ---Purpose: + -- This method is called at the end of each iteration to check if the + -- solution is found. + -- It can be redefined in a sub-class to implement a specific test to + -- stop the iterations. + ---C++:inline + returns Boolean is virtual; IsDone(me) @@ -137,7 +120,6 @@ is is static; - fields Done: Boolean; diff --git a/src/math/math_BrentMinimum.cxx b/src/math/math_BrentMinimum.cxx index a959769dec..406ce26f36 100644 --- a/src/math/math_BrentMinimum.cxx +++ b/src/math/math_BrentMinimum.cxx @@ -23,15 +23,68 @@ #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d) +//======================================================================= +//function : math_BrentMinimum +//purpose : Constructor +//======================================================================= +math_BrentMinimum::math_BrentMinimum(const Standard_Real theTolX, + const Standard_Integer theNbIterations, + const Standard_Real theZEPS) +: a (0.0), + b (0.0), + x (0.0), + fx (0.0), + fv (0.0), + fw (0.0), + XTol(theTolX), + EPSZ(theZEPS), + Done (Standard_False), + iter (0), + Itermax(theNbIterations), + myF (Standard_False) +{ +} + +//======================================================================= +//function : math_BrentMinimum +//purpose : Constructor +//======================================================================= +math_BrentMinimum::math_BrentMinimum(const Standard_Real theTolX, + const Standard_Real theFbx, + const Standard_Integer theNbIterations, + const Standard_Real theZEPS) +: a (0.0), + b (0.0), + x (0.0), + fx (theFbx), + fv (0.0), + fw (0.0), + XTol(theTolX), + EPSZ(theZEPS), + Done (Standard_False), + iter (0), + Itermax(theNbIterations), + myF (Standard_True) +{ +} + +//======================================================================= +//function : ~math_BrentMinimum +//purpose : Destructor +//======================================================================= math_BrentMinimum::~math_BrentMinimum() { } +//======================================================================= +//function : Perform +//purpose : +//======================================================================= void math_BrentMinimum::Perform(math_Function& F, - const Standard_Real ax, - const Standard_Real bx, - const Standard_Real cx) { - + const Standard_Real ax, + const Standard_Real bx, + const Standard_Real cx) +{ Standard_Boolean OK; Standard_Real etemp, fu, p, q, r; Standard_Real tol1, tol2, u, v, w, xm; @@ -104,70 +157,20 @@ void math_BrentMinimum::Perform(math_Function& F, return; } - -math_BrentMinimum::math_BrentMinimum(math_Function& F, - const Standard_Real Ax, - const Standard_Real Bx, - const Standard_Real Cx, - const Standard_Real TolX, - const Standard_Integer NbIterations, - const Standard_Real ZEPS) { - - XTol = TolX; - EPSZ = ZEPS; - Itermax = NbIterations; - myF = Standard_False; - Perform(F, Ax, Bx, Cx); -} - - -// Constructeur d'initialisation des champs. - -math_BrentMinimum::math_BrentMinimum(const Standard_Real TolX, - const Standard_Integer NbIterations, - const Standard_Real ZEPS) { - myF = Standard_False; - XTol = TolX; - EPSZ = ZEPS; - Itermax = NbIterations; -} - -math_BrentMinimum::math_BrentMinimum(const Standard_Real TolX, - const Standard_Real Fbx, - const Standard_Integer NbIterations, - const Standard_Real ZEPS) { - - fx = Fbx; - myF = Standard_True; - XTol = TolX; - EPSZ = ZEPS; - Itermax = NbIterations; -} - - -// Standard_Boolean math_BrentMinimum::IsSolutionReached(math_Function& F) { - Standard_Boolean math_BrentMinimum::IsSolutionReached(math_Function& ) { - -// Standard_Real xm = 0.5 * (a + b); - // modified by NIZHNY-MKK Mon Oct 3 17:45:57 2005.BEGIN -// Standard_Real tol = XTol * fabs(x) + EPSZ; -// return fabs(x - xm) <= 2.0 * tol - 0.5 * (b - a); - Standard_Real TwoTol = 2.0 *(XTol * fabs(x) + EPSZ); - return ((x <= (TwoTol + a)) && (x >= (b - TwoTol))); - // modified by NIZHNY-MKK Mon Oct 3 17:46:00 2005.END +//======================================================================= +//function : Dump +//purpose : +//======================================================================= +void math_BrentMinimum::Dump(Standard_OStream& o) const +{ + o << "math_BrentMinimum "; + if(Done) { + o << " Status = Done \n"; + o << " Location value = " << x <<"\n"; + o << " Minimum value = " << fx << "\n"; + o << " Number of iterations = " << iter <<"\n"; } - - - - void math_BrentMinimum::Dump(Standard_OStream& o) const { - o << "math_BrentMinimum "; - if(Done) { - o << " Status = Done \n"; - o << " Location value = " << x <<"\n"; - o << " Minimum value = " << fx << "\n"; - o << " Number of iterations = " << iter <<"\n"; - } - else { - o << " Status = not Done \n"; - } - } + else { + o << " Status = not Done \n"; + } +} diff --git a/src/math/math_BrentMinimum.lxx b/src/math/math_BrentMinimum.lxx index ea3d734554..948890225a 100644 --- a/src/math/math_BrentMinimum.lxx +++ b/src/math/math_BrentMinimum.lxx @@ -14,7 +14,23 @@ #include -inline Standard_Boolean math_BrentMinimum::IsDone() const { return Done; } +inline void math_BrentMinimum::Delete() const +{ +} + + +inline Standard_Boolean math_BrentMinimum::IsSolutionReached(math_Function&) +{ + const Standard_Real TwoTol = 2.0 * (XTol * fabs(x) + EPSZ); + return (x <= (TwoTol + a)) && (x >= (b - TwoTol)); +} + + +inline Standard_Boolean math_BrentMinimum::IsDone() const +{ + return Done; +} + inline Standard_OStream& operator<< (Standard_OStream& o, const math_BrentMinimum& Br) diff --git a/src/math/math_FRPR.cdl b/src/math/math_FRPR.cdl index 43972e08a2..a2704a1cec 100644 --- a/src/math/math_FRPR.cdl +++ b/src/math/math_FRPR.cdl @@ -30,53 +30,40 @@ raises DimensionError from Standard, is - Create(F: in out MultipleVarFunctionWithGradient; - StartingPoint: Vector; Tolerance: Real; - NbIterations: Integer=200; ZEPS: Real=1.0e-12) - ---Purpose: Computes FRPR minimization function F from input vector - -- StartingPoint. The solution F = Fi is found when 2.0 * - -- abs(Fi - Fi-1) <= Tolerance * (abs(Fi) + abs(Fi-1) + - -- ZEPS). The maximum number of iterations allowed is given - -- by NbIterations. - returns FRPR; - - - Create(F: in out MultipleVarFunctionWithGradient; - Tolerance: Real; - NbIterations: Integer = 200; - ZEPS: Real = 1.0e-12) - ---Purpose: Purpose - -- Initializes the computation of the minimum of F. - -- Warning - -- A call to the Perform method must be made after this - -- initialization to compute the minimum of the function. + Create(theFunction: in MultipleVarFunctionWithGradient; + theTolerance: Real; + theNbIterations: Integer = 200; + theZEPS: Real = 1.0e-12) + ---Purpose: + -- Initializes the computation of the minimum of F. + -- Warning: constructor does not perform computations. returns FRPR; + Delete(me) is static; + ---Purpose: Destructor alias. ---C++: alias " Standard_EXPORT virtual ~math_FRPR();" - Perform(me: in out; F: in out MultipleVarFunctionWithGradient; - StartingPoint: Vector) - ---Purpose: Use this method after a call to the initialization constructor - -- to compute the minimum of function F. - -- Warning - -- The initialization constructor must have been called before - -- the Perform method is called + Perform(me: in out; + theFunction: in out MultipleVarFunctionWithGradient; + theStartingPoint: Vector) + ---Purpose: + -- The solution F = Fi is found when + -- 2.0 * abs(Fi - Fi-1) <= Tolerance * (abs(Fi) + abs(Fi-1) + ZEPS). is static; - IsSolutionReached(me: in out; F: in out MultipleVarFunctionWithGradient) - ---Purpose: - -- The solution F = Fi is found when : - -- 2.0 * abs(Fi - Fi-1) <= Tolerance * (abs(Fi) + abs(Fi-1)) + ZEPS. - -- The maximum number of iterations allowed is given by NbIterations. - - returns Boolean - is virtual; - - + IsSolutionReached(me: in out; theFunction: in out MultipleVarFunctionWithGradient) + ---Purpose: + -- The solution F = Fi is found when: + -- 2.0 * abs(Fi - Fi-1) <= Tolerance * (abs(Fi) + abs(Fi-1)) + ZEPS. + -- The maximum number of iterations allowed is given by NbIterations. + ---C++:inline + returns Boolean is virtual; + + IsDone(me) - ---Purpose: Returns true if the computations are successful, otherwise returns false. + ---Purpose: Returns true if the computations are successful, otherwise returns false. ---C++: inline returns Boolean is static; @@ -158,6 +145,7 @@ is fields + Done: Boolean; TheLocation: Vector is protected; TheGradient: Vector is protected; diff --git a/src/math/math_FRPR.cxx b/src/math/math_FRPR.cxx index a3776e4ab1..333018901b 100644 --- a/src/math/math_FRPR.cxx +++ b/src/math/math_FRPR.cxx @@ -88,8 +88,9 @@ static Standard_Boolean MinimizeDirection(math_Vector& P, math_BracketMinimum Bracket(F, 0.0, 1.0); if(Bracket.IsDone()) { Bracket.Values(ax, xx, bx); - math_BrentMinimum Sol(F, ax, xx, bx, 1.0e-10, 100); - if(Sol.IsDone()) { + math_BrentMinimum Sol(1.e-10); + Sol.Perform(F, ax, xx, bx); + if (Sol.IsDone()) { Standard_Real Scale = Sol.Location(); Result = Sol.Minimum(); Dir.Multiply(Scale); @@ -100,10 +101,45 @@ static Standard_Boolean MinimizeDirection(math_Vector& P, return Standard_False; } - +//======================================================================= +//function : math_FRPR +//purpose : Constructor +//======================================================================= +math_FRPR::math_FRPR(const math_MultipleVarFunctionWithGradient& theFunction, + const Standard_Real theTolerance, + const Standard_Integer theNbIterations, + const Standard_Real theZEPS) + +: TheLocation(1, theFunction.NbVariables()), + TheGradient(1, theFunction.NbVariables()), + TheMinimum (0.0), + PreviousMinimum(0.0), + XTol (theTolerance), + EPSZ (theZEPS), + Done (Standard_False), + Iter (0), + State (0), + TheStatus (math_NotBracketed), + Itermax (theNbIterations) +{ +} + +//======================================================================= +//function : ~math_FRPR +//purpose : Destructor +//======================================================================= +math_FRPR::~math_FRPR() +{ +} + + +//======================================================================= +//function : Perform +//purpose : +//======================================================================= void math_FRPR::Perform(math_MultipleVarFunctionWithGradient& F, - const math_Vector& StartingPoint) { - + const math_Vector& StartingPoint) +{ Standard_Boolean Good; Standard_Integer n = TheLocation.Length(); Standard_Integer j, its; @@ -140,7 +176,7 @@ void math_FRPR::Perform(math_MultipleVarFunctionWithGradient& F, } if(IsSolutionReached(F)) { Done = Standard_True; - State = F.GetStateNumber(); + State = F.GetStateNumber(); TheStatus = math_OK; return; } @@ -175,64 +211,25 @@ void math_FRPR::Perform(math_MultipleVarFunctionWithGradient& F, Done = Standard_False; TheStatus = math_TooManyIterations; return; - - } - - - - Standard_Boolean math_FRPR::IsSolutionReached( -// math_MultipleVarFunctionWithGradient& F) { - math_MultipleVarFunctionWithGradient& ) { - - return (2.0 * fabs(TheMinimum - PreviousMinimum)) <= - XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ); - } - - math_FRPR::math_FRPR(math_MultipleVarFunctionWithGradient& F, - const math_Vector& StartingPoint, - const Standard_Real Tolerance, - const Standard_Integer NbIterations, - const Standard_Real ZEPS) - : TheLocation(1, StartingPoint.Length()), - TheGradient(1, StartingPoint.Length()) { - - XTol = Tolerance; - EPSZ = ZEPS; - Itermax = NbIterations; - Perform(F, StartingPoint); - } - - - math_FRPR::math_FRPR(math_MultipleVarFunctionWithGradient& F, - const Standard_Real Tolerance, - const Standard_Integer NbIterations, - const Standard_Real ZEPS) - : TheLocation(1, F.NbVariables()), - TheGradient(1, F.NbVariables()) { - - XTol = Tolerance; - EPSZ = ZEPS; - Itermax = NbIterations; - } - - - math_FRPR::~math_FRPR() - { - } - - void math_FRPR::Dump(Standard_OStream& o) const { - - o << "math_FRPR "; - if(Done) { - o << " Status = Done \n"; - o << " Location Vector = "<< TheLocation << "\n"; - o << " Minimum value = " << TheMinimum <<"\n"; - o << " Number of iterations = " << Iter <<"\n"; - } - else { - o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n"; - } - } +} + +//======================================================================= +//function : Dump +//purpose : +//======================================================================= +void math_FRPR::Dump(Standard_OStream& o) const +{ + o << "math_FRPR "; + if(Done) { + o << " Status = Done \n"; + o << " Location Vector = "<< TheLocation << "\n"; + o << " Minimum value = " << TheMinimum <<"\n"; + o << " Number of iterations = " << Iter <<"\n"; + } + else { + o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n"; + } +} diff --git a/src/math/math_FRPR.lxx b/src/math/math_FRPR.lxx index 636c0bae89..aa0a9e670f 100644 --- a/src/math/math_FRPR.lxx +++ b/src/math/math_FRPR.lxx @@ -15,8 +15,21 @@ #include #include +inline void math_FRPR::Delete() const +{ +} + +inline Standard_Boolean math_FRPR::IsSolutionReached(math_MultipleVarFunctionWithGradient&) +{ + return 2.0 * fabs(TheMinimum - PreviousMinimum) <= + XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ); +} + +inline Standard_Boolean math_FRPR::IsDone() const +{ + return Done; -inline Standard_Boolean math_FRPR::IsDone() const { return Done; } +} inline Standard_OStream& operator<<(Standard_OStream& o, const math_FRPR& Fr) diff --git a/src/math/math_FunctionRoot.cxx b/src/math/math_FunctionRoot.cxx index 0cd6182ab0..b35e156e6a 100644 --- a/src/math/math_FunctionRoot.cxx +++ b/src/math/math_FunctionRoot.cxx @@ -76,7 +76,8 @@ class math_MyFunctionSetWithDerivatives : public math_FunctionSetWithDerivatives math_MyFunctionSetWithDerivatives Ff(F); V(1)=Guess; Tol(1) = Tolerance; - math_FunctionSetRoot Sol(Ff,V,Tol,NbIterations); + math_FunctionSetRoot Sol(Ff, Tol, NbIterations); + Sol.Perform(Ff, V); Done = Sol.IsDone(); if (Done) { F.GetStateNumber(); @@ -98,7 +99,8 @@ class math_MyFunctionSetWithDerivatives : public math_FunctionSetWithDerivatives Tol(1) = Tolerance; Aa(1)=A; Bb(1)=B; - math_FunctionSetRoot Sol(Ff,V,Tol,Aa,Bb,NbIterations); + math_FunctionSetRoot Sol(Ff, Tol, NbIterations); + Sol.Perform(Ff, V, Aa, Bb); Done = Sol.IsDone(); if (Done) { F.GetStateNumber(); diff --git a/src/math/math_FunctionSetRoot.cdl b/src/math/math_FunctionSetRoot.cdl index 40213802fe..fe09f7405f 100644 --- a/src/math/math_FunctionSetRoot.cdl +++ b/src/math/math_FunctionSetRoot.cdl @@ -57,64 +57,50 @@ is -- constructor. returns FunctionSetRoot from math; - - - - Create(F: in out FunctionSetWithDerivatives; StartingPoint: Vector; - Tolerance: Vector; NbIterations: Integer = 100) - ---Purpose: is used to improve the root of the function F - -- from the initial guess StartingPoint. - -- The maximum number of iterations allowed is given by - -- NbIterations. - -- In this case, the solution is found when: - -- abs(Xi - Xi-1)(j) <= Tolerance(j) for all unknowns. - - returns FunctionSetRoot from math; - - - Create(F: in out FunctionSetWithDerivatives; StartingPoint: Vector; - Tolerance: Vector; infBound, supBound: Vector; - NbIterations: Integer = 100; theStopOnDivergent : Boolean from Standard = Standard_False) - ---Purpose: is used to improve the root of the function F - -- from the initial guess StartingPoint. - -- The maximum number of iterations allowed is given - -- by NbIterations. - -- In this case, the solution is found when: - -- abs(Xi - Xi-1) <= Tolerance for all unknowns. - - returns FunctionSetRoot from math; + Delete(me) is static; + ---Purpose: Destructor alias. ---C++: alias " Standard_EXPORT virtual ~math_FunctionSetRoot();" - SetTolerance(me: in out; Tolerance: Vector) - ---Purpose: Initializes the tolerance values. + SetTolerance(me: in out; Tolerance: Vector) + ---Purpose: Initializes the tolerance values. is static; - - - Perform(me: in out; F: in out FunctionSetWithDerivatives; - StartingPoint: Vector; - infBound, supBound: Vector; theStopOnDivergent : Boolean from Standard = Standard_False) - ---Purpose: Improves the root of function F from the initial guess - -- StartingPoint. infBound and supBound may be given to constrain the solution. - -- Warning - -- This method is called when computation of the solution is - -- not performed by the constructors. + IsSolutionReached(me: in out; F: in out FunctionSetWithDerivatives) + ---Purpose: This routine is called at the end of each iteration + -- to check if the solution was found. It can be redefined + -- in a sub-class to implement a specific test to stop the iterations. + -- In this case, the solution is found when: abs(Xi - Xi-1) <= Tolerance + -- for all unknowns. + ---C++: inline + returns Boolean is virtual; - is static; + Perform(me: in out; + theFunction: in out FunctionSetWithDerivatives; + theStartingPoint: Vector; + theStopOnDivergent: Boolean from Standard = Standard_False) + ---Purpose: + -- Improves the root of function from the initial guess point. + -- The infinum and supremum may be given to constrain the solution. + -- In this case, the solution is found when: abs(Xi - Xi-1)(j) <= Tolerance(j) + -- for all unknowns. + is static; - IsSolutionReached(me: in out; F: in out FunctionSetWithDerivatives) - ---Purpose: This routine is called at the end of each iteration - -- to check if the solution was found. It can be redefined - -- in a sub-class to implement a specific test to stop the - -- iterations. - -- In this case, the solution is found when: - -- abs(Xi - Xi-1) <= Tolerance for all unknowns. - returns Boolean is virtual; + Perform(me: in out; + theFunction: in out FunctionSetWithDerivatives; + theStartingPoint: Vector; + theInfBound, theSupBound: Vector; + theStopOnDivergent: Boolean from Standard = Standard_False) + ---Purpose: + -- Improves the root of function from the initial guess point. + -- The infinum and supremum may be given to constrain the solution. + -- In this case, the solution is found when: abs(Xi - Xi-1) <= Tolerance + -- for all unknowns. + is static; IsDone(me) diff --git a/src/math/math_FunctionSetRoot.cxx b/src/math/math_FunctionSetRoot.cxx index 3facf61631..a57b5177b2 100644 --- a/src/math/math_FunctionSetRoot.cxx +++ b/src/math/math_FunctionSetRoot.cxx @@ -228,8 +228,12 @@ static Standard_Boolean MinimizeDirection(const math_Vector& P0, FSR_DEBUG (" minimisation dans la direction") ax = -1; bx = 0; cx = (P2-P1).Norm()*invnorme; - if (cx < 1.e-2) return Standard_False; - math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d); + if (cx < 1.e-2) + return Standard_False; + + math_BrentMinimum Sol(tol1d, 100, tol1d); + Sol.Perform(F, ax, bx, cx); + if(Sol.IsDone()) { tsol = Sol.Location(); if (Sol.Minimum() < F1) { @@ -322,7 +326,10 @@ static Standard_Boolean MinimizeDirection(const math_Vector& P, ax = 0.0; bx = tsol; cx = 1.0; } FSR_DEBUG(" minimisation dans la direction"); - math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d); + + math_BrentMinimum Sol(tol1d, 100, tol1d); + Sol.Perform(F, ax, bx, cx); + if(Sol.IsDone()) { if (Sol.Minimum() <= Result) { tsol = Sol.Location(); @@ -551,155 +558,109 @@ Standard_Boolean Bounds(const math_Vector& InfBound, - -math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F, - const math_Vector& Tolerance, - const Standard_Integer NbIterations) : -Delta(1, F.NbVariables()), -Sol(1, F.NbVariables()), -DF(1, F.NbEquations(), - 1, F.NbVariables()), - Tol(1,F.NbVariables()), - - InfBound(1, F.NbVariables()), - SupBound(1, F.NbVariables()), - - SolSave(1, F.NbVariables()), - GH(1, F.NbVariables()), - DH(1, F.NbVariables()), - DHSave(1, F.NbVariables()), - FF(1, F.NbEquations()), - PreviousSolution(1, F.NbVariables()), - Save(0, NbIterations), - Constraints(1, F.NbVariables()), - Temp1(1, F.NbVariables()), - Temp2(1, F.NbVariables()), - Temp3(1, F.NbVariables()), - Temp4(1, F.NbEquations()) - -{ - for (Standard_Integer i = 1; i <= Tol.Length(); i++) { - Tol(i) =Tolerance(i); - } - Itermax = NbIterations; -} - -math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F, - const Standard_Integer NbIterations) : -Delta(1, F.NbVariables()), -Sol(1, F.NbVariables()), -DF(1, F.NbEquations(), - 1, F.NbVariables()), - Tol(1, F.NbVariables()), - - InfBound(1, F.NbVariables()), - SupBound(1, F.NbVariables()), - - SolSave(1, F.NbVariables()), - GH(1, F.NbVariables()), - DH(1, F.NbVariables()), - DHSave(1, F.NbVariables()), - FF(1, F.NbEquations()), - PreviousSolution(1, F.NbVariables()), - Save(0, NbIterations), - Constraints(1, F.NbVariables()), - Temp1(1, F.NbVariables()), - Temp2(1, F.NbVariables()), - Temp3(1, F.NbVariables()), - Temp4(1, F.NbEquations()) - +//======================================================================= +//function : math_FunctionSetRoot +//purpose : Constructor +//======================================================================= +math_FunctionSetRoot::math_FunctionSetRoot( + math_FunctionSetWithDerivatives& theFunction, + const math_Vector& theTolerance, + const Standard_Integer theNbIterations) + +: Delta(1, theFunction.NbVariables()), + Sol (1, theFunction.NbVariables()), + DF (1, theFunction.NbEquations() , 1, theFunction.NbVariables()), + Tol (1, theFunction.NbVariables()), + Done (Standard_False), + Kount (0), + State (0), + Itermax (theNbIterations), + InfBound(1, theFunction.NbVariables(), RealFirst()), + SupBound(1, theFunction.NbVariables(), RealLast ()), + SolSave (1, theFunction.NbVariables()), + GH (1, theFunction.NbVariables()), + DH (1, theFunction.NbVariables()), + DHSave (1, theFunction.NbVariables()), + FF (1, theFunction.NbEquations()), + PreviousSolution(1, theFunction.NbVariables()), + Save (0, theNbIterations), + Constraints(1, theFunction.NbVariables()), + Temp1 (1, theFunction.NbVariables()), + Temp2 (1, theFunction.NbVariables()), + Temp3 (1, theFunction.NbVariables()), + Temp4 (1, theFunction.NbEquations()), + myIsDivergent(Standard_False) { - Itermax = NbIterations; + SetTolerance(theTolerance); } - - -math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F, - const math_Vector& StartingPoint, - const math_Vector& Tolerance, - const math_Vector& infBound, - const math_Vector& supBound, - const Standard_Integer NbIterations, - Standard_Boolean theStopOnDivergent) : -Delta(1, F.NbVariables()), -Sol(1, F.NbVariables()), -DF(1, F.NbEquations(), - 1, F.NbVariables()), - Tol(1,F.NbVariables()), - - - InfBound(1, F.NbVariables()), - SupBound(1, F.NbVariables()), - - SolSave(1, F.NbVariables()), - GH(1, F.NbVariables()), - DH(1, F.NbVariables()), - DHSave(1, F.NbVariables()), - FF(1, F.NbEquations()), - PreviousSolution(1, F.NbVariables()), - Save(0, NbIterations), - Constraints(1, F.NbVariables()), - Temp1(1, F.NbVariables()), - Temp2(1, F.NbVariables()), - Temp3(1, F.NbVariables()), - Temp4(1, F.NbEquations()) - +//======================================================================= +//function : math_FunctionSetRoot +//purpose : Constructor +//======================================================================= +math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& theFunction, + const Standard_Integer theNbIterations) + +: Delta(1, theFunction.NbVariables()), + Sol (1, theFunction.NbVariables()), + DF (1, theFunction.NbEquations() , 1, theFunction.NbVariables()), + Tol (1, theFunction.NbVariables()), + Done (Standard_False), + Kount (0), + State (0), + Itermax (theNbIterations), + InfBound(1, theFunction.NbVariables(), RealFirst()), + SupBound(1, theFunction.NbVariables(), RealLast ()), + SolSave (1, theFunction.NbVariables()), + GH (1, theFunction.NbVariables()), + DH (1, theFunction.NbVariables()), + DHSave (1, theFunction.NbVariables()), + FF (1, theFunction.NbEquations()), + PreviousSolution(1, theFunction.NbVariables()), + Save (0, theNbIterations), + Constraints(1, theFunction.NbVariables()), + Temp1 (1, theFunction.NbVariables()), + Temp2 (1, theFunction.NbVariables()), + Temp3 (1, theFunction.NbVariables()), + Temp4 (1, theFunction.NbEquations()), + myIsDivergent(Standard_False) { - - for (Standard_Integer i = 1; i <= Tol.Length(); i++) { - Tol(i) =Tolerance(i); - } - Itermax = NbIterations; - Perform(F, StartingPoint, infBound, supBound, theStopOnDivergent); } -math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F, - const math_Vector& StartingPoint, - const math_Vector& Tolerance, - const Standard_Integer NbIterations) : -Delta(1, F.NbVariables()), -Sol(1, F.NbVariables()), -DF(1, F.NbEquations(), - 1, StartingPoint.Length()), - Tol(1,F.NbVariables()), - - InfBound(1, F.NbVariables()), - SupBound(1, F.NbVariables()), - - SolSave(1, F.NbVariables()), - GH(1, F.NbVariables()), - DH(1, F.NbVariables()), - DHSave(1, F.NbVariables()), - FF(1, F.NbEquations()), - PreviousSolution(1, F.NbVariables()), - Save(0, NbIterations), - Constraints(1, F.NbVariables()), - Temp1(1, F.NbVariables()), - Temp2(1, F.NbVariables()), - Temp3(1, F.NbVariables()), - Temp4(1, F.NbEquations()) - +//======================================================================= +//function : ~math_FunctionSetRoot +//purpose : Destructor +//======================================================================= +math_FunctionSetRoot::~math_FunctionSetRoot() { - for (Standard_Integer i = 1; i <= Tol.Length(); i++) { - Tol(i) = Tolerance(i); - } - Itermax = NbIterations; - InfBound.Init(RealFirst()); - SupBound.Init(RealLast()); - Perform(F, StartingPoint, InfBound, SupBound); + Delete(); } -math_FunctionSetRoot::~math_FunctionSetRoot() +//======================================================================= +//function : SetTolerance +//purpose : +//======================================================================= +void math_FunctionSetRoot::SetTolerance(const math_Vector& theTolerance) { + for (Standard_Integer i = 1; i <= Tol.Length(); ++i) + Tol(i) = theTolerance(i); } -void math_FunctionSetRoot::SetTolerance(const math_Vector& Tolerance) +//======================================================================= +//function : Perform +//purpose : +//======================================================================= +void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& theFunction, + const math_Vector& theStartingPoint, + const Standard_Boolean theStopOnDivergent) { - for (Standard_Integer i = 1; i <= Tol.Length(); i++) { - Tol(i) = Tolerance(i); - } + Perform(theFunction, theStartingPoint, InfBound, SupBound, theStopOnDivergent); } +//======================================================================= +//function : Perform +//purpose : +//======================================================================= void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F, const math_Vector& StartingPoint, const math_Vector& theInfBound, @@ -1232,38 +1193,40 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F, } } - - - -Standard_Boolean math_FunctionSetRoot::IsSolutionReached(math_FunctionSetWithDerivatives& ) { - for(Standard_Integer i = 1; i<= Sol.Length(); i++) { - if(Abs(Delta(i)) > Tol(i)) {return Standard_False;} - } - return Standard_True; -} - - -void math_FunctionSetRoot::Dump(Standard_OStream& o) const { - o<<" math_FunctionSetRoot"; +//======================================================================= +//function : Dump +//purpose : +//======================================================================= +void math_FunctionSetRoot::Dump(Standard_OStream& o) const +{ + o << " math_FunctionSetRoot"; if (Done) { o << " Status = Done\n"; o << " Location value = " << Sol << "\n"; o << " Number of iterations = " << Kount << "\n"; } else { - o<<"Status = Not Done\n"; + o << "Status = Not Done\n"; } } - -void math_FunctionSetRoot::Root(math_Vector& Root) const{ +//======================================================================= +//function : Root +//purpose : +//======================================================================= +void math_FunctionSetRoot::Root(math_Vector& Root) const +{ StdFail_NotDone_Raise_if(!Done, " "); Standard_DimensionError_Raise_if(Root.Length() != Sol.Length(), " "); Root = Sol; } - -void math_FunctionSetRoot::FunctionSetErrors(math_Vector& Err) const{ +//======================================================================= +//function : FunctionSetErrors +//purpose : +//======================================================================= +void math_FunctionSetRoot::FunctionSetErrors(math_Vector& Err) const +{ StdFail_NotDone_Raise_if(!Done, " "); Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " "); Err = Delta; diff --git a/src/math/math_FunctionSetRoot.lxx b/src/math/math_FunctionSetRoot.lxx index adc615c0f7..50937dba46 100644 --- a/src/math/math_FunctionSetRoot.lxx +++ b/src/math/math_FunctionSetRoot.lxx @@ -15,8 +15,26 @@ #include #include +inline void math_FunctionSetRoot::Delete() const +{ +} + + +inline Standard_Boolean math_FunctionSetRoot::IsSolutionReached(math_FunctionSetWithDerivatives&) +{ + for (Standard_Integer i = 1; i <= Sol.Length(); ++i) + if ( Abs(Delta(i)) > Tol(i) ) + return Standard_False; + + return Standard_True; +} + + +inline Standard_Boolean math_FunctionSetRoot::IsDone() const +{ + return Done; +} -inline Standard_Boolean math_FunctionSetRoot::IsDone() const { return Done; } inline Standard_OStream& operator<<(Standard_OStream& o, const math_FunctionSetRoot& F) diff --git a/src/math/math_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx index 0efeedb3ec..ce474511f3 100644 --- a/src/math/math_GlobOptMin.cxx +++ b/src/math/math_GlobOptMin.cxx @@ -247,7 +247,9 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt math_MultipleVarFunctionWithHessian* myTmp = dynamic_cast (myFunc); - math_NewtonMinimum newtonMinimum(*myTmp, thePnt); + math_NewtonMinimum newtonMinimum(*myTmp); + newtonMinimum.Perform(*myTmp, thePnt); + if (newtonMinimum.IsDone()) { newtonMinimum.Location(theOutPnt); @@ -278,7 +280,8 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt for(i = 1; i <= myN; i++) m(1, 1) = 1.0; - math_Powell powell(*myFunc, thePnt, m, 1e-10); + math_Powell powell(*myFunc, 1e-10); + powell.Perform(*myFunc, thePnt, m); if (powell.IsDone()) { diff --git a/src/math/math_NewtonFunctionSetRoot.cdl b/src/math/math_NewtonFunctionSetRoot.cdl index 69a9f27f93..4b8b2e0aa9 100644 --- a/src/math/math_NewtonFunctionSetRoot.cdl +++ b/src/math/math_NewtonFunctionSetRoot.cdl @@ -30,96 +30,66 @@ raises NotDone from StdFail, is - Create(F: in out FunctionSetWithDerivatives; - XTol: Vector; FTol: Real; - NbIterations: Integer = 100) - ---Purpose: - -- This constructor should be used in a sub-class to initialize - -- correctly all the fields of this class. - -- The range (1, F.NbVariables()) must be especially respected for - -- all vectors and matrix declarations. - - returns NewtonFunctionSetRoot; - - - Create(F: in out FunctionSetWithDerivatives; - FTol: Real; NbIterations: Integer = 100) - ---Purpose: - -- This constructor should be used in a sub-class to initialize - -- correctly all the fields of this class. - -- The range (1, F.NbVariables()) must be especially respected for - -- all vectors and matrix declarations. - -- The method SetTolerance must be called before performing the - -- algorithm. - - returns NewtonFunctionSetRoot; - - - - Create(F: in out FunctionSetWithDerivatives; - StartingPoint: Vector; - XTol: Vector; FTol: Real; - NbIterations: Integer = 100) - ---Purpose: - -- The Newton method is done to improve the root of the function F - -- from the initial guess StartingPoint. - -- The tolerance required on the root is given by Tolerance. - -- The solution is found when : - -- abs(Xj - Xj-1)(i) <= XTol(i) and abs(Fi) <= FTol for all i; - -- The maximum number of iterations allowed is given by NbIterations. - - returns NewtonFunctionSetRoot; - - - - Create(F: in out FunctionSetWithDerivatives; - StartingPoint: Vector; - InfBound, SupBound: Vector; - XTol: Vector; FTol: Real; - NbIterations: Integer = 100) - ---Purpose: - -- The Newton method is done to improve the root of the function F - -- from the initial guess StartingPoint. - -- The tolerance required on the root is given by Tolerance. - -- The solution is found when : - -- abs(Xj - Xj-1)(i) <= XTol(i) and abs(Fi) <= FTol for all i; - -- The maximum number of iterations allowed is given by NbIterations. - - returns NewtonFunctionSetRoot; + Create(theFunction: in out FunctionSetWithDerivatives; + theXTolerance: Vector; theFTolerance: Real; + tehNbIterations: Integer = 100) + ---Purpose: + -- Initialize correctly all the fields of this class. + -- The range (1, F.NbVariables()) must be especially respected for + -- all vectors and matrix declarations. + returns NewtonFunctionSetRoot; + + + Create(theFunction: in out FunctionSetWithDerivatives; + theFTolerance: Real; theNbIterations: Integer = 100) + ---Purpose: + -- This constructor should be used in a sub-class to initialize + -- correctly all the fields of this class. + -- The range (1, F.NbVariables()) must be especially respected for + -- all vectors and matrix declarations. + -- The method SetTolerance must be called before performing the algorithm. + returns NewtonFunctionSetRoot; + + + Delete(me) is static; + ---Purpose: Destructor alias. + ---C++: alias " Standard_EXPORT virtual ~math_NewtonFunctionSetRoot();" - ---C++: alias " Standard_EXPORT virtual ~math_NewtonFunctionSetRoot();" - SetTolerance(me: in out; XTol: Vector) - ---Purpose: Initializes the tolerance values for the unknowns. + ---Purpose: Initializes the tolerance values for the unknowns. + is static; + + Perform(me: in out; theFunction: in out FunctionSetWithDerivatives; theStartingPoint: Vector) + ---Purpose: + -- The Newton method is done to improve the root of the function + -- from the initial guess point. The solution is found when: + -- abs(Xj - Xj-1)(i) <= XTol(i) and abs(Fi) <= FTol for all i; is static; - - Perform(me: in out; F:in out FunctionSetWithDerivatives; - StartingPoint: Vector; - InfBound, SupBound: Vector) - ---Purpose: Improves the root of function F from the initial guess - -- StartingPoint. infBound and supBound may be given, to constrain the solution. - -- Warning - -- This method must be called when the solution is not computed by the constructors. + Perform(me: in out; theFunction: in out FunctionSetWithDerivatives; + theStartingPoint: Vector; theInfBound, theSupBound: Vector) + ---Purpose: + -- The Newton method is done to improve the root of the function + -- from the initial guess point. Bounds may be given, to constrain the solution. + -- The solution is found when: + -- abs(Xj - Xj-1)(i) <= XTol(i) and abs(Fi) <= FTol for all i; is static; - IsSolutionReached(me: in out; F: in out FunctionSetWithDerivatives) - ---Purpose: - -- This method is called at the end of each iteration to check if the - -- solution is found. - -- Vectors DeltaX, Fvalues and Jacobian Matrix are consistent with the - -- possible solution Vector Sol and can be inspected to decide whether - -- the solution is reached or not. + ---Purpose: + -- This method is called at the end of each iteration to check if the + -- solution is found. + -- Vectors DeltaX, Fvalues and Jacobian Matrix are consistent with the + -- possible solution Vector Sol and can be inspected to decide whether + -- the solution is reached or not. + ---C++: inline + returns Boolean is virtual; + - returns Boolean - is virtual; - - IsDone(me) ---Purpose: Returns true if the computations are successful, otherwise returns false. ---C++: inline diff --git a/src/math/math_NewtonFunctionSetRoot.cxx b/src/math/math_NewtonFunctionSetRoot.cxx index 57fffeea48..78c5524d43 100644 --- a/src/math/math_NewtonFunctionSetRoot.cxx +++ b/src/math/math_NewtonFunctionSetRoot.cxx @@ -22,128 +22,93 @@ #include #include -Standard_Boolean math_NewtonFunctionSetRoot::IsSolutionReached -// (math_FunctionSetWithDerivatives& F) - (math_FunctionSetWithDerivatives& ) -{ - - for(Standard_Integer i = DeltaX.Lower(); i <= DeltaX.Upper(); i++) { - if(Abs(DeltaX(i)) > TolX(i) || Abs(FValues(i)) > TolF) return Standard_False; - } - return Standard_True; -} - - -// Constructeurs d'initialisation des champs (pour utiliser Perform) +//======================================================================= +//function : math_NewtonFunctionSetRoot +//purpose : Constructor +//======================================================================= math_NewtonFunctionSetRoot::math_NewtonFunctionSetRoot( - math_FunctionSetWithDerivatives& F, - const math_Vector& XTol, - const Standard_Real FTol, - const Standard_Integer NbIterations): - TolX(1, F.NbVariables()), - TolF(FTol), - Indx(1, F.NbVariables()), - Scratch(1, F.NbVariables()), - Sol(1, F.NbVariables()), - DeltaX(1, F.NbVariables()), - FValues(1, F.NbVariables()), - Jacobian(1, F.NbVariables(), - 1, F.NbVariables()), - Itermax(NbIterations) + math_FunctionSetWithDerivatives& theFunction, + const math_Vector& theXTolerance, + const Standard_Real theFTolerance, + const Standard_Integer theNbIterations) + +: TolX (1, theFunction.NbVariables()), + TolF (theFTolerance), + Indx (1, theFunction.NbVariables()), + Scratch (1, theFunction.NbVariables()), + Sol (1, theFunction.NbVariables()), + DeltaX (1, theFunction.NbVariables()), + FValues (1, theFunction.NbVariables()), + Jacobian(1, theFunction.NbVariables(), 1, theFunction.NbVariables()), + Done (Standard_False), + State (0), + Iter (0), + Itermax (theNbIterations) { - for (Standard_Integer i = 1; i <= TolX.Length(); i++) { - TolX(i) = XTol(i); - } + SetTolerance(theXTolerance); } - +//======================================================================= +//function : math_NewtonFunctionSetRoot +//purpose : Constructor +//======================================================================= math_NewtonFunctionSetRoot::math_NewtonFunctionSetRoot( - math_FunctionSetWithDerivatives& F, - const Standard_Real FTol, - const Standard_Integer NbIterations): - TolX(1, F.NbVariables()), - TolF(FTol), - Indx(1, F.NbVariables()), - Scratch(1, F.NbVariables()), - Sol(1, F.NbVariables()), - DeltaX(1, F.NbVariables()), - FValues(1, F.NbVariables()), - Jacobian(1, F.NbVariables(), - 1, F.NbVariables()), - Itermax(NbIterations) + math_FunctionSetWithDerivatives& theFunction, + const Standard_Real theFTolerance, + const Standard_Integer theNbIterations) + +: TolX (1, theFunction.NbVariables()), + TolF (theFTolerance), + Indx (1, theFunction.NbVariables()), + Scratch (1, theFunction.NbVariables()), + Sol (1, theFunction.NbVariables()), + DeltaX (1, theFunction.NbVariables()), + FValues (1, theFunction.NbVariables()), + Jacobian(1, theFunction.NbVariables(), 1, theFunction.NbVariables()), + Done (Standard_False), + State (0), + Iter (0), + Itermax (theNbIterations) { } - -math_NewtonFunctionSetRoot::math_NewtonFunctionSetRoot - (math_FunctionSetWithDerivatives& F, - const math_Vector& StartingPoint, - const math_Vector& XTol, - const Standard_Real FTol, - const Standard_Integer NbIterations) : - - TolX(1, F.NbVariables()), - TolF(FTol), - Indx (1, F.NbVariables()), - Scratch (1, F.NbVariables()), - Sol (1, F.NbVariables()), - DeltaX (1, F.NbVariables()), - FValues (1, F.NbVariables()), - Jacobian(1, F.NbVariables(), - 1, F.NbVariables()), - Itermax(NbIterations) +//======================================================================= +//function : ~math_NewtonFunctionSetRoot +//purpose : Destructor +//======================================================================= +math_NewtonFunctionSetRoot::~math_NewtonFunctionSetRoot() { - for (Standard_Integer i = 1; i <= TolX.Length(); i++) { - TolX(i) = XTol(i); - } - math_Vector UFirst(1, F.NbVariables()), - ULast(1, F.NbVariables()); - UFirst.Init(RealFirst()); - ULast.Init(RealLast()); - Perform(F, StartingPoint, UFirst, ULast); } -math_NewtonFunctionSetRoot::math_NewtonFunctionSetRoot - (math_FunctionSetWithDerivatives& F, - const math_Vector& StartingPoint, - const math_Vector& InfBound, - const math_Vector& SupBound, - const math_Vector& XTol, - const Standard_Real FTol, - const Standard_Integer NbIterations) : - - TolX(1, F.NbVariables()), - TolF(FTol), - Indx (1, F.NbVariables()), - Scratch (1, F.NbVariables()), - Sol (1, F.NbVariables()), - DeltaX (1, F.NbVariables()), - FValues (1, F.NbVariables()), - Jacobian(1, F.NbVariables(), - 1, F.NbVariables()), - Itermax(NbIterations) +//======================================================================= +//function : SetTolerance +//purpose : +//======================================================================= +void math_NewtonFunctionSetRoot::SetTolerance(const math_Vector& theXTolerance) { - for (Standard_Integer i = 1; i <= TolX.Length(); i++) { - TolX(i) = XTol(i); - } - Perform(F, StartingPoint, InfBound, SupBound); + for (Standard_Integer i = 1; i <= TolX.Length(); ++i) + TolX(i) = theXTolerance(i); } -math_NewtonFunctionSetRoot::~math_NewtonFunctionSetRoot() +//======================================================================= +//function : Perform +//purpose : +//======================================================================= +void math_NewtonFunctionSetRoot::Perform( + math_FunctionSetWithDerivatives& theFunction, + const math_Vector& theStartingPoint) { -} + const math_Vector anInf(1, theFunction.NbVariables(), RealFirst()); + const math_Vector aSup (1, theFunction.NbVariables(), RealLast ()); -void math_NewtonFunctionSetRoot::SetTolerance - (const math_Vector& XTol) -{ - for (Standard_Integer i = 1; i <= TolX.Length(); i++) { - TolX(i) = XTol(i); - } + Perform(theFunction, theStartingPoint, anInf, aSup); } - - +//======================================================================= +//function : Perform +//purpose : +//======================================================================= void math_NewtonFunctionSetRoot::Perform( math_FunctionSetWithDerivatives& F, const math_Vector& StartingPoint, @@ -184,6 +149,10 @@ void math_NewtonFunctionSetRoot::Perform( } } +//======================================================================= +//function : Dump +//purpose : +//======================================================================= void math_NewtonFunctionSetRoot::Dump(Standard_OStream& o) const { o <<"math_NewtonFunctionSetRoot "; diff --git a/src/math/math_NewtonFunctionSetRoot.lxx b/src/math/math_NewtonFunctionSetRoot.lxx index 39b66a2bf7..c318555dd5 100644 --- a/src/math/math_NewtonFunctionSetRoot.lxx +++ b/src/math/math_NewtonFunctionSetRoot.lxx @@ -14,8 +14,23 @@ #include -inline Standard_Boolean math_NewtonFunctionSetRoot::IsDone() const - { return Done;} +inline void math_NewtonFunctionSetRoot::Delete() const +{ +} + +inline Standard_Boolean math_NewtonFunctionSetRoot::IsSolutionReached(math_FunctionSetWithDerivatives&) +{ + for (Standard_Integer i = DeltaX.Lower(); i <= DeltaX.Upper(); ++i) + if ( Abs(DeltaX(i)) > TolX(i) || Abs(FValues(i)) > TolF ) + return Standard_False; + + return Standard_True; +} + +inline Standard_Boolean math_NewtonFunctionSetRoot::IsDone() const +{ + return Done; +} inline Standard_OStream& operator<<(Standard_OStream& o, const math_NewtonFunctionSetRoot& N) diff --git a/src/math/math_NewtonMinimum.cdl b/src/math/math_NewtonMinimum.cdl index dc8bac1651..d4a0b886ae 100644 --- a/src/math/math_NewtonMinimum.cdl +++ b/src/math/math_NewtonMinimum.cdl @@ -28,57 +28,41 @@ raises NotDone, DimensionError is - Create(F: in out MultipleVarFunctionWithHessian; - StartingPoint: Vector; - Tolerance: Real=1.0e-7; - NbIterations: Integer=40; - Convexity: Real=1.0e-6; - WithSingularity : Boolean = Standard_True) - ---Purpose: -- Given the starting point StartingPoint, - -- The tolerance required on the solution is given by - -- Tolerance. - -- Iteration are stopped if - -- (!WithSingularity) and H(F(Xi)) is not definite - -- positive (if the smaller eigenvalue of H < Convexity) - -- or IsConverged() returns True for 2 successives Iterations. - -- Warning: Obsolete Constructor (because IsConverged can not be redefined - -- with this. ) - returns NewtonMinimum; - - Create(F: in out MultipleVarFunctionWithHessian; - Tolerance: Real=1.0e-7; - NbIterations: Integer=40; - Convexity: Real=1.0e-6; - WithSingularity : Boolean = Standard_True) + Create(theFunction: in MultipleVarFunctionWithHessian; + theTolerance: Real=1.0e-7; + theNbIterations: Integer=40; + theConvexity: Real=1.0e-6; + theWithSingularity: Boolean = Standard_True) ---Purpose: - -- The tolerance required on the solution is given by - -- Tolerance. - -- Iteration are stopped if - -- (!WithSingularity) and H(F(Xi)) is not definite - -- positive (if the smaller eigenvalue of H < Convexity) - -- or IsConverged() returns True for 2 successives Iterations. - -- Warning: This constructor do not computation + -- The tolerance required on the solution is given by Tolerance. + -- Iteration are stopped if (!WithSingularity) and H(F(Xi)) is not definite + -- positive (if the smaller eigenvalue of H < Convexity) + -- or IsConverged() returns True for 2 successives Iterations. + -- Warning: This constructor does not perform computation. returns NewtonMinimum; - - ---C++: alias " Standard_EXPORT virtual ~math_NewtonMinimum();" - - Perform(me: in out; F: in out MultipleVarFunctionWithHessian; - StartingPoint: Vector) - ---Purpose: Search the solution. + + + Perform(me: in out; theFunction: in out MultipleVarFunctionWithHessian; + theStartingPoint: Vector) + ---Purpose: Search the solution. is static; - + + Delete(me) is static; + ---Purpose: Destructor alias. + ---C++: inline + ---C++: alias " Standard_EXPORT virtual ~math_NewtonMinimum();" + + IsConverged(me) - ---Purpose: This method is called at the end of each - -- iteration to check the convergence : - -- || Xi+1 - Xi || < Tolerance - -- or || F(Xi+1) - F(Xi)|| < Tolerance * || F(Xi) || - -- It can be redefined in a sub-class to implement a specific test. - - returns Boolean - is virtual; - - + ---Purpose: + -- This method is called at the end of each iteration to check the convergence: + -- || Xi+1 - Xi || < Tolerance or || F(Xi+1) - F(Xi)|| < Tolerance * || F(Xi) || + -- It can be redefined in a sub-class to implement a specific test. + ---C++: inline + returns Boolean is virtual; + + IsDone(me) ---Purpose: Tests if an error has occured. ---C++: inline diff --git a/src/math/math_NewtonMinimum.cxx b/src/math/math_NewtonMinimum.cxx index 11f64f6fe1..7d3ebc4ac4 100644 --- a/src/math/math_NewtonMinimum.cxx +++ b/src/math/math_NewtonMinimum.cxx @@ -25,58 +25,49 @@ #include #include -//============================================================================ -math_NewtonMinimum::math_NewtonMinimum(math_MultipleVarFunctionWithHessian& F, - const math_Vector& StartingPoint, - const Standard_Real Tolerance, - const Standard_Integer NbIterations, - const Standard_Real Convexity, - const Standard_Boolean WithSingularity) -//============================================================================ - : TheLocation(1, F.NbVariables()), - TheGradient(1, F.NbVariables()), - TheStep(1, F.NbVariables(), 10*Tolerance), - TheHessian(1, F.NbVariables(), 1, F.NbVariables() ) +//======================================================================= +//function : math_NewtonMinimum +//purpose : Constructor +//======================================================================= +math_NewtonMinimum::math_NewtonMinimum( + const math_MultipleVarFunctionWithHessian& theFunction, + const Standard_Real theTolerance, + const Standard_Integer theNbIterations, + const Standard_Real theConvexity, + const Standard_Boolean theWithSingularity + ) +: TheStatus (math_NotBracketed), + TheLocation(1, theFunction.NbVariables()), + TheGradient(1, theFunction.NbVariables()), + TheStep (1, theFunction.NbVariables(), 10.0 * theTolerance), + TheHessian (1, theFunction.NbVariables(), 1, theFunction.NbVariables()), + PreviousMinimum (0.0), + TheMinimum (0.0), + MinEigenValue (0.0), + XTol (theTolerance), + CTol (theConvexity), + nbiter (0), + NoConvexTreatement(theWithSingularity), + Convex (Standard_True), + Done (Standard_False), + Itermax (theNbIterations) { - XTol = Tolerance; - CTol = Convexity; - Itermax = NbIterations; - NoConvexTreatement = WithSingularity; - Convex = Standard_True; -// Done = Standard_True; -// TheStatus = math_OK; - Perform ( F, StartingPoint); } -//============================================================================ -math_NewtonMinimum::math_NewtonMinimum(math_MultipleVarFunctionWithHessian& F, - const Standard_Real Tolerance, - const Standard_Integer NbIterations, - const Standard_Real Convexity, - const Standard_Boolean WithSingularity) -//============================================================================ - : TheLocation(1, F.NbVariables()), - TheGradient(1, F.NbVariables()), - TheStep(1, F.NbVariables(), 10*Tolerance), - TheHessian(1, F.NbVariables(), 1, F.NbVariables() ) -{ - XTol = Tolerance; - CTol = Convexity; - Itermax = NbIterations; - NoConvexTreatement = WithSingularity; - Convex = Standard_True; - Done = Standard_False; - TheStatus = math_NotBracketed; -} -//============================================================================ +//======================================================================= +//function : ~math_NewtonMinimum +//purpose : Destructor +//======================================================================= math_NewtonMinimum::~math_NewtonMinimum() { } -//============================================================================ +//======================================================================= +//function : Perform +//purpose : +//======================================================================= void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F, - const math_Vector& StartingPoint) -//============================================================================ + const math_Vector& StartingPoint) { math_Vector Point1 (1, F.NbVariables()); Point1 = StartingPoint; @@ -191,26 +182,20 @@ void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F, TheLocation = *precedent; } -//============================================================================ -Standard_Boolean math_NewtonMinimum::IsConverged() const -//============================================================================ -{ - return ( (TheStep.Norm() <= XTol ) || - ( Abs(TheMinimum-PreviousMinimum) <= XTol*Abs(PreviousMinimum) )); -} - -//============================================================================ +//======================================================================= +//function : Dump +//purpose : +//======================================================================= void math_NewtonMinimum::Dump(Standard_OStream& o) const -//============================================================================ { - o<< "math_Newton Optimisation: "; - o << " Done =" << Done << endl; - o << " Status = " << (Standard_Integer)TheStatus << endl; - o <<" Location Vector = " << Location() << endl; - o <<" Minimum value = "<< Minimum()<< endl; - o <<" Previous value = "<< PreviousMinimum << endl; - o <<" Number of iterations = " < #include -inline Standard_Boolean math_Powell::IsDone() const { return Done; } +inline void math_Powell::Delete() const +{ +} + +inline Standard_Boolean math_Powell::IsSolutionReached(math_MultipleVarFunction&) +{ + return 2.0 * fabs(PreviousMinimum - TheMinimum) <= + XTol * (fabs(PreviousMinimum) + fabs(TheMinimum) + EPSZ); +} + +inline Standard_Boolean math_Powell::IsDone() const +{ + return Done; +} inline Standard_OStream& operator<<(Standard_OStream& o, const math_Powell& P)