0027131: [Regression to 6.7] DistShapeShape performance loss

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

Possibility to freeze Lipchitz constant is added to improve performance.

diff --git a/src/Extrema/Extrema_GenExtCC.gxx b/src/Extrema/Extrema_GenExtCC.gxx
index 88ca463c7c7651feeecb2ac6fa4d4168599ddca1..3b564f4c74ce91d934e8fe10f2d23d29c12a7360 100644 (file)
--- a/src/Extrema/Extrema_GenExtCC.gxx
+++ b/src/Extrema/Extrema_GenExtCC.gxx
@@ -201,13 +201,36 @@ void Extrema_GenExtCC::Perform()
   C1.Intervals(anIntervals1, GeomAbs_C2);
   C2.Intervals(anIntervals2, GeomAbs_C2);
 
+  // Lipchitz constant approximation.
+  Standard_Real aLC = 9.0; // Default value.
+  Standard_Boolean isConstLockedFlag = Standard_False;
+  if (C1.GetType() == GeomAbs_Line)
+  {
+    Standard_Real aMaxDer = 1.0 / C2.Resolution(1.0);
+    if (aLC > aMaxDer)
+    {
+      isConstLockedFlag = Standard_True;
+      aLC = aMaxDer;
+    }
+  }
+  if (C2.GetType() == GeomAbs_Line)
+  {
+    Standard_Real aMaxDer = 1.0 / C1.Resolution(1.0);
+    if (aLC > aMaxDer)
+    {
+      isConstLockedFlag = Standard_True;
+      aLC = aMaxDer;
+    }
+  }
+
   Extrema_GlobOptFuncCCC2 aFunc (C1, C2);
-  math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder);
+  math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC);
+  aFinder.SetLipConstState(isConstLockedFlag);
   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 disntance cannot be lower than 0.0
+  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(),
@@ -287,7 +310,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;

diff --git a/src/math/math_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx
index 2560aa252d7953e7b78859f77888bd382ebf4028..bf75a6af713a6c2790514868eae1e3b8ec98a1e2 100644 (file)
--- a/src/math/math_GlobOptMin.cxx
+++ b/src/math/math_GlobOptMin.cxx
@@ -41,6 +41,7 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
   myB(1, myN),
   myGlobA(1, myN),
   myGlobB(1, myN),
+  myIsConstLocked(Standard_False),
   myX(1, myN),
   myTmp(1, myN),
   myV(1, myN),
@@ -53,6 +54,7 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
 
   myFunc = theFunc;
   myC = theC;
+  myInitC = theC;
   myIsFindSingleSolution = Standard_False;
   myFunctionalMinimalValue = -Precision::Infinite();
   myZ = -1;
@@ -85,13 +87,14 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
   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,
@@ -104,6 +107,7 @@ void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
 
   myFunc = theFunc;
   myC = theC;
+  myInitC = theC;
   myZ = -1;
   mySolCount = 0;
 
@@ -131,13 +135,14 @@ void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
   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)
@@ -145,8 +150,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);
@@ -227,8 +230,11 @@ void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
     return;
   }
 
-  // Compute initial values for myF, myY, myC.
-  computeInitialValues();
+  if (!myIsConstLocked)
+  {
+    // Compute initial value for myC.
+    computeInitialValues();
+  }
 
   myE1 = minLength * myTol;
   myE2 = maxLength * myTol;
@@ -236,8 +242,7 @@ void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
   myIsFindSingleSolution = isFindSingleSolution;
   if (isFindSingleSolution)
   {
-    // Run local optimization 
-    // if current value better than optimal.
+    // Run local optimization if current value better than optimal.
     myE3 = 0.0;
   }
   else
@@ -248,13 +253,14 @@ void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution)
       myE3 = - maxLength * myTol * myC / 4.0;
   }
 
-  // Search single solution and current solution in its neighbourhood.
+  // Search single solution and current solution in its neighborhood.
   if (CheckFunctionalStopCriteria())
   {
     myDone = Standard_True;
     return;
   }
 
+  myLastStep = 0.0;
   isFirstCellFilterInvoke = Standard_True;
   computeGlobalExtremum(myN);
 
@@ -338,35 +344,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
-  // 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
+  // Lipchitz const approximation.
   Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
   Standard_Integer aPntNb = 13;
   myFunc->Value(myA, aPrevValDiag);
