COPASI API  4.16.103
COptMethodGASR.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 "COptMethodGASR.h"
20 #include "COptProblem.h"
21 #include "COptItem.h"
22 #include "COptTask.h"
23 
28 
30  COptMethod(CCopasiTask::optimization, CCopasiMethod::GeneticAlgorithmSR, pParent),
31  mGenerations(0),
32  mPopulationSize(0),
33  mpRandom(NULL),
34  mVariableSize(0),
35  mIndividual(0),
36  mCrossOverFalse(0),
37  mCrossOver(0),
38  mValue(0),
39  mpPermutation(NULL),
40  mWins(0),
41  mMutationVarians(0.1),
42  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
43  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
44  mBestIndex(C_INVALID_INDEX),
45  mGeneration(0)
46 
47 {
48  addParameter("Number of Generations", CCopasiParameter::UINT, (unsigned C_INT32) 200);
49  addParameter("Population Size", CCopasiParameter::UINT, (unsigned C_INT32) 20);
50  addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
51  addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
52  addParameter("Pf", CCopasiParameter::DOUBLE, (C_FLOAT64) 0.475); //*****ADDED for SR
53 
54  initObjects();
55 }
56 
58  const CCopasiContainer * pParent):
59  COptMethod(src, pParent),
60  mGenerations(0),
61  mPopulationSize(0),
62  mpRandom(NULL),
63  mVariableSize(0),
64  mIndividual(0),
65  mCrossOverFalse(0),
66  mCrossOver(0),
67  mValue(0),
68  mpPermutation(NULL),
69  mWins(0),
70  mMutationVarians(0.1),
71  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
72  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
73  mBestIndex(C_INVALID_INDEX),
74  mGeneration(0)
75 {initObjects();}
76 
78 {cleanup();}
79 
80 // evaluate the fitness of one individual
81 bool COptMethodGASR::evaluate(const CVector< C_FLOAT64 > & /* individual */)
82 {
83  bool Continue = true;
84 
85  // We do not need to check whether the parametric constraints are fulfilled
86  // since this method allows for parameters outside the bounds
87 
88  // evaluate the fitness
89  Continue = mpOptProblem->calculate();
90 
91  // We do not need to check whether the functional constraints are fulfilled
92  // since this method allows for solutions outside the bounds.
93 
95 
96  return Continue;
97 }
98 
99 bool COptMethodGASR::swap(size_t from, size_t to)
100 {
101  CVector< C_FLOAT64 > * pTmp = mIndividual[to];
102  mIndividual[to] = mIndividual[from];
103  mIndividual[from] = pTmp;
104 
105  //**** Added for swapping individual Phi values also for Stochastic Ranking
106  C_FLOAT64 dTmp = mPhi[to];
107  mPhi[to] = mPhi[from];
108  mPhi[from] = dTmp;
109 
110  dTmp = mValue[to];
111  mValue[to] = mValue[from];
112  mValue[from] = dTmp;
113 
114  size_t iTmp = mWins[to];
115  mWins[to] = mWins[from];
116  mWins[from] = iTmp;
117 
118  return true;
119 }
120 
121 //mutate one individual
123 {
124  size_t j;
125 
126  // mutate the parameters
127  for (j = 0; j < mVariableSize; j++)
128  {
129  C_FLOAT64 & mut = individual[j];
130 
131  // calculate the mutatated parameter
133 
134  // for SR do not force to be within bounds
135 
136  // We need to set the value here so that further checks take
137  // account of the value.
138  (*(*mpSetCalculateVariable)[j])(mut);
139  }
140 
141  return true;
142 }
143 
145  const CVector< C_FLOAT64 > & parent2,
146  CVector< C_FLOAT64 > & child1,
147  CVector< C_FLOAT64 > & child2)
148 {
149  size_t i, crp;
150  size_t nCross = 0;
151 
153 
154  if (mVariableSize > 1)
155  nCross = mpRandom->getRandomU((unsigned C_INT32)(mVariableSize / 2));
156 
157  if (nCross == 0)
158  {
159  // if less than 0 just copy parent to child
160  child1 = parent1;
161  child2 = parent2;
162 
163  return true;
164  }
165 
166  // choose cross over points;
167  // We do not mind if a crossover point gets drawn twice
168  for (i = 0; i < nCross; i++)
169  {
170  crp = mpRandom->getRandomU((unsigned C_INT32)(mVariableSize - 1));
171  mCrossOver[crp] = true;
172  }
173 
174  const CVector< C_FLOAT64 > * pParent1 = & parent1;
175 
176  const CVector< C_FLOAT64 > * pParent2 = & parent2;
177 
178  const CVector< C_FLOAT64 > * pTmp;
179 
180  for (i = 0; i < mVariableSize; i++)
181  {
182  if (mCrossOver[i])
183  {
184  pTmp = pParent1;
185  pParent1 = pParent2;
186  pParent2 = pTmp;
187  }
188 
189  child1[i] = (*pParent1)[i];
190  child2[i] = (*pParent2)[i];
191  }
192 
193  return true;
194 }
195 
197 {
198  size_t i;
199  bool Continue = true;
200 
201  // generate a random order for the parents
203 
204  // reproduce in consecutive pairs
205  for (i = 0; i < mPopulationSize / 2; i++)
208  *mIndividual[mPopulationSize + i * 2],
209  *mIndividual[mPopulationSize + i * 2 + 1]);
210 
211  // check if there is one left over and just copy it
212  if (mPopulationSize % 2 > 0)
213  *mIndividual[2 * mPopulationSize - 1] = *mIndividual[mPopulationSize - 1];
214 
215  // mutate the offspring
216  for (i = mPopulationSize; i < 2 * mPopulationSize && Continue; i++)
217  {
218  mutate(*mIndividual[i]);
219  Continue = evaluate(*mIndividual[i]);
221 
222  /* Calculate the phi value of the individual for SR*/
223  mPhi[i] = phi(i);
224  }
225 
226  return Continue;
227 }
228 
229 // select mPopulationSize individuals
231 {
232  size_t i, j;
233  size_t TotalPopulation = 2 * mPopulationSize;
234  bool wasSwapped;
235  size_t sweepNum = TotalPopulation; // This is default based on paper
236 
237  // Selection Method for Stochastic Ranking
238  // stochastic ranking "bubble sort"
239 
240  for (i = 0; i < sweepNum; i++) // Here sweepNum is optimal number of sweeps from paper
241  {
242  wasSwapped = false;
243 
244  // :TODO: since we are only interested in mPopulationSize highest ranked
245  // individuals the upper limit of the loop can be improved.
246  for (j = 0; j < TotalPopulation - 1; j++) // lambda is number of individuals
247  {
248  if ((mPhi[j] == 0 && mPhi[j + 1] == 0) || // within bounds
249  (mpRandom->getRandomOO() < mPf)) // random chance to compare values outside bounds
250  {
251  // compare obj fcn using mValue alternative code
252  if (mValue[j] > mValue[j + 1])
253  {
254  swap(j, j + 1);
255  wasSwapped = true;
256  }
257  }
258  else //mPhi values are not equal and not within boundary
259  {
260  if (mPhi[j] > mPhi[j + 1]) // j further outside then j+1
261  {
262  swap(j, j + 1);
263  wasSwapped = true;
264  }
265  }
266  }
267 
268  // if no swap then break
269  if (wasSwapped == false) break;
270  }
271 
272  return true;
273 }
274 
275 // evaluate the distance of parameters and constraints to boundaries
277 {
278  C_FLOAT64 phiVal = 0.0;
279  C_FLOAT64 phiCalc;
280 
281  std::vector< COptItem * >::const_iterator it = mpOptItem->begin();
282  std::vector< COptItem * >::const_iterator end = mpOptItem->end();
283  C_FLOAT64 * pValue = mIndividual[indivNum]->array();
284 
285  for (; it != end; ++it, pValue++)
286  {
287  switch ((*it)->checkConstraint())
288  {
289  case - 1:
290  phiCalc = *(*it)->getLowerBoundValue() - *pValue;
291  phiVal += phiCalc * phiCalc;
292  break;
293 
294  case 1:
295  phiCalc = *pValue - *(*it)->getUpperBoundValue();
296  phiVal += phiCalc * phiCalc;
297  break;
298  }
299  }
300 
301  it = mpOptContraints->begin();
302  end = mpOptContraints->end();
303 
304  for (; it != end; ++it)
305  {
306  phiCalc = (*it)->getConstraintViolation();
307 
308  if (phiCalc > 0.0)
309  phiVal += phiCalc * phiCalc;
310  }
311 
312  return phiVal;
313 }
314 
315 // check the best individual at this generation
317 {
318  size_t i, BestIndex = C_INVALID_INDEX;
320 
321  for (i = 0; i < mPopulationSize; i++)
322  if (mValue[i] < BestValue && !(mPhi[i] != 0))
323  {
324  BestIndex = i;
325  BestValue = mValue[i];
326  }
327 
328  return BestIndex;
329 }
330 
331 // initialise the population
332 bool COptMethodGASR::creation(size_t first,
333  size_t last)
334 {
335  size_t Last = std::min< size_t >(last, mPopulationSize);
336 
337  size_t i;
338  size_t j;
339 
340  C_FLOAT64 mn;
341  C_FLOAT64 mx;
342  C_FLOAT64 la;
343 
344  bool Continue = true;
345 
346  for (i = first; i < Last && Continue; i++)
347  {
348  for (j = 0; j < mVariableSize; j++)
349  {
350  // calculate lower and upper bounds
351  COptItem & OptItem = *(*mpOptItem)[j];
352  mn = *OptItem.getLowerBoundValue();
353  mx = *OptItem.getUpperBoundValue();
354 
355  C_FLOAT64 & mut = (*mIndividual[i])[j];
356 
357  try
358  {
359  // determine if linear or log scale
360  if ((mn < 0.0) || (mx <= 0.0))
361  mut = mn + mpRandom->getRandomCC() * (mx - mn);
362  else
363  {
364  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
365 
366  if (la < 1.8)
367  mut = mn + mpRandom->getRandomCC() * (mx - mn);
368  else
369  mut = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC());
370  }
371  }
372 
373  catch (...)
374  {
375  mut = (mx + mn) * 0.5;
376  }
377 
378  // We need to set the value here so that further checks take
379  // account of the value.
380  (*(*mpSetCalculateVariable)[j])(mut);
381  }
382 
383  // calculate its fitness
384  Continue = evaluate(*mIndividual[i]);
386 
387  /* Calculate the phi value of the individual for SR*/
388  mPhi[i] = phi(i);
389  }
390 
391  return Continue;
392 }
393 
395 {
397 }
398 
400 {
401  cleanup();
402 
403  size_t i;
404 
405  if (!COptMethod::initialize()) return false;
406 
407  mGeneration = 0;
408  mGenerations = * getValue("Number of Generations").pUINT;
409 
410  if (mpCallBack)
411  mhGenerations =
412  mpCallBack->addItem("Current Generation",
413  mGeneration,
414  & mGenerations);
415 
416  mGeneration++;
417 
418  mPopulationSize = * getValue("Population Size").pUINT;
419  mPf = *(C_FLOAT64*) getValue("Pf").pDOUBLE;
420 
421  if (mPf < 0.0 || 1.0 < mPf)
422  {
423  mPf = 0.475;
424  setValue("Pf", mPf);
425  }
426 
427  mpRandom =
428  CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
429  * getValue("Seed").pUINT);
430 
431  mVariableSize = mpOptItem->size();
432 
433  mIndividual.resize(2 * mPopulationSize);
435 
436  for (i = 0; i < 2 * mPopulationSize; i++)
438 
440  mCrossOverFalse = false;
442 
443  mValue.resize(2 * mPopulationSize);
444 
445  mpPermutation = new CPermutation(mpRandom, mPopulationSize);
446 
447  mWins.resize(2 * mPopulationSize);
448 
449  // initialise the variance for mutations
450  mMutationVarians = 0.1;
451 
452  return true;
453 }
454 
456 {
457  size_t i;
458 
459  pdelete(mpRandom);
461 
462  for (i = 0; i < mIndividual.size(); i++)
463  pdelete(mIndividual[i]);
464 
465  return true;
466 }
467 
469 {
470  bool Continue = true;
471 
472  if (!initialize())
473  {
474  if (mpCallBack)
476 
477  return false;
478  }
479 
480  // Counters to determine whether the optimization process has stalled
481  // They count the number of generations without advances.
482  size_t Stalled, Stalled10, Stalled30, Stalled50;
483  Stalled = Stalled10 = Stalled30 = Stalled50 = 0;
484 
485  size_t i;
486 
487  // initialise the population
488  // first individual is the initial guess
489  for (i = 0; i < mVariableSize; i++)
490  (*mIndividual[0])[i] = (*mpOptItem)[i]->getStartValue();
491 
492  // calculate the fitness
493  size_t j;
494  std::vector< UpdateMethod *>::const_iterator itMethod = mpSetCalculateVariable->begin();
495 
496  // set the paramter values
497  for (j = 0; j < mVariableSize; j++, ++itMethod)
498  (**itMethod)((*mIndividual[0])[j]);
499 
500  Continue = evaluate(*mIndividual[0]);
502 
503  /* Calculate the phi value of the individual for SR*/
504  mPhi[0] = phi(0);
505 
506  // the others are random
507  Continue = creation(1, mPopulationSize);
508 
509  // get the index of the fittest
510  mBestIndex = fittest();
511 
513  {
514  // and store that value
517 
518  // We found a new best value lets report it.
520  }
521 
522  if (!Continue)
523  {
524  if (mpCallBack)
526 
527  cleanup();
528  return true;
529  }
530 
531  // ITERATE FOR gener GENERATIONS
532  for (mGeneration = 2;
533  mGeneration <= mGenerations && Continue;
534  mGeneration++, Stalled++, Stalled10++, Stalled30++, Stalled50++)
535  {
536  // perturb the population if we have stalled for a while
537  if (Stalled > 50 && Stalled50 > 50)
538  {
539  Continue = creation((size_t)(mPopulationSize * 0.5),
541  Stalled10 = Stalled30 = Stalled50 = 0;
542  }
543  else if (Stalled > 30 && Stalled30 > 30)
544  {
545  Continue = creation((size_t)(mPopulationSize * 0.7),
547  Stalled10 = Stalled30 = 0;
548  }
549  else if (Stalled > 10 && Stalled10 > 10)
550  {
551  Continue = creation((size_t)(mPopulationSize * 0.9),
553  Stalled10 = 0;
554  }
555  // replicate the individuals
556  else
557  Continue = replicate();
558 
559  // select the most fit
560  Continue = select();
561 
562  // get the index of the fittest
563  mBestIndex = fittest();
564 
565  if (mBestIndex != C_INVALID_INDEX &&
567  {
568  Stalled = Stalled10 = Stalled30 = Stalled50 = 0;
570 
572 
573  // We found a new best value lets report it.
575  }
576 
577  if (mpCallBack)
579  }
580 
581  if (mpCallBack)
583 
584  cleanup();
585 
586  return true;
587 }
bool mutate(CVector< C_FLOAT64 > &individual)
COptTask * mpParentTask
Definition: COptMethod.h:58
virtual bool cleanup()
#define pdelete(p)
Definition: copasi.h:215
virtual bool initialize()
Definition: COptMethod.cpp:189
C_FLOAT64 mBestValue
C_FLOAT64 phi(size_t indvNum)
virtual C_FLOAT64 getRandomNormal(const C_FLOAT64 &mean, const C_FLOAT64 &sd)
Definition: CRandom.cpp:292
CVector< C_FLOAT64 > mPhi
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
const size_t & next()
virtual bool optimise()
CRandom * mpRandom
unsigned C_INT32 mPopulationSize
virtual ~COptMethodGASR()
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
const std::vector< COptItem * > * mpOptContraints
Definition: COptMethod.h:75
virtual void output(const COutputInterface::Activity &activity)
#define C_INT32
Definition: copasi.h:90
virtual bool progressItem(const size_t &handle)
CPermutation * mpPermutation
std::vector< CVector< C_FLOAT64 > * > mIndividual
unsigned C_INT32 mGenerations
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)
C_FLOAT64 mEvaluationValue
const std::vector< UpdateMethod * > * mpSetCalculateVariable
Definition: COptMethod.h:65
const Value & getValue() const
void shuffle(const size_t &swaps=C_INVALID_INDEX)
CVector< C_FLOAT64 > mValue
bool setValue(const std::string &name, const CType &value)
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
COptMethodGASR(const COptMethodGASR &src, const CCopasiContainer *pParent=NULL)
unsigned C_INT32 * pUINT
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
virtual bool finishItem(const size_t &handle)
C_FLOAT64 mMutationVarians
CVector< bool > mCrossOver
virtual bool initialize()
CVector< size_t > mWins
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
#define C_FLOAT64
Definition: copasi.h:92
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
bool addParameter(const CCopasiParameter &parameter)
bool evaluate(const CVector< C_FLOAT64 > &individual)
bool crossover(const CVector< C_FLOAT64 > &parent1, const CVector< C_FLOAT64 > &parent2, CVector< C_FLOAT64 > &child1, CVector< C_FLOAT64 > &child2)
CVector< bool > mCrossOverFalse
bool swap(size_t from, size_t to)
const C_FLOAT64 & getCalculateValue() const
unsigned C_INT32 mGeneration
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
bool creation(size_t first, size_t last=std::numeric_limits< size_t >::max())
#define min(a, b)
Definition: f2c.h:175
#define max(a, b)
Definition: f2c.h:176