#include <CHybridMethodODE45.h>

Inheritance diagram for CHybridMethodODE45:

Collaboration diagram for CHybridMethodODE45:

Classes
struct	Data

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

	~CHybridMethodODE45 ()

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

const CVector< C_INT > &	getRoots () const

void	setCurrentState (CState *currentState)

void	setProblem (CTrajectoryProblem *problem)

virtual void	stateChanged ()

	~CTrajectoryMethod ()

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

const CCopasiMethod::SubType &	getSubType () const

const CCopasiTask::Type &	getType () const

virtual void	load (CReadConfig &configBuffer, CReadConfig::Mode mode=CReadConfig::SEARCH)

virtual void	print (std::ostream *ostream) const

virtual void	printResult (std::ostream *ostream) const

virtual bool	setCallBack (CProcessReport *pCallBack)

virtual	~CCopasiMethod ()

Public Member Functions inherited from CCopasiParameterGroup
bool	addGroup (const std::string &name)

bool	addParameter (const CCopasiParameter &parameter)

bool	addParameter (const std::string &name, const CCopasiParameter::Type type)

template<class CType >
bool	addParameter (const std::string &name, const CCopasiParameter::Type type, const CType &value)

void	addParameter (CCopasiParameter *pParameter)

CCopasiParameterGroup *	assertGroup (const std::string &name)

template<class CType >
CCopasiParameter *	assertParameter (const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)

index_iterator	beginIndex () const

name_iterator	beginName () const

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

	CCopasiParameterGroup (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &objectType="ParameterGroup")

void	clear ()

index_iterator	endIndex () const

name_iterator	endName () const

CCopasiParameterGroup *	getGroup (const std::string &name)

const CCopasiParameterGroup *	getGroup (const std::string &name) const

CCopasiParameterGroup *	getGroup (const size_t &index)

const CCopasiParameterGroup *	getGroup (const size_t &index) const

size_t	getIndex (const std::string &name) const

std::string	getKey (const std::string &name) const

std::string	getKey (const size_t &index) const

virtual const std::string &	getName (const size_t &index) const

virtual const CObjectInterface *	getObject (const CCopasiObjectName &cn) const

CCopasiParameter *	getParameter (const std::string &name)

const CCopasiParameter *	getParameter (const std::string &name) const

CCopasiParameter *	getParameter (const size_t &index)

const CCopasiParameter *	getParameter (const size_t &index) const

CCopasiParameter::Type	getType (const std::string &name) const

CCopasiParameter::Type	getType (const size_t &index) const

std::string	getUniqueParameterName (const CCopasiParameter *pParameter) const

const CCopasiParameter::Value &	getValue (const std::string &name) const

const CCopasiParameter::Value &	getValue (const size_t &index) const

CCopasiParameter::Value &	getValue (const std::string &name)

CCopasiParameter::Value &	getValue (const size_t &index)

CCopasiParameterGroup &	operator= (const CCopasiParameterGroup &rhs)

bool	removeParameter (const std::string &name)

bool	removeParameter (const size_t &index)

template<class CType >
bool	setValue (const std::string &name, const CType &value)

template<class CType >
bool	setValue (const size_t &index, const CType &value)

size_t	size () const

bool	swap (const size_t &iFrom, const size_t &iTo)

bool	swap (index_iterator &from, index_iterator &to)

virtual	~CCopasiParameterGroup ()

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

	CCopasiParameter (const std::string &name, const Type &type, const void pValue=NULL, const CCopasiContainer pParent=NULL, const std::string &objectType="Parameter")

virtual CCopasiObjectName	getCN () const

virtual const std::string &	getKey () const

virtual std::string	getObjectDisplayName (bool regular=true, bool richtext=false) const

const CCopasiParameter::Type &	getType () const

const Value &	getValue () const

Value &	getValue ()

virtual void *	getValuePointer () const

CCopasiObject *	getValueReference () const

bool	isValidValue (const C_FLOAT64 &value) const

bool	isValidValue (const C_INT32 &value) const

bool	isValidValue (const unsigned C_INT32 &value) const

bool	isValidValue (const bool &value) const

bool	isValidValue (const std::string &value) const

bool	isValidValue (const CCopasiObjectName &value) const

bool	isValidValue (const std::vector< CCopasiParameter * > &value) const

CCopasiParameter &	operator= (const CCopasiParameter &rhs)

template<class CType >
bool	setValue (const CType &value)

bool	setValue (const std::vector< CCopasiParameter * > &value)

virtual	~CCopasiParameter ()

Public Member Functions inherited from CCopasiContainer
virtual bool	add (CCopasiObject *pObject, const bool &adopt=true)

	CCopasiContainer (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=CCopasiObject::Container)

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

virtual std::string	getChildObjectUnits (const CCopasiObject *pObject) const

virtual const objectMap &	getObjects () const

virtual std::string	getUnits () const

virtual const CCopasiObject *	getValueObject () const

virtual bool	remove (CCopasiObject *pObject)

virtual	~CCopasiContainer ()

Public Member Functions inherited from CCopasiObject
void	addDirectDependency (const CCopasiObject *pObject)

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

void	clearDirectDependencies ()

void	clearRefresh ()

bool	dependsOn (DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const

void	getAllDependencies (DataObjectSet &dependencies, const DataObjectSet &context) const

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

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

CCopasiDataModel *	getObjectDataModel ()

const CCopasiDataModel *	getObjectDataModel () const

const std::string &	getObjectName () const

CCopasiContainer *	getObjectParent () const

const std::string &	getObjectType () const

virtual const CObjectInterface::ObjectSet &	getPrerequisites () const

virtual Refresh *	getRefresh () const

UpdateMethod *	getUpdateMethod () const

bool	hasCircularDependencies (DataObjectSet &candidates, DataObjectSet &verified, const DataObjectSet &context) const

bool	hasUpdateMethod () const

bool	isArray () const

bool	isContainer () const

bool	isDataModel () const

bool	isMatrix () const

bool	isNameVector () const

bool	isNonUniqueName () const

virtual bool	isPrerequisiteForContext (const CObjectInterface *pObject, const CMath::SimulationContextFlag &context, const CObjectInterface::ObjectSet &changedObjects) const

bool	isReference () const

bool	isRoot () const

bool	isSeparator () const

bool	isStaticString () const

bool	isValueBool () const

bool	isValueDbl () const

bool	isValueInt () const

bool	isValueInt64 () const

bool	isValueString () const

bool	isVector () const

virtual bool	mustBeDeleted (const DataObjectSet &deletedObjects) const

void	removeDirectDependency (const CCopasiObject *pObject)

void	setDirectDependencies (const DataObjectSet &directDependencies)

bool	setObjectName (const std::string &name)

virtual bool	setObjectParent (const CCopasiContainer *pParent)

void	setObjectValue (const C_FLOAT64 &value)

void	setObjectValue (const C_INT32 &value)

void	setObjectValue (const bool &value)

template<class CType >
void	setRefresh (CType pType, void(CType::method)(void))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const C_FLOAT64 &))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const C_INT32 &))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const bool &))

virtual	~CCopasiObject ()

Public Member Functions inherited from CObjectInterface
	CObjectInterface ()

virtual	~CObjectInterface ()

Protected Member Functions
void	calculateAmu (size_t rIndex)

C_INT32	checkModel (CModel *model)

	CHybridMethodODE45 (const CCopasiContainer *pParent=NULL)

void	cleanup ()

void	doInverseInterpolation ()

C_FLOAT64	doSingleStep (C_FLOAT64 currentTime, C_FLOAT64 endTime)

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

void	fireReaction (size_t rIndex)

void	fireSlowReaction4Hybrid ()

C_FLOAT64	generateReactionTime (size_t rIndex)

C_FLOAT64	getDefaultAtol (const CModel *pModel) const

size_t	getReactionIndex4Hybrid ()

void	getStochTimeAndIndex (C_FLOAT64 &ds, size_t &rIndex)

void	initMethod (C_FLOAT64 time)

void	integrateDeterministicPart (C_FLOAT64 ds)

void	outputData ()

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

void	outputState (const CState *pS)

void	setupBalances ()

void	setupCalculateSet ()

void	setupDependencyGraph ()

void	setupMetab2React ()

void	setupMetabFlags ()

void	setupMethod ()

void	setupPriorityQueue (C_FLOAT64 startTime=0.0)

void	setupReactAffect ()

void	setupReactionFlags ()

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 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
CVector< C_FLOAT64 >	mAmu

CVector< C_FLOAT64 >	mAmuOld

C_FLOAT64	mAtol

std::set< size_t >	mCalculateSet

Data	mData

bool	mDefaultAtol

CDependencyGraph	mDG

bool	mDoCorrection

CVector< C_FLOAT64 >	mDWork

C_FLOAT64	mEndTime

std::ostringstream	mErrorMsg

size_t	mFirstMetabIndex

bool	mHasAssignments

bool	mHasDetermReaction

bool	mHasStoiReaction

C_INT	mIFlag

CVector< C_INT >	mIWork

std::vector< std::vector < CHybridODE45Balance > >	mLocalBalances

std::vector< std::vector < CHybridODE45Balance > >	mLocalSubstrates

size_t	mMaxBalance

size_t	mMaxSteps

bool	mMaxStepsReached

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

std::vector < CHybridODE45MetabFlag >	mMetabFlags

size_t	mMethod

size_t	mNumReactions

size_t	mNumVariableMetabs

C_INT	mODE45Status

C_FLOAT64	mOldTime

C_FLOAT64 *	mOldY

size_t	mOutputCounter

std::ofstream	mOutputFile

std::string	mOutputFileName

CStateRecord *	mpEventState

CInterpolation *	mpInterpolation

CCopasiVector< CMetab > *	mpMetabolites

CModel *	mpModel

CIndexedPriorityQueue	mPQ

CRandom *	mpRandomGenerator

const CCopasiVectorNS < CReaction > *	mpReactions

CState *	mpState

unsigned C_INT32	mRandomSeed

std::vector< std::set< size_t > >	mReactAffect

CVector< size_t >	mReactionFlags

bool	mReducedModel

C_FLOAT64	mRtol

CVector< C_FLOAT64 >	mStateRecord

CMatrix< C_FLOAT64 >	mStoi

C_FLOAT64	mTime

std::set< size_t >	mUpdateSet

bool	mUseRandomSeed

bool	mUseStateRecord

C_FLOAT64 *	mY

CVector< C_FLOAT64 >	mYdot

CVector< 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
C_INT	fehl_ (pEvalF f, const C_INT neqn, double y, double t, double h__, double yp, double f1, double f2, double f3, double f4, double f5, double s, double yrcd)

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

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

void	initializeParameter ()

C_INT	rkf45_ (pEvalF f, const C_INT neqn, double y, double t, double tout, double relerr, double abserr, C_INT iflag, double work, C_INT iwork, double yrcd)

C_INT	rkfs_ (pEvalF f, const C_INT neqn, double y, double t, double tout, double relerr, double abserr, C_INT iflag, double yp, double h__, double f1, double f2, double f3, double f4, double f5, double savre, double savae, C_INT nfe, C_INT kop, C_INT init, C_INT jflag, C_INT kflag, double yrcd)

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

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

Public Types inherited from CCopasiMethod
enum	SubType { unset = 0, RandomSearch, RandomSearchMaster, SimulatedAnnealing, CoranaWalk, DifferentialEvolution, ScatterSearch, GeneticAlgorithm, EvolutionaryProgram, SteepestDescent, HybridGASA, GeneticAlgorithmSR, HookeJeeves, LevenbergMarquardt, NelderMead, SRES, Statistics, ParticleSwarm, Praxis, TruncatedNewton, Newton, deterministic, LSODAR, directMethod, stochastic, tauLeap, adaptiveSA, hybrid, hybridLSODA, hybridODE45, DsaLsodar, tssILDM, tssILDMModified, tssCSP, mcaMethodReder, scanMethod, lyapWolf, sensMethod, EFMAlgorithm, EFMBitPatternTreeAlgorithm, EFMBitPatternAlgorithm, Householder, crossSectionMethod, linearNoiseApproximation }

Public Types inherited from CCopasiParameterGroup
typedef parameterGroup::iterator	index_iterator

typedef CCopasiContainer::objectMap::iterator	name_iterator

typedef std::vector < CCopasiParameter * >	parameterGroup

Public Types inherited from CCopasiParameter
enum	Type { DOUBLE = 0, UDOUBLE, INT, UINT, BOOL, GROUP, STRING, CN, KEY, FILE, EXPRESSION, INVALID }

