COPASI API  4.16.103
CTrajAdaptiveSA.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 #ifdef WIN32
7 # pragma warning (disable: 4786)
8 # pragma warning (disable: 4243)
9 // warning C4355: 'this' : used in base member initializer list
10 # pragma warning (disable: 4355)
11 #endif // WIN32
12 
13 #include <limits>
14 
15 #include <vector>
16 #include <numeric>
17 #include <limits>
18 #include <set>
19 #include <string>
20 #include <cmath>
21 
22 #include "copasi.h"
23 #include "CTrajAdaptiveSA.h"
25 #include "function/CFunction.h"
27 #include "CTrajectoryProblem.h"
28 #include "model/CState.h"
29 #include "model/CCompartment.h"
30 #include "model/CModel.h"
31 
33  mMethodSpeciesIndex(0),
34  mSpeciesMultiplier(0),
35  mMethodSpecies(0),
36  mModelSpecies(0),
37  mCalculations(),
38  mDependentReactions(0),
39  mSubstrateMultiplier(0),
40  mMethodSubstrates(0),
41  mModelSubstrates(0),
42  mpParticleFlux(NULL)
43 {}
44 
46  mMethodSpeciesIndex(src.mMethodSpeciesIndex),
47  mSpeciesMultiplier(src.mSpeciesMultiplier),
48  mMethodSpecies(src.mMethodSpecies),
49  mModelSpecies(src.mModelSpecies),
50  mCalculations(src.mCalculations),
51  mDependentReactions(src.mDependentReactions),
52  mSubstrateMultiplier(src.mSubstrateMultiplier),
53  mMethodSubstrates(src.mMethodSubstrates),
54  mModelSubstrates(src.mModelSubstrates),
55  mpParticleFlux(src.mpParticleFlux)
56 {}
57 
59 {}
60 
62 {
63  mSpeciesMultiplier = rhs.mSpeciesMultiplier;
64  mMethodSpecies = rhs.mMethodSpecies;
65  mModelSpecies = rhs.mModelSpecies;
66  mCalculations = rhs.mCalculations;
67  mDependentReactions = rhs.mDependentReactions;
68  mSubstrateMultiplier = rhs.mSubstrateMultiplier;
69  mMethodSubstrates = rhs.mMethodSubstrates;
70  mModelSubstrates = rhs.mModelSubstrates;
71  mpParticleFlux = rhs.mpParticleFlux;
72 
73  return * this;
74 }
75 
76 // static
78 {
79  CTrajAdaptiveSA * pMethod = NULL;
80 
81  pMethod = new CTrajAdaptiveSA();
82 
83  return pMethod;
84 }
85 
89  mReactionFiring(0),
91  mAvgDX(0),
92  mSigDX(0),
93  mpMethodSpecies(NULL),
95  mpRandomGenerator(CRandom::createGenerator(CRandom::mt19937)),
96  mpModel(NULL),
97  mNumReactions(0),
98  mNumSpecies(0),
99  mMaxSteps(1000000),
100  mNextReactionTime(0.0),
102  mDoCorrection(true),
103  mAmu(0),
104  mPartitionedAmu(0),
105  mMethodState(),
108  mA0(0.0),
109  mMaxStepsReached(false)
110 {
112 }
113 
115  const CCopasiContainer * pParent):
116  CTrajectoryMethod(src, pParent),
117  mMaxReactionFiring(src.mMaxReactionFiring),
118  mReactionFiring(src.mReactionFiring),
119  mPartitionedReactionFiring(src.mPartitionedReactionFiring),
120  mAvgDX(src.mAvgDX),
121  mSigDX(src.mSigDX),
122  mpMethodSpecies(src.mpMethodSpecies),
123  mSpeciesAfterTau(src.mSpeciesAfterTau),
124  mpRandomGenerator(CRandom::createGenerator(CRandom::mt19937)),
125  mpModel(src.mpModel),
126  mNumReactions(src.mNumReactions),
127  mNumSpecies(src.mNumSpecies),
128  mMaxSteps(src.mMaxSteps),
129  mNextReactionTime(src.mNextReactionTime),
130  mNextReactionIndex(src.mNextReactionIndex),
131  mDoCorrection(src.mDoCorrection),
132  mAmu(src.mAmu),
133  mPartitionedAmu(src.mPartitionedAmu),
134  mMethodState(src.mMethodState),
135  mReactionDependencies(src.mReactionDependencies),
136  mPartitionedDependencies(src.mPartitionedDependencies),
137  mA0(src.mA0),
138  mMaxStepsReached(src.mMaxStepsReached)
139 {
141 }
142 
144 {
146 }
147 
149 {
151  assertParameter("Max Internal Steps", CCopasiParameter::INT, (C_INT32) 1000000);
152  assertParameter("Use Random Seed", CCopasiParameter::BOOL, false);
153  assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) 1);
154 }
155 
157 {
159  return true;
160 }
161 
163 {
164  // do several steps
165  C_FLOAT64 Time = mpCurrentState->getTime();
166  C_FLOAT64 EndTime = Time + deltaT;
167 
168  size_t Steps = 0;
169 
170  while (Time < EndTime)
171  {
172  if (mSSAStepCounter > 0)
173  {
174  Time += doSingleSSAStep(Time, EndTime);
175  mSSAStepCounter--;
176  }
177  else
178  {
179  Time += doSingleTauLeapStep(Time, EndTime);
180  }
181 
182  if (++Steps > mMaxSteps)
183  {
185  }
186  }
187 
189  mpCurrentState->setTime(Time);
190 
191  return NORMAL;
192 }
193 
194 void CTrajAdaptiveSA::start(const CState * initialState)
195 {
196  /* get configuration data */
197  mMaxSteps = * getValue("Max Internal Steps").pINT;
198  mEpsilon = * getValue("Epsilon").pDOUBLE;
199 
200  bool useRandomSeed = * getValue("Use Random Seed").pBOOL;
201  unsigned C_INT32 randomSeed = * getValue("Random Seed").pUINT;
202 
203  if (useRandomSeed) mpRandomGenerator->initialize(randomSeed);
204 
205  //mpCurrentState is initialized. This state is not used internally in the
206  //stochastic solver, but it is used for returning the result after each step.
207  *mpCurrentState = *initialState;
208 
210  assert(mpModel);
211 
213  mDoCorrection = true;
214  else
215  mDoCorrection = false;
216 
217  //initialize the vector of ints that contains the particle numbers
218  //for the discrete simulation. This also floors all particle numbers in the model.
219 
222 
228 
231  mAmu = 0.0;
232 
234 
236 
239 
240  // Initialize the
241  // Create a local copy of the state where the particle number are rounded to integers
243 
244  const CStateTemplate & StateTemplate = mpModel->getStateTemplate();
245 
246  CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
247  CModelEntity *const* endEntity = StateTemplate.endFixed();
249 
251  size_t Index = 1;
252 
253  C_FLOAT64 * pSpeciesAfterTau = mSpeciesAfterTau.array();
254 
255  for (; ppEntity != endEntity; ++ppEntity, ++pValue, ++Index)
256  {
257  if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
258  {
259  *pValue = floor(*pValue + 0.5);
260 
261  if ((*ppEntity)->getStatus() == CModelEntity::REACTIONS &&
262  (*ppEntity)->isUsed())
263  {
264  *pSpeciesAfterTau = *pValue;
265  pSpeciesAfterTau++;
266 
268  {
270  mpMethodSpecies = pValue;
271  }
272  }
273  }
274  }
275 
276  // Update the model state so that the species are all represented by integers.
278  mpModel->updateSimulatedValues(false); //for assignments
279 
280  // TODO handle species of type ASSIGNMENT.
281  // These need to be checked whether they are sufficiently close to an integer
282 
283  C_FLOAT64 * pMethodStateValue = mMethodState.beginIndependent() - 1;
284 
285  // Build the reaction dependencies
286  size_t NumReactions = 0;
287 
290  std::vector< CReactionDependencies >::iterator itDependencies = mReactionDependencies.begin();
291 
292  for (; it != end; ++it)
293  {
294  const CCopasiVector<CChemEqElement> & Balances = (*it)->getChemEq().getBalances();
295  const CCopasiVector<CChemEqElement> & Substrates = (*it)->getChemEq().getSubstrates();
296 
297  // This reactions does not change anything we ignore it
298  if (Balances.size() == 0 && Substrates.size() == 0)
299  {
300  continue;
301  }
302 
303  itDependencies->mpParticleFlux = (C_FLOAT64 *)(*it)->getParticleFluxReference()->getValuePointer();
304 
305  itDependencies->mMethodSpeciesIndex.resize(Balances.size());
306  itDependencies->mSpeciesMultiplier.resize(Balances.size());
307  itDependencies->mMethodSpecies.resize(Balances.size());
308  itDependencies->mModelSpecies.resize(Balances.size());
309 
310  std::set< const CCopasiObject * > changed;
311 
312  // The time is always updated
313  changed.insert(mpModel->getValueReference());
314 
317 
318  size_t Index = 0;
319 
320  for (; itBalance != endBalance; ++itBalance)
321  {
322  const CMetab * pMetab = (*itBalance)->getMetabolite();
323 
324  if (pMetab->getStatus() == CModelEntity::REACTIONS)
325  {
326  itDependencies->mMethodSpeciesIndex[Index] = StateTemplate.getIndex(pMetab) - mFirstReactionSpeciesIndex;
327  itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
328  itDependencies->mMethodSpecies[Index] = pMethodStateValue + StateTemplate.getIndex(pMetab);
329  itDependencies->mModelSpecies[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
330 
331  changed.insert(pMetab->getValueReference());
332 
333  Index++;
334  }
335  }
336 
337  // Correct allocation for metabolites which are not determined by reactions
338  itDependencies->mMethodSpeciesIndex.resize(Index, true);
339  itDependencies->mSpeciesMultiplier.resize(Index, true);
340  itDependencies->mMethodSpecies.resize(Index, true);
341  itDependencies->mModelSpecies.resize(Index, true);
342 
343  itDependencies->mSubstrateMultiplier.resize(Substrates.size());
344  itDependencies->mMethodSubstrates.resize(Substrates.size());
345  itDependencies->mModelSubstrates.resize(Substrates.size());
346 
347  CCopasiVector< CChemEqElement >::const_iterator itSubstrate = Substrates.begin();
348  CCopasiVector< CChemEqElement >::const_iterator endSubstrate = Substrates.end();
349 
350  Index = 0;
351 
352  for (; itSubstrate != endSubstrate; ++itSubstrate, ++Index)
353  {
354  const CMetab * pMetab = (*itSubstrate)->getMetabolite();
355 
356  itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
357  itDependencies->mMethodSubstrates[Index] = pMethodStateValue + StateTemplate.getIndex(pMetab);
358  itDependencies->mModelSubstrates[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
359  }
360 
361  std::set< const CCopasiObject * > dependend;
362 
364  itDependencies->mDependentReactions.resize(mNumReactions);
365 
366  Index = 0;
367  size_t Count = 0;
368 
369  for (; itReaction != end; ++itReaction, ++Index)
370  {
371  if ((*itReaction)->getParticleFluxReference()->dependsOn(changed))
372  {
373  dependend.insert((*itReaction)->getParticleFluxReference());
374  itDependencies->mDependentReactions[Count] = Index;
375 
376  Count++;
377  }
378  }
379 
380  itDependencies->mDependentReactions.resize(Count, true);
381 
382  itDependencies->mCalculations = CCopasiObject::buildUpdateSequence(dependend, changed);
383 
384  ++itDependencies;
385  ++NumReactions;
386  }
387 
388  mNumReactions = NumReactions;
389 
391  mAmu.resize(mNumReactions, true);
392 
393  mA0 = 0;
394  size_t i;
395 
396  for (i = 0; i < mNumReactions; i++)
397  {
398  mA0 += calculateAmu(i);
399  }
400 
401  mMaxStepsReached = false;
404 
405  mSSAStepCounter = 0;
406 
407  return;
408 }
409 
410 // virtual
412 {
413  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
414 
415  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
416 
417  if (pTP->getDuration() < 0.0)
418  {
419  //back integration not possible
421  return false;
422  }
423 
424  if (pTP->getModel()->getTotSteps() < 1)
425  {
426  //at least one reaction necessary
428  return false;
429  }
430 
431  // check for ODEs
432  const CStateTemplate & StateTemplate = pTP->getModel()->getStateTemplate();
433  CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
434  CModelEntity *const* ppEntityEnd = StateTemplate.endIndependent();
435 
436  for (; ppEntity != ppEntityEnd; ++ppEntity)
437  {
438  if ((*ppEntity)->getStatus() == CModelEntity::ODE)
439  {
440  if (dynamic_cast<const CModelValue *>(*ppEntity) != NULL)
441  {
442  // global quantity ode rule found
444  return false;
445  }
446  else if (dynamic_cast<const CCompartment *>(*ppEntity) != NULL)
447  {
448  // compartment ode rule found
450  return false;
451  }
452  else
453  {
454  // species ode rule found
456  return false;
457  }
458  }
459  }
460 
461  //TODO: rewrite CModel::suitableForStochasticSimulation() to use
462  // CCopasiMessage
463  std::string message = pTP->getModel()->suitableForStochasticSimulation();
464 
465  if (message != "")
466  {
467  //model not suitable, message describes the problem
468  CCopasiMessage(CCopasiMessage::EXCEPTION, message.c_str());
469  return false;
470  }
471 
472  if (* getValue("Max Internal Steps").pINT <= 0)
473  {
474  //max steps should be at least 1
476  return false;
477  }
478 
479  return true;
480 }
481 
483 {
484  std::vector< CReactionDependencies >::const_iterator itDependencies = mReactionDependencies.begin();
485  std::vector< CReactionDependencies >::const_iterator endDependencies = mReactionDependencies.end();
486  size_t * pMaxReactionFiring = mMaxReactionFiring.array();
487 
488  for (; itDependencies != endDependencies; ++itDependencies, ++pMaxReactionFiring)
489  {
490  C_FLOAT64 *const* ppModelSpecies = itDependencies->mModelSpecies.array();
491  const C_FLOAT64 * pMultiplicity = itDependencies->mSpeciesMultiplier.array();
492  const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + itDependencies->mSpeciesMultiplier.size();
493 
494  *pMaxReactionFiring = std::numeric_limits< size_t >::max(); // Assigned a maximum value
495 
496  for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, ppModelSpecies++)
497  {
498  if (*pMultiplicity < 0)
499  {
500  size_t TmpMax;
501 
502  if ((TmpMax = (size_t) fabs(**ppModelSpecies / *pMultiplicity)) < *pMaxReactionFiring)
503  {
504  *pMaxReactionFiring = TmpMax;
505  }
506  }
507  }
508  }
509 
510  size_t NonCriticalReactions = 0;
511  size_t InsertCritical = mNumReactions - 1;
512 
513  pMaxReactionFiring = mMaxReactionFiring.array();
514  const C_FLOAT64 * pAmu = mAmu.array();
515 
516  itDependencies = mReactionDependencies.begin();
517  C_FLOAT64 * pReactionFiring = mReactionFiring.array();
518 
519  for (; itDependencies != endDependencies; ++itDependencies, ++pAmu, ++pMaxReactionFiring, ++pReactionFiring)
520  {
521  *pReactionFiring = 0;
522 
523  if (*pAmu == 0 ||
524  *pMaxReactionFiring > UPPER_LIMIT)
525  {
526  mPartitionedDependencies[NonCriticalReactions] = &(*itDependencies);
527  mPartitionedAmu[NonCriticalReactions] = pAmu;
528  mPartitionedReactionFiring[NonCriticalReactions] = pReactionFiring;
529  NonCriticalReactions++;
530  }
531  else
532  {
533  mPartitionedDependencies[InsertCritical] = &(*itDependencies);
534  mPartitionedAmu[InsertCritical] = pAmu;
535  mPartitionedReactionFiring[InsertCritical] = pReactionFiring;
536  InsertCritical--;
537  }
538  }
539 
540  mAvgDX = 0.0;
541  mSigDX = 0.0;
542 
543  const CReactionDependencies **ppOrderedReactions = mPartitionedDependencies.array();
544  const C_FLOAT64 ** ppOrderedAmu = mPartitionedAmu.array();
545  const C_FLOAT64 ** ppOrderedAmuEnd = ppOrderedAmu + NonCriticalReactions;
546 
547  for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedReactions, ++ppOrderedAmu)
548  {
549  const C_FLOAT64 * pMultiplicity = (*ppOrderedReactions)->mSpeciesMultiplier.array();
550  const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + (*ppOrderedReactions)->mSpeciesMultiplier.size();
551  const size_t * pIndex = (*ppOrderedReactions)->mMethodSpeciesIndex.array();
552 
553  for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
554  {
555  mAvgDX[*pIndex] += **ppOrderedAmu * *pMultiplicity;
556  mSigDX[*pIndex] += **ppOrderedAmu * *pMultiplicity * *pMultiplicity;
557  }
558  }
559 
560  C_FLOAT64 MaxTau;
561 
562  if (NonCriticalReactions == 0)
563  {
564  MaxTau = 0;
565  }
566  else
567  {
568  C_FLOAT64 ex, t1, t2;
569  t1 = t2 = std::numeric_limits< C_FLOAT64 >::infinity();
570 
571  for (size_t i = 0; i < mNumReactionSpecies; i++)
572  {
573  C_FLOAT64 t3, t4, t5, t6;
574 
575  ex = (*(mpMethodSpecies + i) * mEpsilon) + 1.0;
576  t3 = fabs(mAvgDX[i]);
577  t4 = mSigDX[i];
578  t5 = ex / t3;
579  t6 = ex * ex / t4;
580 
581  if (t3 != 0 && t5 < t1) t1 = t5;
582 
583  if (t4 != 0 && t6 < t2) t2 = t6;
584  }
585 
586  if (t1 < t2)
587  MaxTau = t1;
588  else
589  MaxTau = t2;
590  }
591 
592  if (MaxTau < (SSA_MULTIPLE / mA0) ||
593  MaxTau == std::numeric_limits< C_FLOAT64 >::infinity())
594  {
598 
599  return 0;
600  }
601 
602  C_FLOAT64 AmuCritical = 0;
603  ppOrderedAmu = mPartitionedAmu.array() + NonCriticalReactions;
604  ppOrderedAmuEnd = mPartitionedAmu.array() + mNumReactions;
605 
606  for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedAmu)
607  {
608  AmuCritical += **ppOrderedAmu;
609  }
610 
611  C_FLOAT64 CriticalReactionTau;
612 
613  if (NonCriticalReactions == mNumReactions || AmuCritical == 0)
614  CriticalReactionTau = std::numeric_limits< C_FLOAT64 >::infinity();
615  else
616  CriticalReactionTau = -log(mpRandomGenerator->getRandomOO()) / AmuCritical;
617 
618  bool isUpdated = false;
619 
620  C_FLOAT64 Tau;
621 
622  while (!isUpdated)
623  {
624  Tau = MaxTau;
625 
626  if (CriticalReactionTau < Tau || Tau == 0) Tau = CriticalReactionTau;
627 
628  if ((endTime - curTime) < Tau) Tau = (endTime - curTime);
629 
630  // Calculate the firing of non critical reactions
631  ppOrderedAmu = mPartitionedAmu.array();
632  ppOrderedAmuEnd = mPartitionedAmu.array() + NonCriticalReactions;
633 
634  C_FLOAT64 **ppOrderedReactionFiring = mPartitionedReactionFiring.array();
635 
636  for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedAmu, ++ppOrderedReactionFiring)
637  {
638  C_FLOAT64 Lambda = **ppOrderedAmu * Tau;
639 
640  if (Lambda < 0.0)
642  else if (Lambda > 2.0e9)
644 
645  **ppOrderedReactionFiring = mpRandomGenerator->getRandomPoisson(Lambda);
646  }
647 
648  size_t CriticalReactionIndex = C_INVALID_INDEX;
649 
650  if (CriticalReactionTau <= Tau)
651  {
652  // Determine the critical reaction which fires.
653  C_FLOAT64 sum = 0;
654  C_FLOAT64 rand = mpRandomGenerator->getRandomOO() * AmuCritical;
655 
656  ppOrderedAmu = mPartitionedAmu.array() + NonCriticalReactions;
657  ppOrderedAmuEnd = mPartitionedAmu.array() + mNumReactions;
658 
659  CriticalReactionIndex = NonCriticalReactions;
660 
661  for (; (sum < rand) && (ppOrderedAmu != ppOrderedAmuEnd); ++ppOrderedAmu, ++CriticalReactionIndex)
662  {
663  sum += **ppOrderedAmu;
664  }
665 
666  CriticalReactionIndex--;
667  }
668 
669  C_FLOAT64 * pSpeciesAfterTau = mSpeciesAfterTau.array();
670  C_FLOAT64 * pSpeciesAfterTauEnd = pSpeciesAfterTau + mNumReactionSpecies;
671  C_FLOAT64 * pMethodSpecies = mpMethodSpecies;
672 
673  for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau, ++pMethodSpecies)
674  {
675  *pSpeciesAfterTau = *pMethodSpecies;
676  }
677 
678  ppOrderedReactions = mPartitionedDependencies.array();
679  const CReactionDependencies **ppOrderedReactionsEnd = ppOrderedReactions + NonCriticalReactions;
680  ppOrderedReactionFiring = mPartitionedReactionFiring.array();
681 
682  for (; ppOrderedReactions != ppOrderedReactionsEnd; ++ppOrderedReactions, ++ppOrderedReactionFiring)
683  {
684  const C_FLOAT64 * pMultiplicity = (*ppOrderedReactions)->mSpeciesMultiplier.array();
685  const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + (*ppOrderedReactions)->mSpeciesMultiplier.size();
686  const size_t *pIndex = (*ppOrderedReactions)->mMethodSpeciesIndex.array();
687 
688  for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
689  {
690  mSpeciesAfterTau[*pIndex] += (*pMultiplicity * **ppOrderedReactionFiring);
691  }
692 
693  **ppOrderedReactionFiring = 0.0;
694  }
695 
696  if (CriticalReactionIndex != C_INVALID_INDEX)
697  {
698  const CReactionDependencies * pCriticalReaction = mPartitionedDependencies[CriticalReactionIndex];
699 
700  const C_FLOAT64 * pMultiplicity = pCriticalReaction->mSpeciesMultiplier.array();
701  const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + pCriticalReaction->mSpeciesMultiplier.size();
702  const size_t *pIndex = pCriticalReaction->mMethodSpeciesIndex.array();
703 
704  for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
705  {
706  mSpeciesAfterTau[*pIndex] += *pMultiplicity;
707  }
708  }
709 
710  pSpeciesAfterTau = mSpeciesAfterTau.array();
711 
712  for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau)
713  {
714  if (*pSpeciesAfterTau < 0)
715  break;
716  }
717 
718  if (pSpeciesAfterTau != pSpeciesAfterTauEnd)
719  {
720  isUpdated = false;
721  MaxTau = MaxTau / 2.0;
722  }
723  else
724  {
725  isUpdated = true;
726 
727  pSpeciesAfterTau = mSpeciesAfterTau.array();
728  pMethodSpecies = mpMethodSpecies;
729 
730  for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau, ++pMethodSpecies)
731  {
732  *pMethodSpecies = *pSpeciesAfterTau;
733  }
734  }
735  }
736 
737  // Update the model time (for explicitly time dependent models)
738  mMethodState.setTime(curTime + Tau);
741 
742  mA0 = 0;
743  size_t i;
744 
745  for (i = 0; i < mNumReactions; i++)
746  {
747  mA0 += calculateAmu(i);
748  }
749 
750  return Tau;
751 }
752 
754 {
756  {
757  if (mA0 == 0)
758  {
759  return endTime - curTime;
760  }
761 
762  // We need to throw an exception if mA0 is NaN
763  if (isnan(mA0))
764  {
766  }
767 
768  mNextReactionTime = curTime - log(mpRandomGenerator->getRandomOO()) / mA0;
769  }
770 
771  if (mNextReactionTime >= endTime)
772  {
773  return endTime - curTime;
774  }
775 
776  // Update the model time (for explicitly time dependent models)
778 
779  // We are sure that we have at least 1 reaction
780  mNextReactionIndex = 0;
781 
782  C_FLOAT64 sum = 0.0;
784 
785  C_FLOAT64 * pAmu = mAmu.array();
786  C_FLOAT64 * endAmu = pAmu + mNumReactions;
787 
788  for (; (sum < rand) && (pAmu != endAmu); ++pAmu, ++mNextReactionIndex)
789  {
790  sum += *pAmu;
791  }
792 
794 
796 
797  // Update the method internal and model species numbers
798  C_FLOAT64 ** ppModelSpecies = Dependencies.mModelSpecies.array();
799  C_FLOAT64 ** ppLocalSpecies = Dependencies.mMethodSpecies.array();
800  C_FLOAT64 * pMultiplicity = Dependencies.mSpeciesMultiplier.array();
801  C_FLOAT64 * pMultiplicityEnd = pMultiplicity + Dependencies.mSpeciesMultiplier.size();
802  //if(mAmu.array()[mNextReactionIndex]==0) //Important check
803  // printf("Error 0\n");
804 
805  for (; pMultiplicity != pMultiplicityEnd; ++pMultiplicity, ++ppLocalSpecies, ++ppModelSpecies)
806  {
807  **ppLocalSpecies += *pMultiplicity;
808  **ppModelSpecies = **ppLocalSpecies;
809  }
810 
811  // Calculate all values which depend on the firing reaction
812  std::vector< Refresh * >::const_iterator itCalcualtion = Dependencies.mCalculations.begin();
813  std::vector< Refresh * >::const_iterator endCalcualtion = Dependencies.mCalculations.end();
814 
815  while (itCalcualtion != endCalcualtion)
816  {
817  (**itCalcualtion++)();
818  }
819 
820  // calculate the propensities which depend on the firing reaction
821  size_t * pDependentReaction = Dependencies.mDependentReactions.array();
822  size_t * endDependentReactions = pDependentReaction + Dependencies.mDependentReactions.size();
823 
824  for (; pDependentReaction != endDependentReactions; ++pDependentReaction)
825  {
826  calculateAmu(*pDependentReaction);
827  }
828 
829  // calculate the total propensity
830  pAmu = mAmu.array();
831 
832  mA0 = 0.0;
833 
834  for (; pAmu != endAmu; ++pAmu)
835  {
836  mA0 += *pAmu;
837  }
838 
840 
841  return mNextReactionTime - curTime;
842 }
843 
844 const C_FLOAT64 & CTrajAdaptiveSA::calculateAmu(const size_t & index)
845 {
846  const CReactionDependencies & Dependencies = mReactionDependencies[index];
847  C_FLOAT64 & Amu = mAmu[index];
848 
849  Amu = *Dependencies.mpParticleFlux;
850 
851  if (Amu < 0.0)
852  {
853  // TODO CRITICAL Create a warning message
854  Amu = 0.0;
855  }
856 
857  if (!mDoCorrection)
858  {
859  return Amu;
860  }
861 
862  C_FLOAT64 SubstrateMultiplier = 1.0;
863  C_FLOAT64 SubstrateDevisor = 1.0;
864  C_FLOAT64 Multiplicity;
865  C_FLOAT64 LowerBound;
866  C_FLOAT64 Number;
867 
868  bool ApplyCorrection = false;
869 
870  const C_FLOAT64 * pMultiplicity = Dependencies.mSubstrateMultiplier.array();
871  const C_FLOAT64 * pEndMultiplicity = pMultiplicity + Dependencies.mSubstrateMultiplier.size();
872  C_FLOAT64 *const* ppLocalSubstrate = Dependencies.mMethodSubstrates.array();
873  C_FLOAT64 *const* ppModelSubstrate = Dependencies.mModelSubstrates.array();
874 
875  for (; pMultiplicity != pEndMultiplicity; ++pMultiplicity, ++ppLocalSubstrate, ++ppModelSubstrate)
876  {
877  Multiplicity = *pMultiplicity;
878 
879  // TODO We should check the error introduced through rounding.
880  **ppLocalSubstrate = floor(**ppModelSubstrate + 0.5);
881 
882  if (Multiplicity > 1.01)
883  {
884  ApplyCorrection = true;
885 
886  Number = **ppLocalSubstrate;
887 
888  LowerBound = Number - Multiplicity;
889  SubstrateDevisor *= pow(Number, Multiplicity - 1.0); //optimization
890  Number -= 1.0;
891 
892  while (Number > LowerBound)
893  {
894  SubstrateMultiplier *= Number;
895  Number -= 1.0;
896  }
897  }
898  }
899 
900  // at least one substrate particle number is zero
901  if (SubstrateMultiplier < 0.5 || SubstrateDevisor < 0.5)
902  {
903  Amu = 0.0;
904  }
905  else if (ApplyCorrection)
906  {
907  Amu *= SubstrateMultiplier / SubstrateDevisor;
908  }
909 
910  return Amu;
911 }
#define UPPER_LIMIT
#define SSA_MULTIPLE
#define pdelete(p)
Definition: copasi.h:215
CVector< size_t > mMaxReactionFiring
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
virtual size_t size() const
virtual bool elevateChildren()
size_t mNumReactionSpecies
std::vector< Refresh * > mCalculations
iterator begin()
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
Definition: CState.h:305
CRandom * mpRandomGenerator
CTrajectoryProblem * mpProblem
CVector< C_FLOAT64 * > mMethodSubstrates
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
size_t mFirstReactionSpeciesIndex
#define C_INVALID_INDEX
Definition: copasi.h:222
CVector< C_FLOAT64 > mAvgDX
size_t getIndex(const CModelEntity *entity) const
Definition: CState.cpp:231
#define C_INT32
Definition: copasi.h:90
CVector< const C_FLOAT64 * > mPartitionedAmu
Definition: CMetab.h:178
virtual bool isValidProblem(const CCopasiProblem *pProblem)
const C_FLOAT64 & calculateAmu(const size_t &index)
CVector< C_FLOAT64 > mReactionFiring
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
size_t getTotSteps() const
Definition: CModel.cpp:1136
C_FLOAT64 mSSAStepCounter
CReactionDependencies & operator=(const CReactionDependencies &rhs)
#define MCTrajectoryMethod
CVector< C_FLOAT64 > mSpeciesAfterTau
iterator end()
std::vector< CType * >::const_iterator const_iterator
Definition: CCopasiVector.h:57
virtual C_FLOAT64 getRandomPoisson(const C_FLOAT64 &mean)
Definition: CRandom.cpp:314
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
std::vector< CReactionDependencies > mReactionDependencies
CModelEntity ** endFixed()
Definition: CState.cpp:213
void setState(const CState &state)
Definition: CModel.cpp:1785
const C_FLOAT64 & getDuration() const
static CTrajAdaptiveSA * createTauLeapMethod()
const Value & getValue() const
C_FLOAT64 doSingleTauLeapStep(const C_FLOAT64 &curTime, const C_FLOAT64 &endTime)
unsigned C_INT32 * pUINT
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
size_t size() const
Definition: CVector.h:100
void setTime(const C_FLOAT64 &time)
Definition: CModel.cpp:1187
unsigned C_INT32 mMaxSteps
virtual Status step(const double &deltaT)
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
CVector< C_FLOAT64 > mAmu
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
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 * getValuePointer() const
#define SSA_UPPER_NUM
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
#define EPS
const CState & getState() const
Definition: CModel.cpp:1771
virtual bool isValidProblem(const CCopasiProblem *pProblem)
const CModelEntity::Status & getStatus() const
CVector< const CReactionDependencies * > mPartitionedDependencies
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
C_FLOAT64 mNextReactionTime
virtual void start(const CState *initialState)
C_FLOAT64 doSingleSSAStep(const C_FLOAT64 &curTime, const C_FLOAT64 &endTime)
CVector< C_FLOAT64 > mSigDX
CTrajAdaptiveSA(const CCopasiContainer *pParent=NULL)
std::string suitableForStochasticSimulation() const
Definition: CModel.cpp:3901
CModel * getModel() const
CModelEntity ** endIndependent()
Definition: CState.cpp:209
CVector< C_FLOAT64 * > mPartitionedReactionFiring
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
CCopasiObject * getValueReference() const
C_FLOAT64 * mpMethodSpecies
#define max(a, b)
Definition: f2c.h:176