0025004: Extrema curve/curve incorrect result

Fixed bug in extrema clustering algorithm.
Tolerances changing is available now.
Testcase with Branin function added.

Test cases for issue CR25004

src/Extrema/Extrema_ExtCC.cxx

@@ -163,6 +163,7 @@ void Extrema_ExtCC::Perform()
   Standard_NullObject_Raise_if (!myC[0] || !myC[1], "Extrema_ExtCC::Perform()")
   myECC.SetParams(*((Adaptor3d_Curve*)myC[0]), 
                   *((Adaptor3d_Curve*)myC[1]), myInf[0], mySup[0], myInf[1], mySup[1]);
+  myECC.SetTolerance(Min(myTol[0], myTol[1]));
   myDone = Standard_False;
   mypoints.Clear();
   mySqDist.Clear();

src/Extrema/Extrema_GenExtCC.cdl

@@ -44,18 +44,21 @@ is
       	---Purpose: It calculates all the distances.
         --          The function F(u,v)=distance(C1(u),C2(v)) has an 
         --          extremum when gradient(f)=0. The algorithm uses
-        --          Evtushenko's global optimization solver..
+        --          Evtushenko's global optimization solver.
 
     Create (C1: Curve1; C2: Curve2; Uinf, Usup, Vinf, Vsup: Real) returns GenExtCC;
     	---Purpose: Calculates all the distances as above
     	--          between Uinf and Usup for C1 and  between Vinf and Vsup 
     	--          for C2.
 
-    SetParams(me: in out; C1: Curve1; C2: Curve2; Uinf, Usup, Vinf, Vsup: Real)
+    SetParams (me: in out; C1: Curve1; C2: Curve2; Uinf, Usup, Vinf, Vsup: Real)
       ---Purpose: Set params in case of empty constructor is usage.
       is static;
 
-    Perform(me: in out) is static;
+    SetTolerance (me: in out; Tol: Real);
+      ---Purpose:
+
+    Perform (me: in out) is static;
     	---Purpose: Performs calculations.
 
 
@@ -93,12 +96,13 @@ is
     	is static;
 
 fields
-	myLowBorder : Vector from math;
+  myCurveMinTol : Real from Standard;
-	myUppBorder : Vector from math;
+  myLowBorder : Vector from math;
-	mySolCount : Integer from Standard;
+  myUppBorder : Vector from math;
-	myPoints1 : SequenceOfReal from TColStd;
+  mySolCount : Integer from Standard;
-	myPoints2 : SequenceOfReal from TColStd;
+  myPoints1 : SequenceOfReal from TColStd;
-	myC : Address from Standard [2];
+  myPoints2 : SequenceOfReal from TColStd;
-    myDone : Boolean;
+  myC : Address from Standard [2];
+  myDone : Boolean;
 
 end GenExtCC;

src/Extrema/Extrema_GenExtCC.gxx

@@ -22,6 +22,10 @@
 #include <TColStd_Array1OfReal.hxx>
 #include <Precision.hxx>
 
+//=======================================================================
+//function : Extrema_GenExtCC
+//purpose  : 
+//=======================================================================
 Extrema_GenExtCC::Extrema_GenExtCC()
 : myLowBorder(1,2),
   myUppBorder(1,2),
@@ -29,8 +33,12 @@ Extrema_GenExtCC::Extrema_GenExtCC()
 {
 }
 
-Extrema_GenExtCC::Extrema_GenExtCC (const Curve1& C1,
+//=======================================================================
-                                    const Curve2& C2)
+//function : Extrema_GenExtCC
+//purpose  : 
+//=======================================================================
+Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
+                                   const Curve2& C2)
 : myLowBorder(1,2),
   myUppBorder(1,2),
   myDone(Standard_False)
@@ -41,15 +49,19 @@ Extrema_GenExtCC::Extrema_GenExtCC (const Curve1& C1,
   myLowBorder(2) = C2.FirstParameter();
   myUppBorder(1) = C1.LastParameter();
   myUppBorder(2) = C2.LastParameter();
+  myCurveMinTol = 1.0e-9;
 }
 
