COPASI API  4.16.103
CTrajectoryMethodDsaLsodar.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2014 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 /**
7  * CTrajectoryMethodDsaLsodar
8  *
9  * This class implements an hybrid algorithm for the simulation of a
10  * biochemical system over time.
11  */
12 
13 /* DEFINE ********************************************************************/
14 
15 #ifdef WIN32
16 #if _MSC_VER < 1600
17 #define min _cpp_min
18 #define max _cpp_max
19 #endif // _MSC_VER
20 #endif // WIN32
21 
22 #include <limits.h>
23 
24 #include "copasi.h"
25 
27 #include "CTrajectoryProblem.h"
28 #include "model/CModel.h"
29 #include "model/CMetab.h"
30 #include "model/CReaction.h"
31 #include "model/CState.h"
32 #include "model/CChemEq.h"
33 #include "model/CChemEqElement.h"
34 #include "model/CCompartment.h"
36 #include "utilities/CMatrix.h"
41 
43  mSpeciesMultiplier(0),
44  mMethodSpecies(0),
45  mModelSpecies(0),
46  mCalculations(),
47  mDependentReactions(0),
48  mSubstrateMultiplier(0),
49  mMethodSubstrates(0),
50  mModelSubstrates(0),
51  mpParticleFlux(NULL),
52  mSpeciesIndex(0)
53 {}
54 
56  mSpeciesMultiplier(src.mSpeciesMultiplier),
57  mMethodSpecies(src.mMethodSpecies),
58  mModelSpecies(src.mModelSpecies),
59  mCalculations(src.mCalculations),
60  mDependentReactions(src.mDependentReactions),
61  mSubstrateMultiplier(src.mSubstrateMultiplier),
62  mMethodSubstrates(src.mMethodSubstrates),
63  mModelSubstrates(src.mModelSubstrates),
64  mpParticleFlux(src.mpParticleFlux),
65  mSpeciesIndex(src.mSpeciesIndex)
66 {}
67 
69 {}
70 
72 {
73  mSpeciesMultiplier = rhs.mSpeciesMultiplier;
74  mMethodSpecies = rhs.mMethodSpecies;
75  mModelSpecies = rhs.mModelSpecies;
76  mCalculations = rhs.mCalculations;
77  mDependentReactions = rhs.mDependentReactions;
78  mSubstrateMultiplier = rhs.mSubstrateMultiplier;
79  mMethodSubstrates = rhs.mMethodSubstrates;
80  mModelSubstrates = rhs.mModelSubstrates;
81  mpParticleFlux = rhs.mpParticleFlux;
82  mSpeciesIndex = rhs.mSpeciesIndex;
83 
84  return * this;
85 }
86 
88  mSpeciesToReactions(),
89  mLowerThreshold(),
90  mUpperThreshold(),
92  mNumReactionSpecies(0),
93  mStochasticReactions(0),
94  mDeterministicReactions(0),
95  mStochasticSpecies(0),
96  mHasStochastic(false),
97  mHasDeterministic(false),
98  mNumLowSpecies(0),
99  mLowSpecies(0)
100 {}
101 
103  mSpeciesToReactions(src.mSpeciesToReactions),
104  mLowerThreshold(src.mLowerThreshold),
105  mUpperThreshold(src.mUpperThreshold),
107  mNumReactionSpecies(src.mNumReactionSpecies),
108  mStochasticReactions(src.mStochasticReactions),
109  mDeterministicReactions(src.mDeterministicReactions),
110  mStochasticSpecies(src.mStochasticSpecies),
111  mHasStochastic(src.mHasStochastic),
112  mHasDeterministic(src.mHasDeterministic),
113  mNumLowSpecies(src.mNumLowSpecies),
114  mLowSpecies(src.mLowSpecies)
115 {}
116 
118 {}
119 
120 void CTrajectoryMethodDsaLsodar::CPartition::intialize(std::vector< CReactionDependencies > & reactions,
121  const speciesToReactionsMap & speciesToReactions,
122  const C_FLOAT64 & lowerThreshold,
123  const C_FLOAT64 & upperThreshold,
124  const CState & state)
125 {
126  mSpeciesToReactions = speciesToReactions;
127  mLowerThreshold = lowerThreshold;
128  mUpperThreshold = upperThreshold;
129 
130  mFirstReactionSpeciesIndex = mSpeciesToReactions.begin()->first;
131  mNumReactionSpecies = mSpeciesToReactions.rbegin()->first - mFirstReactionSpeciesIndex + 1;
132 
133  mStochasticReactions.resize(reactions.size());
134  mStochasticReactions = NULL;
135 
136  mDeterministicReactions.resize(reactions.size());
137  mDeterministicReactions = NULL;
138 
139  mNumLowSpecies.resize(reactions.size());
140  mNumLowSpecies = 0;
141 
142  mStochasticSpecies.resize(mNumReactionSpecies);
143  mStochasticSpecies = false;
144 
145  mLowSpecies.resize(mNumReactionSpecies);
146  mLowSpecies = false;
147 
148  mHasStochastic = false;
149  mHasDeterministic = false;
150 
151  // Create the initial partition by first counting species with low particle number per reaction
152  const C_FLOAT64 * pValue = state.beginIndependent() - 1 + mFirstReactionSpeciesIndex;
153  const C_FLOAT64 * pValueEnd = pValue + mNumReactionSpecies;
154  bool * pLowSpecies = mLowSpecies.array();
155 
156  C_FLOAT64 CriticalValue = (mLowerThreshold + mUpperThreshold) * 0.5;
157 
158  size_t Index = mFirstReactionSpeciesIndex;
159 
160  for (; pValue != pValueEnd; ++pValue, ++Index, ++pLowSpecies)
161  {
162  if (*pValue < CriticalValue)
163  {
164  *pLowSpecies = true;
165  std::pair< speciesToReactionsMap::const_iterator, speciesToReactionsMap::const_iterator > Range
166  = mSpeciesToReactions.equal_range(Index);
167 
168  for (; Range.first != Range.second; ++Range.first)
169  {
170  mNumLowSpecies[Range.first->second]++;
171  }
172  }
173  }
174 
175  // Reactions for which the count of species with low numbers is non zero
176  // are treated stochastically.
177  const size_t * pLow = mNumLowSpecies.array();
178  const size_t * pLowEnd = pLow + mNumLowSpecies.size();
179  std::vector< CReactionDependencies >::iterator itReaction = reactions.begin();
180  CReactionDependencies ** ppStochastic = mStochasticReactions.array();
181  CReactionDependencies ** ppDeterministic = mDeterministicReactions.array();
182 
183  for (; pLow != pLowEnd; ++pLow, ++itReaction, ++ppStochastic, ++ppDeterministic)
184  {
185  if (*pLow != 0)
186  {
187  *ppStochastic = &(*itReaction);
188  mHasStochastic = true;
189  }
190  else
191  {
192  *ppDeterministic = &(*itReaction);
193  mHasDeterministic = true;
194  }
195  }
196 
197  determineStochasticSpecies();
198 }
199 
201 {
202  // Modify the partition if species values move beyond the threshold
203  const C_FLOAT64 * pValue = state.beginIndependent() - 1 + mFirstReactionSpeciesIndex;
204  const C_FLOAT64 * pValueEnd = pValue + mNumReactionSpecies;
205  bool * pLowSpecies = mLowSpecies.array();
206 
207  bool PartitionChanged = false;
208 
209  size_t Index = 0;
210 
211  for (; pValue != pValueEnd; ++pValue, ++Index, ++pLowSpecies)
212  {
213  if (!*pLowSpecies && *pValue < mLowerThreshold)
214  {
215  *pLowSpecies = true;
216  PartitionChanged = true;
217 
218  mStochasticSpecies[Index] = true;
219 
220  std::pair< speciesToReactionsMap::const_iterator, speciesToReactionsMap::const_iterator > Range
221  = mSpeciesToReactions.equal_range(Index + mFirstReactionSpeciesIndex);
222 
223  for (; Range.first != Range.second; ++Range.first)
224  {
225  mNumLowSpecies[Range.first->second]++;
226  }
227  }
228  else if (*pLowSpecies && *pValue > mUpperThreshold)
229  {
230  *pLowSpecies = false;
231  PartitionChanged = true;
232 
233  mStochasticSpecies[Index] = true;
234 
235  std::pair< speciesToReactionsMap::const_iterator, speciesToReactionsMap::const_iterator > Range
236  = mSpeciesToReactions.equal_range(Index);
237 
238  for (; Range.first != Range.second; ++Range.first)
239  {
240  mNumLowSpecies[Range.first->second]--;
241  }
242  }
243  }
244 
245  if (!PartitionChanged)
246  {
247  return false;
248  }
249 
250  PartitionChanged = false;
251 
252  // Reactions for which the count of species with low numbers is non zero
253  // are treated stochastically.
254  const size_t * pLow = mNumLowSpecies.array();
255  const size_t * pLowEnd = pLow + mNumLowSpecies.size();
256  CReactionDependencies ** ppStochastic = mStochasticReactions.array();
257  CReactionDependencies ** ppDeterministic = mDeterministicReactions.array();
258  mHasStochastic = false;
259  mHasDeterministic = false;
260 
261  for (; pLow != pLowEnd; ++pLow, ++ppStochastic, ++ppDeterministic)
262  {
263  if (*pLow != 0)
264  {
265  mHasStochastic = true;
266 
267  if (*ppDeterministic != NULL)
268  {
269  PartitionChanged = true;
270  *ppStochastic = *ppDeterministic;
271  *ppDeterministic = NULL;
272  }
273  }
274  else
275  {
276  mHasDeterministic = true;
277 
278  if (*ppStochastic != NULL)
279  {
280  PartitionChanged = true;
281  *ppDeterministic = *ppStochastic;
282  *ppStochastic = NULL;
283  }
284  }
285  }
286 
287  if (PartitionChanged)
288  {
289  // determineStochasticSpecies();
290  }
291 
292  return PartitionChanged;
293 }
294 
296 {
297  // All species which participate in stochastic reactions are treated stochastically
298  CReactionDependencies **ppReaction = mStochasticReactions.array();
299  CReactionDependencies **ppReactionEnd = ppReaction + mStochasticReactions.size();
300 
301  mStochasticSpecies = false;
302 
303  for (; ppReaction != ppReactionEnd; ++ppReaction)
304  {
305  if (*ppReaction != NULL)
306  {
307  size_t * pSpeciesIndex = (*ppReaction)->mSpeciesIndex.array();
308  size_t * pSpeciesIndexEnd = pSpeciesIndex + (*ppReaction)->mSpeciesIndex.size();
309 
310  for (; pSpeciesIndex != pSpeciesIndexEnd; ++pSpeciesIndex)
311  {
312  mStochasticSpecies[*pSpeciesIndex - mFirstReactionSpeciesIndex] = true;
313  }
314  }
315  }
316 }
317 
318 /**
319  * Default constructor.
320  */
322  const CCopasiContainer * pParent):
323  CLsodaMethod(subType, pParent)
324 {
327 }
328 
330  const CCopasiContainer * pParent):
331  CLsodaMethod(src, pParent)
332 {
335 }
336 
337 /**
338  * Destructor.
339  */
341 {
342  cleanup();
344 }
345 
347 {
348  CCopasiParameter *pParm;
349 
350  mpMaxSteps =
351  assertParameter("Max Internal Steps", CCopasiParameter::UINT, (C_INT32) 1000000)->getValue().pUINT;
352  mpLowerLimit =
354  mpUpperLimit =
357  assertParameter("Partitioning Interval", CCopasiParameter::UINT, (unsigned C_INT32) 1)->getValue().pUINT;
359  assertParameter("Partitioning Stepsize", CCopasiParameter::UDOUBLE, (C_FLOAT64) 0.001)->getValue().pUDOUBLE;
360 
361  // Check whether we have a method with the old parameter names
362  if ((pParm = getParameter("HYBRID.MaxSteps")) != NULL)
363  {
364  setValue("Max Internal Steps", *pParm->getValue().pUINT);
365  removeParameter("HYBRID.MaxSteps");
366 
367  if ((pParm = getParameter("HYBRID.LowerStochLimit")) != NULL)
368  {
369  setValue("Lower Limit", *pParm->getValue().pDOUBLE);
370  removeParameter("HYBRID.LowerStochLimit");
371  }
372 
373  if ((pParm = getParameter("HYBRID.UpperStochLimit")) != NULL)
374  {
375  setValue("Upper Limit", *pParm->getValue().pDOUBLE);
376  removeParameter("HYBRID.UpperStochLimit");
377  }
378 
379  if ((pParm = getParameter("HYBRID.PartitioningInterval")) != NULL)
380  {
381  setValue("Partitioning Interval", *pParm->getValue().pUINT);
382  removeParameter("HYBRID.PartitioningInterval");
383  }
384 
385  if ((pParm = getParameter("UseRandomSeed")) != NULL)
386  {
387  removeParameter("UseRandomSeed");
388  }
389 
390  if ((pParm = getParameter("")) != NULL)
391  {
392  removeParameter("");
393  }
394  }
395 }
396 
397 // virtual
399 {
400  bool success = CLsodaMethod::elevateChildren();
401 
403  return success;
404 }
405 
406 // virtual
408 {
411 }
412 
413 // virtual
415 {
416  // do several steps:
417  C_FLOAT64 Time = mpCurrentState->getTime();
418  C_FLOAT64 EndTime = Time + deltaT;
419 
420  C_FLOAT64 Tolerance = 100.0 * (fabs(EndTime) * std::numeric_limits< C_FLOAT64 >::epsilon() + std::numeric_limits< C_FLOAT64 >::min());
421 
422  size_t Steps = 0;
423 
424  while (fabs(Time - EndTime) > Tolerance)
425  {
426  Time += doSingleStep(Time, EndTime);
427 
428  if (++Steps > *mpMaxSteps)
429  {
431  }
432  }
433 
434  mMethodState.setTime(EndTime); // make sure the end time is exactly the time reach to next interval
436 
437  return NORMAL;
438 }
439 
440 // virtual
442 {
443  C_FLOAT64 DeltaT = 0.0;
444  bool FireReaction = false;
445 
446  // if there are stochastic reactions
447  if (mPartition.mHasStochastic) // there is at least one stochastic reaction
448  {
450  {
451  if (mA0 != 0)
452  {
453  mNextReactionTime = curTime - log(mpRandomGenerator->getRandomOO()) / mA0;
454 
455  // We are sure that we have at least 1 reaction
456  mNextReactionIndex = 0;
457 
458  C_FLOAT64 sum = 0.0;
460 
461  C_FLOAT64 * pAmu = mAmu.array();
462  C_FLOAT64 * endAmu = pAmu + mNumReactions;
463  CReactionDependencies **ppStochastic = mPartition.mStochasticReactions.array();
464 
465  // Only consider stochastic reactions
466  for (; (sum <= rand) && (pAmu != endAmu); ++pAmu, ++mNextReactionIndex, ++ppStochastic)
467  {
468  if (*ppStochastic != NULL)
469  {
470  sum += *pAmu;
471  }
472  }
473 
474  assert(mNextReactionIndex > 0);
475 
476  mNextReactionIndex--;
477  }
478  else
479  {
480  mNextReactionTime = std::numeric_limits< C_FLOAT64 >::infinity();
482  }
483  }
484 
485  if (mNextReactionTime <= endTime)
486  {
487  FireReaction = true;
488  DeltaT = mNextReactionTime - curTime;
489  }
490  else
491  {
492  DeltaT = std::min(*mpPartitioningSteps, endTime - curTime);
493  }
494  }
495  else
496  {
497  DeltaT = std::min(*mpPartitioningSteps, endTime - curTime);
498  }
499 
500  // if there are deterministic reactions
501  if (mPartition.mHasDeterministic) // there is at least one deterministic reaction
502  {
505 
506  // Now the model state and mMethodState are identical and the propensities are
507  // calculated accordingly.
508  }
509 
510  if (FireReaction)
511  {
514 
516 
517  mNextReactionTime = std::numeric_limits< C_FLOAT64 >::infinity();
519 
520  // Now the model state and mMethodState are identical and the propensities are
521  // calculated accordingly.
522  }
523 
525  {
527  {
529  }
530 
532  }
533 
534  return DeltaT;
535 }
536 
537 // virtual
538 void CTrajectoryMethodDsaLsodar::start(const CState * initialState)
539 {
540  CLsodaMethod::start(initialState);
541 
543  mDoCorrection = true;
544  else
545  mDoCorrection = false;
546 
547  //initialize the vector of ints that contains the particle numbers
548  //for the discrete simulation. This also floors all particle numbers in the model.
549 
551 
554  mAmu = 0.0;
555 
556  const CStateTemplate & StateTemplate = mpModel->getStateTemplate();
557 
558  CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
559  CModelEntity *const* endEntity = StateTemplate.endFixed();
561 
563  size_t IndexSpecies = 1;
564 
565  for (; ppEntity != endEntity; ++ppEntity, ++pValue, ++IndexSpecies)
566  {
567  if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
568  {
569  *pValue = floor(*pValue + 0.5);
570 
571  if (mFirstReactionSpeciesIndex == 0 &&
572  (*ppEntity)->getStatus() == CModelEntity::REACTIONS)
573  {
574  mFirstReactionSpeciesIndex = IndexSpecies;
575  }
576  }
577  }
578 
579  // Update the model state so that the species are all represented by integers.
581  mpModel->updateSimulatedValues(false); //for assignments
582 
583  // TODO CRITICAL Handle species of type ASSIGNMENT.
584  // These need to be checked whether they are sufficiently close to an integer
585 
586  // Build the reaction dependencies
587  C_FLOAT64 * pMethodStateValue = mMethodState.beginIndependent() - 1;
588  size_t IndexReactions = 0;
589  std::multimap< size_t, size_t> SpeciesToReactions;
590 
593  std::vector< CReactionDependencies >::iterator itDependencies = mReactionDependencies.begin();
594 
595  for (; it != end; ++it)
596  {
597  const CCopasiVector<CChemEqElement> & Balances = (*it)->getChemEq().getBalances();
598  const CCopasiVector<CChemEqElement> & Substrates = (*it)->getChemEq().getSubstrates();
599 
600  // This reactions does not change anything we ignore it
601  if (Balances.size() == 0 && Substrates.size() == 0)
602  {
603  continue;
604  }
605 
606  itDependencies->mpParticleFlux = (C_FLOAT64 *)(*it)->getParticleFluxReference()->getValuePointer();
607 
608  itDependencies->mSpeciesMultiplier.resize(Balances.size());
609  itDependencies->mMethodSpecies.resize(Balances.size());
610  itDependencies->mModelSpecies.resize(Balances.size());
611 
612  std::set< size_t > SpeciesIndexSet;
613 
614  std::set< const CCopasiObject * > changed;
615 
616  // The time is always updated
617  changed.insert(mpModel->getValueReference());
618 
621 
622  size_t Index = 0;
623 
624  for (; itBalance != endBalance; ++itBalance)
625  {
626  const CMetab * pMetab = (*itBalance)->getMetabolite();
627 
628  if (pMetab->getStatus() == CModelEntity::REACTIONS)
629  {
630  size_t SpeciesIndex = StateTemplate.getIndex(pMetab);
631 
632  itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
633  itDependencies->mMethodSpecies[Index] = pMethodStateValue + SpeciesIndex;
634  itDependencies->mModelSpecies[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
635 
636  changed.insert(pMetab->getValueReference());
637 
638  SpeciesToReactions.insert(std::make_pair(SpeciesIndex, IndexReactions));
639  SpeciesIndexSet.insert(SpeciesIndex);
640 
641  Index++;
642  }
643  }
644 
645  // Correct allocation for metabolites which are not determined by reactions
646  itDependencies->mSpeciesMultiplier.resize(Index, true);
647  itDependencies->mMethodSpecies.resize(Index, true);
648  itDependencies->mModelSpecies.resize(Index, true);
649 
650  itDependencies->mSubstrateMultiplier.resize(Substrates.size());
651  itDependencies->mMethodSubstrates.resize(Substrates.size());
652  itDependencies->mModelSubstrates.resize(Substrates.size());
653 
654  CCopasiVector< CChemEqElement >::const_iterator itSubstrate = Substrates.begin();
655  CCopasiVector< CChemEqElement >::const_iterator endSubstrate = Substrates.end();
656 
657  Index = 0;
658 
659  for (; itSubstrate != endSubstrate; ++itSubstrate, ++Index)
660  {
661  const CMetab * pMetab = (*itSubstrate)->getMetabolite();
662 
663  size_t SpeciesIndex = StateTemplate.getIndex(pMetab);
664 
665  itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
666  itDependencies->mMethodSubstrates[Index] = pMethodStateValue + SpeciesIndex;
667  itDependencies->mModelSubstrates[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
668 
669  if (pMetab->getStatus() == CModelEntity::REACTIONS)
670  {
671  SpeciesIndexSet.insert(SpeciesIndex);
672  }
673  }
674 
675  itDependencies->mSpeciesIndex.resize(SpeciesIndexSet.size());
676  size_t * pSpeciesIndex = itDependencies->mSpeciesIndex.array();
677  std::set< size_t >::const_iterator itSet = SpeciesIndexSet.begin();
678  std::set< size_t >::const_iterator endSet = SpeciesIndexSet.end();
679 
680  for (; itSet != endSet; ++itSet, ++pSpeciesIndex)
681  {
682  *pSpeciesIndex = *itSet;
683  }
684 
685  std::set< const CCopasiObject * > dependend;
686 
688  itDependencies->mDependentReactions.resize(mNumReactions);
689 
690  Index = 0;
691  size_t Count = 0;
692 
693  for (; itReaction != end; ++itReaction, ++Index)
694  {
695  if ((*itReaction)->getParticleFluxReference()->dependsOn(changed))
696  {
697  dependend.insert((*itReaction)->getParticleFluxReference());
698  itDependencies->mDependentReactions[Count] = Index;
699 
700  Count++;
701  }
702  }
703 
704  itDependencies->mDependentReactions.resize(Count, true);
705 
706  itDependencies->mCalculations = CCopasiObject::buildUpdateSequence(dependend, changed);
707 
708  ++itDependencies;
709  ++IndexReactions;
710  }
711 
712  mNumReactions = IndexReactions;
713 
715  mAmu.resize(mNumReactions, true);
716 
717  mPartition.intialize(mReactionDependencies, SpeciesToReactions,
719 
720  mNextReactionTime = std::numeric_limits< C_FLOAT64 >::infinity();
722 
724 
725  return;
726 }
727 
728 // virtual
729 void CTrajectoryMethodDsaLsodar::evalF(const C_FLOAT64 * t, const C_FLOAT64 * /* y */, C_FLOAT64 * ydot)
730 {
731  // If we have no ODEs add a constant one.
732  if (mNoODE)
733  {
734  *ydot = 1.0;
735  return;
736  }
737 
738  mMethodState.setTime(*t);
739 
742 
743  if (*mpReducedModel)
745  else
747 
748  // Mask derivatives of stochastic species;
749  bool * pStochastic = mPartition.mStochasticSpecies.array();
750  bool * pStochasticEnd = pStochastic + mPartition.mStochasticSpecies.size();
751  C_FLOAT64 * pYdot = ydot + mFirstReactionSpeciesIndex - 1;
752 
753  for (; pStochastic != pStochasticEnd; ++pStochastic, ++pYdot)
754  {
755  if (*pStochastic)
756  {
757  *pYdot = 0;
758  }
759  }
760 
761  return;
762 }
763 
764 // virtual
765 void CTrajectoryMethodDsaLsodar::evalR(const C_FLOAT64 * /* t */, const C_FLOAT64 * /* y */, const C_INT * /* nr */, C_FLOAT64 * /* r */)
766 {
767 }
768 
769 /* PROTECTED METHODS *********************************************************/
770 
771 /**
772  * Cleans up memory, etc.
773  */
775 {
776  delete mpRandomGenerator;
777  mpRandomGenerator = NULL;
778  mpModel = NULL;
779  return;
780 }
781 
782 /* DETERMINISTIC STUFF *******************************************************/
783 
785 {
786  CLsodaMethod::step(deltaT);
787 
790 
791  return;
792 }
793 
794 /* STOCHASTIC STUFF **********************************************************/
795 
796 /**
797  * Executes the specified reaction in the system once.
798  *
799  * @param index A size_t specifying the index of the reaction, which
800  * will be fired.
801  * @param time The current time
802  */
804 {
805  CReactionDependencies & Dependencies = mReactionDependencies[index];
806 
807  // Update the method internal and model species numbers
808  C_FLOAT64 ** ppModelSpecies = Dependencies.mModelSpecies.array();
809  C_FLOAT64 ** ppLocalSpecies = Dependencies.mMethodSpecies.array();
810  C_FLOAT64 * pMultiplier = Dependencies.mSpeciesMultiplier.array();
811  C_FLOAT64 * endMultiplier = pMultiplier + Dependencies.mSpeciesMultiplier.size();
812 
813  for (; pMultiplier != endMultiplier; ++pMultiplier, ++ppLocalSpecies, ++ppModelSpecies)
814  {
815  **ppLocalSpecies += *pMultiplier;
816  **ppModelSpecies = **ppLocalSpecies;
817  }
818 
819  // Calculate all values which depend on the firing reaction
820  std::vector< Refresh * >::const_iterator itCalcualtion = Dependencies.mCalculations.begin();
821  std::vector< Refresh * >::const_iterator endCalcualtion = Dependencies.mCalculations.end();
822 
823  while (itCalcualtion != endCalcualtion)
824  {
825  (**itCalcualtion++)();
826  }
827 
828  // We do not need to update the the method state since the only independent state
829  // values are species of type reaction which are all controlled by the method.
830 
831  // calculate the propensities which depend on the firing reaction
832  size_t * pDependentReaction = Dependencies.mDependentReactions.array();
833  size_t * endDependentReactions = pDependentReaction + Dependencies.mDependentReactions.size();
834 
835  // It suffices to recalculate the propensities for stochastic reactions.
836  for (; pDependentReaction != endDependentReactions; ++pDependentReaction)
837  {
838  if (mPartition.mStochasticReactions[*pDependentReaction] != NULL)
839  {
840  calculateAmu(*pDependentReaction);
841  }
842  }
843 
844  // calculate the total propensity
845  C_FLOAT64 * pAmu = mAmu.array();
846  C_FLOAT64 * endAmu = pAmu + mNumReactions;
847 
848  mA0 = 0.0;
849 
850  // We must only consider the propensities for stochastic reactions
851  CReactionDependencies ** ppStochastic = mPartition.mStochasticReactions.array();
852 
853  for (; pAmu != endAmu; ++pAmu, ++ppStochastic)
854  {
855  if (*ppStochastic != NULL)
856  {
857  mA0 += *pAmu;
858  }
859  }
860 
862  mLsodaStatus = 1;
863 
864  destroyRootMask();
865 
866  return;
867 }
868 
870 {
871  CReactionDependencies & Dependencies = mReactionDependencies[index];
872  C_FLOAT64 & Amu = mAmu[index];
873 
874  Amu = *Dependencies.mpParticleFlux;
875 
876  if (Amu < 0.0)
877  {
878  // TODO CRITICAL Create a warning message
879  Amu = 0.0;
880  }
881 
882  if (!mDoCorrection)
883  {
884  return;
885  }
886 
887  C_FLOAT64 SubstrateMultiplier = 1.0;
888  C_FLOAT64 SubstrateDevisor = 1.0;
889  C_FLOAT64 Multiplicity;
890  C_FLOAT64 LowerBound;
891  C_FLOAT64 Number;
892 
893  bool ApplyCorrection = false;
894 
895  C_FLOAT64 * pMultiplier = Dependencies.mSubstrateMultiplier.array();
896  C_FLOAT64 * endMultiplier = pMultiplier + Dependencies.mSubstrateMultiplier.size();
897  C_FLOAT64 ** ppLocalSubstrate = Dependencies.mMethodSubstrates.array();
898  C_FLOAT64 ** ppModelSubstrate = Dependencies.mModelSubstrates.array();
899 
900  for (; pMultiplier != endMultiplier; ++pMultiplier, ++ppLocalSubstrate, ++ppModelSubstrate)
901  {
902  Multiplicity = *pMultiplier;
903 
904  // TODO We should check the error introduced through rounding.
905  **ppLocalSubstrate = floor(**ppModelSubstrate + 0.5);
906 
907  if (Multiplicity > 1.01)
908  {
909  ApplyCorrection = true;
910 
911  Number = **ppLocalSubstrate;
912 
913  LowerBound = Number - Multiplicity;
914  SubstrateDevisor *= pow(Number, Multiplicity - 1.0); //optimization
915  Number -= 1.0;
916 
917  while (Number > LowerBound)
918  {
919  SubstrateMultiplier *= Number;
920  Number -= 1.0;
921  }
922  }
923  }
924 
925  // at least one substrate particle number is zero
926  if (SubstrateMultiplier < 0.5 || SubstrateDevisor < 0.5)
927  {
928  Amu = 0.0;
929  return;
930  }
931 
932  // rate_factor is the rate function divided by substrate_factor.
933  // It would be more efficient if this was generated directly, since in effect we
934  // are multiplying and then dividing by the same thing (substrate_factor)!
935  //C_FLOAT64 rate_factor = mpModel->getReactions()[index]->calculateParticleFlux();
936  if (ApplyCorrection)
937  {
938  Amu *= SubstrateMultiplier / SubstrateDevisor;
939  }
940 
941  return;
942 }
943 
945 {
946  // It suffices to recalculate the propensities for stochastic reactions.
947  CReactionDependencies **ppStochastic = mPartition.mStochasticReactions.array();
948 
949  for (size_t Index = 0; Index != this->mNumReactions; ++Index, ++ppStochastic)
950  {
951  if (*ppStochastic != NULL)
952  {
953  calculateAmu(Index);
954  }
955  }
956 
958  mpModel->updateSimulatedValues(false); //for assignments
959 
960  // calculate the total propensity
961  C_FLOAT64 * pAmu = mAmu.array();
962  C_FLOAT64 * endAmu = pAmu + mNumReactions;
963 
964  mA0 = 0.0;
965 
966  // We must only consider the propensities for stochastic reactions
967  ppStochastic = mPartition.mStochasticReactions.array();
968 
969  for (; pAmu != endAmu; ++pAmu, ++ppStochastic)
970  {
971  if (*ppStochastic != NULL)
972  {
973  mA0 += *pAmu;
974  assert(mA0 >= 0.0);
975  }
976  }
977 
978  return;
979 }
980 
981 /* TESTING THE MODEL AND SETTING UP THINGS ***********************************/
982 
983 //virtual
985 {
986  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
987 
988  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
989 
990  if (pTP->getDuration() < 0.0)
991  {
992  //back integration not possible
994  return false;
995  }
996 
997  //check for rules
998  size_t i, imax = pTP->getModel()->getNumModelValues();
999 
1000  for (i = 0; i < imax; ++i)
1001  {
1002  if (pTP->getModel()->getModelValues()[i]->getStatus() == CModelEntity::ODE)
1003  {
1004  //ode rule found
1006  return false;
1007  }
1008 
1009  /* if (pTP->getModel()->getModelValues()[i]->getStatus()==CModelEntity::ASSIGNMENT)
1010  if (pTP->getModel()->getModelValues()[i]->isUsed())
1011  {
1012  //used assignment found
1013  CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 19);
1014  return false;
1015  }*/
1016  }
1017 
1018  imax = pTP->getModel()->getNumMetabs();
1019 
1020  for (i = 0; i < imax; ++i)
1021  {
1022  if (pTP->getModel()->getMetabolites()[i]->getStatus() == CModelEntity::ODE)
1023  {
1024  //ode rule found
1026  return false;
1027  }
1028  }
1029 
1030  imax = pTP->getModel()->getCompartments().size();
1031 
1032  for (i = 0; i < imax; ++i)
1033  {
1034  if (pTP->getModel()->getCompartments()[i]->getStatus() == CModelEntity::ODE)
1035  {
1036  //ode rule found
1038  return false;
1039  }
1040  }
1041 
1042  //TODO: rewrite CModel::suitableForStochasticSimulation() to use
1043  // CCopasiMessage
1044  std::string message = pTP->getModel()->suitableForStochasticSimulation();
1045 
1046  if (message != "")
1047  {
1048  //model not suitable, message describes the problem
1049  CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
1050  return false;
1051  }
1052 
1053  /* Max Internal Steps */
1054  if (* getValue("Max Internal Steps").pINT <= 0)
1055  {
1056  //max steps should be at least 1
1058  return false;
1059  }
1060 
1061  /* Lower Limit, Upper Limit */
1062  *mpLowerLimit = * getValue("Lower Limit").pDOUBLE;
1063  *mpUpperLimit = * getValue("Upper Limit").pDOUBLE;
1064 
1065  if (*mpLowerLimit > *mpUpperLimit)
1066  {
1068  return false;
1069  }
1070 
1071  /* Partitioning Interval */
1072  // nothing to be done here so far
1073 
1074  /* Use Random Seed */
1075  // should be checked in the widget later on
1076 
1077  /* Random Seed */
1078  // nothing to be done here
1079 
1080  //events are not supported at the moment
1081  if (pTP->getModel()->getEvents().size() > 0)
1082  {
1084  return false;
1085  }
1086 
1087  return true;
1088 }
CTrajectoryMethodDsaLsodar(const CCopasiMethod::SubType &subType=DsaLsodar, const CCopasiContainer *pParent=NULL)
void calculateAmu(const size_t &index)
#define C_INT
Definition: copasi.h:115
C_INT mLsodaStatus
Definition: CLsodaMethod.h:107
CCopasiVectorN< CEvent > & getEvents()
Definition: CModel.cpp:1110
virtual void evalR(const C_FLOAT64 *t, const C_FLOAT64 *y, const C_INT *nr, C_FLOAT64 *r)
const CCopasiVector< CMetab > & getMetabolites() const
Definition: CModel.cpp:1051
const CCopasiVectorN< CModelValue > & getModelValues() const
Definition: CModel.cpp:1060
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
virtual size_t size() const
size_t getNumMetabs() const
Definition: CModel.cpp:1118
iterator begin()
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
virtual void start(const CState *initialState)
virtual Status step(const double &deltaT)
Definition: CState.h:305
CVector< CReactionDependencies * > mStochasticReactions
virtual void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INVALID_INDEX
Definition: copasi.h:222
size_t getIndex(const CModelEntity *entity) const
Definition: CState.cpp:231
CModel * mpModel
Definition: CLsodaMethod.h:159
#define C_INT32
Definition: copasi.h:90
Definition: CMetab.h:178
virtual bool isValidProblem(const CCopasiProblem *pProblem)
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
virtual Status step(const double &deltaT)
CState mMethodState
Definition: CLsodaMethod.h:74
CReactionDependencies & operator=(const CReactionDependencies &rhs)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
bool removeParameter(const std::string &name)
#define MCTrajectoryMethod
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
iterator end()
std::vector< CType * >::const_iterator const_iterator
Definition: CCopasiVector.h:57
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
std::vector< CReactionDependencies > mReactionDependencies
virtual bool elevateChildren()
CModelEntity ** endFixed()
Definition: CState.cpp:213
void setState(const CState &state)
Definition: CModel.cpp:1785
const C_FLOAT64 & getDuration() const
void integrateDeterministicPart(const C_FLOAT64 &deltaT)
bool * mpReducedModel
Definition: CLsodaMethod.h:52
const Value & getValue() const
bool setValue(const std::string &name, const CType &value)
unsigned C_INT32 * pUINT
void intialize(std::vector< CReactionDependencies > &reactions, const speciesToReactionsMap &speciesToReactions, const C_FLOAT64 &lowerThreshold, const C_FLOAT64 &upperThreshold, const CState &state)
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
CCopasiParameter * getParameter(const std::string &name)
void destroyRootMask()
C_FLOAT64 mTime
Definition: CLsodaMethod.h:102
size_t size() const
Definition: CVector.h:100
void setTime(const C_FLOAT64 &time)
Definition: CModel.cpp:1187
void calculateDerivativesX(C_FLOAT64 *derivativesX)
Definition: CModel.cpp:1941
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
size_t getNumModelValues() const
Definition: CModel.cpp:1139
std::multimap< size_t, size_t > speciesToReactionsMap
const ModelType & getModelType() const
Definition: CModel.cpp:2339
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
static std::vector< Refresh * > buildUpdateSequence(const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
virtual void stateChanged()
void calculateDerivatives(C_FLOAT64 *derivatives)
Definition: CModel.cpp:1903
virtual void * getValuePointer() const
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
const CModelEntity::Status & getStatus() const
void fireReaction(const size_t &index)
std::string suitableForStochasticSimulation() const
Definition: CModel.cpp:3901
CModel * getModel() const
#define min(a, b)
Definition: f2c.h:175
virtual void start(const CState *initialState)
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
CCopasiObject * getValueReference() const