#include <CHybridMethodLSODA.h>

Inheritance diagram for CHybridMethodLSODA:

Collaboration diagram for CHybridMethodLSODA:

Classes
struct	Data

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

	~CHybridMethodLSODA ()

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

Static Public Member Functions
static CHybridMethodLSODA *	createHybridMethodLSODA ()

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)

Protected Types
enum	metabStatus { LOW = 0, HIGH }

Protected Types inherited from CCopasiObject
enum	Flag { Container = 0x1, Vector = 0x2, Matrix = 0x4, NameVector = 0x8, Reference = 0x10, ValueBool = 0x20, ValueInt = 0x40, ValueInt64 = 0x80, ValueDbl = 0x100, NonUniqueName = 0x200, StaticString = 0x400, ValueString = 0x800, Separator = 0x1000, ModelEntity = 0x2000, Array = 0x4000, DataModel = 0x8000, Root = 0x10000, Gui = 0x20000 }

Protected Member Functions
void	calculateAmu (size_t rIndex)

void	calculateDerivative (std::vector< C_FLOAT64 > &deriv)

	CHybridMethodLSODA (const CCopasiContainer *pParent=NULL)

void	cleanup ()

virtual C_FLOAT64	doSingleStep (C_FLOAT64 currentTime, C_FLOAT64 endTime)=0

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

void	fireReaction (size_t rIndex)

C_FLOAT64	generateReactionTime (size_t rIndex)

std::set< std::string > *	getAffects (size_t rIndex)

C_FLOAT64	getDefaultAtol (const CModel *pModel) const

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

std::set< size_t > *	getParticipatesIn (size_t rIndex)

void	getState (std::vector< C_FLOAT64 > &target)

void	getStochTimeAndIndex (C_FLOAT64 &ds, size_t &rIndex)

void	initMethod (C_FLOAT64 time)

void	insertDeterministicReaction (size_t rIndex)

void	integrateDeterministicPart (C_FLOAT64 ds)

void	outputData (std::ostream &os, C_INT32 mode)

void	outputDebug (std::ostream &os, size_t level)

void	partitionSystem ()

void	removeDeterministicReaction (size_t rIndex)

void	setState (std::vector< C_FLOAT64 > &source)

void	setupBalances ()

void	setupDependencyGraph ()

void	setupMetab2React ()

void	setupMetab2ReactComplete ()

void	setupMetab2ReactPlusModifier ()

void	setupPartition ()

void	setupPriorityQueue (C_FLOAT64 startTime=0.0)

void	updatePriorityQueue (size_t rIndex, C_FLOAT64 time)

