#include <CTauLeapMethod.h>

Inheritance diagram for CTauLeapMethod:

Collaboration diagram for CTauLeapMethod:

Classes
class	CReactionDependencies

Public Member Functions
	CTauLeapMethod (const CTauLeapMethod &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)

	~CTauLeapMethod ()

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
const C_FLOAT64 &	calculateAmu (const size_t &index)

void	cleanup ()

	CTauLeapMethod (const CCopasiContainer *pParent=NULL)

C_FLOAT64	doSingleStep (C_FLOAT64 ds)

void	updatePropensities ()

bool	updateSystem ()

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

CVector< C_FLOAT64 >	mAvgDX

bool	mDoCorrection

C_FLOAT64	mEpsilon

size_t	mFirstReactionSpeciesIndex

CVector< C_FLOAT64 >	mK

unsigned C_INT32	mMaxSteps

CState	mMethodState

size_t	mNumReactions

size_t	mNumReactionSpecies

CModel *	mpModel

CRandom *	mpRandomGenerator

unsigned C_INT32	mRandomSeed

std::vector < CReactionDependencies >	mReactionDependencies

CVector< C_FLOAT64 >	mSigDX

bool	mUseRandomSeed

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 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 CTauLeapMethod.h.

Constructor & Destructor Documentation

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

Copy constructor

Parameters

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

Definition at line 112 of file CTauLeapMethod.cpp.

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

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

CTauLeapMethod::~CTauLeapMethod ( )

Destructor.

Definition at line 124 of file CTauLeapMethod.cpp.

References cleanup(), and DESTRUCTOR_TRACE.

 {
   cleanup();
   DESTRUCTOR_TRACE;
 }

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

protected

Default constructor.

Parameters

const CCopasiContainer * pParent (default: NULL)

Default constructor.

Definition at line 104 of file CTauLeapMethod.cpp.

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

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

Member Function Documentation

const C_FLOAT64 & CTauLeapMethod::calculateAmu ( const size_t & index )

protected

Calculate one of the propensities

Parameters

const size_t & index

Returns: const C_FLOAT64 & amu

Definition at line 473 of file CTauLeapMethod.cpp.

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

Referenced by updatePropensities().

 {
   const 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 Amu;
     }
 
   C_FLOAT64 SubstrateMultiplier = 1.0;
   C_FLOAT64 SubstrateDevisor = 1.0;
   C_FLOAT64 Multiplicity;
   C_FLOAT64 LowerBound;
   C_FLOAT64 Number;
 
   bool ApplyCorrection = false;
 
   const C_FLOAT64 * pMultiplicity = Dependencies.mSubstrateMultiplier.array();
   const C_FLOAT64 * pEndMultiplicity = pMultiplicity + Dependencies.mSubstrateMultiplier.size();
   C_FLOAT64 *const* ppLocalSubstrate = Dependencies.mMethodSubstrates.array();
   C_FLOAT64 *const* ppModelSubstrate = Dependencies.mModelSubstrates.array();
 
   for (; pMultiplicity != pEndMultiplicity; ++pMultiplicity, ++ppLocalSubstrate, ++ppModelSubstrate)
     {
       Multiplicity = *pMultiplicity;
 
       // 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;
     }
   else if (ApplyCorrection)
     {
       Amu *= SubstrateMultiplier / SubstrateDevisor;
     }
 
   return Amu;
 }

void CTauLeapMethod::cleanup ( )

protected

Cleans up memory, etc.

Initializes the solver and sets the model to be used.

Parameters

model A reference to an instance of a CModel Cleans up memory, etc.

Definition at line 357 of file CTauLeapMethod.cpp.

References mpModel, and mpRandomGenerator.

Referenced by ~CTauLeapMethod().

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

C_FLOAT64 CTauLeapMethod::doSingleStep ( C_FLOAT64 ds )

protected

Simulates the system over the next interval of time. The timestep is given as argument.

Parameters

ds	A C_FLOAT64 specifying the timestep

Definition at line 371 of file CTauLeapMethod.cpp.

References CVectorCore< CType >::array(), CState::beginIndependent(), C_FLOAT64, CCopasiMessage::EXCEPTION, CRandom::getRandomCC(), CRandom::getRandomPoisson(), mAmu, mAvgDX, MCTrajectoryMethod, mEpsilon, mFirstReactionSpeciesIndex, min, mK, mMethodState, mNumReactions, mNumReactionSpecies, mpRandomGenerator, mReactionDependencies, mSigDX, updatePropensities(), and updateSystem().

Referenced by step().

 {
   C_FLOAT64 Lambda, Tmp, Tau, Tau1, Tau2;
 
   updatePropensities();
 
   mAvgDX = 0.0;
   mSigDX = 0.0;
 
   std::vector< CReactionDependencies >::const_iterator itReaction = mReactionDependencies.begin();
   const C_FLOAT64 * pAmu = mAmu.array();
   const C_FLOAT64 * pAmuEnd = pAmu + mNumReactions;
 
   for (; pAmu != pAmuEnd; ++pAmu, ++itReaction)
     {
       const C_FLOAT64 * pMultiplicity = itReaction->mSpeciesMultiplier.array();
       const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + itReaction->mSpeciesMultiplier.size();
       const size_t * pIndex = itReaction->mMethodSpeciesIndex.array();
 
       for (; pMultiplicity != pMultiplicityEnd; ++pMultiplicity, ++pIndex)
         {
           mAvgDX[*pIndex] += *pMultiplicity * *pAmu;
           mSigDX[*pIndex] += *pMultiplicity * *pMultiplicity * *pAmu;
         }
     }
 
   Tau1 = Tau2 = std::numeric_limits< C_FLOAT64 >::infinity();
 
   const C_FLOAT64 * pNumber = mMethodState.beginIndependent() + mFirstReactionSpeciesIndex - 1;
   const C_FLOAT64 * pNumberEnd = pNumber + mNumReactionSpecies;
   C_FLOAT64 * pAvgDX = mAvgDX.array();
   C_FLOAT64 * pSigDX = mSigDX.array();
 
   for (; pNumber != pNumberEnd; ++pNumber, ++pAvgDX, ++pSigDX)
     {
       if ((Tmp = mEpsilon * fabs(*pNumber)) < 1.0)
         Tmp = 1.0;
 
       *pAvgDX = Tmp / fabs(*pAvgDX);
       *pSigDX = (Tmp * Tmp) / fabs(*pSigDX);
 
       if (Tau1 > *pAvgDX)
         Tau1 = *pAvgDX;
 
       if (Tau2 > *pSigDX)
         Tau2 = *pSigDX;
     }
 
   Tau = std::min(Tau1, Tau2);
 
   if (ds < Tau)
     Tau = ds;
 
   pAmu = mAmu.array();
   C_FLOAT64 * pK = mK.array();
   C_FLOAT64 * pKEnd = pK + mNumReactions;
 
   for (; pAmu != pAmuEnd; ++pAmu, ++pK)
     {
       Lambda = *pAmu * Tau;
 
       if (Lambda < 0.0)
         CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 10);
       else if (Lambda > 2.0e9)
         CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 26);
 
       *pK = mpRandomGenerator->getRandomPoisson(Lambda);
     }
 
   while (!updateSystem())
     {
       Tau *= 0.5;
       pK = mK.array();
 
       for (; pK != pKEnd; ++pK)
         {
           *pK *= 0.5;
 
           if (*pK < floor(*pK + 0.75))
             {
               *pK += mpRandomGenerator->getRandomCC() < 0.5 ? - 0.5 : 0.5;
             }
         }
     }
 
 
   return Tau;
 }

bool CTauLeapMethod::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 158 of file CTauLeapMethod.cpp.

References initializeParameter().

 {
   initializeParameter();
   return true;
 }

void CTauLeapMethod::initializeParameter ( )

private

Initialize the method parameter

Definition at line 130 of file CTauLeapMethod.cpp.

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

Referenced by CTauLeapMethod(), and elevateChildren().

 {
   CCopasiParameter *pParm;
 
   assertParameter("Epsilon", CCopasiParameter::DOUBLE, (C_FLOAT64) 0.001);
   assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) 10000);
   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("TAULEAP.Tau")) != NULL)
     {
       removeParameter("TAULEAP.Tau");
 
       if ((pParm = getParameter("TAULEAP.UseRandomSeed")) != NULL)
         {
           setValue("Use Random Seed", *pParm->getValue().pBOOL);
           removeParameter("TAULEAP.UseRandomSeed");
         }
 
       if ((pParm = getParameter("TAULEAP.RandomSeed")) != NULL)
         {
           setValue("Random Seed", *pParm->getValue().pUINT);
           removeParameter("TAULEAP.RandomSeed");
         }
     }
 }

bool CTauLeapMethod::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 585 of file CTauLeapMethod.cpp.

References CCopasiMessage::ERROR, CModel::getCompartments(), CTrajectoryProblem::getDuration(), CModel::getEvents(), CModel::getMetabolites(), CCopasiProblem::getModel(), CModel::getModelValues(), CModel::getNumMetabs(), CModel::getNumModelValues(), CModel::getTotSteps(), CTrajectoryMethod::isValidProblem(), MCTrajectoryMethod, CModelEntity::ODE, 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 ODE 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;
         }
     }
 
   //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;
     }
 
   //events are not supported at the moment
   if (pTP->getModel()->getEvents().size() > 0)
     {
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 23);
       return false;
     }
 
   return true;
 }

void CTauLeapMethod::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 195 of file CTauLeapMethod.cpp.

References CCopasiVector< T >::begin(), CStateTemplate::beginIndependent(), CState::beginIndependent(), C_FLOAT64, C_INT32, CModel::deterministic, CCopasiVector< T >::end(), CStateTemplate::endFixed(), CStateTemplate::getIndex(), CCopasiProblem::getModel(), CModel::getModelType(), CModel::getNumDependentReactionMetabs(), CModel::getNumIndependentReactionMetabs(), CModel::getReactions(), CModel::getState(), CModel::getStateTemplate(), CModelEntity::getStatus(), CCopasiParameter::getValue(), CCopasiObject::getValuePointer(), CModelEntity::getValueReference(), CRandom::initialize(), mAmu, mAvgDX, mDoCorrection, mEpsilon, mFirstReactionSpeciesIndex, mK, mMaxSteps, mMethodState, mNumReactions, mNumReactionSpecies, CTrajectoryMethod::mpCurrentState, mpModel, CTrajectoryMethod::mpProblem, mpRandomGenerator, mRandomSeed, mReactionDependencies, mSigDX, mUseRandomSeed, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pUINT, CModelEntity::REACTIONS, CVector< CType >::resize(), CModel::setState(), CCopasiVector< T >::size(), and CModel::updateSimulatedValues().

 {
   /* get configuration data */
 
   bool useRandomSeed = * getValue("Use Random Seed").pBOOL;
   unsigned C_INT32 randomSeed = * getValue("Random Seed").pUINT;
 
   if (useRandomSeed) mpRandomGenerator->initialize(randomSeed);
 
   mEpsilon = * getValue("Epsilon").pDOUBLE;
   mUseRandomSeed = * getValue("Use Random Seed").pBOOL;
   mRandomSeed = * getValue("Random Seed").pUINT;
   mMaxSteps = * getValue("Max Internal Steps").pUINT;
 
   *mpCurrentState = *initialState;
 
   mpModel = mpProblem->getModel();
   assert(mpModel);
 
   if (mpModel->getModelType() == CModel::deterministic)
     mDoCorrection = true;
   else
     mDoCorrection = false;
 
   // Size the arrays
   mNumReactions = mpModel->getReactions().size();
 
   mReactionDependencies.resize(mNumReactions);
   mK.resize(mNumReactions);
   mAmu.resize(mNumReactions);
   mAmu = 0.0;
 
   mNumReactionSpecies = mpModel->getNumIndependentReactionMetabs() + mpModel->getNumDependentReactionMetabs();
 
   mAvgDX.resize(mNumReactionSpecies);
   mSigDX.resize(mNumReactionSpecies);
 
   // Create a local copy of the state where the particle number are rounded to integers
   mMethodState = mpModel->getState();
 
   const CStateTemplate & StateTemplate = mpModel->getStateTemplate();
 
   CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
   CModelEntity *const* endEntity = StateTemplate.endFixed();
   C_FLOAT64 * pValue = mMethodState.beginIndependent();
 
   mFirstReactionSpeciesIndex = 0;
   size_t Index = 1;
 
   for (; ppEntity != endEntity; ++ppEntity, ++pValue, ++Index)
     {
       if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
         {
           *pValue = floor(*pValue + 0.5);
 
           if (mFirstReactionSpeciesIndex == 0 &&
               (*ppEntity)->getStatus() == CModelEntity::REACTIONS)
             {
               mFirstReactionSpeciesIndex = Index;
             }
         }
     }
 
   // Update the model state so that the species are all represented by integers.
   mpModel->setState(mMethodState);
   mpModel->updateSimulatedValues(false); //for assignments
 
   C_FLOAT64 * pMethodStateValue = mMethodState.beginIndependent() - 1;
 
   // Build the reaction dependencies
   size_t NumReactions = 0;
 
   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->mMethodSpeciesIndex.resize(Balances.size());
       itDependencies->mSpeciesMultiplier.resize(Balances.size());
       itDependencies->mMethodSpecies.resize(Balances.size());
       itDependencies->mModelSpecies.resize(Balances.size());
 
       CCopasiVector< CChemEqElement >::const_iterator itBalance = Balances.begin();
       CCopasiVector< CChemEqElement >::const_iterator endBalance = Balances.end();
 
       Index = 0;
 
       for (; itBalance != endBalance; ++itBalance)
         {
           const CMetab * pMetab = (*itBalance)->getMetabolite();
 
           if (pMetab->getStatus() == CModelEntity::REACTIONS)
             {
               itDependencies->mMethodSpeciesIndex[Index] = StateTemplate.getIndex(pMetab) - mFirstReactionSpeciesIndex;
               itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
               itDependencies->mMethodSpecies[Index] = pMethodStateValue + StateTemplate.getIndex(pMetab);
               itDependencies->mModelSpecies[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
 
               Index++;
             }
         }
 
       // Correct allocation for metabolites which are not determined by reactions
       itDependencies->mMethodSpeciesIndex.resize(Index, true);
       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();
 
           itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
           itDependencies->mMethodSubstrates[Index] = pMethodStateValue + StateTemplate.getIndex(pMetab);
           itDependencies->mModelSubstrates[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
         }
 
       ++itDependencies;
       ++NumReactions;
     }
 
   mNumReactions = NumReactions;
 
   mReactionDependencies.resize(mNumReactions);
   mAmu.resize(mNumReactions, true);
   mK.resize(mNumReactions, true);
 
   return;
 }

CTrajectoryMethod::Status CTauLeapMethod::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 164 of file CTauLeapMethod.cpp.

References C_FLOAT64, doSingleStep(), CCopasiMessage::EXCEPTION, CCopasiProblem::getModel(), CModel::getState(), CState::getTime(), MCTrajectoryMethod, mMaxSteps, mMethodState, CTrajectoryMethod::mpCurrentState, mpModel, CTrajectoryMethod::mpProblem, CTrajectoryMethod::NORMAL, CModel::setState(), CState::setTime(), and CModel::updateSimulatedValues().

 {
   // do several steps
   C_FLOAT64 Time = mpCurrentState->getTime();
   C_FLOAT64 EndTime = Time + deltaT;
 
   size_t Steps = 0;
 
   while (Time < EndTime)
     {
       mMethodState.setTime(Time);
       mpModel->setState(mMethodState);
       mpModel->updateSimulatedValues(false);
 
       // 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.
 
       Time += doSingleStep(EndTime - Time);
 
       if (++Steps > mMaxSteps)
         {
           CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 12);
         }
     }
 
   *mpCurrentState = mpProblem->getModel()->getState();
   mpCurrentState->setTime(Time);
 
   return NORMAL;
 }

void CTauLeapMethod::updatePropensities ( )

protected

Calculate the propensities for all reactions

Definition at line 460 of file CTauLeapMethod.cpp.

References calculateAmu(), mA0, and mNumReactions.

Referenced by doSingleStep().

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

bool CTauLeapMethod::updateSystem ( )

protected

Updates the system according to the probabilistic number of firings mK[i] of each reaction i

Definition at line 546 of file CTauLeapMethod.cpp.

References CVectorCore< CType >::array(), CState::beginIndependent(), C_FLOAT64, mFirstReactionSpeciesIndex, mK, mMethodState, mNumReactions, mNumReactionSpecies, and mReactionDependencies.

Referenced by doSingleStep().

 {
   std::vector< CReactionDependencies >::const_iterator itReaction = mReactionDependencies.begin();
 
   CState OldState(mMethodState);
 
   const C_FLOAT64 * pK = mK.array();
   const C_FLOAT64 * pKEnd = pK + mNumReactions;
 
   for (; pK != pKEnd; ++pK, ++itReaction)
     {
       const C_FLOAT64 * pMultiplicity = itReaction->mSpeciesMultiplier.array();
       const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + itReaction->mSpeciesMultiplier.size();
       C_FLOAT64 * const * pSpecies = itReaction->mMethodSpecies.array();
 
       for (; pMultiplicity != pMultiplicityEnd; ++pMultiplicity, ++pSpecies)
         {
           **pSpecies += *pK * *pMultiplicity;
         }
     }
 
   const C_FLOAT64 * pSpecies = mMethodState.beginIndependent() + mFirstReactionSpeciesIndex - 1;
 
   const C_FLOAT64 * pSpeciesEnd = pSpecies + mNumReactionSpecies;
 
   for (; pSpecies != pSpeciesEnd; ++pSpecies)
     {
       if (*pSpecies < -0.5)
         {
           // We need to undo the changes
           mMethodState = OldState;
           return false;
         }
     }
 
   return true;
 }

Friends And Related Function Documentation

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

friend

Member Data Documentation

C_FLOAT64 CTauLeapMethod::mA0

protected

Total propensity (sum over mAmu[i])

Definition at line 252 of file CTauLeapMethod.h.

Referenced by updatePropensities().

CVector< C_FLOAT64 > CTauLeapMethod::mAmu

protected

A vector of reaction propensities

Definition at line 247 of file CTauLeapMethod.h.

Referenced by calculateAmu(), doSingleStep(), and start().

CVector< C_FLOAT64> CTauLeapMethod::mAvgDX

protected

For tau-selection method

Definition at line 263 of file CTauLeapMethod.h.

Referenced by doSingleStep(), and start().

bool CTauLeapMethod::mDoCorrection

protected

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

Definition at line 299 of file CTauLeapMethod.h.

Referenced by calculateAmu(), and start().

C_FLOAT64 CTauLeapMethod::mEpsilon

protected

The tolerance ratio x(t+t')< eps*x(t)

Definition at line 273 of file CTauLeapMethod.h.

Referenced by doSingleStep(), and start().

size_t CTauLeapMethod::mFirstReactionSpeciesIndex

protected

index of first species in a CState

Definition at line 304 of file CTauLeapMethod.h.

Referenced by doSingleStep(), start(), and updateSystem().

CVector< C_FLOAT64 > CTauLeapMethod::mK

protected

The k-values of the reactions, that is the probabilistic number of firings within one leap.

Definition at line 258 of file CTauLeapMethod.h.

Referenced by doSingleStep(), start(), and updateSystem().

unsigned C_INT32 CTauLeapMethod::mMaxSteps

protected

The maximum number of tau leap steps allowed for an integrations step

Definition at line 278 of file CTauLeapMethod.h.

Referenced by start(), and step().

CState CTauLeapMethod::mMethodState

protected

The method internal state which contains particle rounded particle numbers.

Definition at line 227 of file CTauLeapMethod.h.

Referenced by doSingleStep(), start(), step(), and updateSystem().

size_t CTauLeapMethod::mNumReactions

protected

Number of reactions.

Definition at line 232 of file CTauLeapMethod.h.

Referenced by doSingleStep(), start(), updatePropensities(), and updateSystem().

size_t CTauLeapMethod::mNumReactionSpecies

protected

Number of variable metabolites.

Definition at line 242 of file CTauLeapMethod.h.

Referenced by doSingleStep(), start(), and updateSystem().

CModel* CTauLeapMethod::mpModel

protected

Pointer to the model.

Definition at line 222 of file CTauLeapMethod.h.

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

CRandom* CTauLeapMethod::mpRandomGenerator

protected

The random number generator.

Definition at line 294 of file CTauLeapMethod.h.

Referenced by cleanup(), CTauLeapMethod(), doSingleStep(), and start().

unsigned C_INT32 CTauLeapMethod::mRandomSeed

protected

The random seed to use.

Definition at line 289 of file CTauLeapMethod.h.

Referenced by start().

std::vector< CReactionDependencies > CTauLeapMethod::mReactionDependencies

protected

A vector containing dependency information to minimize the required updates.

Definition at line 237 of file CTauLeapMethod.h.

Referenced by calculateAmu(), doSingleStep(), start(), and updateSystem().

CVector< C_FLOAT64> CTauLeapMethod::mSigDX

protected

For tau-selection method

Definition at line 268 of file CTauLeapMethod.h.

Referenced by doSingleStep(), and start().

bool CTauLeapMethod::mUseRandomSeed

protected

Specifies if the mRandomSeed should be used. otherwise a randomly chosen seed is used.

Definition at line 284 of file CTauLeapMethod.h.

Referenced by start().

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

Classes

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation