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