COPASI API  4.16.103
COptMethodNelderMead.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) 2006 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 // adapted by Stefan Hoops, August 2006, from the original Gepasi file
16 // nelmea.cpp :
17 
18 #include <cmath>
19 
20 #include "copasi.h"
21 
22 #include "COptMethodNelderMead.h"
23 #include "COptProblem.h"
24 #include "COptItem.h"
25 #include "COptTask.h"
26 
29 
31  COptMethod(CCopasiTask::optimization, CCopasiMethod::NelderMead, pParent)
32 {
33  addParameter("Iteration Limit", CCopasiParameter::UINT, (unsigned C_INT32) 200);
34  addParameter("Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.e-005);
36 
37  initObjects();
38 }
39 
41  const CCopasiContainer * pParent):
42  COptMethod(src, pParent)
43 {initObjects();}
44 
46 {cleanup();}
47 
49 {
51 }
52 
53 /*
54  The code for COptMethodNelderMead::optimize() is derived from nelmin.c
55  developed by:
56 
57  Peter & Nigel,
58  Design Software,
59  ----------------------
60  42 Gubberley St,
61  Kenmore, 4069,
62  Australia.
63  */
64 /************************************************/
65 /* */
66 /* CMATH. Copyright (c) 1989 Design Software */
67 /* */
68 /************************************************/
69 
70 /*
71  Purpose ....
72  -------
73  Find the minimum value of a user-specified function.
74 
75  Input ...
76  -----
77  Fnelmin: User specified function to be minimized. NOTE that
78  the function name must be passed.
79  double Fnelmin (n, x)
80  int n;
81  double x[];
82  n : INTEGER. The number of variables over which we are
83  minimizing.
84  xmin : Array [n] of double; contains the coordinates
85  of the starting point. (see note 5)
86  reqmin : The terminating limit for the variance of function values.
87  step : Array [n] of double; determines the
88  size and shape of the initial simplex.
89  The relative magnitudes of its N elements
90  should reflect the units of the N variables.
91  konvge : The convergence check is carried out every KONVGE
92  iterations.
93  kcount : Maximum number of function evaluations.
94  reltol : relative tolerance on minimum check (see note 6)
95  abstol : absolute tolerance on minimum check (see note 6)
96 
97  Output ...
98  ------
99  xmin : Array [id] of double; contains the
100  coordinates of the minimum.
101  ynewlo : The minimum value of the function.
102  icount : Function evaluations performed.
103  numres : Number of restarts.
104  ifault : return status flag ...
105  = 0 No error.
106  = 1 REQMIN, N or KONVGE illegal value.
107  = 2 Termination because KCOUNT exceeded before convergence.
108  = 3 could not allocate memory for workspace
109 
110  This C code written by ... Peter & Nigel,
111  ---------------------- Design Software,
112  42 Gubberley St,
113  Kenmore, 4069,
114  Australia.
115 
116  Version ... 1.0 November 1987
117  ------- 2.0 March 1989 reltol, abstol added
118  zero subscripts implemented
119  2.1 15 April 1989 workspace allocated
120  start[] removed from arguements
121  2.2 28 April 1989 n added to Fnelmin()
122  2.3 7 August 89 fixed true min. check
123  2.4 12 Dec 1989 fixed memory (de)allocation
124 
125  Notes ...
126  -----
127  1. This algorithm is a modified version of ..
128  Algorithm AS 47 Applied Statistics (J. R. Statist. Soc. C),
129  (1971) VOL.20. NO.3
130 
131  2. The data values of RCOEFF (reflection), ECOEFF (extension),
132  and CCOEFF (contraction) can be changed by rewriting the
133  assignment statements below. The values currently set are 1.0,
134  2.0 and 0.5 respectively. These values have been recommended
135  by Nelder and Mead as being best for a general situation.
136 
137  3. The value of REQMIN must be set larger than the rounding
138  error in computing the variance at the minimum. This holds
139  for both Double and Single precision calculations.
140  Chambers and Ertel recommend that REQMIN be set to a value
141  somewhere between the required accuracy of the function at
142  the minimum and the square of this value.
143 
144  4. Note that the elements [0 .. n-1] are utilized in this
145  version of the routine.
146 
147  6. reltol and abstol should be set to zero for well behaved
148  functions where restarts are not a problem.
149 
150  7. References ...
151  Nelder, J.A. and Mead, R. (1965) A simplex method for function
152  minimization. Computer J. 7,308-313
153 
154  O'Neill, R. (1971) Algorithm AS47. Function minimization
155  using a simplex algorithm. Appl. Statist. 20,338-345.
156 
157  Chambers, J.M. and Ertel, J.E. (1974) Remark AS R11.
158  Appl. Statist. 23,250-251.
159 
160  Olsson, D.M. and Nelson, L.S. (1975) The Nelder - Mead
161  simplex procedure for function minimization.
162  Technometrics 17,45-51. (Examples of use.)
163 
164  Benyon, P.R. (1976) Remark AS R15. Appl. Statist. 25,97.
165 
166  Hill, I.D. (1978) Remark AS R28. Appl. Statist. 27,380-382.
167 
168  ----------------------------------------------------------------------
169  */
170 
172 {
173  if (!initialize())
174  {
175  if (mpCallBack)
177 
178  return false;
179  }
180 
181  // set tolerances for local minima test to zero
182  C_FLOAT64 abstol, reltol;
183  abstol = reltol = 0.0;
184 
185  // test convergence every iteration
186  size_t konvge;
187  konvge = 1;
188 
189  size_t iBest, iWorst;
190 
191  // double zero, half, one, delta;
192  C_FLOAT64 ccoeff, rcoeff, ecoeff;
193 
194  /* ---- reflection, extension and contraction coefficients. ---- */
195  /* parameters from Nelder and Mead */
196  //rcoeff = 1.0;
197  //ccoeff = 0.5;
198  //ecoeff = 2.0;
199 
200  /* Parameters from Parkinson and Hutchinson */
201  rcoeff = 2.0;
202  ccoeff = 0.25;
203  ecoeff = 2.5;
204 
205  C_FLOAT64 z, sum, ylo;
206  C_FLOAT64 curmin, del, x, yhi;
207  C_FLOAT64 small;
208  C_FLOAT64 factor = 0.8;
209 
210  size_t jcount, np1, i, j;
211 
212  bool ok, reflok, extnok, contok, quit, found; /* boolean variables */
213 
214  // initial point is first guess but we have to make sure that we
215  // are within the parameter domain
216  for (i = 0; i < mVariableSize; i++)
217  {
218  const COptItem & OptItem = *(*mpOptItem)[i];
219 
220  switch (OptItem.checkConstraint(OptItem.getStartValue()))
221  {
222  case - 1:
223  mCurrent[i] = *OptItem.getLowerBoundValue();
224  break;
225 
226  case 1:
227  mCurrent[i] = *OptItem.getUpperBoundValue();
228  break;
229 
230  case 0:
231  mCurrent[i] = OptItem.getStartValue();
232  break;
233  }
234 
235  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
236 
237  // set the magnitude of each parameter
238  mStep[i] = (*OptItem.getUpperBoundValue() - *OptItem.getLowerBoundValue()) / mScale;
239  }
240 
241  evaluate();
242 
243  if (!isnan(mEvaluationValue))
244  {
245  // and store that value
248 
249  // We found a new best value lets report it.
251  }
252 
253  quit = false;
254  found = false;
255 
256  np1 = mVariableSize + 1;
257  del = 1.0;
258 
259 First:
260 
261  /* ---- Construct the initial simplex. ---- */
262  for (i = 0; i < mVariableSize; ++i)
263  mSimplex[i][mVariableSize] = mCurrent[i];
264 
266 
267  for (j = 0; j < mVariableSize && mContinue; ++j)
268  {
269  x = mCurrent[j];
270 
271  // -- move a little in one dimension only --
272  mCurrent[j] += (mStep[j] * del);
273 
274  // Check constraint
275  const COptItem & OptItem = *(*mpOptItem)[j];
276 
277  if (mCurrent[j] > *OptItem.getUpperBoundValue())
278  mCurrent[j] = x - (mStep[j] * del);
279 
280  // Store the simplex corner
281  for (i = 0; i < mVariableSize; ++i)
282  mSimplex[i][j] = mCurrent[i];
283 
284  // Calculate the value for the corner
285  (*(*mpSetCalculateVariable)[j])(mCurrent[j]);
286  mValue[j] = evaluate();
287 
289  {
290  // and store that value
293 
294  // We found a new best value lets report it.
296  }
297 
298  // Reset
299  mCurrent[j] = x;
300  (*(*mpSetCalculateVariable)[j])(mCurrent[j]);
301  }
302 
303  /* ---- Simplex construction complete; now do some work. ---- */
304 
305  while (!found && !quit && mContinue)
306  {
307  /* ---- take some steps ---- */
308 
309  for (jcount = 1; jcount <= konvge; ++jcount)
310  {
311  /* ---- take a single step ---- */
312 
313  /* ---- find highest and lowest y values. yhi (=y[iWorst]) indicates
314  the vertex of the simplex to be replaced. */
315 
316  ylo = mValue[0];
317  yhi = ylo;
318  iBest = 0;
319  iWorst = 0;
320 
321  for (i = 1; i < np1; ++i)
322  {
323  if (mValue[i] < ylo) /* -- find the lowest value of the objective */
324  {
325  ylo = mValue[i];
326  iBest = i;
327  }
328 
329  if (mValue[i] > yhi) /* -- find the highest value of the objective -- */
330  {
331  yhi = mValue[i];
332  iWorst = i;
333  }
334  }
335 
336  /* ---- Calculate mCentroid, the centroid of the simplex vertices
337  excepting that with the y value yhi. (largest objective value) */
338 
339  for (i = 0; i < mVariableSize; ++i)
340  {
341  /* -- calculate the average in each dimension -- */
342  z = 0.0;
343 
344  for (j = 0; j < np1; ++j)
345  z += mSimplex[i][j];
346 
347  z -= mSimplex[i][iWorst];
348  mCentroid[i] = z / mVariableSize;
349  }
350 
351  /* ---- reflection through the centroid. ---- */
352 
353  for (i = 0; i < mVariableSize; ++i)
354  {
355  mCurrent[i] = (1.0 + rcoeff) * mCentroid[i] - rcoeff * mSimplex[i][iWorst];
356  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
357 
358  // enforce_bounds(&pstar[i], i); /* make sure it is inside */
359  }
360 
361  evaluate();
362  reflok = (mEvaluationValue < ylo); /* -- 1 if we have moved downhill -- */
363 
364  if (reflok) /* then */
365  {
366  ylo = mEvaluationValue;
367 
368  for (i = 0; i < mVariableSize; ++i)
369  mSimplex[i][iWorst] = mCurrent[i];
370 
371  mValue[iWorst] = ylo;
372 
373  // and store that value
375  {
378 
379  // We found a new best value lets report it.
381  }
382 
383  /* ---- successful reflection, so be bold and try extension */
384  for (i = 0; i < mVariableSize; ++i)
385  {
386  mCurrent[i] = ecoeff * mCurrent[i] + (1.0 - ecoeff) * mCentroid[i];
387  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
388 
389  // enforce_bounds(&p2star[i], i); /* make sure it is inside */
390  }
391 
392  evaluate();
393  extnok = (mEvaluationValue < ylo); /* -- 1 if we continue moved downhill -- */
394 
395  if (extnok)
396  {
397  /* ---- retain extension ---- */
398  ylo = mEvaluationValue;
399 
400  for (i = 0; i < mVariableSize; ++i)
401  mSimplex[i][iWorst] = mCurrent[i];
402 
403  mValue[iWorst] = ylo;
404 
405  // and store that value
407  {
410 
411  // We found a new best value lets report it.
413  }
414  }
415  } /* successful reflection */
416  else /* if reflok then else */
417  {
418  /* ---- reflection has not been successful. ystar is the function
419  value at the extended point. now look at the other points
420  of the simplex ...
421  if there are none > ystar then contract to a point within the simplex
422  if there is one > ystar then reduce the size of the extension
423  if there are many > ystar then retain the reflection as is. */
424 
425  size_t L = 0;
426 
427  for (i = 0; i < np1; ++i)
428  if (mValue[i] > mEvaluationValue)
429  L++;
430 
431  if (L == 1) /* then */
432  {
433  /* ---- reflection was not successful but we can try reducing the
434  extension by contracting on the reflection side of the centroid.
435  Record the reflection for the contraction below. */
436  for (i = 0; i < mVariableSize; ++i)
437  mSimplex[i][iWorst] = mCurrent[i];
438 
439  mValue[iWorst] = mEvaluationValue;
440  }
441 
442  if (L == 0 || L == 1) /* then */
443 
444  {
445  /* L == 0 : there are no other points with higher objectives so
446  try a contraction on the y(iWorst) side of the centroid.
447  i.e. within the simplex.
448  L == 1 : there is one other point with a higher objective so
449  try a contraction on the reflection side of the centroid. */
450  for (i = 0; i < mVariableSize; ++i)
451  {
452  mCurrent[i] = ccoeff * mSimplex[i][iWorst] + (1.0 - ccoeff) * mCentroid[i];
453  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
454 
455  // may not need to check boundaries since it is contraction
456  // enforce_bounds(&p2star[i], i); /* make sure it is inside */
457  }
458 
459  evaluate();
460  contok = (mEvaluationValue <= mValue[iWorst]); /* -- 1 if we have not gone uphill -- */
461 
462  if (contok) /* then */
463  {
464  /* ---- retain contraction ---- */
465  for (i = 0; i < mVariableSize; ++i)
466  mSimplex[i][iWorst] = mCurrent[i];
467 
468  mValue[iWorst] = mEvaluationValue;
469  }
470  else
471  {
472  /* ---- contract whole simplex about a point within itself
473  but close to the current minimum. ---- */
474  for (j = 0; j < np1; ++j)
475  {
476  for (i = 0; i < mVariableSize; ++i)
477  {
478  mSimplex[i][j] = (mSimplex[i][j] + mSimplex[i][iBest]) * 0.5;
479  mCurrent[i] = mSimplex[i][j];
480  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
481 
482  //enforce_bounds(&xmin[i], i); /* make sure it is inside */
483  }
484 
485  mValue[j] = evaluate();
486 
488  {
489  /* ---- retain extension ---- */
490  // and store that value
493 
494  // We found a new best value lets report it.
496  }
497  }
498  } /* if (contok) */
499  } /* contraction */
500  else /* if (L == 0) or (L == 1) then else */
501  {
502  /* ---- Retain the first reflection as there are many points in
503  the simplex that have higher objective values than
504  this new point. A change is better than nothing. ---- */
505  for (i = 0; i < mVariableSize; ++i)
506  mSimplex[i][iWorst] = mCurrent[i];
507 
508  mValue[iWorst] = mEvaluationValue;
509  }
510  } /* if Reflok then else .. */
511  } /* ---- procedure TakeAStep ----*/
512 
513  /* are we over the limit for number of function evaluations ? */
514  // quit = (icount > kcount);
515  // are we above the iteration limit?
516  ++mIteration;
517 
518  // signal another iteration
519  if (mpCallBack)
520  mContinue &= mpCallBack->progressItem(mhIteration);
521 
522  quit = (mIteration >= mIterationLimit);
523 
524  if (!quit) /* then */
525  {
526  /* ---- check to see if minimum reached.
527  calculation of the variance must be done in the highest
528  precision available. ---- */
529 
530  /* mean */
531  sum = 0.0;
532  curmin = 0.0;
533 
534  for (i = 0; i < np1; ++i)
535  {
536  sum += mValue[i];
537  curmin += mValue[i] * mValue[i];
538  }
539 
540  sum /= np1;
541  curmin /= np1;
542 
543  /* ---- curmin is the variance of the n+1 Fnelmin values at the
544  vertices. If we haven't reached the minimum to the
545  required accuracy then take a few more steps. */
546  found = (sqrt(curmin - sum * sum) < mTolerance);
547  }
548  } /* while not found and not quit ... */
549 
550  /* **** bail out if necessary **** */
551  if (quit || !mContinue) goto Finish;
552 
553  /* ---- Check around the currently selected point to see if it is a local minimum. */
554  small = fabs(mBestValue) * reltol + abstol;
556 
557  for (i = 0; i < mVariableSize; ++i)
558  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
559 
560  for (i = 0; i < mVariableSize; ++i)
561  {
562  /* ---- check along each dimension ---- */
563  C_FLOAT64 delta = mStep[i] * 1.0e-03;
564 
565  /* -- check along one direction -- */
566  (*(*mpSetCalculateVariable)[i])(mCurrent[i] + delta);
567  evaluate();
568 
569  if ((mEvaluationValue - mBestValue) < -small)
570  break;
571 
572  /* -- now check the other way -- */
573  (*(*mpSetCalculateVariable)[i])(mCurrent[i] - delta);
574  evaluate();
575 
576  if ((mEvaluationValue - mBestValue) < -small)
577  break;
578 
579  /* -- back to start -- */
580  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
581  }
582 
583  ok = (i == mVariableSize);
584 
585  /* ---- return on finding a local minimum to the desired accuracy. ---- */
586 
587  if (ok) /* then */
588  {
589  goto Finish;
590  }
591 
592  /* ---- Reduce the size of the simplex and restart the procedure. ---- */
593 
594  found = 0; /* -- we did not find a 1 minimum -- */
595  del = std::max(del * factor, 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon());
596 
597  goto First;
598 
599 Finish: /* end of procedure */
600 
601  if (mpCallBack)
603 
604  return true;
605 }
606 
608 {
609  return true;
610 }
611 
613 {
614  // We do not need to check whether the parametric constraints are fulfilled
615  // since the parameters are created within the bounds.
616 
619 
620  // when we leave the either the parameter or functional domain
621  // we penalize the objective value by forcing it to be larger
622  // than the best value recorded so far.
627 
628  return mEvaluationValue;
629 }
630 
632 {
633  cleanup();
634 
635  if (!COptMethod::initialize()) return false;
636 
637  mIterationLimit = * getValue("Iteration Limit").pUINT;
638  mTolerance = * getValue("Tolerance").pUDOUBLE;
639  mScale = * getValue("Scale").pUDOUBLE;
640 
641  mIteration = 0;
642 
643  if (mpCallBack)
644  mhIteration =
645  mpCallBack->addItem("Current Iteration",
646  mIteration,
647  & mIterationLimit);
648 
649  mVariableSize = mpOptItem->size();
650 
654 
657 
658  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
659 
660  mContinue = true;
661 
662  return true;
663 }
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
COptTask * mpParentTask
Definition: COptMethod.h:58
CMatrix< C_FLOAT64 > mSimplex
virtual bool initialize()
Definition: COptMethod.cpp:189
unsigned C_INT32 mIteration
COptMethodNelderMead(const COptMethodNelderMead &src, const CCopasiContainer *pParent=NULL)
unsigned C_INT32 mIterationLimit
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
CVector< C_FLOAT64 > mStep
COptProblem * mpOptProblem
Definition: COptMethod.h:56
virtual void output(const COutputInterface::Activity &activity)
CVector< C_FLOAT64 > mValue
#define C_INT32
Definition: copasi.h:90
virtual bool progressItem(const size_t &handle)
CVector< C_FLOAT64 > mCentroid
virtual bool calculate()
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
const C_FLOAT64 & evaluate()
const std::vector< UpdateMethod * > * mpSetCalculateVariable
Definition: COptMethod.h:65
const Value & getValue() const
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)
virtual bool finishItem(const size_t &handle)
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
virtual bool checkFunctionalConstraints()
const C_FLOAT64 & getStartValue() const
Definition: COptItem.cpp:199
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
#define C_FLOAT64
Definition: copasi.h:92
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
bool addParameter(const CCopasiParameter &parameter)
CVector< C_FLOAT64 > mCurrent
virtual bool checkParametricConstraints()
const C_FLOAT64 & getCalculateValue() const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
#define max(a, b)
Definition: f2c.h:176