COPASI API  4.16.103
Public Member Functions | Static Public Member Functions | Protected Types | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Private Member Functions | Friends | List of all members
CHybridMethod Class Referenceabstract

#include <CHybridMethod.h>

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

Public Member Functions

 CHybridMethod (const CHybridMethod &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)
 
 ~CHybridMethod ()
 
- 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 ()
 

Static Public Member Functions

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

Protected Types

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

Protected Member Functions

void calculateAmu (size_t rIndex)
 
void calculateDerivative (std::vector< C_FLOAT64 > &deriv)
 
 CHybridMethod (const CCopasiContainer *pParent=NULL)
 
void cleanup ()
 
virtual C_FLOAT64 doSingleStep (C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
 
void fireReaction (size_t rIndex)
 
C_FLOAT64 generateReactionTime (size_t rIndex)
 
std::set< std::string > * getAffects (size_t rIndex)
 
std::set< std::string > * getDependsOn (size_t rIndex)
 
std::set< size_t > * getParticipatesIn (size_t rIndex)
 
void getState (std::vector< C_FLOAT64 > &target)
 
void getStochTimeAndIndex (C_FLOAT64 &ds, size_t &rIndex)
 
void initMethod (C_FLOAT64 time)
 
void insertDeterministicReaction (size_t rIndex)
 
void integrateDeterministicPart (C_FLOAT64 ds)
 
void integrateDeterministicPartEuler (C_FLOAT64 ds)
 
void outputData (std::ostream &os, C_INT32 mode)
 
void outputDebug (std::ostream &os, size_t level)
 
void partitionSystem ()
 
void removeDeterministicReaction (size_t rIndex)
 
void rungeKutta (C_FLOAT64 dt)
 
void setState (std::vector< C_FLOAT64 > &source)
 
void setupBalances ()
 
void setupDependencyGraph ()
 
void setupMetab2React ()
 
void setupMetab2ReactComplete ()
 
void setupMetab2ReactPlusModifier ()
 
void setupPartition ()
 
void setupPriorityQueue (C_FLOAT64 startTime=0.0)
 
void updatePriorityQueue (size_t rIndex, C_FLOAT64 time)
 
void updateTauMu (size_t rIndex, C_FLOAT64 time)
 
- Protected Member Functions inherited from CTrajectoryMethod
 CTrajectoryMethod (const CCopasiMethod::SubType &subType, const CCopasiContainer *pParent=NULL)
 
- Protected Member Functions inherited from CCopasiMethod
 CCopasiMethod (const CCopasiTask::Type &taskType, const SubType &subType, const CCopasiContainer *pParent=NULL)
 
- Protected Member Functions inherited from CCopasiParameterGroup
 CCopasiParameterGroup ()
 
- Protected Member Functions inherited from CCopasiContainer
template<class CType >
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 C_INT32 checkModel (CModel *model)
 
static bool modelHasAssignments (const CModel *pModel)
 

Protected Attributes

std::vector< C_FLOAT64currentState
 
std::vector< C_FLOAT64k1
 
std::vector< C_FLOAT64k2
 
std::vector< C_FLOAT64k3
 
std::vector< C_FLOAT64k4
 
std::vector< C_FLOAT64mAmu
 
std::vector< C_FLOAT64mAmuOld
 
CDependencyGraph mDG
 
bool mDoCorrection
 
size_t mFirstMetabIndex
 
CHybridStochFlagmFirstReactionFlag
 
bool mHasAssignments
 
std::vector< std::vector
< CHybridBalance > > 
mLocalBalances
 
std::vector< std::vector
< CHybridBalance > > 
mLocalSubstrates
 
C_FLOAT64 mLowerStochLimit
 
size_t mMaxBalance
 
size_t mMaxIntBeforeStep
 
unsigned C_INT32 mMaxSteps
 
bool mMaxStepsReached
 
std::vector< std::set< size_t > > mMetab2React
 
std::vector< metabStatusmMetabFlags
 
size_t mNumVariableMetabs
 
size_t mOutputCounter
 
std::ofstream mOutputFile
 
std::string mOutputFileName
 
unsigned C_INT32 mPartitioningInterval
 
CCopasiVector< CMetab > * mpMetabolites
 
CModelmpModel
 
CIndexedPriorityQueue mPQ
 
CRandommpRandomGenerator
 
const CCopasiVectorNS
< CReaction > * 
mpReactions
 
unsigned C_INT32 mRandomSeed
 
std::vector< CHybridStochFlagmReactionFlags
 
size_t mStepsAfterPartitionSystem
 
C_FLOAT64 mStepsize
 
CMatrix< C_FLOAT64mStoi
 
std::set< size_t > mUpdateSet
 
C_FLOAT64 mUpperStochLimit
 
bool mUseRandomSeed
 
std::vector< 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

void initializeParameter ()
 

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 Attributes inherited from CCopasiMethod
static const std::string SubTypeName []
 
static const char * XMLSubType []
 
- Static Public Attributes inherited from CCopasiParameter
static const std::string TypeName []
 
static const char * XMLType []
 
- Static Public Attributes inherited from CCopasiContainer
static const std::vector
< CCopasiContainer * > 
EmptyList
 
- Static Protected Attributes inherited from CCopasiObject
static CRenameHandlersmpRenameHandler = NULL
 

Detailed Description

Definition at line 113 of file CHybridMethod.h.

Member Enumeration Documentation

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

Enumerator
LOW 
HIGH 

Definition at line 462 of file CHybridMethod.h.

Constructor & Destructor Documentation

CHybridMethod::CHybridMethod ( const CHybridMethod src,
const CCopasiContainer pParent = NULL 
)

Copy constructor

Parameters
constCHybridMethod & src
constCCopasiContainer * pParent (default: NULL)

Definition at line 70 of file CHybridMethod.cpp.

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

71  :
72  CTrajectoryMethod(src, pParent)
73 {
76 }
CRandom * mpRandomGenerator
void initializeParameter()
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
CHybridMethod::~CHybridMethod ( )

Destructor.

Definition at line 81 of file CHybridMethod.cpp.

References cleanup(), and DESTRUCTOR_TRACE.

82 {
83  cleanup();
85 }
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
CHybridMethod::CHybridMethod ( const CCopasiContainer pParent = NULL)
protected

Default constructor.

Parameters
constCCopasiContainer * pParent (default: NULL)

CHybridMethod

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

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

Last change: 14, December 2004

(C) European Media Lab 2003. Default constructor.

Definition at line 63 of file CHybridMethod.cpp.

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

63  :
65 {
68 }
CRandom * mpRandomGenerator
void initializeParameter()
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49

Member Function Documentation

void CHybridMethod::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 691 of file CHybridMethod.cpp.

References C_FLOAT64, C_INT32, C_INT64, CCopasiParameter::getValue(), mAmu, mDoCorrection, mLocalSubstrates, and mpMetabolites.

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

692 {
693  if (!mDoCorrection)
694  {
695  mAmu[rIndex] = (*mpReactions)[rIndex]->calculateParticleFlux();
696  return;
697  }
698 
699  // We need the product of the cmu and hmu for this step.
700  // We calculate this in one go, as there are fewer steps to
701  // perform and we eliminate some possible rounding errors.
702  C_FLOAT64 amu = 1; // initially
703  //size_t total_substrates = 0;
704  C_INT32 num_ident = 0;
705  C_INT64 number = 0;
706  C_INT64 lower_bound;
707  // substrate_factor - The substrates, raised to their multiplicities,
708  // multiplied with one another. If there are, e.g. m substrates of type m,
709  // and n of type N, then substrate_factor = M^m * N^n.
710  C_FLOAT64 substrate_factor = 1;
711  // First, find the reaction associated with this index.
712  // Keep a pointer to this.
713  // Iterate through each substrate in the reaction
714  const std::vector<CHybridBalance> & substrates = mLocalSubstrates[rIndex];
715 
716  int flag = 0;
717 
718  for (size_t i = 0; i < substrates.size(); i++)
719  {
720  num_ident = substrates[i].mMultiplicity;
721 
722  if (num_ident > 1)
723  {
724  flag = 1;
725  number = static_cast<C_INT64>(floor((*mpMetabolites)[substrates[i].mIndex]->getValue()));
726  lower_bound = number - num_ident;
727  substrate_factor = substrate_factor * pow((double) number, (int) num_ident - 1); //optimization
728 
729  number--; // optimization
730 
731  while (number > lower_bound)
732  {
733  amu *= number;
734  number--;
735  }
736  }
737  }
738 
739  if ((amu == 0) || (substrate_factor == 0)) // at least one substrate particle number is zero
740  {
741  mAmu[rIndex] = 0;
742  return;
743  }
744 
745  // rate_factor is the rate function divided by substrate_factor.
746  // It would be more efficient if this was generated directly, since in effect we
747  // are multiplying and then dividing by the same thing (substrate_factor)!
748  C_FLOAT64 rate_factor = (*mpReactions)[rIndex]->calculateParticleFlux();
749 
750  if (flag)
751  {
752  amu *= rate_factor / substrate_factor;;
753  mAmu[rIndex] = amu;
754  }
755  else
756  {
757  mAmu[rIndex] = rate_factor;
758  }
759 
760  return;
761 
762  // a more efficient way to calculate mass action kinetics could be included
763 }
#define C_INT64
Definition: copasi.h:88
#define C_INT32
Definition: copasi.h:90
const Value & getValue() const
CCopasiVector< CMetab > * mpMetabolites
std::vector< C_FLOAT64 > mAmu
long int flag
Definition: f2c.h:52
#define C_FLOAT64
Definition: copasi.h:92
std::vector< std::vector< CHybridBalance > > mLocalSubstrates
void CHybridMethod::calculateDerivative ( std::vector< C_FLOAT64 > &  deriv)
protected

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

Parameters
derivA vector reference of length mNumVariableMetabs, into which the derivative is written

Definition at line 489 of file CHybridMethod.cpp.

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

Referenced by integrateDeterministicPartEuler(), and rungeKutta().

490 {
491  size_t i;
492  C_INT32 bal = 0;
493  CHybridStochFlag * j;
494 
495  // Calculate all the needed kinetic functions of the deterministic reactions
496  for (j = mFirstReactionFlag; j != NULL; j = j->mpNext)
497  {
498  (*mpReactions)[j->mIndex]->calculateParticleFlux();
499  }
500 
501  // For each metabolite add up the contributions of the det. reactions
502  // the number of rows in mStoi equals the number of non-fixed metabolites!
503  // for (i=0; i<mNumVariableMetabs; i++)
504  for (i = 0; i < mStoi.numRows(); i++)
505  {
506  deriv[i] = 0.0;
507 
508  for (j = mFirstReactionFlag; j != NULL; j = j->mpNext)
509  {
510  // juergen: +0.5 to get a rounding out of the static_cast
511  bal = static_cast< C_INT32 >(floor(mStoi[i][j->mIndex] + 0.5));
512  deriv[i] += bal * (*mpReactions)[j->mIndex]->getParticleFlux(); // balance * flux;
513  }
514  }
515 
516  /*
517  for (; i < mNumVariableMetabs; i++) deriv[i] = 0.0; // important to get a correct deriv vector, because mStoi doesn't cover fixed metabolites
518  */
519  return;
520 }
virtual size_t numRows() const
Definition: CMatrix.h:138
CHybridStochFlag * mFirstReactionFlag
CMatrix< C_FLOAT64 > mStoi
#define C_INT32
Definition: copasi.h:90
CHybridStochFlag * mpNext
Definition: CHybridMethod.h:92
C_INT32 CHybridMethod::checkModel ( CModel model)
staticprotected

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

Parameters
modelThe model to be checked
Returns
1, if hybrid simulation is possible; <0, if an error occured.

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

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

Definition at line 804 of file CHybridMethod.cpp.

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

805 {
807  CMatrix <C_FLOAT64> mStoi = model->getStoi();
808  size_t i, numReactions = mpReactions->size();
809  size_t j;
810  C_INT32 multInt;
811 
812  C_FLOAT64 multFloat;
813  // size_t metabSize = mpMetabolites->size();
814 
815  for (i = 0; i < numReactions; i++) // for every reaction
816  {
817  // TEST getCompartmentNumber() == 1
818  if ((*mpReactions)[i]->getCompartmentNumber() != 1) return - 1;
819 
820  // TEST isReversible() == 0
821  if ((*mpReactions)[i]->isReversible() != 0) return - 2;
822 
823  // TEST integer stoichometry
824  // Iterate through each the metabolites
825  // juergen: the number of rows of mStoi equals the number of non-fixed metabs!
826  // for (j=0; i<metabSize; j++)
827  for (j = 0; j < mStoi.numRows(); j++)
828  {
829  multFloat = mStoi[j][i];
830  multInt = static_cast<C_INT32>(floor(multFloat + 0.5)); // +0.5 to get a rounding out of the static_cast to int!
831 
832  if ((multFloat - multInt) > INT_EPSILON) return - 3; // INT_EPSILON in CHybridMethod.h
833  }
834  }
835 
836  return 1; // Model is appropriate for hybrid simulation
837 }
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
const CCopasiVectorNS< CReaction > * mpReactions
#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 CHybridMethod::cleanup ( )
protected

Cleans up memory, etc.

Definition at line 319 of file CHybridMethod.cpp.

References mpModel, and mpRandomGenerator.

Referenced by ~CHybridMethod().

320 {
321  delete mpRandomGenerator;
322  mpRandomGenerator = NULL;
323  mpModel = NULL;
324  return;
325 }
CRandom * mpRandomGenerator
CModel * mpModel
CHybridMethod * CHybridMethod::createHybridMethod ( )
static

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

Definition at line 153 of file CHybridMethod.cpp.

References C_INT32.

Referenced by CTrajectoryMethod::createMethod().

154 {
155  C_INT32 result = 1; // hybrid NextReactionRungeKutta method as default
156  /* if (pProblem && pProblem->getModel())
157  {
158  result = checkModel(pProblem->getModel());
159  }*/
160  CHybridMethod * method = NULL;
161 
162  switch (result)
163  {
164  /* case - 3: // non-integer stoichometry
165  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 1);
166  break;
167  case - 2: // reversible reaction exists
168  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 2);
169  break;
170 
171  case - 1: // more than one compartment involved
172  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 3);
173  break;*/
174  case 1:
175  default:
176  // Everything alright: Hybrid simulation possible
177  method = new CHybridNextReactionRKMethod();
178  break;
179  }
180 
181  return method;
182 }
#define C_INT32
Definition: copasi.h:90
virtual C_FLOAT64 CHybridMethod::doSingleStep ( C_FLOAT64  currentTime,
C_FLOAT64  endTime 
)
protectedpure virtual

Simulates the system over the next interval of time. The current time and the end time of the current step() are given as arguments.

Parameters
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

Implemented in CHybridNextReactionRKMethod.

Referenced by step().

bool CHybridMethod::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 143 of file CHybridMethod.cpp.

References initializeParameter().

144 {
146  return true;
147 }
void initializeParameter()
void CHybridMethod::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.

Executes the specified reaction in the system once.

Parameters
rIndexA size_t specifying the index of the reaction, which will be fired.
timeThe current time

Definition at line 588 of file CHybridMethod.cpp.

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

Referenced by CHybridNextReactionRKMethod::doSingleStep().

589 {
590  // Change the particle numbers according to which step took place.
591  // First, get the vector of balances in the reaction we've got.
592  // (This vector expresses the number change of each metabolite
593  // in the reaction.) Then step through each balance, using its
594  // multiplicity to calculate a new value for the associated
595  // metabolite. Finally, update the metabolite.
596 
597  size_t i;
598  C_FLOAT64 newNumber;
599  CMetab * pMetab;
600 
601  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
602  {
603  pMetab = mLocalBalances[rIndex][i].mpMetabolite;
604  newNumber = pMetab->getValue() + mLocalBalances[rIndex][i].mMultiplicity;
605 
606  pMetab->setValue(newNumber);
607  pMetab->refreshConcentration();
608  }
609 
610  // insert all dependent reactions into the mUpdateSet
611  const std::set <size_t> & dependents = mDG.getDependents(rIndex);
612  std::copy(dependents.begin(), dependents.end(),
613  std::inserter(mUpdateSet, mUpdateSet.begin()));
614 
615  return;
616 }
const std::set< size_t > & getDependents(const size_t &node) const
virtual void setValue(const C_FLOAT64 &value)
Definition: CMetab.h:178
CDependencyGraph mDG
void refreshConcentration()
Definition: CMetab.cpp:281
const C_FLOAT64 & getValue() const
#define C_FLOAT64
Definition: copasi.h:92
std::vector< std::vector< CHybridBalance > > mLocalBalances
std::set< size_t > mUpdateSet
C_FLOAT64 CHybridMethod::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 678 of file CHybridMethod.cpp.

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

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

679 {
680  if (mAmu[rIndex] == 0) return std::numeric_limits<C_FLOAT64>::infinity();
681 
683  return - 1 * log(rand2) / mAmu[rIndex];
684 }
CRandom * mpRandomGenerator
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
std::vector< C_FLOAT64 > mAmu
#define C_FLOAT64
Definition: copasi.h:92
std::set< std::string > * CHybridMethod::getAffects ( size_t  rIndex)
protected

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 1331 of file CHybridMethod.cpp.

References mLocalBalances.

Referenced by setupDependencyGraph().

1332 {
1333  size_t i;
1334  std::set<std::string> *retset = new std::set<std::string>;
1335 
1336  // Get the balances associated with the reaction at this index
1337  // XXX We first get the chemical equation, then the balances, since the getBalances method in CReaction is unimplemented!
1338 
1339  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
1340  {
1341  if (mLocalBalances[rIndex][i].mMultiplicity != 0)
1342  {
1343  retset->insert(mLocalBalances[rIndex][i].mpMetabolite->getKey());
1344  }
1345  }
1346 
1347  return retset;
1348 }
std::vector< std::vector< CHybridBalance > > mLocalBalances
std::set< std::string > * CHybridMethod::getDependsOn ( size_t  rIndex)
protected

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 1298 of file CHybridMethod.cpp.

References mpReactions, and CFunctionParameter::PARAMETER.

Referenced by setupDependencyGraph().

1299 {
1300  std::set<std::string> *retset = new std::set<std::string>;
1301 
1302  size_t i, imax = (*mpReactions)[rIndex]->getFunctionParameters().size();
1303  size_t j, jmax;
1304 
1305  for (i = 0; i < imax; ++i)
1306  {
1307  if ((*mpReactions)[rIndex]->getFunctionParameters()[i]->getUsage() == CFunctionParameter::PARAMETER)
1308  continue;
1309 
1310  //metablist = (*mpReactions)[rIndex]->getParameterMappingMetab(i);
1311  const std::vector <std::string> & metabKeylist =
1312  (*mpReactions)[rIndex]->getParameterMappings()[i];
1313  jmax = metabKeylist.size();
1314 
1315  for (j = 0; j < jmax; ++j)
1316  {
1317  retset->insert(metabKeylist[j]);
1318  }
1319  }
1320 
1321  return retset;
1322 }
const CCopasiVectorNS< CReaction > * mpReactions
std::set< size_t > * CHybridMethod::getParticipatesIn ( size_t  rIndex)
protected

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

Parameters
rIndexThe index of the reaction being executed.
Returns
The set of participating metabolites.

Definition at line 1357 of file CHybridMethod.cpp.

Referenced by setupMetab2ReactPlusModifier().

1358 {
1359  std::set<size_t> *retset = new std::set<size_t>;
1360  return retset;
1361 }
void CHybridMethod::getState ( std::vector< C_FLOAT64 > &  target)
protected

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

Parameters
targetA vector reference of length mNumVariableMetabs, into which the state of the system is written

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

Parameters
targetAn array of C_FLOAT64s with length mNumVariableMetabs, into which the state of the system is written

Definition at line 529 of file CHybridMethod.cpp.

References mNumVariableMetabs.

Referenced by integrateDeterministicPartEuler(), and rungeKutta().

530 {
531  size_t i;
532 
533  for (i = 0; i < mNumVariableMetabs; i++)
534  {
535  target[i] = (*mpMetabolites)[i]->getValue();
536  }
537 
538  return;
539 }
size_t mNumVariableMetabs
void CHybridMethod::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 574 of file CHybridMethod.cpp.

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

Referenced by CHybridNextReactionRKMethod::doSingleStep().

575 {
576  ds = mPQ.topKey();
577  rIndex = mPQ.topIndex();
578  return;
579 }
CIndexedPriorityQueue mPQ
void CHybridMethod::initializeParameter ( )
private

Intialize the method parameter

Definition at line 87 of file CHybridMethod.cpp.

References CCopasiParameterGroup::assertParameter(), CCopasiParameter::BOOL, C_FLOAT64, C_INT32, CCopasiParameter::DOUBLE, CCopasiParameterGroup::getParameter(), CCopasiParameter::getValue(), CCopasiParameter::INT, LOWER_STOCH_LIMIT, MAX_STEPS, PARTITIONING_INTERVAL, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pINT, CCopasiParameter::Value::pUINT, RANDOM_SEED, CCopasiParameterGroup::removeParameter(), RUNGE_KUTTA_STEPSIZE, CCopasiParameterGroup::setValue(), CCopasiParameter::UINT, UPPER_STOCH_LIMIT, and USE_RANDOM_SEED.

Referenced by CHybridMethod(), and elevateChildren().

88 {
89  CCopasiParameter *pParm;
90 
91  assertParameter("Max Internal Steps", CCopasiParameter::INT, (C_INT32) MAX_STEPS);
95  assertParameter("Partitioning Interval", CCopasiParameter::UINT, (unsigned C_INT32) PARTITIONING_INTERVAL);
96  assertParameter("Use Random Seed", CCopasiParameter::BOOL, (bool) USE_RANDOM_SEED);
97  assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) RANDOM_SEED);
98 
99  // Check whether we have a method with the old parameter names
100  if ((pParm = getParameter("HYBRID.MaxSteps")) != NULL)
101  {
102  setValue("Max Internal Steps", *pParm->getValue().pINT);
103  removeParameter("HYBRID.MaxSteps");
104 
105  if ((pParm = getParameter("HYBRID.LowerStochLimit")) != NULL)
106  {
107  setValue("Lower Limit", *pParm->getValue().pDOUBLE);
108  removeParameter("HYBRID.LowerStochLimit");
109  }
110 
111  if ((pParm = getParameter("HYBRID.UpperStochLimit")) != NULL)
112  {
113  setValue("Upper Limit", *pParm->getValue().pDOUBLE);
114  removeParameter("HYBRID.UpperStochLimit");
115  }
116 
117  if ((pParm = getParameter("HYBRID.RungeKuttaStepsize")) != NULL)
118  {
119  setValue("Runge Kutta Stepsize", *pParm->getValue().pDOUBLE);
120  removeParameter("HYBRID.RungeKuttaStepsize");
121  }
122 
123  if ((pParm = getParameter("HYBRID.PartitioningInterval")) != NULL)
124  {
125  setValue("Partitioning Interval", *pParm->getValue().pUINT);
126  removeParameter("HYBRID.PartitioningInterval");
127  }
128 
129  if ((pParm = getParameter("UseRandomSeed")) != NULL)
130  {
131  setValue("Use Random Seed", *pParm->getValue().pBOOL);
132  removeParameter("UseRandomSeed");
133  }
134 
135  if ((pParm = getParameter("")) != NULL)
136  {
137  setValue("Random Seed", *pParm->getValue().pUINT);
138  removeParameter("");
139  }
140  }
141 }
#define USE_RANDOM_SEED
Definition: CHybridMethod.h:67
#define UPPER_STOCH_LIMIT
Definition: CHybridMethod.h:61
#define PARTITIONING_INTERVAL
Definition: CHybridMethod.h:63
#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
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
#define RUNGE_KUTTA_STEPSIZE
Definition: CHybridMethod.h:62
#define LOWER_STOCH_LIMIT
Definition: CHybridMethod.h:60
void CHybridMethod::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 264 of file CHybridMethod.cpp.

