1d33d22b |
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 | |
f82a9555 |
88 | if (aCurrValue < aParticle->Distance) |
1d33d22b |
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 | } |
f82a9555 |
96 | aParticle->Distance = aCurrValue; |
97 | aParticle->BestDistance = aCurrValue; |
1d33d22b |
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); |
1d33d22b |
202 | if(aParticle->Distance < aParticle->BestDistance) |
203 | { |
204 | aParticle->BestDistance = aParticle->Distance; |
205 | for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx) |
206 | aParticle->BestPosition[aDimIdx] = aParticle->Position[aDimIdx]; |
207 | |
208 | if(aParticle->Distance < aBestGlobalDistance) |
209 | { |
210 | aBestGlobalDistance = aParticle->Distance; |
211 | for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx) |
212 | aBestGlobalPosition(aDimIdx + 1) = aParticle->Position[aDimIdx]; |
213 | } |
214 | } |
215 | } |
216 | |
217 | Standard_Boolean isTerminalVelocityReached = Standard_True; |
218 | for(Standard_Integer aDimIdx = 1; aDimIdx <= myN; ++aDimIdx) |
219 | { |
220 | if(aMinimalVelocity(aDimIdx) > aTerminationVelocity(aDimIdx)) |
221 | { |
222 | isTerminalVelocityReached = Standard_False; |
223 | break; |
224 | } |
225 | } |
226 | |
227 | if(isTerminalVelocityReached) |
228 | { |
229 | // Minimum number of steps |
230 | const Standard_Integer aMinSteps = 16; |
231 | |
232 | if (aStep > aMinSteps) |
233 | { |
234 | break; |
235 | } |
236 | |
237 | for (Standard_Integer aPartIdx = 1; aPartIdx <= myNbParticles; ++aPartIdx) |
238 | { |
239 | const Standard_Real aKsi = aRandom.NextReal(); |
240 | |
241 | aParticle = aParticles.GetParticle(aPartIdx); |
242 | |
243 | for(Standard_Integer aDimIdx = 0; aDimIdx < myN; ++aDimIdx) |
244 | { |
245 | if (aParticle->Position[aDimIdx] == aMinUV(aDimIdx + 1) || |
246 | aParticle->Position[aDimIdx] == aMaxUV(aDimIdx + 1)) |
247 | { |
248 | aParticle->Velocity[aDimIdx] = aParticle->Position[aDimIdx] == aMinUV(aDimIdx + 1) ? |
249 | mySteps(aDimIdx + 1) * aKsi : -mySteps(aDimIdx + 1) * aKsi; |
250 | } |
251 | else |
252 | { |
253 | aParticle->Velocity[aDimIdx] = mySteps(aDimIdx + 1) * (aKsi - 0.5) * 2.0; |
254 | } |
255 | } |
256 | } |
257 | } |
258 | } |
259 | theValue = aBestGlobalDistance; |
260 | theOutPnt = aBestGlobalPosition; |
261 | } |