COPASI API  4.16.103
CStochNextReactionMethod.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2002 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <cmath>
16 
17 #include "copasi.h"
19 #include "CTrajectoryMethod.h"
20 #include "model/CCompartment.h"
21 #include "model/CModel.h"
22 
24  : CStochMethod()
25 {}
26 
28 {
29  setupPriorityQueue(start_time);
30 }
31 
33 {
34  C_FLOAT64 steptime = mPQ.topKey();
35 
36  // We need to throw an exception if mA0 is NaN
37  if (isnan(steptime))
38  {
40  }
41 
42  if (steptime >= endTime)
43  {
44  return endTime;
45  }
46  else
47  {
48  size_t reaction_index = mPQ.topIndex();
49  updateSystemState(reaction_index, steptime);
50  updatePriorityQueue(reaction_index, steptime);
51  //printDebugInfo();
52  return steptime;
53  }
54 }
55 
56 //void CStochNextReactionMethod::printDebugInfo()
57 //{
58 // cout << mPQ << endl;
59 //}
60 
62 {
63  C_FLOAT64 time;
64 
65  mPQ.clear();
66 
67  for (size_t i = 0; i < mpModel->getReactions().size(); i++)
68  {
69  time = start_time + generateReactionTime(i);
70  mPQ.pushPair(i, time);
71  }
72 
73  mPQ.buildHeap();
74 }
75 
77 {
78  //first the new time for the currently fired reaction
79  C_FLOAT64 new_time = time + generateReactionTime(reaction_index);
80  mPQ.updateNode(reaction_index, new_time);
81 
82  //now the updates for the other reactions (whose propensities may have changed)
83 
84  //if the model contains assignment we use a less efficient loop over all reactions since
85  // we do not know the exact dependencies
86  if (mHasAssignments)
87  {
88  size_t i;
89 
90  for (i = 0; i < mNumReactions; ++i)
91  {
92  if (i != reaction_index)
93  {
94  if (mAmu[i] == mAmuOld[i])
95  continue;
96 
97  C_FLOAT64 new_time;
98 
99  if (mAmuOld[i] > 0)
100  {
101  new_time = time + (mAmuOld[i] / mAmu[i]) * (mPQ.getKey(i) - time);
102  }
103  else
104  {
105  new_time = time + generateReactionTime(i);
106  }
107 
108  mPQ.updateNode(i, new_time);
109  }
110  }
111  }
112  else
113  {
114  const std::set<size_t> & dep_nodes = mDG.getDependents(reaction_index);
115  std::set<size_t>::const_iterator di;
116 
117  for (di = dep_nodes.begin(); di != dep_nodes.end(); ++di)
118  {
119  if (*di != reaction_index)
120  {
121  size_t index = *di;
122  C_FLOAT64 new_time;
123 
124  if (mAmuOld[index] > 0)
125  {
126  new_time = time + (mAmuOld[index] / mAmu[index]) * (mPQ.getKey(index) - time);
127  }
128  else
129  {
130  new_time = time + generateReactionTime(index);
131  }
132 
133  mPQ.updateNode(index, new_time);
134  }
135  }
136  }
137 }
std::vector< C_FLOAT64 > mAmu
Definition: CStochMethod.h:57
void setupPriorityQueue(C_FLOAT64 start_time=0)
virtual size_t size() const
std::vector< C_FLOAT64 > mAmuOld
Definition: CStochMethod.h:67
CModel * mpModel
Definition: CStochMethod.h:246
const std::set< size_t > & getDependents(const size_t &node) const
void initMethod(C_FLOAT64 start_time)
CDependencyGraph mDG
Definition: CStochMethod.h:252
#define C_UNUSED(p)
Definition: copasi.h:220
C_FLOAT64 getKey(const size_t index) const
C_FLOAT64 doSingleStep(C_FLOAT64 time, C_FLOAT64 endTime)
size_t mNumReactions
Definition: CStochMethod.h:269
#define MCTrajectoryMethod
void updatePriorityQueue(size_t reaction_index, C_FLOAT64 time)
bool mHasAssignments
Definition: CStochMethod.h:79
C_INT32 updateSystemState(size_t reaction_index, const C_FLOAT64 &time)
void updateNode(const size_t index, const C_FLOAT64 key)
#define C_FLOAT64
Definition: copasi.h:92
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
C_FLOAT64 generateReactionTime()
size_t pushPair(const size_t index, const C_FLOAT64 key)