#include <CMathModel.h>

Inheritance diagram for CMathModel:

Collaboration diagram for CMathModel:

Public Member Functions
void	applyInitialValues ()

std::vector< Refresh * >	buildDependendRefreshList (const std::set< const CCopasiObject * > &changedObjects) const

std::vector< Refresh * >	buildRequiredRefreshList (const std::set< const CCopasiObject * > &requiredObjects) const

void	calculateRootDerivatives (CVector< C_FLOAT64 > &rootDerivatives)

	CMathModel (const CCopasiContainer *pParent=NULL)

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

bool	compile (CModel *pModel)

void	evaluateRoots (CVectorCore< C_FLOAT64 > &rootValues, const bool &ignoreDiscrete)

const C_FLOAT64 &	getInitialTime () const

size_t	getNumRoots () const

const CProcessQueue &	getProcessQueue () const

CProcessQueue &	getProcessQueue ()

const C_FLOAT64 &	getProcessQueueExecutionTime () const

const CVector < CMathTrigger::CRootFinder * > &	getRootFinders () const

bool	processQueue (const C_FLOAT64 &time, const bool &equality, CProcessQueue::resolveSimultaneousAssignments pResolveSimultaneousAssignments)

void	processRoots (const C_FLOAT64 &time, const bool &equality, const bool &correct, const CVector< C_INT > &roots)

void	processRoots (const C_FLOAT64 &time, const CVector< C_INT > &roots)

virtual	~CMathModel ()

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

virtual const DataObjectSet &	getDirectDependencies (const DataObjectSet &context=DataObjectSet()) const

virtual const std::string &	getKey () const

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

CCopasiDataModel *	getObjectDataModel ()

const CCopasiDataModel *	getObjectDataModel () const

virtual std::string	getObjectDisplayName (bool regular=true, bool richtext=false) 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

virtual void *	getValuePointer () 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

virtual void	print (std::ostream *ostream) 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
void	calculateRootJacobian (CMatrix< C_FLOAT64 > &jacobian, const CVector< C_FLOAT64 > &rates)

bool	determineInitialRoots (CVector< C_INT > &foundRoots)

Private Attributes
CCopasiVector< CMathEvent >	mEvents

CModel *	mpModel

CProcessQueue	mProcessQueue

CVector< bool >	mRootDiscrete

CVector< CMathEvent * >	mRootIndex2Event

CVector < CMathTrigger::CRootFinder * >	mRootIndex2RootFinder

std::vector< Refresh * >	mRootRefreshes

CVector< C_FLOAT64 * >	mRootValues

Additional Inherited Members
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 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 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 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 CCopasiContainer
objectMap	mObjects

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

Detailed Description

Definition at line 23 of file CMathModel.h.

Constructor & Destructor Documentation

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

Default constructor

Parameters

const CCopasiContainer * pParent (default: NULL)

Definition at line 21 of file CMathModel.cpp.

                                                        :
   CCopasiContainer("MathModel", pParent, "CMathModel"),
   mpModel(NULL),
   mProcessQueue(),
   mEvents("ListOfMathEvents", this),
   mRootValues(),
   mRootDiscrete(),
   mRootRefreshes(),
   mRootIndex2Event(),
   mRootIndex2RootFinder()
 {}

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

Copy constructor

Parameters

const CMathModel &	src
"const	CCopasiContainer * pParent (default: NULL)

Definition at line 33 of file CMathModel.cpp.

References compile(), and mpModel.

                                                          :
   CCopasiContainer(src, pParent),
   mpModel(src.mpModel),
   mProcessQueue(src.mProcessQueue),
   mEvents("ListOfMathEvents", this),
   mRootValues(),
   mRootDiscrete(),
   mRootRefreshes(),
   mRootIndex2Event(),
   mRootIndex2RootFinder()
 {
   // Compile the math model.
   compile(mpModel);
 
   // TODO
   // The problem is that the pointers to events in the process queue point
   // to events in the source model. We can construct a mapping of source event pointers
   // to event pointers and update each process queue entry.
 }

CMathModel::~CMathModel ( )

virtual

Destructor

Definition at line 55 of file CMathModel.cpp.

56 {}

Member Function Documentation

void CMathModel::applyInitialValues ( )

Initialize all values of the math model with their initial values.

Definition at line 372 of file CMathModel.cpp.

