7 # pragma warning (disable: 4786)
8 # pragma warning (disable: 4243)
10 # pragma warning (disable: 4355)
33 mMethodSpeciesIndex(0),
34 mSpeciesMultiplier(0),
38 mDependentReactions(0),
39 mSubstrateMultiplier(0),
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)
117 mMaxReactionFiring(src.mMaxReactionFiring),
118 mReactionFiring(src.mReactionFiring),
119 mPartitionedReactionFiring(src.mPartitionedReactionFiring),
122 mpMethodSpecies(src.mpMethodSpecies),
123 mSpeciesAfterTau(src.mSpeciesAfterTau),
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),
133 mPartitionedAmu(src.mPartitionedAmu),
134 mMethodState(src.mMethodState),
135 mReactionDependencies(src.mReactionDependencies),
136 mPartitionedDependencies(src.mPartitionedDependencies),
138 mMaxStepsReached(src.mMaxStepsReached)
170 while (Time < EndTime)
255 for (; ppEntity != endEntity; ++ppEntity, ++pValue, ++Index)
257 if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
259 *pValue = floor(*pValue + 0.5);
262 (*ppEntity)->isUsed())
264 *pSpeciesAfterTau = *pValue;
286 size_t NumReactions = 0;
292 for (; it != end; ++it)
298 if (Balances.
size() == 0 && Substrates.
size() == 0)
303 itDependencies->mpParticleFlux = (
C_FLOAT64 *)(*it)->getParticleFluxReference()->getValuePointer();
305 itDependencies->mMethodSpeciesIndex.resize(Balances.
size());
306 itDependencies->mSpeciesMultiplier.resize(Balances.
size());
307 itDependencies->mMethodSpecies.resize(Balances.
size());
308 itDependencies->mModelSpecies.resize(Balances.
size());
310 std::set< const CCopasiObject * > changed;
320 for (; itBalance != endBalance; ++itBalance)
322 const CMetab * pMetab = (*itBalance)->getMetabolite();
327 itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
328 itDependencies->mMethodSpecies[Index] = pMethodStateValue + StateTemplate.
getIndex(pMetab);
338 itDependencies->mMethodSpeciesIndex.resize(Index,
true);
339 itDependencies->mSpeciesMultiplier.resize(Index,
true);
340 itDependencies->mMethodSpecies.resize(Index,
true);
341 itDependencies->mModelSpecies.resize(Index,
true);
343 itDependencies->mSubstrateMultiplier.resize(Substrates.
size());
344 itDependencies->mMethodSubstrates.resize(Substrates.
size());
345 itDependencies->mModelSubstrates.resize(Substrates.
size());
352 for (; itSubstrate != endSubstrate; ++itSubstrate, ++Index)
354 const CMetab * pMetab = (*itSubstrate)->getMetabolite();
356 itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
357 itDependencies->mMethodSubstrates[Index] = pMethodStateValue + StateTemplate.
getIndex(pMetab);
361 std::set< const CCopasiObject * > dependend;
369 for (; itReaction != end; ++itReaction, ++Index)
371 if ((*itReaction)->getParticleFluxReference()->dependsOn(changed))
373 dependend.insert((*itReaction)->getParticleFluxReference());
374 itDependencies->mDependentReactions[Count] = Index;
380 itDependencies->mDependentReactions.resize(Count,
true);
436 for (; ppEntity != ppEntityEnd; ++ppEntity)
440 if (dynamic_cast<const CModelValue *>(*ppEntity) != NULL)
446 else if (dynamic_cast<const CCompartment *>(*ppEntity) != NULL)
488 for (; itDependencies != endDependencies; ++itDependencies, ++pMaxReactionFiring)
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();
496 for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, ppModelSpecies++)
498 if (*pMultiplicity < 0)
502 if ((TmpMax = (
size_t) fabs(**ppModelSpecies / *pMultiplicity)) < *pMaxReactionFiring)
504 *pMaxReactionFiring = TmpMax;
510 size_t NonCriticalReactions = 0;
519 for (; itDependencies != endDependencies; ++itDependencies, ++pAmu, ++pMaxReactionFiring, ++pReactionFiring)
521 *pReactionFiring = 0;
529 NonCriticalReactions++;
545 const C_FLOAT64 ** ppOrderedAmuEnd = ppOrderedAmu + NonCriticalReactions;
547 for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedReactions, ++ppOrderedAmu)
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();
553 for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
555 mAvgDX[*pIndex] += **ppOrderedAmu * *pMultiplicity;
556 mSigDX[*pIndex] += **ppOrderedAmu * *pMultiplicity * *pMultiplicity;
562 if (NonCriticalReactions == 0)
569 t1 = t2 = std::numeric_limits< C_FLOAT64 >::infinity();
581 if (t3 != 0 && t5 < t1) t1 = t5;
583 if (t4 != 0 && t6 < t2) t2 = t6;
593 MaxTau == std::numeric_limits< C_FLOAT64 >::infinity())
606 for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedAmu)
608 AmuCritical += **ppOrderedAmu;
613 if (NonCriticalReactions ==
mNumReactions || AmuCritical == 0)
614 CriticalReactionTau = std::numeric_limits< C_FLOAT64 >::infinity();
618 bool isUpdated =
false;
626 if (CriticalReactionTau < Tau || Tau == 0) Tau = CriticalReactionTau;
628 if ((endTime - curTime) < Tau) Tau = (endTime - curTime);
636 for (; ppOrderedAmu != ppOrderedAmuEnd; ++ppOrderedAmu, ++ppOrderedReactionFiring)
642 else if (Lambda > 2.0e9)
650 if (CriticalReactionTau <= Tau)
659 CriticalReactionIndex = NonCriticalReactions;
661 for (; (sum < rand) && (ppOrderedAmu != ppOrderedAmuEnd); ++ppOrderedAmu, ++CriticalReactionIndex)
663 sum += **ppOrderedAmu;
666 CriticalReactionIndex--;
673 for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau, ++pMethodSpecies)
675 *pSpeciesAfterTau = *pMethodSpecies;
682 for (; ppOrderedReactions != ppOrderedReactionsEnd; ++ppOrderedReactions, ++ppOrderedReactionFiring)
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();
688 for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
693 **ppOrderedReactionFiring = 0.0;
704 for (; pMultiplicity != pMultiplicityEnd; pMultiplicity++, pIndex++)
712 for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau)
714 if (*pSpeciesAfterTau < 0)
718 if (pSpeciesAfterTau != pSpeciesAfterTauEnd)
721 MaxTau = MaxTau / 2.0;
730 for (; pSpeciesAfterTau != pSpeciesAfterTauEnd; ++pSpeciesAfterTau, ++pMethodSpecies)
732 *pMethodSpecies = *pSpeciesAfterTau;
759 return endTime - curTime;
773 return endTime - curTime;
805 for (; pMultiplicity != pMultiplicityEnd; ++pMultiplicity, ++ppLocalSpecies, ++ppModelSpecies)
807 **ppLocalSpecies += *pMultiplicity;
808 **ppModelSpecies = **ppLocalSpecies;
812 std::vector< Refresh * >::const_iterator itCalcualtion = Dependencies.
mCalculations.begin();
813 std::vector< Refresh * >::const_iterator endCalcualtion = Dependencies.
mCalculations.end();
815 while (itCalcualtion != endCalcualtion)
817 (**itCalcualtion++)();
824 for (; pDependentReaction != endDependentReactions; ++pDependentReaction)
834 for (; pAmu != endAmu; ++pAmu)
868 bool ApplyCorrection =
false;
875 for (; pMultiplicity != pEndMultiplicity; ++pMultiplicity, ++ppLocalSubstrate, ++ppModelSubstrate)
877 Multiplicity = *pMultiplicity;
880 **ppLocalSubstrate = floor(**ppModelSubstrate + 0.5);
882 if (Multiplicity > 1.01)
884 ApplyCorrection =
true;
886 Number = **ppLocalSubstrate;
888 LowerBound = Number - Multiplicity;
889 SubstrateDevisor *= pow(Number, Multiplicity - 1.0);
892 while (Number > LowerBound)
894 SubstrateMultiplier *= Number;
901 if (SubstrateMultiplier < 0.5 || SubstrateDevisor < 0.5)
905 else if (ApplyCorrection)
907 Amu *= SubstrateMultiplier / SubstrateDevisor;
CVector< size_t > mMaxReactionFiring
virtual C_FLOAT64 getRandomOO()
virtual size_t size() const
virtual bool elevateChildren()
CVector< C_FLOAT64 * > mModelSpecies
size_t mNumReactionSpecies
std::vector< Refresh * > mCalculations
void updateSimulatedValues(const bool &updateMoieties)
CRandom * mpRandomGenerator
CTrajectoryProblem * mpProblem
CVector< C_FLOAT64 * > mMethodSubstrates
void resize(size_t size, const bool ©=false)
size_t mFirstReactionSpeciesIndex
CVector< C_FLOAT64 > mAvgDX
size_t getIndex(const CModelEntity *entity) const
void initializeParameter()
CVector< C_FLOAT64 * > mModelSubstrates
CVector< const C_FLOAT64 * > mPartitionedAmu
virtual bool isValidProblem(const CCopasiProblem *pProblem)
CVector< C_FLOAT64 * > mMethodSpecies
const C_FLOAT64 & calculateAmu(const size_t &index)
CVector< C_FLOAT64 > mReactionFiring
void setTime(const C_FLOAT64 &time)
size_t mNextReactionIndex
size_t getTotSteps() const
C_FLOAT64 mSSAStepCounter
CVector< size_t > mMethodSpeciesIndex
CReactionDependencies & operator=(const CReactionDependencies &rhs)
#define MCTrajectoryMethod
CVector< C_FLOAT64 > mSpeciesAfterTau
std::vector< CType * >::const_iterator const_iterator
virtual C_FLOAT64 getRandomPoisson(const C_FLOAT64 &mean)
size_t getNumDependentReactionMetabs() const
std::vector< CReactionDependencies > mReactionDependencies
CModelEntity ** endFixed()
void setState(const CState &state)
const C_FLOAT64 & getDuration() const
static CTrajAdaptiveSA * createTauLeapMethod()
const Value & getValue() const
C_FLOAT64 doSingleTauLeapStep(const C_FLOAT64 &curTime, const C_FLOAT64 &endTime)
CModelEntity ** beginIndependent()
size_t getNumIndependentReactionMetabs() const
void setTime(const C_FLOAT64 &time)
unsigned C_INT32 mMaxSteps
virtual Status step(const double &deltaT)
const C_FLOAT64 & getTime() const
CVector< C_FLOAT64 > mAmu
CVector< C_FLOAT64 > mSpeciesMultiplier
const CStateTemplate & getStateTemplate() const
const ModelType & getModelType() const
CCopasiVectorNS< CReaction > & getReactions()
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)
CVector< C_FLOAT64 > mSubstrateMultiplier
const CState & getState() const
virtual bool isValidProblem(const CCopasiProblem *pProblem)
const CModelEntity::Status & getStatus() const
CVector< const CReactionDependencies * > mPartitionedDependencies
const CCopasiVector< CMetab > & getMetabolitesX() const
CVector< size_t > mDependentReactions
C_FLOAT64 * mpParticleFlux
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
CModel * getModel() const
CModelEntity ** endIndependent()
CVector< C_FLOAT64 * > mPartitionedReactionFiring
C_FLOAT64 * beginIndependent()
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
CCopasiObject * getValueReference() const
C_FLOAT64 * mpMethodSpecies