61 mSpeciesMultiplier(0),
62 mMethodSpeciesIndex(0),
65 mSubstrateMultiplier(0),
72 mSpeciesMultiplier(src.mSpeciesMultiplier),
73 mMethodSpeciesIndex(src.mMethodSpeciesIndex),
74 mMethodSpecies(src.mMethodSpecies),
75 mModelSpecies(src.mModelSpecies),
76 mSubstrateMultiplier(src.mSubstrateMultiplier),
77 mMethodSubstrates(src.mMethodSubstrates),
78 mModelSubstrates(src.mModelSubstrates),
79 mpParticleFlux(src.mpParticleFlux)
115 mReactionDependencies()
144 if ((pParm =
getParameter(
"TAULEAP.UseRandomSeed")) != NULL)
150 if ((pParm =
getParameter(
"TAULEAP.RandomSeed")) != NULL)
172 while (Time < EndTime)
244 for (; ppEntity != endEntity; ++ppEntity, ++pValue, ++Index)
246 if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
248 *pValue = floor(*pValue + 0.5);
265 size_t NumReactions = 0;
271 for (; it != end; ++it)
277 if (Balances.
size() == 0 && Substrates.
size() == 0)
282 itDependencies->mpParticleFlux = (
C_FLOAT64 *)(*it)->getParticleFluxReference()->getValuePointer();
284 itDependencies->mMethodSpeciesIndex.resize(Balances.
size());
285 itDependencies->mSpeciesMultiplier.resize(Balances.
size());
286 itDependencies->mMethodSpecies.resize(Balances.
size());
287 itDependencies->mModelSpecies.resize(Balances.
size());
294 for (; itBalance != endBalance; ++itBalance)
296 const CMetab * pMetab = (*itBalance)->getMetabolite();
301 itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
302 itDependencies->mMethodSpecies[Index] = pMethodStateValue + StateTemplate.
getIndex(pMetab);
310 itDependencies->mMethodSpeciesIndex.resize(Index,
true);
311 itDependencies->mSpeciesMultiplier.resize(Index,
true);
312 itDependencies->mMethodSpecies.resize(Index,
true);
313 itDependencies->mModelSpecies.resize(Index,
true);
315 itDependencies->mSubstrateMultiplier.resize(Substrates.
size());
316 itDependencies->mMethodSubstrates.resize(Substrates.
size());
317 itDependencies->mModelSubstrates.resize(Substrates.
size());
324 for (; itSubstrate != endSubstrate; ++itSubstrate, ++Index)
326 const CMetab * pMetab = (*itSubstrate)->getMetabolite();
328 itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
329 itDependencies->mMethodSubstrates[Index] = pMethodStateValue + StateTemplate.
getIndex(pMetab);
384 for (; pAmu != pAmuEnd; ++pAmu, ++itReaction)
386 const C_FLOAT64 * pMultiplicity = itReaction->mSpeciesMultiplier.array();
387 const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + itReaction->mSpeciesMultiplier.size();
388 const size_t * pIndex = itReaction->mMethodSpeciesIndex.array();
390 for (; pMultiplicity != pMultiplicityEnd; ++pMultiplicity, ++pIndex)
392 mAvgDX[*pIndex] += *pMultiplicity * *pAmu;
393 mSigDX[*pIndex] += *pMultiplicity * *pMultiplicity * *pAmu;
397 Tau1 = Tau2 = std::numeric_limits< C_FLOAT64 >::infinity();
404 for (; pNumber != pNumberEnd; ++pNumber, ++pAvgDX, ++pSigDX)
406 if ((Tmp =
mEpsilon * fabs(*pNumber)) < 1.0)
409 *pAvgDX = Tmp / fabs(*pAvgDX);
410 *pSigDX = (Tmp * Tmp) / fabs(*pSigDX);
428 for (; pAmu != pAmuEnd; ++pAmu, ++pK)
430 Lambda = *pAmu * Tau;
434 else if (Lambda > 2.0e9)
445 for (; pK != pKEnd; ++pK)
449 if (*pK < floor(*pK + 0.75))
497 bool ApplyCorrection =
false;
504 for (; pMultiplicity != pEndMultiplicity; ++pMultiplicity, ++ppLocalSubstrate, ++ppModelSubstrate)
506 Multiplicity = *pMultiplicity;
509 **ppLocalSubstrate = floor(**ppModelSubstrate + 0.5);
511 if (Multiplicity > 1.01)
513 ApplyCorrection =
true;
515 Number = **ppLocalSubstrate;
517 LowerBound = Number - Multiplicity;
518 SubstrateDevisor *= pow(Number, Multiplicity - 1.0);
521 while (Number > LowerBound)
523 SubstrateMultiplier *= Number;
530 if (SubstrateMultiplier < 0.5 || SubstrateDevisor < 0.5)
534 else if (ApplyCorrection)
536 Amu *= SubstrateMultiplier / SubstrateDevisor;
555 for (; pK != pKEnd; ++pK, ++itReaction)
557 const C_FLOAT64 * pMultiplicity = itReaction->mSpeciesMultiplier.array();
558 const C_FLOAT64 * pMultiplicityEnd = pMultiplicity + itReaction->mSpeciesMultiplier.size();
559 C_FLOAT64 *
const * pSpecies = itReaction->mMethodSpecies.array();
561 for (; pMultiplicity != pMultiplicityEnd; ++pMultiplicity, ++pSpecies)
563 **pSpecies += *pK * *pMultiplicity;
571 for (; pSpecies != pSpeciesEnd; ++pSpecies)
573 if (*pSpecies < -0.5)
608 for (i = 0; i < imax; ++i)
620 for (i = 0; i < imax; ++i)
632 for (i = 0; i < imax; ++i)
CCopasiVectorN< CEvent > & getEvents()
const CCopasiVector< CMetab > & getMetabolites() const
CVector< C_FLOAT64 > mSpeciesMultiplier
CVector< C_FLOAT64 * > mMethodSubstrates
const CCopasiVectorN< CModelValue > & getModelValues() const
virtual size_t size() const
void initializeParameter()
size_t getNumMetabs() const
void updateSimulatedValues(const bool &updateMoieties)
CTrajectoryProblem * mpProblem
void resize(size_t size, const bool ©=false)
void updatePropensities()
size_t mFirstReactionSpeciesIndex
CReactionDependencies & operator=(const CReactionDependencies &rhs)
size_t getIndex(const CModelEntity *entity) const
CVector< C_FLOAT64 * > mMethodSpecies
C_FLOAT64 doSingleStep(C_FLOAT64 ds)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
void setTime(const C_FLOAT64 &time)
const C_FLOAT64 & calculateAmu(const size_t &index)
size_t getTotSteps() const
bool removeParameter(const std::string &name)
#define MCTrajectoryMethod
std::vector< CType * >::const_iterator const_iterator
std::vector< CReactionDependencies > mReactionDependencies
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
virtual C_FLOAT64 getRandomPoisson(const C_FLOAT64 &mean)
virtual C_FLOAT64 getRandomCC()
C_FLOAT64 * mpParticleFlux
size_t getNumDependentReactionMetabs() const
CModelEntity ** endFixed()
void setState(const CState &state)
unsigned C_INT32 mMaxSteps
const C_FLOAT64 & getDuration() const
const Value & getValue() const
CVector< C_FLOAT64 > mAmu
bool setValue(const std::string &name, const CType &value)
virtual Status step(const double &deltaT)
CVector< size_t > mMethodSpeciesIndex
CModelEntity ** beginIndependent()
CCopasiParameter * getParameter(const std::string &name)
CTauLeapMethod(const CTauLeapMethod &src, const CCopasiContainer *pParent=NULL)
size_t getNumIndependentReactionMetabs() const
CRandom * mpRandomGenerator
CVector< C_FLOAT64 > mSigDX
CCopasiVectorNS< CCompartment > & getCompartments()
const C_FLOAT64 & getTime() const
CVector< C_FLOAT64 * > mModelSubstrates
const CStateTemplate & getStateTemplate() const
size_t getNumModelValues() const
const ModelType & getModelType() const
size_t mNumReactionSpecies
CCopasiVectorNS< CReaction > & getReactions()
virtual void * getValuePointer() const
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
CVector< C_FLOAT64 > mAvgDX
unsigned C_INT32 mRandomSeed
const CState & getState() const
const CModelEntity::Status & getStatus() const
virtual bool elevateChildren()
virtual bool isValidProblem(const CCopasiProblem *pProblem)
CVector< C_FLOAT64 > mSubstrateMultiplier
std::string suitableForStochasticSimulation() const
virtual void start(const CState *initialState)
CModel * getModel() const
C_FLOAT64 * beginIndependent()
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
CVector< C_FLOAT64 * > mModelSpecies
CCopasiObject * getValueReference() const