From: aml Date: Tue, 12 Apr 2016 09:15:07 +0000 (+0300) Subject: 0027371: Regression: BRepExtrema works too much slower in 691 (from 670) X-Git-Tag: V7_0_winwerth~74 X-Git-Url: http://git.dev.opencascade.org/gitweb/?p=occt.git;a=commitdiff_plain;h=246c7a7554dd24047ada720fbfce42e9c10b4657;hp=4b5857d330d381ea52951f7b64336b455187a948 0027371: Regression: BRepExtrema works too much slower in 691 (from 670) I Lipschitz constant tuning. Shubert estimation for minimal value is added. II Math_GlobOptMin Refactoring. III Test case is added. class NCollection_CellFilter_Inspector moved into math_GlobOptMin class. --- diff --git a/src/Extrema/Extrema_GenExtCC.gxx b/src/Extrema/Extrema_GenExtCC.gxx index f3604aa333..bb28481f24 100644 --- a/src/Extrema/Extrema_GenExtCC.gxx +++ b/src/Extrema/Extrema_GenExtCC.gxx @@ -210,8 +210,15 @@ void Extrema_GenExtCC::Perform() C1.Intervals(anIntervals1, aContinuity); C2.Intervals(anIntervals2, aContinuity); - // Lipchitz constant approximation. + // Lipchitz constant computation. Standard_Real aLC = 9.0; // Default value. + const Standard_Real aMaxDer1 = 1.0 / C1.Resolution(1.0); + const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0); + const Standard_Real aMaxDer = Max(aMaxDer1, aMaxDer2) * Sqrt(2.0); + if (aLC > aMaxDer) + aLC = aMaxDer; + + // Change constant value according to the concrete curve types. Standard_Boolean isConstLockedFlag = Standard_False; if (C1.GetType() == GeomAbs_Line) { diff --git a/src/math/math_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx index 6c1b23e997..abe10831c9 100644 --- a/src/math/math_GlobOptMin.cxx +++ b/src/math/math_GlobOptMin.cxx @@ -25,6 +25,27 @@ #include #include +//======================================================================= +//function : DistanceToBorder +//purpose : +//======================================================================= +static Standard_Real DistanceToBorder(const math_Vector & theX, + const math_Vector & theMin, + const math_Vector & theMax) +{ + Standard_Real aDist = RealLast(); + + for (Standard_Integer anIdx = theMin.Lower(); anIdx <= theMin.Upper(); ++anIdx) + { + const Standard_Real aDist1 = Abs (theX(anIdx) - theMin(anIdx)); + const Standard_Real aDist2 = Abs (theX(anIdx) - theMax(anIdx)); + + aDist = Min (aDist, Min (aDist1, aDist2)); + } + + return aDist; +} + //======================================================================= //function : math_GlobOptMin @@ -46,7 +67,6 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc, myTmp(1, myN), myV(1, myN), myMaxV(1, myN), - myExpandCoeff(1, myN), myCellSize(0, myN - 1), myFilter(theFunc->NbVariables()), myCont(2) @@ -75,12 +95,6 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc, myMaxV(i) = (myB(i) - myA(i)) / 3.0; } - myExpandCoeff(1) = 1.0; - for(i = 2; i <= myN; i++) - { - myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1)); - } - myTol = theDiscretizationTol; mySameTol = theSameTol; @@ -126,12 +140,6 @@ void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc, myMaxV(i) = (myB(i) - myA(i)) / 3.0; } - myExpandCoeff(1) = 1.0; - for(i = 2; i <= myN; i++) - { - myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1)); - } - myTol = theDiscretizationTol; mySameTol = theSameTol; @@ -162,12 +170,6 @@ void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA, myMaxV(i) = (myB(i) - myA(i)) / 3.0; } - myExpandCoeff(1) = 1.0; - for(i = 2; i <= myN; i++) - { - myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1)); - } - myDone = Standard_False; } @@ -193,14 +195,6 @@ void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol, theSameTol = mySameTol; } -//======================================================================= -//function : ~math_GlobOptMin -//purpose : -//======================================================================= -math_GlobOptMin::~math_GlobOptMin() -{ -} - //======================================================================= //function : Perform //purpose : Compute Global extremum point @@ -208,18 +202,20 @@ math_GlobOptMin::~math_GlobOptMin() // In this algo indexes started from 1, not from 0. void math_GlobOptMin::Perform(const Standard_Boolean isFindSingleSolution) { - Standard_Integer i; + myDone = Standard_False; // Compute parameters range Standard_Real minLength = RealLast(); Standard_Real maxLength = RealFirst(); - for(i = 1; i <= myN; i++) + for(Standard_Integer i = 1; i <= myN; i++) { Standard_Real currentLength = myB(i) - myA(i); if (currentLength < minLength) minLength = currentLength; if (currentLength > maxLength) maxLength = currentLength; + + myV(i) = 0.0; } if (minLength < Precision::PConfusion()) @@ -397,18 +393,16 @@ void math_GlobOptMin::computeInitialValues() void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j) { Standard_Integer i; - Standard_Real d; // Functional in moved point. + Standard_Real d = RealLast(), aPrevVal; // Functional in original and moved points. Standard_Real val = RealLast(); // Local extrema computed in moved point. Standard_Real aStepBestValue = RealLast(); math_Vector aStepBestPoint(1, myN); - Standard_Boolean isInside = Standard_False; - Standard_Real r; - Standard_Boolean isReached = Standard_False; + Standard_Boolean isInside = Standard_False, + isReached = Standard_False; + Standard_Real r1, r2, r; - for(myX(j) = myA(j) + myE1; - (myX(j) < myB(j) + myE1) && (!isReached); - myX(j) += myV(j)) + for(myX(j) = myA(j) + myE1; !isReached; myX(j) += myV(j)) { if (myX(j) > myB(j)) { @@ -422,11 +416,30 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j) if (j == 1) { isInside = Standard_False; + aPrevVal = d; myFunc->Value(myX, d); - r = (d + myZ * myC * myLastStep - myF) * myZ; + r1 = (d + myZ * myC * myLastStep - myF) * myZ; // Evtushenko estimation. + r2 = ((d + aPrevVal - myC * myLastStep) * 0.5 - myF) * myZ; // Shubert / Piyavsky estimation. + r = Min(r1, r2); if(r > myE3) { - isInside = computeLocalExtremum(myX, val, myTmp); + Standard_Real aSaveParam = myX(1); + + // Piyavsky midpoint estimation. + Standard_Real aParam = (2 * myX(1) - myV(1) ) * 0.5 + (aPrevVal - d) * 0.5 / myC; + if (Precision::IsInfinite(aPrevVal)) + aParam = myX(1) - myV(1) * 0.5; // Protection from upper dimension step. + + myX(1) = aParam; + Standard_Real aVal = 0; + myFunc->Value(myX, aVal); + myX(1) = aSaveParam; + + if ( (aVal < d && aVal < aPrevVal) || + DistanceToBorder(myX, myA, myB) < myE1 ) // Condition optimization case near the border. + { + isInside = computeLocalExtremum(myX, val, myTmp); + } } aStepBestValue = (isInside && (val < d))? val : d; aStepBestPoint = (isInside && (val < d))? myTmp : myX; @@ -478,10 +491,9 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j) for(i = 1; i < j; i++) myV(i) = 0.0; } - // Compute step in (j + 1) dimension according to scale. if (j < myN) { - Standard_Real aUpperDimStep = myV(j) * myExpandCoeff(j + 1); + Standard_Real aUpperDimStep = Max(myV(j), myE2); if (myV(j + 1) > aUpperDimStep) { if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step. @@ -573,51 +585,6 @@ Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt) return Standard_False; } -//======================================================================= -//function : NbExtrema -//purpose : -//======================================================================= -Standard_Integer math_GlobOptMin::NbExtrema() -{ - return mySolCount; -} - -//======================================================================= -//function : GetF -//purpose : -//======================================================================= -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 : -//======================================================================= -Standard_Boolean math_GlobOptMin::isDone() -{ - return myDone; -} - //======================================================================= //function : Points //purpose : @@ -693,15 +660,4 @@ void math_GlobOptMin::ComputeInitSol() 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 4e65293a49..e98681181e 100644 --- a/src/math/math_GlobOptMin.hxx +++ b/src/math/math_GlobOptMin.hxx @@ -23,87 +23,42 @@ #include #include -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 &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.
-//! Article: Yu. Evtushenko. Numerical methods for finding global extreme (case of a non-uniform mesh).
+//! This class represents Evtushenko's algorithm of global optimization based on non-uniform mesh. +//! Article: Yu. Evtushenko. Numerical methods for finding global extreme (case of a non-uniform mesh). //! U.S.S.R. Comput. Maths. Math. Phys., Vol. 11, N 6, pp. 38-54. - +//! +//! This method performs search on non-uniform mesh. The search space is a box in R^n space. +//! The default behavior is to find all minimums in that box. Computation of maximums is not supported. +//! +//! The search box can be split into smaller boxes by discontinuity criteria. +//! This functionality is covered by SetGlobalParams and SetLocalParams API. +//! +//! It is possible to set continuity of the local boxes. +//! Such option can forcibly change local extrema search. +//! In other words if theFunc can be casted to the function with Hessian but, continuity is set to 1 +//! Gradient based local optimization method will be used, not Hessian based method. +//! This functionality is covered by SetContinuity and GetContinuity API. +//! +//! It is possible to freeze Lipschitz const to avoid internal modifications on it. +//! This functionality is covered by SetLipConstState and GetLipConstState API. +//! +//! It is possible to perform single solution search. +//! This functionality is covered by first parameter in Perform method. +//! +//! It is possible to set / get minimal value of the functional. +//! It works well together with single solution search. +//! This functionality is covered by SetFunctionalMinimalValue and GetFunctionalMinimalValue API. class math_GlobOptMin { public: + //! Constructor. Perform method is not called from it. + //! @param theFunc - objective functional. + //! @param theLowerBorder - lower corner of the search box. + //! @param theUpperBorder - upper corner of the search box. + //! @param theC - Lipschitz constant. + //! @param theDiscretizationTol - parameter space discretization tolerance. + //! @param theSameTol - functional value space indifference tolerance. Standard_EXPORT math_GlobOptMin(math_MultipleVarFunction* theFunc, const math_Vector& theLowerBorder, const math_Vector& theUpperBorder, @@ -111,6 +66,12 @@ public: const Standard_Real theDiscretizationTol = 1.0e-2, const Standard_Real theSameTol = 1.0e-7); + //! @param theFunc - objective functional. + //! @param theLowerBorder - lower corner of the search box. + //! @param theUpperBorder - upper corner of the search box. + //! @param theC - Lipschitz constant. + //! @param theDiscretizationTol - parameter space discretization tolerance. + //! @param theSameTol - functional value space indifference tolerance. Standard_EXPORT void SetGlobalParams(math_MultipleVarFunction* theFunc, const math_Vector& theLowerBorder, const math_Vector& theUpperBorder, @@ -118,47 +79,129 @@ public: const Standard_Real theDiscretizationTol = 1.0e-2, const Standard_Real theSameTol = 1.0e-7); + //! Method to reduce bounding box. Perform will use this box. + //! @param theLocalA - lower corner of the local box. + //! @param theLocalB - upper corner of the local box. Standard_EXPORT void SetLocalParams(const math_Vector& theLocalA, const math_Vector& theLocalB); + //! Method to set tolerances. + //! @param theDiscretizationTol - parameter space discretization tolerance. + //! @param theSameTol - functional value space indifference tolerance. Standard_EXPORT void SetTol(const Standard_Real theDiscretizationTol, const Standard_Real theSameTol); + //! Method to get tolerances. + //! @param theDiscretizationTol - parameter space discretization tolerance. + //! @param theSameTol - functional value space indifference tolerance. Standard_EXPORT void GetTol(Standard_Real& theDiscretizationTol, Standard_Real& theSameTol); - Standard_EXPORT ~math_GlobOptMin(); - //! @param isFindSingleSolution - defines whether to find single solution or all solutions. Standard_EXPORT void Perform(const Standard_Boolean isFindSingleSolution = Standard_False); + //! Return solution theIndex, 1 <= theIndex <= NbExtrema. + Standard_EXPORT void Points(const Standard_Integer theIndex, math_Vector& theSol); + + //! Set / Get continuity of local borders splits (0 ~ C0, 1 ~ C1, 2 ~ C2). + inline void SetContinuity(const Standard_Integer theCont) { myCont = theCont; } + inline Standard_Integer GetContinuity() const {return myCont; } + + //! Set / Get functional minimal value. + inline void SetFunctionalMinimalValue(const Standard_Real theMinimalValue) + { myFunctionalMinimalValue = theMinimalValue; } + inline Standard_Real GetFunctionalMinimalValue() const {return myFunctionalMinimalValue;} + + //! Set / Get Lipchitz constant modification state. + //! True means that the constant is locked and unlocked otherwise. + inline void SetLipConstState(const Standard_Boolean theFlag) {myIsConstLocked = theFlag; } + inline Standard_Boolean GetLipConstState() const { return myIsConstLocked; } + + //! Return computation state of the algorithm. + inline Standard_Boolean isDone() const {return myDone; } + //! Get best functional value. - Standard_EXPORT Standard_Real GetF(); + inline Standard_Real GetF() const {return myF;} //! Return count of global extremas. - Standard_EXPORT Standard_Integer NbExtrema(); + inline Standard_Integer NbExtrema() const {return mySolCount;} - //! Return solution theIndex, 1 <= theIndex <= NbExtrema. - Standard_EXPORT void Points(const Standard_Integer theIndex, math_Vector& theSol); +private: - Standard_EXPORT Standard_Boolean isDone(); + //! Class for duplicate fast search. For internal usage only. + class NCollection_CellFilter_Inspector + { + public: - //! Set functional minimal value. - Standard_EXPORT void SetFunctionalMinimalValue(const Standard_Real theMinimalValue); + //! Points and target type + typedef math_Vector Point; + typedef math_Vector Target; - //! Lock/Unlock Lipchitz constant for internal modifications. - Standard_EXPORT void SetLipConstState(const Standard_Boolean theFlag); + NCollection_CellFilter_Inspector(const Standard_Integer theDim, + const Standard_Real theTol) + : myCurrent(1, theDim) + { + myTol = theTol * theTol; + myIsFind = Standard_False; + Dimension = theDim; + } - //! Get functional minimal value. - Standard_EXPORT Standard_Real GetFunctionalMinimalValue(); + //! Access to coordinate. + static Standard_Real Coord (int i, const Point &thePnt) + { + return thePnt(i + 1); + } - //! Get continuity of local borders splits. - inline Standard_Integer GetContinuity() const {return myCont; } + //! 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 &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); + } + } - //! Set continuity of local borders splits. - inline void SetContinuity(const Standard_Integer theCont) { myCont = theCont; } + 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; + }; -private: // Compute cell size. void initCellSize(); @@ -173,14 +216,14 @@ private: void computeGlobalExtremum(Standard_Integer theIndex); //! Check possibility to stop computations. - //! Find single solution + in neighbourhood of best possible solution. + //! Find single solution + in neighborhood of best possible solution. Standard_Boolean CheckFunctionalStopCriteria(); //! Computes starting value / approximation: //! myF - initial best value. //! myY - initial best point. //! myC - approximation of Lipschitz constant. - //! to imporve convergence speed. + //! to improve convergence speed. void computeInitialValues(); //! Check that myA <= thePnt <= myB @@ -213,7 +256,7 @@ private: // Algorithm data. Standard_Real myZ; - Standard_Real myE1; // Border coeff. + Standard_Real myE1; // Border coefficient. Standard_Real myE2; // Minimum step size. Standard_Real myE3; // Local extrema starting parameter. @@ -222,7 +265,6 @@ private: 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 myCellSize; Standard_Integer myMinCellFilterSol; diff --git a/tests/bugs/fclasses/bug25635_1 b/tests/bugs/fclasses/bug25635_1 index 2ae860f420..0770219c53 100755 --- a/tests/bugs/fclasses/bug25635_1 +++ b/tests/bugs/fclasses/bug25635_1 @@ -11,7 +11,7 @@ ellipse c2 4 0 2 1 set info [2dextrema c1 c2] -set tol_abs 1.e-5 +set tol_abs 7.e-5 set tol_rel 0.01 #-1 @@ -36,7 +36,7 @@ checkreal "Parameter2" ${Parameter2} ${expected_Parameter2} ${tol_abs} ${tol_rel set expected_OriginX 2. checkreal "OriginX" ${OriginX} ${expected_OriginX} ${tol_abs} ${tol_rel} -set expected_OriginY 7.e-05 +set expected_OriginY 0.0 checkreal "OriginY" ${OriginY} ${expected_OriginY} ${tol_abs} ${tol_rel} set expected_AxisX 1. diff --git a/tests/bugs/fclasses/bug27371 b/tests/bugs/fclasses/bug27371 new file mode 100644 index 0000000000..9b9de46353 --- /dev/null +++ b/tests/bugs/fclasses/bug27371 @@ -0,0 +1,33 @@ +puts "========" +puts "OCC27371" +puts "========" +puts "" +############################################## +# Regression: BRepExtrema works too much slower in 691 (from 670) +############################################## +restore [locate_data_file bug27371.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 + distmini d aShape_2 aShape_1 +} +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.2 +set aDist [dval d_val] +checkreal "Distance value check" $aDist $aDist_Exp $absTol $relTol \ No newline at end of file