COPASI API  4.16.103
Classes | 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
CHybridMethodLSODA Class Referenceabstract

#include <CHybridMethodLSODA.h>

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

Classes

struct  Data
 

Public Member Functions

 CHybridMethodLSODA (const CHybridMethodLSODA &src, const CCopasiContainer *pParent=NULL)
 
virtual bool elevateChildren ()
 
virtual bool isValidProblem (const CCopasiProblem *pProblem)
 
virtual void start (const CState *initialState)
 
virtual Status step (const double &deltaT)
 
 ~CHybridMethodLSODA ()
 
- Public Member Functions inherited from CTrajectoryMethod
 CTrajectoryMethod (const CTrajectoryMethod &src, const CCopasiContainer *pParent=NULL)
 
const CVector< C_INT > & getRoots () const
 
void setCurrentState (CState *currentState)
 
void setProblem (CTrajectoryProblem *problem)
 
virtual void stateChanged ()
 
 ~CTrajectoryMethod ()
 
- Public Member Functions inherited from CCopasiMethod
 CCopasiMethod (const CCopasiMethod &src, const CCopasiContainer *pParent=NULL)
 
const CCopasiMethod::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 CHybridMethodLSODAcreateHybridMethodLSODA ()
 
- 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)
 
 CHybridMethodLSODA (const CCopasiContainer *pParent=NULL)
 
void cleanup ()
 