void	updateTauMu (size_t rIndex, 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 C_INT32	checkModel (CModel *model)

static void	EvalF (const C_INT n, const C_FLOAT64 t, const C_FLOAT64 y, C_FLOAT64 ydot)

static bool	modelHasAssignments (const CModel *pModel)

Protected Attributes
std::vector< C_FLOAT64 >	currentState

std::vector< C_FLOAT64 >	mAmu

std::vector< C_FLOAT64 >	mAmuOld

C_FLOAT64	mAtol

Data	mData

bool	mDefaultAtol

CDependencyGraph	mDG

bool	mDoCorrection

CVector< C_FLOAT64 >	mDWork

C_FLOAT64	mEndt

std::ostringstream	mErrorMsg

size_t	mFirstMetabIndex

CHybridLSODAStochFlag *	mFirstReactionFlag

bool	mHasAssignments

CVector< C_INT >	mIWork

C_INT	mJType

std::vector< std::vector < CHybridLSODABalance > >	mLocalBalances

std::vector< std::vector < CHybridLSODABalance > >	mLocalSubstrates

C_FLOAT64	mLowerStochLimit

CLSODA	mLSODA

C_INT	mLsodaStatus

C_INT32	mMaxBalance

C_INT32	mMaxIntBeforeStep

unsigned C_INT32	mMaxSteps

bool	mMaxStepsReached

std::vector< std::set< size_t > >	mMetab2React

std::vector< metabStatus >	mMetabFlags

size_t	mNumVariableMetabs

size_t	mOutputCounter

std::ofstream	mOutputFile

std::string	mOutputFileName

unsigned C_INT32	mPartitioningInterval

CCopasiVector< CMetab > *	mpMetabolites

CModel *	mpModel

CIndexedPriorityQueue	mPQ

CRandom *	mpRandomGenerator

const CCopasiVectorNS < CReaction > *	mpReactions

CState *	mpState

unsigned C_INT32	mRandomSeed

std::vector < CHybridLSODAStochFlag >	mReactionFlags

bool	mReducedModel

bool	mRestartLSODA

C_FLOAT64	mRtol

C_INT	mState

unsigned C_INT32	mStepsAfterPartitionSystem

C_FLOAT64	mStepsize

CMatrix< C_FLOAT64 >	mStoi

C_FLOAT64	mTime

std::set< size_t >	mUpdateSet

C_FLOAT64	mUpperStochLimit

bool	mUseRandomSeed

C_FLOAT64 *	mY

CVector< C_FLOAT64 >	mYdot

std::vector< C_FLOAT64 >	temp

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

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

Detailed Description

Definition at line 118 of file CHybridMethodLSODA.h.

Member Enumeration Documentation

enum CHybridMethodLSODA::metabStatus

protected

Status of the metabolites. Particle number can be high or low.

Enumerator
LOW
HIGH

Definition at line 476 of file CHybridMethodLSODA.h.

476 {LOW = 0, HIGH};

CHybridMethodLSODA::HIGH

Definition: CHybridMethodLSODA.h:476

CHybridMethodLSODA::LOW

Definition: CHybridMethodLSODA.h:476

Constructor & Destructor Documentation

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

Copy constructor

Parameters

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

Definition at line 74 of file CHybridMethodLSODA.cpp.

References CRandom::createGenerator(), CHybridMethodLSODA::Data::dim, initializeParameter(), mData, mpRandomGenerator, CRandom::mt19937, and CHybridMethodLSODA::Data::pMethod.

                                                                         :
   CTrajectoryMethod(src, pParent)
 {
   assert((void *) &mData == (void *) &mData.dim);
 
   mData.pMethod = this;
 
   mpRandomGenerator = CRandom::createGenerator(CRandom::mt19937);
   initializeParameter();
 }

CHybridMethodLSODA::~CHybridMethodLSODA ( )

Destructor.

Definition at line 89 of file CHybridMethodLSODA.cpp.

References cleanup(), and DESTRUCTOR_TRACE.

 {
   cleanup();
   DESTRUCTOR_TRACE;
 }

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

protected

Default constructor.

Parameters

const CCopasiContainer * pParent (default: NULL)

CHybridMethodLSODA

This class implements an hybrid algorithm for the simulation of a biochemical system over time.

File name: CHybridMethodLSODA.cpp Author: Juergen Pahle Email: juerg.nosp@m.en.p.nosp@m.ahle@.nosp@m.eml-.nosp@m.r.vil.nosp@m.la-b.nosp@m.osch..nosp@m.de

Last change: 09, May 2006

(C) European Media Lab 2003. Default constructor.

Definition at line 63 of file CHybridMethodLSODA.cpp.

References CRandom::createGenerator(), CHybridMethodLSODA::Data::dim, initializeParameter(), mData, mpRandomGenerator, CRandom::mt19937, and CHybridMethodLSODA::Data::pMethod.

                                                                       :
   CTrajectoryMethod(CCopasiMethod::hybridLSODA, pParent)
 {
   assert((void *) &mData == (void *) &mData.dim);
 
   mData.pMethod = this;
 
   mpRandomGenerator = CRandom::createGenerator(CRandom::mt19937);
   initializeParameter();
 }

Member Function Documentation

void CHybridMethodLSODA::calculateAmu ( size_t rIndex )

protected

Calculates an amu value for a given reaction.

Parameters

rIndex A size_t specifying the reaction to be updated

Definition at line 706 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, C_INT32, mAmu, mDoCorrection, and mLocalSubstrates.

Referenced by partitionSystem(), setupPriorityQueue(), and updatePriorityQueue().

 {
   if (!mDoCorrection)
     {
       mAmu[rIndex] = (*mpReactions)[rIndex]->calculateParticleFlux();
       return;
     }
 
   // 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_INT32 number = 0;
   C_INT32 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<CHybridLSODABalance> & substrates = mLocalSubstrates[rIndex];
 
   int flag = 0;
 
   for (size_t i = 0; i < substrates.size(); i++)
     {
       num_ident = substrates[i].mMultiplicity;
 
       if (num_ident > 1)
         {
           flag = 1;
           number = static_cast<C_INT32>((*mpMetabolites)[substrates[i].mIndex]->getValue());
           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[rIndex] = 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 = (*mpReactions)[rIndex]->calculateParticleFlux();
 
   if (flag)
     {
       amu *= rate_factor / substrate_factor;;
       mAmu[rIndex] = amu;
     }
   else
     {
       mAmu[rIndex] = rate_factor;
     }
 
   return;
 
   // a more efficient way to calculate mass action kinetics could be included
 }

void CHybridMethodLSODA::calculateDerivative ( std::vector< C_FLOAT64 > & deriv )

protected

Calculates the derivative of the system and writes it into the vector deriv. Length of deriv must be mNumVariableMetabs. CAUTION: Only deterministic reactions are taken into account. That is, this is only the derivative of the deterministic part of the system.

Parameters

deriv A vector reference of length mNumVariableMetabs, into which the derivative is written

Definition at line 542 of file CHybridMethodLSODA.cpp.

References C_INT32, mFirstReactionFlag, CHybridLSODAStochFlag::mIndex, CHybridLSODAStochFlag::mpNext, mStoi, and CMatrix< CType >::numRows().

Referenced by evalF().

 {
   size_t i;
   C_INT32 bal = 0;
   CHybridLSODAStochFlag * j;
 
   // Calculate all the needed kinetic functions of the deterministic reactions
   for (j = mFirstReactionFlag; j != NULL; j = j->mpNext)
     {
       (*mpReactions)[j->mIndex]->calculateParticleFlux();
     }
 
   // For each metabolite add up the contributions of the det. reactions
   // the number of rows in mStoi equals the number of non-fixed metabolites!
   //  for (i=0; i<mNumVariableMetabs; i++)
   for (i = 0; i < mStoi.numRows(); i++)
     {
       deriv[i] = 0.0;
 
       for (j = mFirstReactionFlag; j != NULL; j = j->mpNext)
         {
           // juergen: +0.5 to get a rounding out of the static_cast
           bal = static_cast<C_INT32>(floor(mStoi[i][j->mIndex] + 0.5));
           deriv[i] += bal * (*mpReactions)[j->mIndex]->getParticleFlux(); //  balance * flux;
         }
     }
 
   /*
     for (; i < mNumVariableMetabs; i++) deriv[i] = 0.0; // important to get a correct deriv vector, because mStoi doesn't cover fixed metabolites
   */
   return;
 }

C_INT32 CHybridMethodLSODA::checkModel ( CModel * model )

staticprotected

Test the model if it is proper to perform stochastic simulations on. Several properties are tested (e.g. integer stoichometry, all reactions take place in one compartment only, irreversibility...).

Parameters

model The model to be checked

Returns: 1, if hybrid simulation is possible; <0, if an error occured.

Test the model if it is proper to perform stochastic simulations on. Several properties are tested (e.g. integer stoichometry, all reactions take place in one compartment only, irreversibility...).

Returns: 0, if everything is ok; <0, if an error occured.

Definition at line 819 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, C_INT32, CModel::getReactions(), CModel::getStoi(), INT_EPSILON, mpReactions, mStoi, CMatrix< CType >::numRows(), and CCopasiVector< T >::size().

 {
   CCopasiVectorNS <CReaction> * mpReactions = &model->getReactions();
   CMatrix <C_FLOAT64> mStoi = model->getStoi();
   C_INT32 multInt;
   size_t i, j, numReactions = mpReactions->size();
   C_FLOAT64 multFloat;
   //  size_t metabSize = mpMetabolites->size();
 
   for (i = 0; i < numReactions; i++) // for every reaction
     {
       // TEST getCompartmentNumber() == 1
       if ((*mpReactions)[i]->getCompartmentNumber() != 1) return - 1;
 
       // TEST isReversible() == 0
       if ((*mpReactions)[i]->isReversible() != 0) return - 2;
 
       // TEST integer stoichometry
       // Iterate through each the metabolites
       // juergen: the number of rows of mStoi equals the number of non-fixed metabs!
       //  for (j=0; i<metabSize; j++)
       for (j = 0; j < mStoi.numRows(); j++)
         {
           multFloat = mStoi[j][i];
           multInt = static_cast<C_INT32>(floor(multFloat + 0.5)); // +0.5 to get a rounding out of the static_cast to int!
 
           if ((multFloat - multInt) > INT_EPSILON) return - 3; // INT_EPSILON in CHybridMethod.h
         }
     }
 
   return 1; // Model is appropriate for hybrid simulation
 }

void CHybridMethodLSODA::cleanup ( )

protected

Cleans up memory, etc.

Definition at line 382 of file CHybridMethodLSODA.cpp.

References mpModel, and mpRandomGenerator.

Referenced by ~CHybridMethodLSODA().

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

CHybridMethodLSODA * CHybridMethodLSODA::createHybridMethodLSODA ( )

static

Creates a HybridMethod adequate for the problem. (only CHybridNextReactionLSODAMethod so far)

Creates a HybridMethodLSODA adequate for the problem.

Definition at line 175 of file CHybridMethodLSODA.cpp.

References C_INT32.

Referenced by CTrajectoryMethod::createMethod().

 {
   C_INT32 result = 1;
   CHybridMethodLSODA * method = NULL;
 
   switch (result)
     {
         /*    case - 3:        // non-integer stoichometry
         CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 1);
         break;
         case - 2:        // reversible reaction exists
         CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 2);
         break;
 
         case - 1:        // more than one compartment involved
         CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 3);
         break;*/
       case 1:
       default:
         // Everything alright: Hybrid simulation possible
         method = new CHybridNextReactionLSODAMethod();
         break;
     }
 
   return method;
 }

virtual C_FLOAT64 CHybridMethodLSODA::doSingleStep	(	C_FLOAT64	currentTime,
		C_FLOAT64	endTime
	)

protectedpure 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

Implemented in CHybridNextReactionLSODAMethod.

Referenced by step().

bool CHybridMethodLSODA::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 166 of file CHybridMethodLSODA.cpp.

References initializeParameter().

 {
   initializeParameter();
   return true;
 }

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

staticprotected

Definition at line 392 of file CHybridMethodLSODA.cpp.

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

Referenced by integrateDeterministicPart().

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

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

protected

This evaluates the derivatives for the complete model

Definition at line 395 of file CHybridMethodLSODA.cpp.

References calculateDerivative(), CCopasiProblem::getModel(), mNumVariableMetabs, CTrajectoryMethod::mpProblem, mpState, CModel::setState(), CState::setTime(), temp, and CModel::updateSimulatedValues().

Referenced by EvalF().

 {
   size_t i;
 
   CModel * pModel = mpProblem->getModel();
 
   mpState->setTime(*t);
 
   pModel->setState(*mpState);
 
   pModel->updateSimulatedValues(false);
 
   calculateDerivative(temp);
 
   for (i = 0; i < mNumVariableMetabs; i++)
     ydot[i] = temp[i];
 
   return;
 }

void CHybridMethodLSODA::fireReaction ( size_t rIndex )

protected

Executes the specified reaction in the system once.

Parameters

rIndex A size_t specifying the index of the reaction, which will be fired.

Executes the specified reaction in the system once.

Parameters

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

Definition at line 600 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, CDependencyGraph::getDependents(), CModelEntity::getValue(), mDG, mLocalBalances, mRestartLSODA, mUpdateSet, CMetab::refreshConcentration(), and CModelEntity::setValue().

Referenced by CHybridNextReactionLSODAMethod::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.
 
   size_t i;
   C_FLOAT64 newNumber;
 
   CMetab * pMetab;
 
   for (i = 0; i < mLocalBalances[rIndex].size(); i++)
     {
       pMetab = mLocalBalances[rIndex][i].mpMetabolite;
       newNumber = pMetab->getValue() + mLocalBalances[rIndex][i].mMultiplicity;
 
       pMetab->setValue(newNumber);
       pMetab->refreshConcentration();
     }
 
   // insert all dependent reactions into the mUpdateSet
   const std::set <size_t> & dependents = mDG.getDependents(rIndex);
   std::copy(dependents.begin(), dependents.end(),
             std::inserter(mUpdateSet, mUpdateSet.begin()));
 
   mRestartLSODA = true;
 
   return;
 }

C_FLOAT64 CHybridMethodLSODA::generateReactionTime ( size_t rIndex )

protected

Generates a putative reaction time for the given reaction.

Parameters

rIndex A size_t specifying the index of the reaction

Returns: A C_FLOAT64 holding the calculated reaction time

Definition at line 693 of file CHybridMethodLSODA.cpp.

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

Referenced by partitionSystem(), setupPriorityQueue(), updatePriorityQueue(), and updateTauMu().

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

std::set< std::string > * CHybridMethodLSODA::getAffects ( size_t rIndex )

protected

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

Parameters

rIndex The index of the reaction being executed.

Returns: The set of affected metabolites.

Definition at line 1349 of file CHybridMethodLSODA.cpp.

References mLocalBalances.

Referenced by setupDependencyGraph().

 {
   size_t i;
   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!
 
   for (i = 0; i < mLocalBalances[rIndex].size(); i++)
     {
       if (mLocalBalances[rIndex][i].mMultiplicity != 0)
         {
           retset->insert(mLocalBalances[rIndex][i].mpMetabolite->getKey());
         }
     }
 
   return retset;
 }

C_FLOAT64 CHybridMethodLSODA::getDefaultAtol ( const CModel * pModel ) const

protected

Calculate the default absolute tolerance

Parameters

const CModel * pModel

Returns: C_FLOAT64 defaultAtol

Definition at line 361 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, CModel::getCompartments(), CModel::getQuantity2NumberFactor(), CCopasiParameter::getValue(), max, and CCopasiVector< T >::size().

Referenced by initMethod().

 {
 
   if (!pModel) return 1.0e009;
 
   const CCopasiVectorNS< CCompartment > & Compartment = pModel->getCompartments();
   size_t i, imax;
 
   C_FLOAT64 Volume = std::numeric_limits< C_FLOAT64 >::max();
 
   for (i = 0, imax = Compartment.size(); i < imax; i++)
     if (Compartment[i]->getValue() < Volume) Volume = Compartment[i]->getValue();
 
   if (Volume == std::numeric_limits< C_FLOAT64 >::max()) return 1.0e009;
 
   return Volume * pModel->getQuantity2NumberFactor() * 1.e-12;
 }

std::set< std::string > * CHybridMethodLSODA::getDependsOn ( size_t rIndex )

protected

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

Parameters

rIndex The index of the reaction being executed.

Returns: The set of metabolites depended on.

Definition at line 1316 of file CHybridMethodLSODA.cpp.

References mpReactions, and CFunctionParameter::PARAMETER.

Referenced by setupDependencyGraph().

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

std::set< size_t > * CHybridMethodLSODA::getParticipatesIn ( size_t rIndex )

protected

Gets the set of metabolites, which participate in the given reaction either as substrate, product or modifier.

Parameters

rIndex The index of the reaction being executed.

Returns: The set of participating metabolites.

Definition at line 1375 of file CHybridMethodLSODA.cpp.

Referenced by setupMetab2ReactPlusModifier().

 {
   std::set<size_t> *retset = new std::set<size_t>;
   return retset;
 }

void CHybridMethodLSODA::getState ( std::vector< C_FLOAT64 > & target )

protected

Gathers the state of the system into the array target. Later on CState should be used for this. Length of the array target must be mNumVariableMetabs.

Parameters

target A vector reference of length mNumVariableMetabs, into which the state of the system is written

void CHybridMethodLSODA::getStochTimeAndIndex	(	C_FLOAT64 &	ds,
		size_t &	rIndex
	)

protected

Find the reaction index and the reaction time of the stochastic (!) reaction with the lowest reaction time.

Parameters

ds	A reference to a C_FLOAT64. The putative reaction time for the first stochastic reaction is written into this variable.
rIndex	A reference to a size_t. The index of the first stochastic reaction is written into this variable.

Definition at line 586 of file CHybridMethodLSODA.cpp.

References mPQ, CIndexedPriorityQueue::topIndex(), and CIndexedPriorityQueue::topKey().

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

 {
   ds = mPQ.topKey();
   rIndex = mPQ.topIndex();
   return;
 }

void CHybridMethodLSODA::initializeParameter ( )

private

Intialize the method parameter

Definition at line 95 of file CHybridMethodLSODA.cpp.

Referenced by CHybridMethodLSODA(), and elevateChildren().

 {
   CCopasiParameter *pParm;
 
   assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) MAX_STEPS);
   assertParameter("Lower Limit", CCopasiParameter::DOUBLE, (C_FLOAT64) LOWER_STOCH_LIMIT);
   assertParameter("Upper Limit", CCopasiParameter::DOUBLE, (C_FLOAT64) UPPER_STOCH_LIMIT);
 
   assertParameter("Partitioning Interval", CCopasiParameter::UINT, (unsigned C_INT32) PARTITIONING_INTERVAL);
   assertParameter("Use Random Seed", CCopasiParameter::BOOL, (bool) USE_RANDOM_SEED);
   assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) RANDOM_SEED);
 
   // 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)
         {
           setValue("Use Random Seed", *pParm->getValue().pBOOL);
           removeParameter("UseRandomSeed");
         }
 
       if ((pParm = getParameter("")) != NULL)
         {
           setValue("Random Seed", *pParm->getValue().pUINT);
           removeParameter("");
         }
     }
 
   /* LSODA ****************************************************************************/
 
   addParameter("Partitioning Stepsize",
                CCopasiParameter::DOUBLE, (C_FLOAT64) PARTITIONING_STEPSIZE);
 
   addParameter("Relative Tolerance",
                CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-006);
 
   addParameter("Use Default Absolute Tolerance",
                CCopasiParameter::BOOL, (bool) true);
 
   addParameter("Absolute Tolerance",
                CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e009);
 
   addParameter("Max Internal Steps (LSODA)",
                CCopasiParameter::UINT, (unsigned C_INT32) 10000);
 
   // These parameters are no longer supported.
   removeParameter("Adams Max Order");
   removeParameter("BDF Max Order");
 }

