#include <CStochMethod.h>

Inheritance diagram for CStochMethod:

Collaboration diagram for CStochMethod:

Public Member Functions
	CStochMethod (const CStochMethod &src, const CCopasiContainer *pParent=NULL)

virtual bool	elevateChildren ()

virtual bool	isValidProblem (const CCopasiProblem *pProblem)

virtual void	start (const CState *initialState)

virtual Status	step (const double &deltaT)

	~CStochMethod ()

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)

virtual void	stateChanged ()

	~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
C_INT32	calculateAmu (size_t reaction_index)

	CStochMethod (const CCopasiContainer *pParent=NULL)

virtual C_FLOAT64	doSingleStep (C_FLOAT64 time, C_FLOAT64 endtime)=0

size_t	generateReactionIndex ()

C_FLOAT64	generateReactionTime ()

C_FLOAT64	generateReactionTime (size_t reaction_index)

virtual void	initMethod (C_FLOAT64 start_time)=0

void	setupDependencyGraphAndBalances ()

C_INT32	updatePropensities ()

C_INT32	updateSystemState (size_t reaction_index, const C_FLOAT64 &time)

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)

Static Protected Member Functions
static bool	modelHasAssignments (const CModel *pModel)

Protected Attributes
C_FLOAT64	mA0

std::vector< C_FLOAT64 >	mAmu

std::vector< C_FLOAT64 >	mAmuOld

CDependencyGraph	mDG

bool	mDoCorrection

size_t	mFirstMetabIndex

bool	mHasAssignments

std::vector< std::vector < CStochBalance > >	mLocalBalances

std::vector< std::vector < CStochBalance > >	mLocalSubstrates

std::vector< C_INT64 >	mNumbers

size_t	mNumReactions

CModel *	mpModel

CRandom *	mpRandomGenerator

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
std::set< std::string > *	getAffects (size_t reaction_index)

std::set< std::string > *	getDependsOn (size_t reaction_index)

void	initializeParameter ()

Static Private Member Functions
static C_INT32	checkModel (CModel *pmodel)

Private Attributes
C_INT32	mMaxBalance

C_INT64	mMaxIntBeforeStep

size_t	mMaxSteps

bool	mMaxStepsReached

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 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 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 47 of file CStochMethod.h.

Constructor & Destructor Documentation

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

protected

Default constructor.

Parameters

const CCopasiContainer * pParent (default: NULL)

Definition at line 48 of file CStochMethod.cpp.

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

                                                           :
   CTrajectoryMethod(CCopasiMethod::stochastic, pParent)
 {
   initializeParameter();
   mpRandomGenerator = CRandom::createGenerator(CRandom::mt19937);
 }

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

Copy constructor.

Parameters

const	CStochMethod & src,
const	CCopasiContainer * pParent (Default: NULL)

Definition at line 55 of file CStochMethod.cpp.

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

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

CStochMethod::~CStochMethod ( )

Destructor.

Definition at line 63 of file CStochMethod.cpp.

References mpRandomGenerator.

 {
   delete mpRandomGenerator;
   mpRandomGenerator = NULL;
 }

Member Function Documentation

C_INT32 CStochMethod::calculateAmu ( size_t reaction_index )

protected

Calculate one of the propensities

Returns: mFail

See Also: mFail

Definition at line 237 of file CStochMethod.cpp.

References C_FLOAT64, C_INT32, C_INT64, CModel::getReactions(), mAmu, mDoCorrection, mLocalSubstrates, mNumbers, and mpModel.

Referenced by updatePropensities(), and updateSystemState().

 {
   C_FLOAT64 rate_factor = mpModel->getReactions()[index]->calculateParticleFlux();
 
   if (rate_factor < 0.0)
     {
       // TODO CRITICAL Create a warning message
       rate_factor = 0.0;
     }
 
   if (!mDoCorrection)
     {
       mAmu[index] = rate_factor;
       return 0;
     }
 
   // We need the product of the cmu and hmu for this step.
   // We calculate this in one go, as there are fewer steps to
   // perform and we eliminate some possible rounding errors.
   C_FLOAT64 amu = 1; // initially
   //size_t total_substrates = 0;
   C_INT32 num_ident = 0;
   C_INT64 number = 0;
   C_INT64 lower_bound;
   // substrate_factor - The substrates, raised to their multiplicities,
   // multiplied with one another. If there are, e.g. m substrates of type m,
   // and n of type N, then substrate_factor = M^m * N^n.
   C_FLOAT64 substrate_factor = 1;
   // First, find the reaction associated with this index.
   // Keep a pointer to this.
   // Iterate through each substrate in the reaction
   const std::vector<CStochBalance> & substrates = mLocalSubstrates[index];
 
   int flag = 0;
 
   for (size_t i = 0; i < substrates.size(); i++)
     {
       num_ident = substrates[i].mMultiplicity;
 
       if (num_ident > 1)
         {
           flag = 1;
           number = mNumbers[substrates[i].mIndex];
           lower_bound = number - num_ident;
           substrate_factor = substrate_factor * pow((double) number, (int)(num_ident - 1));  //optimization
 
           number--; //optimization
 
           while (number > lower_bound)
             {
               amu *= number;
               number--;
             }
         }
     }
 
   if ((amu == 0) || (substrate_factor == 0))  // at least one substrate particle number is zero
     {
       mAmu[index] = 0;
       return 0;
     }
 
   // 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)!
 
   if (flag)
     {
       //cout << "Rate factor = " << rate_factor << endl;
       amu *= rate_factor / substrate_factor;
       mAmu[index] = amu;
     }
   else
     {
       mAmu[index] = rate_factor;
     }
 
   return 0;
 }

C_INT32 CStochMethod::checkModel ( CModel * pmodel )

staticprivate

This checks if a model is suitable for stochastic simulation. It returns a suggestion which method to use

Parameters

model The model to check

Returns: 1: direct method, 2: next reaction method, -1: no stochastic simulation possible

Definition at line 42 of file CStochMethod.cpp.

 {
   // Here several checks will be performed to validate the model
   return 2; // suggest next reaction method
 }

virtual C_FLOAT64 CStochMethod::doSingleStep	(	C_FLOAT64	time,
		C_FLOAT64	endtime
	)

protectedpure virtual

Do one iteration of the simulation

Returns: Current simulation time or -1 if error.

Implemented in CStochNextReactionMethod.

Referenced by step().

bool CStochMethod::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 104 of file CStochMethod.cpp.

References initializeParameter().

 {
   initializeParameter();
   return true;
 }

size_t CStochMethod::generateReactionIndex ( )

protected

Generate the index of a putative reaction.

Returns: The reaction index

Definition at line 389 of file CStochMethod.cpp.

References C_FLOAT64, CRandom::getRandomOO(), CModel::getReactions(), mA0, mAmu, mpModel, mpRandomGenerator, and CCopasiVector< T >::size().

 {
   C_FLOAT64 rand1 = mpRandomGenerator->getRandomOO() * mA0;
   C_FLOAT64 sum = 0;
   size_t index = 0;
 
   while (index < (mpModel->getReactions().size() - 1))
     {
       sum += mAmu[index] /* /mA0 */;
 
       if (rand1 <= sum)
         {return index;}
 
       index++;
     }
 
   return index;
 }

C_FLOAT64 CStochMethod::generateReactionTime ( )

protected

Generate the putative time taken before any reaction takes place

Returns: The time before the reaction

Definition at line 408 of file CStochMethod.cpp.

References C_FLOAT64, CRandom::getRandomOO(), mA0, and mpRandomGenerator.

Referenced by CStochNextReactionMethod::setupPriorityQueue(), and CStochNextReactionMethod::updatePriorityQueue().

 {
   if (mA0 == 0) return std::numeric_limits<C_FLOAT64>::infinity();
 
   C_FLOAT64 rand2 = mpRandomGenerator->getRandomOO();
   return - 1 * log(rand2) / mA0;
 }

C_FLOAT64 CStochMethod::generateReactionTime ( size_t reaction_index )

protected

Generate the putative time taken before a special reaction takes place

Returns: The time before the reaction

Definition at line 416 of file CStochMethod.cpp.

References C_FLOAT64, CRandom::getRandomOO(), mAmu, and mpRandomGenerator.

 {
   if (mAmu[reaction_index] == 0) return std::numeric_limits<C_FLOAT64>::infinity();
 
   C_FLOAT64 rand2 = mpRandomGenerator->getRandomOO();
   return - 1 * log(rand2) / mAmu[reaction_index];
 }

std::set< std::string > * CStochMethod::getAffects ( size_t reaction_index )

private

Get the set of metabolites which change number when a given reaction is executed.

Parameters

reaction_index The index of the reaction being executed.

Returns: The set of affected metabolites.

Definition at line 560 of file CStochMethod.cpp.

References CModelEntity::FIXED, CCopasiParameter::getKey(), CModel::getReactions(), mpModel, and CCopasiVector< T >::size().

Referenced by setupDependencyGraphAndBalances().

 {
   std::set<std::string> *retset = new std::set<std::string>;
 
   // Get the balances  associated with the reaction at this index
   // XXX We first get the chemical equation, then the balances, since the getBalances method in CReaction is unimplemented!
   const CCopasiVector<CChemEqElement> & balances = mpModel->getReactions()[reaction_index]->getChemEq().getBalances();
 
   for (size_t i = 0; i < balances.size(); i++)
     {
       if (!balances[i]->getMetabolite()) continue;
 
       if (fabs(balances[i]->getMultiplicity()) >= 0.1)
         if (balances[i]->getMetabolite()->getStatus() != CModelEntity::FIXED)
           {
             retset->insert(balances[i]->getMetabolite()->getKey());
           }
     }
 
   return retset;
 }

std::set< std::string > * CStochMethod::getDependsOn ( size_t reaction_index )

private

Get the set of metabolites on which a given reaction depends.

Parameters

reaction_index The index of the reaction being executed.

Returns: The set of metabolites depended on.

Definition at line 534 of file CStochMethod.cpp.

References CModel::getReactions(), mpModel, CFunctionParameter::PARAMETER, and CCopasiVector< T >::size().

Referenced by setupDependencyGraphAndBalances().

 {
   std::set<std::string> *retset = new std::set<std::string>;
 
   size_t i, imax = mpModel->getReactions()[reaction_index]->getFunctionParameters().size();
   size_t j, jmax;
 
   for (i = 0; i < imax; ++i)
     {
       if (mpModel->getReactions()[reaction_index]->getFunctionParameters()[i]->getUsage() == CFunctionParameter::PARAMETER)
         continue;
 
       //metablist = mpModel->getReactions()[reaction_index]->getParameterMappingMetab(i);
       const std::vector <std::string> & metabKeylist =
         mpModel->getReactions()[reaction_index]->getParameterMappings()[i];
       jmax = metabKeylist.size();
 
       for (j = 0; j < jmax; ++j)
         {
           retset->insert(metabKeylist[j]);
         }
     }
 
   return retset;
 }

void CStochMethod::initializeParameter ( )

private

Intialize the method parameter

Definition at line 69 of file CStochMethod.cpp.

References CCopasiParameterGroup::assertParameter(), CCopasiParameter::BOOL, C_INT32, CCopasiParameterGroup::getParameter(), CCopasiParameter::getValue(), CCopasiParameter::INT, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pINT, CCopasiParameter::Value::pUINT, CCopasiParameterGroup::removeParameter(), CCopasiParameterGroup::setValue(), and CCopasiParameter::UINT.

Referenced by CStochMethod(), and elevateChildren().

 {
   CCopasiParameter *pParm;
 
   assertParameter("Max Internal Steps", CCopasiParameter::INT, (C_INT32) 1000000);
   assertParameter("Subtype", CCopasiParameter::UINT, (unsigned C_INT32) 2);
   assertParameter("Use Random Seed", CCopasiParameter::BOOL, false);
   assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) 1);
 
   // Check whether we have a method with the old parameter names
   if ((pParm = getParameter("STOCH.MaxSteps")) != NULL)
     {
       setValue("Max Internal Steps", *pParm->getValue().pINT);
       removeParameter("STOCH.MaxSteps");
 
       if ((pParm = getParameter("STOCH.Subtype")) != NULL)
         {
           setValue("Subtype", *pParm->getValue().pUINT);
           removeParameter("STOCH.Subtype");
         }
 
       if ((pParm = getParameter("STOCH.UseRandomSeed")) != NULL)
         {
           setValue("Use Random Seed", *pParm->getValue().pBOOL);
           removeParameter("STOCH.UseRandomSeed");
         }
 
       if ((pParm = getParameter("STOCH.RandomSeed")) != NULL)
         {
           setValue("Random Seed", *pParm->getValue().pUINT);
           removeParameter("STOCH.RandomSeed");
         }
     }
 }

virtual void CStochMethod::initMethod ( C_FLOAT64 start_time )

protectedpure virtual

Initialization.

Implemented in CStochNextReactionMethod.

Referenced by start().

bool CStochMethod::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 583 of file CStochMethod.cpp.

References CCopasiMessage::ERROR, CModel::getCompartments(), CTrajectoryProblem::getDuration(), CModel::getEvents(), CModel::getMetabolites(), CCopasiProblem::getModel(), CModel::getModelValues(), CModel::getNumMetabs(), CModel::getNumModelValues(), CModel::getTotSteps(), CCopasiParameter::getValue(), CTrajectoryMethod::isValidProblem(), MCTrajectoryMethod, CModelEntity::ODE, 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;
     }
 
   if (pTP->getModel()->getTotSteps() < 1)
     {
       //at least one reaction necessary
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 17);
       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;
         }
     }
 
   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;
         }
     }
 
   //   if (modelHasAssignments(pTP->getModel()))
   //     {
   //       CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 19);
   //       return false;
   //}
   //  this test is disabled since we have preliminary support for assignments
 
   //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;
     }
 
   if (* getValue("Max Internal Steps").pINT <= 0)
     {
       //max steps should be at least 1
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 15);
       return false;
     }
 
   //events are not supported at the moment
   if (pTP->getModel()->getEvents().size() > 0)
     {
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 23);
       return false;
     }
 
   return true;
 }

bool CStochMethod::modelHasAssignments ( const CModel * pModel )

staticprotected

tests if the model contains a global value with an assignment rule that is used in calculations

Definition at line 676 of file CStochMethod.cpp.

References CModelEntity::ASSIGNMENT, CModel::getCompartments(), CModel::getMetabolites(), CModel::getModelValues(), CModel::getNumMetabs(), CModel::getNumModelValues(), and CCopasiVector< T >::size().

Referenced by start().

 {
   size_t i, imax = pModel->getNumModelValues();
 
   for (i = 0; i < imax; ++i)
     {
       if (pModel->getModelValues()[i]->getStatus() == CModelEntity::ASSIGNMENT)
         if (pModel->getModelValues()[i]->isUsed())
           {
             //used assignment found
             return true;
           }
     }
 
   imax = pModel->getNumMetabs();
 
   for (i = 0; i < imax; ++i)
     {
       if (pModel->getMetabolites()[i]->getStatus() == CModelEntity::ASSIGNMENT)
         if (pModel->getMetabolites()[i]->isUsed())
           {
             //used assignment found
             return true;
           }
     }
 
   imax = pModel->getCompartments().size();
 
   for (i = 0; i < imax; ++i)
     {
       if (pModel->getCompartments()[i]->getStatus() == CModelEntity::ASSIGNMENT)
         if (pModel->getCompartments()[i]->isUsed())
           {
             //used assignment found
             return true;
           }
     }
 
   return false;
 }

void CStochMethod::setupDependencyGraphAndBalances ( )

protected

Set up the dependency graph and the balances

Definition at line 424 of file CStochMethod.cpp.

References CDependencyGraph::addDependent(), C_INT32, CDependencyGraph::clear(), CModelEntity::FIXED, getAffects(), getDependsOn(), CCopasiVector< T >::getIndex(), CModel::getMetabolitesX(), CModel::getReactions(), max, mDG, CStochBalance::mIndex, mLocalBalances, mLocalSubstrates, mMaxBalance, mMaxIntBeforeStep, mMaxSteps, CStochBalance::mMultiplicity, mNumReactions, mpModel, CDependencyGraph::resize(), and CCopasiVector< T >::size().

Referenced by start().

 {
   mDG.clear();
   std::vector< std::set<std::string>* > DependsOn;
   std::vector< std::set<std::string>* > Affects;
   size_t i, j;
   // Do for each reaction:
 
   for (i = 0; i < mNumReactions; i++)
     {
       // Get the set of metabolites  which affect the value of amu for this
       // reaction i.e. the set on which amu depends. This may be  more than
       // the set of substrates, since the kinetics can involve other
       // reactants, e.g. catalysts. We thus need to step through the
       // rate function and pick out every reactant which can vary.
       DependsOn.push_back(getDependsOn(i));
       // Get the set of metabolites which are affected when this reaction takes place
       Affects.push_back(getAffects(i));
     }
 
   // For each possible pair of reactions i and j, if the intersection of
   // Affects(i) with DependsOn(j) is non-empty, add a dependency edge from i to j.
   mDG.resize(mNumReactions);
 
   for (i = 0; i < mNumReactions; i++)
     {
       for (j = 0; j < mNumReactions; j++)
         {
           // Determine whether the intersection of these two sets is non-empty
           // Could also do this with set_intersection generic algorithm, but that
           // would require operator<() to be defined on the set elements.
 
           std::set<std::string>::iterator iter = Affects[i]->begin();
 
           for (; iter != Affects[i]->end(); iter++)
             {
               if (DependsOn[j]->count(*iter))
                 {
                   // The set intersection is non-empty
                   mDG.addDependent(i, j);
                   break;
                 }
             }
         }
 
       // Ensure that self edges are included
       //mDG.addDependent(i, i);
     }
 
   // Create local copy of balances and substrates list
   CStochBalance bb;
   C_INT32 maxBalance = 0;
 
   mLocalBalances.clear();
   mLocalBalances.resize(mNumReactions);
   mLocalSubstrates.clear();
   mLocalSubstrates.resize(mNumReactions);
 
   for (i = 0; i < mNumReactions; i++)
     {
       const CCopasiVector<CChemEqElement> * bbb;
 
       bbb = &mpModel->getReactions()[i]->getChemEq().getBalances();
 
       //TODO clear old local balances and substrates
       for (j = 0; j < bbb->size(); j++)
         {
           //bb.mIndex = mpModel->getMetabolites().getIndex((*bbb)[j]->getMetabolite().getObjectName(());
           //bb.mIndex = mpModel->findMetabByKey((*bbb)[j]->getMetaboliteKey());
           assert((*bbb)[j]->getMetabolite());
           bb.mIndex = mpModel->getMetabolitesX().getIndex((*bbb)[j]->getMetabolite());
           bb.mMultiplicity = static_cast<C_INT32>(floor((*bbb)[j]->getMultiplicity() + 0.5));
 
           if (((*bbb)[j]->getMetabolite()->getStatus()) != CModelEntity::FIXED)
             {
               if (bb.mMultiplicity > maxBalance) maxBalance = bb.mMultiplicity;
 
               mLocalBalances[i].push_back(bb);
             }
         }
 
       bbb = &mpModel->getReactions()[i]->getChemEq().getSubstrates();
 
       for (j = 0; j < bbb->size(); j++)
         {
           //bb.mIndex = mpModel->getMetabolites().getIndex((*bbb)[j]->getMetabolite().getObjectName(());
           //bb.mIndex = mpModel->findMetabByKey((*bbb)[j]->getMetaboliteKey());
           assert((*bbb)[j]->getMetabolite());
           bb.mIndex = mpModel->getMetabolitesX().getIndex((*bbb)[j]->getMetabolite());
           bb.mMultiplicity = static_cast<C_INT32>(floor((*bbb)[j]->getMultiplicity() + 0.5));
 
           if (1)
             {
               mLocalSubstrates[i].push_back(bb);
             }
         }
     }
 
   mMaxBalance = maxBalance;
   mMaxIntBeforeStep = std::numeric_limits< C_INT64 >::max() - 1 - mMaxSteps * mMaxBalance;
 
   // Delete the memory allocated in getDependsOn() and getAffects()
   // since this is allocated in other functions.
   for (i = 0; i < mNumReactions; i++)
     {
       delete DependsOn[i];
       delete Affects[i];
     }
 }

void CStochMethod::start ( const CState * initialState )

virtual

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

Parameters

const CState * initialState

Reimplemented from CTrajectoryMethod.

Definition at line 160 of file CStochMethod.cpp.

References C_FLOAT64, C_INT32, C_INT64, CModel::deterministic, CStateTemplate::getIndex(), CModel::getMetabolitesX(), CCopasiProblem::getModel(), CModel::getModelType(), CModel::getReactions(), CModel::getStateTemplate(), CState::getTime(), CCopasiParameter::getValue(), CRandom::initialize(), initMethod(), mAmu, mAmuOld, mDoCorrection, mFirstMetabIndex, mHasAssignments, mMaxSteps, mMaxStepsReached, mNumbers, mNumReactions, modelHasAssignments(), CTrajectoryMethod::mpCurrentState, mpModel, CTrajectoryMethod::mpProblem, mpRandomGenerator, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pINT, CCopasiParameter::Value::pUINT, setupDependencyGraphAndBalances(), CCopasiVector< T >::size(), updatePropensities(), and CModel::updateSimulatedValues().

 {
   /* get configuration data */
   mMaxSteps = * getValue("Max Internal Steps").pINT;
 
   bool useRandomSeed = * getValue("Use Random Seed").pBOOL;
   unsigned C_INT32 randomSeed = * getValue("Random Seed").pUINT;
 
   if (useRandomSeed) mpRandomGenerator->initialize(randomSeed);
 
   //mpCurrentState is initialized. This state is not used internally in the
   //stochastic solver, but it is used for returning the result after each step.
   *mpCurrentState = *initialState;
 
   mpModel = mpProblem->getModel();
   assert(mpModel);
 
   if (mpModel->getModelType() == CModel::deterministic)
     mDoCorrection = true;
   else
     mDoCorrection = false;
 
   mHasAssignments = modelHasAssignments(mpModel);
 
   size_t i;
 
   //initialize the vector of ints that contains the particle numbers
   //for the discrete simulation. This also floors all particle numbers in the model.
   mNumbers.resize(mpModel->getMetabolitesX().size());
 
   for (i = 0; i < mNumbers.size(); ++i)
     {
       mNumbers[i] = (C_INT64) mpModel->getMetabolitesX()[i]->getValue();
       mpModel->getMetabolitesX()[i]->setValue((C_FLOAT64) mNumbers[i]);
       mpModel->getMetabolitesX()[i]->refreshConcentration();
     }
 
   mFirstMetabIndex = mpModel->getStateTemplate().getIndex(mpModel->getMetabolitesX()[0]);
 
   mpModel->updateSimulatedValues(false); //for assignments
   //mpModel->updateNonSimulatedValues(); //for assignments
 
   mNumReactions = mpModel->getReactions().size();
 
   mAmu.clear(); mAmuOld.clear();
 
   for (i = 0; i < mNumReactions; i++)
     {
       mAmu.push_back(0);
       mAmuOld.push_back(0);
     }
 
   setupDependencyGraphAndBalances();
   updatePropensities();
 
   // call init of the specific simulation method
   initMethod(mpCurrentState->getTime());
 
   mMaxStepsReached = false;
   return;
 }

CTrajectoryMethod::Status CStochMethod::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

Returns: Status status

Reimplemented from CTrajectoryMethod.

Definition at line 110 of file CStochMethod.cpp.

References CState::beginIndependent(), C_FLOAT64, doSingleStep(), CCopasiMessage::EXCEPTION, CModel::getMetabolitesX(), CCopasiProblem::getModel(), CModel::getNumVariableMetabs(), CModel::getState(), CState::getTime(), MCTrajectoryMethod, mFirstMetabIndex, mMaxIntBeforeStep, mMaxSteps, CTrajectoryMethod::mpCurrentState, CTrajectoryMethod::mpProblem, CTrajectoryMethod::NORMAL, and CState::setTime().

 {
   // write the current state to the model:
   //mpProblem->getModel()->setState(mpCurrentState); //?
 
   // check for possible overflows:
   size_t i;
   size_t imax;
 
   // :TODO: Bug 774: This assumes that the number of variable metabs is the number
   // of metabs determined by reaction. In addition they are expected at the beginning of the
   // MetabolitesX which is not the case if we have metabolites of type ODE.
   for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++)
     if (mpProblem->getModel()->getMetabolitesX()[i]->getValue() >= mMaxIntBeforeStep)
       {
         CCopasiMessage(CCopasiMessage::EXCEPTION, "at least one particle number got to big.");
         // TODO:throw exception or something like that
       }
 
   // do several steps:
   C_FLOAT64 time = mpCurrentState->getTime();
   C_FLOAT64 endtime = time + deltaT;
 
   size_t Steps = 0;
 
   while (time < endtime)
     {
       time = doSingleStep(time, endtime);
       mpCurrentState->setTime(time);
 
       if (++Steps > mMaxSteps)
         {
           CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 12);
         }
     }
 
   *mpCurrentState = mpProblem->getModel()->getState();
   mpCurrentState->setTime(time);
 
   // get back the particle numbers:
 
   /* Set the variable Metabolites */
   C_FLOAT64 * Dbl = mpCurrentState->beginIndependent() + mFirstMetabIndex - 1;
 
   for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, Dbl++)
     *Dbl = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
 
   return NORMAL;
 }

C_INT32 CStochMethod::updatePropensities ( )

protected

Calculate the propensities for all reactions

Returns: mFail

See Also: mFail

Definition at line 222 of file CStochMethod.cpp.

References calculateAmu(), mA0, mAmu, mAmuOld, and mNumReactions.

Referenced by start(), and updateSystemState().

 {
   //mA0Old = mA0;
   mA0 = 0;
 
   for (size_t i = 0; i < mNumReactions; i++)
     {
       mAmuOld[i] = mAmu[i];
       calculateAmu(i);
       mA0 += mAmu[i];
     }
 
   return 0;
 }

C_INT32 CStochMethod::updateSystemState	(	size_t	reaction_index,
		const C_FLOAT64 &	time
	)

protected

Update the particle numbers according to which reaction ocurred

Returns: mFail

See Also: mFail

Definition at line 317 of file CStochMethod.cpp.

References CModelEntity::ASSIGNMENT, C_FLOAT64, C_INT64, calculateAmu(), CDependencyGraph::getDependents(), CModel::getMetabolitesX(), mA0, mAmu, mAmuOld, mDG, mHasAssignments, mLocalBalances, mNumbers, mpModel, CMetab::refreshConcentration(), CModel::setTime(), CModelEntity::setValue(), updatePropensities(), and CModel::updateSimulatedValues().

Referenced by CStochNextReactionMethod::doSingleStep().

 {
   // Change the particle numbers according to which step took place.
   // First, get the vector of balances in the reaction we've got.
   // (This vector expresses the number change of each metabolite
   // in the reaction.) Then step through each balance, using its
   // multiplicity to calculate a new value for the associated
   // metabolite. Finally, update the metabolite.
 
   const std::vector<CStochBalance> & bals = mLocalBalances[rxn];
 
   std::vector<CStochBalance>::const_iterator bi;
   CMetab* pTmpMetab;
 
   for (bi = bals.begin(); bi != bals.end(); bi++)
     {
       mNumbers[bi->mIndex] = mNumbers[bi->mIndex] + bi->mMultiplicity;
       pTmpMetab = mpModel->getMetabolitesX()[bi->mIndex];
       pTmpMetab->setValue((C_FLOAT64) mNumbers[bi->mIndex]);
       pTmpMetab->refreshConcentration();
     }
 
   if (mHasAssignments)
     {
       // this is less efficient but can deal with assignments.
       //TODO: handle dependencies for assignments also.
       mpModel->setTime(time);
       mpModel->updateSimulatedValues(false);
 
       //now potentially species with assignments have non integer
       //particle numbers. This needs to be rounded. Also the updated
       //particle numbers need to be copied to the vector of integers.
       //(the integer values may be used to calculate the propensities for
       //higher order kinetics).
       size_t i, imax = mNumbers.size();
 
       for (i = 0; i < imax; ++i)
         {
           if (mpModel->getMetabolitesX()[i]->getStatus() == CModelEntity::ASSIGNMENT)
             {
               mNumbers[i] = (C_INT64) mpModel->getMetabolitesX()[i]->getValue();
               mpModel->getMetabolitesX()[i]->setValue((C_FLOAT64) mNumbers[i]);
               mpModel->getMetabolitesX()[i]->refreshConcentration();
             }
         }
 
       //now the propensities can be updated
       updatePropensities();
     }
   else
     {
       const std::set<size_t> & dep_nodes = mDG.getDependents(rxn);
 
       std::set<size_t>::const_iterator it;
       size_t ii;
 
       for (it = dep_nodes.begin(); it != dep_nodes.end(); it++)
         {
           ii = *it;
           mAmuOld[ii] = mAmu[ii];
           calculateAmu(ii);
         }
 
       //mA0Old = mA0;
 
       mA0 = 0;
       mA0 = std::accumulate(mAmu.begin(), mAmu.end(), mA0);
     }
 
   return 0;
 }

Friends And Related Function Documentation

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

friend

Member Data Documentation

C_FLOAT64 CStochMethod::mA0

protected

The sum of the propensities

Definition at line 62 of file CStochMethod.h.

Referenced by generateReactionIndex(), generateReactionTime(), updatePropensities(), and updateSystemState().

std::vector<C_FLOAT64> CStochMethod::mAmu

protected

The propensities for the reactions

Definition at line 57 of file CStochMethod.h.

Referenced by calculateAmu(), generateReactionIndex(), generateReactionTime(), start(), CStochNextReactionMethod::updatePriorityQueue(), updatePropensities(), and updateSystemState().

std::vector<C_FLOAT64> CStochMethod::mAmuOld

protected

The stored propensities for the reactions before the last update

Definition at line 67 of file CStochMethod.h.

Referenced by start(), CStochNextReactionMethod::updatePriorityQueue(), updatePropensities(), and updateSystemState().

CDependencyGraph CStochMethod::mDG

protected

The graph of reactions and their dependent reactions. When a reaction is executed, the propensities for each of its dependents must be updated.

Definition at line 252 of file CStochMethod.h.

Referenced by setupDependencyGraphAndBalances(), CStochNextReactionMethod::updatePriorityQueue(), and updateSystemState().

bool CStochMethod::mDoCorrection

protected

indicates if the correction N^2 -> N*(N-1) should be performed

Definition at line 72 of file CStochMethod.h.

Referenced by calculateAmu(), and start().

size_t CStochMethod::mFirstMetabIndex

protected

index of first metab in a CState

Definition at line 274 of file CStochMethod.h.

Referenced by start(), and step().

bool CStochMethod::mHasAssignments

protected

Indicates whether the model has global quantities with assignment rules. If it has, we will use a less efficient way to update the model state to handle this.

Definition at line 79 of file CStochMethod.h.

Referenced by start(), CStochNextReactionMethod::updatePriorityQueue(), and updateSystemState().

std::vector< std::vector <CStochBalance> > CStochMethod::mLocalBalances

protected

The balances of the reactions as integers

Definition at line 257 of file CStochMethod.h.

Referenced by setupDependencyGraphAndBalances(), and updateSystemState().

std::vector< std::vector <CStochBalance> > CStochMethod::mLocalSubstrates

protected

the substrates of each reaction with their multiplicities

Definition at line 262 of file CStochMethod.h.

Referenced by calculateAmu(), and setupDependencyGraphAndBalances().