COPASI API  4.16.103
CHybridMethodLSODA.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2006 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * CHybridMethodLSODA
17  *
18  * This class implements an hybrid algorithm for the simulation of a
19  * biochemical system over time.
20  *
21  * File name: CHybridMethodLSODA.cpp
22  * Author: Juergen Pahle
23  * Email: juergen.pahle@eml-r.villa-bosch.de
24  *
25  * Last change: 09, May 2006
26  *
27  * (C) European Media Lab 2003.
28  */
29 
30 /* DEFINE ********************************************************************/
31 
32 #ifdef WIN32
33 #if _MSC_VER < 1600
34 #define min _cpp_min
35 #define max _cpp_max
36 #endif // _MSC_VER
37 #endif // WIN32
38 
39 #include <limits.h>
40 #include <iterator>
41 
42 #include "copasi.h"
43 
44 #include "CHybridMethodLSODA.h"
45 #include "CTrajectoryProblem.h"
46 #include "model/CModel.h"
47 #include "model/CMetab.h"
48 #include "model/CReaction.h"
49 #include "model/CState.h"
50 #include "model/CChemEq.h"
51 #include "model/CChemEqElement.h"
52 #include "model/CCompartment.h"
54 #include "utilities/CMatrix.h"
57 #include "utilities/CVersion.h"
59 
60 /**
61  * Default constructor.
62  */
64  CTrajectoryMethod(CCopasiMethod::hybridLSODA, pParent)
65 {
66  assert((void *) &mData == (void *) &mData.dim);
67 
68  mData.pMethod = this;
69 
72 }
73 
75  const CCopasiContainer * pParent):
76  CTrajectoryMethod(src, pParent)
77 {
78  assert((void *) &mData == (void *) &mData.dim);
79 
80  mData.pMethod = this;
81 
84 }
85 
86 /**
87  * Destructor.
88  */
90 {
91  cleanup();
93 }
94 
96 {
97  CCopasiParameter *pParm;
98 
99  assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) MAX_STEPS);
102 
103  assertParameter("Partitioning Interval", CCopasiParameter::UINT, (unsigned C_INT32) PARTITIONING_INTERVAL);
104  assertParameter("Use Random Seed", CCopasiParameter::BOOL, (bool) USE_RANDOM_SEED);
105  assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) RANDOM_SEED);
106 
107  // Check whether we have a method with the old parameter names
108  if ((pParm = getParameter("HYBRID.MaxSteps")) != NULL)
109  {
110  setValue("Max Internal Steps", *pParm->getValue().pUINT);
111  removeParameter("HYBRID.MaxSteps");
112 
113  if ((pParm = getParameter("HYBRID.LowerStochLimit")) != NULL)
114  {
115  setValue("Lower Limit", *pParm->getValue().pDOUBLE);
116  removeParameter("HYBRID.LowerStochLimit");
117  }
118 
119  if ((pParm = getParameter("HYBRID.UpperStochLimit")) != NULL)
120  {
121  setValue("Upper Limit", *pParm->getValue().pDOUBLE);
122  removeParameter("HYBRID.UpperStochLimit");
123  }
124 
125  if ((pParm = getParameter("HYBRID.PartitioningInterval")) != NULL)
126  {
127  setValue("Partitioning Interval", *pParm->getValue().pUINT);
128  removeParameter("HYBRID.PartitioningInterval");
129  }
130 
131  if ((pParm = getParameter("UseRandomSeed")) != NULL)
132  {
133  setValue("Use Random Seed", *pParm->getValue().pBOOL);
134  removeParameter("UseRandomSeed");
135  }
136 
137  if ((pParm = getParameter("")) != NULL)
138  {
139  setValue("Random Seed", *pParm->getValue().pUINT);
140  removeParameter("");
141  }
142  }
143 
144  /* LSODA ****************************************************************************/
145 
146  addParameter("Partitioning Stepsize",
148 
149  addParameter("Relative Tolerance",
151 
152  addParameter("Use Default Absolute Tolerance",
153  CCopasiParameter::BOOL, (bool) true);
154 
155  addParameter("Absolute Tolerance",
157 
158  addParameter("Max Internal Steps (LSODA)",
159  CCopasiParameter::UINT, (unsigned C_INT32) 10000);
160 
161  // These parameters are no longer supported.
162  removeParameter("Adams Max Order");
163  removeParameter("BDF Max Order");
164 }
165 
167 {
169  return true;
170 }
171 
172 /**
173  * Creates a HybridMethodLSODA adequate for the problem.
174  */
176 {
177  C_INT32 result = 1;
178  CHybridMethodLSODA * method = NULL;
179 
180  switch (result)
181  {
182  /* case - 3: // non-integer stoichometry
183  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 1);
184  break;
185  case - 2: // reversible reaction exists
186  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 2);
187  break;
188 
189  case - 1: // more than one compartment involved
190  CCopasiMessage(CCopasiMessage::ERROR, MCTrajectoryMethod + 3);
191  break;*/
192  case 1:
193  default:
194  // Everything alright: Hybrid simulation possible
195  method = new CHybridNextReactionLSODAMethod();
196  break;
197  }
198 
199  return method;
200 }
201 
203 {
204  // write the current state to the model
205  // mpProblem->getModel()->setState(mpCurrentState); // is that correct?
206 
207  // check for possible overflows
208  size_t i;
209  size_t imax;
210 
211  // :TODO: Bug 774: This assumes that the number of variable metabs is the number
212  // of metabs determined by reaction. In addition they are expected at the beginning of the
213  // MetabolitesX which is not the case if we have metabolites of type ODE.
214  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++)
215  if (mpProblem->getModel()->getMetabolitesX()[i]->getValue() >= mMaxIntBeforeStep)
216  {
217  // throw exception or something like that
218  }
219 
220  // do several steps
221  C_FLOAT64 time = mpProblem->getModel()->getTime();
222  C_FLOAT64 endTime = time + deltaT;
223 
224  for (i = 0; ((i < mMaxSteps) && (time < endTime)); i++)
225  {
226  time = doSingleStep(time, endTime);
227  }
228 
229  //mpCurrentState->setTime(time);
230 
231  if ((i >= mMaxSteps) && (!mMaxStepsReached))
232  {
233  mMaxStepsReached = true; //only report this message once
234  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.");
235  }
236 
237  // get back the particle numbers
238 
239  /* Set the variable metabolites */
241 
242  for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, Dbl++)
243  *Dbl = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
244 
245  //the task expects the result in mpCurrentState
247 
248  mpCurrentState->setTime(time);
249 
250  return NORMAL;
251 }
252 
253 void CHybridMethodLSODA::start(const CState * initialState)
254 {
255  //*mpCurrentState = *initialState;
256 
257  mpProblem->getModel()->setState(*initialState);
258  mRestartLSODA = true;
259 
261 
263  assert(mpModel);
264 
265  if (mpModel->getModelType() == CModel::deterministic)
266  mDoCorrection = true;
267  else
268  mDoCorrection = false;
269 
271 
272  mFirstMetabIndex = mpModel->getStateTemplate().getIndex(mpModel->getMetabolitesX()[0]);
273 
274 //mpProblem->getModel()->setState(*mpCurrentState);
275 
276  // call init of the simulation method, can be overloaded in derived classes
277  // mDim[1] = (C_INT) (void *) this;
279 
280  return;
281 }
282 
283 /* PROTECTED METHODS *********************************************************/
284 
285 /**
286  * Initializes the solver and sets the model to be used.
287  *
288  * @param model A reference to an instance of a CModel
289  */
291 {
293  mAmu.clear();
294  mAmu.resize(mpReactions->size());
295  mAmuOld.clear();
296  mAmuOld.resize(mpReactions->size());
298  //mNumVariableMetabs = mpModel->getNumVariableMetabs(); // ind + dep metabs, without fixed metabs
300 
301  temp.clear();
302  temp.resize(mNumVariableMetabs);
303  //currentState.clear();
304  //currentState.resize(mNumVariableMetabs);
305 
306  /* get configuration data */
307  mMaxSteps = * getValue("Max Internal Steps").pUINT;
308  mLowerStochLimit = * getValue("Lower Limit").pDOUBLE;
309  mUpperStochLimit = * getValue("Upper Limit").pDOUBLE;
310  mStepsize = * getValue("Partitioning Stepsize").pDOUBLE;
311  mPartitioningInterval = * getValue("Partitioning Interval").pUINT;
312  mUseRandomSeed = * getValue("Use Random Seed").pBOOL;
313  mRandomSeed = * getValue("Random Seed").pUINT;
314 
316 
317  mStoi = mpModel->getStoi();
319  mUpdateSet.clear();
320 
321  setupBalances(); // initialize mLocalBalances and mLocalSubstrates (has to be called first!)
322  setupDependencyGraph(); // initialize mDG
323  setupMetab2React(); // initialize mMetab2React
324  setupPartition(); // initialize mReactionFlags
325  setupPriorityQueue(start_time); // initialize mPQ
326 
327  mMaxStepsReached = false;
328 
329  /* CONFIGURE LSODA ***********************************************************************/
330 
332 
334 
335  mRtol = * getValue("Relative Tolerance").pUDOUBLE;
336 
337  mDefaultAtol = * getValue("Use Default Absolute Tolerance").pBOOL;
338 
339  if (mDefaultAtol)
340  {
342  setValue("Absolute Tolerance", mAtol);
343  }
344  else
345  mAtol = * getValue("Absolute Tolerance").pUDOUBLE;
346 
347  mDWork.resize(22 + mData.dim * std::max<C_INT>(16, mData.dim + 9));
348  mDWork[4] = mDWork[5] = mDWork[6] = mDWork[7] = mDWork[8] = mDWork[9] = 0.0;
349  mIWork.resize(20 + mData.dim);
350 
351  mIWork[4] = mIWork[6] = mIWork[9] = 0;
352  mIWork[5] = * getValue("Max Internal Steps (LSODA)").pUINT;
353  mIWork[7] = 12;
354  mIWork[8] = 5;
355 
356  /***********************************************************************************/
357 
358  return;
359 }
360 
362 {
363 
364  if (!pModel) return 1.0e009;
365 
366  const CCopasiVectorNS< CCompartment > & Compartment = pModel->getCompartments();
367  size_t i, imax;
368 
370 
371  for (i = 0, imax = Compartment.size(); i < imax; i++)
372  if (Compartment[i]->getValue() < Volume) Volume = Compartment[i]->getValue();
373 
374  if (Volume == std::numeric_limits< C_FLOAT64 >::max()) return 1.0e009;
375 
376  return Volume * pModel->getQuantity2NumberFactor() * 1.e-12;
377 }
378 
379 /**
380  * Cleans up memory, etc.
381  */
383 {
384  delete mpRandomGenerator;
385  mpRandomGenerator = NULL;
386  mpModel = NULL;
387  return;
388 }
389 
390 /* DETERMINISTIC STUFF *******************************************************/
391 
392 void CHybridMethodLSODA::EvalF(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot)
393 {static_cast<Data *>((void *) n)->pMethod->evalF(t, y, ydot);}
394 
395 void CHybridMethodLSODA::evalF(const C_FLOAT64 * t, const C_FLOAT64 * /* y */, C_FLOAT64 * ydot)
396 {
397  size_t i;
398 
399  CModel * pModel = mpProblem->getModel();
400 
401  mpState->setTime(*t);
402 
403  pModel->setState(*mpState);
404 
405  pModel->updateSimulatedValues(false);
406 
408 
409  for (i = 0; i < mNumVariableMetabs; i++)
410  ydot[i] = temp[i];
411 
412  return;
413 }
414 
415 /**
416  * Integrates the deterministic reactions of the system over the specified
417  * time interval.
418  *
419  * @param ds A C_FLOAT64 specifying the stepsize.
420  */
421 
423 {
424  CHybridLSODAStochFlag * react = NULL;
425 
426  /* RESET LSODA **************************************************************/
427 
428  if (mRestartLSODA)
429  {
430  mRestartLSODA = false;
431  mLsodaStatus = 1;
432  }
433 
434  mState = 1;
435  mJType = 2;
436 
437  mErrorMsg.str("");
439 
440  C_INT one = 1;
441  C_INT DSize = (C_INT) mDWork.size();
442  C_INT ISize = (C_INT) mIWork.size();
443 
444  //mpState = new CState (mpProblem->getModel()->getState());
446 
447  mY = mpState->beginIndependent();// + mFirstMetabIndex - 1;
448  //for (i = 0, imax = mpProblem->getModel()->getNumVariableMetabs(); i < imax; i++, mY++)
449  // *mY = mpProblem->getModel()->getMetabolitesX()[i]->getValue();
450 
451  mTime = mpState->getTime();
452 
453  C_FLOAT64 EndTime = mTime + deltaT;
454 
455  C_FLOAT64 tdist , d__1, d__2, w0;
456  tdist = fabs(deltaT);
457  d__1 = fabs(mTime), d__2 = fabs(EndTime);
458  w0 = std::max(d__1, d__2);
459 
460  if (tdist < std::numeric_limits< C_FLOAT64 >::epsilon() * 2. * w0) //just do nothing
461  {
462  //mTime = mTime + deltaT;
463  //mpState->setTime(mTime);
464  //*mpCurrentState = *mpState;
465  mpProblem->getModel()->setTime(EndTime);
466 
467  return;
468  }
469 
470  if (!mData.dim) //just do nothing if there are no variables
471  {
472  //mTime = mTime + deltaT;
473  mpProblem->getModel()->setTime(EndTime);
474  //mpState->setTime(mTime);
475  //*mpCurrentState = *mpState;
476 
477  return;
478  }
479 
480  mLSODA(&EvalF , // 1. evaluate F
481  &mData.dim, // 2. number of variables
482  mY , // 3. the array of current concentrations
483  &mTime , // 4. the current time
484  &EndTime , // 5. the final time
485  &one , // 6. scalar error control
486  &mRtol , // 7. relative tolerance array
487  &mAtol , // 8. absolute tolerance array
488  &mState , // 9. output by overshoot & interpolatation
489  &mLsodaStatus , // 10. the state control variable
490  &one , // 11. futher options (one)
491  mDWork.array() , // 12. the double work array
492  &DSize , // 13. the double work array size
493  mIWork.array() , // 14. the int work array
494  &ISize , // 15. the int work array size
495  NULL , // 16. evaluate J (not given)
496  &mJType); // 17. the type of jacobian calculate (2)
497 
498  if (mLsodaStatus == -1) mLsodaStatus = 2;
499 
500  if ((mLsodaStatus != 1) && (mLsodaStatus != 2) && (mLsodaStatus != -1))
501  {
503  }
504 
506  //*mpCurrentState = *mpState;
508 
509  /* size_t i;
510 
511  for (i = 0; i < mNumVariableMetabs; i++)
512  {
513  (*mpMetabolites)[i]->setValue(mY[i]);
514  (*mpMetabolites)[i]->refreshConcentration();
515  } */
516 
517  mpModel->updateSimulatedValues(false); //for assignments
518 
519  //pdelete(mpState);
520 
521  // find the set union of all reactions, which depend on one of the deterministic reactions. The propensities of the stochastic reactions in this set union will be updated later in the method updatePriorityQueue().
522 
523  for (react = mFirstReactionFlag; react != NULL; react = react->mpNext)
524  {
525  const std::set <size_t> & dependents = mDG.getDependents(react->mIndex);
526  std::copy(dependents.begin(), dependents.end(),
527  std::inserter(mUpdateSet, mUpdateSet.begin()));
528  }
529 
530  return;
531 }
532 
533 /**
534  * Calculates the derivative of the system and writes it into the vector
535  * deriv. Length of deriv must be mNumVariableMetabs.
536  * CAUTION: Only deterministic reactions are taken into account. That is,
537  * this is only the derivative of the deterministic part of the system.
538  *
539  * @param deriv A vector reference of length mNumVariableMetabs, into
540  * which the derivative is written
541  */
542 void CHybridMethodLSODA::calculateDerivative(std::vector <C_FLOAT64> & deriv)
543 {
544  size_t i;
545  C_INT32 bal = 0;
547 
548  // Calculate all the needed kinetic functions of the deterministic reactions
549  for (j = mFirstReactionFlag; j != NULL; j = j->mpNext)
550  {
551  (*mpReactions)[j->mIndex]->calculateParticleFlux();
552  }
553 
554  // For each metabolite add up the contributions of the det. reactions
555  // the number of rows in mStoi equals the number of non-fixed metabolites!
556  // for (i=0; i<mNumVariableMetabs; i++)
557  for (i = 0; i < mStoi.numRows(); i++)
558  {
559  deriv[i] = 0.0;
560 
561  for (j = mFirstReactionFlag; j != NULL; j = j->mpNext)
562  {
563  // juergen: +0.5 to get a rounding out of the static_cast
564  bal = static_cast<C_INT32>(floor(mStoi[i][j->mIndex] + 0.5));
565  deriv[i] += bal * (*mpReactions)[j->mIndex]->getParticleFlux(); // balance * flux;
566  }
567  }
568 
569  /*
570  for (; i < mNumVariableMetabs; i++) deriv[i] = 0.0; // important to get a correct deriv vector, because mStoi doesn't cover fixed metabolites
571  */
572  return;
573 }
574 
575 /* STOCHASTIC STUFF **********************************************************/
576 
577 /**
578  * Find the reaction index and the reaction time of the stochastic (!)
579  * reaction with the lowest reaction time.
580  *
581  * @param ds A reference to a C_FLOAT64. The putative reaction time for the
582  * first stochastic reaction is written into this variable.
583  * @param rIndex A reference to a size_t. The index of the first
584  * stochastic reaction is written into this variable.
585  */
587 {
588  ds = mPQ.topKey();
589  rIndex = mPQ.topIndex();
590  return;
591 }
592 
593 /**
594  * Executes the specified reaction in the system once.
595  *
596  * @param rIndex A size_t specifying the index of the reaction, which
597  * will be fired.
598  * @param time The current time
599  */
601 {
602  // Change the particle numbers according to which step took place.
603  // First, get the vector of balances in the reaction we've got.
604  // (This vector expresses the number change of each metabolite
605  // in the reaction.) Then step through each balance, using its
606  // multiplicity to calculate a new value for the associated
607  // metabolite. Finally, update the metabolite.
608 
609  size_t i;
610  C_FLOAT64 newNumber;
611 
612  CMetab * pMetab;
613 
614  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
615  {
616  pMetab = mLocalBalances[rIndex][i].mpMetabolite;
617  newNumber = pMetab->getValue() + mLocalBalances[rIndex][i].mMultiplicity;
618 
619  pMetab->setValue(newNumber);
620  pMetab->refreshConcentration();
621  }
622 
623  // insert all dependent reactions into the mUpdateSet
624  const std::set <size_t> & dependents = mDG.getDependents(rIndex);
625  std::copy(dependents.begin(), dependents.end(),
626  std::inserter(mUpdateSet, mUpdateSet.begin()));
627 
628  mRestartLSODA = true;
629 
630  return;
631 }
632 
633 /**
634  * Updates the priority queue.
635  *
636  * @param rIndex A size_t giving the index of the fired reaction (-1, if no
637  * stochastic reaction has fired)
638  * @param time A C_FLOAT64 holding the current time
639  */
640 
642 {
643  C_FLOAT64 newTime;
644  size_t index;
645  std::set <size_t>::iterator iter, iterEnd;
646 
647  //if the model contains assignments we use a less efficient loop over all (stochastic) reactions to capture all changes
648  // we do not know the exact dependencies. TODO: this should be changed later in order to get a more efficient update scheme
649  if (mHasAssignments)
650  {
652 
653  for (index = 0; index < mpReactions->size(); index++)
654  {
655  if (mReactionFlags[index].mpPrev == NULL) // Reaction is stochastic!
656  {
657  mAmuOld[index] = mAmu[index];
658  calculateAmu(index);
659 
660  if (mAmuOld[index] != mAmu[index])
661  if (index != rIndex) updateTauMu(index, time);
662  }
663  }
664  }
665  else
666  {
667  // iterate through the set of affected reactions and update the stochastic ones in the priority queue
668  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
669  {
670  if (mReactionFlags[*iter].mpPrev == NULL) // reaction is stochastic!
671  {
672  index = *iter;
673  mAmuOld[index] = mAmu[index];
674  calculateAmu(index);
675 
676  if (*iter != rIndex) updateTauMu(index, time);
677  }
678  }
679  }
680 
681  // draw new random number and update the reaction just fired
682  if ((rIndex != C_INVALID_INDEX) && (mReactionFlags[rIndex].mpPrev == NULL))
683  {
684  // reaction is stochastic
685  newTime = time + generateReactionTime(rIndex);
686  mPQ.updateNode(rIndex, newTime);
687  }
688 
689  // empty the mUpdateSet
690  mUpdateSet.clear();
691  return;
692 }
694 {
695  if (mAmu[rIndex] == 0) return std::numeric_limits<C_FLOAT64>::infinity();
696 
698  return - 1 * log(rand2) / mAmu[rIndex];
699 }
700 
701 /**
702  * Calculates an amu value for a given reaction.
703  *
704  * @param rIndex A size_t specifying the reaction to be updated
705  */
707 {
708  if (!mDoCorrection)
709  {
710  mAmu[rIndex] = (*mpReactions)[rIndex]->calculateParticleFlux();
711  return;
712  }
713 
714  // We need the product of the cmu and hmu for this step.
715  // We calculate this in one go, as there are fewer steps to
716  // perform and we eliminate some possible rounding errors.
717  C_FLOAT64 amu = 1; // initially
718  //size_t total_substrates = 0;
719  C_INT32 num_ident = 0;
720  C_INT32 number = 0;
721  C_INT32 lower_bound;
722  // substrate_factor - The substrates, raised to their multiplicities,
723  // multiplied with one another. If there are, e.g. m substrates of type m,
724  // and n of type N, then substrate_factor = M^m * N^n.
725  C_FLOAT64 substrate_factor = 1;
726  // First, find the reaction associated with this index.
727  // Keep a pointer to this.
728  // Iterate through each substrate in the reaction
729  const std::vector<CHybridLSODABalance> & substrates = mLocalSubstrates[rIndex];
730 
731  int flag = 0;
732 
733  for (size_t i = 0; i < substrates.size(); i++)
734  {
735  num_ident = substrates[i].mMultiplicity;
736 
737  if (num_ident > 1)
738  {
739  flag = 1;
740  number = static_cast<C_INT32>((*mpMetabolites)[substrates[i].mIndex]->getValue());
741  lower_bound = number - num_ident;
742  substrate_factor = substrate_factor * pow((double) number, (int) num_ident - 1); //optimization
743 
744  number--; // optimization
745 
746  while (number > lower_bound)
747  {
748  amu *= number;
749  number--;
750  }
751  }
752  }
753 
754  if ((amu == 0) || (substrate_factor == 0)) // at least one substrate particle number is zero
755  {
756  mAmu[rIndex] = 0;
757  return;
758  }
759 
760  // rate_factor is the rate function divided by substrate_factor.
761  // It would be more efficient if this was generated directly, since in effect we
762  // are multiplying and then dividing by the same thing (substrate_factor)!
763  C_FLOAT64 rate_factor = (*mpReactions)[rIndex]->calculateParticleFlux();
764 
765  if (flag)
766  {
767  amu *= rate_factor / substrate_factor;;
768  mAmu[rIndex] = amu;
769  }
770  else
771  {
772  mAmu[rIndex] = rate_factor;
773  }
774 
775  return;
776 
777  // a more efficient way to calculate mass action kinetics could be included
778 }
779 
780 /**
781  * Updates the putative reaction time of a stochastic reaction in the
782  * priority queue. The corresponding amu and amu_old must be set prior to
783  * the call of this method.
784  *
785  * @param rIndex A size_t specifying the index of the reaction
786  */
788 {
789  C_FLOAT64 newTime;
790 
791  // One must make sure that the calculation yields reasonable results even in the cases where mAmu=0 or mAmuOld=0 or both =0. Therefore mAmuOld=0 is checked. If mAmuOld equals 0, then a new random number has to be drawn, because tau equals inf and the stoch. information is lost. If both values equal 0, then tau should remain inf and the update of the queue can be skipped!
792 
793  if (mAmuOld[rIndex] == 0.0)
794  {
795  if (mAmu[rIndex] != 0.0)
796  {
797  newTime = time + generateReactionTime(rIndex);
798  mPQ.updateNode(rIndex, newTime);
799  }
800  }
801  else
802  {
803  newTime = time + (mAmuOld[rIndex] / mAmu[rIndex]) * (mPQ.getKey(rIndex) - time);
804  mPQ.updateNode(rIndex, newTime);
805  }
806 
807  return;
808 }
809 
810 /* TESTING THE MODEL AND SETTING UP THINGS ***********************************/
811 
812 /**
813  * Test the model if it is proper to perform stochastic simulations on.
814  * Several properties are tested (e.g. integer stoichometry, all reactions
815  * take place in one compartment only, irreversibility...).
816  *
817  * @return 0, if everything is ok; <0, if an error occured.
818  */
820 {
822  CMatrix <C_FLOAT64> mStoi = model->getStoi();
823  C_INT32 multInt;
824  size_t i, j, numReactions = mpReactions->size();
825  C_FLOAT64 multFloat;
826  // size_t metabSize = mpMetabolites->size();
827 
828  for (i = 0; i < numReactions; i++) // for every reaction
829  {
830  // TEST getCompartmentNumber() == 1
831  if ((*mpReactions)[i]->getCompartmentNumber() != 1) return - 1;
832 
833  // TEST isReversible() == 0
834  if ((*mpReactions)[i]->isReversible() != 0) return - 2;
835 
836  // TEST integer stoichometry
837  // Iterate through each the metabolites
838  // juergen: the number of rows of mStoi equals the number of non-fixed metabs!
839  // for (j=0; i<metabSize; j++)
840  for (j = 0; j < mStoi.numRows(); j++)
841  {
842  multFloat = mStoi[j][i];
843  multInt = static_cast<C_INT32>(floor(multFloat + 0.5)); // +0.5 to get a rounding out of the static_cast to int!
844 
845  if ((multFloat - multInt) > INT_EPSILON) return - 3; // INT_EPSILON in CHybridMethod.h
846  }
847  }
848 
849  return 1; // Model is appropriate for hybrid simulation
850 }
851 
852 /**
853  * Sets up an internal representation of the balances for each reaction.
854  * This is done in order to be able to deal with fixed metabolites and
855  * to avoid a time consuming search for the indices of metabolites in the
856  * model.
857  */
859 {
860  size_t i, j;
861  CHybridLSODABalance newElement;
862  C_INT32 maxBalance = 0;
863  size_t numReactions;
864 
865  numReactions = mpReactions->size();
866  mLocalBalances.clear();
867  mLocalBalances.resize(numReactions);
868  mLocalSubstrates.clear();
869  mLocalSubstrates.resize(numReactions);
870 
871  for (i = 0; i < numReactions; i++)
872  {
873  const CCopasiVector <CChemEqElement> * balances =
874  &(*mpReactions)[i]->getChemEq().getBalances();
875 
876  for (j = 0; j < balances->size(); j++)
877  {
878  newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
879  newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
880  // + 0.5 to get a rounding out of the static_cast to C_INT32!
881  newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
882 
883  if ((newElement.mpMetabolite->getStatus()) != CModelEntity::FIXED)
884  {
885  if (newElement.mMultiplicity > maxBalance) maxBalance = newElement.mMultiplicity;
886 
887  mLocalBalances[i].push_back(newElement); // element is copied for the push_back
888  }
889  }
890 
891  balances = &(*mpReactions)[i]->getChemEq().getSubstrates();
892 
893  for (j = 0; j < balances->size(); j++)
894  {
895  newElement.mpMetabolite = const_cast < CMetab* >((*balances)[j]->getMetabolite());
896  newElement.mIndex = mpModel->getMetabolitesX().getIndex(newElement.mpMetabolite);
897  // + 0.5 to get a rounding out of the static_cast to C_INT32!
898  newElement.mMultiplicity = static_cast<C_INT32>(floor((*balances)[j]->getMultiplicity() + 0.5));
899 
900  mLocalSubstrates[i].push_back(newElement); // element is copied for the push_back
901  }
902  }
903 
904  mMaxBalance = maxBalance;
905  mMaxIntBeforeStep = INT_MAX - 1 - mMaxSteps * mMaxBalance;
906 
907  return;
908 }
909 
910 /**
911  * Sets up the dependency graph.
912  */
914 {
915  mDG.clear();
916  std::vector< std::set<std::string>* > DependsOn;
917  std::vector< std::set<std::string>* > Affects;
918  size_t numReactions = mpReactions->size();
919  size_t i, j;
920 
921  // Do for each reaction:
922  for (i = 0; i < numReactions; i++)
923  {
924  // Get the set of metabolites which affect the value of amu for this
925  // reaction i.e. the set on which amu depends. This may be more than
926  // the set of substrates, since the kinetics can involve other
927  // reactants, e.g. catalysts. We thus need to step through the
928  // rate function and pick out every reactant which can vary.
929  DependsOn.push_back(getDependsOn(i));
930  // Get the set of metabolites which are affected when this reaction takes place
931  Affects.push_back(getAffects(i));
932  }
933 
934  // For each possible pair of reactions i and j, if the intersection of
935  // Affects(i) with DependsOn(j) is non-empty, add a dependency edge from i to j.
936  for (i = 0; i < numReactions; i++)
937  {
938  for (j = 0; j < numReactions; j++)
939  {
940  // Determine whether the intersection of these two sets is non-empty
941  // Could also do this with set_intersection generic algorithm, but that
942  // would require operator<() to be defined on the set elements.
943 
944  std::set<std::string>::iterator iter = Affects[i]->begin();
945 
946  for (; iter != Affects[i]->end(); iter++)
947  {
948  if (DependsOn[j]->count(*iter))
949  {
950  // The set intersection is non-empty
951  mDG.addDependent(i, j);
952  break;
953  }
954  }
955  }
956 
957  // Ensure that self edges are included
958  //mDG.addDependent(i, i);
959  }
960 
961  // Delete the memory allocated in getDependsOn() and getAffects()
962  // since this is allocated in other functions.
963  for (i = 0; i < numReactions; i++)
964  {
965  delete DependsOn[i];
966  delete Affects[i];
967  }
968 
969  return;
970 }
971 
972 /**
973  * Creates for each metabolite a set of reaction indices. If the metabolite
974  * participates in a reaction as substrate or product (that means:
975  * balance != 0) this reaction is added to the corresponding set.
976  */
978 {
979  size_t i, j;
980  size_t metaboliteIndex;
981 
982  // Resize mMetab2React and create an initial set for each metabolite
983  mMetab2React.clear();
984  mMetab2React.resize(mpMetabolites->size());
985 
986  // Iterate over all reactions
987  for (i = 0; i < mLocalBalances.size(); i++)
988  {
989  // Get the set of metabolites which take part in this reaction
990  for (j = 0; j < mLocalBalances[i].size(); j++)
991  {
992  // find metaboliteIndex and insert the reaction into the set
993  metaboliteIndex = mLocalBalances[i][j].mIndex;
994  mMetab2React[metaboliteIndex].insert(i);
995  }
996  }
997 
998  return;
999 }
1000 
1001 /**
1002  * Creates for each metabolite a set of reaction indices. If the metabolite
1003  * participates in a reaction as substrate, product or modifier this
1004  * reaction is added to the corresponding set.
1005  */
1007 {
1008  std::vector< std::set<size_t>* > participatesIn;
1009  size_t numReactions = mpReactions->size();
1010  size_t i;
1011 
1012  // Resize mMetab2React and create an initial set for each metabolite
1013  mMetab2React.resize(mpMetabolites->size());
1014 
1015  // Do for each reaction:
1016  for (i = 0; i < numReactions; i++)
1017  {
1018  participatesIn.push_back(getParticipatesIn(i));
1019  }
1020 
1021  // Iterate over all reactions
1022  for (i = 0; i < numReactions; i++)
1023  {
1024  // Get the set of metabolites which take part in this reaction
1025  std::set<size_t>::iterator iter = participatesIn[i]->begin();
1026 
1027  for (; iter != participatesIn[i]->end(); iter++)
1028  mMetab2React[*iter].insert(i);
1029  }
1030 
1031  for (i = 0; i < numReactions; i++)
1032  {
1033  delete participatesIn[i];
1034  }
1035 
1036  return;
1037 }
1038 
1039 /**
1040  * Creates for each metabolite a set of reaction indices. Each reaction is
1041  * dependent on each metabolite resulting in a complete switch.
1042  */
1044 {
1045  size_t i, j;
1046 
1047  // Resize mMetab2React and create an initial set for each metabolite
1048  mMetab2React.resize(mpMetabolites->size());
1049 
1050  // Iterate over all metabolites
1051  for (i = 0; i < mpMetabolites->size(); i++)
1052  {
1053  // Iterate over all reactions
1054  for (j = 0; j < mpReactions->size(); j++)
1055  {
1056  mMetab2React[i].insert(j);
1057  }
1058  }
1059 
1060  return;
1061 }
1062 
1063 /**
1064  * Creates an initial partitioning of the system. Deterministic and
1065  * stochastic reactions are determined. The vector mReactionFlags and
1066  * the vector mMetabFlags are initialized.
1067  */
1069 {
1070  size_t i, j;
1071  CHybridLSODAStochFlag * prevFlag;
1072  C_FLOAT64 averageStochLimit = (mUpperStochLimit + mLowerStochLimit) / 2;
1073 
1074  // initialize vector mMetabFlags
1075  mMetabFlags.clear();
1077 
1078  for (i = 0; i < mMetabFlags.size(); i++)
1079  {
1080  if ((*mpMetabolites)[i]->getValue() < averageStochLimit)
1081  {
1082  mMetabFlags[i] = LOW;
1083  (*mpMetabolites)[i]->setValue(floor((*mpMetabolites)[i]->getValue()));
1084  (*mpMetabolites)[i]->refreshConcentration();
1085  }
1086  else
1087  mMetabFlags[i] = HIGH;
1088  }
1089 
1090  // initialize vector mReactionFlags
1091  mReactionFlags.clear();
1092  mReactionFlags.resize(mLocalBalances.size());
1093 
1094  for (i = 0; i < mLocalBalances.size(); i++)
1095  {
1096  mReactionFlags[i].mIndex = i;
1097  mReactionFlags[i].mValue = 0;
1098 
1099  for (j = 0; j < mLocalBalances[i].size(); j++)
1100  {
1101  if (mMetabFlags[mLocalBalances[i][j].mIndex] == LOW)
1102  {
1103  mReactionFlags[i].mValue++;
1104  }
1105  }
1106  }
1107 
1108  mFirstReactionFlag = NULL;
1109  prevFlag = NULL;
1110 
1111  for (i = 0; i < mLocalBalances.size(); i++)
1112  {
1113  if (mReactionFlags[i].mValue == 0)
1114  {
1115  if (mFirstReactionFlag != NULL)
1116  {
1117  prevFlag->mpNext = &mReactionFlags[i];
1118  mReactionFlags[i].mpPrev = prevFlag;
1119  prevFlag = &mReactionFlags[i];
1120  }
1121  else
1122  {
1124  mReactionFlags[i].mpPrev = &mReactionFlags[i]; // Important to distinguish between stochastic (prev == NULL) and deterministic (prev != NULL) reactions
1125  prevFlag = &mReactionFlags[i];
1126  }
1127  }
1128  else
1129  {
1130  mReactionFlags[i].mpPrev = NULL;
1131  mReactionFlags[i].mpNext = NULL;
1132  }
1133  }
1134 
1135  if (prevFlag != NULL)
1136  {
1137  prevFlag->mpNext = NULL;
1138  }
1139 
1140  return;
1141 }
1142 
1143 /**
1144  * Sets up the priority queue.
1145  *
1146  * @param startTime The time at which the simulation starts.
1147  */
1149 {
1150  size_t i;
1151  C_FLOAT64 time;
1152 
1153  mPQ.clear();
1155 
1156  for (i = 0; i < mpReactions->size(); i++)
1157  {
1158  if (mReactionFlags[i].mpPrev == NULL) // Reaction is stochastic!
1159  {
1160  calculateAmu(i);
1161  time = startTime + generateReactionTime(i);
1162  mPQ.insertStochReaction(i, time);
1163  }
1164  }
1165 
1166  return;
1167 }
1168 
1169 /* HELPER METHODS ************************************************************/
1170 
1171 /**
1172  * Updates the partitioning of the system depending on the particle
1173  * numbers present.
1174  */
1176 {
1177  size_t i;
1178  std::set <size_t>::iterator iter, iterEnd;
1179  C_FLOAT64 key;
1180 
1181  for (i = 0; i < mNumVariableMetabs; i++)
1182  {
1183  if ((mMetabFlags[i] == LOW) && ((*mpMetabolites)[i]->getValue() >= mUpperStochLimit))
1184  {
1185  mRestartLSODA = true;
1186  mMetabFlags[i] = HIGH;
1187 
1188  // go through all corresponding reactions and update flags
1189  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; iter++)
1190  {
1191  mReactionFlags[*iter].mValue--;
1192 
1193  // if reaction gets deterministic, insert it into the linked list of deterministic reactions
1194  if (mReactionFlags[*iter].mValue == 0)
1195  {
1197  mPQ.removeStochReaction(*iter);
1198  }
1199  }
1200  }
1201 
1202  if ((mMetabFlags[i] == HIGH) && ((*mpMetabolites)[i]->getValue() < mLowerStochLimit))
1203  {
1204  mRestartLSODA = true;
1205  mMetabFlags[i] = LOW;
1206  (*mpMetabolites)[i]->setValue(floor((*mpMetabolites)[i]->getValue()));
1207  (*mpMetabolites)[i]->refreshConcentration();
1208 
1209  // go through all corresponding reactions and update flags
1210  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; iter++)
1211  {
1212  if (mReactionFlags[*iter].mValue == 0)
1213  {
1215  /*
1216  mPQ.insertStochReaction(*iter, 1234567.8); // juergen: have to beautify this, number has to be the biggest C_FLOAT64 !!!
1217  */
1218  calculateAmu(*iter);
1219  mAmuOld[*iter] = mAmu[*iter];
1220  //key = mpCurrentState->getTime() + generateReactionTime(*iter); NATALIA
1221 
1222  key = mpProblem->getModel()->getTime() + generateReactionTime(*iter);
1223 
1224  mPQ.insertStochReaction(*iter, key);
1225  }
1226 
1227  mReactionFlags[*iter].mValue++;
1228  }
1229  }
1230  }
1231 
1232  return;
1233 }
1234 
1235 /**
1236  * Inserts a new deterministic reaction into the linked list in the
1237  * array mReactionFlags.
1238  *
1239  * @param rIndex A size_t giving the index of the reaction to be inserted
1240  * into the list of deterministic reactions.
1241  */
1243 {
1244  if (mReactionFlags[rIndex].mpPrev == NULL)
1245  // reaction is stochastic (avoids double insertions)
1246  {
1247  if (mFirstReactionFlag != NULL)
1248  // there are deterministic reactions already
1249  {
1251  mReactionFlags[rIndex].mpNext = mFirstReactionFlag;
1254  }
1255  else
1256  {
1257  // there are no deterministic reactions
1258  // Important to distinguish between stochastic (prev == NULL) and deterministic (prev != NULL) reactions
1259  mReactionFlags[rIndex].mpPrev = &mReactionFlags[rIndex];
1261  }
1262 
1263  mAmu[rIndex] = 0.0;
1264  mAmuOld[rIndex] = 0.0;
1265  }
1266 
1267  return;
1268 }
1269 
1270 /**
1271  * Removes a deterministic reaction from the linked list in the
1272  * array mReactionFlags.
1273  *
1274  * @param rIndex A size_t giving the index of the reaction to be removed
1275  * from the list of deterministic reactions.
1276  */
1278 {
1279  if (mReactionFlags[rIndex].mpPrev != NULL)
1280  // reaction is deterministic
1281  {
1282  if (&mReactionFlags[rIndex] != mFirstReactionFlag)
1283  // reactionFlag is not the first in the linked list
1284  {
1285  mReactionFlags[rIndex].mpPrev->mpNext = mReactionFlags[rIndex].mpNext;
1286 
1287  if (mReactionFlags[rIndex].mpNext != NULL)
1288  mReactionFlags[rIndex].mpNext->mpPrev = mReactionFlags[rIndex].mpPrev;
1289  }
1290  else
1291  // reactionFlag is the first in the linked list
1292  {
1293  if (mReactionFlags[rIndex].mpNext != NULL) // reactionFlag is not the only one in the linked list
1294  {
1295  mFirstReactionFlag = mReactionFlags[rIndex].mpNext;
1297  }
1298  else // reactionFlag is the only one in the linked list
1299  {
1300  mFirstReactionFlag = NULL;
1301  }
1302  }
1303  }
1304 
1305  mReactionFlags[rIndex].mpPrev = NULL;
1306  mReactionFlags[rIndex].mpNext = NULL;
1307  return;
1308 }
1309 
1310 /**
1311  * Gets the set of metabolites on which a given reaction depends.
1312  *
1313  * @param rIndex The index of the reaction being executed.
1314  * @return The set of metabolites depended on.
1315  */
1316 std::set<std::string> *CHybridMethodLSODA::getDependsOn(size_t rIndex)
1317 {
1318  std::set<std::string> *retset = new std::set<std::string>;
1319 
1320  size_t i, imax = (*mpReactions)[rIndex]->getFunctionParameters().size();
1321  size_t j, jmax;
1322 
1323  for (i = 0; i < imax; ++i)
1324  {
1325  if ((*mpReactions)[rIndex]->getFunctionParameters()[i]->getUsage() == CFunctionParameter::PARAMETER)
1326  continue;
1327 
1328  //metablist = (*mpReactions)[rIndex]->getParameterMappingMetab(i);
1329  const std::vector <std::string> & metabKeylist =
1330  (*mpReactions)[rIndex]->getParameterMappings()[i];
1331  jmax = metabKeylist.size();
1332 
1333  for (j = 0; j < jmax; ++j)
1334  {
1335  retset->insert(metabKeylist[j]);
1336  }
1337  }
1338 
1339  return retset;
1340 }
1341 
1342 /**
1343  * Gets the set of metabolites which change number when a given
1344  * reaction is executed.
1345  *
1346  * @param rIndex The index of the reaction being executed.
1347  * @return The set of affected metabolites.
1348  */
1349 std::set<std::string> *CHybridMethodLSODA::getAffects(size_t rIndex)
1350 {
1351  size_t i;
1352  std::set<std::string> *retset = new std::set<std::string>;
1353 
1354  // Get the balances associated with the reaction at this index
1355  // XXX We first get the chemical equation, then the balances, since the getBalances method in CReaction is unimplemented!
1356 
1357  for (i = 0; i < mLocalBalances[rIndex].size(); i++)
1358  {
1359  if (mLocalBalances[rIndex][i].mMultiplicity != 0)
1360  {
1361  retset->insert(mLocalBalances[rIndex][i].mpMetabolite->getKey());
1362  }
1363  }
1364 
1365  return retset;
1366 }
1367 
1368 /**
1369  * Gets the set of metabolites, which participate in the given
1370  * reaction either as substrate, product or modifier.
1371  *
1372  * @param rIndex The index of the reaction being executed.
1373  * @return The set of participating metabolites.
1374  */
1375 std::set<size_t> *CHybridMethodLSODA::getParticipatesIn(size_t /* rIndex */)
1376 {
1377  std::set<size_t> *retset = new std::set<size_t>;
1378  return retset;
1379 }
1380 
1381 /**
1382  * Prints out data on standard output. Deprecated.
1383  */
1384 void CHybridMethodLSODA::outputData(std::ostream & os, C_INT32 mode)
1385 {
1386  static unsigned C_INT32 counter = 0;
1387  size_t i;
1388 
1389  switch (mode)
1390  {
1391  case 0:
1392 
1393  if (mOutputCounter == (counter++))
1394  {
1395  counter = 0;
1396  os << mpCurrentState->getTime() << " : ";
1397 
1398  for (i = 0; i < mpMetabolites->size(); i++)
1399  {
1400  os << (*mpMetabolites)[i]->getValue() << " ";
1401  }
1402 
1403  os << std::endl;
1404  }
1405 
1406  break;
1407 
1408  case 1:
1409  os << mpCurrentState->getTime() << " : ";
1410 
1411  for (i = 0; i < mpMetabolites->size(); i++)
1412  {
1413  os << (*mpMetabolites)[i]->getValue() << " ";
1414  }
1415 
1416  os << std::endl;
1417  break;
1418 
1419  default:
1420  ;
1421  }
1422 
1423  return;
1424 }
1425 
1426 /**
1427  * Prints out various data on standard output for debugging purposes.
1428  */
1429 void CHybridMethodLSODA::outputDebug(std::ostream & os, size_t level)
1430 {
1431  size_t i, j;
1432  std::set< size_t >::iterator iter, iterEnd;
1433 
1434  os << "outputDebug(" << level << ") *********************************************** BEGIN" << std::endl;
1435 
1436  switch (level)
1437  {
1438  case 0: // Everything !!!
1439  os << " Name: " << CCopasiParameter::getObjectName() << std::endl;
1440  os << "current time: " << mpCurrentState->getTime() << std::endl;
1441  os << "mNumVariableMetabs: " << mNumVariableMetabs << std::endl;
1442  os << "mMaxSteps: " << mMaxSteps << std::endl;
1443  os << "mMaxBalance: " << mMaxBalance << std::endl;
1444  os << "mMaxIntBeforeStep: " << mMaxIntBeforeStep << std::endl;
1445  os << "mpReactions.size(): " << mpReactions->size() << std::endl;
1446 
1447  for (i = 0; i < mpReactions->size(); i++)
1448  os << *(*mpReactions)[i] << std::endl;
1449 
1450  os << "mpMetabolites.size(): " << mpMetabolites->size() << std::endl;
1451 
1452  for (i = 0; i < mpMetabolites->size(); i++)
1453  os << *(*mpMetabolites)[i] << std::endl;
1454 
1455  os << "mStoi: " << std::endl;
1456 
1457  for (i = 0; i < mStoi.numRows(); i++)
1458  {
1459  for (j = 0; j < mStoi.numCols(); j++)
1460  os << mStoi[i][j] << " ";
1461 
1462  os << std::endl;
1463  }
1464 
1465  os << "temp: ";
1466 
1467  for (i = 0; i < mNumVariableMetabs; i++)
1468  os << temp[i] << " ";
1469 
1470  os << std::endl;
1471 
1472  os << "mReactionFlags: " << std::endl;
1473 
1474  for (i = 0; i < mLocalBalances.size(); i++)
1475  os << mReactionFlags[i];
1476 
1477  os << "mFirstReactionFlag: " << std::endl;
1478  if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
1479 
1480  os << "mMetabFlags: " << std::endl;
1481 
1482  for (i = 0; i < mMetabFlags.size(); i++)
1483  {
1484  if (mMetabFlags[i] == LOW)
1485  os << "LOW ";
1486  else
1487  os << "HIGH ";
1488  }
1489 
1490  os << std::endl;
1491  os << "mLocalBalances: " << std::endl;
1492 
1493  for (i = 0; i < mLocalBalances.size(); i++)
1494  {
1495  for (j = 0; j < mLocalBalances[i].size(); j++)
1496  os << mLocalBalances[i][j];
1497 
1498  os << std::endl;
1499  }
1500 
1501  os << "mLocalSubstrates: " << std::endl;
1502 
1503  for (i = 0; i < mLocalSubstrates.size(); i++)
1504  {
1505  for (j = 0; j < mLocalSubstrates[i].size(); j++)
1506  os << mLocalSubstrates[i][j];
1507 
1508  os << std::endl;
1509  }
1510 
1511  os << "mLowerStochLimit: " << mLowerStochLimit << std::endl;
1512  os << "mUpperStochLimit: " << mUpperStochLimit << std::endl;
1513  //deprecated: os << "mOutputCounter: " << mOutputCounter << endl;
1514  os << "mPartitioningInterval: " << mPartitioningInterval << std::endl;
1515  os << "mStepsAfterPartitionSystem: " << mStepsAfterPartitionSystem << std::endl;
1516  os << "mStepsize: " << mStepsize << std::endl;
1517  os << "mMetab2React: " << std::endl;
1518 
1519  for (i = 0; i < mMetab2React.size(); i++)
1520  {
1521  os << i << ": ";
1522 
1523  for (iter = mMetab2React[i].begin(), iterEnd = mMetab2React[i].end(); iter != iterEnd; ++iter)
1524  os << *iter << " ";
1525 
1526  os << std::endl;
1527  }
1528 
1529  os << "mAmu: " << std::endl;
1530 
1531  for (i = 0; i < mpReactions->size(); i++)
1532  os << mAmu[i] << " ";
1533 
1534  os << std::endl;
1535  os << "mAmuOld: " << std::endl;
1536 
1537  for (i = 0; i < mpReactions->size(); i++)
1538  os << mAmuOld[i] << " ";
1539 
1540  os << std::endl;
1541  os << "mUpdateSet: " << std::endl;
1542 
1543  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
1544  os << *iter;
1545 
1546  os << std::endl;
1547  os << "mpRandomGenerator: " << mpRandomGenerator << std::endl;
1548  os << "mDG: " << std::endl << mDG;
1549  os << "mPQ: " << std::endl << mPQ;
1550  os << "Particle numbers: " << std::endl;
1551 
1552  for (i = 0; i < mpMetabolites->size(); i++)
1553  {
1554  os << (*mpMetabolites)[i]->getValue() << " ";
1555  }
1556 
1557  os << std::endl;
1558  break;
1559 
1560  case 1: // Variable values only
1561  os << "current time: " << mpCurrentState->getTime() << std::endl;
1562  /*
1563  case 1:
1564  os << "mTime: " << mpCurrentState->getTime() << std::endl;
1565  os << "oldState: ";
1566  for (i = 0; i < mDim; i++)
1567  os << oldState[i] << " ";
1568  os << std::endl;
1569  os << "x: ";
1570  for (i = 0; i < mDim; i++)
1571  os << x[i] << " ";
1572  os << std::endl;
1573  os << "y: ";
1574  for (i = 0; i < mDim; i++)
1575  os << y[i] << " ";
1576  os << std::endl;
1577  os << "increment: ";
1578  for (i = 0; i < mDim; i++)
1579  os << increment[i] << " ";
1580  os << std::endl;*/
1581  os << "temp: ";
1582 
1583  for (i = 0; i < mNumVariableMetabs; i++)
1584  os << temp[i] << " ";
1585 
1586  os << std::endl;
1587 
1588  os << "mReactionFlags: " << std::endl;
1589 
1590  for (i = 0; i < mLocalBalances.size(); i++)
1591  os << mReactionFlags[i];
1592 
1593  os << "mFirstReactionFlag: " << std::endl;
1594  if (mFirstReactionFlag == NULL) os << "NULL" << std::endl; else os << *mFirstReactionFlag;
1595 
1596  os << "mMetabFlags: " << std::endl;
1597 
1598  for (i = 0; i < mMetabFlags.size(); i++)
1599  {
1600  if (mMetabFlags[i] == LOW)
1601  os << "LOW ";
1602  else
1603  os << "HIGH ";
1604  }
1605 
1606  os << std::endl;
1607  os << "mAmu: " << std::endl;
1608 
1609  for (i = 0; i < mpReactions->size(); i++)
1610  os << mAmu[i] << " ";
1611 
1612  os << std::endl;
1613  os << "mAmuOld: " << std::endl;
1614 
1615  for (i = 0; i < mpReactions->size(); i++)
1616  os << mAmuOld[i] << " ";
1617 
1618  os << std::endl;
1619  os << "mUpdateSet: " << std::endl;
1620 
1621  for (iter = mUpdateSet.begin(), iterEnd = mUpdateSet.end(); iter != iterEnd; iter++)
1622  os << *iter;
1623 
1624  os << std::endl;
1625  os << "mPQ: " << std::endl << mPQ;
1626  os << "Particle numbers: " << std::endl;
1627 
1628  for (i = 0; i < mpMetabolites->size(); i++)
1629  {
1630  os << (*mpMetabolites)[i]->getValue() << " ";
1631  }
1632 
1633  os << std::endl;
1634  break;
1635 
1636  case 2:
1637  break;
1638 
1639  default:
1640  ;
1641  }
1642 
1643  os << "outputDebug(" << level << ") ************************************************* END" << std::endl;
1644  return;
1645 }
1646 
1647 std::ostream & operator<<(std::ostream & os, const CHybridLSODAStochFlag & d)
1648 {
1649  os << "CHybridLSODAStochFlag " << std::endl;
1650  os << " mIndex: " << d.mIndex << " mValue: " << d.mValue << std::endl;
1651 
1652  if (d.mpPrev != NULL)
1653  os << " prevIndex: " << d.mpPrev->mIndex << " prevPointer: " << d.mpPrev << std::endl;
1654  else
1655  os << " prevPointer: NULL" << std::endl;
1656 
1657  if (d.mpNext != NULL)
1658  os << " nextIndex: " << d.mpNext->mIndex << " nextPointer: " << d.mpNext << std::endl;
1659  else
1660  os << " nextPointer: NULL" << std::endl;
1661 
1662  return os;
1663 }
1664 
1665 std::ostream & operator<<(std::ostream & os, const CHybridLSODABalance & d)
1666 {
1667  os << "CHybridLSODABalance" << std::endl;
1668  os << " mIndex: " << d.mIndex << " mMultiplicity: " << d.mMultiplicity
1669  << " mpMetabolite: " << d.mpMetabolite << std::endl;
1670  return os;
1671 }
1672 
1673 //virtual
1675 {
1676  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
1677 
1678  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
1679 
1680  if (pTP->getDuration() < 0.0)
1681  {
1682  //back integration not possible
1684  return false;
1685  }
1686 
1687  //check for rules
1688  size_t i, imax = pTP->getModel()->getNumModelValues();
1689 
1690  for (i = 0; i < imax; ++i)
1691  {
1692  if (pTP->getModel()->getModelValues()[i]->getStatus() == CModelEntity::ODE)
1693  {
1694  //ode rule found
1696  return false;
1697  }
1698  }
1699 
1700  imax = pTP->getModel()->getNumMetabs();
1701 
1702  for (i = 0; i < imax; ++i)
1703  {
1704  if (pTP->getModel()->getMetabolites()[i]->getStatus() == CModelEntity::ODE)
1705  {
1706  //ode rule found
1708  return false;
1709  }
1710  }
1711 
1712  imax = pTP->getModel()->getCompartments().size();
1713 
1714  for (i = 0; i < imax; ++i)
1715  {
1716  if (pTP->getModel()->getCompartments()[i]->getStatus() == CModelEntity::ODE)
1717  {
1718  //ode rule found
1720  return false;
1721  }
1722  }
1723 
1724  std::string message = pTP->getModel()->suitableForStochasticSimulation();
1725 
1726  if (message != "")
1727  {
1728  //model not suitable, message describes the problem
1729  CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
1730  return false;
1731  }
1732 
1733  mLowerStochLimit = * getValue("Lower Limit").pDOUBLE;
1734  mUpperStochLimit = * getValue("Upper Limit").pDOUBLE;
1735 
1738 
1739  //events are not supported at the moment
1740  if (pTP->getModel()->getEvents().size() > 0)
1741  {
1743  return false;
1744  }
1745 
1746  return true;
1747 }
1748 
1749 //static
1751 {
1752  size_t i, imax = pModel->getNumModelValues();
1753 
1754  for (i = 0; i < imax; ++i)
1755  {
1756  if (pModel->getModelValues()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1757  if (pModel->getModelValues()[i]->isUsed())
1758  {
1759  //used assignment found
1760  return true;
1761  }
1762  }
1763 
1764  imax = pModel->getNumMetabs();
1765 
1766  for (i = 0; i < imax; ++i)
1767  {
1768  if (pModel->getMetabolites()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1769  if (pModel->getMetabolites()[i]->isUsed())
1770  {
1771  //used assignment found
1772  return true;
1773  }
1774  }
1775 
1776  imax = pModel->getCompartments().size();
1777 
1778  for (i = 0; i < imax; ++i)
1779  {
1780  if (pModel->getCompartments()[i]->getStatus() == CModelEntity::ASSIGNMENT)
1781  if (pModel->getCompartments()[i]->isUsed())
1782  {
1783  //used assignment found
1784  return true;
1785  }
1786  }
1787 
1788  return false;
1789 }
C_FLOAT64 getDefaultAtol(const CModel *pModel) const
#define C_INT
Definition: copasi.h:115
CVector< C_FLOAT64 > mDWork
void outputData(std::ostream &os, C_INT32 mode)
void getStochTimeAndIndex(C_FLOAT64 &ds, size_t &rIndex)
CCopasiVectorN< CEvent > & getEvents()
Definition: CModel.cpp:1110
unsigned C_INT32 mPartitioningInterval
std::vector< std::set< size_t > > mMetab2React
std::vector< C_FLOAT64 > mAmu
const CCopasiVector< CMetab > & getMetabolites() const
Definition: CModel.cpp:1051
const std::string & getObjectName() const
unsigned C_INT32 mMaxSteps
CIndexedPriorityQueue mPQ
const CCopasiVectorN< CModelValue > & getModelValues() const
Definition: CModel.cpp:1060
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
virtual size_t size() const
void addDependent(const size_t &node, const size_t &dependent)
static C_INT32 checkModel(CModel *model)
void updateTauMu(size_t rIndex, C_FLOAT64 time)
size_t getNumMetabs() const
Definition: CModel.cpp:1118
unsigned C_INT32 mStepsAfterPartitionSystem
void setOstream(std::ostream &os)
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
virtual bool elevateChildren()
virtual bool isValidProblem(const CCopasiProblem *pProblem)
#define UPPER_STOCH_LIMIT
Definition: CHybridMethod.h:61
static bool modelHasAssignments(const CModel *pModel)
#define INT_EPSILON
Definition: CHybridMethod.h:59
Definition: CState.h:305
const std::set< size_t > & getDependents(const size_t &node) const
CTrajectoryProblem * mpProblem
#define PARTITIONING_INTERVAL
Definition: CHybridMethod.h:63
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
C_FLOAT64 generateReactionTime(size_t rIndex)
CHybridLSODAStochFlag * mpNext
#define C_INT32
Definition: copasi.h:90
C_FLOAT64 getKey(const size_t index) const
std::ostringstream mErrorMsg
Definition: CMetab.h:178
virtual bool isValidProblem(const CCopasiProblem *pProblem)
CVector< C_INT > mIWork
void updatePriorityQueue(size_t rIndex, C_FLOAT64 time)
unsigned C_INT32 mRandomSeed
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
#define MAX_STEPS
Definition: CHybridMethod.h:58
CCopasiVector< CMetab > * mpMetabolites
CHybridMethodLSODA * pMethod
CMatrix< C_FLOAT64 > mStoi
static CHybridMethodLSODA * createHybridMethodLSODA()
void setupPriorityQueue(C_FLOAT64 startTime=0.0)
bool removeParameter(const std::string &name)
#define MCTrajectoryMethod
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
virtual Status step(const double &deltaT)
void integrateDeterministicPart(C_FLOAT64 ds)
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 getNumVariableMetabs() const
Definition: CModel.cpp:1121
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
void setState(const CState &state)
Definition: CModel.cpp:1785
size_t insertStochReaction(const size_t index, const C_FLOAT64 key)
const C_FLOAT64 & getDuration() const
CVector< C_FLOAT64 > mYdot
std::vector< std::vector< CHybridLSODABalance > > mLocalBalances
const Value & getValue() const
#define PARTITIONING_STEPSIZE
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)=0
bool setValue(const std::string &name, const CType &value)
unsigned C_INT32 * pUINT
std::vector< C_FLOAT64 > temp
CCopasiParameter * getParameter(const std::string &name)
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
CHybridMethodLSODA(const CHybridMethodLSODA &src, const CCopasiContainer *pParent=NULL)
std::set< size_t > * getParticipatesIn(size_t rIndex)
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
std::vector< CHybridLSODAStochFlag > mReactionFlags
size_t removeStochReaction(const size_t index)
const C_FLOAT64 & getTime() const
Definition: CModel.cpp:1190
void removeDeterministicReaction(size_t rIndex)
size_t size() const
Definition: CVector.h:100
long int flag
Definition: f2c.h:52
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
void setTime(const C_FLOAT64 &time)
Definition: CModel.cpp:1187
void refreshConcentration()
Definition: CMetab.cpp:281
std::set< std::string > * getDependsOn(size_t rIndex)
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
void updateNode(const size_t index, const C_FLOAT64 key)
const C_FLOAT64 & getValue() const
virtual void start(const CState *initialState)
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
virtual size_t getIndex(const CCopasiObject *pObject) const
std::ostream & operator<<(std::ostream &os, const CHybridLSODAStochFlag &d)
void initMethod(C_FLOAT64 time)
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
void outputDebug(std::ostream &os, size_t level)
CHybridLSODAStochFlag * mpPrev
size_t getNumModelValues() const
Definition: CModel.cpp:1139
bool addParameter(const CCopasiParameter &parameter)
CDependencyGraph mDG
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
const CCopasiVectorNS< CReaction > * mpReactions
std::set< std::string > * getAffects(size_t rIndex)
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
Definition: CModel.h:50
std::vector< metabStatus > mMetabFlags
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
virtual size_t numCols() const
Definition: CMatrix.h:144
std::string suitableForStochasticSimulation() const
Definition: CModel.cpp:3901
CModel * getModel() const
void insertDeterministicReaction(size_t rIndex)
std::vector< std::vector< CHybridLSODABalance > > mLocalSubstrates
std::vector< C_FLOAT64 > mAmuOld
void calculateDerivative(std::vector< C_FLOAT64 > &deriv)
CHybridLSODAStochFlag * mFirstReactionFlag
void initializeIndexPointer(const size_t numberOfReactions)
#define LOWER_STOCH_LIMIT
Definition: CHybridMethod.h:60
void fireReaction(size_t rIndex)
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
std::set< size_t > mUpdateSet
void calculateAmu(size_t rIndex)
#define max(a, b)
Definition: f2c.h:176