void CHybridMethodLSODA::initMethod ( C_FLOAT64 start_time )

protected

Initializes the solver.

Parameters

time	the current time

Initializes the solver and sets the model to be used.

Parameters

model A reference to an instance of a CModel

Definition at line 290 of file CHybridMethodLSODA.cpp.

References C_INT, CHybridMethodLSODA::Data::dim, getDefaultAtol(), CModel::getMetabolitesX(), CCopasiProblem::getModel(), CModel::getNumDependentReactionMetabs(), CModel::getNumIndependentReactionMetabs(), CModel::getReactions(), CModel::getStoi(), CCopasiParameter::getValue(), CRandom::initialize(), mAmu, mAmuOld, mAtol, mData, mDefaultAtol, mDWork, mIWork, mLowerStochLimit, mMaxSteps, mMaxStepsReached, mNumVariableMetabs, mPartitioningInterval, mpMetabolites, mpModel, CTrajectoryMethod::mpProblem, mpRandomGenerator, mpReactions, mRandomSeed, mRtol, mStepsAfterPartitionSystem, mStepsize, mStoi, mUpdateSet, mUpperStochLimit, mUseRandomSeed, mYdot, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pUDOUBLE, CCopasiParameter::Value::pUINT, CVector< CType >::resize(), setupBalances(), setupDependencyGraph(), setupMetab2React(), setupPartition(), setupPriorityQueue(), CCopasiParameterGroup::setValue(), CCopasiVector< T >::size(), and temp.

Referenced by start().

 {
   mpReactions = &mpModel->getReactions();
   mAmu.clear();
   mAmu.resize(mpReactions->size());
   mAmuOld.clear();
   mAmuOld.resize(mpReactions->size());
   mpMetabolites = &(const_cast < CCopasiVector < CMetab > & >(mpModel->getMetabolitesX()));
   //mNumVariableMetabs = mpModel->getNumVariableMetabs(); // ind + dep metabs, without fixed metabs
   mNumVariableMetabs = mpModel->getNumIndependentReactionMetabs() + mpModel->getNumDependentReactionMetabs();
 
   temp.clear();
   temp.resize(mNumVariableMetabs);
   //currentState.clear();
   //currentState.resize(mNumVariableMetabs);
 
   /* get configuration data */
   mMaxSteps = * getValue("Max Internal Steps").pUINT;
   mLowerStochLimit = * getValue("Lower Limit").pDOUBLE;
   mUpperStochLimit = * getValue("Upper Limit").pDOUBLE;
   mStepsize = * getValue("Partitioning Stepsize").pDOUBLE;
   mPartitioningInterval = * getValue("Partitioning Interval").pUINT;
   mUseRandomSeed = * getValue("Use Random Seed").pBOOL;
   mRandomSeed = * getValue("Random Seed").pUINT;
 
   if (mUseRandomSeed) mpRandomGenerator->initialize(mRandomSeed);
 
   mStoi = mpModel->getStoi();
   mStepsAfterPartitionSystem = 0;
   mUpdateSet.clear();
 
   setupBalances(); // initialize mLocalBalances and mLocalSubstrates (has to be called first!)
   setupDependencyGraph(); // initialize mDG
   setupMetab2React(); // initialize mMetab2React
   setupPartition(); // initialize mReactionFlags
   setupPriorityQueue(start_time); // initialize mPQ
 
   mMaxStepsReached = false;
 
   /* CONFIGURE LSODA ***********************************************************************/
 
   mData.dim = (C_INT) mNumVariableMetabs;
 
   mYdot.resize(mData.dim);
 
   mRtol = * getValue("Relative Tolerance").pUDOUBLE;
 
   mDefaultAtol = * getValue("Use Default Absolute Tolerance").pBOOL;
 
   if (mDefaultAtol)
     {
       mAtol = getDefaultAtol(mpProblem->getModel());
       setValue("Absolute Tolerance", mAtol);
     }
   else
     mAtol = * getValue("Absolute Tolerance").pUDOUBLE;
 
   mDWork.resize(22 + mData.dim * std::max<C_INT>(16, mData.dim + 9));
   mDWork[4] = mDWork[5] = mDWork[6] = mDWork[7] = mDWork[8] = mDWork[9] = 0.0;
   mIWork.resize(20 + mData.dim);
 
   mIWork[4] = mIWork[6] = mIWork[9] = 0;
   mIWork[5] = * getValue("Max Internal Steps (LSODA)").pUINT;
   mIWork[7] = 12;
   mIWork[8] = 5;
 
   /***********************************************************************************/
 
   return;
 }

void CHybridMethodLSODA::insertDeterministicReaction ( size_t rIndex )

protected

Inserts a new deterministic reaction into the linked list in the vector mReactionFlags.

Parameters

rIndex A size_t giving the index of the reaction to be inserted into the list of deterministic reactions.

Inserts a new deterministic reaction into the linked list in the array mReactionFlags.

Parameters

rIndex A size_t giving the index of the reaction to be inserted into the list of deterministic reactions.

Definition at line 1242 of file CHybridMethodLSODA.cpp.

References mAmu, mAmuOld, mFirstReactionFlag, CHybridLSODAStochFlag::mpPrev, and mReactionFlags.

