25 mEvents(
"ListOfMathEvents", this),
30 mRootIndex2RootFinder()
37 mProcessQueue(src.mProcessQueue),
38 mEvents(
"ListOfMathEvents", this),
43 mRootIndex2RootFinder()
83 std::vector< CCopasiContainer * > Context;
89 size_t RootFinderCount = 0;
91 for (; itEvent != endEvent; ++itEvent)
95 success &= pEvent->
compile(*itEvent, Context);
119 std::set< const CCopasiObject * > StateVariables;
123 for (; ppEntity != ppEntityEnd; ++ppEntity)
125 StateVariables.insert((*ppEntity)->getValueReference());
128 std::set< const CCopasiObject * > RootValuesDependencies;
134 for (; itMathEvent != endMathEvent; ++itMathEvent)
137 (*itMathEvent)->getMathTrigger().getRootFinders().
begin();
139 (*itMathEvent)->getMathTrigger().getRootFinders().
end();
142 for (; itRootFinder != endRootFinder;
143 ++itRootFinder, ++ppRootValue, ++pRootDiscrete, ++ppEvent, ++ppRootFinder)
146 *ppRootValue = (*itRootFinder)->getRootValuePtr();
149 (*itRootFinder)->determineDiscrete(StateVariables);
150 *pRootDiscrete = (*itRootFinder)->isDiscrete();
153 *ppEvent = *itMathEvent;
156 *ppRootFinder = *itRootFinder;
159 RootValuesDependencies.insert(*itRootFinder);
170 const bool & ignoreDiscrete)
173 std::vector< Refresh * >::const_iterator itRefresh =
mRootRefreshes.begin();
174 std::vector< Refresh * >::const_iterator endRefresh =
mRootRefreshes.end();
176 while (itRefresh != endRefresh)
190 for (; pSrc != pSrcEnd; ++pSrc, ++pDiscrete, ++pTarget)
192 *pTarget = (*pDiscrete) ? 1.0 : **pSrc;
197 for (; pSrc != pSrcEnd; ++pSrc, ++pTarget)
212 const bool & equality,
219 const bool & equality,
220 const bool & correct,
224 std::vector< Refresh * >::const_iterator itRefresh =
mRootRefreshes.begin();
225 std::vector< Refresh * >::const_iterator endRefresh =
mRootRefreshes.end();
227 while (itRefresh != endRefresh)
236 const C_INT *pFoundRootEnd = pFoundRoot + foundRoots.
size();
245 while (pFoundRoot != pFoundRootEnd)
248 if (*pFoundRoot < 1 && (*ppRootFinder)->isEquality() == equality)
253 ++pFoundRoot; ++ppRootFinder;
257 pFoundRoot = foundRoots.
array();
262 bool TriggerBefore =
true;
264 while (pFoundRoot != pFoundRootEnd)
266 pProcessEvent = *ppEvent;
271 while (pFoundRoot != pFoundRootEnd &&
272 *ppEvent == pProcessEvent)
277 (*ppRootFinder)->toggle(time, equality, correct);
280 ++pFoundRoot; ++ppEvent; ++ppRootFinder;
286 if (TriggerAfter ==
true &&
287 TriggerBefore ==
false)
300 std::vector< Refresh * >::const_iterator itRefresh =
mRootRefreshes.begin();
301 std::vector< Refresh * >::const_iterator endRefresh =
mRootRefreshes.end();
303 while (itRefresh != endRefresh)
312 const C_INT *pFoundRootEnd = pFoundRoot + foundRoots.
size();
321 bool TriggerBefore =
true;
323 while (pFoundRoot != pFoundRootEnd)
325 pProcessEvent = *ppEvent;
330 while (pFoundRoot != pFoundRootEnd &&
331 *ppEvent == pProcessEvent)
336 (*ppRootFinder)->toggle(time);
339 ++pFoundRoot; ++ppEvent; ++ppRootFinder;
345 if (TriggerAfter ==
true &&
346 TriggerBefore ==
false)
379 std::vector< Refresh * >::const_iterator itRefresh =
mRootRefreshes.begin();
380 std::vector< Refresh * >::const_iterator endRefresh =
mRootRefreshes.end();
382 while (itRefresh != endRefresh)
390 for (; itMathEvent != endMathEvent; ++itMathEvent)
392 (*itMathEvent)->getMathTrigger().applyInitialValues();
443 dgemm_(&T, &T, &M, &N, &K, &Alpha, Rates.
array(), &M,
444 Jacobian.
array(), &
K, &Beta, pDerivative, &M);
456 std::set< const CCopasiObject * > UpToDate;
464 return std::vector< Refresh * >();
478 std::set< const CCopasiObject * > RequiredObjects;
494 for (; ppEntity != ppEndEntity; ++ppEntity)
496 switch ((*ppEntity)->getStatus())
502 if ((*ppEntity)->getRateReference()->dependsOn(changedObjects, changedObjects))
504 RequiredObjects.insert((*ppEntity)->getRateReference());
510 pSpecies =
dynamic_cast< const CMetab *
>(*ppEntity);
512 if (pSpecies != NULL &&
538 for (; itMoiety != endMoiety; ++itMoiety)
540 if ((*itMoiety)->getValueReference()->dependsOn(changedObjects, changedObjects))
542 RequiredObjects.insert((*itMoiety)->getValueReference());
546 std::set< const CCopasiObject * > UpToDate;
564 for (; ppRootFinder != ppEndRootFinder;
565 ++ppRootFinder, ++pRootValue, ++pDerivative, ++pFoundRoot)
567 if (**pRootValue == 0.0)
569 if (!(*ppRootFinder)->isEquality() &&
578 if ((*ppRootFinder)->isEquality() &&
581 (*ppRootFinder)->toggle(std::numeric_limits< C_FLOAT64 >::quiet_NaN(),
true,
false);
604 jacobian.
resize(NumRows, NumCols);
621 for (; pX != pXEnd; ++pX, ++Col, ++pRate)
625 if (fabs(*pRate) < 1e4 * std::numeric_limits< C_FLOAT64 >::epsilon() * fabs(Store) ||
643 InvDelta = 500.0 / Store;
648 X1 = Store - 0.001 * *pRate;
649 X2 = Store + 0.001 * *pRate;
650 InvDelta = 500.0 / *pRate;
663 pJacobian = jacobian.
array() + Col;
667 for (; pJacobian < pJacobianEnd; pJacobian += NumCols, ++pY1, ++pY2)
668 * pJacobian = (*pY2 - *pY1) * InvDelta;
void initialize(CMathModel *pMathModel)
const C_FLOAT64 & getProcessQueueExecutionTime() const
CCopasiVectorN< CEvent > & getEvents()
bool compile(const CEvent *pEvent, std::vector< CCopasiContainer * > listOfContainer)
size_t getNumIndependent() const
CVector< C_FLOAT64 * > mRootValues
const C_FLOAT64 & getProcessQueueExecutionTime() const
void calculateRootJacobian(CMatrix< C_FLOAT64 > &jacobian, const CVector< C_FLOAT64 > &rates)
void updateSimulatedValues(const bool &updateMoieties)
std::vector< Refresh * > buildRequiredRefreshList(const std::set< const CCopasiObject * > &requiredObjects) const
void fire(const C_FLOAT64 &time, const bool &equality, CProcessQueue &processQueue)
void evaluateRoots(CVectorCore< C_FLOAT64 > &rootValues, const bool &ignoreDiscrete)
void resize(size_t size, const bool ©=false)
void applyInitialValues()
std::vector< Refresh * > buildDependendRefreshList(const std::set< const CCopasiObject * > &changedObjects) const
CModelEntity ** endDependent()
const C_FLOAT64 & getInitialTime() const
const size_t & size() const
CProcessQueue mProcessQueue
const CCopasiVector< CMoiety > & getMoieties() const
CMathTrigger & getMathTrigger()
size_t getNumDependentReactionMetabs() const
void setState(const CState &state)
CMathModel(const CCopasiContainer *pParent=NULL)
virtual bool add(const CType &src)
const C_FLOAT64 & getInitialTime() const
size_t getNumIndependent() const
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
virtual void resize(size_t rows, size_t cols, const bool ©=false)
bool process(const C_FLOAT64 &time, const bool &priorToOutput, resolveSimultaneousAssignments pResolveSimultaneousAssignments)
bool dependsOn(DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const
std::vector< Refresh * > mRootRefreshes
iterator(* resolveSimultaneousAssignments)(const std::multimap< CKey, CAction > &, const C_FLOAT64 &, const bool &, const size_t &)
CVector< bool > mRootDiscrete
const CStateTemplate & getStateTemplate() const
virtual size_t size() const
static std::vector< Refresh * > buildUpdateSequence(const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())
void calculateDerivatives(C_FLOAT64 *derivatives)
void calculateTrueValue()
CVector< CMathTrigger::CRootFinder * > mRootIndex2RootFinder
CCopasiVector< CMathEvent > mEvents
const CState & getState() const
void calculateRootDerivatives(CVector< C_FLOAT64 > &rootDerivatives)
const CProcessQueue & getProcessQueue() const
CVector< CMathEvent * > mRootIndex2Event
size_t getNumRoots() const
bool determineInitialRoots(CVector< C_INT > &foundRoots)
const CVector< CMathTrigger::CRootFinder * > & getRootFinders() const
void processRoots(const C_FLOAT64 &time, const bool &equality, const bool &correct, const CVector< C_INT > &roots)
bool compile(CModel *pModel)
C_FLOAT64 * beginIndependent()
CModelEntity ** getEntities()
CCopasiVector< CRootFinder > & getRootFinders()
bool processQueue(const C_FLOAT64 &time, const bool &equality, CProcessQueue::resolveSimultaneousAssignments pResolveSimultaneousAssignments)
CCopasiObject * getValueReference() const