References currentState, CModel::getMetabolitesX(), CModel::getNumDependentReactionMetabs(), CModel::getNumIndependentReactionMetabs(), CModel::getReactions(), CModel::getStoi(), CCopasiParameter::getValue(), CRandom::initialize(), k1, k2, k3, k4, mAmu, mAmuOld, mLowerStochLimit, mMaxSteps, mMaxStepsReached, mNumVariableMetabs, mPartitioningInterval, mpMetabolites, mpModel, mpRandomGenerator, mpReactions, mRandomSeed, mStepsAfterPartitionSystem, mStepsize, mStoi, mUpdateSet, mUpperStochLimit, mUseRandomSeed, CCopasiParameter::Value::pBOOL, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pINT, CCopasiParameter::Value::pUINT, setupBalances(), setupDependencyGraph(), setupMetab2React(), setupPartition(), setupPriorityQueue(), CCopasiVector< T >::size(), and temp.

Referenced by start().

265 {
267  mAmu.clear();
268  mAmu.resize(mpReactions->size());
269  mAmuOld.clear();
270  mAmuOld.resize(mpReactions->size());
272  //mNumVariableMetabs = mpModel->getNumVariableMetabs(); // ind + dep metabs, without fixed metabs
273  mNumVariableMetabs = mpModel->getNumIndependentReactionMetabs() + mpModel->getNumDependentReactionMetabs(); // ind + dep metabs, without fixed metabs
274  // mNumVariableMetabs = mpCurrentState->getNumVariable(); // mpBeginFixed - mpBeginIndependent
275 
276  temp.clear();
277  temp.resize(mNumVariableMetabs);
278  currentState.clear();
280 
281  k1.clear();
282  k1.resize(mNumVariableMetabs);
283  k2.clear();
284  k2.resize(mNumVariableMetabs);
285  k3.clear();
286  k3.resize(mNumVariableMetabs);
287  k4.clear();
288  k4.resize(mNumVariableMetabs);
289 
290  /* get configuration data */
291  mMaxSteps = * getValue("Max Internal Steps").pINT;
292  mLowerStochLimit = * getValue("Lower Limit").pDOUBLE;
293  mUpperStochLimit = * getValue("Upper Limit").pDOUBLE;
294  mStepsize = * getValue("Runge Kutta Stepsize").pDOUBLE;
295  mPartitioningInterval = * getValue("Partitioning Interval").pUINT;
296  mUseRandomSeed = * getValue("Use Random Seed").pBOOL;
297  mRandomSeed = * getValue("Random Seed").pUINT;
298 
300 
301  mStoi = mpModel->getStoi();
303  mUpdateSet.clear();
304 
305  setupBalances(); // initialize mLocalBalances and mLocalSubstrates (has to be called first!)
306  setupDependencyGraph(); // initialize mDG
307  setupMetab2React(); // initialize mMetab2React
308  setupPartition(); // initialize mReactionFlags
309  setupPriorityQueue(start_time); // initialize mPQ
310 
311  mMaxStepsReached = false;
312 
313  return;
314 }
CRandom * mpRandomGenerator
C_FLOAT64 mLowerStochLimit
unsigned C_INT32 mMaxSteps
CModel * mpModel
virtual size_t size() const
std::vector< C_FLOAT64 > k3
CMatrix< C_FLOAT64 > mStoi
const CCopasiVectorNS< CReaction > * mpReactions
std::vector< C_FLOAT64 > k4
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
const Value & getValue() const
unsigned C_INT32 * pUINT
void setupMetab2React()
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
CCopasiVector< CMetab > * mpMetabolites
std::vector< C_FLOAT64 > mAmu
size_t mStepsAfterPartitionSystem
C_FLOAT64 mStepsize
unsigned C_INT32 mPartitioningInterval
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
std::vector< C_FLOAT64 > currentState
std::vector< C_FLOAT64 > k2
unsigned C_INT32 mRandomSeed
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
C_FLOAT64 mUpperStochLimit
std::vector< C_FLOAT64 > k1
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
std::set< size_t > mUpdateSet
std::vector< C_FLOAT64 > mAmuOld
std::vector< C_FLOAT64 > temp
size_t mNumVariableMetabs
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
void setupDependencyGraph()
void CHybridMethod::insertDeterministicReaction ( size_t  rIndex)
protected

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

