COPASI API  4.16.103
Classes | Public Member Functions | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Private Member Functions | Friends | List of all members
CHybridMethodODE45 Class Reference

#include <CHybridMethodODE45.h>

Inheritance diagram for CHybridMethodODE45:
Inheritance graph
[legend]
Collaboration diagram for CHybridMethodODE45:
Collaboration graph
[legend]

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::SubTypegetSubType () const
 
const CCopasiTask::TypegetType () 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)
 
CCopasiParameterGroupassertGroup (const std::string &name)
 
template<class CType >
CCopasiParameterassertParameter (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
 
CCopasiParameterGroupgetGroup (const std::string &name)
 
const CCopasiParameterGroupgetGroup (const std::string &name) const
 
CCopasiParameterGroupgetGroup (const size_t &index)
 
const CCopasiParameterGroupgetGroup (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 CObjectInterfacegetObject (const CCopasiObjectName &cn) const
 
CCopasiParametergetParameter (const std::string &name)
 
const CCopasiParametergetParameter (const std::string &name) const
 
CCopasiParametergetParameter (const size_t &index)
 
const CCopasiParametergetParameter (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::ValuegetValue (const std::string &name) const
 
const CCopasiParameter::ValuegetValue (const size_t &index) const
 
CCopasiParameter::ValuegetValue (const std::string &name)
 
CCopasiParameter::ValuegetValue (const size_t &index)
 
CCopasiParameterGroupoperator= (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::TypegetType () const
 
const ValuegetValue () const
 
ValuegetValue ()
 
virtual voidgetValuePointer () const
 
CCopasiObjectgetValueReference () 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
 
CCopasiParameteroperator= (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 objectMapgetObjects () const
 
virtual std::string getUnits () const
 
virtual const CCopasiObjectgetValueObject () 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 DataObjectSetgetDirectDependencies (const DataObjectSet &context=DataObjectSet()) const
 
CCopasiContainergetObjectAncestor (const std::string &type) const
 
CCopasiDataModelgetObjectDataModel ()
 
const CCopasiDataModelgetObjectDataModel () const
 
const std::string & getObjectName () const
 
CCopasiContainergetObjectParent () const
 
const std::string & getObjectType () const
 
virtual const
CObjectInterface::ObjectSet
getPrerequisites () const
 
virtual RefreshgetRefresh () const
 
UpdateMethodgetUpdateMethod () 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 >
CCopasiObjectaddMatrixReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
 
template<class CType >
CCopasiObjectaddObjectReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
 
template<class CType >
CCopasiObjectaddVectorReference (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_FLOAT64mAmu
 
CVector< C_FLOAT64mAmuOld
 
C_FLOAT64 mAtol
 
std::set< size_t > mCalculateSet
 
Data mData
 
bool mDefaultAtol
 
CDependencyGraph mDG
 
bool mDoCorrection
 
CVector< C_FLOAT64mDWork
 
C_FLOAT64 mEndTime
 
std::ostringstream mErrorMsg
 
size_t mFirstMetabIndex
 
bool mHasAssignments
 
bool mHasDetermReaction
 
bool mHasStoiReaction
 
C_INT mIFlag
 
CVector< C_INTmIWork
 
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_FLOAT64mOldY
 
size_t mOutputCounter
 
std::ofstream mOutputFile
 
std::string mOutputFileName
 
CStateRecordmpEventState
 
CInterpolationmpInterpolation
 
CCopasiVector< CMetab > * mpMetabolites
 
CModelmpModel
 
CIndexedPriorityQueue mPQ
 
CRandommpRandomGenerator
 
const CCopasiVectorNS
< CReaction > * 
mpReactions
 
CStatempState
 
unsigned C_INT32 mRandomSeed
 
std::vector< std::set< size_t > > mReactAffect
 
CVector< size_t > mReactionFlags
 
bool mReducedModel
 
C_FLOAT64 mRtol
 
CVector< C_FLOAT64mStateRecord
 
CMatrix< C_FLOAT64mStoi
 
C_FLOAT64 mTime
 
std::set< size_t > mUpdateSet
 
bool mUseRandomSeed
 
bool mUseStateRecord
 
C_FLOAT64mY
 
CVector< C_FLOAT64mYdot
 
CVector< C_FLOAT64temp
 
- Protected Attributes inherited from CTrajectoryMethod
CStatempCurrentState
 
CTrajectoryProblemmpProblem
 
CVector< C_INTmRoots
 
- Protected Attributes inherited from CCopasiMethod
CProcessReportmpCallBack
 
- Protected Attributes inherited from CCopasiParameter
std::string mKey
 
CCopasiObjectmpValueReference
 
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

CTrajectoryMethodCTrajectoryMethod::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 CTrajectoryMethodcreateMethod (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 CRenameHandlersmpRenameHandler = 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.

61  :
63 {
64  assert((void *) &mData == (void *) &mData.dim);
65  mData.pMethod = this;
66 
69 }
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
CHybridMethodODE45 * pMethod
CHybridMethodODE45::CHybridMethodODE45 ( const CHybridMethodODE45 src,
const CCopasiContainer pParent = NULL 
)

Copy constructor

Parameters
constCHybridMethodODE45 & src
constCCopasiContainer * 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.

75  :
76  CTrajectoryMethod(src, pParent)
77 {
78  assert((void *) &mData == (void *) &mData.dim);
79  mData.pMethod = this;
80 
83 }
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
CHybridMethodODE45 * pMethod
CHybridMethodODE45::~CHybridMethodODE45 ( )

Destructor.

Definition at line 88 of file CHybridMethodODE45.cpp.

References cleanup(), and DESTRUCTOR_TRACE.

89 {
90  cleanup();
92 }
#define DESTRUCTOR_TRACE
Definition: copasi.h:206

Member Function Documentation

void CHybridMethodODE45::calculateAmu ( size_t  rIndex)
protected

Calculates an amu value for a given reaction.

Parameters
rIndexA 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().

1155 {
1156  if (!mDoCorrection)
1157  {
1158  mAmu[rIndex] = (*mpReactions)[rIndex]->calculateParticleFlux();
1159  return;
1160  }
1161 
1162  // We need the product of the cmu and hmu for this step.
1163  // We calculate this in one go, as there are fewer steps to
1164  // perform and we eliminate some possible rounding errors.
1165  C_FLOAT64 amu = 1; // initially
1166  //size_t total_substrates = 0;
1167  C_INT32 num_ident = 0;
1168  C_INT32 number = 0;
1169  C_INT32 lower_bound;
1170  // substrate_factor - The substrates, raised to their multiplicities,
1171  // multiplied with one another. If there are, e.g. m substrates of type m,
1172  // and n of type N, then substrate_factor = M^m * N^n.
1173  C_FLOAT64 substrate_factor = 1;
1174  // First, find the reaction associated with this index.
1175  // Keep a pointer to this.
1176  // Iterate through each substrate in the reaction
1177  const std::vector<CHybridODE45Balance> & substrates = mLocalSubstrates[rIndex];
1178 
1179  int flag = 0;
1180 
1181  for (size_t i = 0; i < substrates.size(); i++)
1182  {
1183  num_ident = substrates[i].mMultiplicity;
1184 
1185  if (num_ident > 1)
1186  {
1187  flag = 1;
1188  number = static_cast<C_INT32>((*mpMetabolites)[substrates[i].mIndex]->getValue());
1189  lower_bound = number - num_ident;
1190  substrate_factor = substrate_factor
1191  * pow((double) number, (int) num_ident - 1); //optimization
1192 
1193  number--; // optimization
1194 
1195  while (number > lower_bound)
1196  {
1197  amu *= number;
1198  number--;
1199  }
1200  }
1201  }
1202 
1203  if ((amu == 0) || (substrate_factor == 0)) // at least one substrate particle number is zero
1204  {
1205  mAmu[rIndex] = 0;
1206  return;
1207  }
1208 
1209  // rate_factor is the rate function divided by substrate_factor.
1210  // It would be more efficient if this was generated directly, since in effect we
1211  // are multiplying and then dividing by the same thing (substrate_factor)!
1212  C_FLOAT64 rate_factor = (*mpReactions)[rIndex]->calculateParticleFlux();
1213 
1214  if (flag)
1215  {
1216  amu *= rate_factor / substrate_factor;;
1217  mAmu[rIndex] = amu;
1218  }
1219  else
1220  {
1221  mAmu[rIndex] = rate_factor;
1222  }
1223 
1224  return;
1225 
1226  // a more efficient way to calculate mass action kinetics could be included
1227 }
#define C_INT32
Definition: copasi.h:90
std::vector< std::vector< CHybridODE45Balance > > mLocalSubstrates
CVector< C_FLOAT64 > mAmu
long int flag
Definition: f2c.h:52
#define C_FLOAT64
Definition: copasi.h:92
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().

2441 {
2443  CMatrix <C_FLOAT64> mStoi = model->getStoi();
2444  C_INT32 multInt;
2445  size_t i, j, numReactions = mpReactions->size();
2446  C_FLOAT64 multFloat;
2447  // size_t metabSize = mpMetabolites->size();
2448 
2449  for (i = 0; i < numReactions; i++) // for every reaction
2450  {
2451  // TEST getCompartmentNumber() == 1
2452  if ((*mpReactions)[i]->getCompartmentNumber() != 1) return - 1;
2453 
2454  // TEST isReversible() == 0
2455  if ((*mpReactions)[i]->isReversible() != 0) return - 2;
2456 
2457  // TEST integer stoichometry
2458  // Iterate through each the metabolites
2459  // juergen: the number of rows of mStoi equals the number of non-fixed metabs!
2460  // for (j=0; i<metabSize; j++)
2461  for (j = 0; j < mStoi.numRows(); j++)
2462  {
2463  multFloat = mStoi[j][i];
2464  multInt = static_cast<C_INT32>(floor(multFloat + 0.5)); // +0.5 to get a rounding out of the static_cast to int!
2465 
2466  if ((multFloat - multInt) > INT_EPSILON) return - 3; // INT_EPSILON in CHybridMethod.h
2467  }
2468  }
2469 
2470  return 1; // Model is appropriate for hybrid simulation
2471 }
const CCopasiVectorNS< CReaction > * mpReactions
virtual size_t size() const
virtual size_t numRows() const
Definition: CMatrix.h:138
#define INT_EPSILON
Definition: CHybridMethod.h:59
CMatrix< C_FLOAT64 > mStoi
#define C_INT32
Definition: copasi.h:90
#define C_FLOAT64
Definition: copasi.h:92
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
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().

381 {
382  delete mpRandomGenerator;
383  mpRandomGenerator = NULL;
384  mpModel = NULL;
385 
386  delete mpState;
387  mpState = NULL;
388 
389  delete mpInterpolation;
390  mpInterpolation = NULL;
391 
392  delete [] mY;
393  mY = NULL;
394 
395  delete [] mOldY;
396  mOldY = NULL;
397 
398  return;
399 }
CInterpolation * mpInterpolation
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().

978 {
979 
980  //==(1)==for one-step method, reset record each time when do interpolation
982 
983  //==(2)==set record in class interpolation
985  size_t offset;
986 
987  if (mUseStateRecord)
988  {
989  for (size_t i = 0; i < INTERP_RECORD_NUM - 2; i++) //record the middle 4 states
990  {
991  offset = i * (mData.dim + 1);
993  mStateRecord.array() + offset + 1);
994  }
995  }
996 
998 
999  //==(3)==do interpolation
1001 
1002  //==(4)==record the state
1003  mTime = mpEventState->getTime();
1004  C_FLOAT64 * tmpY = mpEventState->getArray() + 1;
1005  C_FLOAT64 * stateY = mpState->beginIndependent();
1006 
1007  for (size_t i = 0; i < mData.dim - 1; i++)
1008  stateY[i] = tmpY[i]; //write result into mpState
1009 
1010  mpState->setTime(mTime);
1012  mpModel->updateSimulatedValues(false); //for assignments?????????
1013 
1014  return;
1015 }
void recordState(const C_FLOAT64 time, const C_FLOAT64 *y)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
#define INTERP_RECORD_NUM
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
void setState(const CState &state)
Definition: CModel.cpp:1785
C_FLOAT64 * getArray() const
const C_FLOAT64 & getTime() const
CInterpolation * mpInterpolation
CStateRecord * getInterpolationState()
CVector< C_FLOAT64 > mStateRecord
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
CStateRecord * mpEventState
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
currentTimeA C_FLOAT64 specifying the current time
endTimeA 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
currentTimeA C_FLOAT64 specifying the current time
endTimeA 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().

754 {
755  C_FLOAT64 ds = 0.0;
756  size_t rIndex = 0;
757 
758  //1----pure SSA method
759  if (mMethod == STOCHASTIC) //has only stochastic reactions
760  {
761  getStochTimeAndIndex(ds, rIndex);
762 
763  if (ds > endTime)
764  ds = endTime;
765  else
766  {
767  fireReaction(rIndex);
768  updatePriorityQueue(rIndex, ds);
769  }
770 
771  //population has been changed and record to the mCurrentState
772  //in fireReaction(). So here just store time
773  *mpState = mpModel->getState();
774  mpState->setTime(ds);
776 
777  //outputState(mpCurrentState);
778  //getchar();
779  }
780  //2----Method with Deterministic Part
781  else if (mMethod == DETERMINISTIC) //has only deterministic reactions
782  {
783  integrateDeterministicPart(endTime - currentTime);
784  ds = mpState->getTime();
786  // Till now, state has been recorded
787  }
788  else if (mMethod == HYBRID)//Hybrid Method
789  {
790  integrateDeterministicPart(endTime - currentTime);
791  ds = mpState->getTime();
792 
793  if (mODE45Status == HAS_EVENT) //fire slow reaction
794  {
796 
799  }
800  else
802  }
803 
804  return ds;
805 }
#define CONTINUE
#define HYBRID
#define HAS_EVENT
void getStochTimeAndIndex(C_FLOAT64 &ds, size_t &rIndex)
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
void updatePriorityQueue(size_t rIndex, C_FLOAT64 time)
#define NEW_STEP
void fireReaction(size_t rIndex)
void setState(const CState &state)
Definition: CModel.cpp:1785
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
#define STOCHASTIC
const CState & getState() const
Definition: CModel.cpp:1771
void integrateDeterministicPart(C_FLOAT64 ds)
#define DETERMINISTIC
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().

97 {
99  return true;
100 }
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().

1027 {
1028  size_t i;
1029 
1030  //========No need to record State, right?========
1031  //maybe here I just want the propensities amu???
1032  //calculateDerivative(temp);
1033 
1034  //1====calculate propensities
1035  //just care about reactions involving fast metab
1036 
1037  //(1)Put current time *t and independent values *y into model.
1038  // This step seemes necessary, since propensity calculation process
1039  // requires functions called from the model class.
1040  mpState->setTime(*t);
1041  C_FLOAT64 * tmpY = mpState->beginIndependent();//mpState is a local copy
1042 
1043  for (i = 0; i < mData.dim - 1; ++i)
1044  tmpY[i] = y[i];
1045 
1047  mpModel->updateSimulatedValues(false); //really?
1048 
1049  //(2)Calculate propensities.
1050  std::set <size_t>::iterator reactIt = mCalculateSet.begin();
1051  size_t reactID;
1052 
1053  for (; reactIt != mCalculateSet.end(); reactIt++)
1054  {
1055  reactID = *reactIt;
1056  calculateAmu(reactID);
1057  }
1058 
1059  //2====calculate derivative
1060  //(1)initialize
1061  for (i = 0; i < mData.dim; i++)
1062  ydot[i] = 0.0;
1063 
1064  //(2)go through all the reactions and
1065  //update derivatives
1066  std::vector <CHybridODE45Balance>::iterator metabIt;
1067  std::vector <CHybridODE45Balance>::iterator metabEndIt;
1068  size_t metabIndex;
1069 
1070  for (i = 0; i < mNumReactions; i++)
1071  {
1072  if (mReactionFlags[i] == SLOW) //slow reaction
1073  ydot[mData.dim - 1] += mAmu[i];
1074  else //fast reaction
1075  {
1076  metabIt = mLocalBalances[i].begin();
1077  metabEndIt = mLocalBalances[i].end();
1078 
1079  for (; metabIt != metabEndIt; metabIt++)
1080  {
1081  metabIndex = metabIt->mIndex;
1082  ydot[metabIndex] += metabIt->mMultiplicity * mAmu[i];
1083  }
1084  }
1085  }
1086 
1087  return;
1088 }
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
void setState(const CState &state)
Definition: CModel.cpp:1785
CVector< C_FLOAT64 > mAmu
#define SLOW
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
CVector< size_t > mReactionFlags
void calculateAmu(size_t rIndex)
std::set< size_t > mCalculateSet
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_().

2272 {
2273  /* System generated locals */
2274  C_INT i__1;
2275  double d__1;
2276 
2277  /* Local variables */
2278  static C_INT k;
2279  static double ch;
2280  static C_INT base;
2281 
2282  /* fehlberg fourth-fifth order runge-kutta method */
2283 
2284  /* fehl integrates a system of neqn first order */
2285  /* ordinary differential equations of the form */
2286  /* dy(i)/dt=f(t,y(1),---,y(neqn)) */
2287  /* where the initial values y(i) and the initial derivatives */
2288  /* yp(i) are specified at the starting point t. fehl advances */
2289  /* the solution over the fixed step h and returns */
2290  /* the fifth order (sixth order accurate locally) solution */
2291  /* approximation at t+h in array s(i). */
2292  /* f1,---,f5 are arrays of dimension neqn which are needed */
2293  /* for internal storage. */
2294  /* the formulas have been grouped to control loss of significance. */
2295  /* fehl should be called with an h not smaller than 13 units of */
2296  /* roundoff in t so that the various independent arguments can be */
2297  /* distinguished. */
2298 
2299  /* Parameter adjustments */
2300  --yrcd;
2301  --s;
2302  --f5;
2303  --f4;
2304  --f3;
2305  --f2;
2306  --f1;
2307  --yp;
2308  --y;
2309 
2310  /* Function Body */
2311 
2312  //~~~~~~~~(1)~~~~~~~~
2313  ch = *h__ / 4.;
2314  i__1 = *neqn;
2315 
2316  for (k = 1; k <= i__1; ++k)
2317  {
2318  f5[k] = y[k] + ch * yp[k];
2319  }
2320 
2321  d__1 = *t + ch;
2322  (*f)(neqn, &d__1, &f5[1], &f1[1]);
2323 
2324  /* Record trcd and yrcd at t+h/4 */
2325  base = 1;
2326  yrcd[base] = *t + ch;
2327  i__1 = *neqn;
2328 
2329  for (k = 1; k <= i__1; ++k)
2330  {
2331  yrcd[base + k] = f5[k];
2332  }
2333 
2334  //~~~~~~~~(2)~~~~~~~~
2335  ch = *h__ * 3. / 32.;
2336  i__1 = *neqn;
2337 
2338  for (k = 1; k <= i__1; ++k)
2339  {
2340  f5[k] = y[k] + ch * (yp[k] + f1[k] * 3.);
2341  }
2342 
2343  d__1 = *t + *h__ * 3. / 8.;
2344  (*f)(neqn, &d__1, &f5[1], &f2[1]);
2345 
2346  /* Record trcd and yrcd at t+3h/8 */
2347  base = *neqn + 2;
2348  yrcd[base] = *t + *h__ * 3. / 8.;
2349  i__1 = *neqn;
2350 
2351  for (k = 1; k <= i__1; ++k)
2352  {
2353  yrcd[base + k] = f5[k];
2354  }
2355 
2356  //~~~~~~~~(3)~~~~~~~~
2357  ch = *h__ / 2197.;
2358  i__1 = *neqn;
2359 
2360  for (k = 1; k <= i__1; ++k)
2361  {
2362  f5[k] = y[k] + ch * (yp[k] * 1932. + (f2[k] * 7296.
2363  - f1[k] * 7200.));
2364  }
2365 
2366  d__1 = *t + *h__ * 12. / 13.;
2367  (*f)(neqn, &d__1, &f5[1], &f3[1]);
2368 
2369  /* Record trcd and yrcd at t+12h/13 */
2370  base = (*neqn + 1) * 3 + 1;
2371  yrcd[base] = *t + *h__ * 12. / 13.;
2372  i__1 = *neqn;
2373 
2374  for (k = 1; k <= i__1; ++k)
2375  {
2376  yrcd[base + k] = f5[k];
2377  }
2378 
2379  //~~~~~~~~(4)~~~~~~~~
2380  ch = *h__ / 4104.;
2381  i__1 = *neqn;
2382 
2383  for (k = 1; k <= i__1; ++k)
2384  {
2385  f5[k] = y[k] + ch * (yp[k] * 8341. - f3[k] * 845.
2386  + (f2[k] * 29440. - f1[k] * 32832.));
2387  }
2388 
2389  d__1 = *t + *h__;
2390  (*f)(neqn, &d__1, &f5[1], &f4[1]);
2391 
2392  //~~~~~~~~(5)~~~~~~~~
2393  ch = *h__ / 20520.;
2394  i__1 = *neqn;
2395 
2396  for (k = 1; k <= i__1; ++k)
2397  {
2398  f1[k] = y[k] + ch * (yp[k] * -6080.
2399  + (f3[k] * 9295. - f4[k] * 5643.)
2400  + (f1[k] * 41040. - f2[k] * 28352.));
2401  }
2402 
2403  d__1 = *t + *h__ / 2.;
2404  (*f)(neqn, &d__1, &f1[1], &f5[1]);
2405 
2406  /* Record trcd and yrcd at t+h/2 */
2407  //base = (*neqn + 1 << 1) + 1;
2408  base = (*neqn + 1) * 2 + 1;
2409  yrcd[base] = *t + *h__ / 2.;
2410  i__1 = *neqn;
2411 
2412  for (k = 1; k <= i__1; ++k)
2413  {
2414  yrcd[base + k] = f1[k];
2415  }
2416 
2417  //~~~~~~~~(6)~~~~~~~~
2418  /* compute approximate solution at t+h */
2419  ch = *h__ / 7618050.;
2420  i__1 = *neqn;
2421 
2422  for (k = 1; k <= i__1; ++k)
2423  {
2424  s[k] = y[k] + ch * (yp[k] * 902880.
2425  + (f3[k] * 3855735. - f4[k] * 1371249.)
2426  + (f2[k] * 3953664. + f5[k] * 277020.));
2427  }
2428 
2429  return 0;
2430 } /* fehl_ */
#define C_INT
Definition: copasi.h:115
void CHybridMethodODE45::fireReaction ( size_t  rIndex)
protected

Executes the specified reaction in the system once.

Parameters
rIndexA 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().

1411 {
1412  // Change the particle numbers according to which step took place.
1413  // First, get the vector of balances in the reaction we've got.
1414  // (This vector expresses the number change of each metabolite
1415  // in the reaction.) Then step through each balance, using its
1416  // multiplicity to calculate a new value for the associated
1417  // metabolite. Finally, update the metabolite.
1418 
1419  size_t i;
1420  C_FLOAT64 newNumber;
1421 
1422  CMetab * pMetab;
1423 
1424  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
1425  {
1426  pMetab = mLocalBalances[rIndex][i].mpMetabolite;
1427  newNumber = pMetab->getValue() + mLocalBalances[rIndex][i].mMultiplicity;
1428 
1429  pMetab->setValue(newNumber);
1430  pMetab->refreshConcentration();
1431  }
1432 
1433  return;
1434 }
virtual void setValue(const C_FLOAT64 &value)
Definition: CMetab.h:178
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
void refreshConcentration()
Definition: CMetab.cpp:281
const C_FLOAT64 & getValue() const
#define C_FLOAT64
Definition: copasi.h:92
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().

808 {
809  //First Update Propensity
810  std::set <size_t>::iterator reactIt = mCalculateSet.begin();
811  size_t reactID;
812 
813  for (; reactIt != mCalculateSet.end(); reactIt++)
814  {
815  reactID = *reactIt;
816  calculateAmu(reactID);
817  }
818 
819  reactID = getReactionIndex4Hybrid();
820  fireReaction(reactID); //Directly update current state in global view
821  *mpState = mpModel->getState();
822 
823  //update corresponding propensities
824  std::set <size_t>::iterator updateIt
825  = mReactAffect[reactID].begin();
826  const std::set <size_t>::iterator updateEndIt
827  = mReactAffect[reactID].end();
828 
829  for (; updateIt != updateEndIt; updateIt++)
830  {
831  reactID = *updateIt;
832  calculateAmu(reactID);
833  }
834 
835  return;
836 }
void fireReaction(size_t rIndex)
std::vector< std::set< size_t > > mReactAffect
const CState & getState() const
Definition: CModel.cpp:1771
void calculateAmu(size_t rIndex)
std::set< size_t > mCalculateSet
C_FLOAT64 CHybridMethodODE45::generateReactionTime ( size_t  rIndex)
protected

Generates a putative reaction time for the given reaction.

Parameters
rIndexA 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().

1108 {
1109  if (mAmu[rIndex] == 0) return std::numeric_limits<C_FLOAT64>::infinity();
1110 
1112  return - 1 * log(rand2) / mAmu[rIndex];
1113 }
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
CVector< C_FLOAT64 > mAmu
#define C_FLOAT64
Definition: copasi.h:92
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
rIndexThe 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().

1267 {
1268  size_t i;
1269  std::set<std::string> *retset = new std::set<std::string>;
1270 
1271  // Get the balances associated with the reaction at this index
1272  // XXX We first get the chemical equation, then the balances, since the getBalances method in CReaction is unimplemented!
1273 
1274  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
1275  {
1276  if (mLocalBalances[rIndex][i].mMultiplicity != 0)
1277  retset->insert(mLocalBalances[rIndex][i].mpMetabolite->getKey());
1278  }
1279 
1280  return retset;
1281 }
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
C_FLOAT64 CHybridMethodODE45::getDefaultAtol ( const CModel pModel) const
protected

Calculate the default absolute tolerance

Parameters
constCModel * 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().

840 {
841  if (!pModel)
842  return 1.0e009;
843 
844  const CCopasiVectorNS< CCompartment > & Compartment =
845  pModel->getCompartments();
846 
847  size_t i, imax;
848 
850 
851  for (i = 0, imax = Compartment.size(); i < imax; i++)
852  {
853  if (Compartment[i]->getValue() < Volume)
854  Volume = Compartment[i]->getValue();
855  }
856 
858  return 1.0e009;
859 
860  return Volume * pModel->getQuantity2NumberFactor() * 1.e-12;
861 }
virtual size_t size() const
const C_FLOAT64 & getQuantity2NumberFactor() const
Definition: CModel.cpp:2354
const Value & getValue() const
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
#define C_FLOAT64
Definition: copasi.h:92
#define max(a, b)
Definition: f2c.h:176
std::set< std::string > * CHybridMethodODE45::getDependsOn ( size_t  rIndex)
private

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

Parameters
rIndexThe 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().

1236 {
1237  std::set<std::string> *retset = new std::set<std::string>;
1238 
1239  size_t i, imax = (*mpReactions)[rIndex]->getFunctionParameters().size();
1240  size_t j, jmax;
1241 
1242  for (i = 0; i < imax; ++i)
1243  {
1244  if ((*mpReactions)[rIndex]->getFunctionParameters()[i]->getUsage()
1246  continue;
1247 
1248  const std::vector <std::string> & metabKeylist =
1249  (*mpReactions)[rIndex]->getParameterMappings()[i];
1250  jmax = metabKeylist.size();
1251 
1252  for (j = 0; j < jmax; ++j)
1253  retset->insert(metabKeylist[j]);
1254  }
1255 
1256  return retset;
1257 }
const CCopasiVectorNS< CReaction > * mpReactions
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().

1442 {
1443  //calculate sum of amu
1444  C_FLOAT64 mAmuSum = 0.0;
1445 
1446  for (int i = 0; i < mNumReactions; i++)
1447  {
1448  if (mReactionFlags[i] == SLOW)
1449  mAmuSum += mAmu[i];
1450  }
1451 
1452  //get the threshold
1454  C_FLOAT64 threshold = mAmuSum * rand2;
1455 
1456  //get the reaction index
1457  size_t rIndex;
1458  C_FLOAT64 tmp = 0.0;
1459 
1460  //is there some algorithm that can get a log() complex?
1461  for (rIndex = 0; rIndex < mNumReactions; rIndex++)
1462  {
1463  if (mReactionFlags[rIndex] == SLOW)
1464  {
1465  tmp += mAmu[rIndex];
1466 
1467  if (tmp >= threshold)
1468  return rIndex;
1469  }
1470  }
1471  return rIndex;
1472 }
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
CVector< C_FLOAT64 > mAmu
#define SLOW
#define C_FLOAT64
Definition: copasi.h:92
CVector< size_t > mReactionFlags
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
dsA reference to a C_FLOAT64. The putative reaction time for the first stochastic reaction is written into this variable.
rIndexA 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().

1101 {
1102  ds = mPQ.topKey();
1103  rIndex = mPQ.topIndex();
1104  return;
1105 }
CIndexedPriorityQueue mPQ
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().

176 {
177  CCopasiParameter *pParm;
178 
179  //-----------------------------------------------------------------------
180  // This part will be dealt with later, when considering the partition strategy.
181  //Now, under the condition of using statistic partition given by users,
182  //we won't adopt such parameters listed here.
183  assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) MAX_STEPS);
184  //assertParameter("Lower Limit", CCopasiParameter::DOUBLE, (C_FAT64) LOWER_STOCH_LIMIT);
185  //assertParameter("Upper Limit", CCopasiParameter::DOUBLE, (C_FLOAT64) UPPER_STOCH_LIMIT);
186 
187  //assertParameter("Partitioning Interval", CCopasiParameter::UINT, (unsigned C_INT32) PARTITIONING_INTERVAL);
188  //-----------------------------------------------------------------------
189 
190  assertParameter("Use Random Seed", CCopasiParameter::BOOL, (bool) USE_RANDOM_SEED);
191  assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) RANDOM_SEED);
192 
193  // Check whether we have a method with the old parameter names
194  if ((pParm = getParameter("HYBRID.MaxSteps")) != NULL)
195  {
196  //-------------------------------------------------------------------
197  setValue("Max Internal Steps", *pParm->getValue().pUINT);
198  removeParameter("HYBRID.MaxSteps");
199 
200  /*
201  if ((pParm = getParameter("HYBRID.LowerStochLimit")) != NULL)
202  {
203  setValue("Lower Limit", *pParm->getValue().pDOUBLE);
204  removeParameter("HYBRID.LowerStochLimit");
205  }
206 
207  if ((pParm = getParameter("HYBRID.UpperStochLimit")) != NULL)
208  {
209  setValue("Upper Limit", *pParm->getValue().pDOUBLE);
210  removeParameter("HYBRID.UpperStochLimit");
211  }
212 
213  if ((pParm = getParameter("HYBRID.PartitioningInterval")) != NULL)//need this part?
214  {
215  setValue("Partitioning Interval", *pParm->getValue().pUINT);
216  removeParameter("HYBRID.PartitioningInterval");
217  }
218  */
219  //-------------------------------------------------------------------
220 
221  if ((pParm = getParameter("UseRandomSeed")) != NULL)
222  {
223  setValue("Use Random Seed", *pParm->getValue().pBOOL);
224  removeParameter("UseRandomSeed");
225  }
226 
227  if ((pParm = getParameter("")) != NULL)
228  {
229  setValue("Random Seed", *pParm->getValue().pUINT);
230  removeParameter("");
231  }
232  }
233 
234  /* ODE45 ****************************************************************************/
235  //----------------------------------------------------------------------
236  //addParameter("Partitioning Stepsize",
237  // CCopasiParameter::DOUBLE, (C_FLOAT64) PARTITIONING_STEPSIZE);
238 
239  //addParameter("Max Internal Steps (LSOA)",
240  // CCopasiParameter::UINT, (unsigned C_INT32) 10000);
241  //-----------------------------------------------------------------------
242 
243  addParameter("Relative Tolerance",
245 
246  addParameter("Use Default Absolute Tolerance",
247  CCopasiParameter::BOOL, (bool) true);
248 
249  addParameter("Absolute Tolerance",
251 
252  //????????????????????????????????????????????????????????
253  // These parameters are no longer supported.
254  //removeParameter("Adams Max Order");
255  //removeParameter("BDF Max Order");
256  //????????????????????????????????????????????????????????
257 }
#define USE_RANDOM_SEED
Definition: CHybridMethod.h:67
#define RANDOM_SEED
Definition: CHybridMethod.h:68
#define C_INT32
Definition: copasi.h:90
#define MAX_STEPS
Definition: CHybridMethod.h:58
bool removeParameter(const std::string &name)
const Value & getValue() const
bool setValue(const std::string &name, const CType &value)
unsigned C_INT32 * pUINT
CCopasiParameter * getParameter(const std::string &name)
#define C_FLOAT64
Definition: copasi.h:92
bool addParameter(const CCopasiParameter &parameter)
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
void CHybridMethodODE45::initMethod ( C_FLOAT64  start_time)
protected

Initializes the solver.

Parameters
timethe current time

Initializes the solver and sets the model to be used.

Parameters
modelA 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().

296 {
297  //(1)----set attributes related with REACTIONS
300 
303 
305 
306  //(2)----set attributes related with METABS
308 
311 
312  setupBalances(); //initialize mLocalBalances and mLocalSubstrates
313  setupMetab2React(); //initialize mMetab2React
314  setupMetabFlags();
315 
316  //(3)----set attributes related with STATE
318 
319  //(4)----set attributes related with SYSTEM
320  mMaxSteps = * getValue("Max Internal Steps").pUINT;
321  mMaxStepsReached = false;
322  setupMethod();
323 
324  //(5)----set attributes related with STOCHASTIC
325  mUseRandomSeed = * getValue("Use Random Seed").pBOOL;
326  mRandomSeed = * getValue("Random Seed").pUINT;
327 
328  if (mUseRandomSeed)
330 
331  mStoi = mpModel->getStoi();
332  mUpdateSet.clear();// how to set this parameter
333  setupCalculateSet(); //should be done after setupBalances()
335 
336  setupDependencyGraph(); // initialize mDG
337  setupPriorityQueue(start_time); // initialize mPQ
338 
339  //(6)----set attributes for INTERPOLATION
342 
343  //(7)----set attributes for ODE45
345  mData.dim = (C_INT)(mNumVariableMetabs + 1); //one more for sum of propensities
347 
348  mRtol = * getValue("Relative Tolerance").pUDOUBLE;
349  mDefaultAtol = * getValue("Use Default Absolute Tolerance").pBOOL;
350 
351  if (mDefaultAtol)
352  {
354  setValue("Absolute Tolerance", mAtol);
355  }
356  else
357  mAtol = * getValue("Absolute Tolerance").pUDOUBLE;
358 
359  mDWork.resize(3 + 6 * mData.dim);
360  mIWork.resize(5);
361 
362  //setupUpdateSet();
363 
364  // we don't want to directly record new results into mpState, since
365  // the sum of slow reaction propensities is also recorded in mY
366  mY = new C_FLOAT64[mData.dim];
367  mOldY = new C_FLOAT64[mData.dim];
368 
369  //first calculate propensities, in the next integration process
370  //system will just update the propensities record in mCalculateSet
371  for (size_t i = 0; i < mNumReactions; i++)
372  calculateAmu(i);
373 
374  return;
375 }
std::set< size_t > mUpdateSet
#define C_INT
Definition: copasi.h:115
CVector< C_FLOAT64 > temp
const CCopasiVectorNS< CReaction > * mpReactions
virtual size_t size() const
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
CVector< C_FLOAT64 > mAmuOld
CTrajectoryProblem * mpProblem
CMatrix< C_FLOAT64 > mStoi
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
unsigned C_INT32 mRandomSeed
#define INTERP_RECORD_NUM
C_FLOAT64 getDefaultAtol(const CModel *pModel) const
CVector< C_FLOAT64 > mYdot
#define NEW_STEP
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
CVector< C_FLOAT64 > mAmu
const Value & getValue() const
bool setValue(const std::string &name, const CType &value)
unsigned C_INT32 * pUINT
CInterpolation * mpInterpolation
CVector< C_INT > mIWork
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
CVector< C_FLOAT64 > mStateRecord
CVector< C_FLOAT64 > mDWork
CCopasiVector< CMetab > * mpMetabolites
#define C_FLOAT64
Definition: copasi.h:92
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
CModel * getModel() const
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
void calculateAmu(size_t rIndex)
void CHybridMethodODE45::integrateDeterministicPart ( C_FLOAT64  deltaT)
protected

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

Parameters
dsA 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().

871 {
872 
873  //1----Set Parameters for ODE45 solver
874  //=(1)= solver error message
875  mErrorMsg.str("");
876 
877  //=(2)= set state
878  *mpState = mpModel->getState();
879 
880  //=(3)= set time and old time
881  mTime = mpState->getTime();
882  C_FLOAT64 EndTime = mTime + deltaT;
883 
884  mOldTime = mTime;
885 
886  //=(4)= set y and old_y
887  C_FLOAT64 * stateY = mpState->beginIndependent();
888 
889  for (size_t i = 0; i < mData.dim - 1; i++, stateY++)
890  {
891  mOldY[i] = *stateY;
892  mY[i] = *stateY;
893  }
894 
895  if (mODE45Status == NEW_STEP)
896  {
898  mY[mData.dim - 1] = log(rand2);
899  }
900 
901  mOldY[mData.dim - 1] = mY[mData.dim - 1];
902 
903  //2----Reset ODE45 Solver
904  mIFlag = -1; //integration by one step
905 
906  //3----If time increment is too small, do nothing
907  C_FLOAT64 tdist , d__1, d__2, w0;
908  tdist = fabs(deltaT); //absolute time increment
909  d__1 = fabs(mTime), d__2 = fabs(EndTime);
910  w0 = std::max(d__1, d__2);
911 
912  if (tdist < std::numeric_limits< C_FLOAT64 >::epsilon() * 2. * w0) //just do nothing
913  {
914  mpModel->setTime(EndTime);
915  return;
916  }
917 
918  //4----just do nothing if there are no variables
919  if (!mData.dim)
920  {
921  mpModel->setTime(EndTime);
922  return;
923  }
924 
925  while ((mODE45Status == CONTINUE)
926  || (mODE45Status == NEW_STEP))
927  {
928  //(1)----ODE Solver
929  rkf45_(&EvalF , //1. evaluate F
930  &mData.dim, //2. number of variables
931  mY , //3. the array of current concentrations
932  &mTime , //4. the current time
933  &EndTime , //5. the final time
934  &mRtol , //6. relative tolerance array
935  &mAtol , //7. absolute tolerance array
936  &mIFlag , //8. the state for input and output
937  mDWork.array() , //9. the double work array
938  mIWork.array() , //10. the int work array
939  mStateRecord.array()); //11. the array to record temp state
940 
941  //(2)----set mODE45Status
942  if (mIFlag == 2)
943  {
944  if (mY[mData.dim - 1] >= 0.0)
946  else
947  mODE45Status = REACH_END_TIME; //reach EndTime
948  }
949  else if (mIFlag == -2) //finish one step integration
950  {
951  if (mY[mData.dim - 1] >= 0.0) //has an event
953  else
955  }
956  else //error happens
957  {
959  std::cout << "Error Type: " << mIFlag << std::endl;
961  }//end if
962  }//end while
963 
964  //5----Record State
965  stateY = mpState->beginIndependent();
966 
967  for (size_t i = 0; i < mData.dim - 1; i++)
968  stateY[i] = mY[i]; //write result into mpState
969 
972  mpModel->updateSimulatedValues(false); //for assignments?????????
973 
974  return;
975 }
#define CONTINUE
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
#define HAS_EVENT
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
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)
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
#define NEW_STEP
#define MCTrajectoryMethod
void setState(const CState &state)
Definition: CModel.cpp:1785
#define HAS_ERR
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
CVector< C_INT > mIWork
CVector< C_FLOAT64 > mStateRecord
void setTime(const C_FLOAT64 &time)
Definition: CModel.cpp:1187
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
CVector< C_FLOAT64 > mDWork
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
std::ostringstream mErrorMsg
const CState & getState() const
Definition: CModel.cpp:1771
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
#define REACH_END_TIME
#define max(a, b)
Definition: f2c.h:176
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().

104 {
105  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
106 
107  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
108 
109  if (pTP->getDuration() < 0.0)
110  {
111  //back integration not possible
113  return false;
114  }
115 
116  //check for rules
117  size_t i, imax = pTP->getModel()->getNumModelValues();
118 
119  for (i = 0; i < imax; ++i)
120  {
121  if (pTP->getModel()->getModelValues()[i]->getStatus() == CModelEntity::ODE)
122  {
123  //ode rule found
125  return false;
126  }
127  }
128 
129  imax = pTP->getModel()->getNumMetabs();
130 
131  for (i = 0; i < imax; ++i)
132  {
133  if (pTP->getModel()->getMetabolites()[i]->getStatus() == CModelEntity::ODE)
134  {
135  //ode rule found
137  return false;
138  }
139  }
140 
141  imax = pTP->getModel()->getCompartments().size();
142 
143  for (i = 0; i < imax; ++i)
144  {
145  if (pTP->getModel()->getCompartments()[i]->getStatus() == CModelEntity::ODE)
146  {
147  //ode rule found
149  return false;
150  }
151  }
152 
153 /*
154  std::string message = pTP->getModel()->suitableForStochasticSimulation();
155 
156  if (message != "")
157  {
158  //model not suitable, message describes the problem
159  CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
160  return false;
161  }
162 */
163 
164  //events are not supported at the moment
165  if (pTP->getModel()->getEvents().size() > 0)
166  {
168  return false;
169  }
170 
171 
172  return true;
173 }
CCopasiVectorN< CEvent > & getEvents()
Definition: CModel.cpp:1110
const CCopasiVector< CMetab > & getMetabolites() const
Definition: CModel.cpp:1051
const CCopasiVectorN< CModelValue > & getModelValues() const
Definition: CModel.cpp:1060
virtual size_t size() const
size_t getNumMetabs() const
Definition: CModel.cpp:1118
virtual bool isValidProblem(const CCopasiProblem *pProblem)
#define MCTrajectoryMethod
const C_FLOAT64 & getDuration() const
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
size_t getNumModelValues() const
Definition: CModel.cpp:1139
CModel * getModel() const
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().

2870 {
2871  size_t i, imax = pModel->getNumModelValues();
2872 
2873  for (i = 0; i < imax; ++i)
2874  {
2875  if (pModel->getModelValues()[i]->getStatus() == CModelEntity::ASSIGNMENT)
2876  if (pModel->getModelValues()[i]->isUsed())
2877  {
2878  //used assignment found
2879  return true;
2880  }
2881  }
2882 
2883  imax = pModel->getNumMetabs();
2884 
2885  for (i = 0; i < imax; ++i)
2886  {
2887  if (pModel->getMetabolites()[i]->getStatus() == CModelEntity::ASSIGNMENT)
2888  if (pModel->getMetabolites()[i]->isUsed())
2889  {
2890  //used assignment found
2891  return true;
2892  }
2893  }
2894 
2895  imax = pModel->getCompartments().size();
2896 
2897  for (i = 0; i < imax; ++i)
2898  {
2899  if (pModel->getCompartments()[i]->getStatus() == CModelEntity::ASSIGNMENT)
2900  if (pModel->getCompartments()[i]->isUsed())
2901  {
2902  //used assignment found
2903  return true;
2904  }
2905  }
2906 
2907  return false;
2908 }
const CCopasiVector< CMetab > & getMetabolites() const
Definition: CModel.cpp:1051
const CCopasiVectorN< CModelValue > & getModelValues() const
Definition: CModel.cpp:1060
virtual size_t size() const
size_t getNumMetabs() const
Definition: CModel.cpp:1118
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
size_t getNumModelValues() const
Definition: CModel.cpp:1139
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().

2699 {
2700  std::cout << "============Boolean============" << std::endl;
2701 
2702  if (mDoCorrection)
2703  std::cout << "mDoCorrection: Yes" << std::endl;
2704  else
2705  std::cout << "mDoCorrection: No" << std::endl;
2706 
2707  if (mReducedModel)
2708  std::cout << "mReducedModel: Yes" << std::endl;
2709  else
2710  std::cout << "mReducedModel: No" << std::endl;
2711 
2712  std::cout << std::endl;
2713 
2714  std::cout << "============Metab============" << std::endl;
2715 
2716  std::cout << "~~~~mNumVariableMetabs:~~~~ " << mNumVariableMetabs << std::endl;
2717  std::cout << "~~~~mFirstMetabIndex:~~~~ " << mFirstMetabIndex << std::endl;
2718  std::cout << "~~~~mpMetabolites:~~~~" << std::endl;
2719 
2720  for (size_t i = 0; i < mpMetabolites->size(); i++)
2721  {
2722  std::cout << "metab #" << i << " name: " << (*mpMetabolites)[i]->getObjectDisplayName() << std::endl;
2723  std::cout << "value pointer: " << (*mpMetabolites)[i]->getValuePointer() << std::endl;
2724  std::cout << "value: " << *((double *)(*mpMetabolites)[i]->getValuePointer()) << std::endl;
2725  }
2726 
2727  std::cout << std::endl;
2728 
2729  std::cout << "~~~~mMetabFlags:~~~~ " << std::endl;
2730 
2731  for (size_t i = 0; i < mMetabFlags.size(); ++i)
2732  {
2733  std::cout << "metab #" << i << std::endl;
2734  std::cout << "mFlag: " << mMetabFlags[i].mFlag << std::endl;
2735  std::cout << "mFastReactions: ";
2736  std::set<size_t>::iterator it = mMetabFlags[i].mFastReactions.begin();
2737  std::set<size_t>::iterator endIt = mMetabFlags[i].mFastReactions.end();
2738 
2739  for (; it != endIt; it++)
2740  std::cout << *it << " ";
2741 
2742  std::cout << std::endl;
2743  }
2744 
2745  std::cout << std::endl;
2746 
2747  std::cout << "============Reaction============" << std::endl;
2748 
2749  std::cout << "~~~~mNumReactions:~~~~ " << mNumReactions << std::endl;
2750 
2751  for (size_t i = 0; i < mNumReactions; ++i)
2752  {
2753  std::cout << "Reaction #: " << i << " Flag: " << mReactionFlags[i] << std::endl;
2754  }
2755 
2756  std::cout << "~~~~mLocalBalances:~~~~ " << std::endl;
2757 
2758  for (size_t i = 0; i < mLocalBalances.size(); ++i)
2759  {
2760  std::cout << "Reaction: " << i << std::endl;
2761 
2762  for (size_t j = 0; j < mLocalBalances[i].size(); ++j)
2763  {
2764  std::cout << "Index: " << mLocalBalances[i][j].mIndex << std::endl;
2765  std::cout << "mMultiplicity: " << mLocalBalances[i][j].mMultiplicity << std::endl;
2766  std::cout << "mpMetablite: " << mLocalBalances[i][j].mpMetabolite << std::endl;
2767  }
2768  }
2769 
2770  std::cout << "~~~~mLocalSubstrates:~~~~ " << std::endl;
2771 
2772  for (size_t i = 0; i < mLocalSubstrates.size(); ++i)
2773  {
2774  std::cout << "Reaction: " << i << std::endl;
2775 
2776  for (size_t j = 0; j < mLocalSubstrates[i].size(); ++j)
2777  {
2778  std::cout << "Index: " << mLocalSubstrates[i][j].mIndex << std::endl;
2779  std::cout << "mMultiplicity: " << mLocalSubstrates[i][j].mMultiplicity << std::endl;
2780  std::cout << "mpMetablite: " << mLocalSubstrates[i][j].mpMetabolite << std::endl;
2781  }
2782  }
2783 
2784  std::cout << "~~~~mMetab2React:~~~~ " << std::endl;
2785 
2786  for (size_t i = 0; i < mMetab2React.size(); i++)
2787  {
2788  std::cout << "metab #: " << i << std::endl;
2789  std::set<size_t>::iterator it = mMetab2React[i].begin();
2790  std::set<size_t>::iterator endIt = mMetab2React[i].end();
2791  std::cout << "React: ";
2792 
2793  for (; it != endIt; it++)
2794  std::cout << *it << " ";
2795 
2796  std::cout << std::endl;
2797  }
2798 
2799  std::cout << std::endl;
2800  std::cout << "~~~~mReactAffect:~~~~ " << std::endl;
2801 
2802  for (size_t i = 0; i < mReactAffect.size(); i++)
2803  {
2804  std::cout << "react #: " << i << std::endl;
2805  std::set<size_t>::iterator it = mReactAffect[i].begin();
2806  std::set<size_t>::iterator endIt = mReactAffect[i].end();
2807  std::cout << "affect: ";
2808 
2809  for (; it != endIt; it++)
2810  std::cout << *it << " ";
2811 
2812  std::cout << std::endl;
2813  }
2814 
2815  if (mHasStoiReaction)
2816  std::cout << "mHasStoiReaction: Yes" << std::endl;
2817  else
2818  std::cout << "mHasStoiReaction: No" << std::endl;
2819 
2820  if (mHasDetermReaction)
2821  std::cout << "mHasDetermReaction: Yes" << std::endl;
2822  else
2823  std::cout << "mHasDetermReaction: No" << std::endl;
2824 
2825  getchar();
2826  return;
2827 }
virtual size_t size() const
std::vector< std::set< size_t > > mMetab2React
std::vector< CHybridODE45MetabFlag > mMetabFlags
std::vector< std::vector< CHybridODE45Balance > > mLocalSubstrates
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
CCopasiVector< CMetab > * mpMetabolites
std::vector< std::set< size_t > > mReactAffect
CVector< size_t > mReactionFlags
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.

2477 {
2478  size_t i, j;
2479  std::set< size_t >::iterator iter, iterEnd;
2480 
2481  os << "outputDebug(" << level << ") *********************************************** BEGIN" << std::endl;
2482 
2483  switch (level)
2484  {
2485  case 0: // Everything !!!
2486  os << " Name: " << CCopasiParameter::getObjectName() << std::endl;
2487  os << "current time: " << mpCurrentState->getTime() << std::endl;
2488  os << "mNumVariableMetabs: " << mNumVariableMetabs << std::endl;
2489  os << "mMaxSteps: " << mMaxSteps << std::endl;
2490  os << "mMaxBalance: " << mMaxBalance << std::endl;
2491  //os << "mMaxIntBeforeStep: " << mMaxIntBeforeStep << std::endl;
2492  os << "mpReactions.size(): " << mpReactions->size() << std::endl;
2493 
2494  for (i = 0; i < mpReactions->size(); i++)
2495  os << *(*mpReactions)[i] << std::endl;
2496 
2497  os << "mpMetabolites.size(): " << mpMetabolites->size() << std::endl;
2498 
2499  for (i = 0; i < mpMetabolites->size(); i++)
2500  os << *(*mpMetabolites)[i] << std::endl;
2501 
2502  os << "mStoi: " << std::endl;
2503 
2504  for (i = 0; i < mStoi.numRows(); i++)
2505  {
2506  for (j = 0; j < mStoi.numCols(); j++)
2507  os << mStoi[i][j] << " ";
2508 
2509  os << std::endl;
2510  }
2511 
2512  os << "temp: ";
2513 
2514  for (i = 0; i < mNumVariableMetabs; i++)
2515  os << temp[i] << " ";
2516 
2517  os << std::endl;
2518 
2519  os << "mReactionFlags: " << std::endl;
2520 
2521  for (i = 0; i < mLocalBalances.size(); i++)
2522  os << mReactionFlags[i];
2523 
2524  //os << "mFirstReactionFlag: " << std::endl;
2525  // if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
2526 
2527  os << "mMetabFlags: " << std::endl;
2528 
2529  for (i = 0; i < mMetabFlags.size(); i++)
2530  {
2531  if (mMetabFlags[i].mFlag == SLOW)
2532  os << "SLOW ";
2533  else
2534  os << "FAST ";
2535  }
2536 
2537  os << std::endl;
2538  os << "mLocalBalances: " << std::endl;
2539 
2540  /*
2541  for (i = 0; i < mLocalBalances.size(); i++)
2542  {
2543  for (j = 0; j < mLocalBalances[i].size(); j++)
2544  os << mLocalBalances[i][j];
2545 
2546  os << std::endl;
2547  }
2548 
2549  os << "mLocalSubstrates: " << std::endl;
2550 
2551  for (i = 0; i < mLocalSubstrates.size(); i++)
2552  {
2553  for (j = 0; j < mLocalSubstrates[i].size(); j++)
2554  os << mLocalSubstrates[i][j];
2555 
2556  os << std::endl;
2557  }
2558  */
2559  //os << "mLowerStochLimit: " << mLowerStochLimit << std::endl;
2560  //os << "mUpperStochLimit: " << mUpperStochLimit << std::endl;
2561  //deprecated: os << "mOutputCounter: " << mOutputCounter << endl;
2562  //os << "mPartitioningInterval: " << mPartitioningInterval << std::endl;
2563  //os << "mStepsAfterPartitionSystem: " << mStepsAfterPartitionSystem << std::endl;
2564  //os << "mStepsize: " << mStepsize << std::endl;
2565  os << "mMetab2React: " << std::endl;
2566 
2567  for (i = 0; i < mMetab2React.size(); i++)
2568  {
2569  os << i << ": ";
2570 
2571  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; ++iter)
2572  os << *iter << " ";
2573 
2574  os << std::endl;
2575  }
2576 
2577  os << "mAmu: " << std::endl;
2578 
2579  for (i = 0; i < mpReactions->size(); i++)
2580  os << mAmu[i] << " ";
2581 
2582  os << std::endl;
2583  os << "mAmuOld: " << std::endl;
2584 
2585  for (i = 0; i < mpReactions->size(); i++)
2586  os << mAmuOld[i] << " ";
2587 
2588  os << std::endl;
2589  os << "mUpdateSet: " << std::endl;
2590 
2591  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
2592  os << *iter;
2593 
2594  os << std::endl;
2595  os << "mpRandomGenerator: " << mpRandomGenerator << std::endl;
2596  os << "mDG: " << std::endl << mDG;
2597  os << "mPQ: " << std::endl << mPQ;
2598  os << "Particle numbers: " << std::endl;
2599 
2600  for (i = 0; i < mpMetabolites->size(); i++)
2601  {
2602  os << (*mpMetabolites)[i]->getValue() << " ";
2603  }
2604 
2605  os << std::endl;
2606  break;
2607 
2608  case 1: // Variable values only
2609  os << "current time: " << mpCurrentState->getTime() << std::endl;
2610  /*
2611  case 1:
2612  os << "mTime: " << mpCurrentState->getTime() << std::endl;
2613  os << "oldState: ";
2614  for (i = 0; i < mDim; i++)
2615  os << oldState[i] << " ";
2616  os << std::endl;
2617  os << "x: ";
2618  for (i = 0; i < mDim; i++)
2619  os << x[i] << " ";
2620  os << std::endl;
2621  os << "y: ";
2622  for (i = 0; i < mDim; i++)
2623  os << y[i] << " ";
2624  os << std::endl;
2625  os << "increment: ";
2626  for (i = 0; i < mDim; i++)
2627  os << increment[i] << " ";
2628  os << std::endl;*/
2629  os << "temp: ";
2630 
2631  for (i = 0; i < mNumVariableMetabs; i++)
2632  os << temp[i] << " ";
2633 
2634  os << std::endl;
2635 
2636  os << "mReactionFlags: " << std::endl;
2637 
2638  for (i = 0; i < mLocalBalances.size(); i++)
2639  os << mReactionFlags[i];
2640 
2641  // os << "mFirstReactionFlag: " << std::endl;
2642  //if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
2643 
2644  os << "mMetabFlags: " << std::endl;
2645 
2646  for (i = 0; i < mMetabFlags.size(); i++)
2647  {
2648  if (mMetabFlags[i].mFlag == SLOW)
2649  os << "SLOW ";
2650  else
2651  os << "FAST ";
2652  }
2653 
2654  os << std::endl;
2655  os << "mAmu: " << std::endl;
2656 
2657  for (i = 0; i < mpReactions->size(); i++)
2658  os << mAmu[i] << " ";
2659 
2660  os << std::endl;
2661  os << "mAmuOld: " << std::endl;
2662 
2663  for (i = 0; i < mpReactions->size(); i++)
2664  os << mAmuOld[i] << " ";
2665 
2666  os << std::endl;
2667  os << "mUpdateSet: " << std::endl;
2668 
2669  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
2670  os << *iter;
2671 
2672  os << std::endl;
2673  os << "mPQ: " << std::endl << mPQ;
2674  os << "Particle numbers: " << std::endl;
2675 
2676  for (i = 0; i < mpMetabolites->size(); i++)
2677  {
2678  os << (*mpMetabolites)[i]->getValue() << " ";
2679  }
2680 
2681  os << std::endl;
2682  break;
2683 
2684  case 2:
2685  break;
2686 
2687  default:
2688  ;
2689  }
2690 
2691  os << "outputDebug(" << level << ") ************************************************* END" << std::endl;
2692  return;
2693 }
std::set< size_t > mUpdateSet
CVector< C_FLOAT64 > temp
const CCopasiVectorNS< CReaction > * mpReactions
const std::string & getObjectName() const
virtual size_t size() const
CDependencyGraph mDG
std::vector< std::set< size_t > > mMetab2React
virtual size_t numRows() const
Definition: CMatrix.h:138
CVector< C_FLOAT64 > mAmuOld
CMatrix< C_FLOAT64 > mStoi
std::vector< CHybridODE45MetabFlag > mMetabFlags
CVector< C_FLOAT64 > mAmu
#define SLOW
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
CIndexedPriorityQueue mPQ
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
CCopasiVector< CMetab > * mpMetabolites
virtual size_t numCols() const
Definition: CMatrix.h:144
CVector< size_t > mReactionFlags
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().

2833 {
2834  const C_FLOAT64 time = mpState->getTime();
2835  std::cout << "**State Output**" << std::endl;
2836  std::cout << "Time: " << time << std::endl;
2837 
2838  std::cout << "Indep #: " << mpState->getNumIndependent() << " Id: ";
2839  const C_FLOAT64 *pIt = mpState->beginIndependent();
2840  const C_FLOAT64 *pEnd = mpState->endIndependent();
2841 
2842  for (; pIt != pEnd; pIt++)
2843  std::cout << *pIt << " ";
2844 
2845  std::cout << std::endl;
2846 
2847  std::cout << "Dep #: " << mpState->getNumDependent() << " Id: ";
2848  pIt = mpState->beginDependent();
2849  pEnd = mpState->endDependent();
2850 
2851  for (; pIt != pEnd; pIt++)
2852  std::cout << *pIt << " ";
2853 
2854  std::cout << std::endl;
2855 
2856  std::cout << "Fix #: " << mpState->getNumIndependent() << " Id: ";
2857  pIt = mpState->beginFixed();
2858  pEnd = mpState->endFixed();
2859 
2860  for (; pIt != pEnd; pIt++)
2861  std::cout << *pIt << " ";
2862 
2863  std::cout << std::endl;
2864 
2865  getchar();
2866  return;
2867 }
C_FLOAT64 * endDependent()
Definition: CState.cpp:331
C_FLOAT64 * endIndependent()
Definition: CState.cpp:329
C_FLOAT64 * endFixed()
Definition: CState.cpp:333
size_t getNumDependent() const
Definition: CState.cpp:344
size_t getNumIndependent() const
Definition: CState.cpp:342
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 * beginFixed()
Definition: CState.cpp:332
C_FLOAT64 * beginDependent()
Definition: CState.cpp:330
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
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().

1480 {
1481  static C_INT k1, k2, k3, k4, k5, k6, k1m;
1482 
1483  /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
1484  /* This is a modefied fehlberg fourth-fifth order runge-kutta method */
1485  /* based on h.a.watts and l.f.shapine's fortran ODE solver rkf45.f */
1486  /* with one more input double precision vector recording time */
1487  /* points between t and t+h in each steap and */
1488  /* approximations of y in such corresponding times. */
1489  /* yrcd is a vector of length 4*(neqn+1), where "4" is related with */
1490  /* four middle time points: */
1491  /* c1=t+h/4, c2=t+3h/8, c3=t+h/2, c4=t+12h/13. */
1492  /* yrcd can be seperated into four parts, each of which is a double */
1493  /* precision vector of length neqn+1. The first member is assigned */
1494  /* to time record, and the last neqn members are for y at the */
1495  /* corresponding time. */
1496  /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
1497 
1498  /* fehlberg fourth-fifth order runge-kutta method */
1499 
1500  /* written by h.a.watts and l.f.shampine */
1501  /* sandia laboratories */
1502  /* albuquerque,new mexico */
1503 
1504  /* rkf45 is primarily designed to solve non-stiff and mildly stiff */
1505  /* differential equations when derivative evaluations are inexpensive. */
1506  /* rkf45 should generally not be used when the user is demanding */
1507  /* high accuracy. */
1508 
1509  /* abstract */
1510 
1511  /* subroutine rkf45 integrates a system of neqn first order */
1512  /* ordinary differential equations of the form */
1513  /* dy(i)/dt = f(t,y(1),y(2),...,y(neqn)) */
1514  /* where the y(i) are given at t . */
1515  /* typically the subroutine is used to integrate from t to tout but it */
1516  /* can be used as a one-step integrator to advance the solution a */
1517  /* single step in the direction of tout. on return the parameters in */
1518  /* the call list are set for continuing the integration. the user has */
1519  /* only to call rkf45 again (and perhaps define a new value for tout). */
1520  /* actually, rkf45 is an interfacing routine which calls subroutine */
1521  /* rkfs for the solution. rkfs in turn calls subroutine fehl which */
1522  /* computes an approximate solution over one step. */
1523 
1524  /* rkf45 uses the runge-kutta-fehlberg (4,5) method described */
1525  /* in the reference */
1526  /* e.fehlberg , low-order classical runge-kutta formulas with stepsize */
1527  /* control , nasa tr r-315 */
1528 
1529  /* the performance of rkf45 is illustrated in the reference */
1530  /* l.f.shampine,h.a.watts,s.davenport, solving non-stiff ordinary */
1531  /* differential equations-the state of the art , */
1532  /* sandia laboratories report sand75-0182 , */
1533  /* to appear in siam review. */
1534 
1535  /* the parameters represent- */
1536  /* f -- subroutine f(t,y,yp) to evaluate derivatives yp(i)=dy(i)/dt */
1537  /* neqn -- number of equations to be integrated */
1538  /* y(*) -- solution vector at t */
1539  /* t -- independent variable */
1540  /* tout -- output point at which solution is desired */
1541  /* relerr,abserr -- relative and absolute error tolerances for local */
1542  /* error test. at each step the code requires that */
1543  /* abs(local error) .le. relerr*abs(y) + abserr */
1544  /* for each component of the local error and solution vectors */
1545  /* iflag -- indicator for status of integration */
1546  /* work(*) -- array to hold information internal to rkf45 which is */
1547  /* necessary for subsequent calls. must be dimensioned */
1548  /* at least 3+6*neqn */
1549  /* iwork(*) -- C_INT array used to hold information internal to */
1550  /* rkf45 which is necessary for subsequent calls. must be */
1551  /* dimensioned at least 5 */
1552 
1553  /* first call to rkf45 */
1554 
1555  /* the user must provide storage in his calling program for the arrays */
1556  /* in the call list - y(neqn) , work(3+6*neqn) , iwork(5) , */
1557  /* declare f in an external statement, supply subroutine f(t,y,yp) and */
1558  /* initialize the following parameters- */
1559 
1560  /* neqn -- number of equations to be integrated. (neqn .ge. 1) */
1561  /* y(*) -- vector of initial conditions */
1562  /* t -- starting point of integration , must be a variable */
1563  /* tout -- output point at which solution is desired. */
1564  /* t=tout is allowed on the first call only, in which case */
1565  /* rkf45 returns with iflag=2 if continuation is possible. */
1566  /* relerr,abserr -- relative and absolute local error tolerances */
1567  /* which must be non-negative. relerr must be a variable while */
1568  /* abserr may be a constant. the code should normally not be */
1569  /* used with relative error control smaller than about 1.e-8 . */
1570  /* to avoid limiting precision difficulties the code requires */
1571  /* relerr to be larger than an internally computed relative */
1572  /* error parameter which is machine dependent. in particular, */
1573  /* pure absolute error is not permitted. if a smaller than */
1574  /* allowable value of relerr is attempted, rkf45 increases */
1575  /* relerr appropriately and returns control to the user before */
1576  /* continuing the integration. */
1577  /* iflag -- +1,-1 indicator to initialize the code for each new */
1578  /* problem. normal input is +1. the user should set iflag=-1 */
1579  /* only when one-step integrator control is essential. in this */
1580  /* case, rkf45 attempts to advance the solution a single step */
1581  /* in the direction of tout each time it is called. since this */
1582  /* mode of operation results in extra computing overhead, it */
1583  /* should be avoided unless needed. */
1584 
1585  /* output from rkf45 */
1586 
1587  /* y(*) -- solution at t */
1588  /* t -- last point reached in integration. */
1589  /* iflag = 2 -- integration reached tout. indicates successful retur */
1590  /* and is the normal mode for continuing integration. */
1591  /* =-2 -- a single successful step in the direction of tout */
1592  /* has been taken. normal mode for continuing */
1593  /* integration one step at a time. */
1594  /* = 3 -- integration was not completed because relative error */
1595  /* tolerance was too small. relerr has been increased */
1596  /* appropriately for continuing. */
1597  /* = 4 -- integration was not completed because more than */
1598  /* 3000 derivative evaluations were needed. this */
1599  /* is approximately 500 steps. */
1600  /* = 5 -- integration was not completed because solution */
1601  /* vanished making a pure relative error test */
1602  /* impossible. must use non-zero abserr to continue. */
1603  /* using the one-step integration mode for one step */
1604  /* is a good way to proceed. */
1605  /* = 6 -- integration was not completed because requested */
1606  /* accuracy could not be achieved using smallest */
1607  /* allowable stepsize. user must increase the error */
1608  /* tolerance before continued integration can be */
1609  /* attempted. */
1610  /* = 7 -- it is likely that rkf45 is inefficient for solving */
1611  /* this problem. too much output is restricting the */
1612  /* natural stepsize choice. use the one-step integrator */
1613  /* mode. */
1614  /* = 8 -- invalid input parameters */
1615  /* this indicator occurs if any of the following is */
1616  /* satisfied - neqn .le. 0 */
1617  /* t=tout and iflag .ne. +1 or -1 */
1618  /* relerr or abserr .lt. 0. */
1619  /* iflag .eq. 0 or .lt. -2 or .gt. 8 */
1620  /* work(*),iwork(*) -- information which is usually of no interest */
1621  /* to the user but necessary for subsequent calls. */
1622  /* work(1),...,work(neqn) contain the first derivatives */
1623  /* of the solution vector y at t. work(neqn+1) contains */
1624  /* the stepsize h to be attempted on the next step. */
1625  /* iwork(1) contains the derivative evaluation counter. */
1626 
1627  /* subsequent calls to rkf45 */
1628 
1629  /* subroutine rkf45 returns with all information needed to continue */
1630  /* the integration. if the integration reached tout, the user need onl */
1631  /* define a new tout and call rkf45 again. in the one-step integrator */
1632  /* mode (iflag=-2) the user must keep in mind that each step taken is */
1633  /* in the direction of the current tout. upon reaching tout (indicated */
1634  /* by changing iflag to 2),the user must then define a new tout and */
1635  /* reset iflag to -2 to continue in the one-step integrator mode. */
1636 
1637  /* if the integration was not completed but the user still wants to */
1638  /* continue (iflag=3,4 cases), he just calls rkf45 again. with iflag=3 */
1639  /* the relerr parameter has been adjusted appropriately for continuing */
1640  /* the integration. in the case of iflag=4 the function counter will */
1641  /* be reset to 0 and another 3000 function evaluations are allowed. */
1642 
1643  /* however,in the case iflag=5, the user must first alter the error */
1644  /* criterion to use a positive value of abserr before integration can */
1645  /* proceed. if he does not,execution is terminated. */
1646 
1647  /* also,in the case iflag=6, it is necessary for the user to reset */
1648  /* iflag to 2 (or -2 when the one-step integration mode is being used) */
1649  /* as well as increasing either abserr,relerr or both before the */
1650  /* integration can be continued. if this is not done, execution will */
1651  /* be terminated. the occurrence of iflag=6 indicates a trouble spot */
1652  /* (solution is changing rapidly,singularity may be present) and it */
1653  /* often is inadvisable to continue. */
1654 
1655  /* if iflag=7 is encountered, the user should use the one-step */
1656  /* integration mode with the stepsize determined by the code or */
1657  /* consider switching to the adams codes de/step,intrp. if the user */
1658  /* insists upon continuing the integration with rkf45, he must reset */
1659  /* iflag to 2 before calling rkf45 again. otherwise,execution will be */
1660  /* terminated. */
1661 
1662  /* if iflag=8 is obtained, integration can not be continued unless */
1663  /* the invalid input parameters are corrected. */
1664 
1665  /* it should be noted that the arrays work,iwork contain information */
1666  /* required for subsequent integration. accordingly, work and iwork */
1667  /* should not be altered. */
1668 
1669  /* compute indices for the splitting of the work array */
1670 
1671  /* Parameter adjustments */
1672  --yrcd;
1673  --y;
1674  --work;
1675  --iwork;
1676 
1677  /* Function Body */
1678  k1m = *neqn + 1;
1679  //k1m = *neqn;
1680  k1 = k1m + 1;
1681  k2 = k1 + *neqn;
1682  k3 = k2 + *neqn;
1683  k4 = k3 + *neqn;
1684  k5 = k4 + *neqn;
1685  k6 = k5 + *neqn;
1686 
1687  /* this interfacing routine merely relieves the user of a long */
1688  /* calling list via the splitting apart of two working storage */
1689  /* arrays. if this is not compatible with the users compiler, */
1690  /* he must use rkfs directly. */
1691 
1692  rkfs_((pEvalF)f, neqn, &y[1], t, tout, relerr, abserr, iflag,
1693  &work[1], &work[k1m], &work[k1], &work[k2], &work[k3],
1694  &work[k4], &work[k5], &work[k6], &work[k6 + 1], &iwork[1],
1695  &iwork[2], &iwork[3], &iwork[4], &iwork[5], &yrcd[1]);
1696 
1697  return 0;
1698 } /* rkf45_ */
#define C_INT
Definition: copasi.h:115
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)
void(* pEvalF)(const C_INT *, const double *, const double *, double *)
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_().

1706 {
1707  /* Initialized data */
1708 
1709  static double remin = 1e-12;
1710  static C_INT maxnfe = 3000;
1711 
1712  /* System generated locals */
1713  C_INT i__1;
1714  double d__1, d__2, d__3, d__4;
1715 
1716  /* Builtin functions */
1717  /* Subroutine */
1718  //int s_stop(char *, ftnlen);
1719  //double pow_dd(double *, double *);
1720  //C_INT i_sign(C_INT *, C_INT *);
1721  //double d_sign(double *, double *);
1722 
1723  /* Local variables */
1724  static double a;
1725  static C_INT k;
1726  static double s, ae, ee, dt, et, u26, rer, tol, ypk;
1727  static double hmin, toln;
1728  static C_INT mflag;
1729  static double scale, eeoet;
1730  static int hfaild;
1731  static double esttol, twoeps;
1732  static int output;
1733 
1734  /* fehlberg fourth-fifth order runge-kutta method */
1735 
1736  /* rkfs integrates a system of first order ordinary differential */
1737  /* equations as described in the comments for rkf45 . */
1738  /* the arrays yp,f1,f2,f3,f4,and f5 (of dimension at least neqn) and */
1739  /* the variables h,savre,savae,nfe,kop,init,jflag,and kflag are used */
1740  /* internally by the code and appear in the call list to eliminate */
1741  /* local retention of variables between calls. accordingly, they */
1742  /* should not be altered. items of possible interest are */
1743  /* yp - derivative of solution vector at t */
1744  /* h - an appropriate stepsize to be used for the next step */
1745  /* nfe- counter on the number of derivative function evaluations */
1746 
1747  /* remin is the minimum acceptable value of relerr. attempts */
1748  /* to obtain higher accuracy with this subroutine are usually */
1749  /* very expensive and often unsuccessful. */
1750 
1751  /* Parameter adjustments */
1752  --yrcd;
1753  --f5;
1754  --f4;
1755  --f3;
1756  --f2;
1757  --f1;
1758  --yp;
1759  --y;
1760 
1761  /* Function Body */
1762 
1763  /* the expense is controlled by restricting the number */
1764  /* of function evaluations to be approximately maxnfe. */
1765  /* as set, this corresponds to about 500 steps. */
1766 
1767  /* here two constants emboding the machine epsilon is present */
1768  /* twoesp is set to twice the machine epsilon while u26 is set */
1769  /* to 26 times the machine epsilon */
1770 
1771  /* data twoeps, u26/4.4d-16, 5.72d-15/ *** */
1772  twoeps = 2.f * std::numeric_limits< C_FLOAT64 >::epsilon();
1773  //twoeps = 4.4e-16;
1774  u26 = twoeps * 13.f;
1775 
1776  /* check input parameters */
1777  if (*neqn < 1)
1778  {
1779  goto L10;
1780  }
1781 
1782  if (*relerr < 0. || *abserr < 0.)
1783  {
1784  goto L10;
1785  }
1786 
1787  mflag = abs(*iflag);
1788 
1789  if (mflag >= 1 && mflag <= 8)
1790  {
1791  goto L20;
1792  }
1793 
1794  /* invalid input */
1795 L10:
1796  *iflag = 8;
1797  return 0;
1798 
1799  /* is this the first call */
1800 L20:
1801 
1802  if (mflag == 1)
1803  {
1804  goto L50;
1805  }
1806 
1807  /* check continuation possibilities */
1808 
1809  if (*t == *tout && *kflag != 3)
1810  {
1811  goto L10;
1812  }
1813 
1814  if (mflag != 2)
1815  {
1816  goto L25;
1817  }
1818 
1819  /* iflag = +2 or -2 */
1820  if (*kflag == 3)
1821  {
1822  goto L45;
1823  }
1824 
1825  if (*init == 0)
1826  {
1827  goto L45;
1828  }
1829 
1830  if (*kflag == 4)
1831  {
1832  goto L40;
1833  }
1834 
1835  if (*kflag == 5 && *abserr == 0.)
1836  {
1837  goto L30;
1838  }
1839 
1840  if (*kflag == 6 && *relerr <= *savre && *abserr <= *savae)
1841  {
1842  goto L30;
1843  }
1844 
1845  goto L50;
1846 
1847  /* iflag = 3,4,5,6,7 or 8 */
1848 L25:
1849 
1850  if (*iflag == 3)
1851  {
1852  goto L45;
1853  }
1854 
1855  if (*iflag == 4)
1856  {
1857  goto L40;
1858  }
1859 
1860  if (*iflag == 5 && *abserr > 0.)
1861  {
1862  goto L45;
1863  }
1864 
1865  /* integration cannot be continued since user did not respond to */
1866  /* the instructions pertaining to iflag=5,6,7 or 8 */
1867 L30:
1868  //s_stop("", (ftnlen)0);
1869  std::cout << "STOP 0 statement executed" << std::endl;
1870  return -1;
1871 
1872  /* reset function evaluation counter */
1873 L40:
1874  *nfe = 0;
1875 
1876  if (mflag == 2)
1877  {
1878  goto L50;
1879  }
1880 
1881  /* reset flag value from previous call */
1882 L45:
1883  *iflag = *jflag;
1884 
1885  if (*kflag == 3)
1886  {
1887  mflag = abs(*iflag);
1888  }
1889 
1890  /* save input iflag and set continuation flag value for subsequent */
1891  /* input checking */
1892 L50:
1893  *jflag = *iflag;
1894  *kflag = 0;
1895 
1896  /* save relerr and abserr for checking input on subsequent calls */
1897  *savre = *relerr;
1898  *savae = *abserr;
1899 
1900  /* restrict relative error tolerance to be at least as large as */
1901  /* 2*eps+remin to avoid limiting precision difficulties arising */
1902  /* from impossible accuracy requests */
1903 
1904  rer = twoeps + remin;
1905 
1906  if (*relerr >= rer)
1907  {
1908  goto L55;
1909  }
1910 
1911  /* relative error tolerance too small */
1912  *relerr = rer;
1913  *iflag = 3;
1914  *kflag = 3;
1915  return 0;
1916 
1917 L55:
1918  dt = *tout - *t;
1919 
1920  if (mflag == 1)
1921  {
1922  goto L60;
1923  }
1924 
1925  if (*init == 0)
1926  {
1927  goto L65;
1928  }
1929 
1930  goto L80;
1931 
1932  /* initialization -- */
1933  /* set initialization completion indicator,init */
1934  /* set indicator for too many output points,kop */
1935  /* evaluate initial derivatives */
1936  /* set counter for function evaluations,nfe */
1937  /* evaluate initial derivatives */
1938  /* set counter for function evaluations,nfe */
1939  /* estimate starting stepsize */
1940 
1941 L60:
1942  *init = 0;
1943  *kop = 0;
1944 
1945  a = *t;
1946  (*f)(neqn, &a, &y[1], &yp[1]);
1947  *nfe = 1;
1948 
1949  if (*t != *tout)
1950  {
1951  goto L65;
1952  }
1953 
1954  *iflag = 2;
1955  return 0;
1956 
1957 L65:
1958  *init = 1;
1959  *h__ = fabs(dt);
1960  toln = 0.f;
1961  i__1 = *neqn;
1962 
1963  for (k = 1; k <= i__1; ++k)
1964  {
1965  tol = *relerr * (d__1 = y[k], fabs(d__1)) + *abserr;
1966 
1967  if (tol <= 0.f)
1968  {
1969  goto L70;
1970  }
1971 
1972  toln = tol;
1973  ypk = (d__1 = yp[k], fabs(d__1));
1974  /* Computing 5th power */
1975  d__1 = *h__, d__2 = d__1, d__1 *= d__1;
1976 
1977  if (ypk * (d__2 * (d__1 * d__1)) > tol)
1978  {
1979  d__3 = tol / ypk;
1980  *h__ = pow(d__3, .2);
1981  }
1982 
1983 L70:
1984  ;
1985  }
1986 
1987  if (toln <= 0.)
1988  {
1989  *h__ = 0.;
1990  }
1991 
1992  /* Computing MAX */
1993  /* Computing MAX */
1994  d__3 = fabs(*t), d__4 = fabs(dt);
1995  d__1 = *h__, d__2 = u26 * std::max(d__3, d__4);
1996  *h__ = std::max(d__1, d__2);
1997  //*jflag = i_sign(2, iflag);
1998  *jflag = *iflag >= 0 ? 2 : -2;
1999 
2000  /* set stepsize for integration in the direction from t to tout */
2001 
2002 L80:
2003  //*h__ = d_sign(h__, &dt);
2004  *h__ = dt >= 0 ? fabs(*h__) : -fabs(*h__);
2005 
2006  /* test to see if rkf45 is being severely impacted by too many */
2007  /* output points */
2008 
2009  if (fabs(*h__) >= fabs(dt) * 2.)
2010  {
2011  ++(*kop);
2012  }
2013 
2014  if (*kop != 100)
2015  {
2016  goto L85;
2017  }
2018 
2019  /* unnecessary frequency of output */
2020  *kop = 0;
2021  *iflag = 7;
2022  return 0;
2023 
2024 L85:
2025 
2026  if (fabs(dt) > u26 * fabs(*t))
2027  {
2028  goto L95;
2029  }
2030 
2031  /* if too close to output point,extrapolate and return */
2032  i__1 = *neqn;
2033 
2034  for (k = 1; k <= i__1; ++k)
2035  {
2036  /* L90: */
2037  y[k] += dt * yp[k];
2038  }
2039 
2040  a = *tout;
2041  (*f)(neqn, &a, &y[1], &yp[1]);
2042  ++(*nfe);
2043 
2044  std::cout << "Step too Small" << std::endl;
2045  mUseStateRecord = false;
2046  goto L300;
2047 
2048  /* initialize output point indicator */
2049 
2050 L95:
2051  output = 0;
2052 
2053  /* to avoid premature underflow in the error tolerance function, */
2054  /* scale the error tolerances */
2055 
2056  scale = 2. / *relerr;
2057  ae = scale * *abserr;
2058 
2059  /* step by step integration */
2060 
2061 L100:
2062  hfaild = 0;
2063 
2064  /* set smallest allowable stepsize */
2065  hmin = u26 * fabs(*t);
2066 
2067  /* adjust stepsize if necessary to hit the output point. */
2068  /* look ahead two steps to avoid drastic changes in the stepsize and */
2069  /* thus lessen the impact of output points on the code. */
2070 
2071  dt = *tout - *t;
2072 
2073  if (fabs(dt) >= fabs(*h__) * 2.)
2074  {
2075  goto L200;
2076  }
2077 
2078  if (fabs(dt) > fabs(*h__))
2079  {
2080  goto L150;
2081  }
2082 
2083  /* the next successful step will complete the integration to the */
2084  /* output point */
2085 
2086  output = 1;
2087  *h__ = dt;
2088  goto L200;
2089 
2090 L150:
2091  *h__ = dt * .5;
2092 
2093  /* core integrator for taking a single step */
2094 
2095  /* the tolerances have been scaled to avoid premature underflow in */
2096  /* computing the error tolerance function et. */
2097  /* to avoid problems with zero crossings,relative error is measured */
2098  /* using the average of the magnitudes of the solution at the */
2099  /* beginning and end of a step. */
2100  /* the error estimate formula has been grouped to control loss of */
2101  /* significance. */
2102  /* to distinguish the various arguments, h is not permitted */
2103  /* to become smaller than 26 units of roundoff in t. */
2104  /* practical limits on the change in the stepsize are enforced to */
2105  /* smooth the stepsize selection process and to avoid excessive */
2106  /* chattering on problems having discontinuities. */
2107  /* to prevent unnecessary failures, the code uses 9/10 the stepsize */
2108  /* it estimates will succeed. */
2109  /* after a step failure, the stepsize is not allowed to increase for */
2110  /* the next attempted step. this makes the code more efficient on */
2111  /* problems having discontinuities and more effective in general */
2112  /* since local extrapolation is being used and extra caution seems */
2113  /* warranted. */
2114 
2115  /* test number of derivative function evaluations. */
2116  /* if okay,try to advance the integration from t to t+h */
2117 
2118 L200:
2119 
2120  if (*nfe <= maxnfe)
2121  {
2122  goto L220;
2123  }
2124 
2125  /* too much work */
2126  *iflag = 4;
2127  *kflag = 4;
2128  return 0;
2129 
2130  /* advance an approximate solution over one step of length h */
2131 
2132 L220:
2133  fehl_((pEvalF)f, neqn, &y[1], t, h__, &yp[1], &f1[1], &f2[1],
2134  &f3[1], &f4[1], &f5[1], &f1[1], &yrcd[1]);
2135  *nfe += 5;
2136  mUseStateRecord = true;
2137 
2138  /* compute and test allowable tolerances versus local error estimates*/
2139  /* and remove scaling of tolerances. note that relative error is */
2140  /* measured with respect to the average of the magnitudes of the */
2141  /* solution at the beginning and end of the step. */
2142 
2143  eeoet = 0.;
2144  i__1 = *neqn;
2145 
2146  for (k = 1; k <= i__1; ++k)
2147  {
2148  et = (d__1 = y[k], fabs(d__1)) + (d__2 = f1[k], fabs(d__2)) + ae;
2149 
2150  if (et > 0.)
2151  {
2152  goto L240;
2153  }
2154 
2155  /* inappropriate error tolerance */
2156  *iflag = 5;
2157  return 0;
2158 
2159 L240:
2160  ee = (d__1 = yp[k] * -2090. + (f3[k] * 21970. - f4[k] * 15048.) + (f2[k] * 22528. - f5[k] * 27360.), fabs(d__1));
2161  /* L250: */
2162  /* Computing MAX */
2163  d__1 = eeoet, d__2 = ee / et;
2164  eeoet = std::max(d__1, d__2);
2165  }
2166 
2167  esttol = fabs(*h__) * eeoet * scale / 752400.;
2168 
2169  if (esttol <= 1.)
2170  {
2171  goto L260;
2172  }
2173 
2174  /* unsuccessful step */
2175  /* reduce the stepsize , try again */
2176  /* the decrease is limited to a factor of 1/10 */
2177 
2178  hfaild = 1;
2179  output = 0;
2180  s = .1;
2181 
2182  if (esttol < 59049.)
2183  {
2184  s = .9 / pow(esttol, .2);
2185  }
2186 
2187  *h__ = s * *h__;
2188 
2189  if (fabs(*h__) > hmin)
2190  {
2191  goto L200;
2192  }
2193 
2194  /* requested error unattainable at smallest allowable stepsize */
2195  *iflag = 6;
2196  *kflag = 6;
2197  return 0;
2198 
2199  /* successful step */
2200  /* store solution at t+h */
2201  /* and evaluate derivatives there */
2202 
2203 L260:
2204  *t += *h__;
2205  i__1 = *neqn;
2206 
2207  for (k = 1; k <= i__1; ++k)
2208  {
2209  /* L270: */
2210  y[k] = f1[k];
2211  }
2212 
2213  a = *t;
2214  (*f)(neqn, &a, &y[1], &yp[1]);
2215  ++(*nfe);
2216 
2217  /* choose next stepsize */
2218  /* the increase is limited to a factor of 5 */
2219  /* if step failure has just occurred, next */
2220  /* stepsize is not allowed to increase */
2221 
2222  s = 5.;
2223 
2224  if (esttol > 1.889568e-4)
2225  {
2226  s = .9 / pow(esttol, .2);
2227  }
2228 
2229  if (hfaild)
2230  {
2231  s = std::min(s, 1.);
2232  }
2233 
2234  /* Computing MAX */
2235  d__2 = s * fabs(*h__);
2236  d__1 = std::max(d__2, hmin);
2237  //*h__ = d_sign(&d__1, h__);
2238  *h__ = *h__ >= 0 ? fabs(d__1) : -fabs(d__1);
2239 
2240  /* end of core integrator */
2241 
2242  /* should we take another step */
2243  if (output)
2244  {
2245  goto L300;
2246  }
2247 
2248  if (*iflag > 0)
2249  {
2250  goto L100;
2251  }
2252 
2253  /* integration successfully completed */
2254 
2255  /* one-step mode */
2256  *iflag = -2;
2257  return 0;
2258 
2259  /* interval mode */
2260 L300:
2261  *t = *tout;
2262  *iflag = 2;
2263  return 0;
2264 } /* rkfs_ */
#define C_INT
Definition: copasi.h:115
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)
void(* pEvalF)(const C_INT *, const double *, const double *, double *)
#define abs(x)
Definition: f2c.h:173
#define min(a, b)
Definition: f2c.h:175
#define max(a, b)
Definition: f2c.h:176
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().

485 {
486  size_t i, j;
487  CHybridODE45Balance newElement;
488  size_t maxBalance = 0;
489 
490  mLocalBalances.clear();
492  mLocalSubstrates.clear();
494 
495  for (i = 0; i < mNumReactions; i++)
496  {
497  const CCopasiVector <CChemEqElement> * balances =
498  &(*mpReactions)[i]->getChemEq().getBalances();
499 
500  for (j = 0; j < balances->size(); j++)
501  {
502  newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
503  newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
504  // + 0.5 to get a rounding out of the static_cast to C_INT32!
505  newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity()
506  + 0.5));
507 
508  if ((newElement.mpMetabolite->getStatus()) != CModelEntity::FIXED)
509  {
510  if (newElement.mMultiplicity > maxBalance) maxBalance = newElement.mMultiplicity;
511 
512  mLocalBalances[i].push_back(newElement); // element is copied for the push_back
513  }
514  }
515 
516  balances = &(*mpReactions)[i]->getChemEq().getSubstrates();
517 
518  for (j = 0; j < balances->size(); j++)
519  {
520  newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
521  newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
522  // + 0.5 to get a rounding out of the static_cast to C_INT32!
523  newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity()
524  + 0.5));
525 
526  mLocalSubstrates[i].push_back(newElement); // element is copied for the push_back
527  }
528  }
529 
530  mMaxBalance = maxBalance;
531 
532  return;
533 }
virtual size_t size() const
#define C_INT32
Definition: copasi.h:90
Definition: CMetab.h:178
std::vector< std::vector< CHybridODE45Balance > > mLocalSubstrates
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
virtual size_t getIndex(const CCopasiObject *pObject) const
const CModelEntity::Status & getStatus() const
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
void CHybridMethodODE45::setupCalculateSet ( )
protected

