0025058: Regression of performance of BRepExtrema_ExtCC (1000 times slower)

Added initial values approximation to improve performance.
Local optimization start coefficient fixed.

Test case for issue CR25058

src/math/math_GlobOptMin.cxx

@@ -175,6 +175,9 @@ void math_GlobOptMin::Perform()
 {
   Standard_Integer i;
 
+  // Compute initial values for myF, myY, myC.
+  computeInitialValues();
+
   // Compute parameters range
   Standard_Real minLength = RealLast();
   Standard_Real maxLength = RealFirst();
@@ -187,36 +190,12 @@ void math_GlobOptMin::Perform()
       maxLength = currentLength;
   }
 
-  myE1 = minLength * myTol       / myC;
-  myE2 = maxLength * myTol * 2.0 / myC;
-  myE3 = - maxLength * myTol / 4.0;
-
-  // Compute start point.
-  math_Vector aPnt(1,myN);
-  for(i = 1; i <= myN; i++)
-  {
-    Standard_Real currCentral = (myA(i) + myB(i)) / 2.0;
-    aPnt(i) = currCentral;
-  }
-
-  myFunc->Value(aPnt, myF);
-
-  math_Vector aExtremumPoint(1,myN);
-  Standard_Real aExtremumValue = RealLast();
-  if (computeLocalExtremum(aPnt, aExtremumValue, aExtremumPoint))
-  {
-    // Local Extremum finds better solution than midpoint.
-    if (aExtremumValue < myF)
-    {
-      myF = aExtremumValue;
-      aPnt = aExtremumPoint;
-    }
-  }
-
-  myY.Clear();
-  for(i = 1; i <= myN; i++)
-    myY.Append(aPnt(i));
-  mySolCount++;
+  myE1 = minLength * myTol;
+  myE2 = maxLength * myTol;
+  if (myC > 1.0)
+    myE3 = - maxLength * myTol / 4.0;
+  else
+    myE3 = - maxLength * myTol * myC / 4.0;
 
   computeGlobalExtremum(myN);
 
@@ -285,6 +264,72 @@ Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt
     return Standard_False;
 }
 
+//=======================================================================
+//function : computeInitialValues
+//purpose  : 
+//=======================================================================
+void math_GlobOptMin::computeInitialValues()
+{
+  Standard_Integer i;
+  math_Vector aCurrPnt(1, myN);
+  math_Vector aBestPnt(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
+  Standard_Real aLipConst = 0.0, aPrevVal;
+  Standard_Integer aPntNb = 13;
+  myFunc->Value(myA, aPrevVal);
+  Standard_Real aStep = (myB - myA).Norm() / aPntNb;
+  for(i = 1; i <= aPntNb; i++)
+  {
+    aCurrPnt = myA + (myB - myA) * i / (aPntNb - 1);
+    myFunc->Value(aCurrPnt, aCurrVal);
+
+    if(Abs(aCurrVal - aPrevVal) / aStep > aLipConst)
+      aLipConst = Abs(aCurrVal - aPrevVal) / aStep;
+
+    aPrevVal = aCurrVal;
+  }
+  aLipConst *= Sqrt(myN);
+
+  if (aLipConst < myC * 0.1)
+  {
+    myC = Max(aLipConst * 0.1, 0.01);
+  }
+  else if (aLipConst > myC * 10)
+  {
+    myC = Min(myC * 2, 30.0);
+  }
+}
+
 //=======================================================================
 //function : ComputeGlobalExtremum
 //purpose  :

src/math/math_GlobOptMin.hxx

@@ -58,7 +58,7 @@ public:
   //! Get best functional value.
   Standard_EXPORT Standard_Real GetF();
 
-  //! Return count of global extremas. NbExtrema <= MAX_SOLUTIONS.
+  //! Return count of global extremas.
   Standard_EXPORT Standard_Integer NbExtrema();
 
   //! Return solution i, 1 <= i <= NbExtrema.
@@ -74,6 +74,13 @@ private:
 
   void computeGlobalExtremum(Standard_Integer theIndex);
 
+  //! Computes starting value / approximation:
+  // myF - initial best value.
+  // myY - initial best point.
+  // myC - approximation of Lipschitz constant.
+  // to imporve convergence speed.
+  void computeInitialValues();
+
   //! Check that myA <= pnt <= myB
   Standard_Boolean isInside(const math_Vector& thePnt);