COPASI API  4.16.103
CStochDirectMethod.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) 2002 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #ifdef WIN32
16 # pragma warning (disable: 4786)
17 # pragma warning (disable: 4243)
18 // warning C4355: 'this' : used in base member initializer list
19 # pragma warning (disable: 4355)
20 #endif // WIN32
21 
22 #include <limits.h>
23 
24 #include <vector>
25 #include <numeric>
26 #include <limits>
27 #include <set>
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <string.h>
31 #include <cmath>
32 
33 #include "copasi.h"
34 #include "CStochDirectMethod.h"
36 #include "function/CFunction.h"
38 #include "CTrajectoryMethod.h"
39 #include "CTrajectoryProblem.h"
40 #include "model/CState.h"
41 #include "model/CCompartment.h"
42 #include "model/CModel.h"
43 
45  mSpeciesMultiplier(0),
46  mMethodSpecies(0),
47  mModelSpecies(0),
48  mCalculations(),
49  mDependentReactions(0),
50  mSubstrateMultiplier(0),
51  mMethodSubstrates(0),
52  mModelSubstrates(0),
53  mpParticleFlux(NULL)
54 {}
55 
57  mSpeciesMultiplier(src.mSpeciesMultiplier),
58  mMethodSpecies(src.mMethodSpecies),
59  mModelSpecies(src.mModelSpecies),
60  mCalculations(src.mCalculations),
61  mDependentReactions(src.mDependentReactions),
62  mSubstrateMultiplier(src.mSubstrateMultiplier),
63  mMethodSubstrates(src.mMethodSubstrates),
64  mModelSubstrates(src.mModelSubstrates),
65  mpParticleFlux(src.mpParticleFlux)
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 
83  return * this;
84 }
85 
88  mpRandomGenerator(CRandom::createGenerator(CRandom::mt19937)),
89  mpModel(NULL),
90  mNumReactions(0),
91  mMaxSteps(1000000),
92  mNextReactionTime(0.0),
94  mDoCorrection(true),
95  mAmu(0),
96  mA0(0.0),
97  mMethodState(),
99  mMaxStepsReached(false)
100 {
102 }
103 
105  const CCopasiContainer * pParent):
106  CTrajectoryMethod(src, pParent),
107  mpRandomGenerator(CRandom::createGenerator(CRandom::mt19937)),
108  mpModel(src.mpModel),
109  mNumReactions(src.mNumReactions),
110  mMaxSteps(src.mMaxSteps),
111  mNextReactionTime(src.mNextReactionTime),
112  mNextReactionIndex(src.mNextReactionIndex),
113  mDoCorrection(src.mDoCorrection),
114  mAmu(src.mAmu),
115  mA0(src.mA0),
116  mMethodState(src.mMethodState),
117  mReactionDependencies(src.mReactionDependencies),
118  mMaxStepsReached(src.mMaxStepsReached)
119 {
121 }
122 
124 {
126 }
127 
129 {
130  assertParameter("Max Internal Steps", CCopasiParameter::INT, (C_INT32) 1000000);
131  assertParameter("Use Random Seed", CCopasiParameter::BOOL, false);
132  assertParameter("Random Seed", CCopasiParameter::UINT, (unsigned C_INT32) 1);
133 }
134 
136 {
138  return true;
139 }
140 
142 {
143  // do several steps
144  C_FLOAT64 Time = mpCurrentState->getTime();
145  C_FLOAT64 EndTime = Time + deltaT;
146 
147  size_t Steps = 0;
148 
149  while (Time < EndTime)
150  {
151  mMethodState.setTime(Time);
152  Time += doSingleStep(Time, EndTime);
153 
154  if (++Steps > mMaxSteps)
155  {
157  }
158  }
159 
161  mpCurrentState->setTime(Time);
162 
163  return NORMAL;
164 }
165 
166 void CStochDirectMethod::start(const CState * initialState)
167 {
168  /* get configuration data */
169  mMaxSteps = * getValue("Max Internal Steps").pINT;
170 
171  bool useRandomSeed = * getValue("Use Random Seed").pBOOL;
172  unsigned C_INT32 randomSeed = * getValue("Random Seed").pUINT;
173 
174  if (useRandomSeed) mpRandomGenerator->initialize(randomSeed);
175 
176  //mpCurrentState is initialized. This state is not used internally in the
177  //stochastic solver, but it is used for returning the result after each step.
178  *mpCurrentState = *initialState;
179 
181  assert(mpModel);
182 
184  mDoCorrection = true;
185  else
186  mDoCorrection = false;
187 
188  //initialize the vector of ints that contains the particle numbers
189  //for the discrete simulation. This also floors all particle numbers in the model.
190 
192 
195  mAmu = 0.0;
196 
197  // Create a local copy of the state where the particle number are rounded to integers
199 
200  const CStateTemplate & StateTemplate = mpModel->getStateTemplate();
201 
202  CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
203  CModelEntity *const* endEntity = StateTemplate.endFixed();
205 
206  for (; ppEntity != endEntity; ++ppEntity, ++pValue)
207  {
208  if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
209  {
210  *pValue = floor(*pValue + 0.5);
211  }
212  }
213 
214  // Update the model state so that the species are all represented by integers.
216  mpModel->updateSimulatedValues(false); //for assignments
217 
218  // TODO CRITICAL Handle species of type ASSIGNMENT.
219  // These need to be checked whether they are sufficiently close to an integer
220 
221  C_FLOAT64 * pMethodStateValue = mMethodState.beginIndependent() - 1;
222 
223  // Build the reaction dependencies
224  size_t NumReactions = 0;
225 
228  std::vector< CReactionDependencies >::iterator itDependencies = mReactionDependencies.begin();
229 
230  for (; it != end; ++it)
231  {
232  const CCopasiVector<CChemEqElement> & Balances = (*it)->getChemEq().getBalances();
233  const CCopasiVector<CChemEqElement> & Substrates = (*it)->getChemEq().getSubstrates();
234 
235  // This reactions does not change anything we ignore it
236  if (Balances.size() == 0 && Substrates.size() == 0)
237  {
238  continue;
239  }
240 
241  itDependencies->mpParticleFlux = (C_FLOAT64 *)(*it)->getParticleFluxReference()->getValuePointer();
242 
243  itDependencies->mSpeciesMultiplier.resize(Balances.size());
244  itDependencies->mMethodSpecies.resize(Balances.size());
245  itDependencies->mModelSpecies.resize(Balances.size());
246 
247  std::set< const CCopasiObject * > changed;
248 
249  // The time is always updated
250  changed.insert(mpModel->getValueReference());
251 
254 
255  size_t Index = 0;
256 
257  for (; itBalance != endBalance; ++itBalance)
258  {
259  const CMetab * pMetab = (*itBalance)->getMetabolite();
260 
261  if (pMetab->getStatus() == CModelEntity::REACTIONS)
262  {
263  itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
264  itDependencies->mMethodSpecies[Index] = pMethodStateValue + StateTemplate.getIndex(pMetab);
265  itDependencies->mModelSpecies[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
266 
267  changed.insert(pMetab->getValueReference());
268 
269  Index++;
270  }
271  }
272 
273  // Correct allocation for metabolites which are not determined by reactions
274  itDependencies->mSpeciesMultiplier.resize(Index, true);
275  itDependencies->mMethodSpecies.resize(Index, true);
276  itDependencies->mModelSpecies.resize(Index, true);
277 
278  itDependencies->mSubstrateMultiplier.resize(Substrates.size());
279  itDependencies->mMethodSubstrates.resize(Substrates.size());
280  itDependencies->mModelSubstrates.resize(Substrates.size());
281 
282  CCopasiVector< CChemEqElement >::const_iterator itSubstrate = Substrates.begin();
283  CCopasiVector< CChemEqElement >::const_iterator endSubstrate = Substrates.end();
284 
285  Index = 0;
286 
287  for (; itSubstrate != endSubstrate; ++itSubstrate, ++Index)
288  {
289  const CMetab * pMetab = (*itSubstrate)->getMetabolite();
290 
291  itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
292  itDependencies->mMethodSubstrates[Index] = pMethodStateValue + StateTemplate.getIndex(pMetab);
293  itDependencies->mModelSubstrates[Index] = (C_FLOAT64 *) pMetab->getValueReference()->getValuePointer();
294  }
295 
296  std::set< const CCopasiObject * > dependend;
297 
299  itDependencies->mDependentReactions.resize(mNumReactions);
300 
301  Index = 0;
302  size_t Count = 0;
303 
304  for (; itReaction != end; ++itReaction, ++Index)
305  {
306  if ((*itReaction)->getParticleFluxReference()->dependsOn(changed))
307  {
308  dependend.insert((*itReaction)->getParticleFluxReference());
309  itDependencies->mDependentReactions[Count] = Index;
310 
311  Count++;
312  }
313  }
314 
315  itDependencies->mDependentReactions.resize(Count, true);
316 
317  itDependencies->mCalculations = CCopasiObject::buildUpdateSequence(dependend, changed);
318 
319  ++itDependencies;
320  ++NumReactions;
321  }
322 
323  mNumReactions = NumReactions;
324 
326  mAmu.resize(mNumReactions, true);
327 
328  size_t i;
329 
330  for (i = 0; i < mNumReactions; i++)
331  {
332  calculateAmu(i);
333  }
334 
335  mA0 = 0;
336 
337  for (i = 0; i < mNumReactions; i++)
338  {
339  mA0 += mAmu[i];
340  }
341 
342  mMaxStepsReached = false;
343 
346 
347  return;
348 }
349 
350 // virtual
352 {
353  if (!CTrajectoryMethod::isValidProblem(pProblem)) return false;
354 
355  const CTrajectoryProblem * pTP = dynamic_cast<const CTrajectoryProblem *>(pProblem);
356 
357  if (pTP->getDuration() < 0.0)
358  {
359  //back integration not possible
361  return false;
362  }
363 
364  if (pTP->getModel()->getTotSteps() < 1)
365  {
366  //at least one reaction necessary
368  return false;
369  }
370 
371  // check for ODEs
372  const CStateTemplate & StateTemplate = pTP->getModel()->getStateTemplate();
373  CModelEntity *const* ppEntity = StateTemplate.beginIndependent();
374  CModelEntity *const* ppEntityEnd = StateTemplate.endIndependent();
375 
376  for (; ppEntity != ppEntityEnd; ++ppEntity)
377  {
378  if ((*ppEntity)->getStatus() == CModelEntity::ODE)
379  {
380  if (dynamic_cast<const CModelValue *>(*ppEntity) != NULL)
381  {
382  // global quantity ode rule found
384  return false;
385  }
386  else if (dynamic_cast<const CCompartment *>(*ppEntity) != NULL)
387  {
388  // compartment ode rule found
390  return false;
391  }
392  else
393  {
394  // species ode rule found
396  return false;
397  }
398  }
399  }
400 
401  //TODO: rewrite CModel::suitableForStochasticSimulation() to use
402  // CCopasiMessage
403  std::string message = pTP->getModel()->suitableForStochasticSimulation();
404 
405  if (message != "")
406  {
407  //model not suitable, message describes the problem
408  CCopasiMessage(CCopasiMessage::ERROR, message.c_str());
409  return false;
410  }
411 
412  if (* getValue("Max Internal Steps").pINT <= 0)
413  {
414  //max steps should be at least 1
416  return false;
417  }
418 
419  //events are not supported at the moment
420  if (pTP->getModel()->getEvents().size() > 0)
421  {
423  return false;
424  }
425 
426  return true;
427 }
428 
430 {
432  {
433  if (mA0 == 0)
434  {
435  return endTime - curTime;
436  }
437 
438  // We need to throw an exception if mA0 is NaN
439  if (isnan(mA0))
440  {
442  }
443 
444  mNextReactionTime = curTime - log(mpRandomGenerator->getRandomOO()) / mA0;
445 
446  // We are sure that we have at least 1 reaction
447  mNextReactionIndex = 0;
448 
449  C_FLOAT64 sum = 0.0;
451 
452  C_FLOAT64 * pAmu = mAmu.array();
453  C_FLOAT64 * endAmu = pAmu + mNumReactions;
454 
455  for (; (sum < rand) && (pAmu != endAmu); ++pAmu, ++mNextReactionIndex)
456  {
457  sum += *pAmu;
458  }
459 
460  mNextReactionIndex--;
461  }
462 
463  if (mNextReactionTime >= endTime)
464  {
465  return endTime - curTime;
466  }
467 
468  // Update the model time (for explicitly time dependent models)
470 
472 
473  // Update the method internal and model species numbers
474  C_FLOAT64 ** ppModelSpecies = Dependencies.mModelSpecies.array();
475  C_FLOAT64 ** ppLocalSpecies = Dependencies.mMethodSpecies.array();
476  C_FLOAT64 * pMultiplier = Dependencies.mSpeciesMultiplier.array();
477  C_FLOAT64 * endMultiplier = pMultiplier + Dependencies.mSpeciesMultiplier.size();
478 
479  for (; pMultiplier != endMultiplier; ++pMultiplier, ++ppLocalSpecies, ++ppModelSpecies)
480  {
481  **ppLocalSpecies += *pMultiplier;
482  **ppModelSpecies = **ppLocalSpecies;
483  }
484 
485  // Calculate all values which depend on the firing reaction
486  std::vector< Refresh * >::const_iterator itCalcualtion = Dependencies.mCalculations.begin();
487  std::vector< Refresh * >::const_iterator endCalcualtion = Dependencies.mCalculations.end();
488 
489  while (itCalcualtion != endCalcualtion)
490  {
491  (**itCalcualtion++)();
492  }
493 
494  // We do not need to update the the method state since the only independent state
495  // values are species of type reaction which are all controlled by the method.
496 
497  // calculate the propensities which depend on the firing reaction
498  size_t * pDependentReaction = Dependencies.mDependentReactions.array();
499  size_t * endDependentReactions = pDependentReaction + Dependencies.mDependentReactions.size();
500 
501  for (; pDependentReaction != endDependentReactions; ++pDependentReaction)
502  {
503  calculateAmu(*pDependentReaction);
504  }
505 
506  // calculate the total propensity
507  C_FLOAT64 * pAmu = mAmu.array();
508  C_FLOAT64 * endAmu = pAmu + mNumReactions;
509 
510  mA0 = 0.0;
511 
512  for (; pAmu != endAmu; ++pAmu)
513  {
514  mA0 += *pAmu;
515  }
516 
518 
519  return mNextReactionTime - curTime;
520 }
521 
522 void CStochDirectMethod::calculateAmu(const size_t & index)
523 {
524  CReactionDependencies & Dependencies = mReactionDependencies[index];
525  C_FLOAT64 & Amu = mAmu[index];
526 
527  Amu = *Dependencies.mpParticleFlux;
528 
529  if (Amu < 0.0)
530  {
531  // TODO CRITICAL Create a warning message
532  Amu = 0.0;
533  }
534 
535  if (!mDoCorrection)
536  {
537  return;
538  }
539 
540  C_FLOAT64 SubstrateMultiplier = 1.0;
541  C_FLOAT64 SubstrateDevisor = 1.0;
542  C_FLOAT64 Multiplicity;
543  C_FLOAT64 LowerBound;
544  C_FLOAT64 Number;
545 
546  bool ApplyCorrection = false;
547 
548  C_FLOAT64 * pMultiplier = Dependencies.mSubstrateMultiplier.array();
549  C_FLOAT64 * endMultiplier = pMultiplier + Dependencies.mSubstrateMultiplier.size();
550  C_FLOAT64 ** ppLocalSubstrate = Dependencies.mMethodSubstrates.array();
551  C_FLOAT64 ** ppModelSubstrate = Dependencies.mModelSubstrates.array();
552 
553  for (; pMultiplier != endMultiplier; ++pMultiplier, ++ppLocalSubstrate, ++ppModelSubstrate)
554  {
555  Multiplicity = *pMultiplier;
556 
557  // TODO We should check the error introduced through rounding.
558  **ppLocalSubstrate = floor(**ppModelSubstrate + 0.5);
559 
560  if (Multiplicity > 1.01)
561  {
562  ApplyCorrection = true;
563 
564  Number = **ppLocalSubstrate;
565 
566  LowerBound = Number - Multiplicity;
567  SubstrateDevisor *= pow(Number, Multiplicity - 1.0); //optimization
568  Number -= 1.0;
569 
570  while (Number > LowerBound)
571  {
572  SubstrateMultiplier *= Number;
573  Number -= 1.0;
574  }
575  }
576  }
577 
578  // at least one substrate particle number is zero
579  if (SubstrateMultiplier < 0.5 || SubstrateDevisor < 0.5)
580  {
581  Amu = 0.0;
582  return;
583  }
584 
585  // rate_factor is the rate function divided by substrate_factor.
586  // It would be more efficient if this was generated directly, since in effect we
587  // are multiplying and then dividing by the same thing (substrate_factor)!
588  //C_FLOAT64 rate_factor = mpModel->getReactions()[index]->calculateParticleFlux();
589  if (ApplyCorrection)
590  {
591  Amu *= SubstrateMultiplier / SubstrateDevisor;
592  }
593 
594  return;
595 }
virtual bool isValidProblem(const CCopasiProblem *pProblem)
CReactionDependencies & operator=(const CReactionDependencies &rhs)
#define pdelete(p)
Definition: copasi.h:215
CCopasiVectorN< CEvent > & getEvents()
Definition: CModel.cpp:1110
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
virtual size_t size() const
void calculateAmu(const size_t &index)
CVector< C_FLOAT64 > mAmu
iterator begin()
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
Definition: CState.h:305
std::vector< CReactionDependencies > mReactionDependencies
CTrajectoryProblem * mpProblem
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
#define C_INT32
Definition: copasi.h:90
Definition: CMetab.h:178
virtual bool isValidProblem(const CCopasiProblem *pProblem)
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
size_t getTotSteps() const
Definition: CModel.cpp:1136
#define MCTrajectoryMethod
iterator end()
std::vector< CType * >::const_iterator const_iterator
Definition: CCopasiVector.h:57
CModelEntity ** endFixed()
Definition: CState.cpp:213
void setState(const CState &state)
Definition: CModel.cpp:1785
const C_FLOAT64 & getDuration() const
const Value & getValue() const
virtual void start(const CState *initialState)
unsigned C_INT32 * pUINT
C_FLOAT64 doSingleStep(const C_FLOAT64 &curTime, const C_FLOAT64 &endTime)
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
unsigned C_INT32 mMaxSteps
size_t size() const
Definition: CVector.h:100
CStochDirectMethod(const CCopasiContainer *pParent=NULL)
void setTime(const C_FLOAT64 &time)
Definition: CModel.cpp:1187
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
virtual bool elevateChildren()
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
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
const CState & getState() const
Definition: CModel.cpp:1771
const CModelEntity::Status & getStatus() const
std::string suitableForStochasticSimulation() const
Definition: CModel.cpp:3901
virtual Status step(const double &deltaT)
CModel * getModel() const
CModelEntity ** endIndependent()
Definition: CState.cpp:209
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
CCopasiObject * getValueReference() const