Public Types inherited from CCopasiContainer
typedef std::multimap < std::string, CCopasiObject * >	objectMap

Public Types inherited from CCopasiObject
typedef std::set< const CCopasiObject * >	DataObjectSet

typedef std::vector< Refresh * >	DataUpdateSequence

Public Types inherited from CObjectInterface
typedef std::set< const CObjectInterface * >	ObjectSet

typedef std::vector < CObjectInterface * >	UpdateSequence

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

Static Public Member Functions inherited from CCopasiObject
static std::vector< Refresh * >	buildUpdateSequence (const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())

static void	setRenameHandler (CRenameHandler *rh)

Static Public Attributes inherited from CCopasiMethod
static const std::string	SubTypeName []

static const char *	XMLSubType []

Static Public Attributes inherited from CCopasiParameter
static const std::string	TypeName []

static const char *	XMLType []

Static Public Attributes inherited from CCopasiContainer
static const std::vector < CCopasiContainer * >	EmptyList

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

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

Detailed Description

Definition at line 133 of file CHybridMethodODE45.h.

Constructor & Destructor Documentation

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

protected

Default Constructor

CHybridMethodODE45

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

File name: CHybridMethodODE45.cpp Author: Shuo Wang Email: shuow.nosp@m.ang..nosp@m.learn.nosp@m.er@g.nosp@m.mail..nosp@m.com

Last change: 26, Aug 2013 Default constructor.

Definition at line 61 of file CHybridMethodODE45.cpp.

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

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

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

Copy constructor

Parameters

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

Copy Constructor

Definition at line 74 of file CHybridMethodODE45.cpp.

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

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

CHybridMethodODE45::~CHybridMethodODE45 ( )

Destructor.

Definition at line 88 of file CHybridMethodODE45.cpp.

References cleanup(), and DESTRUCTOR_TRACE.

 {
   cleanup();
   DESTRUCTOR_TRACE;
 }

Member Function Documentation

void CHybridMethodODE45::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 1154 of file CHybridMethodODE45.cpp.

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

Referenced by evalF(), fireSlowReaction4Hybrid(), initMethod(), 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<CHybridODE45Balance> & 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
 }

C_INT32 CHybridMethodODE45::checkModel ( CModel * model )

protected

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 2440 of file CHybridMethodODE45.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 CHybridMethodODE45::cleanup ( )

protected

Cleans up memory, etc.

Definition at line 380 of file CHybridMethodODE45.cpp.

References mOldY, mpInterpolation, mpModel, mpRandomGenerator, mpState, and mY.

Referenced by ~CHybridMethodODE45().

 {
   delete mpRandomGenerator;
   mpRandomGenerator = NULL;
   mpModel = NULL;
 
   delete mpState;
   mpState = NULL;
 
   delete mpInterpolation;
   mpInterpolation = NULL;
 
   delete [] mY;
   mY = NULL;
 
   delete [] mOldY;
   mOldY = NULL;
 
   return;
 }

void CHybridMethodODE45::doInverseInterpolation ( )

protected

Do inverse interpolation to find the state when a slow reaction is fired.

Definition at line 977 of file CHybridMethodODE45.cpp.

References CVectorCore< CType >::array(), CState::beginIndependent(), C_FLOAT64, CHybridMethodODE45::Data::dim, CStateRecord::getArray(), CInterpolation::getInterpolationState(), CStateRecord::getTime(), INTERP_RECORD_NUM, mData, mOldTime, mOldY, mpEventState, mpInterpolation, mpModel, mpState, mStateRecord, mTime, mUseStateRecord, mY, CInterpolation::recordReset(), CInterpolation::recordState(), CModel::setState(), CState::setTime(), and CModel::updateSimulatedValues().

Referenced by doSingleStep().

 {
 
   //==(1)==for one-step method, reset record each time when do interpolation
   mpInterpolation->recordReset();
 
   //==(2)==set record in class interpolation
   mpInterpolation->recordState(mOldTime, mOldY);
   size_t offset;
 
   if (mUseStateRecord)
     {
       for (size_t i = 0; i < INTERP_RECORD_NUM - 2; i++) //record the middle 4 states
         {
           offset = i * (mData.dim + 1);
           mpInterpolation->recordState(mStateRecord[offset],
                                        mStateRecord.array() + offset + 1);
         }
     }
 
   mpInterpolation->recordState(mTime, mY);
 
   //==(3)==do interpolation
   mpEventState = mpInterpolation->getInterpolationState();
 
   //==(4)==record the state
   mTime = mpEventState->getTime();
   C_FLOAT64 * tmpY   = mpEventState->getArray() + 1;
   C_FLOAT64 * stateY = mpState->beginIndependent();
 
   for (size_t i = 0; i < mData.dim - 1; i++)
     stateY[i] = tmpY[i]; //write result into mpState
 
   mpState->setTime(mTime);
   mpModel->setState(*mpState);
   mpModel->updateSimulatedValues(false); //for assignments?????????
 
   return;
 }

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

protected

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

Simulates the system over the next interval of time. The new time after this step is returned.

Parameters

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

Returns: A C_FLOAT giving the new time

Definition at line 753 of file CHybridMethodODE45.cpp.

References C_FLOAT64, CONTINUE, DETERMINISTIC, doInverseInterpolation(), fireReaction(), fireSlowReaction4Hybrid(), CModel::getState(), getStochTimeAndIndex(), CState::getTime(), HAS_EVENT, HYBRID, integrateDeterministicPart(), mMethod, mODE45Status, mpModel, mpState, NEW_STEP, CModel::setState(), CState::setTime(), STOCHASTIC, and updatePriorityQueue().

Referenced by step().

 {
   C_FLOAT64 ds = 0.0;
   size_t rIndex = 0;
 
   //1----pure SSA method
   if (mMethod == STOCHASTIC) //has only stochastic reactions
     {
       getStochTimeAndIndex(ds, rIndex);
 
       if (ds > endTime)
         ds = endTime;
       else
         {
           fireReaction(rIndex);
           updatePriorityQueue(rIndex, ds);
         }
 
       //population has been changed and record to the mCurrentState
       //in fireReaction(). So here just store time
       *mpState = mpModel->getState();
       mpState->setTime(ds);
       mpModel->setState(*mpState);
 
       //outputState(mpCurrentState);
       //getchar();
     }
   //2----Method with Deterministic Part
   else if (mMethod == DETERMINISTIC) //has only deterministic reactions
     {
       integrateDeterministicPart(endTime - currentTime);
       ds = mpState->getTime();
       mODE45Status = CONTINUE;
       // Till now, state has been recorded
     }
   else if (mMethod == HYBRID)//Hybrid Method
     {
       integrateDeterministicPart(endTime - currentTime);
       ds = mpState->getTime();
 
       if (mODE45Status == HAS_EVENT) //fire slow reaction
         {
           doInverseInterpolation();
 
           fireSlowReaction4Hybrid();
           mODE45Status = NEW_STEP;
         }
       else
         mODE45Status = CONTINUE;
     }
 
   return ds;
 }

bool CHybridMethodODE45::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 96 of file CHybridMethodODE45.cpp.

References initializeParameter().

 {
   initializeParameter();
   return true;
 }

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

staticprotected

Dummy Function for calculating derivative of ODE systems

Dummy f function for calculating derivative of y

Definition at line 1020 of file CHybridMethodODE45.cpp.

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

Referenced by integrateDeterministicPart().

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

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

protected

This evaluates the derivatives for the complete model

Derivative Calculation Function

Definition at line 1026 of file CHybridMethodODE45.cpp.

References CState::beginIndependent(), C_FLOAT64, calculateAmu(), CHybridMethodODE45::Data::dim, mAmu, mCalculateSet, mData, mLocalBalances, mNumReactions, mpModel, mpState, mReactionFlags, CModel::setState(), CState::setTime(), SLOW, and CModel::updateSimulatedValues().

Referenced by EvalF().

 {
   size_t i;
 
   //========No need to record State, right?========
   //maybe here I just want the propensities amu???
   //calculateDerivative(temp);
 
   //1====calculate propensities
   //just care about reactions involving fast metab
 
   //(1)Put current time *t and independent values *y into model.
   // This step seemes necessary, since propensity calculation process
   // requires functions called from the model class.
   mpState->setTime(*t);
   C_FLOAT64 * tmpY = mpState->beginIndependent();//mpState is a local copy
 
   for (i = 0; i < mData.dim - 1; ++i)
     tmpY[i] = y[i];
 
   mpModel->setState(*mpState);
   mpModel->updateSimulatedValues(false); //really?
 
   //(2)Calculate propensities.
   std::set <size_t>::iterator reactIt = mCalculateSet.begin();
   size_t reactID;
 
   for (; reactIt != mCalculateSet.end(); reactIt++)
     {
       reactID = *reactIt;
       calculateAmu(reactID);
     }
 
   //2====calculate derivative
   //(1)initialize
   for (i = 0; i < mData.dim; i++)
     ydot[i] = 0.0;
 
   //(2)go through all the reactions and
   //update derivatives
   std::vector <CHybridODE45Balance>::iterator metabIt;
   std::vector <CHybridODE45Balance>::iterator metabEndIt;
   size_t metabIndex;
 
   for (i = 0; i < mNumReactions; i++)
     {
       if (mReactionFlags[i] == SLOW) //slow reaction
         ydot[mData.dim - 1] += mAmu[i];
       else //fast reaction
         {
           metabIt    = mLocalBalances[i].begin();
           metabEndIt = mLocalBalances[i].end();
 
           for (; metabIt != metabEndIt; metabIt++)
             {
               metabIndex = metabIt->mIndex;
               ydot[metabIndex] += metabIt->mMultiplicity * mAmu[i];
             }
         }
     }
 
   return;
 }

C_INT CHybridMethodODE45::fehl_	(	pEvalF	f,
		const C_INT *	neqn,
		double *	y,
		double *	t,
		double *	h__,
		double *	yp,
		double *	f1,
		double *	f2,
		double *	f3,
		double *	f4,
		double *	f5,
		double *	s,
		double *	yrcd
	)

private

Definition at line 2266 of file CHybridMethodODE45.cpp.

References C_INT.

Referenced by rkfs_().

 {
   /* System generated locals */
   C_INT i__1;
   double d__1;
 
   /* Local variables */
   static C_INT k;
   static double ch;
   static C_INT base;
 
   /* fehlberg fourth-fifth order runge-kutta method */
 
   /* fehl integrates a system of neqn first order */
   /* ordinary differential equations of the form */
   /*             dy(i)/dt=f(t,y(1),---,y(neqn)) */
   /* where the initial values y(i) and the initial derivatives */
   /* yp(i) are specified at the starting point t. fehl advances */
   /* the solution over the fixed step h and returns */
   /* the fifth order (sixth order accurate locally) solution */
   /* approximation at t+h in array s(i). */
   /* f1,---,f5 are arrays of dimension neqn which are needed */
   /* for internal storage. */
   /* the formulas have been grouped to control loss of significance. */
   /* fehl should be called with an h not smaller than 13 units of */
   /* roundoff in t so that the various independent arguments can be */
   /* distinguished. */
 
   /* Parameter adjustments */
   --yrcd;
   --s;
   --f5;
   --f4;
   --f3;
   --f2;
   --f1;
   --yp;
   --y;
 
   /* Function Body */
 
   //~~~~~~~~(1)~~~~~~~~
   ch = *h__ / 4.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       f5[k] = y[k] + ch * yp[k];
     }
 
   d__1 = *t + ch;
   (*f)(neqn, &d__1, &f5[1], &f1[1]);
 
   /* Record trcd and yrcd at t+h/4 */
   base = 1;
   yrcd[base] = *t + ch;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       yrcd[base + k] = f5[k];
     }
 
   //~~~~~~~~(2)~~~~~~~~
   ch = *h__ * 3. / 32.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       f5[k] = y[k] + ch * (yp[k] + f1[k] * 3.);
     }
 
   d__1 = *t + *h__ * 3. / 8.;
   (*f)(neqn, &d__1, &f5[1], &f2[1]);
 
   /* Record trcd and yrcd at t+3h/8 */
   base = *neqn + 2;
   yrcd[base] = *t + *h__ * 3. / 8.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       yrcd[base + k] = f5[k];
     }
 
   //~~~~~~~~(3)~~~~~~~~
   ch = *h__ / 2197.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       f5[k] = y[k] + ch * (yp[k] * 1932. + (f2[k] * 7296.
                                             - f1[k] * 7200.));
     }
 
   d__1 = *t + *h__ * 12. / 13.;
   (*f)(neqn, &d__1, &f5[1], &f3[1]);
 
   /* Record trcd and yrcd at t+12h/13 */
   base = (*neqn + 1) * 3 + 1;
   yrcd[base] = *t + *h__ * 12. / 13.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       yrcd[base + k] = f5[k];
     }
 
   //~~~~~~~~(4)~~~~~~~~
   ch = *h__ / 4104.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       f5[k] = y[k] + ch * (yp[k] * 8341. - f3[k] * 845.
                            + (f2[k] * 29440. - f1[k] * 32832.));
     }
 
   d__1 = *t + *h__;
   (*f)(neqn, &d__1, &f5[1], &f4[1]);
 
   //~~~~~~~~(5)~~~~~~~~
   ch = *h__ / 20520.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       f1[k] = y[k] + ch * (yp[k] * -6080.
                            + (f3[k] * 9295. - f4[k] * 5643.)
                            + (f1[k] * 41040. - f2[k] * 28352.));
     }
 
   d__1 = *t + *h__ / 2.;
   (*f)(neqn, &d__1, &f1[1], &f5[1]);
 
   /* Record trcd and yrcd at t+h/2 */
   //base = (*neqn + 1 << 1) + 1;
   base = (*neqn + 1) * 2 + 1;
   yrcd[base] = *t + *h__ / 2.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       yrcd[base + k] = f1[k];
     }
 
   //~~~~~~~~(6)~~~~~~~~
   /* compute approximate solution at t+h */
   ch = *h__ / 7618050.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       s[k] = y[k] + ch * (yp[k] * 902880.
                           + (f3[k] * 3855735. - f4[k] * 1371249.)
                           + (f2[k] * 3953664. + f5[k] * 277020.));
     }
 
   return 0;
 } /* fehl_ */

