COPASI API  4.16.103
CHybridMethodODE45.h
Go to the documentation of this file.
1 // Copyright (C) 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 /**
7  * CHybridMethodODE45
8  *
9  * This class implements an hybrid algorithm to simulate a biochemical
10  * system over time.
11  *
12  * File name: CHybridMethodODE45.h
13  * Author: Shuo Wang
14  * Email: shuowang.learner@gmail.com
15  *
16  * Last change: 26, Aug 2013
17  *
18  */
19 /**
20  * Partition the system into a deterministic part and a stochastic part.
21  * That is, every reaction is either classified deterministic or
22  * stochastic. Deterministic reactions involve only those metabolites (on
23  * substrate and product side), which have a high particle number.
24  * Therefore it is legal to integrate this part of the system with e.g. a
25  * numerical integrator. The concentration and relative particle number
26  * change should be low enough, so that the probabilities of all the
27  * reactions in the system show only little changes. The stochastic
28  * reactions must be simulated with an exact stochastic method (e.g. next
29  * reaction method (Gibson)), because their firing changes the reaction
30  * probabilities in the system significantly.
31  */
32 
33 #ifndef COPASI_CHybridMethodODE45
34 #define COPASI_CHybridMethodODE45
35 
36 /* INCLUDES ******************************************************************/
37 #include <set>
38 #include <vector>
39 #include <sstream>
40 #include <fstream>
41 #include <cstdlib>
42 #include <new>
43 
45 
51 #include "CInterpolation.h"
52 
53 /* DEFINE ********************************************************************/
54 #define MAX_STEPS 1000000
55 #define INT_EPSILON 0.1
56 //#define LOWER_STOCH_LIMIT 800 //800
57 //#define UPPER_STOCH_LIMIT 1000 //1000
58 //#define RUNGE_KUTTA_STEPSIZE 0.001
59 
60 //#define PARTITIONING_STEPSIZE 0.001
61 //#define PARTITIONING_INTERVAL 1
62 //#define OUTPUT_COUNTER 100
63 
64 //#define DEFAULT_OUTPUT_FILE "hybrid_lsoda.output"
65 #define SUBTYPE 1
66 #define USE_RANDOM_SEED false
67 #define RANDOM_SEED 1
68 
69 //Simulation Part
70 #define SLOW 0
71 #define FAST 1
72 
73 //Method Used
74 #define STOCHASTIC 0
75 #define DETERMINISTIC 1
76 #define HYBRID 2
77 
78 
79 //Interpolation Part
80 #define INTERP_RECORD_NUM 6
81 #define HAS_ERR -1
82 #define REACH_END_TIME 0
83 #define CONTINUE 1
84 #define HAS_EVENT 2
85 #define NEW_STEP 3
86 
87 /* Function Pointer **********************************************************/
88 typedef void (*pEvalF)(const C_INT*, const double*, const double*, double*);
89 
90 /* CLASSES *******************************************************************/
91 
92 class CMetab;
93 class CTrajectoryProblem;
94 class CState;
95 class CModel;
96 class CVersion;
97 class CReaction;
98 class CRandom;
100 class CDependencyGraph;
101 
102 class CStateRecord;
103 class CInterpolation;
104 
105 /**
106  * Internal representation of the balances of each reaction.
107  * The index of each metabolite in the reaction is provided.
108  */
110 {
111 public:
112  size_t mIndex;
115 
116  // insert operator
117  //friend std::ostream & operator<<(std::ostream & os, const CHybridODE45Balance & d);
118 };
119 
120 /**
121  * A class to record whether a metab is slow or fast
122  * @param fastReactions is applied to store which reactions this
123  * metab participates
124  * @param flag, if set is empty -> false, else, -> true
125  */
127 {
128 public:
129  std::set<size_t> mFastReactions;
130  size_t mFlag;
131 };
132 
134 {
135 
136  friend CTrajectoryMethod *
138 
139 public:
140  struct Data
141  {
144  };
145 
146  //================Function for Class================
147 protected:
148  /**
149  * Default Constructor
150  */
151  CHybridMethodODE45(const CCopasiContainer * pParent = NULL);
152 
153 public:
154  /**
155  * Copy constructor
156  * @param const CHybridMethodODE45 & src
157  * @param const CCopasiContainer * pParent (default: NULL)
158  */
160  const CCopasiContainer * pParent = NULL);
161  /**
162  * Destructor.
163  */
165 
166  //================Function for System================
167 public:
168 
169  /**
170  * Check if the method is suitable for this problem
171  * @return bool suitability of the method
172  */
173  virtual bool isValidProblem(const CCopasiProblem * pProblem);
174 
175  /**
176  * This methods must be called to elevate subgroups to
177  * derived objects. The default implementation does nothing.
178  * @return bool success
179  */
180  virtual bool elevateChildren();
181 
182  /**
183  * This instructs the method to prepare for integration
184  * starting with the initialState given.
185  * @param "const CState *" initialState
186  */
187  virtual void start(const CState * initialState);
188 
189 protected:
190  /**
191  * Initializes the solver.
192  * @param time the current time
193  */
194  void initMethod(C_FLOAT64 time);
195 
196  /**
197  * Cleans up memory, etc.
198  */
199  void cleanup();
200 
201 private:
202  /**
203  * Intialize the method parameter
204  */
205  void initializeParameter();
206 
207  //================Function for Model================
208 public:
209 
210 protected:
211  /**
212  * setup mMetabFlags
213  */
214  void setupMetabFlags();
215 
216  /**
217  * setup mReactionFlags
218  */
219  void setupReactionFlags();
220 
221  /**
222  * setup mMethod
223  */
224  void setupMethod();
225 
226  /**
227  *
228  */
229  void setupCalculateSet();
230 
231  /**
232  * Sets up an internal representation of the balances for each reaction.
233  * This is done in order to be able to deal with fixed metabolites and
234  * to avoid a time consuming search for the indices of metabolites in the
235  * model.
236  */
237  void setupBalances();
238 
239  /**
240  * Sets up an internal representation of the balances for each reaction.
241  * This is done in order to be able to deal with fixed metabolites and
242  * to avoid a time consuming search for the indices of metabolites in the
243  * model.
244  */
245  void setupMetab2React();
246 
247  void setupReactAffect();
248 
249  //================Function for ODE45================
250 public:
251 
252 protected:
253  /**
254  * Calculate the default absolute tolerance
255  * @param const CModel * pModel
256  * @return C_FLOAT64 defaultAtol
257  */
258  C_FLOAT64 getDefaultAtol(const CModel * pModel) const;
259 
260  /**
261  * Integrates the deterministic reactions of the system over the
262  * specified time interval.
263  *
264  * @param ds A C_FLOAT64 specifying the stepsize.
265  */
267 
268  /**
269  * Dummy Function for calculating derivative of ODE systems
270  */
271  static void EvalF(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot);
272 
273  /**
274  * This evaluates the derivatives for the complete model
275  */
276  void evalF(const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot);
277 
278  //================Function for Simulation================
279 public:
280  /**
281  * This instructs the method to calculate a time step of deltaT
282  * starting with the current state, i.e., the result of the previous
283  * step.
284  * The new state (after deltaT) is expected in the current state.
285  * The return value is the actual timestep taken.
286  * @param "const double &" deltaT
287  * @return Status status
288  */
289  virtual Status step(const double & deltaT);
290 
291 protected:
292  /**
293  * Simulates the system over the next interval of time. The current time
294  * and the end time of the current step() are given as arguments.
295  *
296  * @param currentTime A C_FLOAT64 specifying the current time
297  * @param endTime A C_FLOAT64 specifying the end time of the step()
298  * @return A C_FLOAT giving the new time
299  */
300  C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime);
301 
302  /**
303  * Calculates an amu value for a given reaction.
304  *
305  * @param rIndex A size_t specifying the reaction to be updated
306  */
307  void calculateAmu(size_t rIndex);
308 
309  /**
310  * Do inverse interpolation to find the state when a slow reaction
311  * is fired.
312  */
313  void doInverseInterpolation();
314 
315  /**
316  * Fire slow reaction and update populations and propensities
317  * when Hybrid Method is used
318  */
320 
321  //================Function for Stoichastic Part================
322 protected:
323 
324  /**
325  * Sets up the dependency graph
326  */
327  void setupDependencyGraph();
328 
329  /**
330  * Sets up the priority queue.
331  *
332  * @param startTime The time at which the simulation starts.
333  */
334  void setupPriorityQueue(C_FLOAT64 startTime = 0.0);
335 
336  /**
337  * Updates the priority queue.
338  *
339  * @param rIndex A size_t giving the index of the fired reaction
340  * @param time A C_FLOAT64 holding the time taken by this reaction
341  */
342  void updatePriorityQueue(size_t rIndex, C_FLOAT64 time);
343 
344  /**
345  * Executes the specified reaction in the system once.
346  *
347  * @param rIndex A size_t specifying the index of the reaction, which
348  * will be fired.
349  */
350  void fireReaction(size_t rIndex);
351 
352  /**
353  * Function return a reaction index which is a slow reaction firing
354  * at the event time.
355  */
356  size_t getReactionIndex4Hybrid();
357 
358  /**
359  * Find the reaction index and the reaction time of the stochastic (!)
360  * reaction with the lowest reaction time.
361  *
362  * @param ds A reference to a C_FLOAT64. The putative reaction time for
363  * the first stochastic reaction is written into this variable.
364  * @param rIndex A reference to a size_t. The index of the first
365  * stochastic reaction is written into this variable.
366  */
367  void getStochTimeAndIndex(C_FLOAT64 & ds, size_t & rIndex);
368 
369  /**
370  * Generates a putative reaction time for the given reaction.
371  *
372  * @param rIndex A size_t specifying the index of the reaction
373  * @return A C_FLOAT64 holding the calculated reaction time
374  */
375  C_FLOAT64 generateReactionTime(size_t rIndex);
376 
377  /**
378  * Updates the putative reaction time of a stochastic reaction in the
379  * priority queue. The corresponding amu and amu_old must be set prior to
380  * the call of this method.
381  *
382  * @param rIndex A size_t specifying the index of the reaction
383  * @param time A C_FLOAT64 specifying the current time
384  */
385  void updateTauMu(size_t rIndex, C_FLOAT64 time);
386 
387 private:
388  /**
389  * Gets the set of metabolites on which a given reaction depends.
390  *
391  * @param rIndex The index of the reaction being executed.
392  * @return The set of metabolites depended on.
393  */
394  std::set <std::string> *getDependsOn(size_t rIndex);
395 
396  /**
397  * Gets the set of metabolites which change number when a given
398  * reaction is executed.
399  *
400  * @param rIndex The index of the reaction being executed.
401  * @return The set of affected metabolites.
402  */
403  std::set <std::string> *getAffects(size_t rIndex);
404 
405  //================Function for C Code from f2c================
406 private:
407  C_INT rkf45_(pEvalF f, const C_INT *neqn, double *y, double *t,
408  double *tout, double *relerr, double *abserr,
409  C_INT *iflag, double *work, C_INT *iwork,
410  double *yrcd);
411 
412  C_INT rkfs_(pEvalF f, const C_INT *neqn, double *y, double *
413  t, double *tout, double *relerr,
414  double *abserr, C_INT *iflag, double *yp,
415  double *h__, double *f1, double *f2,
416  double *f3, double *f4, double *f5,
417  double *savre, double *savae, C_INT *nfe,
418  C_INT *kop, C_INT *init, C_INT *jflag,
419  C_INT *kflag, double *yrcd);
420 
421  C_INT fehl_(pEvalF f, const C_INT *neqn, double *y, double *t,
422  double *h__, double *yp, double *f1,
423  double *f2, double *f3, double *f4,
424  double *f5, double *s, double *yrcd);
425 
426  //================Help Functions================
427 protected:
428  /**
429  * Test the model if it is proper to perform stochastic simulations on.
430  * Several properties are tested (e.g. integer stoichometry, all reactions
431  * take place in one compartment only, irreversibility...).
432  *
433  * @return 0, if everything is ok; <0, if an error occured.
434  */
435  C_INT32 checkModel(CModel * model);
436 
437  /**
438  * Prints out data on standard output.
439  */
440  //void outputData(std::ostream & os, C_INT32 mode);
441  void outputData();
442  void outputState(const CState * pS);
443 
444  /**
445  * Prints out various data on standard output for debugging purposes.
446  */
447  void outputDebug(std::ostream & os, size_t level);
448 
449  /**
450  * tests if the model contains a global value with an assignment rule that is
451  * used in calculations
452  */
453  static bool modelHasAssignments(const CModel* pModel);
454 
455 //Attributes:
456  //================Model Related================
457 
458  //~~~~~~~~Model Describtion~~~~~~~~
459  /**
460  * Pointer to the model
461  */
463 
464  /**
465  * The stoichometry matrix of the model.
466  */
468 
469  /**
470  * indicates if the correction N^2 -> N*(N-1) should be performed
471  */
473 
474  /**
475  *
476  */
478 
479  //~~~~~~~~Metabs Related~~~~~~~~
480  /**
481  * Dimension of the system. Total number of metabolites.
482  */
484 
485  /**
486  * A pointer to the metabolites of the model.
487  */
489 
490  /**
491  * index of the first metab in CState
492  */
494 
495  /**
496  * Vector holding information on the status of metabolites. They can
497  * be FAST or SLOW.
498  */
499  std::vector<CHybridODE45MetabFlag> mMetabFlags;
500 
501  //~~~~~~~~Reaction Related~~~~~~~~
502  /**
503  * Number of Reactions
504  */
506 
507  /**
508  * Fast reactions are set FAST and
509  * slow ones are SLOW
510  */
512 
513  /**
514  * A pointer to the reactions of the model.
515  */
517 
518  /**
519  * Internal representation of the balances of each reaction. The index of
520  * each metabolite in the reaction is provided.
521  */
522  std::vector <std::vector <CHybridODE45Balance> > mLocalBalances;
523 
524  /**
525  * Internal representation of the substrates of each reaction. The index
526  * of each substrate-metabolite is provided.
527  */
528  std::vector <std::vector <CHybridODE45Balance> > mLocalSubstrates;
529 
530  /**
531  * Vector of relations between metabolites with reactions.
532  */
533  std::vector <std::set <size_t> > mMetab2React;
534 
535  /**
536  * Vector of sets that the indeces of propencities which should be
537  * updated after one reaction has been fired
538  */
539  std::vector <std::set<size_t> > mReactAffect;
540 
541  /**
542  * Bool value
543  */
545 
546  /**
547  *
548  */
550 
551  /**
552  * An Integer showes the method now CHybridMethod used
553  * 0 == Stochastic
554  * 1 == Deterministic
555  * 2 == Hybrid
556  */
557  size_t mMethod;
558 
559  //================Attributes for State================
560  /**
561  * A pointer to the current state in complete model view.
562  */
564 
565  /**
566  * Vectors to hold the system state and intermediate results
567  */
569 
570  //=================Attributes for ODE45 Solver================
571 
572  /**
573  * Max number of doSingleStep() per step()
574  */
575  size_t mMaxSteps;
577 
578  /** Is this useful in my method?
579  * maximal increase of a particle number in one step.
580  */
581  size_t mMaxBalance;
582 
583  /**
584  * mData.dim is the dimension of the ODE system.
585  * mData.pMethod contains CLsodaMethod * this to be used
586  * in the static method EvalF
587  */
589 
590  /**
591  * Pointer to the array with left hand side values.
592  */
594 
595  /**
596  * Vector containig the derivatives after calling eval
597  */
599 
600  /**
601  * Current time.
602  */
604 
605  /**
606  * Requested end time.
607  */
609 
610  /**
611  * ODE45 state, corresponding to iflag in rkf45, an ODE45
612  * solver argument
613  */
615 
616  /**
617  * state of ODE45, indicating what to do next in the step part
618  * -1, error emerge
619  * 0, finish all integration
620  * 1, finish one step integration
621  * 2, has event
622  */
624 
625  /**
626  * Relative tolerance.
627  */
629 
630  /**
631  *
632  */
634 
635  /**
636  * Absolute tolerance.
637  */
639 
642 
643  /**
644  * A boolean variable to show whether mStateRecord is
645  * used or not. If t is too close to tout, rkfs45() won't
646  * go through into fehl() and uses Forward Euler method
647  */
649 
650  /**
651  * The propensities of the stochastic reactions.
652  */
655 
656  /**
657  * Set of the reactions, which must be updated.
658  * Reactions have at lease one fast metab as subtract with
659  * non-zero balance
660  */
661  std::set <size_t> mCalculateSet;
662 
663  /**
664  * Set of the reactions, which must update after one
665  * slow reaction fires
666  */
667  std::set <size_t> mUpdateSet;
668 
669  //================Attributes for Interpolation================
670  /**
671  * Record of y and time of the previous step
672  */
675 
676  /**
677  * A Pointer to Object of CInterpolation
678  */
680 
681  /**
682  * Pointer to the State at Event Happens
683  */
685 
686  /**
687  * A vector to record middle state for interpolation
688  */
690 
691  //================Stochastic Related================
692  /**
693  * The random number generator.
694  */
696 
697  /**
698  * Specifies if the mRandomSeed should be used.
699  * otherwise a randomly chosen seed is used.
700  */
702 
703  /**
704  * The random seed to use.
705  */
707 
708  /**
709  * The graph of reactions and their dependent reactions.
710  * When a reaction is executed, the propensities for each of
711  * its dependents must be updated.
712  */
714 
715  /**
716  * The set of putative stochastic (!) reactions and associated times at
717  * which each reaction occurs. This is represented as a priority queue,
718  * sorted by the reaction time. This heap changes dynamically as
719  * stochastic reactions become deterministic (delete this reaction from
720  * the queue) or vice versa (insert a new reaction into the queue).
721  */
723 
724  //========System Related========
725  /**
726  * File output stream to write data.
727  */
728  std::ofstream mOutputFile;
729 
730  /**
731  * Output filename.
732  */
733  std::string mOutputFileName;
734 
735  /**
736  * Output counter.
737  */
739 
740  /**
741  * Indicates whether the model has global quantities
742  * with assignment rules.
743  * If it has, we will use a less efficient way to update the model
744  * state to handle this.
745  */
747 
748  /**
749  *
750  */
751  std::ostringstream mErrorMsg;
752 };
753 
754 #endif // COPASI_CHybridMethodODE45
std::set< size_t > mUpdateSet
#define C_INT
Definition: copasi.h:115
CVector< C_FLOAT64 > temp
const CCopasiVectorNS< CReaction > * mpReactions
C_INT fehl_(pEvalF f, const C_INT *neqn, double *y, double *t, double *h__, double *yp, double *f1, double *f2, double *f3, double *f4, double *f5, double *s, double *yrcd)
virtual void start(const CState *initialState)
CDependencyGraph mDG
std::vector< std::set< size_t > > mMetab2React
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
void getStochTimeAndIndex(C_FLOAT64 &ds, size_t &rIndex)
void outputState(const CState *pS)
C_INT rkfs_(pEvalF f, const C_INT *neqn, double *y, double *t, double *tout, double *relerr, double *abserr, C_INT *iflag, double *yp, double *h__, double *f1, double *f2, double *f3, double *f4, double *f5, double *savre, double *savae, C_INT *nfe, C_INT *kop, C_INT *init, C_INT *jflag, C_INT *kflag, double *yrcd)
Definition: CState.h:305
CVector< C_FLOAT64 > mAmuOld
C_FLOAT64 generateReactionTime(size_t rIndex)
CMatrix< C_FLOAT64 > mStoi
unsigned C_INT32 mRandomSeed
C_INT32 checkModel(CModel *model)
C_INT rkf45_(pEvalF f, const C_INT *neqn, double *y, double *t, double *tout, double *relerr, double *abserr, C_INT *iflag, double *work, C_INT *iwork, double *yrcd)
std::vector< CHybridODE45MetabFlag > mMetabFlags
#define C_INT32
Definition: copasi.h:90
Definition: CMetab.h:178
C_FLOAT64 getDefaultAtol(const CModel *pModel) const
void updatePriorityQueue(size_t rIndex, C_FLOAT64 time)
CVector< C_FLOAT64 > mYdot
virtual bool elevateChildren()
std::vector< std::vector< CHybridODE45Balance > > mLocalSubstrates
void fireReaction(size_t rIndex)
std::set< std::string > * getDependsOn(size_t rIndex)
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
CVector< C_FLOAT64 > mAmu
std::ofstream mOutputFile
std::set< std::string > * getAffects(size_t rIndex)
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
CInterpolation * mpInterpolation
CIndexedPriorityQueue mPQ
CVector< C_INT > mIWork
virtual Status step(const double &deltaT)
static CTrajectoryMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::deterministic)
CHybridMethodODE45 * pMethod
CVector< C_FLOAT64 > mStateRecord
std::set< size_t > mFastReactions
static bool modelHasAssignments(const CModel *pModel)
CVector< C_FLOAT64 > mDWork
void updateTauMu(size_t rIndex, C_FLOAT64 time)
CCopasiVector< CMetab > * mpMetabolites
#define C_FLOAT64
Definition: copasi.h:92
void(* pEvalF)(const C_INT *, const double *, const double *, double *)
std::vector< std::set< size_t > > mReactAffect
C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)
() void(yyvaluep))
std::ostringstream mErrorMsg
Definition: CModel.h:50
void integrateDeterministicPart(C_FLOAT64 ds)
CHybridMethodODE45(const CCopasiContainer *pParent=NULL)
CVector< size_t > mReactionFlags
void calculateAmu(size_t rIndex)
void initMethod(C_FLOAT64 time)
std::set< size_t > mCalculateSet
CStateRecord * mpEventState
void outputDebug(std::ostream &os, size_t level)
virtual bool isValidProblem(const CCopasiProblem *pProblem)