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