From: aml Date: Tue, 28 Jul 2015 09:18:04 +0000 (+0300) Subject: 0027131: [Regression to 6.7] DistShapeShape performance loss X-Git-Url: http://git.dev.opencascade.org/gitweb/?a=commitdiff_plain;h=38c459922e9ee1a56f3b794add1c97bfdd33b28f;p=occt-copy.git 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. --- diff --git a/src/BRepMesh/BRepMesh_CircleTool.hxx b/src/BRepMesh/BRepMesh_CircleTool.hxx index c74a2a4fe6..a81e739b63 100644 --- a/src/BRepMesh/BRepMesh_CircleTool.hxx +++ b/src/BRepMesh/BRepMesh_CircleTool.hxx @@ -25,6 +25,7 @@ #include #include #include +#include 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 aCellSize(aCellSizeC[0], 1, 2); myCellFilter.Reset(aCellSize, myAllocator); } diff --git a/src/BRepMesh/BRepMesh_VertexTool.hxx b/src/BRepMesh/BRepMesh_VertexTool.hxx index 3fe895b3c8..1c632e88d2 100644 --- a/src/BRepMesh/BRepMesh_VertexTool.hxx +++ b/src/BRepMesh/BRepMesh_VertexTool.hxx @@ -14,6 +14,7 @@ #ifndef _BRepMesh_VertexTool_HeaderFile #define _BRepMesh_VertexTool_HeaderFile +#include #include #include #include @@ -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 aCellSize(aCellSizeC[0], 1, 2); myCellFilter.Reset(aCellSize, myAllocator); mySelector.Clear(); } diff --git a/src/Extrema/Extrema_GenExtCC.gxx b/src/Extrema/Extrema_GenExtCC.gxx index 6acba89bd9..a342a61b9a 100644 --- a/src/Extrema/Extrema_GenExtCC.gxx +++ b/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; diff --git a/src/NCollection/NCollection_CellFilter.hxx b/src/NCollection/NCollection_CellFilter.hxx index 2d49f117b6..9c6dfac4d1 100644 --- a/src/NCollection/NCollection_CellFilter.hxx +++ b/src/NCollection/NCollection_CellFilter.hxx @@ -17,6 +17,8 @@ #define NCollection_CellFilter_HeaderFile #include +#include +#include #include #include #include @@ -110,37 +112,39 @@ enum NCollection_CellFilter_Action * Note that method Inspect() can be const and/or virtual. */ -template -class NCollection_CellFilter -{ +template class NCollection_CellFilter +{ public: typedef TYPENAME Inspector::Target Target; typedef TYPENAME Inspector::Point Point; 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 + //! Note: the cell size must be ensured to be greater than //! maximal co-ordinate of the involved points divided by INT_MAX, //! in order to avoid integer overflow of cell index. //! //! 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, - const Handle(NCollection_IncAllocator)& theAlloc=0) + //! Constructor when dimension count is unknown at compilation time. + 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. - //! Note: the cell size in each dimension must be ensured to be greater than - //! maximal co-ordinate of the involved points by this dimension divided by INT_MAX, - //! in order to avoid integer overflow of cell index. - NCollection_CellFilter (Standard_Real theCellSize[], - const Handle(NCollection_IncAllocator)& theAlloc=0) + //! Constructor when dimenstion count is known at compilation time. + NCollection_CellFilter (const Standard_Real theCellSize = 0, + const Handle(NCollection_IncAllocator)& theAlloc = 0) + : myCellSize(0, Inspector::Dimension - 1) { + 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++) - myCellSize[i] = theCellSize; + for (int i=0; i < myDim; i++) + 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 theCellSize, const Handle(NCollection_IncAllocator)& theAlloc=0) { - for (int i=0; i < Inspector::Dimension; i++) - myCellSize[i] = theCellSize[i]; + myCellSize = theCellSize; 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[]) - : Objects(0) + Cell (const Point& thePnt, + const NCollection_Array1& 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++) - index[i] = theOther.index[i]; + Standard_Integer myDim = Standard_Integer(theOther.index.Size()); + 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 index; ListNode *Objects; }; @@ -452,9 +469,10 @@ protected: } protected: + Standard_Integer myDim; Handle(NCollection_BaseAllocator) myAllocator; NCollection_Map myCells; - Standard_Real myCellSize [Inspector::Dimension]; + NCollection_Array1 myCellSize; }; /** @@ -504,4 +522,3 @@ struct NCollection_CellFilter_InspectorXY }; #endif - diff --git a/src/NCollection/NCollection_LocalArray.hxx b/src/NCollection/NCollection_LocalArray.hxx index ae2bbf91ca..20580cc31c 100644 --- a/src/NCollection/NCollection_LocalArray.hxx +++ b/src/NCollection/NCollection_LocalArray.hxx @@ -20,13 +20,10 @@ //! Auxiliary class optimizing creation of array buffer //! (using stack allocation for small arrays). -template class NCollection_LocalArray +template 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; }; diff --git a/src/math/math_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx index ef8b51ea88..b8f0869c6b 100644 --- a/src/math/math_GlobOptMin.cxx +++ b/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(myFunc)) { - math_MultipleVarFunctionWithHessian* myTmp = + math_MultipleVarFunctionWithHessian* aTmp = dynamic_cast (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(myFunc)) { - math_MultipleVarFunctionWithGradient* myTmp = + math_MultipleVarFunctionWithGradient* aTmp = dynamic_cast (myFunc); - math_BFGS bfgs(myTmp->NbVariables()); - bfgs.Perform(*myTmp, thePnt); + math_BFGS bfgs(aTmp->NbVariables()); + 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 - // 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); @@ -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,7 +399,8 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j) Standard_Real r; Standard_Boolean isReached = Standard_False; - for(myX(j) = myA(j) + myE1; + + 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) + { + for(i = 0; i < mySolCount; i++) + { + isSame = Standard_True; + for(j = 1; j <= myN; j++) + { + if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j)) + { + isSame = Standard_False; + break; + } + } + if (isSame == Standard_True) + return Standard_True; + } + } + else { - isSame = Standard_True; - for(j = 1; j <= myN; j++) + NCollection_CellFilter_Inspector anInspector(myN, Precision::PConfusion()); + if (isFirstCellFilterInvoke) { - if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j)) + myFilter.Reset(myCellSize); + + // Copy initial data into cell filter. + for(Standard_Integer aSolIdx = 0; aSolIdx < mySolCount; aSolIdx++) { - isSame = Standard_False; - break; + 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; +} diff --git a/src/math/math_GlobOptMin.hxx b/src/math/math_GlobOptMin.hxx index 40a6d49445..5bf006ac05 100644 --- a/src/math/math_GlobOptMin.hxx +++ b/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 +#include +#include #include #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).
//! 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 myCellSize; + Standard_Integer myMinCellFilterSol; + Standard_Boolean isFirstCellFilterInvoke; + NCollection_CellFilter myFilter; + Standard_Real myF; // Current value of Global optimum. }; 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