References CCopasiVector< T >::begin(), determineInitialRoots(), CCopasiVector< T >::end(), CModel::getInitialTime(), CProcessQueue::initialize(), mEvents, mpModel, mProcessQueue, mRootRefreshes, and processRoots().

Referenced by CModel::applyInitialValues().

 {
   // Clear the process queue.
   mProcessQueue.initialize(this);
 
   // Prepare the roots, i.e., evaluate them
   // Apply all needed refresh calls to calculate the current root values.
   std::vector< Refresh * >::const_iterator itRefresh = mRootRefreshes.begin();
   std::vector< Refresh * >::const_iterator endRefresh = mRootRefreshes.end();
 
   while (itRefresh != endRefresh)
     (**itRefresh++)();
 
   // Now calculate the current root activities
   CCopasiVector< CMathEvent >::const_iterator itMathEvent = mEvents.begin();
   CCopasiVector< CMathEvent >::const_iterator endMathEvent = mEvents.end();
 
   // for each event
   for (; itMathEvent != endMathEvent; ++itMathEvent)
     {
       (*itMathEvent)->getMathTrigger().applyInitialValues();
     }
 
   // We need to schedule events which fire at t > t_0
 
   // The roots which are checked for inequality with a current root value
   // zero are the candidates. We should toggle them if their time derivative
   // is positive.
 
   // The time derivative of the roots can be calculated as:
   // Dr_i/Dt = sum_j (dr_i/dx_j dx_j/dt) + dr_i/dt
   // dx_j/dt are the rates of the state variables
   CVector< C_INT > InitialRoots;
 
   if (determineInitialRoots(InitialRoots))
     {
       processRoots(mpModel->getInitialTime(), false, false, InitialRoots);
     }
 }

std::vector< Refresh * > CMathModel::buildDependendRefreshList ( const std::set< const CCopasiObject * > & changedObjects ) const

Build a list of refresh calls needed to assure that all objects depending on the changed objects are up to date.

Parameters

const std::set< const CCopasiObject * > & changedObjects

Returns: std::vector< Refresh * > refreshList

Definition at line 461 of file CMathModel.cpp.

References CModelEntity::ASSIGNMENT, CCopasiVector< T >::begin(), CCopasiObject::buildUpdateSequence(), CCopasiObject::dependsOn(), CCopasiVector< T >::end(), CModelEntity::FIXED, CStateTemplate::getEntities(), CModel::getMoieties(), CModel::getStateTemplate(), CModelEntity::getValueReference(), mpModel, CModelEntity::ODE, CModelEntity::REACTIONS, CStateTemplate::size(), and CModelEntity::TIME.

