COPASI API  4.16.103
CHybridMethod.h
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2012 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) 2003 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * CHybridMethod
17  *
18  * This class implements an hybrid algorithm to simulate a biochemical
19  * system over time.
20  *
21  * File name: CHybridMethod.h
22  * Author: Juergen Pahle
23  * Email: juergen.pahle@eml-r.villa-bosch.de
24  *
25  * Last change: 15, December 2004
26  *
27  * (C) European Media Lab 2003.
28  */
29 /**
30  * Partition the system into a deterministic part and a stochastic part.
31  * That is, every reaction is either classified deterministic or
32  * stochastic. Deterministic reactions involve only those metabolites (on
33  * substrate and product side), which have a high particle number.
34  * Therefore it is legal to integrate this part of the system with e.g. a
35  * numerical integrator. The concentration and relative particle number
36  * change should be low enough, so that the probabilities of all the
37  * reactions in the system show only little changes. The stochastic
38  * reactions must be simulated with an exact stochastic method (e.g. next
39  * reaction method (Gibson)), because their firing changes the reaction
40  * probabilities in the system significantly.
41  */
42 #ifndef COPASI_CHybridMethod
43 #define COPASI_CHybridMethod
44 
45 /* INCLUDES ******************************************************************/
47 #include <set>
48 #include <vector>
49 #include <iostream>
50 #include <fstream>
51 #include "utilities/CVersion.h"
52 #include "utilities/CMatrix.h"
56 
57 /* DEFINE ********************************************************************/
58 #define MAX_STEPS 1000000
59 #define INT_EPSILON 0.1
60 #define LOWER_STOCH_LIMIT 800 //800
61 #define UPPER_STOCH_LIMIT 1000 //1000
62 #define RUNGE_KUTTA_STEPSIZE 0.001
63 #define PARTITIONING_INTERVAL 1
64 #define OUTPUT_COUNTER 100
65 //#define DEFAULT_OUTPUT_FILE "hybrid.output"
66 #define SUBTYPE 1
67 #define USE_RANDOM_SEED false
68 #define RANDOM_SEED 1
69 
70 /* CLASSES *******************************************************************/
71 
72 class CMetab;
73 class CTrajectoryProblem;
74 class CState;
75 class CModel;
76 class CVersion;
77 class CReaction;
78 class CRandom;
80 class CDependencyGraph;
81 
82 /**
83  * Element of the array mReactionFlags. The deterministic reactions are
84  * organized in a linked list to be able to iterate over them quickly.
85  */
87 {
88 public:
89  size_t mIndex;
90  size_t mValue;
93 
94  // insert operator
95  friend std::ostream & operator<<(std::ostream & os, const CHybridStochFlag & d);
96 };
97 
98 /**
99  * Internal representation of the balances of each reaction. The index of
100  * each metabolite in the reaction is provided.
101  */
103 {
104 public:
105  size_t mIndex;
108 
109  // insert operator
110  friend std::ostream & operator<<(std::ostream & os, const CHybridBalance & d);
111 };
112 
114 {
115  friend CTrajectoryMethod *
117 
118  /* PUBLIC METHODS **********************************************************/
119 
120 public:
121  /**
122  * Copy constructor
123  * @param const CHybridMethod & src
124  * @param const CCopasiContainer * pParent (default: NULL)
125  */
126  CHybridMethod(const CHybridMethod & src,
127  const CCopasiContainer * pParent = NULL);
128 
129  /**
130  * Destructor.
131  */
132  ~CHybridMethod();
133 
134  /**
135  * This methods must be called to elevate subgroups to
136  * derived objects. The default implementation does nothing.
137  * @return bool success
138  */
139  virtual bool elevateChildren();
140 
141  /**
142  * Creates a HybridMethod adequate for the problem.
143  * (only CHybridNextReactionRKMethod so far)
144  */
146 
147  /**
148  * This instructs the method to calculate a time step of deltaT
149  * starting with the current state, i.e., the result of the previous
150  * step.
151  * The new state (after deltaT) is expected in the current state.
152  * The return value is the actual timestep taken.
153  * @param "const double &" deltaT
154  * @return Status status
155  */
156  virtual Status step(const double & deltaT);
157 
158  /**
159  * This instructs the method to prepare for integration
160  * starting with the initialState given.
161  * @param "const CState *" initialState
162  */
163  virtual void start(const CState * initialState);
164 
165  /**
166  * Check if the method is suitable for this problem
167  * @return bool suitability of the method
168  */
169  virtual bool isValidProblem(const CCopasiProblem * pProblem);
170 
171  /* PROTECTED METHODS *******************************************************/
172 
173 protected:
174  /**
175  * Default constructor.
176  * @param const CCopasiContainer * pParent (default: NULL)
177  */
178  CHybridMethod(const CCopasiContainer * pParent = NULL);
179 
180  /**
181  * Initializes the solver.
182  * @param time the current time
183  */
184  void initMethod(C_FLOAT64 time);
185 
186  /**
187  * Cleans up memory, etc.
188  */
189  void cleanup();
190 
191  /* VIRTUAL METHODS *********************************************************/
192 
193  /**
194  * Simulates the system over the next interval of time. The current time
195  * and the end time of the current step() are given as arguments.
196  *
197  * @param currentTime A C_FLOAT64 specifying the current time
198  * @param endTime A C_FLOAT64 specifying the end time of the step()
199  * @return A C_FLOAT giving the new time
200  */
201  virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime) = 0;
202 
203  /* DETERMINISTIC STUFF *****************************************************/
204 
205  /**
206  * Integrates the deterministic reactions of the system over the
207  * specified time interval.
208  *
209  * @param ds A C_FLOAT64 specifying the stepsize.
210  */
212 
213  /**
214  * Integrates the deterministic reactions of the system over the
215  * specified time interval.
216  *
217  * @param ds A C_FLOAT64 specifying the stepsize.
218  */
220 
221  /**
222  * Does one 4th order RungeKutta step to integrate the system
223  * numerically.
224  *
225  * @param dt A C_FLOAT64 specifying the stepsize
226  * @param result A reference to a vector, into which the result, that is
227  * the increment vector, will be written
228  */
229  void rungeKutta(C_FLOAT64 dt);
230 
231  /**
232  * Calculates the derivative of the system and writes it into the vector
233  * deriv. Length of deriv must be mNumVariableMetabs.
234  * CAUTION: Only deterministic reactions are taken into account. That is,
235  * this is only the derivative of the deterministic part of the system.
236  *
237  * @param deriv A vector reference of length mNumVariableMetabs, into
238  * which the derivative is written
239  */
240  void calculateDerivative(std::vector <C_FLOAT64> & deriv);
241 
242  /**
243  * Gathers the state of the system into the array target. Later on CState
244  * should be used for this. Length of the array target must be
245  * mNumVariableMetabs.
246  *
247  * @param target A vector reference of length mNumVariableMetabs, into
248  * which the state of the system is written
249  */
250  void getState(std::vector <C_FLOAT64> & target);
251 
252  /**
253  * Writes the state specified in the array source into the model.
254  * Length of the vector source must be mNumVariableMetabs.
255  * (Number of non-fixed metabolites in the model).
256  *
257  * @param source A vector reference with length mNumVariableMetabs,
258  * holding the state of the system to be set in the model
259  */
260  void setState(std::vector <C_FLOAT64> & source);
261 
262  /* STOCHASTIC STUFF ********************************************************/
263 
264  /**
265  * Find the reaction index and the reaction time of the stochastic (!)
266  * reaction with the lowest reaction time.
267  *
268  * @param ds A reference to a C_FLOAT64. The putative reaction time for
269  * the first stochastic reaction is written into this variable.
270  * @param rIndex A reference to a size_t. The index of the first
271  * stochastic reaction is written into this variable.
272  */
273  void getStochTimeAndIndex(C_FLOAT64 & ds, size_t & rIndex);
274 
275  /**
276  * Executes the specified reaction in the system once.
277  *
278  * @param rIndex A size_t specifying the index of the reaction, which
279  * will be fired.
280  */
281  void fireReaction(size_t rIndex);
282 
283  /**
284  * Updates the priority queue.
285  *
286  * @param rIndex A size_t giving the index of the fired reaction
287  * @param time A C_FLOAT64 holding the time taken by this reaction
288  */
289  void updatePriorityQueue(size_t rIndex, C_FLOAT64 time);
290 
291  /**
292  * Generates a putative reaction time for the given reaction.
293  *
294  * @param rIndex A size_t specifying the index of the reaction
295  * @return A C_FLOAT64 holding the calculated reaction time
296  */
297  C_FLOAT64 generateReactionTime(size_t rIndex);
298 
299  /**
300  * Calculates an amu value for a given reaction.
301  *
302  * @param rIndex A size_t specifying the reaction to be updated
303  */
304  void calculateAmu(size_t rIndex);
305 
306  /**
307  * Updates the putative reaction time of a stochastic reaction in the
308  * priority queue. The corresponding amu and amu_old must be set prior to
309  * the call of this method.
310  *
311  * @param rIndex A size_t specifying the index of the reaction
312  * @param time A C_FLOAT64 specifying the current time
313  */
314  void updateTauMu(size_t rIndex, C_FLOAT64 time);
315 
316  /* TESTING THE MODEL AND SETTING UP THINGS *********************************/
317 
318  /**
319  * Test the model if it is proper to perform stochastic simulations on.
320  * Several properties are tested (e.g. integer stoichometry, all
321  * reactions take place in one compartment only, irreversibility...).
322  *
323  * @param model The model to be checked
324  * @return 1, if hybrid simulation is possible; <0, if an error occured.
325  */
326  static C_INT32 checkModel(CModel * model);
327 
328  /**
329  * tests if the model contains a global value with an assignment rule that is
330  * used in calculations
331  */
332  static bool modelHasAssignments(const CModel* pModel);
333 
334  /**
335  * Sets up an internal representation of the balances for each reaction.
336  * This is done in order to be able to deal with fixed metabolites and
337  * to avoid a time consuming search for the indices of metabolites in the
338  * model.
339  */
340  void setupBalances();
341 
342  /**
343  * Sets up the dependency graph
344  */
345  void setupDependencyGraph();
346 
347  /**
348  * Creates for each metabolite a set of reaction indices. If the
349  * metabolite participates in a reaction as substrate or product this
350  * reaction is added to the corresponding set.
351  */
352  void setupMetab2React();
353 
354  /**
355  * Creates for each metabolite a set of reaction indices. If the
356  * metabolite participates in a reaction as substrate, product or
357  * modifier this reaction is added to the corresponding set.
358  */
360 
361  /**
362  * Creates for each metabolite a set of reaction indices. Each reaction
363  * is dependent on each metabolite resulting in a complete switch.
364  */
366 
367  /**
368  * Creates an initial partitioning of the system. Deterministic and
369  * stochastic reactions are determined. The array mStochReactions is
370  * initialized.
371  */
372  void setupPartition();
373 
374  /**
375  * Sets up the priority queue.
376  *
377  * @param startTime The time at which the simulation starts.
378  */
379  void setupPriorityQueue(C_FLOAT64 startTime = 0.0);
380 
381  /* HELPER METHODS **********************************************************/
382 
383  /**
384  * Updates the partitioning of the system depending on the particle
385  * numbers present.
386  */
387  void partitionSystem();
388 
389  /**
390  * Inserts a new deterministic reaction into the linked list in the
391  * vector mReactionFlags.
392  *
393  * @param rIndex A size_t giving the index of the reaction to be
394  * inserted into the list of deterministic reactions.
395  */
396  void insertDeterministicReaction(size_t rIndex);
397 
398  /**
399  * Removes a deterministic reaction from the linked list in the
400  * vector mReactionFlags.
401  *
402  * @param rIndex A size_t giving the index of the reaction to be
403  * removed from the list of deterministic reactions.
404  */
405  void removeDeterministicReaction(size_t rIndex);
406 
407  /**
408  * Gets the set of metabolites on which a given reaction depends.
409  *
410  * @param rIndex The index of the reaction being executed.
411  * @return The set of metabolites depended on.
412  */
413  std::set <std::string> *getDependsOn(size_t rIndex);
414 
415  /**
416  * Gets the set of metabolites which change number when a given
417  * reaction is executed.
418  *
419  * @param rIndex The index of the reaction being executed.
420  * @return The set of affected metabolites.
421  */
422  std::set <std::string> *getAffects(size_t rIndex);
423 
424  /**
425  * Gets the set of metabolites, which participate in the given
426  * reaction either as substrate, product or modifier.
427  *
428  * @param rIndex The index of the reaction being executed.
429  * @return The set of participating metabolites.
430  */
431  std::set <size_t> *getParticipatesIn(size_t rIndex);
432 
433  /**
434  * Prints out data on standard output.
435  */
436  void outputData(std::ostream & os, C_INT32 mode);
437 
438  /**
439  * Prints out various data on standard output for debugging purposes.
440  */
441  void outputDebug(std::ostream & os, size_t level);
442 
443  /* PRIVATE METHODS *********************************************************/
444 
445 private:
446  /**
447  * Intialize the method parameter
448  */
449  void initializeParameter();
450 
451 private:
452 
453  /* VARIABLES ***************************************************************/
454 
455 public:
456 
457 protected:
458 
459  /**
460  * Status of the metabolites. Particle number can be high or low.
461  */
462  enum metabStatus {LOW = 0, HIGH};
463 
464  /**
465  * Pointer to the model
466  */
468 
469  /**
470  * Dimension of the system. Total number of metabolites.
471  */
473 
474  /**
475  * index of the first metab in CState
476  */
478 
479  /**
480  * Max number of doSingleStep() per step()
481  */
482  unsigned C_INT32 mMaxSteps;
483 
485 
486  /**
487  * Specifies if the mRandomSeed should be used.
488  * otherwise a randomly chosen seed is used.
489  */
491 
492  /**
493  * The random seed to use.
494  */
496 
497  /**
498  * maximal increase of a particle number in one step.
499  */
500  size_t mMaxBalance;
501 
502  /**
503  * This is set to maxint - mMaxSteps*mMaxBalance
504  */
506 
507  /**
508  * indicates if the correction N^2 -> N*(N-1) should be performed
509  */
511 
512  /**
513  * A pointer to the reactions of the model.
514  */
516 
517  /**
518  * A pointer to the metabolites of the model.
519  */
521 
522  /**
523  * The stoichometry matrix of the model.
524  */
526 
527  /**
528  * Vectors to hold the system state and intermediate results
529  */
530  std::vector <C_FLOAT64> temp;
531  std::vector <C_FLOAT64> currentState;
532  std::vector <C_FLOAT64> k1;
533  std::vector <C_FLOAT64> k2;
534  std::vector <C_FLOAT64> k3;
535  std::vector <C_FLOAT64> k4;
536 
537  /**
538  * Vector to hold information about how many metabolites of a reaction
539  * have low particle numbers. If no metabolite of one reaction has
540  * low particle numbers this reaction will be simulated det.
541  */
542  std::vector <CHybridStochFlag> mReactionFlags;
544 
545  /**
546  * Vector holding information on the status of metabolites. They can
547  * have low or high particle numbers.
548  */
549  std::vector <metabStatus> mMetabFlags;
550 
551  /**
552  * Internal representation of the balances of each reaction. The index of
553  * each metabolite in the reaction is provided.
554  */
555  std::vector <std::vector <CHybridBalance> > mLocalBalances;
556 
557  /**
558  * Internal representation of the substrates of each reaction. The index
559  * of each substrate-metabolite is provided.
560  */
561  std::vector <std::vector <CHybridBalance> > mLocalSubstrates;
562 
563  /**
564  * Limit for particle numbers. If a particle number is < StochLimit the
565  * corresponding reactions must be simulated stochastically.
566  */
569 
570  /**
571  * The system gets repartitioned after this number of elementary steps.
572  */
574 
575  /**
576  * Number of elementary steps after the last partitioning.
577  */
579 
580  /**
581  * Stepsize for the rungeKutta steps of the numerical integrator.
582  */
584 
585  /**
586  * Vector of relations between metabolites to reactions.
587  */
588  std::vector <std::set <size_t> > mMetab2React;
589 
590  /**
591  * The propensities of the stochastic reactions.
592  */
593  std::vector <C_FLOAT64> mAmu;
594  std::vector <C_FLOAT64> mAmuOld;
595 
596  /**
597  * Set of the reactions, which must be updated.
598  */
599  std::set <size_t> mUpdateSet;
600 
601  /**
602  * The random number generator.
603  */
605 
606  /**
607  * The graph of reactions and their dependent reactions. When a reaction is
608  * executed, the propensities for each of its dependents must be updated.
609  */
611 
612  /**
613  * The set of putative stochastic (!) reactions and associated times at
614  * which each reaction occurs. This is represented as a priority queue,
615  * sorted by the reaction time. This heap changes dynamically as
616  * stochastic reactions become deterministic (delete this reaction from
617  * the queue) or vice versa (insert a new reaction into the queue).
618  */
620 
621  // DEPRECATED:
622  /**
623  * File output stream to write data.
624  */
625  std::ofstream mOutputFile;
626 
627  /**
628  * Output filename.
629  */
630  std::string mOutputFileName;
631 
632  /**
633  * Output counter.
634  */
636 
637  /**
638  * Indicates whether the model has global quantities with assignment rules.
639  * If it has, we will use a less efficient way to update the model
640  * state to handle this.
641  */
643 
644 private:
645 };
646 
648 
649 #endif // COPASI_CHybridMethod
CRandom * mpRandomGenerator
C_FLOAT64 mLowerStochLimit
unsigned C_INT32 mMaxSteps
CModel * mpModel
void outputData(std::ostream &os, C_INT32 mode)
void getState(std::vector< C_FLOAT64 > &target)
void getStochTimeAndIndex(C_FLOAT64 &ds, size_t &rIndex)
size_t mMaxIntBeforeStep
void initMethod(C_FLOAT64 time)
void calculateDerivative(std::vector< C_FLOAT64 > &deriv)
static C_INT32 checkModel(CModel *model)
std::vector< C_FLOAT64 > k3
std::vector< metabStatus > mMetabFlags
CHybridStochFlag * mFirstReactionFlag
void setState(std::vector< C_FLOAT64 > &source)
Definition: CState.h:305
size_t mMaxBalance
std::set< size_t > * getParticipatesIn(size_t rIndex)
void setupMetab2ReactComplete()
CMatrix< C_FLOAT64 > mStoi
#define C_INT32
Definition: copasi.h:90
CIndexedPriorityQueue mPQ
void setupMetab2ReactPlusModifier()
Definition: CMetab.h:178
virtual Status step(const double &deltaT)
size_t mFirstMetabIndex
CHybridStochFlag * mpPrev
Definition: CHybridMethod.h:91
void initializeParameter()
void updateTauMu(size_t rIndex, C_FLOAT64 time)
void integrateDeterministicPart(C_FLOAT64 ds)
void removeDeterministicReaction(size_t rIndex)
void rungeKutta(C_FLOAT64 dt)
const CCopasiVectorNS< CReaction > * mpReactions
std::vector< C_FLOAT64 > k4
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
void outputDebug(std::ostream &os, size_t level)
CDependencyGraph mDG
CHybridMethod(const CHybridMethod &src, const CCopasiContainer *pParent=NULL)
virtual bool elevateChildren()
C_FLOAT64 generateReactionTime(size_t rIndex)
std::string mOutputFileName
friend std::ostream & operator<<(std::ostream &os, const CHybridStochFlag &d)
std::set< std::string > * getAffects(size_t rIndex)
void setupMetab2React()
CHybridStochFlag * mpNext
Definition: CHybridMethod.h:92
CCopasiVector< CMetab > * mpMetabolites
void calculateAmu(size_t rIndex)
std::vector< C_FLOAT64 > mAmu
size_t mStepsAfterPartitionSystem
static CTrajectoryMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::deterministic)
C_FLOAT64 mStepsize
void updatePriorityQueue(size_t rIndex, C_FLOAT64 time)
unsigned C_INT32 mPartitioningInterval
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
#define C_FLOAT64
Definition: copasi.h:92
size_t mOutputCounter
std::vector< C_FLOAT64 > currentState
std::vector< C_FLOAT64 > k2
void insertDeterministicReaction(size_t rIndex)
std::vector< std::vector< CHybridBalance > > mLocalSubstrates
unsigned C_INT32 mRandomSeed
static CHybridMethod * createHybridMethod()
std::vector< std::vector< CHybridBalance > > mLocalBalances
std::vector< std::set< size_t > > mMetab2React
void integrateDeterministicPartEuler(C_FLOAT64 ds)
std::ofstream mOutputFile
Definition: CModel.h:50
virtual void start(const CState *initialState)
CMetab * mpMetabolite
C_FLOAT64 mUpperStochLimit
friend std::ostream & operator<<(std::ostream &os, const CHybridBalance &d)
std::vector< C_FLOAT64 > k1
C_INT32 mMultiplicity
std::set< size_t > mUpdateSet
std::vector< C_FLOAT64 > mAmuOld
std::vector< C_FLOAT64 > temp
void fireReaction(size_t rIndex)
std::vector< CHybridStochFlag > mReactionFlags
std::set< std::string > * getDependsOn(size_t rIndex)
size_t mNumVariableMetabs
static bool modelHasAssignments(const CModel *pModel)
void setupDependencyGraph()