Parameters
rIndexA size_t giving the index of the reaction to be inserted into the list of deterministic reactions.

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

Parameters
rIndexA size_t giving the index of the reaction to be inserted into the list of deterministic reactions.

Definition at line 1224 of file CHybridMethod.cpp.

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

Referenced by partitionSystem().

1225 {
1226  if (mReactionFlags[rIndex].mpPrev == NULL)
1227  // reaction is stochastic (avoids double insertions)
1228  {
1229  if (mFirstReactionFlag != NULL)
1230  // there are deterministic reactions already
1231  {
1233  mReactionFlags[rIndex].mpNext = mFirstReactionFlag;
1236  }
1237  else
1238  {
1239  // there are no deterministic reactions
1240  // Important to distinguish between stochastic (prev == NULL) and deterministic (prev != NULL) reactions
1241  mReactionFlags[rIndex].mpPrev = &mReactionFlags[rIndex];
1243  }
1244 
1245  mAmu[rIndex] = 0.0;
1246  mAmuOld[rIndex] = 0.0;
1247  }
1248 
1249  return;
1250 }
CHybridStochFlag * mFirstReactionFlag
CHybridStochFlag * mpPrev
Definition: CHybridMethod.h:91
std::vector< C_FLOAT64 > mAmu
std::vector< C_FLOAT64 > mAmuOld
std::vector< CHybridStochFlag > mReactionFlags
void CHybridMethod::integrateDeterministicPart ( C_FLOAT64  dt)
protected

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

Parameters
dsA C_FLOAT64 specifying the stepsize.

Definition at line 335 of file CHybridMethod.cpp.

References C_FLOAT64, CDependencyGraph::getDependents(), mDG, mFirstReactionFlag, CHybridStochFlag::mIndex, CHybridStochFlag::mpNext, mStepsize, mUpdateSet, and rungeKutta().

Referenced by CHybridNextReactionRKMethod::doSingleStep().

336 {
337  C_FLOAT64 integrationTime = 0.0;
338  CHybridStochFlag * react = NULL;
339 
340  // This method uses a 4th order RungeKutta-method to integrate the deterministic part of the system. Maybe a better numerical method (adaptive stepsize, lsoda, ...) should be introduced here later on
341 
342  while ((dt - integrationTime) > mStepsize)
343  {
344  rungeKutta(mStepsize); // for the deterministic part of the system
345  integrationTime += mStepsize;
346  }
347 
348  rungeKutta(dt - integrationTime);
349 
350  // find the set union of all reactions, which depend on one of the deterministic reactions. The propensities of the stochastic reactions in this set union will be updated later in the method updatePriorityQueue().
351  for (react = mFirstReactionFlag; react != NULL; react = react->mpNext)
352  {
353  const std::set <size_t> & dependents = mDG.getDependents(react->mIndex);
354  std::copy(dependents.begin(), dependents.end(),
355  std::inserter(mUpdateSet, mUpdateSet.begin()));
356  }
357 
358  return;
359 }
CHybridStochFlag * mFirstReactionFlag
const std::set< size_t > & getDependents(const size_t &node) const
void rungeKutta(C_FLOAT64 dt)
CDependencyGraph mDG
CHybridStochFlag * mpNext
Definition: CHybridMethod.h:92
C_FLOAT64 mStepsize
#define C_FLOAT64
Definition: copasi.h:92
std::set< size_t > mUpdateSet
void CHybridMethod::integrateDeterministicPartEuler ( C_FLOAT64  dt)
protected

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

Parameters
dsA C_FLOAT64 specifying the stepsize.

Definition at line 367 of file CHybridMethod.cpp.

References C_FLOAT64, calculateDerivative(), currentState, CDependencyGraph::getDependents(), getState(), mDG, mFirstReactionFlag, CHybridStochFlag::mIndex, mNumVariableMetabs, CHybridStochFlag::mpNext, mStepsize, mUpdateSet, setState(), and temp.

368 {
369  C_FLOAT64 integrationTime = 0.0;
370  CHybridStochFlag * react = NULL;
371  size_t i;
372 
373  while ((dt - integrationTime) > mStepsize)
374  {
377 
378  for (i = 0; i < mNumVariableMetabs; i++)
379  temp[i] = currentState[i] + (temp[i] * mStepsize);
380 
381  setState(temp);
382  integrationTime += mStepsize;
383  }
384 
387 
388  for (i = 0; i < mNumVariableMetabs; i++)
389  temp[i] = currentState[i] + (temp[i] * (dt - integrationTime));
390 
391  setState(temp);
392 
393  // find the set union of all reactions, which depend on one of the deterministic reactions. The propensities of the stochastic reactions in this set union will be updated later in the method updatePriorityQueue().
394  for (react = mFirstReactionFlag; react != NULL; react = react->mpNext)
395  {
396  const std::set <size_t> & dependents = mDG.getDependents(react->mIndex);
397  std::copy(dependents.begin(), dependents.end(),
398  std::inserter(mUpdateSet, mUpdateSet.begin()));
399  }
400 
401  return;
402 }
void getState(std::vector< C_FLOAT64 > &target)
void calculateDerivative(std::vector< C_FLOAT64 > &deriv)
CHybridStochFlag * mFirstReactionFlag
void setState(std::vector< C_FLOAT64 > &source)
const std::set< size_t > & getDependents(const size_t &node) const
CDependencyGraph mDG
CHybridStochFlag * mpNext
Definition: CHybridMethod.h:92
C_FLOAT64 mStepsize
#define C_FLOAT64
Definition: copasi.h:92
std::vector< C_FLOAT64 > currentState
std::set< size_t > mUpdateSet
std::vector< C_FLOAT64 > temp
size_t mNumVariableMetabs
bool CHybridMethod::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 1714 of file CHybridMethod.cpp.

References CModelEntity::ASSIGNMENT, CCopasiMessage::ERROR, CModel::getCompartments(), CTrajectoryProblem::getDuration(), CModel::getEvents(), CModel::getMetabolites(), CCopasiProblem::getModel(), CModel::getModelValues(), CModel::getNumMetabs(), CModel::getNumModelValues(), CModel::getTotSteps(), CCopasiParameter::getValue(), CTrajectoryMethod::isValidProblem(), MCTrajectoryMethod, mLowerStochLimit, mUpperStochLimit, CModelEntity::ODE, CCopasiParameter::Value::pDOUBLE, CCopasiParameter::Value::pINT, CCopasiVector< T >::size(), and CModel::suitableForStochasticSimulation().

1715 {
1716  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
1717 
1718  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
1719 
1720  if (pTP->getDuration() < 0.0)
1721  {
1722  //back integration not possible
1724  return false;
1725  }
1726 
1727  if (pTP->getModel()->getTotSteps() < 1)
1728  {
1729  //at least one reaction necessary
1731  return false;
1732  }
1733 
1734  //check for rules
1735  size_t i, imax = pTP->getModel()->getNumModelValues();
1736 
1737  for (i = 0; i < imax; ++i)
1738  {
1739  if (pTP->getModel()->getModelValues()[i]->getStatus() == CModelEntity::ODE)
1740  {
1741  //ode rule found
1743  return false;
1744  }
1745 
1746  /* if (pTP->getModel()->getModelValues()[i]->getStatus()==CModelEntity::ASSIGNMENT)
1747  if (pTP->getModel()->getModelValues()[i]->isUsed())
1748  {
1749  //used assignment found
1750  CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 19);
1751  return false;
1752  }*/
1753  }
1754 
1755  imax = pTP->getModel()->getNumMetabs();
1756 
1757  for (i = 0; i < imax; ++i)
1758  {
1759  if (pTP->getModel()->getMetabolites()[i]->getStatus() == CModelEntity::ODE)
1760  {
1761  //ode rule found
1763  return false;
1764  }
1765 
1766  if (pTP->getModel()->getMetabolites()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1767  if (pTP->getModel()->getMetabolites()[i]->isUsed())
1768  {
1769  //used assignment for species found
1771  return false;
1772  }
1773  }
1774 
1775  imax = pTP->getModel()->getCompartments().size();
1776 
1777  for (i = 0; i < imax; ++i)
1778  {
1779  if (pTP->getModel()->getCompartments()[i]->getStatus() == CModelEntity::ODE)
1780  {
1781  //ode rule found
1783  return false;
1784  }
1785  }
1786 
1787  //TODO: rewrite CModel::suitableForStochasticSimulation() to use
1788  // CCopasiMessage
1789  std::string message = pTP->getModel()->suitableForStochasticSimulation();
1790 
1791  if (message != "")
1792  {
1793  //model not suitable, message describes the problem
1794  CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
1795  return false;
1796  }
1797 
1798  /* Max Internal Steps */
1799  if (* getValue("Max Internal Steps").pINT <= 0)
1800  {
1801  //max steps should be at least 1
1803  return false;
1804  }
1805 
1806  /* Lower Limit, Upper Limit */
1807  mLowerStochLimit = * getValue("Lower Limit").pDOUBLE;
1808  mUpperStochLimit = * getValue("Upper Limit").pDOUBLE;
1809 
1811  {
1813  return false;
1814  }
1815 
1816  /* Runge Kutta Stepsize */
1817  if (* getValue("Runge Kutta Stepsize").pDOUBLE <= 0.0)
1818  {
1819  // Runge Kutta Stepsize must be positive
1821  return false;
1822  }
1823 
1824  /* Partitioning Interval */
1825  // nothing to be done here so far
1826 
1827  /* Use Random Seed */
1828  // should be checked in the widget later on
1829 
1830  /* Random Seed */
1831  // nothing to be done here
1832 
1833  //events are not supported at the moment
1834  if (pTP->getModel()->getEvents().size() > 0)
1835  {
1837  return false;
1838  }
1839 
1840  return true;
1841 }
C_FLOAT64 mLowerStochLimit
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)
size_t getTotSteps() const
Definition: CModel.cpp:1136
#define MCTrajectoryMethod
const C_FLOAT64 & getDuration() const
const Value & getValue() const
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
size_t getNumModelValues() const
Definition: CModel.cpp:1139
C_FLOAT64 mUpperStochLimit
std::string suitableForStochasticSimulation() const
Definition: CModel.cpp:3901
CModel * getModel() const
bool CHybridMethod::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 1844 of file CHybridMethod.cpp.

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

Referenced by start().

