COPASI API  4.16.103
COptMethodDE.cpp
Go to the documentation of this file.
1 // Copyright (C) 2012 - 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 #include <limits>
7 #include <string>
8 #include <cmath>
9 
10 #include "copasi.h"
11 #include "COptMethodDE.h"
12 #include "COptProblem.h"
13 #include "COptItem.h"
14 #include "COptTask.h"
15 
19 #include "utilities/CSort.h"
21 
23  COptMethod(CCopasiTask::optimization, CCopasiMethod::DifferentialEvolution, pParent),
24  mGenerations(0),
25  mPopulationSize(0),
26  mpRandom(NULL),
27  mVariableSize(0),
28  mIndividual(0),
29  mpPermutation(NULL),
30  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
31  mValue(0),
32  mMutationVarians(0.1),
33  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
34  mBestIndex(C_INVALID_INDEX),
35  mGeneration(0)
36 
37 {
38  addParameter("Number of Generations", CCopasiParameter::UINT, (unsigned C_INT32) 2000);
39  addParameter("Population Size", CCopasiParameter::UINT, (unsigned C_INT32) 10);
40  addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
41  addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
42 
43  initObjects();
44 }
45 
47  const CCopasiContainer * pParent):
48  COptMethod(src, pParent),
49  mGenerations(0),
50  mPopulationSize(0),
51  mpRandom(NULL),
52  mVariableSize(0),
53  mIndividual(0),
54  mpPermutation(NULL),
55  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
56  mValue(0),
57  mMutationVarians(0.1),
58  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
59  mBestIndex(C_INVALID_INDEX),
60  mGeneration(0)
61 {initObjects();}
62 
64 {cleanup();}
65 
66 // evaluate the fitness of one individual
67 bool COptMethodDE::evaluate(const CVector< C_FLOAT64 > & /* individual */)
68 {
69  bool Continue = true;
70 
71  // We do not need to check whether the parametric constraints are fulfilled
72  // since the parameters are created within the bounds.
73 
74  // evaluate the fitness
75  Continue &= mpOptProblem->calculate();
76 
77  // check whether the functional constraints are fulfilled
79  mEvaluationValue = std::numeric_limits<C_FLOAT64>::infinity();
80  else
82 
83  return Continue;
84 }
85 
87 {
88  size_t i, j;
89  size_t a, b, c;
90 
91  bool Continue = true;
92 
93  for (i = mPopulationSize; i < 2 * mPopulationSize && Continue; i++)
94  {
96 
97  // MUTATION a, b, c in [0, mPopulationSize - 1]
98  a = mpPermutation->pick();
99  // b is guaranteed to be different from a
100  b = mpPermutation->next();
101  // c is guaranteed to be different from a and b
102  c = mpPermutation->next();
103 
104  // MUTATE CURRENT GENERATION
105  for (j = 0; j < mVariableSize; j++)
106  {
107  COptItem & OptItem = *(*mpOptItem)[j];
108  C_FLOAT64 & mut = (*mIndividual[i])[j];
109 
110  mut = (*mIndividual[c])[j] + 2 * ((*mIndividual[a])[j] - (*mIndividual[b])[j]);
111 
112  // force it to be within the bounds
113  switch (OptItem.checkConstraint(mut))
114  {
115  case - 1:
116  mut = *OptItem.getLowerBoundValue();
117  break;
118 
119  case 1:
120  mut = *OptItem.getUpperBoundValue();
121  break;
122  }
123 
124  // We need to set the value here so that further checks take
125  // account of the value.
126  (*(*mpSetCalculateVariable)[j])(mut);
127  }
128 
129  Continue &= evaluate(*mIndividual[i]);
131  }
132 
133  //CROSSOVER MUTATED GENERATION WITH THE CURRENT ONE
134  for (i = 2 * mPopulationSize; i < 3 * mPopulationSize && Continue; i++)
135  {
136  for (j = 0; j < mVariableSize; j++)
137  {
138 
139  COptItem & OptItem = *(*mpOptItem)[j];
140  C_FLOAT64 & mut = (*mIndividual[i])[j];
141 
142  size_t r = mpRandom->getRandomU(mPopulationSize - 1);
143 
144  if (r < 0.6 * mPopulationSize)
145  {
146  mut = (*mIndividual[i - mPopulationSize])[j] *
148  }
149 
150  else
151  mut = (*mIndividual[i - 2 * mPopulationSize])[j];
152 
153  // force it to be within the bounds
154  switch (OptItem.checkConstraint(mut))
155  {
156  case - 1:
157  mut = *OptItem.getLowerBoundValue();
158  break;
159 
160  case 1:
161  mut = *OptItem.getUpperBoundValue();
162  break;
163  }
164 
165  (*(*mpSetCalculateVariable)[j])(mut);
166  }
167 
168  Continue &= evaluate(*mIndividual[i]);
170  }
171 
172  //SELECT NEXT GENERATION
173  for (i = 2 * mPopulationSize; i < 3 * mPopulationSize && Continue; i++)
174  {
175  if (mValue[i - mPopulationSize] > mValue[i] && mValue[i - 2 * mPopulationSize] > mValue[i])
176  {
177  *mIndividual[i - 2 * mPopulationSize] = *mIndividual[i];
178  mValue[i - 2 * mPopulationSize] = mValue[i];
179  }
180  else if (mValue[i - mPopulationSize] < mValue[i - 2 * mPopulationSize])
181  {
184  }
185  else if (mBestIndex != i && mBestIndex == C_INVALID_INDEX)
186  {
187  for (j = 0; j < mVariableSize; j++)
188  {
189  COptItem & OptItem = *(*mpOptItem)[j];
190  C_FLOAT64 & mut = (*mIndividual[i - 2 * mPopulationSize])[j];
191 
192  size_t r = mpRandom->getRandomU(mPopulationSize - 1);
193 
194  if (r < 0.6 * mPopulationSize)
196 
197  // force it to be within the bounds
198  switch (OptItem.checkConstraint(mut))
199  {
200  case - 1:
201  mut = *OptItem.getLowerBoundValue();
202  break;
203 
204  case 1:
205  mut = *OptItem.getUpperBoundValue();
206  break;
207  }
208 
209  (*(*mpSetCalculateVariable)[j])(mut);
210  }
211 
212  Continue &= evaluate(*mIndividual[i - 2 * mPopulationSize]);
214  }
215  }
216 
217  return Continue;
218 }
219 
220 // check the best individual at this generation
222 {
223  size_t i, BestIndex = C_INVALID_INDEX;
225 
226  for (i = 0; i < mPopulationSize; i++)
227  if (mValue[i] < BestValue)
228  {
229  BestIndex = i;
230  BestValue = mValue[i];
231  }
232 
233  return BestIndex;
234 }
235 
236 // check the best individual at this generation
238 {
239  size_t i;
240 
241  for (i = 2 * mPopulationSize; i < 3 * mPopulationSize; i++)
242  if (mValue[i] < mValue[i - 2 * mPopulationSize])
243  {
244  *mIndividual[i - 2 * mPopulationSize] = *mIndividual[i];
245  mValue[i - 2 * mPopulationSize] = mValue[i];
246  }
247 }
248 
249 // Initialize the population
250 bool COptMethodDE::creation(size_t first, size_t last)
251 {
252  size_t Last = std::min(last, (size_t) mPopulationSize);
253  size_t i, j;
254 
255  C_FLOAT64 mn;
256  C_FLOAT64 mx;
257  C_FLOAT64 la;
258 
259  bool Continue = true;
260 
261  for (i = first; i < Last && Continue; i++)
262  {
263  // We do not want to loose the best individual;
264  if (mBestIndex != i)
265  for (j = 0; j < mVariableSize; j++)
266  {
267  // calculate lower and upper bounds
268  COptItem & OptItem = *(*mpOptItem)[j];
269  mn = *OptItem.getLowerBoundValue();
270  mx = *OptItem.getUpperBoundValue();
271 
272  C_FLOAT64 & mut = (*mIndividual[i])[j];
273 
274  try
275  {
276  // determine if linear or log scale
277  if ((mn < 0.0) || (mx <= 0.0))
278  mut = mn + mpRandom->getRandomCC() * (mx - mn);
279  else
280  {
281  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
282 
283  if (la < 1.8)
284  mut = mn + mpRandom->getRandomCC() * (mx - mn);
285  else
286  mut = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC());
287  }
288  }
289 
290  catch (...)
291  {
292  mut = (mx + mn) * 0.5;
293  }
294 
295  // force it to be within the bounds
296  switch (OptItem.checkConstraint(mut))
297  {
298  case - 1:
299  mut = *OptItem.getLowerBoundValue();
300  break;
301 
302  case 1:
303  mut = *OptItem.getUpperBoundValue();
304  break;
305  }
306 
307  // We need to set the value here so that further checks take
308  // account of the value.
309  (*(*mpSetCalculateVariable)[j])(mut);
310  }
311 
312  // calculate its fitness
313  Continue &= evaluate(*mIndividual[i]);
315  }
316 
317  return Continue;
318 }
319 
321 {
323 }
324 
326 {
327  cleanup();
328 
329  size_t i;
330 
331  if (!COptMethod::initialize())
332  {
333  if (mpCallBack)
335 
336  return false;
337  }
338 
339  mGenerations = * getValue("Number of Generations").pUINT;
340  mGeneration = 0;
341 
342  if (mpCallBack)
343  mhGenerations =
344  mpCallBack->addItem("Current Generation",
345  mGeneration,
346  & mGenerations);
347 
348  mGeneration++;
349 
350  mPopulationSize = * getValue("Population Size").pUINT;
351 
352  if (mPopulationSize < 4)
353  {
354  mPopulationSize = 4;
355  setValue("Population Size", mPopulationSize);
356  }
357 
358  mpRandom =
359  CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
360  * getValue("Seed").pUINT);
361 
363 
364  mVariableSize = mpOptItem->size();
365 
366  mIndividual.resize(3 * mPopulationSize);
367 
368  for (i = 0; i < 3 * mPopulationSize; i++)
370 
371  mValue.resize(3 * mPopulationSize);
372  mValue = std::numeric_limits<C_FLOAT64>::infinity();
373  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
374 
375  mMutationVarians = 0.1;
376 
377  return true;
378 }
379 
381 {
382  size_t i;
383 
384  pdelete(mpRandom);
386 
387  for (i = 0; i < mIndividual.size(); i++)
388  pdelete(mIndividual[i]);
389 
390  return true;
391 }
392 
394 {
395  bool Continue = true;
396 
397  if (!initialize())
398  {
399  if (mpCallBack)
401 
402  return false;
403  }
404 
405  size_t i;
406 
407  // initialise the population
408  // first individual is the initial guess
409  for (i = 0; i < mVariableSize; i++)
410  {
411  C_FLOAT64 & mut = (*mIndividual[0])[i];
412  COptItem & OptItem = *(*mpOptItem)[i];
413 
414  mut = OptItem.getStartValue();
415 
416  // force it to be within the bounds
417  switch (OptItem.checkConstraint(mut))
418  {
419  case - 1:
420  mut = *OptItem.getLowerBoundValue();
421  break;
422 
423  case 1:
424  mut = *OptItem.getUpperBoundValue();
425  break;
426  }
427 
428  // We need to set the value here so that further checks take
429  // account of the value.
430  (*(*mpSetCalculateVariable)[i])(mut);
431  }
432 
433  Continue &= evaluate(*mIndividual[0]);
435 
436  if (!isnan(mEvaluationValue))
437  {
438  // and store that value
439  mBestValue = mValue[0];
440  Continue &= mpOptProblem->setSolution(mBestValue, *mIndividual[0]);
441 
442  // We found a new best value lets report it.
444  }
445 
446  // the others are random
447  Continue &= creation(1, mPopulationSize);
448 
449  mBestIndex = fittest();
450 
451  if (mBestIndex != C_INVALID_INDEX &&
453  {
454  // and store that value
457 
458  // We found a new best value lets report it.
460  }
461 
462  if (!Continue)
463  {
464  if (mpCallBack)
466 
467  cleanup();
468  return true;
469  }
470 
471  size_t Stalled = 0;
472 
473  // ITERATE FOR gener GENERATIONS
474  for (mGeneration = 2;
475  mGeneration <= mGenerations && Continue;
476  mGeneration++, Stalled++)
477  {
478  if (Stalled > 10)
479  {
480  Continue &= creation((size_t) 0.4 * mPopulationSize, (size_t) 0.8 * mPopulationSize);
481  }
482 
483  // perturb the population a bit
484  Continue &= creation((size_t)(mPopulationSize * 0.9), mPopulationSize);
485 
486  Continue &= replicate();
487 
488  // get the index of the fittest
489  mBestIndex = fittest();
490 
491  if (mBestIndex != C_INVALID_INDEX &&
493  {
494  Stalled = 0;
496 
498 
499  // We found a new best value lets report it.
500  //if (mpReport) mpReport->printBody();
502  }
503 
504  if (mpCallBack)
505  Continue &= mpCallBack->progressItem(mhGenerations);
506  }
507 
508  if (mpCallBack)
510 
511  cleanup();
512 
513  return true;
514 }
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
COptTask * mpParentTask
Definition: COptMethod.h:58
virtual bool optimise()
virtual bool cleanup()
#define pdelete(p)
Definition: copasi.h:215
virtual bool initialize()
Definition: COptMethod.cpp:189
unsigned C_INT32 mPopulationSize
Definition: COptMethodDE.h:114
virtual C_FLOAT64 getRandomNormal(const C_FLOAT64 &mean, const C_FLOAT64 &sd)
Definition: CRandom.cpp:292
const size_t & next()
CPermutation * mpPermutation
Definition: COptMethodDE.h:134
size_t mBestIndex
Definition: COptMethodDE.h:152
C_FLOAT64 mBestValue
Definition: COptMethodDE.h:151
bool evaluate(const CVector< C_FLOAT64 > &individual)
bool creation(size_t first, size_t last=std::numeric_limits< size_t >::max())
virtual ~COptMethodDE()
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)
#define C_INT32
Definition: copasi.h:90
const size_t & pick()
virtual bool progressItem(const size_t &handle)
size_t mVariableSize
Definition: COptMethodDE.h:124
unsigned C_INT32 mGenerations
Definition: COptMethodDE.h:104
size_t mhGenerations
Definition: COptMethodDE.h:109
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
virtual unsigned C_INT32 getRandomU()
Definition: CRandom.cpp:173
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
bool replicate()
C_FLOAT64 mEvaluationValue
Definition: COptMethodDE.h:139
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
virtual bool initialize()
unsigned C_INT32 * pUINT
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
virtual bool finishItem(const size_t &handle)
virtual bool checkFunctionalConstraints()
const C_FLOAT64 & getStartValue() const
Definition: COptItem.cpp:199
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
void initObjects()
#define C_FLOAT64
Definition: copasi.h:92
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
bool addParameter(const CCopasiParameter &parameter)
CRandom * mpRandom
Definition: COptMethodDE.h:119
std::vector< CVector< C_FLOAT64 > * > mIndividual
Definition: COptMethodDE.h:129
C_FLOAT64 mMutationVarians
Definition: COptMethodDE.h:149
const C_FLOAT64 & getCalculateValue() const
unsigned C_INT32 mGeneration
Definition: COptMethodDE.h:153
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
COptMethodDE(const COptMethodDE &src, const CCopasiContainer *pParent=NULL)
#define min(a, b)
Definition: f2c.h:175
size_t fittest()
CVector< C_FLOAT64 > mValue
Definition: COptMethodDE.h:144
#define max(a, b)
Definition: f2c.h:176