src/math/math_PSO.cxx

   1 // Created on: 2014-07-18
   2 // Created by: Alexander Malyshev
   3 // Copyright (c) 2014-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16 #include <math_PSO.hxx>
  17 #include <math_PSOParticlesPool.hxx>
  18
  19 const Standard_Real aBorderDivisor = 1.0e+4;
  20
  21 //=======================================================================
  22 //function : math_PSO
  23 //purpose  : Constructor
  24 //=======================================================================
  25 math_PSO::math_PSO(math_MultipleVarFunction* theFunc,
  26                    const math_Vector& theLowBorder,
  27                    const math_Vector& theUppBorder,
  28                    const math_Vector& theSteps,
  29                    const Standard_Integer theNbParticles,
  30                    const Standard_Integer theNbIter)
  31 : myLowBorder(1, theFunc->NbVariables()),
  32   myUppBorder(1, theFunc->NbVariables()),
  33   mySteps(1, theFunc->NbVariables())
  34 {
  35   myN = theFunc->NbVariables();
  36   myNbParticles = theNbParticles;
  37   myNbIter = theNbIter;
  38   myFunc = theFunc;
  39
  40   myLowBorder = theLowBorder;
  41   myUppBorder = theUppBorder;
  42   mySteps     = theSteps;
  43 }
  44
  45 //=======================================================================
  46 //function : math_PSO
  47 //purpose  : Perform
  48 //=======================================================================
  49 void math_PSO::Perform(math_PSOParticlesPool& theParticles,
  50                        Standard_Integer theNbParticles,
  51                        Standard_Real& theValue,
  52                        math_Vector& theOutPnt,
  53                        const Standard_Integer theNbIter)
  54 {
  55   performPSOWithGivenParticles(theParticles, theNbParticles, theValue, theOutPnt, theNbIter);
  56 }
  57
  58 //=======================================================================
  59 //function : math_PSO
  60 //purpose  : Perform
  61 //=======================================================================
  62 void math_PSO::Perform(const math_Vector& theSteps,
  63                        Standard_Real& theValue,
  64                        math_Vector& theOutPnt,
  65                        const Standard_Integer theNbIter)
  66 {
  67   // Initialization.
  68   math_Vector aMinUV(1, myN), aMaxUV(1, myN);
  69   aMinUV = myLowBorder + (myUppBorder - myLowBorder) / aBorderDivisor;
  70   aMaxUV = myUppBorder - (myUppBorder - myLowBorder) / aBorderDivisor;
  71   myNbIter = theNbIter;
  72   mySteps = theSteps;
  73
  74   // To generate initial distribution it is necessary to have grid steps.
  75   math_PSOParticlesPool aPool(myNbParticles, myN);
  76
  77   // Generate initial particles distribution.
  78   Standard_Boolean isRegularGridFinished = Standard_False;
  79   Standard_Real aCurrValue;
  80   math_Vector   aCurrPoint(1, myN);
  81
  82   PSO_Particle* aParticle = aPool.GetWorstParticle();
  83   aCurrPoint = aMinUV;
  84   do
  85   {
  86     myFunc->Value(aCurrPoint, aCurrValue);
  87
  88     if (aCurrValue * aCurrValue < aParticle->Distance)
  89     {
  90       Standard_Integer aDimIdx;
  91       for(aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
  92       {
  93         aParticle->Position[aDimIdx] = aCurrPoint(aDimIdx + 1);
  94         aParticle->BestPosition[aDimIdx] = aCurrPoint(aDimIdx + 1);
  95       }
  96       aParticle->Distance     = aCurrValue * aCurrValue;
  97       aParticle->BestDistance = aCurrValue * aCurrValue;
  98
  99       aParticle = aPool.GetWorstParticle();
 100     }
 101
 102     // Step.
 103     aCurrPoint(1) += Max(mySteps(1), 1.0e-15); // Avoid too small step
 104     for(Standard_Integer aDimIdx = 1; aDimIdx < myN; ++aDimIdx)
 105     {
 106       if(aCurrPoint(aDimIdx) > aMaxUV(aDimIdx))
 107       {
 108         aCurrPoint(aDimIdx) = aMinUV(aDimIdx);
 109         aCurrPoint(aDimIdx + 1) += mySteps(aDimIdx + 1);
 110       }
 111       else
 112         break;
 113     }
 114
 115     // Stop criteria.
 116     if(aCurrPoint(myN) > aMaxUV(myN))
 117       isRegularGridFinished = Standard_True;
 118   }
 119   while(!isRegularGridFinished);
 120
 121   performPSOWithGivenParticles(aPool, myNbParticles, theValue, theOutPnt, theNbIter);
 122 }
 123
 124 //=======================================================================
 125 //function : math_PSO
 126 //purpose  : Perform
 127 //=======================================================================
 128 void math_PSO::performPSOWithGivenParticles(math_PSOParticlesPool& theParticles,
 129                                             Standard_Integer theNbParticles,
 130                                             Standard_Real& theValue,
 131                                             math_Vector& theOutPnt,
 132                                             const Standard_Integer theNbIter)
 133 {
 134   math_Vector aMinUV(1, myN), aMaxUV(1, myN);
 135   aMinUV = myLowBorder + (myUppBorder - myLowBorder) / aBorderDivisor;
 136   aMaxUV = myUppBorder - (myUppBorder - myLowBorder) / aBorderDivisor;
 137   myNbIter = theNbIter;
 138   myNbParticles = theNbParticles;
 139   math_PSOParticlesPool& aParticles = theParticles;
 140
 141   // Current particle data.
 142   math_Vector aCurrPoint(1, myN);
 143   math_Vector aBestGlobalPosition(1, myN);
 144   PSO_Particle* aParticle = 0;
 145
 146
 147   // Generate initial particle velocities.
 148   math_BullardGenerator aRandom;
 149   for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx)
 150   {
 151     aParticle = aParticles.GetParticle(aPartIdx);
 152     for(Standard_Integer aDimIdx = 1; aDimIdx <= myN; ++aDimIdx)
 153     {
 154       const Standard_Real aKsi = aRandom.NextReal();
 155       aParticle->Velocity[aDimIdx - 1] = mySteps(aDimIdx) * (aKsi - 0.5) * 2.0;
 156     }
 157   }
 158
 159   aParticle = aParticles.GetBestParticle();
 160   for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
 161     aBestGlobalPosition(aDimIdx + 1) = aParticle->Position[aDimIdx];
 162   Standard_Real aBestGlobalDistance = aParticle->Distance;
 163
 164   // This velocity is used for detecting stagnation state.
 165   math_Vector aTerminationVelocity(1, myN);
 166   aTerminationVelocity = mySteps / 2048.0;
 167   // This velocity shows minimal velocity on current step.
 168   math_Vector aMinimalVelocity(1, myN);
 169
 170   // Run PSO iterations
 171   for (Standard_Integer aStep = 1; aStep < myNbIter; ++aStep)
 172   {
 173     aMinimalVelocity.Init(RealLast());
 174
 175     for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx)
 176     {
 177       const Standard_Real aKsi1 = aRandom.NextReal();
 178       const Standard_Real aKsi2 = aRandom.NextReal();
 179
 180       aParticle = aParticles.GetParticle(aPartIdx);
 181
 182       const Standard_Real aRetentWeight = 0.72900;
 183       const Standard_Real aPersonWeight = 1.49445;
 184       const Standard_Real aSocialWeight = 1.49445;
 185
 186       for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
 187       {
 188         aParticle->Velocity[aDimIdx] = aParticle->Velocity[aDimIdx] * aRetentWeight
 189           + (aParticle->BestPosition[aDimIdx]  - aParticle->Position[aDimIdx]) * (aPersonWeight * aKsi1)
 190           + (aBestGlobalPosition(aDimIdx + 1) - aParticle->Position[aDimIdx]) * (aSocialWeight * aKsi2);
 191
 192         aParticle->Position[aDimIdx] += aParticle->Velocity[aDimIdx];
 193         aParticle->Position[aDimIdx] = Min(Max(aParticle->Position[aDimIdx], aMinUV(aDimIdx + 1)),
 194                                           aMaxUV(aDimIdx + 1));
 195         aCurrPoint(aDimIdx + 1) = aParticle->Position[aDimIdx];
 196
 197         aMinimalVelocity(aDimIdx + 1) = Min(Abs(aParticle->Velocity[aDimIdx]),
 198                                             aMinimalVelocity(aDimIdx + 1));
 199       }
 200
 201       myFunc->Value(aCurrPoint, aParticle->Distance);
 202       aParticle->Distance *= aParticle->Distance;
 203       if(aParticle->Distance < aParticle->BestDistance)
 204       {
 205         aParticle->BestDistance = aParticle->Distance;
 206         for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
 207           aParticle->BestPosition[aDimIdx] = aParticle->Position[aDimIdx];
 208
 209         if(aParticle->Distance < aBestGlobalDistance)
 210         {
 211           aBestGlobalDistance = aParticle->Distance;
 212           for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
 213             aBestGlobalPosition(aDimIdx + 1) = aParticle->Position[aDimIdx];
 214         }
 215       }
 216     }
 217
 218     Standard_Boolean isTerminalVelocityReached = Standard_True;
 219     for(Standard_Integer aDimIdx = 1; aDimIdx <= myN; ++aDimIdx)
 220     {
 221       if(aMinimalVelocity(aDimIdx) > aTerminationVelocity(aDimIdx))
 222       {
 223         isTerminalVelocityReached = Standard_False;
 224         break;
 225       }
 226     }
 227
 228     if(isTerminalVelocityReached)
 229     {
 230       // Minimum number of steps
 231       const Standard_Integer aMinSteps = 16;
 232
 233       if (aStep > aMinSteps)
 234       {
 235         break;
 236       }
 237
 238       for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx)
 239       {
 240         const Standard_Real aKsi = aRandom.NextReal();
 241
 242         aParticle = aParticles.GetParticle(aPartIdx);
 243
 244         for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
 245         {
 246           if (aParticle->Position[aDimIdx] == aMinUV(aDimIdx + 1) ||
 247               aParticle->Position[aDimIdx] == aMaxUV(aDimIdx + 1))
 248           {
 249             aParticle->Velocity[aDimIdx] = aParticle->Position[aDimIdx] == aMinUV(aDimIdx + 1) ?
 250               mySteps(aDimIdx + 1) * aKsi : -mySteps(aDimIdx + 1) * aKsi;
 251           }
 252           else
 253           {
 254             aParticle->Velocity[aDimIdx] = mySteps(aDimIdx + 1) * (aKsi - 0.5) * 2.0;
 255           }
 256         }
 257       }
 258     }
 259   }
 260   theValue  = aBestGlobalDistance;
 261   theOutPnt = aBestGlobalPosition;
 262 }