1845 {
1846  size_t i, imax = pModel->getNumModelValues();
1847 
1848  for (i = 0; i < imax; ++i)
1849  {
1850  if (pModel->getModelValues()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1851  if (pModel->getModelValues()[i]->isUsed())
1852  {
1853  //used assignment found
1854  return true;
1855  }
1856  }
1857 
1858  imax = pModel->getNumMetabs();
1859 
1860  for (i = 0; i < imax; ++i)
1861  {
1862  if (pModel->getMetabolites()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1863  if (pModel->getMetabolites()[i]->isUsed())
1864  {
1865  //used assignment found
1866  return true;
1867  }
1868  }
1869 
1870  imax = pModel->getCompartments().size();
1871 
1872  for (i = 0; i < imax; ++i)
1873  {
1874  if (pModel->getCompartments()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1875  if (pModel->getCompartments()[i]->isUsed())
1876  {
1877  //used assignment found
1878  return true;
1879  }
1880  }
1881 
1882  return false;
1883 }
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 CHybridMethod::outputData ( std::ostream &  os,
C_INT32  mode 
)
protected

Prints out data on standard output.

Prints out data on standard output. Deprecated.

Definition at line 1366 of file CHybridMethod.cpp.

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

1367 {
1368  static unsigned C_INT32 counter = 0;
1369  size_t i;
1370 
1371  switch (mode)
1372  {
1373  case 0:
1374 
1375  if (mOutputCounter == (counter++))
1376  {
1377  counter = 0;
1378  os << mpCurrentState->getTime() << " : ";
1379 
1380  for (i = 0; i < mpMetabolites->size(); i++)
1381  {
1382  os << (*mpMetabolites)[i]->getValue() << " ";
1383  }
1384 
1385  os << std::endl;
1386  }
1387 
1388  break;
1389 
1390  case 1:
1391  os << mpCurrentState->getTime() << " : ";
1392 
1393  for (i = 0; i < mpMetabolites->size(); i++)
1394  {
1395  os << (*mpMetabolites)[i]->getValue() << " ";
1396  }
1397 
1398  os << std::endl;
1399  break;
1400 
1401  default:
1402  ;
1403  }
1404 
1405  return;
1406 }
virtual size_t size() const
#define C_INT32
Definition: copasi.h:90
CCopasiVector< CMetab > * mpMetabolites
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
size_t mOutputCounter
void CHybridMethod::outputDebug ( std::ostream &  os,
size_t  level 
)
protected

Prints out various data on standard output for debugging purposes.

Definition at line 1411 of file CHybridMethod.cpp.

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

1412 {
1413  size_t i, j;
1414  std::set <size_t>::iterator iter, iterEnd;
1415 
1416  os << "outputDebug(" << level << ") *********************************************** BEGIN" << std::endl;
1417 
1418  switch (level)
1419  {
1420  case 0: // Everything !!!
1421  os << " Name: " << CCopasiParameter::getObjectName() << std::endl;
1422  os << "current time: " << mpCurrentState->getTime() << std::endl;
1423  os << "mNumVariableMetabs: " << mNumVariableMetabs << std::endl;
1424  os << "mMaxSteps: " << mMaxSteps << std::endl;
1425  os << "mMaxBalance: " << mMaxBalance << std::endl;
1426  os << "mMaxIntBeforeStep: " << mMaxIntBeforeStep << std::endl;
1427  os << "mpReactions.size(): " << mpReactions->size() << std::endl;
1428 
1429  for (i = 0; i < mpReactions->size(); i++)
1430  os << *(*mpReactions)[i] << std::endl;
1431 
1432  os << "mpMetabolites.size(): " << mpMetabolites->size() << std::endl;
1433 
1434  for (i = 0; i < mpMetabolites->size(); i++)
1435  os << *(*mpMetabolites)[i] << std::endl;
1436 
1437  os << "mStoi: " << std::endl;
1438 
1439  for (i = 0; i < mStoi.numRows(); i++)
1440  {
1441  for (j = 0; j < mStoi.numCols(); j++)
1442  os << mStoi[i][j] << " ";
1443 
1444  os << std::endl;
1445  }
1446 
1447  os << "temp: ";
1448 
1449  for (i = 0; i < mNumVariableMetabs; i++)
1450  os << temp[i] << " ";
1451 
1452  os << std::endl;
1453  os << "curentState: ";
1454 
1455  for (i = 0; i < mNumVariableMetabs; i++)
1456  os << currentState[i] << " ";
1457 
1458  os << std::endl;
1459  os << "k1: ";
1460 
1461  for (i = 0; i < mNumVariableMetabs; i++)
1462  os << k1[i] << " ";
1463 
1464  os << std::endl;
1465  os << "k2: ";
1466 
1467  for (i = 0; i < mNumVariableMetabs; i++)
1468  os << k2[i] << " ";
1469 
1470  os << std::endl;
1471  os << "k3: ";
1472 
1473  for (i = 0; i < mNumVariableMetabs; i++)
1474  os << k3[i] << " ";
1475 
1476  os << std::endl;
1477  os << "k4: ";
1478 
1479  for (i = 0; i < mNumVariableMetabs; i++)
1480  os << k4[i] << " ";
1481 
1482  os << std::endl;
1483  os << "mReactionFlags: " << std::endl;
1484 
1485  for (i = 0; i < mLocalBalances.size(); i++)
1486  os << mReactionFlags[i];
1487 
1488  os << "mFirstReactionFlag: " << std::endl;
1489  if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
1490 
1491  os << "mMetabFlags: " << std::endl;
1492 
1493  for (i = 0; i < mMetabFlags.size(); i++)
1494  {
1495  if (mMetabFlags[i] == LOW)
1496  os << "LOW ";
1497  else
1498  os << "HIGH ";
1499  }
1500 
1501  os << std::endl;
1502  os << "mLocalBalances: " << std::endl;
1503 
1504  for (i = 0; i < mLocalBalances.size(); i++)
1505  {
1506  for (j = 0; j < mLocalBalances[i].size(); j++)
1507  os << mLocalBalances[i][j];
1508 
1509  os << std::endl;
1510  }
1511 
1512  os << "mLocalSubstrates: " << std::endl;
1513 
1514  for (i = 0; i < mLocalSubstrates.size(); i++)
1515  {
1516  for (j = 0; j < mLocalSubstrates[i].size(); j++)
1517  os << mLocalSubstrates[i][j];
1518 
1519  os << std::endl;
1520  }
1521 
1522  os << "mLowerStochLimit: " << mLowerStochLimit << std::endl;
1523  os << "mUpperStochLimit: " << mUpperStochLimit << std::endl;
1524  //deprecated: os << "mOutputCounter: " << mOutputCounter << endl;
1525  os << "mPartitioningInterval: " << mPartitioningInterval << std::endl;
1526  os << "mStepsAfterPartitionSystem: " << mStepsAfterPartitionSystem << std::endl;
1527  os << "mStepsize: " << mStepsize << std::endl;
1528  os << "mMetab2React: " << std::endl;
1529 
1530  for (i = 0; i < mMetab2React.size(); i++)
1531  {
1532  os << i << ": ";
1533 
1534  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; ++iter)
1535  os << *iter << " ";
1536 
1537  os << std::endl;
1538  }
1539 
1540  os << "mAmu: " << std::endl;
1541 
1542  for (i = 0; i < mpReactions->size(); i++)
1543  os << mAmu[i] << " ";
1544 
1545  os << std::endl;
1546  os << "mAmuOld: " << std::endl;
1547 
1548  for (i = 0; i < mpReactions->size(); i++)
1549  os << mAmuOld[i] << " ";
1550 
1551  os << std::endl;
1552  os << "mUpdateSet: " << std::endl;
1553 
1554  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
1555  os << *iter;
1556 
1557  os << std::endl;
1558  os << "mpRandomGenerator: " << mpRandomGenerator << std::endl;
1559  os << "mDG: " << std::endl << mDG;
1560  os << "mPQ: " << std::endl << mPQ;
1561  os << "Particle numbers: " << std::endl;
1562 
1563  for (i = 0; i < mpMetabolites->size(); i++)
1564  {
1565  os << (*mpMetabolites)[i]->getValue() << " ";
1566  }
1567 
1568  os << std::endl;
1569  break;
1570 
1571  case 1: // Variable values only
1572  os << "current time: " << mpCurrentState->getTime() << std::endl;
1573  /*
1574  case 1:
1575  os << "mTime: " << mpCurrentState->getTime() << std::endl;
1576  os << "oldState: ";
1577  for (i = 0; i < mDim; i++)
1578  os << oldState[i] << " ";
1579  os << std::endl;
1580  os << "x: ";
1581  for (i = 0; i < mDim; i++)
1582  os << x[i] << " ";
1583  os << std::endl;
1584  os << "y: ";
1585  for (i = 0; i < mDim; i++)
1586  os << y[i] << " ";
1587  os << std::endl;
1588  os << "increment: ";
1589  for (i = 0; i < mDim; i++)
1590  os << increment[i] << " ";
1591  os << std::endl;*/
1592  os << "temp: ";
1593 
1594  for (i = 0; i < mNumVariableMetabs; i++)
1595  os << temp[i] << " ";
1596 
1597  os << std::endl;
1598  os << "currentState: ";
1599 
1600  for (i = 0; i < mNumVariableMetabs; i++)
1601  os << currentState[i] << " ";
1602 
1603  os << std::endl;
1604  os << "k1: ";
1605 
1606  for (i = 0; i < mNumVariableMetabs; i++)
1607  os << k1[i] << " ";
1608 
1609  os << std::endl;
1610  os << "k2: ";
1611 
1612  for (i = 0; i < mNumVariableMetabs; i++)
1613  os << k2[i] << " ";
1614 
1615  os << std::endl;
1616  os << "k3: ";
1617 
1618  for (i = 0; i < mNumVariableMetabs; i++)
1619  os << k3[i] << " ";
1620 
1621  os << std::endl;
1622  os << "k4: ";
1623 
1624  for (i = 0; i < mNumVariableMetabs; i++)
1625  os << k4[i] << " ";
1626 
1627  os << std::endl;
1628  os << "mReactionFlags: " << std::endl;
1629 
1630  for (i = 0; i < mLocalBalances.size(); i++)
1631  os << mReactionFlags[i];
1632 
1633  os << "mFirstReactionFlag: " << std::endl;
1634  if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
1635 
1636  os << "mMetabFlags: " << std::endl;
1637 
1638  for (i = 0; i < mMetabFlags.size(); i++)
1639  {
1640  if (mMetabFlags[i] == LOW)
1641  os << "LOW ";
1642  else
1643  os << "HIGH ";
1644  }
1645 
1646  os << std::endl;
1647  os << "mAmu: " << std::endl;
1648 
1649  for (i = 0; i < mpReactions->size(); i++)
1650  os << mAmu[i] << " ";
1651 
1652  os << std::endl;
1653  os << "mAmuOld: " << std::endl;
1654 
1655  for (i = 0; i < mpReactions->size(); i++)
1656  os << mAmuOld[i] << " ";
1657 
1658  os << std::endl;
1659  os << "mUpdateSet: " << std::endl;
1660 
1661  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
1662  os << *iter;
1663 
1664  os << std::endl;
1665  os << "mPQ: " << std::endl << mPQ;
1666  os << "Particle numbers: " << std::endl;
1667 
1668  for (i = 0; i < mpMetabolites->size(); i++)
1669  {
1670  os << (*mpMetabolites)[i]->getValue() << " ";
1671  }
1672 
1673  os << std::endl;
1674  break;
1675 
1676  case 2:
1677  break;
1678 
1679  default:
1680  ;
1681  }
1682 
1683  os << "outputDebug(" << level << ") ************************************************* END" << std::endl;
1684  return;
1685 }
CRandom * mpRandomGenerator
C_FLOAT64 mLowerStochLimit
unsigned C_INT32 mMaxSteps
const std::string & getObjectName() const
size_t mMaxIntBeforeStep
virtual size_t size() const
std::vector< C_FLOAT64 > k3
std::vector< metabStatus > mMetabFlags
virtual size_t numRows() const
Definition: CMatrix.h:138
CHybridStochFlag * mFirstReactionFlag
size_t mMaxBalance
CMatrix< C_FLOAT64 > mStoi
CIndexedPriorityQueue mPQ
const CCopasiVectorNS< CReaction > * mpReactions
std::vector< C_FLOAT64 > k4
CDependencyGraph mDG
CCopasiVector< CMetab > * mpMetabolites
std::vector< C_FLOAT64 > mAmu
size_t mStepsAfterPartitionSystem
C_FLOAT64 mStepsize
unsigned C_INT32 mPartitioningInterval
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
std::vector< C_FLOAT64 > currentState
std::vector< C_FLOAT64 > k2
std::vector< std::vector< CHybridBalance > > mLocalSubstrates
std::vector< std::vector< CHybridBalance > > mLocalBalances
std::vector< std::set< size_t > > mMetab2React
C_FLOAT64 mUpperStochLimit
std::vector< C_FLOAT64 > k1
std::set< size_t > mUpdateSet
std::vector< C_FLOAT64 > mAmuOld
virtual size_t numCols() const
Definition: CMatrix.h:144
std::vector< C_FLOAT64 > temp
std::vector< CHybridStochFlag > mReactionFlags
size_t mNumVariableMetabs
void CHybridMethod::partitionSystem ( )
protected

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

Definition at line 1162 of file CHybridMethod.cpp.

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

Referenced by CHybridNextReactionRKMethod::doSingleStep().

1163 {
1164  size_t i;
1165  std::set <size_t>::iterator iter, iterEnd;
1166  C_FLOAT64 key;
1167 
1168  for (i = 0; i < mNumVariableMetabs; i++)
1169  {
1170  if ((mMetabFlags[i] == LOW) && ((*mpMetabolites)[i]->getValue() >= mUpperStochLimit))
1171  {
1172  mMetabFlags[i] = HIGH;
1173 
1174  // go through all corresponding reactions and update flags
1175  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; iter++)
1176  {
1177  mReactionFlags[*iter].mValue--;
1178 
1179  // if reaction gets deterministic, insert it into the linked list of deterministic reactions
1180  if (mReactionFlags[*iter].mValue == 0)
1181  {
1183  mPQ.removeStochReaction(*iter);
1184  }
1185  }
1186  }
1187 
1188  if ((mMetabFlags[i] == HIGH) && ((*mpMetabolites)[i]->getValue() < mLowerStochLimit))
1189  {
1190  mMetabFlags[i] = LOW;
1191  (*mpMetabolites)[i]->setValue(floor((*mpMetabolites)[i]->getValue()));
1192  (*mpMetabolites)[i]->refreshConcentration();
1193 
1194  // go through all corresponding reactions and update flags
1195  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; iter++)
1196  {
1197  if (mReactionFlags[*iter].mValue == 0)
1198  {
1200  /*
1201  mPQ.insertStochReaction(*iter, 1234567.8); // juergen: have to beautify this, number has to be the biggest C_FLOAT64 !!!
1202  */
1203  calculateAmu(*iter);
1204  mAmuOld[*iter] = mAmu[*iter];
1205  key = mpCurrentState->getTime() + generateReactionTime(*iter);
1206  mPQ.insertStochReaction(*iter, key);
1207  }
1208 
1209  mReactionFlags[*iter].mValue++;
1210  }
1211  }
1212  }
1213 
1214  return;
1215 }
C_FLOAT64 mLowerStochLimit
std::vector< metabStatus > mMetabFlags
CIndexedPriorityQueue mPQ
void removeDeterministicReaction(size_t rIndex)
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
const Value & getValue() const
C_FLOAT64 generateReactionTime(size_t rIndex)
size_t removeStochReaction(const size_t index)
CCopasiVector< CMetab > * mpMetabolites
void calculateAmu(size_t rIndex)
std::vector< C_FLOAT64 > mAmu
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
void insertDeterministicReaction(size_t rIndex)
std::vector< std::set< size_t > > mMetab2React
C_FLOAT64 mUpperStochLimit
std::vector< C_FLOAT64 > mAmuOld
std::vector< CHybridStochFlag > mReactionFlags
size_t mNumVariableMetabs
void CHybridMethod::removeDeterministicReaction ( size_t  rIndex)
protected

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

