COPASI API  4.16.103
CHybridMethodODE45.cpp
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 for the simulation of a
10  * biochemical system over time.
11  *
12  * File name: CHybridMethodODE45.cpp
13  * Author: Shuo Wang
14  * Email: shuowang.learner@gmail.com
15  *
16  * Last change: 26, Aug 2013
17  *
18  */
19 
20 /* DEFINE ********************************************************************/
21 
22 #ifdef WIN32
23 #if _MSC_VER < 1600
24 #define min _cpp_min
25 #define max _cpp_max
26 #endif // _MSC_VER
27 #endif // WIN32
28 
29 #include <stdio.h>
30 #include <limits.h>
31 #include <iterator>
32 #include <iostream>
33 #include <algorithm>
34 #include <new>
35 #include <cmath>
36 
37 #include "copasi.h"
38 
39 #include "CHybridMethodODE45.h"
40 #include "CTrajectoryProblem.h"
41 #include "model/CModel.h"
42 #include "model/CMetab.h"
43 #include "model/CReaction.h"
44 #include "model/CState.h"
45 #include "model/CChemEq.h"
46 #include "model/CChemEqElement.h"
47 #include "model/CCompartment.h"
49 #include "utilities/CMatrix.h"
52 #include "utilities/CVersion.h"
54 
55 #include "CInterpolation.h"
56 
57 //================Function for Class================
58 /**
59  * Default constructor.
60  */
62  CTrajectoryMethod(CCopasiMethod::hybridODE45, pParent)
63 {
64  assert((void *) &mData == (void *) &mData.dim);
65  mData.pMethod = this;
66 
69 }
70 
71 /**
72  * Copy Constructor
73  */
75  const CCopasiContainer * pParent):
76  CTrajectoryMethod(src, pParent)
77 {
78  assert((void *) &mData == (void *) &mData.dim);
79  mData.pMethod = this;
80 
83 }
84 
85 /**
86  * Destructor.
87  */
89 {
90  cleanup();
92 }
93 
94 //================Function for System================
95 
97 {
99  return true;
100 }
101 
102 //virtual
104 {
105  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
106 
107  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
108 
109  if (pTP->getDuration() < 0.0)
110  {
111  //back integration not possible
113  return false;
114  }
115 
116  //check for rules
117  size_t i, imax = pTP->getModel()->getNumModelValues();
118 
119  for (i = 0; i < imax; ++i)
120  {
121  if (pTP->getModel()->getModelValues()[i]->getStatus() == CModelEntity::ODE)
122  {
123  //ode rule found
125  return false;
126  }
127  }
128 
129  imax = pTP->getModel()->getNumMetabs();
130 
131  for (i = 0; i < imax; ++i)
132  {
133  if (pTP->getModel()->getMetabolites()[i]->getStatus() == CModelEntity::ODE)
134  {
135  //ode rule found
137  return false;
138  }
139  }
140 
141  imax = pTP->getModel()->getCompartments().size();
142 
143  for (i = 0; i < imax; ++i)
144  {
145  if (pTP->getModel()->getCompartments()[i]->getStatus() == CModelEntity::ODE)
146  {
147  //ode rule found
149  return false;
150  }
151  }
152 
153 /*
154  std::string message = pTP->getModel()->suitableForStochasticSimulation();
155 
156  if (message != "")
157  {
158  //model not suitable, message describes the problem
159  CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
160  return false;
161  }
162 */
163 
164  //events are not supported at the moment
165  if (pTP->getModel()->getEvents().size() > 0)
166  {
168  return false;
169  }
170 
171 
172  return true;
173 }
174 
176 {
177  CCopasiParameter *pParm;
178 
179  //-----------------------------------------------------------------------
180  // This part will be dealt with later, when considering the partition strategy.
181  //Now, under the condition of using statistic partition given by users,
182  //we won't adopt such parameters listed here.
183  assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) MAX_STEPS);
184  //assertParameter("Lower Limit", CCopasiParameter::DOUBLE, (C_FAT64) LOWER_STOCH_LIMIT);
185  //assertParameter("Upper Limit", CCopasiParameter::DOUBLE, (C_FLOAT64) UPPER_STOCH_LIMIT);
186 
187  //assertParameter("Partitioning Interval", CCopasiParameter::UINT, (unsigned C_INT32) PARTITIONING_INTERVAL);
188  //-----------------------------------------------------------------------
189 
190  assertParameter("Use Random Seed", CCopasiParameter::BOOL, (bool) USE_RANDOM_SEED);
191  assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) RANDOM_SEED);
192 
193  // Check whether we have a method with the old parameter names
194  if ((pParm = getParameter("HYBRID.MaxSteps")) != NULL)
195  {
196  //-------------------------------------------------------------------
197  setValue("Max Internal Steps", *pParm->getValue().pUINT);
198  removeParameter("HYBRID.MaxSteps");
199 
200  /*
201  if ((pParm = getParameter("HYBRID.LowerStochLimit")) != NULL)
202  {
203  setValue("Lower Limit", *pParm->getValue().pDOUBLE);
204  removeParameter("HYBRID.LowerStochLimit");
205  }
206 
207  if ((pParm = getParameter("HYBRID.UpperStochLimit")) != NULL)
208  {
209  setValue("Upper Limit", *pParm->getValue().pDOUBLE);
210  removeParameter("HYBRID.UpperStochLimit");
211  }
212 
213  if ((pParm = getParameter("HYBRID.PartitioningInterval")) != NULL)//need this part?
214  {
215  setValue("Partitioning Interval", *pParm->getValue().pUINT);
216  removeParameter("HYBRID.PartitioningInterval");
217  }
218  */
219  //-------------------------------------------------------------------
220 
221  if ((pParm = getParameter("UseRandomSeed")) != NULL)
222  {
223  setValue("Use Random Seed", *pParm->getValue().pBOOL);
224  removeParameter("UseRandomSeed");
225  }
226 
227  if ((pParm = getParameter("")) != NULL)
228  {
229  setValue("Random Seed", *pParm->getValue().pUINT);
230  removeParameter("");
231  }
232  }
233 
234  /* ODE45 ****************************************************************************/
235  //----------------------------------------------------------------------
236  //addParameter("Partitioning Stepsize",
237  // CCopasiParameter::DOUBLE, (C_FLOAT64) PARTITIONING_STEPSIZE);
238 
239  //addParameter("Max Internal Steps (LSOA)",
240  // CCopasiParameter::UINT, (unsigned C_INT32) 10000);
241  //-----------------------------------------------------------------------
242 
243  addParameter("Relative Tolerance",
245 
246  addParameter("Use Default Absolute Tolerance",
247  CCopasiParameter::BOOL, (bool) true);
248 
249  addParameter("Absolute Tolerance",
251 
252  //????????????????????????????????????????????????????????
253  // These parameters are no longer supported.
254  //removeParameter("Adams Max Order");
255  //removeParameter("BDF Max Order");
256  //????????????????????????????????????????????????????????
257 }
258 
259 void CHybridMethodODE45::start(const CState * initialState)
260 {
261  //set inital state
262  mpProblem->getModel()->setState(*initialState);
263 
264  //set mpState
266 
267  //set the mpModel
269  assert(mpModel);
270 
271  //set mDoCorrection
272  if (mpModel->getModelType() == CModel::deterministic)
273  mDoCorrection = true;
274  else
275  mDoCorrection = false;
276 
277  //set mHasAssignments
279 
280  //set mFirstMetabIndex
281  mFirstMetabIndex = mpModel->getStateTemplate().getIndex(mpModel->getMetabolitesX()[0]);
282 
283  // Call initMethod function
285 
286  //output data to check init status
287  //outputData();
288  return;
289 }
290 
291 /**
292  * Initializes the solver and sets the model to be used.
293  * @param model A reference to an instance of a CModel
294  */
296 {
297  //(1)----set attributes related with REACTIONS
300 
303 
305 
306  //(2)----set attributes related with METABS
308 
311 
312  setupBalances(); //initialize mLocalBalances and mLocalSubstrates
313  setupMetab2React(); //initialize mMetab2React
314  setupMetabFlags();
315 
316  //(3)----set attributes related with STATE
318 
319  //(4)----set attributes related with SYSTEM
320  mMaxSteps = * getValue("Max Internal Steps").pUINT;
321  mMaxStepsReached = false;
322  setupMethod();
323 
324  //(5)----set attributes related with STOCHASTIC
325  mUseRandomSeed = * getValue("Use Random Seed").pBOOL;
326  mRandomSeed = * getValue("Random Seed").pUINT;
327 
328  if (mUseRandomSeed)
330 
331  mStoi = mpModel->getStoi();
332  mUpdateSet.clear();// how to set this parameter
333  setupCalculateSet(); //should be done after setupBalances()
335 
336  setupDependencyGraph(); // initialize mDG
337  setupPriorityQueue(start_time); // initialize mPQ
338 
339  //(6)----set attributes for INTERPOLATION
342 
343  //(7)----set attributes for ODE45
345  mData.dim = (C_INT)(mNumVariableMetabs + 1); //one more for sum of propensities
347 
348  mRtol = * getValue("Relative Tolerance").pUDOUBLE;
349  mDefaultAtol = * getValue("Use Default Absolute Tolerance").pBOOL;
350 
351  if (mDefaultAtol)
352  {
354  setValue("Absolute Tolerance", mAtol);
355  }
356  else
357  mAtol = * getValue("Absolute Tolerance").pUDOUBLE;
358 
359  mDWork.resize(3 + 6 * mData.dim);
360  mIWork.resize(5);
361 
362  //setupUpdateSet();
363 
364  // we don't want to directly record new results into mpState, since
365  // the sum of slow reaction propensities is also recorded in mY
366  mY = new C_FLOAT64[mData.dim];
367  mOldY = new C_FLOAT64[mData.dim];
368 
369  //first calculate propensities, in the next integration process
370  //system will just update the propensities record in mCalculateSet
371  for (size_t i = 0; i < mNumReactions; i++)
372  calculateAmu(i);
373 
374  return;
375 }
376 
377 /**
378  * Cleans up memory, etc.
379  */
381 {
382  delete mpRandomGenerator;
383  mpRandomGenerator = NULL;
384  mpModel = NULL;
385 
386  delete mpState;
387  mpState = NULL;
388 
389  delete mpInterpolation;
390  mpInterpolation = NULL;
391 
392  delete [] mY;
393  mY = NULL;
394 
395  delete [] mOldY;
396  mOldY = NULL;
397 
398  return;
399 }
400 
401 //================Function for Model================
403 {
404  //size_t rctIndex;
406 
407  std::vector<CHybridODE45MetabFlag>::iterator metabIt
408  = mMetabFlags.begin();
409  const std::vector<CHybridODE45MetabFlag>::iterator metabEndIt
410  = mMetabFlags.end();
411 
412  // initialization
413  for (; metabIt != metabEndIt; ++metabIt)
414  {
415  metabIt->mFlag = SLOW;
416  metabIt->mFastReactions.clear();
417  }
418 
419  std::vector<CHybridODE45Balance>::iterator itBalance;
420  std::vector<CHybridODE45Balance>::iterator itEndBalance;
421 
422  // go over all mLocalBalances, if balance != 0,
423  // and reaction is FAST, insert into set
424  for (size_t rct = 0; rct < mNumReactions; ++rct)
425  {
426  if (mReactionFlags[rct] == FAST)
427  {
428  itBalance = mLocalBalances[rct].begin();
429  itEndBalance = mLocalBalances[rct].end();
430 
431  for (; itBalance != itEndBalance; ++itBalance)
432  {
433  size_t metab = itBalance->mIndex;
434  mMetabFlags[metab].mFastReactions.insert(rct);
435  mMetabFlags[metab].mFlag = FAST;
436  }
437  }
438  }
439 
440  return;
441 }
442 
444 {
445  mHasStoiReaction = false;
446  mHasDetermReaction = false;
448 
449  for (size_t rct = 0; rct < mNumReactions; rct++)
450  {
451  if ((*mpReactions)[rct]->isFast())
452  {
453  mReactionFlags[rct] = FAST;
454  mHasDetermReaction = true;
455  }
456  else
457  {
458  mReactionFlags[rct] = SLOW;
459  mHasStoiReaction = true;
460  }
461  }
462 
463  return;
464 }
465 
467 {
472  else
473  mMethod = HYBRID;
474 
475  return;
476 }
477 
478 /**
479  * Sets up an internal representation of the balances for each reaction.
480  * This is done in order to be able to deal with fixed metabolites and
481  * to avoid a time consuming search for the indices of metabolites in the
482  * model.
483  */
485 {
486  size_t i, j;
487  CHybridODE45Balance newElement;
488  size_t maxBalance = 0;
489 
490  mLocalBalances.clear();
492  mLocalSubstrates.clear();
494 
495  for (i = 0; i < mNumReactions; i++)
496  {
497  const CCopasiVector <CChemEqElement> * balances =
498  &(*mpReactions)[i]->getChemEq().getBalances();
499 
500  for (j = 0; j < balances->size(); j++)
501  {
502  newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
503  newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
504  // + 0.5 to get a rounding out of the static_cast to C_INT32!
505  newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity()
506  + 0.5));
507 
508  if ((newElement.mpMetabolite->getStatus()) != CModelEntity::FIXED)
509  {
510  if (newElement.mMultiplicity > maxBalance) maxBalance = newElement.mMultiplicity;
511 
512  mLocalBalances[i].push_back(newElement); // element is copied for the push_back
513  }
514  }
515 
516  balances = &(*mpReactions)[i]->getChemEq().getSubstrates();
517 
518  for (j = 0; j < balances->size(); j++)
519  {
520  newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
521  newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
522  // + 0.5 to get a rounding out of the static_cast to C_INT32!
523  newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity()
524  + 0.5));
525 
526  mLocalSubstrates[i].push_back(newElement); // element is copied for the push_back
527  }
528  }
529 
530  mMaxBalance = maxBalance;
531 
532  return;
533 }
534 
535 /**
536  * Sets up the priority queue, for pure SSA.
537  *
538  * @param startTime The time at which the simulation starts.
539  */
541 {
542  size_t i;
543  C_FLOAT64 time;
544 
545  mPQ.clear();
547 
548  for (i = 0; i < mNumReactions; i++)
549  {
550  calculateAmu(i);
551  time = startTime + generateReactionTime(i);
552  mPQ.insertStochReaction(i, time);
553  }
554 
555  return;
556 }
557 
558 /**
559  * Creates for each metabolite a set of reaction indices.
560  * If the metabolite participates in a reaction as substrate
561  * or modifiers, this reaction is added to the corresponding set.
562  */
564 {
565  mMetab2React.clear();
567 
568  CMetab * pMetab;
569  CReaction * pReaction;
570 
571  //Check wether metabs play as substrates and modifiers
572  for (size_t rct = 0; rct < mNumReactions; ++rct)
573  {
574  size_t index;
575 
576  //deal with substrates
578  (*mpReactions)[rct]->getChemEq().getSubstrates();
579 
580  for (size_t i = 0; i < metab.size(); i++)
581  {
582  pMetab = const_cast < CMetab* >
583  (metab[i]->getMetabolite());
584  index = mpModel->getMetabolitesX().getIndex(pMetab);
585 
586  if ((pMetab->getStatus()) != CModelEntity::FIXED)
587  mMetab2React[index].insert(rct);
588  }
589 
590  //deal with modifiers
591  metab = (*mpReactions)[rct]->getChemEq().getModifiers();
592 
593  for (size_t i = 0; i < metab.size(); i++)
594  {
595  pMetab = const_cast < CMetab* >
596  (metab[i]->getMetabolite());
597  index = mpModel->getMetabolitesX().getIndex(pMetab);
598 
599  if ((pMetab->getStatus()) != CModelEntity::FIXED)
600  mMetab2React[index].insert(rct);
601  }
602  }
603 
604  //Check wether metabs appear in reaction laws
605  std::set< const CCopasiObject * > changed;
606 
607  for (size_t metab = 0; metab < mNumVariableMetabs; metab++)
608  {
609  pMetab = (*mpMetabolites)[metab];
610 
611  if (pMetab->getStatus() != CModelEntity::FIXED)
612  {
613  changed.clear();
614  changed.insert(mpModel->getValueReference());
615  changed.insert(pMetab->getValueReference());
616 
617  for (size_t react = 0; react < mNumReactions; react++)
618  {
619  pReaction = (*mpReactions)[react];
620 
621  if (pReaction->getParticleFluxReference()->dependsOn(changed))
622  mMetab2React[metab].insert(react);
623  }
624  }
625  }
626 
627  return;
628 }
629 
630 // A fast metab appearring in a reaction as a substrate or
631 // a modifier makes a change of propensity.
633 {
634  mCalculateSet.clear();
635 
636  std::set<size_t>::iterator reactIt;
637 
638  std::vector<CHybridODE45MetabFlag>::iterator metabIt
639  = mMetabFlags.begin();
640  const std::vector<CHybridODE45MetabFlag>::iterator metabEndIt
641  = mMetabFlags.end();
642 
643  std::vector< std::set<size_t> >::iterator m2rIt =
644  mMetab2React.begin();
645 
646  for (; metabIt != metabEndIt; ++metabIt, ++m2rIt)
647  {
648  //for fast metabs, insert reaction into mCalculateSet
649  if (metabIt->mFlag == FAST) //fast metab
650  {
651  for (reactIt = m2rIt->begin();
652  reactIt != m2rIt->end();
653  reactIt++) // check all reactions that involves metab
654  mCalculateSet.insert(*reactIt);
655  }
656  }
657 
658  return;
659 }
660 
662 {
663  mReactAffect.clear();
664  mReactAffect.resize(mNumReactions);
665 
666  std::vector<std::set<size_t> >::iterator rctIt =
667  mReactAffect.begin();
668  const std::vector<std::set<size_t> >::iterator rctEndIt =
669  mReactAffect.end();
670 
671  for (size_t rct = 0; rctIt != rctEndIt; ++rct, ++rctIt)
672  {
673  rctIt->clear();
674 
675  std::vector<CHybridODE45Balance>::iterator balIt =
676  mLocalBalances[rct].begin();
677  const std::vector<CHybridODE45Balance>::iterator balEndIt =
678  mLocalBalances[rct].end();
679 
680  for (; balIt != balEndIt; ++balIt)
681  {
682  size_t metabId = balIt->mIndex;
683 
684  std::set<size_t>::iterator it =
685  mMetab2React[metabId].begin();
686  const std::set<size_t>::iterator endIt =
687  mMetab2React[metabId].end();
688 
689  for (; it != endIt; ++it)
690  rctIt->insert(*it);
691  }
692  }
693 
694  return;
695 }
696 
697 //========Function for Simulation========
699 {
700  // write the current state to the model???
701 
702  // check for possible overflows
703  size_t i;
704  //size_t imax;
705 
706  // :TODO: Bug 774: This assumes that the number of variable metabs is the number
707  // of metabs determined by reaction. In addition they are expected at the beginning of the
708  // MetabolitesX which is not the case if we have metabolites of type ODE.
709  //for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++)
710  // if (mpProblem->getModel()->getMetabolitesX()[i]->getValue() >= mMaxIntBeforeStep)
711  // {
712  // throw exception or something like that
713  //}
714 
715  // do several steps
716  C_FLOAT64 time = mpModel->getTime();
717  C_FLOAT64 endTime = time + deltaT;
718 
719  for (i = 0; ((i < mMaxSteps) && (time < endTime)); i++)
720  time = doSingleStep(time, endTime);
721 
722  // Warning Message
723  if ((i >= mMaxSteps) && (!mMaxStepsReached))
724  {
725  mMaxStepsReached = true; //only report this message once
726  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.");
727  }
728 
729  // get back the particle numbers
730 
731  /* Set the variable metabolites ????*/
732  //C_FLOAT64 * Dbl = mpCurrentState->beginIndependent() + mFirstMetabIndex - 1;
733 
734  //for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, Dbl++)
735  // *Dbl = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
736 
737  //the task expects the result in mpCurrentState
738  //*mpCurrentState = mpProblem->getModel()->getState();
740  //mpCurrentState->setTime(time);
741 
742  return NORMAL;
743 }
744 
745 /**
746  * Simulates the system over the next interval of time. The new time after
747  * this step is returned.
748  *
749  * @param currentTime A C_FLOAT64 specifying the current time
750  * @param endTime A C_FLOAT64 specifying the endTime of the current step()
751  * @return A C_FLOAT giving the new time
752  */
754 {
755  C_FLOAT64 ds = 0.0;
756  size_t rIndex = 0;
757 
758  //1----pure SSA method
759  if (mMethod == STOCHASTIC) //has only stochastic reactions
760  {
761  getStochTimeAndIndex(ds, rIndex);
762 
763  if (ds > endTime)
764  ds = endTime;
765  else
766  {
767  fireReaction(rIndex);
768  updatePriorityQueue(rIndex, ds);
769  }
770 
771  //population has been changed and record to the mCurrentState
772  //in fireReaction(). So here just store time
773  *mpState = mpModel->getState();
774  mpState->setTime(ds);
776 
777  //outputState(mpCurrentState);
778  //getchar();
779  }
780  //2----Method with Deterministic Part
781  else if (mMethod == DETERMINISTIC) //has only deterministic reactions
782  {
783  integrateDeterministicPart(endTime - currentTime);
784  ds = mpState->getTime();
786  // Till now, state has been recorded
787  }
788  else if (mMethod == HYBRID)//Hybrid Method
789  {
790  integrateDeterministicPart(endTime - currentTime);
791  ds = mpState->getTime();
792 
793  if (mODE45Status == HAS_EVENT) //fire slow reaction
794  {
796 
799  }
800  else
802  }
803 
804  return ds;
805 }
806 
808 {
809  //First Update Propensity
810  std::set <size_t>::iterator reactIt = mCalculateSet.begin();
811  size_t reactID;
812 
813  for (; reactIt != mCalculateSet.end(); reactIt++)
814  {
815  reactID = *reactIt;
816  calculateAmu(reactID);
817  }
818 
819  reactID = getReactionIndex4Hybrid();
820  fireReaction(reactID); //Directly update current state in global view
821  *mpState = mpModel->getState();
822 
823  //update corresponding propensities
824  std::set <size_t>::iterator updateIt
825  = mReactAffect[reactID].begin();
826  const std::set <size_t>::iterator updateEndIt
827  = mReactAffect[reactID].end();
828 
829  for (; updateIt != updateEndIt; updateIt++)
830  {
831  reactID = *updateIt;
832  calculateAmu(reactID);
833  }
834 
835  return;
836 }
837 
838 //========Function for ODE45========
840 {
841  if (!pModel)
842  return 1.0e009;
843 
844  const CCopasiVectorNS< CCompartment > & Compartment =
845  pModel->getCompartments();
846 
847  size_t i, imax;
848 
850 
851  for (i = 0, imax = Compartment.size(); i < imax; i++)
852  {
853  if (Compartment[i]->getValue() < Volume)
854  Volume = Compartment[i]->getValue();
855  }
856 
858  return 1.0e009;
859 
860  return Volume * pModel->getQuantity2NumberFactor() * 1.e-12;
861 }
862 
863 /**
864  * Integrates the deterministic reactions of the system over
865  * the specified time interval.
866  *
867  * @param ds A C_FLOAT64 specifying the stepsize.
868  */
869 
871 {
872 
873  //1----Set Parameters for ODE45 solver
874  //=(1)= solver error message
875  mErrorMsg.str("");
876 
877  //=(2)= set state
878  *mpState = mpModel->getState();
879 
880  //=(3)= set time and old time
881  mTime = mpState->getTime();
882  C_FLOAT64 EndTime = mTime + deltaT;
883 
884  mOldTime = mTime;
885 
886  //=(4)= set y and old_y
887  C_FLOAT64 * stateY = mpState->beginIndependent();
888 
889  for (size_t i = 0; i < mData.dim - 1; i++, stateY++)
890  {
891  mOldY[i] = *stateY;
892  mY[i] = *stateY;
893  }
894 
895  if (mODE45Status == NEW_STEP)
896  {
898  mY[mData.dim - 1] = log(rand2);
899  }
900 
901  mOldY[mData.dim - 1] = mY[mData.dim - 1];
902 
903  //2----Reset ODE45 Solver
904  mIFlag = -1; //integration by one step
905 
906  //3----If time increment is too small, do nothing
907  C_FLOAT64 tdist , d__1, d__2, w0;
908  tdist = fabs(deltaT); //absolute time increment
909  d__1 = fabs(mTime), d__2 = fabs(EndTime);
910  w0 = std::max(d__1, d__2);
911 
912  if (tdist < std::numeric_limits< C_FLOAT64 >::epsilon() * 2. * w0) //just do nothing
913  {
914  mpModel->setTime(EndTime);
915  return;
916  }
917 
918  //4----just do nothing if there are no variables
919  if (!mData.dim)
920  {
921  mpModel->setTime(EndTime);
922  return;
923  }
924 
925  while ((mODE45Status == CONTINUE)
926  || (mODE45Status == NEW_STEP))
927  {
928  //(1)----ODE Solver
929  rkf45_(&EvalF , //1. evaluate F
930  &mData.dim, //2. number of variables
931  mY , //3. the array of current concentrations
932  &mTime , //4. the current time
933  &EndTime , //5. the final time
934  &mRtol , //6. relative tolerance array
935  &mAtol , //7. absolute tolerance array
936  &mIFlag , //8. the state for input and output
937  mDWork.array() , //9. the double work array
938  mIWork.array() , //10. the int work array
939  mStateRecord.array()); //11. the array to record temp state
940 
941  //(2)----set mODE45Status
942  if (mIFlag == 2)
943  {
944  if (mY[mData.dim - 1] >= 0.0)
946  else
947  mODE45Status = REACH_END_TIME; //reach EndTime
948  }
949  else if (mIFlag == -2) //finish one step integration
950  {
951  if (mY[mData.dim - 1] >= 0.0) //has an event
953  else
955  }
956  else //error happens
957  {
959  std::cout << "Error Type: " << mIFlag << std::endl;
961  }//end if
962  }//end while
963 
964  //5----Record State
965  stateY = mpState->beginIndependent();
966 
967  for (size_t i = 0; i < mData.dim - 1; i++)
968  stateY[i] = mY[i]; //write result into mpState
969 
972  mpModel->updateSimulatedValues(false); //for assignments?????????
973 
974  return;
975 }
976 
978 {
979 
980  //==(1)==for one-step method, reset record each time when do interpolation
982 
983  //==(2)==set record in class interpolation
985  size_t offset;
986 
987  if (mUseStateRecord)
988  {
989  for (size_t i = 0; i < INTERP_RECORD_NUM - 2; i++) //record the middle 4 states
990  {
991  offset = i * (mData.dim + 1);
993  mStateRecord.array() + offset + 1);
994  }
995  }
996 
998 
999  //==(3)==do interpolation
1001 
1002  //==(4)==record the state
1003  mTime = mpEventState->getTime();
1004  C_FLOAT64 * tmpY = mpEventState->getArray() + 1;
1005  C_FLOAT64 * stateY = mpState->beginIndependent();
1006 
1007  for (size_t i = 0; i < mData.dim - 1; i++)
1008  stateY[i] = tmpY[i]; //write result into mpState
1009 
1010  mpState->setTime(mTime);
1012  mpModel->updateSimulatedValues(false); //for assignments?????????
1013 
1014  return;
1015 }
1016 
1017 /**
1018  * Dummy f function for calculating derivative of y
1019  */
1020 void CHybridMethodODE45::EvalF(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot)
1021 {static_cast<Data *>((void *) n)->pMethod->evalF(t, y, ydot);}
1022 
1023 /**
1024  * Derivative Calculation Function
1025  */
1026 void CHybridMethodODE45::evalF(const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot)
1027 {
1028  size_t i;
1029 
1030  //========No need to record State, right?========
1031  //maybe here I just want the propensities amu???
1032  //calculateDerivative(temp);
1033 
1034  //1====calculate propensities
1035  //just care about reactions involving fast metab
1036 
1037  //(1)Put current time *t and independent values *y into model.
1038  // This step seemes necessary, since propensity calculation process
1039  // requires functions called from the model class.
1040  mpState->setTime(*t);
1041  C_FLOAT64 * tmpY = mpState->beginIndependent();//mpState is a local copy
1042 
1043  for (i = 0; i < mData.dim - 1; ++i)
1044  tmpY[i] = y[i];
1045 
1047  mpModel->updateSimulatedValues(false); //really?
1048 
1049  //(2)Calculate propensities.
1050  std::set <size_t>::iterator reactIt = mCalculateSet.begin();
1051  size_t reactID;
1052 
1053  for (; reactIt != mCalculateSet.end(); reactIt++)
1054  {
1055  reactID = *reactIt;
1056  calculateAmu(reactID);
1057  }
1058 
1059  //2====calculate derivative
1060  //(1)initialize
1061  for (i = 0; i < mData.dim; i++)
1062  ydot[i] = 0.0;
1063 
1064  //(2)go through all the reactions and
1065  //update derivatives
1066  std::vector <CHybridODE45Balance>::iterator metabIt;
1067  std::vector <CHybridODE45Balance>::iterator metabEndIt;
1068  size_t metabIndex;
1069 
1070  for (i = 0; i < mNumReactions; i++)
1071  {
1072  if (mReactionFlags[i] == SLOW) //slow reaction
1073  ydot[mData.dim - 1] += mAmu[i];
1074  else //fast reaction
1075  {
1076  metabIt = mLocalBalances[i].begin();
1077  metabEndIt = mLocalBalances[i].end();
1078 
1079  for (; metabIt != metabEndIt; metabIt++)
1080  {
1081  metabIndex = metabIt->mIndex;
1082  ydot[metabIndex] += metabIt->mMultiplicity * mAmu[i];
1083  }
1084  }
1085  }
1086 
1087  return;
1088 }
1089 
1090 //========Function for Stoichastic========
1091 /**
1092  * Find the reaction index and the reaction time of the stochastic (!)
1093  * reaction with the lowest reaction time.
1094  *
1095  * @param ds A reference to a C_FLOAT64. The putative reaction time for the
1096  * first stochastic reaction is written into this variable.
1097  * @param rIndex A reference to a size_t. The index of the first
1098  * stochastic reaction is written into this variable.
1099  */
1101 {
1102  ds = mPQ.topKey();
1103  rIndex = mPQ.topIndex();
1104  return;
1105 }
1106 
1108 {
1109  if (mAmu[rIndex] == 0) return std::numeric_limits<C_FLOAT64>::infinity();
1110 
1112  return - 1 * log(rand2) / mAmu[rIndex];
1113 }
1114 
1115 /**
1116  * Updates the putative reaction time of a stochastic reaction in the
1117  * priority queue. The corresponding amu and amu_old must be set prior to
1118  * the call of this method.
1119  *
1120  * @param rIndex A size_t specifying the index of the reaction
1121  */
1123 {
1124  C_FLOAT64 newTime;
1125 
1126  // One must make sure that the calculation yields reasonable results even in the cases
1127  // where mAmu=0 or mAmuOld=0 or both=0. Therefore mAmuOld=0 is checked. If mAmuOld equals 0,
1128  // then a new random number has to be drawn, because tau equals inf and the stoch.
1129  // information is lost.
1130  // If both values equal 0, then tau should remain inf and the update of the queue can
1131  // be skipped!
1132  if (mAmuOld[rIndex] == 0.0)
1133  {
1134  if (mAmu[rIndex] != 0.0)
1135  {
1136  newTime = time + generateReactionTime(rIndex);
1137  mPQ.updateNode(rIndex, newTime);
1138  }
1139  }
1140  else
1141  {
1142  newTime = time + (mAmuOld[rIndex] / mAmu[rIndex]) * (mPQ.getKey(rIndex) - time);
1143  mPQ.updateNode(rIndex, newTime);
1144  }
1145 
1146  return;
1147 }
1148 
1149 /**
1150  * Calculates an amu value for a given reaction.
1151  *
1152  * @param rIndex A size_t specifying the reaction to be updated
1153  */
1155 {
1156  if (!mDoCorrection)
1157  {
1158  mAmu[rIndex] = (*mpReactions)[rIndex]->calculateParticleFlux();
1159  return;
1160  }
1161 
1162  // We need the product of the cmu and hmu for this step.
1163  // We calculate this in one go, as there are fewer steps to
1164  // perform and we eliminate some possible rounding errors.
1165  C_FLOAT64 amu = 1; // initially
1166  //size_t total_substrates = 0;
1167  C_INT32 num_ident = 0;
1168  C_INT32 number = 0;
1169  C_INT32 lower_bound;
1170  // substrate_factor - The substrates, raised to their multiplicities,
1171  // multiplied with one another. If there are, e.g. m substrates of type m,
1172  // and n of type N, then substrate_factor = M^m * N^n.
1173  C_FLOAT64 substrate_factor = 1;
1174  // First, find the reaction associated with this index.
1175  // Keep a pointer to this.
1176  // Iterate through each substrate in the reaction
1177  const std::vector<CHybridODE45Balance> & substrates = mLocalSubstrates[rIndex];
1178 
1179  int flag = 0;
1180 
1181  for (size_t i = 0; i < substrates.size(); i++)
1182  {
1183  num_ident = substrates[i].mMultiplicity;
1184 
1185  if (num_ident > 1)
1186  {
1187  flag = 1;
1188  number = static_cast<C_INT32>((*mpMetabolites)[substrates[i].mIndex]->getValue());
1189  lower_bound = number - num_ident;
1190  substrate_factor = substrate_factor
1191  * pow((double) number, (int) num_ident - 1); //optimization
1192 
1193  number--; // optimization
1194 
1195  while (number > lower_bound)
1196  {
1197  amu *= number;
1198  number--;
1199  }
1200  }
1201  }
1202 
1203  if ((amu == 0) || (substrate_factor == 0)) // at least one substrate particle number is zero
1204  {
1205  mAmu[rIndex] = 0;
1206  return;
1207  }
1208 
1209  // rate_factor is the rate function divided by substrate_factor.
1210  // It would be more efficient if this was generated directly, since in effect we
1211  // are multiplying and then dividing by the same thing (substrate_factor)!
1212  C_FLOAT64 rate_factor = (*mpReactions)[rIndex]->calculateParticleFlux();
1213 
1214  if (flag)
1215  {
1216  amu *= rate_factor / substrate_factor;;
1217  mAmu[rIndex] = amu;
1218  }
1219  else
1220  {
1221  mAmu[rIndex] = rate_factor;
1222  }
1223 
1224  return;
1225 
1226  // a more efficient way to calculate mass action kinetics could be included
1227 }
1228 
1229 /**
1230  * Gets the set of metabolites on which a given reaction depends.
1231  *
1232  * @param rIndex The index of the reaction being executed.
1233  * @return The set of metabolites depended on.
1234  */
1235 std::set<std::string> *CHybridMethodODE45::getDependsOn(size_t rIndex)
1236 {
1237  std::set<std::string> *retset = new std::set<std::string>;
1238 
1239  size_t i, imax = (*mpReactions)[rIndex]->getFunctionParameters().size();
1240  size_t j, jmax;
1241 
1242  for (i = 0; i < imax; ++i)
1243  {
1244  if ((*mpReactions)[rIndex]->getFunctionParameters()[i]->getUsage()
1246  continue;
1247 
1248  const std::vector <std::string> & metabKeylist =
1249  (*mpReactions)[rIndex]->getParameterMappings()[i];
1250  jmax = metabKeylist.size();
1251 
1252  for (j = 0; j < jmax; ++j)
1253  retset->insert(metabKeylist[j]);
1254  }
1255 
1256  return retset;
1257 }
1258 
1259 /**
1260  * Gets the set of metabolites which change number when a given
1261  * reaction is executed.
1262  *
1263  * @param rIndex The index of the reaction being executed.
1264  * @return The set of affected metabolites.
1265  */
1266 std::set<std::string> *CHybridMethodODE45::getAffects(size_t rIndex)
1267 {
1268  size_t i;
1269  std::set<std::string> *retset = new std::set<std::string>;
1270 
1271  // Get the balances associated with the reaction at this index
1272  // XXX We first get the chemical equation, then the balances, since the getBalances method in CReaction is unimplemented!
1273 
1274  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
1275  {
1276  if (mLocalBalances[rIndex][i].mMultiplicity != 0)
1277  retset->insert(mLocalBalances[rIndex][i].mpMetabolite->getKey());
1278  }
1279 
1280  return retset;
1281 }
1282 
1283 /**
1284  * Sets up the dependency graph.
1285  */
1287 {
1288  mDG.clear();
1289  std::vector< std::set<std::string>* > DependsOn;
1290  std::vector< std::set<std::string>* > Affects;
1291  size_t numReactions = mpReactions->size();
1292  size_t i, j;
1293 
1294  // Do for each reaction:
1295  for (i = 0; i < mNumReactions; i++)
1296  {
1297  // Get the set of metabolites which affect the value of amu for this
1298  // reaction i.e. the set on which amu depends. This may be more than
1299  // the set of substrates, since the kinetics can involve other
1300  // reactants, e.g. catalysts. We thus need to step through the
1301  // rate function and pick out every reactant which can vary.
1302  DependsOn.push_back(getDependsOn(i));
1303  // Get the set of metabolites which are affected when this reaction takes place
1304  Affects.push_back(getAffects(i));
1305  }
1306 
1307  // For each possible pair of reactions i and j, if the intersection of
1308  // Affects(i) with DependsOn(j) is non-empty, add a dependency edge from i to j.
1309  for (i = 0; i < mNumReactions; i++)
1310  {
1311  for (j = 0; j < mNumReactions; j++)
1312  {
1313  // Determine whether the intersection of these two sets is non-empty
1314  // Could also do this with set_intersection generic algorithm, but that
1315  // would require operator<() to be defined on the set elements.
1316 
1317  std::set<std::string>::iterator iter = Affects[i]->begin();
1318 
1319  for (; iter != Affects[i]->end(); iter++)
1320  {
1321  if (DependsOn[j]->count(*iter))
1322  {
1323  // The set intersection is non-empty
1324  mDG.addDependent(i, j);
1325  break;
1326  }
1327  }
1328  }
1329 
1330  // Ensure that self edges are included
1331  //mDG.addDependent(i, i);
1332  }
1333 
1334  // Delete the memory allocated in getDependsOn() and getAffects()
1335  // since this is allocated in other functions.
1336  for (i = 0; i < numReactions; i++)
1337  {
1338  delete DependsOn[i];
1339  delete Affects[i];
1340  }
1341 
1342  return;
1343 }
1344 
1345 /**
1346  * Updates the priority queue.
1347  *
1348  * @param rIndex A size_t giving the index of the fired reaction (-1, if no
1349  * stochastic reaction has fired)
1350  * @param time A C_FLOAT64 holding the current time
1351  */
1352 
1354 {
1355  C_FLOAT64 newTime;
1356  size_t index;
1357  std::set <size_t>::iterator iter, iterEnd;
1358 
1359  //if the model contains assignments we use a less efficient loop over all (stochastic)
1360  //reactions to capture all changes
1361  //we do not know the exact dependencies.
1362  //TODO: this should be changed later in order to get a more efficient update scheme
1363  if (mHasAssignments)
1364  {
1366 
1367  for (index = 0; index < mNumReactions; index++)
1368  {
1369  mAmuOld[index] = mAmu[index];
1370  calculateAmu(index);
1371 
1372  if (mAmuOld[index] != mAmu[index])
1373  if (index != rIndex) updateTauMu(index, time);
1374  }
1375  }
1376  else
1377  {
1378  // iterate through the set of affected reactions and update the stochastic
1379  // ones in the priority queue
1380  iter = mReactAffect[rIndex].begin();
1381  iterEnd = mReactAffect[rIndex].end();
1382 
1383  for (; iter != iterEnd; iter++)
1384  {
1385  index = *iter;
1386  mAmuOld[index] = mAmu[index];
1387  calculateAmu(index);
1388 
1389  if (*iter != rIndex) updateTauMu(index, time);
1390  }
1391  }
1392 
1393  // draw new random number and update the reaction just fired
1394  if (rIndex != C_INVALID_INDEX)
1395  {
1396  // reaction is stochastic
1397  newTime = time + generateReactionTime(rIndex);
1398  mPQ.updateNode(rIndex, newTime);
1399  }
1400 
1401  return;
1402 }
1403 
1404 /**
1405  * Executes the specified reaction in the system once.
1406  *
1407  * @param rIndex A size_t specifying the index of the reaction, which
1408  * will be fired.
1409  */
1411 {
1412  // Change the particle numbers according to which step took place.
1413  // First, get the vector of balances in the reaction we've got.
1414  // (This vector expresses the number change of each metabolite
1415  // in the reaction.) Then step through each balance, using its
1416  // multiplicity to calculate a new value for the associated
1417  // metabolite. Finally, update the metabolite.
1418 
1419  size_t i;
1420  C_FLOAT64 newNumber;
1421 
1422  CMetab * pMetab;
1423 
1424  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
1425  {
1426  pMetab = mLocalBalances[rIndex][i].mpMetabolite;
1427  newNumber = pMetab->getValue() + mLocalBalances[rIndex][i].mMultiplicity;
1428 
1429  pMetab->setValue(newNumber);
1430  pMetab->refreshConcentration();
1431  }
1432 
1433  return;
1434 }
1435 
1436 /**
1437  *
1438  *
1439  */
1440 
1442 {
1443  //calculate sum of amu
1444  C_FLOAT64 mAmuSum = 0.0;
1445 
1446  for (int i = 0; i < mNumReactions; i++)
1447  {
1448  if (mReactionFlags[i] == SLOW)
1449  mAmuSum += mAmu[i];
1450  }
1451 
1452  //get the threshold
1454  C_FLOAT64 threshold = mAmuSum * rand2;
1455 
1456  //get the reaction index
1457  size_t rIndex;
1458  C_FLOAT64 tmp = 0.0;
1459 
1460  //is there some algorithm that can get a log() complex?
1461  for (rIndex = 0; rIndex < mNumReactions; rIndex++)
1462  {
1463  if (mReactionFlags[rIndex] == SLOW)
1464  {
1465  tmp += mAmu[rIndex];
1466 
1467  if (tmp >= threshold)
1468  return rIndex;
1469  }
1470  }
1471  return rIndex;
1472 }
1473 
1474 /*========Functions for C Code from f2c========*/
1475 C_INT CHybridMethodODE45::rkf45_(pEvalF f, const C_INT *neqn, double *y,
1476  double *t, double *tout,
1477  double *relerr, double *abserr,
1478  C_INT *iflag, double *work,
1479  C_INT *iwork, double *yrcd)
1480 {
1481  static C_INT k1, k2, k3, k4, k5, k6, k1m;
1482 
1483  /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
1484  /* This is a modefied fehlberg fourth-fifth order runge-kutta method */
1485  /* based on h.a.watts and l.f.shapine's fortran ODE solver rkf45.f */
1486  /* with one more input double precision vector recording time */
1487  /* points between t and t+h in each steap and */
1488  /* approximations of y in such corresponding times. */
1489  /* yrcd is a vector of length 4*(neqn+1), where "4" is related with */
1490  /* four middle time points: */
1491  /* c1=t+h/4, c2=t+3h/8, c3=t+h/2, c4=t+12h/13. */
1492  /* yrcd can be seperated into four parts, each of which is a double */
1493  /* precision vector of length neqn+1. The first member is assigned */
1494  /* to time record, and the last neqn members are for y at the */
1495  /* corresponding time. */
1496  /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
1497 
1498  /* fehlberg fourth-fifth order runge-kutta method */
1499 
1500  /* written by h.a.watts and l.f.shampine */
1501  /* sandia laboratories */
1502  /* albuquerque,new mexico */
1503 
1504  /* rkf45 is primarily designed to solve non-stiff and mildly stiff */
1505  /* differential equations when derivative evaluations are inexpensive. */
1506  /* rkf45 should generally not be used when the user is demanding */
1507  /* high accuracy. */
1508 
1509  /* abstract */
1510 
1511  /* subroutine rkf45 integrates a system of neqn first order */
1512  /* ordinary differential equations of the form */
1513  /* dy(i)/dt = f(t,y(1),y(2),...,y(neqn)) */
1514  /* where the y(i) are given at t . */
1515  /* typically the subroutine is used to integrate from t to tout but it */
1516  /* can be used as a one-step integrator to advance the solution a */
1517  /* single step in the direction of tout. on return the parameters in */
1518  /* the call list are set for continuing the integration. the user has */
1519  /* only to call rkf45 again (and perhaps define a new value for tout). */
1520  /* actually, rkf45 is an interfacing routine which calls subroutine */
1521  /* rkfs for the solution. rkfs in turn calls subroutine fehl which */
1522  /* computes an approximate solution over one step. */
1523 
1524  /* rkf45 uses the runge-kutta-fehlberg (4,5) method described */
1525  /* in the reference */
1526  /* e.fehlberg , low-order classical runge-kutta formulas with stepsize */
1527  /* control , nasa tr r-315 */
1528 
1529  /* the performance of rkf45 is illustrated in the reference */
1530  /* l.f.shampine,h.a.watts,s.davenport, solving non-stiff ordinary */
1531  /* differential equations-the state of the art , */
1532  /* sandia laboratories report sand75-0182 , */
1533  /* to appear in siam review. */
1534 
1535  /* the parameters represent- */
1536  /* f -- subroutine f(t,y,yp) to evaluate derivatives yp(i)=dy(i)/dt */
1537  /* neqn -- number of equations to be integrated */
1538  /* y(*) -- solution vector at t */
1539  /* t -- independent variable */
1540  /* tout -- output point at which solution is desired */
1541  /* relerr,abserr -- relative and absolute error tolerances for local */
1542  /* error test. at each step the code requires that */
1543  /* abs(local error) .le. relerr*abs(y) + abserr */
1544  /* for each component of the local error and solution vectors */
1545  /* iflag -- indicator for status of integration */
1546  /* work(*) -- array to hold information internal to rkf45 which is */
1547  /* necessary for subsequent calls. must be dimensioned */
1548  /* at least 3+6*neqn */
1549  /* iwork(*) -- C_INT array used to hold information internal to */
1550  /* rkf45 which is necessary for subsequent calls. must be */
1551  /* dimensioned at least 5 */
1552 
1553  /* first call to rkf45 */
1554 
1555  /* the user must provide storage in his calling program for the arrays */
1556  /* in the call list - y(neqn) , work(3+6*neqn) , iwork(5) , */
1557  /* declare f in an external statement, supply subroutine f(t,y,yp) and */
1558  /* initialize the following parameters- */
1559 
1560  /* neqn -- number of equations to be integrated. (neqn .ge. 1) */
1561  /* y(*) -- vector of initial conditions */
1562  /* t -- starting point of integration , must be a variable */
1563  /* tout -- output point at which solution is desired. */
1564  /* t=tout is allowed on the first call only, in which case */
1565  /* rkf45 returns with iflag=2 if continuation is possible. */
1566  /* relerr,abserr -- relative and absolute local error tolerances */
1567  /* which must be non-negative. relerr must be a variable while */
1568  /* abserr may be a constant. the code should normally not be */
1569  /* used with relative error control smaller than about 1.e-8 . */
1570  /* to avoid limiting precision difficulties the code requires */
1571  /* relerr to be larger than an internally computed relative */
1572  /* error parameter which is machine dependent. in particular, */
1573  /* pure absolute error is not permitted. if a smaller than */
1574  /* allowable value of relerr is attempted, rkf45 increases */
1575  /* relerr appropriately and returns control to the user before */
1576  /* continuing the integration. */
1577  /* iflag -- +1,-1 indicator to initialize the code for each new */
1578  /* problem. normal input is +1. the user should set iflag=-1 */
1579  /* only when one-step integrator control is essential. in this */
1580  /* case, rkf45 attempts to advance the solution a single step */
1581  /* in the direction of tout each time it is called. since this */
1582  /* mode of operation results in extra computing overhead, it */
1583  /* should be avoided unless needed. */
1584 
1585  /* output from rkf45 */
1586 
1587  /* y(*) -- solution at t */
1588  /* t -- last point reached in integration. */
1589  /* iflag = 2 -- integration reached tout. indicates successful retur */
1590  /* and is the normal mode for continuing integration. */
1591  /* =-2 -- a single successful step in the direction of tout */
1592  /* has been taken. normal mode for continuing */
1593  /* integration one step at a time. */
1594  /* = 3 -- integration was not completed because relative error */
1595  /* tolerance was too small. relerr has been increased */
1596  /* appropriately for continuing. */
1597  /* = 4 -- integration was not completed because more than */
1598  /* 3000 derivative evaluations were needed. this */
1599  /* is approximately 500 steps. */
1600  /* = 5 -- integration was not completed because solution */
1601  /* vanished making a pure relative error test */
1602  /* impossible. must use non-zero abserr to continue. */
1603  /* using the one-step integration mode for one step */
1604  /* is a good way to proceed. */
1605  /* = 6 -- integration was not completed because requested */
1606  /* accuracy could not be achieved using smallest */
1607  /* allowable stepsize. user must increase the error */
1608  /* tolerance before continued integration can be */
1609  /* attempted. */
1610  /* = 7 -- it is likely that rkf45 is inefficient for solving */
1611  /* this problem. too much output is restricting the */
1612  /* natural stepsize choice. use the one-step integrator */
1613  /* mode. */
1614  /* = 8 -- invalid input parameters */
1615  /* this indicator occurs if any of the following is */
1616  /* satisfied - neqn .le. 0 */
1617  /* t=tout and iflag .ne. +1 or -1 */
1618  /* relerr or abserr .lt. 0. */
1619  /* iflag .eq. 0 or .lt. -2 or .gt. 8 */
1620  /* work(*),iwork(*) -- information which is usually of no interest */
1621  /* to the user but necessary for subsequent calls. */
1622  /* work(1),...,work(neqn) contain the first derivatives */
1623  /* of the solution vector y at t. work(neqn+1) contains */
1624  /* the stepsize h to be attempted on the next step. */
1625  /* iwork(1) contains the derivative evaluation counter. */
1626 
1627  /* subsequent calls to rkf45 */
1628 
1629  /* subroutine rkf45 returns with all information needed to continue */
1630  /* the integration. if the integration reached tout, the user need onl */
1631  /* define a new tout and call rkf45 again. in the one-step integrator */
1632  /* mode (iflag=-2) the user must keep in mind that each step taken is */
1633  /* in the direction of the current tout. upon reaching tout (indicated */
1634  /* by changing iflag to 2),the user must then define a new tout and */
1635  /* reset iflag to -2 to continue in the one-step integrator mode. */
1636 
1637  /* if the integration was not completed but the user still wants to */
1638  /* continue (iflag=3,4 cases), he just calls rkf45 again. with iflag=3 */
1639  /* the relerr parameter has been adjusted appropriately for continuing */
1640  /* the integration. in the case of iflag=4 the function counter will */
1641  /* be reset to 0 and another 3000 function evaluations are allowed. */
1642 
1643  /* however,in the case iflag=5, the user must first alter the error */
1644  /* criterion to use a positive value of abserr before integration can */
1645  /* proceed. if he does not,execution is terminated. */
1646 
1647  /* also,in the case iflag=6, it is necessary for the user to reset */
1648  /* iflag to 2 (or -2 when the one-step integration mode is being used) */
1649  /* as well as increasing either abserr,relerr or both before the */
1650  /* integration can be continued. if this is not done, execution will */
1651  /* be terminated. the occurrence of iflag=6 indicates a trouble spot */
1652  /* (solution is changing rapidly,singularity may be present) and it */
1653  /* often is inadvisable to continue. */
1654 
1655  /* if iflag=7 is encountered, the user should use the one-step */
1656  /* integration mode with the stepsize determined by the code or */
1657  /* consider switching to the adams codes de/step,intrp. if the user */
1658  /* insists upon continuing the integration with rkf45, he must reset */
1659  /* iflag to 2 before calling rkf45 again. otherwise,execution will be */
1660  /* terminated. */
1661 
1662  /* if iflag=8 is obtained, integration can not be continued unless */
1663  /* the invalid input parameters are corrected. */
1664 
1665  /* it should be noted that the arrays work,iwork contain information */
1666  /* required for subsequent integration. accordingly, work and iwork */
1667  /* should not be altered. */
1668 
1669  /* compute indices for the splitting of the work array */
1670 
1671  /* Parameter adjustments */
1672  --yrcd;
1673  --y;
1674  --work;
1675  --iwork;
1676 
1677  /* Function Body */
1678  k1m = *neqn + 1;
1679  //k1m = *neqn;
1680  k1 = k1m + 1;
1681  k2 = k1 + *neqn;
1682  k3 = k2 + *neqn;
1683  k4 = k3 + *neqn;
1684  k5 = k4 + *neqn;
1685  k6 = k5 + *neqn;
1686 
1687  /* this interfacing routine merely relieves the user of a long */
1688  /* calling list via the splitting apart of two working storage */
1689  /* arrays. if this is not compatible with the users compiler, */
1690  /* he must use rkfs directly. */
1691 
1692  rkfs_((pEvalF)f, neqn, &y[1], t, tout, relerr, abserr, iflag,
1693  &work[1], &work[k1m], &work[k1], &work[k2], &work[k3],
1694  &work[k4], &work[k5], &work[k6], &work[k6 + 1], &iwork[1],
1695  &iwork[2], &iwork[3], &iwork[4], &iwork[5], &yrcd[1]);
1696 
1697  return 0;
1698 } /* rkf45_ */
1699 
1700 C_INT CHybridMethodODE45::rkfs_(pEvalF f, const C_INT *neqn, double *y, double *
1701  t, double *tout, double *relerr, double *abserr, C_INT *
1702  iflag, double *yp, double *h__, double *f1, double *
1703  f2, double *f3, double *f4, double *f5, double *savre,
1704  double *savae, C_INT *nfe, C_INT *kop, C_INT *init,
1705  C_INT *jflag, C_INT *kflag, double *yrcd)
1706 {
1707  /* Initialized data */
1708 
1709  static double remin = 1e-12;
1710  static C_INT maxnfe = 3000;
1711 
1712  /* System generated locals */
1713  C_INT i__1;
1714  double d__1, d__2, d__3, d__4;
1715 
1716  /* Builtin functions */
1717  /* Subroutine */
1718  //int s_stop(char *, ftnlen);
1719  //double pow_dd(double *, double *);
1720  //C_INT i_sign(C_INT *, C_INT *);
1721  //double d_sign(double *, double *);
1722 
1723  /* Local variables */
1724  static double a;
1725  static C_INT k;
1726  static double s, ae, ee, dt, et, u26, rer, tol, ypk;
1727  static double hmin, toln;
1728  static C_INT mflag;
1729  static double scale, eeoet;
1730  static int hfaild;
1731  static double esttol, twoeps;
1732  static int output;
1733 
1734  /* fehlberg fourth-fifth order runge-kutta method */
1735 
1736  /* rkfs integrates a system of first order ordinary differential */
1737  /* equations as described in the comments for rkf45 . */
1738  /* the arrays yp,f1,f2,f3,f4,and f5 (of dimension at least neqn) and */
1739  /* the variables h,savre,savae,nfe,kop,init,jflag,and kflag are used */
1740  /* internally by the code and appear in the call list to eliminate */
1741  /* local retention of variables between calls. accordingly, they */
1742  /* should not be altered. items of possible interest are */
1743  /* yp - derivative of solution vector at t */
1744  /* h - an appropriate stepsize to be used for the next step */
1745  /* nfe- counter on the number of derivative function evaluations */
1746 
1747  /* remin is the minimum acceptable value of relerr. attempts */
1748  /* to obtain higher accuracy with this subroutine are usually */
1749  /* very expensive and often unsuccessful. */
1750 
1751  /* Parameter adjustments */
1752  --yrcd;
1753  --f5;
1754  --f4;
1755  --f3;
1756  --f2;
1757  --f1;
1758  --yp;
1759  --y;
1760 
1761  /* Function Body */
1762 
1763  /* the expense is controlled by restricting the number */
1764  /* of function evaluations to be approximately maxnfe. */
1765  /* as set, this corresponds to about 500 steps. */
1766 
1767  /* here two constants emboding the machine epsilon is present */
1768  /* twoesp is set to twice the machine epsilon while u26 is set */
1769  /* to 26 times the machine epsilon */
1770 
1771  /* data twoeps, u26/4.4d-16, 5.72d-15/ *** */
1772  twoeps = 2.f * std::numeric_limits< C_FLOAT64 >::epsilon();
1773  //twoeps = 4.4e-16;
1774  u26 = twoeps * 13.f;
1775 
1776  /* check input parameters */
1777  if (*neqn < 1)
1778  {
1779  goto L10;
1780  }
1781 
1782  if (*relerr < 0. || *abserr < 0.)
1783  {
1784  goto L10;
1785  }
1786 
1787  mflag = abs(*iflag);
1788 
1789  if (mflag >= 1 && mflag <= 8)
1790  {
1791  goto L20;
1792  }
1793 
1794  /* invalid input */
1795 L10:
1796  *iflag = 8;
1797  return 0;
1798 
1799  /* is this the first call */
1800 L20:
1801 
1802  if (mflag == 1)
1803  {
1804  goto L50;
1805  }
1806 
1807  /* check continuation possibilities */
1808 
1809  if (*t == *tout && *kflag != 3)
1810  {
1811  goto L10;
1812  }
1813 
1814  if (mflag != 2)
1815  {
1816  goto L25;
1817  }
1818 
1819  /* iflag = +2 or -2 */
1820  if (*kflag == 3)
1821  {
1822  goto L45;
1823  }
1824 
1825  if (*init == 0)
1826  {
1827  goto L45;
1828  }
1829 
1830  if (*kflag == 4)
1831  {
1832  goto L40;
1833  }
1834 
1835  if (*kflag == 5 && *abserr == 0.)
1836  {
1837  goto L30;
1838  }
1839 
1840  if (*kflag == 6 && *relerr <= *savre && *abserr <= *savae)
1841  {
1842  goto L30;
1843  }
1844 
1845  goto L50;
1846 
1847  /* iflag = 3,4,5,6,7 or 8 */
1848 L25:
1849 
1850  if (*iflag == 3)
1851  {
1852  goto L45;
1853  }
1854 
1855  if (*iflag == 4)
1856  {
1857  goto L40;
1858  }
1859 
1860  if (*iflag == 5 && *abserr > 0.)
1861  {
1862  goto L45;
1863  }
1864 
1865  /* integration cannot be continued since user did not respond to */
1866  /* the instructions pertaining to iflag=5,6,7 or 8 */
1867 L30:
1868  //s_stop("", (ftnlen)0);
1869  std::cout << "STOP 0 statement executed" << std::endl;
1870  return -1;
1871 
1872  /* reset function evaluation counter */
1873 L40:
1874  *nfe = 0;
1875 
1876  if (mflag == 2)
1877  {
1878  goto L50;
1879  }
1880 
1881  /* reset flag value from previous call */
1882 L45:
1883  *iflag = *jflag;
1884 
1885  if (*kflag == 3)
1886  {
1887  mflag = abs(*iflag);
1888  }
1889 
1890  /* save input iflag and set continuation flag value for subsequent */
1891  /* input checking */
1892 L50:
1893  *jflag = *iflag;
1894  *kflag = 0;
1895 
1896  /* save relerr and abserr for checking input on subsequent calls */
1897  *savre = *relerr;
1898  *savae = *abserr;
1899 
1900  /* restrict relative error tolerance to be at least as large as */
1901  /* 2*eps+remin to avoid limiting precision difficulties arising */
1902  /* from impossible accuracy requests */
1903 
1904  rer = twoeps + remin;
1905 
1906  if (*relerr >= rer)
1907  {
1908  goto L55;
1909  }
1910 
1911  /* relative error tolerance too small */
1912  *relerr = rer;
1913  *iflag = 3;
1914  *kflag = 3;
1915  return 0;
1916 
1917 L55:
1918  dt = *tout - *t;
1919 
1920  if (mflag == 1)
1921  {
1922  goto L60;
1923  }
1924 
1925  if (*init == 0)
1926  {
1927  goto L65;
1928  }
1929 
1930  goto L80;
1931 
1932  /* initialization -- */
1933  /* set initialization completion indicator,init */
1934  /* set indicator for too many output points,kop */
1935  /* evaluate initial derivatives */
1936  /* set counter for function evaluations,nfe */
1937  /* evaluate initial derivatives */
1938  /* set counter for function evaluations,nfe */
1939  /* estimate starting stepsize */
1940 
1941 L60:
1942  *init = 0;
1943  *kop = 0;
1944 
1945  a = *t;
1946  (*f)(neqn, &a, &y[1], &yp[1]);
1947  *nfe = 1;
1948 
1949  if (*t != *tout)
1950  {
1951  goto L65;
1952  }
1953 
1954  *iflag = 2;
1955  return 0;
1956 
1957 L65:
1958  *init = 1;
1959  *h__ = fabs(dt);
1960  toln = 0.f;
1961  i__1 = *neqn;
1962 
1963  for (k = 1; k <= i__1; ++k)
1964  {
1965  tol = *relerr * (d__1 = y[k], fabs(d__1)) + *abserr;
1966 
1967  if (tol <= 0.f)
1968  {
1969  goto L70;
1970  }
1971 
1972  toln = tol;
1973  ypk = (d__1 = yp[k], fabs(d__1));
1974  /* Computing 5th power */
1975  d__1 = *h__, d__2 = d__1, d__1 *= d__1;
1976 
1977  if (ypk * (d__2 * (d__1 * d__1)) > tol)
1978  {
1979  d__3 = tol / ypk;
1980  *h__ = pow(d__3, .2);
1981  }
1982 
1983 L70:
1984  ;
1985  }
1986 
1987  if (toln <= 0.)
1988  {
1989  *h__ = 0.;
1990  }
1991 
1992  /* Computing MAX */
1993  /* Computing MAX */
1994  d__3 = fabs(*t), d__4 = fabs(dt);
1995  d__1 = *h__, d__2 = u26 * std::max(d__3, d__4);
1996  *h__ = std::max(d__1, d__2);
1997  //*jflag = i_sign(2, iflag);
1998  *jflag = *iflag >= 0 ? 2 : -2;
1999 
2000  /* set stepsize for integration in the direction from t to tout */
2001 
2002 L80:
2003  //*h__ = d_sign(h__, &dt);
2004  *h__ = dt >= 0 ? fabs(*h__) : -fabs(*h__);
2005 
2006  /* test to see if rkf45 is being severely impacted by too many */
2007  /* output points */
2008 
2009  if (fabs(*h__) >= fabs(dt) * 2.)
2010  {
2011  ++(*kop);
2012  }
2013 
2014  if (*kop != 100)
2015  {
2016  goto L85;
2017  }
2018 
2019  /* unnecessary frequency of output */
2020  *kop = 0;
2021  *iflag = 7;
2022  return 0;
2023 
2024 L85:
2025 
2026  if (fabs(dt) > u26 * fabs(*t))
2027  {
2028  goto L95;
2029  }
2030 
2031  /* if too close to output point,extrapolate and return */
2032  i__1 = *neqn;
2033 
2034  for (k = 1; k <= i__1; ++k)
2035  {
2036  /* L90: */
2037  y[k] += dt * yp[k];
2038  }
2039 
2040  a = *tout;
2041  (*f)(neqn, &a, &y[1], &yp[1]);
2042  ++(*nfe);
2043 
2044  std::cout << "Step too Small" << std::endl;
2045  mUseStateRecord = false;
2046  goto L300;
2047 
2048  /* initialize output point indicator */
2049 
2050 L95:
2051  output = 0;
2052 
2053  /* to avoid premature underflow in the error tolerance function, */
2054  /* scale the error tolerances */
2055 
2056  scale = 2. / *relerr;
2057  ae = scale * *abserr;
2058 
2059  /* step by step integration */
2060 
2061 L100:
2062  hfaild = 0;
2063 
2064  /* set smallest allowable stepsize */
2065  hmin = u26 * fabs(*t);
2066 
2067  /* adjust stepsize if necessary to hit the output point. */
2068  /* look ahead two steps to avoid drastic changes in the stepsize and */
2069  /* thus lessen the impact of output points on the code. */
2070 
2071  dt = *tout - *t;
2072 
2073  if (fabs(dt) >= fabs(*h__) * 2.)
2074  {
2075  goto L200;
2076  }
2077 
2078  if (fabs(dt) > fabs(*h__))
2079  {
2080  goto L150;
2081  }
2082 
2083  /* the next successful step will complete the integration to the */
2084  /* output point */
2085 
2086  output = 1;
2087  *h__ = dt;
2088  goto L200;
2089 
2090 L150:
2091  *h__ = dt * .5;
2092 
2093  /* core integrator for taking a single step */
2094 
2095  /* the tolerances have been scaled to avoid premature underflow in */
2096  /* computing the error tolerance function et. */
2097  /* to avoid problems with zero crossings,relative error is measured */
2098  /* using the average of the magnitudes of the solution at the */
2099  /* beginning and end of a step. */
2100  /* the error estimate formula has been grouped to control loss of */
2101  /* significance. */
2102  /* to distinguish the various arguments, h is not permitted */
2103  /* to become smaller than 26 units of roundoff in t. */
2104  /* practical limits on the change in the stepsize are enforced to */
2105  /* smooth the stepsize selection process and to avoid excessive */
2106  /* chattering on problems having discontinuities. */
2107  /* to prevent unnecessary failures, the code uses 9/10 the stepsize */
2108  /* it estimates will succeed. */
2109  /* after a step failure, the stepsize is not allowed to increase for */
2110  /* the next attempted step. this makes the code more efficient on */
2111  /* problems having discontinuities and more effective in general */
2112  /* since local extrapolation is being used and extra caution seems */
2113  /* warranted. */
2114 
2115  /* test number of derivative function evaluations. */
2116  /* if okay,try to advance the integration from t to t+h */
2117 
2118 L200:
2119 
2120  if (*nfe <= maxnfe)
2121  {
2122  goto L220;
2123  }
2124 
2125  /* too much work */
2126  *iflag = 4;
2127  *kflag = 4;
2128  return 0;
2129 
2130  /* advance an approximate solution over one step of length h */
2131 
2132 L220:
2133  fehl_((pEvalF)f, neqn, &y[1], t, h__, &yp[1], &f1[1], &f2[1],
2134  &f3[1], &f4[1], &f5[1], &f1[1], &yrcd[1]);
2135  *nfe += 5;
2136  mUseStateRecord = true;
2137 
2138  /* compute and test allowable tolerances versus local error estimates*/
2139  /* and remove scaling of tolerances. note that relative error is */
2140  /* measured with respect to the average of the magnitudes of the */
2141  /* solution at the beginning and end of the step. */
2142 
2143  eeoet = 0.;
2144  i__1 = *neqn;
2145 
2146  for (k = 1; k <= i__1; ++k)
2147  {
2148  et = (d__1 = y[k], fabs(d__1)) + (d__2 = f1[k], fabs(d__2)) + ae;
2149 
2150  if (et > 0.)
2151  {
2152  goto L240;
2153  }
2154 
2155  /* inappropriate error tolerance */
2156  *iflag = 5;
2157  return 0;
2158 
2159 L240:
2160  ee = (d__1 = yp[k] * -2090. + (f3[k] * 21970. - f4[k] * 15048.) + (f2[k] * 22528. - f5[k] * 27360.), fabs(d__1));
2161  /* L250: */
2162  /* Computing MAX */
2163  d__1 = eeoet, d__2 = ee / et;
2164  eeoet = std::max(d__1, d__2);
2165  }
2166 
2167  esttol = fabs(*h__) * eeoet * scale / 752400.;
2168 
2169  if (esttol <= 1.)
2170  {
2171  goto L260;
2172  }
2173 
2174  /* unsuccessful step */
2175  /* reduce the stepsize , try again */
2176  /* the decrease is limited to a factor of 1/10 */
2177 
2178  hfaild = 1;
2179  output = 0;
2180  s = .1;
2181 
2182  if (esttol < 59049.)
2183  {
2184  s = .9 / pow(esttol, .2);
2185  }
2186 
2187  *h__ = s * *h__;
2188 
2189  if (fabs(*h__) > hmin)
2190  {
2191  goto L200;
2192  }
2193 
2194  /* requested error unattainable at smallest allowable stepsize */
2195  *iflag = 6;
2196  *kflag = 6;
2197  return 0;
2198 
2199  /* successful step */
2200  /* store solution at t+h */
2201  /* and evaluate derivatives there */
2202 
2203 L260:
2204  *t += *h__;
2205  i__1 = *neqn;
2206 
2207  for (k = 1; k <= i__1; ++k)
2208  {
2209  /* L270: */
2210  y[k] = f1[k];
2211  }
2212 
2213  a = *t;
2214  (*f)(neqn, &a, &y[1], &yp[1]);
2215  ++(*nfe);
2216 
2217  /* choose next stepsize */
2218  /* the increase is limited to a factor of 5 */
2219  /* if step failure has just occurred, next */
2220  /* stepsize is not allowed to increase */
2221 
2222  s = 5.;
2223 
2224  if (esttol > 1.889568e-4)
2225  {
2226  s = .9 / pow(esttol, .2);
2227  }
2228 
2229  if (hfaild)
2230  {
2231  s = std::min(s, 1.);
2232  }
2233 
2234  /* Computing MAX */
2235  d__2 = s * fabs(*h__);
2236  d__1 = std::max(d__2, hmin);
2237  //*h__ = d_sign(&d__1, h__);
2238  *h__ = *h__ >= 0 ? fabs(d__1) : -fabs(d__1);
2239 
2240  /* end of core integrator */
2241 
2242  /* should we take another step */
2243  if (output)
2244  {
2245  goto L300;
2246  }
2247 
2248  if (*iflag > 0)
2249  {
2250  goto L100;
2251  }
2252 
2253  /* integration successfully completed */
2254 
2255  /* one-step mode */
2256  *iflag = -2;
2257  return 0;
2258 
2259  /* interval mode */
2260 L300:
2261  *t = *tout;
2262  *iflag = 2;
2263  return 0;
2264 } /* rkfs_ */
2265 
2266 C_INT CHybridMethodODE45::fehl_(pEvalF f, const C_INT *neqn, double *y,
2267  double *t, double *h__,
2268  double *yp, double *f1,
2269  double *f2, double *f3,
2270  double *f4, double *f5,
2271  double *s, double *yrcd)
2272 {
2273  /* System generated locals */
2274  C_INT i__1;
2275  double d__1;
2276 
2277  /* Local variables */
2278  static C_INT k;
2279  static double ch;
2280  static C_INT base;
2281 
2282  /* fehlberg fourth-fifth order runge-kutta method */
2283 
2284  /* fehl integrates a system of neqn first order */
2285  /* ordinary differential equations of the form */
2286  /* dy(i)/dt=f(t,y(1),---,y(neqn)) */
2287  /* where the initial values y(i) and the initial derivatives */
2288  /* yp(i) are specified at the starting point t. fehl advances */
2289  /* the solution over the fixed step h and returns */
2290  /* the fifth order (sixth order accurate locally) solution */
2291  /* approximation at t+h in array s(i). */
2292  /* f1,---,f5 are arrays of dimension neqn which are needed */
2293  /* for internal storage. */
2294  /* the formulas have been grouped to control loss of significance. */
2295  /* fehl should be called with an h not smaller than 13 units of */
2296  /* roundoff in t so that the various independent arguments can be */
2297  /* distinguished. */
2298 
2299  /* Parameter adjustments */
2300  --yrcd;
2301  --s;
2302  --f5;
2303  --f4;
2304  --f3;
2305  --f2;
2306  --f1;
2307  --yp;
2308  --y;
2309 
2310  /* Function Body */
2311 
2312  //~~~~~~~~(1)~~~~~~~~
2313  ch = *h__ / 4.;
2314  i__1 = *neqn;
2315 
2316  for (k = 1; k <= i__1; ++k)
2317  {
2318  f5[k] = y[k] + ch * yp[k];
2319  }
2320 
2321  d__1 = *t + ch;
2322  (*f)(neqn, &d__1, &f5[1], &f1[1]);
2323 
2324  /* Record trcd and yrcd at t+h/4 */
2325  base = 1;
2326  yrcd[base] = *t + ch;
2327  i__1 = *neqn;
2328 
2329  for (k = 1; k <= i__1; ++k)
2330  {
2331  yrcd[base + k] = f5[k];
2332  }
2333 
2334  //~~~~~~~~(2)~~~~~~~~
2335  ch = *h__ * 3. / 32.;
2336  i__1 = *neqn;
2337 
2338  for (k = 1; k <= i__1; ++k)
2339  {
2340  f5[k] = y[k] + ch * (yp[k] + f1[k] * 3.);
2341  }
2342 
2343  d__1 = *t + *h__ * 3. / 8.;
2344  (*f)(neqn, &d__1, &f5[1], &f2[1]);
2345 
2346  /* Record trcd and yrcd at t+3h/8 */
2347  base = *neqn + 2;
2348  yrcd[base] = *t + *h__ * 3. / 8.;
2349  i__1 = *neqn;
2350 
2351  for (k = 1; k <= i__1; ++k)
2352  {
2353  yrcd[base + k] = f5[k];
2354  }
2355 
2356  //~~~~~~~~(3)~~~~~~~~
2357  ch = *h__ / 2197.;
2358  i__1 = *neqn;
2359 
2360  for (k = 1; k <= i__1; ++k)
2361  {
2362  f5[k] = y[k] + ch * (yp[k] * 1932. + (f2[k] * 7296.
2363  - f1[k] * 7200.));
2364  }
2365 
2366  d__1 = *t + *h__ * 12. / 13.;
2367  (*f)(neqn, &d__1, &f5[1], &f3[1]);
2368 
2369  /* Record trcd and yrcd at t+12h/13 */
2370  base = (*neqn + 1) * 3 + 1;
2371  yrcd[base] = *t + *h__ * 12. / 13.;
2372  i__1 = *neqn;
2373 
2374  for (k = 1; k <= i__1; ++k)
2375  {
2376  yrcd[base + k] = f5[k];
2377  }
2378 
2379  //~~~~~~~~(4)~~~~~~~~
2380  ch = *h__ / 4104.;
2381  i__1 = *neqn;
2382 
2383  for (k = 1; k <= i__1; ++k)
2384  {
2385  f5[k] = y[k] + ch * (yp[k] * 8341. - f3[k] * 845.
2386  + (f2[k] * 29440. - f1[k] * 32832.));
2387  }
2388 
2389  d__1 = *t + *h__;
2390  (*f)(neqn, &d__1, &f5[1], &f4[1]);
2391 
2392  //~~~~~~~~(5)~~~~~~~~
2393  ch = *h__ / 20520.;
2394  i__1 = *neqn;
2395 
2396  for (k = 1; k <= i__1; ++k)
2397  {
2398  f1[k] = y[k] + ch * (yp[k] * -6080.
2399  + (f3[k] * 9295. - f4[k] * 5643.)
2400  + (f1[k] * 41040. - f2[k] * 28352.));
2401  }
2402 
2403  d__1 = *t + *h__ / 2.;
2404  (*f)(neqn, &d__1, &f1[1], &f5[1]);
2405 
2406  /* Record trcd and yrcd at t+h/2 */
2407  //base = (*neqn + 1 << 1) + 1;
2408  base = (*neqn + 1) * 2 + 1;
2409  yrcd[base] = *t + *h__ / 2.;
2410  i__1 = *neqn;
2411 
2412  for (k = 1; k <= i__1; ++k)
2413  {
2414  yrcd[base + k] = f1[k];
2415  }
2416 
2417  //~~~~~~~~(6)~~~~~~~~
2418  /* compute approximate solution at t+h */
2419  ch = *h__ / 7618050.;
2420  i__1 = *neqn;
2421 
2422  for (k = 1; k <= i__1; ++k)
2423  {
2424  s[k] = y[k] + ch * (yp[k] * 902880.
2425  + (f3[k] * 3855735. - f4[k] * 1371249.)
2426  + (f2[k] * 3953664. + f5[k] * 277020.));
2427  }
2428 
2429  return 0;
2430 } /* fehl_ */
2431 
2432 //========Help Function========
2433 /**
2434  * Test the model if it is proper to perform stochastic simulations on.
2435  * Several properties are tested (e.g. integer stoichometry, all reactions
2436  * take place in one compartment only, irreversibility...).
2437  *
2438  * @return 0, if everything is ok; <0, if an error occured.
2439  */
2441 {
2443  CMatrix <C_FLOAT64> mStoi = model->getStoi();
2444  C_INT32 multInt;
2445  size_t i, j, numReactions = mpReactions->size();
2446  C_FLOAT64 multFloat;
2447  // size_t metabSize = mpMetabolites->size();
2448 
2449  for (i = 0; i < numReactions; i++) // for every reaction
2450  {
2451  // TEST getCompartmentNumber() == 1
2452  if ((*mpReactions)[i]->getCompartmentNumber() != 1) return - 1;
2453 
2454  // TEST isReversible() == 0
2455  if ((*mpReactions)[i]->isReversible() != 0) return - 2;
2456 
2457  // TEST integer stoichometry
2458  // Iterate through each the metabolites
2459  // juergen: the number of rows of mStoi equals the number of non-fixed metabs!
2460  // for (j=0; i<metabSize; j++)
2461  for (j = 0; j < mStoi.numRows(); j++)
2462  {
2463  multFloat = mStoi[j][i];
2464  multInt = static_cast<C_INT32>(floor(multFloat + 0.5)); // +0.5 to get a rounding out of the static_cast to int!
2465 
2466  if ((multFloat - multInt) > INT_EPSILON) return - 3; // INT_EPSILON in CHybridMethod.h
2467  }
2468  }
2469 
2470  return 1; // Model is appropriate for hybrid simulation
2471 }
2472 
2473 /**
2474  * Prints out various data on standard output for debugging purposes.
2475  */
2476 void CHybridMethodODE45::outputDebug(std::ostream & os, size_t level)
2477 {
2478  size_t i, j;
2479  std::set< size_t >::iterator iter, iterEnd;
2480 
2481  os << "outputDebug(" << level << ") *********************************************** BEGIN" << std::endl;
2482 
2483  switch (level)
2484  {
2485  case 0: // Everything !!!
2486  os << " Name: " << CCopasiParameter::getObjectName() << std::endl;
2487  os << "current time: " << mpCurrentState->getTime() << std::endl;
2488  os << "mNumVariableMetabs: " << mNumVariableMetabs << std::endl;
2489  os << "mMaxSteps: " << mMaxSteps << std::endl;
2490  os << "mMaxBalance: " << mMaxBalance << std::endl;
2491  //os << "mMaxIntBeforeStep: " << mMaxIntBeforeStep << std::endl;
2492  os << "mpReactions.size(): " << mpReactions->size() << std::endl;
2493 
2494  for (i = 0; i < mpReactions->size(); i++)
2495  os << *(*mpReactions)[i] << std::endl;
2496 
2497  os << "mpMetabolites.size(): " << mpMetabolites->size() << std::endl;
2498 
2499  for (i = 0; i < mpMetabolites->size(); i++)
2500  os << *(*mpMetabolites)[i] << std::endl;
2501 
2502  os << "mStoi: " << std::endl;
2503 
2504  for (i = 0; i < mStoi.numRows(); i++)
2505  {
2506  for (j = 0; j < mStoi.numCols(); j++)
2507  os << mStoi[i][j] << " ";
2508 
2509  os << std::endl;
2510  }
2511 
2512  os << "temp: ";
2513 
2514  for (i = 0; i < mNumVariableMetabs; i++)
2515  os << temp[i] << " ";
2516 
2517  os << std::endl;
2518 
2519  os << "mReactionFlags: " << std::endl;
2520 
2521  for (i = 0; i < mLocalBalances.size(); i++)
2522  os << mReactionFlags[i];
2523 
2524  //os << "mFirstReactionFlag: " << std::endl;
2525  // if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
2526 
2527  os << "mMetabFlags: " << std::endl;
2528 
2529  for (i = 0; i < mMetabFlags.size(); i++)
2530  {
2531  if (mMetabFlags[i].mFlag == SLOW)
2532  os << "SLOW ";
2533  else
2534  os << "FAST ";
2535  }
2536 
2537  os << std::endl;
2538  os << "mLocalBalances: " << std::endl;
2539 
2540  /*
2541  for (i = 0; i < mLocalBalances.size(); i++)
2542  {
2543  for (j = 0; j < mLocalBalances[i].size(); j++)
2544  os << mLocalBalances[i][j];
2545 
2546  os << std::endl;
2547  }
2548 
2549  os << "mLocalSubstrates: " << std::endl;
2550 
2551  for (i = 0; i < mLocalSubstrates.size(); i++)
2552  {
2553  for (j = 0; j < mLocalSubstrates[i].size(); j++)
2554  os << mLocalSubstrates[i][j];
2555 
2556  os << std::endl;
2557  }
2558  */
2559  //os << "mLowerStochLimit: " << mLowerStochLimit << std::endl;
2560  //os << "mUpperStochLimit: " << mUpperStochLimit << std::endl;
2561  //deprecated: os << "mOutputCounter: " << mOutputCounter << endl;
2562  //os << "mPartitioningInterval: " << mPartitioningInterval << std::endl;
2563  //os << "mStepsAfterPartitionSystem: " << mStepsAfterPartitionSystem << std::endl;
2564  //os << "mStepsize: " << mStepsize << std::endl;
2565  os << "mMetab2React: " << std::endl;
2566 
2567  for (i = 0; i < mMetab2React.size(); i++)
2568  {
2569  os << i << ": ";
2570 
2571  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; ++iter)
2572  os << *iter << " ";
2573 
2574  os << std::endl;
2575  }
2576 
2577  os << "mAmu: " << std::endl;
2578 
2579  for (i = 0; i < mpReactions->size(); i++)
2580  os << mAmu[i] << " ";
2581 
2582  os << std::endl;
2583  os << "mAmuOld: " << std::endl;
2584 
2585  for (i = 0; i < mpReactions->size(); i++)
2586  os << mAmuOld[i] << " ";
2587 
2588  os << std::endl;
2589  os << "mUpdateSet: " << std::endl;
2590 
2591  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
2592  os << *iter;
2593 
2594  os << std::endl;
2595  os << "mpRandomGenerator: " << mpRandomGenerator << std::endl;
2596  os << "mDG: " << std::endl << mDG;
2597  os << "mPQ: " << std::endl << mPQ;
2598  os << "Particle numbers: " << std::endl;
2599 
2600  for (i = 0; i < mpMetabolites->size(); i++)
2601  {
2602  os << (*mpMetabolites)[i]->getValue() << " ";
2603  }
2604 
2605  os << std::endl;
2606  break;
2607 
2608  case 1: // Variable values only
2609  os << "current time: " << mpCurrentState->getTime() << std::endl;
2610  /*
2611  case 1:
2612  os << "mTime: " << mpCurrentState->getTime() << std::endl;
2613  os << "oldState: ";
2614  for (i = 0; i < mDim; i++)
2615  os << oldState[i] << " ";
2616  os << std::endl;
2617  os << "x: ";
2618  for (i = 0; i < mDim; i++)
2619  os << x[i] << " ";
2620  os << std::endl;
2621  os << "y: ";
2622  for (i = 0; i < mDim; i++)
2623  os << y[i] << " ";
2624  os << std::endl;
2625  os << "increment: ";
2626  for (i = 0; i < mDim; i++)
2627  os << increment[i] << " ";
2628  os << std::endl;*/
2629  os << "temp: ";
2630 
2631  for (i = 0; i < mNumVariableMetabs; i++)
2632  os << temp[i] << " ";
2633 
2634  os << std::endl;
2635 
2636  os << "mReactionFlags: " << std::endl;
2637 
2638  for (i = 0; i < mLocalBalances.size(); i++)
2639  os << mReactionFlags[i];
2640 
2641  // os << "mFirstReactionFlag: " << std::endl;
2642  //if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
2643 
2644  os << "mMetabFlags: " << std::endl;
2645 
2646  for (i = 0; i < mMetabFlags.size(); i++)
2647  {
2648  if (mMetabFlags[i].mFlag == SLOW)
2649  os << "SLOW ";
2650  else
2651  os << "FAST ";
2652  }
2653 
2654  os << std::endl;
2655  os << "mAmu: " << std::endl;
2656 
2657  for (i = 0; i < mpReactions->size(); i++)
2658  os << mAmu[i] << " ";
2659 
2660  os << std::endl;
2661  os << "mAmuOld: " << std::endl;
2662 
2663  for (i = 0; i < mpReactions->size(); i++)
2664  os << mAmuOld[i] << " ";
2665 
2666  os << std::endl;
2667  os << "mUpdateSet: " << std::endl;
2668 
2669  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
2670  os << *iter;
2671 
2672  os << std::endl;
2673  os << "mPQ: " << std::endl << mPQ;
2674  os << "Particle numbers: " << std::endl;
2675 
2676  for (i = 0; i < mpMetabolites->size(); i++)
2677  {
2678  os << (*mpMetabolites)[i]->getValue() << " ";
2679  }
2680 
2681  os << std::endl;
2682  break;
2683 
2684  case 2:
2685  break;
2686 
2687  default:
2688  ;
2689  }
2690 
2691  os << "outputDebug(" << level << ") ************************************************* END" << std::endl;
2692  return;
2693 }
2694 
2695 /**
2696  * Prints out data on standard output. Deprecated.
2697  */
2699 {
2700  std::cout << "============Boolean============" << std::endl;
2701 
2702  if (mDoCorrection)
2703  std::cout << "mDoCorrection: Yes" << std::endl;
2704  else
2705  std::cout << "mDoCorrection: No" << std::endl;
2706 
2707  if (mReducedModel)
2708  std::cout << "mReducedModel: Yes" << std::endl;
2709  else
2710  std::cout << "mReducedModel: No" << std::endl;
2711 
2712  std::cout << std::endl;
2713 
2714  std::cout << "============Metab============" << std::endl;
2715 
2716  std::cout << "~~~~mNumVariableMetabs:~~~~ " << mNumVariableMetabs << std::endl;
2717  std::cout << "~~~~mFirstMetabIndex:~~~~ " << mFirstMetabIndex << std::endl;
2718  std::cout << "~~~~mpMetabolites:~~~~" << std::endl;
2719 
2720  for (size_t i = 0; i < mpMetabolites->size(); i++)
2721  {
2722  std::cout << "metab #" << i << " name: " << (*mpMetabolites)[i]->getObjectDisplayName() << std::endl;
2723  std::cout << "value pointer: " << (*mpMetabolites)[i]->getValuePointer() << std::endl;
2724  std::cout << "value: " << *((double *)(*mpMetabolites)[i]->getValuePointer()) << std::endl;
2725  }
2726 
2727  std::cout << std::endl;
2728 
2729  std::cout << "~~~~mMetabFlags:~~~~ " << std::endl;
2730 
2731  for (size_t i = 0; i < mMetabFlags.size(); ++i)
2732  {
2733  std::cout << "metab #" << i << std::endl;
2734  std::cout << "mFlag: " << mMetabFlags[i].mFlag << std::endl;
2735  std::cout << "mFastReactions: ";
2736  std::set<size_t>::iterator it = mMetabFlags[i].mFastReactions.begin();
2737  std::set<size_t>::iterator endIt = mMetabFlags[i].mFastReactions.end();
2738 
2739  for (; it != endIt; it++)
2740  std::cout << *it << " ";
2741 
2742  std::cout << std::endl;
2743  }
2744 
2745  std::cout << std::endl;
2746 
2747  std::cout << "============Reaction============" << std::endl;
2748 
2749  std::cout << "~~~~mNumReactions:~~~~ " << mNumReactions << std::endl;
2750 
2751  for (size_t i = 0; i < mNumReactions; ++i)
2752  {
2753  std::cout << "Reaction #: " << i << " Flag: " << mReactionFlags[i] << std::endl;
2754  }
2755 
2756  std::cout << "~~~~mLocalBalances:~~~~ " << std::endl;
2757 
2758  for (size_t i = 0; i < mLocalBalances.size(); ++i)
2759  {
2760  std::cout << "Reaction: " << i << std::endl;
2761 
2762  for (size_t j = 0; j < mLocalBalances[i].size(); ++j)
2763  {
2764  std::cout << "Index: " << mLocalBalances[i][j].mIndex << std::endl;
2765  std::cout << "mMultiplicity: " << mLocalBalances[i][j].mMultiplicity << std::endl;
2766  std::cout << "mpMetablite: " << mLocalBalances[i][j].mpMetabolite << std::endl;
2767  }
2768  }
2769 
2770  std::cout << "~~~~mLocalSubstrates:~~~~ " << std::endl;
2771 
2772  for (size_t i = 0; i < mLocalSubstrates.size(); ++i)
2773  {
2774  std::cout << "Reaction: " << i << std::endl;
2775 
2776  for (size_t j = 0; j < mLocalSubstrates[i].size(); ++j)
2777  {
2778  std::cout << "Index: " << mLocalSubstrates[i][j].mIndex << std::endl;
2779  std::cout << "mMultiplicity: " << mLocalSubstrates[i][j].mMultiplicity << std::endl;
2780  std::cout << "mpMetablite: " << mLocalSubstrates[i][j].mpMetabolite << std::endl;
2781  }
2782  }
2783 
2784  std::cout << "~~~~mMetab2React:~~~~ " << std::endl;
2785 
2786  for (size_t i = 0; i < mMetab2React.size(); i++)
2787  {
2788  std::cout << "metab #: " << i << std::endl;
2789  std::set<size_t>::iterator it = mMetab2React[i].begin();
2790  std::set<size_t>::iterator endIt = mMetab2React[i].end();
2791  std::cout << "React: ";
2792 
2793  for (; it != endIt; it++)
2794  std::cout << *it << " ";
2795 
2796  std::cout << std::endl;
2797  }
2798 
2799  std::cout << std::endl;
2800  std::cout << "~~~~mReactAffect:~~~~ " << std::endl;
2801 
2802  for (size_t i = 0; i < mReactAffect.size(); i++)
2803  {
2804  std::cout << "react #: " << i << std::endl;
2805  std::set<size_t>::iterator it = mReactAffect[i].begin();
2806  std::set<size_t>::iterator endIt = mReactAffect[i].end();
2807  std::cout << "affect: ";
2808 
2809  for (; it != endIt; it++)
2810  std::cout << *it << " ";
2811 
2812  std::cout << std::endl;
2813  }
2814 
2815  if (mHasStoiReaction)
2816  std::cout << "mHasStoiReaction: Yes" << std::endl;
2817  else
2818  std::cout << "mHasStoiReaction: No" << std::endl;
2819 
2820  if (mHasDetermReaction)
2821  std::cout << "mHasDetermReaction: Yes" << std::endl;
2822  else
2823  std::cout << "mHasDetermReaction: No" << std::endl;
2824 
2825  getchar();
2826  return;
2827 }
2828 
2829 /**
2830  * Print State Data for Debug
2831  */
2833 {
2834  const C_FLOAT64 time = mpState->getTime();
2835  std::cout << "**State Output**" << std::endl;
2836  std::cout << "Time: " << time << std::endl;
2837 
2838  std::cout << "Indep #: " << mpState->getNumIndependent() << " Id: ";
2839  const C_FLOAT64 *pIt = mpState->beginIndependent();
2840  const C_FLOAT64 *pEnd = mpState->endIndependent();
2841 
2842  for (; pIt != pEnd; pIt++)
2843  std::cout << *pIt << " ";
2844 
2845  std::cout << std::endl;
2846 
2847  std::cout << "Dep #: " << mpState->getNumDependent() << " Id: ";
2848  pIt = mpState->beginDependent();
2849  pEnd = mpState->endDependent();
2850 
2851  for (; pIt != pEnd; pIt++)
2852  std::cout << *pIt << " ";
2853 
2854  std::cout << std::endl;
2855 
2856  std::cout << "Fix #: " << mpState->getNumIndependent() << " Id: ";
2857  pIt = mpState->beginFixed();
2858  pEnd = mpState->endFixed();
2859 
2860  for (; pIt != pEnd; pIt++)
2861  std::cout << *pIt << " ";
2862 
2863  std::cout << std::endl;
2864 
2865  getchar();
2866  return;
2867 }
2868 
2870 {
2871  size_t i, imax = pModel->getNumModelValues();
2872 
2873  for (i = 0; i < imax; ++i)
2874  {
2875  if (pModel->getModelValues()[i]->getStatus() == CModelEntity::ASSIGNMENT)
2876  if (pModel->getModelValues()[i]->isUsed())
2877  {
2878  //used assignment found
2879  return true;
2880  }
2881  }
2882 
2883  imax = pModel->getNumMetabs();
2884 
2885  for (i = 0; i < imax; ++i)
2886  {
2887  if (pModel->getMetabolites()[i]->getStatus() == CModelEntity::ASSIGNMENT)
2888  if (pModel->getMetabolites()[i]->isUsed())
2889  {
2890  //used assignment found
2891  return true;
2892  }
2893  }
2894 
2895  imax = pModel->getCompartments().size();
2896 
2897  for (i = 0; i < imax; ++i)
2898  {
2899  if (pModel->getCompartments()[i]->getStatus() == CModelEntity::ASSIGNMENT)
2900  if (pModel->getCompartments()[i]->isUsed())
2901  {
2902  //used assignment found
2903  return true;
2904  }
2905  }
2906 
2907  return false;
2908 }
#define CONTINUE
std::set< size_t > mUpdateSet
#define C_INT
Definition: copasi.h:115
CVector< C_FLOAT64 > temp
CCopasiObject * getParticleFluxReference()
Definition: CReaction.cpp:207
CCopasiVectorN< CEvent > & getEvents()
Definition: CModel.cpp:1110
const CCopasiVectorNS< CReaction > * mpReactions
const CCopasiVector< CMetab > & getMetabolites() const
Definition: CModel.cpp:1051
const std::string & getObjectName() const
const CCopasiVectorN< CModelValue > & getModelValues() const
Definition: CModel.cpp:1060
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)
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
#define HYBRID
virtual size_t size() const
void addDependent(const size_t &node, const size_t &dependent)
CDependencyGraph mDG
#define FAST
void recordState(const C_FLOAT64 time, const C_FLOAT64 *y)
#define HAS_EVENT
size_t getNumMetabs() const
Definition: CModel.cpp:1118
std::vector< std::set< size_t > > mMetab2React
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
virtual size_t numRows() const
Definition: CMatrix.h:138
#define USE_RANDOM_SEED
Definition: CHybridMethod.h:67
void getStochTimeAndIndex(C_FLOAT64 &ds, size_t &rIndex)
void outputState(const CState *pS)
#define INT_EPSILON
Definition: CHybridMethod.h:59
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)
CTrajectoryProblem * mpProblem
CMatrix< C_FLOAT64 > mStoi
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
virtual void setValue(const C_FLOAT64 &value)
#define RANDOM_SEED
Definition: CHybridMethod.h:68
#define C_INVALID_INDEX
Definition: copasi.h:222
unsigned C_INT32 mRandomSeed
C_INT32 checkModel(CModel *model)
#define INTERP_RECORD_NUM
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
C_FLOAT64 getKey(const size_t index) const
Definition: CMetab.h:178
virtual bool isValidProblem(const CCopasiProblem *pProblem)
C_FLOAT64 getDefaultAtol(const CModel *pModel) const
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
#define MAX_STEPS
Definition: CHybridMethod.h:58
void updatePriorityQueue(size_t rIndex, C_FLOAT64 time)
CVector< C_FLOAT64 > mYdot
virtual bool elevateChildren()
C_FLOAT64 * endDependent()
Definition: CState.cpp:331
C_FLOAT64 * endIndependent()
Definition: CState.cpp:329
std::vector< std::vector< CHybridODE45Balance > > mLocalSubstrates
C_FLOAT64 * endFixed()
Definition: CState.cpp:333
bool removeParameter(const std::string &name)
#define NEW_STEP
void fireReaction(size_t rIndex)
size_t getNumDependent() const
Definition: CState.cpp:344
#define MCTrajectoryMethod
std::set< std::string > * getDependsOn(size_t rIndex)
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
const C_FLOAT64 & getQuantity2NumberFactor() const
Definition: CModel.cpp:2354
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
void setState(const CState &state)
Definition: CModel.cpp:1785
CVector< C_FLOAT64 > mAmu
#define HAS_ERR
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
const C_FLOAT64 & getDuration() const
C_FLOAT64 * getArray() const
const Value & getValue() const
std::set< std::string > * getAffects(size_t rIndex)
#define SLOW
bool setValue(const std::string &name, const CType &value)
unsigned C_INT32 * pUINT
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
size_t getNumIndependent() const
Definition: CState.cpp:342
std::vector< std::vector< CHybridODE45Balance > > mLocalBalances
const C_FLOAT64 & getTime() const
CInterpolation * mpInterpolation
CCopasiParameter * getParameter(const std::string &name)
CIndexedPriorityQueue mPQ
CVector< C_INT > mIWork
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
virtual Status step(const double &deltaT)
const C_FLOAT64 & getTime() const
Definition: CModel.cpp:1190
CStateRecord * getInterpolationState()
bool dependsOn(DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const
CHybridMethodODE45 * pMethod
CVector< C_FLOAT64 > mStateRecord
long int flag
Definition: f2c.h:52
void setTime(const C_FLOAT64 &time)
Definition: CModel.cpp:1187
void refreshConcentration()
Definition: CMetab.cpp:281
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
void updateNode(const size_t index, const C_FLOAT64 key)
static bool modelHasAssignments(const CModel *pModel)
const C_FLOAT64 & getValue() const
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
virtual size_t getIndex(const CCopasiObject *pObject) const
CVector< C_FLOAT64 > mDWork
void updateTauMu(size_t rIndex, C_FLOAT64 time)
CCopasiVector< CMetab > * mpMetabolites
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
void(* pEvalF)(const C_INT *, const double *, const double *, double *)
#define STOCHASTIC
size_t getNumModelValues() const
Definition: CModel.cpp:1139
bool addParameter(const CCopasiParameter &parameter)
#define abs(x)
Definition: f2c.h:173
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
std::vector< std::set< size_t > > mReactAffect
C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
std::ostringstream mErrorMsg
Definition: CModel.h:50
C_FLOAT64 * beginFixed()
Definition: CState.cpp:332
const CState & getState() const
Definition: CModel.cpp:1771
const CModelEntity::Status & getStatus() const
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
void integrateDeterministicPart(C_FLOAT64 ds)
#define DETERMINISTIC
virtual size_t numCols() const
Definition: CMatrix.h:144
C_FLOAT64 * beginDependent()
Definition: CState.cpp:330
CHybridMethodODE45(const CCopasiContainer *pParent=NULL)
CModel * getModel() const
#define min(a, b)
Definition: f2c.h:175
void initializeIndexPointer(const size_t numberOfReactions)
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
CVector< size_t > mReactionFlags
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
void calculateAmu(size_t rIndex)
void initMethod(C_FLOAT64 time)
#define REACH_END_TIME
std::set< size_t > mCalculateSet
CStateRecord * mpEventState
void outputDebug(std::ostream &os, size_t level)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
CCopasiObject * getValueReference() const
#define max(a, b)
Definition: f2c.h:176