1 // Created on: 2014-07-18
2 // Created by: Alexander Malyshev
3 // Copyright (c) 2014-2014 OPEN CASCADE SAS
5 // This file is part of Open CASCADE Technology software library.
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.
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
16 #include <math_PSO.hxx>
17 #include <math_PSOParticlesPool.hxx>
19 const Standard_Real aBorderDivisor = 1.0e+4;
21 //=======================================================================
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())
35 myN = theFunc->NbVariables();
36 myNbParticles = theNbParticles;
40 myLowBorder = theLowBorder;
41 myUppBorder = theUppBorder;
45 //=======================================================================
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)
55 performPSOWithGivenParticles(theParticles, theNbParticles, theValue, theOutPnt, theNbIter);
58 //=======================================================================
61 //=======================================================================
62 void math_PSO::Perform(const math_Vector& theSteps,
63 Standard_Real& theValue,
64 math_Vector& theOutPnt,
65 const Standard_Integer theNbIter)
68 math_Vector aMinUV(1, myN), aMaxUV(1, myN);
69 aMinUV = myLowBorder + (myUppBorder - myLowBorder) / aBorderDivisor;
70 aMaxUV = myUppBorder - (myUppBorder - myLowBorder) / aBorderDivisor;
74 // To generate initial distribution it is necessary to have grid steps.
75 math_PSOParticlesPool aPool(myNbParticles, myN);
77 // Generate initial particles distribution.
78 Standard_Boolean isRegularGridFinished = Standard_False;
79 Standard_Real aCurrValue;
80 math_Vector aCurrPoint(1, myN);
82 PSO_Particle* aParticle = aPool.GetWorstParticle();
86 myFunc->Value(aCurrPoint, aCurrValue);
88 if (aCurrValue < aParticle->Distance)
90 Standard_Integer aDimIdx;
91 for(aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
93 aParticle->Position[aDimIdx] = aCurrPoint(aDimIdx + 1);
94 aParticle->BestPosition[aDimIdx] = aCurrPoint(aDimIdx + 1);
96 aParticle->Distance = aCurrValue;
97 aParticle->BestDistance = aCurrValue;
99 aParticle = aPool.GetWorstParticle();
103 aCurrPoint(1) += Max(mySteps(1), 1.0e-15); // Avoid too small step
104 for(Standard_Integer aDimIdx = 1; aDimIdx < myN; ++aDimIdx)
106 if(aCurrPoint(aDimIdx) > aMaxUV(aDimIdx))
108 aCurrPoint(aDimIdx) = aMinUV(aDimIdx);
109 aCurrPoint(aDimIdx + 1) += mySteps(aDimIdx + 1);
116 if(aCurrPoint(myN) > aMaxUV(myN))
117 isRegularGridFinished = Standard_True;
119 while(!isRegularGridFinished);
121 performPSOWithGivenParticles(aPool, myNbParticles, theValue, theOutPnt, theNbIter);
124 //=======================================================================
125 //function : math_PSO
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)
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;
141 // Current particle data.
142 math_Vector aCurrPoint(1, myN);
143 math_Vector aBestGlobalPosition(1, myN);
144 PSO_Particle* aParticle = 0;
147 // Generate initial particle velocities.
148 math_BullardGenerator aRandom;
149 for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx)
151 aParticle = aParticles.GetParticle(aPartIdx);
152 for(Standard_Integer aDimIdx = 1; aDimIdx <= myN; ++aDimIdx)
154 const Standard_Real aKsi = aRandom.NextReal();
155 aParticle->Velocity[aDimIdx - 1] = mySteps(aDimIdx) * (aKsi - 0.5) * 2.0;
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;
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);
170 // Run PSO iterations
171 for (Standard_Integer aStep = 1; aStep < myNbIter; ++aStep)
173 aMinimalVelocity.Init(RealLast());
175 for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx)
177 const Standard_Real aKsi1 = aRandom.NextReal();
178 const Standard_Real aKsi2 = aRandom.NextReal();
180 aParticle = aParticles.GetParticle(aPartIdx);
182 const Standard_Real aRetentWeight = 0.72900;
183 const Standard_Real aPersonWeight = 1.49445;
184 const Standard_Real aSocialWeight = 1.49445;
186 for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
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);
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];
197 aMinimalVelocity(aDimIdx + 1) = Min(Abs(aParticle->Velocity[aDimIdx]),
198 aMinimalVelocity(aDimIdx + 1));
201 myFunc->Value(aCurrPoint, aParticle->Distance);
202 if(aParticle->Distance < aParticle->BestDistance)
204 aParticle->BestDistance = aParticle->Distance;
205 for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
206 aParticle->BestPosition[aDimIdx] = aParticle->Position[aDimIdx];
208 if(aParticle->Distance < aBestGlobalDistance)
210 aBestGlobalDistance = aParticle->Distance;
211 for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
212 aBestGlobalPosition(aDimIdx + 1) = aParticle->Position[aDimIdx];
217 Standard_Boolean isTerminalVelocityReached = Standard_True;
218 for(Standard_Integer aDimIdx = 1; aDimIdx <= myN; ++aDimIdx)
220 if(aMinimalVelocity(aDimIdx) > aTerminationVelocity(aDimIdx))
222 isTerminalVelocityReached = Standard_False;
227 if(isTerminalVelocityReached)
229 // Minimum number of steps
230 const Standard_Integer aMinSteps = 16;
232 if (aStep > aMinSteps)
237 for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx)
239 const Standard_Real aKsi = aRandom.NextReal();
241 aParticle = aParticles.GetParticle(aPartIdx);
243 for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx)
245 if (aParticle->Position[aDimIdx] == aMinUV(aDimIdx + 1) ||
246 aParticle->Position[aDimIdx] == aMaxUV(aDimIdx + 1))
248 aParticle->Velocity[aDimIdx] = aParticle->Position[aDimIdx] == aMinUV(aDimIdx + 1) ?
249 mySteps(aDimIdx + 1) * aKsi : -mySteps(aDimIdx + 1) * aKsi;
253 aParticle->Velocity[aDimIdx] = mySteps(aDimIdx + 1) * (aKsi - 0.5) * 2.0;
259 theValue = aBestGlobalDistance;
260 theOutPnt = aBestGlobalPosition;