Parameters
rIndexA size_t giving the index of the reaction to be removed from the list of deterministic reactions.

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

Parameters
rIndexA size_t giving the index of the reaction to be removed from the list of deterministic reactions.

Definition at line 1259 of file CHybridMethod.cpp.

References mFirstReactionFlag, CHybridStochFlag::mpPrev, and mReactionFlags.

Referenced by partitionSystem().

1260 {
1261  if (mReactionFlags[rIndex].mpPrev != NULL)
1262  // reaction is deterministic
1263  {
1264  if (&mReactionFlags[rIndex] != mFirstReactionFlag)
1265  // reactionFlag is not the first in the linked list
1266  {
1267  mReactionFlags[rIndex].mpPrev->mpNext = mReactionFlags[rIndex].mpNext;
1268 
1269  if (mReactionFlags[rIndex].mpNext != NULL)
1270  mReactionFlags[rIndex].mpNext->mpPrev = mReactionFlags[rIndex].mpPrev;
1271  }
1272  else
1273  // reactionFlag is the first in the linked list
1274  {
1275  if (mReactionFlags[rIndex].mpNext != NULL) // reactionFlag is not the only one in the linked list
1276  {
1277  mFirstReactionFlag = mReactionFlags[rIndex].mpNext;
1279  }
1280  else // reactionFlag is the only one in the linked list
1281  {
1282  mFirstReactionFlag = NULL;
1283  }
1284  }
1285  }
1286 
1287  mReactionFlags[rIndex].mpPrev = NULL;
1288  mReactionFlags[rIndex].mpNext = NULL;
1289  return;
1290 }
CHybridStochFlag * mFirstReactionFlag
CHybridStochFlag * mpPrev
Definition: CHybridMethod.h:91
std::vector< CHybridStochFlag > mReactionFlags
void CHybridMethod::rungeKutta ( C_FLOAT64  dt)
protected

Does one 4th order RungeKutta step to integrate the system numerically.

Parameters
dtA C_FLOAT64 specifying the stepsize
resultA reference to a vector, into which the result, that is the increment vector, will be written

Definition at line 412 of file CHybridMethod.cpp.

References calculateDerivative(), currentState, getState(), k1, k2, k3, k4, mNumVariableMetabs, setState(), and temp.

Referenced by integrateDeterministicPart().

