COPASI API  4.16.103
COptMethodPS.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2006 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <cmath>
16 
17 #include "copasi.h"
18 
19 #include "COptMethodPS.h"
20 #include "COptProblem.h"
21 #include "COptItem.h"
22 #include "COptTask.h"
23 
27 #include "utilities/CSort.h"
29 
31  COptMethod(CCopasiTask::optimization, CCopasiMethod::ParticleSwarm, pParent),
32  mIterationLimit(0),
33  mSwarmSize(0),
34  mVariance(0.0),
35  mpRandom(NULL),
36  mIteration(0),
37  mhIteration(C_INVALID_INDEX),
38  mVariableSize(0),
39  mIndividuals(),
40  mValues(),
41  mVelocities(),
42  mBestValues(),
43  mBestPositions(),
44  mpPermutation(NULL),
45  mInformants(),
46  mNumInformedMin(0),
47  mNumInformed(0),
48  mBestIndex(0),
49  mEvaluationValue(0),
50  mContinue(true)
51 {
52  addParameter("Iteration Limit", CCopasiParameter::UINT, (unsigned C_INT32) 2000);
53  addParameter("Swarm Size", CCopasiParameter::UINT, (unsigned C_INT32) 50);
54  addParameter("Std. Deviation", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-6);
55  addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
56  addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
57 
58  initObjects();
59 }
60 
62  const CCopasiContainer * pParent):
63  COptMethod(src, pParent),
64  mIterationLimit(0),
65  mSwarmSize(0),
66  mVariance(0.0),
67  mpRandom(NULL),
68  mIteration(0),
69  mhIteration(C_INVALID_INDEX),
70  mVariableSize(0),
71  mIndividuals(),
72  mValues(),
73  mVelocities(),
74  mBestValues(),
75  mBestPositions(),
76  mpPermutation(NULL),
77  mInformants(),
78  mNumInformedMin(0),
79  mNumInformed(0),
80  mBestIndex(0),
81  mEvaluationValue(0),
82  mContinue(true)
83 {initObjects();}
84 
86 {cleanup();}
87 
88 // evaluate the fitness of one individual
90 {
91  // We do not need to check whether the parametric constraints are fulfilled
92  // since the parameters are created within the bounds.
93 
94  // evaluate the fitness
96 
97  // check wheter the functional constraints are fulfilled
99  mEvaluationValue = std::numeric_limits<C_FLOAT64>::infinity();
100  else
102 
103  return mEvaluationValue;
104 }
105 
106 // move an individual
107 bool COptMethodPS::move(const size_t & index)
108 {
109  const C_FLOAT64 w = 1 / (2 * log(2.0));
110  const C_FLOAT64 c = 0.5 + log(2.0);
111 
112  bool Improved = false;
113 
114  C_FLOAT64 * pIndividual = mIndividuals[index].array();
115  C_FLOAT64 * pEnd = pIndividual + mVariableSize;
116  C_FLOAT64 * pVelocity = mVelocities[index];
117  C_FLOAT64 * pBestPosition = mBestPositions[index];
118  std::vector< COptItem * >::const_iterator itOptItem = mpOptItem->begin();
119  std::vector< UpdateMethod * >::const_iterator itSetCalculateVariable = mpSetCalculateVariable->begin();
120 
121  C_FLOAT64 * pBestInformantPosition = mBestPositions[index];
122  C_FLOAT64 BestInformantValue = mBestValues[index];
123 
124  std::set< size_t >::const_iterator itInformant = mInformants[index].begin();
125  std::set< size_t >::const_iterator endInformant = mInformants[index].end();
126 
127  size_t i = mNumInformed + mNumInformedMin;
128 
129  for (; i && itInformant != endInformant; --i, ++itInformant)
130  if (mBestValues[*itInformant] < BestInformantValue)
131  {
132  BestInformantValue = mBestValues[*itInformant];
133  pBestInformantPosition = mBestPositions[*itInformant];
134  }
135 
136  for (; pIndividual != pEnd;
137  ++pIndividual, ++pVelocity, ++pBestPosition, ++itOptItem, ++itSetCalculateVariable, ++pBestInformantPosition)
138  {
139  *pVelocity *= w;
140  *pVelocity += c * mpRandom->getRandomCC() * (*pBestPosition - *pIndividual);
141  *pVelocity += c * mpRandom->getRandomCC() * (*pBestInformantPosition - *pIndividual);
142 
143  *pIndividual += *pVelocity;
144 
145  COptItem & OptItem = **itOptItem;
146 
147  // force it to be within the bounds
148  switch (OptItem.checkConstraint(*pIndividual))
149  {
150  case - 1:
151  *pIndividual = *OptItem.getLowerBoundValue();
152  *pVelocity = 0.0;
153  break;
154 
155  case 1:
156  *pIndividual = *OptItem.getUpperBoundValue();
157  *pVelocity = 0.0;
158  break;
159  }
160 
161  // We need to set the value here so that further checks take
162  // account of the value.
163  (**itSetCalculateVariable)(*pIndividual);
164  }
165 
166  // calculate its fitness
167  mValues[index] = evaluate();
168 
169  // Check if we improved individually
170  if (mEvaluationValue < mBestValues[index])
171  {
172  Improved = true;
173 
174  // Save the individually best value;
175  mBestValues[index] = mEvaluationValue;
176  memcpy(mBestPositions[index], mIndividuals[index].array(), sizeof(C_FLOAT64) * mVariableSize);
177 
178  // Check if we improved globally
180  {
181  // and store that value
182  mBestIndex = index;
184 
185  // We found a new best value lets report it.
187  }
188  }
189 
190  return Improved;
191 }
192 
193 // initialise an individual
194 bool COptMethodPS::create(const size_t & index)
195 {
196  C_FLOAT64 * pIndividual = mIndividuals[index].array();
197  C_FLOAT64 * pEnd = pIndividual + mVariableSize;
198  C_FLOAT64 * pVelocity = mVelocities[index];
199  C_FLOAT64 * pBestPosition = mBestPositions[index];
200  std::vector< COptItem * >::const_iterator itOptItem = mpOptItem->begin();
201  std::vector< UpdateMethod * >::const_iterator itSetCalculateVariable = mpSetCalculateVariable->begin();
202 
203  C_FLOAT64 mn, mx, la;
204 
205  for (; pIndividual != pEnd;
206  ++pIndividual, ++pVelocity, ++pBestPosition, ++itOptItem, ++itSetCalculateVariable)
207  {
208  COptItem & OptItem = **itOptItem;
209 
210  // calculate lower and upper bounds
211  mn = *OptItem.getLowerBoundValue();
212  mx = *OptItem.getUpperBoundValue();
213 
214  try
215  {
216  // First determine the location of the intervall
217  // Secondly determine whether to distribute the parameter linearly or not
218  // depending on the location and act uppon it.
219  if (0.0 <= mn) // the interval [mn, mx) is in [0, inf)
220  {
221  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
222 
223  if (la < 1.8 || !(mn > 0.0)) // linear
224  {
225  *pIndividual = mn + mpRandom->getRandomCC() * (mx - mn);
226  *pVelocity = mn + mpRandom->getRandomCC() * (mx - mn) - *pIndividual;
227  }
228  else
229  {
230  *pIndividual = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC());
231  *pVelocity =
232  pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC()) - *pIndividual;
233  }
234  }
235  else if (mx > 0) // 0 is in the interval (mn, mx)
236  {
237  la = log10(mx) + log10(-mn);
238 
239  if (la < 3.6) // linear
240  {
241  *pIndividual = mn + mpRandom->getRandomCC() * (mx - mn);
242  *pVelocity = mn + mpRandom->getRandomCC() * (mx - mn) - *pIndividual;
243  }
244  else
245  {
246  C_FLOAT64 mean = (mx + mn) * 0.5;
247  C_FLOAT64 sigma = std::min(std::numeric_limits< C_FLOAT64 >::max(), mx - mn) / 3.0;
248 
249  do
250  {
251  *pIndividual = mpRandom->getRandomNormal(mean, sigma);
252  }
253  while ((*pIndividual < mn) || (*pIndividual > mx));
254 
255  *pVelocity = mpRandom->getRandomNormal(mean, sigma) - *pIndividual;
256  }
257  }
258  else // the interval (mn, mx] is in (-inf, 0]
259  {
260  // Switch lower and upper bound and change sign, i.e.,
261  // we can treat it similarly as location 1:
262  mx = - *OptItem.getLowerBoundValue();
263  mn = - *OptItem.getUpperBoundValue();
264 
265  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
266 
267  if (la < 1.8 || !(mn > 0.0)) // linear
268  {
269  *pIndividual = - (mn + mpRandom->getRandomCC() * (mx - mn));
270  *pVelocity = - (mn + mpRandom->getRandomCC() * (mx - mn)) - *pIndividual;
271  }
272  else
273  {
274  *pIndividual = - pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC());
275  *pVelocity =
276  - pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC()) - *pIndividual;
277  }
278  }
279  }
280 
281  catch (...)
282  {
283  *pIndividual = (mx + mn) * 0.5;
284  *pVelocity = 0.0;
285  }
286 
287  // force it to be within the bounds
288  switch (OptItem.checkConstraint(*pIndividual))
289  {
290  case - 1:
291  *pIndividual = *OptItem.getLowerBoundValue();
292  break;
293 
294  case 1:
295  *pIndividual = *OptItem.getUpperBoundValue();
296  break;
297  }
298 
299  *pBestPosition = *pIndividual;
300 
301  // We need to set the value here so that further checks take
302  // account of the value.
303  (**itSetCalculateVariable)(*pIndividual);
304  }
305 
306  // calculate its fitness
307  mBestValues[index] = mValues[index] = evaluate();
308 
309  if (mBestValues[index] < mBestValues[mBestIndex])
310  {
311  // and store that value
312  mBestIndex = index;
314 
315  // We found a new best value lets report it.
317  }
318 
319  return mContinue;
320 }
321 
323 {
325 }
326 
328 {
329  cleanup();
330 
331  if (!COptMethod::initialize()) return false;
332 
333  mIterationLimit = * getValue("Iteration Limit").pUINT;
334  mIteration = 0;
335 
336  if (mpCallBack)
337  mhIteration =
338  mpCallBack->addItem("Iteration Limit",
339  mIteration,
340  & mIterationLimit);
341 
342  mSwarmSize = * getValue("Swarm Size").pUINT;
343 
344  if (mSwarmSize < 5)
345  {
346  mSwarmSize = 5;
347  setValue("Swarm Size", mSwarmSize);
348  }
349 
350  mVariance = *getValue("Std. Deviation").pUDOUBLE;
351  mVariance *= mVariance;
352 
353  mpRandom =
354  CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
355  * getValue("Seed").pUINT);
356 
357  mVariableSize = mpOptItem->size();
358 
360 
361  size_t i;
362 
363  for (i = 0; i < mSwarmSize; i++)
364  mIndividuals[i].resize(mVariableSize);
365 
366  mValues.resize(mSwarmSize);
367  mVelocities.resize(mSwarmSize, mVariableSize);
368  mBestValues.resize(mSwarmSize);
369  mBestPositions.resize(mSwarmSize, mVariableSize);
370 
371  mNumInformedMin = std::max<size_t>(mSwarmSize / 10, 5) - 1;
373 
374  mpPermutation = new CPermutation(mpRandom, mSwarmSize);
375 
376  mContinue = true;
377 
378  return mContinue;
379 }
380 
382 {
383  pdelete(mpRandom);
385 
386  return true;
387 }
388 
390 {
391  if (mNumInformed < mSwarmSize)
392  mNumInformed++;
393  else
394  return;
395 
396  mInformants.clear();
397  mInformants.resize(mSwarmSize);
399 
400  size_t i, j;
401  size_t Informant;
402 
403  for (i = 0; i < mSwarmSize; i++)
404  {
405  mInformants[i].insert(i);
406 
407  Informant = mpPermutation->pick();
408 
409  for (j = 1; j < mNumInformed; j++, Informant = mpPermutation->next())
410  {
411  if (Informant == i)
412  {
413  Informant = mpPermutation->next();
414  }
415 
416  mInformants[Informant].insert(i);
417  }
418  }
419 
420  return;
421 }
422 
424 {
425  if (mNumInformed > mNumInformedMin + 1)
426  mNumInformed--;
427 
428  // Check whether the swarm has settled
429  C_FLOAT64 * pValue = mValues.array();
430  C_FLOAT64 * pEnd = pValue + mSwarmSize;
431 
432  C_FLOAT64 Delta;
433 
434  C_FLOAT64 Mean = 0.0;
435  C_FLOAT64 Variance = 0.0;
436  size_t N = 0;
437 
438  for (; pValue != pEnd; ++pValue)
439  {
440  // We need to deal with infinity values since they indicate failure
441  if (*pValue == std::numeric_limits<C_FLOAT64>::infinity())
442  return false;
443 
444  Delta = *pValue - Mean;
445  Mean += Delta / ++N;
446  // This uses the new mean, i.e., not Delta * Delta
447  Variance += Delta * (*pValue - Mean);
448  }
449 
450  Variance /= (N - 1);
451 
452  if (Variance > mVariance)
453  return false;
454 
455  // The variance of the function value is smaller than required. We now
456  // Check the variance of the flock positions.
457  CVector< C_FLOAT64 > FirstMoments(mVariableSize);
458  CVector< C_FLOAT64 > SecondMoments(mVariableSize);
459  FirstMoments = 0.0;
460  SecondMoments = 0.0;
461 
462  CVector< C_FLOAT64 > * pIndividual = mIndividuals.array();
463  CVector< C_FLOAT64 > * pIndividualEnd = pIndividual + mSwarmSize;
464 
465  C_FLOAT64 * pFirstMoment;
466  C_FLOAT64 * pSecondMoment;
467  pEnd = FirstMoments.array() + mVariableSize;
468 
469  for (; pIndividual != pIndividualEnd; ++pIndividual)
470  {
471  pFirstMoment = FirstMoments.array();
472  pSecondMoment = SecondMoments.array();
473  pValue = pIndividual->array();
474 
475  for (; pFirstMoment != pEnd; ++pFirstMoment, ++pSecondMoment, ++pValue)
476  {
477  *pFirstMoment += *pValue;
478  *pSecondMoment += *pValue * *pValue;
479  }
480  }
481 
482  pFirstMoment = FirstMoments.array();
483  pSecondMoment = SecondMoments.array();
484 
485  for (; pFirstMoment != pEnd; ++pFirstMoment, ++pSecondMoment)
486  {
487  Variance = (*pSecondMoment - *pFirstMoment * *pFirstMoment / mSwarmSize) / (mSwarmSize - 1);
488 
489  if (Variance > mVariance) return false;
490  }
491 
492  return true;
493 }
494 
496 {
497  size_t i;
498 
499  if (!initialize())
500  {
501  if (mpCallBack)
503 
504  return false;
505  }
506 
507  C_FLOAT64 * pIndividual = mIndividuals[0].array();
508  C_FLOAT64 * pEnd = pIndividual + mVariableSize;
509  C_FLOAT64 * pVelocity = mVelocities[0];
510  C_FLOAT64 * pBestPosition = mBestPositions[0];
511  std::vector< COptItem * >::const_iterator itOptItem = mpOptItem->begin();
512  std::vector< UpdateMethod * >::const_iterator itSetCalculateVariable = mpSetCalculateVariable->begin();
513 
514  // initialise the population
515  // first individual is the initial guess
516  for (; pIndividual != pEnd;
517  ++pIndividual, ++pVelocity, ++pBestPosition, ++itOptItem, ++itSetCalculateVariable)
518  {
519  COptItem & OptItem = **itOptItem;
520 
521  *pIndividual = OptItem.getStartValue();
522 
523  // force it to be within the bounds
524  switch (OptItem.checkConstraint(*pIndividual))
525  {
526  case - 1:
527  *pIndividual = *OptItem.getLowerBoundValue();
528  break;
529 
530  case 1:
531  *pIndividual = *OptItem.getUpperBoundValue();
532  break;
533  }
534 
535  *pBestPosition = *pIndividual;
536  *pVelocity = 0.0;
537 
538  // We need to set the value here so that further checks take
539  // account of the value.
540  (**itSetCalculateVariable)(*pIndividual);
541  }
542 
543  // calculate its fitness
544  mBestValues[0] = mValues[0] = evaluate();
545 
546  // and store that value
547  mBestIndex = 0;
549 
550  // We found a new best value lets report it.
552 
553  // the others are random
554  for (i = 1; i < mSwarmSize && mContinue; i++)
555  create(i);
556 
557  // create the informant list
558  buildInformants();
559 
560  bool Improved;
561 
563  {
564  Improved = false;
565 
566  for (i = 0; i < mSwarmSize && mContinue; i++)
567  Improved |= move(i);
568 
569  if (!Improved)
570  buildInformants();
571  else if (reachedStdDeviation())
572  break;
573 
574  if (mpCallBack)
575  mContinue &= mpCallBack->progressItem(mhIteration);
576  }
577 
578  if (mpCallBack)
580 
581  cleanup();
582 
583  return true;
584 }
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
unsigned C_INT32 mIteration
Definition: COptMethodPS.h:132
COptTask * mpParentTask
Definition: COptMethod.h:58
const C_FLOAT64 & evaluate()
C_FLOAT64 mEvaluationValue
Definition: COptMethodPS.h:197
#define pdelete(p)
Definition: copasi.h:215
size_t mNumInformedMin
Definition: COptMethodPS.h:182
virtual bool initialize()
Definition: COptMethod.cpp:189
bool move(const size_t &index)
size_t mVariableSize
Definition: COptMethodPS.h:142
virtual C_FLOAT64 getRandomNormal(const C_FLOAT64 &mean, const C_FLOAT64 &sd)
Definition: CRandom.cpp:292
void initObjects()
const size_t & next()
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INVALID_INDEX
Definition: copasi.h:222
COptProblem * mpOptProblem
Definition: COptMethod.h:56
virtual void output(const COutputInterface::Activity &activity)
bool create(const size_t &index)
#define C_INT32
Definition: copasi.h:90
const size_t & pick()
CMatrix< C_FLOAT64 > mBestPositions
Definition: COptMethodPS.h:167
virtual bool progressItem(const size_t &handle)
unsigned C_INT32 mIterationLimit
Definition: COptMethodPS.h:112
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
virtual bool calculate()
virtual C_FLOAT64 getRandomCC()
Definition: CRandom.cpp:235
CVector< C_FLOAT64 > mValues
Definition: COptMethodPS.h:152
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
size_t mNumInformed
Definition: COptMethodPS.h:187
const std::vector< UpdateMethod * > * mpSetCalculateVariable
Definition: COptMethod.h:65
const Value & getValue() const
void shuffle(const size_t &swaps=C_INVALID_INDEX)
bool setValue(const std::string &name, const CType &value)
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
unsigned C_INT32 * pUINT
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
CPermutation * mpPermutation
Definition: COptMethodPS.h:172
void buildInformants()
virtual bool finishItem(const size_t &handle)
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
virtual bool checkFunctionalConstraints()
bool reachedStdDeviation()
const C_FLOAT64 & getStartValue() const
Definition: COptItem.cpp:199
unsigned C_INT32 mSwarmSize
Definition: COptMethodPS.h:117
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
size_t mBestIndex
Definition: COptMethodPS.h:192
#define C_FLOAT64
Definition: copasi.h:92
std::vector< std::set< size_t > > mInformants
Definition: COptMethodPS.h:177
COptMethodPS(const COptMethodPS &src, const CCopasiContainer *pParent=NULL)
CVector< CVector< C_FLOAT64 > > mIndividuals
Definition: COptMethodPS.h:147
CType * array()
Definition: CVector.h:139
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
C_FLOAT64 mVariance
Definition: COptMethodPS.h:122
bool addParameter(const CCopasiParameter &parameter)
virtual bool optimise()
virtual ~COptMethodPS()
CRandom * mpRandom
Definition: COptMethodPS.h:127
size_t mhIteration
Definition: COptMethodPS.h:137
CMatrix< C_FLOAT64 > mVelocities
Definition: COptMethodPS.h:157
const C_FLOAT64 & getCalculateValue() const
virtual bool cleanup()
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
#define min(a, b)
Definition: f2c.h:175
virtual bool initialize()
CVector< C_FLOAT64 > mBestValues
Definition: COptMethodPS.h:162
#define max(a, b)
Definition: f2c.h:176