COPASI API  4.16.103
COptMethodGA.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) 2004 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <cmath>
16 #include <limits>
17 #include <string>
18 #ifdef DEBUG_OPT
19 # include <iostream>
20 # include <fstream>
21 #endif
22 
23 #include "copasi.h"
24 
25 #include "COptMethodGA.h"
26 #include "COptProblem.h"
27 #include "COptItem.h"
28 #include "COptTask.h"
29 
33 #include "utilities/CSort.h"
35 
37  COptMethod(CCopasiTask::optimization, CCopasiMethod::GeneticAlgorithm, pParent),
38  mGenerations(0),
39  mPopulationSize(0),
40  mpRandom(NULL),
41  mVariableSize(0),
42  mIndividual(0),
43  mCrossOverFalse(0),
44  mCrossOver(0),
45  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
46  mValue(0),
47  mpPermutation(NULL),
48  mLosses(0),
49  mPivot(0),
50  mMutationVarians(0.1),
51  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
52  mBestIndex(C_INVALID_INDEX),
53  mGeneration(0)
54 
55 {
56  addParameter("Number of Generations", CCopasiParameter::UINT, (unsigned C_INT32) 200);
57  addParameter("Population Size", CCopasiParameter::UINT, (unsigned C_INT32) 20);
58  addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
59  addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
60 
61  initObjects();
62 }
63 
65  const CCopasiContainer * pParent):
66  COptMethod(src, pParent),
67  mGenerations(0),
68  mPopulationSize(0),
69  mpRandom(NULL),
70  mVariableSize(0),
71  mIndividual(0),
72  mCrossOverFalse(0),
73  mCrossOver(0),
74  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
75  mValue(0),
76  mpPermutation(NULL),
77  mLosses(0),
78  mPivot(0),
79  mMutationVarians(0.1),
80  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
81  mBestIndex(C_INVALID_INDEX),
82  mGeneration(0)
83 {initObjects();}
84 
86 {cleanup();}
87 
88 // evaluate the fitness of one individual
89 bool COptMethodGA::evaluate(const CVector< C_FLOAT64 > & /* individual */)
90 {
91  bool Continue = true;
92 
93  // We do not need to check whether the parametric constraints are fulfilled
94  // since the parameters are created within the bounds.
95 
96  // evaluate the fitness
97  Continue &= mpOptProblem->calculate();
98 
99  // check whether the functional constraints are fulfilled
101  mEvaluationValue = std::numeric_limits<C_FLOAT64>::infinity();
102  else
104 
105  return Continue;
106 }
107 
108 bool COptMethodGA::swap(size_t from, size_t to)
109 {
110  CVector< C_FLOAT64 > * pTmp = mIndividual[to];
111  mIndividual[to] = mIndividual[from];
112  mIndividual[from] = pTmp;
113 
114  C_FLOAT64 dTmp = mValue[to];
115  mValue[to] = mValue[from];
116  mValue[from] = dTmp;
117 
118  size_t iTmp = mLosses[to];
119  mLosses[to] = mLosses[from];
120  mLosses[from] = iTmp;
121 
122  return true;
123 }
124 
125 //mutate one individual
127 {
128  size_t j;
129 
130  // mutate the parameters
131  for (j = 0; j < mVariableSize; j++)
132  {
133  COptItem & OptItem = *(*mpOptItem)[j];
134  C_FLOAT64 & mut = individual[j];
135 
136  // calculate the mutatated parameter
138 
139  // force it to be within the bounds
140  switch (OptItem.checkConstraint(mut))
141  {
142  case - 1:
143  mut = *OptItem.getLowerBoundValue();
144  break;
145 
146  case 1:
147  mut = *OptItem.getUpperBoundValue();
148  break;
149  }
150 
151  // We need to set the value here so that further checks take
152  // account of the value.
153  (*(*mpSetCalculateVariable)[j])(mut);
154  }
155 
156  return true;
157 }
158 
160  const CVector< C_FLOAT64 > & parent2,
161  CVector< C_FLOAT64 > & child1,
162  CVector< C_FLOAT64 > & child2)
163 {
164  size_t i, crp;
165  size_t nCross = 0;
166 
168 
169  if (mVariableSize > 1)
170  nCross = mpRandom->getRandomU((unsigned C_INT32)(mVariableSize / 2));
171 
172  if (nCross == 0)
173  {
174  // if less than 0 just copy parent to child
175  child1 = parent1;
176  child2 = parent2;
177 
178  return true;
179  }
180 
181  // choose cross over points;
182  // We do not mind if a crossover point gets drawn twice
183  for (i = 0; i < nCross; i++)
184  {
185  crp = mpRandom->getRandomU((unsigned C_INT32)(mVariableSize - 1));
186  mCrossOver[crp] = true;
187  }
188 
189  const CVector< C_FLOAT64 > * pParent1 = & parent1;
190 
191  const CVector< C_FLOAT64 > * pParent2 = & parent2;
192 
193  const CVector< C_FLOAT64 > * pTmp;
194 
195  for (i = 0; i < mVariableSize; i++)
196  {
197  if (mCrossOver[i])
198  {
199  pTmp = pParent1;
200  pParent1 = pParent2;
201  pParent2 = pTmp;
202  }
203 
204  child1[i] = (*pParent1)[i];
205  child2[i] = (*pParent2)[i];
206  }
207 
208  return true;
209 }
210 
212 {
213  size_t i;
214  bool Continue = true;
215 
216  // generate a random order for the parents
218 
219  // reproduce in consecutive pairs
220  for (i = 0; i < mPopulationSize / 2; i++)
223  *mIndividual[mPopulationSize + i * 2],
224  *mIndividual[mPopulationSize + i * 2 + 1]);
225 
226  // check if there is one left over and just copy it
227  if (mPopulationSize % 2 > 0)
228  *mIndividual[2 * mPopulationSize - 1] = *mIndividual[mPopulationSize - 1];
229 
230  // mutate the offspring
231  for (i = mPopulationSize; i < 2 * mPopulationSize && Continue; i++)
232  {
233  mutate(*mIndividual[i]);
234  Continue &= evaluate(*mIndividual[i]);
236  }
237 
238  return Continue;
239 }
240 
241 // select mPopulationSize individuals
243 {
244  size_t i, j, nopp, opp;
245  size_t TotalPopulation = 2 * mPopulationSize;
246 
247  // tournament competition
248  mLosses = 0; // Set all wins to 0.
249 
250  // compete with ~ 20% of the TotalPopulation
251  nopp = std::max<size_t>(1, mPopulationSize / 5);
252 
253  // parents and offspring are all in competition
254  for (i = 0; i < TotalPopulation; i++)
255  for (j = 0; j < nopp; j++)
256  {
257  // get random opponent
258  do
259  {
260  opp = mpRandom->getRandomU((unsigned C_INT32)(TotalPopulation - 1));
261  }
262  while (i == opp);
263 
264  if (mValue[i] < mValue[opp])
265  mLosses[opp]++;
266  else
267  mLosses[i]++;
268  }
269 
270  // selection of top mPopulationSize winners
273  mLosses.array() + TotalPopulation,
274  mPivot);
275 
278 
279  return true;
280 }
281 
282 // check the best individual at this generation
284 {
285  size_t i, BestIndex = C_INVALID_INDEX;
287 
288  for (i = 0; i < mPopulationSize && !mLosses[i]; i++)
289  if (mValue[i] < BestValue)
290  {
291  BestIndex = i;
292  BestValue = mValue[i];
293  }
294 
295  return BestIndex;
296 }
297 
298 #ifdef DEBUG_OPT
299 // serialize the population to a file for debugging purposes
300 bool COptMethodGA::serializepop(size_t first, size_t last)
301 {
302  size_t Last = std::min(last, (size_t) mPopulationSize);
303 
304  size_t i;
305  size_t j;
306 
307  std::ofstream ofile;
308 //std::ifstream test("gapop.txt");
309 //if( !test.good() ) return false;
310 
311  // open the file for output, in append mode
312  ofile.open( "gapop.txt", std::ios::out | std::ios::app );
313  if ( ! ofile.is_open() )
314  {
315  std::cerr << "error opening file \'gapop.txt\'" << std::endl;
316  return false;
317  }
318 
319  for (i = first; i < Last; i++)
320  {
321  for (j = 0; j < mVariableSize; j++)
322  {
323  C_FLOAT64 & mut = (*mIndividual[i])[j];
324  // print this parameter
325  ofile << mut << "\t";
326  }
327  // print the fitness of the individual
328  ofile << mValue[i] << std::endl;
329  }
330  ofile << std::endl;
331  ofile.close();
332  return true;
333 }
334 #endif
335 
336 // initialise the population
337 bool COptMethodGA::creation(size_t first,
338  size_t last)
339 {
340  size_t Last = std::min(last, (size_t) mPopulationSize);
341 
342  size_t i;
343  size_t j;
344 
345  C_FLOAT64 mn;
346  C_FLOAT64 mx;
347  C_FLOAT64 la;
348 
349  bool Continue = true;
350 
351  for (i = first; i < Last && Continue; i++)
352  {
353  for (j = 0; j < mVariableSize; j++)
354  {
355  // calculate lower and upper bounds
356  COptItem & OptItem = *(*mpOptItem)[j];
357  mn = *OptItem.getLowerBoundValue();
358  mx = *OptItem.getUpperBoundValue();
359 
360  C_FLOAT64 & mut = (*mIndividual[i])[j];
361 
362  try
363  {
364  // determine if linear or log scale
365  if ((mn < 0.0) || (mx <= 0.0))
366  mut = mn + mpRandom->getRandomCC() * (mx - mn);
367  else
368  {
369  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
370 
371  if (la < 1.8)
372  mut = mn + mpRandom->getRandomCC() * (mx - mn);
373  else
374  mut = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC());
375  }
376  }
377 
378  catch (...)
379  {
380  mut = (mx + mn) * 0.5;
381  }
382 
383  // force it to be within the bounds
384  switch (OptItem.checkConstraint(mut))
385  {
386  case - 1:
387  mut = *OptItem.getLowerBoundValue();
388  break;
389 
390  case 1:
391  mut = *OptItem.getUpperBoundValue();
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 
400  // calculate its fitness
401  Continue &= evaluate(*mIndividual[i]);
403  }
404 
405  return Continue;
406 }
407 
409 {
411 }
412 
414 {
415  cleanup();
416 
417  size_t i;
418 
419  if (!COptMethod::initialize())
420  {
421  if (mpCallBack)
423 
424  return false;
425  }
426 
427  mGenerations = * getValue("Number of Generations").pUINT;
428  mGeneration = 0;
429 
430  if (mpCallBack)
431  mhGenerations =
432  mpCallBack->addItem("Current Generation",
433  mGeneration,
434  & mGenerations);
435 
436  mGeneration++;
437 
438  mPopulationSize = * getValue("Population Size").pUINT;
439  mpRandom =
440  CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
441  * getValue("Seed").pUINT);
442 
443  mVariableSize = mpOptItem->size();
444 
445  mIndividual.resize(2 * mPopulationSize);
446 
447  for (i = 0; i < 2 * mPopulationSize; i++)
449 
451  mCrossOverFalse = false;
453 
454  mValue.resize(2 * mPopulationSize);
455  mValue = std::numeric_limits<C_FLOAT64>::infinity();
456  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
457 
458  mpPermutation = new CPermutation(mpRandom, mPopulationSize);
459 
460  mLosses.resize(2 * mPopulationSize);
461 
462  // Initialize the variance for mutations
463  mMutationVarians = 0.1;
464 
465  return true;
466 }
467 
469 {
470  size_t i;
471 
472  pdelete(mpRandom);
474 
475  for (i = 0; i < mIndividual.size(); i++)
476  pdelete(mIndividual[i]);
477 
478  return true;
479 }
480 
482 {
483  bool Continue = true;
484 
485  if (!initialize())
486  {
487  // initialisation failed, we exit
488  if (mpCallBack)
490  return false;
491  }
492 
493  // Counters to determine whether the optimization process has stalled
494  // They count the number of generations without advances.
495  size_t Stalled, Stalled10, Stalled30, Stalled50;
496  Stalled = Stalled10 = Stalled30 = Stalled50 = 0;
497 
498  size_t i;
499 
500  // initialise the population
501  // first individual is the initial guess
502  for (i = 0; i < mVariableSize; i++)
503  {
504  C_FLOAT64 & mut = (*mIndividual[0])[i];
505  COptItem & OptItem = *(*mpOptItem)[i];
506 
507  mut = OptItem.getStartValue();
508 
509  // force it to be within the bounds
510  switch (OptItem.checkConstraint(mut))
511  {
512  case - 1:
513  mut = *OptItem.getLowerBoundValue();
514  break;
515 
516  case 1:
517  mut = *OptItem.getUpperBoundValue();
518  break;
519  }
520 
521  // We need to set the value here so that further checks take
522  // account of the value.
523  (*(*mpSetCalculateVariable)[i])(mut);
524  }
525 
526  Continue &= evaluate(*mIndividual[0]);
528 
529  if (!isnan(mEvaluationValue))
530  {
531  // and store that value
532  mBestValue = mValue[0];
533  Continue &= mpOptProblem->setSolution(mBestValue, *mIndividual[0]);
534 
535  // We found a new best value lets report it.
537  }
538 
539  // the others are random
540  Continue &= creation(1, mPopulationSize);
541 
542 #ifdef DEBUG_OPT
544 #endif
545 
546  Continue &= select();
547  mBestIndex = fittest();
548 
549  if (mBestIndex != C_INVALID_INDEX &&
551  {
552  // and store that value
555 
556  // We found a new best value lets report it.
558  }
559 
560  // test if the user wants to stop, and do so if needed
561  if (!Continue)
562  {
563  if (mpCallBack)
565  cleanup();
566  return true;
567  }
568 
569  // ITERATE FOR gener GENERATIONS
570  for (mGeneration = 2;
571  mGeneration <= mGenerations && Continue;
572  mGeneration++, Stalled++, Stalled10++, Stalled30++, Stalled50++)
573  {
574  // perturb the population if we have stalled for a while
575  if (Stalled > 50 && Stalled50 > 50)
576  {
577  Continue &= creation((size_t)(mPopulationSize / 2),
579  Stalled10 = Stalled30 = Stalled50 = 0;
580  }
581  else if (Stalled > 30 && Stalled30 > 30)
582  {
583  Continue &= creation((size_t)(mPopulationSize * 0.7),
585  Stalled10 = Stalled30 = 0;
586  }
587  else if (Stalled > 10 && Stalled10 > 10)
588  {
589  Continue &= creation((size_t)(mPopulationSize * 0.9),
591  Stalled10 = 0;
592  }
593  // replicate the individuals
594  else
595  Continue &= replicate();
596 
597  // select the most fit
598  Continue &= select();
599 
600  // get the index of the fittest
601  mBestIndex = fittest();
602 
603  if (mBestIndex != C_INVALID_INDEX &&
605  {
606  // reset the stalled counters, since we made progress this time
607  Stalled = Stalled10 = Stalled30 = Stalled50 = 0;
608  // keep best value
610  // pass the current best value upstream
612  // We found a new best value lets report it.
614  }
615 
616  if (mpCallBack)
617  Continue &= mpCallBack->progressItem(mhGenerations);
618 #ifdef DEBUG_OPT
620 #endif
621  }
622 
623  if (mpCallBack)
625 
626  cleanup();
627  return true;
628 }
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
virtual bool initialize()
COptTask * mpParentTask
Definition: COptMethod.h:58
#define pdelete(p)
Definition: copasi.h:215
virtual bool initialize()
Definition: COptMethod.cpp:189
C_FLOAT64 mMutationVarians
Definition: COptMethodGA.h:217
virtual C_FLOAT64 getRandomNormal(const C_FLOAT64 &mean, const C_FLOAT64 &sd)
Definition: CRandom.cpp:292
bool swap(size_t from, size_t to)
unsigned C_INT32 mGenerations
Definition: COptMethodGA.h:152
const size_t & next()
CVector< size_t > mPivot
Definition: COptMethodGA.h:212
std::vector< CVector< C_FLOAT64 > * > mIndividual
Definition: COptMethodGA.h:177
CVector< bool > mCrossOverFalse
Definition: COptMethodGA.h:182
bool crossover(const CVector< C_FLOAT64 > &parent1, const CVector< C_FLOAT64 > &parent2, CVector< C_FLOAT64 > &child1, CVector< C_FLOAT64 > &child2)
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
C_FLOAT64 mEvaluationValue
Definition: COptMethodGA.h:192
virtual void output(const COutputInterface::Activity &activity)
bool applyPartialPivot(const CVector< size_t > &pivot, const size_t &ordered, SwapMethod swap)
Definition: CSort.h:375
#define C_INT32
Definition: copasi.h:90
size_t fittest()
CRandom * mpRandom
Definition: COptMethodGA.h:167
virtual bool progressItem(const size_t &handle)
CVector< size_t > mLosses
Definition: COptMethodGA.h:207
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
virtual bool calculate()
CVector< bool > mCrossOver
Definition: COptMethodGA.h:187
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
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
void initObjects()
const Value & getValue() const
size_t mVariableSize
Definition: COptMethodGA.h:172
size_t mBestIndex
Definition: COptMethodGA.h:220
void shuffle(const size_t &swaps=C_INVALID_INDEX)
unsigned C_INT32 mGeneration
Definition: COptMethodGA.h:221
unsigned C_INT32 mPopulationSize
Definition: COptMethodGA.h:162
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
bool serializepop(size_t first, size_t last)
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()
virtual ~COptMethodGA()
bool mutate(CVector< C_FLOAT64 > &individual)
const C_FLOAT64 & getStartValue() const
Definition: COptItem.cpp:199
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
COptMethodGA(const COptMethodGA &src, const CCopasiContainer *pParent=NULL)
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 mBestValue
Definition: COptMethodGA.h:219
CType * array()
Definition: CVector.h:139
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
virtual bool cleanup()
bool addParameter(const CCopasiParameter &parameter)
bool evaluate(const CVector< C_FLOAT64 > &individual)
bool creation(size_t first, size_t last=std::numeric_limits< size_t >::max())
CVector< C_FLOAT64 > mValue
Definition: COptMethodGA.h:197
size_t mhGenerations
Definition: COptMethodGA.h:157
const C_FLOAT64 & getCalculateValue() const
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 optimise()
CPermutation * mpPermutation
Definition: COptMethodGA.h:202
#define max(a, b)
Definition: f2c.h:176