413 {
414  size_t i;
415 
416  /* save current state */
418 
419  /* k1 step: k1 = dt*f(x(n)) */
420  calculateDerivative(temp); // systemState == x(n)
421 
422  for (i = 0; i < mNumVariableMetabs; i++)
423  {
424  k1[i] = temp[i] * dt;
425  }
426 
427  /* k2 step: k2 = dt*f(x(n) + k1/2) */
428  for (i = 0; i < mNumVariableMetabs; i++)
429  {
430  temp[i] = k1[i] / 2.0 + currentState[i];
431  }
432 
433  setState(temp);
434  calculateDerivative(temp); // systemState == x(n) + k1/2
435 
436  for (i = 0; i < mNumVariableMetabs; i++)
437  {
438  k2[i] = temp[i] * dt;
439  }
440 
441  /* k3 step: k3 = dt*f(x(n) + k2/2) */
442  for (i = 0; i < mNumVariableMetabs; i++)
443  {
444  temp[i] = k2[i] / 2.0 + currentState[i];
445  }
446 
447  setState(temp);
448  calculateDerivative(temp); // systemState == x(n) + k2/2
449 
450  for (i = 0; i < mNumVariableMetabs; i++)
451  {
452  k3[i] = temp[i] * dt;
453  }
454 
455  /* k4 step: k4 = dt*f(x(n) + k3); */
456  for (i = 0; i < mNumVariableMetabs; i++)
457  {
458  temp[i] = k3[i] + currentState[i];
459  }
460 
461  setState(temp);
462  calculateDerivative(temp); // systemState == x(n) + k3
463 
464  for (i = 0; i < mNumVariableMetabs; i++)
465  {
466  k4[i] = temp[i] * dt;
467  }
468 
469  /* Find next position: x(n+1) = x(n) + 1/6*(k1 + 2*k2 + 2*k3 + k4) */
470  for (i = 0; i < mNumVariableMetabs; i++)
471  {
472  temp[i] = currentState[i] + (1.0 / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
473  }
474 
475  setState(temp);
476 
477  return;
478 }
void getState(std::vector< C_FLOAT64 > &target)
void calculateDerivative(std::vector< C_FLOAT64 > &deriv)
std::vector< C_FLOAT64 > k3
void setState(std::vector< C_FLOAT64 > &source)
std::vector< C_FLOAT64 > k4
std::vector< C_FLOAT64 > currentState
std::vector< C_FLOAT64 > k2
std::vector< C_FLOAT64 > k1
std::vector< C_FLOAT64 > temp
size_t mNumVariableMetabs
void CHybridMethod::setState ( std::vector< C_FLOAT64 > &  source)
protected

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

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

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

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

Definition at line 549 of file CHybridMethod.cpp.

References mNumVariableMetabs, mpModel, and CModel::updateSimulatedValues().

Referenced by integrateDeterministicPartEuler(), and rungeKutta().

550 {
551  size_t i;
552 
553  for (i = 0; i < mNumVariableMetabs; i++)
554  {
555  (*mpMetabolites)[i]->setValue(source[i]);
556  (*mpMetabolites)[i]->refreshConcentration();
557  }
558 
559  mpModel->updateSimulatedValues(false); //for assignments
560  return;
561 }
CModel * mpModel
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
size_t mNumVariableMetabs
void CHybridMethod::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 845 of file CHybridMethod.cpp.

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

Referenced by initMethod().

846 {
847  size_t i, j;
848  CHybridBalance newElement;
849  C_INT32 maxBalance = 0;
850  size_t numReactions;
851 
852  numReactions = mpReactions->size();
853  mLocalBalances.clear();
854  mLocalBalances.resize(numReactions);
855  mLocalSubstrates.clear();
856  mLocalSubstrates.resize(numReactions);
857 
858  for (i = 0; i < numReactions; i++)
859  {
860  const CCopasiVector <CChemEqElement> * balances =
861  &(*mpReactions)[i]->getChemEq().getBalances();
862 
863  for (j = 0; j < balances->size(); j++)
864  {
865  newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
866  newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
867  // + 0.5 to get a rounding out of the static_cast to size_t!
868  newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
869 
870  if ((newElement.mpMetabolite->getStatus()) != CModelEntity::FIXED)
871  {
872  if (newElement.mMultiplicity > maxBalance) maxBalance = newElement.mMultiplicity;
873 
874  mLocalBalances[i].push_back(newElement); // element is copied for the push_back
875  }
876  }
877 
878  balances = &(*mpReactions)[i]->getChemEq().getSubstrates();
879 
880  for (j = 0; j < balances->size(); j++)
881  {
882  newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
883  newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
884  // + 0.5 to get a rounding out of the static_cast to size_t!
885  newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
886 
887  mLocalSubstrates[i].push_back(newElement); // element is copied for the push_back
888  }
889  }
890 
891  mMaxBalance = maxBalance;
892  mMaxIntBeforeStep = INT_MAX - 1 - mMaxSteps * mMaxBalance;
893 
894  return;
895 }
unsigned C_INT32 mMaxSteps
CModel * mpModel
size_t mMaxIntBeforeStep
virtual size_t size() const
size_t mMaxBalance
#define C_INT32
Definition: copasi.h:90
Definition: CMetab.h:178
const CCopasiVectorNS< CReaction > * mpReactions
virtual size_t getIndex(const CCopasiObject *pObject) const
std::vector< std::vector< CHybridBalance > > mLocalSubstrates
std::vector< std::vector< CHybridBalance > > mLocalBalances
CMetab * mpMetabolite
const CModelEntity::Status & getStatus() const
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
C_INT32 mMultiplicity
void CHybridMethod::setupDependencyGraph ( )
protected

Sets up the dependency graph

Sets up the dependency graph.

Definition at line 900 of file CHybridMethod.cpp.

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

Referenced by initMethod().

901 {
902  mDG.clear();
903  std::vector< std::set<std::string>* > DependsOn;
904  std::vector< std::set<std::string>* > Affects;
905  size_t numReactions = mpReactions->size();
906  size_t i, j;
907 
908  // Do for each reaction:
909  for (i = 0; i < numReactions; i++)
910  {
911  // Get the set of metabolites which affect the value of amu for this
912  // reaction i.e. the set on which amu depends. This may be more than
913  // the set of substrates, since the kinetics can involve other
914  // reactants, e.g. catalysts. We thus need to step through the
915  // rate function and pick out every reactant which can vary.
916  DependsOn.push_back(getDependsOn(i));
917  // Get the set of metabolites which are affected when this reaction takes place
918  Affects.push_back(getAffects(i));
919  }
920 
921  // For each possible pair of reactions i and j, if the intersection of
922  // Affects(i) with DependsOn(j) is non-empty, add a dependency edge from i to j.
923  for (i = 0; i < numReactions; i++)
924  {
925  for (j = 0; j < numReactions; j++)
926  {
927  // Determine whether the intersection of these two sets is non-empty
928  // Could also do this with set_intersection generic algorithm, but that
929  // would require operator<() to be defined on the set elements.
930 
931  std::set<std::string>::iterator iter = Affects[i]->begin();
932 
933  for (; iter != Affects[i]->end(); iter++)
934  {
935  if (DependsOn[j]->count(*iter))
936  {
937  // The set intersection is non-empty
938  mDG.addDependent(i, j);
939  break;
940  }
941  }
942  }
943 
944  // Ensure that self edges are included
945  //mDG.addDependent(i, i);
946  }
947 
948  // Delete the memory allocated in getDependsOn() and getAffects()
949  // since this is allocated in other functions.
950  for (i = 0; i < numReactions; i++)
951  {
952  delete DependsOn[i];
953  delete Affects[i];
954  }
955 
956  return;
957 }
virtual size_t size() const
void addDependent(const size_t &node, const size_t &dependent)
const CCopasiVectorNS< CReaction > * mpReactions
CDependencyGraph mDG
std::set< std::string > * getAffects(size_t rIndex)
std::set< std::string > * getDependsOn(size_t rIndex)
void CHybridMethod::setupMetab2React ( )
protected

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

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

Definition at line 964 of file CHybridMethod.cpp.

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

Referenced by initMethod().

965 {
966  size_t i, j;
967  size_t metaboliteIndex;
968 
969  // Resize mMetab2React and create an initial set for each metabolite
970  mMetab2React.clear();
971  mMetab2React.resize(mpMetabolites->size());
972 
973  // Iterate over all reactions
974  for (i = 0; i < mLocalBalances.size(); i++)
975  {
976  // Get the set of metabolites which take part in this reaction
977  for (j = 0; j < mLocalBalances[i].size(); j++)
978  {
979  // find metaboliteIndex and insert the reaction into the set
980  metaboliteIndex = mLocalBalances[i][j].mIndex;
981  mMetab2React[metaboliteIndex].insert(i);
982  }
983  }
984 
985  return;
986 }
virtual size_t size() const
CCopasiVector< CMetab > * mpMetabolites
std::vector< std::vector< CHybridBalance > > mLocalBalances
std::vector< std::set< size_t > > mMetab2React
void CHybridMethod::setupMetab2ReactComplete ( )
protected

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

Definition at line 1030 of file CHybridMethod.cpp.

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

1031 {
1032  size_t i, j;
1033 
1034  // Resize mMetab2React and create an initial set for each metabolite
1035  mMetab2React.resize(mpMetabolites->size());
1036 
1037  // Iterate over all metabolites
1038  for (i = 0; i < mpMetabolites->size(); i++)
1039  {
1040  // Iterate over all reactions
1041  for (j = 0; j < mpReactions->size(); j++)
1042  {
1043  mMetab2React[i].insert(j);
1044  }
1045  }
1046 
1047  return;
1048 }
virtual size_t size() const
const CCopasiVectorNS< CReaction > * mpReactions
CCopasiVector< CMetab > * mpMetabolites
std::vector< std::set< size_t > > mMetab2React
void CHybridMethod::setupMetab2ReactPlusModifier ( )
protected

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

Definition at line 993 of file CHybridMethod.cpp.

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

994 {
995  std::vector< std::set<size_t>* > participatesIn;
996  size_t numReactions = mpReactions->size();
997  size_t i;
998 
999  // Resize mMetab2React and create an initial set for each metabolite
1000  mMetab2React.resize(mpMetabolites->size());
1001 
1002  // Do for each reaction:
1003  for (i = 0; i < numReactions; i++)
1004  {
1005  participatesIn.push_back(getParticipatesIn(i));
1006  }
1007 
1008  // Iterate over all reactions
1009  for (i = 0; i < numReactions; i++)
1010  {
1011  // Get the set of metabolites which take part in this reaction
1012  std::set<size_t>::iterator iter = participatesIn[i]->begin();
1013 
1014  for (; iter != participatesIn[i]->end(); iter++)
1015  mMetab2React[*iter].insert(i);
1016  }
1017 
1018  for (i = 0; i < numReactions; i++)
1019  {
1020  delete participatesIn[i];
1021  }
1022 
1023  return;
1024 }
virtual size_t size() const
std::set< size_t > * getParticipatesIn(size_t rIndex)
const CCopasiVectorNS< CReaction > * mpReactions
CCopasiVector< CMetab > * mpMetabolites
std::vector< std::set< size_t > > mMetab2React
void CHybridMethod::setupPartition ( )
protected

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

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

Definition at line 1055 of file CHybridMethod.cpp.

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

Referenced by initMethod().

1056 {
1057  size_t i, j;
1058  CHybridStochFlag * prevFlag;
1059  C_FLOAT64 averageStochLimit = (mUpperStochLimit + mLowerStochLimit) / 2;
1060 
1061  // initialize vector mMetabFlags
1062  mMetabFlags.clear();
1064 
1065  for (i = 0; i < mMetabFlags.size(); i++)
1066  {
1067  if ((*mpMetabolites)[i]->getValue() < averageStochLimit)
1068  {
1069  mMetabFlags[i] = LOW;
1070  (*mpMetabolites)[i]->setValue(floor((*mpMetabolites)[i]->getValue()));
1071  (*mpMetabolites)[i]->refreshConcentration();
1072  }
1073  else
1074  mMetabFlags[i] = HIGH;
1075  }
1076 
1077  // initialize vector mReactionFlags
1078  mReactionFlags.clear();
1079  mReactionFlags.resize(mLocalBalances.size());
1080 
1081  for (i = 0; i < mLocalBalances.size(); i++)
1082  {
1083  mReactionFlags[i].mIndex = i;
1084  mReactionFlags[i].mValue = 0;
1085 
1086  for (j = 0; j < mLocalBalances[i].size(); j++)
1087  {
1088  if (mMetabFlags[mLocalBalances[i][j].mIndex] == LOW)
1089  {
1090  mReactionFlags[i].mValue++;
1091  }
1092  }
1093  }
1094 
1095  mFirstReactionFlag = NULL;
1096  prevFlag = NULL;
1097 
1098  for (i = 0; i < mLocalBalances.size(); i++)
1099  {
1100  if (mReactionFlags[i].mValue == 0)
1101  {
1102  if (mFirstReactionFlag != NULL)
1103  {
1104  prevFlag->mpNext = &mReactionFlags[i];
1105  mReactionFlags[i].mpPrev = prevFlag;
1106  prevFlag = &mReactionFlags[i];
1107  }
1108  else
1109  {
1111  mReactionFlags[i].mpPrev = &mReactionFlags[i]; // Important to distinguish between stochastic (prev == NULL) and deterministic (prev != NULL) reactions
1112  prevFlag = &mReactionFlags[i];
1113  }
1114  }
1115  else
1116  {
1117  mReactionFlags[i].mpPrev = NULL;
1118  mReactionFlags[i].mpNext = NULL;
1119  }
1120  }
1121 
1122  if (prevFlag != NULL)
1123  {
1124  prevFlag->mpNext = NULL;
1125  }
1126 
1127  return;
1128 }
C_FLOAT64 mLowerStochLimit
std::vector< metabStatus > mMetabFlags
CHybridStochFlag * mFirstReactionFlag
const Value & getValue() const
CHybridStochFlag * mpNext
Definition: CHybridMethod.h:92
CCopasiVector< CMetab > * mpMetabolites
#define C_FLOAT64
Definition: copasi.h:92
std::vector< std::vector< CHybridBalance > > mLocalBalances
C_FLOAT64 mUpperStochLimit
std::vector< CHybridStochFlag > mReactionFlags
size_t mNumVariableMetabs
void CHybridMethod::setupPriorityQueue ( C_FLOAT64  startTime = 0.0)
protected

Sets up the priority queue.

Parameters
startTimeThe time at which the simulation starts.

Definition at line 1135 of file CHybridMethod.cpp.

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

Referenced by initMethod().

1136 {
1137  size_t i;
1138  C_FLOAT64 time;
1139 
1140  mPQ.clear();
1142 
1143  for (i = 0; i < mpReactions->size(); i++)
1144  {
1145  if (mReactionFlags[i].mpPrev == NULL) // Reaction is stochastic!
1146  {
1147  calculateAmu(i);
1148  time = startTime + generateReactionTime(i);
1149  mPQ.insertStochReaction(i, time);
1150  }
1151  }
1152 
1153  return;
1154 }
virtual size_t size() const
CIndexedPriorityQueue mPQ
const CCopasiVectorNS< CReaction > * mpReactions
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
C_FLOAT64 generateReactionTime(size_t rIndex)
void calculateAmu(size_t rIndex)
#define C_FLOAT64
Definition: copasi.h:92
std::vector< CHybridStochFlag > mReactionFlags
void initializeIndexPointer(const size_t numberOfReactions)
void CHybridMethod::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 230 of file CHybridMethod.cpp.

References CModel::deterministic, CStateTemplate::getIndex(), CModel::getMetabolitesX(), CCopasiProblem::getModel(), CModel::getModelType(), CModel::getStateTemplate(), CState::getTime(), initMethod(), mDoCorrection, mFirstMetabIndex, mHasAssignments, modelHasAssignments(), CTrajectoryMethod::mpCurrentState, mpModel, CTrajectoryMethod::mpProblem, CModel::setState(), and CModel::updateSimulatedValues().

231 {
232  *mpCurrentState = *initialState;
233 
235  assert(mpModel);
236 
238  mDoCorrection = true;
239  else
240  mDoCorrection = false;
241 
243 
245 
247 
248  mpModel->updateSimulatedValues(false); //for assignments
249  //mpModel->updateNonSimulatedValues(); //for assignments
250 
251  // call init of the simulation method, can be overloaded in derived classes
253 
254  return;
255 }
CModel * mpModel
void initMethod(C_FLOAT64 time)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
CTrajectoryProblem * mpProblem
size_t getIndex(const CModelEntity *entity) const
Definition: CState.cpp:231
size_t mFirstMetabIndex
void setState(const CState &state)
Definition: CModel.cpp:1785
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
const ModelType & getModelType() const
Definition: CModel.cpp:2339
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
CModel * getModel() const
static bool modelHasAssignments(const CModel *pModel)
CTrajectoryMethod::Status CHybridMethod::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 184 of file CHybridMethod.cpp.

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

185 {
186  // write the current state to the model
187  // mpProblem->getModel()->setState(mpCurrentState); // is that correct?
188 
189  // check for possible overflows
190  size_t i;
191  size_t imax;
192 
193  // :TODO: Bug 774: This assumes that the number of variable metabs is the number
194  // of metabs determined by reaction. In addition they are expected at the beginning of the
195  // MetabolitesX which is not the case if we have metabolites of type ODE.
196  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++)
197  if (mpProblem->getModel()->getMetabolitesX()[i]->getValue() >= mMaxIntBeforeStep)
198  {
199  // throw exception or something like that
200  }
201 
202  // do several steps
203  C_FLOAT64 time = mpCurrentState->getTime();
204  C_FLOAT64 endTime = time + deltaT;
205 
206  for (i = 0; ((i < mMaxSteps) && (time < endTime)); i++)
207  {
208  time = doSingleStep(time, endTime);
209  }
210 
211  mpCurrentState->setTime(time);
212 
213  if ((i >= mMaxSteps) && (!mMaxStepsReached))
214  {
215  mMaxStepsReached = true; //only report this message once
216  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.");
217  }
218 
219  // get back the particle numbers
220 
221  /* Set the variable metabolites */
223 
224  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, Dbl++)
225  *Dbl = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
226 
227  return NORMAL;
228 }
unsigned C_INT32 mMaxSteps
size_t mMaxIntBeforeStep
CTrajectoryProblem * mpProblem
size_t mFirstMetabIndex
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
size_t getNumVariableMetabs() const
Definition: CModel.cpp:1121
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
CModel * getModel() const
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
void CHybridMethod::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 625 of file CHybridMethod.cpp.

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

