COPASI API  4.16.103
Public Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
COptMethodLevenbergMarquardt Class Reference

#include <COptMethodLevenbergMarquardt.h>

Inheritance diagram for COptMethodLevenbergMarquardt:
Inheritance graph
[legend]
Collaboration diagram for COptMethodLevenbergMarquardt:
Collaboration graph
[legend]

Public Member Functions

 COptMethodLevenbergMarquardt (const COptMethodLevenbergMarquardt &src, const CCopasiContainer *pParent=NULL)
 
virtual bool optimise ()
 
virtual ~COptMethodLevenbergMarquardt ()
 
- Public Member Functions inherited from COptMethod
 COptMethod (const COptMethod &src, const CCopasiContainer *pParent=NULL)
 
bool isBounded (void)
 
virtual bool isValidProblem (const CCopasiProblem *pProblem)
 
void setProblem (COptProblem *problem)
 
virtual ~COptMethod ()
 
- Public Member Functions inherited from CCopasiMethod
 CCopasiMethod (const CCopasiMethod &src, const CCopasiContainer *pParent=NULL)
 
const CCopasiMethod::SubTypegetSubType () const
 
const CCopasiTask::TypegetType () const
 
virtual void load (CReadConfig &configBuffer, CReadConfig::Mode mode=CReadConfig::SEARCH)
 
virtual void print (std::ostream *ostream) const
 
virtual void printResult (std::ostream *ostream) const
 
virtual bool setCallBack (CProcessReport *pCallBack)
 
virtual ~CCopasiMethod ()
 
- Public Member Functions inherited from CCopasiParameterGroup
bool addGroup (const std::string &name)
 
bool addParameter (const CCopasiParameter &parameter)
 
bool addParameter (const std::string &name, const CCopasiParameter::Type type)
 
template<class CType >
bool addParameter (const std::string &name, const CCopasiParameter::Type type, const CType &value)
 
void addParameter (CCopasiParameter *pParameter)
 
CCopasiParameterGroupassertGroup (const std::string &name)
 