Definition at line 632 of file CHybridMethodODE45.cpp.

References FAST, mCalculateSet, mMetab2React, and mMetabFlags.

Referenced by initMethod().

633 {
634  mCalculateSet.clear();
635 
636  std::set<size_t>::iterator reactIt;
637 
638  std::vector<CHybridODE45MetabFlag>::iterator metabIt
639  = mMetabFlags.begin();
640  const std::vector<CHybridODE45MetabFlag>::iterator metabEndIt
641  = mMetabFlags.end();
642 
643  std::vector< std::set<size_t> >::iterator m2rIt =
644  mMetab2React.begin();
645 
646  for (; metabIt != metabEndIt; ++metabIt, ++m2rIt)
647  {
648  //for fast metabs, insert reaction into mCalculateSet
649  if (metabIt->mFlag == FAST) //fast metab
650  {
651  for (reactIt = m2rIt->begin();
652  reactIt != m2rIt->end();
653  reactIt++) // check all reactions that involves metab
654  mCalculateSet.insert(*reactIt);
655  }
656  }
657 
658  return;
659 }
#define FAST
std::vector< std::set< size_t > > mMetab2React
std::vector< CHybridODE45MetabFlag > mMetabFlags
std::set< size_t > mCalculateSet
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().

1287 {
1288  mDG.clear();
1289  std::vector< std::set<std::string>* > DependsOn;
1290  std::vector< std::set<std::string>* > Affects;
1291  size_t numReactions = mpReactions->size();
1292  size_t i, j;
1293 
1294  // Do for each reaction:
1295  for (i = 0; i < mNumReactions; i++)
1296  {
1297  // Get the set of metabolites which affect the value of amu for this
1298  // reaction i.e. the set on which amu depends. This may be more than
1299  // the set of substrates, since the kinetics can involve other
1300  // reactants, e.g. catalysts. We thus need to step through the
1301  // rate function and pick out every reactant which can vary.
1302  DependsOn.push_back(getDependsOn(i));
1303  // Get the set of metabolites which are affected when this reaction takes place
1304  Affects.push_back(getAffects(i));
1305  }
1306 
1307  // For each possible pair of reactions i and j, if the intersection of
1308  // Affects(i) with DependsOn(j) is non-empty, add a dependency edge from i to j.
1309  for (i = 0; i < mNumReactions; i++)
1310  {
1311  for (j = 0; j < mNumReactions; j++)
1312  {
1313  // Determine whether the intersection of these two sets is non-empty
1314  // Could also do this with set_intersection generic algorithm, but that
1315  // would require operator<() to be defined on the set elements.
1316 
1317  std::set<std::string>::iterator iter = Affects[i]->begin();
1318 
1319  for (; iter != Affects[i]->end(); iter++)
1320  {
1321  if (DependsOn[j]->count(*iter))
1322  {
1323  // The set intersection is non-empty
1324  mDG.addDependent(i, j);
1325  break;
1326  }
1327  }
1328  }
1329 
1330  // Ensure that self edges are included
1331  //mDG.addDependent(i, i);
1332  }
1333 
1334  // Delete the memory allocated in getDependsOn() and getAffects()
1335  // since this is allocated in other functions.
1336  for (i = 0; i < numReactions; i++)
1337  {
1338  delete DependsOn[i];
1339  delete Affects[i];
1340  }
1341 
1342  return;
1343 }
const CCopasiVectorNS< CReaction > * mpReactions
virtual size_t size() const
void addDependent(const size_t &node, const size_t &dependent)
CDependencyGraph mDG
std::set< std::string > * getDependsOn(size_t rIndex)
std::set< std::string > * getAffects(size_t rIndex)
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().