Referenced by CMathEvent::compile().

 {
   if (mpModel == NULL)
     return std::vector< Refresh * >();
 
   // We cannot use the default dependencies for a species which concentration was changed
   // by an event as the concentration is the dependent value. This must not be done for
   // species which where not changed by the event, i.e., the solution must be context specific.
   // This implies that the methods getDirectDependencies and its callers must be aware of the species
   // concentrations in the list of the changed objects.
 
   // CCopasiObject::mDependencies is private this means the only place where can be accessed is
   // CCopasiObject itself. We need to replace all read access the mDirectDependencies with
   // getDirectDependencies in CCopasiObject and than make this method context sensitive. This new
   // method will be virtual and must be overloaded for the concentration and particle number in species.
 
   // We build a list of all simulated objects which depend on the changed values.
   std::set< const CCopasiObject * > RequiredObjects;
 
   // The values which need to be updated are all values needed for simulation.
   // These are all objects which are required to calculate the reaction rates and
   // the right hand side of ODEs in other words all rates of entities of type reaction
   // or ODE.
 
   // Due to the context sensitivity, i.e., the changed concentration is no longer dependent
   // on the particle number and compartment volume, we will miss all particle numbers of changed
   // species which do not directly appear in any calculation.
 
   CModelEntity *const* ppEntity = mpModel->getStateTemplate().getEntities();
   CModelEntity *const* ppEndEntity = ppEntity + mpModel->getStateTemplate().size();
 
   const CMetab * pSpecies;
 
   for (; ppEntity != ppEndEntity; ++ppEntity)
     {
       switch ((*ppEntity)->getStatus())
         {
             // First rates of all entities depending on ODEs and species depending on reactions.
           case CModelEntity::ODE:
           case CModelEntity::REACTIONS:
 
             if ((*ppEntity)->getRateReference()->dependsOn(changedObjects, changedObjects))
               {
                 RequiredObjects.insert((*ppEntity)->getRateReference());
               }
 
             // The break statement is intentionally missing since we need to check
             // the particle values of species for the above and for fixed.
           case CModelEntity::FIXED:
             pSpecies = dynamic_cast< const CMetab * >(*ppEntity);
 
             if (pSpecies != NULL &&
                 pSpecies->getValueReference()->dependsOn(changedObjects, changedObjects))
               {
                 RequiredObjects.insert(pSpecies->getValueReference());
               }
 
             break;
 
           case CModelEntity::ASSIGNMENT:
             // Assignments may have to be recalculated but only if an ODE or a
             // reaction flux depends on it.
             break;
 
           case CModelEntity::TIME:
             // The model time will not advance during an event. For discrete time
             // simulation the users should add a global quantity representing the time.
             break;
         }
     }
 
   // Since we have discontinuous change we have to evaluate the total amount of the moieties
   // which participant changed concentration.
 
   CCopasiVector< CMoiety >::const_iterator itMoiety = mpModel->getMoieties().begin();
   CCopasiVector< CMoiety >::const_iterator endMoiety = mpModel->getMoieties().end();
 
   for (; itMoiety != endMoiety; ++itMoiety)
     {
       if ((*itMoiety)->getValueReference()->dependsOn(changedObjects, changedObjects))
         {
           RequiredObjects.insert((*itMoiety)->getValueReference());
         }
     }
 
   std::set< const CCopasiObject * > UpToDate;
   return CCopasiObject::buildUpdateSequence(RequiredObjects, UpToDate, changedObjects);
 }

std::vector< Refresh * > CMathModel::buildRequiredRefreshList ( const std::set< const CCopasiObject * > & requiredObjects ) const

Build a list of refresh calls needed to assure that required objects are up to date.

Parameters

const std::set< const CCopasiObject * > & requiredObjects

Returns: std::vector< Refresh * > refreshList

Definition at line 452 of file CMathModel.cpp.

References CCopasiObject::buildUpdateSequence().

Referenced by compile(), and CMathEvent::compile().

 {
   // We assume that all values determined by ODEs and reactions are current.
   // This means for reduced models moieties must have been applied.
   std::set< const CCopasiObject * > UpToDate;
 
   return CCopasiObject::buildUpdateSequence(requiredObjects, UpToDate);
 }

void CMathModel::calculateRootDerivatives ( CVector< C_FLOAT64 > & rootDerivatives )

Calculate the time derivative of all roots

Parameters

CVector< C_FLOAT64 > & rootDerivatives

Definition at line 417 of file CMathModel.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), C_FLOAT64, C_INT, CModel::calculateDerivatives(), calculateRootJacobian(), dgemm_(), CModel::getNumDependentReactionMetabs(), CStateTemplate::getNumIndependent(), CModel::getStateTemplate(), K, mpModel, mRootValues, CVector< CType >::resize(), CVectorCore< CType >::size(), and CModel::updateSimulatedValues().

Referenced by CModel::calculateRootDerivatives(), and determineInitialRoots().

 {
   size_t NumCols =
     mpModel->getStateTemplate().getNumIndependent() + mpModel->getNumDependentReactionMetabs() + 1;
 
   CVector< C_FLOAT64 > Rates(NumCols);
   C_FLOAT64 * pRate = Rates.array();
   *pRate = 1.0; // for time
 
   mpModel->updateSimulatedValues(false);
   mpModel->calculateDerivatives(pRate + 1);
 
   CMatrix< C_FLOAT64 > Jacobian;
   calculateRootJacobian(Jacobian, Rates);
 
   rootDerivatives.resize(mRootValues.size());
   C_FLOAT64 * pDerivative = rootDerivatives.array();
 
   // Now calculate derivatives of all metabolites determined by reactions
   char T = 'N';
   C_INT M = 1;
   C_INT N = (C_INT) mRootValues.size();
   C_INT K = (C_INT) NumCols;
   C_FLOAT64 Alpha = 1.0;
   C_FLOAT64 Beta = 0.0;
 
   dgemm_(&T, &T, &M, &N, &K, &Alpha, Rates.array(), &M,
          Jacobian.array(), &K, &Beta, pDerivative, &M);
 }

void CMathModel::calculateRootJacobian	(	CMatrix< C_FLOAT64 > &	jacobian,
		const CVector< C_FLOAT64 > &	rates
	)

private

Calculate the Jacobian for the roots.

Parameters

CMatrix<	C_FLOAT64 > & jacobian
const	CVector< C_FLOAT64 > & rates

Definition at line 591 of file CMathModel.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), CState::beginIndependent(), C_FLOAT64, evaluateRoots(), CModel::getNumDependentReactionMetabs(), CState::getNumIndependent(), CModel::getState(), min, mpModel, mRootValues, CMatrix< CType >::resize(), CModel::setState(), CVectorCore< CType >::size(), and CMatrix< CType >::size().

Referenced by calculateRootDerivatives().

 {
   CState State = mpModel->getState();
 
   size_t NumRows = mRootValues.size();
 
   // Partial derivatives with respect to time and all variables determined by ODEs and reactions.
   size_t NumCols =
     State.getNumIndependent() + mpModel->getNumDependentReactionMetabs() + 1;
 
   size_t Col = 0;
 
   jacobian.resize(NumRows, NumCols);
 
   C_FLOAT64 X1 = 0.0;
   C_FLOAT64 X2 = 0.0;
   C_FLOAT64 InvDelta = 0.0;
 
   CVector< C_FLOAT64 > Y1(NumRows);
   CVector< C_FLOAT64 > Y2(NumRows);
 
   C_FLOAT64 * pX = State.beginIndependent() - 1;
   C_FLOAT64 * pXEnd = pX + NumCols;
 
   C_FLOAT64 * pJacobian = jacobian.array();
   C_FLOAT64 * pJacobianEnd = pJacobian + jacobian.size();
 
   const C_FLOAT64 * pRate = rates.array();
 
   for (; pX != pXEnd; ++pX, ++Col, ++pRate)
     {
       C_FLOAT64 Store = *pX;
 
       if (fabs(*pRate) < 1e4 * std::numeric_limits< C_FLOAT64 >::epsilon() * fabs(Store) ||
           fabs(*pRate) < 1e4 * std::numeric_limits< C_FLOAT64 >::min())
         {
           if (fabs(Store) < 100.0 * std::numeric_limits< C_FLOAT64 >::min())
             {
               X1 = 0.0;
 
               if (Store < 0.0)
                 X2 = -200.0 * std::numeric_limits< C_FLOAT64 >::min();
               else
                 X2 = 200.0 * std::numeric_limits< C_FLOAT64 >::min();
 
               InvDelta = X2;
             }
           else
             {
               X1 = 0.999 * Store;
               X2 = 1.001 * Store;
               InvDelta = 500.0 / Store;
             }
         }
       else
         {
           X1 = Store - 0.001 * *pRate;
           X2 = Store + 0.001 * *pRate;
           InvDelta = 500.0 / *pRate;
         }
 
       *pX = X1;
       mpModel->setState(State);
       evaluateRoots(Y1, true);
 
       *pX = X2;
       mpModel->setState(State);
       evaluateRoots(Y2, true);
 
       *pX = Store;
 
       pJacobian = jacobian.array() + Col;
       C_FLOAT64 * pY1 = Y1.array();
       C_FLOAT64 * pY2 = Y2.array();
 
       for (; pJacobian < pJacobianEnd; pJacobian += NumCols, ++pY1, ++pY2)
         * pJacobian = (*pY2 - *pY1) * InvDelta;
     }
 
   // Undo the changes from evaluateRoots.
   mpModel->setState(State);
   evaluateRoots(Y2, true);
 }

bool CMathModel::compile ( CModel * pModel )

Compile the mathematical model based on given the model.

Parameters

CModel * pModel

Returns: bool success

Definition at line 58 of file CMathModel.cpp.

References CCopasiVector< T >::add(), CVectorCore< CType >::array(), CCopasiVector< T >::begin(), buildRequiredRefreshList(), C_FLOAT64, CCopasiVector< T >::clear(), CMathEvent::compile(), CCopasiVector< T >::end(), CStateTemplate::endDependent(), CStateTemplate::getEntities(), CModel::getEvents(), CMathEvent::getMathTrigger(), CMathTrigger::getRootFinders(), CModel::getStateTemplate(), mEvents, mpModel, mRootDiscrete, mRootIndex2Event, mRootIndex2RootFinder, mRootRefreshes, mRootValues, and CVector< CType >::resize().

Referenced by CMathModel(), and CModel::compile().

 {
   // For each CEvent we need to create a CMathEvent
   mEvents.clear();
 
   // We clear the sequence of root refreshes to avoid accessing invalid pointers.
   mRootRefreshes.clear();
 
   // We clear the map to avoid accessing no longer existing events.
   mRootIndex2Event.resize(0);
 
   // We clear the map to avoid accessing no longer existing root finders.
   mRootIndex2RootFinder.resize(0);
 
   mpModel = pModel;
 
   if (mpModel == NULL)
     return false;
 
   // Now we are ready to start the compilations.
   bool success = true;
 
   // The context in which this model is compiled.
   // TODO When the math model is the only thing required for simulation this should be changed
   // to this.
   std::vector< CCopasiContainer * > Context;
   Context.push_back(mpModel);
 
   CCopasiVectorN< CEvent >::const_iterator itEvent = mpModel->getEvents().begin();
   CCopasiVectorN< CEvent >::const_iterator endEvent = mpModel->getEvents().end();
 
   size_t RootFinderCount = 0;
 
   for (; itEvent != endEvent; ++itEvent)
     {
       CMathEvent * pEvent = new CMathEvent();
       mEvents.add(pEvent, true);
       success &= pEvent->compile(*itEvent, Context);
 
       RootFinderCount += pEvent->getMathTrigger().getRootFinders().size();
     }
 
   // We need to create a CVector of pointers to the current root finder values
   mRootValues.resize(RootFinderCount);
   C_FLOAT64 ** ppRootValue = mRootValues.array();
 
   // We need to create a CVector indicating whether the roots only change
   // during discrete events.
   mRootDiscrete.resize(RootFinderCount);
   bool * pRootDiscrete = mRootDiscrete.array();
 
   // We need create a map of root indexes to events.
   mRootIndex2Event.resize(RootFinderCount);
   CMathEvent ** ppEvent = mRootIndex2Event.array();
 
   // We need create a map of root indexes to root finders.
   mRootIndex2RootFinder.resize(RootFinderCount);
   CMathTrigger::CRootFinder ** ppRootFinder = mRootIndex2RootFinder.array();
 
   // We need to list all state variables to find out whether the roots
   // only change during discrete events.
   std::set< const CCopasiObject * > StateVariables;
   CModelEntity *const*ppEntity = mpModel->getStateTemplate().getEntities();
   CModelEntity *const*ppEntityEnd = mpModel->getStateTemplate().endDependent();
 
   for (; ppEntity != ppEntityEnd; ++ppEntity)
     {
       StateVariables.insert((*ppEntity)->getValueReference());
     }
 
   std::set< const CCopasiObject * > RootValuesDependencies;
 
   CCopasiVector< CMathEvent >::const_iterator itMathEvent = mEvents.begin();
   CCopasiVector< CMathEvent >::const_iterator endMathEvent = mEvents.end();
 
   // for each event
   for (; itMathEvent != endMathEvent; ++itMathEvent)
     {
       CCopasiVector< CMathTrigger::CRootFinder >::const_iterator itRootFinder =
         (*itMathEvent)->getMathTrigger().getRootFinders().begin();
       CCopasiVector< CMathTrigger::CRootFinder >::const_iterator endRootFinder =
         (*itMathEvent)->getMathTrigger().getRootFinders().end();
 
       // for each root finder
       for (; itRootFinder != endRootFinder;
            ++itRootFinder, ++ppRootValue, ++pRootDiscrete, ++ppEvent, ++ppRootFinder)
         {
           // Update the vector of pointers to current root values.
           *ppRootValue = (*itRootFinder)->getRootValuePtr();
 
           // Store whether a root changes only during discrete events.
           (*itRootFinder)->determineDiscrete(StateVariables);
           *pRootDiscrete = (*itRootFinder)->isDiscrete();
 
           // Build the mapping from root values indexes to CMathEvents
           *ppEvent = *itMathEvent;
 
           // Build the mapping from root values indexes to root finders
           *ppRootFinder = *itRootFinder;
 
           // The root finder needs to be up to date
           RootValuesDependencies.insert(*itRootFinder);
         }
     }
 
   // We need to compile a refresh sequence which calculates all root values.
   mRootRefreshes = buildRequiredRefreshList(RootValuesDependencies);
 
   return success;
 }