template<class CType >
CCopasiParameterassertParameter (const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
 
index_iterator beginIndex () const
 
name_iterator beginName () const
 
 CCopasiParameterGroup (const CCopasiParameterGroup &src, const CCopasiContainer *pParent=NULL)
 
 CCopasiParameterGroup (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &objectType="ParameterGroup")
 
void clear ()
 
virtual bool elevateChildren ()
 
index_iterator endIndex () const
 
name_iterator endName () const
 
CCopasiParameterGroupgetGroup (const std::string &name)
 
const CCopasiParameterGroupgetGroup (const std::string &name) const
 
CCopasiParameterGroupgetGroup (const size_t &index)
 
const CCopasiParameterGroupgetGroup (const size_t &index) const
 
size_t getIndex (const std::string &name) const
 
std::string getKey (const std::string &name) const
 
std::string getKey (const size_t &index) const
 
virtual const std::string & getName (const size_t &index) const
 
virtual const CObjectInterfacegetObject (const CCopasiObjectName &cn) const
 
CCopasiParametergetParameter (const std::string &name)
 
const CCopasiParametergetParameter (const std::string &name) const
 
CCopasiParametergetParameter (const size_t &index)
 
const CCopasiParametergetParameter (const size_t &index) const
 
CCopasiParameter::Type getType (const std::string &name) const
 
CCopasiParameter::Type getType (const size_t &index) const
 
std::string getUniqueParameterName (const CCopasiParameter *pParameter) const
 
const CCopasiParameter::ValuegetValue (const std::string &name) const
 
const CCopasiParameter::ValuegetValue (const size_t &index) const
 
CCopasiParameter::ValuegetValue (const std::string &name)
 
CCopasiParameter::ValuegetValue (const size_t &index)
 
CCopasiParameterGroupoperator= (const CCopasiParameterGroup &rhs)
 
bool removeParameter (const std::string &name)
 
bool removeParameter (const size_t &index)
 
template<class CType >
bool setValue (const std::string &name, const CType &value)
 
template<class CType >
bool setValue (const size_t &index, const CType &value)
 
size_t size () const
 
bool swap (const size_t &iFrom, const size_t &iTo)
 
bool swap (index_iterator &from, index_iterator &to)
 
virtual ~CCopasiParameterGroup ()
 
- Public Member Functions inherited from CCopasiParameter
 CCopasiParameter (const CCopasiParameter &src, const CCopasiContainer *pParent=NULL)
 
 CCopasiParameter (const std::string &name, const Type &type, const void *pValue=NULL, const CCopasiContainer *pParent=NULL, const std::string &objectType="Parameter")
 
virtual CCopasiObjectName getCN () const
 
virtual const std::string & getKey () const
 
virtual std::string getObjectDisplayName (bool regular=true, bool richtext=false) const
 
const CCopasiParameter::TypegetType () const
 
const ValuegetValue () const
 
ValuegetValue ()
 
virtual voidgetValuePointer () const
 
CCopasiObjectgetValueReference () const
 
bool isValidValue (const C_FLOAT64 &value) const
 
bool isValidValue (const C_INT32 &value) const
 
bool isValidValue (const unsigned C_INT32 &value) const
 
bool isValidValue (const bool &value) const
 
bool isValidValue (const std::string &value) const
 
bool isValidValue (const CCopasiObjectName &value) const
 
bool isValidValue (const std::vector< CCopasiParameter * > &value) const
 
CCopasiParameteroperator= (const CCopasiParameter &rhs)
 
template<class CType >
bool setValue (const CType &value)
 
bool setValue (const std::vector< CCopasiParameter * > &value)
 
virtual ~CCopasiParameter ()
 
- Public Member Functions inherited from CCopasiContainer
virtual bool add (CCopasiObject *pObject, const bool &adopt=true)
 
 CCopasiContainer (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=CCopasiObject::Container)
 
 CCopasiContainer (const CCopasiContainer &src, const CCopasiContainer *pParent=NULL)
 
virtual std::string getChildObjectUnits (const CCopasiObject *pObject) const
 
virtual const objectMapgetObjects () const
 
virtual std::string getUnits () const
 
virtual const CCopasiObjectgetValueObject () const
 
virtual bool remove (CCopasiObject *pObject)
 
virtual ~CCopasiContainer ()
 
- Public Member Functions inherited from CCopasiObject
void addDirectDependency (const CCopasiObject *pObject)
 
 CCopasiObject (const CCopasiObject &src, const CCopasiContainer *pParent=NULL)
 
void clearDirectDependencies ()
 
void clearRefresh ()
 
bool dependsOn (DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const
 
void getAllDependencies (DataObjectSet &dependencies, const DataObjectSet &context) const
 
virtual const DataObjectSetgetDirectDependencies (const DataObjectSet &context=DataObjectSet()) const
 
CCopasiContainergetObjectAncestor (const std::string &type) const
 
CCopasiDataModelgetObjectDataModel ()
 
const CCopasiDataModelgetObjectDataModel () const
 
const std::string & getObjectName () const
 
CCopasiContainergetObjectParent () const
 
const std::string & getObjectType () const
 
virtual const
CObjectInterface::ObjectSet
getPrerequisites () const
 
virtual RefreshgetRefresh () const
 
UpdateMethodgetUpdateMethod () const
 
bool hasCircularDependencies (DataObjectSet &candidates, DataObjectSet &verified, const DataObjectSet &context) const
 
bool hasUpdateMethod () const
 
bool isArray () const
 
bool isContainer () const
 
bool isDataModel () const
 
bool isMatrix () const
 
bool isNameVector () const
 
bool isNonUniqueName () const
 
virtual bool isPrerequisiteForContext (const CObjectInterface *pObject, const CMath::SimulationContextFlag &context, const CObjectInterface::ObjectSet &changedObjects) const
 
bool isReference () const
 
bool isRoot () const
 
bool isSeparator () const
 
bool isStaticString () const
 
bool isValueBool () const
 
bool isValueDbl () const
 
bool isValueInt () const
 
bool isValueInt64 () const
 
bool isValueString () const
 
bool isVector () const
 
virtual bool mustBeDeleted (const DataObjectSet &deletedObjects) const
 
void removeDirectDependency (const CCopasiObject *pObject)
 
void setDirectDependencies (const DataObjectSet &directDependencies)
 
bool setObjectName (const std::string &name)
 
virtual bool setObjectParent (const CCopasiContainer *pParent)
 
void setObjectValue (const C_FLOAT64 &value)
 
void setObjectValue (const C_INT32 &value)
 
void setObjectValue (const bool &value)
 
template<class CType >
void setRefresh (CType *pType, void(CType::*method)(void))
 
template<class CType >
void setUpdateMethod (CType *pType, void(CType::*method)(const C_FLOAT64 &))
 
template<class CType >
void setUpdateMethod (CType *pType, void(CType::*method)(const C_INT32 &))
 
template<class CType >
void setUpdateMethod (CType *pType, void(CType::*method)(const bool &))
 
virtual ~CCopasiObject ()
 
- Public Member Functions inherited from CObjectInterface
 CObjectInterface ()
 
virtual ~CObjectInterface ()
 

Private Member Functions

virtual bool cleanup ()
 
 COptMethodLevenbergMarquardt (const CCopasiContainer *pParent=NULL)
 
const C_FLOAT64evaluate ()
 
void gradient ()
 
void hessian ()
 
virtual bool initialize ()
 
void initObjects ()
 

Private Attributes

CVector< C_FLOAT64mBest
 
C_FLOAT64 mBestValue
 
bool mContinue
 
CVector< C_FLOAT64mCurrent
 
C_FLOAT64 mEvaluationValue
 
CVector< C_FLOAT64mGradient
 
bool mHaveResiduals
 
CMatrix< C_FLOAT64mHessian
 
CMatrix< C_FLOAT64mHessianLM
 
size_t mhIteration
 
unsigned C_INT32 mIteration
 
unsigned C_INT32 mIterationLimit
 
C_FLOAT64 mModulation
 
CMatrix< C_FLOAT64mResidualJacobianT
 
CVector< C_FLOAT64mStep
 
CVector< C_FLOAT64mTemp
 
C_FLOAT64 mTolerance
 
size_t mVariableSize
 

Friends

COptMethodCOptMethod::createMethod (CCopasiMethod::SubType subType)
 

Additional Inherited Members

- Public Types inherited from CCopasiMethod
enum  SubType {
  unset = 0, RandomSearch, RandomSearchMaster, SimulatedAnnealing,
  CoranaWalk, DifferentialEvolution, ScatterSearch, GeneticAlgorithm,
  EvolutionaryProgram, SteepestDescent, HybridGASA, GeneticAlgorithmSR,
  HookeJeeves, LevenbergMarquardt, NelderMead, SRES,
  Statistics, ParticleSwarm, Praxis, TruncatedNewton,
  Newton, deterministic, LSODAR, directMethod,
  stochastic, tauLeap, adaptiveSA, hybrid,
  hybridLSODA, hybridODE45, DsaLsodar, tssILDM,
  tssILDMModified, tssCSP, mcaMethodReder, scanMethod,
  lyapWolf, sensMethod, EFMAlgorithm, EFMBitPatternTreeAlgorithm,
  EFMBitPatternAlgorithm, Householder, crossSectionMethod, linearNoiseApproximation
}
 
- Public Types inherited from CCopasiParameterGroup
typedef parameterGroup::iterator index_iterator
 
typedef
CCopasiContainer::objectMap::iterator 
name_iterator
 
typedef std::vector
< CCopasiParameter * > 
parameterGroup
 
- Public Types inherited from CCopasiParameter
enum  Type {
  DOUBLE = 0, UDOUBLE, INT, UINT,
  BOOL, GROUP, STRING, CN,
  KEY, FILE, EXPRESSION, INVALID
}
 
- Public Types inherited from CCopasiContainer
typedef std::multimap
< std::string, CCopasiObject * > 
objectMap
 
- Public Types inherited from CCopasiObject
typedef std::set< const
CCopasiObject * > 
DataObjectSet
 
typedef std::vector< Refresh * > DataUpdateSequence
 
- Public Types inherited from CObjectInterface
typedef std::set< const
CObjectInterface * > 
ObjectSet
 
typedef std::vector
< CObjectInterface * > 
UpdateSequence
 
- Static Public Member Functions inherited from COptMethod
static COptMethodcreateMethod (CCopasiMethod::SubType subType=CCopasiMethod::RandomSearch)
 
- Static Public Member Functions inherited from CCopasiObject
static std::vector< Refresh * > buildUpdateSequence (const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
 
static void setRenameHandler (CRenameHandler *rh)
 
- Static Public Attributes inherited from CCopasiMethod
static const std::string SubTypeName []
 
static const char * XMLSubType []
 
- Static Public Attributes inherited from CCopasiParameter
static const std::string TypeName []
 
static const char * XMLType []
 
- Static Public Attributes inherited from CCopasiContainer
static const std::vector
< CCopasiContainer * > 
EmptyList
 
- Protected Types inherited from CCopasiObject
enum  Flag {
  Container = 0x1, Vector = 0x2, Matrix = 0x4, NameVector = 0x8,
  Reference = 0x10, ValueBool = 0x20, ValueInt = 0x40, ValueInt64 = 0x80,
  ValueDbl = 0x100, NonUniqueName = 0x200, StaticString = 0x400, ValueString = 0x800,
  Separator = 0x1000, ModelEntity = 0x2000, Array = 0x4000, DataModel = 0x8000,
  Root = 0x10000, Gui = 0x20000
}
 
- Protected Member Functions inherited from COptMethod
 COptMethod (const CCopasiTask::Type &taskType, const SubType &subType, const CCopasiContainer *pParent=NULL)
 
- Protected Member Functions inherited from CCopasiMethod
 CCopasiMethod (const CCopasiTask::Type &taskType, const SubType &subType, const CCopasiContainer *pParent=NULL)
 
- Protected Member Functions inherited from CCopasiParameterGroup
 CCopasiParameterGroup ()
 
- Protected Member Functions inherited from CCopasiContainer
template<class CType >
CCopasiObjectaddMatrixReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
 
template<class CType >
CCopasiObjectaddObjectReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
 
template<class CType >
CCopasiObjectaddVectorReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
 
void initObjects ()
 
- Protected Member Functions inherited from CCopasiObject
 CCopasiObject ()
 
 CCopasiObject (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=0)
 
- Protected Attributes inherited from COptMethod
const bool mBounds
 
const std::vector< COptItem * > * mpOptContraints
 
const std::vector< COptItem * > * mpOptItem
 
COptProblemmpOptProblem
 
COptTaskmpParentTask
 
const std::vector
< UpdateMethod * > * 
mpSetCalculateVariable
 
- Protected Attributes inherited from CCopasiMethod
CProcessReportmpCallBack
 
- Protected Attributes inherited from CCopasiParameter
std::string mKey
 
CCopasiObjectmpValueReference
 
size_t mSize
 
Value mValue
 
- Protected Attributes inherited from CCopasiContainer
objectMap mObjects
 
- Static Protected Attributes inherited from CCopasiObject
static CRenameHandlersmpRenameHandler = NULL
 

Detailed Description

Definition at line 31 of file COptMethodLevenbergMarquardt.h.

Constructor & Destructor Documentation

COptMethodLevenbergMarquardt::COptMethodLevenbergMarquardt ( const COptMethodLevenbergMarquardt src,
const CCopasiContainer pParent = NULL 
)

Copy Constructor

Parameters
constCOptMethodLevenbergMarquardt & src
constCCopasiContainer * pParent (default: NULL)

Definition at line 65 of file COptMethodLevenbergMarquardt.cpp.

References initObjects().

66  :
67  COptMethod(src, pParent),
71  mIteration(0),
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),
86 {initObjects();}
#define C_INVALID_INDEX
Definition: copasi.h:222
COptMethodLevenbergMarquardt::~COptMethodLevenbergMarquardt ( )
virtual

Destructor

Definition at line 88 of file COptMethodLevenbergMarquardt.cpp.

References cleanup().

COptMethodLevenbergMarquardt::COptMethodLevenbergMarquardt ( const CCopasiContainer pParent = NULL)
private

Default Constructor

Parameters
constCCopasiContainer * pParent (default: NULL)

Definition at line 34 of file COptMethodLevenbergMarquardt.cpp.

References CCopasiParameterGroup::addParameter(), C_FLOAT64, C_INT32, CCopasiParameter::DOUBLE, initObjects(), and CCopasiParameter::UINT.

34  :
36  mIterationLimit(2000),
37  mTolerance(1.e-006),
38  mModulation(1.e-006),
39  mIteration(0),
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),
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 }
#define C_INVALID_INDEX
Definition: copasi.h:222
#define C_INT32
Definition: copasi.h:90
#define C_FLOAT64
Definition: copasi.h:92
bool addParameter(const CCopasiParameter &parameter)