void CHybridMethodODE45::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.

Definition at line 1410 of file CHybridMethodODE45.cpp.

References C_FLOAT64, CModelEntity::getValue(), mLocalBalances, CMetab::refreshConcentration(), and CModelEntity::setValue().

Referenced by doSingleStep(), and fireSlowReaction4Hybrid().

 {
   // 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();
     }
 
   return;
 }

void CHybridMethodODE45::fireSlowReaction4Hybrid ( )

protected

Fire slow reaction and update populations and propensities when Hybrid Method is used

Definition at line 807 of file CHybridMethodODE45.cpp.

References calculateAmu(), fireReaction(), getReactionIndex4Hybrid(), CModel::getState(), mCalculateSet, mpModel, mpState, and mReactAffect.

Referenced by doSingleStep().

 {
   //First Update Propensity
   std::set <size_t>::iterator reactIt = mCalculateSet.begin();
   size_t reactID;
 
   for (; reactIt != mCalculateSet.end(); reactIt++)
     {
       reactID = *reactIt;
       calculateAmu(reactID);
     }
 
   reactID = getReactionIndex4Hybrid();
   fireReaction(reactID); //Directly update current state in global view
   *mpState = mpModel->getState();
 
   //update corresponding propensities
   std::set <size_t>::iterator updateIt
     = mReactAffect[reactID].begin();
   const std::set <size_t>::iterator updateEndIt
     = mReactAffect[reactID].end();
 
   for (; updateIt != updateEndIt; updateIt++)
     {
       reactID = *updateIt;
       calculateAmu(reactID);
     }
 
   return;
 }

C_FLOAT64 CHybridMethodODE45::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 1107 of file CHybridMethodODE45.cpp.

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

Referenced by 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 > * CHybridMethodODE45::getAffects ( size_t rIndex )

private

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 1266 of file CHybridMethodODE45.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 CHybridMethodODE45::getDefaultAtol ( const CModel * pModel ) const

protected

Calculate the default absolute tolerance

Parameters

const CModel * pModel

Returns: C_FLOAT64 defaultAtol

Definition at line 839 of file CHybridMethodODE45.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 > * CHybridMethodODE45::getDependsOn ( size_t rIndex )

private

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 1235 of file CHybridMethodODE45.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;
 
       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;
 }

size_t CHybridMethodODE45::getReactionIndex4Hybrid ( )

protected

Function return a reaction index which is a slow reaction firing at the event time.

Definition at line 1441 of file CHybridMethodODE45.cpp.

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

Referenced by fireSlowReaction4Hybrid().

 {
   //calculate sum of amu
   C_FLOAT64 mAmuSum = 0.0;
 
   for (int i = 0; i < mNumReactions; i++)
     {
       if (mReactionFlags[i] == SLOW)
         mAmuSum += mAmu[i];
     }
 
   //get the threshold
   C_FLOAT64 rand2 = mpRandomGenerator->getRandomOO();
   C_FLOAT64 threshold = mAmuSum * rand2;
 
   //get the reaction index
   size_t rIndex;
   C_FLOAT64 tmp = 0.0;
 
   //is there some algorithm that can get a log() complex?
   for (rIndex = 0; rIndex < mNumReactions; rIndex++)
     {
       if (mReactionFlags[rIndex] == SLOW)
         {
           tmp += mAmu[rIndex];
 
           if (tmp >= threshold)
             return rIndex;
         }
     }
   return rIndex;
 }

void CHybridMethodODE45::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 1100 of file CHybridMethodODE45.cpp.

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

Referenced by doSingleStep().

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

void CHybridMethodODE45::initializeParameter ( )

private

Intialize the method parameter

Definition at line 175 of file CHybridMethodODE45.cpp.

References CCopasiParameterGroup::addParameter(), CCopasiParameterGroup::assertParameter(), CCopasiParameter::BOOL, C_FLOAT64, C_INT32, CCopasiParameterGroup::getParameter(), CCopasiParameter::getValue(), MAX_STEPS, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pUINT, RANDOM_SEED, CCopasiParameterGroup::removeParameter(), CCopasiParameterGroup::setValue(), CCopasiParameter::UDOUBLE, CCopasiParameter::UINT, and USE_RANDOM_SEED.

Referenced by CHybridMethodODE45(), and elevateChildren().

 {
   CCopasiParameter *pParm;
 
   //-----------------------------------------------------------------------
   // This part will be dealt with later, when considering the partition strategy.
   //Now, under the condition of using statistic partition given by users,
   //we won't adopt such parameters listed here.
   assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) MAX_STEPS);
   //assertParameter("Lower Limit", CCopasiParameter::DOUBLE, (C_FAT64) 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)//need this part?
         {
           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("");
         }
     }
 
   /* ODE45 ****************************************************************************/
   //----------------------------------------------------------------------
   //addParameter("Partitioning Stepsize",
   //             CCopasiParameter::DOUBLE, (C_FLOAT64) PARTITIONING_STEPSIZE);
 
   //addParameter("Max Internal Steps (LSOA)",
   //             CCopasiParameter::UINT, (unsigned C_INT32) 10000);
   //-----------------------------------------------------------------------
 
   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);
 
   //????????????????????????????????????????????????????????
   // These parameters are no longer supported.
   //removeParameter("Adams Max Order");
   //removeParameter("BDF Max Order");
   //????????????????????????????????????????????????????????
 }

void CHybridMethodODE45::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 295 of file CHybridMethodODE45.cpp.

References C_FLOAT64, C_INT, calculateAmu(), CHybridMethodODE45::Data::dim, getDefaultAtol(), CModel::getMetabolitesX(), CCopasiProblem::getModel(), CModel::getNumDependentReactionMetabs(), CModel::getNumIndependentReactionMetabs(), CModel::getReactions(), CModel::getStoi(), CCopasiParameter::getValue(), CRandom::initialize(), INTERP_RECORD_NUM, mAmu, mAmuOld, mAtol, mData, mDefaultAtol, mDWork, mIWork, mMaxSteps, mMaxStepsReached, mNumReactions, mNumVariableMetabs, mODE45Status, mOldY, mpInterpolation, mpMetabolites, mpModel, CTrajectoryMethod::mpProblem, mpRandomGenerator, mpReactions, mRandomSeed, mRtol, mStateRecord, mStoi, mUpdateSet, mUseRandomSeed, mY, mYdot, NEW_STEP, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pUDOUBLE, CCopasiParameter::Value::pUINT, CVector< CType >::resize(), setupBalances(), setupCalculateSet(), setupDependencyGraph(), setupMetab2React(), setupMetabFlags(), setupMethod(), setupPriorityQueue(), setupReactAffect(), setupReactionFlags(), CCopasiParameterGroup::setValue(), CCopasiVector< T >::size(), and temp.

Referenced by start().

 {
   //(1)----set attributes related with REACTIONS
   mpReactions   = &mpModel->getReactions();
   mNumReactions = mpReactions->size();
 
   mAmu.resize(mNumReactions);
   mAmuOld.resize(mNumReactions);
 
   setupReactionFlags();
 
   //(2)----set attributes related with METABS
   mpMetabolites = &(const_cast < CCopasiVector < CMetab > & >(mpModel->getMetabolitesX()));
 
   mNumVariableMetabs = mpModel->getNumIndependentReactionMetabs()
                        + mpModel->getNumDependentReactionMetabs();
 
   setupBalances();    //initialize mLocalBalances and mLocalSubstrates
   setupMetab2React(); //initialize mMetab2React
   setupMetabFlags();
 
   //(3)----set attributes related with STATE
   temp.resize(mNumVariableMetabs + 1);
 
   //(4)----set attributes related with SYSTEM
   mMaxSteps = * getValue("Max Internal Steps").pUINT;
   mMaxStepsReached = false;
   setupMethod();
 
   //(5)----set attributes related with STOCHASTIC
   mUseRandomSeed = * getValue("Use Random Seed").pBOOL;
   mRandomSeed = * getValue("Random Seed").pUINT;
 
   if (mUseRandomSeed)
     mpRandomGenerator->initialize(mRandomSeed);
 
   mStoi = mpModel->getStoi();
   mUpdateSet.clear();// how to set this parameter
   setupCalculateSet(); //should be done after setupBalances()
   setupReactAffect();
 
   setupDependencyGraph(); // initialize mDG
   setupPriorityQueue(start_time); // initialize mPQ
 
   //(6)----set attributes for INTERPOLATION
   mpInterpolation = new CInterpolation(INTERP_RECORD_NUM, mNumVariableMetabs);
   mStateRecord.resize((INTERP_RECORD_NUM - 2) * (2 + mNumVariableMetabs));
 
   //(7)----set attributes for ODE45
   mODE45Status = NEW_STEP;
   mData.dim = (C_INT)(mNumVariableMetabs + 1);  //one more for sum of propensities
   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(3 + 6 * mData.dim);
   mIWork.resize(5);
 
   //setupUpdateSet();
 
   // we don't want to directly record new results into mpState, since
   // the sum of slow reaction propensities is also recorded in mY
   mY    = new C_FLOAT64[mData.dim];
   mOldY = new C_FLOAT64[mData.dim];
 
   //first calculate propensities, in the next integration process
   //system will just update the propensities record in mCalculateSet
   for (size_t i = 0; i < mNumReactions; i++)
     calculateAmu(i);
 
   return;
 }

void CHybridMethodODE45::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 870 of file CHybridMethodODE45.cpp.

References CVectorCore< CType >::array(), CState::beginIndependent(), C_FLOAT64, CONTINUE, CHybridMethodODE45::Data::dim, EvalF(), CCopasiMessage::EXCEPTION, CRandom::getRandomOO(), CModel::getState(), CState::getTime(), HAS_ERR, HAS_EVENT, mAtol, max, MCTrajectoryMethod, mData, mDWork, mErrorMsg, mIFlag, mIWork, mODE45Status, mOldTime, mOldY, mpModel, mpRandomGenerator, mpState, mRtol, mStateRecord, mTime, mY, NEW_STEP, REACH_END_TIME, rkf45_(), CModel::setState(), CState::setTime(), CModel::setTime(), and CModel::updateSimulatedValues().

Referenced by doSingleStep().

 {
 
   //1----Set Parameters for ODE45 solver
   //=(1)= solver error message
   mErrorMsg.str("");
 
   //=(2)= set state
   *mpState = mpModel->getState();
 
   //=(3)= set time and old time
   mTime = mpState->getTime();
   C_FLOAT64 EndTime = mTime + deltaT;
 
   mOldTime = mTime;
 
   //=(4)= set y and old_y
   C_FLOAT64 * stateY = mpState->beginIndependent();
 
   for (size_t i = 0; i < mData.dim - 1; i++, stateY++)
     {
       mOldY[i] = *stateY;
       mY[i]    = *stateY;
     }
 
   if (mODE45Status == NEW_STEP)
     {
       C_FLOAT64 rand2 = mpRandomGenerator->getRandomOO();
       mY[mData.dim - 1] = log(rand2);
     }
 
   mOldY[mData.dim - 1] = mY[mData.dim - 1];
 
   //2----Reset ODE45 Solver
   mIFlag = -1; //integration by one step
 
   //3----If time increment is too small, do nothing
   C_FLOAT64 tdist , d__1, d__2, w0;
   tdist = fabs(deltaT); //absolute time increment
   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
     {
       mpModel->setTime(EndTime);
       return;
     }
 
   //4----just do nothing if there are no variables
   if (!mData.dim)
     {
       mpModel->setTime(EndTime);
       return;
     }
 
   while ((mODE45Status == CONTINUE)
          || (mODE45Status == NEW_STEP))
     {
       //(1)----ODE Solver
       rkf45_(&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
              &mRtol ,                 //6. relative tolerance array
              &mAtol ,                 //7. absolute tolerance array
              &mIFlag ,                //8. the state for input and output
              mDWork.array() ,         //9. the double work array
              mIWork.array() ,         //10. the int work array
              mStateRecord.array());   //11. the array to record temp state
 
       //(2)----set mODE45Status
       if (mIFlag == 2)
         {
           if (mY[mData.dim - 1] >= 0.0)
             mODE45Status = HAS_EVENT;
           else
             mODE45Status = REACH_END_TIME; //reach EndTime
         }
       else if (mIFlag == -2) //finish one step integration
         {
           if (mY[mData.dim - 1] >= 0.0) //has an event
             mODE45Status = HAS_EVENT;
           else
             mODE45Status = CONTINUE;
         }
       else //error happens
         {
           mODE45Status = HAS_ERR;
           std::cout << "Error Type: " << mIFlag << std::endl;
           CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 6, mErrorMsg.str().c_str());
         }//end if
     }//end while
 
   //5----Record State
   stateY = mpState->beginIndependent();
 
   for (size_t i = 0; i < mData.dim - 1; i++)
     stateY[i] = mY[i]; //write result into mpState
 
   mpState->setTime(mTime);
   mpModel->setState(*mpState);
   mpModel->updateSimulatedValues(false); //for assignments?????????
 
   return;
 }

bool CHybridMethodODE45::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 103 of file CHybridMethodODE45.cpp.

References CCopasiMessage::ERROR, CCopasiMessage::EXCEPTION, CModel::getCompartments(), CTrajectoryProblem::getDuration(), CModel::getEvents(), CModel::getMetabolites(), CCopasiProblem::getModel(), CModel::getModelValues(), CModel::getNumMetabs(), CModel::getNumModelValues(), CTrajectoryMethod::isValidProblem(), MCTrajectoryMethod, CModelEntity::ODE, and CCopasiVector< T >::size().

 {
   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;
     }
 */
   
   //events are not supported at the moment
   if (pTP->getModel()->getEvents().size() > 0)
     {
       CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 23);
       return false;
     }
   
 
   return true;
 }

bool CHybridMethodODE45::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 2869 of file CHybridMethodODE45.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 CHybridMethodODE45::outputData ( )

protected

Prints out data on standard output.

Prints out data on standard output. Deprecated.

Definition at line 2698 of file CHybridMethodODE45.cpp.

References mDoCorrection, mFirstMetabIndex, mHasDetermReaction, mHasStoiReaction, mLocalBalances, mLocalSubstrates, mMetab2React, mMetabFlags, mNumReactions, mNumVariableMetabs, mpMetabolites, mReactAffect, mReactionFlags, mReducedModel, and CCopasiVector< T >::size().

 {
   std::cout << "============Boolean============" << std::endl;
 
   if (mDoCorrection)
     std::cout << "mDoCorrection: Yes" << std::endl;
   else
     std::cout << "mDoCorrection: No" << std::endl;
 
   if (mReducedModel)
     std::cout << "mReducedModel: Yes" << std::endl;
   else
     std::cout << "mReducedModel: No" << std::endl;
 
   std::cout << std::endl;
 
   std::cout << "============Metab============" << std::endl;
 
   std::cout << "~~~~mNumVariableMetabs:~~~~ " << mNumVariableMetabs << std::endl;
   std::cout << "~~~~mFirstMetabIndex:~~~~ " << mFirstMetabIndex << std::endl;
   std::cout << "~~~~mpMetabolites:~~~~" << std::endl;
 
   for (size_t i = 0; i < mpMetabolites->size(); i++)
     {
       std::cout << "metab #" << i << " name: " << (*mpMetabolites)[i]->getObjectDisplayName() << std::endl;
       std::cout << "value pointer: " << (*mpMetabolites)[i]->getValuePointer() << std::endl;
       std::cout << "value: " << *((double *)(*mpMetabolites)[i]->getValuePointer()) << std::endl;
     }
 
   std::cout << std::endl;
 
   std::cout << "~~~~mMetabFlags:~~~~ " << std::endl;
 
   for (size_t i = 0; i < mMetabFlags.size(); ++i)
     {
       std::cout << "metab #" << i << std::endl;
       std::cout << "mFlag: " << mMetabFlags[i].mFlag << std::endl;
       std::cout << "mFastReactions: ";
       std::set<size_t>::iterator it = mMetabFlags[i].mFastReactions.begin();
       std::set<size_t>::iterator endIt = mMetabFlags[i].mFastReactions.end();
 
       for (; it != endIt; it++)
         std::cout << *it << " ";
 
       std::cout << std::endl;
     }
 
   std::cout << std::endl;
 
   std::cout << "============Reaction============" << std::endl;
 
   std::cout << "~~~~mNumReactions:~~~~ " << mNumReactions << std::endl;
 
   for (size_t i = 0; i < mNumReactions; ++i)
     {
       std::cout << "Reaction #: " << i << " Flag: " << mReactionFlags[i] << std::endl;
     }
 
   std::cout << "~~~~mLocalBalances:~~~~ " << std::endl;
 
   for (size_t i = 0; i < mLocalBalances.size(); ++i)
     {
       std::cout << "Reaction: " << i << std::endl;
 
       for (size_t j = 0; j < mLocalBalances[i].size(); ++j)
         {
           std::cout << "Index: " << mLocalBalances[i][j].mIndex << std::endl;
           std::cout << "mMultiplicity: " << mLocalBalances[i][j].mMultiplicity << std::endl;
           std::cout << "mpMetablite: " << mLocalBalances[i][j].mpMetabolite << std::endl;
         }
     }
 
   std::cout << "~~~~mLocalSubstrates:~~~~ " << std::endl;
 
   for (size_t i = 0; i < mLocalSubstrates.size(); ++i)
     {
       std::cout << "Reaction: " << i << std::endl;
 
       for (size_t j = 0; j < mLocalSubstrates[i].size(); ++j)
         {
           std::cout << "Index: " << mLocalSubstrates[i][j].mIndex << std::endl;
           std::cout << "mMultiplicity: " << mLocalSubstrates[i][j].mMultiplicity << std::endl;
           std::cout << "mpMetablite: " << mLocalSubstrates[i][j].mpMetabolite << std::endl;
         }
     }
 
   std::cout << "~~~~mMetab2React:~~~~ " << std::endl;
 
   for (size_t i = 0; i < mMetab2React.size(); i++)
     {
       std::cout << "metab #: " << i << std::endl;
       std::set<size_t>::iterator it = mMetab2React[i].begin();
       std::set<size_t>::iterator endIt = mMetab2React[i].end();
       std::cout << "React: ";
 
       for (; it != endIt; it++)
         std::cout << *it << " ";
 
       std::cout << std::endl;
     }
 
   std::cout << std::endl;
   std::cout << "~~~~mReactAffect:~~~~ " << std::endl;
 
   for (size_t i = 0; i < mReactAffect.size(); i++)
     {
       std::cout << "react #: " << i << std::endl;
       std::set<size_t>::iterator it = mReactAffect[i].begin();
       std::set<size_t>::iterator endIt = mReactAffect[i].end();
       std::cout << "affect: ";
 
       for (; it != endIt; it++)
         std::cout << *it << " ";
 
       std::cout << std::endl;
     }
 
   if (mHasStoiReaction)
     std::cout << "mHasStoiReaction: Yes" << std::endl;
   else
     std::cout << "mHasStoiReaction: No" << std::endl;
 
   if (mHasDetermReaction)
     std::cout << "mHasDetermReaction: Yes" << std::endl;
   else
     std::cout << "mHasDetermReaction: No" << std::endl;
 
   getchar();
   return;
 }

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

protected

Prints out various data on standard output for debugging purposes.

Definition at line 2476 of file CHybridMethodODE45.cpp.

References CCopasiObject::getObjectName(), CState::getTime(), mAmu, mAmuOld, mDG, mLocalBalances, mMaxBalance, mMaxSteps, mMetab2React, mMetabFlags, mNumVariableMetabs, CTrajectoryMethod::mpCurrentState, mpMetabolites, mPQ, mpRandomGenerator, mpReactions, mReactionFlags, mStoi, mUpdateSet, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CCopasiVector< T >::size(), SLOW, 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].mFlag == SLOW)
               os << "SLOW ";
             else
               os << "FAST ";
           }
 
         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].mFlag == SLOW)
               os << "SLOW ";
             else
               os << "FAST ";
           }
 
         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 CHybridMethodODE45::outputState ( const CState * mpState )

protected

Print State Data for Debug

Definition at line 2832 of file CHybridMethodODE45.cpp.

References CState::beginDependent(), CState::beginFixed(), CState::beginIndependent(), C_FLOAT64, CState::endDependent(), CState::endFixed(), CState::endIndependent(), CState::getNumDependent(), CState::getNumIndependent(), and CState::getTime().

 {
   const C_FLOAT64 time = mpState->getTime();
   std::cout << "**State Output**" << std::endl;
   std::cout << "Time: " << time << std::endl;
 
   std::cout << "Indep #: " << mpState->getNumIndependent() << " Id: ";
   const C_FLOAT64 *pIt = mpState->beginIndependent();
   const C_FLOAT64 *pEnd = mpState->endIndependent();
 
   for (; pIt != pEnd; pIt++)
     std::cout << *pIt << " ";
 
   std::cout << std::endl;
 
   std::cout << "Dep #: " << mpState->getNumDependent() << " Id: ";
   pIt = mpState->beginDependent();
   pEnd = mpState->endDependent();
 
   for (; pIt != pEnd; pIt++)
     std::cout << *pIt << " ";
 
   std::cout << std::endl;
 
   std::cout << "Fix #: " << mpState->getNumIndependent() << " Id: ";
   pIt = mpState->beginFixed();
   pEnd = mpState->endFixed();
 
   for (; pIt != pEnd; pIt++)
     std::cout << *pIt << " ";
 
   std::cout << std::endl;
 
   getchar();
   return;
 }

C_INT CHybridMethodODE45::rkf45_	(	pEvalF	f,
		const C_INT *	neqn,
		double *	y,
		double *	t,
		double *	tout,
		double *	relerr,
		double *	abserr,
		C_INT *	iflag,
		double *	work,
		C_INT *	iwork,
		double *	yrcd
	)

private

Definition at line 1475 of file CHybridMethodODE45.cpp.

References C_INT, and rkfs_().