virtual C_FLOAT64 doSingleStep (C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
 
void evalF (const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
 
void fireReaction (size_t rIndex)
 
C_FLOAT64 generateReactionTime (size_t rIndex)
 
std::set< std::string > * getAffects (size_t rIndex)
 
C_FLOAT64 getDefaultAtol (const CModel *pModel) const
 
std::set< std::string > * getDependsOn (size_t rIndex)
 
std::set< size_t > * getParticipatesIn (size_t rIndex)
 
void getState (std::vector< C_FLOAT64 > &target)
 
void getStochTimeAndIndex (C_FLOAT64 &ds, size_t &rIndex)
 
void initMethod (C_FLOAT64 time)
 
void insertDeterministicReaction (size_t rIndex)
 
void integrateDeterministicPart (C_FLOAT64 ds)
 
void outputData (std::ostream &os, C_INT32 mode)
 
void outputDebug (std::ostream &os, size_t level)
 
void partitionSystem ()
 
void removeDeterministicReaction (size_t rIndex)
 
void setState (std::vector< C_FLOAT64 > &source)
 
void setupBalances ()
 
void setupDependencyGraph ()
 
void setupMetab2React ()
 
void setupMetab2ReactComplete ()
 
void setupMetab2ReactPlusModifier ()
 
void setupPartition ()
 
void setupPriorityQueue (C_FLOAT64 startTime=0.0)
 
void updatePriorityQueue (size_t rIndex, C_FLOAT64 time)
 
void updateTauMu (size_t rIndex, C_FLOAT64 time)
 
- Protected Member Functions inherited from CTrajectoryMethod
 CTrajectoryMethod (const CCopasiMethod::SubType &subType, const CCopasiContainer *pParent=NULL)
 
- Protected Member Functions inherited from CCopasiMethod
 CCopasiMethod (const CCopasiTask::Type &taskType, const SubType &subType, const CCopasiContainer *pParent=NULL)
 
- Protected Member Functions inherited from CCopasiParameterGroup
 CCopasiParameterGroup ()
 
- Protected Member Functions inherited from CCopasiContainer
template<class CType >
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 void EvalF (const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
 
static bool modelHasAssignments (const CModel *pModel)
 

Protected Attributes

std::vector< C_FLOAT64currentState
 
std::vector< C_FLOAT64mAmu
 
std::vector< C_FLOAT64mAmuOld
 
C_FLOAT64 mAtol
 
Data mData
 
bool mDefaultAtol
 
CDependencyGraph mDG
 
bool mDoCorrection
 
CVector< C_FLOAT64mDWork
 
C_FLOAT64 mEndt
 
std::ostringstream mErrorMsg
 
size_t mFirstMetabIndex
 
CHybridLSODAStochFlagmFirstReactionFlag
 
bool mHasAssignments
 
CVector< C_INTmIWork
 
C_INT mJType
 
std::vector< std::vector
< CHybridLSODABalance > > 
mLocalBalances
 
std::vector< std::vector
< CHybridLSODABalance > > 
mLocalSubstrates
 
C_FLOAT64 mLowerStochLimit
 
CLSODA mLSODA
 
C_INT mLsodaStatus
 
C_INT32 mMaxBalance
 
C_INT32 mMaxIntBeforeStep
 
unsigned C_INT32 mMaxSteps
 
bool mMaxStepsReached
 
std::vector< std::set< size_t > > mMetab2React
 
std::vector< 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
 
CStatempState
 
unsigned C_INT32 mRandomSeed
 
std::vector
< CHybridLSODAStochFlag
mReactionFlags
 
bool mReducedModel
 
bool mRestartLSODA
 
C_FLOAT64 mRtol
 
C_INT mState
 
unsigned C_INT32 mStepsAfterPartitionSystem
 
C_FLOAT64 mStepsize
 
CMatrix< C_FLOAT64mStoi
 
C_FLOAT64 mTime
 
std::set< size_t > mUpdateSet
 
C_FLOAT64 mUpperStochLimit
 
bool mUseRandomSeed
 
C_FLOAT64mY
 
CVector< C_FLOAT64mYdot
 
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 118 of file CHybridMethodLSODA.h.

Member Enumeration Documentation

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

Enumerator
LOW 
HIGH 

Definition at line 476 of file CHybridMethodLSODA.h.

Constructor & Destructor Documentation

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

Copy constructor

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

Definition at line 74 of file CHybridMethodLSODA.cpp.

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

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

Destructor.

Definition at line 89 of file CHybridMethodLSODA.cpp.

References cleanup(), and DESTRUCTOR_TRACE.

90 {
91  cleanup();
93 }
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
CHybridMethodLSODA::CHybridMethodLSODA ( const CCopasiContainer pParent = NULL)
protected

Default constructor.

Parameters
constCCopasiContainer * pParent (default: NULL)

CHybridMethodLSODA

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

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

Last change: 09, May 2006

(C) European Media Lab 2003. Default constructor.

Definition at line 63 of file CHybridMethodLSODA.cpp.

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

63  :
65 {
66  assert((void *) &mData == (void *) &mData.dim);
67 
68  mData.pMethod = this;
69 
72 }
CHybridMethodLSODA * pMethod
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49

Member Function Documentation

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

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

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

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

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

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

Definition at line 542 of file CHybridMethodLSODA.cpp.

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

Referenced by evalF().

543 {
544  size_t i;
545  C_INT32 bal = 0;
547 
548  // Calculate all the needed kinetic functions of the deterministic reactions
549  for (j = mFirstReactionFlag; j != NULL; j = j->mpNext)
550  {
551  (*mpReactions)[j->mIndex]->calculateParticleFlux();
552  }
553 
554  // For each metabolite add up the contributions of the det. reactions
555  // the number of rows in mStoi equals the number of non-fixed metabolites!
556  // for (i=0; i<mNumVariableMetabs; i++)
557  for (i = 0; i < mStoi.numRows(); i++)
558  {
559  deriv[i] = 0.0;
560 
561  for (j = mFirstReactionFlag; j != NULL; j = j->mpNext)
562  {
563  // juergen: +0.5 to get a rounding out of the static_cast
564  bal = static_cast<C_INT32>(floor(mStoi[i][j->mIndex] + 0.5));
565  deriv[i] += bal * (*mpReactions)[j->mIndex]->getParticleFlux(); // balance * flux;
566  }
567  }
568 
569  /*
570  for (; i < mNumVariableMetabs; i++) deriv[i] = 0.0; // important to get a correct deriv vector, because mStoi doesn't cover fixed metabolites
571  */
572  return;
573 }
virtual size_t numRows() const
Definition: CMatrix.h:138
CHybridLSODAStochFlag * mpNext
#define C_INT32
Definition: copasi.h:90
CMatrix< C_FLOAT64 > mStoi
CHybridLSODAStochFlag * mFirstReactionFlag
C_INT32 CHybridMethodLSODA::checkModel ( CModel model)
staticprotected

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

Parameters
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 819 of file CHybridMethodLSODA.cpp.

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

820 {
822  CMatrix <C_FLOAT64> mStoi = model->getStoi();
823  C_INT32 multInt;
824  size_t i, j, numReactions = mpReactions->size();
825  C_FLOAT64 multFloat;
826  // size_t metabSize = mpMetabolites->size();
827 
828  for (i = 0; i < numReactions; i++) // for every reaction
829  {
830  // TEST getCompartmentNumber() == 1
831  if ((*mpReactions)[i]->getCompartmentNumber() != 1) return - 1;
832 
833  // TEST isReversible() == 0
834  if ((*mpReactions)[i]->isReversible() != 0) return - 2;
835 
836  // TEST integer stoichometry
837  // Iterate through each the metabolites
838  // juergen: the number of rows of mStoi equals the number of non-fixed metabs!
839  // for (j=0; i<metabSize; j++)
840  for (j = 0; j < mStoi.numRows(); j++)
841  {
842  multFloat = mStoi[j][i];
843  multInt = static_cast<C_INT32>(floor(multFloat + 0.5)); // +0.5 to get a rounding out of the static_cast to int!
844 
845  if ((multFloat - multInt) > INT_EPSILON) return - 3; // INT_EPSILON in CHybridMethod.h
846  }
847  }
848 
849  return 1; // Model is appropriate for hybrid simulation
850 }
virtual size_t size() const
virtual size_t numRows() const
Definition: CMatrix.h:138
#define INT_EPSILON
Definition: CHybridMethod.h:59
#define C_INT32
Definition: copasi.h:90
CMatrix< C_FLOAT64 > mStoi
#define C_FLOAT64
Definition: copasi.h:92
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
const CCopasiVectorNS< CReaction > * mpReactions
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
void CHybridMethodLSODA::cleanup ( )
protected

Cleans up memory, etc.

Definition at line 382 of file CHybridMethodLSODA.cpp.

References mpModel, and mpRandomGenerator.

Referenced by ~CHybridMethodLSODA().

383 {
384  delete mpRandomGenerator;
385  mpRandomGenerator = NULL;
386  mpModel = NULL;
387  return;
388 }
CHybridMethodLSODA * CHybridMethodLSODA::createHybridMethodLSODA ( )
static

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

Creates a HybridMethodLSODA adequate for the problem.

Definition at line 175 of file CHybridMethodLSODA.cpp.

References C_INT32.

Referenced by CTrajectoryMethod::createMethod().

176 {
177  C_INT32 result = 1;
178  CHybridMethodLSODA * method = NULL;
179 
180  switch (result)
181  {
182  /* case - 3: // non-integer stoichometry
183  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 1);
184  break;
185  case - 2: // reversible reaction exists
186  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 2);
187  break;
188 
189  case - 1: // more than one compartment involved
190  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 3);
191  break;*/
192  case 1:
193  default:
194  // Everything alright: Hybrid simulation possible
195  method = new CHybridNextReactionLSODAMethod();
196  break;
197  }
198 
199  return method;
200 }
#define C_INT32
Definition: copasi.h:90
virtual C_FLOAT64 CHybridMethodLSODA::doSingleStep ( C_FLOAT64  currentTime,
C_FLOAT64  endTime 
)
protectedpure virtual

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

Parameters
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 CHybridNextReactionLSODAMethod.

Referenced by step().

bool CHybridMethodLSODA::elevateChildren ( )
virtual

This methods must be called to elevate subgroups to derived objects. The default implementation does nothing.

Returns
bool success

Reimplemented from CCopasiParameterGroup.

Definition at line 166 of file CHybridMethodLSODA.cpp.

References initializeParameter().

167 {
169  return true;
170 }
void CHybridMethodLSODA::EvalF ( const C_INT n,
const C_FLOAT64 t,
const C_FLOAT64 y,
C_FLOAT64 ydot 
)
staticprotected

Definition at line 392 of file CHybridMethodLSODA.cpp.

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

Referenced by integrateDeterministicPart().

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

This evaluates the derivatives for the complete model

Definition at line 395 of file CHybridMethodLSODA.cpp.

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

Referenced by EvalF().

396 {
397  size_t i;
398 
399  CModel * pModel = mpProblem->getModel();
400 
401  mpState->setTime(*t);
402 
403  pModel->setState(*mpState);
404 
405  pModel->updateSimulatedValues(false);
406 
408 
409  for (i = 0; i < mNumVariableMetabs; i++)
410  ydot[i] = temp[i];
411 
412  return;
413 }
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
CTrajectoryProblem * mpProblem
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
void setState(const CState &state)
Definition: CModel.cpp:1785
std::vector< C_FLOAT64 > temp
Definition: CModel.h:50
CModel * getModel() const
void calculateDerivative(std::vector< C_FLOAT64 > &deriv)
void CHybridMethodLSODA::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 600 of file CHybridMethodLSODA.cpp.

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

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

601 {
602  // Change the particle numbers according to which step took place.
603  // First, get the vector of balances in the reaction we've got.
604  // (This vector expresses the number change of each metabolite
605  // in the reaction.) Then step through each balance, using its
606  // multiplicity to calculate a new value for the associated
607  // metabolite. Finally, update the metabolite.
608 
609  size_t i;
610  C_FLOAT64 newNumber;
611 
612  CMetab * pMetab;
613 
614  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
615  {
616  pMetab = mLocalBalances[rIndex][i].mpMetabolite;
617  newNumber = pMetab->getValue() + mLocalBalances[rIndex][i].mMultiplicity;
618 
619  pMetab->setValue(newNumber);
620  pMetab->refreshConcentration();
621  }
622 
623  // insert all dependent reactions into the mUpdateSet
624  const std::set <size_t> & dependents = mDG.getDependents(rIndex);
625  std::copy(dependents.begin(), dependents.end(),
626  std::inserter(mUpdateSet, mUpdateSet.begin()));
627 
628  mRestartLSODA = true;
629 
630  return;
631 }
const std::set< size_t > & getDependents(const size_t &node) const
virtual void setValue(const C_FLOAT64 &value)
Definition: CMetab.h:178
std::vector< std::vector< CHybridLSODABalance > > mLocalBalances
void refreshConcentration()
Definition: CMetab.cpp:281
const C_FLOAT64 & getValue() const
#define C_FLOAT64
Definition: copasi.h:92
CDependencyGraph mDG
std::set< size_t > mUpdateSet
C_FLOAT64 CHybridMethodLSODA::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 693 of file CHybridMethodLSODA.cpp.

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

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

694 {
695  if (mAmu[rIndex] == 0) return std::numeric_limits<C_FLOAT64>::infinity();
696 
698  return - 1 * log(rand2) / mAmu[rIndex];
699 }
std::vector< C_FLOAT64 > mAmu
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
#define C_FLOAT64
Definition: copasi.h:92
std::set< std::string > * CHybridMethodLSODA::getAffects ( size_t  rIndex)
protected

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

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

Definition at line 1349 of file CHybridMethodLSODA.cpp.

References mLocalBalances.

Referenced by setupDependencyGraph().

1350 {
1351  size_t i;
1352  std::set<std::string> *retset = new std::set<std::string>;
1353 
1354  // Get the balances associated with the reaction at this index
1355  // XXX We first get the chemical equation, then the balances, since the getBalances method in CReaction is unimplemented!
1356 
1357  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
1358  {
1359  if (mLocalBalances[rIndex][i].mMultiplicity != 0)
1360  {
1361  retset->insert(mLocalBalances[rIndex][i].mpMetabolite->getKey());
1362  }
1363  }
1364 
1365  return retset;
1366 }
std::vector< std::vector< CHybridLSODABalance > > mLocalBalances
C_FLOAT64 CHybridMethodLSODA::getDefaultAtol ( const CModel pModel) const
protected

Calculate the default absolute tolerance

Parameters
constCModel * pModel
Returns
C_FLOAT64 defaultAtol

Definition at line 361 of file CHybridMethodLSODA.cpp.

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

Referenced by initMethod().

362 {
363 
364  if (!pModel) return 1.0e009;
365 
366  const CCopasiVectorNS< CCompartment > & Compartment = pModel->getCompartments();
367  size_t i, imax;
368 
370 
371  for (i = 0, imax = Compartment.size(); i < imax; i++)
372  if (Compartment[i]->getValue() < Volume) Volume = Compartment[i]->getValue();
373 
374  if (Volume == std::numeric_limits< C_FLOAT64 >::max()) return 1.0e009;
375 
376  return Volume * pModel->getQuantity2NumberFactor() * 1.e-12;
377 }
virtual size_t size() const
const C_FLOAT64 & getQuantity2NumberFactor() const
Definition: CModel.cpp:2354
const Value & getValue() const
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
#define C_FLOAT64
Definition: copasi.h:92
#define max(a, b)
Definition: f2c.h:176
std::set< std::string > * CHybridMethodLSODA::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 1316 of file CHybridMethodLSODA.cpp.

References mpReactions, and CFunctionParameter::PARAMETER.

Referenced by setupDependencyGraph().

1317 {
1318  std::set<std::string> *retset = new std::set<std::string>;
1319 
1320  size_t i, imax = (*mpReactions)[rIndex]->getFunctionParameters().size();
1321  size_t j, jmax;
1322 
1323  for (i = 0; i < imax; ++i)
1324  {
1325  if ((*mpReactions)[rIndex]->getFunctionParameters()[i]->getUsage() == CFunctionParameter::PARAMETER)
1326  continue;
1327 
1328  //metablist = (*mpReactions)[rIndex]->getParameterMappingMetab(i);
1329  const std::vector <std::string> & metabKeylist =
1330  (*mpReactions)[rIndex]->getParameterMappings()[i];
1331  jmax = metabKeylist.size();
1332 
1333  for (j = 0; j < jmax; ++j)
1334  {
1335  retset->insert(metabKeylist[j]);
1336  }
1337  }
1338 
1339  return retset;
1340 }
const CCopasiVectorNS< CReaction > * mpReactions
std::set< size_t > * CHybridMethodLSODA::getParticipatesIn ( size_t  rIndex)
protected

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

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

Definition at line 1375 of file CHybridMethodLSODA.cpp.

Referenced by setupMetab2ReactPlusModifier().

1376 {
1377  std::set<size_t> *retset = new std::set<size_t>;
1378  return retset;
1379 }
void CHybridMethodLSODA::getState ( std::vector< C_FLOAT64 > &  target)
protected

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

Parameters
targetA vector reference of length mNumVariableMetabs, into which the state of the system is written
void CHybridMethodLSODA::getStochTimeAndIndex ( C_FLOAT64 ds,
size_t &  rIndex 
)
protected

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

Parameters
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 586 of file CHybridMethodLSODA.cpp.

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

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

587 {
588  ds = mPQ.topKey();
589  rIndex = mPQ.topIndex();
590  return;
591 }
CIndexedPriorityQueue mPQ
void CHybridMethodLSODA::initializeParameter ( )
private

Intialize the method parameter

Definition at line 95 of file CHybridMethodLSODA.cpp.

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

Referenced by CHybridMethodLSODA(), and elevateChildren().

96 {
97  CCopasiParameter *pParm;
98 
99  assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) MAX_STEPS);
102 
103  assertParameter("Partitioning Interval", CCopasiParameter::UINT, (unsigned C_INT32) PARTITIONING_INTERVAL);
104  assertParameter("Use Random Seed", CCopasiParameter::BOOL, (bool) USE_RANDOM_SEED);
105  assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) RANDOM_SEED);
106 
107  // Check whether we have a method with the old parameter names
108  if ((pParm = getParameter("HYBRID.MaxSteps")) != NULL)
109  {
110  setValue("Max Internal Steps", *pParm->getValue().pUINT);
111  removeParameter("HYBRID.MaxSteps");
112 
113  if ((pParm = getParameter("HYBRID.LowerStochLimit")) != NULL)
114  {
115  setValue("Lower Limit", *pParm->getValue().pDOUBLE);
116  removeParameter("HYBRID.LowerStochLimit");
117  }
118 
119  if ((pParm = getParameter("HYBRID.UpperStochLimit")) != NULL)
120  {
121  setValue("Upper Limit", *pParm->getValue().pDOUBLE);
122  removeParameter("HYBRID.UpperStochLimit");
123  }
124 
125  if ((pParm = getParameter("HYBRID.PartitioningInterval")) != NULL)
126  {
127  setValue("Partitioning Interval", *pParm->getValue().pUINT);
128  removeParameter("HYBRID.PartitioningInterval");
129  }
130 
131  if ((pParm = getParameter("UseRandomSeed")) != NULL)
132  {
133  setValue("Use Random Seed", *pParm->getValue().pBOOL);
134  removeParameter("UseRandomSeed");
135  }
136 
137  if ((pParm = getParameter("")) != NULL)
138  {
139  setValue("Random Seed", *pParm->getValue().pUINT);
140  removeParameter("");
141  }
142  }
143 
144  /* LSODA ****************************************************************************/
145 
146  addParameter("Partitioning Stepsize",
148 
149  addParameter("Relative Tolerance",
151 
152  addParameter("Use Default Absolute Tolerance",
153  CCopasiParameter::BOOL, (bool) true);
154 
155  addParameter("Absolute Tolerance",
157 
158  addParameter("Max Internal Steps (LSODA)",
159  CCopasiParameter::UINT, (unsigned C_INT32) 10000);
160 
161  // These parameters are no longer supported.
162  removeParameter("Adams Max Order");
163  removeParameter("BDF Max Order");
164 }
#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
#define PARTITIONING_STEPSIZE
bool setValue(const std::string &name, const CType &value)
unsigned C_INT32 * pUINT
CCopasiParameter * getParameter(const std::string &name)
#define C_FLOAT64
Definition: copasi.h:92
bool addParameter(const CCopasiParameter &parameter)
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
#define LOWER_STOCH_LIMIT
Definition: CHybridMethod.h:60
void CHybridMethodLSODA::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 290 of file CHybridMethodLSODA.cpp.

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

Referenced by start().

291 {
293  mAmu.clear();
294  mAmu.resize(mpReactions->size());
295  mAmuOld.clear();
296  mAmuOld.resize(mpReactions->size());
298  //mNumVariableMetabs = mpModel->getNumVariableMetabs(); // ind + dep metabs, without fixed metabs
300 
301  temp.clear();
302  temp.resize(mNumVariableMetabs);
303  //currentState.clear();
304  //currentState.resize(mNumVariableMetabs);
305 
306  /* get configuration data */
307  mMaxSteps = * getValue("Max Internal Steps").pUINT;
308  mLowerStochLimit = * getValue("Lower Limit").pDOUBLE;
309  mUpperStochLimit = * getValue("Upper Limit").pDOUBLE;
310  mStepsize = * getValue("Partitioning Stepsize").pDOUBLE;
311  mPartitioningInterval = * getValue("Partitioning Interval").pUINT;
312  mUseRandomSeed = * getValue("Use Random Seed").pBOOL;
313  mRandomSeed = * getValue("Random Seed").pUINT;
314 
316 
317  mStoi = mpModel->getStoi();
319  mUpdateSet.clear();
320 
321  setupBalances(); // initialize mLocalBalances and mLocalSubstrates (has to be called first!)
322  setupDependencyGraph(); // initialize mDG
323  setupMetab2React(); // initialize mMetab2React
324  setupPartition(); // initialize mReactionFlags
325  setupPriorityQueue(start_time); // initialize mPQ
326 
327  mMaxStepsReached = false;
328 
329  /* CONFIGURE LSODA ***********************************************************************/
330 
332 
334 
335  mRtol = * getValue("Relative Tolerance").pUDOUBLE;
336 
337  mDefaultAtol = * getValue("Use Default Absolute Tolerance").pBOOL;
338 
339  if (mDefaultAtol)
340  {
342  setValue("Absolute Tolerance", mAtol);
343  }
344  else
345  mAtol = * getValue("Absolute Tolerance").pUDOUBLE;
346 
347  mDWork.resize(22 + mData.dim * std::max<C_INT>(16, mData.dim + 9));
348  mDWork[4] = mDWork[5] = mDWork[6] = mDWork[7] = mDWork[8] = mDWork[9] = 0.0;
349  mIWork.resize(20 + mData.dim);
350 
351  mIWork[4] = mIWork[6] = mIWork[9] = 0;
352  mIWork[5] = * getValue("Max Internal Steps (LSODA)").pUINT;
353  mIWork[7] = 12;
354  mIWork[8] = 5;
355 
356  /***********************************************************************************/
357 
358  return;
359 }
C_FLOAT64 getDefaultAtol(const CModel *pModel) const
#define C_INT
Definition: copasi.h:115
CVector< C_FLOAT64 > mDWork
unsigned C_INT32 mPartitioningInterval
std::vector< C_FLOAT64 > mAmu
unsigned C_INT32 mMaxSteps
virtual size_t size() const
unsigned C_INT32 mStepsAfterPartitionSystem
CTrajectoryProblem * mpProblem
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
CVector< C_INT > mIWork
unsigned C_INT32 mRandomSeed
CCopasiVector< CMetab > * mpMetabolites
CMatrix< C_FLOAT64 > mStoi
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
CVector< C_FLOAT64 > mYdot
const Value & getValue() const
bool setValue(const std::string &name, const CType &value)
unsigned C_INT32 * pUINT
std::vector< C_FLOAT64 > temp
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
const CCopasiVectorNS< CReaction > * mpReactions
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
CModel * getModel() const
std::vector< C_FLOAT64 > mAmuOld
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
std::set< size_t > mUpdateSet
void CHybridMethodLSODA::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 1242 of file CHybridMethodLSODA.cpp.

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

Referenced by partitionSystem().

1243 {
1244  if (mReactionFlags[rIndex].mpPrev == NULL)
1245  // reaction is stochastic (avoids double insertions)
1246  {
1247  if (mFirstReactionFlag != NULL)
1248  // there are deterministic reactions already
1249  {
1251  mReactionFlags[rIndex].mpNext = mFirstReactionFlag;
1254  }
1255  else
1256  {
1257  // there are no deterministic reactions
1258  // Important to distinguish between stochastic (prev == NULL) and deterministic (prev != NULL) reactions
1259  mReactionFlags[rIndex].mpPrev = &mReactionFlags[rIndex];
1261  }
1262 
1263  mAmu[rIndex] = 0.0;
1264  mAmuOld[rIndex] = 0.0;
1265  }
1266 
1267  return;
1268 }
std::vector< C_FLOAT64 > mAmu
std::vector< CHybridLSODAStochFlag > mReactionFlags
CHybridLSODAStochFlag * mpPrev
std::vector< C_FLOAT64 > mAmuOld
CHybridLSODAStochFlag * mFirstReactionFlag
void CHybridMethodLSODA::integrateDeterministicPart ( C_FLOAT64  deltaT)
protected

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

Parameters
dsA C_FLOAT64 specifying the stepsize.

Definition at line 422 of file CHybridMethodLSODA.cpp.

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

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

423 {
424  CHybridLSODAStochFlag * react = NULL;
425 
426  /* RESET LSODA **************************************************************/
427 
428  if (mRestartLSODA)
429  {
430  mRestartLSODA = false;
431  mLsodaStatus = 1;
432  }
433 
434  mState = 1;
435  mJType = 2;
436 
437  mErrorMsg.str("");
439 
440  C_INT one = 1;
441  C_INT DSize = (C_INT) mDWork.size();
442  C_INT ISize = (C_INT) mIWork.size();
443 
444  //mpState = new CState (mpProblem->getModel()->getState());
446 
447  mY = mpState->beginIndependent();// + mFirstMetabIndex - 1;
448  //for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, mY++)
449  // *mY = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
450 
451  mTime = mpState->getTime();
452 
453  C_FLOAT64 EndTime = mTime + deltaT;
454 
455  C_FLOAT64 tdist , d__1, d__2, w0;
456  tdist = fabs(deltaT);
457  d__1 = fabs(mTime), d__2 = fabs(EndTime);
458  w0 = std::max(d__1, d__2);
459 
460  if (tdist < std::numeric_limits< C_FLOAT64 >::epsilon() * 2. * w0) //just do nothing
461  {
462  //mTime = mTime + deltaT;
463  //mpState->setTime(mTime);
464  //*mpCurrentState = *mpState;
465  mpProblem->getModel()->setTime(EndTime);
466 
467  return;
468  }
469 
470  if (!mData.dim) //just do nothing if there are no variables
471  {
472  //mTime = mTime + deltaT;
473  mpProblem->getModel()->setTime(EndTime);
474  //mpState->setTime(mTime);
475  //*mpCurrentState = *mpState;
476 
477  return;
478  }
479 
480  mLSODA(&EvalF , // 1. evaluate F
481  &mData.dim, // 2. number of variables
482  mY , // 3. the array of current concentrations
483  &mTime , // 4. the current time
484  &EndTime , // 5. the final time
485  &one , // 6. scalar error control
486  &mRtol , // 7. relative tolerance array
487  &mAtol , // 8. absolute tolerance array
488  &mState , // 9. output by overshoot & interpolatation
489  &mLsodaStatus , // 10. the state control variable
490  &one , // 11. futher options (one)
491  mDWork.array() , // 12. the double work array
492  &DSize , // 13. the double work array size
493  mIWork.array() , // 14. the int work array
494  &ISize , // 15. the int work array size
495  NULL , // 16. evaluate J (not given)
496  &mJType); // 17. the type of jacobian calculate (2)
497 
498  if (mLsodaStatus == -1) mLsodaStatus = 2;
499 
500  if ((mLsodaStatus != 1) && (mLsodaStatus != 2) && (mLsodaStatus != -1))
501  {
503  }
504 
506  //*mpCurrentState = *mpState;
508 
509  /* size_t i;
510 
511  for (i = 0; i < mNumVariableMetabs; i++)
512  {
513  (*mpMetabolites)[i]->setValue(mY[i]);
514  (*mpMetabolites)[i]->refreshConcentration();
515  } */
516 
517  mpModel->updateSimulatedValues(false); //for assignments
518 
519  //pdelete(mpState);
520 
521  // 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().
522 
523  for (react = mFirstReactionFlag; react != NULL; react = react->mpNext)
524  {
525  const std::set <size_t> & dependents = mDG.getDependents(react->mIndex);
526  std::copy(dependents.begin(), dependents.end(),
527  std::inserter(mUpdateSet, mUpdateSet.begin()));
528  }
529 
530  return;
531 }
#define C_INT
Definition: copasi.h:115
CVector< C_FLOAT64 > mDWork
void setOstream(std::ostream &os)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
const std::set< size_t > & getDependents(const size_t &node) const
CTrajectoryProblem * mpProblem
CHybridLSODAStochFlag * mpNext
std::ostringstream mErrorMsg
CVector< C_INT > mIWork
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
#define MCTrajectoryMethod
void setState(const CState &state)
Definition: CModel.cpp:1785
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
size_t size() const
Definition: CVector.h:100
void setTime(const C_FLOAT64 &time)
Definition: CModel.cpp:1187
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
CDependencyGraph mDG
const CState & getState() const
Definition: CModel.cpp:1771
CModel * getModel() const
CHybridLSODAStochFlag * mFirstReactionFlag
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
std::set< size_t > mUpdateSet
#define max(a, b)
Definition: f2c.h:176
bool CHybridMethodLSODA::isValidProblem ( const CCopasiProblem pProblem)
virtual

Check if the method is suitable for this problem

Returns
bool suitability of the method

Reimplemented from CTrajectoryMethod.

Definition at line 1674 of file CHybridMethodLSODA.cpp.

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

1675 {
1676  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
1677 
1678  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
1679 
1680  if (pTP->getDuration() < 0.0)
1681  {
1682  //back integration not possible
1684  return false;
1685  }
1686 
1687  //check for rules
1688  size_t i, imax = pTP->getModel()->getNumModelValues();
1689 
1690  for (i = 0; i < imax; ++i)
1691  {
1692  if (pTP->getModel()->getModelValues()[i]->getStatus() == CModelEntity::ODE)
1693  {
1694  //ode rule found
1696  return false;
1697  }
1698  }
1699 
1700  imax = pTP->getModel()->getNumMetabs();
1701 
1702  for (i = 0; i < imax; ++i)
1703  {
1704  if (pTP->getModel()->getMetabolites()[i]->getStatus() == CModelEntity::ODE)
1705  {
1706  //ode rule found
1708  return false;
1709  }
1710  }
1711 
1712  imax = pTP->getModel()->getCompartments().size();
1713 
1714  for (i = 0; i < imax; ++i)
1715  {
1716  if (pTP->getModel()->getCompartments()[i]->getStatus() == CModelEntity::ODE)
1717  {
1718  //ode rule found
1720  return false;
1721  }
1722  }
1723 
1724  std::string message = pTP->getModel()->suitableForStochasticSimulation();
1725 
1726  if (message != "")
1727  {
1728  //model not suitable, message describes the problem
1729  CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
1730  return false;
1731  }
1732 
1733  mLowerStochLimit = * getValue("Lower Limit").pDOUBLE;
1734  mUpperStochLimit = * getValue("Upper Limit").pDOUBLE;
1735 
1738 
1739  //events are not supported at the moment
1740  if (pTP->getModel()->getEvents().size() > 0)
1741  {
1743  return false;
1744  }
1745 
1746  return true;
1747 }
CCopasiVectorN< CEvent > & getEvents()
Definition: CModel.cpp:1110
const CCopasiVector< CMetab > & getMetabolites() const
Definition: CModel.cpp:1051
const CCopasiVectorN< CModelValue > & getModelValues() const
Definition: CModel.cpp:1060
virtual size_t size() const
size_t getNumMetabs() const
Definition: CModel.cpp:1118
virtual bool isValidProblem(const CCopasiProblem *pProblem)
#define MCTrajectoryMethod
const C_FLOAT64 & getDuration() const
const Value & getValue() const
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
size_t getNumModelValues() const
Definition: CModel.cpp:1139
std::string suitableForStochasticSimulation() const
Definition: CModel.cpp:3901
CModel * getModel() const
bool CHybridMethodLSODA::modelHasAssignments ( const CModel pModel)
staticprotected

tests if the model contains a global value with an assignment rule that is used in calculations

Definition at line 1750 of file CHybridMethodLSODA.cpp.

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

Referenced by start().

1751 {
1752  size_t i, imax = pModel->getNumModelValues();
1753 
1754  for (i = 0; i < imax; ++i)
1755  {
1756  if (pModel->getModelValues()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1757  if (pModel->getModelValues()[i]->isUsed())
1758  {
1759  //used assignment found
1760  return true;
1761  }
1762  }
1763 
1764  imax = pModel->getNumMetabs();
1765 
1766  for (i = 0; i < imax; ++i)
1767  {
1768  if (pModel->getMetabolites()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1769  if (pModel->getMetabolites()[i]->isUsed())
1770  {
1771  //used assignment found
1772  return true;
1773  }
1774  }
1775 
1776  imax = pModel->getCompartments().size();
1777 
1778  for (i = 0; i < imax; ++i)
1779  {
1780  if (pModel->getCompartments()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1781  if (pModel->getCompartments()[i]->isUsed())
1782  {
1783  //used assignment found
1784  return true;
1785  }
1786  }
1787 
1788  return false;
1789 }
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 CHybridMethodLSODA::outputData ( std::ostream &  os,
C_INT32  mode 
)
protected

Prints out data on standard output.

Prints out data on standard output. Deprecated.

Definition at line 1384 of file CHybridMethodLSODA.cpp.

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

1385 {
1386  static unsigned C_INT32 counter = 0;
1387  size_t i;
1388 
1389  switch (mode)
1390  {
1391  case 0:
1392 
1393  if (mOutputCounter == (counter++))
1394  {
1395  counter = 0;
1396  os << mpCurrentState->getTime() << " : ";
1397 
1398  for (i = 0; i < mpMetabolites->size(); i++)
1399  {
1400  os << (*mpMetabolites)[i]->getValue() << " ";
1401  }
1402 
1403  os << std::endl;
1404  }
1405 
1406  break;
1407 
1408  case 1:
1409  os << mpCurrentState->getTime() << " : ";
1410 
1411  for (i = 0; i < mpMetabolites->size(); i++)
1412  {
1413  os << (*mpMetabolites)[i]->getValue() << " ";
1414  }
1415 
1416  os << std::endl;
1417  break;
1418 
1419  default:
1420  ;
1421  }
1422 
1423  return;
1424 }
virtual size_t size() const
#define C_INT32
Definition: copasi.h:90
CCopasiVector< CMetab > * mpMetabolites
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
void CHybridMethodLSODA::outputDebug ( std::ostream &  os,
size_t  level 
)
protected

Prints out various data on standard output for debugging purposes.

Definition at line 1429 of file CHybridMethodLSODA.cpp.

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

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

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

Definition at line 1175 of file CHybridMethodLSODA.cpp.

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

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

1176 {
1177  size_t i;
1178  std::set <size_t>::iterator iter, iterEnd;
1179  C_FLOAT64 key;
1180 
1181  for (i = 0; i < mNumVariableMetabs; i++)
1182  {
1183  if ((mMetabFlags[i] == LOW) && ((*mpMetabolites)[i]->getValue() >= mUpperStochLimit))
1184  {
1185  mRestartLSODA = true;
1186  mMetabFlags[i] = HIGH;
1187 
1188  // go through all corresponding reactions and update flags
1189  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; iter++)
1190  {
1191  mReactionFlags[*iter].mValue--;
1192 
1193  // if reaction gets deterministic, insert it into the linked list of deterministic reactions
1194  if (mReactionFlags[*iter].mValue == 0)
1195  {
1197  mPQ.removeStochReaction(*iter);
1198  }
1199  }
1200  }
1201 
1202  if ((mMetabFlags[i] == HIGH) && ((*mpMetabolites)[i]->getValue() < mLowerStochLimit))
1203  {
1204  mRestartLSODA = true;
1205  mMetabFlags[i] = LOW;
1206  (*mpMetabolites)[i]->setValue(floor((*mpMetabolites)[i]->getValue()));
1207  (*mpMetabolites)[i]->refreshConcentration();
1208 
1209  // go through all corresponding reactions and update flags
1210  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; iter++)
1211  {
1212  if (mReactionFlags[*iter].mValue == 0)
1213  {
1215  /*
1216  mPQ.insertStochReaction(*iter, 1234567.8); // juergen: have to beautify this, number has to be the biggest C_FLOAT64 !!!
1217  */
1218  calculateAmu(*iter);
1219  mAmuOld[*iter] = mAmu[*iter];
1220  //key = mpCurrentState->getTime() + generateReactionTime(*iter); NATALIA
1221 
1222  key = mpProblem->getModel()->getTime() + generateReactionTime(*iter);
1223 
1224  mPQ.insertStochReaction(*iter, key);
1225  }
1226 
1227  mReactionFlags[*iter].mValue++;
1228  }
1229  }
1230  }
1231 
1232  return;
1233 }
std::vector< std::set< size_t > > mMetab2React
std::vector< C_FLOAT64 > mAmu
CIndexedPriorityQueue mPQ
CTrajectoryProblem * mpProblem
C_FLOAT64 generateReactionTime(size_t rIndex)
CCopasiVector< CMetab > * mpMetabolites
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
const Value & getValue() const
std::vector< CHybridLSODAStochFlag > mReactionFlags
size_t removeStochReaction(const size_t index)
const C_FLOAT64 & getTime() const
Definition: CModel.cpp:1190
void removeDeterministicReaction(size_t rIndex)
#define C_FLOAT64
Definition: copasi.h:92
std::vector< metabStatus > mMetabFlags
CModel * getModel() const
void insertDeterministicReaction(size_t rIndex)
std::vector< C_FLOAT64 > mAmuOld
void calculateAmu(size_t rIndex)
void CHybridMethodLSODA::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 1277 of file CHybridMethodLSODA.cpp.

References mFirstReactionFlag, CHybridLSODAStochFlag::mpPrev, and mReactionFlags.

Referenced by partitionSystem().

1278 {
1279  if (mReactionFlags[rIndex].mpPrev != NULL)
1280  // reaction is deterministic
1281  {
1282  if (&mReactionFlags[rIndex] != mFirstReactionFlag)
1283  // reactionFlag is not the first in the linked list
1284  {
1285  mReactionFlags[rIndex].mpPrev->mpNext = mReactionFlags[rIndex].mpNext;
1286 
1287  if (mReactionFlags[rIndex].mpNext != NULL)
1288  mReactionFlags[rIndex].mpNext->mpPrev = mReactionFlags[rIndex].mpPrev;
1289  }
1290  else
1291  // reactionFlag is the first in the linked list
1292  {
1293  if (mReactionFlags[rIndex].mpNext != NULL) // reactionFlag is not the only one in the linked list
1294  {
1295  mFirstReactionFlag = mReactionFlags[rIndex].mpNext;
1297  }
1298  else // reactionFlag is the only one in the linked list
1299  {
1300  mFirstReactionFlag = NULL;
1301  }
1302  }
1303  }
1304 
1305  mReactionFlags[rIndex].mpPrev = NULL;
1306  mReactionFlags[rIndex].mpNext = NULL;
1307  return;
1308 }
std::vector< CHybridLSODAStochFlag > mReactionFlags
CHybridLSODAStochFlag * mpPrev
CHybridLSODAStochFlag * mFirstReactionFlag
void CHybridMethodLSODA::setState ( std::vector< C_FLOAT64 > &  source)
protected

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

Parameters
sourceA vector reference with length mNumVariableMetabs, holding the state of the system to be set in the model
void CHybridMethodLSODA::setupBalances ( )
protected

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

Definition at line 858 of file CHybridMethodLSODA.cpp.

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

Referenced by initMethod().

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

Sets up the dependency graph

Sets up the dependency graph.

Definition at line 913 of file CHybridMethodLSODA.cpp.

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

Referenced by initMethod().

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

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

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

Definition at line 977 of file CHybridMethodLSODA.cpp.

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

Referenced by initMethod().

978 {
979  size_t i, j;
980  size_t metaboliteIndex;
981 
982  // Resize mMetab2React and create an initial set for each metabolite
983  mMetab2React.clear();
984  mMetab2React.resize(mpMetabolites->size());
985 
986  // Iterate over all reactions
987  for (i = 0; i < mLocalBalances.size(); i++)
988  {
989  // Get the set of metabolites which take part in this reaction
990  for (j = 0; j < mLocalBalances[i].size(); j++)
991  {
992  // find metaboliteIndex and insert the reaction into the set
993  metaboliteIndex = mLocalBalances[i][j].mIndex;
994  mMetab2React[metaboliteIndex].insert(i);
995  }
996  }
997 
998  return;
999 }
std::vector< std::set< size_t > > mMetab2React
virtual size_t size() const
CCopasiVector< CMetab > * mpMetabolites
std::vector< std::vector< CHybridLSODABalance > > mLocalBalances
void CHybridMethodLSODA::setupMetab2ReactComplete ( )
protected

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

Definition at line 1043 of file CHybridMethodLSODA.cpp.

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

1044 {
1045  size_t i, j;
1046 
1047  // Resize mMetab2React and create an initial set for each metabolite
1048  mMetab2React.resize(mpMetabolites->size());
1049 
1050  // Iterate over all metabolites
1051  for (i = 0; i < mpMetabolites->size(); i++)
1052  {
1053  // Iterate over all reactions
1054  for (j = 0; j < mpReactions->size(); j++)
1055  {
1056  mMetab2React[i].insert(j);
1057  }
1058  }
1059 
1060  return;
1061 }
std::vector< std::set< size_t > > mMetab2React
virtual size_t size() const
CCopasiVector< CMetab > * mpMetabolites
const CCopasiVectorNS< CReaction > * mpReactions
void CHybridMethodLSODA::setupMetab2ReactPlusModifier ( )
protected

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

Definition at line 1006 of file CHybridMethodLSODA.cpp.

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

1007 {
1008  std::vector< std::set<size_t>* > participatesIn;
1009  size_t numReactions = mpReactions->size();
1010  size_t i;
1011 
1012  // Resize mMetab2React and create an initial set for each metabolite
1013  mMetab2React.resize(mpMetabolites->size());
1014 
1015  // Do for each reaction:
1016  for (i = 0; i < numReactions; i++)
1017  {
1018  participatesIn.push_back(getParticipatesIn(i));
1019  }
1020 
1021  // Iterate over all reactions
1022  for (i = 0; i < numReactions; i++)
1023  {
1024  // Get the set of metabolites which take part in this reaction
1025  std::set<size_t>::iterator iter = participatesIn[i]->begin();
1026 
1027  for (; iter != participatesIn[i]->end(); iter++)
1028  mMetab2React[*iter].insert(i);
1029  }
1030 
1031  for (i = 0; i < numReactions; i++)
1032  {
1033  delete participatesIn[i];
1034  }
1035 
1036  return;
1037 }
std::vector< std::set< size_t > > mMetab2React
virtual size_t size() const
CCopasiVector< CMetab > * mpMetabolites
std::set< size_t > * getParticipatesIn(size_t rIndex)
const CCopasiVectorNS< CReaction > * mpReactions
void CHybridMethodLSODA::setupPartition ( )
protected

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

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

Definition at line 1068 of file CHybridMethodLSODA.cpp.

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

Referenced by initMethod().

1069 {
1070  size_t i, j;
1071  CHybridLSODAStochFlag * prevFlag;
1072  C_FLOAT64 averageStochLimit = (mUpperStochLimit + mLowerStochLimit) / 2;
1073 
1074  // initialize vector mMetabFlags
1075  mMetabFlags.clear();
1077 
1078  for (i = 0; i < mMetabFlags.size(); i++)
1079  {
1080  if ((*mpMetabolites)[i]->getValue() < averageStochLimit)
1081  {
1082  mMetabFlags[i] = LOW;
1083  (*mpMetabolites)[i]->setValue(floor((*mpMetabolites)[i]->getValue()));
1084  (*mpMetabolites)[i]->refreshConcentration();
1085  }
1086  else
1087  mMetabFlags[i] = HIGH;
1088  }
1089 
1090  // initialize vector mReactionFlags
1091  mReactionFlags.clear();
1092  mReactionFlags.resize(mLocalBalances.size());
1093 
1094  for (i = 0; i < mLocalBalances.size(); i++)
1095  {
1096  mReactionFlags[i].mIndex = i;
1097  mReactionFlags[i].mValue = 0;
1098 
1099  for (j = 0; j < mLocalBalances[i].size(); j++)
1100  {
1101  if (mMetabFlags[mLocalBalances[i][j].mIndex] == LOW)
1102  {
1103  mReactionFlags[i].mValue++;
1104  }
1105  }
1106  }
1107 
1108  mFirstReactionFlag = NULL;
1109  prevFlag = NULL;
1110 
1111  for (i = 0; i < mLocalBalances.size(); i++)
1112  {
1113  if (mReactionFlags[i].mValue == 0)
1114  {
1115  if (mFirstReactionFlag != NULL)
1116  {
1117  prevFlag->mpNext = &mReactionFlags[i];
1118  mReactionFlags[i].mpPrev = prevFlag;
1119  prevFlag = &mReactionFlags[i];
1120  }
1121  else
1122  {
1124  mReactionFlags[i].mpPrev = &mReactionFlags[i]; // Important to distinguish between stochastic (prev == NULL) and deterministic (prev != NULL) reactions
1125  prevFlag = &mReactionFlags[i];
1126  }
1127  }
1128  else
1129  {
1130  mReactionFlags[i].mpPrev = NULL;
1131  mReactionFlags[i].mpNext = NULL;
1132  }
1133  }
1134 
1135  if (prevFlag != NULL)
1136  {
1137  prevFlag->mpNext = NULL;
1138  }
1139 
1140  return;
1141 }
CHybridLSODAStochFlag * mpNext
CCopasiVector< CMetab > * mpMetabolites
std::vector< std::vector< CHybridLSODABalance > > mLocalBalances
const Value & getValue() const
std::vector< CHybridLSODAStochFlag > mReactionFlags
#define C_FLOAT64
Definition: copasi.h:92
std::vector< metabStatus > mMetabFlags
CHybridLSODAStochFlag * mFirstReactionFlag
void CHybridMethodLSODA::setupPriorityQueue ( C_FLOAT64  startTime = 0.0)
protected

Sets up the priority queue.

Parameters
startTimeThe time at which the simulation starts.

Definition at line 1148 of file CHybridMethodLSODA.cpp.

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

Referenced by initMethod().

1149 {
1150  size_t i;
1151  C_FLOAT64 time;
1152 
1153  mPQ.clear();
1155 
1156  for (i = 0; i < mpReactions->size(); i++)
1157  {
1158  if (mReactionFlags[i].mpPrev == NULL) // Reaction is stochastic!
1159  {
1160  calculateAmu(i);
1161  time = startTime + generateReactionTime(i);
1162  mPQ.insertStochReaction(i, time);
1163  }
1164  }
1165 
1166  return;
1167 }
CIndexedPriorityQueue mPQ
virtual size_t size() const
C_FLOAT64 generateReactionTime(size_t rIndex)
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
std::vector< CHybridLSODAStochFlag > mReactionFlags
#define C_FLOAT64
Definition: copasi.h:92
const CCopasiVectorNS< CReaction > * mpReactions
void initializeIndexPointer(const size_t numberOfReactions)
void calculateAmu(size_t rIndex)
void CHybridMethodLSODA::start ( const CState initialState)
virtual

This instructs the method to prepare for integration starting with the initialState given.

Parameters
const CState *initialState

Reimplemented from CTrajectoryMethod.

Definition at line 253 of file CHybridMethodLSODA.cpp.

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

254 {
255  //*mpCurrentState = *initialState;
256 
257  mpProblem->getModel()->setState(*initialState);
258  mRestartLSODA = true;
259 
261 
263  assert(mpModel);
264 
265  if (mpModel->getModelType() == CModel::deterministic)
266  mDoCorrection = true;
267  else
268  mDoCorrection = false;
269 
271 
272  mFirstMetabIndex = mpModel->getStateTemplate().getIndex(mpModel->getMetabolitesX()[0]);
273 
274 //mpProblem->getModel()->setState(*mpCurrentState);
275 
276  // call init of the simulation method, can be overloaded in derived classes
277  // mDim[1] = (C_INT) (void *) this;
279 
280  return;
281 }
static bool modelHasAssignments(const CModel *pModel)
Definition: CState.h:305
CTrajectoryProblem * mpProblem
void setState(const CState &state)
Definition: CModel.cpp:1785
const C_FLOAT64 & getTime() const
Definition: CModel.cpp:1190
void initMethod(C_FLOAT64 time)
const CState & getState() const
Definition: CModel.cpp:1771
CModel * getModel() const
CTrajectoryMethod::Status CHybridMethodLSODA::step ( const double &  deltaT)
virtual

This instructs the method to calculate a time step of deltaT starting with the current state, i.e., the result of the previous step. The new state (after deltaT) is expected in the current state. The return value is the actual timestep taken.

Parameters
const double &deltaT
Returns
Status status

Reimplemented from CTrajectoryMethod.

Definition at line 202 of file CHybridMethodLSODA.cpp.

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

203 {
204  // write the current state to the model
205  // mpProblem->getModel()->setState(mpCurrentState); // is that correct?
206 
207  // check for possible overflows
208  size_t i;
209  size_t imax;
210 
211  // :TODO: Bug 774: This assumes that the number of variable metabs is the number
212  // of metabs determined by reaction. In addition they are expected at the beginning of the
213  // MetabolitesX which is not the case if we have metabolites of type ODE.
214  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++)
215  if (mpProblem->getModel()->getMetabolitesX()[i]->getValue() >= mMaxIntBeforeStep)
216  {
217  // throw exception or something like that
218  }
219 
220  // do several steps
221  C_FLOAT64 time = mpProblem->getModel()->getTime();
222  C_FLOAT64 endTime = time + deltaT;
223 
224  for (i = 0; ((i < mMaxSteps) && (time < endTime)); i++)
225  {
226  time = doSingleStep(time, endTime);
227  }
228 
229  //mpCurrentState->setTime(time);
230 
231  if ((i >= mMaxSteps) && (!mMaxStepsReached))
232  {
233  mMaxStepsReached = true; //only report this message once
234  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.");
235  }
236 
237  // get back the particle numbers
238 
239  /* Set the variable metabolites */
241 
242  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, Dbl++)
243  *Dbl = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
244 
245  //the task expects the result in mpCurrentState
247 
248  mpCurrentState->setTime(time);
249 
250  return NORMAL;
251 }
unsigned C_INT32 mMaxSteps
CTrajectoryProblem * mpProblem
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
size_t getNumVariableMetabs() const
Definition: CModel.cpp:1121
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
const C_FLOAT64 & getTime() const
Definition: CModel.cpp:1190
#define C_FLOAT64
Definition: copasi.h:92
const CState & getState() const
Definition: CModel.cpp:1771
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
CModel * getModel() const
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
void CHybridMethodLSODA::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 641 of file CHybridMethodLSODA.cpp.

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

Referenced by CHybridNextReactionLSODAMethod::doSingleStep().

642 {
643  C_FLOAT64 newTime;
644  size_t index;
645  std::set <size_t>::iterator iter, iterEnd;
646 
647  //if the model contains assignments we use a less efficient loop over all (stochastic) reactions to capture all changes
648  // we do not know the exact dependencies. TODO: this should be changed later in order to get a more efficient update scheme
649  if (mHasAssignments)
650  {
652 
653  for (index = 0; index < mpReactions->size(); index++)
654  {
655  if (mReactionFlags[index].mpPrev == NULL) // Reaction is stochastic!
656  {
657  mAmuOld[index] = mAmu[index];
658  calculateAmu(index);
659 
660  if (mAmuOld[index] != mAmu[index])
661  if (index != rIndex) updateTauMu(index, time);
662  }
663  }
664  }
665  else
666  {
667  // iterate through the set of affected reactions and update the stochastic ones in the priority queue
668  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
669  {
670  if (mReactionFlags[*iter].mpPrev == NULL) // reaction is stochastic!
671  {
672  index = *iter;
673  mAmuOld[index] = mAmu[index];
674  calculateAmu(index);
675 
676  if (*iter != rIndex) updateTauMu(index, time);
677  }
678  }
679  }
680 
681  // draw new random number and update the reaction just fired
682  if ((rIndex != C_INVALID_INDEX) && (mReactionFlags[rIndex].mpPrev == NULL))
683  {
684  // reaction is stochastic
685  newTime = time + generateReactionTime(rIndex);
686  mPQ.updateNode(rIndex, newTime);
687  }
688 
689  // empty the mUpdateSet
690  mUpdateSet.clear();
691  return;
692 }
std::vector< C_FLOAT64 > mAmu
CIndexedPriorityQueue mPQ
virtual size_t size() const
void updateTauMu(size_t rIndex, C_FLOAT64 time)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
#define C_INVALID_INDEX
Definition: copasi.h:222
C_FLOAT64 generateReactionTime(size_t rIndex)
std::vector< CHybridLSODAStochFlag > mReactionFlags
void updateNode(const size_t index, const C_FLOAT64 key)
#define C_FLOAT64
Definition: copasi.h:92
const CCopasiVectorNS< CReaction > * mpReactions
std::vector< C_FLOAT64 > mAmuOld
std::set< size_t > mUpdateSet
void calculateAmu(size_t rIndex)
void CHybridMethodLSODA::updateTauMu ( size_t  rIndex,
C_FLOAT64  time 
)
protected

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

Parameters
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 787 of file CHybridMethodLSODA.cpp.

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

Referenced by updatePriorityQueue().

788 {
789  C_FLOAT64 newTime;
790 
791  // 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!
792 
793  if (mAmuOld[rIndex] == 0.0)
794  {
795  if (mAmu[rIndex] != 0.0)
796  {
797  newTime = time + generateReactionTime(rIndex);
798  mPQ.updateNode(rIndex, newTime);
799  }
800  }
801  else
802  {
803  newTime = time + (mAmuOld[rIndex] / mAmu[rIndex]) * (mPQ.getKey(rIndex) - time);
804  mPQ.updateNode(rIndex, newTime);
805  }
806 
807  return;
808 }
std::vector< C_FLOAT64 > mAmu
CIndexedPriorityQueue mPQ
C_FLOAT64 generateReactionTime(size_t rIndex)
C_FLOAT64 getKey(const size_t index) const
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> CHybridMethodLSODA::currentState
protected

Definition at line 545 of file CHybridMethodLSODA.h.

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

The propensities of the stochastic reactions.

Definition at line 603 of file CHybridMethodLSODA.h.

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

std::vector<C_FLOAT64> CHybridMethodLSODA::mAmuOld
protected
C_FLOAT64 CHybridMethodLSODA::mAtol
protected

Absolute tolerance.

Definition at line 693 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and integrateDeterministicPart().

Data CHybridMethodLSODA::mData
protected

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

Definition at line 648 of file CHybridMethodLSODA.h.

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

bool CHybridMethodLSODA::mDefaultAtol
protected

Definition at line 688 of file CHybridMethodLSODA.h.

Referenced by initMethod().

CDependencyGraph CHybridMethodLSODA::mDG
protected

The graph of reactions and their dependent reactions. When a reaction is executed, the propensities for each of its dependents must be updated.

Definition at line 620 of file CHybridMethodLSODA.h.

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

bool CHybridMethodLSODA::mDoCorrection
protected

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

Definition at line 524 of file CHybridMethodLSODA.h.

Referenced by calculateAmu(), and start().

CVector< C_FLOAT64 > CHybridMethodLSODA::mDWork
protected

Definition at line 702 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and integrateDeterministicPart().

C_FLOAT64 CHybridMethodLSODA::mEndt
protected

Requested end time.

Definition at line 668 of file CHybridMethodLSODA.h.

std::ostringstream CHybridMethodLSODA::mErrorMsg
protected

Definition at line 698 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

size_t CHybridMethodLSODA::mFirstMetabIndex
protected

index of the first metab in CState

Definition at line 491 of file CHybridMethodLSODA.h.

Referenced by start(), and step().

CHybridLSODAStochFlag* CHybridMethodLSODA::mFirstReactionFlag
protected
bool CHybridMethodLSODA::mHasAssignments
protected

Indicates whether the model has global quantities with assignment rules. If it has, we will use a less efficient way to update the model state to handle this.

Definition at line 727 of file CHybridMethodLSODA.h.

Referenced by start(), and updatePriorityQueue().

CVector< C_INT > CHybridMethodLSODA::mIWork
protected

Definition at line 703 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and integrateDeterministicPart().

C_INT CHybridMethodLSODA::mJType
protected

Definition at line 704 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

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

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

Definition at line 565 of file CHybridMethodLSODA.h.

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

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

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

Definition at line 571 of file CHybridMethodLSODA.h.

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

C_FLOAT64 CHybridMethodLSODA::mLowerStochLimit
protected

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

Definition at line 577 of file CHybridMethodLSODA.h.

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

CLSODA CHybridMethodLSODA::mLSODA
protected

Definition at line 700 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

C_INT CHybridMethodLSODA::mLsodaStatus
protected

LSODA state.

Definition at line 673 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

C_INT32 CHybridMethodLSODA::mMaxBalance
protected

maximal increase of a particle number in one step.

Definition at line 514 of file CHybridMethodLSODA.h.

Referenced by outputDebug(), and setupBalances().

C_INT32 CHybridMethodLSODA::mMaxIntBeforeStep
protected

This is set to maxint - mMaxSteps*mMaxBalance

Definition at line 519 of file CHybridMethodLSODA.h.

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

unsigned C_INT32 CHybridMethodLSODA::mMaxSteps
protected

Max number of doSingleStep() per step()

Definition at line 496 of file CHybridMethodLSODA.h.

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

bool CHybridMethodLSODA::mMaxStepsReached
protected

Definition at line 498 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and step().

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

Vector of relations between metabolites to reactions.

Definition at line 598 of file CHybridMethodLSODA.h.

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

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

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

Definition at line 559 of file CHybridMethodLSODA.h.

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

size_t CHybridMethodLSODA::mNumVariableMetabs
protected

Dimension of the system. Total number of metabolites.

Definition at line 486 of file CHybridMethodLSODA.h.

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

size_t CHybridMethodLSODA::mOutputCounter
protected

Output counter.

Definition at line 720 of file CHybridMethodLSODA.h.

Referenced by outputData().

std::ofstream CHybridMethodLSODA::mOutputFile
protected

File output stream to write data.

Definition at line 710 of file CHybridMethodLSODA.h.

std::string CHybridMethodLSODA::mOutputFileName
protected

Output filename.

Definition at line 715 of file CHybridMethodLSODA.h.

unsigned C_INT32 CHybridMethodLSODA::mPartitioningInterval
protected

The system gets repartitioned after this number of elementary steps.

Definition at line 583 of file CHybridMethodLSODA.h.

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

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

Pointer to the model

Definition at line 481 of file CHybridMethodLSODA.h.

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

CIndexedPriorityQueue CHybridMethodLSODA::mPQ
protected

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

Definition at line 629 of file CHybridMethodLSODA.h.

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

CRandom* CHybridMethodLSODA::mpRandomGenerator
protected

The random number generator.

Definition at line 614 of file CHybridMethodLSODA.h.

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

const CCopasiVectorNS<CReaction>* CHybridMethodLSODA::mpReactions
protected
CState* CHybridMethodLSODA::mpState
protected

A pointer to the current state in complete model view.

Definition at line 642 of file CHybridMethodLSODA.h.

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

unsigned C_INT32 CHybridMethodLSODA::mRandomSeed
protected

The random seed to use.

Definition at line 509 of file CHybridMethodLSODA.h.

Referenced by initMethod().

std::vector<CHybridLSODAStochFlag> CHybridMethodLSODA::mReactionFlags
protected

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

Definition at line 552 of file CHybridMethodLSODA.h.

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

bool CHybridMethodLSODA::mReducedModel
protected

Definition at line 678 of file CHybridMethodLSODA.h.

bool CHybridMethodLSODA::mRestartLSODA
protected

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

Definition at line 637 of file CHybridMethodLSODA.h.

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

C_FLOAT64 CHybridMethodLSODA::mRtol
protected

Relative tolerance.

Definition at line 683 of file CHybridMethodLSODA.h.

Referenced by initMethod(), and integrateDeterministicPart().

C_INT CHybridMethodLSODA::mState
protected

Definition at line 701 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

unsigned C_INT32 CHybridMethodLSODA::mStepsAfterPartitionSystem
protected

Number of elementary steps after the last partitioning.

Definition at line 588 of file CHybridMethodLSODA.h.

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

C_FLOAT64 CHybridMethodLSODA::mStepsize
protected

Partitioning stepsize

Definition at line 593 of file CHybridMethodLSODA.h.

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

CMatrix<C_FLOAT64> CHybridMethodLSODA::mStoi
protected

The stoichometry matrix of the model.

Definition at line 539 of file CHybridMethodLSODA.h.

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

C_FLOAT64 CHybridMethodLSODA::mTime
protected

Current time.

Definition at line 663 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

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

Set of the reactions, which must be updated.

Definition at line 609 of file CHybridMethodLSODA.h.

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

C_FLOAT64 CHybridMethodLSODA::mUpperStochLimit
protected
bool CHybridMethodLSODA::mUseRandomSeed
protected

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

Definition at line 504 of file CHybridMethodLSODA.h.

Referenced by initMethod().

C_FLOAT64* CHybridMethodLSODA::mY
protected

Pointer to the array with left hand side values.

Definition at line 653 of file CHybridMethodLSODA.h.

Referenced by integrateDeterministicPart().

CVector< C_FLOAT64 > CHybridMethodLSODA::mYdot
protected

Vector containig the derivatives after calling eval

Definition at line 658 of file CHybridMethodLSODA.h.

Referenced by initMethod().

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

Vectors to hold the system state and intermediate results

Definition at line 544 of file CHybridMethodLSODA.h.

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


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