diff --git a/src/Extrema/Extrema_GenExtCS.cdl b/src/Extrema/Extrema_GenExtCS.cdl index ac855cdbc6..f36de818bb 100644 --- a/src/Extrema/Extrema_GenExtCS.cdl +++ b/src/Extrema/Extrema_GenExtCS.cdl @@ -22,7 +22,7 @@ class GenExtCS from Extrema uses POnSurf from Extrema, POnCurv from Extrema, - Pnt from gp, + Pnt from gp, HArray1OfPnt from TColgp, HArray2OfPnt from TColgp, FuncExtCS from Extrema, @@ -139,17 +139,15 @@ is fields myDone : Boolean; myInit : Boolean; - mytmin : Real; - mytsup : Real; - myumin : Real; - myusup : Real; - myvmin : Real; - myvsup : Real; + mytmin : Real; + mytsup : Real; + myumin : Real; + myusup : Real; + myvmin : Real; + myvsup : Real; mytsample : Integer; myusample : Integer; myvsample : Integer; - mypoints1 : HArray1OfPnt from TColgp; - mypoints2 : HArray2OfPnt from TColgp; mytol1 : Real; mytol2 : Real; myF : FuncExtCS from Extrema; diff --git a/src/Extrema/Extrema_GenExtCS.cxx b/src/Extrema/Extrema_GenExtCS.cxx index 812c15aa16..ce61b3182b 100644 --- a/src/Extrema/Extrema_GenExtCS.cxx +++ b/src/Extrema/Extrema_GenExtCS.cxx @@ -16,27 +16,27 @@ #include -#include -#include -#include - -#include #include -#include - #include -#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include -#include -#include -#include +const Standard_Real aMaxParamVal = 1.0e+10; +const Standard_Real aBorderDivisor = 1.0e+4; //======================================================================= //function : Extrema_GenExtCS //purpose : //======================================================================= - Extrema_GenExtCS::Extrema_GenExtCS() { myDone = Standard_False; @@ -47,14 +47,13 @@ Extrema_GenExtCS::Extrema_GenExtCS() //function : Extrema_GenExtCS //purpose : //======================================================================= - -Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C, - const Adaptor3d_Surface& S, - const Standard_Integer NbT, - const Standard_Integer NbU, - const Standard_Integer NbV, - const Standard_Real Tol1, - const Standard_Real Tol2) +Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C, + const Adaptor3d_Surface& S, + const Standard_Integer NbT, + const Standard_Integer NbU, + const Standard_Integer NbV, + const Standard_Real Tol1, + const Standard_Real Tol2) { Initialize(S, NbU, NbV, Tol2); Perform(C, NbT, Tol1); @@ -65,18 +64,18 @@ Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C, //purpose : //======================================================================= -Extrema_GenExtCS::Extrema_GenExtCS (const Adaptor3d_Curve& C, - const Adaptor3d_Surface& S, - const Standard_Integer NbT, - const Standard_Integer NbU, - const Standard_Integer NbV, - const Standard_Real tmin, - const Standard_Real tsup, - const Standard_Real Umin, +Extrema_GenExtCS::Extrema_GenExtCS (const Adaptor3d_Curve& C, + const Adaptor3d_Surface& S, + const Standard_Integer NbT, + const Standard_Integer NbU, + const Standard_Integer NbV, + const Standard_Real tmin, + const Standard_Real tsup, + const Standard_Real Umin, const Standard_Real Usup, - const Standard_Real Vmin, - const Standard_Real Vsup, - const Standard_Real Tol1, + const Standard_Real Vmin, + const Standard_Real Vsup, + const Standard_Real Tol1, const Standard_Real Tol2) { Initialize(S, NbU, NbV, Umin,Usup,Vmin,Vsup,Tol2); @@ -87,10 +86,9 @@ Extrema_GenExtCS::Extrema_GenExtCS (const Adaptor3d_Curve& C, //function : Initialize //purpose : //======================================================================= - -void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, - const Standard_Integer NbU, - const Standard_Integer NbV, +void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, + const Standard_Integer NbU, + const Standard_Integer NbV, const Standard_Real Tol2) { myumin = S.FirstUParameter(); @@ -104,7 +102,6 @@ void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, //function : Initialize //purpose : //======================================================================= - void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, const Standard_Integer NbU, const Standard_Integer NbV, @@ -123,15 +120,15 @@ void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, myvsup = Vsup; mytol2 = Tol2; - const Standard_Real aTrimMaxU = Precision::IsInfinite (myusup) ? 1.0e+10 : myusup; - const Standard_Real aTrimMinU = Precision::IsInfinite (myumin) ? -1.0e+10 : myumin; - const Standard_Real aTrimMaxV = Precision::IsInfinite (myvsup) ? 1.0e+10 : myvsup; - const Standard_Real aTrimMinV = Precision::IsInfinite (myvmin) ? -1.0e+10 : myvmin; + const Standard_Real aTrimMaxU = Precision::IsInfinite (myusup) ? aMaxParamVal : myusup; + const Standard_Real aTrimMinU = Precision::IsInfinite (myumin) ? -aMaxParamVal : myumin; + const Standard_Real aTrimMaxV = Precision::IsInfinite (myvsup) ? aMaxParamVal : myvsup; + const Standard_Real aTrimMinV = Precision::IsInfinite (myvmin) ? -aMaxParamVal : myvmin; - const Standard_Real aMinU = aTrimMinU + (aTrimMaxU - aTrimMinU) / 1.0e4; - const Standard_Real aMinV = aTrimMinV + (aTrimMaxV - aTrimMinV) / 1.0e4; - const Standard_Real aMaxU = aTrimMaxU - (aTrimMaxU - aTrimMinU) / 1.0e4; - const Standard_Real aMaxV = aTrimMaxV - (aTrimMaxV - aTrimMinV) / 1.0e4; + const Standard_Real aMinU = aTrimMinU + (aTrimMaxU - aTrimMinU) / aBorderDivisor; + const Standard_Real aMinV = aTrimMinV + (aTrimMaxV - aTrimMinV) / aBorderDivisor; + const Standard_Real aMaxU = aTrimMaxU - (aTrimMaxU - aTrimMinU) / aBorderDivisor; + const Standard_Real aMaxV = aTrimMaxV - (aTrimMaxV - aTrimMinV) / aBorderDivisor; const Standard_Real aStepSU = (aMaxU - aMinU) / myusample; const Standard_Real aStepSV = (aMaxV - aMinV) / myvsample; @@ -153,7 +150,6 @@ void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, //function : Perform //purpose : //======================================================================= - void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C, const Standard_Integer NbT, const Standard_Real Tol1) @@ -163,86 +159,10 @@ void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C, Perform(C, NbT, mytmin, mytsup,Tol1); } -namespace -{ - //! Vector type for storing curve and surface parameters. - typedef NCollection_Vec3 Extrema_Vec3d; - - //! Describes particle for using in PSO algorithm. - struct Extrema_Particle - { - //! Creates unitialized particle. - Extrema_Particle() - : Distance (DBL_MAX) - { - // - } - - //! Creates particle with the given position and distance. - Extrema_Particle (const Extrema_Vec3d& thePosition, Standard_Real theDistance) - : Position (thePosition), - Distance (theDistance), - BestPosition (thePosition), - BestDistance (theDistance) - { - // - } - - //! Compares the particles according to their distances. - bool operator< (const Extrema_Particle& thePnt) const - { - return Distance < thePnt.Distance; - } - - Extrema_Vec3d Position; - Extrema_Vec3d Velocity; - Standard_Real Distance; - - Extrema_Vec3d BestPosition; - Standard_Real BestDistance; - }; - - //! Fast random number generator (the algorithm proposed by Ian C. Bullard). - class Extrema_Random - { - private: - - unsigned int myStateHi; - unsigned int myStateLo; - - public: - - //! Creates new Xorshift 64-bit RNG. - Extrema_Random (unsigned int theSeed = 1) - : myStateHi (theSeed) - { - myStateLo = theSeed ^ 0x49616E42; - } - - //! Generates new 64-bit integer value. - unsigned int NextInt() - { - myStateHi = (myStateHi >> 2) + (myStateHi << 2); - - myStateHi += myStateLo; - myStateLo += myStateHi; - - return myStateHi; - } - - //! Generates new floating-point value. - Standard_Real NextReal() - { - return NextInt() / static_cast (0xFFFFFFFFu); - } - }; -} - //======================================================================= //function : Perform //purpose : //======================================================================= - void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C, const Standard_Integer NbT, const Standard_Real tmin, @@ -259,30 +179,30 @@ void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C, Standard_Real trimusup = myusup, trimumin = myumin,trimvsup = myvsup,trimvmin = myvmin; if (Precision::IsInfinite(trimusup)){ - trimusup = 1.0e+10; + trimusup = aMaxParamVal; } if (Precision::IsInfinite(trimvsup)){ - trimvsup = 1.0e+10; + trimvsup = aMaxParamVal; } if (Precision::IsInfinite(trimumin)){ - trimumin = -1.0e+10; + trimumin = -aMaxParamVal; } if (Precision::IsInfinite(trimvmin)){ - trimvmin = -1.0e+10; + trimvmin = -aMaxParamVal; } // math_Vector Tol(1, 3); Tol(1) = mytol1; Tol(2) = mytol2; Tol(3) = mytol2; - math_Vector UV(1,3), UVinf(1,3), UVsup(1,3); - UVinf(1) = mytmin; - UVinf(2) = trimumin; - UVinf(3) = trimvmin; + math_Vector TUV(1,3), TUVinf(1,3), TUVsup(1,3); + TUVinf(1) = mytmin; + TUVinf(2) = trimumin; + TUVinf(3) = trimvmin; - UVsup(1) = mytsup; - UVsup(2) = trimusup; - UVsup(3) = trimvsup; + TUVsup(1) = mytsup; + TUVsup(2) = trimusup; + TUVsup(3) = trimvsup; // 18/02/02 akm vvv : (OCC163) bad extremas - extrusion surface versus the line. // Try to preset the initial solution as extrema between // extrusion direction and the curve. @@ -304,7 +224,7 @@ void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C, { aLocator.Points (iExt, aP1, aP2); // Parameter on curve - UV(1) = aP1.Parameter(); + TUV(1) = aP1.Parameter(); // To find parameters on surf, try ExtPS Extrema_ExtPS aPreciser (aP1.Value(), *myS, mytol2, mytol2); if (aPreciser.IsDone()) @@ -313,43 +233,39 @@ void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C, Standard_Integer iPExt; for (iPExt=1; iPExt<=aPreciser.NbExt(); iPExt++) { - aPreciser.Point(iPExt).Parameter(UV(2),UV(3)); - math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup); + aPreciser.Point(iPExt).Parameter(TUV(2),TUV(3)); + math_FunctionSetRoot S1 (myF,TUV,Tol,TUVinf,TUVsup); } } else { // Failed... try the point on iso line - UV(2) = dfUFirst; - UV(3) = aP2.Parameter(); - math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup); + TUV(2) = dfUFirst; + TUV(3) = aP2.Parameter(); + math_FunctionSetRoot S1 (myF,TUV,Tol,TUVinf,TUVsup); } } // for (iExt=1; iExt<=aLocator.NbExt(); iExt++) } // if (aLocator.IsDone() && aLocator.NbExt()>0) } // if (myS.Type() == GeomAbs_ExtrusionSurface) else { - // Number of particles used in PSO algorithm (particle swarm optimization) + // Number of particles used in PSO algorithm (particle swarm optimization). const Standard_Integer aNbParticles = 32; - NCollection_Array1 aParticles (1, aNbParticles); + math_PSOParticlesPool aParticles(aNbParticles, 3); - const Extrema_Vec3d aMinUV (UVinf(1) + (UVsup(1) - UVinf(1)) / 1.0e4, - UVinf(2) + (UVsup(2) - UVinf(2)) / 1.0e4, - UVinf(3) + (UVsup(3) - UVinf(3)) / 1.0e4); + math_Vector aMinTUV(1,3); + aMinTUV = TUVinf + (TUVsup - TUVinf) / aBorderDivisor; - const Extrema_Vec3d aMaxUV (UVsup(1) - (UVsup(1) - UVinf(1)) / 1.0e4, - UVsup(2) - (UVsup(2) - UVinf(2)) / 1.0e4, - UVsup(3) - (UVsup(3) - UVinf(3)) / 1.0e4); + math_Vector aMaxTUV(1,3); + aMaxTUV = TUVsup - (TUVsup - TUVinf) / aBorderDivisor; - Standard_Real aStepCU = (aMaxUV.x() - aMinUV.x()) / mytsample; - Standard_Real aStepSU = (aMaxUV.y() - aMinUV.y()) / myusample; - Standard_Real aStepSV = (aMaxUV.z() - aMinUV.z()) / myvsample; + Standard_Real aStepCU = (aMaxTUV(1) - aMinTUV(1)) / mytsample; + Standard_Real aStepSU = (aMaxTUV(2) - aMinTUV(2)) / myusample; + Standard_Real aStepSV = (aMaxTUV(3) - aMinTUV(3)) / myvsample; // Correct number of curve samples in case of low resolution - Standard_Real aScaleFactor = 5.0; - Standard_Real aResolutionCU = aStepCU / C.Resolution (1.0); Standard_Real aMinResolution = aScaleFactor * Min (aResolutionCU, @@ -361,174 +277,64 @@ void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C, { const Standard_Integer aMaxNbNodes = 50; - mytsample = Min (aMaxNbNodes, RealToInt ( - mytsample * aResolutionCU / aMinResolution)); + mytsample = Min(aMaxNbNodes, + RealToInt(mytsample * aResolutionCU / aMinResolution)); - aStepCU = (aMaxUV.x() - aMinUV.x()) / mytsample; + aStepCU = (aMaxTUV(1) - aMinTUV(1)) / mytsample; } } - // Pre-compute curve sample points + // Pre-compute curve sample points. TColgp_HArray1OfPnt aCurvPnts (0, mytsample); - Standard_Real aCU = aMinUV.x(); + Standard_Real aCU = aMinTUV(1); for (Standard_Integer aCUI = 0; aCUI <= mytsample; aCUI++, aCU += aStepCU) - { aCurvPnts.SetValue (aCUI, C.Value (aCU)); - } - - NCollection_Array1::iterator aWorstPnt = aParticles.begin(); + PSO_Particle* aParticle = aParticles.GetWorstParticle(); // Select specified number of particles from pre-computed set of samples - Standard_Real aSU = aMinUV.y(); + Standard_Real aSU = aMinTUV(2); for (Standard_Integer aSUI = 0; aSUI <= myusample; aSUI++, aSU += aStepSU) { - Standard_Real aSV = aMinUV.z(); + Standard_Real aSV = aMinTUV(3); for (Standard_Integer aSVI = 0; aSVI <= myvsample; aSVI++, aSV += aStepSV) { - Standard_Real aCU = aMinUV.x(); + Standard_Real aCU = aMinTUV(1); for (Standard_Integer aCUI = 0; aCUI <= mytsample; aCUI++, aCU += aStepCU) { - Standard_Real aSqDist = mySurfPnts->Value (aSUI, aSVI).SquareDistance (aCurvPnts.Value (aCUI)); + Standard_Real aSqDist = mySurfPnts->Value(aSUI, aSVI).SquareDistance(aCurvPnts.Value(aCUI)); - if (aSqDist < aWorstPnt->Distance) + if (aSqDist < aParticle->Distance) { - *aWorstPnt = Extrema_Particle (Extrema_Vec3d (aCU, aSU, aSV), aSqDist); + aParticle->Position[0] = aCU; + aParticle->Position[1] = aSU; + aParticle->Position[2] = aSV; - aWorstPnt = std::max_element (aParticles.begin(), aParticles.end()); + aParticle->BestPosition[0] = aCU; + aParticle->BestPosition[1] = aSU; + aParticle->BestPosition[2] = aSV; + + aParticle->Distance = aSqDist; + aParticle->BestDistance = aSqDist; + + aParticle = aParticles.GetWorstParticle(); } } } } - Extrema_Random aRandom; - - // Generate initial particle velocities - for (Standard_Integer aPartIdx = 1; aPartIdx <= aNbParticles; ++aPartIdx) - { - const Standard_Real aKsi1 = aRandom.NextReal(); - const Standard_Real aKsi2 = aRandom.NextReal(); - const Standard_Real aKsi3 = aRandom.NextReal(); - - aParticles.ChangeValue (aPartIdx).Velocity.x() = aStepCU * (aKsi1 - 0.5) * 2.0; - aParticles.ChangeValue (aPartIdx).Velocity.y() = aStepSU * (aKsi2 - 0.5) * 2.0; - aParticles.ChangeValue (aPartIdx).Velocity.z() = aStepSV * (aKsi3 - 0.5) * 2.0; - } - - NCollection_Array1::iterator aBestPnt = - std::min_element (aParticles.begin(), aParticles.end()); - - Extrema_Vec3d aBestGlobalPosition = aBestPnt->Position; - Standard_Real aBestGlobalDistance = aBestPnt->Distance; - - // This velocity is used for detecting stagnation state - Extrema_Vec3d aTerminationVelocity (aStepCU / 2048, aStepSU / 2048, aStepSV / 2048); - - // Run PSO iterations - for (Standard_Integer aStep = 1, aMaxNbSteps = 100; aStep < aMaxNbSteps; ++aStep) - { - Extrema_Vec3d aMinimalVelocity (DBL_MAX, DBL_MAX, DBL_MAX); - - for (Standard_Integer aPartIdx = 1; aPartIdx <= aParticles.Size(); ++aPartIdx) - { - const Standard_Real aKsi1 = aRandom.NextReal(); - const Standard_Real aKsi2 = aRandom.NextReal(); - - Extrema_Particle& aPoint = aParticles.ChangeValue (aPartIdx); - - const Standard_Real aRetentWeight = 0.72900; - const Standard_Real aPersonWeight = 1.49445; - const Standard_Real aSocialWeight = 1.49445; - - aPoint.Velocity = aPoint.Velocity * aRetentWeight + - (aPoint.BestPosition - aPoint.Position) * (aPersonWeight * aKsi1) + - (aBestGlobalPosition - aPoint.Position) * (aSocialWeight * aKsi2); - - aPoint.Position += aPoint.Velocity; - - aPoint.Position.x() = Min (Max (aPoint.Position.x(), aMinUV.x()), aMaxUV.x()); - aPoint.Position.y() = Min (Max (aPoint.Position.y(), aMinUV.y()), aMaxUV.y()); - aPoint.Position.z() = Min (Max (aPoint.Position.z(), aMinUV.z()), aMaxUV.z()); - - aPoint.Distance = myS->Value (aPoint.Position.y(), - aPoint.Position.z()).SquareDistance (C.Value (aPoint.Position.x())); - - if (aPoint.Distance < aPoint.BestDistance) - { - aPoint.BestDistance = aPoint.Distance; - aPoint.BestPosition = aPoint.Position; - - if (aPoint.Distance < aBestGlobalDistance) - { - aBestGlobalDistance = aPoint.Distance; - aBestGlobalPosition = aPoint.Position; - } - } - - aMinimalVelocity.x() = Min (Abs (aPoint.Velocity.x()), aMinimalVelocity.x()); - aMinimalVelocity.y() = Min (Abs (aPoint.Velocity.y()), aMinimalVelocity.y()); - aMinimalVelocity.z() = Min (Abs (aPoint.Velocity.x()), aMinimalVelocity.z()); - } - - if (aMinimalVelocity.x() < aTerminationVelocity.x() - && aMinimalVelocity.y() < aTerminationVelocity.y() - && aMinimalVelocity.z() < aTerminationVelocity.z()) - { - // Minimum number of steps - const Standard_Integer aMinSteps = 16; - - if (aStep > aMinSteps) - { - break; - } - - for (Standard_Integer aPartIdx = 1; aPartIdx <= aNbParticles; ++aPartIdx) - { - const Standard_Real aKsi1 = aRandom.NextReal(); - const Standard_Real aKsi2 = aRandom.NextReal(); - const Standard_Real aKsi3 = aRandom.NextReal(); - - Extrema_Particle& aPoint = aParticles.ChangeValue (aPartIdx); - - if (aPoint.Position.x() == aMinUV.x() || aPoint.Position.x() == aMaxUV.x()) - { - aPoint.Velocity.x() = aPoint.Position.x() == aMinUV.x() ? - aStepCU * aKsi1 : -aStepCU * aKsi1; - } - else - { - aPoint.Velocity.x() = aStepCU * (aKsi1 - 0.5) * 2.0; - } - - if (aPoint.Position.y() == aMinUV.y() || aPoint.Position.y() == aMaxUV.y()) - { - aPoint.Velocity.y() = aPoint.Position.y() == aMinUV.y() ? - aStepSU * aKsi2 : -aStepSU * aKsi2; - } - else - { - aPoint.Velocity.y() = aStepSU * (aKsi2 - 0.5) * 2.0; - } - - if (aPoint.Position.z() == aMinUV.z() || aPoint.Position.z() == aMaxUV.z()) - { - aPoint.Velocity.z() = aPoint.Position.z() == aMinUV.z() ? - aStepSV * aKsi3 : -aStepSV * aKsi3; - } - else - { - aPoint.Velocity.z() = aStepSV * (aKsi3 - 0.5) * 2.0; - } - } - } - } + math_Vector aStep(1,3); + aStep(1) = aStepCU; + aStep(2) = aStepSU; + aStep(3) = aStepSV; // Find min approximation - UV(1) = aBestGlobalPosition.x(); - UV(2) = aBestGlobalPosition.y(); - UV(3) = aBestGlobalPosition.z(); - - math_FunctionSetRoot anA (myF, UV, Tol, UVinf, UVsup, 100, Standard_False); + Standard_Real aValue; + Extrema_GlobOptFuncCS aFunc(&C, myS); + math_PSO aPSO(&aFunc, TUVinf, TUVsup, aStep); + aPSO.Perform(aParticles, aNbParticles, aValue, TUV); + + math_FunctionSetRoot anA (myF, TUV, Tol, TUVinf, TUVsup, 100, Standard_False); } myDone = Standard_True; @@ -538,7 +344,6 @@ void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C, //function : IsDone //purpose : //======================================================================= - Standard_Boolean Extrema_GenExtCS::IsDone() const { return myDone; @@ -548,7 +353,6 @@ Standard_Boolean Extrema_GenExtCS::IsDone() const //function : NbExt //purpose : //======================================================================= - Standard_Integer Extrema_GenExtCS::NbExt() const { if (!IsDone()) { StdFail_NotDone::Raise(); } @@ -559,7 +363,6 @@ Standard_Integer Extrema_GenExtCS::NbExt() const //function : Value //purpose : //======================================================================= - Standard_Real Extrema_GenExtCS::SquareDistance(const Standard_Integer N) const { if (!IsDone()) { StdFail_NotDone::Raise(); } @@ -570,7 +373,6 @@ Standard_Real Extrema_GenExtCS::SquareDistance(const Standard_Integer N) const //function : PointOnCurve //purpose : //======================================================================= - const Extrema_POnCurv& Extrema_GenExtCS::PointOnCurve(const Standard_Integer N) const { if (!IsDone()) { StdFail_NotDone::Raise(); } @@ -581,7 +383,6 @@ const Extrema_POnCurv& Extrema_GenExtCS::PointOnCurve(const Standard_Integer N) //function : PointOnSurface //purpose : //======================================================================= - const Extrema_POnSurf& Extrema_GenExtCS::PointOnSurface(const Standard_Integer N) const { if (!IsDone()) { StdFail_NotDone::Raise(); } @@ -592,7 +393,6 @@ const Extrema_POnSurf& Extrema_GenExtCS::PointOnSurface(const Standard_Integer N //function : BidonSurface //purpose : //======================================================================= - Adaptor3d_SurfacePtr Extrema_GenExtCS::BidonSurface() const { return (Adaptor3d_SurfacePtr)0L; @@ -602,7 +402,6 @@ Adaptor3d_SurfacePtr Extrema_GenExtCS::BidonSurface() const //function : BidonCurve //purpose : //======================================================================= - Adaptor3d_CurvePtr Extrema_GenExtCS::BidonCurve() const { return (Adaptor3d_CurvePtr)0L; diff --git a/src/Extrema/Extrema_GlobOptFuncCS.cxx b/src/Extrema/Extrema_GlobOptFuncCS.cxx new file mode 100644 index 0000000000..53982401c4 --- /dev/null +++ b/src/Extrema/Extrema_GlobOptFuncCS.cxx @@ -0,0 +1,227 @@ +// Created on: 2014-06-23 +// Created by: Alexander Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement + +#include + +#include +#include +#include +#include +#include + +//! F = Sqrt( (x1(cu) - x2(su,sv))^2 +// + (y1(cu) - y2(su,sv))^2 +// + (z1(cu) - z2(su,sv))^2 ) + + +//======================================================================= +//function : value +//purpose : +//======================================================================= +void Extrema_GlobOptFuncCS::value(Standard_Real cu, + Standard_Real su, + Standard_Real sv, + Standard_Real &F) +{ + F = myC->Value(cu).Distance(myS->Value(su, sv)); +} + +//======================================================================= +//function : gradient +//purpose : +//======================================================================= +void Extrema_GlobOptFuncCS::gradient(Standard_Real cu, + Standard_Real su, + Standard_Real sv, + math_Vector &G) +{ + gp_Pnt CD0, SD0; + gp_Vec CD1, SD1U, SD1V; + + myC->D1(cu, CD0, CD1); + myS->D1(su, sv, SD0, SD1U, SD1V); + + G(1) = + (CD0.X() - SD0.X()) * CD1.X() + + (CD0.Y() - SD0.Y()) * CD1.Y() + + (CD0.Z() - SD0.Z()) * CD1.Z(); + G(2) = - (CD0.X() - SD0.X()) * SD1U.X() + - (CD0.Y() - SD0.Y()) * SD1U.Y() + - (CD0.Z() - SD0.Z()) * SD1U.Z(); + G(3) = - (CD0.X() - SD0.X()) * SD1V.X() + - (CD0.Y() - SD0.Y()) * SD1V.Y() + - (CD0.Z() - SD0.Z()) * SD1V.Z(); +} + +//======================================================================= +//function : hessian +//purpose : +//======================================================================= +void Extrema_GlobOptFuncCS::hessian(Standard_Real cu, + Standard_Real su, + Standard_Real sv, + math_Matrix &H) +{ + gp_Pnt CD0, SD0; + gp_Vec CD1, SD1U, SD1V, CD2, SD2UU, SD2UV, SD2VV; + + myC->D2(cu, CD0, CD1, CD2); + myS->D2(su, sv, SD0, SD1U, SD1V, SD2UU, SD2VV, SD2UV); + + H(1,1) = + CD1.X() * CD1.X() + + CD1.Y() * CD1.Y() + + CD1.Z() * CD1.Z() + + (CD0.X() - SD0.X()) * CD2.X() + + (CD0.Y() - SD0.Y()) * CD2.Y() + + (CD0.Z() - SD0.Z()) * CD2.Z(); + H(1,2) = - CD1.X() * SD1U.X() + - CD1.Y() * SD1U.Y() + - CD1.Z() * SD1U.Z(); + H(1,3) = - CD1.X() * SD1V.X() + - CD1.Y() * SD1V.Y() + - CD1.Z() * SD1V.Z(); + H(2,1) = H(1,2); + H(2,2) = + SD1U.X() * SD1U.X() + + SD1U.Y() * SD1U.Y() + + SD1U.Z() * SD1U.Z() + - (CD0.X() - SD0.X()) * SD2UU.X() + - (CD0.X() - SD0.X()) * SD2UU.Y() + - (CD0.X() - SD0.X()) * SD2UU.Z(); + H(2,3) = + SD1U.X() * SD1V.X() + + SD1U.Y() * SD1V.Y() + + SD1U.Z() * SD1V.Z() + - (CD0.X() - SD0.X()) * SD2UV.X() + - (CD0.X() - SD0.X()) * SD2UV.Y() + - (CD0.X() - SD0.X()) * SD2UV.Z(); + H(3,1) = H(1,3); + H(3,2) = H(2,3); + H(3,3) = + SD1V.X() * SD1V.X() + + SD1V.Y() * SD1V.Y() + + SD1V.Z() * SD1V.Z() + - (CD0.X() - SD0.X()) * SD2VV.X() + - (CD0.X() - SD0.X()) * SD2VV.Y() + - (CD0.X() - SD0.X()) * SD2VV.Z(); +} + +//======================================================================= +//function : checkInputData +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCS::checkInputData(const math_Vector &X, + Standard_Real &cu, + Standard_Real &su, + Standard_Real &sv) +{ + Standard_Integer aStartIndex = X.Lower(); + cu = X(aStartIndex); + su = X(aStartIndex + 1); + sv = X(aStartIndex + 2); + + if (cu < myC->FirstParameter() || + cu > myC->LastParameter() || + su < myS->FirstUParameter() || + su > myS->LastUParameter() || + sv < myS->FirstVParameter() || + sv > myS->LastVParameter()) + { + return Standard_False; + } + return Standard_True; +} + +//======================================================================= +//function : Extrema_GlobOptFuncCS +//purpose : Constructor +//======================================================================= +Extrema_GlobOptFuncCS::Extrema_GlobOptFuncCS(const Adaptor3d_Curve *C, + const Adaptor3d_Surface *S) +: myC(C), + myS(S) +{ +} + +//======================================================================= +//function : NbVariables +//purpose : +//======================================================================= +Standard_Integer Extrema_GlobOptFuncCS::NbVariables() const +{ + return 3; +} + +//======================================================================= +//function : Value +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCS::Value(const math_Vector &X, + Standard_Real &F) +{ + Standard_Real cu, su, sv; + if (!checkInputData(X, cu, su, sv)) + return Standard_False; + + value(cu, su, sv, F); + return Standard_True; +} + +//======================================================================= +//function : Gradient +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCS::Gradient(const math_Vector &X, + math_Vector &G) +{ + Standard_Real cu, su, sv; + if (!checkInputData(X, cu, su, sv)) + return Standard_False; + + gradient(cu, su, sv, G); + return Standard_True; +} + +//======================================================================= +//function : Values +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCS::Values(const math_Vector &X, + Standard_Real &F, + math_Vector &G) +{ + Standard_Real cu, su, sv; + if (!checkInputData(X, cu, su, sv)) + return Standard_False; + + value(cu, su, sv, F); + gradient(cu, su, sv, G); + return Standard_True; +} + +//======================================================================= +//function : Values +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCS::Values(const math_Vector &X, + Standard_Real &F, + math_Vector &G, + math_Matrix &H) +{ + Standard_Real cu, su, sv; + if (!checkInputData(X, cu, su, sv)) + return Standard_False; + + value(cu, su, sv, F); + gradient(cu, su, sv, G); + hessian(cu, su, sv, H); + return Standard_True; +} \ No newline at end of file diff --git a/src/Extrema/Extrema_GlobOptFuncCS.hxx b/src/Extrema/Extrema_GlobOptFuncCS.hxx new file mode 100644 index 0000000000..228b4f3111 --- /dev/null +++ b/src/Extrema/Extrema_GlobOptFuncCS.hxx @@ -0,0 +1,81 @@ +// Created on: 2014-06-23 +// Created by: Alexander Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement + +#ifndef _Extrema_GlobOptFuncCS_HeaderFile +#define _Extrema_GlobOptFuncCS_HeaderFile + + +#include +#include +#include +#include +#include + +//! This class implements function which calculate Eucluidean distance +//! between point on curve and point on surface in case of continuity is C2. +class Extrema_GlobOptFuncCS : public math_MultipleVarFunctionWithHessian +{ +public: + + //! Curve and surface should exist during all the lifetime of Extrema_GlobOptFuncCS. + Standard_EXPORT Extrema_GlobOptFuncCS(const Adaptor3d_Curve *C, + const Adaptor3d_Surface *S); + + Standard_EXPORT virtual Standard_Integer NbVariables() const; + + Standard_EXPORT virtual Standard_Boolean Value(const math_Vector &theX, + Standard_Real &theF); + + Standard_EXPORT virtual Standard_Boolean Gradient(const math_Vector &theX, + math_Vector &theG); + + Standard_EXPORT virtual Standard_Boolean Values(const math_Vector &theX, + Standard_Real &theF, + math_Vector &theG); + + Standard_EXPORT virtual Standard_Boolean Values(const math_Vector &theX, + Standard_Real &theF, + math_Vector &theG, + math_Matrix &theH); + +private: + + Standard_Boolean checkInputData(const math_Vector &X, + Standard_Real &cu, + Standard_Real &su, + Standard_Real &sv); + + void value(Standard_Real cu, + Standard_Real su, + Standard_Real sv, + Standard_Real &F); + + void gradient(Standard_Real cu, + Standard_Real su, + Standard_Real sv, + math_Vector &G); + + void hessian(Standard_Real cu, + Standard_Real su, + Standard_Real sv, + math_Matrix &H); + + Extrema_GlobOptFuncCS & operator = (const Extrema_GlobOptFuncCS & theOther); + + const Adaptor3d_Curve *myC; + const Adaptor3d_Surface *myS; +}; + +#endif diff --git a/src/Extrema/FILES b/src/Extrema/FILES index f091b03b59..84b25ea88f 100644 --- a/src/Extrema/FILES +++ b/src/Extrema/FILES @@ -1,3 +1,5 @@ Extrema_HUBTreeOfSphere.hxx Extrema_GlobOptFuncCC.hxx Extrema_GlobOptFuncCC.cxx +Extrema_GlobOptFuncCS.hxx +Extrema_GlobOptFuncCS.cxx diff --git a/src/math/FILES b/src/math/FILES index 60b777e2d9..e608228a2a 100755 --- a/src/math/FILES +++ b/src/math/FILES @@ -10,4 +10,9 @@ math_IntegerVector.hxx math_IntegerVector.cxx math_SingleTab.hxx math_GlobOptMin.hxx -math_GlobOptMin.cxx \ No newline at end of file +math_GlobOptMin.cxx +math_PSO.hxx +math_PSO.cxx +math_BullardGenerator.hxx +math_PSOParticlesPool.hxx +math_PSOParticlesPool.cxx \ No newline at end of file diff --git a/src/math/math.cdl b/src/math/math.cdl index ea190ffe9d..7fe5438b52 100644 --- a/src/math/math.cdl +++ b/src/math/math.cdl @@ -52,6 +52,9 @@ is deferred class FunctionSet; deferred class FunctionSetWithDerivatives; imported GlobOptMin; + imported PSO; + imported PSOParticlesPool; + imported BullardGenerator; class IntegerRandom; class Gauss; class GaussLeastSquare; diff --git a/src/math/math_BullardGenerator.hxx b/src/math/math_BullardGenerator.hxx new file mode 100644 index 0000000000..ff604ec720 --- /dev/null +++ b/src/math/math_BullardGenerator.hxx @@ -0,0 +1,57 @@ +// Created on: 2014-07-18 +// Created by: Alexander Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _math_BullardGenerator_HeaderFile +#define _math_BullardGenerator_HeaderFile + +#include + +//! Fast random number generator (the algorithm proposed by Ian C. Bullard). +class math_BullardGenerator +{ +public: + + //! Creates new Xorshift 64-bit RNG. + Standard_EXPORT math_BullardGenerator(unsigned int theSeed = 1) + : myStateHi (theSeed) + { + myStateLo = theSeed ^ 0x49616E42; + } + + //! Generates new 64-bit integer value. + Standard_EXPORT unsigned int NextInt() + { + myStateHi = (myStateHi >> 2) + (myStateHi << 2); + + myStateHi += myStateLo; + myStateLo += myStateHi; + + return myStateHi; + } + + //! Generates new floating-point value. + Standard_EXPORT Standard_Real NextReal() + { + return NextInt() / static_cast (0xFFFFFFFFu); + } + +private: + + unsigned int myStateHi; + unsigned int myStateLo; + +}; + +#endif \ No newline at end of file diff --git a/src/math/math_PSO.cxx b/src/math/math_PSO.cxx new file mode 100644 index 0000000000..93fecf9acf --- /dev/null +++ b/src/math/math_PSO.cxx @@ -0,0 +1,262 @@ +// Created on: 2014-07-18 +// Created by: Alexander Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include + +const Standard_Real aBorderDivisor = 1.0e+4; + +//======================================================================= +//function : math_PSO +//purpose : Constructor +//======================================================================= +math_PSO::math_PSO(math_MultipleVarFunction* theFunc, + const math_Vector& theLowBorder, + const math_Vector& theUppBorder, + const math_Vector& theSteps, + const Standard_Integer theNbParticles, + const Standard_Integer theNbIter) +: myLowBorder(1, theFunc->NbVariables()), + myUppBorder(1, theFunc->NbVariables()), + mySteps(1, theFunc->NbVariables()) +{ + myN = theFunc->NbVariables(); + myNbParticles = theNbParticles; + myNbIter = theNbIter; + myFunc = theFunc; + + myLowBorder = theLowBorder; + myUppBorder = theUppBorder; + mySteps = theSteps; +} + +//======================================================================= +//function : math_PSO +//purpose : Perform +//======================================================================= +void math_PSO::Perform(math_PSOParticlesPool& theParticles, + Standard_Integer theNbParticles, + Standard_Real& theValue, + math_Vector& theOutPnt, + const Standard_Integer theNbIter) +{ + performPSOWithGivenParticles(theParticles, theNbParticles, theValue, theOutPnt, theNbIter); +} + +//======================================================================= +//function : math_PSO +//purpose : Perform +//======================================================================= +void math_PSO::Perform(const math_Vector& theSteps, + Standard_Real& theValue, + math_Vector& theOutPnt, + const Standard_Integer theNbIter) +{ + // Initialization. + math_Vector aMinUV(1, myN), aMaxUV(1, myN); + aMinUV = myLowBorder + (myUppBorder - myLowBorder) / aBorderDivisor; + aMaxUV = myUppBorder - (myUppBorder - myLowBorder) / aBorderDivisor; + myNbIter = theNbIter; + mySteps = theSteps; + + // To generate initial distribution it is necessary to have grid steps. + math_PSOParticlesPool aPool(myNbParticles, myN); + + // Generate initial particles distribution. + Standard_Boolean isRegularGridFinished = Standard_False; + Standard_Real aCurrValue; + math_Vector aCurrPoint(1, myN); + + PSO_Particle* aParticle = aPool.GetWorstParticle(); + aCurrPoint = aMinUV; + do + { + myFunc->Value(aCurrPoint, aCurrValue); + + if (aCurrValue * aCurrValue < aParticle->Distance) + { + Standard_Integer aDimIdx; + for(aDimIdx = 0; aDimIdx < myN; ++aDimIdx) + { + aParticle->Position[aDimIdx] = aCurrPoint(aDimIdx + 1); + aParticle->BestPosition[aDimIdx] = aCurrPoint(aDimIdx + 1); + } + aParticle->Distance = aCurrValue * aCurrValue; + aParticle->BestDistance = aCurrValue * aCurrValue; + + aParticle = aPool.GetWorstParticle(); + } + + // Step. + aCurrPoint(1) += Max(mySteps(1), 1.0e-15); // Avoid too small step + for(Standard_Integer aDimIdx = 1; aDimIdx < myN; ++aDimIdx) + { + if(aCurrPoint(aDimIdx) > aMaxUV(aDimIdx)) + { + aCurrPoint(aDimIdx) = aMinUV(aDimIdx); + aCurrPoint(aDimIdx + 1) += mySteps(aDimIdx + 1); + } + else + break; + } + + // Stop criteria. + if(aCurrPoint(myN) > aMaxUV(myN)) + isRegularGridFinished = Standard_True; + } + while(!isRegularGridFinished); + + performPSOWithGivenParticles(aPool, myNbParticles, theValue, theOutPnt, theNbIter); +} + +//======================================================================= +//function : math_PSO +//purpose : Perform +//======================================================================= +void math_PSO::performPSOWithGivenParticles(math_PSOParticlesPool& theParticles, + Standard_Integer theNbParticles, + Standard_Real& theValue, + math_Vector& theOutPnt, + const Standard_Integer theNbIter) +{ + math_Vector aMinUV(1, myN), aMaxUV(1, myN); + aMinUV = myLowBorder + (myUppBorder - myLowBorder) / aBorderDivisor; + aMaxUV = myUppBorder - (myUppBorder - myLowBorder) / aBorderDivisor; + myNbIter = theNbIter; + myNbParticles = theNbParticles; + math_PSOParticlesPool& aParticles = theParticles; + + // Current particle data. + math_Vector aCurrPoint(1, myN); + math_Vector aBestGlobalPosition(1, myN); + PSO_Particle* aParticle = 0; + + + // Generate initial particle velocities. + math_BullardGenerator aRandom; + for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx) + { + aParticle = aParticles.GetParticle(aPartIdx); + for(Standard_Integer aDimIdx = 1; aDimIdx <= myN; ++aDimIdx) + { + const Standard_Real aKsi = aRandom.NextReal(); + aParticle->Velocity[aDimIdx - 1] = mySteps(aDimIdx) * (aKsi - 0.5) * 2.0; + } + } + + aParticle = aParticles.GetBestParticle(); + for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx) + aBestGlobalPosition(aDimIdx + 1) = aParticle->Position[aDimIdx]; + Standard_Real aBestGlobalDistance = aParticle->Distance; + + // This velocity is used for detecting stagnation state. + math_Vector aTerminationVelocity(1, myN); + aTerminationVelocity = mySteps / 2048.0; + // This velocity shows minimal velocity on current step. + math_Vector aMinimalVelocity(1, myN); + + // Run PSO iterations + for (Standard_Integer aStep = 1; aStep < myNbIter; ++aStep) + { + aMinimalVelocity.Init(RealLast()); + + for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx) + { + const Standard_Real aKsi1 = aRandom.NextReal(); + const Standard_Real aKsi2 = aRandom.NextReal(); + + aParticle = aParticles.GetParticle(aPartIdx); + + const Standard_Real aRetentWeight = 0.72900; + const Standard_Real aPersonWeight = 1.49445; + const Standard_Real aSocialWeight = 1.49445; + + for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx) + { + aParticle->Velocity[aDimIdx] = aParticle->Velocity[aDimIdx] * aRetentWeight + + (aParticle->BestPosition[aDimIdx] - aParticle->Position[aDimIdx]) * (aPersonWeight * aKsi1) + + (aBestGlobalPosition(aDimIdx + 1) - aParticle->Position[aDimIdx]) * (aSocialWeight * aKsi2); + + aParticle->Position[aDimIdx] += aParticle->Velocity[aDimIdx]; + aParticle->Position[aDimIdx] = Min(Max(aParticle->Position[aDimIdx], aMinUV(aDimIdx + 1)), + aMaxUV(aDimIdx + 1)); + aCurrPoint(aDimIdx + 1) = aParticle->Position[aDimIdx]; + + aMinimalVelocity(aDimIdx + 1) = Min(Abs(aParticle->Velocity[aDimIdx]), + aMinimalVelocity(aDimIdx + 1)); + } + + myFunc->Value(aCurrPoint, aParticle->Distance); + aParticle->Distance *= aParticle->Distance; + if(aParticle->Distance < aParticle->BestDistance) + { + aParticle->BestDistance = aParticle->Distance; + for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx) + aParticle->BestPosition[aDimIdx] = aParticle->Position[aDimIdx]; + + if(aParticle->Distance < aBestGlobalDistance) + { + aBestGlobalDistance = aParticle->Distance; + for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx) + aBestGlobalPosition(aDimIdx + 1) = aParticle->Position[aDimIdx]; + } + } + } + + Standard_Boolean isTerminalVelocityReached = Standard_True; + for(Standard_Integer aDimIdx = 1; aDimIdx <= myN; ++aDimIdx) + { + if(aMinimalVelocity(aDimIdx) > aTerminationVelocity(aDimIdx)) + { + isTerminalVelocityReached = Standard_False; + break; + } + } + + if(isTerminalVelocityReached) + { + // Minimum number of steps + const Standard_Integer aMinSteps = 16; + + if (aStep > aMinSteps) + { + break; + } + + for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx) + { + const Standard_Real aKsi = aRandom.NextReal(); + + aParticle = aParticles.GetParticle(aPartIdx); + + for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx) + { + if (aParticle->Position[aDimIdx] == aMinUV(aDimIdx + 1) || + aParticle->Position[aDimIdx] == aMaxUV(aDimIdx + 1)) + { + aParticle->Velocity[aDimIdx] = aParticle->Position[aDimIdx] == aMinUV(aDimIdx + 1) ? + mySteps(aDimIdx + 1) * aKsi : -mySteps(aDimIdx + 1) * aKsi; + } + else + { + aParticle->Velocity[aDimIdx] = mySteps(aDimIdx + 1) * (aKsi - 0.5) * 2.0; + } + } + } + } + } + theValue = aBestGlobalDistance; + theOutPnt = aBestGlobalPosition; +} diff --git a/src/math/math_PSO.hxx b/src/math/math_PSO.hxx new file mode 100644 index 0000000000..308c70b70d --- /dev/null +++ b/src/math/math_PSO.hxx @@ -0,0 +1,80 @@ +// Created on: 2014-07-18 +// Created by: Alexander Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _math_PSO_HeaderFile +#define _math_PSO_HeaderFile + +#include +#include +#include + +class math_PSOParticlesPool; + +//! In this class implemented variation of Particle Swarm Optimization (PSO) method. +//! A. Ismael F. Vaz, L. N. Vicente +//! "A particle swarm pattern search method for bound constrained global optimization" +class math_PSO +{ +public: + + /** + * Constructor. + * + * @param theFunc defines the objective function. It should exist during all lifetime of class instance. + * @param theLowBorder defines lower border of search space. + * @param theUppBorder defines upper border of search space. + * @param theSteps defines steps of regular grid, used for particle generation. + This parameter used to define stop condition (TerminalVelocity). + * @param theNbParticles defines number of particles. + * @param theNbIter defines maximum number of iterations. + */ + Standard_EXPORT math_PSO(math_MultipleVarFunction* theFunc, + const math_Vector& theLowBorder, + const math_Vector& theUppBorder, + const math_Vector& theSteps, + const Standard_Integer theNbParticles = 32, + const Standard_Integer theNbIter = 100); + + //! Perform computations, particles array is constructed inside of this function. + Standard_EXPORT void Perform(const math_Vector& theSteps, + Standard_Real& theValue, + math_Vector& theOutPnt, + const Standard_Integer theNbIter = 100); + + //! Perform computations with given particles array. + Standard_EXPORT void Perform(math_PSOParticlesPool& theParticles, + Standard_Integer theNbParticles, + Standard_Real& theValue, + math_Vector& theOutPnt, + const Standard_Integer theNbIter = 100); + +private: + + void performPSOWithGivenParticles(math_PSOParticlesPool& theParticles, + Standard_Integer theNbParticles, + Standard_Real& theValue, + math_Vector& theOutPnt, + const Standard_Integer theNbIter = 100); + + math_MultipleVarFunction *myFunc; + math_Vector myLowBorder; // Lower border. + math_Vector myUppBorder; // Upper border. + math_Vector mySteps; // steps used in PSO algorithm. + Standard_Integer myN; // Dimension count. + Standard_Integer myNbParticles; // Particles number. + Standard_Integer myNbIter; +}; + +#endif diff --git a/src/math/math_PSOParticlesPool.cxx b/src/math/math_PSOParticlesPool.cxx new file mode 100644 index 0000000000..f3b7bf6ca6 --- /dev/null +++ b/src/math/math_PSOParticlesPool.cxx @@ -0,0 +1,79 @@ +// Created on: 2014-07-18 +// Created by: Alexander Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include + + +//======================================================================= +//function : math_PSOParticlesPool +//purpose : Constructor +//======================================================================= +math_PSOParticlesPool::math_PSOParticlesPool(const Standard_Integer theParticlesCount, + const Standard_Integer theDimensionCount) +: myParticlesPool(1, theParticlesCount), + myMemory(0, theParticlesCount * (theDimensionCount // Position + + theDimensionCount // Velocity + + theDimensionCount) // BestPosition + - 1) +{ + myParticlesCount = theParticlesCount; + myDimensionCount = theDimensionCount; + + // Pointers adjusting. + Standard_Integer aParIdx, aShiftIdx; + for(aParIdx = 1; aParIdx <= myParticlesCount; ++aParIdx) + { + aShiftIdx = (theDimensionCount * 3) * (aParIdx - 1); + myParticlesPool(aParIdx).Position = &myMemory(aShiftIdx); + myParticlesPool(aParIdx).Velocity = &myMemory(aShiftIdx + theDimensionCount); + myParticlesPool(aParIdx).BestPosition = &myMemory(aShiftIdx + 2 * theDimensionCount); + } +} + +//======================================================================= +//function : ~math_PSOParticlesPool +//purpose : Destructor +//======================================================================= +math_PSOParticlesPool::~math_PSOParticlesPool() +{ +} + +//======================================================================= +//function : GetParticle +//purpose : +//======================================================================= +PSO_Particle* math_PSOParticlesPool::GetParticle(const Standard_Integer theIdx) +{ + return &myParticlesPool(theIdx); +} + +//======================================================================= +//function : GetBestParticle +//purpose : +//======================================================================= +PSO_Particle* math_PSOParticlesPool::GetBestParticle() +{ + return &*std::min_element(myParticlesPool.begin(), myParticlesPool.end()); +} + +//======================================================================= +//function : GetWorstParticle +//purpose : +//======================================================================= +PSO_Particle* math_PSOParticlesPool::GetWorstParticle() +{ + return &*std::max_element(myParticlesPool.begin(), myParticlesPool.end()); +} \ No newline at end of file diff --git a/src/math/math_PSOParticlesPool.hxx b/src/math/math_PSOParticlesPool.hxx new file mode 100644 index 0000000000..7ddfd1de0a --- /dev/null +++ b/src/math/math_PSOParticlesPool.hxx @@ -0,0 +1,73 @@ +// Created on: 2014-07-18 +// Created by: Alexander Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _math_PSOParticlesPool_HeaderFile +#define _math_PSOParticlesPool_HeaderFile + +#include + +//! Describes particle pool for using in PSO algorithm. +//! Indexes: +//! 0 <= aDimidx <= myDimensionCount - 1 +struct PSO_Particle +{ + Standard_Real* Position; // Data for pointers allocated within PSOParticlesPool instance. + Standard_Real* Velocity; // Not need to delete it manually. + Standard_Real* BestPosition; + Standard_Real Distance; + Standard_Real BestDistance; + + PSO_Particle() + { + Distance = RealLast(); + BestDistance = RealLast(); + Position = 0; + Velocity = 0; + BestPosition = 0; + } + + //! Compares the particles according to their distances. + bool operator< (const PSO_Particle& thePnt) const + { + return Distance < thePnt.Distance; + } +}; + +// Indexes: +// 1 <= aParticleIdx <= myParticlesCount +class math_PSOParticlesPool +{ +public: + + Standard_EXPORT math_PSOParticlesPool(const Standard_Integer theParticlesCount, + const Standard_Integer theDimensionCount); + + Standard_EXPORT PSO_Particle* GetParticle(const Standard_Integer theIdx); + + Standard_EXPORT PSO_Particle* GetBestParticle(); + + Standard_EXPORT PSO_Particle* GetWorstParticle(); + + Standard_EXPORT ~math_PSOParticlesPool(); + +private: + + NCollection_Array1 myParticlesPool; + NCollection_Array1 myMemory; // Stores particles vector data. + Standard_Integer myParticlesCount; + Standard_Integer myDimensionCount; +}; + +#endif