Member Function Documentation

bool COptMethodLevenbergMarquardt::cleanup ( )
privatevirtual

Cleanup arrays and pointers.

Returns
bool success

Reimplemented from COptMethod.

Definition at line 374 of file COptMethodLevenbergMarquardt.cpp.

Referenced by initialize(), and ~COptMethodLevenbergMarquardt().

375 {
376  return true;
377 }
const C_FLOAT64 & COptMethodLevenbergMarquardt::evaluate ( )
private

Evaluate the objective function

Returns
bool continue

Definition at line 379 of file COptMethodLevenbergMarquardt.cpp.

References COptProblem::calculate(), COptProblem::checkFunctionalConstraints(), COptProblem::checkParametricConstraints(), COptProblem::getCalculateValue(), mBestValue, mContinue, mEvaluationValue, and COptMethod::mpOptProblem.

Referenced by gradient(), hessian(), and optimise().

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 }
COptProblem * mpOptProblem
Definition: COptMethod.h:56
virtual bool calculate()
virtual bool checkFunctionalConstraints()
virtual bool checkParametricConstraints()
const C_FLOAT64 & getCalculateValue() const
void COptMethodLevenbergMarquardt::gradient ( )
private

Calculate the gradient of the objective at the current parameters

Definition at line 448 of file COptMethodLevenbergMarquardt.cpp.

