COPASI API  4.16.103
COptMethodLevenbergMarquardt.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 Pedro Mendes, April 2005, from the original Gepasi file
16 // levmarq.cpp : Restricted step Newton method (Marquardt iteration)
17 
18 #include <cmath>
19 #include "copasi.h"
20 
22 #include "COptProblem.h"
23 #include "COptItem.h"
24 #include "COptTask.h"
25 
28 
29 #include "lapack/lapackwrap.h"
30 #include "lapack/blaswrap.h"
31 
32 #define LAMBDA_MAX 1e80
33 
35  COptMethod(CCopasiTask::optimization, CCopasiMethod::LevenbergMarquardt, pParent),
36  mIterationLimit(2000),
37  mTolerance(1.e-006),
38  mModulation(1.e-006),
39  mIteration(0),
40  mhIteration(C_INVALID_INDEX),
41  mVariableSize(0),
42  mCurrent(),
43  mBest(),
44  mGradient(),
45  mStep(),
46  mHessian(),
47  mHessianLM(),
48  mTemp(),
49  mBestValue(std::numeric_limits< C_FLOAT64 >::infinity()),
50  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::infinity()),
51  mContinue(true),
52  mHaveResiduals(false),
53  mResidualJacobianT()
54 {
55  addParameter("Iteration Limit", CCopasiParameter::UINT, (unsigned C_INT32) 2000);
56  addParameter("Tolerance", CCopasiParameter::DOUBLE, (C_FLOAT64) 1.e-006);
57 
58 #ifdef COPASI_DEBUG
59  addParameter("Modulation", CCopasiParameter::DOUBLE, (C_FLOAT64) 1.e-006);
60 #endif // COPASI_DEBUG
61 
62  initObjects();
63 }
64 
66  const CCopasiContainer * pParent):
67  COptMethod(src, pParent),
68  mIterationLimit(src.mIterationLimit),
69  mTolerance(src.mTolerance),
70  mModulation(src.mModulation),
71  mIteration(0),
72  mhIteration(C_INVALID_INDEX),
73  mVariableSize(0),
74  mCurrent(),
75  mBest(),
76  mGradient(),
77  mStep(),
78  mHessian(),
79  mHessianLM(),
80  mTemp(),
81  mBestValue(std::numeric_limits< C_FLOAT64 >::infinity()),
82  mEvaluationValue(std::numeric_limits< C_FLOAT64 >::infinity()),
83  mContinue(true),
84  mHaveResiduals(false),
85  mResidualJacobianT()
86 {initObjects();}
87 
89 {cleanup();}
90 
92 {
94 
95 #ifndef COPASI_DEBUG
96  removeParameter("Modulation");
97 #endif
98 }
99 
101 {
102  if (!initialize())
103  {
104  if (mpCallBack)
106 
107  return false;
108  }
109 
110  C_INT dim, starts, info, nrhs;
111  C_INT one = 1;
112 
113  size_t i;
114  C_FLOAT64 LM_lambda, nu, convp, convx;
115  bool calc_hess;
116  nrhs = 1;
117 
118  dim = (C_INT) mVariableSize;
119 
120 #ifdef XXXX
121  CVector< C_INT > Pivot(dim);
122  CVector< C_FLOAT64 > Work(1);
123  C_INT LWork = -1;
124 
125  // SUBROUTINE DSYTRF(UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
126  dsytrf_("L", &dim, mHessianLM.array(), &dim, Pivot.array(),
127  Work.array(), &LWork, &info);
128  LWork = Work[0];
129  Work.resize(LWork);
130 #endif // XXXX
131 
132  // initial point is first guess but we have to make sure that we
133  // are within the parameter domain
134  for (i = 0; i < mVariableSize; i++)
135  {
136  const COptItem & OptItem = *(*mpOptItem)[i];
137 
138  switch (OptItem.checkConstraint(OptItem.getStartValue()))
139  {
140  case -1:
141  mCurrent[i] = *OptItem.getLowerBoundValue();
142  break;
143 
144  case 1:
145  mCurrent[i] = *OptItem.getUpperBoundValue();
146  break;
147 
148  case 0:
149  mCurrent[i] = OptItem.getStartValue();
150  break;
151  }
152 
153  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
154  }
155 
156  // keep the current parameter for later
157  mBest = mCurrent;
158 
159  evaluate();
160 
161  if (!isnan(mEvaluationValue))
162  {
163  // and store that value
166 
167  // We found a new best value lets report it.
169  }
170 
171  // Initialize LM_lambda
172  LM_lambda = 1.0;
173  nu = 2.0;
174  calc_hess = true;
175  starts = 1;
176 
177  for (mIteration = 0; (mIteration < mIterationLimit) && (nu != 0.0) && mContinue; mIteration++)
178  {
179  // calculate gradient and Hessian
180  if (calc_hess) hessian();
181 
182  calc_hess = true;
183 
184  // mHessianLM will be used for Cholesky decomposition
186 
187  // add Marquardt lambda to Hessian
188  // put -gradient in h
189  for (i = 0; i < mVariableSize; i++)
190  {
191  mHessianLM[i][i] *= 1.0 + LM_lambda; // Improved
192  // mHessianLM[i][i] += LM_lambda; // Original
193  mStep[i] = - mGradient[i];
194  }
195 
196  // SUBROUTINE DSYTRF(UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
197  // dsytrf_("L", &dim, mHessianLM.array(), &dim, Pivot.array(),
198  // Work.array(), &LWork, &info);
199 
200  // SUBROUTINE DPOTRF(UPLO, N, A, LDA, INFO)
201  // Cholesky factorization
202  char UPLO = 'L';
203  dpotrf_(&UPLO, &dim, mHessianLM.array(), &dim, &info);
204 
205  // if Hessian is positive definite solve Hess * h = -grad
206  if (info == 0)
207  {
208  // SUBROUTINE DPOTRS(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
209  dpotrs_(&UPLO, &dim, &one, mHessianLM.array(), &dim, mStep.array(), &dim, &info);
210 
211  // SUBROUTINE DSYTRS(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO);
212  // dsytrs_("L", &dim, &one, mHessianLM.array(), &dim, Pivot.array(), mStep.array(),
213  // &dim, &info);
214  }
215  else
216  {
217  // We are in a concave region. Thus the current step is an over estimation.
218  // We reduce it by dividing by lambda
219  for (i = 0; i < mVariableSize; i++)
220  {
221  mStep[i] /= LM_lambda;
222  }
223  }
224 
225 //REVIEW:START
226 // This code is different between Gepasi and COPASI
227 // Gepasi goes along the direction until it hits the boundary
228 // COPASI moves in a different direction; this could be a problem
229 
230  // Force the parameters to stay within the defined boundaries.
231  // Forcing the parameters gives better solution than forcing the steps.
232  // It gives same results with Gepasi.
233  for (i = 0; i < mVariableSize; i++)
234  {
235  mCurrent[i] = mBest[i] + mStep[i];
236 
237  const COptItem & OptItem = *(*mpOptItem)[i];
238 
239  switch (OptItem.checkConstraint(mCurrent[i]))
240  {
241  case - 1:
242  mCurrent[i] = *OptItem.getLowerBoundValue() + std::numeric_limits< C_FLOAT64 >::epsilon();
243  break;
244 
245  case 1:
246  mCurrent[i] = *OptItem.getUpperBoundValue() - std::numeric_limits< C_FLOAT64 >::epsilon();
247  break;
248 
249  case 0:
250  break;
251  }
252  }
253 
254 // This is the Gepasi code, which would do the truncation along the search line
255 
256  // Assure that the step will stay within the bounds but is
257  // in its original direction.
258  /* C_FLOAT64 Factor = 1.0;
259  while (true)
260  {
261  convp = 0.0;
262 
263  for (i = 0; i < mVariableSize; i++)
264  {
265  mStep[i] *= Factor;
266  mCurrent[i] = mBest[i] + mStep[i];
267  }
268 
269  Factor = 1.0;
270 
271  for (i = 0; i < mVariableSize; i++)
272  {
273  const COptItem & OptItem = *(*mpOptItem)[i];
274 
275  switch (OptItem.checkConstraint(mCurrent[i]))
276  {
277  case - 1:
278  Factor =
279  std::min(Factor, (*OptItem.getLowerBoundValue() - mBest[i]) / mStep[i]);
280  break;
281 
282  case 1:
283  Factor =
284  std::min(Factor, (*OptItem.getUpperBoundValue() - mBest[i]) / mStep[i]);
285  break;
286 
287  case 0:
288  break;
289  }
290  }
291 
292  if (Factor == 1.0) break;
293  }
294  */
295 //REVIEW:END
296 
297  // calculate the relative change in each parameter
298  for (convp = 0.0, i = 0; i < mVariableSize; i++)
299  {
300  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
301  convp += fabs((mCurrent[i] - mBest[i]) / mBest[i]);
302  }
303 
304  // calculate the objective function value
305  evaluate();
306 
307  // set the convergence check amplitudes
308  // convx has the relative change in objective function value
309  convx = (mBestValue - mEvaluationValue) / mBestValue;
310  // convp has the average relative change in parameter values
311  convp /= mVariableSize;
312 
314  {
315  // keep this value
317 
318  // store the new parameter set
319  mBest = mCurrent;
320 
321  // Inform the problem about the new solution.
323 
324  // We found a new best value lets report it.
326 
327  // decrease LM_lambda
328  LM_lambda /= nu;
329 
330  if ((convp < mTolerance) && (convx < mTolerance))
331  {
332  if (starts < 3)
333  {
334  // let's restart with lambda=1
335  LM_lambda = 1.0;
336  starts++;
337  }
338  else
339  // signal the end
340  nu = 0.0;
341  }
342  }
343  else
344  {
345  // restore the old parameter values
346  mCurrent = mBest;
347 
348  for (i = 0; i < mVariableSize; i++)
349  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
350 
351  // if lambda too high terminate
352  if (LM_lambda > LAMBDA_MAX) nu = 0.0;
353  else
354  {
355  // increase lambda
356  LM_lambda *= nu * 2;
357  // don't recalculate the Hessian
358  calc_hess = false;
359  // correct the number of iterations
360  mIteration--;
361  }
362  }
363 
364  if (mpCallBack)
366  }
367 
368  if (mpCallBack)
370 
371  return true;
372 }
373 
375 {
376  return true;
377 }
378 
380 {
381  // We do not need to check whether the parametric constraints are fulfilled
382  // since the parameters are created within the bounds.
383 
386 
387  // when we leave either the parameter or functional domain
388  // we penalize the objective value by forcing it to be larger
389  // than the best value recorded so far.
394 
395  return mEvaluationValue;
396 }
397 
399 {
400  cleanup();
401 
402  if (!COptMethod::initialize()) return false;
403 
404  mModulation = 0.001;
405  mIterationLimit = * getValue("Iteration Limit").pUINT;
406  mTolerance = * getValue("Tolerance").pDOUBLE;
407 
408 #ifdef COPASI_DEBUG
409  mModulation = * getValue("Modulation").pDOUBLE;
410 #endif // COPASI_DEBUG
411 
412  mIteration = 0;
413 
414  if (mpCallBack)
415  mhIteration =
416  mpCallBack->addItem("Current Iteration",
417  mIteration,
418  & mIterationLimit);
419 
420  mVariableSize = mpOptItem->size();
421 
427 
428  mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
429 
430  mContinue = true;
431 
432  CFitProblem * pFitProblem = dynamic_cast< CFitProblem * >(mpOptProblem);
433 
434  if (pFitProblem != NULL)
435  // if (false)
436  {
437  mHaveResiduals = true;
438  pFitProblem->setResidualsRequired(true);
440  }
441  else
442  mHaveResiduals = false;
443 
444  return true;
445 }
446 
447 // evaluate the value of the gradient by forward finite differences
449 {
450  size_t i;
451 
452  C_FLOAT64 y;
453  C_FLOAT64 x;
454  C_FLOAT64 mod1;
455 
456  mod1 = 1.0 + mModulation;
457 
458  y = evaluate();
459 
460  for (i = 0; i < mVariableSize && mContinue; i++)
461  {
462 //REVIEW:START
463  if ((x = mCurrent[i]) != 0.0)
464  {
465  (*(*mpSetCalculateVariable)[i])(x * mod1);
466  mGradient[i] = (evaluate() - y) / (x * mModulation);
467  }
468 
469  else
470  {
471  (*(*mpSetCalculateVariable)[i])(mModulation);
472  mGradient[i] = (evaluate() - y) / mModulation;
473  }
474 
475 //REVIEW:END
476  (*(*mpSetCalculateVariable)[i])(x);
477  }
478 }
479 
480 //evaluate the Hessian
482 {
483  size_t i, j;
484  C_FLOAT64 mod1;
485 
486  mod1 = 1.0 + mModulation;
487 
488  if (mHaveResiduals)
489  {
490  evaluate();
491 
492  const CVector< C_FLOAT64 > & Residuals =
493  static_cast<CFitProblem *>(mpOptProblem)->getResiduals();
494 
495  const CVector< C_FLOAT64 > CurrentResiduals = Residuals;
496 
497  size_t ResidualSize = Residuals.size();
498 
499  C_FLOAT64 * pJacobianT = mResidualJacobianT.array();
500 
501  const C_FLOAT64 * pCurrentResiduals;
502  const C_FLOAT64 * pEnd = CurrentResiduals.array() + ResidualSize;
503  const C_FLOAT64 * pResiduals;
504 
505  C_FLOAT64 Delta;
506  C_FLOAT64 x;
507 
508  for (i = 0; i < mVariableSize && mContinue; i++)
509  {
510 //REVIEW:START
511  if ((x = mCurrent[i]) != 0.0)
512  {
513  Delta = 1.0 / (x * mModulation);
514  (*(*mpSetCalculateVariable)[i])(x * mod1);
515  }
516 
517  else
518  {
519  Delta = 1.0 / mModulation;
520  (*(*mpSetCalculateVariable)[i])(mModulation);
521 //REVIEW:END
522  }
523 
524  // evaluate another column of the Jacobian
525  evaluate();
526  pCurrentResiduals = CurrentResiduals.array();
527  pResiduals = Residuals.array();
528 
529  for (; pCurrentResiduals != pEnd; pCurrentResiduals++, pResiduals++, pJacobianT++)
530  *pJacobianT = (*pResiduals - *pCurrentResiduals) * Delta;
531 
532  (*(*mpSetCalculateVariable)[i])(x);
533  }
534 
535 #ifdef XXXX
536  // calculate the gradient
537  C_INT m = 1;
538  C_INT n = mGradient.size();
539  C_INT k = mResidualJacobianT.numCols(); /* == CurrentResiduals.size() */
540 
541  char op = 'N';
542 
543  C_FLOAT64 Alpha = 1.0;
544  C_FLOAT64 Beta = 0.0;
545 
546  dgemm_("N", "T", &m, &n, &k, &Alpha,
547  const_cast<C_FLOAT64 *>(CurrentResiduals.array()), &m,
548  mResidualJacobianT.array(), &n, &Beta,
549  const_cast<C_FLOAT64 *>(mGradient.array()), &m);
550 #endif //XXXX
551 
552  C_FLOAT64 * pGradient = mGradient.array();
553  pJacobianT = mResidualJacobianT.array();
554 
555  for (i = 0; i < mVariableSize; i++, pGradient++)
556  {
557  *pGradient = 0.0;
558  pCurrentResiduals = CurrentResiduals.array();
559 
560  for (; pCurrentResiduals != pEnd; pCurrentResiduals++, pJacobianT++)
561  *pGradient += *pJacobianT * *pCurrentResiduals;
562 
563  // This is formally correct but cancels out with factor 2 below
564  // *pGradient *= 2.0;
565  }
566 
567  // calculate the Hessian
568  C_FLOAT64 * pHessian;
569  C_FLOAT64 * pJacobian;
570 
571  for (i = 0; i < mVariableSize; i++)
572  {
573  pHessian = mHessian[i];
574 
575  for (j = 0; j <= i; j++, pHessian++)
576  {
577  *pHessian = 0.0;
578  pJacobianT = mResidualJacobianT[i];
579  pEnd = pJacobianT + ResidualSize;
580  pJacobian = mResidualJacobianT[j];
581 
582  for (; pJacobianT != pEnd; pJacobianT++, pJacobian++)
583  *pHessian += *pJacobianT * *pJacobian;
584 
585  // This is formally correct but cancels out with factor 2 above
586  // *pHessian *= 2.0;
587  }
588  }
589  }
590  else
591  {
592  C_FLOAT64 Delta;
593  C_FLOAT64 x;
594 
595  // calculate the gradient
596  gradient();
597 
598  // and store it
599  mTemp = mGradient;
600 
601  // calculate rows of the Hessian
602  for (i = 0; i < mVariableSize; i++)
603  {
604 //REVIEW:START
605  if ((x = mCurrent[i]) != 0.0)
606  {
607  mCurrent[i] = x * mod1;
608  Delta = 1.0 / (x * mModulation);
609  }
610  else
611  {
612  mCurrent[i] = mModulation;
613  Delta = 1.0 / mModulation;
614 //REVIEW:END
615  }
616 
617  (*(*mpSetCalculateVariable)[i])(mCurrent[i]);
618  gradient();
619 
620  for (j = 0; j <= i; j++)
621  mHessian[i][j] = (mGradient[j] - mTemp[j]) * Delta;
622 
623  // restore the original parameter value
624  mCurrent[i] = x;
625  (*(*mpSetCalculateVariable)[i])(x);
626  }
627 
628  // restore the gradient
629  mGradient = mTemp;
630  }
631 
632  // make the matrix symmetric
633  // not really necessary
634  for (i = 0; i < mVariableSize; i++)
635  for (j = i + 1; j < mVariableSize; j++)
636  mHessian[i][j] = mHessian[j][i];
637 }
#define C_INT
Definition: copasi.h:115
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
COptTask * mpParentTask
Definition: COptMethod.h:58
virtual bool initialize()
Definition: COptMethod.cpp:189
bool setResidualsRequired(const bool &required)
int dpotrf_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info)
int dsytrf_(char *uplo, integer *n, doublereal *a, integer *lda, integer *ipiv, doublereal *work, integer *lwork, integer *info)
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 CVector< C_FLOAT64 > & getResiduals() const
virtual void output(const COutputInterface::Activity &activity)
#define C_INT32
Definition: copasi.h:90
virtual bool progressItem(const size_t &handle)
bool removeParameter(const std::string &name)
virtual bool calculate()
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
const std::vector< UpdateMethod * > * mpSetCalculateVariable
Definition: COptMethod.h:65
const Value & getValue() const
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
unsigned C_INT32 * pUINT
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
virtual bool finishItem(const size_t &handle)
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
virtual bool checkFunctionalConstraints()
size_t size() const
Definition: CVector.h:100
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
CType * array()
Definition: CVector.h:139
COptMethodLevenbergMarquardt(const COptMethodLevenbergMarquardt &src, const CCopasiContainer *pParent=NULL)
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
bool addParameter(const CCopasiParameter &parameter)
#define LAMBDA_MAX
virtual bool checkParametricConstraints()
virtual size_t numCols() const
Definition: CMatrix.h:144
int dpotrs_(char *uplo, integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *info)
const C_FLOAT64 & getCalculateValue() const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CProcessReport * mpCallBack
virtual CType * array()
Definition: CMatrix.h:337