COPASI API  4.16.103
COptMethodEP2.cpp
Go to the documentation of this file.
1 // Begin CVS Header
2 // $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/optimization/COptMethodEP2.cpp,v $
3 // $Revision: 1.10 $
4 // $Name: $
5 // $Author: shoops $
6 // $Date: 2012/04/23 21:11:20 $
7 // End CVS Header
8 
9 // Copyright (C) 2012 by Pedro Mendes, Virginia Tech Intellectual
10 // Properties, Inc., University of Heidelberg, and The University
11 // of Manchester.
12 // All rights reserved.
13 
14 // Copyright (C) 2008 by Pedro Mendes, Virginia Tech Intellectual
15 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
16 // and The University of Manchester.
17 // All rights reserved.
18 
19 // Copyright (C) 2001 - 2007 by Pedro Mendes, Virginia Tech Intellectual
20 // Properties, Inc. and EML Research, gGmbH.
21 // All rights reserved.
22 
23 /***************************************************************************
24  COptMethodEP2.cpp - Evolutionary Program Optimizer
25  -------------------------------------------------
26 
27  Implemented by Dingjun Chen
28 
29  Starting Date: Dec. 2003
30 
31  COPASI project group
32  ***************************************************************************/
33 
34 /***************************************************************************
35  * This is the implementation of the Eovlutionary Program for Optimization. The
36  * class is inherited from the COptAlgorithm class
37  ***************************************************************************/
38 
39 #include <vector>
40 
41 #include "copasi.h"
42 #include "COptMethod.h"
43 #include "CRealProblem.h"
45 
47  COptMethod(src)
48 {}
49 
51 {
52  addParameter("EvolutionaryProgram2.Iterations",
54  (unsigned C_INT32) 10000);
55  addParameter("EvolutionaryProgram2.PopulationSize",
57  (unsigned C_INT32) 600);
58  addParameter("EvolutionaryProgram2.RandomGenerator.Type",
60  (C_INT32)CRandom::mt19937);
61  addParameter("EvolutionaryProgram2.RandomGenerator.Seed",
63  (C_INT32) 0);
64 }
65 
66 /**
67  * Destructor
68  */
70 {}
71 
72 /**
73  * GA Optimizer Function:
74  * Returns: nothing
75  */
77 {
78  NumGeneration = (C_INT32) getValue("EvolutionaryProgram2.Iterations");
79  PopulationSize = (C_INT32) getValue("EvolutionaryProgram2.PopulationSize");
81 
82  /* Create a random number generator */
84  Type = (CRandom::Type)(C_INT32) getValue("EvolutionaryProgram2.RandomGenerator.Type");
85  C_INT32 Seed;
86  Seed = (C_INT32) getValue("EvolutionaryProgram2.RandomGenerator.Seed");
87  CRandom * pRand = CRandom::createGenerator(Type, Seed);
88 
89  assert(pRand);
90 
91  const double ** Minimum = mpOptProblem->getParameterMin().array();
92  const double ** Maximum = mpOptProblem->getParameterMax().array();
93 
94  std::vector< UpdateMethod * > & Parameter = mpOptProblem->getCalculateVariableUpdateMethods();
95 
96  double current_best_value, la;
97  int i, j, last_update, u10, u30, u50;
98  bool linear;
99 
100  // create array for function value of candidate
101  CandidateValue = new double[2 * PopulationSize];
102 
103  // create array for function value rate of candidate
104  CandidateValueRate = new double[PopulationSize];
105 
106  // create array for tournament competition strategy
107  WinScore = new int[2 * PopulationSize];
108 
109  // create array for crossover points
110  CrossPoint = new int[NumParameter];
111 
112  // create the population array
114 
115  // create the individuals
116  for (i = 0; i < 2*PopulationSize; i++)
117  individual[i].resize(NumParameter);
118 
119  // Prepare inital population
120  //Generate the initial population at random
121  for (i = 0; i < PopulationSize; i++)
122  {
123  for (j = 0; j < NumParameter; j++)
124  {
125  try
126  {
127  // determine if linear or log scale
128  linear = FALSE; la = 1.0;
129 
130  if ((*Maximum[j] <= 0.0) || (*Minimum[j] < 0.0)) linear = TRUE;
131  else
132  {
133  la = log10(*Maximum[j]) - log10(std::min(*Minimum[j], std::numeric_limits< C_FLOAT64 >::epsilon()));
134 
135  if (la < 1.8) linear = TRUE;
136  }
137 
138  // set it to a random value within the interval
139  if (linear)
140  individual[i][j] = *Minimum[j] + pRand->getRandomCC() * (*Maximum[j] - *Minimum[j]);
141  else
142  individual[i][j] = *Minimum[j] * pow(10.0, la * pRand->getRandomCC());
143  }
144  catch (int)
145  {
146  individual[i][j] = (*Maximum[j] - *Minimum[j]) * 0.5 + *Minimum[j];
147  }
148  }
149 
150  try
151  {
152  // calculate its fitness value
153  for (int kk = 0; kk < NumParameter; kk++)
154  (*Parameter[kk])(individual[i][kk]);
155 
157  }
158  catch (int)
159  {
161  }
162  } // initialization ends
163 
164  //Set fitness map rate
165  // double q=0.8;
166  // double constant1=100;
167  double constant1 = 100;
168  // double constant2=20;
169  double constant2 = 10;
170  /* for (i=0; i<PopulationSize; i++)
171  {
172  CandidateValueRate[i]=1-q*pow((1-q),i);
173  }
174  */
175  // get the index of the fittest
177 
178  // initialise the update registers
179  last_update = 0;
180  u10 = u30 = u50 = 0;
181 
182  // store that value
183  current_best_value = Get_BestFoundSoFar_candidate();
184 
185  std::ofstream finalout("debugopt.dat");
186 
187  if (!finalout)
188  {
189  std::cout << "debugopt.dat cannot be opened!" << std::endl;
190  exit(1);
191  }
192 
193  finalout << "----------------------------- the best result at each generation---------------------" << std::endl;
194  finalout << "Generation\t" << "Best candidate value for object function\t" << "Display " << NumParameter << " parameters" << std::endl;
195  finalout << std::endl;
196  finalout.close();
197 
198  srand(time(NULL)); rand();
199 
200  // Iterations of GA
201  for (i = 0; i < NumGeneration; i++)
202  {
203  std::cout << std::endl;
204  std::cout << "EP2 is processing at generation " << i << std::endl;
205 
206  TrackDataFile(i);
207 
208  select(5);
209 
210  //Mutating all parents
211  for (int nn = PopulationSize; nn < 2*PopulationSize; nn++)
212  {
213  double mut;
214 
215  // mutate the parameters
216  for (int j = 0; j < NumParameter; j++)
217  {
218  if (floor(10*rand() / RAND_MAX) > 5)
219  {
220  // if(i<NumGeneration*0.5) mut = individual[nn-PopulationSize][j]*(1 + pRand->getRandomCC());
221  if (i < NumGeneration*0.5) mut = individual[nn - PopulationSize][j] * (1 + sqrt(constant1 * CandidateValueRate[nn] + constant2) * pRand->getRandomCC());
222  else
223  {
224  if (i < NumGeneration*0.6) mut = individual[nn - PopulationSize][j] * (1 + sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.5 * pRand->getRandomCC());
225  else
226  {
227  if (i < NumGeneration*0.7) mut = individual[nn - PopulationSize][j] * (1 + sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.25 * pRand->getRandomCC());
228  else
229  {
230  if (i < NumGeneration*0.8) mut = individual[nn - PopulationSize][j] * (1 + sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.1 * pRand->getRandomCC());
231  else
232  {
233  if (i < NumGeneration*0.9) mut = individual[nn - PopulationSize][j] * (1 + sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.01 * pRand->getRandomCC());
234  else
235  {
236  if (i < NumGeneration*0.95) mut = individual[nn - PopulationSize][j] * (1 + sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.001 * pRand->getRandomCC());
237  else mut = individual[nn - PopulationSize][j] * (1 + sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.0001 * pRand->getRandomCC());
238  }
239  }
240  }
241  }
242  }
243  }
244  else
245  {
246  if (i < NumGeneration*0.5) mut = individual[nn - PopulationSize][j] * (1 - sqrt(constant1 * CandidateValueRate[nn] + constant2) * pRand->getRandomCC());
247  else
248  {
249  if (i < NumGeneration*0.6) mut = individual[nn - PopulationSize][j] * (1 - sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.5 * pRand->getRandomCC());
250  else
251  {
252  if (i < NumGeneration*0.7) mut = individual[nn - PopulationSize][j] * (1 - sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.25 * pRand->getRandomCC());
253  else
254  {
255  if (i < NumGeneration*0.8) mut = individual[nn - PopulationSize][j] * (1 - sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.1 * pRand->getRandomCC());
256  else
257  {
258  if (i < NumGeneration*0.9) mut = individual[nn - PopulationSize][j] * (1 - sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.01 * pRand->getRandomCC());
259  else
260  {
261  if (i < NumGeneration*0.95) mut = individual[nn - PopulationSize][j] * (1 - sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.001 * pRand->getRandomCC());
262  else mut = individual[nn - PopulationSize][j] * (1 - sqrt(constant1 * CandidateValueRate[nn] + constant2) * 0.0001 * pRand->getRandomCC());
263  }
264  }
265  }
266  }
267  }
268  }
269 
270  // check boundary and force it to be within the bounds
271  if (mut <= *Minimum[j]) mut = *Minimum[j] + std::numeric_limits< C_FLOAT64 >::epsilon();
272  else
273  {
274  if (mut < *Minimum[j]) mut = *Minimum[j];
275  }
276 
277  if (mut >= *Maximum[j]) mut = *Maximum[j] - std::numeric_limits< C_FLOAT64 >::epsilon();
278  else
279  {
280  if (mut > *Maximum[j]) mut = *Maximum[j];
281  }
282 
283  // store it
284  individual[nn][j] = mut;
285  }
286 
287  // evaluate the fitness
288  for (int kk = 0; kk < NumParameter; kk++)
289  (*Parameter[kk])(individual[nn][kk]);
290 
292  }
293 
294  // select approprate individuals as new population
295  select(2);
296 
297  // get the index of the fittest
299 
300  if (Get_BestFoundSoFar_candidate() != current_best_value)
301  {
302  last_update = i;
303  current_best_value = Get_BestFoundSoFar_candidate();
304  }
305 
306  if (u10) u10--;
307 
308  if (u30) u30--;
309 
310  if (u50) u50--;
311 
312  if ((u50 == 0) && (i - last_update > 50))
313  {
314  for (int mm = PopulationSize / 2; mm < PopulationSize; mm++)
315  {
316  for (int jj = 0; jj < NumParameter; jj++)
317  {
318  try
319  {
320  // determine if linear or log scale
321  linear = FALSE; la = 1.0;
322 
323  if ((*Maximum[jj] <= 0.0) || (*Minimum[jj] < 0.0)) linear = TRUE;
324  else
325  {
326  la = log10(*Maximum[jj]) - log10(std::min(*Minimum[jj], std::numeric_limits< C_FLOAT64 >::epsilon()));
327 
328  if (la < 1.8) linear = TRUE;
329  }
330 
331  // set it to a random value within the interval
332  if (linear)
333  individual[mm][jj] = *Minimum[jj] + pRand->getRandomCC() * (*Maximum[jj] - *Minimum[jj]);
334  else
335  individual[mm][jj] = *Minimum[jj] * pow(10.0, la * pRand->getRandomCC());
336  }
337  catch (int)
338  {
339  individual[mm][jj] = (*Maximum[jj] - *Minimum[jj]) * 0.5 + *Minimum[jj];
340  }
341  }
342 
343  try
344  {
345  // calculate its fitness
346  for (int kk = 0; kk < NumParameter; kk++)
347  (*Parameter[kk])(individual[mm][kk]);
348 
350  }
351  catch (int)
352  {
354  }
355  }
356 
357  //end external for loop
359  u50 = 50; u30 = 30; u10 = 10;
360  }
361  else
362  {
363  if ((u30 == 0) && (i - last_update > 30))
364  {
365  for (int mm = (int)floor(PopulationSize * 0.7); mm < PopulationSize; mm++)
366  {
367  for (int jj = 0; jj < NumParameter; jj++)
368  {
369  try
370  {
371  // determine if linear or log scale
372  linear = FALSE; la = 1.0;
373 
374  if ((*Maximum[jj] <= 0.0) || (*Minimum[jj] < 0.0)) linear = TRUE;
375  else
376  {
377  la = log10(*Maximum[jj]) - log10(std::min(*Minimum[jj], std::numeric_limits< C_FLOAT64 >::epsilon()));
378 
379  if (la < 1.8) linear = TRUE;
380  }
381 
382  // set it to a random value within the interval
383  if (linear)
384  individual[mm][jj] = *Minimum[jj] + pRand->getRandomCC() * (*Maximum[jj] - *Minimum[jj]);
385  else
386  individual[mm][jj] = *Minimum[jj] * pow(10.0, la * pRand->getRandomCC());
387  }
388  catch (int)
389  {
390  individual[mm][jj] = (*Maximum[jj] - *Minimum[jj]) * 0.5 + *Minimum[jj];
391  }
392  }
393 
394  try
395  {
396  // calculate its fitness
397  for (int kk = 0; kk < NumParameter; kk++)
398  (*Parameter[kk])(individual[mm][kk]);
399 
401  }
402  catch (int)
403  {
405  }
406  }
407 
408  //end external for loop
410  u30 = 30; u10 = 10;
411  }
412  else
413  {
414  if ((u10 == 0) && (i - last_update > 10))
415  {
416  for (int mm = (int) floor(PopulationSize * 0.9); mm < PopulationSize; mm++)
417  {
418  for (int jj = 0; jj < NumParameter; jj++)
419  {
420  try
421  {
422  // determine if linear or log scale
423  linear = FALSE; la = 1.0;
424 
425  if ((*Maximum[jj] <= 0.0) || (*Minimum[jj] < 0.0)) linear = TRUE;
426  else
427  {
428  la = log10(*Maximum[jj]) - log10(std::min(*Minimum[jj], std::numeric_limits< C_FLOAT64 >::epsilon()));
429 
430  if (la < 1.8) linear = TRUE;
431  }
432 
433  // set it to a random value within the interval
434  if (linear)
435  individual[mm][jj] = *Minimum[jj] + pRand->getRandomCC() * (*Maximum[jj] - *Minimum[jj]);
436  else
437  individual[mm][jj] = *Minimum[jj] * pow(10.0, la * pRand->getRandomCC());
438  }
439  catch (int)
440  {
441  individual[mm][jj] = (*Maximum[jj] - *Minimum[jj]) * 0.5 + *Minimum[jj];
442  }
443  }
444 
445  try
446  {
447  // calculate its fitness
448  for (int kk = 0; kk < NumParameter; kk++)
449  (*Parameter[kk])(individual[mm][kk]);
450 
452  }
453  catch (int)
454  {
456  }
457  }
458 
459  //end external for loop
461  u10 = 10;
462  //u10=50;
463  }
464  }
465  }
466  } // end iteration of generations
467 
468  //store the combination of the BestFoundSoFar parameter values found so far
469  mpOptProblem->setSolutionVariables(individual[BestFoundSoFar]);
470 
471  //set the BestFoundSoFar function value
472  mpOptProblem->setSolutionValue(CandidateValue[BestFoundSoFar]);
473 
474  //free memory space
475  delete CrossPoint;
476  delete WinScore;
477  delete CandidateValue;
478 
479  std::cout << std::endl;
480  std::cout << "GA has successfully done!" << std::endl;
481 
482  return 0;
483 }
484 
485 // copy individual o to position d
486 void COptMethodEP2::copy(int o, int d)
487 {
488  int i;
489 
490  for (i = 0; i < NumParameter; i++)
491  {
492  individual[d][i] = individual[o][i];
493  }
494 
496 }
497 
498 // swap individuals o and d
499 void COptMethodEP2::swap(int o, int d)
500 {
501  int i;
502  double tmp;
503 
504  for (i = 0; i < NumParameter; i++)
505  {
506  tmp = individual[d][i];
507  individual[d][i] = individual[o][i];
508  individual[o][i] = tmp;
509  }
510 
511  tmp = CandidateValue[d];
513  CandidateValue[o] = tmp;
514 
515  i = WinScore[d];
516  WinScore[d] = WinScore[o];
517  WinScore[o] = i;
518 }
519 
520 // exchange individuals o and d
521 void COptMethodEP2::exchange(int o, int d)
522 {
523  int i;
524  double tmp;
525 
526  for (i = 0; i < NumParameter; i++)
527  {
528  tmp = individual[d][i];
529  individual[d][i] = individual[o][i];
530  individual[o][i] = tmp;
531  }
532 
533  tmp = CandidateValue[d];
535  CandidateValue[o] = tmp;
536 }
537 
538 // select PopulationSize individuals
539 void COptMethodEP2::select(int SelectionStrategy)
540 {
541  int i, j, TournamentSize, RandomRival;
542  int RandomIndividual;
543  //srand(time(NULL)); rand();
544 
545  switch (SelectionStrategy)
546  {
547  case 1: // parent-offspring competition
548 
549  for (i = PopulationSize; i < 2*PopulationSize; i++)
550  {
551  // if offspring is fitter keep it
552  for (j = 0; j < PopulationSize; j++)
553  {
554  if (CandidateValue[i] < CandidateValue[j]) exchange(i, j);
555  }
556  }
557 
558  break;
559  case 2: // tournament competition
560  // compete with 20% of the population
561  TournamentSize = PopulationSize / 5;
562 
563  // but at least one
564  if (TournamentSize < 1) TournamentSize = 1;
565 
566  // parents and offspring are all in competition
567  for (i = 0; i < 2*PopulationSize; i++)
568  {
569  WinScore[i] = 0;
570 
571  for (j = 0; j < TournamentSize; j++)
572  {
573  // get random rival
574  RandomRival = (int)fabs(floor((PopulationSize * 2 - 1) * rand() / RAND_MAX));
575 
576  if (CandidateValue[i] <= CandidateValue[RandomRival]) WinScore[i]++;
577  }
578  }
579 
580  // selection of top PopulationSize winners
581  for (i = 0; i < PopulationSize; i++)
582  {
583  for (j = i + 1; j < 2*PopulationSize; j++)
584  {if (WinScore[i] < WinScore[j]) swap(i, j);}
585  }
586 
587  break;
588 
589  // ranking strategy without proportionate-fitness
590  case 3:
591 
592  for (i = 0; i < PopulationSize; i++)
593  {
594  for (j = i + 1; j < 2*PopulationSize; j++)
595  {
596  if (CandidateValue[i] > CandidateValue[j]) exchange(i, j);
597  }
598  }
599 
600  break;
601  // Randomly select P individuals from population of 2P
602  case 4:
603 
604  for (i = 0; i < PopulationSize; i++)
605  {
606  RandomIndividual = (int)fabs(floor((PopulationSize * 2 - 1) * rand() / RAND_MAX));
607  exchange(i, RandomIndividual);
608  }
609 
610  break;
611 
612  // ranking strategy for EP2
613  case 5:
614 
615  for (i = 0; i < PopulationSize; i++)
616  {
617  for (j = i + 1; j < PopulationSize; j++)
618  {
619  if (CandidateValue[i] > CandidateValue[j]) exchange(i, j);
620  }
621  }
622 
623  for (i = 0; i < PopulationSize; i++)
624  {
625  CandidateValueRate[i] = (CandidateValue[i] - CandidateValue[0]) / (CandidateValue[PopulationSize - 1] - CandidateValue[0]);
626  }
627 
628  break;
629  }
630 }
631 
632 // check the best individual at this generation
634 {
635  int i, b;
636  double f;
637  f = CandidateValue[0];
638  b = 0;
639 
640  for (i = 1; i < PopulationSize; i++)
641  {
642  if (CandidateValue[i] < f)
643  {
644  b = i;
645  f = CandidateValue[i];
646  }
647  }
648 
649  return b;
650 }
651 
653 {
654  std::ofstream finalout("debugopt.dat", std::ios::app);
655 
656  if (!finalout)
657  {
658  std::cout << "debugopt.dat cannot be opened!" << std::endl;
659  exit(1);
660  }
661 
662  finalout << "#" << i << "\t" << std::setprecision(8) << CandidateValue[BestFoundSoFar] << std::endl;
663 
664  for (int j = 0; j < NumParameter; j++)
665  {
666  finalout << std::setprecision(8) << individual[BestFoundSoFar][j] << "\t";
667  }
668 
669  finalout << std::endl;
670  finalout << std::endl;
671 
672  finalout.close();
673 }
674 
675 /**********************************the end **************************************/
size_t getVariableSize() const
virtual void exchange(int o, int d)
virtual int fittest(void)
#define TRUE
Definition: CGA.h:25
qreal linear(qreal a, qreal b, qreal t)
double * CandidateValue
Definition: COptMethodEP2.h:61
virtual void swap(int o, int d)
const std::vector< UpdateMethod * > & getCalculateVariableUpdateMethods() const
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
COptProblem * mpOptProblem
Definition: COptMethod.h:56
#define C_INT32
Definition: copasi.h:90
virtual ~COptMethodEP2()
CVector< CVector< C_FLOAT64 > > individual
Definition: COptMethodEP2.h:60
#define FALSE
Definition: CGA.h:26
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
const Value & getValue() const
virtual void select(int method)
virtual bool optimise()
virtual void TrackDataFile(int i)
virtual void copy(int o, int d)
double Get_BestFoundSoFar_candidate()
bool addParameter(const CCopasiParameter &parameter)
double * CandidateValueRate
Definition: COptMethodEP2.h:62
#define min(a, b)
Definition: f2c.h:175
#define max(a, b)
Definition: f2c.h:176