43 mSpeciesMultiplier(0),
47 mDependentReactions(0),
48 mSubstrateMultiplier(0),
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)
88 mSpeciesToReactions(),
92 mNumReactionSpecies(0),
93 mStochasticReactions(0),
94 mDeterministicReactions(0),
95 mStochasticSpecies(0),
96 mHasStochastic(false),
97 mHasDeterministic(false),
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)
126 mSpeciesToReactions = speciesToReactions;
127 mLowerThreshold = lowerThreshold;
128 mUpperThreshold = upperThreshold;
133 mStochasticReactions.resize(reactions.size());
134 mStochasticReactions = NULL;
136 mDeterministicReactions.resize(reactions.size());
137 mDeterministicReactions = NULL;
139 mNumLowSpecies.resize(reactions.size());
142 mStochasticSpecies.resize(mNumReactionSpecies);
143 mStochasticSpecies =
false;
145 mLowSpecies.resize(mNumReactionSpecies);
148 mHasStochastic =
false;
149 mHasDeterministic =
false;
153 const C_FLOAT64 * pValueEnd = pValue + mNumReactionSpecies;
154 bool * pLowSpecies = mLowSpecies.array();
156 C_FLOAT64 CriticalValue = (mLowerThreshold + mUpperThreshold) * 0.5;
160 for (; pValue != pValueEnd; ++pValue, ++Index, ++pLowSpecies)
162 if (*pValue < CriticalValue)
165 std::pair< speciesToReactionsMap::const_iterator, speciesToReactionsMap::const_iterator > Range
166 = mSpeciesToReactions.equal_range(Index);
168 for (; Range.first != Range.second; ++Range.first)
170 mNumLowSpecies[Range.first->second]++;
177 const size_t * pLow = mNumLowSpecies.array();
178 const size_t * pLowEnd = pLow + mNumLowSpecies.size();
179 std::vector< CReactionDependencies >::iterator itReaction = reactions.begin();
183 for (; pLow != pLowEnd; ++pLow, ++itReaction, ++ppStochastic, ++ppDeterministic)
187 *ppStochastic = &(*itReaction);
188 mHasStochastic =
true;
192 *ppDeterministic = &(*itReaction);
193 mHasDeterministic =
true;
197 determineStochasticSpecies();
204 const C_FLOAT64 * pValueEnd = pValue + mNumReactionSpecies;
205 bool * pLowSpecies = mLowSpecies.array();
207 bool PartitionChanged =
false;
211 for (; pValue != pValueEnd; ++pValue, ++Index, ++pLowSpecies)
213 if (!*pLowSpecies && *pValue < mLowerThreshold)
216 PartitionChanged =
true;
218 mStochasticSpecies[Index] =
true;
220 std::pair< speciesToReactionsMap::const_iterator, speciesToReactionsMap::const_iterator > Range
223 for (; Range.first != Range.second; ++Range.first)
225 mNumLowSpecies[Range.first->second]++;
228 else if (*pLowSpecies && *pValue > mUpperThreshold)
230 *pLowSpecies =
false;
231 PartitionChanged =
true;
233 mStochasticSpecies[Index] =
true;
235 std::pair< speciesToReactionsMap::const_iterator, speciesToReactionsMap::const_iterator > Range
236 = mSpeciesToReactions.equal_range(Index);
238 for (; Range.first != Range.second; ++Range.first)
240 mNumLowSpecies[Range.first->second]--;
245 if (!PartitionChanged)
250 PartitionChanged =
false;
254 const size_t * pLow = mNumLowSpecies.array();
255 const size_t * pLowEnd = pLow + mNumLowSpecies.size();
258 mHasStochastic =
false;
259 mHasDeterministic =
false;
261 for (; pLow != pLowEnd; ++pLow, ++ppStochastic, ++ppDeterministic)
265 mHasStochastic =
true;
267 if (*ppDeterministic != NULL)
269 PartitionChanged =
true;
270 *ppStochastic = *ppDeterministic;
271 *ppDeterministic = NULL;
276 mHasDeterministic =
true;
278 if (*ppStochastic != NULL)
280 PartitionChanged =
true;
281 *ppDeterministic = *ppStochastic;
282 *ppStochastic = NULL;
287 if (PartitionChanged)
292 return PartitionChanged;
301 mStochasticSpecies =
false;
303 for (; ppReaction != ppReactionEnd; ++ppReaction)
305 if (*ppReaction != NULL)
308 size_t * pSpeciesIndexEnd = pSpeciesIndex + (*ppReaction)->mSpeciesIndex.size();
310 for (; pSpeciesIndex != pSpeciesIndexEnd; ++pSpeciesIndex)
367 if ((pParm =
getParameter(
"HYBRID.LowerStochLimit")) != NULL)
373 if ((pParm =
getParameter(
"HYBRID.UpperStochLimit")) != NULL)
379 if ((pParm =
getParameter(
"HYBRID.PartitioningInterval")) != NULL)
424 while (fabs(Time - EndTime) > Tolerance)
444 bool FireReaction =
false;
466 for (; (sum <= rand) && (pAmu != endAmu); ++pAmu, ++
mNextReactionIndex, ++ppStochastic)
468 if (*ppStochastic != NULL)
474 assert(mNextReactionIndex > 0);
476 mNextReactionIndex--;
563 size_t IndexSpecies = 1;
565 for (; ppEntity != endEntity; ++ppEntity, ++pValue, ++IndexSpecies)
567 if (dynamic_cast< const CMetab * >(*ppEntity) != NULL)
569 *pValue = floor(*pValue + 0.5);
588 size_t IndexReactions = 0;
589 std::multimap< size_t, size_t> SpeciesToReactions;
595 for (; it != end; ++it)
601 if (Balances.
size() == 0 && Substrates.
size() == 0)
606 itDependencies->mpParticleFlux = (
C_FLOAT64 *)(*it)->getParticleFluxReference()->getValuePointer();
608 itDependencies->mSpeciesMultiplier.resize(Balances.
size());
609 itDependencies->mMethodSpecies.resize(Balances.
size());
610 itDependencies->mModelSpecies.resize(Balances.
size());
612 std::set< size_t > SpeciesIndexSet;
614 std::set< const CCopasiObject * > changed;
624 for (; itBalance != endBalance; ++itBalance)
626 const CMetab * pMetab = (*itBalance)->getMetabolite();
630 size_t SpeciesIndex = StateTemplate.
getIndex(pMetab);
632 itDependencies->mSpeciesMultiplier[Index] = floor((*itBalance)->getMultiplicity() + 0.5);
633 itDependencies->mMethodSpecies[Index] = pMethodStateValue + SpeciesIndex;
638 SpeciesToReactions.insert(std::make_pair(SpeciesIndex, IndexReactions));
639 SpeciesIndexSet.insert(SpeciesIndex);
646 itDependencies->mSpeciesMultiplier.resize(Index,
true);
647 itDependencies->mMethodSpecies.resize(Index,
true);
648 itDependencies->mModelSpecies.resize(Index,
true);
650 itDependencies->mSubstrateMultiplier.resize(Substrates.
size());
651 itDependencies->mMethodSubstrates.resize(Substrates.
size());
652 itDependencies->mModelSubstrates.resize(Substrates.
size());
659 for (; itSubstrate != endSubstrate; ++itSubstrate, ++Index)
661 const CMetab * pMetab = (*itSubstrate)->getMetabolite();
663 size_t SpeciesIndex = StateTemplate.
getIndex(pMetab);
665 itDependencies->mSubstrateMultiplier[Index] = floor((*itSubstrate)->getMultiplicity() + 0.5);
666 itDependencies->mMethodSubstrates[Index] = pMethodStateValue + SpeciesIndex;
671 SpeciesIndexSet.insert(SpeciesIndex);
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();
680 for (; itSet != endSet; ++itSet, ++pSpeciesIndex)
682 *pSpeciesIndex = *itSet;
685 std::set< const CCopasiObject * > dependend;
693 for (; itReaction != end; ++itReaction, ++Index)
695 if ((*itReaction)->getParticleFluxReference()->dependsOn(changed))
697 dependend.insert((*itReaction)->getParticleFluxReference());
698 itDependencies->mDependentReactions[Count] = Index;
704 itDependencies->mDependentReactions.resize(Count,
true);
753 for (; pStochastic != pStochasticEnd; ++pStochastic, ++pYdot)
813 for (; pMultiplier != endMultiplier; ++pMultiplier, ++ppLocalSpecies, ++ppModelSpecies)
815 **ppLocalSpecies += *pMultiplier;
816 **ppModelSpecies = **ppLocalSpecies;
820 std::vector< Refresh * >::const_iterator itCalcualtion = Dependencies.
mCalculations.begin();
821 std::vector< Refresh * >::const_iterator endCalcualtion = Dependencies.
mCalculations.end();
823 while (itCalcualtion != endCalcualtion)
825 (**itCalcualtion++)();
836 for (; pDependentReaction != endDependentReactions; ++pDependentReaction)
853 for (; pAmu != endAmu; ++pAmu, ++ppStochastic)
855 if (*ppStochastic != NULL)
893 bool ApplyCorrection =
false;
900 for (; pMultiplier != endMultiplier; ++pMultiplier, ++ppLocalSubstrate, ++ppModelSubstrate)
902 Multiplicity = *pMultiplier;
905 **ppLocalSubstrate = floor(**ppModelSubstrate + 0.5);
907 if (Multiplicity > 1.01)
909 ApplyCorrection =
true;
911 Number = **ppLocalSubstrate;
913 LowerBound = Number - Multiplicity;
914 SubstrateDevisor *= pow(Number, Multiplicity - 1.0);
917 while (Number > LowerBound)
919 SubstrateMultiplier *= Number;
926 if (SubstrateMultiplier < 0.5 || SubstrateDevisor < 0.5)
938 Amu *= SubstrateMultiplier / SubstrateDevisor;
949 for (
size_t Index = 0; Index != this->
mNumReactions; ++Index, ++ppStochastic)
951 if (*ppStochastic != NULL)
969 for (; pAmu != endAmu; ++pAmu, ++ppStochastic)
971 if (*ppStochastic != NULL)
1000 for (i = 0; i < imax; ++i)
1020 for (i = 0; i < imax; ++i)
1032 for (i = 0; i < imax; ++i)
CTrajectoryMethodDsaLsodar(const CCopasiMethod::SubType &subType=DsaLsodar, const CCopasiContainer *pParent=NULL)
void calculateAmu(const size_t &index)
virtual bool elevateChildren()
CVector< size_t > mDependentReactions
CVector< size_t > mSpeciesIndex
CCopasiVectorN< CEvent > & getEvents()
virtual void evalR(const C_FLOAT64 *t, const C_FLOAT64 *y, const C_INT *nr, C_FLOAT64 *r)
const CCopasiVector< CMetab > & getMetabolites() const
const CCopasiVectorN< CModelValue > & getModelValues() const
virtual C_FLOAT64 getRandomOO()
C_FLOAT64 mNextReactionTime
void initializeParameter()
virtual size_t size() const
CVector< C_FLOAT64 > mSpeciesMultiplier
size_t getNumMetabs() const
void updateSimulatedValues(const bool &updateMoieties)
virtual void start(const CState *initialState)
virtual Status step(const double &deltaT)
CVector< CReactionDependencies * > mStochasticReactions
CVector< C_FLOAT64 * > mMethodSpecies
C_FLOAT64 * mpPartitioningSteps
virtual void stateChanged()
virtual void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
void resize(size_t size, const bool ©=false)
size_t getIndex(const CModelEntity *entity) const
virtual bool isValidProblem(const CCopasiProblem *pProblem)
unsigned C_INT32 * mpMaxSteps
virtual C_FLOAT64 doSingleStep(C_FLOAT64 currentTime, C_FLOAT64 endTime)
void setTime(const C_FLOAT64 &time)
virtual Status step(const double &deltaT)
size_t mNextReactionIndex
CReactionDependencies & operator=(const CReactionDependencies &rhs)
CVector< C_FLOAT64 > mAmu
virtual bool isValidProblem(const CCopasiProblem *pProblem)
bool removeParameter(const std::string &name)
void calculatePropensities()
#define MCTrajectoryMethod
std::vector< CType * >::const_iterator const_iterator
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
C_FLOAT64 * mpParticleFlux
size_t mFirstReactionSpeciesIndex
std::vector< CReactionDependencies > mReactionDependencies
virtual bool elevateChildren()
CModelEntity ** endFixed()
void setState(const CState &state)
const C_FLOAT64 & getDuration() const
void integrateDeterministicPart(const C_FLOAT64 &deltaT)
const Value & getValue() const
bool setValue(const std::string &name, const CType &value)
void intialize(std::vector< CReactionDependencies > &reactions, const speciesToReactionsMap &speciesToReactions, const C_FLOAT64 &lowerThreshold, const C_FLOAT64 &upperThreshold, const CState &state)
CModelEntity ** beginIndependent()
CCopasiParameter * getParameter(const std::string &name)
CVector< C_FLOAT64 > mSubstrateMultiplier
void setTime(const C_FLOAT64 &time)
void calculateDerivativesX(C_FLOAT64 *derivativesX)
CCopasiVectorNS< CCompartment > & getCompartments()
CVector< C_FLOAT64 * > mModelSpecies
const C_FLOAT64 & getTime() const
bool rePartition(const CState &state)
CVector< C_FLOAT64 * > mMethodSubstrates
CVector< bool > mStochasticSpecies
const CStateTemplate & getStateTemplate() const
size_t getNumModelValues() const
std::multimap< size_t, size_t > speciesToReactionsMap
const ModelType & getModelType() const
CCopasiVectorNS< CReaction > & getReactions()
void determineStochasticSpecies()
static std::vector< Refresh * > buildUpdateSequence(const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
virtual void stateChanged()
void calculateDerivatives(C_FLOAT64 *derivatives)
virtual void * getValuePointer() const
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
unsigned C_INT32 * mpPartitioningInterval
const CModelEntity::Status & getStatus() const
CVector< C_FLOAT64 * > mModelSubstrates
void fireReaction(const size_t &index)
std::vector< Refresh * > mCalculations
size_t mStepsAfterPartitionSystem
~CTrajectoryMethodDsaLsodar()
std::string suitableForStochasticSimulation() const
CModel * getModel() const
CRandom * mpRandomGenerator
virtual void start(const CState *initialState)
C_FLOAT64 * beginIndependent()
CCopasiObject * getValueReference() const