bool CMathModel::determineInitialRoots ( CVector< C_INT > & foundRoots )

private

Determine whether there are any initial roots and mark them.

Parameters

CVector< C_INT > & foundRoots

Returns: bool found

Definition at line 550 of file CMathModel.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, C_INT, calculateRootDerivatives(), mRootIndex2RootFinder, mRootValues, CVector< CType >::resize(), and CVectorCore< CType >::size().

Referenced by applyInitialValues().

 {
   bool Found = false;
 
   CVector< C_FLOAT64 > RootDerivatives;
   calculateRootDerivatives(RootDerivatives);
 
   foundRoots.resize(mRootValues.size());
   C_INT * pFoundRoot = foundRoots.array();
   C_FLOAT64 ** pRootValue = mRootValues.array();
   C_FLOAT64 * pDerivative = RootDerivatives.array();
   CMathTrigger::CRootFinder ** ppRootFinder = mRootIndex2RootFinder.array();
   CMathTrigger::CRootFinder ** ppEndRootFinder = ppRootFinder + mRootIndex2RootFinder.size();
 
   for (; ppRootFinder != ppEndRootFinder;
        ++ppRootFinder, ++pRootValue, ++pDerivative, ++pFoundRoot)
     {
       if (**pRootValue == 0.0)
         {
           if (!(*ppRootFinder)->isEquality() &&
               *pDerivative > 0.0)
             {
               *pFoundRoot = 1;
               Found = true;
 
               continue;
             }
 
           if ((*ppRootFinder)->isEquality() &&
               *pDerivative < 0.0)
             {
               (*ppRootFinder)->toggle(std::numeric_limits< C_FLOAT64 >::quiet_NaN(), true, false);
             }
         }
 
       *pFoundRoot = 0;
     }
 
   return Found;
 }

void CMathModel::evaluateRoots	(	CVectorCore< C_FLOAT64 > &	rootValues,
		const bool &	ignoreDiscrete
	)

Evaluate all root values for the current state of the model. If ignoreDiscrete is true discrete roots evaluate to 1.0.

Parameters

CVectorCore<	double > & rootValues
const	bool & ignoreDiscrete

Definition at line 169 of file CMathModel.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, mRootDiscrete, mRootRefreshes, mRootValues, and CVectorCore< CType >::size().

Referenced by calculateRootJacobian(), CModel::evaluateRoots(), CProcessQueue::process(), and CProcessQueue::rootsFound().

 {
   // Apply all needed refresh calls to calculate the current root values.
   std::vector< Refresh * >::const_iterator itRefresh = mRootRefreshes.begin();
   std::vector< Refresh * >::const_iterator endRefresh = mRootRefreshes.end();
 
   while (itRefresh != endRefresh)
     (**itRefresh++)();
 
   // Copy the current values to the output vector.
   assert(rootValues.size() == mRootValues.size());
 
   C_FLOAT64 *pTarget = rootValues.array();
   C_FLOAT64 **pSrc = mRootValues.array();
   C_FLOAT64 **pSrcEnd = pSrc + mRootValues.size();
 
   if (ignoreDiscrete)
     {
       bool *pDiscrete = mRootDiscrete.array();
 
       for (; pSrc != pSrcEnd; ++pSrc, ++pDiscrete, ++pTarget)
         {
           *pTarget = (*pDiscrete) ? 1.0 : **pSrc;
         }
     }
   else
     {
       for (; pSrc != pSrcEnd; ++pSrc, ++pTarget)
         {
           *pTarget = **pSrc;
         }
     }
 
   return;
 }

const C_FLOAT64 & CMathModel::getInitialTime ( ) const

Retreive the initial time

Returns: const C_FLOAT64 & initialTime

Definition at line 206 of file CMathModel.cpp.

References CModel::getInitialTime(), and mpModel.

Referenced by CProcessQueue::initialize().

 {
   return mpModel->getInitialTime();
 }

size_t CMathModel::getNumRoots ( ) const

Retrieve the number of roots used in checking for discontinuities.

Returns: size_t numRoots

Definition at line 412 of file CMathModel.cpp.

References mRootValues, and CVectorCore< CType >::size().

Referenced by CModel::getNumRoots(), and CProcessQueue::initialize().

 {
   return mRootValues.size();
 }

const CProcessQueue & CMathModel::getProcessQueue ( ) const

Definition at line 362 of file CMathModel.cpp.

References mProcessQueue.

Referenced by CCrossSectionTask::finish(), CTrajectoryTask::initialize(), CCrossSectionTask::process(), and CCrossSectionTask::restore().

 {
   return mProcessQueue;
 }

CProcessQueue & CMathModel::getProcessQueue ( )

Definition at line 367 of file CMathModel.cpp.

References mProcessQueue.

 {
   return mProcessQueue;
 }

const C_FLOAT64 & CMathModel::getProcessQueueExecutionTime ( ) const

Retrieve the next execution time scheduled in the process queue

Returns: const C_FLOAT64 & processQueueExecutionTime

Definition at line 357 of file CMathModel.cpp.

References CProcessQueue::getProcessQueueExecutionTime(), and mProcessQueue.

Referenced by CModel::getProcessQueueExecutionTime().

 {
   return mProcessQueue.getProcessQueueExecutionTime();
 }

const CVector< CMathTrigger::CRootFinder * > & CMathModel::getRootFinders ( ) const

Retrieve a vector of root finders

Returns: const CVector< CMathTrigger::CRootFinder * > & rootFinders

Definition at line 447 of file CMathModel.cpp.

References mRootIndex2RootFinder.

Referenced by CModel::getRootFinders(), and CProcessQueue::rootsFound().

 {
   return mRootIndex2RootFinder;
 }

bool CMathModel::processQueue	(	const C_FLOAT64 &	time,
		const bool &	equality,
		CProcessQueue::resolveSimultaneousAssignments	pResolveSimultaneousAssignments
	)

Process events scheduled at the given which a are checked for equality or not

Parameters

const	C_FLOAT64 & time
const	bool & equality
CProcessQueue::resolveSimultaneousAssignments	pResolveSimultaneousAssignments

Returns: bool stateChanged

Definition at line 211 of file CMathModel.cpp.

References mProcessQueue, and CProcessQueue::process().

Referenced by CModel::processQueue().

 {
   return mProcessQueue.process(time, equality, pResolveSimultaneousAssignments);
 }

void CMathModel::processRoots	(	const C_FLOAT64 &	time,
		const bool &	equality,
		const bool &	correct,
		const CVector< C_INT > &	roots
	)

Check whether the roots which have value 1 lead to firing of events and schedule them if needed. This method is intended to be called from integration methods, which can distinguish between equality and inequality.

Parameters

const	C_FLOAT64 & time
const	bool & equality
const	bool & correct
const	CVector< C_INT > & roots

Definition at line 218 of file CMathModel.cpp.

References CVectorCore< CType >::array(), C_INT, CMathTrigger::calculate(), CMathTrigger::CRootFinder::calculateTrueValue(), CMathEvent::fire(), CMathEvent::getMathTrigger(), mProcessQueue, mRootIndex2Event, mRootIndex2RootFinder, mRootRefreshes, and CVectorCore< CType >::size().

Referenced by applyInitialValues(), CProcessQueue::process(), and CModel::processRoots().

 {
   // Apply all needed refresh calls to calculate the current root values.
   std::vector< Refresh * >::const_iterator itRefresh = mRootRefreshes.begin();
   std::vector< Refresh * >::const_iterator endRefresh = mRootRefreshes.end();
 
   while (itRefresh != endRefresh)
     (**itRefresh++)();
 
   assert(foundRoots.size() == mRootIndex2Event.size());
 
   // All events associated with the found roots need to be evaluated whether they fire.
   // In case one fires the corresponding event needs to be scheduled in the process queue.
 
   const C_INT *pFoundRoot = foundRoots.array();
   const C_INT *pFoundRootEnd = pFoundRoot + foundRoots.size();
 
   CMathEvent ** ppEvent = mRootIndex2Event.array();
   CMathEvent * pProcessEvent = NULL;
 
   CMathTrigger::CRootFinder **ppRootFinder = mRootIndex2RootFinder.array();
 
   if (correct)
     {
       while (pFoundRoot != pFoundRootEnd)
         {
           // We reevaluate the state of the non found roots, which should be save.
           if (*pFoundRoot < 1 && (*ppRootFinder)->isEquality() == equality)
             {
               (*ppRootFinder)->calculateTrueValue();
             }
 
           ++pFoundRoot; ++ppRootFinder;
         }
     }
 
   pFoundRoot = foundRoots.array();
   ppRootFinder = mRootIndex2RootFinder.array();
 
   // We go through the list of roots and process the events
   // which need to be checked whether they fire.
   bool TriggerBefore = true;
 
   while (pFoundRoot != pFoundRootEnd)
     {
       pProcessEvent = *ppEvent;
       TriggerBefore = pProcessEvent->getMathTrigger().calculate();
 
       // Process the events for which we have found a root.
       // A found root is indicated by roots[i] = 1 or 0 otherwise.
       while (pFoundRoot != pFoundRootEnd &&
              *ppEvent == pProcessEvent)
         {
           // We must only toggle the roots which are marked.
           if (*pFoundRoot > 0)
             {
               (*ppRootFinder)->toggle(time, equality, correct);
             }
 
           ++pFoundRoot; ++ppEvent; ++ppRootFinder;
         }
 
       bool TriggerAfter = pProcessEvent->getMathTrigger().calculate();
 
       // Check whether the event fires
       if (TriggerAfter == true &&
           TriggerBefore == false)
         {
           pProcessEvent->fire(time, equality, mProcessQueue);
         }
     }
 
   return;
 }

void CMathModel::processRoots	(	const C_FLOAT64 &	time,
		const CVector< C_INT > &	roots
	)

Check whether the roots which have value 1 lead to firing of events and schedule them if needed. This method is intended to be called during discrete event processing e.g. while processing the event queue.

Parameters

const	C_FLOAT64 & time
const	bool & equality
const	bool & correct
const	CVector< C_INT > & roots

Definition at line 296 of file CMathModel.cpp.

References CVectorCore< CType >::array(), C_INT, CMathTrigger::calculate(), CMathEvent::fire(), CMathEvent::getMathTrigger(), mProcessQueue, mRootIndex2Event, mRootIndex2RootFinder, mRootRefreshes, and CVectorCore< CType >::size().

 {
   // Apply all needed refresh calls to calculate the current root values.
   std::vector< Refresh * >::const_iterator itRefresh = mRootRefreshes.begin();
   std::vector< Refresh * >::const_iterator endRefresh = mRootRefreshes.end();
 
   while (itRefresh != endRefresh)
     (**itRefresh++)();
 
   assert(foundRoots.size() == mRootIndex2Event.size());
 
   // All events associated with the found roots need to be evaluated whether they fire.
   // In case one fires the corresponding event needs to be scheduled in the process queue.
 
   const C_INT *pFoundRoot = foundRoots.array();
   const C_INT *pFoundRootEnd = pFoundRoot + foundRoots.size();
 
   CMathEvent ** ppEvent = mRootIndex2Event.array();
   CMathEvent * pProcessEvent = NULL;
 
   CMathTrigger::CRootFinder **ppRootFinder = mRootIndex2RootFinder.array();
 
   // We go through the list of roots and process the events
   // which need to be checked whether they fire.
   bool TriggerBefore = true;
 
   while (pFoundRoot != pFoundRootEnd)
     {
       pProcessEvent = *ppEvent;
       TriggerBefore = pProcessEvent->getMathTrigger().calculate();
 
       // Process the events for which we have found a root.
       // A found root is indicated by roots[i] = 1 or 0 otherwise.
       while (pFoundRoot != pFoundRootEnd &&
              *ppEvent == pProcessEvent)
         {
           // We must only toggle the roots which are marked.
           if (*pFoundRoot > 0)
             {
               (*ppRootFinder)->toggle(time);
             }
 
           ++pFoundRoot; ++ppEvent; ++ppRootFinder;
         }
 
       bool TriggerAfter = pProcessEvent->getMathTrigger().calculate();
 
       // Check whether the event fires
       if (TriggerAfter == true &&
           TriggerBefore == false)
         {
           // We are dealing with discrete root processing and can not distinguish
           // between equality and inequality. Thus we just pick equality = true
           pProcessEvent->fire(time, true, mProcessQueue);
         }
     }
 
   return;
 }