COPASI API  4.16.103
COptMethodEP.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) 2005 - 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 "COptMethodEP.h"
20 #include "COptProblem.h"
21 #include "COptItem.h"
22 #include "COptTask.h"
23 
26 #include "utilities/CSort.h"
27 
29  COptMethod(CCopasiTask::optimization, CCopasiMethod::EvolutionaryProgram, pParent),
30  mGenerations(0),
31  mGeneration(0),
32  mPopulationSize(0),
33  mpRandom(NULL),
34  mBestIndex(C_INVALID_INDEX),
35  mLosses(0),
36  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
37  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
38  mValue(0),
39  mVariableSize(0),
40  mIndividual(0),
41  mVariance(0)
42 {
43  addParameter("Number of Generations", CCopasiParameter::UINT, (unsigned C_INT32) 200);
44  addParameter("Population Size", CCopasiParameter::UINT, (unsigned C_INT32) 20);
45  addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
46  addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
47 
48  initObjects();
49 }
50 
52  const CCopasiContainer * pParent): COptMethod(src, pParent),
53  mGenerations(0),
54  mGeneration(0),
55  mPopulationSize(0),
56  mpRandom(NULL),
57  mBestIndex(C_INVALID_INDEX),
58  mLosses(0),
59  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
60  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
61  mValue(0),
62  mVariableSize(0),
63  mIndividual(0),
64  mVariance(0)
65 {initObjects();}
66 
68 {
69  cleanup();
70 }
71 
73 {
74  if (!initialize())
75  {
76  if (mpCallBack)
78 
79  return false;
80  }
81 
82  bool Continue = true;
83 
84  // Initialize the population
85  Continue = creation();
86 
87  // get the index of the fittest
88  mBestIndex = fittest();
89 
91  {
92  // and store that value
95 
96  // We found a new best value lets report it.
98  }
99 
100  if (!Continue)
101  {
102  if (mpCallBack)
104 
105  cleanup();
106  return true;
107  }
108 
109  // iterate over Generations
110  for (mGeneration = 2; mGeneration <= mGenerations && Continue; mGeneration++)
111  {
112  // replicate the individuals
113  Continue = replicate();
114 
115  // select the most fit
116  Continue = select();
117 
118  // get the index of the fittest
119  mBestIndex = fittest();
120 
121  if (mBestIndex != C_INVALID_INDEX &&
123  {
125 
127 
128  // We found a new best value lets report it.
129  //if (mpReport) mpReport->printBody();
131  }
132 
133  if (mpCallBack)
135  }
136 
137  if (mpCallBack)
139 
140  cleanup();
141 
142  return true;
143 }
144 
146 {
147  size_t i;
148  // pdelete all used variables
149  pdelete(mpRandom);
150 
151  for (i = 0; i < mIndividual.size(); i++)
152  {
153  pdelete(mIndividual[i]);
154  pdelete(mVariance[i]);
155  }
156 
157  return true;
158 }
159 
161 {
162  cleanup();
163 
164  size_t i;
165 
166  if (!COptMethod::initialize()) return false;
167 
168  mGenerations = * getValue("Number of Generations").pUINT;
169  mGeneration = 0;
170 
171  if (mpCallBack)
172  mhGenerations =
173  mpCallBack->addItem("Current Generation",
174  mGeneration,
175  & mGenerations);
176 
177  mGeneration++;
178 
179  mPopulationSize = * getValue("Population Size").pUINT;
180  mpRandom =
181  CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
182  * getValue("Seed").pUINT);
183 
184  mVariableSize = mpOptItem->size();
185 
186  mIndividual.resize(2 * mPopulationSize);
187  mVariance.resize(2 * mPopulationSize);
188 
189  for (i = 0; i < 2 * mPopulationSize; i++)
190  {
193  }
194 
195  mValue.resize(2 * mPopulationSize);
196  mValue = std::numeric_limits< C_FLOAT64 >::infinity();
197 
198  mLosses.resize(2 * mPopulationSize);
199  mLosses = 0;
200 
201  // Initialize the parameters to update the variances
202  tau1 = 1.0 / sqrt(2 * double(mVariableSize));
203  tau2 = 1.0 / sqrt(2 * sqrt(double(mVariableSize)));
204 
205  return true;
206 }
207 
208 // evaluate the fitness of one individual
209 bool COptMethodEP::evaluate(const CVector< C_FLOAT64 > & /* individual */)
210 {
211  bool Continue = true;
212 
213  // We do not need to check whether the parametric constraints are fulfilled
214  // since the parameters are created within the bounds.
215 
216  // evaluate the fitness
217  Continue = mpOptProblem->calculate();
218 
219  // check whether the functional constraints are fulfilled
221  mEvaluationValue = std::numeric_limits<C_FLOAT64>::infinity();
222  else
224 
225  return Continue;
226 }
227 
229 {
231 }
232 
234 {
235  size_t i;
236  size_t j;
237 
238  C_FLOAT64 mn;
239  C_FLOAT64 mx;
240  C_FLOAT64 la;
241 
242  bool Continue = true;
243 
244  // set the first individual to the initial guess
245 
246  for (i = 0; i < mVariableSize; i++)
247  {
248  C_FLOAT64 & mut = (*mIndividual[0])[i];
249  COptItem & OptItem = *(*mpOptItem)[i];
250 
251  mut = OptItem.getStartValue();
252 
253  // force it to be within the bounds
254  switch (OptItem.checkConstraint(mut))
255  {
256  case - 1:
257  mut = *OptItem.getLowerBoundValue();
258 
259  if (!OptItem.checkLowerBound(mut)) // Inequality
260  {
261  if (mut == 0.0)
263  else
264  mut += mut * std::numeric_limits< C_FLOAT64 >::epsilon();
265  }
266 
267  break;
268 
269  case 1:
270  mut = *OptItem.getUpperBoundValue();
271 
272  if (!OptItem.checkUpperBound(mut)) // Inequality
273  {
274  if (mut == 0.0)
276  else
277  mut -= mut * std::numeric_limits< C_FLOAT64 >::epsilon();
278  }
279 
280  break;
281  }
282 
283  // We need to set the value here so that further checks take
284  // account of the value.
285  (*(*mpSetCalculateVariable)[i])(mut);
286 
287  // Set the variance for this parameter.
288  (*mVariance[0])[i] = fabs(mut) * 0.5;
289  }
290 
291  Continue = evaluate(*mIndividual[0]);
293  //candx[0] = evaluate(0);
294 
295  // and copy it to the rest
296  // for(i=1; i<half; i++)
297  // candx[i] = candx[0];
298  // set the other half to random values within the boundaries
299  // for(i=half; i<mVariableSizeze; i++)
300 
301  for (i = 1; i < mPopulationSize; i++)
302  {
303  for (j = 0; j < mVariableSize; j++)
304  {
305  C_FLOAT64 & mut = (*mIndividual[i])[j];
306  COptItem & OptItem = *(*mpOptItem)[j];
307 
308  // calculate lower and upper bounds
309  mn = *OptItem.getLowerBoundValue();
310  mx = *OptItem.getUpperBoundValue();
311 
312  try
313  {
314  // First determine the location of the interval
315  // Secondly determine whether to distribute the parameter linearly or not
316  // depending on the location and act uppon it.
317  if (0.0 <= mn) // the interval [mn, mx) is in [0, inf)
318  {
319  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
320 
321  if (la < 1.8 || !(mn > 0.0)) // linear
322  mut = mn + mpRandom->getRandomCC() * (mx - mn);
323  else
324  mut = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC());
325  }
326  else if (mx > 0) // 0 is in the interval (mn, mx)
327  {
328  la = log10(mx) + log10(-mn);
329 
330  if (la < 3.6) // linear
331  mut = mn + mpRandom->getRandomCC() * (mx - mn);
332  else
333  {
334  C_FLOAT64 mean = (mx + mn) * 0.5;
335  C_FLOAT64 sigma = mean * 0.01;
336 
337  do
338  {
339  mut = mpRandom->getRandomNormal(mean, sigma);
340  }
341  while ((mut < mn) || (mut > mx));
342  }
343  }
344  else // the interval (mn, mx] is in (-inf, 0]
345  {
346  // Switch lower and upper bound and change sign, i.e.,
347  // we can treat it similarly as location 1:
348  mx = - *OptItem.getLowerBoundValue();
349  mn = - *OptItem.getUpperBoundValue();
350 
351  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
352 
353  if (la < 1.8 || !(mn > 0.0)) // linear
354  mut = - (mn + mpRandom->getRandomCC() * (mx - mn));
355  else
356  mut = - pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC());
357  }
358  }
359 
360  catch (...)
361  {
362  mut = (mx + mn) * 0.5;
363  }
364 
365  // force it to be within the bounds
366  switch (OptItem.checkConstraint(mut))
367  {
368  case - 1:
369  mut = *OptItem.getLowerBoundValue();
370 
371  if (!OptItem.checkLowerBound(mut)) // Inequality
372  {
373  if (mut == 0.0)
375  else
376  mut += mut * std::numeric_limits< C_FLOAT64 >::epsilon();
377  }
378 
379  break;
380 
381  case 1:
382  mut = *OptItem.getUpperBoundValue();
383 
384  if (!OptItem.checkUpperBound(mut)) // Inequality
385  {
386  if (mut == 0.0)
388  else
389  mut -= mut * std::numeric_limits< C_FLOAT64 >::epsilon();
390  }
391 
392  break;
393  }
394 
395  // We need to set the value here so that further checks take
396  // account of the value.
397  (*(*mpSetCalculateVariable)[j])(mut);
398 
399  // Set the variance for this parameter.
400  (*mVariance[i])[j] = fabs(mut) * 0.5;
401  }
402 
403  // calculate its fitness
404  Continue = evaluate(*mIndividual[i]);
406  }
407 
408  return Continue;
409 }
410 
412 {
413  size_t i, j, nopp, opp;
414  size_t TotalPopulation = 2 * mPopulationSize;
415 
416  // tournament competition
417  mLosses = 0; // Set all losses to 0.
418 
419  // compete with ~ 20% of the TotalPopulation
420  nopp = std::max<size_t>(1, mPopulationSize / 5);
421 
422  // parents and offspring are all in competition
423  for (i = 0; i < TotalPopulation; i++)
424  for (j = 0; j < nopp; j++)
425  {
426  // get random opponent
427  do
428  {
429  opp = mpRandom->getRandomU((unsigned C_INT32)(TotalPopulation - 1));
430  }
431  while (i == opp);
432 
433  if (mValue[i] < mValue[opp])
434  mLosses[opp]++;
435  else
436  mLosses[i]++;
437  }
438 
441  mLosses.array() + TotalPopulation,
442  mPivot);
443 
446 
447  return true;
448 }
449 
450 bool COptMethodEP::swap(size_t from, size_t to)
451 {
452  CVector< C_FLOAT64 > * pTmp = mIndividual[to];
453  mIndividual[to] = mIndividual[from];
454  mIndividual[from] = pTmp;
455 
456  pTmp = mVariance[to];
457  mVariance[to] = mVariance[from];
458  mVariance[from] = pTmp;
459 
460  C_FLOAT64 dTmp = mValue[to];
461  mValue[to] = mValue[from];
462  mValue[from] = dTmp;
463 
464  size_t iTmp = mLosses[to];
465  mLosses[to] = mLosses[from];
466  mLosses[from] = iTmp;
467 
468  return true;
469 }
470 
472 {
473  size_t i, BestIndex = 0;
474  C_FLOAT64 BestValue = mValue[0];
475 
476  for (i = 1; i < mPopulationSize && !mLosses[i]; i++)
477  if (mValue[i] < BestValue)
478  {
479  BestIndex = i;
480  BestValue = mValue[i];
481  }
482 
483  return BestIndex;
484 }
485 
487 {
488  size_t i;
489  size_t j;
490  bool Continue = true;
491 
492  // iterate over parents
493  for (i = 0; i < mPopulationSize && Continue; i++)
494  {
495  // replicate them
496  for (j = 0; j < mVariableSize; j++)
497  {
498  (*mIndividual[mPopulationSize + i])[j] = (*mIndividual[i])[j];
499  (*mVariance[mPopulationSize + i])[j] = (*mVariance[i])[j];
500  }
501 
502  mValue[mPopulationSize + i] = mValue[i];
503 
504  // possibly mutate the offspring
505  Continue = mutate(mPopulationSize + i);
506  }
507 
508  return Continue;
509 }
510 
511 bool COptMethodEP::mutate(size_t i)
512 {
513  size_t j;
514  C_FLOAT64 v1;
515 
516  CVector<C_FLOAT64> & Individual = *mIndividual[i];
517  CVector<C_FLOAT64> & Variance = *mVariance[i];
518 
519  v1 = mpRandom->getRandomNormal01();
520 
521  // update the variances
522  for (j = 0; j < mVariableSize; j++)
523  {
524  C_FLOAT64 & mut = Individual[j];
525  COptItem & OptItem = *(*mpOptItem)[j];
526 
527  try
528  {
529  // update the parameter for the variances
530  Variance[j] =
531  std::max(Variance[j] * exp(tau1 * v1 + tau2 * mpRandom->getRandomNormal01()), 1e-8);
532 
533  // calculate the mutated parameter
534  mut += Variance[j] * mpRandom->getRandomNormal01();
535  }
536 
537  catch (...)
538  {
539  mut = (*OptItem.getUpperBoundValue() + *OptItem.getLowerBoundValue()) * 0.5;
540  }
541 
542  // force it to be within the bounds
543  switch (OptItem.checkConstraint(mut))
544  {
545  case - 1:
546  mut = *OptItem.getLowerBoundValue();
547  break;
548 
549  case 1:
550  mut = *OptItem.getUpperBoundValue();
551  break;
552  }
553 
554  // We need to set the value here so that further checks take
555  // account of the value.
556  (*(*mpSetCalculateVariable)[j])(mut);
557  }
558 
559  // calculate its fitness
560  bool Continue = evaluate(Individual);
562 
563  return Continue;
564 }
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
COptTask * mpParentTask
Definition: COptMethod.h:58
#define pdelete(p)
Definition: copasi.h:215
virtual bool initialize()
Definition: COptMethod.cpp:189
size_t fittest()
bool checkLowerBound(const C_FLOAT64 &value) const
Definition: COptItem.cpp:435
virtual C_FLOAT64 getRandomNormal(const C_FLOAT64 &mean, const C_FLOAT64 &sd)
Definition: CRandom.cpp:292
C_FLOAT64 mEvaluationValue
Definition: COptMethodEP.h:170
size_t mVariableSize
Definition: COptMethodEP.h:180
unsigned C_INT32 mPopulationSize
Definition: COptMethodEP.h:144
CRandom * mpRandom
Definition: COptMethodEP.h:149
virtual bool cleanup()
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
size_t mhGenerations
Definition: COptMethodEP.h:139
#define C_INVALID_INDEX
Definition: copasi.h:222
COptProblem * mpOptProblem
Definition: COptMethod.h:56
virtual void output(const COutputInterface::Activity &activity)
bool applyPartialPivot(const CVector< size_t > &pivot, const size_t &ordered, SwapMethod swap)
Definition: CSort.h:375
virtual bool optimise()
#define C_INT32
Definition: copasi.h:90
CVector< size_t > mPivot
Definition: COptMethodEP.h:200
virtual bool progressItem(const size_t &handle)
virtual bool initialize()
CVector< size_t > mLosses
Definition: COptMethodEP.h:160
bool checkUpperBound(const C_FLOAT64 &value) const
Definition: COptItem.cpp:438
virtual C_FLOAT64 getRandomNormal01()
Definition: CRandom.cpp:261
bool swap(size_t from, size_t to)
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
virtual bool calculate()
bool evaluate(const CVector< C_FLOAT64 > &individual)
virtual C_FLOAT64 getRandomCC()
Definition: CRandom.cpp:235
void partialSortWithPivot(RandomAccessIterator first, RandomAccessIterator middle, RandomAccessIterator last, CVector< size_t > &pivot)
Definition: CSort.h:149
virtual unsigned C_INT32 getRandomU()
Definition: CRandom.cpp:173
virtual ~COptMethodEP()
COptMethodEP(const CCopasiContainer *pParent=NULL)
bool mutate(size_t i)
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
const Value & getValue() const
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)
std::vector< CVector< C_FLOAT64 > * > mVariance
Definition: COptMethodEP.h:190
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
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
std::vector< CVector< C_FLOAT64 > * > mIndividual
Definition: COptMethodEP.h:185
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
bool addParameter(const CCopasiParameter &parameter)
C_FLOAT64 mBestValue
Definition: COptMethodEP.h:165
void initObjects()
size_t mBestIndex
Definition: COptMethodEP.h:151
unsigned C_INT32 mGeneration
Definition: COptMethodEP.h:134
CVector< C_FLOAT64 > mValue
Definition: COptMethodEP.h:175
const C_FLOAT64 & getCalculateValue() const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
unsigned C_INT32 mGenerations
Definition: COptMethodEP.h:129
#define min(a, b)
Definition: f2c.h:175
#define max(a, b)
Definition: f2c.h:176