Porting Particle Swarm Optimization (PSO) to math package.
uses POnSurf from Extrema,
POnCurv from Extrema,
- Pnt from gp,
+ Pnt from gp,
HArray1OfPnt from TColgp,
HArray2OfPnt from TColgp,
FuncExtCS from Extrema,
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;
#include <Extrema_GenExtCS.ixx>
-#include <math_Vector.hxx>
-#include <math_FunctionSetRoot.hxx>
-#include <Precision.hxx>
-
-#include <Geom_Line.hxx>
#include <Adaptor3d_HCurve.hxx>
-#include <GeomAdaptor_Curve.hxx>
-
#include <Extrema_ExtCC.hxx>
-#include <Extrema_POnCurv.hxx>
#include <Extrema_ExtPS.hxx>
+#include <Extrema_GlobOptFuncCS.hxx>
+#include <Extrema_POnCurv.hxx>
+#include <Geom_Line.hxx>
+#include <GeomAdaptor_Curve.hxx>
+#include <math_FunctionSetRoot.hxx>
+#include <math_PSO.hxx>
+#include <math_PSOParticlesPool.hxx>
+#include <math_Vector.hxx>
+#include <Precision.hxx>
+#include <TColgp_HArray1OfPnt.hxx>
-#include <algorithm>
-#include <NCollection_Vec3.hxx>
-#include <NCollection_Array1.hxx>
+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;
//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);
//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);
//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();
//function : Initialize
//purpose :
//=======================================================================
-
void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S,
const Standard_Integer NbU,
const Standard_Integer NbV,
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;
//function : Perform
//purpose :
//=======================================================================
-
void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C,
const Standard_Integer NbT,
const Standard_Real Tol1)
Perform(C, NbT, mytmin, mytsup,Tol1);
}
-namespace
-{
- //! Vector type for storing curve and surface parameters.
- typedef NCollection_Vec3<Standard_Real> 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<Standard_Real> (0xFFFFFFFFu);
- }
- };
-}
-
//=======================================================================
//function : Perform
//purpose :
//=======================================================================
-
void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C,
const Standard_Integer NbT,
const Standard_Real tmin,
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;
-
- UVsup(1) = mytsup;
- UVsup(2) = trimusup;
- UVsup(3) = trimvsup;
+ math_Vector TUV(1,3), TUVinf(1,3), TUVsup(1,3);
+ TUVinf(1) = mytmin;
+ TUVinf(2) = trimumin;
+ TUVinf(3) = trimvmin;
+
+ 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.
{
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())
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<Extrema_Particle> 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,
{
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<Extrema_Particle>::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;
- 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<Extrema_Particle>::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();
+ aParticle->Distance = aSqDist;
+ aParticle->BestDistance = aSqDist;
- 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;
+ aParticle = aParticles.GetWorstParticle();
}
}
}
}
+ 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;
//function : IsDone
//purpose :
//=======================================================================
-
Standard_Boolean Extrema_GenExtCS::IsDone() const
{
return myDone;
//function : NbExt
//purpose :
//=======================================================================
-
Standard_Integer Extrema_GenExtCS::NbExt() const
{
if (!IsDone()) { StdFail_NotDone::Raise(); }
//function : Value
//purpose :
//=======================================================================
-
Standard_Real Extrema_GenExtCS::SquareDistance(const Standard_Integer N) const
{
if (!IsDone()) { StdFail_NotDone::Raise(); }
//function : PointOnCurve
//purpose :
//=======================================================================
-
const Extrema_POnCurv& Extrema_GenExtCS::PointOnCurve(const Standard_Integer N) const
{
if (!IsDone()) { StdFail_NotDone::Raise(); }
//function : PointOnSurface
//purpose :
//=======================================================================
-
const Extrema_POnSurf& Extrema_GenExtCS::PointOnSurface(const Standard_Integer N) const
{
if (!IsDone()) { StdFail_NotDone::Raise(); }
//function : BidonSurface
//purpose :
//=======================================================================
-
Adaptor3d_SurfacePtr Extrema_GenExtCS::BidonSurface() const
{
return (Adaptor3d_SurfacePtr)0L;
//function : BidonCurve
//purpose :
//=======================================================================
-
Adaptor3d_CurvePtr Extrema_GenExtCS::BidonCurve() const
{
return (Adaptor3d_CurvePtr)0L;
--- /dev/null
+// 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 <Extrema_GlobOptFuncCS.hxx>
+
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+#include <math_Vector.hxx>
+#include <Standard_Integer.hxx>
+#include <Standard_OutOfRange.hxx>
+
+//! 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
--- /dev/null
+// 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 <Adaptor3d_Curve.hxx>
+#include <Adaptor3d_Surface.hxx>
+#include <math_Matrix.hxx>
+#include <math_Vector.hxx>
+#include <math_MultipleVarFunctionWithHessian.hxx>
+
+//! 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
Extrema_HUBTreeOfSphere.hxx
Extrema_GlobOptFuncCC.hxx
Extrema_GlobOptFuncCC.cxx
+Extrema_GlobOptFuncCS.hxx
+Extrema_GlobOptFuncCS.cxx
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
deferred class FunctionSet;
deferred class FunctionSetWithDerivatives;
imported GlobOptMin;
+ imported PSO;
+ imported PSOParticlesPool;
+ imported BullardGenerator;
class IntegerRandom;
class Gauss;
class GaussLeastSquare;
--- /dev/null
+// 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 <Standard_Real.hxx>
+
+//! 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<Standard_Real> (0xFFFFFFFFu);
+ }
+
+private:
+
+ unsigned int myStateHi;
+ unsigned int myStateLo;
+
+};
+
+#endif
\ No newline at end of file
--- /dev/null
+// 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 <math_PSO.hxx>
+#include <math_PSOParticlesPool.hxx>
+
+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;
+}
--- /dev/null
+// 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 <math_BullardGenerator.hxx>
+#include <math_MultipleVarFunction.hxx>
+#include <math_Vector.hxx>
+
+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
--- /dev/null
+// 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 <math_PSOParticlesPool.hxx>
+#include <algorithm>
+
+
+//=======================================================================
+//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
--- /dev/null
+// 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 <NCollection_Array1.hxx>
+
+//! 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<PSO_Particle> myParticlesPool;
+ NCollection_Array1<Standard_Real> myMemory; // Stores particles vector data.
+ Standard_Integer myParticlesCount;
+ Standard_Integer myDimensionCount;
+};
+
+#endif
projects / occt-copy.git / commitdiff