References C_FLOAT64, evaluate(), mContinue, mCurrent, mGradient, mModulation, and mVariableSize.

Referenced by hessian().

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 }
#define C_FLOAT64
Definition: copasi.h:92
void COptMethodLevenbergMarquardt::hessian ( )
private

Calculate the Hessian of the objective at the current point

Definition at line 481 of file COptMethodLevenbergMarquardt.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), C_FLOAT64, C_INT, dgemm_(), evaluate(), gradient(), mContinue, mCurrent, mGradient, mHaveResiduals, mHessian, mModulation, COptMethod::mpOptProblem, mResidualJacobianT, mTemp, mVariableSize, CMatrix< CType >::numCols(), and CVectorCore< CType >::size().

Referenced by optimise().

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
COptProblem * mpOptProblem
Definition: COptMethod.h:56
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)
size_t size() const
Definition: CVector.h:100
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
virtual size_t numCols() const
Definition: CMatrix.h:144
virtual CType * array()
Definition: CMatrix.h:337
bool COptMethodLevenbergMarquardt::initialize ( )
privatevirtual

Initialize arrays and pointer.

Returns
bool success

Reimplemented from COptMethod.

Definition at line 398 of file COptMethodLevenbergMarquardt.cpp.

References CProcessReport::addItem(), cleanup(), CFitProblem::getResiduals(), CCopasiParameter::getValue(), COptMethod::initialize(), mBest, mBestValue, mContinue, mCurrent, mGradient, mHaveResiduals, mHessian, mhIteration, mIteration, mIterationLimit, mModulation, CCopasiMethod::mpCallBack, COptMethod::mpOptItem, COptMethod::mpOptProblem, mResidualJacobianT, mStep, mTolerance, mVariableSize, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pUINT, CMatrix< CType >::resize(), CVector< CType >::resize(), CFitProblem::setResidualsRequired(), and CVectorCore< CType >::size().