Referenced by integrateDeterministicPart().

 {
   static C_INT k1, k2, k3, k4, k5, k6, k1m;
 
   /*     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
   /*     This is a modefied fehlberg fourth-fifth order runge-kutta method */
   /*     based on h.a.watts and l.f.shapine's fortran ODE solver rkf45.f */
   /*     with one more input double precision vector recording time */
   /*     points between t and t+h in each steap and */
   /*     approximations of y in such corresponding times. */
   /*     yrcd is a vector of length 4*(neqn+1), where "4" is related with */
   /*     four middle time points: */
   /*     c1=t+h/4, c2=t+3h/8, c3=t+h/2, c4=t+12h/13. */
   /*     yrcd can be seperated into four parts, each of which is a double */
   /*     precision vector of length neqn+1. The first member is assigned */
   /*     to time record, and the last neqn members are for y at the */
   /*     corresponding time. */
   /*     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
 
   /*     fehlberg fourth-fifth order runge-kutta method */
 
   /*     written by h.a.watts and l.f.shampine */
   /*                   sandia laboratories */
   /*                  albuquerque,new mexico */
 
   /*    rkf45 is primarily designed to solve non-stiff and mildly stiff */
   /*    differential equations when derivative evaluations are inexpensive. */
   /*    rkf45 should generally not be used when the user is demanding */
   /*    high accuracy. */
 
   /* abstract */
 
   /*    subroutine  rkf45  integrates a system of neqn first order */
   /*    ordinary differential equations of the form */
   /*             dy(i)/dt = f(t,y(1),y(2),...,y(neqn)) */
   /*              where the y(i) are given at t . */
   /*    typically the subroutine is used to integrate from t to tout but it */
   /*    can be used as a one-step integrator to advance the solution a */
   /*    single step in the direction of tout.  on return the parameters in */
   /*    the call list are set for continuing the integration. the user has */
   /*    only to call rkf45 again (and perhaps define a new value for tout). */
   /*    actually, rkf45 is an interfacing routine which calls subroutine */
   /*    rkfs for the solution.  rkfs in turn calls subroutine  fehl which */
   /*    computes an approximate solution over one step. */
 
   /*    rkf45  uses the runge-kutta-fehlberg (4,5)  method described */
   /*    in the reference */
   /*    e.fehlberg , low-order classical runge-kutta formulas with stepsize */
   /*                 control , nasa tr r-315 */
 
   /*    the performance of rkf45 is illustrated in the reference */
   /*    l.f.shampine,h.a.watts,s.davenport, solving non-stiff ordinary */
   /*                 differential equations-the state of the art , */
   /*                 sandia laboratories report sand75-0182 , */
   /*                 to appear in siam review. */
 
   /*    the parameters represent- */
   /*      f -- subroutine f(t,y,yp) to evaluate derivatives yp(i)=dy(i)/dt */
   /*      neqn -- number of equations to be integrated */
   /*      y(*) -- solution vector at t */
   /*      t -- independent variable */
   /*      tout -- output point at which solution is desired */
   /*      relerr,abserr -- relative and absolute error tolerances for local */
   /*            error test. at each step the code requires that */
   /*                 abs(local error) .le. relerr*abs(y) + abserr */
   /*            for each component of the local error and solution vectors */
   /*      iflag -- indicator for status of integration */
   /*      work(*) -- array to hold information internal to rkf45 which is */
   /*            necessary for subsequent calls. must be dimensioned */
   /*            at least  3+6*neqn */
   /*      iwork(*) -- C_INT array used to hold information internal to */
   /*            rkf45 which is necessary for subsequent calls. must be */
   /*            dimensioned at least  5 */
 
   /*  first call to rkf45 */
 
   /*    the user must provide storage in his calling program for the arrays */
   /*    in the call list  -      y(neqn) , work(3+6*neqn) , iwork(5)  , */
   /*    declare f in an external statement, supply subroutine f(t,y,yp) and */
   /*    initialize the following parameters- */
 
   /*      neqn -- number of equations to be integrated.  (neqn .ge. 1) */
   /*      y(*) -- vector of initial conditions */
   /*      t -- starting point of integration , must be a variable */
   /*      tout -- output point at which solution is desired. */
   /*            t=tout is allowed on the first call only, in which case */
   /*            rkf45 returns with iflag=2 if continuation is possible. */
   /*      relerr,abserr -- relative and absolute local error tolerances */
   /*            which must be non-negative. relerr must be a variable while */
   /*            abserr may be a constant. the code should normally not be */
   /*            used with relative error control smaller than about 1.e-8 . */
   /*            to avoid limiting precision difficulties the code requires */
   /*            relerr to be larger than an internally computed relative */
   /*            error parameter which is machine dependent. in particular, */
   /*            pure absolute error is not permitted. if a smaller than */
   /*            allowable value of relerr is attempted, rkf45 increases */
   /*            relerr appropriately and returns control to the user before */
   /*            continuing the integration. */
   /*      iflag -- +1,-1  indicator to initialize the code for each new */
   /*            problem. normal input is +1. the user should set iflag=-1 */
   /*            only when one-step integrator control is essential. in this */
   /*            case, rkf45 attempts to advance the solution a single step */
   /*            in the direction of tout each time it is called. since this */
   /*            mode of operation results in extra computing overhead, it */
   /*            should be avoided unless needed. */
 
   /*  output from rkf45 */
 
   /*      y(*) -- solution at t */
   /*      t -- last point reached in integration. */
   /*      iflag = 2 -- integration reached tout. indicates successful retur */
   /*                   and is the normal mode for continuing integration. */
   /*            =-2 -- a single successful step in the direction of tout */
   /*                   has been taken. normal mode for continuing */
   /*                   integration one step at a time. */
   /*            = 3 -- integration was not completed because relative error */
   /*                   tolerance was too small. relerr has been increased */
   /*                   appropriately for continuing. */
   /*            = 4 -- integration was not completed because more than */
   /*                   3000 derivative evaluations were needed. this */
   /*                   is approximately 500 steps. */
   /*            = 5 -- integration was not completed because solution */
   /*                   vanished making a pure relative error test */
   /*                   impossible. must use non-zero abserr to continue. */
   /*                   using the one-step integration mode for one step */
   /*                   is a good way to proceed. */
   /*            = 6 -- integration was not completed because requested */
   /*                   accuracy could not be achieved using smallest */
   /*                   allowable stepsize. user must increase the error */
   /*                   tolerance before continued integration can be */
   /*                   attempted. */
   /*            = 7 -- it is likely that rkf45 is inefficient for solving */
   /*                   this problem. too much output is restricting the */
   /*                   natural stepsize choice. use the one-step integrator */
   /*                   mode. */
   /*            = 8 -- invalid input parameters */
   /*                   this indicator occurs if any of the following is */
   /*                   satisfied -   neqn .le. 0 */
   /*                                 t=tout  and  iflag .ne. +1 or -1 */
   /*                                 relerr or abserr .lt. 0. */
   /*                                 iflag .eq. 0  or  .lt. -2  or  .gt. 8 */
   /*      work(*),iwork(*) -- information which is usually of no interest */
   /*                   to the user but necessary for subsequent calls. */
   /*                   work(1),...,work(neqn) contain the first derivatives */
   /*                   of the solution vector y at t. work(neqn+1) contains */
   /*                   the stepsize h to be attempted on the next step. */
   /*                   iwork(1) contains the derivative evaluation counter. */
 
   /*  subsequent calls to rkf45 */
 
   /*    subroutine rkf45 returns with all information needed to continue */
   /*    the integration. if the integration reached tout, the user need onl */
   /*    define a new tout and call rkf45 again. in the one-step integrator */
   /*    mode (iflag=-2) the user must keep in mind that each step taken is */
   /*    in the direction of the current tout. upon reaching tout (indicated */
   /*    by changing iflag to 2),the user must then define a new tout and */
   /*    reset iflag to -2 to continue in the one-step integrator mode. */
 
   /*    if the integration was not completed but the user still wants to */
   /*    continue (iflag=3,4 cases), he just calls rkf45 again. with iflag=3 */
   /*    the relerr parameter has been adjusted appropriately for continuing */
   /*    the integration. in the case of iflag=4 the function counter will */
   /*    be reset to 0 and another 3000 function evaluations are allowed. */
 
   /*    however,in the case iflag=5, the user must first alter the error */
   /*    criterion to use a positive value of abserr before integration can */
   /*    proceed. if he does not,execution is terminated. */
 
   /*    also,in the case iflag=6, it is necessary for the user to reset */
   /*    iflag to 2 (or -2 when the one-step integration mode is being used) */
   /*    as well as increasing either abserr,relerr or both before the */
   /*    integration can be continued. if this is not done, execution will */
   /*    be terminated. the occurrence of iflag=6 indicates a trouble spot */
   /*    (solution is changing rapidly,singularity may be present) and it */
   /*    often is inadvisable to continue. */
 
   /*    if iflag=7 is encountered, the user should use the one-step */
   /*    integration mode with the stepsize determined by the code or */
   /*    consider switching to the adams codes de/step,intrp. if the user */
   /*    insists upon continuing the integration with rkf45, he must reset */
   /*    iflag to 2 before calling rkf45 again. otherwise,execution will be */
   /*    terminated. */
 
   /*    if iflag=8 is obtained, integration can not be continued unless */
   /*    the invalid input parameters are corrected. */
 
   /*    it should be noted that the arrays work,iwork contain information */
   /*    required for subsequent integration. accordingly, work and iwork */
   /*    should not be altered. */
 
   /*     compute indices for the splitting of the work array */
 
   /* Parameter adjustments */
   --yrcd;
   --y;
   --work;
   --iwork;
 
   /* Function Body */
   k1m = *neqn + 1;
   //k1m = *neqn;
   k1 = k1m + 1;
   k2 = k1 + *neqn;
   k3 = k2 + *neqn;
   k4 = k3 + *neqn;
   k5 = k4 + *neqn;
   k6 = k5 + *neqn;
 
   /*     this interfacing routine merely relieves the user of a long */
   /*     calling list via the splitting apart of two working storage */
   /*     arrays. if this is not compatible with the users compiler, */
   /*     he must use rkfs directly. */
 
   rkfs_((pEvalF)f, neqn, &y[1], t, tout, relerr, abserr, iflag,
         &work[1], &work[k1m], &work[k1], &work[k2], &work[k3],
         &work[k4], &work[k5], &work[k6], &work[k6 + 1], &iwork[1],
         &iwork[2], &iwork[3], &iwork[4], &iwork[5], &yrcd[1]);
 
   return 0;
 } /* rkf45_ */

C_INT CHybridMethodODE45::rkfs_	(	pEvalF	f,
		const C_INT *	neqn,
		double *	y,
		double *	t,
		double *	tout,
		double *	relerr,
		double *	abserr,
		C_INT *	iflag,
		double *	yp,
		double *	h__,
		double *	f1,
		double *	f2,
		double *	f3,
		double *	f4,
		double *	f5,
		double *	savre,
		double *	savae,
		C_INT *	nfe,
		C_INT *	kop,
		C_INT *	init,
		C_INT *	jflag,
		C_INT *	kflag,
		double *	yrcd
	)

private

Definition at line 1700 of file CHybridMethodODE45.cpp.

References abs, C_INT, fehl_(), max, min, and mUseStateRecord.

Referenced by rkf45_().

 {
   /* Initialized data */
 
   static double remin = 1e-12;
   static C_INT maxnfe = 3000;
 
   /* System generated locals */
   C_INT i__1;
   double d__1, d__2, d__3, d__4;
 
   /* Builtin functions */
   /* Subroutine */
   //int s_stop(char *, ftnlen);
   //double pow_dd(double *, double *);
   //C_INT i_sign(C_INT *, C_INT *);
   //double d_sign(double *, double *);
 
   /* Local variables */
   static double a;
   static C_INT k;
   static double s, ae, ee, dt, et, u26, rer, tol, ypk;
   static double hmin, toln;
   static C_INT mflag;
   static double scale, eeoet;
   static int hfaild;
   static double esttol, twoeps;
   static int output;
 
   /* fehlberg fourth-fifth order runge-kutta method */
 
   /* rkfs integrates a system of first order ordinary differential */
   /* equations as described in the comments for rkf45 . */
   /* the arrays yp,f1,f2,f3,f4,and f5 (of dimension at least neqn) and */
   /* the variables h,savre,savae,nfe,kop,init,jflag,and kflag are used */
   /* internally by the code and appear in the call list to eliminate */
   /* local retention of variables between calls. accordingly, they */
   /* should not be altered. items of possible interest are */
   /*     yp - derivative of solution vector at t */
   /*     h  - an appropriate stepsize to be used for the next step */
   /*     nfe- counter on the number of derivative function evaluations */
 
   /*  remin is the minimum acceptable value of relerr.  attempts */
   /*  to obtain higher accuracy with this subroutine are usually */
   /*  very expensive and often unsuccessful. */
 
   /* Parameter adjustments */
   --yrcd;
   --f5;
   --f4;
   --f3;
   --f2;
   --f1;
   --yp;
   --y;
 
   /* Function Body */
 
   /* the expense is controlled by restricting the number */
   /* of function evaluations to be approximately maxnfe. */
   /* as set, this corresponds to about 500 steps. */
 
   /* here two constants emboding the machine epsilon is present */
   /* twoesp is set to twice the machine epsilon while u26 is set */
   /* to 26 times the machine epsilon */
 
   /* data twoeps, u26/4.4d-16, 5.72d-15/  *** */
   twoeps = 2.f * std::numeric_limits< C_FLOAT64 >::epsilon();
   //twoeps = 4.4e-16;
   u26    = twoeps * 13.f;
 
   /* check input parameters */
   if (*neqn < 1)
     {
       goto L10;
     }
 
   if (*relerr < 0. || *abserr < 0.)
     {
       goto L10;
     }
 
   mflag = abs(*iflag);
 
   if (mflag >= 1 && mflag <= 8)
     {
       goto L20;
     }
 
   /* invalid input */
 L10:
   *iflag = 8;
   return 0;
 
   /* is this the first call */
 L20:
 
   if (mflag == 1)
     {
       goto L50;
     }
 
   /* check continuation possibilities */
 
   if (*t == *tout && *kflag != 3)
     {
       goto L10;
     }
 
   if (mflag != 2)
     {
       goto L25;
     }
 
   /* iflag = +2 or -2 */
   if (*kflag == 3)
     {
       goto L45;
     }
 
   if (*init == 0)
     {
       goto L45;
     }
 
   if (*kflag == 4)
     {
       goto L40;
     }
 
   if (*kflag == 5 && *abserr == 0.)
     {
       goto L30;
     }
 
   if (*kflag == 6 && *relerr <= *savre && *abserr <= *savae)
     {
       goto L30;
     }
 
   goto L50;
 
   /* iflag = 3,4,5,6,7 or 8 */
 L25:
 
   if (*iflag == 3)
     {
       goto L45;
     }
 
   if (*iflag == 4)
     {
       goto L40;
     }
 
   if (*iflag == 5 && *abserr > 0.)
     {
       goto L45;
     }
 
   /* integration cannot be continued since user did not respond to */
   /* the instructions pertaining to iflag=5,6,7 or 8 */
 L30:
   //s_stop("", (ftnlen)0);
   std::cout << "STOP 0 statement executed" << std::endl;
   return -1;
 
   /* reset function evaluation counter */
 L40:
   *nfe = 0;
 
   if (mflag == 2)
     {
       goto L50;
     }
 
   /* reset flag value from previous call */
 L45:
   *iflag = *jflag;
 
   if (*kflag == 3)
     {
       mflag = abs(*iflag);
     }
 
   /* save input iflag and set continuation flag value for subsequent */
   /* input checking */
 L50:
   *jflag = *iflag;
   *kflag = 0;
 
   /* save relerr and abserr for checking input on subsequent calls */
   *savre = *relerr;
   *savae = *abserr;
 
   /* restrict relative error tolerance to be at least as large as */
   /* 2*eps+remin to avoid limiting precision difficulties arising */
   /* from impossible accuracy requests */
 
   rer = twoeps + remin;
 
   if (*relerr >= rer)
     {
       goto L55;
     }
 
   /* relative error tolerance too small */
   *relerr = rer;
   *iflag = 3;
   *kflag = 3;
   return 0;
 
 L55:
   dt = *tout - *t;
 
   if (mflag == 1)
     {
       goto L60;
     }
 
   if (*init == 0)
     {
       goto L65;
     }
 
   goto L80;
 
   /* initialization -- */
   /* set initialization completion indicator,init */
   /* set indicator for too many output points,kop */
   /* evaluate initial derivatives */
   /* set counter for function evaluations,nfe */
   /* evaluate initial derivatives */
   /* set counter for function evaluations,nfe */
   /* estimate starting stepsize */
 
 L60:
   *init = 0;
   *kop = 0;
 
   a = *t;
   (*f)(neqn, &a, &y[1], &yp[1]);
   *nfe = 1;
 
   if (*t != *tout)
     {
       goto L65;
     }
 
   *iflag = 2;
   return 0;
 
 L65:
   *init = 1;
   *h__ = fabs(dt);
   toln = 0.f;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       tol = *relerr * (d__1 = y[k], fabs(d__1)) + *abserr;
 
       if (tol <= 0.f)
         {
           goto L70;
         }
 
       toln = tol;
       ypk = (d__1 = yp[k], fabs(d__1));
       /* Computing 5th power */
       d__1 = *h__, d__2 = d__1, d__1 *= d__1;
 
       if (ypk * (d__2 * (d__1 * d__1)) > tol)
         {
           d__3 = tol / ypk;
           *h__ = pow(d__3, .2);
         }
 
 L70:
       ;
     }
 
   if (toln <= 0.)
     {
       *h__ = 0.;
     }
 
   /* Computing MAX */
   /* Computing MAX */
   d__3 = fabs(*t), d__4 = fabs(dt);
   d__1 = *h__, d__2 = u26 * std::max(d__3, d__4);
   *h__ = std::max(d__1, d__2);
   //*jflag = i_sign(2, iflag);
   *jflag = *iflag >= 0 ? 2 : -2;
 
   /* set stepsize for integration in the direction from t to tout */
 
 L80:
   //*h__ = d_sign(h__, &dt);
   *h__ = dt >= 0 ? fabs(*h__) : -fabs(*h__);
 
   /* test to see if rkf45 is being severely impacted by too many */
   /* output points */
 
   if (fabs(*h__) >= fabs(dt) * 2.)
     {
       ++(*kop);
     }
 
   if (*kop != 100)
     {
       goto L85;
     }
 
   /* unnecessary frequency of output */
   *kop = 0;
   *iflag = 7;
   return 0;
 
 L85:
 
   if (fabs(dt) > u26 * fabs(*t))
     {
       goto L95;
     }
 
   /* if too close to output point,extrapolate and return */
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       /* L90: */
       y[k] += dt * yp[k];
     }
 
   a = *tout;
   (*f)(neqn, &a, &y[1], &yp[1]);
   ++(*nfe);
 
   std::cout << "Step too Small" << std::endl;
   mUseStateRecord = false;
   goto L300;
 
   /* initialize output point indicator */
 
 L95:
   output = 0;
 
   /* to avoid premature underflow in the error tolerance function, */
   /* scale the error tolerances */
 
   scale = 2. / *relerr;
   ae = scale * *abserr;
 
   /* step by step integration */
 
 L100:
   hfaild = 0;
 
   /* set smallest allowable stepsize */
   hmin = u26 * fabs(*t);
 
   /* adjust stepsize if necessary to hit the output point. */
   /* look ahead two steps to avoid drastic changes in the stepsize and */
   /* thus lessen the impact of output points on the code. */
 
   dt = *tout - *t;
 
   if (fabs(dt) >= fabs(*h__) * 2.)
     {
       goto L200;
     }
 
   if (fabs(dt) > fabs(*h__))
     {
       goto L150;
     }
 
   /* the next successful step will complete the integration to the */
   /* output point */
 
   output = 1;
   *h__ = dt;
   goto L200;
 
 L150:
   *h__ = dt * .5;
 
   /* core integrator for taking a single step */
 
   /* the tolerances have been scaled to avoid premature underflow in */
   /* computing the error tolerance function et. */
   /* to avoid problems with zero crossings,relative error is measured */
   /* using the average of the magnitudes of the solution at the */
   /* beginning and end of a step. */
   /* the error estimate formula has been grouped to control loss of */
   /* significance. */
   /* to distinguish the various arguments, h is not permitted */
   /* to become smaller than 26 units of roundoff in t. */
   /* practical limits on the change in the stepsize are enforced to */
   /* smooth the stepsize selection process and to avoid excessive */
   /* chattering on problems having discontinuities. */
   /* to prevent unnecessary failures, the code uses 9/10 the stepsize */
   /* it estimates will succeed. */
   /* after a step failure, the stepsize is not allowed to increase for */
   /* the next attempted step. this makes the code more efficient on */
   /* problems having discontinuities and more effective in general */
   /* since local extrapolation is being used and extra caution seems */
   /* warranted. */
 
   /* test number of derivative function evaluations. */
   /* if okay,try to advance the integration from t to t+h */
 
 L200:
 
   if (*nfe <= maxnfe)
     {
       goto L220;
     }
 
   /* too much work */
   *iflag = 4;
   *kflag = 4;
   return 0;
 
   /* advance an approximate solution over one step of length h */
 
 L220:
   fehl_((pEvalF)f, neqn, &y[1], t, h__, &yp[1], &f1[1], &f2[1],
         &f3[1], &f4[1], &f5[1], &f1[1], &yrcd[1]);
   *nfe += 5;
   mUseStateRecord = true;
 
   /* compute and test allowable tolerances versus local error estimates*/
   /* and remove scaling of tolerances. note that relative error is */
   /* measured with respect to the average of the magnitudes of the */
   /* solution at the beginning and end of the step. */
 
   eeoet = 0.;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       et = (d__1 = y[k], fabs(d__1)) + (d__2 = f1[k], fabs(d__2)) + ae;
 
       if (et > 0.)
         {
           goto L240;
         }
 
       /* inappropriate error tolerance */
       *iflag = 5;
       return 0;
 
 L240:
       ee = (d__1 = yp[k] * -2090. + (f3[k] * 21970. - f4[k] * 15048.) + (f2[k] * 22528. - f5[k] * 27360.), fabs(d__1));
       /* L250: */
       /* Computing MAX */
       d__1 = eeoet, d__2 = ee / et;
       eeoet = std::max(d__1, d__2);
     }
 
   esttol = fabs(*h__) * eeoet * scale / 752400.;
 
   if (esttol <= 1.)
     {
       goto L260;
     }
 
   /* unsuccessful step */
   /* reduce the stepsize , try again */
   /* the decrease is limited to a factor of 1/10 */
 
   hfaild = 1;
   output = 0;
   s = .1;
 
   if (esttol < 59049.)
     {
       s = .9 / pow(esttol, .2);
     }
 
   *h__ = s * *h__;
 
   if (fabs(*h__) > hmin)
     {
       goto L200;
     }
 
   /* requested error unattainable at smallest allowable stepsize */
   *iflag = 6;
   *kflag = 6;
   return 0;
 
   /* successful step */
   /* store solution at t+h */
   /* and evaluate derivatives there */
 
 L260:
   *t += *h__;
   i__1 = *neqn;
 
   for (k = 1; k <= i__1; ++k)
     {
       /* L270: */
       y[k] = f1[k];
     }
 
   a = *t;
   (*f)(neqn, &a, &y[1], &yp[1]);
   ++(*nfe);
 
   /* choose next stepsize */
   /* the increase is limited to a factor of 5 */
   /* if step failure has just occurred, next */
   /* stepsize is not allowed to increase */
 
   s = 5.;
 
   if (esttol > 1.889568e-4)
     {
       s = .9 / pow(esttol, .2);
     }
 
   if (hfaild)
     {
       s = std::min(s, 1.);
     }
 
   /* Computing MAX */
   d__2 = s * fabs(*h__);
   d__1 = std::max(d__2, hmin);
   //*h__ = d_sign(&d__1, h__);
   *h__ = *h__ >= 0 ? fabs(d__1) : -fabs(d__1);
 
   /* end of core integrator */
 
   /* should we take another step */
   if (output)
     {
       goto L300;
     }
 
   if (*iflag > 0)
     {
       goto L100;
     }
 
   /* integration successfully completed */
 
   /* one-step mode */
   *iflag = -2;
   return 0;
 
   /* interval mode */
 L300:
   *t = *tout;
   *iflag = 2;
   return 0;
 } /* rkfs_ */

void CHybridMethodODE45::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 484 of file CHybridMethodODE45.cpp.

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

