diff --git a/src/Extrema/Extrema_GenExtCC.gxx b/src/Extrema/Extrema_GenExtCC.gxx index 88ca463c7c..3b564f4c74 100644 --- 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 2560aa252d..bf75a6af71 100644 --- 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 1d42ab91e3..4c4e52096f 100644 --- 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 myCellSize; diff --git a/tests/bugs/fclasses/bug27131 b/tests/bugs/fclasses/bug27131 new file mode 100644 index 0000000000..7f833ff7bf --- /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 516eaee789..8cd60bb33e 100755 --- 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 }