Referenced by CHybridNextReactionRKMethod::doSingleStep().

626 {
627  C_FLOAT64 newTime;
628  size_t index;
629  std::set <size_t>::iterator iter, iterEnd;
630 
631  //if the model contains assignments we use a less efficient loop over all (stochastic) reactions to capture all changes
632  // we do not know the exact dependencies. TODO: this should be changed later in order to get a more efficient update scheme
633  if (mHasAssignments)
634  {
636 
637  for (index = 0; index < mpReactions->size(); index++)
638  {
639  if (mReactionFlags[index].mpPrev == NULL) // Reaction is stochastic!
640  {
641  mAmuOld[index] = mAmu[index];
642  calculateAmu(index);
643 
644  if (mAmuOld[index] != mAmu[index])
645  if (index != rIndex) updateTauMu(index, time);
646  }
647  }
648  }
649  else
650  {
651  // iterate through the set of affected reactions and update the stochastic ones in the priority queue
652  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
653  {
654  if (mReactionFlags[*iter].mpPrev == NULL) // reaction is stochastic!
655  {
656  index = *iter;
657  mAmuOld[index] = mAmu[index];
658  calculateAmu(index);
659 
660  if (*iter != rIndex) updateTauMu(index, time);
661  }
662  }
663  }
664 
665  // draw new random number and update the reaction just fired
666  if ((rIndex != C_INVALID_INDEX) && (mReactionFlags[rIndex].mpPrev == NULL))
667  {
668  // reaction is stochastic
669  newTime = time + generateReactionTime(rIndex);
670  mPQ.updateNode(rIndex, newTime);
671  }
672 
673  // empty the mUpdateSet
674  mUpdateSet.clear();
675  return;
676 }
CModel * mpModel
virtual size_t size() const
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
#define C_INVALID_INDEX
Definition: copasi.h:222
CIndexedPriorityQueue mPQ
void updateTauMu(size_t rIndex, C_FLOAT64 time)
const CCopasiVectorNS< CReaction > * mpReactions
C_FLOAT64 generateReactionTime(size_t rIndex)
void calculateAmu(size_t rIndex)
std::vector< C_FLOAT64 > mAmu
void updateNode(const size_t index, const C_FLOAT64 key)
#define C_FLOAT64
Definition: copasi.h:92
std::set< size_t > mUpdateSet
std::vector< C_FLOAT64 > mAmuOld
std::vector< CHybridStochFlag > mReactionFlags
void CHybridMethod::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 772 of file CHybridMethod.cpp.

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

Referenced by updatePriorityQueue().

773 {
774  C_FLOAT64 newTime;
775 
776  // One must make sure that the calculation yields reasonable results even in the cases where mAmu=0 or mAmuOld=0 or both =0. Therefore mAmuOld=0 is checked. If mAmuOld equals 0, then a new random number has to be drawn, because tau equals inf and the stoch. information is lost. If both values equal 0, then tau should remain inf and the update of the queue can be skipped!
777 
778  if (mAmuOld[rIndex] == 0.0)
779  {
780  if (mAmu[rIndex] != 0.0)
781  {
782  newTime = time + generateReactionTime(rIndex);
783  mPQ.updateNode(rIndex, newTime);
784  }
785  }
786  else
787  {
788  newTime = time + (mAmuOld[rIndex] / mAmu[rIndex]) * (mPQ.getKey(rIndex) - time);
789  mPQ.updateNode(rIndex, newTime);
790  }
791 
792  return;
793 }
CIndexedPriorityQueue mPQ
C_FLOAT64 getKey(const size_t index) const
C_FLOAT64 generateReactionTime(size_t rIndex)
std::vector< C_FLOAT64 > mAmu
void updateNode(const size_t index, const C_FLOAT64 key)
#define C_FLOAT64
Definition: copasi.h:92
std::vector< C_FLOAT64 > mAmuOld

Friends And Related Function Documentation

Member Data Documentation

std::vector<C_FLOAT64> CHybridMethod::currentState
protected
std::vector<C_FLOAT64> CHybridMethod::k1
protected

Definition at line 532 of file CHybridMethod.h.

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

std::vector<C_FLOAT64> CHybridMethod::k2
protected

Definition at line 533 of file CHybridMethod.h.

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

std::vector<C_FLOAT64> CHybridMethod::k3
protected

Definition at line 534 of file CHybridMethod.h.

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

std::vector<C_FLOAT64> CHybridMethod::k4
protected

Definition at line 535 of file CHybridMethod.h.

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

std::vector<C_FLOAT64> CHybridMethod::mAmu
protected

The propensities of the stochastic reactions.

Definition at line 593 of file CHybridMethod.h.

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

std::vector<C_FLOAT64> CHybridMethod::mAmuOld
protected
CDependencyGraph CHybridMethod::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 610 of file CHybridMethod.h.

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

bool CHybridMethod::mDoCorrection
protected

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

Definition at line 510 of file CHybridMethod.h.

Referenced by calculateAmu(), and start().

size_t CHybridMethod::mFirstMetabIndex
protected

index of the first metab in CState

Definition at line 477 of file CHybridMethod.h.

Referenced by start(), and step().

CHybridStochFlag* CHybridMethod::mFirstReactionFlag
protected
bool CHybridMethod::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 642 of file CHybridMethod.h.

Referenced by start(), and updatePriorityQueue().

std::vector<std::vector <CHybridBalance> > CHybridMethod::mLocalBalances
protected

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

Definition at line 555 of file CHybridMethod.h.

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

std::vector<std::vector <CHybridBalance> > CHybridMethod::mLocalSubstrates
protected

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

Definition at line 561 of file CHybridMethod.h.

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

C_FLOAT64 CHybridMethod::mLowerStochLimit
protected

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

Definition at line 567 of file CHybridMethod.h.

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

size_t CHybridMethod::mMaxBalance
protected

maximal increase of a particle number in one step.

Definition at line 500 of file CHybridMethod.h.

Referenced by outputDebug(), and setupBalances().

size_t CHybridMethod::mMaxIntBeforeStep
protected

This is set to maxint - mMaxSteps*mMaxBalance

Definition at line 505 of file CHybridMethod.h.

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

unsigned C_INT32 CHybridMethod::mMaxSteps
protected

Max number of doSingleStep() per step()

Definition at line 482 of file CHybridMethod.h.

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

bool CHybridMethod::mMaxStepsReached
protected

Definition at line 484 of file CHybridMethod.h.

Referenced by initMethod(), and step().

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

Vector of relations between metabolites to reactions.

Definition at line 588 of file CHybridMethod.h.

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

std::vector<metabStatus> CHybridMethod::mMetabFlags
protected

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

Definition at line 549 of file CHybridMethod.h.

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

size_t CHybridMethod::mNumVariableMetabs
protected

Dimension of the system. Total number of metabolites.

Definition at line 472 of file CHybridMethod.h.

Referenced by getState(), initMethod(), integrateDeterministicPartEuler(), outputDebug(), partitionSystem(), rungeKutta(), setState(), and setupPartition().

size_t CHybridMethod::mOutputCounter
protected

Output counter.

Definition at line 635 of file CHybridMethod.h.

Referenced by outputData().

std::ofstream CHybridMethod::mOutputFile
protected

File output stream to write data.

Definition at line 625 of file CHybridMethod.h.

std::string CHybridMethod::mOutputFileName
protected

Output filename.

Definition at line 630 of file CHybridMethod.h.

unsigned C_INT32 CHybridMethod::mPartitioningInterval
protected

The system gets repartitioned after this number of elementary steps.

Definition at line 573 of file CHybridMethod.h.

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

CCopasiVector<CMetab>* CHybridMethod::mpMetabolites
protected
CModel* CHybridMethod::mpModel
protected

Pointer to the model

Definition at line 467 of file CHybridMethod.h.

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

CIndexedPriorityQueue CHybridMethod::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 619 of file CHybridMethod.h.

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

CRandom* CHybridMethod::mpRandomGenerator
protected

The random number generator.

Definition at line 604 of file CHybridMethod.h.

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

const CCopasiVectorNS<CReaction>* CHybridMethod::mpReactions
protected
unsigned C_INT32 CHybridMethod::mRandomSeed
protected

The random seed to use.

Definition at line 495 of file CHybridMethod.h.

Referenced by initMethod().

std::vector<CHybridStochFlag> CHybridMethod::mReactionFlags
protected

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

Definition at line 542 of file CHybridMethod.h.

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

size_t CHybridMethod::mStepsAfterPartitionSystem
protected

Number of elementary steps after the last partitioning.

Definition at line 578 of file CHybridMethod.h.

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

C_FLOAT64 CHybridMethod::mStepsize
protected

Stepsize for the rungeKutta steps of the numerical integrator.

Definition at line 583 of file CHybridMethod.h.

Referenced by CHybridNextReactionRKMethod::doSingleStep(), initMethod(), integrateDeterministicPart(), integrateDeterministicPartEuler(), and outputDebug().

CMatrix<C_FLOAT64> CHybridMethod::mStoi
protected

The stoichometry matrix of the model.

Definition at line 525 of file CHybridMethod.h.

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

std::set<size_t> CHybridMethod::mUpdateSet
protected

Set of the reactions, which must be updated.

Definition at line 599 of file CHybridMethod.h.

Referenced by fireReaction(), initMethod(), integrateDeterministicPart(), integrateDeterministicPartEuler(), outputDebug(), and updatePriorityQueue().

C_FLOAT64 CHybridMethod::mUpperStochLimit
protected
bool CHybridMethod::mUseRandomSeed
protected

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

Definition at line 490 of file CHybridMethod.h.

Referenced by initMethod().

std::vector<C_FLOAT64> CHybridMethod::temp
protected

Vectors to hold the system state and intermediate results

Definition at line 530 of file CHybridMethod.h.

Referenced by initMethod(), integrateDeterministicPartEuler(), outputDebug(), and rungeKutta().


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