Referenced by optimise().

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 }
virtual bool initialize()
Definition: COptMethod.cpp:189
bool setResidualsRequired(const bool &required)
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
COptProblem * mpOptProblem
Definition: COptMethod.h:56
const CVector< C_FLOAT64 > & getResiduals() const
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
const Value & getValue() const
unsigned C_INT32 * pUINT
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
size_t size() const
Definition: CVector.h:100
const std::vector< COptItem * > * mpOptItem
Definition: COptMethod.h:70
CProcessReport * mpCallBack
void COptMethodLevenbergMarquardt::initObjects ( )
private

Initialize contained objects.

Definition at line 91 of file COptMethodLevenbergMarquardt.cpp.

References CCopasiContainer::addObjectReference(), mIteration, CCopasiParameterGroup::removeParameter(), and CCopasiObject::ValueInt.

Referenced by COptMethodLevenbergMarquardt().

92 {
94 
95 #ifndef COPASI_DEBUG
96  removeParameter("Modulation");
97 #endif
98 }
bool removeParameter(const std::string &name)
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
bool COptMethodLevenbergMarquardt::optimise ( void  )
virtual

Execute the optimization algorithm calling simulation routine when needed. It is noted that this procedure can give feedback of its progress by the callback function set with SetCallback. @ return success;

Reimplemented from COptMethod.

Definition at line 100 of file COptMethodLevenbergMarquardt.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), C_FLOAT64, C_INT, COptItem::checkConstraint(), dpotrf_(), dpotrs_(), dsytrf_(), COutputInterface::DURING, evaluate(), CProcessReport::finishItem(), COptItem::getLowerBoundValue(), COptItem::getStartValue(), COptItem::getUpperBoundValue(), hessian(), initialize(), LAMBDA_MAX, mBest, mBestValue, mContinue, mCurrent, mEvaluationValue, mGradient, mHessian, mHessianLM, mhIteration, mIteration, mIterationLimit, CCopasiMethod::mpCallBack, COptMethod::mpOptProblem, COptMethod::mpParentTask, COptMethod::mpSetCalculateVariable, mStep, mTolerance, mVariableSize, CCopasiTask::output(), CProcessReport::progressItem(), CVector< CType >::resize(), and COptProblem::setSolution().

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 }
#define C_INT
Definition: copasi.h:115
virtual C_INT32 checkConstraint() const
Definition: COptItem.cpp:401
COptTask * mpParentTask
Definition: COptMethod.h:58
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)
COptProblem * mpOptProblem
Definition: COptMethod.h:56
virtual void output(const COutputInterface::Activity &activity)
virtual bool progressItem(const size_t &handle)
const std::vector< UpdateMethod * > * mpSetCalculateVariable
Definition: COptMethod.h:65
const C_FLOAT64 * getLowerBoundValue() const
Definition: COptItem.h:191
virtual bool setSolution(const C_FLOAT64 &value, const CVector< C_FLOAT64 > &variables)
virtual bool finishItem(const size_t &handle)
const C_FLOAT64 & getStartValue() const
Definition: COptItem.cpp:199
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
const C_FLOAT64 * getUpperBoundValue() const
Definition: COptItem.h:198
#define LAMBDA_MAX
int dpotrs_(char *uplo, integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *info)
CProcessReport * mpCallBack
virtual CType * array()
Definition: CMatrix.h:337

Friends And Related Function Documentation

Member Data Documentation

CVector< C_FLOAT64 > COptMethodLevenbergMarquardt::mBest
private

The last individual

Definition at line 139 of file COptMethodLevenbergMarquardt.h.

Referenced by initialize(), and optimise().

C_FLOAT64 COptMethodLevenbergMarquardt::mBestValue
private

The best value found so far

Definition at line 169 of file COptMethodLevenbergMarquardt.h.

Referenced by evaluate(), initialize(), and optimise().

bool COptMethodLevenbergMarquardt::mContinue
private

Flag indicating whether the computation shall continue

Definition at line 179 of file COptMethodLevenbergMarquardt.h.

Referenced by evaluate(), gradient(), hessian(), initialize(), and optimise().

CVector< C_FLOAT64 > COptMethodLevenbergMarquardt::mCurrent
private

The current solution guess

Definition at line 134 of file COptMethodLevenbergMarquardt.h.

Referenced by gradient(), hessian(), initialize(), and optimise().

C_FLOAT64 COptMethodLevenbergMarquardt::mEvaluationValue
private

The result of a function evaluation

Definition at line 174 of file COptMethodLevenbergMarquardt.h.

Referenced by evaluate(), and optimise().

CVector< C_FLOAT64 > COptMethodLevenbergMarquardt::mGradient
private

The gradient

Definition at line 144 of file COptMethodLevenbergMarquardt.h.

Referenced by gradient(), hessian(), initialize(), and optimise().

bool COptMethodLevenbergMarquardt::mHaveResiduals
private

Indicate whether we have access to the residuals, i.e., it is a fitting problem we are working on.

Definition at line 185 of file COptMethodLevenbergMarquardt.h.

Referenced by hessian(), and initialize().

CMatrix<C_FLOAT64> COptMethodLevenbergMarquardt::mHessian
private

The Hessian matrix

Definition at line 154 of file COptMethodLevenbergMarquardt.h.

Referenced by hessian(), initialize(), and optimise().

CMatrix<C_FLOAT64> COptMethodLevenbergMarquardt::mHessianLM
private

The Hessian matrix modified according to Levenberg-Marquardt

Definition at line 159 of file COptMethodLevenbergMarquardt.h.

Referenced by optimise().

size_t COptMethodLevenbergMarquardt::mhIteration
private

Handle to the process report item "Current Iteration"

Definition at line 124 of file COptMethodLevenbergMarquardt.h.

Referenced by initialize(), and optimise().

unsigned C_INT32 COptMethodLevenbergMarquardt::mIteration
private

The number of iterations

Definition at line 119 of file COptMethodLevenbergMarquardt.h.

Referenced by initialize(), initObjects(), and optimise().

unsigned C_INT32 COptMethodLevenbergMarquardt::mIterationLimit
private

The maximum number of iterations

Definition at line 104 of file COptMethodLevenbergMarquardt.h.

Referenced by initialize(), and optimise().

C_FLOAT64 COptMethodLevenbergMarquardt::mModulation
private

The modulation factor

Definition at line 114 of file COptMethodLevenbergMarquardt.h.

Referenced by gradient(), hessian(), and initialize().

CMatrix< C_FLOAT64 > COptMethodLevenbergMarquardt::mResidualJacobianT
private

The transpose jacobian of the residuals.

Definition at line 190 of file COptMethodLevenbergMarquardt.h.

Referenced by hessian(), and initialize().

CVector< C_FLOAT64 > COptMethodLevenbergMarquardt::mStep
private

The step taken

Definition at line 149 of file COptMethodLevenbergMarquardt.h.

Referenced by initialize(), and optimise().

CVector< C_FLOAT64 > COptMethodLevenbergMarquardt::mTemp
private

Vector for temporary values

Definition at line 164 of file COptMethodLevenbergMarquardt.h.

Referenced by hessian().

C_FLOAT64 COptMethodLevenbergMarquardt::mTolerance
private

The tolerance

Definition at line 109 of file COptMethodLevenbergMarquardt.h.

Referenced by initialize(), and optimise().

size_t COptMethodLevenbergMarquardt::mVariableSize
private

number of parameters

Definition at line 129 of file COptMethodLevenbergMarquardt.h.

Referenced by gradient(), hessian(), initialize(), and optimise().


The documentation for this class was generated from the following files: