0025086: Implementation PSO in package math

Porting Particle Swarm Optimization (PSO) to math package.

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;
+    mytmin     : Real;
-    mytsup    : Real;
+    mytsup     : Real;
-    myumin    : Real;
+    myumin     : Real;
-    myusup    : Real;
+    myusup     : Real;
-    myvmin    : Real;
+    myvmin     : Real;
-    myvsup    : 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;

src/Extrema/Extrema_GenExtCS.cxx

@@ -16,27 +16,27 @@
 
 #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>
+const Standard_Real aMaxParamVal = 1.0e+10;
-#include <NCollection_Vec3.hxx>
+const Standard_Real aBorderDivisor = 1.0e+4;
-#include <NCollection_Array1.hxx>
 
 //=======================================================================
 //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,
-Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C, 
+                                   const Adaptor3d_Surface& S,
-  const Adaptor3d_Surface& S, 
+                                   const Standard_Integer NbT,
-  const Standard_Integer NbT, 
+                                   const Standard_Integer NbU,
-  const Standard_Integer NbU, 
+                                   const Standard_Integer NbV,
-  const Standard_Integer NbV, 
+                                   const Standard_Real Tol1,
-  const Standard_Real Tol1, 
+                                   const Standard_Real Tol2)
-  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, 
+Extrema_GenExtCS::Extrema_GenExtCS (const Adaptor3d_Curve& C,
-                                    const Adaptor3d_Surface& S, 
+                                    const Adaptor3d_Surface& S,
-                                    const Standard_Integer NbT, 
+                                    const Standard_Integer NbT,
-                                    const Standard_Integer NbU, 
+                                    const Standard_Integer NbU,
-                                    const Standard_Integer NbV, 
+                                    const Standard_Integer NbV,
-                                    const Standard_Real tmin, 
+                                    const Standard_Real tmin,
-                                    const Standard_Real tsup, 
+                                    const Standard_Real tsup,
-                                    const Standard_Real Umin, 
+                                    const Standard_Real Umin,
                                     const Standard_Real Usup,
-                                    const Standard_Real Vmin, 
+                                    const Standard_Real Vmin,
-                                    const Standard_Real Vsup, 
+                                    const Standard_Real Vsup,
-                                    const Standard_Real Tol1, 
+                                    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,
-void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S, 
+                                   const Standard_Integer NbU,
-                                   const Standard_Integer NbU, 
+                                   const Standard_Integer NbV,
-                                   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 aTrimMaxU = Precision::IsInfinite (myusup) ?  aMaxParamVal : myusup;
-  const Standard_Real aTrimMinU = Precision::IsInfinite (myumin) ? -1.0e+10 : myumin;
+  const Standard_Real aTrimMinU = Precision::IsInfinite (myumin) ? -aMaxParamVal : myumin;
-  const Standard_Real aTrimMaxV = Precision::IsInfinite (myvsup) ?  1.0e+10 : myvsup;
+  const Standard_Real aTrimMaxV = Precision::IsInfinite (myvsup) ?  aMaxParamVal : myvsup;
-  const Standard_Real aTrimMinV = Precision::IsInfinite (myvmin) ? -1.0e+10 : myvmin;
+  const Standard_Real aTrimMinV = Precision::IsInfinite (myvmin) ? -aMaxParamVal : myvmin;
 
-  const Standard_Real aMinU = aTrimMinU + (aTrimMaxU - aTrimMinU) / 1.0e4;
+  const Standard_Real aMinU = aTrimMinU + (aTrimMaxU - aTrimMinU) / aBorderDivisor;
-  const Standard_Real aMinV = aTrimMinV + (aTrimMaxV - aTrimMinV) / 1.0e4;
+  const Standard_Real aMinV = aTrimMinV + (aTrimMaxV - aTrimMinV) / aBorderDivisor;
-  const Standard_Real aMaxU = aTrimMaxU - (aTrimMaxU - aTrimMinU) / 1.0e4;
+  const Standard_Real aMaxU = aTrimMaxU - (aTrimMaxU - aTrimMinU) / aBorderDivisor;
-  const Standard_Real aMaxV = aTrimMaxV - (aTrimMaxV - aTrimMinV) / 1.0e4;
+  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<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,
@@ -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);
+  math_Vector TUV(1,3), TUVinf(1,3), TUVsup(1,3);
-  UVinf(1) = mytmin;
+  TUVinf(1) = mytmin;
-  UVinf(2) = trimumin;
+  TUVinf(2) = trimumin;
-  UVinf(3) = trimvmin;
+  TUVinf(3) = trimvmin;
 
-  UVsup(1) = mytsup;
+  TUVsup(1) = mytsup;
-  UVsup(2) = trimusup;
+  TUVsup(2) = trimusup;
-  UVsup(3) = trimvsup;
+  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));
+            aPreciser.Point(iPExt).Parameter(TUV(2),TUV(3));
-            math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
+            math_FunctionSetRoot S1 (myF,TUV,Tol,TUVinf,TUVsup);
           }
         }
         else
         {
           // Failed... try the point on iso line
-          UV(2) = dfUFirst;
+          TUV(2) = dfUFirst;
-          UV(3) = aP2.Parameter();
+          TUV(3) = aP2.Parameter();
-          math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
+          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,
+    math_Vector aMinTUV(1,3);
-                                UVinf(2) + (UVsup(2) - UVinf(2)) / 1.0e4,
+    aMinTUV = TUVinf + (TUVsup - TUVinf) / aBorderDivisor;
-                                UVinf(3) + (UVsup(3) - UVinf(3)) / 1.0e4);
 
-    const Extrema_Vec3d aMaxUV (UVsup(1) - (UVsup(1) - UVinf(1)) / 1.0e4,
+    math_Vector aMaxTUV(1,3); 
-                                UVsup(2) - (UVsup(2) - UVinf(2)) / 1.0e4,
+    aMaxTUV = TUVsup - (TUVsup - TUVinf) / aBorderDivisor;
-                                UVsup(3) - (UVsup(3) - UVinf(3)) / 1.0e4);
 
-    Standard_Real aStepCU = (aMaxUV.x() - aMinUV.x()) / mytsample;
+    Standard_Real aStepCU = (aMaxTUV(1) - aMinTUV(1)) / mytsample;
-    Standard_Real aStepSU = (aMaxUV.y() - aMinUV.y()) / myusample;
+    Standard_Real aStepSU = (aMaxTUV(2) - aMinTUV(2)) / myusample;
-    Standard_Real aStepSV = (aMaxUV.z() - aMinUV.z()) / myvsample;
+    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 = Min(aMaxNbNodes,
-          mytsample * aResolutionCU / aMinResolution));
+                        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;
+
+            aParticle->Distance     = aSqDist;
+            aParticle->BestDistance = aSqDist;
+
+            aParticle = aParticles.GetWorstParticle();
           }
         }
       }
     }
 
-    Extrema_Random aRandom;
+    math_Vector aStep(1,3);
-
+    aStep(1) = aStepCU;
-    // Generate initial particle velocities
+    aStep(2) = aStepSU;
-    for (Standard_Integer aPartIdx = 1; aPartIdx <= aNbParticles; ++aPartIdx)
+    aStep(3) = aStepSV;
-    {
-      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();
-
-        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;
-          }
-        }
-      }
-    }
 
     // Find min approximation
-    UV(1) = aBestGlobalPosition.x();
+    Standard_Real aValue;
-    UV(2) = aBestGlobalPosition.y();
+    Extrema_GlobOptFuncCS aFunc(&C, myS);
-    UV(3) = aBestGlobalPosition.z();
+    math_PSO aPSO(&aFunc, TUVinf, TUVsup, aStep);
-    
+    aPSO.Perform(aParticles, aNbParticles, aValue, TUV);
-    math_FunctionSetRoot anA (myF, UV, Tol, UVinf, UVsup, 100, Standard_False);
+
+    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;

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 <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;
+}

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 <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

src/Extrema/FILES

@@ -1,3 +1,5 @@
 Extrema_HUBTreeOfSphere.hxx
 Extrema_GlobOptFuncCC.hxx
 Extrema_GlobOptFuncCC.cxx
+Extrema_GlobOptFuncCS.hxx
+Extrema_GlobOptFuncCS.cxx

src/math/FILES

@@ -10,4 +10,9 @@ math_IntegerVector.hxx
 math_IntegerVector.cxx
 math_SingleTab.hxx
 math_GlobOptMin.hxx
-math_GlobOptMin.cxx
+math_GlobOptMin.cxx
+math_PSO.hxx
+math_PSO.cxx
+math_BullardGenerator.hxx
+math_PSOParticlesPool.hxx
+math_PSOParticlesPool.cxx

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;

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 <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

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 <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;
+}

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 <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

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 <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());
+}

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 <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