564 {
565  mMetab2React.clear();
567 
568  CMetab * pMetab;
569  CReaction * pReaction;
570 
571  //Check wether metabs play as substrates and modifiers
572  for (size_t rct = 0; rct < mNumReactions; ++rct)
573  {
574  size_t index;
575 
576  //deal with substrates
578  (*mpReactions)[rct]->getChemEq().getSubstrates();
579 
580  for (size_t i = 0; i < metab.size(); i++)
581  {
582  pMetab = const_cast < CMetab* >
583  (metab[i]->getMetabolite());
584  index = mpModel->getMetabolitesX().getIndex(pMetab);
585 
586  if ((pMetab->getStatus()) != CModelEntity::FIXED)
587  mMetab2React[index].insert(rct);
588  }
589 
590  //deal with modifiers
591  metab = (*mpReactions)[rct]->getChemEq().getModifiers();
592 
593  for (size_t i = 0; i < metab.size(); i++)
594  {
595  pMetab = const_cast < CMetab* >
596  (metab[i]->getMetabolite());
597  index = mpModel->getMetabolitesX().getIndex(pMetab);
598 
599  if ((pMetab->getStatus()) != CModelEntity::FIXED)
600  mMetab2React[index].insert(rct);
601  }
602  }
603 
604  //Check wether metabs appear in reaction laws
605  std::set< const CCopasiObject * > changed;
606 
607  for (size_t metab = 0; metab < mNumVariableMetabs; metab++)
608  {
609  pMetab = (*mpMetabolites)[metab];
610 
611  if (pMetab->getStatus() != CModelEntity::FIXED)
612  {
613  changed.clear();
614  changed.insert(mpModel->getValueReference());
615  changed.insert(pMetab->getValueReference());
616 
617  for (size_t react = 0; react < mNumReactions; react++)
618  {
619  pReaction = (*mpReactions)[react];
620 
621  if (pReaction->getParticleFluxReference()->dependsOn(changed))
622  mMetab2React[metab].insert(react);
623  }
624  }
625  }
626 
627  return;
628 }
CCopasiObject * getParticleFluxReference()
Definition: CReaction.cpp:207
virtual size_t size() const
std::vector< std::set< size_t > > mMetab2React
Definition: CMetab.h:178
bool dependsOn(DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const
virtual size_t getIndex(const CCopasiObject *pObject) const
const CModelEntity::Status & getStatus() const
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
CCopasiObject * getValueReference() const
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().

403 {
404  //size_t rctIndex;
406 
407  std::vector<CHybridODE45MetabFlag>::iterator metabIt
408  = mMetabFlags.begin();
409  const std::vector<CHybridODE45MetabFlag>::iterator metabEndIt
410  = mMetabFlags.end();
411 
412  // initialization
413  for (; metabIt != metabEndIt; ++metabIt)
414  {
415  metabIt->mFlag = SLOW;
416  metabIt->mFastReactions.clear();
417  }
418 
419  std::vector<CHybridODE45Balance>::iterator itBalance;
420  std::vector<CHybridODE45Balance>::iterator itEndBalance;
421 
422  // go over all mLocalBalances, if balance != 0,
423  // and reaction is FAST, insert into set
424  for (size_t rct = 0; rct < mNumReactions; ++rct)
425  {
426  if (mReactionFlags[rct] == FAST)
427  {
428  itBalance = mLocalBalances[rct].begin();
429  itEndBalance = mLocalBalances[rct].end();
430 
431  for (; itBalance != itEndBalance; ++itBalance)
432  {
433  size_t metab = itBalance->mIndex;
434  mMetabFlags[metab].mFastReactions.insert(rct);
435  mMetabFlags[metab].mFlag = FAST;
436  }
437  }
438  }
439 
440  return;
441 }
#define FAST
std::vector< CHybridODE45MetabFlag > mMetabFlags
#define SLOW
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
CVector< size_t > mReactionFlags
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().

467 {
472  else
473  mMethod = HYBRID;
474 
475  return;
476 }
#define HYBRID
#define STOCHASTIC
#define DETERMINISTIC
void CHybridMethodODE45::setupPriorityQueue ( C_FLOAT64  startTime = 0.0)
protected

Sets up the priority queue.

Parameters
startTimeThe time at which the simulation starts.

Sets up the priority queue, for pure SSA.

Parameters
startTimeThe 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().

541 {
542  size_t i;
543  C_FLOAT64 time;
544 
545  mPQ.clear();
547 
548  for (i = 0; i < mNumReactions; i++)
549  {
550  calculateAmu(i);
551  time = startTime + generateReactionTime(i);
552  mPQ.insertStochReaction(i, time);
553  }
554 
555  return;
556 }
C_FLOAT64 generateReactionTime(size_t rIndex)
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
CIndexedPriorityQueue mPQ
#define C_FLOAT64
Definition: copasi.h:92
void initializeIndexPointer(const size_t numberOfReactions)
void calculateAmu(size_t rIndex)
void CHybridMethodODE45::setupReactAffect ( )
protected

Definition at line 661 of file CHybridMethodODE45.cpp.

References mLocalBalances, mMetab2React, mNumReactions, and mReactAffect.

Referenced by initMethod().

662 {
663  mReactAffect.clear();
664  mReactAffect.resize(mNumReactions);
665 
666  std::vector<std::set<size_t> >::iterator rctIt =
667  mReactAffect.begin();
668  const std::vector<std::set<size_t> >::iterator rctEndIt =
669  mReactAffect.end();
670 
671  for (size_t rct = 0; rctIt != rctEndIt; ++rct, ++rctIt)
672  {
673  rctIt->clear();
674 
675  std::vector<CHybridODE45Balance>::iterator balIt =
676  mLocalBalances[rct].begin();
677  const std::vector<CHybridODE45Balance>::iterator balEndIt =
678  mLocalBalances[rct].end();
679 
680  for (; balIt != balEndIt; ++balIt)
681  {
682  size_t metabId = balIt->mIndex;
683 
684  std::set<size_t>::iterator it =
685  mMetab2React[metabId].begin();
686  const std::set<size_t>::iterator endIt =
687  mMetab2React[metabId].end();
688 
689  for (; it != endIt; ++it)
690  rctIt->insert(*it);
691  }
692  }
693 
694  return;
695 }
std::vector< std::set< size_t > > mMetab2React
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
std::vector< std::set< size_t > > mReactAffect
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().

444 {
445  mHasStoiReaction = false;
446  mHasDetermReaction = false;
448 
449  for (size_t rct = 0; rct < mNumReactions; rct++)
450  {
451  if ((*mpReactions)[rct]->isFast())
452  {
453  mReactionFlags[rct] = FAST;
454  mHasDetermReaction = true;
455  }
456  else
457  {
458  mReactionFlags[rct] = SLOW;
459  mHasStoiReaction = true;
460  }
461  }
462 
463  return;
464 }
const CCopasiVectorNS< CReaction > * mpReactions
#define FAST
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define SLOW
CVector< size_t > mReactionFlags
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().

260 {
261  //set inital state
262  mpProblem->getModel()->setState(*initialState);
263 
264  //set mpState
266 
267  //set the mpModel
269  assert(mpModel);
270 
271  //set mDoCorrection
272  if (mpModel->getModelType() == CModel::deterministic)
273  mDoCorrection = true;
274  else
275  mDoCorrection = false;
276 
277  //set mHasAssignments
279 
280  //set mFirstMetabIndex
281  mFirstMetabIndex = mpModel->getStateTemplate().getIndex(mpModel->getMetabolitesX()[0]);
282 
283  // Call initMethod function
285 
286  //output data to check init status
287  //outputData();
288  return;
289 }
Definition: CState.h:305
CTrajectoryProblem * mpProblem
void setState(const CState &state)
Definition: CModel.cpp:1785
const C_FLOAT64 & getTime() const
Definition: CModel.cpp:1190
static bool modelHasAssignments(const CModel *pModel)
const CState & getState() const
Definition: CModel.cpp:1771
CModel * getModel() const
void initMethod(C_FLOAT64 time)
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.

699 {
700  // write the current state to the model???
701 
702  // check for possible overflows
703  size_t i;
704  //size_t imax;
705 
706  // :TODO: Bug 774: This assumes that the number of variable metabs is the number
707  // of metabs determined by reaction. In addition they are expected at the beginning of the
708  // MetabolitesX which is not the case if we have metabolites of type ODE.
709  //for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++)
710  // if (mpProblem->getModel()->getMetabolitesX()[i]->getValue() >= mMaxIntBeforeStep)
711  // {
712  // throw exception or something like that
713  //}
714 
715  // do several steps
716  C_FLOAT64 time = mpModel->getTime();
717  C_FLOAT64 endTime = time + deltaT;
718 
719  for (i = 0; ((i < mMaxSteps) && (time < endTime)); i++)
720  time = doSingleStep(time, endTime);
721 
722  // Warning Message
723  if ((i >= mMaxSteps) && (!mMaxStepsReached))
724  {
725  mMaxStepsReached = true; //only report this message once
726  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.");
727  }
728 
729  // get back the particle numbers
730 
731  /* Set the variable metabolites ????*/
732  //C_FLOAT64 * Dbl = mpCurrentState->beginIndependent() + mFirstMetabIndex - 1;
733 
734  //for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, Dbl++)
735  // *Dbl = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
736 
737  //the task expects the result in mpCurrentState
738  //*mpCurrentState = mpProblem->getModel()->getState();
740  //mpCurrentState->setTime(time);
741 
742  return NORMAL;
743 }
const C_FLOAT64 & getTime() const
Definition: CModel.cpp:1190
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)
void CHybridMethodODE45::updatePriorityQueue ( size_t  rIndex,
C_FLOAT64  time 
)
protected

