#include <CLyapWolfMethod.h>

Inheritance diagram for CLyapWolfMethod:

Collaboration diagram for CLyapWolfMethod:

Classes
struct	Data

Public Member Functions
bool	calculate ()

	CLyapWolfMethod (const CLyapWolfMethod &src, const CCopasiContainer *pParent=NULL)

virtual bool	elevateChildren ()

void	evalF (const C_FLOAT64 t, const C_FLOAT64 y, C_FLOAT64 *ydot)

virtual bool	isValidProblem (const CCopasiProblem *pProblem)

virtual void	start ()

virtual double	step (const double &deltaT)

	~CLyapWolfMethod ()

Public Member Functions inherited from CLyapMethod
	CLyapMethod (const CLyapMethod &src, const CCopasiContainer *pParent=NULL)

void	setProblem (CLyapProblem *problem)

virtual void	start (const CState *initialState)

	~CLyapMethod ()

Public Member Functions inherited from CCopasiMethod
	CCopasiMethod (const CCopasiMethod &src, const CCopasiContainer *pParent=NULL)

const CCopasiMethod::SubType &	getSubType () const

const CCopasiTask::Type &	getType () 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)

CCopasiParameterGroup *	assertGroup (const std::string &name)

template<class CType >
CCopasiParameter *	assertParameter (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 ()

index_iterator	endIndex () const

name_iterator	endName () const

CCopasiParameterGroup *	getGroup (const std::string &name)

const CCopasiParameterGroup *	getGroup (const std::string &name) const

CCopasiParameterGroup *	getGroup (const size_t &index)

const CCopasiParameterGroup *	getGroup (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 CObjectInterface *	getObject (const CCopasiObjectName &cn) const

CCopasiParameter *	getParameter (const std::string &name)

const CCopasiParameter *	getParameter (const std::string &name) const

CCopasiParameter *	getParameter (const size_t &index)

const CCopasiParameter *	getParameter (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::Value &	getValue (const std::string &name) const

const CCopasiParameter::Value &	getValue (const size_t &index) const

CCopasiParameter::Value &	getValue (const std::string &name)

CCopasiParameter::Value &	getValue (const size_t &index)

CCopasiParameterGroup &	operator= (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::Type &	getType () const

const Value &	getValue () const

Value &	getValue ()

virtual void *	getValuePointer () const

CCopasiObject *	getValueReference () 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

CCopasiParameter &	operator= (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 objectMap &	getObjects () const

virtual std::string	getUnits () const

virtual const CCopasiObject *	getValueObject () 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 DataObjectSet &	getDirectDependencies (const DataObjectSet &context=DataObjectSet()) const

CCopasiContainer *	getObjectAncestor (const std::string &type) const

CCopasiDataModel *	getObjectDataModel ()

const CCopasiDataModel *	getObjectDataModel () const

const std::string &	getObjectName () const

CCopasiContainer *	getObjectParent () const

const std::string &	getObjectType () const

virtual const CObjectInterface::ObjectSet &	getPrerequisites () const

virtual Refresh *	getRefresh () const

UpdateMethod *	getUpdateMethod () 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 ()

Static Public Member Functions
static void	EvalF (const C_INT n, const C_FLOAT64 t, const C_FLOAT64 y, C_FLOAT64 ydot)

Static Public Member Functions inherited from CLyapMethod
static CLyapMethod *	createMethod (CCopasiMethod::SubType subType=CCopasiMethod::lyapWolf)

static bool	isValidSubType (const CCopasiMethod::SubType &subType)

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)

Private Member Functions
	CLyapWolfMethod (const CCopasiContainer *pParent=NULL)

void	initializeParameter ()

void	orthonormalize ()

Static Private Member Functions
static void	add (C_FLOAT64 dbl1, const C_FLOAT64 dbl1End, const C_FLOAT64 &f, const C_FLOAT64 *dbl2)

static C_FLOAT64	norm (const C_FLOAT64 dbl1, const C_FLOAT64 dbl2)

static C_FLOAT64	product (const C_FLOAT64 dbl1, const C_FLOAT64 dbl1End, const C_FLOAT64 *dbl2)

static void	scalarmult (C_FLOAT64 dbl1, const C_FLOAT64 dbl2, const C_FLOAT64 &f)

Private Attributes
CVector< C_FLOAT64 >	mAtol

Data	mData

bool	mDoDivergence

CVector< C_FLOAT64 >	mDWork

C_FLOAT64	mEndt

std::ostringstream	mErrorMsg

CVector< C_INT >	mIWork

CMatrix< C_FLOAT64 >	mJacobian

C_INT	mJType

CLSODA	mLSODA

C_INT	mLsodaStatus

CVector< C_FLOAT64 >	mNorms

unsigned C_INT32	mNumExp

CState *	mpState

CLyapTask *	mpTask

bool	mReducedModel

C_FLOAT64	mRtol

C_INT	mState

C_FLOAT64	mSumDivergence

CVector< C_FLOAT64 >	mSumExponents

size_t	mSystemSize

C_FLOAT64	mTime

CVector< C_FLOAT64 >	mVariables

Friends
CLyapMethod *	CLyapMethod::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 Attributes inherited from CLyapMethod
static const CLyapValidSubTypes	ValidSubTypes

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 CLyapMethod
	CLyapMethod (const CCopasiMethod::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 >
CCopasiObject *	addMatrixReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

template<class CType >
CCopasiObject *	addObjectReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

template<class CType >
CCopasiObject *	addVectorReference (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 CLyapMethod
CLyapProblem *	mpProblem

Protected Attributes inherited from CCopasiMethod
CProcessReport *	mpCallBack

Protected Attributes inherited from CCopasiParameter
std::string	mKey

CCopasiObject *	mpValueReference

size_t	mSize

Value	mValue

Protected Attributes inherited from CCopasiContainer
objectMap	mObjects

Static Protected Attributes inherited from CCopasiObject
static CRenameHandler *	smpRenameHandler = NULL

Detailed Description

Definition at line 36 of file CLyapWolfMethod.h.

Constructor & Destructor Documentation

CLyapWolfMethod::CLyapWolfMethod ( const CCopasiContainer * pParent = NULL )

private

Default constructor.

Parameters

const CCopasiContainer * pParent (default: NULL)

Definition at line 29 of file CLyapWolfMethod.cpp.

References CLyapWolfMethod::Data::dim, initializeParameter(), mData, and CLyapWolfMethod::Data::pMethod.

                                                                 :
   CLyapMethod(CCopasiMethod::lyapWolf, pParent),
   mpState(NULL)
 {
   assert((void *) &mData == (void *) &mData.dim);
 
   mData.pMethod = this;
   initializeParameter();
 }

CLyapWolfMethod::CLyapWolfMethod	(	const CLyapWolfMethod &	src,
		const CCopasiContainer *	pParent = `NULL`
	)

Copy constructor.

Parameters

const CLsodaMethod &	src
const	CCopasiContainer * pParent (default: NULL)

Definition at line 39 of file CLyapWolfMethod.cpp.

References CLyapWolfMethod::Data::dim, initializeParameter(), mData, and CLyapWolfMethod::Data::pMethod.

                                                                   :
   CLyapMethod(src, pParent),
   mpState(NULL)
 {
   assert((void *) &mData == (void *) &mData.dim);
 
   mData.pMethod = this;
   initializeParameter();
 }

CLyapWolfMethod::~CLyapWolfMethod ( )

Destructor.

Definition at line 50 of file CLyapWolfMethod.cpp.

References mpState, and pdelete.

 {
   pdelete(mpState);
 }

Member Function Documentation

void CLyapWolfMethod::add	(	C_FLOAT64 *	dbl1,
		const C_FLOAT64 *	dbl1End,
		const C_FLOAT64 &	f,
		const C_FLOAT64 *	dbl2
	)

staticprivate

Definition at line 494 of file CLyapWolfMethod.cpp.

Referenced by orthonormalize().

 {
   //calculate v1 = v1 + f * v2
   for (; dbl1 != dbl1End; ++dbl1, ++dbl2)
     *dbl1 += f * *dbl2;
 }

bool CLyapWolfMethod::calculate ( void )

virtual

handles the complete calculation

Implements CLyapMethod.

Definition at line 326 of file CLyapWolfMethod.cpp.

References CVectorCore< CType >::array(), CState::beginIndependent(), C_FLOAT64, CCopasiProblem::getModel(), CCopasiObject::getObjectParent(), CLyapProblem::getTransientTime(), CCopasiParameter::getValue(), CLyapTask::mAverageDivergence, mDoDivergence, CLyapTask::methodCallback(), CLyapTask::mExponents, CLyapTask::mIntervalDivergence, CLyapTask::mLocalExponents, mLsodaStatus, mNorms, mNumExp, CLyapMethod::mpProblem, mpState, mpTask, mReducedModel, mSumDivergence, mSumExponents, mSystemSize, mTime, mVariables, orthonormalize(), CCopasiParameter::Value::pUDOUBLE, CModel::setState(), CState::setTime(), CVectorCore< CType >::size(), start(), step(), and CModel::updateSimulatedValues().

 {
   //get the pointer to the parent task. It's needed for output
   mpTask = dynamic_cast<CLyapTask*>(getObjectParent());
   assert(mpTask);
 
   //initialize LSODA
   start();
 
   C_FLOAT64 stepSize = *getValue("Orthonormalization Interval").pUDOUBLE;
   C_FLOAT64 transientTime = mpProblem->getTransientTime() + mTime;
   C_FLOAT64 endTime = mTime + *getValue("Overall time").pUDOUBLE;
   C_FLOAT64 startTime = mTime;
 
   bool flagProceed = true;
   C_FLOAT64 handlerFactor = 100.0 / (endTime - startTime);
 
   //** do the transient **
   C_FLOAT64 CompareTime = transientTime - 100.0 * fabs(transientTime) * std::numeric_limits< C_FLOAT64 >::epsilon();
 
   if (mTime < CompareTime)
     {
       do
         {
           step(transientTime - mTime);
 
           if (mTime > CompareTime) break;
 
           /* Here we will do conditional event processing */
 
           /* Currently this is correct since no events are processed. */
           //CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 12);
 
           flagProceed &= mpTask->methodCallback((mTime - startTime) * handlerFactor, true);
         }
       while (flagProceed);
 
       //mpLyapProblem->getModel()->setState(*mpCurrentState);
       //mpLyapProblem->getModel()->updateSimulatedValues();
     }
   else
     {
     }
 
   //copy working array to state
   memcpy(mpState->beginIndependent(), mVariables.array(), mSystemSize * sizeof(C_FLOAT64));
   mpState->setTime(mTime);
 
   //copy state to model and do output
   mpProblem->getModel()->setState(*mpState);
   mpProblem->getModel()->updateSimulatedValues(mReducedModel);
   mpTask->methodCallback((mTime - startTime) * handlerFactor, false);
   //********
 
   orthonormalize();
 
   if (mDoDivergence)
     *(mVariables.array() + mVariables.size() - 1) = 0; //divergence
 
   mLsodaStatus = 1; //the state has changed, we need to restart lsoda
 
   size_t i;
 
   C_FLOAT64 realStepSize;
 
   do
     {
       realStepSize = step(stepSize);
 
       orthonormalize();
       mLsodaStatus = 1; //the state has changed, we need to restart lsoda
 
       //process results of orthonormalisation
       for (i = 0; i < mNumExp; ++i)
         {
           mpTask->mLocalExponents[i] = log(mNorms[i]);
           mSumExponents[i] += mpTask->mLocalExponents[i];
           mpTask->mLocalExponents[i] = mpTask->mLocalExponents[i] / realStepSize;
           mpTask->mExponents[i] = mSumExponents[i] / (mTime - transientTime);
         }
 
       //process result of divergence integration
       if (mDoDivergence)
         {
           mSumDivergence += *(mVariables.array() + mVariables.size() - 1);
           mpTask->mIntervalDivergence = *(mVariables.array() + mVariables.size() - 1) / realStepSize;
           *(mVariables.array() + mVariables.size() - 1) = 0;
           mpTask->mAverageDivergence = mSumDivergence / (mTime - transientTime);
         }
 
       //copy working array to state
       memcpy(mpState->beginIndependent(), mVariables.array(), mSystemSize * sizeof(C_FLOAT64));
       mpState->setTime(mTime);
       //copy state to model and do output
       mpProblem->getModel()->setState(*mpState);
       mpProblem->getModel()->updateSimulatedValues(mReducedModel);
       flagProceed &= mpTask->methodCallback((mTime - startTime) * handlerFactor, false);
     }
   while ((mTime < endTime) && flagProceed);
 
   return flagProceed;
 }

bool CLyapWolfMethod::elevateChildren ( )

virtual

This methods must be called to elevate subgroups to derived objects. The default implementation does nothing.

Returns: bool success

Reimplemented from CCopasiParameterGroup.

Definition at line 89 of file CLyapWolfMethod.cpp.

References initializeParameter().

 {
   initializeParameter();
   return true;
 }

void CLyapWolfMethod::EvalF	(	const C_INT *	n,
		const C_FLOAT64 *	t,
		const C_FLOAT64 *	y,
		C_FLOAT64 *	ydot
	)

static

Calculate the default absolute tolerance

Parameters

const CModel * pModel

Returns: C_FLOAT64 defaultAtol

Definition at line 239 of file CLyapWolfMethod.cpp.

References evalF(), and CLyapWolfMethod::Data::pMethod.

Referenced by step().

240 {static_cast<Data *>((void *) n)->pMethod->evalF(t, y, ydot);}

void CLyapWolfMethod::evalF	(	const C_FLOAT64 *	t,
		const C_FLOAT64 *	y,
		C_FLOAT64 *	ydot
	)

This evaluates the derivatives for the complete model

Definition at line 242 of file CLyapWolfMethod.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), CState::beginIndependent(), C_FLOAT64, CModel::calculateDerivatives(), CModel::calculateDerivativesX(), CModel::calculateJacobianX(), CCopasiProblem::getModel(), mDoDivergence, mJacobian, mNumExp, CLyapMethod::mpProblem, mpState, mReducedModel, mSystemSize, mVariables, CModel::setState(), CState::setTime(), and CModel::updateSimulatedValues().

Referenced by EvalF().

 {
   assert(y == mVariables.array());
 
   //set time in model
   mpState->setTime(*t);
 
   //copy working array to state
   memcpy(mpState->beginIndependent(), mVariables.array(), mSystemSize * sizeof(C_FLOAT64));
 
   //copy state to model
   CModel * pModel = mpProblem->getModel();
   pModel->setState(*mpState);
   pModel->updateSimulatedValues(mReducedModel);
 
   //the model
   if (mReducedModel)
     pModel->calculateDerivativesX(ydot);
   else
     pModel->calculateDerivatives(ydot);
 
   //the linearized model
 
   //get jacobian
   pModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
 
   //empty dummy entries... to be removed later
   //C_FLOAT64 *dbl, *dblEnd = ydot + mData.dim;
   //for (dbl = ydot + mSystemSize; dbl != dblEnd; ++dbl)
   //  *dbl = 0;
 
   //   char T = 'N';
   //   C_INT M = mSystemSize;
   //   C_INT N = mSystemSize;
   //   C_INT K = 1;
   //   //C_INT LD = mUnscaledElasticities.numCols();
   //
   //   C_FLOAT64 Alpha = 1.0;
   //   C_FLOAT64 Beta = 0.0;
   //
   //   dgemm_(&T, &T, &K, &N, &M, &Alpha,
   //          const_cast<C_FLOAT64*>(y) + mSystemSize, &K,
   //          mJacobian.array(), &M,
   //          &Beta,
   //          ydot + mSystemSize, &M);
 
   C_FLOAT64 *dbl1;
   const C_FLOAT64 *dbl2, *dbl3, *dbl1end, *dbl3end;
 
   dbl1 = ydot + mSystemSize;
 
   size_t i;
 
   for (i = 1; i <= mNumExp; ++i)
     {
       //dbl1 += mSystemSize;
       dbl1end = dbl1 + mSystemSize;
       dbl2 = mJacobian.array();
 
       for (; dbl1 != dbl1end; ++dbl1)
         {
           *dbl1 = 0.0;
 
           dbl3 = y + i * mSystemSize;
           dbl3end = dbl3 + mSystemSize;
 
           for (; dbl3 != dbl3end; ++dbl3, ++dbl2)
             *dbl1 += *dbl2 * *dbl3;
         }
     }
 
   if (!mDoDivergence) return;
 
   //divergence; trace of jacobian
   *dbl1 = 0;
   dbl2 = mJacobian.array();
 
   for (i = 0; i < mSystemSize; ++i, dbl2 += (mSystemSize + 1))
     * dbl1 += *dbl2;
 
   return;
 }

void CLyapWolfMethod::initializeParameter ( )

private

Intialize the method parameter

Definition at line 55 of file CLyapWolfMethod.cpp.

References CCopasiParameterGroup::assertParameter(), C_FLOAT64, C_INT32, CCopasiParameterGroup::getParameter(), CCopasiParameter::getValue(), CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pUDOUBLE, CCopasiParameterGroup::removeParameter(), CCopasiParameterGroup::setValue(), CCopasiParameter::UDOUBLE, and CCopasiParameter::UINT.

Referenced by CLyapWolfMethod(), and elevateChildren().

 {
   assertParameter("Orthonormalization Interval", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0);
   assertParameter("Overall time", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1000.0);
   assertParameter("Relative Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-6);
   assertParameter("Absolute Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-12);
   assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) 10000);
 
   // Check whether we have an (obsolete) parameter "Use Default Absolute Tolerance"
   CCopasiParameter *pParm;
 
   if ((pParm = getParameter("Use Default Absolute Tolerance")) != NULL)
     {
       C_FLOAT64 NewValue;
 
       if (*pParm->getValue().pBOOL)
         {
           // The default
           NewValue = 1.e-12;
         }
       else
         {
           NewValue = *getValue("Absolute Tolerance").pUDOUBLE;
         }
 
       setValue("Absolute Tolerance", NewValue);
       removeParameter("Use Default Absolute Tolerance");
     }
 
   // These parameters are no longer supported.
   removeParameter("Adams Max Order");
   removeParameter("BDF Max Order");
 }

bool CLyapWolfMethod::isValidProblem ( const CCopasiProblem * pProblem )

virtual

Check if the method is suitable for this problem

Returns: bool suitability of the method

Reimplemented from CLyapMethod.

Definition at line 503 of file CLyapWolfMethod.cpp.

References C_FLOAT64, CCopasiMessage::EXCEPTION, CLyapProblem::getTransientTime(), CCopasiParameter::getValue(), CLyapMethod::isValidProblem(), MCLyap, and CCopasiParameter::Value::pUDOUBLE.

 {
   if (!CLyapMethod::isValidProblem(pProblem)) return false;
 
   const CLyapProblem * pLP = dynamic_cast<const CLyapProblem *>(pProblem);
   assert(pLP);
 
   C_FLOAT64 stepSize = *getValue("Orthonormalization Interval").pUDOUBLE;
   C_FLOAT64 transientTime = pLP->getTransientTime();
   C_FLOAT64 endTime = *getValue("Overall time").pUDOUBLE;
 
   if (transientTime >= endTime)
     {
       //
       CCopasiMessage(CCopasiMessage::EXCEPTION, MCLyap + 4);
       return false;
     }
 
   if (stepSize > (endTime - transientTime))
     {
       //
       CCopasiMessage(CCopasiMessage::EXCEPTION, MCLyap + 5);
       return false;
     }
 
   return true;
 }

C_FLOAT64 CLyapWolfMethod::norm	(	const C_FLOAT64 *	dbl1,
		const C_FLOAT64 *	dbl2
	)

staticprivate

Definition at line 463 of file CLyapWolfMethod.cpp.

References C_FLOAT64.

Referenced by orthonormalize().

 {
   C_FLOAT64 sum = 0;
 
   for (; dbl1 != dbl2; ++dbl1)
     sum += *dbl1 * *dbl1;
 
   return sqrt(sum);
 }

void CLyapWolfMethod::orthonormalize ( )

private

Definition at line 429 of file CLyapWolfMethod.cpp.

References add(), CVectorCore< CType >::array(), C_FLOAT64, mNorms, mNumExp, mSystemSize, mVariables, norm(), product(), and scalarmult().

Referenced by calculate(), and start().

 {
   if (mNumExp < 1) return;
 
   //TODO generalize
   C_FLOAT64 *dbl, *dblEnd;
 
   dbl = mVariables.array() + mSystemSize;
   dblEnd = dbl + mSystemSize;
   mNorms[0] = norm(dbl, dblEnd);
   scalarmult(dbl, dblEnd, 1 / mNorms[0]);
 
   size_t i, j;
 
   for (i = 1; i < mNumExp; ++i)
     {
       dbl += mSystemSize;
       dblEnd = dbl + mSystemSize;
 
       //orthogonalisation
       for (j = 0; j < i; ++j)
         {
           add(dbl, dblEnd,
               -product(dbl, dblEnd, mVariables.array() + (j + 1)*mSystemSize),
               mVariables.array() + (j + 1)*mSystemSize);
         }
 
       //normalisation
       mNorms[i] = norm(dbl, dblEnd);
       scalarmult(dbl, dblEnd, 1 / mNorms[i]);
     }
 }

C_FLOAT64 CLyapWolfMethod::product	(	const C_FLOAT64 *	dbl1,
		const C_FLOAT64 *	dbl1End,
		const C_FLOAT64 *	dbl2
	)

staticprivate

Definition at line 482 of file CLyapWolfMethod.cpp.

References C_FLOAT64.

Referenced by orthonormalize().

 {
   C_FLOAT64 sum = 0;
 
   for (; dbl1 != dbl1End; ++dbl1, ++dbl2)
     sum += *dbl1 * *dbl2;
 
   return sum;
 }

void CLyapWolfMethod::scalarmult	(	C_FLOAT64 *	dbl1,
		const C_FLOAT64 *	dbl2,
		const C_FLOAT64 &	f
	)

staticprivate

Definition at line 474 of file CLyapWolfMethod.cpp.

Referenced by orthonormalize().

 {
   for (; dbl1 != dbl2; ++dbl1)
     *dbl1 *= f;
 }

void CLyapWolfMethod::start ( )

virtual

This instructs the method to prepare for integration starting with the initialState given.

Parameters

const CState * initialState

Definition at line 140 of file CLyapWolfMethod.cpp.

References CVectorCore< CType >::array(), CState::beginIndependent(), C_FLOAT64, C_INT, CLyapWolfMethod::Data::dim, CLyapProblem::divergenceRequested(), CLyapProblem::getExponentNumber(), CCopasiProblem::getModel(), CModel::getNumDependentReactionMetabs(), CState::getNumIndependent(), CModel::getState(), CState::getTime(), CCopasiParameter::getValue(), CModel::initializeAtolVector(), mAtol, CLyapTask::mAverageDivergence, mData, mDoDivergence, mDWork, mErrorMsg, CLyapTask::mExponents, CLyapTask::mIntervalDivergence, mIWork, mJacobian, mJType, CLyapTask::mLocalExponents, mLSODA, mLsodaStatus, mNorms, mNumExp, CLyapMethod::mpProblem, mpState, mpTask, mReducedModel, mRtol, mState, mSumDivergence, mSumExponents, mSystemSize, mTime, mVariables, orthonormalize(), pdelete, CCopasiParameter::Value::pUDOUBLE, CCopasiParameter::Value::pUINT, CMatrix< CType >::resize(), CVector< CType >::resize(), and CInternalSolver::setOstream().

Referenced by calculate().

 {
   /* Reset lsoda */
   mLsodaStatus = 1;
   mState = 1;
   mJType = 2;
   mErrorMsg.str("");
   mLSODA.setOstream(mErrorMsg);
 
   /* Release previous state and make the initialState the current */
   pdelete(mpState);
   mpState = new CState(mpProblem->getModel()->getState());
 
   //mY = mpState->beginIndependent();
   mTime = mpState->getTime();
 
   mReducedModel = true; /* *getValue("Integrate Reduced Model").pBOOL;*/
 
   if (mReducedModel)
     mSystemSize = mpState->getNumIndependent();
   else
     mSystemSize = mpState->getNumIndependent() + mpProblem->getModel()->getNumDependentReactionMetabs();
 
   mNumExp = mpProblem->getExponentNumber();
   mDoDivergence = mpProblem->divergenceRequested();
 
   //calculate the number of variables for lsoda integration
   if (mDoDivergence)
     mData.dim = (C_INT)(mSystemSize * (1 + mNumExp) + 1);
   else
     mData.dim = (C_INT)(mSystemSize * (1 + mNumExp));
 
   //reserve space for exponents. The vectors in the task are resized by the task because they
   //need to have a minimum size defined in the task
   //pTask->mLocalExponents.resize(mNumExp);
   mSumExponents.resize(mNumExp);
   mNorms.resize(mNumExp);
   //pTask->mExponents.resize(mNumExp);
 
   //initialize the vector on which lsoda will work
   mVariables.resize(mData.dim);
   memcpy(mVariables.array(), mpState->beginIndependent(), mSystemSize * sizeof(C_FLOAT64));
 
   //generate base vectors. Just define some arbitrary starting vectors that are not too specific and orthonormalize
   // first fill the array with 0.1
   C_FLOAT64 *dbl, *dblEnd = mVariables.array() + mData.dim;
 
   for (dbl = mVariables.array() + mSystemSize; dbl != dblEnd; ++dbl)
     *dbl = 0.01;
 
   //now add 1.0
   if (mNumExp > 0)
     for (dbl = mVariables.array() + mSystemSize; dbl < dblEnd; dbl += (mSystemSize + 1))
       * dbl = 1.0;
 
   orthonormalize();
 
   //reserve space for jacobian
   mJacobian.resize(mSystemSize, mSystemSize);
 
   size_t i, imax = mNumExp;
 
   for (i = 0; i < imax; ++i)
     {
       mpTask->mLocalExponents[i] = 0;
       mSumExponents[i] = 0;
       mpTask->mExponents[i] = 0;
     }
 
   mpTask->mIntervalDivergence = 0;
   mSumDivergence = 0;
   mpTask->mAverageDivergence = 0;
 
   /* Configure lsoda */
   mRtol = * getValue("Relative Tolerance").pUDOUBLE;
 
   C_FLOAT64 * pTolerance = getValue("Absolute Tolerance").pUDOUBLE;
   CVector< C_FLOAT64 > tmpAtol = mpProblem->getModel()->initializeAtolVector(*pTolerance, mReducedModel);
 
   mAtol.resize(mData.dim);
 
   for (i = 0; i < mSystemSize; ++i)
     mAtol[i] = tmpAtol[i];
 
   for (i = mSystemSize; (int)i < mData.dim; ++i)
     mAtol[i] = 1e-20;
 
   mDWork.resize(22 + mData.dim * std::max<C_INT>(16, mData.dim + 9));
   mDWork[4] = mDWork[5] = mDWork[6] = mDWork[7] = mDWork[8] = mDWork[9] = 0.0;
   mIWork.resize(20 + mData.dim);
   mIWork[4] = mIWork[6] = mIWork[9] = 0;
 
   mIWork[5] = * getValue("Max Internal Steps").pUINT;
   mIWork[7] = 12;
   mIWork[8] = 5;
 
   return;
 }

double CLyapWolfMethod::step ( const double & deltaT )

virtual

This instructs the method to calculate a time step of deltaT starting with the current state, i.e., the result of the previous step. The new state (after deltaT) is expected in the current state. The return value is the actual timestep taken.

Parameters

const double & deltaT

Reimplemented from CLyapMethod.

Definition at line 95 of file CLyapWolfMethod.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, C_INT, CLyapWolfMethod::Data::dim, EvalF(), CCopasiMessage::EXCEPTION, mAtol, MCTrajectoryMethod, mData, mDWork, mErrorMsg, mIWork, mJType, mLSODA, mLsodaStatus, mpState, mRtol, mState, mTime, mVariables, CState::setTime(), and CVectorCore< CType >::size().

Referenced by calculate().

 {
   if (!mData.dim) //just do nothing if there are no variables
     {
       mTime = mTime + deltaT;
       mpState->setTime(mTime);
 
       return deltaT;
     }
 
   C_FLOAT64 startTime = mTime;
   C_FLOAT64 EndTime = mTime + deltaT;
   C_INT one = 1;
   C_INT two = 2;
   C_INT DSize = (C_INT) mDWork.size();
   C_INT ISize = (C_INT) mIWork.size();
 
   mLSODA(&EvalF , //  1. evaluate F
          &mData.dim , //  2. number of variables
          mVariables.array(), //  3. the array of current concentrations
          &mTime , //  4. the current time
          &EndTime , //  5. the final time
          &two , //  6. vector absolute error, scalar relative error
          &mRtol , //  7. relative tolerance array
          mAtol.array() , //  8. absolute tolerance array
          &mState , //  9. output by overshoot & interpolation
          &mLsodaStatus , // 10. the state control variable
          &one , // 11. futher options (one)
          mDWork.array() , // 12. the double work array
          &DSize , // 13. the double work array size
          mIWork.array() , // 14. the int work array
          &ISize , // 15. the int work array size
          NULL , // 16. evaluate J (not given)
          &mJType);           // 17. the type of Jacobian calculate (2)
 
   if (mLsodaStatus == -1) mLsodaStatus = 2;
 
   if ((mLsodaStatus != 1) && (mLsodaStatus != 2) && (mLsodaStatus != -1))
     {
       CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 6, mErrorMsg.str().c_str());
     }
 
   return mTime - startTime;
 }

Friends And Related Function Documentation

CLyapMethod* CLyapMethod::createMethod ( CCopasiMethod::SubType subType )

friend

Member Data Documentation

CVector<C_FLOAT64> CLyapWolfMethod::mAtol

private

Absolute tolerance.

Definition at line 124 of file CLyapWolfMethod.h.

Referenced by start(), and step().

Data CLyapWolfMethod::mData

private

mData.dim is the dimension of the ODE system. mData.pMethod contains CLsodaMethod * this to be used in the static method EvalF

Definition at line 65 of file CLyapWolfMethod.h.

Referenced by CLyapWolfMethod(), start(), and step().

bool CLyapWolfMethod::mDoDivergence

private

Definition at line 77 of file CLyapWolfMethod.h.

Referenced by calculate(), evalF(), and start().

CVector< C_FLOAT64 > CLyapWolfMethod::mDWork

private

Definition at line 133 of file CLyapWolfMethod.h.

Referenced by start(), and step().

C_FLOAT64 CLyapWolfMethod::mEndt

private

Requested end time.

Definition at line 99 of file CLyapWolfMethod.h.

std::ostringstream CLyapWolfMethod::mErrorMsg

private

Definition at line 129 of file CLyapWolfMethod.h.

Referenced by start(), and step().

CVector< C_INT > CLyapWolfMethod::mIWork

private

Definition at line 134 of file CLyapWolfMethod.h.

Referenced by start(), and step().

CMatrix<C_FLOAT64> CLyapWolfMethod::mJacobian

private

Definition at line 79 of file CLyapWolfMethod.h.

Referenced by evalF(), and start().

C_INT CLyapWolfMethod::mJType

private

Definition at line 141 of file CLyapWolfMethod.h.

Referenced by start(), and step().

CLSODA CLyapWolfMethod::mLSODA

private

Definition at line 131 of file CLyapWolfMethod.h.

Referenced by start(), and step().

C_INT CLyapWolfMethod::mLsodaStatus

private

LSODA state.

Definition at line 104 of file CLyapWolfMethod.h.

Referenced by calculate(), start(), and step().

CVector<C_FLOAT64> CLyapWolfMethod::mNorms

private

Definition at line 84 of file CLyapWolfMethod.h.

Referenced by calculate(), orthonormalize(), and start().

unsigned C_INT32 CLyapWolfMethod::mNumExp

private

Number of exponents to calculate

Definition at line 75 of file CLyapWolfMethod.h.

Referenced by calculate(), evalF(), orthonormalize(), and start().

CState* CLyapWolfMethod::mpState

private

A pointer to the current state in complete model view.

Definition at line 59 of file CLyapWolfMethod.h.

Referenced by calculate(), evalF(), start(), step(), and ~CLyapWolfMethod().

CLyapTask* CLyapWolfMethod::mpTask

private

Definition at line 143 of file CLyapWolfMethod.h.

Referenced by calculate(), and start().

bool CLyapWolfMethod::mReducedModel

private

Definition at line 109 of file CLyapWolfMethod.h.

Referenced by calculate(), evalF(), and start().

C_FLOAT64 CLyapWolfMethod::mRtol

private

Relative tolerance.

Definition at line 114 of file CLyapWolfMethod.h.

Referenced by start(), and step().

C_INT CLyapWolfMethod::mState

private

Definition at line 132 of file CLyapWolfMethod.h.

Referenced by start(), and step().

C_FLOAT64 CLyapWolfMethod::mSumDivergence

private

integrated divergence

Definition at line 89 of file CLyapWolfMethod.h.

Referenced by calculate(), and start().

CVector<C_FLOAT64> CLyapWolfMethod::mSumExponents

private

Definition at line 82 of file CLyapWolfMethod.h.

Referenced by calculate(), and start().

size_t CLyapWolfMethod::mSystemSize

private

Number of variables in the model

Definition at line 70 of file CLyapWolfMethod.h.

Referenced by calculate(), evalF(), orthonormalize(), and start().

C_FLOAT64 CLyapWolfMethod::mTime

private

Current time.

Definition at line 94 of file CLyapWolfMethod.h.

Referenced by calculate(), start(), and step().

CVector< C_FLOAT64 > CLyapWolfMethod::mVariables

private

the vector that contains the extended system

Definition at line 139 of file CLyapWolfMethod.h.

Referenced by calculate(), evalF(), orthonormalize(), start(), and step().

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

Classes

Public Member Functions

Static Public Member Functions

Private Member Functions

Static Private Member Functions

Private Attributes

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation