COPASI API  4.16.103
COptMethodSS.cpp
Go to the documentation of this file.
1 // Copyright (C) 2013 - 2014 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // DEBUG_OPT allows doing debug output only on this class and is preferable
7 // to using COPASI_DEBUG, which would make this output on for any debug version
8 #undef DEBUG_OPT
9 
10 #include <limits>
11 #include <string>
12 #include <cmath>
13 #ifdef DEBUG_OPT
14 # include <iostream>
15 # include <fstream>
16 #endif
17 
18 #include "copasi.h"
19 #include "COptMethodSS.h"
20 #include "COptProblem.h"
22 #include "COptItem.h"
23 #include "COptTask.h"
24 
27 #include "utilities/CSort.h"
29 
31  COptMethod(CCopasiTask::optimization, CCopasiMethod::ScatterSearch, pParent),
32  mIterations(0),
33  mPopulationSize(0),
34  mVariableSize(0),
35  mRefSet(0),
36  mRefSetVal(0),
37  mPool(0),
38  mPoolVal(0),
39  mPoolSize(0),
40  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
41  mIteration(0),
42  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
43  mBestIndex(C_INVALID_INDEX),
44  mpRandom(NULL),
45  mpOptProblemLocal(NULL),
46  mpLocalMinimizer(NULL)
47 {
48  addParameter("Number of Iterations", CCopasiParameter::UINT, (unsigned C_INT32) 200);
49 // we no longer give the user choice of rng, we use the mersenne twister!
50 // but in DEBUG versions we should still have access to it
51 #ifdef COPASI_DEBUG
52  addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
53  addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
54 #endif
55 
56  initObjects();
57 }
58 
60  const CCopasiContainer * pParent):
61  COptMethod(src, pParent),
62  mIterations(0),
63  mPopulationSize(0),
64  mVariableSize(0),
65  mRefSet(0),
66  mRefSetVal(0),
67  mPool(0),
68  mPoolVal(0),
69  mPoolSize(0),
70  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
71  mIteration(0),
72  mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
73  mBestIndex(C_INVALID_INDEX),
74  mpRandom(NULL),
75  mpOptProblemLocal(NULL),
76  mpLocalMinimizer(NULL)
77 {
78  // remove eventual existing parameters from the release version.
79  initObjects();
80 }
81 
83 {cleanup();}
84 
85 // virtual
87 {
88 #ifndef COPASI_DEBUG
89  removeParameter("Random Number Generator");
90  removeParameter("Seed");
91 #endif
92  return true;
93 }
94 
96 {
98 }
99 
101 {
102  cleanup();
103 
104  size_t i;
105 
106  if (!COptMethod::initialize())
107  {
108  if (mpCallBack)
110 
111  return false;
112  }
113 
114  // get total number of iterations
115  mIterations = * getValue("Number of Iterations").pUINT;
116  // set current iteration to zero
117  mIteration = 0;
118 
119  if (mpCallBack)
120  mhIterations =
121  mpCallBack->addItem("Current Iteration",
122  mIteration,
123  & mIterations);
124 
125  mIteration++;
126 
127 #ifdef COPASI_DEBUG
128  mpRandom =
129  CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
130  * getValue("Seed").pUINT);
131 #else
132  // Use the default random number generator which is the Mersenne Twister
134 #endif
135 
136  mCloseValue = 0.001;
137 
138  // frequency of local optimiser; this is hardcoded as 20 per Jose Egea,
139  mLocalFreq = 20;
140 
141  // get number of variables in the problem
142  mVariableSize = mpOptItem->size();
143 
144  // calculate the number of individuals in population
145  mPopulationSize = (C_INT32) ceil(1.0 + sqrt(1.0 + 40.0 * mVariableSize) / 2.0);
146 
147  if (mPopulationSize % 2 != 0) mPopulationSize++;
148 
149  CFitProblem * pFitProblem = dynamic_cast< CFitProblem * >(mpOptProblem);
150 
151  if (pFitProblem != NULL)
152  {
153  // this is a least squares problem (param estimation)
154  // let's use our favorite lsq method
156  // the intermediate local minimizations use a rather relaxed tolerance
157  mpLocalMinimizer->setValue("Tolerance", (C_FLOAT64) 1.e-003);
158  // TODO: not sure if we should let this one go that long...
159  mpLocalMinimizer->setValue("Iteration Limit", (C_INT32) 2000);
160  }
161  else
162  {
163  // this is a generic optimisation problem
164  // let's use Hooke and Jeeves
166  // with a rather relaxed tolerance (1e-3) for intermediate minimizations
167  mpLocalMinimizer->setValue("Tolerance", (C_FLOAT64) 1.e-003);
168  mpLocalMinimizer->setValue("Iteration Limit", (C_INT32) 50);
169  mpLocalMinimizer->setValue("Rho", (C_FLOAT64) 0.2);
170  }
171 
172  // local minimization problem (starts as a copy of the current problem)
173  if (pFitProblem != NULL)
174  {
175  // this is a least squares problem (param estimation)
176  mpOptProblemLocal = new CFitProblem(*pFitProblem, getObjectParent());
177  }
178  else
179  {
180  // this is a generic optimisation problem
182  }
183 
184  // the local optimization method should not have a callback
186 
187  // set object parent (this is needed or else initialize() will fail)
189  // we also have to initialize the subtask
191  // initialize it
193  // no statistics to be calculated in the local problems
195  // do not randomize the initial values
198 
199  // create matrix for the RefSet (population)
200  mRefSet.resize(mPopulationSize);
201 
202  for (i = 0; i < mPopulationSize; i++)
203  mRefSet[i] = new CVector< C_FLOAT64 >(mVariableSize);
204 
205  // create vector for function values (of RefSet members)
206  mRefSetVal.resize(mPopulationSize);
207  mRefSetVal = std::numeric_limits<C_FLOAT64>::infinity();
208 
209  // create vector for counting stuck iterations
210  mStuck.resize(mPopulationSize);
211  mStuck = 0;
212 
213  // create matrix for the RefSet children
214  mChild.resize(mPopulationSize);
215 
216  for (i = 0; i < mPopulationSize; i++)
217  mChild[i] = new CVector< C_FLOAT64 >(mVariableSize);
218 
219  // create vector for function values (of child members)
220  mChildVal.resize(mPopulationSize);
221  mChildVal = std::numeric_limits<C_FLOAT64>::infinity();
222 
223  // we have not generated any children yet
224  mChildrenGenerated = false;
225 
226  // create matrix for the pool of diverse solutions
227  // this will also be used to store the initial and
228  // final positions of local optimizations
229  if (10 * mVariableSize > 2 * mIterations / mLocalFreq)
230  {
231  mPoolSize = 10 * mVariableSize;
232  }
233  else
234  {
236  }
237 
238  mPool.resize(mPoolSize);
239 
240  for (i = 0; i < mPoolSize; i++)
241  mPool[i] = new CVector< C_FLOAT64 >(mVariableSize);
242 
243  mPoolVal.resize(mPoolSize);
244  mPoolVal = std::numeric_limits<C_FLOAT64>::infinity();
245 
246  // best is infinity, so anything will improve it
247  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
248 
249  // array for frequencies
250  // initialized at 1
251  mFreq.resize(mVariableSize);
252 
253  for (i = 0; i < mVariableSize; i++)
254  {
255  mFreq[i] = new CVector< C_INT32 >(4);
256  *mFreq[i] = 1;
257  }
258 
259  // vector for probabilities
260  // initialized at 0
261  mProb.resize(4);
262  mProb = 0.0;
263 
264  return true;
265 }
266 
268 {
269  size_t i;
270 
271  pdelete(mpRandom);
272 
274 
276 
277  for (i = 0; i < mRefSet.size(); i++)
278  pdelete(mRefSet[i]);
279 
280  for (i = 0; i < mChild.size(); i++)
281  pdelete(mChild[i]);
282 
283  for (i = 0; i < mPool.size(); i++)
284  pdelete(mPool[i]);
285 
286  for (i = 0; i < mFreq.size(); i++)
287  pdelete(mFreq[i]);
288 
289  return true;
290 }
291 
292 // Find a local minimum
293 // solution has initial guess on entry, and solution on exit
294 // fval has value of objective function on exit
296 {
297  bool Running = true;
298  unsigned C_INT32 i;
299 
301 
302  // first we set up the problem
303  // (optmethod and optproblem already setup in initialization)
304  // let's get the list of parameters
305  std::vector<COptItem *> optitem = mpOptProblemLocal->getOptItemList();
306 
307  // and set them to the values passed in solution
308  for (i = 0; i < mVariableSize; i++)
309  {
310  optitem[i]->setStartValue(solution[i]);
311  }
312 
313  // reset the function counter of the local minimizer
315 
316  // run it
317  Running &= mpLocalMinimizer->optimise();
318  // add the function evaluations taken in local to the global problem
320  // pass the results on to the calling parameters
322 
323  for (i = 0; i < mVariableSize; i++)
324  {
325  solution[i] = mpOptProblemLocal->getSolutionVariables()[i];
326  }
327 
328  return Running;
329 }
330 
331 // evaluate the fitness of one individual
332 bool COptMethodSS::evaluate(const CVector< C_FLOAT64 > & /* individual */)
333 {
334  bool Running = true;
335 
336  // We do not need to check whether the parametric constraints are fulfilled
337  // since the parameters are created within the bounds.
338 
339  // evaluate the fitness
340  Running &= mpOptProblem->calculate();
341 
342  // check whether the functional constraints are fulfilled
344  mEvaluationValue = std::numeric_limits<C_FLOAT64>::infinity();
345  else
347 
348  return Running;
349 }
350 
351 // randomize the values of RefSet[i]
353 {
354  C_FLOAT64 mn, mx, la; // for boundaries of rnd
355  bool Running = true; // flag for invalid values
356 
357  for (C_INT32 j = 0; j < mVariableSize; j++)
358  {
359  // get pointers to appropriate elements (easier reading of code)
360  COptItem & OptItem = *(*mpOptItem)[j];
361  C_FLOAT64 & Sol = (*mRefSet[i])[j];
362  // calculate lower and upper bounds for this variable
363  mn = *OptItem.getLowerBoundValue();
364  mx = *OptItem.getUpperBoundValue();
365 
366  try
367  {
368  // calculate orders of magnitude of the interval
369  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
370 
371  // determine if linear or log scale
372  if ((mn < 0.0) || (mx <= 0.0))
373  Sol = mn + mpRandom->getRandomCC() * (mx - mn);
374  else
375  {
376  if (la < 1.8)
377  Sol = mn + mpRandom->getRandomCC() * (mx - mn);
378  else
379  Sol = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()))
380  + la * mpRandom->getRandomCC());
381  }
382  }
383  catch (...)
384  {
385  // if there were errors, let's just stay with the midpoint
386  Sol = (mx + mn) * 0.5;
387  }
388 
389  // force it to be within the bounds
390  switch (OptItem.checkConstraint(Sol))
391  {
392  case - 1:
393  Sol = *OptItem.getLowerBoundValue();
394  break;
395 
396  case 1:
397  Sol = *OptItem.getUpperBoundValue();
398  break;
399  }
400 
401  // We need to set the value here so that further checks take
402  // account of the value.
403  (*(*mpSetCalculateVariable)[j])(Sol);
404  }
405 
406  // calculate its fitness
407  Running = evaluate(*mRefSet[i]);
409  // reset the stuck flag
410  mStuck[i] = 1;
411  return Running;
412 }
413 
414 // initialise the population using the
415 // Diversification Generation Method
417 {
418  C_INT32 i, j, k, l; // counters
419  C_FLOAT64 mn, mx, la; // for boundaries of rnd
420  C_FLOAT64 sum; // to calculate a summation
421  C_FLOAT64 a; // to hold a random number
422  bool Running = true; // flag for invalid values
423 
424  // initialize all entries of the Pool.
425  // first 4 candidates as a latin hypercube
426  for (i = 0; (i < 4) && Running; i++)
427  {
428  for (j = 0; j < mVariableSize; j++)
429  {
430  // get pointers to appropriate elements (easier reading of code)
431  COptItem & OptItem = *(*mpOptItem)[j];
432  C_FLOAT64 & Sol = (*mPool[i])[j];
433  // calculate lower and upper bounds for this variable
434  mn = *OptItem.getLowerBoundValue();
435  mx = *OptItem.getUpperBoundValue();
436 
437  try
438  {
439  // calculate orders of magnitude of the interval
440  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
441 
442  // determine if linear or log scale
443  if ((mn < 0.0) || (mx <= 0.0))
444  Sol = mn + (mpRandom->getRandomCC() + (C_FLOAT64) i) * (mx - mn) * 0.25;
445  else
446  {
447  if (la < 1.8)
448  Sol = mn + (mpRandom->getRandomCC() + (C_FLOAT64) i) * (mx - mn) * 0.25;
449  else
450  Sol = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()))
451  + la * 0.25 * (mpRandom->getRandomCC() + (C_FLOAT64) i));
452  }
453  }
454  catch (...)
455  {
456  // TODO: this sounds a bit daft in this context, what else could be done, though?
457  Sol = (mx + mn) * 0.5;
458  }
459 
460  // force it to be within the bounds
461  switch (OptItem.checkConstraint(Sol))
462  {
463  case - 1:
464  Sol = *OptItem.getLowerBoundValue();
465  break;
466 
467  case 1:
468  Sol = *OptItem.getUpperBoundValue();
469  break;
470  }
471 
472  // We need to set the value here so that further checks take
473  // account of the value.
474  (*(*mpSetCalculateVariable)[j])(Sol);
475  }
476 
477  // calculate its fitness
478  Running &= evaluate(*mPool[i]);
480  }
481 
482  // next we add the initial guess from the user
483  for (j = 0; j < mVariableSize; j++)
484  {
485  COptItem & OptItem = *(*mpOptItem)[j];
486  C_FLOAT64 & Sol = (*mPool[i])[j];
487 
488  // get the vector of initial value
489  Sol = OptItem.getStartValue();
490 
491  // force it to be within the bounds
492  switch (OptItem.checkConstraint(Sol))
493  {
494  case - 1:
495  Sol = *OptItem.getLowerBoundValue();
496  break;
497 
498  case 1:
499  Sol = *OptItem.getUpperBoundValue();
500  break;
501  }
502 
503  // We need to set the value here so that further checks take
504  // account of the value.
505  (*(*mpSetCalculateVariable)[j])(Sol);
506  }
507 
508  // calculate its fitness
509  Running &= evaluate(*mPool[i]);
511  i++;
512 
513  // the remaining entries depend on probabilities
514  for (; (i < mPoolSize) && Running; i++)
515  {
516  for (j = 0; j < mVariableSize; j++)
517  {
518  // get pointers to appropriate elements (easier reading of code)
519  COptItem & OptItem = *(*mpOptItem)[j];
520  C_FLOAT64 & Sol = (*mPool[i])[j];
521  // calculate lower and upper bounds for this variable
522  mn = *OptItem.getLowerBoundValue();
523  mx = *OptItem.getUpperBoundValue();
524 
525  for (k = 0; k < 4; k++)
526  {
527  for (l = 0, sum = 0.0 ; l < 4; l++)
528  sum += 1.0 / (*mFreq[j])[l];
529 
530  mProb[k] = 1.0 / (*mFreq[j])[k] / sum;
531 
532  // we only want the cumulative probabilities
533  if (k > 0) mProb[k] += mProb[k - 1];
534  }
535 
536  a = mpRandom->getRandomCC();
537 
538  for (k = 0; k < 4; k++)
539  {
540  // note that the original is <= but numerically < is essentially the same and faster
541  if (a < mProb[k])
542  {
543  try
544  {
545  // calculate orders of magnitude of the interval
546  la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
547 
548  // determine if linear or log scale
549  if ((mn < 0.0) || (mx <= 0.0))
550  Sol = mn + (mpRandom->getRandomCC() + (C_FLOAT64) k) * (mx - mn) * 0.25;
551  else
552  {
553  if (la < 1.8)
554  Sol = mn + (mpRandom->getRandomCC() + (C_FLOAT64) k) * (mx - mn) * 0.25;
555  else
556  Sol = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()))
557  + la * 0.25 * (mpRandom->getRandomCC() + (C_FLOAT64) k));
558  }
559  }
560  catch (...)
561  {
562  // TODO: this sounds a bit daft in this context, what else could be done, though?
563  Sol = (mx + mn) * 0.5;
564  }
565 
566  // force it to be within the bounds
567  switch (OptItem.checkConstraint(Sol))
568  {
569  case - 1:
570  Sol = *OptItem.getLowerBoundValue();
571  break;
572 
573  case 1:
574  Sol = *OptItem.getUpperBoundValue();
575  break;
576  }
577 
578  // We need to set the value here so that further checks take
579  // account of the value.
580  (*(*mpSetCalculateVariable)[j])(Sol);
581  // increase the frequency
582  (*mFreq[j])[k] += 1;
583  break;
584  }
585  }
586  }
587 
588  // calculate its fitness
589  Running &= evaluate(*mPool[i]);
591  }
592 
593  // at this point the pool is formed
594  // now we partially sort the pool for the h = b/2 top elements,
595  // where b is mPopulationSize, This is done in place with heap sort
596  CVector< C_FLOAT64 > *tempvec;
597  C_FLOAT64 tempval;
598  C_INT32 parent, child;
599  C_INT32 h = mPopulationSize / 2;
600 
601  // top h are the heap, we first have to sort it
602  for (i = 1; i < h; i++)
603  {
604  // bubble-up
605  child = i;
606 
607  for (;;)
608  {
609  if (child == 0) break;
610 
611  parent = floor((double)(child - 1) / 2);
612 
613  if (mPoolVal[child] < mPoolVal[parent])
614  {
615  // swap with parent
616  tempval = mPoolVal[child];
617  mPoolVal[child] = mPoolVal[parent];
618  mPoolVal[parent] = tempval;
619  tempvec = mPool[child];
620  mPool[child] = mPool[parent];
621  mPool[parent] = tempvec;
622  // make parent the new child
623  child = parent;
624  }
625  else break;
626  }
627  }
628 
629  for (i = h; i < mPoolSize; i++)
630  {
631  child = 0;
632 
633  // check if this element is smaller than any of the leafs
634  for (size_t leaf = h / 2; leaf < h; leaf++)
635  {
636  if (mPoolVal[i] < mPoolVal[leaf])
637  {
638  // keep if this leaf is worse than previous
639  if (mPoolVal[child] < mPoolVal[leaf])
640  child = leaf;
641  }
642  }
643 
644  if (child == 0) continue;
645 
646  if (mPoolVal[i] < mPoolVal[child])
647  {
648  // swap i and h+j
649  tempval = mPoolVal[child];
650  mPoolVal[child] = mPoolVal[i];
651  mPoolVal[i] = tempval;
652  tempvec = mPool[child];
653  mPool[child] = mPool[i];
654  mPool[i] = tempvec;
655 
656  // now bubble-up
657  for (;;)
658  {
659  if (child == 0) break;
660 
661  parent = floor((double)(child - 1) / 2);
662 
663  if (mPoolVal[child] < mPoolVal[parent])
664  {
665  // swap with parent
666  tempval = mPoolVal[child];
667  mPoolVal[child] = mPoolVal[parent];
668  mPoolVal[parent] = tempval;
669  tempvec = mPool[child];
670  mPool[child] = mPool[parent];
671  mPool[parent] = tempvec;
672  // make parent the new child
673  child = parent;
674  }
675  else break;
676  }
677  }
678  }
679 
680  // since some leafs are not in order in the array, we now do
681  // a bubble sort (note: this is best case for bubble sort)
682  // we use j because we do not want to change the value of h
683  j = h;
684 
685  do
686  {
687  k = 0;
688 
689  for (i = 0; i <= j; i++)
690  {
691  if (mPoolVal[i] > mPoolVal[i + 1])
692  {
693  tempval = mPoolVal[i];
694  mPoolVal[i] = mPoolVal[i + 1];
695  mPoolVal[i + 1] = tempval;
696  tempvec = mPool[i];
697  mPool[i] = mPool[i + 1];
698  mPool[i + 1] = tempvec;
699  k = i;
700  }
701  }
702 
703  j = k;
704  }
705  while (j > 1);
706 
707  // at this point the pool is formed and partially sorted
708  // now we build the RefSet:
709  // 1st b/2 are the best ones in the Pool (sorted already)
710  // the 2nd b/2 are random (or later can be made diverse by
711  // maximising the Euclidean distances)
712  for (i = 0; i < mPopulationSize; i++)
713  {
714  (*mRefSet[i]) = (*mPool[i]); // copy the whole vector
715  mRefSetVal[i] = mPoolVal[i]; // keep the value
716  mStuck[i] = 1; // initialize the mStuck counter
717  }
718 
719  // RefSet needs to be fully sorted. Note that top half is sorted
720  // and we are garanteed that it contains the best values
721  // so it is only bottom half that needs sorting.
722  sortRefSet(h, mPopulationSize);
723  // we're done
724  return Running;
725 }
726 
727 // sort the RefSet and associated RefSetVal & Stuck
728 // between positions lower and upper.
730 {
731  C_INT32 i, j, k;
732  C_INT32 parent, child;
733  CVector< C_FLOAT64 > *tempvec;
734  C_FLOAT64 tempval;
735 
736  // Use heap sort
737  for (i = lower + 1; i < upper; i++)
738  {
739  // bubble-up
740  child = i;
741 
742  for (;;)
743  {
744  if (child == 0) break;
745 
746  parent = floor((double)(child - 1) / 2);
747 
748  if (mRefSetVal[child] < mRefSetVal[parent])
749  {
750  // swap with parent
751  tempval = mRefSetVal[child];
752  mRefSetVal[child] = mRefSetVal[parent];
753  mRefSetVal[parent] = tempval;
754  tempval = mStuck[child];
755  mStuck[child] = mStuck[parent];
756  mStuck[parent] = tempval;
757  tempvec = mRefSet[child];
758  mRefSet[child] = mRefSet[parent];
759  mRefSet[parent] = tempvec;
760  // make parent the new child
761  child = parent;
762  }
763  else break;
764  }
765  }
766 
767  // since some leafs are not in order in the array, we now do
768  // a bubble sort (note: this is best case for bubble sort)
769  j = upper - 1;
770 
771  do
772  {
773  k = lower;
774 
775  for (i = lower; i < j; i++)
776  {
777  if (mRefSetVal[i] > mRefSetVal[i + 1])
778  {
779  tempval = mRefSetVal[i];
780  mRefSetVal[i] = mRefSetVal[i + 1];
781  mRefSetVal[i + 1] = tempval;
782  tempval = mStuck[i];
783  mStuck[i] = mStuck[i + 1];
784  mStuck[i + 1] = tempval;
785  tempvec = mRefSet[i];
786  mRefSet[i] = mRefSet[i + 1];
787  mRefSet[i + 1] = tempvec;
788  k = i;
789  }
790  }
791 
792  j = k;
793  }
794  while (j > lower);
795 }
796 
797 // check if all the indexes of a Child member are closer to
798 // the indexes of a Pool member than a certain relative distance
800 {
801  C_FLOAT64 mx;
802 
803  for (C_INT32 k = 0; k < mVariableSize; k++)
804  {
805  mx = (fabs((*mChild[i])[k]) + fabs((*mPool[j])[k])) / 2.0;
806 
807  if (fabs((*mChild[i])[k] - (*mPool[j])[k]) / mx > dist) return false;
808  }
809 
810  return true;
811 }
812 
813 // check if all the indexes of two refset members
814 // are closer than a certain relative distance
816 {
817  C_FLOAT64 mx;
818 
819  for (C_INT32 k = 0; k < mVariableSize; k++)
820  {
821  mx = (fabs((*mRefSet[i])[k]) + fabs((*mRefSet[j])[k])) / 2.0;
822 
823  if (fabs((*mRefSet[i])[k] - (*mRefSet[j])[k]) / mx > dist) return false;
824  }
825 
826  return true;
827 }
828 
829 // combine individuals in the RefSet two by two
830 // this is a sort of (1+1)-ES strategy
832 {
833  C_INT32 i, j, k, l; // counters
834  C_FLOAT64 mn, mx; // for bounds on parameters
835  C_FLOAT64 beta; // bias
836  C_FLOAT64 la; // for orders of magnitude
837  C_FLOAT64 alpha; // 1 or -1
838  C_FLOAT64 bm2; // b-2
839  C_FLOAT64 omatb; // (1-alpha*beta)*0.5
840  C_FLOAT64 dd; // (x(i) - x(j) ) / 2 * (1-alpha*beta)
841  C_FLOAT64 c1, c2; // coordinates of the edges of hyperectangle
842  CVector< C_FLOAT64 > xnew, xpr;
843  C_FLOAT64 xnewval, xprval; // to hold temp value of "parent" in go-beyond strategy
844  C_FLOAT64 lambda = 1.0; // damping factor for go-beyond strategy
845  C_INT32 improvement; // count iterations with improvement in go-beyond strategy
846  bool Running = true; // flag for invalid values
847 
848  // make xnew large enough
849  xnew.resize(mVariableSize);
850  xpr.resize(mVariableSize);
851  // calculate a constant
852  bm2 = C_FLOAT64(mPopulationSize) - 2.0;
853  // signal no children yet
854  mChildrenGenerated = false;
855 
856  // generate children for each member of the population
857  for (i = 0; (i < mPopulationSize) && Running; i++)
858  {
859  // keep the parent value in childval[i] so that we only accept better than that
860  mChildVal[i] = mRefSetVal[i];
861 
862  for (j = 0; j < mPopulationSize; j++)
863  {
864  // no self-reproduction...
865  if (i != j)
866  {
867  if (i < j) alpha = 1.0; else alpha = -1.0;
868 
869  beta = (fabs(C_FLOAT64(j) - C_FLOAT64(i)) - 1.0) / bm2;
870  omatb = (1.0 + alpha * beta) * 0.5;
871 
872  // generate a child
873  for (k = 0; k < mVariableSize; k++)
874  {
875  // get the bounds of this parameter
876  COptItem & OptItem = *(*mpOptItem)[k];
877  mn = *OptItem.getLowerBoundValue();
878  mx = *OptItem.getUpperBoundValue();
879 
880  try
881  {
882  // calculate orders of magnitude of the interval
883  if (((*mRefSet[i])[k] > 0.0) && ((*mRefSet[j])[k] > 0.0))
884  {
885  la = log10((*mRefSet[i])[k]) - log10(std::max((*mRefSet[j])[k], std::numeric_limits< C_FLOAT64 >::min()));
886  }
887  else
888  la = 1.0;
889 
890  dd = ((*mRefSet[i])[k] - (*mRefSet[j])[k]) * omatb;
891  // one of the box limits
892  c1 = (*mRefSet[i])[k] - dd;
893 
894  // force it to be within the bounds
895  switch (OptItem.checkConstraint(c1))
896  {
897  case -1:
898  c1 = mn;
899  break;
900 
901  case 1:
902  c1 = mx;
903  break;
904  }
905 
906  // the other box limit
907  c2 = (*mRefSet[i])[k] + dd;
908 
909  // force it to be within the bounds
910  switch (OptItem.checkConstraint(c2))
911  {
912  case -1:
913  c2 = mn;
914  break;
915 
916  case 1:
917  c2 = mx;
918  break;
919  }
920 
921  xnew[k] = c1 + (c2 - c1) * mpRandom->getRandomCC();
922  }
923  catch (...)
924  {
925  // if something failed leave the value intact
926  xnew[k] = (*mRefSet[i])[k];
927  }
928 
929  // We need to set the value here so that further checks take
930  // account of the value.
931  (*(*mpSetCalculateVariable)[k])(xnew[k]);
932  }
933 
934  // calculate the child's fitness
935  Running &= evaluate(xnew);
936 
937  // keep it if it is better than the previous one
938  if (mEvaluationValue < mChildVal[i])
939  {
940  // keep a copy of this vector in mChild
941  (*mChild[i]) = xnew;
942  // and the fitness value
944  // signal that child is better than parent
945  mStuck[i] = 0;
946  // signal we have generated a child (improvement)
947  mChildrenGenerated = true;
948  }
949  }
950  }
951 
952  // now we apply the go-beyond strategy for
953  // cases where the child was an improvement
954  if (mStuck[i] == 0)
955  {
956  // copy the parent
957  xpr = (*mRefSet[i]);
958  xprval = mRefSetVal[i];
959  lambda = 1.0; // this is really 1/lambda so we can use mult rather than div
960  improvement = 1;
961 
962  // while newval < childval
963  for (; ;)
964  {
965  for (k = 0; k < mVariableSize; k++)
966  {
967  dd = (xpr[i] - (*mChild[i])[k]) * lambda;
968  xnew[k] = (*mChild[i])[k] + dd * mpRandom->getRandomCC();
969  // get the bounds of this parameter
970  COptItem & OptItem = *(*mpOptItem)[k];
971 
972  // put it on the bounds if it had exceeded them
973  switch (OptItem.checkConstraint(xnew[k]))
974  {
975  case -1:
976  xnew[k] = *OptItem.getLowerBoundValue();
977  break;
978 
979  case 1:
980  xnew[k] = *OptItem.getUpperBoundValue();
981  break;
982  }
983 
984  // We need to set the value here so that further checks take
985  // account of the value.
986  (*(*mpSetCalculateVariable)[k])(xnew[k]);
987  }
988 
989  // calculate the child's fitness
990  Running &= evaluate(xnew);
991  xnewval = mEvaluationValue;
992 
993  // if there was no improvement we finish here => exit for(;;)
994  if (mChildVal[i] <= xnewval) break;
995 
996  // old child becomes parent
997  xpr = (*mChild[i]);
998  xprval = mChildVal[i];
999  // new child becomes child
1000  (*mChild[i]) = xnew;
1001  mChildVal[i] = xnewval;
1002  // mark improvement
1003  improvement++;
1004 
1005  if (improvement == 2)
1006  {
1007  lambda *= 0.5;
1008  improvement = 0;
1009  }
1010  }
1011  }
1012  }
1013 
1014  return Running;
1015 }
1016 
1018 {
1019  C_INT32 i, best;
1020  C_FLOAT64 bestVal, fvalmin;
1021  bool Running = true;
1022 
1023  // signal nothing found yet
1024  best = -1;
1025  bestVal = std::numeric_limits<C_FLOAT64>::infinity();
1026 
1027  // find the best child
1028  for (i = 0; i < mPopulationSize; i++)
1029  {
1030  if ((mStuck[i] == 0) && (mChildVal[i] < bestVal))
1031  {
1032  bestVal = mChildVal[i];
1033  best = i;
1034  }
1035  }
1036 
1037  // no child in this iteration? exit now
1038  if (best == -1) return true;
1039 
1040  // check if this child is not close to previous ones
1041  for (i = 0; i < mLocalStored; i++)
1042  {
1043  // is the other one like me?
1044  if (closerChild(best, i, mCloseValue))
1045  {
1046  // it is too close, exit now
1047  return true;
1048  }
1049  }
1050 
1051  // store the initial position
1052  *(mPool[mLocalStored]) = *(mChild[best]);
1053  mPoolVal[mLocalStored] = mChildVal[best]; //bestVal;
1054  mLocalStored++;
1055 
1056  // do local minimization on it
1057  Running &= localmin(*(mChild[best]), mChildVal[best]);
1058 
1059  // store the result
1060  *(mPool[mLocalStored]) = *(mChild[best]);
1061  mPoolVal[mLocalStored] = mChildVal[best];
1062  mLocalStored++;
1063 
1064  // clear the local optimisation counter
1065  mLocalIter = 1;
1066 
1067  return Running;
1068 }
1069 
1071 {
1072  bool Running = true;
1073  bool needsort;
1074  size_t i, j;
1075  C_FLOAT64 mx, mn, la;
1076 
1077  if (!initialize())
1078  {
1079  // initialisation failed, we exit
1080  if (mpCallBack)
1082 
1083  return false;
1084  }
1085 
1086  mIteration = 0;
1087 
1088  // create the Pool of diverse candidate solutions
1089  Running &= creation();
1090 
1091  // best value is (always) at position zero
1092  // store that value
1093  mBestValue = mRefSetVal[0];
1094  // set it upstream
1095  Running &= mpOptProblem->setSolution(mBestValue, *mRefSet[0]);
1096  // We found a new best value let's report it.
1098 
1099  // test if the user wants to stop, and do so if needed
1100  if (!Running)
1101  {
1102  if (mpCallBack)
1104 
1105  cleanup();
1106  return true;
1107  }
1108 
1109  // mPool is now going to be used to keep track of initial and final
1110  // points of local minimizations (to avoid running them more than once)
1112  // reset the number of stored minimizations
1113  mLocalStored = 0;
1114  // reset the counter for local minimisation
1115  mLocalIter = 1;
1116 
1117  // run the mIterations (and count the creation as being the first)
1118  for (mIteration = 1; mIteration < mIterations && Running; mIteration++)
1119  {
1120  // check for stagnation or similarity
1121  needsort = false;
1122 
1123  for (i = 0; i < mPopulationSize; i++)
1124  {
1125  // are we stuck? (20 iterations)
1126  if (mStuck[i] == 19)
1127  {
1128  // substitute this one by a random guess
1129  Running &= randomize(i);
1130  needsort = true;
1131  mStuck[i] = 1;
1132  }
1133  else
1134  {
1135  // check if another RefSet member is similar to us (relative dist 0.1%)
1136  for (j = i + 1; j < mPopulationSize; j++)
1137  {
1138  // is the other one like me?
1139  if (closerRefSet(i, j, mCloseValue))
1140  //if (distRefSet(i, j) < 0.01)
1141  {
1142  // randomize the other one because it has worse value
1143  Running &= randomize(j);
1144  needsort = true;
1145  mStuck[j] = 1;
1146  }
1147  }
1148  }
1149  }
1150 
1151  // sort the RefSet if needed
1152  if (needsort) sortRefSet(0, mPopulationSize);
1153 
1154  // create children by combination
1155  Running &= combination();
1156 
1157  // check if we have to run a local search
1159  {
1160  // carry out a local search
1161  Running &= childLocalMin();
1162  }
1163  else
1164  {
1165  // count this
1166  mLocalIter++;
1167  }
1168 
1169  // substitute the parents for children or increment stuck counter
1170  needsort = false;
1171 
1172  for (i = 0; i < mPopulationSize; i++)
1173  {
1174  // check if child was better than parent
1175  if (mStuck[i] == 0)
1176  {
1177  // copy the child into the population
1178  (*mRefSet[i]) = (*mChild[i]);
1179  // keep its obj funct value
1180  mRefSetVal[i] = mChildVal[i];
1181  // and reset the stuck counter
1182  mStuck[i] = 1;
1183  needsort = true;
1184  }
1185  else
1186  {
1187  // nothing to do, increment parent's stuck counters
1188  mStuck[i]++;
1189  }
1190  }
1191 
1192  // sort the RefSet if needed
1193  if (needsort) sortRefSet(0, mPopulationSize);
1194 
1195  // have we made any progress?
1196  if (mRefSetVal[0] < mBestValue)
1197  {
1198  // and store that value
1199  mBestValue = mRefSetVal[0];
1200  Running &= mpOptProblem->setSolution(mBestValue, *mRefSet[0]);
1201  // We found a new best value lets report it.
1203  }
1204 
1205  if (mpCallBack)
1206  Running &= mpCallBack->progressItem(mhIterations);
1207  }
1208 
1209  // end of loop for iterations
1210 
1211  // the best ever might not be what is on position 0, so bring it back
1213 
1214  // now let's do a final local minimisation with a tighter tolerance
1215 
1216  if (Running)
1217  {
1218  mpLocalMinimizer->setValue("Tolerance", (C_FLOAT64) 1.e-006);
1219  Running &= localmin(*(mRefSet[0]), mRefSetVal[0]);
1220  }
1221 
1222  // has it improved?
1223  if (mRefSetVal[0] < mBestValue)
1224  {
1225  // and store that value
1226  mBestValue = mRefSetVal[0];
1227  Running &= mpOptProblem->setSolution(mBestValue, *mRefSet[0]);
1228  // We found a new best value lets report it.
1230  }
1231 
1232  if (mpCallBack)
1234 
1235  cleanup();
1236  return true;
1237 }
1238 
1239 #ifdef DEBUG_OPT
1240 bool COptMethodSS::serializeprob(void)
1241 {
1242  std::ofstream ofile;
1243  // open the file for output, in append mode
1244  ofile.open("ssprob.txt", std::ios::out); // | std::ios::app );
1245 
1246  if (! ofile.is_open())
1247  {
1248  std::cerr << "error opening file \'ssprob.txt\'" << std::endl;
1249  return false;
1250  }
1251 
1252  ofile << mProb << std::endl;
1253  ofile.close();
1254  return true;
1255 }
1256 // serialize the pool to a file for debugging purposes
1257 bool COptMethodSS::serializepool(size_t first, size_t last)
1258 {
1259  size_t Last = std::min(last, (size_t) mPoolSize);
1260 
1261  size_t i;
1262  size_t j;
1263 
1264  std::ofstream ofile;
1265 
1266  // open the file for output, in append mode
1267  ofile.open("sspool.txt", std::ios::out | std::ios::app);
1268 
1269  if (! ofile.is_open())
1270  {
1271  std::cerr << "error opening file \'sspool.txt\'" << std::endl;
1272  return false;
1273  }
1274 
1275  for (i = first; i < Last; i++)
1276  {
1277  for (j = 0; j < mVariableSize; j++)
1278  {
1279  C_FLOAT64 & mut = (*mPool[i])[j];
1280  // print this parameter
1281  ofile << mut << "\t";
1282  }
1283 
1284  // print the fitness of the individual
1285  ofile << mPoolVal[i] << std::endl;
1286  }
1287 
1288  ofile << std::endl;
1289  ofile.close();
1290  return true;
1291 }
1292 // write informative messages about the progress of the refset
1293 bool COptMethodSS::inforefset(C_INT32 type, C_INT32 element)
1294 {
1295  std::ofstream ofile;
1296 
1297  // open the file for output, in append mode
1298  ofile.open("ssrefset.txt", std::ios::out | std::ios::app);
1299 
1300  if (! ofile.is_open())
1301  {
1302  std::cerr << "error opening file \'ssrefset.txt\'" << std::endl;
1303  return false;
1304  }
1305 
1306  switch (type)
1307  {
1308  case 1: ofile << "element " << element << " improved in combination" << std::endl; break;
1309 
1310  case 2: ofile << "element " << element << " improved in go-beyond" << std::endl; break;
1311 
1312  case 3: ofile << "No element improved in iteration " << element << std::endl ; break;
1313 
1314  case 4: ofile << "element " << element << " randomized, too close to another" << std::endl; break;
1315 
1316  case 5: ofile << "element " << element << " randomized, was stuck" << std::endl; break;
1317 
1318  case 6: ofile << "child" << element << " too close to previous solution, no local min" << std::endl; break;
1319 
1320  case 7: ofile << "c1 is NaN (element" << element << ")" << std::endl; break;
1321 
1322  case 8: ofile << "c2 is NaN (element" << element << ")" << std::endl; break;
1323 
1324  case 9: ofile << "xnew[k] is NaN (element" << element << ")" << std::endl; break;
1325 
1326  case 10: ofile << "Children of " << element << std::endl; break;
1327 
1328  case 11: ofile << "Local minimization at value " << element << std::endl; break;
1329  }
1330 
1331  ofile.close();
1332  return true;
1333 }
1334 
1335 // serialize the population to a file for debugging purposes
1337 {
1338  C_INT32 Last = std::min(last, mPopulationSize);
1339  C_INT32 i, j;
1340  std::ofstream ofile;
1341 
1342  // open the file for output, in append mode
1343  ofile.open("ssrefset.txt", std::ios::out | std::ios::app);
1344 
1345  if (! ofile.is_open())
1346  {
1347  std::cerr << "error opening file \'ssrefset.txt\'" << std::endl;
1348  return false;
1349  }
1350 
1351  ofile << std::endl << "Refset" << std::endl;
1352 
1353  for (i = first; i < Last; i++)
1354  {
1355  for (j = 0; j < mVariableSize; j++)
1356  {
1357  C_FLOAT64 & mut = (*mRefSet[i])[j];
1358  // print this parameter
1359  ofile << mut << "\t";
1360  }
1361 
1362  // print the fitness of the individual
1363  ofile << mRefSetVal[i] << std::endl;
1364  }
1365 
1366  ofile << std::endl;
1367  ofile.close();
1368  return true;
1369 }
1370 // serialize the population to a file for debugging purposes
1372 {
1373  C_INT32 i;
1374  std::ofstream ofile;
1375 
1376  // open the file for output, in append mode
1377  ofile.open("ssrefset.txt", std::ios::out | std::ios::app);
1378 
1379  if (! ofile.is_open())
1380  {
1381  std::cerr << "error opening file \'ssrefset.txt\'" << std::endl;
1382  return false;
1383  }
1384 
1385  for (i = 0; i < mVariableSize; i++)
1386  {
1387  // print this parameter
1388  ofile << x[i] << "\t";
1389  }
1390 
1391  // print the fitness of the individual
1392  ofile << xval << std::endl;
1393 
1394  ofile.close();
1395  return true;
1396 }
1397 // serialize the population to a file for debugging purposes
1398 bool COptMethodSS::serializechildren(size_t first, size_t last)
1399 {
1400  size_t Last = std::min(last, (size_t) mPopulationSize);
1401 
1402  size_t i;
1403  size_t j;
1404 
1405  std::ofstream ofile;
1406 
1407  // open the file for output, in append mode
1408  ofile.open("sschild.txt", std::ios::out | std::ios::app);
1409 
1410  if (! ofile.is_open())
1411  {
1412  std::cerr << "error opening file \'sschild.txt\'" << std::endl;
1413  return false;
1414  }
1415 
1416  for (i = first; i < Last; i++)
1417  {
1418  for (j = 0; j < mVariableSize; j++)
1419  {
1420  C_FLOAT64 & mut = (*mChild[i])[j];
1421  // print this parameter
1422  ofile << mut << "\t";
1423  }
1424 
1425  // print the fitness of the individual
1426  ofile << mChildVal[i] << std::endl;
1427  }
1428 
1429  ofile << std::endl;
1430  ofile.close();
1431  return true;
1432 }
1433 #endif
virtual bool setObjectParent(const CCopasiContainer *pParent)
C_FLOAT64 mCloseValue
Definition: COptMethodSS.h:289
std::vector< CVector< C_INT32 > * > mFreq
Definition: COptMethodSS.h:259
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
COptTask * mpParentTask
Definition: COptMethod.h:58
#define pdelete(p)
Definition: copasi.h:215
void incrementEvaluations(unsigned C_INT32 increment)
virtual bool initialize()
Definition: COptMethod.cpp:189
bool closerRefSet(C_INT32 i, C_INT32 j, C_FLOAT64 dist)
unsigned C_INT32 mIterations
Definition: COptMethodSS.h:182
bool serializechildren(size_t first, size_t last)
CVector< C_FLOAT64 > mProb
Definition: COptMethodSS.h:264
void setProblem(COptProblem *problem)
Definition: COptMethod.cpp:170
virtual bool elevateChildren()
static COptMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::RandomSearch)
Definition: COptMethod.cpp:50
bool creation(void)
bool serializeprob(void)
virtual bool initialize()
CRandom * mpRandom
Definition: COptMethodSS.h:299
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
void sortRefSet(C_INT32 lower, C_INT32 upper)
virtual void output(const COutputInterface::Activity &activity)
CVector< C_FLOAT64 > mChildVal
Definition: COptMethodSS.h:238
std::vector< CVector< C_FLOAT64 > * > mChild
Definition: COptMethodSS.h:233
#define C_INT32
Definition: copasi.h:90
bool inforefset(C_INT32 type, C_INT32 element)
bool mChildrenGenerated
Definition: COptMethodSS.h:212
CVector< C_INT32 > mStuck
Definition: COptMethodSS.h:228
virtual bool progressItem(const size_t &handle)
virtual bool cleanup()
void resetEvaluations()
C_INT32 mPopulationSize
Definition: COptMethodSS.h:187
bool evaluate(const CVector< C_FLOAT64 > &individual)
virtual bool setCallBack(CProcessReport *pCallBack)
bool removeParameter(const std::string &name)
std::vector< CVector< C_FLOAT64 > * > mPool
Definition: COptMethodSS.h:244
const C_FLOAT64 & getSolutionValue() const
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
virtual bool calculate()
virtual ~COptMethodSS()
virtual C_FLOAT64 getRandomCC()
Definition: CRandom.cpp:235
bool childLocalMin(void)
const std::vector< COptItem * > & getOptItemList() const
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
COptMethod * mpLocalMinimizer
Definition: COptMethodSS.h:309
unsigned C_INT32 mVariableSize
Definition: COptMethodSS.h:192
size_t mhIterations
Definition: COptMethodSS.h:294
const Value & getValue() const
bool localmin(CVector< C_FLOAT64 > &solution, C_FLOAT64 &fval)
COptMethodSS(const COptMethodSS &src, const CCopasiContainer *pParent=NULL)
COptProblem * mpOptProblemLocal
Definition: COptMethodSS.h:304
bool setValue(const std::string &name, const CType &value)
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
const CVector< C_FLOAT64 > & getSolutionVariables() const
unsigned C_INT32 * pUINT
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
bool serializevector(CVector< C_FLOAT64 > x, C_FLOAT64 xval)
virtual bool finishItem(const size_t &handle)
virtual bool initializeSubtaskBeforeOutput()
CVector< C_FLOAT64 > mPoolVal
Definition: COptMethodSS.h:249
bool randomize(C_INT32 i)
unsigned C_INT32 mLocalStored
Definition: COptMethodSS.h:207
virtual bool checkFunctionalConstraints()
bool closerChild(C_INT32 i, C_INT32 j, C_FLOAT64 dist)
virtual bool optimise()
Definition: COptMethod.cpp:184
const C_FLOAT64 & getStartValue() const
Definition: COptItem.cpp:199
void setRandomizeStartValues(const bool &randomize)
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
bool serializerefset(C_INT32 first, C_INT32 last)
unsigned C_INT32 mLocalIter
Definition: COptMethodSS.h:202
#define C_FLOAT64
Definition: copasi.h:92
CVector< C_FLOAT64 > mRefSetVal
Definition: COptMethodSS.h:222
size_t mPoolSize
Definition: COptMethodSS.h:254
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
bool serializepool(size_t first, size_t last)
virtual bool initialize()
bool combination(void)
bool addParameter(const CCopasiParameter &parameter)
C_FLOAT64 mEvaluationValue
Definition: COptMethodSS.h:269
unsigned C_INT32 mIteration
Definition: COptMethodSS.h:274
std::vector< CVector< C_FLOAT64 > * > mRefSet
Definition: COptMethodSS.h:217
void initObjects()
const unsigned C_INT32 & getFunctionEvaluations() const
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
C_FLOAT64 mBestValue
Definition: COptMethodSS.h:279
void setCalculateStatistics(const bool &calculate)
virtual bool optimise()
CCopasiContainer * getObjectParent() const
unsigned C_INT32 mLocalFreq
Definition: COptMethodSS.h:197
#define max(a, b)
Definition: f2c.h:176