0027131: [Regression to 6.7] DistShapeShape performance loss

Lipchitz constant approximation and fixes in global optimization algorithm added to improve performance.
Test case added.

Fix backporting:

0026593: Coding rules - revert compatibility of NCollection_CellFilter constructor with old code

Restored old constructor and old behavior where possible.

Minor correction.

0026395: Merge clasees NCollection_CellFilter_NDim and NCollection_CellFilter

Deleted exceed class CellFilterNDim.
Now dimension count used as input parameter in NCollection_CellFilter.

minor corrections.

src/BRepMesh/BRepMesh_CircleTool.hxx

@@ -25,6 +25,7 @@
 #include <Standard_Integer.hxx>
 #include <Standard_Boolean.hxx>
 #include <BRepMesh.hxx>
+#include <NCollection_Array1.hxx>
 
 class gp_Circ2d;
 
@@ -67,7 +68,8 @@ public:
   inline void SetCellSize(const Standard_Real theSizeX,
                           const Standard_Real theSizeY)
   {
-    Standard_Real aCellSize[2] = { theSizeX, theSizeY };
+    Standard_Real aCellSizeC[2] = { theSizeX, theSizeY };
+    NCollection_Array1<Standard_Real> aCellSize(aCellSizeC[0], 1, 2);
     myCellFilter.Reset(aCellSize, myAllocator);
   }
 

src/BRepMesh/BRepMesh_VertexTool.hxx

@@ -14,6 +14,7 @@
 #ifndef _BRepMesh_VertexTool_HeaderFile
 #define _BRepMesh_VertexTool_HeaderFile
 
+#include <NCollection_Array1.hxx>
 #include <Standard.hxx>
 #include <Standard_DefineAlloc.hxx>
 #include <Standard_Macro.hxx>
@@ -54,7 +55,8 @@ public:
   Standard_EXPORT void SetCellSize(const Standard_Real theSizeX,
                                    const Standard_Real theSizeY)
   {
-    Standard_Real aCellSize[2] = { theSizeX, theSizeY };
+    Standard_Real aCellSizeC[2] = { theSizeX, theSizeY };
+    NCollection_Array1<Standard_Real> aCellSize(aCellSizeC[0], 1, 2);
     myCellFilter.Reset(aCellSize, myAllocator);
     mySelector.Clear();
   }

src/Extrema/Extrema_GenExtCC.gxx

@@ -198,12 +198,20 @@ void Extrema_GenExtCC::Perform()
   C1.Intervals(anIntervals1, GeomAbs_C2);
   C2.Intervals(anIntervals2, GeomAbs_C2);
 
+  // Lipchitz constant approximation.
+  Standard_Real aLC = 9.0; // Default value.
+  if (C1.GetType() == GeomAbs_Line)
+    aLC = Min(aLC, 1.0 / C2.Resolution(1.0));
+  if (C2.GetType() == GeomAbs_Line)
+    aLC = Min(aLC, 1.0 / C1.Resolution(1.0));
+
   Extrema_GlobOptFuncCCC2 aFunc (C1, C2);
-  math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder);
+  math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC);
   Standard_Real aDiscTol = 1.0e-2;
   Standard_Real aValueTol = 1.0e-2;
   Standard_Real aSameTol = myCurveMinTol / (aDiscTol);
   aFinder.SetTol(aDiscTol, aSameTol);
+  aFinder.SetFunctionalMinimalValue(0.0); // Best distance cannot be lower than 0.0.
 
   // Size computed to have cell index inside of int32 value.
   const Standard_Real aCellSize = Max(anIntervals1.Upper() - anIntervals1.Lower(),
@@ -283,7 +291,7 @@ void Extrema_GenExtCC::Perform()
   // Check for infinity solutions case, for this:
   // Sort points lexicographically and check midpoint between each two neighboring points.
   // If all midpoints functional value is acceptable
-  // then set myParallel flag to true and return one soulution.
+  // then set myParallel flag to true and return one solution.
   std::sort(aPnts.begin(), aPnts.end(), comp);
   Standard_Boolean isParallel = Standard_False;
   Standard_Real aVal = 0.0;

src/NCollection/NCollection_CellFilter.hxx

@@ -17,6 +17,8 @@
 #define NCollection_CellFilter_HeaderFile
 
 #include <Standard_Real.hxx>
+#include <NCollection_LocalArray.hxx>
+#include <NCollection_Array1.hxx>
 #include <NCollection_List.hxx>
 #include <NCollection_Map.hxx>
 #include <NCollection_DataMap.hxx>
@@ -110,8 +112,7 @@ enum NCollection_CellFilter_Action
  *   Note that method Inspect() can be const and/or virtual.
  */
 
-template <class Inspector> 
+template <class Inspector> class NCollection_CellFilter
-class NCollection_CellFilter
 {
 public:
   typedef TYPENAME Inspector::Target Target;
@@ -119,7 +120,7 @@ public:
 
 public:
 
-  //! Constructor; initialized by cell size.
+  //! Constructor; initialized by dimension count and cell size.
   //!
   //! Note: the cell size must be ensured to be greater than
   //! maximal co-ordinate of the involved points divided by INT_MAX,
@@ -128,19 +129,22 @@ public:
   //! By default cell size is 0, which is invalid; thus if default
   //! constructor is used, the tool must be initialized later with
   //! appropriate cell size by call to Reset()
-  NCollection_CellFilter (Standard_Real theCellSize=0,
+  //! Constructor when dimension count is unknown at compilation time.
-                          const Handle(NCollection_IncAllocator)& theAlloc=0) 
+  NCollection_CellFilter (const Standard_Integer theDim,
+                          const Standard_Real theCellSize = 0,
+                          const Handle(NCollection_IncAllocator)& theAlloc = 0)
+  : myCellSize(0, theDim - 1)
   {
+    myDim = theDim;
     Reset (theCellSize, theAlloc);
   }
 
-  //! Constructor; initialized by cell sizes along each dimension.
+  //! Constructor when dimenstion count is known at compilation time.
-  //! Note: the cell size in each dimension must be ensured to be greater than 
+  NCollection_CellFilter (const Standard_Real theCellSize = 0,
-  //! maximal co-ordinate of the involved points by this dimension divided by INT_MAX,
+                          const Handle(NCollection_IncAllocator)& theAlloc = 0)
-  //! in order to avoid integer overflow of cell index.
+  : myCellSize(0, Inspector::Dimension - 1)
-  NCollection_CellFilter (Standard_Real theCellSize[],
-                          const Handle(NCollection_IncAllocator)& theAlloc=0) 
   {
+    myDim = Inspector::Dimension;
     Reset (theCellSize, theAlloc);
   }
 
@@ -148,17 +152,16 @@ public:
   void Reset (Standard_Real theCellSize, 
               const Handle(NCollection_IncAllocator)& theAlloc=0)
   {
-    for (int i=0; i < Inspector::Dimension; i++)
+    for (int i=0; i < myDim; i++)
-      myCellSize[i] = theCellSize;    
+      myCellSize(i) = theCellSize;
     resetAllocator ( theAlloc );
   }
 
   //! Clear the data structures and set new cell sizes and allocator
-  void Reset (Standard_Real theCellSize[], 
+  void Reset (NCollection_Array1<Standard_Real> theCellSize, 
               const Handle(NCollection_IncAllocator)& theAlloc=0)
   {
-    for (int i=0; i < Inspector::Dimension; i++)
+    myCellSize = theCellSize;
-      myCellSize[i] = theCellSize[i];
     resetAllocator ( theAlloc );
   }
   
@@ -180,7 +183,7 @@ public:
     Cell aCellMax (thePntMax, myCellSize);
     Cell aCell = aCellMin;
     // add object recursively into all cells in range
-    iterateAdd (Inspector::Dimension-1, aCell, aCellMin, aCellMax, theTarget);
+    iterateAdd (myDim-1, aCell, aCellMin, aCellMax, theTarget);
   }
 
   //! Find a target object at a point and remove it from the structures.
@@ -204,7 +207,7 @@ public:
     Cell aCellMax (thePntMax, myCellSize);
     Cell aCell = aCellMin;
     // remove object recursively from all cells in range
-    iterateRemove (Inspector::Dimension-1, aCell, aCellMin, aCellMax, theTarget);
+    iterateRemove (myDim-1, aCell, aCellMin, aCellMax, theTarget);
   }
 
   //! Inspect all targets in the cell corresponding to the given point
@@ -225,7 +228,7 @@ public:
     Cell aCellMax (thePntMax, myCellSize);
     Cell aCell = aCellMin;
     // inspect object recursively into all cells in range
-    iterateInspect (Inspector::Dimension-1, aCell, 
+    iterateInspect (myDim-1, aCell, 
                     aCellMin, aCellMax, theInspector);
   }
 
@@ -238,7 +241,14 @@ protected:
   /**
    * Auxiliary class for storing points belonging to the cell as the list
    */
-  struct ListNode {
+  struct ListNode
+  {
+    ListNode()
+    {
+      // Empty constructor is forbidden.
+      Standard_NoSuchObject::Raise("NCollection_CellFilter::ListNode()");
+    }
+
     Target Object;
     ListNode *Next;
   };
@@ -251,17 +261,16 @@ protected:
   struct Cell
   {
   public:
-    //! Empty constructor -- required only for NCollection_Map,
-    //! therefore does not initialize index (avoid cycle)
-    Cell () : Objects(0) {}
 
     //! Constructor; computes cell indices
-    Cell (const Point& thePnt, const Standard_Real theCellSize[])
+    Cell (const Point& thePnt, 
-      : Objects(0)
+          const NCollection_Array1<Standard_Real>& theCellSize)
+      : index(theCellSize.Size()),
+        Objects(0)
     {
-      for (int i=0; i < Inspector::Dimension; i++)
+      for (int i = 0; i < theCellSize.Size(); i++)
       {
-          Standard_Real val = (Standard_Real)(Inspector::Coord(i, thePnt) / theCellSize[i]);
+          Standard_Real val = (Standard_Real)(Inspector::Coord(i, thePnt) / theCellSize(theCellSize.Lower() + i));
           //If the value of index is greater than
           //INT_MAX it is decreased correspondingly for the value of INT_MAX. If the value
           //of index is less than INT_MIN it is increased correspondingly for the absolute
@@ -273,13 +282,19 @@ protected:
     }
 
     //! Copy constructor: ensure that list is not deleted twice
-    Cell (const Cell& theOther) { (*this) = theOther; }
+    Cell (const Cell& theOther)
+      : index(theOther.index.Size())
+    {
+      (*this) = theOther;
+    }
 
     //! Assignment operator: ensure that list is not deleted twice
     void operator = (const Cell& theOther)
     {
-      for (int i=0; i < Inspector::Dimension; i++)
+      Standard_Integer myDim = Standard_Integer(theOther.index.Size());
-        index[i] = theOther.index[i];
+      for(Standard_Integer anIdx = 0; anIdx < myDim; anIdx++)
+        index[anIdx] = theOther.index[anIdx];
+
       Objects = theOther.Objects;
       ((Cell&)theOther).Objects = 0;
     }
@@ -296,7 +311,8 @@ protected:
     //! Compare cell with other one
     Standard_Boolean IsEqual (const Cell& theOther) const
     {
-      for (int i=0; i < Inspector::Dimension; i++) 
+      Standard_Integer myDim = Standard_Integer(theOther.index.Size());
+      for (int i=0; i < myDim; i++) 
         if ( index[i] != theOther.index[i] ) return Standard_False;
       return Standard_True;
     }
@@ -305,15 +321,16 @@ protected:
     Standard_Integer HashCode (const Standard_Integer theUpper) const
     {
       // number of bits per each dimension in the hash code
-      const Standard_Size aShiftBits = (BITS(long)-1) / Inspector::Dimension;
+      Standard_Integer myDim = Standard_Integer(index.Size());
+      const Standard_Size aShiftBits = (BITS(long)-1) / myDim;
       long aCode=0;
-      for (int i=0; i < Inspector::Dimension; i++)
+      for (int i=0; i < myDim; i++)
         aCode = ( aCode << aShiftBits ) ^ index[i];
       return (unsigned)aCode % theUpper;
     }
 
   public:
-    long index[Inspector::Dimension];
+    NCollection_LocalArray<long, 10> index;
     ListNode *Objects;
   };
 
@@ -452,9 +469,10 @@ protected:
   }
 
 protected:
+  Standard_Integer myDim;
   Handle(NCollection_BaseAllocator) myAllocator;
   NCollection_Map<Cell>             myCells;
-  Standard_Real                     myCellSize [Inspector::Dimension];
+  NCollection_Array1<Standard_Real> myCellSize;
 };
 
 /**
@@ -504,4 +522,3 @@ struct NCollection_CellFilter_InspectorXY
 };
 
 #endif
-

src/NCollection/NCollection_LocalArray.hxx

@@ -20,13 +20,10 @@
 
 //! Auxiliary class optimizing creation of array buffer 
 //! (using stack allocation for small arrays).
-template<class theItem> class NCollection_LocalArray
+template<class theItem, Standard_Integer MAX_ARRAY_SIZE = 1024> class NCollection_LocalArray
 {
 public:
 
-  // 1K * sizeof (theItem)
-  static const size_t MAX_ARRAY_SIZE = 1024;
-
   NCollection_LocalArray (const size_t theSize)
   : myPtr (myBuffer)
   {
@@ -34,7 +31,7 @@ public:
   }
 
   NCollection_LocalArray ()
-  : myPtr (myBuffer) {}
+  : myPtr (myBuffer), mySize(0) {}
 
   ~NCollection_LocalArray()
   {
@@ -48,6 +45,13 @@ public:
       myPtr = (theItem*)Standard::Allocate (theSize * sizeof(theItem));
     else
       myPtr = myBuffer;
+
+    mySize = theSize;
+  }
+
+  size_t Size() const
+  {
+    return mySize;
   }
 
   operator theItem*() const
@@ -72,6 +76,7 @@ protected:
 
   theItem  myBuffer[MAX_ARRAY_SIZE];
   theItem* myPtr;
+  size_t mySize;
 
 };
 

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,13 +45,17 @@ 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;
 
   myFunc = theFunc;
   myC = theC;
+  myInitC = theC;
   myIsFindSingleSolution = Standard_False;
+  myFunctionalMinimalValue = -Precision::Infinite();
   myZ = -1;
   mySolCount = 0;
 
@@ -78,12 +82,18 @@ 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();
+  ComputeInitSol();
+
   myDone = Standard_False;
 }
 
 //=======================================================================
 //function : SetGlobalParams
-//purpose  : Set params without memory allocation.
+//purpose  : Set parameters without memory allocation.
 //=======================================================================
 void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
                                       const math_Vector& theA,
@@ -96,6 +106,7 @@ void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
 
   myFunc = theFunc;
   myC = theC;
+  myInitC = theC;
   myZ = -1;
   mySolCount = 0;
 
@@ -122,12 +133,15 @@ void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
   myTol = theDiscretizationTol;
   mySameTol = theSameTol;
 
+  initCellSize();
+  ComputeInitSol();
+
   myDone = Standard_False;
 }
 
 //=======================================================================
 //function : SetLocalParams
-//purpose  : Set params without memory allocation.
+//purpose  : Set parameters without memory allocation.
 //=======================================================================
 void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
                                      const math_Vector& theLocalB)
@@ -135,8 +149,6 @@ void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
   Standard_Integer i;
 
   myZ = -1;
-  mySolCount = 0;
-
   for(i = 1; i <= myN; i++)
   {
     myA(i) = theLocalA(i);
@@ -217,7 +229,7 @@ void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
     return;
   }
 
-  // Compute initial values for myF, myY, myC.
+  // Compute initial value for myC.
   computeInitialValues();
 
   myE1 = minLength * myTol;
@@ -238,6 +250,14 @@ void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
       myE3 = - maxLength * myTol * myC / 4.0;
   }
 
+  // Search single solution and current solution in its neighborhood.
+  if (CheckFunctionalStopCriteria())
+  {
+    myDone = Standard_True;
+    return;
+  }
+
+  isFirstCellFilterInvoke = Standard_True;
   computeGlobalExtremum(myN);
 
   myDone = Standard_True;
@@ -256,11 +276,11 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt
   //Newton method
   if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
   {
-    math_MultipleVarFunctionWithHessian* myTmp = 
+    math_MultipleVarFunctionWithHessian* aTmp = 
       dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
-    math_NewtonMinimum newtonMinimum(*myTmp);
+    math_NewtonMinimum newtonMinimum(*aTmp);
     newtonMinimum.SetBoundary(myGlobA, myGlobB);
-    newtonMinimum.Perform(*myTmp, thePnt);
+    newtonMinimum.Perform(*aTmp, thePnt);
 
     if (newtonMinimum.IsDone())
     {
@@ -273,10 +293,10 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt
   // BFGS method used.
   if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
   {
-    math_MultipleVarFunctionWithGradient* myTmp = 
+    math_MultipleVarFunctionWithGradient* aTmp =
       dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
-    math_BFGS bfgs(myTmp->NbVariables());
+    math_BFGS bfgs(aTmp->NbVariables());
-    bfgs.Perform(*myTmp, thePnt);
+    bfgs.Perform(*aTmp, thePnt);
     if (bfgs.IsDone())
     {
       bfgs.Location(theOutPnt);
@@ -320,35 +340,8 @@ void math_GlobOptMin::computeInitialValues()
   math_Vector aBestPnt(1, myN);
   math_Vector aParamStep(1, myN);
   Standard_Real aCurrVal = RealLast();
-  Standard_Real aBestVal = RealLast();
 
-  // Check functional value in midpoint, low and upp point border and
+  // Lipchitz const approximation.
-  // in each point try to perform local optimization.
-  aBestPnt = (myA + myB) * 0.5;
-  myFunc->Value(aBestPnt, aBestVal);
-
-  for(i = 1; i <= 3; i++)
-  {
-    aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
-
-    if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
-    {
-      // Local Extremum finds better solution than current point.
-      if (aCurrVal < aBestVal)
-      {
-        aBestVal = aCurrVal;
-        aBestPnt = aCurrPnt;
-      }
-    }
-  }
-
-  myF = aBestVal;
-  myY.Clear();
-  for(i = 1; i <= myN; i++)
-    myY.Append(aBestPnt(i));
-  mySolCount++;
-
-  // Lipschitz const approximation
   Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
   Standard_Integer aPntNb = 13;
   myFunc->Value(myA, aPrevValDiag);
@@ -371,16 +364,23 @@ void math_GlobOptMin::computeInitialValues()
     aPrevValProj = aCurrVal;
   }
 
+  myC = myInitC;
   aLipConst *= Sqrt(myN) / aStep;
-
   if (aLipConst < myC * 0.1)
-  {
     myC = Max(aLipConst * 0.1, 0.01);
-  }
   else if (aLipConst > myC * 10)
-  {
     myC = Min(myC * 2, 30.0);
+
+  // Clear all solutions except one.
+  if (myY.Size() != myN)
+  {
+    for(i = 1; i <= myN; i++)
+      aBestPnt(i) = myY(i);
+    myY.Clear();
+    for(i = 1; i <= myN; i++)
+      myY.Append(aBestPnt(i));
   }
+  mySolCount = 1;
 }
 
 //=======================================================================
@@ -399,6 +399,7 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
   Standard_Real r;
   Standard_Boolean isReached = Standard_False;
 
+
   for(myX(j) = myA(j) + myE1;
      (myX(j) < myB(j) + myE1) && (!isReached);
       myX(j) += myV(j))
@@ -409,11 +410,14 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
       isReached = Standard_True;
     }
 
+    if (CheckFunctionalStopCriteria())
+      return; // Best possible value is obtained.
+
     if (j == 1)
     {
       isInside = Standard_False;
       myFunc->Value(myX, d);
-      r = (d + myZ * myC * aRealStep - myF) * myZ;
+      r = (d + myZ * myC * myV(1) - myF) * myZ;
       if(r > myE3)
       {
         isInside = computeLocalExtremum(myX, val, myTmp);
@@ -449,8 +453,13 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
         for(i = 1; i <= myN; i++)
           myY.Append(aStepBestPoint(i));
         mySolCount++;
+
+        isFirstCellFilterInvoke = Standard_True;
       }
 
+      if (CheckFunctionalStopCriteria())
+        return; // Best possible value is obtained.
+
       aRealStep = myE2 + Abs(myF - d) / myC;
       myV(1) = Min(aRealStep, myMaxV(1));
     }
@@ -505,20 +514,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(i = 0; i < mySolCount; i++)
-    for(j = 1; j <= myN; j++)
     {
-      if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
+      isSame = Standard_True;
+      for(j = 1; j <= myN; j++)
       {
-        isSame = Standard_False;
+        if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
-        break;
+        {
+          isSame = Standard_False;
+          break;
+        }
+      }
+      if (isSame == Standard_True)
+        return Standard_True;
+    }
+  }
+  else
+  {
+    NCollection_CellFilter_Inspector 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;
 }
@@ -541,6 +585,24 @@ Standard_Real math_GlobOptMin::GetF()
   return myF;
 }
 
+//=======================================================================
+//function : SetFunctionalMinimalValue
+//purpose  :
+//=======================================================================
+void math_GlobOptMin::SetFunctionalMinimalValue(const Standard_Real theMinimalValue)
+{
+  myFunctionalMinimalValue = theMinimalValue;
+}
+
+//=======================================================================
+//function : GetFunctionalMinimalValue
+//purpose  :
+//=======================================================================
+Standard_Real math_GlobOptMin::GetFunctionalMinimalValue()
+{
+  return myFunctionalMinimalValue;
+}
+
 //=======================================================================
 //function : IsDone
 //purpose  :
@@ -561,3 +623,70 @@ 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));
+  }
+}
+
+//=======================================================================
+//function : CheckFunctionalStopCriteria
+//purpose  :
+//=======================================================================
+Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
+{
+  // Search single solution and current solution in its neighborhood.
+  if (myIsFindSingleSolution &&
+      Abs (myF - myFunctionalMinimalValue) < mySameTol * 0.01)
+    return Standard_True;
+
+  return Standard_False;
+}
+
+//=======================================================================
+//function : ComputeInitSol
+//purpose  :
+//=======================================================================
+void math_GlobOptMin::ComputeInitSol()
+{
+  Standard_Real aCurrVal, aBestVal;
+  math_Vector aCurrPnt(1, myN);
+  math_Vector aBestPnt(1, myN);
+  math_Vector aParamStep(1, myN);
+  // Check functional value in midpoint, lower and upper border points and
+  // in each point try to perform local optimization.
+  aBestPnt = (myGlobA + myGlobB) * 0.5;
+  myFunc->Value(aBestPnt, aBestVal);
+
+  Standard_Integer i;
+  for(i = 1; i <= 3; i++)
+  {
+    aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
+
+    if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
+    {
+      // Local search tries to find better solution than current point.
+      if (aCurrVal < aBestVal)
+      {
+        aBestVal = aCurrVal;
+        aBestPnt = aCurrPnt;
+      }
+    }
+  }
+
+  myF = aBestVal;
+  myY.Clear();
+  for(i = 1; i <= myN; i++)
+    myY.Append(aBestPnt(i));
+  mySolCount = 1;
+
+  myDone = Standard_False;
+}

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_CellFilter.hxx>
 #include <math_MultipleVarFunction.hxx>
 #include <NCollection_Sequence.hxx>
 #include <Standard_Type.hxx>
 
+class NCollection_CellFilter_Inspector
+{
+public:
+
+  //! Points and target type
+  typedef math_Vector Point;
+  typedef math_Vector Target;
+
+  NCollection_CellFilter_Inspector(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.
@@ -65,16 +141,32 @@ public:
   //! Return solution theIndex, 1 <= theIndex <= NbExtrema.
   Standard_EXPORT void Points(const Standard_Integer theIndex, math_Vector& theSol);
 
-  Standard_Boolean isDone();
+  Standard_EXPORT Standard_Boolean isDone();
+
+  //! Set functional minimal value.
+  Standard_EXPORT void SetFunctionalMinimalValue(const Standard_Real theMinimalValue);
+
+  //! Get functional minimal value.
+  Standard_EXPORT Standard_Real GetFunctionalMinimalValue();
 
 private:
 
+  // Compute cell size.
+  void initCellSize();
+
+  // Compute initial solution
+  void ComputeInitSol();
+
   math_GlobOptMin & operator = (const math_GlobOptMin & theOther);
 
   Standard_Boolean computeLocalExtremum(const math_Vector& thePnt, Standard_Real& theVal, math_Vector& theOutPnt);
 
   void computeGlobalExtremum(Standard_Integer theIndex);
 
+  //! Check possibility to stop computations.
+  //! Find single solution + in neighbourhood of best possible solution.
+  Standard_Boolean CheckFunctionalStopCriteria();
+
   //! Computes starting value / approximation:
   //! myF - initial best value.
   //! myY - initial best point.
@@ -99,8 +191,10 @@ private:
   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
+  Standard_Real myC; //Lipchitz constant, default 9
+  Standard_Real myInitC; // Lipchitz constant initial value.
   Standard_Boolean myIsFindSingleSolution; // Default value is false.
+  Standard_Real myFunctionalMinimalValue; // Default value is -Precision::Infinite
 
   // Output.
   Standard_Boolean myDone;
@@ -119,6 +213,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_CellFilter<NCollection_CellFilter_Inspector> myFilter;
+
   Standard_Real myF; // Current value of Global optimum.
 };
 

tests/bugs/fclasses/bug27131

@@ -0,0 +1,32 @@
+puts "========"
+puts "OCC27131"
+puts "========"
+puts ""
+##############################################
+# DistShapeShape works slow on attached shapes
+##############################################
+restore [locate_data_file bug27131.brep] aShape
+explode aShape
+
+cpulimit 20
+
+# Check computation time
+chrono h reset; chrono h start
+for { set i 1 } { $i <= 100 } { incr i } {
+  distmini d aShape_1 aShape_2
+}
+chrono h stop; chrono h show
+
+regexp {CPU user time: (\d*)} [dchrono h show] dummy sec
+if {$sec > 1} {
+  puts "Error: too long computation time $sec seconds"
+} else {
+  puts "Computation time is OK"
+}
+
+# Check result of distance distance
+set absTol 1.0e-10
+set relTol 0.001
+set aDist_Exp 0.0029087110153708622
+set aDist [dval d_val]
+checkreal "Distance value check" $aDist $aDist_Exp $absTol $relTol