-
+//=======================================================================
-Extrema_GenExtCC::Extrema_GenExtCC (const Curve1& C1,
+//function : Extrema_GenExtCC
-                                    const Curve2& C2,
+//purpose  : 
-                                    const Standard_Real Uinf,
+//=======================================================================
-                                    const Standard_Real Usup,
+Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
-                                    const Standard_Real Vinf,
+                                   const Curve2& C2,
-                                    const Standard_Real Vsup)
+                                   const Standard_Real Uinf,
+                                   const Standard_Real Usup,
+                                   const Standard_Real Vinf,
+                                   const Standard_Real Vsup)
 : myLowBorder(1,2),
   myUppBorder(1,2),
   myDone(Standard_False)
@@ -60,14 +72,19 @@ Extrema_GenExtCC::Extrema_GenExtCC (const Curve1& C1,
   myLowBorder(2) = Vinf;
   myUppBorder(1) = Usup;
   myUppBorder(2) = Vsup;
+  myCurveMinTol = 1.0e-9;
 }
 
-void Extrema_GenExtCC::SetParams (const Curve1& C1,
+//=======================================================================
-                                  const Curve2& C2,
+//function : SetParams
-                                  const Standard_Real Uinf,
+//purpose  : 
-                                  const Standard_Real Usup,
+//=======================================================================
-                                  const Standard_Real Vinf,
+void Extrema_GenExtCC::SetParams(const Curve1& C1,
-                                  const Standard_Real Vsup)
+                                 const Curve2& C2,
+                                 const Standard_Real Uinf,
+                                 const Standard_Real Usup,
+                                 const Standard_Real Vinf,
+                                 const Standard_Real Vsup)
 {
   myC[0] = (Standard_Address)&C1;
   myC[1] = (Standard_Address)&C2;
@@ -77,8 +94,20 @@ void Extrema_GenExtCC::SetParams (const Curve1& C1,
   myUppBorder(2) = Vsup;
 }
 
-//=============================================================================
+//=======================================================================
-void Extrema_GenExtCC::Perform ()
+//function : SetTolerance
+//purpose  : 
+//=======================================================================
+void Extrema_GenExtCC::SetTolerance(Standard_Real theTol)
+{
+  myCurveMinTol = theTol;
+}
+
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void Extrema_GenExtCC::Perform()
 {
   myDone = Standard_False;
 
@@ -95,6 +124,10 @@ void Extrema_GenExtCC::Perform ()
 
   math_MultipleVarFunction *aFunc = new Extrema_GlobOptFuncCCC2(C1, C2);
   math_GlobOptMin aFinder(aFunc, myLowBorder, myUppBorder);
+  Standard_Real aDiscTol = 1.0e-2;
+  Standard_Real aValueTol = 1.0e-2;
+  Standard_Real aSameTol = myCurveMinTol / (aDiscTol);
+  aFinder.SetTol(aDiscTol, aSameTol);
 
   Standard_Integer i,j,k;
   math_Vector aFirstBorderInterval(1,2);
@@ -114,12 +147,15 @@ void Extrema_GenExtCC::Perform ()
       aFinder.Perform();
 
       aCurrF = aFinder.GetF();
-      if (aCurrF < aF + Precision::Confusion())
+      if (aCurrF < aF + aSameTol * aValueTol)
       {
-        if (aCurrF > Abs(aF - Precision::Confusion()) || (aCurrF < 1.0e-15 && aF < 1.0e-15))
+        if (aCurrF > aF - aSameTol * aValueTol)
         {
-          Standard_Integer myTmpSolCount = aFinder.NbExtrema();
+          if (aCurrF < aF)
+            aF = aCurrF;
+
           math_Vector sol(1,2);
+          Standard_Integer myTmpSolCount = aFinder.NbExtrema();
           for(k = 1; k <= myTmpSolCount; k++)
           {
             aFinder.Points(k, sol);
@@ -127,7 +163,7 @@ void Extrema_GenExtCC::Perform ()
             myPoints2.Append(sol(2));
           }
           mySolCount += myTmpSolCount;
-        } // if (aCurrF > aF - Precision::Confusion())
+        } // if (aCurrF > aF - aSameTol * aValueTol)
         else
         {
           aF = aCurrF;
@@ -142,7 +178,7 @@ void Extrema_GenExtCC::Perform ()
             myPoints2.Append(sol(2));
           }
         } // else
-      } //if (aCurrF < aF + Precision::Confusion())
+      } //if (aCurrF < aF + aSameTol * aValueTol)
     }
   }
 
@@ -151,10 +187,10 @@ void Extrema_GenExtCC::Perform ()
   {
     for(j = i + 1; j <= mySolCount; j++)
     {
-      if ((myPoints1(i) - myPoints1(j)) < (myUppBorder(1) - myLowBorder(1)) * Precision::Confusion() &&
+      if (Abs(myPoints1(i) - myPoints1(j)) < (myUppBorder(1) - myLowBorder(1)) * Precision::Confusion() &&
-          (myPoints2(i) - myPoints2(j)) < (myUppBorder(2) - myLowBorder(2)) * Precision::Confusion())
+          Abs(myPoints2(i) - myPoints2(j)) < (myUppBorder(2) - myLowBorder(2)) * Precision::Confusion())
       {
-        // Points with indexes i and j is in same cluster, delete j point from it.
+        // Points with indexes i and j is in same cluster, delete j point from extrema array.
         myPoints1.Remove(j);
         myPoints2.Remove(j);
         j--;
@@ -166,34 +202,46 @@ void Extrema_GenExtCC::Perform ()
   delete aFunc;
   myDone = Standard_True;
 }
-//=============================================================================
 
-Standard_Boolean Extrema_GenExtCC::IsDone () const 
+//=======================================================================
+//function : IsDone
+//purpose  : 
+//=======================================================================
+Standard_Boolean Extrema_GenExtCC::IsDone() const 
 {
   return myDone; 
 }
-//=============================================================================
 
-Standard_Integer Extrema_GenExtCC::NbExt () const
+//=======================================================================
+//function : NbExt
+//purpose  : 
+//=======================================================================
+Standard_Integer Extrema_GenExtCC::NbExt() const
 {
   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()")
 
   return mySolCount;
 }
-//=============================================================================
 
-Standard_Real Extrema_GenExtCC::SquareDistance (const Standard_Integer N) const
+//=======================================================================
+//function : SquareDistance
+//purpose  : 
+//=======================================================================
+Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const
 {
   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()")
   Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()")
 
   return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
 }
-//=============================================================================
 
-void Extrema_GenExtCC::Points (const Standard_Integer N,
+//=======================================================================
-                               POnC& P1,
+//function : Points
-                               POnC& P2) const
+//purpose  : 
+//=======================================================================
+void Extrema_GenExtCC::Points(const Standard_Integer N,
+                              POnC& P1,
+                              POnC& P2) const
 {
   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()")
   Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()")

src/QABugs/QABugs_19.cxx

@@ -2393,6 +2393,154 @@ static Standard_Integer OCC24889 (Draw_Interpretor& theDI,
   return 0;
 }
 
+#include <math_GlobOptMin.hxx>
+#include <math_MultipleVarFunctionWithHessian.hxx>
+//=======================================================================
+//function : OCC25004
+//purpose  : Check extremaCC on Branin function.
+//=======================================================================
+// Function is:
+// f(u,v) = a*(v - b*u^2 + c*u-r)^2+s(1-t)*cos(u)+s
+// Standard borders are:
+// -5 <= u <= 10
+//  0 <= v <= 15
+class BraninFunction : public math_MultipleVarFunctionWithHessian
+{
+public:
+  BraninFunction()
+  {
+    a = 1.0;
+    b = 5.1 / (4.0 * M_PI * M_PI);
+    c = 5.0 / M_PI;
+    r = 6.0;
+    s = 10.0;
+    t = 1.0 / (8.0 *  M_PI);
+  }
+  virtual Standard_Integer NbVariables() const
+  {
+    return 2;
+  }
+  virtual Standard_Boolean Value(const math_Vector& X,Standard_Real& F)
+  {
+    Standard_Real u = X(1);
+    Standard_Real v = X(2);
+
+    Standard_Real aSqPt = (v - b * u * u + c * u - r); // Square Part of function.
+    Standard_Real aLnPt = s * (1 - t) * cos(u); // Linear part of funcrtion.
+    F = a * aSqPt * aSqPt + aLnPt + s;
+    return Standard_True;
+  }
+  virtual Standard_Boolean Gradient(const math_Vector& X,math_Vector& G)
+  {
+    Standard_Real u = X(1);
+    Standard_Real v = X(2);
+
+    Standard_Real aSqPt = (v - b * u * u + c * u - r); // Square Part of function.
+    G(1) = 2 * a * aSqPt * (c - 2 * b * u) - s * (1 - t) * sin(u);
+    G(2) = 2 * a * aSqPt;
+
+    return Standard_True;
+  }
+  virtual Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G)
+  {
+    Value(X,F);
+    Gradient(X,G);
+
+    return Standard_True;
+  }
+  virtual Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G,math_Matrix& H)
+  {
+    Value(X,F);
+    Gradient(X,G);
+
+    Standard_Real u = X(1);
+    Standard_Real v = X(2);
+
+    Standard_Real aSqPt = (v - b * u * u + c * u - r); // Square Part of function.
+    Standard_Real aTmpPt = c - 2 * b *u; // Tmp part.
+    H(1,1) = 2 * a * aTmpPt * aTmpPt - 4 * a * b * aSqPt - s * (1 - t) * cos(u);
+    H(1,2) = 2 * a * aTmpPt;
+    H(2,1) = H(1,2);
+    H(2,2) = 2 * a;
+
+    return Standard_True;
+  }
+
+private:
+  // Standard parameters.
+  Standard_Real a, b, c, r, s, t;
+};
+
+static Standard_Integer OCC25004 (Draw_Interpretor& theDI,
+                                  Standard_Integer /*theNArg*/,
+                                  const char** /*theArgs*/)
+{
+  math_MultipleVarFunction* aFunc = new BraninFunction();
+
+  math_Vector aLower(1,2), aUpper(1,2);
+  aLower(1) = -5;
+  aLower(2) =  0;
+  aUpper(1) = 10;
+  aUpper(2) = 15;
+
+  Standard_Integer aGridOrder = 16;
+  math_Vector aFuncValues(1, aGridOrder * aGridOrder);
+
+  Standard_Real aLipConst = 0;
+  math_Vector aCurrPnt1(1, 2), aCurrPnt2(1, 2);
+
+  // Get Lipshitz constant estimation on regular grid.
+  Standard_Integer i, j, idx = 1;
+  for(i = 1; i <= aGridOrder; i++)
+  {
+    for(j = 1; j <= aGridOrder; j++)
+    {
+      aCurrPnt1(1) = aLower(1) + (aUpper(1) - aLower(1)) * (i - 1) / (aGridOrder - 1.0);
+      aCurrPnt1(2) = aLower(2) + (aUpper(2) - aLower(2)) * (j - 1) / (aGridOrder - 1.0);
+
+      aFunc->Value(aCurrPnt1, aFuncValues(idx));
+      idx++;
+    }
+  }
+
+  Standard_Integer k, l;
+  Standard_Integer idx1, idx2;
+  for(i = 1; i <= aGridOrder; i++)
+  for(j = 1; j <= aGridOrder; j++)
+  for(k = 1; k <= aGridOrder; k++)
+  for(l = 1; l <= aGridOrder; l++)
+    {
+      if (i == k && j == l) 
+        continue;
+
+      aCurrPnt1(1) = aLower(1) + (aUpper(1) - aLower(1)) * (i - 1) / (aGridOrder - 1.0);
+      aCurrPnt1(2) = aLower(2) + (aUpper(2) - aLower(2)) * (j - 1) / (aGridOrder - 1.0);
+      idx1 = (i - 1) * aGridOrder + j;
+
+      aCurrPnt2(1) = aLower(1) + (aUpper(1) - aLower(1)) * (k - 1) / (aGridOrder - 1.0);
+      aCurrPnt2(2) = aLower(2) + (aUpper(2) - aLower(2)) * (l - 1) / (aGridOrder - 1.0);
+      idx2 = (k - 1) * aGridOrder + l;
+
+      aCurrPnt1.Add(-aCurrPnt2);
+      Standard_Real dist = aCurrPnt1.Norm();
+
+      Standard_Real C = Abs(aFuncValues(idx1) - aFuncValues(idx2)) / dist;
+      if (C > aLipConst)
+        aLipConst = C;
+    }
+
+  math_GlobOptMin aFinder(aFunc, aLower, aUpper, aLipConst);
+  aFinder.Perform();
+  //(-pi , 12.275), (pi , 2.275), (9.42478, 2.475)
+
+  Standard_Real anExtValue = aFinder.GetF();
+  theDI << "F = " << anExtValue << "\n";
+
+  Standard_Integer aNbExt = aFinder.NbExtrema();
+  theDI << "NbExtrema = " << aNbExt << "\n";
+
+  return 0;
+}
 
 
 void QABugs::Commands_19(Draw_Interpretor& theCommands) {
@@ -2441,5 +2589,6 @@ void QABugs::Commands_19(Draw_Interpretor& theCommands) {
   theCommands.Add ("OCC24931", "OCC24931", __FILE__, OCC24931, group);
   theCommands.Add ("OCC24945", "OCC24945", __FILE__, OCC24945, group);
   theCommands.Add ("OCC23950", "OCC23950", __FILE__, OCC23950, group);
+  theCommands.Add ("OCC25004", "OCC25004", __FILE__, OCC25004, group);
   return;
 }

src/math/math_GlobOptMin.cxx

@@ -21,7 +21,6 @@
 #include <math_MultipleVarFunctionWithHessian.hxx>
 #include <math_NewtonMinimum.hxx>
 #include <math_Powell.hxx>
-#include <Precision.hxx>
 #include <Standard_Integer.hxx>
 #include <Standard_Real.hxx>
 
@@ -38,7 +37,9 @@ const Handle(Standard_Type)& STANDARD_TYPE(math_GlobOptMin)
 math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
                                  const math_Vector& theA,
                                  const math_Vector& theB,
-                                 Standard_Real theC)
+                                 const Standard_Real theC,
+                                 const Standard_Real theDiscretizationTol,
+                                 const Standard_Real theSameTol)
 : myN(theFunc->NbVariables()),
   myA(1, myN),
   myB(1, myN),
@@ -46,7 +47,8 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
   myGlobB(1, myN),
   myX(1, myN),
   myTmp(1, myN),
-  myV(1, myN)
+  myV(1, myN),
+  myMaxV(1, myN)
 {
   Standard_Integer i;
 
@@ -64,6 +66,14 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
     myB(i) = theB(i);
   }
 
+  for(i = 1; i <= myN; i++)
+  {
+    myMaxV(i) = (myB(i) - myA(i)) / 3.0;
+  }
+
+  myTol = theDiscretizationTol;
+  mySameTol = theSameTol;
+
   myDone = Standard_False;
 }
 
@@ -74,7 +84,9 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
                                       const math_Vector& theA,
                                       const math_Vector& theB,
-                                      Standard_Real theC)
+                                      const Standard_Real theC,
+                                      const Standard_Real theDiscretizationTol,
+                                      const Standard_Real theSameTol)
 {
   Standard_Integer i;
 
@@ -92,6 +104,9 @@ void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
     myB(i) = theB(i);
   }
 
+  myTol = theDiscretizationTol;
+  mySameTol = theSameTol;
+
   myDone = Standard_False;
 }
 
@@ -113,9 +128,36 @@ void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
     myB(i) = theLocalB(i);
   }
 
+  for(i = 1; i <= myN; i++)
+  {
+    myMaxV(i) = (myB(i) - myA(i)) / 3.0;
+  }
+
   myDone = Standard_False;
 }
 
+//=======================================================================
+//function : SetTol
+//purpose  : Set algorithm tolerances.
+//=======================================================================
+void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
+                             const Standard_Real theSameTol)
+{
+  myTol = theDiscretizationTol;
+  mySameTol = theSameTol;
+}
+
+//=======================================================================
+//function : GetTol
+//purpose  : Get algorithm tolerances.
+//=======================================================================
+void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
+                             Standard_Real& theSameTol)
+{
+  theDiscretizationTol = myTol;
+  theSameTol = mySameTol;
+}
+
 //=======================================================================
 //function : ~math_GlobOptMin
 //purpose  : 
@@ -145,13 +187,11 @@ void math_GlobOptMin::Perform()
       maxLength = currentLength;
   }
 
-  Standard_Real myTol = 1e-2;
+  myE1 = minLength * myTol       / myC;
+  myE2 = maxLength * myTol * 2.0 / myC;
+  myE3 = - maxLength * myTol / 4.0;
 
-  myE1 = minLength * myTol / myC;
+  // Compute start point.
-  myE2 = 2.0 * myTol / myC * maxLength;
-  myE3 = - maxLength * myTol / 4;
-
-  // Compure start point.
   math_Vector aPnt(1,myN);
   for(i = 1; i <= myN; i++)
   {
@@ -255,10 +295,12 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
   Standard_Real d; // Functional in moved point.
   Standard_Real val = RealLast(); // Local extrema computed in moved point.
   Standard_Real stepBestValue = RealLast();
-  math_Vector stepBestPoint(1,2);
+  Standard_Real realStep = RealLast();
+  math_Vector stepBestPoint(1, myN);
   Standard_Boolean isInside = Standard_False;
   Standard_Real r;
 
+
   for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j))
   {
     if (myX(j) > myB(j))
@@ -277,7 +319,7 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
       stepBestPoint = (isInside && (val < d))? myTmp : myX;
 
       // Solutions are close to each other.
-      if (Abs(stepBestValue - myF) < Precision::Confusion() * 0.01)
+      if (Abs(stepBestValue - myF) < mySameTol * 0.01)
       {
         if (!isStored(stepBestPoint))
         {
@@ -290,7 +332,7 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
       }
 
       // New best solution.
-      if ((stepBestValue - myF) * myZ > Precision::Confusion() * 0.01)
+      if ((stepBestValue - myF) * myZ > mySameTol * 0.01)
       {
         mySolCount = 0;
         myF = stepBestValue;
@@ -300,19 +342,20 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
         mySolCount++;
       }
 
-      myV(1) = myE2 + Abs(myF - d) / myC;
+      realStep = myE2 + Abs(myF - d) / myC;
+      myV(1) = Min(realStep, myMaxV(1));
     }
     else
     {
       myV(j) = RealLast() / 2.0;
       computeGlobalExtremum(j - 1);
     }
-    if ((j < myN) && (myV(j + 1) > myV(j)))
+    if ((j < myN) && (myV(j + 1) > realStep))
     {
-      if (myV(j) > (myB(j + 1) - myA(j + 1)) / 3.0) // Case of too big step.
+      if (realStep > myMaxV(j + 1)) // Case of too big step.
-        myV(j + 1) = (myB(j + 1) - myA(j + 1)) / 3.0; 
+        myV(j + 1) = myMaxV(j + 1); 
       else
-        myV(j + 1) = myV(j);
+        myV(j + 1) = realStep;
     }
   }
 }
@@ -347,7 +390,7 @@ Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
     isSame = Standard_True;
     for(j = 1; j <= myN; j++)
     {
-      if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) -  myA(j)) * Precision::Confusion())
+      if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) -  myA(j)) * mySameTol)
       {
         isSame = Standard_False;
         break;

src/math/math_GlobOptMin.hxx

@@ -31,16 +31,26 @@ public:
   Standard_EXPORT math_GlobOptMin(math_MultipleVarFunction* theFunc,
                                  const math_Vector& theA,
                                  const math_Vector& theB,
-                                 Standard_Real theC = 9);
+                                 const Standard_Real theC = 9,
+                                 const Standard_Real theDiscretizationTol = 1.0e-2,
+                                 const Standard_Real theSameTol = 1.0e-7);
 
   Standard_EXPORT void SetGlobalParams(math_MultipleVarFunction* theFunc,
                                        const math_Vector& theA,
                                        const math_Vector& theB,
-                                       Standard_Real theC = 9);
+                                       const Standard_Real theC = 9,
+                                       const Standard_Real theDiscretizationTol = 1.0e-2,
+                                       const Standard_Real theSameTol = 1.0e-7);
 
   Standard_EXPORT void SetLocalParams(const math_Vector& theLocalA,
                                       const math_Vector& theLocalB);
 
+  Standard_EXPORT void SetTol(const Standard_Real theDiscretizationTol,
+                              const Standard_Real theSameTol);
+
+  Standard_EXPORT void GetTol(Standard_Real& theDiscretizationTol,
+                              Standard_Real& theSameTol);
+
   Standard_EXPORT ~math_GlobOptMin();
 
   Standard_EXPORT void Perform();
@@ -54,6 +64,8 @@ public:
   //! Return solution i, 1 <= i <= NbExtrema.
   Standard_EXPORT void Points(const Standard_Integer theIndex, math_Vector& theSol);
 
+  Standard_Boolean isDone();
+
 private:
 
   math_GlobOptMin & operator = (const math_GlobOptMin & theOther);
@@ -66,8 +78,6 @@ private:
   Standard_Boolean isInside(const math_Vector& thePnt);
 
   Standard_Boolean isStored(const math_Vector &thePnt);
-  
-  Standard_Boolean isDone();
 
   // Input.
   math_MultipleVarFunction* myFunc;
@@ -76,6 +86,11 @@ private:
   math_Vector myB; // Right border on current C2 interval.
   math_Vector myGlobA; // Global left border.
   math_Vector myGlobB; // Global right border.
+  Standard_Real myTol; // Discretization tolerance, default 1.0e-2.
+  Standard_Real mySameTol; // points with ||p1 - p2|| < mySameTol is equal,
+                           // function values |val1 - val2| * 0.01 < mySameTol is equal,
+                           // default value is 1.0e-7.
+  Standard_Real myC; //Lipschitz constant, default 9
 
   // Output.
   Standard_Boolean myDone;
@@ -84,14 +99,14 @@ private:
 
   // Algorithm data.
   Standard_Real myZ;
-  Standard_Real myC; //Lipschitz constant
   Standard_Real myE1; // Border coeff.
   Standard_Real myE2; // Minimum step size.
   Standard_Real myE3; // Local extrema starting parameter.
 
-  math_Vector myX; // Current modified solution
+  math_Vector myX; // Current modified solution.
-  math_Vector myTmp; // Current modified solution
+  math_Vector myTmp; // Current modified solution.
   math_Vector myV; // Steps array.
+  math_Vector myMaxV; // Max Steps array.
 
   Standard_Real myF; // Current value of Global optimum.
 };

tests/bugs/modalg_5/bug23706_10

@@ -10,19 +10,12 @@ cpulimit 1500
 
 bsplinecurve r3 2 6 1 3 2 1 3 1 4 1 5 1 6 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 4 8 3 1 5 9 3 1 9 7 3 1
 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 3 1 3 9 3 1 5 7 3 1 7 4 3 1
-extrema r3 r4
 
-cvalue ext_1 0 x y z
+set info [extrema r3 r4] 
-set info [dump x] 
-regexp "(\[-0-9.+eE\])" $info full xx
-set info [dump y] 
-regexp "(\[-0-9.+eE\])" $info full yy
-set info [dump z] 
-regexp "(\[-0-9.+eE\])" $info full zz
 
-if { sqrt(($xx - 4.0) * ($xx - 4.0) + 
+regexp {Extrema 1 is point : } $info full pp
-	      ($yy - 8.0) * ($yy - 8.0) + 
+
-	      ($zz - 3.0) * ($zz - 3.0)) > 1.0e-7} {
+if { $pp == "4 8 3"} {
     puts "Error : Point of extrema is wrong"
 } else {
     puts "OK: Point of extrema is valid"

tests/bugs/modalg_5/bug23706_11

@@ -10,7 +10,7 @@ bsplinecurve r1 2 5 1 3 2 1 3 1 4 1 5 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 5 9 3 1
 bsplinecurve r2 2 5 2 3 2.5 1 3 1 3.5 1 4 3 -1 2 3 1 1 11 3 1 3 9 3 1 3 9 3 1 3 9 3 1 5 7 3 1 7 4 3 1
 set info [extrema r1 r2]
 
-if { [llength $info] != 1 } {
+if { [llength $info] != 3 } {
     puts "Error : Extrema is wrong"
 } else {
     puts "OK: Extrema is valid"

tests/bugs/modalg_5/bug25004_1

@@ -0,0 +1,22 @@
+puts "============"
+puts "OCC25004"
+puts "============"
+puts ""
+##########################################################################################################
+# Extrema_ExtCC::Extrema curve/curve incorrect result
+##########################################################################################################
+
+pload QAcommands
+
+set info [OCC25004]
+regexp {F += +([-0-9.+eE]+)} $info full aF
+regexp {NbExtrema += +([-0-9.+eE]+)} $info full aNb
+
+set expected_F 0.39788735772
+set expected_aNb 3
+set tol_abs_dist 1.0e-12
+set tol_rel_dist 0.1
+
+checkreal "value F" ${aF} ${expected_F} ${tol_abs_dist} ${tol_rel_dist}
+checkreal "Nbextrema" ${aNb} ${expected_aNb} ${tol_abs_dist} ${tol_rel_dist}
+

tests/bugs/modalg_5/bug25004_2

@@ -0,0 +1,22 @@
+puts "=========="
+puts "OCC25004  "
+puts "=========="
+puts ""
+##################################################
+## Extrema curve/curve incorrect result
+##################################################
+
+restore [locate_data_file bug25004_c1.draw] c1
+restore [locate_data_file bug25004_c2.draw] c2
+
+set list [extrema c1 c2]
+set nb [llength ${list}]
+if { ${nb} == 2} {
+    puts "OK : command extrema works properly"
+} else {
+    puts "Error : command extrema does NOT work properly"
+}
+
+smallview
+fit
+set only_screen_axo 1

tests/bugs/modalg_5/bug25004_3

@@ -0,0 +1,22 @@
+puts "=========="
+puts "OCC25004  "
+puts "=========="
+puts ""
+##################################################
+## Extrema curve/curve incorrect result
+##################################################
+
+restore [locate_data_file bug25004_bc1.draw] c1
+restore [locate_data_file bug25004_bc2.draw] c2
+
+set list [extrema c1 c2]
+set nb [llength ${list}]
+if { ${nb} == 2} {
+    puts "OK : command extrema works properly"
+} else {
+    puts "Error : command extrema does NOT work properly"
+}
+
+smallview
+fit
+set only_screen_axo 1