COPASI API  4.16.103
CHybridMethodLSODA.h
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) 2006 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * CHybridMethodLSODA
17  *
18  * This class implements an hybrid algorithm to simulate a biochemical
19  * system over time.
20  *
21  * File name: CHybridMethodLSODA.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_CHybridMethodLSODA
43 #define COPASI_CHybridMethodLSODA
44 
45 /* INCLUDES ******************************************************************/
46 #include <set>
47 #include <vector>
48 #include <sstream>
49 #include <fstream>
50 
52 
59 
60 /* DEFINE ********************************************************************/
61 #define MAX_STEPS 1000000
62 #define INT_EPSILON 0.1
63 #define LOWER_STOCH_LIMIT 800 //800
64 #define UPPER_STOCH_LIMIT 1000 //1000
65 #define RUNGE_KUTTA_STEPSIZE 0.001
66 #define PARTITIONING_STEPSIZE 0.001
67 
68 #define PARTITIONING_INTERVAL 1
69 #define OUTPUT_COUNTER 100
70 //#define DEFAULT_OUTPUT_FILE "hybrid_lsoda.output"
71 #define SUBTYPE 1
72 #define USE_RANDOM_SEED false
73 #define RANDOM_SEED 1
74 
75 /* CLASSES *******************************************************************/
76 
77 class CMetab;
78 class CTrajectoryProblem;
79 class CState;
80 class CModel;
81 class CVersion;
82 class CReaction;
83 class CRandom;
85 class CDependencyGraph;
86 
87 /**
88  * Element of the array mReactionFlags. The deterministic reactions are
89  * organized in a linked list to be able to iterate over them quickly.
90  */
92 {
93 public:
94  size_t mIndex;
98 
99  // insert operator
100  friend std::ostream & operator<<(std::ostream & os, const CHybridLSODAStochFlag & d);
101 };
102 
103 /**
104  * Internal representation of the balances of each reaction. The index of
105  * each metabolite in the reaction is provided.
106  */
108 {
109 public:
110  size_t mIndex;
113 
114  // insert operator
115  friend std::ostream & operator<<(std::ostream & os, const CHybridLSODABalance & d);
116 };
117 
119 {
120 
121  friend CTrajectoryMethod *
123 
124  /* PUBLIC METHODS **********************************************************/
125 
126 public:
127  struct Data
128  {
131  };
132 
133 public:
134  /**
135  * Copy constructor
136  * @param const CHybridMethodLSODA & src
137  * @param const CCopasiContainer * pParent (default: NULL)
138  */
140  const CCopasiContainer * pParent = NULL);
141 
142  /**
143  * Destructor.
144  */
146 
147  /**
148  * This methods must be called to elevate subgroups to
149  * derived objects. The default implementation does nothing.
150  * @return bool success
151  */
152  virtual bool elevateChildren();
153 
154  /**
155  * Creates a HybridMethod adequate for the problem.
156  * (only CHybridNextReactionLSODAMethod so far)
157  */
158  static CHybridMethodLSODA *
160 
161  /**
162  * This instructs the method to calculate a time step of deltaT
163  * starting with the current state, i.e., the result of the previous
164  * step.
165  * The new state (after deltaT) is expected in the current state.
166  * The return value is the actual timestep taken.
167  * @param "const double &" deltaT
168  * @return Status status
169  */
170  virtual Status step(const double & deltaT);
171 
172  /**
173  * This instructs the method to prepare for integration
174  * starting with the initialState given.
175  * @param "const CState *" initialState
176  */
177  virtual void start(const CState * initialState);
178 
179  /**
180  * Check if the method is suitable for this problem
181  * @return bool suitability of the method
182  */
183  virtual bool isValidProblem(const CCopasiProblem * pProblem);
184 
185  /* PROTECTED METHODS *******************************************************/
186 
187 protected:
188  /**
189  * Default constructor.
190  * @param const CCopasiContainer * pParent (default: NULL)
191  */
192  CHybridMethodLSODA(const CCopasiContainer * pParent = NULL);
193 
194  /**
195  * Initializes the solver.
196  * @param time the current time
197  */
198  void initMethod(C_FLOAT64 time);
199 
200  /**
201  * Cleans up memory, etc.
202  */
203  void cleanup();
204 
205  /* VIRTUAL METHODS *********************************************************/
206 
207  /**
208  * Simulates the system over the next interval of time. The current time
209  * and the end time of the current step() are given as arguments.
210  *
211  * @param currentTime A C_FLOAT64 specifying the current time
212  * @param endTime A C_FLOAT64 specifying the end time of the step()
213  * @return A C_FLOAT giving the new time
214  */
215  virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime) = 0;
216 
217  /* DETERMINISTIC STUFF *****************************************************/
218 
219  /**
220  * Integrates the deterministic reactions of the system over the
221  * specified time interval.
222  *
223  * @param ds A C_FLOAT64 specifying the stepsize.
224  */
226 
227  /**
228  * Calculates the derivative of the system and writes it into the vector
229  * deriv. Length of deriv must be mNumVariableMetabs.
230  * CAUTION: Only deterministic reactions are taken into account. That is,
231  * this is only the derivative of the deterministic part of the system.
232  *
233  * @param deriv A vector reference of length mNumVariableMetabs, into
234  * which the derivative is written
235  */
236 
237  void calculateDerivative(std::vector <C_FLOAT64> & deriv);
238 
239  /**
240  * Gathers the state of the system into the array target. Later on CState
241  * should be used for this. Length of the array target must be
242  * mNumVariableMetabs.
243  *
244  * @param target A vector reference of length mNumVariableMetabs, into
245  * which the state of the system is written
246  */
247  void getState(std::vector <C_FLOAT64> & target);
248 
249  /**
250  * Writes the state specified in the array source into the model.
251  * Length of the vector source must be mNumVariableMetabs.
252  * (Number of non-fixed metabolites in the model).
253  *
254  * @param source A vector reference with length mNumVariableMetabs,
255  * holding the state of the system to be set in the model
256  */
257  void setState(std::vector <C_FLOAT64> & source);
258 
259  /* specific LSODA STAFF */
260 
261  /**
262  * Calculate the default absolute tolerance
263  * @param const CModel * pModel
264  * @return C_FLOAT64 defaultAtol
265  */
266  C_FLOAT64 getDefaultAtol(const CModel * pModel) const;
267 
268  static void EvalF(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot);
269 
270  /**
271  * This evaluates the derivatives for the complete model
272  */
273 
274  void evalF(const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot);
275 
276  /* STOCHASTIC STUFF ********************************************************/
277 
278  /**
279  * Find the reaction index and the reaction time of the stochastic (!)
280  * reaction with the lowest reaction time.
281  *
282  * @param ds A reference to a C_FLOAT64. The putative reaction time for
283  * the first stochastic reaction is written into this variable.
284  * @param rIndex A reference to a size_t. The index of the first
285  * stochastic reaction is written into this variable.
286  */
287  void getStochTimeAndIndex(C_FLOAT64 & ds, size_t & rIndex);
288 
289  /**
290  * Executes the specified reaction in the system once.
291  *
292  * @param rIndex A size_t specifying the index of the reaction, which
293  * will be fired.
294  */
295  void fireReaction(size_t rIndex);
296 
297  /**
298  * Updates the priority queue.
299  *
300  * @param rIndex A size_t giving the index of the fired reaction
301  * @param time A C_FLOAT64 holding the time taken by this reaction
302  */
303  void updatePriorityQueue(size_t rIndex, C_FLOAT64 time);
304 
305  /**
306  * Generates a putative reaction time for the given reaction.
307  *
308  * @param rIndex A size_t specifying the index of the reaction
309  * @return A C_FLOAT64 holding the calculated reaction time
310  */
311  C_FLOAT64 generateReactionTime(size_t rIndex);
312 
313  /**
314  * Calculates an amu value for a given reaction.
315  *
316  * @param rIndex A size_t specifying the reaction to be updated
317  */
318  void calculateAmu(size_t rIndex);
319 
320  /**
321  * Updates the putative reaction time of a stochastic reaction in the
322  * priority queue. The corresponding amu and amu_old must be set prior to
323  * the call of this method.
324  *
325  * @param rIndex A size_t specifying the index of the reaction
326  * @param time A C_FLOAT64 specifying the current time
327  */
328  void updateTauMu(size_t rIndex, C_FLOAT64 time);
329 
330  /* TESTING THE MODEL AND SETTING UP THINGS *********************************/
331 
332  /**
333  * Test the model if it is proper to perform stochastic simulations on.
334  * Several properties are tested (e.g. integer stoichometry, all
335  * reactions take place in one compartment only, irreversibility...).
336  *
337  * @param model The model to be checked
338  * @return 1, if hybrid simulation is possible; <0, if an error occured.
339  */
340  static C_INT32 checkModel(CModel * model);
341 
342  /**
343  * tests if the model contains a global value with an assignment rule that is
344  * used in calculations
345  */
346  static bool modelHasAssignments(const CModel* pModel);
347 
348  /**
349  * Sets up an internal representation of the balances for each reaction.
350  * This is done in order to be able to deal with fixed metabolites and
351  * to avoid a time consuming search for the indices of metabolites in the
352  * model.
353  */
354  void setupBalances();
355 
356  /**
357  * Sets up the dependency graph
358  */
359  void setupDependencyGraph();
360 
361  /**
362  * Creates for each metabolite a set of reaction indices. If the
363  * metabolite participates in a reaction as substrate or product this
364  * reaction is added to the corresponding set.
365  */
366  void setupMetab2React();
367 
368  /**
369  * Creates for each metabolite a set of reaction indices. If the
370  * metabolite participates in a reaction as substrate, product or
371  * modifier this reaction is added to the corresponding set.
372  */
374 
375  /**
376  * Creates for each metabolite a set of reaction indices. Each reaction
377  * is dependent on each metabolite resulting in a complete switch.
378  */
380 
381  /**
382  * Creates an initial partitioning of the system. Deterministic and
383  * stochastic reactions are determined. The array mStochReactions is
384  * initialized.
385  */
386  void setupPartition();
387 
388  /**
389  * Sets up the priority queue.
390  *
391  * @param startTime The time at which the simulation starts.
392  */
393  void setupPriorityQueue(C_FLOAT64 startTime = 0.0);
394 
395  /* HELPER METHODS **********************************************************/
396 
397  /**
398  * Updates the partitioning of the system depending on the particle
399  * numbers present.
400  */
401  void partitionSystem();
402 
403  /**
404  * Inserts a new deterministic reaction into the linked list in the
405  * vector mReactionFlags.
406  *
407  * @param rIndex A size_t giving the index of the reaction to be
408  * inserted into the list of deterministic reactions.
409  */
410  void insertDeterministicReaction(size_t rIndex);
411 
412  /**
413  * Removes a deterministic reaction from the linked list in the
414  * vector mReactionFlags.
415  *
416  * @param rIndex A size_t giving the index of the reaction to be
417  * removed from the list of deterministic reactions.
418  */
419  void removeDeterministicReaction(size_t rIndex);
420 
421  /**
422  * Gets the set of metabolites on which a given reaction depends.
423  *
424  * @param rIndex The index of the reaction being executed.
425  * @return The set of metabolites depended on.
426  */
427  std::set <std::string> *getDependsOn(size_t rIndex);
428 
429  /**
430  * Gets the set of metabolites which change number when a given
431  * reaction is executed.
432  *
433  * @param rIndex The index of the reaction being executed.
434  * @return The set of affected metabolites.
435  */
436  std::set <std::string> *getAffects(size_t rIndex);
437 
438  /**
439  * Gets the set of metabolites, which participate in the given
440  * reaction either as substrate, product or modifier.
441  *
442  * @param rIndex The index of the reaction being executed.
443  * @return The set of participating metabolites.
444  */
445  std::set <size_t> *getParticipatesIn(size_t rIndex);
446 
447  /**
448  * Prints out data on standard output.
449  */
450  void outputData(std::ostream & os, C_INT32 mode);
451 
452  /**
453  * Prints out various data on standard output for debugging purposes.
454  */
455  void outputDebug(std::ostream & os, size_t level);
456 
457  /* PRIVATE METHODS *********************************************************/
458 
459 private:
460  /**
461  * Intialize the method parameter
462  */
463  void initializeParameter();
464 
465 private:
466 
467  /* VARIABLES ***************************************************************/
468 
469 public:
470 
471 protected:
472 
473  /**
474  * Status of the metabolites. Particle number can be high or low.
475  */
476  enum metabStatus {LOW = 0, HIGH};
477 
478  /**
479  * Pointer to the model
480  */
482 
483  /**
484  * Dimension of the system. Total number of metabolites.
485  */
487 
488  /**
489  * index of the first metab in CState
490  */
492 
493  /**
494  * Max number of doSingleStep() per step()
495  */
496  unsigned C_INT32 mMaxSteps;
497 
499 
500  /**
501  * Specifies if the mRandomSeed should be used.
502  * otherwise a randomly chosen seed is used.
503  */
505 
506  /**
507  * The random seed to use.
508  */
510 
511  /**
512  * maximal increase of a particle number in one step.
513  */
515 
516  /**
517  * This is set to maxint - mMaxSteps*mMaxBalance
518  */
520 
521  /**
522  * indicates if the correction N^2 -> N*(N-1) should be performed
523  */
525 
526  /**
527  * A pointer to the reactions of the model.
528  */
530 
531  /**
532  * A pointer to the metabolites of the model.
533  */
535 
536  /**
537  * The stoichometry matrix of the model.
538  */
540 
541  /**
542  * Vectors to hold the system state and intermediate results
543  */
544  std::vector <C_FLOAT64> temp;
545  std::vector <C_FLOAT64> currentState;
546 
547  /**
548  * Vector to hold information about how many metabolites of a reaction
549  * have low particle numbers. If no metabolite of one reaction has
550  * low particle numbers this reaction will be simulated det.
551  */
552  std::vector <CHybridLSODAStochFlag> mReactionFlags;
554 
555  /**
556  * Vector holding information on the status of metabolites. They can
557  * have low or high particle numbers.
558  */
559  std::vector <metabStatus> mMetabFlags;
560 
561  /**
562  * Internal representation of the balances of each reaction. The index of
563  * each metabolite in the reaction is provided.
564  */
565  std::vector <std::vector <CHybridLSODABalance> > mLocalBalances;
566 
567  /**
568  * Internal representation of the substrates of each reaction. The index
569  * of each substrate-metabolite is provided.
570  */
571  std::vector <std::vector <CHybridLSODABalance> > mLocalSubstrates;
572 
573  /**
574  * Limit for particle numbers. If a particle number is < StochLimit the
575  * corresponding reactions must be simulated stochastically.
576  */
579 
580  /**
581  * The system gets repartitioned after this number of elementary steps.
582  */
584 
585  /**
586  * Number of elementary steps after the last partitioning.
587  */
589 
590  /**
591  * Partitioning stepsize
592  */
594 
595  /**
596  * Vector of relations between metabolites to reactions.
597  */
598  std::vector <std::set <size_t> > mMetab2React;
599 
600  /**
601  * The propensities of the stochastic reactions.
602  */
603  std::vector <C_FLOAT64> mAmu;
604  std::vector <C_FLOAT64> mAmuOld;
605 
606  /**
607  * Set of the reactions, which must be updated.
608  */
609  std::set <size_t> mUpdateSet;
610 
611  /**
612  * The random number generator.
613  */
615 
616  /**
617  * The graph of reactions and their dependent reactions. When a reaction is
618  * executed, the propensities for each of its dependents must be updated.
619  */
621 
622  /**
623  * The set of putative stochastic (!) reactions and associated times at
624  * which each reaction occurs. This is represented as a priority queue,
625  * sorted by the reaction time. This heap changes dynamically as
626  * stochastic reactions become deterministic (delete this reaction from
627  * the queue) or vice versa (insert a new reaction into the queue).
628  */
630 
631  /* specific LSODA STAFF */
632 
633  /**
634  * this flag indicates whether the next LSODA call needs to reinitialize the
635  * integrater (if the partitioning has changed or a stochastic reaction was fired)
636  */
638 
639  /**
640  * A pointer to the current state in complete model view.
641  */
643 
644  /**
645  * mData.dim is the dimension of the ODE system.
646  * mData.pMethod contains CLsodaMethod * this to be used in the static method EvalF
647  */
649 
650  /**
651  * Pointer to the array with left hand side values.
652  */
654 
655  /**
656  * Vector containig the derivatives after calling eval
657  */
659 
660  /**
661  * Current time.
662  */
664 
665  /**
666  * Requested end time.
667  */
669 
670  /**
671  * LSODA state.
672  */
674 
675  /**
676  *
677  */
679 
680  /**
681  * Relative tolerance.
682  */
684 
685  /**
686  *
687  */
689 
690  /**
691  * Absolute tolerance.
692  */
694 
695  /**
696  *
697  */
698  std::ostringstream mErrorMsg;
699 
705 
706  // DEPRECATED:
707  /**
708  * File output stream to write data.
709  */
710  std::ofstream mOutputFile;
711 
712  /**
713  * Output filename.
714  */
715  std::string mOutputFileName;
716 
717  /**
718  * Output counter.
719  */
721 
722  /**
723  * Indicates whether the model has global quantities with assignment rules.
724  * If it has, we will use a less efficient way to update the model
725  * state to handle this.
726  */
728 };
729 
731 
732 #endif // COPASI_CHybridMethodLSODA
C_FLOAT64 getDefaultAtol(const CModel *pModel) const
#define C_INT
Definition: copasi.h:115
void setState(std::vector< C_FLOAT64 > &source)
CVector< C_FLOAT64 > mDWork
void outputData(std::ostream &os, C_INT32 mode)
void getStochTimeAndIndex(C_FLOAT64 &ds, size_t &rIndex)
unsigned C_INT32 mPartitioningInterval
std::vector< std::set< size_t > > mMetab2React
std::vector< C_FLOAT64 > mAmu
unsigned C_INT32 mMaxSteps
CIndexedPriorityQueue mPQ
static C_INT32 checkModel(CModel *model)
void updateTauMu(size_t rIndex, C_FLOAT64 time)
unsigned C_INT32 mStepsAfterPartitionSystem
virtual bool elevateChildren()
virtual bool isValidProblem(const CCopasiProblem *pProblem)
static bool modelHasAssignments(const CModel *pModel)
Definition: CState.h:305
C_FLOAT64 generateReactionTime(size_t rIndex)
CHybridLSODAStochFlag * mpNext
#define C_INT32
Definition: copasi.h:90
std::ostringstream mErrorMsg
Definition: CMetab.h:178
CVector< C_INT > mIWork
void updatePriorityQueue(size_t rIndex, C_FLOAT64 time)
unsigned C_INT32 mRandomSeed
friend std::ostream & operator<<(std::ostream &os, const CHybridLSODABalance &d)
CCopasiVector< CMetab > * mpMetabolites
CHybridMethodLSODA * pMethod
CMatrix< C_FLOAT64 > mStoi
static CHybridMethodLSODA * createHybridMethodLSODA()
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
virtual Status step(const double &deltaT)
void integrateDeterministicPart(C_FLOAT64 ds)
Definition: CLSODA.h:24
CVector< C_FLOAT64 > mYdot
std::vector< std::vector< CHybridLSODABalance > > mLocalBalances
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
std::vector< C_FLOAT64 > temp
std::ofstream mOutputFile
void getState(std::vector< C_FLOAT64 > &target)
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
std::vector< C_FLOAT64 > currentState
CHybridMethodLSODA(const CHybridMethodLSODA &src, const CCopasiContainer *pParent=NULL)
friend std::ostream & operator<<(std::ostream &os, const CHybridLSODAStochFlag &d)
std::set< size_t > * getParticipatesIn(size_t rIndex)
std::vector< CHybridLSODAStochFlag > mReactionFlags
void removeDeterministicReaction(size_t rIndex)
static CTrajectoryMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::deterministic)
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
std::set< std::string > * getDependsOn(size_t rIndex)
virtual void start(const CState *initialState)
void initMethod(C_FLOAT64 time)
#define C_FLOAT64
Definition: copasi.h:92
void outputDebug(std::ostream &os, size_t level)
CHybridLSODAStochFlag * mpPrev
CDependencyGraph mDG
const CCopasiVectorNS< CReaction > * mpReactions
std::set< std::string > * getAffects(size_t rIndex)
Definition: CModel.h:50
std::vector< metabStatus > mMetabFlags
void insertDeterministicReaction(size_t rIndex)
std::vector< C_FLOAT64 > mAmuOld
std::vector< std::vector< CHybridLSODABalance > > mLocalSubstrates
void calculateDerivative(std::vector< C_FLOAT64 > &deriv)
CHybridLSODAStochFlag * mFirstReactionFlag
void fireReaction(size_t rIndex)
std::set< size_t > mUpdateSet
void calculateAmu(size_t rIndex)