Updates the priority queue.

Parameters
rIndexA size_t giving the index of the fired reaction
timeA C_FLOAT64 holding the time taken by this reaction

Updates the priority queue.

Parameters
rIndexA size_t giving the index of the fired reaction (-1, if no stochastic reaction has fired)
timeA 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().

1354 {
1355  C_FLOAT64 newTime;
1356  size_t index;
1357  std::set <size_t>::iterator iter, iterEnd;
1358 
1359  //if the model contains assignments we use a less efficient loop over all (stochastic)
1360  //reactions to capture all changes
1361  //we do not know the exact dependencies.
1362  //TODO: this should be changed later in order to get a more efficient update scheme
1363  if (mHasAssignments)
1364  {
1366 
1367  for (index = 0; index < mNumReactions; index++)
1368  {
1369  mAmuOld[index] = mAmu[index];
1370  calculateAmu(index);
1371 
1372  if (mAmuOld[index] != mAmu[index])
1373  if (index != rIndex) updateTauMu(index, time);
1374  }
1375  }
1376  else
1377  {
1378  // iterate through the set of affected reactions and update the stochastic
1379  // ones in the priority queue
1380  iter = mReactAffect[rIndex].begin();
1381  iterEnd = mReactAffect[rIndex].end();
1382 
1383  for (; iter != iterEnd; iter++)
1384  {
1385  index = *iter;
1386  mAmuOld[index] = mAmu[index];
1387  calculateAmu(index);
1388 
1389  if (*iter != rIndex) updateTauMu(index, time);
1390  }
1391  }
1392 
1393  // draw new random number and update the reaction just fired
1394  if (rIndex != C_INVALID_INDEX)
1395  {
1396  // reaction is stochastic
1397  newTime = time + generateReactionTime(rIndex);
1398  mPQ.updateNode(rIndex, newTime);
1399  }
1400 
1401  return;
1402 }
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
CVector< C_FLOAT64 > mAmuOld
C_FLOAT64 generateReactionTime(size_t rIndex)
#define C_INVALID_INDEX
Definition: copasi.h:222
CVector< C_FLOAT64 > mAmu
CIndexedPriorityQueue mPQ
void updateNode(const size_t index, const C_FLOAT64 key)
void updateTauMu(size_t rIndex, C_FLOAT64 time)
#define C_FLOAT64
Definition: copasi.h:92
std::vector< std::set< size_t > > mReactAffect
void calculateAmu(size_t rIndex)
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
rIndexA size_t specifying the index of the reaction
timeA 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
rIndexA 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().

1123 {
1124  C_FLOAT64 newTime;
1125 
1126  // One must make sure that the calculation yields reasonable results even in the cases
1127  // where mAmu=0 or mAmuOld=0 or both=0. Therefore mAmuOld=0 is checked. If mAmuOld equals 0,
1128  // then a new random number has to be drawn, because tau equals inf and the stoch.
1129  // information is lost.
1130  // If both values equal 0, then tau should remain inf and the update of the queue can
1131  // be skipped!
1132  if (mAmuOld[rIndex] == 0.0)
1133  {
1134  if (mAmu[rIndex] != 0.0)
1135  {
1136  newTime = time + generateReactionTime(rIndex);
1137  mPQ.updateNode(rIndex, newTime);
1138  }
1139  }
1140  else
1141  {
1142  newTime = time + (mAmuOld[rIndex] / mAmu[rIndex]) * (mPQ.getKey(rIndex) - time);
1143  mPQ.updateNode(rIndex, newTime);
1144  }
1145 
1146  return;
1147 }
CVector< C_FLOAT64 > mAmuOld
C_FLOAT64 generateReactionTime(size_t rIndex)
C_FLOAT64 getKey(const size_t index) const
CVector< C_FLOAT64 > mAmu
CIndexedPriorityQueue mPQ
void updateNode(const size_t index, const C_FLOAT64 key)
#define C_FLOAT64
Definition: copasi.h:92

Friends And Related Function Documentation

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