#include <CTrajectoryMethodDsaLsodar.h>

Inheritance diagram for CTrajectoryMethodDsaLsodar:

Collaboration diagram for CTrajectoryMethodDsaLsodar:

Classes
class	CPartition

class	CReactionDependencies

Public Member Functions
void	calculateAmu (const size_t &index)

void	calculatePropensities ()

void	cleanup ()

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

virtual C_FLOAT64	doSingleStep (C_FLOAT64 currentTime, C_FLOAT64 endTime)

virtual bool	elevateChildren ()

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

virtual void	evalR (const C_FLOAT64 t, const C_FLOAT64 y, const C_INT nr, C_FLOAT64 r)

void	fireReaction (const size_t &index)

void	integrateDeterministicPart (const C_FLOAT64 &deltaT)

virtual bool	isValidProblem (const CCopasiProblem *pProblem)

virtual void	start (const CState *initialState)

virtual void	stateChanged ()

virtual Status	step (const double &deltaT)

	~CTrajectoryMethodDsaLsodar ()

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

virtual void	evalJ (const C_FLOAT64 t, const C_FLOAT64 y, const C_INT ml, const C_INT mu, C_FLOAT64 pd, const C_INT nRowPD)

	~CLsodaMethod ()

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

const CVector< C_INT > &	getRoots () const

void	setCurrentState (CState *currentState)

void	setProblem (CTrajectoryProblem *problem)

	~CTrajectoryMethod ()

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 ()

Protected Member Functions
	CTrajectoryMethodDsaLsodar (const CCopasiMethod::SubType &subType=DsaLsodar, const CCopasiContainer *pParent=NULL)

Protected Member Functions inherited from CLsodaMethod
	CLsodaMethod (const CCopasiMethod::SubType &subType=deterministic, const CCopasiContainer *pParent=NULL)

void	destroyRootMask ()

Protected Member Functions inherited from CTrajectoryMethod
	CTrajectoryMethod (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
C_FLOAT64	mA0

CVector< C_FLOAT64 >	mAmu

bool	mDoCorrection

size_t	mFirstReactionSpeciesIndex

bool	mMaxStepsReached

size_t	mNextReactionIndex

C_FLOAT64	mNextReactionTime

size_t	mNumReactions

CPartition	mPartition

C_FLOAT64 *	mpLowerLimit

unsigned C_INT32 *	mpMaxSteps

unsigned C_INT32 *	mpPartitioningInterval

C_FLOAT64 *	mpPartitioningSteps

CRandom *	mpRandomGenerator

C_FLOAT64 *	mpUpperLimit

std::vector < CReactionDependencies >	mReactionDependencies

size_t	mStepsAfterPartitionSystem

Protected Attributes inherited from CLsodaMethod
C_INT	mLsodaStatus

CState	mMethodState

bool	mNoODE

CModel *	mpModel

bool *	mpReducedModel

RootMasking	mRootMasking

C_FLOAT64	mTime

Protected Attributes inherited from CTrajectoryMethod
CState *	mpCurrentState

CTrajectoryProblem *	mpProblem

CVector< C_INT >	mRoots

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

Private Member Functions
void	initializeParameter ()

Friends
CTrajectoryMethod *	CTrajectoryMethod::createMethod (CCopasiMethod::SubType subType)

Additional Inherited Members
Public Types inherited from CTrajectoryMethod
enum	Status { FAILURE = -1, NORMAL = 0, ROOT = 1 }

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 CLsodaMethod
static void	EvalF (const C_INT n, const C_FLOAT64 t, const C_FLOAT64 y, C_FLOAT64 ydot)

static void	EvalJ (const C_INT n, const C_FLOAT64 t, const C_FLOAT64 y, const C_INT ml, const C_INT mu, C_FLOAT64 pd, const C_INT *nRowPD)

static void	EvalR (const C_INT n, const C_FLOAT64 t, const C_FLOAT64 y, const C_INT nr, C_FLOAT64 *r)

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

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 CLsodaMethod
enum	RootMasking { NONE = 0, ALL, DISCRETE }

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 }

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

Detailed Description

Definition at line 52 of file CTrajectoryMethodDsaLsodar.h.

Constructor & Destructor Documentation

CTrajectoryMethodDsaLsodar::CTrajectoryMethodDsaLsodar	(	const CCopasiMethod::SubType &	subType = `DsaLsodar`,
		const CCopasiContainer *	pParent = `NULL`
	)

protected

Default constructor.

Parameters

const	CCopasiMethod::SubType & subType (default: DsaLsodar)
const	CCopasiContainer * pParent (default: NULL)

Default constructor.

Definition at line 321 of file CTrajectoryMethodDsaLsodar.cpp.

References CRandom::createGenerator(), initializeParameter(), mpRandomGenerator, and CRandom::mt19937.

                                      :
   CLsodaMethod(subType, pParent)
 {
   mpRandomGenerator = CRandom::createGenerator(CRandom::mt19937);
   initializeParameter();
 }

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

Copy constructor

Parameters

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

Definition at line 329 of file CTrajectoryMethodDsaLsodar.cpp.

References CRandom::createGenerator(), initializeParameter(), mpRandomGenerator, and CRandom::mt19937.

                                      :
   CLsodaMethod(src, pParent)
 {
   mpRandomGenerator = CRandom::createGenerator(CRandom::mt19937);
   initializeParameter();
 }

CTrajectoryMethodDsaLsodar::~CTrajectoryMethodDsaLsodar ( )

Destructor.

Definition at line 340 of file CTrajectoryMethodDsaLsodar.cpp.

References cleanup(), and DESTRUCTOR_TRACE.

 {
   cleanup();
   DESTRUCTOR_TRACE;
 }

Member Function Documentation

void CTrajectoryMethodDsaLsodar::calculateAmu ( const size_t & index )

Calculates an amu value for a given reaction.

Parameters

const size_t & index

Definition at line 869 of file CTrajectoryMethodDsaLsodar.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, mAmu, mDoCorrection, CTrajectoryMethodDsaLsodar::CReactionDependencies::mMethodSubstrates, CTrajectoryMethodDsaLsodar::CReactionDependencies::mModelSubstrates, CTrajectoryMethodDsaLsodar::CReactionDependencies::mpParticleFlux, mReactionDependencies, CTrajectoryMethodDsaLsodar::CReactionDependencies::mSubstrateMultiplier, and CVectorCore< CType >::size().

Referenced by calculatePropensities(), and fireReaction().

 {
   CReactionDependencies & Dependencies = mReactionDependencies[index];
   C_FLOAT64 & Amu = mAmu[index];
 
   Amu = *Dependencies.mpParticleFlux;
 
   if (Amu < 0.0)
     {
       // TODO CRITICAL Create a warning message
       Amu = 0.0;
     }
 
   if (!mDoCorrection)
     {
       return;
     }
 
   C_FLOAT64 SubstrateMultiplier = 1.0;
   C_FLOAT64 SubstrateDevisor = 1.0;
   C_FLOAT64 Multiplicity;
   C_FLOAT64 LowerBound;
   C_FLOAT64 Number;
 
   bool ApplyCorrection = false;
 
   C_FLOAT64 * pMultiplier = Dependencies.mSubstrateMultiplier.array();
   C_FLOAT64 * endMultiplier = pMultiplier + Dependencies.mSubstrateMultiplier.size();
   C_FLOAT64 ** ppLocalSubstrate = Dependencies.mMethodSubstrates.array();
   C_FLOAT64 ** ppModelSubstrate = Dependencies.mModelSubstrates.array();
 
   for (; pMultiplier != endMultiplier; ++pMultiplier, ++ppLocalSubstrate, ++ppModelSubstrate)
     {
       Multiplicity = *pMultiplier;
 
       // TODO We should check the error introduced through rounding.
       **ppLocalSubstrate = floor(**ppModelSubstrate + 0.5);
 
       if (Multiplicity > 1.01)
         {
           ApplyCorrection = true;
 
           Number = **ppLocalSubstrate;
 
           LowerBound = Number - Multiplicity;
           SubstrateDevisor *= pow(Number, Multiplicity - 1.0);  //optimization
           Number -= 1.0;
 
           while (Number > LowerBound)
             {
               SubstrateMultiplier *= Number;
               Number -= 1.0;
             }
         }
     }
 
   // at least one substrate particle number is zero
   if (SubstrateMultiplier < 0.5 || SubstrateDevisor < 0.5)
     {
       Amu = 0.0;
       return;
     }
 
   // rate_factor is the rate function divided by substrate_factor.
   // It would be more efficient if this was generated directly, since in effect we
   // are multiplying and then dividing by the same thing (substrate_factor)!
   //C_FLOAT64 rate_factor = mpModel->getReactions()[index]->calculateParticleFlux();
   if (ApplyCorrection)
     {
       Amu *= SubstrateMultiplier / SubstrateDevisor;
     }
 
   return;
 }

void CTrajectoryMethodDsaLsodar::calculatePropensities ( )

Calculate the propensities of all stochastic reactions

Definition at line 944 of file CTrajectoryMethodDsaLsodar.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, calculateAmu(), mA0, mAmu, CLsodaMethod::mMethodState, mNumReactions, mPartition, CLsodaMethod::mpModel, CTrajectoryMethodDsaLsodar::CPartition::mStochasticReactions, CModel::setState(), and CModel::updateSimulatedValues().

Referenced by doSingleStep(), and start().

 {
   // It suffices to recalculate the propensities for stochastic reactions.
   CReactionDependencies **ppStochastic = mPartition.mStochasticReactions.array();
 
   for (size_t Index = 0; Index != this->mNumReactions; ++Index, ++ppStochastic)
     {
       if (*ppStochastic != NULL)
         {
           calculateAmu(Index);
         }
     }
 
   mpModel->setState(mMethodState);
   mpModel->updateSimulatedValues(false); //for assignments
 
   // calculate the total propensity
   C_FLOAT64 * pAmu = mAmu.array();
   C_FLOAT64 * endAmu = pAmu + mNumReactions;
 
   mA0 = 0.0;
 
   // We must only consider the propensities for stochastic reactions
   ppStochastic = mPartition.mStochasticReactions.array();
 
   for (; pAmu != endAmu; ++pAmu, ++ppStochastic)
     {
       if (*ppStochastic != NULL)
         {
           mA0 += *pAmu;
           assert(mA0 >= 0.0);
         }
     }
 
   return;
 }

void CTrajectoryMethodDsaLsodar::cleanup ( )

Cleans up memory, etc.

Definition at line 774 of file CTrajectoryMethodDsaLsodar.cpp.

References CLsodaMethod::mpModel, and mpRandomGenerator.

Referenced by ~CTrajectoryMethodDsaLsodar().

 {
   delete mpRandomGenerator;
   mpRandomGenerator = NULL;
   mpModel = NULL;
   return;
 }

C_FLOAT64 CTrajectoryMethodDsaLsodar::doSingleStep	(	C_FLOAT64	currentTime,
		C_FLOAT64	endTime
	)

virtual

Simulates the system over the next interval of time. The current time and the end time of the current step() are given as arguments.

Parameters

currentTime	A C_FLOAT64 specifying the current time
endTime	A C_FLOAT64 specifying the end time of the step()

Returns: A C_FLOAT giving the new time

Definition at line 441 of file CTrajectoryMethodDsaLsodar.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, C_INVALID_INDEX, calculatePropensities(), fireReaction(), CRandom::getRandomOO(), integrateDeterministicPart(), mA0, mAmu, CTrajectoryMethodDsaLsodar::CPartition::mHasDeterministic, CTrajectoryMethodDsaLsodar::CPartition::mHasStochastic, min, CLsodaMethod::mMethodState, mNextReactionIndex, mNextReactionTime, mNumReactions, mPartition, CTrajectoryMethod::mpCurrentState, CLsodaMethod::mpModel, mpPartitioningInterval, mpPartitioningSteps, mpRandomGenerator, mStepsAfterPartitionSystem, CTrajectoryMethodDsaLsodar::CPartition::mStochasticReactions, CTrajectoryMethodDsaLsodar::CPartition::rePartition(), CState::setTime(), and CModel::setTime().

Referenced by step().

 {
   C_FLOAT64 DeltaT = 0.0;
   bool FireReaction = false;
 
   // if there are stochastic reactions
   if (mPartition.mHasStochastic) // there is at least one stochastic reaction
     {
       if (mNextReactionIndex == C_INVALID_INDEX)
         {
           if (mA0 != 0)
             {
               mNextReactionTime = curTime - log(mpRandomGenerator->getRandomOO()) / mA0;
 
               // We are sure that we have at least 1 reaction
               mNextReactionIndex = 0;
 
               C_FLOAT64 sum = 0.0;
               C_FLOAT64 rand = mpRandomGenerator->getRandomOO() * mA0;
 
               C_FLOAT64 * pAmu = mAmu.array();
               C_FLOAT64 * endAmu = pAmu + mNumReactions;
               CReactionDependencies **ppStochastic = mPartition.mStochasticReactions.array();
 
               // Only consider stochastic reactions
               for (; (sum <= rand) && (pAmu != endAmu); ++pAmu, ++mNextReactionIndex, ++ppStochastic)
                 {
                   if (*ppStochastic != NULL)
                     {
                       sum += *pAmu;
                     }
                 }
 
               assert(mNextReactionIndex > 0);
 
               mNextReactionIndex--;
             }
           else
             {
               mNextReactionTime = std::numeric_limits< C_FLOAT64 >::infinity();
               mNextReactionIndex = C_INVALID_INDEX;
             }
         }
 
       if (mNextReactionTime <= endTime)
         {
           FireReaction = true;
           DeltaT = mNextReactionTime - curTime;
         }
       else
         {
           DeltaT = std::min(*mpPartitioningSteps, endTime - curTime);
         }
     }
   else
     {
       DeltaT = std::min(*mpPartitioningSteps, endTime - curTime);
     }
 
   // if there are deterministic reactions
   if (mPartition.mHasDeterministic) // there is at least one deterministic reaction
     {
       integrateDeterministicPart(DeltaT);
       calculatePropensities();
 
       // Now the model state and mMethodState are identical and the propensities are
       // calculated accordingly.
     }
 
   if (FireReaction)
     {
       mpModel->setTime(mNextReactionTime);
       mMethodState.setTime(mNextReactionTime);
 
       fireReaction(mNextReactionIndex);
 
       mNextReactionTime = std::numeric_limits< C_FLOAT64 >::infinity();
       mNextReactionIndex = C_INVALID_INDEX;
 
       // Now the model state and mMethodState are identical and the propensities are
       // calculated accordingly.
     }
 
   if (++mStepsAfterPartitionSystem >= *mpPartitioningInterval)
     {
       if (mPartition.rePartition(*mpCurrentState))
         {
           calculatePropensities();
         }
 
       mStepsAfterPartitionSystem = 0;
     }
 
   return DeltaT;
 }

bool CTrajectoryMethodDsaLsodar::elevateChildren ( )

virtual

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

Returns: bool success

Reimplemented from CLsodaMethod.

Definition at line 398 of file CTrajectoryMethodDsaLsodar.cpp.

References CLsodaMethod::elevateChildren(), and initializeParameter().

 {
   bool success = CLsodaMethod::elevateChildren();
 
   initializeParameter();
   return success;
 }

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

virtual

This evaluates the derivatives

Reimplemented from CLsodaMethod.

Definition at line 729 of file CTrajectoryMethodDsaLsodar.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, CModel::calculateDerivatives(), CModel::calculateDerivativesX(), mFirstReactionSpeciesIndex, CLsodaMethod::mMethodState, CLsodaMethod::mNoODE, mPartition, CLsodaMethod::mpModel, CLsodaMethod::mpReducedModel, CTrajectoryMethodDsaLsodar::CPartition::mStochasticSpecies, CModel::setState(), CState::setTime(), CVectorCore< CType >::size(), and CModel::updateSimulatedValues().

 {
   // If we have no ODEs add a constant one.
   if (mNoODE)
     {
       *ydot = 1.0;
       return;
     }
 
   mMethodState.setTime(*t);
 
   mpModel->setState(mMethodState);
   mpModel->updateSimulatedValues(*mpReducedModel);
 
   if (*mpReducedModel)
     mpModel->calculateDerivativesX(ydot);
   else
     mpModel->calculateDerivatives(ydot);
 
   // Mask derivatives of stochastic species;
   bool * pStochastic = mPartition.mStochasticSpecies.array();
   bool * pStochasticEnd = pStochastic + mPartition.mStochasticSpecies.size();
   C_FLOAT64 * pYdot = ydot + mFirstReactionSpeciesIndex - 1;
 
   for (; pStochastic != pStochasticEnd; ++pStochastic, ++pYdot)
     {
       if (*pStochastic)
         {
           *pYdot = 0;
         }
     }
 
   return;
 }

void CTrajectoryMethodDsaLsodar::evalR	(	const C_FLOAT64 *	t,
		const C_FLOAT64 *	y,
		const C_INT *	nr,
		C_FLOAT64 *	r
	)

virtual

This evaluates the roots

Reimplemented from CLsodaMethod.

Definition at line 765 of file CTrajectoryMethodDsaLsodar.cpp.

766 {

767 }

void CTrajectoryMethodDsaLsodar::fireReaction ( const size_t & index )

Executes the specified reaction in the system once.

Parameters

const size_t & index

Executes the specified reaction in the system once.

Parameters

index	A size_t specifying the index of the reaction, which will be fired.
time	The current time

Definition at line 803 of file CTrajectoryMethodDsaLsodar.cpp.

References CVectorCore< CType >::array(), C_FLOAT64, calculateAmu(), CLsodaMethod::destroyRootMask(), CState::getTime(), mA0, mAmu, CTrajectoryMethodDsaLsodar::CReactionDependencies::mCalculations, CTrajectoryMethodDsaLsodar::CReactionDependencies::mDependentReactions, CLsodaMethod::mLsodaStatus, CTrajectoryMethodDsaLsodar::CReactionDependencies::mMethodSpecies, CLsodaMethod::mMethodState, CTrajectoryMethodDsaLsodar::CReactionDependencies::mModelSpecies, mNumReactions, mPartition, mReactionDependencies, CTrajectoryMethodDsaLsodar::CReactionDependencies::mSpeciesMultiplier, CTrajectoryMethodDsaLsodar::CPartition::mStochasticReactions, CLsodaMethod::mTime, and CVectorCore< CType >::size().

Referenced by doSingleStep().

 {
   CReactionDependencies & Dependencies = mReactionDependencies[index];
 
   // Update the method internal and model species numbers
   C_FLOAT64 ** ppModelSpecies = Dependencies.mModelSpecies.array();
   C_FLOAT64 ** ppLocalSpecies = Dependencies.mMethodSpecies.array();
   C_FLOAT64 * pMultiplier = Dependencies.mSpeciesMultiplier.array();
   C_FLOAT64 * endMultiplier = pMultiplier + Dependencies.mSpeciesMultiplier.size();
 
   for (; pMultiplier != endMultiplier; ++pMultiplier, ++ppLocalSpecies, ++ppModelSpecies)
     {
       **ppLocalSpecies += *pMultiplier;
       **ppModelSpecies = **ppLocalSpecies;
     }
 
   // Calculate all values which depend on the firing reaction
   std::vector< Refresh * >::const_iterator itCalcualtion =  Dependencies.mCalculations.begin();
   std::vector< Refresh * >::const_iterator endCalcualtion =  Dependencies.mCalculations.end();
 
   while (itCalcualtion != endCalcualtion)
     {
       (**itCalcualtion++)();
     }
 
   // We do not need to update the the method state since the only independent state
   // values are species of type reaction which are all controlled by the method.
 
   // calculate the propensities which depend on the firing reaction
   size_t * pDependentReaction = Dependencies.mDependentReactions.array();
   size_t * endDependentReactions = pDependentReaction + Dependencies.mDependentReactions.size();
 
   // It suffices to recalculate the propensities for stochastic reactions.
   for (; pDependentReaction != endDependentReactions; ++pDependentReaction)
     {
       if (mPartition.mStochasticReactions[*pDependentReaction] != NULL)
         {
           calculateAmu(*pDependentReaction);
         }
     }
 
   // calculate the total propensity
   C_FLOAT64 * pAmu = mAmu.array();
   C_FLOAT64 * endAmu = pAmu + mNumReactions;
 
   mA0 = 0.0;
 
   // We must only consider the propensities for stochastic reactions
   CReactionDependencies ** ppStochastic = mPartition.mStochasticReactions.array();
 
   for (; pAmu != endAmu; ++pAmu, ++ppStochastic)
     {
       if (*ppStochastic != NULL)
         {
           mA0 += *pAmu;
         }
     }
 
   mTime = mMethodState.getTime();
   mLsodaStatus = 1;
 
   destroyRootMask();
 
   return;
 }

void CTrajectoryMethodDsaLsodar::initializeParameter ( )

private

Initialize the method parameter

Definition at line 346 of file CTrajectoryMethodDsaLsodar.cpp.

References CCopasiParameterGroup::assertParameter(), C_FLOAT64, C_INT32, CCopasiParameterGroup::getParameter(), CCopasiParameter::getValue(), mpLowerLimit, mpMaxSteps, mpPartitioningInterval, mpPartitioningSteps, mpUpperLimit, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pUDOUBLE, CCopasiParameter::Value::pUINT, CCopasiParameterGroup::removeParameter(), CCopasiParameterGroup::setValue(), CCopasiParameter::UDOUBLE, and CCopasiParameter::UINT.

Referenced by CTrajectoryMethodDsaLsodar(), and elevateChildren().

 {
   CCopasiParameter *pParm;
 
   mpMaxSteps =
     assertParameter("Max Internal Steps", CCopasiParameter::UINT, (C_INT32) 1000000)->getValue().pUINT;
   mpLowerLimit =
     assertParameter("Lower Limit", CCopasiParameter::UDOUBLE, (C_FLOAT64) 800.0)->getValue().pUDOUBLE;
   mpUpperLimit =
     assertParameter("Upper Limit", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1000.0)->getValue().pUDOUBLE;
   mpPartitioningInterval =
     assertParameter("Partitioning Interval", CCopasiParameter::UINT, (unsigned C_INT32) 1)->getValue().pUINT;
   mpPartitioningSteps =
     assertParameter("Partitioning Stepsize", CCopasiParameter::UDOUBLE, (C_FLOAT64) 0.001)->getValue().pUDOUBLE;
 
   // Check whether we have a method with the old parameter names
   if ((pParm = getParameter("HYBRID.MaxSteps")) != NULL)
     {
       setValue("Max Internal Steps", *pParm->getValue().pUINT);
       removeParameter("HYBRID.MaxSteps");
 
       if ((pParm = getParameter("HYBRID.LowerStochLimit")) != NULL)
         {
           setValue("Lower Limit", *pParm->getValue().pDOUBLE);
           removeParameter("HYBRID.LowerStochLimit");
         }
 
       if ((pParm = getParameter("HYBRID.UpperStochLimit")) != NULL)
         {
           setValue("Upper Limit", *pParm->getValue().pDOUBLE);
           removeParameter("HYBRID.UpperStochLimit");
         }
 
       if ((pParm = getParameter("HYBRID.PartitioningInterval")) != NULL)
         {
           setValue("Partitioning Interval", *pParm->getValue().pUINT);
           removeParameter("HYBRID.PartitioningInterval");
         }
 
       if ((pParm = getParameter("UseRandomSeed")) != NULL)
         {
           removeParameter("UseRandomSeed");
         }
 
       if ((pParm = getParameter("")) != NULL)
         {
           removeParameter("");
         }
     }
 }

void CTrajectoryMethodDsaLsodar::integrateDeterministicPart ( const C_FLOAT64 & deltaT )

Integrates the deterministic reactions of the system over the specified time interval.

Parameters

const C_FLOAT64 & deltaT.

Definition at line 784 of file CTrajectoryMethodDsaLsodar.cpp.

References CLsodaMethod::mMethodState, CLsodaMethod::mpModel, CLsodaMethod::mpReducedModel, CModel::setState(), CLsodaMethod::step(), and CModel::updateSimulatedValues().

Referenced by doSingleStep().

 {
   CLsodaMethod::step(deltaT);
 
   mpModel->setState(mMethodState);
   mpModel->updateSimulatedValues(*mpReducedModel);
 
   return;
 }

bool CTrajectoryMethodDsaLsodar::isValidProblem ( const CCopasiProblem * pProblem )

virtual

Check if the method is suitable for this problem

Returns: bool suitability of the method

Reimplemented from CTrajectoryMethod.

Definition at line 984 of file CTrajectoryMethodDsaLsodar.cpp.

References CCopasiMessage::ERROR, CModel::getCompartments(), CTrajectoryProblem::getDuration(), CModel::getEvents(), CModel::getMetabolites(), CCopasiProblem::getModel(), CModel::getModelValues(), CModel::getNumMetabs(), CModel::getNumModelValues(), CCopasiParameter::getValue(), CTrajectoryMethod::isValidProblem(), MCTrajectoryMethod, mpLowerLimit, mpUpperLimit, CModelEntity::ODE, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pINT, CCopasiVector< T >::size(), and CModel::suitableForStochasticSimulation().

 {
   if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
 
   const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
 
   if (pTP->getDuration() < 0.0)
     {
       //back integration not possible
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 9);
       return false;
     }
 
   //check for rules
   size_t i, imax = pTP->getModel()->getNumModelValues();
 
   for (i = 0; i < imax; ++i)
     {
       if (pTP->getModel()->getModelValues()[i]->getStatus() == CModelEntity::ODE)
         {
           //ode rule found
           CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 18);
           return false;
         }
 
       /*      if (pTP->getModel()->getModelValues()[i]->getStatus()==CModelEntity::ASSIGNMENT)
               if (pTP->getModel()->getModelValues()[i]->isUsed())
                 {
                   //used assignment found
                   CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 19);
                   return false;
                 }*/
     }
 
   imax = pTP->getModel()->getNumMetabs();
 
   for (i = 0; i < imax; ++i)
     {
       if (pTP->getModel()->getMetabolites()[i]->getStatus() == CModelEntity::ODE)
         {
           //ode rule found
           CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 20);
           return false;
         }
     }
 
   imax = pTP->getModel()->getCompartments().size();
 
   for (i = 0; i < imax; ++i)
     {
       if (pTP->getModel()->getCompartments()[i]->getStatus() == CModelEntity::ODE)
         {
           //ode rule found
           CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 21);
           return false;
         }
     }
 
   //TODO: rewrite CModel::suitableForStochasticSimulation() to use
   //      CCopasiMessage
   std::string message = pTP->getModel()->suitableForStochasticSimulation();
 
   if (message != "")
     {
       //model not suitable, message describes the problem
       CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
       return false;
     }
 
   /* Max Internal Steps */
   if (* getValue("Max Internal Steps").pINT <= 0)
     {
       //max steps should be at least 1
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 15);
       return false;
     }
 
   /* Lower Limit, Upper Limit */
   *mpLowerLimit = * getValue("Lower Limit").pDOUBLE;
   *mpUpperLimit = * getValue("Upper Limit").pDOUBLE;
 
   if (*mpLowerLimit > *mpUpperLimit)
     {
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 4, *mpLowerLimit, *mpUpperLimit);
       return false;
     }
 
   /* Partitioning Interval */
   // nothing to be done here so far
 
   /* Use Random Seed */
   // should be checked in the widget later on
 
   /* Random Seed */
   // nothing to be done here
 
   //events are not supported at the moment
   if (pTP->getModel()->getEvents().size() > 0)
     {
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 23);
       return false;
     }
 
   return true;
 }

void CTrajectoryMethodDsaLsodar::start ( const CState * initialState )

virtual

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

Parameters

const CState * initialState

Reimplemented from CLsodaMethod.

Definition at line 538 of file CTrajectoryMethodDsaLsodar.cpp.

References CCopasiVector< T >::begin(), CStateTemplate::beginIndependent(), CState::beginIndependent(), CCopasiObject::buildUpdateSequence(), C_FLOAT64, C_INVALID_INDEX, calculatePropensities(), CModel::deterministic, CCopasiVector< T >::end(), CStateTemplate::endFixed(), CStateTemplate::getIndex(), CModel::getModelType(), CModel::getReactions(), CModel::getStateTemplate(), CModelEntity::getStatus(), CCopasiObject::getValuePointer(), CModelEntity::getValueReference(), CTrajectoryMethodDsaLsodar::CPartition::intialize(), mAmu, mDoCorrection, mFirstReactionSpeciesIndex, CLsodaMethod::mMethodState, mNextReactionIndex, mNextReactionTime, mNumReactions, mPartition, mpLowerLimit, CLsodaMethod::mpModel, mpUpperLimit, mReactionDependencies, CModelEntity::REACTIONS, CVector< CType >::resize(), CModel::setState(), CCopasiVector< T >::size(), CLsodaMethod::start(), and CModel::updateSimulatedValues().

 {
   CLsodaMethod::start(initialState);
 
   if (mpModel->getModelType() == CModel::deterministic)
     mDoCorrection = true;
   else
     mDoCorrection = false;
 
   //initialize the vector of ints that contains the particle numbers
   //for the discrete simulation. This also floors all particle numbers in the model.
 
   mNumReactions = mpModel->getReactions().size();
 
   mReactionDependencies.resize(mNumReactions);
   mAmu.resize(mNumReactions);
   mAmu = 0.0;
 
   const CStateTemplate & StateTemplate = mpModel->getStateTemplate();
 
   CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
   CModelEntity *const* endEntity = StateTemplate.endFixed();
   C_FLOAT64 * pValue = mMethodState.beginIndependent();
 
   mFirstReactionSpeciesIndex = 0;
   size_t IndexSpecies = 1;
 
   for (; ppEntity != endEntity; ++ppEntity, ++pValue, ++IndexSpecies)
     {
       if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
         {
           *pValue = floor(*pValue + 0.5);
 
           if (mFirstReactionSpeciesIndex == 0 &&
               (*ppEntity)->getStatus() == CModelEntity::REACTIONS)
             {
               mFirstReactionSpeciesIndex = IndexSpecies;
             }
         }
     }
 
   // Update the model state so that the species are all represented by integers.
   mpModel->setState(mMethodState);
   mpModel->updateSimulatedValues(false); //for assignments
 
   // TODO CRITICAL Handle species of type ASSIGNMENT.
   // These need to be checked whether they are sufficiently close to an integer
 
   // Build the reaction dependencies
   C_FLOAT64 * pMethodStateValue = mMethodState.beginIndependent() - 1;
   size_t IndexReactions = 0;
   std::multimap< size_t, size_t> SpeciesToReactions;
 
   CCopasiVector< CReaction >::const_iterator it = mpModel->getReactions().begin();
   CCopasiVector< CReaction >::const_iterator end = mpModel->getReactions().end();
   std::vector< CReactionDependencies >::iterator itDependencies = mReactionDependencies.begin();
 
   for (; it  != end; ++it)
     {
       const CCopasiVector<CChemEqElement> & Balances = (*it)->getChemEq().getBalances();
       const CCopasiVector<CChemEqElement> & Substrates = (*it)->getChemEq().getSubstrates();
 
       // This reactions does not change anything we ignore it
       if (Balances.size() == 0 && Substrates.size() == 0)
         {
           continue;
         }
 
       itDependencies->mpParticleFlux = (C_FLOAT64 *)(*it)->getParticleFluxReference()->getValuePointer();
 
       itDependencies->mSpeciesMultiplier.resize(Balances.size());
       itDependencies->mMethodSpecies.resize(Balances.size());
       itDependencies->mModelSpecies.resize(Balances.size());
 
       std::set< size_t > SpeciesIndexSet;
 
       std::set< const CCopasiObject * > changed;
 
       // The time is always updated
       changed.insert(mpModel->getValueReference());
 
       CCopasiVector< CChemEqElement >::const_iterator itBalance = Balances.begin();
       CCopasiVector< CChemEqElement >::const_iterator endBalance = Balances.end();
 
       size_t Index = 0;
 
       for (; itBalance != endBalance; ++itBalance)
         {
           const CMetab * pMetab = (*itBalance)->getMetabolite();
 
           if (pMetab->getStatus() == CModelEntity::REACTIONS)
             {
               size_t SpeciesIndex = StateTemplate.getIndex(pMetab);
 
               itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
               itDependencies->mMethodSpecies[Index] = pMethodStateValue + SpeciesIndex;
               itDependencies->mModelSpecies[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
 
               changed.insert(pMetab->getValueReference());
 
               SpeciesToReactions.insert(std::make_pair(SpeciesIndex, IndexReactions));
               SpeciesIndexSet.insert(SpeciesIndex);
 
               Index++;
             }
         }
 
       // Correct allocation for metabolites which are not determined by reactions
       itDependencies->mSpeciesMultiplier.resize(Index, true);
       itDependencies->mMethodSpecies.resize(Index, true);
       itDependencies->mModelSpecies.resize(Index, true);
 
       itDependencies->mSubstrateMultiplier.resize(Substrates.size());
       itDependencies->mMethodSubstrates.resize(Substrates.size());
       itDependencies->mModelSubstrates.resize(Substrates.size());
 
       CCopasiVector< CChemEqElement >::const_iterator itSubstrate = Substrates.begin();
       CCopasiVector< CChemEqElement >::const_iterator endSubstrate = Substrates.end();
 
       Index = 0;
 
       for (; itSubstrate != endSubstrate; ++itSubstrate, ++Index)
         {
           const CMetab * pMetab = (*itSubstrate)->getMetabolite();
 
           size_t SpeciesIndex = StateTemplate.getIndex(pMetab);
 
           itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
           itDependencies->mMethodSubstrates[Index] = pMethodStateValue + SpeciesIndex;
           itDependencies->mModelSubstrates[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
 
           if (pMetab->getStatus() == CModelEntity::REACTIONS)
             {
               SpeciesIndexSet.insert(SpeciesIndex);
             }
         }
 
       itDependencies->mSpeciesIndex.resize(SpeciesIndexSet.size());
       size_t * pSpeciesIndex = itDependencies->mSpeciesIndex.array();
       std::set< size_t >::const_iterator itSet = SpeciesIndexSet.begin();
       std::set< size_t >::const_iterator endSet = SpeciesIndexSet.end();
 
       for (; itSet != endSet; ++itSet, ++pSpeciesIndex)
         {
           *pSpeciesIndex = *itSet;
         }
 
       std::set< const CCopasiObject * > dependend;
 
       CCopasiVector< CReaction >::const_iterator itReaction = mpModel->getReactions().begin();
       itDependencies->mDependentReactions.resize(mNumReactions);
 
       Index = 0;
       size_t Count = 0;
 
       for (; itReaction != end; ++itReaction, ++Index)
         {
           if ((*itReaction)->getParticleFluxReference()->dependsOn(changed))
             {
               dependend.insert((*itReaction)->getParticleFluxReference());
               itDependencies->mDependentReactions[Count] = Index;
 
               Count++;
             }
         }
 
       itDependencies->mDependentReactions.resize(Count, true);
 
       itDependencies->mCalculations = CCopasiObject::buildUpdateSequence(dependend, changed);
 
       ++itDependencies;
       ++IndexReactions;
     }
 
   mNumReactions = IndexReactions;
 
   mReactionDependencies.resize(mNumReactions);
   mAmu.resize(mNumReactions, true);
 
   mPartition.intialize(mReactionDependencies, SpeciesToReactions,
                        *mpLowerLimit, *mpUpperLimit, mMethodState);
 
   mNextReactionTime = std::numeric_limits< C_FLOAT64 >::infinity();
   mNextReactionIndex = C_INVALID_INDEX;
 
   calculatePropensities();
 
   return;
 }

void CTrajectoryMethodDsaLsodar::stateChanged ( )

virtual

Inform the trajectory method that the state has changed outside its control

Reimplemented from CLsodaMethod.

Definition at line 407 of file CTrajectoryMethodDsaLsodar.cpp.

References CLsodaMethod::mMethodState, CTrajectoryMethod::mpCurrentState, and CLsodaMethod::stateChanged().

 {
   *mpCurrentState = mMethodState;
   CLsodaMethod::stateChanged();
 }

CTrajectoryMethod::Status CTrajectoryMethodDsaLsodar::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 time step taken.

Parameters

const double & deltaT

Returns: Status status

Reimplemented from CLsodaMethod.

Definition at line 414 of file CTrajectoryMethodDsaLsodar.cpp.

References C_FLOAT64, doSingleStep(), CCopasiMessage::EXCEPTION, CState::getTime(), MCTrajectoryMethod, min, CLsodaMethod::mMethodState, CTrajectoryMethod::mpCurrentState, mpMaxSteps, CTrajectoryMethod::NORMAL, and CState::setTime().

 {
   // do several steps:
   C_FLOAT64 Time = mpCurrentState->getTime();
   C_FLOAT64 EndTime = Time + deltaT;
 
   C_FLOAT64 Tolerance = 100.0 * (fabs(EndTime) * std::numeric_limits< C_FLOAT64 >::epsilon() + std::numeric_limits< C_FLOAT64 >::min());
 
   size_t Steps = 0;
 
   while (fabs(Time - EndTime) > Tolerance)
     {
       Time += doSingleStep(Time, EndTime);
 
       if (++Steps > *mpMaxSteps)
         {
           CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 12);
         }
     }
 
   mMethodState.setTime(EndTime); // make sure the end time is exactly the time reach to next interval
   *mpCurrentState = mMethodState;
 
   return NORMAL;
 }