Referenced by initMethod().

 {
   size_t i, j;
   CHybridODE45Balance newElement;
   size_t maxBalance = 0;
 
   mLocalBalances.clear();
   mLocalBalances.resize(mNumReactions);
   mLocalSubstrates.clear();
   mLocalSubstrates.resize(mNumReactions);
 
   for (i = 0; i < mNumReactions; 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;
 
   return;
 }

void CHybridMethodODE45::setupCalculateSet ( )

protected

Definition at line 632 of file CHybridMethodODE45.cpp.

References FAST, mCalculateSet, mMetab2React, and mMetabFlags.

Referenced by initMethod().

 {
   mCalculateSet.clear();
 
   std::set<size_t>::iterator reactIt;
 
   std::vector<CHybridODE45MetabFlag>::iterator metabIt
     = mMetabFlags.begin();
   const std::vector<CHybridODE45MetabFlag>::iterator metabEndIt
     = mMetabFlags.end();
 
   std::vector< std::set<size_t> >::iterator m2rIt =
     mMetab2React.begin();
 
   for (; metabIt != metabEndIt; ++metabIt, ++m2rIt)
     {
       //for fast metabs, insert reaction into mCalculateSet
       if (metabIt->mFlag == FAST) //fast metab
         {
           for (reactIt  = m2rIt->begin();
                reactIt != m2rIt->end();
                reactIt++) // check all reactions that involves metab
             mCalculateSet.insert(*reactIt);
         }
     }
 
   return;
 }

void CHybridMethodODE45::setupDependencyGraph ( )

protected

Sets up the dependency graph

Sets up the dependency graph.

Definition at line 1286 of file CHybridMethodODE45.cpp.

References CDependencyGraph::addDependent(), CDependencyGraph::clear(), getAffects(), getDependsOn(), mDG, mNumReactions, 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 < mNumReactions; i++)
     {
       // Get the set of metabolites  which affect the value of amu for this
       // reaction i.e. the set on which amu depends. This may be  more than
       // the set of substrates, since the kinetics can involve other
       // reactants, e.g. catalysts. We thus need to step through the
       // rate function and pick out every reactant which can vary.
       DependsOn.push_back(getDependsOn(i));
       // Get the set of metabolites which are affected when this reaction takes place
       Affects.push_back(getAffects(i));
     }
 
   // For each possible pair of reactions i and j, if the intersection of
   // Affects(i) with DependsOn(j) is non-empty, add a dependency edge from i to j.
   for (i = 0; i < mNumReactions; i++)
     {
       for (j = 0; j < mNumReactions; j++)
         {
           // Determine whether the intersection of these two sets is non-empty
           // Could also do this with set_intersection generic algorithm, but that
           // would require operator<() to be defined on the set elements.
 
           std::set<std::string>::iterator iter = Affects[i]->begin();
 
           for (; iter != Affects[i]->end(); iter++)
             {
               if (DependsOn[j]->count(*iter))
                 {
                   // The set intersection is non-empty
                   mDG.addDependent(i, j);
                   break;
                 }
             }
         }
 
       // Ensure that self edges are included
       //mDG.addDependent(i, i);
     }
 
   // 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 CHybridMethodODE45::setupMetab2React ( )

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.

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

Definition at line 563 of file CHybridMethodODE45.cpp.

References CCopasiObject::dependsOn(), CModelEntity::FIXED, CCopasiVector< T >::getIndex(), CModel::getMetabolitesX(), CReaction::getParticleFluxReference(), CModelEntity::getStatus(), CModelEntity::getValueReference(), mMetab2React, mNumReactions, mNumVariableMetabs, mpModel, and CCopasiVector< T >::size().

Referenced by initMethod().

 {
   mMetab2React.clear();
   mMetab2React.resize(mNumVariableMetabs);
 
   CMetab * pMetab;
   CReaction * pReaction;
 
   //Check wether metabs play as substrates and modifiers
   for (size_t rct = 0; rct < mNumReactions; ++rct)
     {
       size_t index;
 
       //deal with substrates
       CCopasiVector <CChemEqElement> metab =
         (*mpReactions)[rct]->getChemEq().getSubstrates();
 
       for (size_t i = 0; i < metab.size(); i++)
         {
           pMetab = const_cast < CMetab* >
                    (metab[i]->getMetabolite());
           index = mpModel->getMetabolitesX().getIndex(pMetab);
 
           if ((pMetab->getStatus()) != CModelEntity::FIXED)
             mMetab2React[index].insert(rct);
         }
 
       //deal with modifiers
       metab = (*mpReactions)[rct]->getChemEq().getModifiers();
 
       for (size_t i = 0; i < metab.size(); i++)
         {
           pMetab = const_cast < CMetab* >
                    (metab[i]->getMetabolite());
           index = mpModel->getMetabolitesX().getIndex(pMetab);
 
           if ((pMetab->getStatus()) != CModelEntity::FIXED)
             mMetab2React[index].insert(rct);
         }
     }
 
   //Check wether metabs appear in reaction laws
   std::set< const CCopasiObject * > changed;
 
   for (size_t metab = 0; metab < mNumVariableMetabs; metab++)
     {
       pMetab = (*mpMetabolites)[metab];
 
       if (pMetab->getStatus() != CModelEntity::FIXED)
         {
           changed.clear();
           changed.insert(mpModel->getValueReference());
           changed.insert(pMetab->getValueReference());
 
           for (size_t react = 0; react < mNumReactions; react++)
             {
               pReaction = (*mpReactions)[react];
 
               if (pReaction->getParticleFluxReference()->dependsOn(changed))
                 mMetab2React[metab].insert(react);
             }
         }
     }
 
   return;
 }

void CHybridMethodODE45::setupMetabFlags ( )

protected

setup mMetabFlags

Definition at line 402 of file CHybridMethodODE45.cpp.

References FAST, mLocalBalances, mMetabFlags, mNumReactions, mNumVariableMetabs, mReactionFlags, and SLOW.

Referenced by initMethod().

 {
   //size_t rctIndex;
   mMetabFlags.resize(mNumVariableMetabs);
 
   std::vector<CHybridODE45MetabFlag>::iterator metabIt
     = mMetabFlags.begin();
   const std::vector<CHybridODE45MetabFlag>::iterator metabEndIt
     = mMetabFlags.end();
 
   // initialization
   for (; metabIt != metabEndIt; ++metabIt)
     {
       metabIt->mFlag = SLOW;
       metabIt->mFastReactions.clear();
     }
 
   std::vector<CHybridODE45Balance>::iterator  itBalance;
   std::vector<CHybridODE45Balance>::iterator itEndBalance;
 
   // go over all mLocalBalances, if balance != 0,
   // and reaction is FAST, insert into set
   for (size_t rct = 0; rct < mNumReactions; ++rct)
     {
       if (mReactionFlags[rct] == FAST)
         {
           itBalance = mLocalBalances[rct].begin();
           itEndBalance = mLocalBalances[rct].end();
 
           for (; itBalance != itEndBalance; ++itBalance)
             {
               size_t metab = itBalance->mIndex;
               mMetabFlags[metab].mFastReactions.insert(rct);
               mMetabFlags[metab].mFlag = FAST;
             }
         }
     }
 
   return;
 }

void CHybridMethodODE45::setupMethod ( )

protected

setup mMethod

Definition at line 466 of file CHybridMethodODE45.cpp.

References DETERMINISTIC, HYBRID, mHasDetermReaction, mHasStoiReaction, mMethod, and STOCHASTIC.

Referenced by initMethod().

 {
   if (mHasStoiReaction && !mHasDetermReaction)
     mMethod = STOCHASTIC;
   else if (!mHasStoiReaction && mHasDetermReaction)
     mMethod = DETERMINISTIC;
   else
     mMethod = HYBRID;
 
   return;
 }

void CHybridMethodODE45::setupPriorityQueue ( C_FLOAT64 startTime = 0.0 )

protected

Sets up the priority queue.

Parameters

startTime The time at which the simulation starts.

Sets up the priority queue, for pure SSA.

Parameters

startTime The time at which the simulation starts.

Definition at line 540 of file CHybridMethodODE45.cpp.

References C_FLOAT64, calculateAmu(), CIndexedPriorityQueue::clear(), generateReactionTime(), CIndexedPriorityQueue::initializeIndexPointer(), CIndexedPriorityQueue::insertStochReaction(), mNumReactions, and mPQ.

Referenced by initMethod().

 {
   size_t i;
   C_FLOAT64 time;
 
   mPQ.clear();
   mPQ.initializeIndexPointer(mNumReactions);
 
   for (i = 0; i < mNumReactions; i++)
     {
       calculateAmu(i);
       time = startTime + generateReactionTime(i);
       mPQ.insertStochReaction(i, time);
     }
 
   return;
 }

void CHybridMethodODE45::setupReactAffect ( )

protected

Definition at line 661 of file CHybridMethodODE45.cpp.

References mLocalBalances, mMetab2React, mNumReactions, and mReactAffect.

Referenced by initMethod().

 {
   mReactAffect.clear();
   mReactAffect.resize(mNumReactions);
 
   std::vector<std::set<size_t> >::iterator rctIt =
     mReactAffect.begin();
   const std::vector<std::set<size_t> >::iterator rctEndIt =
     mReactAffect.end();
 
   for (size_t rct = 0; rctIt != rctEndIt; ++rct, ++rctIt)
     {
       rctIt->clear();
 
       std::vector<CHybridODE45Balance>::iterator balIt =
         mLocalBalances[rct].begin();
       const std::vector<CHybridODE45Balance>::iterator balEndIt =
         mLocalBalances[rct].end();
 
       for (; balIt != balEndIt; ++balIt)
         {
           size_t metabId = balIt->mIndex;
 
           std::set<size_t>::iterator it =
             mMetab2React[metabId].begin();
           const std::set<size_t>::iterator endIt =
             mMetab2React[metabId].end();
 
           for (; it != endIt; ++it)
             rctIt->insert(*it);
         }
     }
 
   return;
 }

void CHybridMethodODE45::setupReactionFlags ( )

protected

setup mReactionFlags

Definition at line 443 of file CHybridMethodODE45.cpp.

References FAST, mHasDetermReaction, mHasStoiReaction, mNumReactions, mpReactions, mReactionFlags, CVector< CType >::resize(), and SLOW.

Referenced by initMethod().

 {
   mHasStoiReaction = false;
   mHasDetermReaction = false;
   mReactionFlags.resize(mNumReactions);
 
   for (size_t rct = 0; rct < mNumReactions; rct++)
     {
       if ((*mpReactions)[rct]->isFast())
         {
           mReactionFlags[rct] = FAST;
           mHasDetermReaction = true;
         }
       else
         {
           mReactionFlags[rct] = SLOW;
           mHasStoiReaction = true;
         }
     }
 
   return;
 }

void CHybridMethodODE45::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 259 of file CHybridMethodODE45.cpp.

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

 {
   //set inital state
   mpProblem->getModel()->setState(*initialState);
 
   //set mpState
   mpState = new CState(mpProblem->getModel()->getState());
 
   //set the mpModel
   mpModel = mpProblem->getModel();
   assert(mpModel);
 
   //set mDoCorrection
   if (mpModel->getModelType() == CModel::deterministic)
     mDoCorrection = true;
   else
     mDoCorrection = false;
 
   //set mHasAssignments
   mHasAssignments = modelHasAssignments(mpModel);
 
   //set mFirstMetabIndex
   mFirstMetabIndex = mpModel->getStateTemplate().getIndex(mpModel->getMetabolitesX()[0]);
 
   // Call initMethod function
   initMethod(mpProblem->getModel()->getTime());
 
   //output data to check init status
   //outputData();
   return;
 }

CTrajectoryMethod::Status CHybridMethodODE45::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 698 of file CHybridMethodODE45.cpp.

References C_FLOAT64, doSingleStep(), CModel::getTime(), mMaxSteps, mMaxStepsReached, CTrajectoryMethod::mpCurrentState, mpModel, mpState, CTrajectoryMethod::NORMAL, and CCopasiMessage::WARNING.

 {
   // write the current state to the model???
 
   // 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 = mpModel->getTime();
   C_FLOAT64 endTime = time + deltaT;
 
   for (i = 0; ((i < mMaxSteps) && (time < endTime)); i++)
     time = doSingleStep(time, endTime);
 
   // Warning Message
   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 = *mpState;
   //mpCurrentState->setTime(time);
 
   return NORMAL;
 }

void CHybridMethodODE45::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 1353 of file CHybridMethodODE45.cpp.

References C_FLOAT64, C_INVALID_INDEX, calculateAmu(), generateReactionTime(), mAmu, mAmuOld, mHasAssignments, mNumReactions, mpModel, mPQ, mReactAffect, CIndexedPriorityQueue::updateNode(), CModel::updateSimulatedValues(), and updateTauMu().

Referenced by 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 < mNumReactions; index++)
         {
           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
       iter = mReactAffect[rIndex].begin();
       iterEnd = mReactAffect[rIndex].end();
 
       for (; iter != iterEnd; iter++)
         {
           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)
     {
       // reaction is stochastic
       newTime = time + generateReactionTime(rIndex);
       mPQ.updateNode(rIndex, newTime);
     }
 
   return;
 }

void CHybridMethodODE45::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 1122 of file CHybridMethodODE45.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

CVector<C_FLOAT64> CHybridMethodODE45::mAmu

protected

The propensities of the stochastic reactions.

Definition at line 653 of file CHybridMethodODE45.h.

Referenced by calculateAmu(), evalF(), generateReactionTime(), getReactionIndex4Hybrid(), initMethod(), outputDebug(), updatePriorityQueue(), and updateTauMu().

CVector<C_FLOAT64> CHybridMethodODE45::mAmuOld

protected

Definition at line 654 of file CHybridMethodODE45.h.

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

C_FLOAT64 CHybridMethodODE45::mAtol

protected

Absolute tolerance.

Definition at line 638 of file CHybridMethodODE45.h.

Referenced by initMethod(), and integrateDeterministicPart().

std::set<size_t> CHybridMethodODE45::mCalculateSet

protected

Set of the reactions, which must be updated. Reactions have at lease one fast metab as subtract with non-zero balance

Definition at line 661 of file CHybridMethodODE45.h.

Referenced by evalF(), fireSlowReaction4Hybrid(), and setupCalculateSet().

Data CHybridMethodODE45::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 588 of file CHybridMethodODE45.h.

Referenced by CHybridMethodODE45(), doInverseInterpolation(), evalF(), initMethod(), and integrateDeterministicPart().

bool CHybridMethodODE45::mDefaultAtol

protected

Definition at line 633 of file CHybridMethodODE45.h.

Referenced by initMethod().

CDependencyGraph CHybridMethodODE45::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 713 of file CHybridMethodODE45.h.

Referenced by outputDebug(), and setupDependencyGraph().

bool CHybridMethodODE45::mDoCorrection

protected

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

Definition at line 472 of file CHybridMethodODE45.h.

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

CVector< C_FLOAT64 > CHybridMethodODE45::mDWork

protected

Definition at line 640 of file CHybridMethodODE45.h.

Referenced by initMethod(), and integrateDeterministicPart().

C_FLOAT64 CHybridMethodODE45::mEndTime

protected

Requested end time.

Definition at line 608 of file CHybridMethodODE45.h.

std::ostringstream CHybridMethodODE45::mErrorMsg

protected

Definition at line 751 of file CHybridMethodODE45.h.

Referenced by integrateDeterministicPart().

size_t CHybridMethodODE45::mFirstMetabIndex

protected

index of the first metab in CState

Definition at line 493 of file CHybridMethodODE45.h.

Referenced by outputData(), and start().

bool CHybridMethodODE45::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 746 of file CHybridMethodODE45.h.

Referenced by start(), and updatePriorityQueue().

bool CHybridMethodODE45::mHasDetermReaction

protected

Definition at line 549 of file CHybridMethodODE45.h.

Referenced by outputData(), setupMethod(), and setupReactionFlags().

bool CHybridMethodODE45::mHasStoiReaction

protected

Bool value

Definition at line 544 of file CHybridMethodODE45.h.

Referenced by outputData(), setupMethod(), and setupReactionFlags().

C_INT CHybridMethodODE45::mIFlag

protected

ODE45 state, corresponding to iflag in rkf45, an ODE45 solver argument

Definition at line 614 of file CHybridMethodODE45.h.

Referenced by integrateDeterministicPart().

CVector< C_INT > CHybridMethodODE45::mIWork

protected

Definition at line 641 of file CHybridMethodODE45.h.

Referenced by initMethod(), and integrateDeterministicPart().

std::vector<std::vector <CHybridODE45Balance> > CHybridMethodODE45::mLocalBalances

protected

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

Definition at line 522 of file CHybridMethodODE45.h.

Referenced by evalF(), fireReaction(), getAffects(), outputData(), outputDebug(), setupBalances(), setupMetabFlags(), and setupReactAffect().

std::vector<std::vector <CHybridODE45Balance> > CHybridMethodODE45::mLocalSubstrates

protected

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

Definition at line 528 of file CHybridMethodODE45.h.

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

size_t CHybridMethodODE45::mMaxBalance

protected

Is this useful in my method? maximal increase of a particle number in one step.

Definition at line 581 of file CHybridMethodODE45.h.

Referenced by outputDebug(), and setupBalances().

size_t CHybridMethodODE45::mMaxSteps

protected

Max number of doSingleStep() per step()

Definition at line 575 of file CHybridMethodODE45.h.

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

bool CHybridMethodODE45::mMaxStepsReached

protected

Definition at line 576 of file CHybridMethodODE45.h.

Referenced by initMethod(), and step().

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

protected

Vector of relations between metabolites with reactions.

Definition at line 533 of file CHybridMethodODE45.h.

Referenced by outputData(), outputDebug(), setupCalculateSet(), setupMetab2React(), and setupReactAffect().

std::vector<CHybridODE45MetabFlag> CHybridMethodODE45::mMetabFlags

protected

Vector holding information on the status of metabolites. They can be FAST or SLOW.

Definition at line 499 of file CHybridMethodODE45.h.

Referenced by outputData(), outputDebug(), setupCalculateSet(), and setupMetabFlags().

size_t CHybridMethodODE45::mMethod

protected

An Integer showes the method now CHybridMethod used 0 == Stochastic 1 == Deterministic 2 == Hybrid

Definition at line 557 of file CHybridMethodODE45.h.

Referenced by doSingleStep(), and setupMethod().

size_t CHybridMethodODE45::mNumReactions

protected

Number of Reactions

Definition at line 505 of file CHybridMethodODE45.h.

Referenced by evalF(), getReactionIndex4Hybrid(), initMethod(), outputData(), setupBalances(), setupDependencyGraph(), setupMetab2React(), setupMetabFlags(), setupPriorityQueue(), setupReactAffect(), setupReactionFlags(), and updatePriorityQueue().

size_t CHybridMethodODE45::mNumVariableMetabs

protected

Dimension of the system. Total number of metabolites.

Definition at line 483 of file CHybridMethodODE45.h.

Referenced by initMethod(), outputData(), outputDebug(), setupMetab2React(), and setupMetabFlags().

C_INT CHybridMethodODE45::mODE45Status

protected

state of ODE45, indicating what to do next in the step part -1, error emerge 0, finish all integration 1, finish one step integration 2, has event

Definition at line 623 of file CHybridMethodODE45.h.

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

C_FLOAT64 CHybridMethodODE45::mOldTime

protected

Definition at line 674 of file CHybridMethodODE45.h.

Referenced by doInverseInterpolation(), and integrateDeterministicPart().

C_FLOAT64* CHybridMethodODE45::mOldY

protected

Record of y and time of the previous step

Definition at line 673 of file CHybridMethodODE45.h.

Referenced by cleanup(), doInverseInterpolation(), initMethod(), and integrateDeterministicPart().

size_t CHybridMethodODE45::mOutputCounter

protected

Output counter.

Definition at line 738 of file CHybridMethodODE45.h.

std::ofstream CHybridMethodODE45::mOutputFile

protected

File output stream to write data.

Definition at line 728 of file CHybridMethodODE45.h.

std::string CHybridMethodODE45::mOutputFileName

protected

Output filename.

Definition at line 733 of file CHybridMethodODE45.h.

CStateRecord* CHybridMethodODE45::mpEventState

protected

Pointer to the State at Event Happens

Definition at line 684 of file CHybridMethodODE45.h.

Referenced by doInverseInterpolation().

CInterpolation* CHybridMethodODE45::mpInterpolation

protected

A Pointer to Object of CInterpolation

Definition at line 679 of file CHybridMethodODE45.h.

Referenced by cleanup(), doInverseInterpolation(), and initMethod().

CCopasiVector<CMetab>* CHybridMethodODE45::mpMetabolites

protected

A pointer to the metabolites of the model.

Definition at line 488 of file CHybridMethodODE45.h.

Referenced by initMethod(), outputData(), and outputDebug().

CModel* CHybridMethodODE45::mpModel

protected

Pointer to the model

Definition at line 462 of file CHybridMethodODE45.h.

Referenced by cleanup(), doInverseInterpolation(), doSingleStep(), evalF(), fireSlowReaction4Hybrid(), initMethod(), integrateDeterministicPart(), setupBalances(), setupMetab2React(), start(), step(), and updatePriorityQueue().

CIndexedPriorityQueue CHybridMethodODE45::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 722 of file CHybridMethodODE45.h.

Referenced by getStochTimeAndIndex(), outputDebug(), setupPriorityQueue(), updatePriorityQueue(), and updateTauMu().

CRandom* CHybridMethodODE45::mpRandomGenerator

protected

The random number generator.

Definition at line 695 of file CHybridMethodODE45.h.

Referenced by CHybridMethodODE45(), cleanup(), generateReactionTime(), getReactionIndex4Hybrid(), initMethod(), integrateDeterministicPart(), and outputDebug().

const CCopasiVectorNS<CReaction>* CHybridMethodODE45::mpReactions

protected

A pointer to the reactions of the model.

Definition at line 516 of file CHybridMethodODE45.h.

Referenced by checkModel(), getDependsOn(), initMethod(), outputDebug(), setupDependencyGraph(), and setupReactionFlags().

CState* CHybridMethodODE45::mpState

protected

A pointer to the current state in complete model view.

Definition at line 563 of file CHybridMethodODE45.h.

Referenced by cleanup(), doInverseInterpolation(), doSingleStep(), evalF(), fireSlowReaction4Hybrid(), integrateDeterministicPart(), start(), and step().

unsigned C_INT32 CHybridMethodODE45::mRandomSeed

protected

The random seed to use.

Definition at line 706 of file CHybridMethodODE45.h.

Referenced by initMethod().

std::vector<std::set<size_t> > CHybridMethodODE45::mReactAffect

protected

Vector of sets that the indeces of propencities which should be updated after one reaction has been fired

Definition at line 539 of file CHybridMethodODE45.h.

Referenced by fireSlowReaction4Hybrid(), outputData(), setupReactAffect(), and updatePriorityQueue().

CVector<size_t> CHybridMethodODE45::mReactionFlags

protected

Fast reactions are set FAST and slow ones are SLOW

Definition at line 511 of file CHybridMethodODE45.h.

Referenced by evalF(), getReactionIndex4Hybrid(), outputData(), outputDebug(), setupMetabFlags(), and setupReactionFlags().

bool CHybridMethodODE45::mReducedModel

protected

Definition at line 477 of file CHybridMethodODE45.h.

Referenced by outputData().

C_FLOAT64 CHybridMethodODE45::mRtol

protected

Relative tolerance.

Definition at line 628 of file CHybridMethodODE45.h.

Referenced by initMethod(), and integrateDeterministicPart().

CVector<C_FLOAT64> CHybridMethodODE45::mStateRecord

protected

A vector to record middle state for interpolation

Definition at line 689 of file CHybridMethodODE45.h.

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

CMatrix<C_FLOAT64> CHybridMethodODE45::mStoi

protected

The stoichometry matrix of the model.

Definition at line 467 of file CHybridMethodODE45.h.

Referenced by checkModel(), initMethod(), and outputDebug().

C_FLOAT64 CHybridMethodODE45::mTime

protected

Current time.

Definition at line 603 of file CHybridMethodODE45.h.

Referenced by doInverseInterpolation(), and integrateDeterministicPart().

std::set<size_t> CHybridMethodODE45::mUpdateSet

protected

Set of the reactions, which must update after one slow reaction fires

Definition at line 667 of file CHybridMethodODE45.h.

Referenced by initMethod(), and outputDebug().

bool CHybridMethodODE45::mUseRandomSeed

protected

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

Definition at line 701 of file CHybridMethodODE45.h.

Referenced by initMethod().

bool CHybridMethodODE45::mUseStateRecord

protected

A boolean variable to show whether mStateRecord is used or not. If t is too close to tout, rkfs45() won't go through into fehl() and uses Forward Euler method

Definition at line 648 of file CHybridMethodODE45.h.

Referenced by doInverseInterpolation(), and rkfs_().

C_FLOAT64* CHybridMethodODE45::mY

protected

Pointer to the array with left hand side values.

Definition at line 593 of file CHybridMethodODE45.h.

Referenced by cleanup(), doInverseInterpolation(), initMethod(), and integrateDeterministicPart().

CVector< C_FLOAT64 > CHybridMethodODE45::mYdot

protected

Vector containig the derivatives after calling eval

Definition at line 598 of file CHybridMethodODE45.h.

Referenced by initMethod().

CVector<C_FLOAT64> CHybridMethodODE45::temp

protected

Vectors to hold the system state and intermediate results

Definition at line 568 of file CHybridMethodODE45.h.

Referenced by initMethod(), and outputDebug().

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

Classes

Public Member Functions

Protected Member Functions

Static Protected Member Functions

Protected Attributes

Private Member Functions

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation