COPASI API  4.16.103
CProcessQueue.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2014 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) 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 #include <limits>
12 
13 #include "copasi.h"
14 
15 #include "CProcessQueue.h"
16 #include "CMathModel.h"
17 
18 #include "function/CExpression.h"
19 
21  mExecutionTime(0.0),
22  mCascadingLevel(0),
23  mEquality(false)
24 {}
25 
27  mExecutionTime(src.mExecutionTime),
29  mEquality(src.mEquality)
30 {}
31 
32 CProcessQueue::CKey::CKey(const C_FLOAT64 & executionTime,
33  const bool & equality,
34  const size_t & cascadingLevel) :
35  mExecutionTime(executionTime),
36  mCascadingLevel(cascadingLevel),
37  mEquality(equality)
38 {}
39 
41 {}
42 
44 {
45  if (mExecutionTime != rhs.mExecutionTime)
46  return mExecutionTime < rhs.mExecutionTime;
47 
49  return mCascadingLevel > rhs.mCascadingLevel;
50 
51  return mEquality > rhs.mEquality;
52 }
53 
54 //*********************************************************
55 
57  mType(),
58  mValues(),
59  mpEvent(NULL),
60  mpProcessQueue(NULL)
61 {}
62 
64  mType(src.mType),
65  mValues(src.mValues),
66  mpEvent(src.mpEvent),
67  mpProcessQueue(src.mpProcessQueue)
68 {}
69 
71  CProcessQueue * pProcessQueue) :
72  mType(CProcessQueue::CAction::Calculation),
73  mValues(),
74  mpEvent(pEvent),
75  mpProcessQueue(pProcessQueue)
76 {}
77 
79  CMathEvent * pEvent,
80  CProcessQueue * pProcessQueue) :
81  mType(),
82  mValues(values),
83  mpEvent(pEvent),
84  mpProcessQueue(pProcessQueue)
85 {
86  switch (pEvent->getType())
87  {
88  case CEvent::Assignment:
89  mType = Assignment;
90  break;
91 
92  case CEvent::CutPlane:
93  mType = Callback;
94  break;
95  }
96 }
97 
99 {}
100 
102 {
103  // Assume that nothing is changed
104  bool StateChanged = false;
105 
106  switch (mType)
107  {
108  case Calculation:
109 
110  if (mpEvent->delayAssignment())
111  {
112  mpProcessQueue->addAssignment(mpEvent->getAssignmentTime(mpProcessQueue->mTime),
113  mpProcessQueue->mEquality,
114  mpEvent->getTargetValues(),
115  mpEvent);
116  }
117  else
118  {
119  StateChanged = mpEvent->executeAssignment();
120  }
121 
122  break;
123 
124  case Assignment:
125  StateChanged = mpEvent->setTargetValues(mValues);
126  break;
127 
128  case Callback:
129 
130  if (mpProcessQueue->mpEventCallBack)
131  {
132  (*mpProcessQueue->mpEventCallBack)(mpProcessQueue->mpCallbackTask, mpEvent->getType());
133  }
134 
135  break;
136  }
137 
138  return StateChanged;
139 }
140 
142 {
143  return mType;
144 }
145 
146 //*********************************************************
147 
149  mActions(),
150  mExecutionLimit(10000),
152  mTime(0.0),
153  mEquality(true),
154  mCascadingLevel(0),
156  mEventIdSet(),
157  mpMathModel(NULL),
158  mRootsFound(0),
159  mRootValues1(0),
160  mRootValues2(0),
165  mpCallbackTask(NULL),
166  mpEventCallBack(NULL)
167 {}
168 
170  mActions(src.mActions),
171  mExecutionLimit(src.mExecutionLimit),
172  mExecutionCounter(src.mExecutionCounter),
173  mTime(src.mTime),
174  mEquality(src.mEquality),
175  mCascadingLevel(src.mCascadingLevel),
176  mSimultaneousAssignmentsFound(src.mSimultaneousAssignmentsFound),
177  mEventIdSet(src.mEventIdSet),
178  mpMathModel(src.mpMathModel),
179  mRootsFound(src.mRootsFound),
180  mRootValues1(src.mRootValues1),
181  mRootValues2(src.mRootValues2),
182  mpRootValuesBefore(&src.mRootValues1 == src.mpRootValuesBefore ? &mRootValues1 : &mRootValues2),
183  mpRootValuesAfter(&src.mRootValues1 == src.mpRootValuesAfter ? &mRootValues1 : &mRootValues2),
184  mpResolveSimultaneousAssignments(src.mpResolveSimultaneousAssignments),
185  mContinueSimultaneousEvents(src.mContinueSimultaneousEvents),
186  mpCallbackTask(src.mpCallbackTask),
187  mpEventCallBack(src.mpEventCallBack)
188 {}
189 
191 {}
192 
193 bool CProcessQueue::addAssignment(const C_FLOAT64 & executionTime,
194  const bool & equality,
195  const CVector< C_FLOAT64 > & values,
196  CMathEvent * pEvent)
197 {
198  // It is not possible to proceed backwards in time.
199  if (executionTime < mTime) return false;
200 
201  size_t CascadingLevel = mCascadingLevel;
202 
203  if (executionTime > mTime)
204  CascadingLevel = 0;
205 
206  mActions.insert(std::make_pair(CKey(executionTime, equality, CascadingLevel),
207  CAction(values, pEvent, this)));
208 
209  return true;
210 }
211 
212 bool CProcessQueue::addCalculation(const C_FLOAT64 & executionTime,
213  const bool & equality,
214  CMathEvent * pEvent)
215 {
216  // It is not possible to proceed backwards in time.
217  if (executionTime < mTime) return false;
218 
219  size_t CascadingLevel = mCascadingLevel;
220 
221  if (executionTime > mTime)
222  CascadingLevel = 0;
223 
224  mActions.insert(std::make_pair(CKey(executionTime, equality, CascadingLevel),
225  CAction(pEvent, this)));
226 
227  return true;
228 }
229 
231 {
232  mpMathModel = pMathModel;
233  assert(mpMathModel != NULL);
234 
236 
237  mActions.clear();
238  mEventIdSet.clear();
240 
241  size_t NumRoots = mpMathModel->getNumRoots();
242  mRootsFound.resize(NumRoots);
243  mRootsFound = 0;
244  mRootValues1.resize(NumRoots);
245  mRootValues2.resize(NumRoots);
248 
249  return;
250 }
251 
253  const bool & priorToOutput,
254  resolveSimultaneousAssignments pResolveSimultaneousAssignments)
255 {
256  if (getProcessQueueExecutionTime() > time)
257  return false;
258 
259  mTime = time;
260  mEquality = priorToOutput;
261  mpResolveSimultaneousAssignments = pResolveSimultaneousAssignments;
262  mExecutionCounter = 0;
263  mCascadingLevel = 0;
264 
265  bool success = true;
266  bool stateChanged = false;
267 
269 
270  iterator itAction = getAction();
271 
272  // The algorithm below will work properly for user ordered events
273  // as the queue enforces the proper ordering.
274  while (success &&
275  itAction != mActions.end() &&
277  {
278  if (!executeAction(itAction))
279  {
280  itAction = getAction();
281  continue;
282  }
283 
284  stateChanged = true;
285 
286  // We switch to the next cascading level so that events triggered by the
287  // execution of assignments are properly scheduled.
288  mCascadingLevel++;
289 
290  // We need to compare the roots before the execution and after
291  // to determine which roots need to be charged.
292  if (rootsFound())
293  {
294  // We have to deal with both types of found roots.
296  }
297 
298  // Note, applying the events may have added new events to the queue.
299  // The setting of the equality flag for these events may be either true
300  // or false.
301 
302  // Independent from the setting priorToOutput we must only consider equality
303  mEquality = true;
304 
305  // Retrieve the pending calculations.
306  itAction = getAction();
307 
308  while (itAction == mActions.end() &&
309  mCascadingLevel > 0)
310  {
311  // If we are here we have no more actions for this level.
312  mCascadingLevel--;
313 
314  if (mCascadingLevel == 0)
315  {
316  mEquality = priorToOutput;
317  }
318 
319  itAction = getAction();
320  }
321  }
322 
323  return stateChanged;
324 }
325 
327 {
328  CKey Pending(mTime, mEquality, mCascadingLevel);
329  range PendingActions = mActions.equal_range(Pending);
330 
331  if (PendingActions.first == PendingActions.second)
332  {
333  return mActions.end();
334  }
335 
336  iterator nextAction = PendingActions.first;
337  ++nextAction;
338 
339  if (nextAction != PendingActions.second)
340  {
342  {
343 
344  // The resolution of simultaneous events is algorithm dependent.
345  // The simulation routine should provide a call back function.
347  {
349  }
350 
352  }
353  else if (mSimultaneousAssignmentsFound == false)
354  {
356  CCopasiMessage(CCopasiMessage::WARNING_FILTERED, "CMathModel (1): Simultaneous event assignments encountered.");
357  }
358  }
359 
360  return PendingActions.first;
361 }
362 
364 {
365  bool StateChanged = itAction->second.process();
366 
367  if (StateChanged)
368  {
370  }
371 
372  mActions.erase(itAction);
373 
374  return StateChanged;
375 }
376 
378 {
379  bool rootsFound = false;
380 
381  // Calculate the current root values
383 
384  // Compare the root values before and after;
385  C_INT * pRootFound = mRootsFound.array();
386  C_INT * pRootEnd = pRootFound + mRootsFound.size();
387  C_FLOAT64 * pValueBefore = mpRootValuesBefore->array();
388  C_FLOAT64 * pValueAfter = mpRootValuesAfter->array();
389  CMathTrigger::CRootFinder *const* ppRootFinder = mpMathModel->getRootFinders().array();
390 
391  for (; pRootFound != pRootEnd; ++pRootFound, ++pValueBefore, ++pValueAfter, ++ppRootFinder)
392  {
393  // Root values which did not change are not found
394  if (*pValueBefore == *pValueAfter)
395  {
396  *pRootFound = 0;
397  continue;
398  }
399 
400  // Handle equality
401  if ((*ppRootFinder)->isEquality())
402  {
403  if ((*ppRootFinder)->isTrue())
404  {
405  if (*pValueAfter >= 0.0 || *pValueAfter > *pValueBefore)
406  {
407  *pRootFound = 0;
408  }
409  else
410  {
411  *pRootFound = 1;
412  rootsFound = true;
413  }
414  }
415  else
416  {
417  if (*pValueAfter < 0.0 || *pValueAfter < *pValueBefore)
418  {
419  *pRootFound = 0;
420  }
421  else
422  {
423  *pRootFound = 1;
424  rootsFound = true;
425  }
426  }
427  }
428  else
429  {
430  if ((*ppRootFinder)->isTrue())
431  {
432  if (*pValueAfter > 0.0 || *pValueAfter > *pValueBefore)
433  {
434  *pRootFound = 0;
435  }
436  else
437  {
438  *pRootFound = 1;
439  rootsFound = true;
440  }
441  }
442  else
443  {
444  if (*pValueAfter <= 0.0 || *pValueAfter < *pValueBefore)
445  {
446  *pRootFound = 0;
447  }
448  else
449  {
450  *pRootFound = 1;
451  rootsFound = true;
452  }
453  }
454  }
455  }
456 
457  // Swap before and after.
460  mpRootValuesAfter = pTmp;
461 
462  return rootsFound;
463 }
464 
465 // static
467 {
468  return range.first != range.second;
469 }
470 
471 void CProcessQueue::destroyEventId(const size_t & eventId)
472 {
473  mEventIdSet.erase(eventId);
474 }
475 
477 {
478  static C_FLOAT64 Infinity =
479  std::numeric_limits< C_FLOAT64 >::infinity();
480 
481  if (mActions.empty())
482  {
483  return Infinity;
484  }
485 
486  return mActions.begin()->first.getExecutionTime();
487 }
488 
489 void CProcessQueue::setContinueSimultaneousEvents(const bool & continueSimultaneousEvents)
490 {
491  mContinueSimultaneousEvents = continueSimultaneousEvents;
492 }
493 
495 {
497 }
498 
499 void CProcessQueue::setEventCallBack(void* pTask, EventCallBack ecb)
500 {
502  mpEventCallBack = ecb;
503 }
504 
505 std::ostream &operator<<(std::ostream &os, const CProcessQueue & o)
506 {
507  os << "Process Queue" << std::endl;
508  std::multimap< CProcessQueue::CKey, CProcessQueue::CAction >::const_iterator it;
509 
510  if (o.mActions.size()) os << " Actions:" << std::endl;
511 
512  for (it = o.mActions.begin(); it != o.mActions.end(); ++it)
513  {
514  os << "exec time " << it->first.mExecutionTime
515  << ", cascading lvl " << it->first.mCascadingLevel
516  << ", " << (it->first.mEquality ? "equality, " : "inequality");
517 
518  os << std::endl;
519 
520  CMathEvent * pEvent = it->second.mpEvent;
521 
522  os << "pEvent: 0x" << pEvent << ", Action: ";
523 
524  switch (it->second.mType)
525  {
527 
528  if (pEvent->delayAssignment())
529  {
530  os << "Calculation";
531  }
532  else
533  {
534  os << "Calculation & Assignment";
535  }
536 
537  break;
538 
540  os << "Assignment";
541  break;
542 
544  os << "Callback";
545  break;
546  }
547 
548  os << std::endl << std::endl;
549  }
550 
551  return os;
552 }
Header file of class CExpression.
void initialize(CMathModel *pMathModel)
C_FLOAT64 mTime
#define C_INT
Definition: copasi.h:115
iterator getAction()
std::ostream & operator<<(std::ostream &os, const CProcessQueue &o)
CVector< C_FLOAT64 > * mpRootValuesAfter
const C_FLOAT64 & getProcessQueueExecutionTime() const
C_FLOAT64 mExecutionTime
Definition: CProcessQueue.h:79
CVector< C_FLOAT64 > mRootValues1
static bool notEmpty(const range &range)
void evaluateRoots(CVectorCore< C_FLOAT64 > &rootValues, const bool &ignoreDiscrete)
Definition: CMathModel.cpp:169
const Type & getType() const
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
bool addAssignment(const C_FLOAT64 &executionTime, const bool &equality, const CVector< C_FLOAT64 > &values, CMathEvent *pEvent)
size_t mExecutionLimit
std::pair< std::multimap< CKey, CAction >::iterator, std::multimap< CKey, CAction >::iterator > range
std::set< size_t > mEventIdSet
const C_FLOAT64 & getInitialTime() const
Definition: CMathModel.cpp:206
const CEvent::Type & getType() const
size_t mCascadingLevel
#define MCMathModel
EventCallBack mpEventCallBack
CTSSATask * pTask
bool operator<(const CKey &rhs) const
void setContinueSimultaneousEvents(const bool &continueSimultaneousEvents)
void setEventCallBack(void *pTask, EventCallBack ecb)
CMathModel * mpMathModel
CVector< C_INT > mRootsFound
bool executeAction(CProcessQueue::iterator itAction)
bool process(const C_FLOAT64 &time, const bool &priorToOutput, resolveSimultaneousAssignments pResolveSimultaneousAssignments)
CVector< C_FLOAT64 > mRootValues2
CVector< C_FLOAT64 > * mpRootValuesBefore
void destroyEventId(const size_t &eventId)
const bool & getContinueSimultaneousEvents() const
size_t size() const
Definition: CVector.h:100
bool mContinueSimultaneousEvents
#define C_FLOAT64
Definition: copasi.h:92
void * mpCallbackTask
size_t mExecutionCounter
CType * array()
Definition: CVector.h:139
bool addCalculation(const C_FLOAT64 &executionTime, const bool &equality, CMathEvent *pEvent)
bool mSimultaneousAssignmentsFound
resolveSimultaneousAssignments mpResolveSimultaneousAssignments
size_t getNumRoots() const
Definition: CMathModel.cpp:412
std::multimap< CKey, CAction > mActions
const CVector< CMathTrigger::CRootFinder * > & getRootFinders() const
Definition: CMathModel.cpp:447
void processRoots(const C_FLOAT64 &time, const bool &equality, const bool &correct, const CVector< C_INT > &roots)
Definition: CMathModel.cpp:218
std::multimap< CKey, CAction >::iterator iterator
#define max(a, b)
Definition: f2c.h:176