0026184: GeomAPI_ExtremaCurveCurve hangs on parallel b-spline curves

Class CellFilterNDim added.
Class CellFulterNDim used in GlobOptMin to improve performance in case of many solutions.
Memory leak eliminated.

Test cases for issue CR26184

Small correction of test cases for issue CR26184

src/math/math_GlobOptMin.cxx

@@ -1,6 +1,6 @@
 // Created on: 2014-01-20
 // Created by: Alexaner Malyshev
-// Copyright (c) 2014-2014 OPEN CASCADE SAS
+// Copyright (c) 2014-2015 OPEN CASCADE SAS
 //
 // This file is part of Open CASCADE Technology software library.
 //
@@ -45,7 +45,9 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
   myTmp(1, myN),
   myV(1, myN),
   myMaxV(1, myN),
-  myExpandCoeff(1, myN)
+  myExpandCoeff(1, myN),
+  myCellSize(0, myN - 1),
+  myFilter(theFunc->NbVariables())
 {
   Standard_Integer i;
 
@@ -78,6 +80,11 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
   myTol = theDiscretizationTol;
   mySameTol = theSameTol;
 
+  const Standard_Integer aMaxSquareSearchSol = 200;
+  Standard_Integer aSolNb = Standard_Integer(Pow(3.0, Standard_Real(myN)));
+  myMinCellFilterSol = Max(2 * aSolNb, aMaxSquareSearchSol);
+  initCellSize();
+
   myDone = Standard_False;
 }
 
@@ -122,6 +129,8 @@ void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
   myTol = theDiscretizationTol;
   mySameTol = theSameTol;
 
+  initCellSize();
+
   myDone = Standard_False;
 }
 
@@ -238,6 +247,7 @@ void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
       myE3 = - maxLength * myTol * myC / 4.0;
   }
 
+  isFirstCellFilterInvoke = Standard_True;
   computeGlobalExtremum(myN);
 
   myDone = Standard_True;
@@ -398,7 +408,6 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
   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))
@@ -444,6 +453,8 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
         for(i = 1; i <= myN; i++)
           myY.Append(aStepBestPoint(i));
         mySolCount++;
+
+        isFirstCellFilterInvoke = Standard_True;
       }
 
       aRealStep = myE2 + Abs(myF - d) / myC;
@@ -500,20 +511,55 @@ Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
   math_Vector aTol(1, myN);
   aTol = (myB -  myA) * mySameTol;
 
-  for(i = 0; i < mySolCount; i++)
+  // C1 * n^2 = C2 * 3^dim * n
+  if (mySolCount < myMinCellFilterSol)
   {
-    isSame = Standard_True;
-    for(j = 1; j <= myN; j++)
+    for(i = 0; i < mySolCount; i++)
     {
-      if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
+      isSame = Standard_True;
+      for(j = 1; j <= myN; j++)
       {
-        isSame = Standard_False;
-        break;
+        if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
+        {
+          isSame = Standard_False;
+          break;
+        }
+      }
+      if (isSame == Standard_True)
+        return Standard_True;
+    }
+  }
+  else
+  {
+    NCollection_CellFilter_NDimInspector anInspector(myN, Precision::PConfusion());
+    if (isFirstCellFilterInvoke)
+    {
+      myFilter.Reset(myCellSize);
+
+      // Copy initial data into cell filter.
+      for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++)
+      {
+        math_Vector aVec(1, myN);
+        for(Standard_Integer aSolDim = 1; aSolDim <= myN; aSolDim++)
+          aVec(aSolDim) = myY(aSolIdx * myN + aSolDim);
+
+        myFilter.Add(aVec, aVec);
       }
     }
-    if (isSame == Standard_True)
-      return Standard_True;
 
+    isFirstCellFilterInvoke = Standard_False;
+
+    math_Vector aLow(1, myN), anUp(1, myN);
+    anInspector.Shift(thePnt, myCellSize, aLow, anUp);
+
+    anInspector.ClearFind();
+    anInspector.SetCurrent(thePnt);
+    myFilter.Inspect(aLow, anUp, anInspector);
+    if (!anInspector.isFind())
+    {
+      // Point is out of close cells, add new one.
+      myFilter.Add(thePnt, thePnt);
+    }
   }
   return Standard_False;
 }
@@ -556,3 +602,16 @@ void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSo
   for(j = 1; j <= myN; j++)
     theSol(j) = myY((theIndex - 1) * myN + j);
 }
+
+//=======================================================================
+//function : initCellSize
+//purpose  :
+//=======================================================================
+void math_GlobOptMin::initCellSize()
+{
+  for(Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
+  {
+    myCellSize(anIdx - 1) = (myGlobB(anIdx) - myGlobA(anIdx))
+      * Precision::PConfusion() / (2.0 * Sqrt(2.0));
+  }
+}

src/math/math_GlobOptMin.hxx

@@ -1,6 +1,6 @@
 // Created on: 2014-01-20
 // Created by: Alexaner Malyshev
-// Copyright (c) 2014-2014 OPEN CASCADE SAS
+// Copyright (c) 2014-2015 OPEN CASCADE SAS
 //
 // This file is part of Open CASCADE Technology software library.
 //
@@ -16,10 +16,86 @@
 #ifndef _math_GlobOptMin_HeaderFile
 #define _math_GlobOptMin_HeaderFile
 
+#include <gp_Pnt.hxx>
+#include <gp_Pnt2d.hxx>
+#include <NCollection_CellFilterNDim.hxx>
 #include <math_MultipleVarFunction.hxx>
 #include <NCollection_Sequence.hxx>
 #include <Standard_Type.hxx>
 
+class NCollection_CellFilter_NDimInspector
+{
+public:
+
+  //! Points and target type
+  typedef math_Vector Point;
+  typedef math_Vector Target;
+
+  NCollection_CellFilter_NDimInspector(const Standard_Integer theDim,
+                                       const Standard_Real theTol)
+  : myCurrent(1, theDim)
+  {
+    myTol = theTol * theTol;
+    myIsFind = Standard_False;
+    Dimension = theDim;
+  }
+
+  //! Access to co-ordinate
+  static Standard_Real Coord (int i, const Point &thePnt)
+  {
+    return thePnt(i + 1);
+  }
+
+  //! Auxiliary method to shift point by each coordinate on given value;
+  //! useful for preparing a points range for Inspect with tolerance
+  void Shift (const Point& thePnt,
+              const NCollection_Array1<Standard_Real> &theTol,
+              Point& theLowPnt,
+              Point& theUppPnt) const
+  {
+    for(Standard_Integer anIdx = 1; anIdx <= Dimension; anIdx++)
+    {
+      theLowPnt(anIdx) = thePnt(anIdx) - theTol(anIdx - 1);
+      theUppPnt(anIdx) = thePnt(anIdx) + theTol(anIdx - 1);
+    }
+  }
+
+  void ClearFind()
+  {
+    myIsFind = Standard_False;
+  }
+
+  Standard_Boolean isFind()
+  {
+    return myIsFind;
+  }
+
+  //! Set current point to search for coincidence
+  void SetCurrent (const math_Vector& theCurPnt)
+  { 
+    myCurrent = theCurPnt;
+  }
+
+  //! Implementation of inspection method
+  NCollection_CellFilter_Action Inspect (const Target& theObject)
+  {
+    Standard_Real aSqDist = (myCurrent - theObject).Norm2();
+
+    if(aSqDist < myTol)
+    {
+      myIsFind = Standard_True;
+    }
+
+    return CellFilter_Keep;
+  }
+
+private:
+  Standard_Real myTol;
+  math_Vector myCurrent;
+  Standard_Boolean myIsFind;
+  Standard_Integer Dimension;
+};
+
 //! This class represents Evtushenko's algorithm of global optimization based on nonuniform mesh.<br>
 //! Article: Yu. Evtushenko. Numerical methods for finding global extreme (case of a non-uniform mesh). <br>
 //! U.S.S.R. Comput. Maths. Math. Phys., Vol. 11, N 6, pp. 38-54.
@@ -69,6 +145,9 @@ public:
 
 private:
 
+  // Compute cell size.
+  void initCellSize();
+
   math_GlobOptMin & operator = (const math_GlobOptMin & theOther);
 
   Standard_Boolean computeLocalExtremum(const math_Vector& thePnt, Standard_Real& theVal, math_Vector& theOutPnt);
@@ -119,6 +198,11 @@ private:
   math_Vector myMaxV; // Max Steps array.
   math_Vector myExpandCoeff; // Define expand coefficient between neighboring indiced dimensions.
 
+  NCollection_Array1<Standard_Real> myCellSize;
+  Standard_Integer myMinCellFilterSol;
+  Standard_Boolean isFirstCellFilterInvoke;
+  NCollection_CellFilterNDim<NCollection_CellFilter_NDimInspector> myFilter;
+
   Standard_Real myF; // Current value of Global optimum.
 };