Referenced by partitionSystem().

 {
   if (mReactionFlags[rIndex].mpPrev == NULL)
     // reaction is stochastic (avoids double insertions)
     {
       if (mFirstReactionFlag != NULL)
         // there are deterministic reactions already
         {
           mFirstReactionFlag->mpPrev = &mReactionFlags[rIndex];
           mReactionFlags[rIndex].mpNext = mFirstReactionFlag;
           mFirstReactionFlag = &mReactionFlags[rIndex];
           mFirstReactionFlag->mpPrev = mFirstReactionFlag;
         }
       else
         {
           // there are no deterministic reactions
           // Important to distinguish between stochastic (prev == NULL) and deterministic (prev != NULL) reactions
           mReactionFlags[rIndex].mpPrev = &mReactionFlags[rIndex];
           mFirstReactionFlag = &mReactionFlags[rIndex];
         }
 
       mAmu[rIndex] = 0.0;
       mAmuOld[rIndex] = 0.0;
     }
 
   return;
 }

void CHybridMethodLSODA::integrateDeterministicPart ( C_FLOAT64 deltaT )

protected

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

Parameters

ds	A C_FLOAT64 specifying the stepsize.

Definition at line 422 of file CHybridMethodLSODA.cpp.

References CVectorCore< CType >::array(), CState::beginIndependent(), C_FLOAT64, C_INT, CHybridMethodLSODA::Data::dim, EvalF(), CCopasiMessage::EXCEPTION, CDependencyGraph::getDependents(), CCopasiProblem::getModel(), CModel::getState(), CState::getTime(), mAtol, max, MCTrajectoryMethod, mData, mDG, mDWork, mErrorMsg, mFirstReactionFlag, CHybridLSODAStochFlag::mIndex, mIWork, mJType, mLSODA, mLsodaStatus, mpModel, CHybridLSODAStochFlag::mpNext, CTrajectoryMethod::mpProblem, mpState, mRestartLSODA, mRtol, mState, mTime, mUpdateSet, mY, CInternalSolver::setOstream(), CModel::setState(), CState::setTime(), CModel::setTime(), CVectorCore< CType >::size(), and CModel::updateSimulatedValues().

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

 {
   CHybridLSODAStochFlag * react = NULL;
 
   /* RESET LSODA **************************************************************/
 
   if (mRestartLSODA)
     {
       mRestartLSODA = false;
       mLsodaStatus = 1;
     }
 
   mState = 1;
   mJType = 2;
 
   mErrorMsg.str("");
   mLSODA.setOstream(mErrorMsg);
 
   C_INT one = 1;
   C_INT DSize = (C_INT) mDWork.size();
   C_INT ISize = (C_INT) mIWork.size();
 
   //mpState = new CState (mpProblem->getModel()->getState());
   *mpState = mpProblem->getModel()->getState();
 
   mY = mpState->beginIndependent();// + mFirstMetabIndex - 1;
   //for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, mY++)
   //  *mY = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
 
   mTime = mpState->getTime();
 
   C_FLOAT64 EndTime = mTime + deltaT;
 
   C_FLOAT64 tdist , d__1, d__2, w0;
   tdist = fabs(deltaT);
   d__1 = fabs(mTime), d__2 = fabs(EndTime);
   w0 = std::max(d__1, d__2);
 
   if (tdist < std::numeric_limits< C_FLOAT64 >::epsilon() * 2. * w0) //just do nothing
     {
       //mTime = mTime + deltaT;
       //mpState->setTime(mTime);
       //*mpCurrentState = *mpState;
       mpProblem->getModel()->setTime(EndTime);
 
       return;
     }
 
   if (!mData.dim) //just do nothing if there are no variables
     {
       //mTime = mTime + deltaT;
       mpProblem->getModel()->setTime(EndTime);
       //mpState->setTime(mTime);
       //*mpCurrentState = *mpState;
 
       return;
     }
 
   mLSODA(&EvalF , //  1. evaluate F
          &mData.dim, //  2. number of variables
          mY , //  3. the array of current concentrations
          &mTime , //  4. the current time
          &EndTime , //  5. the final time
          &one , //  6. scalar error control
          &mRtol , //  7. relative tolerance array
          &mAtol , //  8. absolute tolerance array
          &mState , //  9. output by overshoot & interpolatation
          &mLsodaStatus , // 10. the state control variable
          &one , // 11. futher options (one)
          mDWork.array() , // 12. the double work array
          &DSize , // 13. the double work array size
          mIWork.array() , // 14. the int work array
          &ISize , // 15. the int work array size
          NULL , // 16. evaluate J (not given)
          &mJType);        // 17. the type of jacobian calculate (2)
 
   if (mLsodaStatus == -1) mLsodaStatus = 2;
 
   if ((mLsodaStatus != 1) && (mLsodaStatus != 2) && (mLsodaStatus != -1))
     {
       CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 6, mErrorMsg.str().c_str());
     }
 
   mpState->setTime(mTime);
   //*mpCurrentState = *mpState;
   mpProblem->getModel()->setState(*mpState);
 
   /* size_t i;
 
   for (i = 0; i < mNumVariableMetabs; i++)
     {
       (*mpMetabolites)[i]->setValue(mY[i]);
       (*mpMetabolites)[i]->refreshConcentration();
     } */
 
   mpModel->updateSimulatedValues(false); //for assignments
 
   //pdelete(mpState);
 
   // find the set union of all reactions, which depend on one of the deterministic reactions. The propensities of the stochastic reactions in this set union will be updated later in the method updatePriorityQueue().
 
   for (react = mFirstReactionFlag; react != NULL; react = react->mpNext)
     {
       const std::set <size_t> & dependents = mDG.getDependents(react->mIndex);
       std::copy(dependents.begin(), dependents.end(),
                 std::inserter(mUpdateSet, mUpdateSet.begin()));
     }
 
   return;
 }

bool CHybridMethodLSODA::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 1674 of file CHybridMethodLSODA.cpp.

References CCopasiMessage::ERROR, CCopasiMessage::EXCEPTION, CModel::getCompartments(), CTrajectoryProblem::getDuration(), CModel::getEvents(), CModel::getMetabolites(), CCopasiProblem::getModel(), CModel::getModelValues(), CModel::getNumMetabs(), CModel::getNumModelValues(), CCopasiParameter::getValue(), CTrajectoryMethod::isValidProblem(), MCTrajectoryMethod, mLowerStochLimit, mUpperStochLimit, CModelEntity::ODE, CCopasiParameter::Value::pDOUBLE, 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::EXCEPTION, 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;
         }
     }
 
   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;
         }
     }
 
   std::string message = pTP->getModel()->suitableForStochasticSimulation();
 
   if (message != "")
     {
       //model not suitable, message describes the problem
       CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
       return false;
     }
 
   mLowerStochLimit = * getValue("Lower Limit").pDOUBLE;
   mUpperStochLimit = * getValue("Upper Limit").pDOUBLE;
 
   if (mLowerStochLimit > mUpperStochLimit)
     CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 4, mLowerStochLimit, mUpperStochLimit);
 
   //events are not supported at the moment
   if (pTP->getModel()->getEvents().size() > 0)
     {
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 23);
       return false;
     }
 
   return true;
 }

bool CHybridMethodLSODA::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 1750 of file CHybridMethodLSODA.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 CHybridMethodLSODA::outputData	(	std::ostream &	os,
		C_INT32	mode
	)

protected

Prints out data on standard output.

Prints out data on standard output. Deprecated.

Definition at line 1384 of file CHybridMethodLSODA.cpp.

References C_INT32, CState::getTime(), mOutputCounter, CTrajectoryMethod::mpCurrentState, mpMetabolites, and CCopasiVector< T >::size().

 {
   static unsigned C_INT32 counter = 0;
   size_t i;
 
   switch (mode)
     {
       case 0:
 
         if (mOutputCounter == (counter++))
           {
             counter = 0;
             os << mpCurrentState->getTime() << " : ";
 
             for (i = 0; i < mpMetabolites->size(); i++)
               {
                 os << (*mpMetabolites)[i]->getValue() << " ";
               }
 
             os << std::endl;
           }
 
         break;
 
       case 1:
         os << mpCurrentState->getTime() << " : ";
 
         for (i = 0; i < mpMetabolites->size(); i++)
           {
             os << (*mpMetabolites)[i]->getValue() << " ";
           }
 
         os << std::endl;
         break;
 
       default:
         ;
     }
 
   return;
 }

void CHybridMethodLSODA::outputDebug	(	std::ostream &	os,
		size_t	level
	)

protected

Prints out various data on standard output for debugging purposes.

Definition at line 1429 of file CHybridMethodLSODA.cpp.

References CCopasiObject::getObjectName(), CState::getTime(), LOW, mAmu, mAmuOld, mDG, mFirstReactionFlag, mLocalBalances, mLocalSubstrates, mLowerStochLimit, mMaxBalance, mMaxIntBeforeStep, mMaxSteps, mMetab2React, mMetabFlags, mNumVariableMetabs, mPartitioningInterval, CTrajectoryMethod::mpCurrentState, mpMetabolites, mPQ, mpRandomGenerator, mpReactions, mReactionFlags, mStepsAfterPartitionSystem, mStepsize, mStoi, mUpdateSet, mUpperStochLimit, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CCopasiVector< T >::size(), and temp.

 {
   size_t i, j;
   std::set< size_t >::iterator iter, iterEnd;
 
   os << "outputDebug(" << level << ") *********************************************** BEGIN" << std::endl;
 
   switch (level)
     {
       case 0:                              // Everything !!!
         os << " Name: " << CCopasiParameter::getObjectName() << std::endl;
         os << "current time: " << mpCurrentState->getTime() << std::endl;
         os << "mNumVariableMetabs: " << mNumVariableMetabs << std::endl;
         os << "mMaxSteps: " << mMaxSteps << std::endl;
         os << "mMaxBalance: " << mMaxBalance << std::endl;
         os << "mMaxIntBeforeStep: " << mMaxIntBeforeStep << std::endl;
         os << "mpReactions.size(): " << mpReactions->size() << std::endl;
 
         for (i = 0; i < mpReactions->size(); i++)
           os << *(*mpReactions)[i] << std::endl;
 
         os << "mpMetabolites.size(): " << mpMetabolites->size() << std::endl;
 
         for (i = 0; i < mpMetabolites->size(); i++)
           os << *(*mpMetabolites)[i] << std::endl;
 
         os << "mStoi: " << std::endl;
 
         for (i = 0; i < mStoi.numRows(); i++)
           {
             for (j = 0; j < mStoi.numCols(); j++)
               os << mStoi[i][j] << " ";
 
             os << std::endl;
           }
 
         os << "temp: ";
 
         for (i = 0; i < mNumVariableMetabs; i++)
           os << temp[i] << " ";
 
         os << std::endl;
 
         os << "mReactionFlags: " << std::endl;
 
         for (i = 0; i < mLocalBalances.size(); i++)
           os << mReactionFlags[i];
 
         os << "mFirstReactionFlag: " << std::endl;
         if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
 
         os << "mMetabFlags: " << std::endl;
 
         for (i = 0; i < mMetabFlags.size(); i++)
           {
             if (mMetabFlags[i] == LOW)
               os << "LOW ";
             else
               os << "HIGH ";
           }
 
         os << std::endl;
         os << "mLocalBalances: " << std::endl;
 
         for (i = 0; i < mLocalBalances.size(); i++)
           {
             for (j = 0; j < mLocalBalances[i].size(); j++)
               os << mLocalBalances[i][j];
 
             os << std::endl;
           }
 
         os << "mLocalSubstrates: " << std::endl;
 
         for (i = 0; i < mLocalSubstrates.size(); i++)
           {
             for (j = 0; j < mLocalSubstrates[i].size(); j++)
               os << mLocalSubstrates[i][j];
 
             os << std::endl;
           }
 
         os << "mLowerStochLimit: " << mLowerStochLimit << std::endl;
         os << "mUpperStochLimit: " << mUpperStochLimit << std::endl;
         //deprecated:      os << "mOutputCounter: " << mOutputCounter << endl;
         os << "mPartitioningInterval: " << mPartitioningInterval << std::endl;
         os << "mStepsAfterPartitionSystem: " << mStepsAfterPartitionSystem << std::endl;
         os << "mStepsize: " << mStepsize << std::endl;
         os << "mMetab2React: " << std::endl;
 
         for (i = 0; i < mMetab2React.size(); i++)
           {
             os << i << ": ";
 
             for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; ++iter)
               os << *iter << " ";
 
             os << std::endl;
           }
 
         os << "mAmu: " << std::endl;
 
         for (i = 0; i < mpReactions->size(); i++)
           os << mAmu[i] << " ";
 
         os << std::endl;
         os << "mAmuOld: " << std::endl;
 
         for (i = 0; i < mpReactions->size(); i++)
           os << mAmuOld[i] << " ";
 
         os << std::endl;
         os << "mUpdateSet: " << std::endl;
 
         for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
           os << *iter;
 
         os << std::endl;
         os << "mpRandomGenerator: " << mpRandomGenerator << std::endl;
         os << "mDG: " << std::endl << mDG;
         os << "mPQ: " << std::endl << mPQ;
         os << "Particle numbers: " << std::endl;
 
         for (i = 0; i < mpMetabolites->size(); i++)
           {
             os << (*mpMetabolites)[i]->getValue() << " ";
           }
 
         os << std::endl;
         break;
 
       case 1:                               // Variable values only
         os << "current time: " << mpCurrentState->getTime() << std::endl;
         /*
         case 1:
         os << "mTime: " << mpCurrentState->getTime() << std::endl;
         os << "oldState: ";
         for (i = 0; i < mDim; i++)
           os << oldState[i] << " ";
         os << std::endl;
         os << "x: ";
         for (i = 0; i < mDim; i++)
           os << x[i] << " ";
         os << std::endl;
         os << "y: ";
         for (i = 0; i < mDim; i++)
           os << y[i] << " ";
         os << std::endl;
         os << "increment: ";
         for (i = 0; i < mDim; i++)
           os << increment[i] << " ";
         os << std::endl;*/
         os << "temp: ";
 
         for (i = 0; i < mNumVariableMetabs; i++)
           os << temp[i] << " ";
 
         os << std::endl;
 
         os << "mReactionFlags: " << std::endl;
 
         for (i = 0; i < mLocalBalances.size(); i++)
           os << mReactionFlags[i];
 
         os << "mFirstReactionFlag: " << std::endl;
         if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
 
         os << "mMetabFlags: " << std::endl;
 
         for (i = 0; i < mMetabFlags.size(); i++)
           {
             if (mMetabFlags[i] == LOW)
               os << "LOW ";
             else
               os << "HIGH ";
           }
 
         os << std::endl;
         os << "mAmu: " << std::endl;
 
         for (i = 0; i < mpReactions->size(); i++)
           os << mAmu[i] << " ";
 
         os << std::endl;
         os << "mAmuOld: " << std::endl;
 
         for (i = 0; i < mpReactions->size(); i++)
           os << mAmuOld[i] << " ";
 
         os << std::endl;
         os << "mUpdateSet: " << std::endl;
 
         for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
           os << *iter;
 
         os << std::endl;
         os << "mPQ: " << std::endl << mPQ;
         os << "Particle numbers: " << std::endl;
 
         for (i = 0; i < mpMetabolites->size(); i++)
           {
             os << (*mpMetabolites)[i]->getValue() << " ";
           }
 
         os << std::endl;
         break;
 
       case 2:
         break;
 
       default:
         ;
     }
 
   os << "outputDebug(" << level << ") ************************************************* END" << std::endl;
   return;
 }

void CHybridMethodLSODA::partitionSystem ( )

protected

Updates the partitioning of the system depending on the particle numbers present.

Definition at line 1175 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, calculateAmu(), generateReactionTime(), CCopasiProblem::getModel(), CModel::getTime(), CCopasiParameter::getValue(), HIGH, insertDeterministicReaction(), CIndexedPriorityQueue::insertStochReaction(), LOW, mAmu, mAmuOld, mLowerStochLimit, mMetab2React, mMetabFlags, mNumVariableMetabs, mpMetabolites, CTrajectoryMethod::mpProblem, mPQ, mReactionFlags, mRestartLSODA, mUpperStochLimit, CCopasiParameter::mValue, removeDeterministicReaction(), and CIndexedPriorityQueue::removeStochReaction().

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

 {
   size_t i;
   std::set <size_t>::iterator iter, iterEnd;
   C_FLOAT64 key;
 
   for (i = 0; i < mNumVariableMetabs; i++)
     {
       if ((mMetabFlags[i] == LOW) && ((*mpMetabolites)[i]->getValue() >= mUpperStochLimit))
         {
           mRestartLSODA = true;
           mMetabFlags[i] = HIGH;
 
           // go through all corresponding reactions and update flags
           for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; iter++)
             {
               mReactionFlags[*iter].mValue--;
 
               // if reaction gets deterministic, insert it into the linked list of deterministic reactions
               if (mReactionFlags[*iter].mValue == 0)
                 {
                   insertDeterministicReaction(*iter);
                   mPQ.removeStochReaction(*iter);
                 }
             }
         }
 
       if ((mMetabFlags[i] == HIGH) && ((*mpMetabolites)[i]->getValue() < mLowerStochLimit))
         {
           mRestartLSODA = true;
           mMetabFlags[i] = LOW;
           (*mpMetabolites)[i]->setValue(floor((*mpMetabolites)[i]->getValue()));
           (*mpMetabolites)[i]->refreshConcentration();
 
           // go through all corresponding reactions and update flags
           for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; iter++)
             {
               if (mReactionFlags[*iter].mValue == 0)
                 {
                   removeDeterministicReaction(*iter);
                   /*
                     mPQ.insertStochReaction(*iter, 1234567.8);  // juergen: have to beautify this, number has to be the biggest C_FLOAT64 !!!
                   */
                   calculateAmu(*iter);
                   mAmuOld[*iter] = mAmu[*iter];
                   //key = mpCurrentState->getTime() + generateReactionTime(*iter); NATALIA
 
                   key = mpProblem->getModel()->getTime() + generateReactionTime(*iter);
 
                   mPQ.insertStochReaction(*iter, key);
                 }
 
               mReactionFlags[*iter].mValue++;
             }
         }
     }
 
   return;
 }

void CHybridMethodLSODA::removeDeterministicReaction ( size_t rIndex )

protected

Removes a deterministic reaction from the linked list in the vector mReactionFlags.

Parameters

rIndex A size_t giving the index of the reaction to be removed from the list of deterministic reactions.

Removes a deterministic reaction from the linked list in the array mReactionFlags.

Parameters

rIndex A size_t giving the index of the reaction to be removed from the list of deterministic reactions.

Definition at line 1277 of file CHybridMethodLSODA.cpp.

References mFirstReactionFlag, CHybridLSODAStochFlag::mpPrev, and mReactionFlags.