@@ -389,16 +368,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)
+  else if (aLipConst > myC * 5)
+    myC = Min(myC * 5, 50.0);
+
+  // Clear all solutions except one.
+  if (myY.Size() != myN)
   {
-    myC = Min(myC * 2, 30.0);
+    for(i = 1; i <= myN; i++)
+      aBestPnt(i) = myY(i);
+    myY.Clear();
+    for(i = 1; i <= myN; i++)
+      myY.Append(aBestPnt(i));
   }
+  mySolCount = 1;
 }
 
 //=======================================================================
@@ -411,12 +397,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 aStepBestValue = RealLast();
-  Standard_Real aRealStep = 0.0;
   math_Vector aStepBestPoint(1, myN);
   Standard_Boolean isInside = Standard_False;
   Standard_Real r;
   Standard_Boolean isReached = Standard_False;
 
+
   for(myX(j) = myA(j) + myE1;
      (myX(j) < myB(j) + myE1) && (!isReached);
       myX(j) += myV(j))
@@ -434,7 +420,7 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
     {
       isInside = Standard_False;
       myFunc->Value(myX, d);
-      r = (d + myZ * myC * aRealStep - myF) * myZ;
+      r = (d + myZ * myC * myLastStep - myF) * myZ;
       if(r > myE3)
       {
         isInside = computeLocalExtremum(myX, val, myTmp);
@@ -477,8 +463,8 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
       if (CheckFunctionalStopCriteria())
         return; // Best possible value is obtained.
 
-      aRealStep = myE2 + Abs(myF - d) / myC;
-      myV(1) = Min(aRealStep, myMaxV(1));
+      myV(1) = Min(myE2 + Abs(myF - d) / myC, myMaxV(1));
+      myLastStep = myV(1);
     }
     else
     {
@@ -660,10 +646,59 @@ void math_GlobOptMin::initCellSize()
 //=======================================================================
 Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
 {
-  // Search single solution and current solution in its neighbourhood.
+  // 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;
+}
+
+//=======================================================================
+//function : SetLipConstState
+//purpose  :
+//=======================================================================
+void math_GlobOptMin::SetLipConstState(const Standard_Boolean theFlag)
+{
+  myIsConstLocked = theFlag;
+}
\ No newline at end of file

diff --git a/src/math/math_GlobOptMin.hxx b/src/math/math_GlobOptMin.hxx
index 1d42ab91e36fd58961425c300b5cc1e073077ec3..4c4e52096fb4eede3622aab559d0bf57dcd7e41f 100644 (file)
--- a/src/math/math_GlobOptMin.hxx
+++ b/src/math/math_GlobOptMin.hxx
@@ -146,6 +146,9 @@ public:
   //! Set functional minimal value.
   Standard_EXPORT void SetFunctionalMinimalValue(const Standard_Real theMinimalValue);
 
+  //! Lock/Unlock Lipchitz constant for internal modifications.
+  Standard_EXPORT void SetLipConstState(const Standard_Boolean theFlag);
+
   //! Get functional minimal value.
   Standard_EXPORT Standard_Real GetFunctionalMinimalValue();
 
@@ -154,6 +157,9 @@ 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);
@@ -188,9 +194,11 @@ 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
+  Standard_Boolean myIsConstLocked;  // Is constant locked for modifications.
 
   // Output.
   Standard_Boolean myDone;
@@ -207,6 +215,7 @@ private:
   math_Vector myTmp; // Current modified solution.
   math_Vector myV; // Steps array.
   math_Vector myMaxV; // Max Steps array.
+  Standard_Real myLastStep; // Last step.
   math_Vector myExpandCoeff; // Define expand coefficient between neighboring indiced dimensions.
 
   NCollection_Array1<Standard_Real> myCellSize;

diff --git a/tests/bugs/fclasses/bug27131 b/tests/bugs/fclasses/bug27131
new file mode 100644 (file)
index 0000000..7f833ff
--- /dev/null
+++ b/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
\ No newline at end of file

diff --git a/tests/bugs/modalg_5/bug23706_14 b/tests/bugs/modalg_5/bug23706_14
index 516eaee789da328c126b6f1321b4a7dbe6518d20..8cd60bb33e958d30eccd66f1847637536adccbe1 100755 (executable)
--- a/tests/bugs/modalg_5/bug23706_14
+++ b/tests/bugs/modalg_5/bug23706_14
@@ -13,7 +13,7 @@ set info [2dextrema b9 b10]
 set status 0
 for { set i 1 } { $i <= 1 } { incr i 1 } {
     regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp
-    if { abs($pp - 3.6712618987696386) > 1.0e-7 } {
+    if { abs($pp - 3.6710534171261284) > 1.0e-7 } {
        puts "Error : Extrema is wrong on dist $i"
        set status 1
     }