COPASI API  4.16.103
CCrossSectionTask.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 /**
7  * CCrossSectionTask class.
8  *
9  * This class implements a cross section task.
10  */
11 
12 #include <string>
13 #include <sstream>
14 #include <cmath>
15 
16 #include "copasi.h"
17 #include <model/CEvent.h>
18 
19 #include "CCrossSectionTask.h"
20 #include "CCrossSectionProblem.h"
21 //#include "CCrossSectionMethod.h"
22 #include "model/CModel.h"
23 #include "model/CMathModel.h"
24 #include "model/CState.h"
25 #include "report/CKeyFactory.h"
26 #include "report/CReport.h"
30 
31 const unsigned int CCrossSectionTask::ValidMethods[] =
32 {
35 };
36 
38  CCopasiTask(CCopasiTask::crosssection, pParent),
39  mTimeSeriesRequested(true),
40  mTimeSeries(),
41  mpCrossSectionProblem(NULL),
42  mpTrajectoryMethod(NULL),
43  mUpdateMoieties(false),
44  mpCurrentState(NULL),
45  mpCurrentTime(NULL),
46  mOutputStartTime(0.0),
47  mStartTime(0.0),
48  mNumCrossings(0),
49  mOutputStartNumCrossings(0),
50  mMaxNumCrossings(std::numeric_limits< size_t >::max()),
51  mhProgress(C_INVALID_INDEX),
52  mProgressMax(1),
53  mProgressValue(0),
54  mProgressFactor(1),
55  mpEvent(NULL),
56  mState(CCrossSectionTask::TRANSIENT),
57  mStatesRing(),
58  mStatesRingCounter(0),
59  mPreviousCrossingTime(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
60  mPeriod(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
61  mAveragePeriod(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
62  mLastPeriod(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
63  mPeriodicity(-1),
64  mLastFreq(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
65  mFreq(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
66  mAverageFreq(std::numeric_limits< C_FLOAT64 >::quiet_NaN())
67 {
68  initObjects();
69  mpProblem = new CCrossSectionProblem(this);
71  this->add(mpMethod, true);
72 }
73 
75  const CCopasiContainer * pParent):
76  CCopasiTask(src, pParent),
77  mTimeSeriesRequested(src.mTimeSeriesRequested),
78  mTimeSeries(),
79  mpCrossSectionProblem(NULL),
80  mpTrajectoryMethod(NULL),
81  mUpdateMoieties(false),
82  mpCurrentState(NULL),
83  mpCurrentTime(NULL),
84  mOutputStartTime(0.0),
85  mStartTime(0.0),
86  mNumCrossings(0),
87  mOutputStartNumCrossings(0),
88  mMaxNumCrossings(std::numeric_limits< size_t >::max()),
89  mhProgress(C_INVALID_INDEX),
90  mProgressMax(1),
91  mProgressValue(0),
92  mProgressFactor(1),
93  mpEvent(NULL),
94  mState(CCrossSectionTask::TRANSIENT),
95  mStatesRing(),
96  mStatesRingCounter(0),
97  mPreviousCrossingTime(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
98  mPeriod(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
99  mAveragePeriod(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
100  mLastPeriod(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
101  mPeriodicity(-1),
102  mLastFreq(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
103  mFreq(std::numeric_limits< C_FLOAT64 >::quiet_NaN()),
104  mAverageFreq(std::numeric_limits< C_FLOAT64 >::quiet_NaN())
105 {
106  initObjects();
107  mpProblem =
108  new CCrossSectionProblem(*static_cast< CCrossSectionProblem * >(src.mpProblem), this);
109 
111  * mpMethod = * src.mpMethod;
112 
113  mpMethod->elevateChildren();
114 
115  this->add(mpMethod, true);
116 }
117 
119 {
120  cleanup();
121 }
122 
124 {
126 }
127 
129 {
137 }
138 
139 #define RING_SIZE 16
140 
142  COutputHandler * pOutputHandler,
143  std::ostream * pOstream)
144 {
145  assert(mpProblem && mpMethod);
146 
148  assert(mpCrossSectionProblem);
149 
150  mpTrajectoryMethod = dynamic_cast<CTrajectoryMethod *>(mpMethod);
151  assert(mpTrajectoryMethod);
152 
154 
155  bool success = mpMethod->isValidProblem(mpProblem);
156 
157  CCopasiParameter * pParameter = mpMethod->getParameter("Integrate Reduced Model");
158 
159  if (pParameter != NULL)
160  mUpdateMoieties = *pParameter->getValue().pBOOL;
161  else
162  mUpdateMoieties = false;
163 
167 
168  //init the ring buffer for the states
169  mStatesRing.resize(RING_SIZE);
170  mStatesRingCounter = 0;
171 
172  // Handle the time series as a regular output.
173  mTimeSeriesRequested = true;//mpCrossSectionProblem->timeSeriesRequested();
174 
175  if ((pOutputHandler != NULL) &&
178  {
179  mTimeSeries.allocate(20);
180  pOutputHandler->addInterface(&mTimeSeries);
181  }
182  else
183  {
184  mTimeSeries.clear();
185  }
186 
187  createEvent();
188 
189  if (!CCopasiTask::initialize(of, pOutputHandler, pOstream)) success = false;
190 
191  return success;
192 }
193 
195 {
196  if (mpEvent != NULL) return;
197 
199 
201  {
202  int count = 0;
203  std::string name = "__cutplane";
204 
205  while (pModel->getEvents().getIndex(name) != C_INVALID_INDEX)
206  {
207  std::stringstream str;
208  str << "__cutplane" << ++count;
209  name = str.str();
210  }
211 
212  mpEvent = pModel->createEvent(name);
214 
215  std::stringstream expression;
216  expression << "<" << mpCrossSectionProblem->getSingleObjectCN() << "> "
217  << (mpCrossSectionProblem->isPositiveDirection() ? std::string(" > ") : std::string(" < "))
219 
220  mpEvent ->setTriggerExpression(expression.str());
221 
222  pModel->compileIfNecessary(NULL);
223  }
224 }
225 
227 {
228  // TODO: remove event
229  if (mpEvent != NULL)
230  {
232  mpEvent = NULL;
233  }
234 }
235 
236 bool CCrossSectionTask::process(const bool & useInitialValues)
237 {
238  processStart(useInitialValues);
239 
240  //this instructs the process queue to call back whenever an event is
241  //executed
243 
244  mPreviousCrossingTime = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
245  mPeriod = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
246  mAveragePeriod = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
247  mLastPeriod = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
248  mPeriodicity = -1;
249  mLastFreq = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
250  mFreq = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
251  mAverageFreq = std::numeric_limits< C_FLOAT64 >::quiet_NaN();
252 
254 
255  //the output starts only after "outputStartTime" has passed
257  {
259  MaxDuration += mpCrossSectionProblem->getOutputStartTime();
260  }
261  else
262  {
264  }
265 
266  const C_FLOAT64 EndTime = *mpCurrentTime + MaxDuration;
267 
269 
270  // It suffices to reach the end time within machine precision
271  C_FLOAT64 CompareEndTime = mOutputStartTime - 100.0 * (fabs(EndTime) * std::numeric_limits< C_FLOAT64 >::epsilon() + std::numeric_limits< C_FLOAT64 >::min());
272 
275  else
276  mMaxNumCrossings = 0;
277 
280  else
282 
284 
285  bool flagProceed = true;
286  mProgressFactor = 100.0 / (MaxDuration + mpCrossSectionProblem->getOutputStartTime());
287  mProgressValue = 0;
288 
289  if (mpCallBack != NULL)
290  {
291  mpCallBack->setName("performing simulation...");
292  mProgressMax = 100.0;
293  mhProgress = mpCallBack->addItem("Completion",
295  &mProgressMax);
296  }
297 
298  mState = TRANSIENT;
299  mStatesRingCounter = 0;
300 
301  mNumCrossings = 0;
302 
303  try
304  {
305  do
306  {
307  flagProceed &= processStep(EndTime);
308  }
309  while ((*mpCurrentTime < CompareEndTime) && flagProceed);
310  }
311 
312  catch (int)
313  {
316  finish();
318  }
319 
320  catch (CCopasiException & Exception)
321  {
324  finish();
325  throw CCopasiException(Exception.getMessage());
326  }
327 
328  finish();
329 
330  return true;
331 }
332 
333 void CCrossSectionTask::processStart(const bool & useInitialValues)
334 {
335  if (useInitialValues)
337 
339 
342 
343  return;
344 }
345 
347 {
348  CModel * pModel = mpCrossSectionProblem->getModel();
349  bool StateChanged = false;
350  bool proceed = true;
351 
352  C_FLOAT64 Tolerance = 100.0 * (fabs(endTime) * std::numeric_limits< C_FLOAT64 >::epsilon() + std::numeric_limits< C_FLOAT64 >::min());
353  C_FLOAT64 NextTime = endTime;
354 
355  while (proceed)
356  {
357  // TODO Provide a call back method for resolving simultaneous assignments.
358  //pModel->getMathModel()->getProcessQueue().printDebug();
359 
360  //execute events for inequalities
361  StateChanged |= pModel->processQueue(*mpCurrentTime, false, NULL);
362 
363  if (StateChanged)
364  {
365  *mpCurrentState = pModel->getState();
367  StateChanged = false;
368  }
369 
370  // std::min suffices since events are only supported in forward integration.
371  NextTime = std::min(endTime, pModel->getProcessQueueExecutionTime());
372 
373  switch (mpTrajectoryMethod->step(NextTime - *mpCurrentTime))
374  {
376  pModel->setState(*mpCurrentState);
378 
379  // TODO Provide a call back method for resolving simultaneous assignments.
380 
381  //execute events for equalities
382  StateChanged |= pModel->processQueue(*mpCurrentTime, true, NULL);
383 
384  // If the state change happens to coincide with end of the step we have to return and
385  // inform the integrator of eventual state changes.
386  if (fabs(*mpCurrentTime - endTime) < Tolerance)
387  {
388  if (StateChanged)
389  {
390  *mpCurrentState = pModel->getState();
392  StateChanged = false;
393  }
394 
395  return true;
396  }
397 
398  break;
399 
401  //we arrive here whenever the integrator finds a root
402  //this does not necessarily mean an event fires
403  pModel->setState(*mpCurrentState);
405 
406  //this checks whether equality events are triggered
407  pModel->processRoots(*mpCurrentTime, true, true, mpTrajectoryMethod->getRoots());
408 
409  // TODO Provide a call back method for resolving simultaneous assignments.
410 
411  //execute scheduled events for equalities
412  StateChanged |= pModel->processQueue(*mpCurrentTime, true, NULL);
413 
414  //this checks whether inequality events are triggered
415  pModel->processRoots(*mpCurrentTime, false, true, mpTrajectoryMethod->getRoots());
416 
417  // If the root happens to coincide with end of the step we have to return and
418  // inform the integrator of eventual state changes.
419  if (fabs(*mpCurrentTime - endTime) < Tolerance)
420  {
421  if (StateChanged)
422  {
423  *mpCurrentState = pModel->getState();
425  StateChanged = false;
426  }
427 
428  return true;
429  }
430 
431  break;
432 
434  finish();
436  return false;
437  break;
438  }
439 
440  proceed = mpCallBack == NULL || mpCallBack->proceed();
441 
442  if (mState == FINISH)
443  proceed = false;
444  }
445 
446  return proceed;
447 }
448 
450 {
451  //reset call back
453 
455 
457 }
458 
460 {
461  bool success = CCopasiTask::restore();
462  CModel * pModel = mpProblem->getModel();
463 
464  removeEvent();
465 
466  if (mUpdateModel)
467  {
468 
469  pModel->setState(*mpCurrentState);
471  pModel->setInitialState(pModel->getState());
472  pModel->updateInitialValues();
473  }
474 
475  //reset call back
477 
478  return success;
479 }
480 
481 bool CCrossSectionTask::setMethodType(const int & type)
482 {
484 
485  if (!isValidMethod(Type, ValidMethods)) return false;
486 
487  if (mpMethod->getSubType() == Type) return true;
488 
489  pdelete(mpMethod);
490  mpMethod = createMethod(Type);
491  this->add(mpMethod, true);
492 
493  CCopasiParameter * pParameter = mpMethod->getParameter("Integrate Reduced Model");
494 
495  if (pParameter != NULL)
496  mUpdateMoieties = *pParameter->getValue().pBOOL;
497  else
498  mUpdateMoieties = false;
499 
500  return true;
501 }
502 
503 // virtual
505 {
507 
508  return CTrajectoryMethod::createMethod(Type);
509 }
510 
512 {return mpCurrentState;}
513 
515 {return mTimeSeries;}
516 
517 //static
519 {static_cast<CCrossSectionTask *>(pCSTask)->eventCallBack(type);}
520 
522 {
523  // std::cout << "event call back: " << type << std::endl;
524 
525  //do progress reporting
526  if (mpCallBack != NULL)
527  {
529 
531  {
532  mState = FINISH;
533  }
534  }
535 
536  //do nothing else if the event is not representing a cut plane
537  if (type != CEvent::CutPlane)
538  return;
539 
542 
543  // count the crossings
544  ++mNumCrossings;
545 
546  // now check if we can transition to the main state
547  if (mState == TRANSIENT)
548  {
549  // if output is not delayed by time and not delayed by number of crossings
550  // and also not delayed by convergence
551  // output is started immediately
555  {
556  mState = MAIN;
557  }
559  {
560  mState = MAIN;
561  }
563  {
564  mState = MAIN;
565  }
567  mStatesRingCounter > 0)
568  {
569  int i;
570 
571  for (i = mStatesRingCounter - 1; i >= 0 && i >= ((int) mStatesRingCounter) - RING_SIZE; --i)
572  {
575 
576  if (tmp < mpCrossSectionProblem->getConvergenceOutTolerance())
577  {
579  mPeriod = *mpCurrentTime - mStatesRing[i % RING_SIZE].getTime();
580  mFreq = 1 / mPeriod;
583 
584  mState = MAIN;
585  break;
586  }
587 
588  //std::cout << i << " " << tmp << std::endl;
589  }
590  }
591 
592  if (mState == MAIN)
593  {
594  mStatesRingCounter = 0;
595  mNumCrossings = 1;
596  }
597  }
598 
599  if (mState == MAIN)
600  {
601  //check for convergence (actually for the occurrence of similar states)
602  if (mStatesRingCounter > 0)
603  {
604  int i;
605 
606  for (i = mStatesRingCounter - 1; i >= 0 && i >= ((int) mStatesRingCounter) - RING_SIZE; --i)
607  {
610 
611  if (tmp < mpCrossSectionProblem->getConvergenceTolerance())
612  {
614  mPeriod = *mpCurrentTime - mStatesRing[i % RING_SIZE].getTime();
615  mFreq = 1 / mPeriod;
618 
620  mState = FINISH;
621 
622  break;
623  }
624 
625  //std::cout << "MAIN" << i << " " << tmp << std::endl;
626  }
627  }
628 
629  if (!isnan(mPreviousCrossingTime))
630  {
632  }
633 
634  mLastFreq = 1 / mLastPeriod;
635 
637 
638  //check if the conditions for stopping are met
639  //we don't have to check for maximum duration, this is done elsewhere
641  mState = FINISH;
642  }
643 
644  //store state in ring buffer
648 }
649 
650 //static
652 {
653  if (!s1 || !s2)
654  {
655  return std::numeric_limits< C_FLOAT64 >::quiet_NaN();
656  }
657 
658  if (s1->endIndependent() - s1->beginIndependent()
659  != s2->endIndependent() - s2->beginIndependent())
660  {
661  return std::numeric_limits< C_FLOAT64 >::quiet_NaN();
662  }
663 
664  C_FLOAT64 ret = 0;
665 
666  const C_FLOAT64 * p1 = s1->beginIndependent();
667  const C_FLOAT64 * p1End = s1->endIndependent();
668  const C_FLOAT64 * p2 = s2->beginIndependent();
669 
670  for (; p1 != p1End; ++p1, ++p2)
671  {
672  ret += (*p1 != *p2) ? pow((*p1 - *p2) / (fabs(*p1) + fabs(*p2)), 2) : 0.0;
673  }
674 
675  return 2.0 * sqrt(ret);
676 }
virtual bool restore()
C_FLOAT64 mPreviousCrossingTime
CCopasiMethod * mpMethod
Definition: CCopasiTask.h:239
#define pdelete(p)
Definition: copasi.h:215
const std::string & getSingleObjectCN() const
CEvent * createEvent(const std::string &name)
Definition: CModel.cpp:2928
CCopasiVectorN< CEvent > & getEvents()
Definition: CModel.cpp:1110
void clear()
Definition: CTimeSeries.cpp:96
void setProblem(CTrajectoryProblem *problem)
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
bool getFlagLimitCrossings() const
virtual void stateChanged()
const CCopasiMessage & getMessage() const
virtual bool setName(const std::string &name)
void setInitialState(const CState &state)
Definition: CModel.cpp:1774
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
CProcessReport * mpCallBack
Definition: CCopasiTask.h:249
void processStart(const bool &useInitialValues)
virtual bool setMethodType(const int &type)
void allocate(const size_t &steps)
Definition: CTimeSeries.cpp:73
Definition: CState.h:305
Type
Definition: CEvent.h:155
const CMathModel * getMathModel() const
Definition: CModel.cpp:4722
#define C_INVALID_INDEX
Definition: copasi.h:222
virtual void output(const COutputInterface::Activity &activity)
void setCurrentState(CState *currentState)
virtual size_t getIndex(const std::string &name) const
const CCopasiMethod::SubType & getSubType() const
virtual bool process(const bool &useInitialValues)
void applyInitialValues()
Definition: CModel.cpp:1236
virtual bool progressItem(const size_t &handle)
bool getFlagLimitConvergence() const
virtual bool isValidProblem(const CCopasiProblem *pProblem)
C_FLOAT64 * endIndependent()
Definition: CState.cpp:329
CCopasiProblem * mpProblem
Definition: CCopasiTask.h:234
const CVector< C_INT > & getRoots() const
virtual bool proceed()
virtual void addInterface(COutputInterface *pInterface)
void eventCallBack(CEvent::Type type)
#define MCTrajectoryMethod
const C_FLOAT64 * mpCurrentTime
void setState(const CState &state)
Definition: CModel.cpp:1785
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
const C_FLOAT64 & getDuration() const
void setEventCallBack(void *pTask, EventCallBack ecb)
void processRoots(const C_FLOAT64 &time, const bool &equality, const bool &correct, const CVector< C_INT > &roots)
Definition: CModel.cpp:4692
bool setTriggerExpression(const std::string &expression)
Definition: CEvent.cpp:474
CCrossSectionTask(const CCopasiContainer *pParent=NULL)
const Value & getValue() const
std::vector< CState > mStatesRing
const C_FLOAT64 & getThreshold() const
virtual void start(const CState *initialState)
virtual bool finishItem(const size_t &handle)
CTimeSeries mTimeSeries
virtual CCopasiMethod * createMethod(const int &type) const
CCopasiParameter * getParameter(const std::string &name)
bool mUpdateModel
Definition: CCopasiTask.h:223
const unsigned C_INT32 & getCrossingsLimit() const
CTrajectoryMethod * mpTrajectoryMethod
bool processQueue(const C_FLOAT64 &time, const bool &equality, CProcessQueue::resolveSimultaneousAssignments pResolveSimultaneousAssignments)
Definition: CModel.cpp:4685
static CTrajectoryMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::deterministic)
const unsigned C_INT32 & getOutCrossingsLimit() const
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
const C_FLOAT64 & getOutputStartTime() const
bool getFlagLimitOutCrossings() const
#define C_FLOAT64
Definition: copasi.h:92
bool compileIfNecessary(CProcessReport *pProcessReport)
Definition: CModel.cpp:612
C_FLOAT64 mOutputStartTime
static void EventCallBack(void *pCSTask, CEvent::Type type)
void setType(const Type &type)
Definition: CEvent.cpp:704
static bool isValidMethod(const unsigned int &method, const unsigned int *validMethods)
Definition: CCopasiTask.cpp:89
bool processStep(const C_FLOAT64 &nextTime)
virtual bool add(CCopasiObject *pObject, const bool &adopt=true)
const C_FLOAT64 & getProcessQueueExecutionTime() const
Definition: CModel.cpp:4702
bool updateInitialValues()
Definition: CModel.cpp:1461
bool getFlagLimitOutConvergence() const
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
Definition: CModel.h:50
Header file of class CEvent.
const CState & getState() const
Definition: CModel.cpp:1771
virtual bool restore()
const CProcessQueue & getProcessQueue() const
Definition: CMathModel.cpp:362
virtual Status step(const double &deltaT)
const CTimeSeries & getTimeSeries() const
#define RING_SIZE
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CCrossSectionProblem * mpCrossSectionProblem
CModel * getModel() const
#define min(a, b)
Definition: f2c.h:175
static C_FLOAT64 relativeDifferenceOfStates(CState *s1, CState *s2)
static const unsigned int ValidMethods[]
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
bool removeEvent(const size_t index, const bool &recursive=true)
Definition: CModel.cpp:2945
#define max(a, b)
Definition: f2c.h:176