Referenced by partitionSystem().

 {
   if (mReactionFlags[rIndex].mpPrev != NULL)
     // reaction is deterministic
     {
       if (&mReactionFlags[rIndex] != mFirstReactionFlag)
         // reactionFlag is not the first in the linked list
         {
           mReactionFlags[rIndex].mpPrev->mpNext = mReactionFlags[rIndex].mpNext;
 
           if (mReactionFlags[rIndex].mpNext != NULL)
             mReactionFlags[rIndex].mpNext->mpPrev = mReactionFlags[rIndex].mpPrev;
         }
       else
         // reactionFlag is the first in the linked list
         {
           if (mReactionFlags[rIndex].mpNext != NULL) // reactionFlag is not the only one in the linked list
             {
               mFirstReactionFlag = mReactionFlags[rIndex].mpNext;
               mFirstReactionFlag->mpPrev = mFirstReactionFlag;
             }
           else // reactionFlag is the only one in the linked list
             {
               mFirstReactionFlag = NULL;
             }
         }
     }
 
   mReactionFlags[rIndex].mpPrev = NULL;
   mReactionFlags[rIndex].mpNext = NULL;
   return;
 }

void CHybridMethodLSODA::setState ( std::vector< C_FLOAT64 > & source )

protected

Writes the state specified in the array source into the model. Length of the vector source must be mNumVariableMetabs. (Number of non-fixed metabolites in the model).

Parameters

source A vector reference with length mNumVariableMetabs, holding the state of the system to be set in the model

void CHybridMethodLSODA::setupBalances ( )

protected

Sets up an internal representation of the balances for each reaction. This is done in order to be able to deal with fixed metabolites and to avoid a time consuming search for the indices of metabolites in the model.

Definition at line 858 of file CHybridMethodLSODA.cpp.

References C_INT32, CModelEntity::FIXED, CCopasiVector< T >::getIndex(), CModel::getMetabolitesX(), CModelEntity::getStatus(), CHybridLSODABalance::mIndex, mLocalBalances, mLocalSubstrates, mMaxBalance, mMaxIntBeforeStep, mMaxSteps, CHybridLSODABalance::mMultiplicity, CHybridLSODABalance::mpMetabolite, mpModel, mpReactions, and CCopasiVector< T >::size().

Referenced by initMethod().

 {
   size_t i, j;
   CHybridLSODABalance newElement;
   C_INT32 maxBalance = 0;
   size_t numReactions;
 
   numReactions = mpReactions->size();
   mLocalBalances.clear();
   mLocalBalances.resize(numReactions);
   mLocalSubstrates.clear();
   mLocalSubstrates.resize(numReactions);
 
   for (i = 0; i < numReactions; i++)
     {
       const CCopasiVector <CChemEqElement> * balances =
         &(*mpReactions)[i]->getChemEq().getBalances();
 
       for (j = 0; j < balances->size(); j++)
         {
           newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
           newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
           // + 0.5 to get a rounding out of the static_cast to C_INT32!
           newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
 
           if ((newElement.mpMetabolite->getStatus()) != CModelEntity::FIXED)
             {
               if (newElement.mMultiplicity > maxBalance) maxBalance = newElement.mMultiplicity;
 
               mLocalBalances[i].push_back(newElement); // element is copied for the push_back
             }
         }
 
       balances = &(*mpReactions)[i]->getChemEq().getSubstrates();
 
       for (j = 0; j < balances->size(); j++)
         {
           newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
           newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
           // + 0.5 to get a rounding out of the static_cast to C_INT32!
           newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
 
           mLocalSubstrates[i].push_back(newElement); // element is copied for the push_back
         }
     }
 
   mMaxBalance = maxBalance;
   mMaxIntBeforeStep = INT_MAX - 1 - mMaxSteps * mMaxBalance;
 
   return;
 }

void CHybridMethodLSODA::setupDependencyGraph ( )

protected

Sets up the dependency graph

Sets up the dependency graph.

Definition at line 913 of file CHybridMethodLSODA.cpp.

References CDependencyGraph::addDependent(), CDependencyGraph::clear(), getAffects(), getDependsOn(), mDG, mpReactions, and CCopasiVector< T >::size().

Referenced by initMethod().

 {
   mDG.clear();
   std::vector< std::set<std::string>* > DependsOn;
   std::vector< std::set<std::string>* > Affects;
   size_t numReactions = mpReactions->size();
   size_t i, j;
 
   // Do for each reaction:
   for (i = 0; i < numReactions; 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.
   for (i = 0; i < numReactions; i++)
     {
       for (j = 0; j < numReactions; 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);
     }
 
   // Delete the memory allocated in getDependsOn() and getAffects()
   // since this is allocated in other functions.
   for (i = 0; i < numReactions; i++)
     {
       delete DependsOn[i];
       delete Affects[i];
     }
 
   return;
 }

void CHybridMethodLSODA::setupMetab2React ( )

protected

Creates for each metabolite a set of reaction indices. If the metabolite participates in a reaction as substrate or product this reaction is added to the corresponding set.

Creates for each metabolite a set of reaction indices. If the metabolite participates in a reaction as substrate or product (that means: balance != 0) this reaction is added to the corresponding set.

Definition at line 977 of file CHybridMethodLSODA.cpp.

References mLocalBalances, mMetab2React, mpMetabolites, and CCopasiVector< T >::size().

Referenced by initMethod().

 {
   size_t i, j;
   size_t metaboliteIndex;
 
   // Resize mMetab2React and create an initial set for each metabolite
   mMetab2React.clear();
   mMetab2React.resize(mpMetabolites->size());
 
   // Iterate over all reactions
   for (i = 0; i < mLocalBalances.size(); i++)
     {
       // Get the set of metabolites which take part in this reaction
       for (j = 0; j < mLocalBalances[i].size(); j++)
         {
           // find metaboliteIndex and insert the reaction into the set
           metaboliteIndex = mLocalBalances[i][j].mIndex;
           mMetab2React[metaboliteIndex].insert(i);
         }
     }
 
   return;
 }

void CHybridMethodLSODA::setupMetab2ReactComplete ( )

protected

Creates for each metabolite a set of reaction indices. Each reaction is dependent on each metabolite resulting in a complete switch.

Definition at line 1043 of file CHybridMethodLSODA.cpp.

References mMetab2React, mpMetabolites, mpReactions, and CCopasiVector< T >::size().

 {
   size_t i, j;
 
   // Resize mMetab2React and create an initial set for each metabolite
   mMetab2React.resize(mpMetabolites->size());
 
   // Iterate over all metabolites
   for (i = 0; i < mpMetabolites->size(); i++)
     {
       // Iterate over all reactions
       for (j = 0; j < mpReactions->size(); j++)
         {
           mMetab2React[i].insert(j);
         }
     }
 
   return;
 }

void CHybridMethodLSODA::setupMetab2ReactPlusModifier ( )

protected

Creates for each metabolite a set of reaction indices. If the metabolite participates in a reaction as substrate, product or modifier this reaction is added to the corresponding set.

Definition at line 1006 of file CHybridMethodLSODA.cpp.

References getParticipatesIn(), mMetab2React, mpMetabolites, mpReactions, and CCopasiVector< T >::size().

 {
   std::vector< std::set<size_t>* > participatesIn;
   size_t numReactions = mpReactions->size();
   size_t i;
 
   // Resize mMetab2React and create an initial set for each metabolite
   mMetab2React.resize(mpMetabolites->size());
 
   // Do for each reaction:
   for (i = 0; i < numReactions; i++)
     {
       participatesIn.push_back(getParticipatesIn(i));
     }
 
   // Iterate over all reactions
   for (i = 0; i < numReactions; i++)
     {
       // Get the set of metabolites which take part in this reaction
       std::set<size_t>::iterator iter = participatesIn[i]->begin();
 
       for (; iter != participatesIn[i]->end(); iter++)
         mMetab2React[*iter].insert(i);
     }
 
   for (i = 0; i < numReactions; i++)
     {
       delete participatesIn[i];
     }
 
   return;
 }

void CHybridMethodLSODA::setupPartition ( )

protected

Creates an initial partitioning of the system. Deterministic and stochastic reactions are determined. The array mStochReactions is initialized.

Creates an initial partitioning of the system. Deterministic and stochastic reactions are determined. The vector mReactionFlags and the vector mMetabFlags are initialized.

Definition at line 1068 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, CCopasiParameter::getValue(), HIGH, LOW, mFirstReactionFlag, mLocalBalances, mLowerStochLimit, mMetabFlags, mNumVariableMetabs, mpMetabolites, CHybridLSODAStochFlag::mpNext, mReactionFlags, mUpperStochLimit, and CCopasiParameter::mValue.

Referenced by initMethod().

 {
   size_t i, j;
   CHybridLSODAStochFlag * prevFlag;
   C_FLOAT64 averageStochLimit = (mUpperStochLimit + mLowerStochLimit) / 2;
 
   // initialize vector mMetabFlags
   mMetabFlags.clear();
   mMetabFlags.resize(mNumVariableMetabs);
 
   for (i = 0; i < mMetabFlags.size(); i++)
     {
       if ((*mpMetabolites)[i]->getValue() < averageStochLimit)
         {
           mMetabFlags[i] = LOW;
           (*mpMetabolites)[i]->setValue(floor((*mpMetabolites)[i]->getValue()));
           (*mpMetabolites)[i]->refreshConcentration();
         }
       else
         mMetabFlags[i] = HIGH;
     }
 
   // initialize vector mReactionFlags
   mReactionFlags.clear();
   mReactionFlags.resize(mLocalBalances.size());
 
   for (i = 0; i < mLocalBalances.size(); i++)
     {
       mReactionFlags[i].mIndex = i;
       mReactionFlags[i].mValue = 0;
 
       for (j = 0; j < mLocalBalances[i].size(); j++)
         {
           if (mMetabFlags[mLocalBalances[i][j].mIndex] == LOW)
             {
               mReactionFlags[i].mValue++;
             }
         }
     }
 
   mFirstReactionFlag = NULL;
   prevFlag = NULL;
 
   for (i = 0; i < mLocalBalances.size(); i++)
     {
       if (mReactionFlags[i].mValue == 0)
         {
           if (mFirstReactionFlag != NULL)
             {
               prevFlag->mpNext = &mReactionFlags[i];
               mReactionFlags[i].mpPrev = prevFlag;
               prevFlag = &mReactionFlags[i];
             }
           else
             {
               mFirstReactionFlag = &mReactionFlags[i];
               mReactionFlags[i].mpPrev = &mReactionFlags[i]; // Important to distinguish between stochastic (prev == NULL) and deterministic (prev != NULL) reactions
               prevFlag = &mReactionFlags[i];
             }
         }
       else
         {
           mReactionFlags[i].mpPrev = NULL;
           mReactionFlags[i].mpNext = NULL;
         }
     }
 
   if (prevFlag != NULL)
     {
       prevFlag->mpNext = NULL;
     }
 
   return;
 }

void CHybridMethodLSODA::setupPriorityQueue ( C_FLOAT64 startTime = 0.0 )

protected

Sets up the priority queue.

Parameters

startTime The time at which the simulation starts.

Definition at line 1148 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, calculateAmu(), CIndexedPriorityQueue::clear(), generateReactionTime(), CIndexedPriorityQueue::initializeIndexPointer(), CIndexedPriorityQueue::insertStochReaction(), mPQ, mpReactions, mReactionFlags, and CCopasiVector< T >::size().

Referenced by initMethod().

 {
   size_t i;
   C_FLOAT64 time;
 
   mPQ.clear();
   mPQ.initializeIndexPointer(mpReactions->size());
 
   for (i = 0; i < mpReactions->size(); i++)
     {
       if (mReactionFlags[i].mpPrev == NULL) // Reaction is stochastic!
         {
           calculateAmu(i);
           time = startTime + generateReactionTime(i);
           mPQ.insertStochReaction(i, time);
         }
     }
 
   return;
 }

void CHybridMethodLSODA::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 253 of file CHybridMethodLSODA.cpp.

References CModel::deterministic, CCopasiProblem::getModel(), CModel::getState(), CModel::getTime(), initMethod(), mDoCorrection, mFirstMetabIndex, mHasAssignments, modelHasAssignments(), mpModel, CTrajectoryMethod::mpProblem, mpState, mRestartLSODA, and CModel::setState().

 {
   //*mpCurrentState = *initialState;
 
   mpProblem->getModel()->setState(*initialState);
   mRestartLSODA = true;
 
   mpState = new CState(mpProblem->getModel()->getState());
 
   mpModel = mpProblem->getModel();
   assert(mpModel);
 
   if (mpModel->getModelType() == CModel::deterministic)
     mDoCorrection = true;
   else
     mDoCorrection = false;
 
   mHasAssignments = modelHasAssignments(mpModel);
 
   mFirstMetabIndex = mpModel->getStateTemplate().getIndex(mpModel->getMetabolitesX()[0]);
 
 //mpProblem->getModel()->setState(*mpCurrentState);
 
   // call init of the simulation method, can be overloaded in derived classes
   // mDim[1] = (C_INT) (void *) this;
   initMethod(mpProblem->getModel()->getTime());
 
   return;
 }

CTrajectoryMethod::Status CHybridMethodLSODA::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 202 of file CHybridMethodLSODA.cpp.

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

 {
   // write the current state to the model
   //  mpProblem->getModel()->setState(mpCurrentState); // is that correct?
 
   // 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)
       {
         // throw exception or something like that
       }
 
   // do several steps
   C_FLOAT64 time = mpProblem->getModel()->getTime();
   C_FLOAT64 endTime = time + deltaT;
 
   for (i = 0; ((i < mMaxSteps) && (time < endTime)); i++)
     {
       time = doSingleStep(time, endTime);
     }
 
   //mpCurrentState->setTime(time);
 
   if ((i >= mMaxSteps) && (!mMaxStepsReached))
     {
       mMaxStepsReached = true; //only report this message once
       CCopasiMessage(CCopasiMessage::WARNING, "maximum number of reaction events was reached in at least one simulation step.\nThat means time intervals in the output may not be what you requested.");
     }
 
   // 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();
 
   //the task expects the result in mpCurrentState
   *mpCurrentState = mpProblem->getModel()->getState();
 
   mpCurrentState->setTime(time);
 
   return NORMAL;
 }

void CHybridMethodLSODA::updatePriorityQueue	(	size_t	rIndex,
		C_FLOAT64	time
	)

protected

Updates the priority queue.

Parameters

rIndex	A size_t giving the index of the fired reaction
time	A C_FLOAT64 holding the time taken by this reaction

Updates the priority queue.

Parameters

rIndex	A size_t giving the index of the fired reaction (-1, if no stochastic reaction has fired)
time	A C_FLOAT64 holding the current time

Definition at line 641 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, C_INVALID_INDEX, calculateAmu(), generateReactionTime(), mAmu, mAmuOld, mHasAssignments, mpModel, mPQ, mpReactions, mReactionFlags, mUpdateSet, CCopasiVector< T >::size(), CIndexedPriorityQueue::updateNode(), CModel::updateSimulatedValues(), and updateTauMu().

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

 {
   C_FLOAT64 newTime;
   size_t index;
   std::set <size_t>::iterator iter, iterEnd;
 
   //if the model contains assignments we use a less efficient loop over all (stochastic) reactions to capture all changes
   // we do not know the exact dependencies. TODO: this should be changed later in order to get a more efficient update scheme
   if (mHasAssignments)
     {
       mpModel->updateSimulatedValues(false);
 
       for (index = 0; index < mpReactions->size(); index++)
         {
           if (mReactionFlags[index].mpPrev == NULL) // Reaction is stochastic!
             {
               mAmuOld[index] = mAmu[index];
               calculateAmu(index);
 
               if (mAmuOld[index] != mAmu[index])
                 if (index != rIndex) updateTauMu(index, time);
             }
         }
     }
   else
     {
       // iterate through the set of affected reactions and update the stochastic ones in the priority queue
       for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
         {
           if (mReactionFlags[*iter].mpPrev == NULL) // reaction is stochastic!
             {
               index = *iter;
               mAmuOld[index] = mAmu[index];
               calculateAmu(index);
 
               if (*iter != rIndex) updateTauMu(index, time);
             }
         }
     }
 
   // draw new random number and update the reaction just fired
   if ((rIndex != C_INVALID_INDEX) && (mReactionFlags[rIndex].mpPrev == NULL))
     {
       // reaction is stochastic
       newTime = time + generateReactionTime(rIndex);
       mPQ.updateNode(rIndex, newTime);
     }
 
   // empty the mUpdateSet
   mUpdateSet.clear();
   return;
 }

void CHybridMethodLSODA::updateTauMu	(	size_t	rIndex,
		C_FLOAT64	time
	)

protected

Updates the putative reaction time of a stochastic reaction in the priority queue. The corresponding amu and amu_old must be set prior to the call of this method.

Parameters

rIndex	A size_t specifying the index of the reaction
time	A C_FLOAT64 specifying the current time

Updates the putative reaction time of a stochastic reaction in the priority queue. The corresponding amu and amu_old must be set prior to the call of this method.

Parameters

rIndex A size_t specifying the index of the reaction

Definition at line 787 of file CHybridMethodLSODA.cpp.

References C_FLOAT64, generateReactionTime(), CIndexedPriorityQueue::getKey(), mAmu, mAmuOld, mPQ, and CIndexedPriorityQueue::updateNode().

Referenced by updatePriorityQueue().

 {
   C_FLOAT64 newTime;
 
   // One must make sure that the calculation yields reasonable results even in the cases where mAmu=0 or mAmuOld=0 or both =0. Therefore mAmuOld=0 is checked. If mAmuOld equals 0, then a new random number has to be drawn, because tau equals inf and the stoch. information is lost. If both values equal 0, then tau should remain inf and the update of the queue can be skipped!
 
   if (mAmuOld[rIndex] == 0.0)
     {
       if (mAmu[rIndex] != 0.0)
         {
           newTime = time + generateReactionTime(rIndex);
           mPQ.updateNode(rIndex, newTime);
         }
     }
   else
     {
       newTime = time + (mAmuOld[rIndex] / mAmu[rIndex]) * (mPQ.getKey(rIndex) - time);
       mPQ.updateNode(rIndex, newTime);
     }
 
   return;
 }

Friends And Related Function Documentation

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

friend

Member Data Documentation

std::vector<C_FLOAT64> CHybridMethodLSODA::currentState

protected

Definition at line 545 of file CHybridMethodLSODA.h.

std::vector<C_FLOAT64> CHybridMethodLSODA::mAmu

protected

The propensities of the stochastic reactions.

Definition at line 603 of file CHybridMethodLSODA.h.

Referenced by calculateAmu(), generateReactionTime(), initMethod(), insertDeterministicReaction(), outputDebug(), partitionSystem(), updatePriorityQueue(), and updateTauMu().

std::vector<C_FLOAT64> CHybridMethodLSODA::mAmuOld

protected

Definition at line 604 of file CHybridMethodLSODA.h.

Referenced by initMethod(), insertDeterministicReaction(), outputDebug(), partitionSystem(), updatePriorityQueue(), and updateTauMu().

C_FLOAT64 CHybridMethodLSODA::mAtol

protected

Absolute tolerance.

Definition at line 693 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and integrateDeterministicPart().

Data CHybridMethodLSODA::mData

protected

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

Definition at line 648 of file CHybridMethodLSODA.h.

Referenced by CHybridMethodLSODA(), initMethod(), and integrateDeterministicPart().

bool CHybridMethodLSODA::mDefaultAtol

protected

Definition at line 688 of file CHybridMethodLSODA.h.

Referenced by initMethod().

CDependencyGraph CHybridMethodLSODA::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 620 of file CHybridMethodLSODA.h.

Referenced by fireReaction(), integrateDeterministicPart(), outputDebug(), and setupDependencyGraph().

bool CHybridMethodLSODA::mDoCorrection

protected

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

Definition at line 524 of file CHybridMethodLSODA.h.

Referenced by calculateAmu(), and start().

CVector< C_FLOAT64 > CHybridMethodLSODA::mDWork

protected

Definition at line 702 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and integrateDeterministicPart().

C_FLOAT64 CHybridMethodLSODA::mEndt

protected

Requested end time.

Definition at line 668 of file CHybridMethodLSODA.h.

std::ostringstream CHybridMethodLSODA::mErrorMsg

protected

Definition at line 698 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

size_t CHybridMethodLSODA::mFirstMetabIndex

protected

index of the first metab in CState

Definition at line 491 of file CHybridMethodLSODA.h.

Referenced by start(), and step().

CHybridLSODAStochFlag* CHybridMethodLSODA::mFirstReactionFlag

protected

Definition at line 553 of file CHybridMethodLSODA.h.

Referenced by calculateDerivative(), CHybridNextReactionLSODAMethod::doSingleStep(), insertDeterministicReaction(), integrateDeterministicPart(), outputDebug(), removeDeterministicReaction(), and setupPartition().

bool CHybridMethodLSODA::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 727 of file CHybridMethodLSODA.h.

Referenced by start(), and updatePriorityQueue().

CVector< C_INT > CHybridMethodLSODA::mIWork

protected

Definition at line 703 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and integrateDeterministicPart().

C_INT CHybridMethodLSODA::mJType

protected

Definition at line 704 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

std::vector<std::vector <CHybridLSODABalance> > CHybridMethodLSODA::mLocalBalances

protected

Internal representation of the balances of each reaction. The index of each metabolite in the reaction is provided.

Definition at line 565 of file CHybridMethodLSODA.h.

Referenced by fireReaction(), getAffects(), outputDebug(), setupBalances(), setupMetab2React(), and setupPartition().

std::vector<std::vector <CHybridLSODABalance> > CHybridMethodLSODA::mLocalSubstrates

protected

Internal representation of the substrates of each reaction. The index of each substrate-metabolite is provided.

Definition at line 571 of file CHybridMethodLSODA.h.

Referenced by calculateAmu(), outputDebug(), and setupBalances().

C_FLOAT64 CHybridMethodLSODA::mLowerStochLimit

protected

Limit for particle numbers. If a particle number is < StochLimit the corresponding reactions must be simulated stochastically.

Definition at line 577 of file CHybridMethodLSODA.h.

Referenced by initMethod(), isValidProblem(), outputDebug(), partitionSystem(), and setupPartition().

CLSODA CHybridMethodLSODA::mLSODA

protected

Definition at line 700 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

C_INT CHybridMethodLSODA::mLsodaStatus

protected

LSODA state.

Definition at line 673 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

C_INT32 CHybridMethodLSODA::mMaxBalance

protected

maximal increase of a particle number in one step.

Definition at line 514 of file CHybridMethodLSODA.h.

Referenced by outputDebug(), and setupBalances().

C_INT32 CHybridMethodLSODA::mMaxIntBeforeStep

protected

This is set to maxint - mMaxSteps*mMaxBalance

Definition at line 519 of file CHybridMethodLSODA.h.

Referenced by outputDebug(), setupBalances(), and step().

unsigned C_INT32 CHybridMethodLSODA::mMaxSteps

protected

Max number of doSingleStep() per step()

Definition at line 496 of file CHybridMethodLSODA.h.

Referenced by initMethod(), outputDebug(), setupBalances(), and step().

bool CHybridMethodLSODA::mMaxStepsReached

protected

Definition at line 498 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and step().

std::vector<std::set <size_t> > CHybridMethodLSODA::mMetab2React

protected

Vector of relations between metabolites to reactions.

Definition at line 598 of file CHybridMethodLSODA.h.

Referenced by outputDebug(), partitionSystem(), setupMetab2React(), setupMetab2ReactComplete(), and setupMetab2ReactPlusModifier().

std::vector<metabStatus> CHybridMethodLSODA::mMetabFlags

protected

Vector holding information on the status of metabolites. They can have low or high particle numbers.

Definition at line 559 of file CHybridMethodLSODA.h.

Referenced by outputDebug(), partitionSystem(), and setupPartition().

size_t CHybridMethodLSODA::mNumVariableMetabs

protected

Dimension of the system. Total number of metabolites.

Definition at line 486 of file CHybridMethodLSODA.h.

Referenced by evalF(), initMethod(), outputDebug(), partitionSystem(), and setupPartition().

size_t CHybridMethodLSODA::mOutputCounter

protected

Output counter.

Definition at line 720 of file CHybridMethodLSODA.h.

Referenced by outputData().

std::ofstream CHybridMethodLSODA::mOutputFile

protected

File output stream to write data.

Definition at line 710 of file CHybridMethodLSODA.h.

std::string CHybridMethodLSODA::mOutputFileName

protected

Output filename.

Definition at line 715 of file CHybridMethodLSODA.h.

unsigned C_INT32 CHybridMethodLSODA::mPartitioningInterval

protected

The system gets repartitioned after this number of elementary steps.

Definition at line 583 of file CHybridMethodLSODA.h.

Referenced by CHybridNextReactionLSODAMethod::doSingleStep(), initMethod(), and outputDebug().

CCopasiVector<CMetab>* CHybridMethodLSODA::mpMetabolites

protected

A pointer to the metabolites of the model.

Definition at line 534 of file CHybridMethodLSODA.h.

Referenced by initMethod(), outputData(), outputDebug(), partitionSystem(), setupMetab2React(), setupMetab2ReactComplete(), setupMetab2ReactPlusModifier(), and setupPartition().

CModel* CHybridMethodLSODA::mpModel

protected

Pointer to the model

Definition at line 481 of file CHybridMethodLSODA.h.

Referenced by cleanup(), initMethod(), integrateDeterministicPart(), setupBalances(), start(), and updatePriorityQueue().

CIndexedPriorityQueue CHybridMethodLSODA::mPQ

protected

The set of putative stochastic (!) reactions and associated times at which each reaction occurs. This is represented as a priority queue, sorted by the reaction time. This heap changes dynamically as stochastic reactions become deterministic (delete this reaction from the queue) or vice versa (insert a new reaction into the queue).

Definition at line 629 of file CHybridMethodLSODA.h.

Referenced by CHybridNextReactionLSODAMethod::doSingleStep(), getStochTimeAndIndex(), outputDebug(), partitionSystem(), setupPriorityQueue(), updatePriorityQueue(), and updateTauMu().

CRandom* CHybridMethodLSODA::mpRandomGenerator

protected

The random number generator.

Definition at line 614 of file CHybridMethodLSODA.h.

Referenced by CHybridMethodLSODA(), cleanup(), generateReactionTime(), initMethod(), and outputDebug().

const CCopasiVectorNS<CReaction>* CHybridMethodLSODA::mpReactions

protected

A pointer to the reactions of the model.

Definition at line 529 of file CHybridMethodLSODA.h.

Referenced by checkModel(), getDependsOn(), initMethod(), outputDebug(), setupBalances(), setupDependencyGraph(), setupMetab2ReactComplete(), setupMetab2ReactPlusModifier(), setupPriorityQueue(), and updatePriorityQueue().

CState* CHybridMethodLSODA::mpState

protected

A pointer to the current state in complete model view.

Definition at line 642 of file CHybridMethodLSODA.h.

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

unsigned C_INT32 CHybridMethodLSODA::mRandomSeed

protected

The random seed to use.

Definition at line 509 of file CHybridMethodLSODA.h.

Referenced by initMethod().

std::vector<CHybridLSODAStochFlag> CHybridMethodLSODA::mReactionFlags

protected

Vector to hold information about how many metabolites of a reaction have low particle numbers. If no metabolite of one reaction has low particle numbers this reaction will be simulated det.

Definition at line 552 of file CHybridMethodLSODA.h.

Referenced by insertDeterministicReaction(), outputDebug(), partitionSystem(), removeDeterministicReaction(), setupPartition(), setupPriorityQueue(), and updatePriorityQueue().

bool CHybridMethodLSODA::mReducedModel

protected

Definition at line 678 of file CHybridMethodLSODA.h.

bool CHybridMethodLSODA::mRestartLSODA

protected

this flag indicates whether the next LSODA call needs to reinitialize the integrater (if the partitioning has changed or a stochastic reaction was fired)

Definition at line 637 of file CHybridMethodLSODA.h.

Referenced by fireReaction(), integrateDeterministicPart(), partitionSystem(), and start().