COPASI API  4.16.103
CTrajectoryTask.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) 2008 - 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 // Copyright (C) 2002 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * CTrajectoryTask class.
17  *
18  * This class implements a trajectory task which is comprised of a
19  * of a problem and a method. Additionally calls to the reporting
20  * methods are done when initialized.
21  *
22  * Created for COPASI by Stefan Hoops 2002
23  */
24 
25 #include <string>
26 #include <cmath>
27 
28 #include "copasi.h"
29 
30 #include "CTrajectoryTask.h"
31 #include "CTrajectoryProblem.h"
32 #include "CTrajectoryMethod.h"
33 #include "model/CModel.h"
34 #include "model/CMathModel.h"
35 #include "model/CModel.h"
36 #include "model/CState.h"
37 #include "report/CKeyFactory.h"
38 #include "report/CReport.h"
42 
43 #define XXXX_Reporting
44 
45 bool fle(const C_FLOAT64 & d1, const C_FLOAT64 & d2)
46 {return (d1 <= d2);}
47 
48 bool fl(const C_FLOAT64 & d1, const C_FLOAT64 & d2)
49 {return (d1 < d2);}
50 
51 bool ble(const C_FLOAT64 & d1, const C_FLOAT64 & d2)
52 {return (d1 >= d2);}
53 
54 bool bl(const C_FLOAT64 & d1, const C_FLOAT64 & d2)
55 {return (d1 > d2);}
56 
57 const unsigned int CTrajectoryTask::ValidMethods[] =
58 {
66 #ifdef COPASI_DEBUG
69 #endif // COPASI_DEBUG
71 };
72 
74  CCopasiTask(CCopasiTask::timeCourse, pParent),
75  mTimeSeriesRequested(true),
76  mTimeSeries(),
77  mpTrajectoryProblem(NULL),
78  mpTrajectoryMethod(NULL),
79  mUpdateMoieties(false),
80  mpCurrentState(NULL),
81  mpCurrentTime(NULL),
82  mOutputStartTime(0.0),
83  mpLessOrEqual(&fle),
84  mpLess(&fl)
85 {
86  mpProblem = new CTrajectoryProblem(this);
88  this->add(mpMethod, true);
89 
90  CCopasiParameter * pParameter = mpMethod->getParameter("Integrate Reduced Model");
91 
92  if (pParameter != NULL)
93  mUpdateMoieties = *pParameter->getValue().pBOOL;
94  else
95  mUpdateMoieties = false;
96 }
97 
99  const CCopasiContainer * pParent):
100  CCopasiTask(src, pParent),
101  mTimeSeriesRequested(src.mTimeSeriesRequested),
102  mTimeSeries(),
103  mpTrajectoryProblem(NULL),
104  mpTrajectoryMethod(NULL),
105  mUpdateMoieties(false),
106  mpCurrentState(NULL),
107  mpCurrentTime(NULL),
108  mOutputStartTime(0.0),
109  mpLessOrEqual(src.mpLessOrEqual),
110  mpLess(src.mpLess)
111 {
112  mpProblem =
113  new CTrajectoryProblem(*static_cast< CTrajectoryProblem * >(src.mpProblem), this);
114 
116  * mpMethod = * src.mpMethod;
117 
118  mpMethod->elevateChildren();
119 
120  this->add(mpMethod, true);
121 
122  CCopasiParameter * pParameter = mpMethod->getParameter("Integrate Reduced Model");
123 
124  if (pParameter != NULL)
125  mUpdateMoieties = *pParameter->getValue().pBOOL;
126  else
127  mUpdateMoieties = false;
128 }
129 
131 {
132  cleanup();
133 }
134 
136 {
138 }
139 
140 void CTrajectoryTask::load(CReadConfig & configBuffer)
141 {
142  configBuffer.getVariable("Dynamics", "bool", &mScheduled,
144 
146  mpProblem = new CTrajectoryProblem(this);
147  ((CTrajectoryProblem *) mpProblem)->load(configBuffer);
148 
149  pdelete(mpMethod);
151  this->add(mpMethod, true);
152 
153  CCopasiParameter * pParameter = mpMethod->getParameter("Integrate Reduced Model");
154 
155  if (pParameter != NULL)
156  mUpdateMoieties = *pParameter->getValue().pBOOL;
157 
159 }
160 
162  COutputHandler * pOutputHandler,
163  std::ostream * pOstream)
164 {
165  assert(mpProblem && mpMethod);
166 
168  assert(mpTrajectoryProblem);
169 
170  mpTrajectoryMethod = dynamic_cast<CTrajectoryMethod *>(mpMethod);
171  assert(mpTrajectoryMethod);
172 
174 
175  bool success = mpMethod->isValidProblem(mpProblem);
176 
177  CCopasiParameter * pParameter = mpMethod->getParameter("Integrate Reduced Model");
178 
179  if (pParameter != NULL)
180  mUpdateMoieties = *pParameter->getValue().pBOOL;
181  else
182  mUpdateMoieties = false;
183 
187 
188  //mpOutputStartTime = & mpTrajectoryProblem->getOutputStartTime();
189 
190  // Handle the time series as a regular output.
192 
193  if ((pOutputHandler != NULL) &&
196  {
198  pOutputHandler->addInterface(&mTimeSeries);
199  }
200  else
201  {
202  mTimeSeries.clear();
203  }
204 
206 
207  if (!CCopasiTask::initialize(of, pOutputHandler, pOstream)) success = false;
208 
209  return success;
210 }
211 
212 bool CTrajectoryTask::process(const bool & useInitialValues)
213 {
214  //*****
215 
216  processStart(useInitialValues);
217 
218  //*****
219 
220  //size_t FailCounter = 0;
221 
224  C_FLOAT64 StepNumber = fabs(Duration) / StepSize;
225 
226  if (isnan(StepNumber) || StepNumber < 1.0)
227  StepNumber = 1.0;
228 
229  //the output starts only after "outputStartTime" has passed
230  if (useInitialValues)
232  else
234 
235  C_FLOAT64 NextTimeToReport;
236 
237  const C_FLOAT64 EndTime = *mpCurrentTime + Duration;
238  const C_FLOAT64 StartTime = *mpCurrentTime;
239  C_FLOAT64 CompareEndTime;
240 
241  if (StepSize < 0.0)
242  {
243  mpLessOrEqual = &ble;
244  mpLess = &bl;
245 
246  // It suffices to reach the end time within machine precision
247  CompareEndTime = EndTime + 100.0 * (fabs(EndTime) * std::numeric_limits< C_FLOAT64 >::epsilon() + std::numeric_limits< C_FLOAT64 >::min());
248  }
249  else
250  {
251  mpLessOrEqual = &fle;
252  mpLess = &fl;
253 
254  // It suffices to reach the end time within machine precision
255  CompareEndTime = EndTime - 100.0 * (fabs(EndTime) * std::numeric_limits< C_FLOAT64 >::epsilon() + std::numeric_limits< C_FLOAT64 >::min());
256  }
257 
258  unsigned C_INT32 StepCounter = 1;
259 
260  if (StepSize == 0.0 && Duration != 0.0)
261  {
263  return false;
264  }
265 
267 
268  bool flagProceed = true;
269  C_FLOAT64 handlerFactor = 100.0 / Duration;
270 
271  C_FLOAT64 Percentage = 0;
272  size_t hProcess;
273 
274  if (mpCallBack != NULL && StepNumber > 1.0)
275  {
276  mpCallBack->setName("performing simulation...");
277  C_FLOAT64 hundred = 100;
278  hProcess = mpCallBack->addItem("Completion",
279  Percentage,
280  &hundred);
281  }
282 
284 
285  try
286  {
287  do
288  {
289  // This is numerically more stable then adding
290  // mpTrajectoryProblem->getStepSize().
291  NextTimeToReport =
292  StartTime + (EndTime - StartTime) * StepCounter++ / StepNumber;
293 
294  flagProceed &= processStep(NextTimeToReport);
295 
296  if (mpCallBack != NULL && StepNumber > 1.0)
297  {
298  Percentage = (*mpCurrentTime - StartTime) * handlerFactor;
299  flagProceed &= mpCallBack->progressItem(hProcess);
300  }
301 
303  {
305  }
306  }
307  while ((*mpLess)(*mpCurrentTime, CompareEndTime) && flagProceed);
308  }
309 
310  catch (int)
311  {
314 
316  {
318  }
319 
320  if (mpCallBack != NULL && StepNumber > 1.0) mpCallBack->finishItem(hProcess);
321 
323 
325  }
326 
327  catch (CCopasiException & Exception)
328  {
331 
333  {
335  }
336 
337  if (mpCallBack != NULL && StepNumber > 1.0) mpCallBack->finishItem(hProcess);
338 
340 
341  throw CCopasiException(Exception.getMessage());
342  }
343 
344  if (mpCallBack != NULL && StepNumber > 1.0) mpCallBack->finishItem(hProcess);
345 
347 
348  return true;
349 }
350 
351 void CTrajectoryTask::processStart(const bool & useInitialValues)
352 {
353  if (useInitialValues)
355 
357 
360 
361  return;
362 }
363 
365 {
366  CModel * pModel = mpTrajectoryProblem->getModel();
367  bool StateChanged = false;
368  bool proceed = true;
369 
370  C_FLOAT64 Tolerance = 100.0 * (fabs(endTime) * std::numeric_limits< C_FLOAT64 >::epsilon() + std::numeric_limits< C_FLOAT64 >::min());
371  C_FLOAT64 NextTime = endTime;
372 
373  while (proceed)
374  {
375  // TODO Provide a call back method for resolving simultaneous assignments.
376  StateChanged = pModel->processQueue(*mpCurrentTime, false, NULL);
377 
378  if (StateChanged)
379  {
382  {
384  }
385 
386  *mpCurrentState = pModel->getState();
388  StateChanged = false;
389  }
390 
391  // std::min suffices since events are only supported in forward integration.
392  NextTime = std::min(endTime, pModel->getProcessQueueExecutionTime());
393 
394  switch (mpTrajectoryMethod->step(NextTime - *mpCurrentTime))
395  {
397  pModel->setState(*mpCurrentState);
399 
403  {
405  }
406 
407  // TODO Provide a call back method for resolving simultaneous assignments.
408  StateChanged = pModel->processQueue(*mpCurrentTime, true, NULL);
409 
410  if (fabs(*mpCurrentTime - endTime) < Tolerance)
411  {
412  if (StateChanged)
413  {
414  *mpCurrentState = pModel->getState();
416  StateChanged = false;
417  }
418 
419  return true;
420  }
421 
422  if (StateChanged)
423  {
426  {
428  }
429 
430  *mpCurrentState = pModel->getState();
432  StateChanged = false;
433  }
434 
435  break;
436 
438  pModel->setState(*mpCurrentState);
440 
441  pModel->processRoots(*mpCurrentTime, true, true, mpTrajectoryMethod->getRoots());
442 
446  {
448  }
449 
450  // TODO Provide a call back method for resolving simultaneous assignments.
451  StateChanged = pModel->processQueue(*mpCurrentTime, true, NULL);
452 
453  pModel->processRoots(*mpCurrentTime, false, true, mpTrajectoryMethod->getRoots());
454 
455  // If the root happens to coincide with end of the step we have to return
456  if (fabs(*mpCurrentTime - endTime) < Tolerance)
457  {
458  if (StateChanged)
459  {
460  *mpCurrentState = pModel->getState();
462  StateChanged = false;
463  }
464 
465  return true;
466  }
467 
469  (StateChanged || *mpCurrentTime == pModel->getProcessQueueExecutionTime()) &&
471  {
473  }
474 
475  if (StateChanged)
476  {
477  *mpCurrentState = pModel->getState();
479  StateChanged = false;
480  }
481 
482  break;
483 
486 
487  return false;
488  break;
489  }
490 
491  proceed = mpCallBack == NULL || mpCallBack->proceed();
492  }
493 
494  return proceed;
495 }
496 
498 {
499  bool success = CCopasiTask::restore();
500 
501  if (mUpdateModel)
502  {
503  CModel * pModel = mpProblem->getModel();
504 
505  pModel->setState(*mpCurrentState);
507  pModel->setInitialState(pModel->getState());
508  pModel->updateInitialValues();
509  }
510 
511  return success;
512 }
513 
514 bool CTrajectoryTask::setMethodType(const int & type)
515 {
517 
518  if (!isValidMethod(Type, ValidMethods)) return false;
519 
520  if (mpMethod->getSubType() == Type) return true;
521 
522  pdelete(mpMethod);
523  mpMethod = createMethod(Type);
524  this->add(mpMethod, true);
525 
526  CCopasiParameter * pParameter = mpMethod->getParameter("Integrate Reduced Model");
527 
528  if (pParameter != NULL)
529  mUpdateMoieties = *pParameter->getValue().pBOOL;
530  else
531  mUpdateMoieties = false;
532 
533  return true;
534 }
535 
536 // virtual
538 {
540 
541  return CTrajectoryMethod::createMethod(Type);
542 }
543 
545 {return mpCurrentState;}
546 
548 {return mTimeSeries;}
CCopasiMethod * mpMethod
Definition: CCopasiTask.h:239
#define pdelete(p)
Definition: copasi.h:215
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
CTimeSeries mTimeSeries
static const unsigned int ValidMethods[]
void clear()
Definition: CTimeSeries.cpp:96
void setProblem(CTrajectoryProblem *problem)
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
const unsigned C_INT32 & getStepNumber() const
virtual void stateChanged()
bool mScheduled
Definition: CCopasiTask.h:217
const CCopasiMessage & getMessage() const
virtual bool setName(const std::string &name)
void setInitialState(const CState &state)
Definition: CModel.cpp:1774
bool processStep(const C_FLOAT64 &nextTime)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
CProcessReport * mpCallBack
Definition: CCopasiTask.h:249
CTrajectoryProblem * mpTrajectoryProblem
void allocate(const size_t &steps)
Definition: CTimeSeries.cpp:73
Definition: CState.h:305
virtual bool process(const bool &useInitialValues)
const CMathModel * getMathModel() const
Definition: CModel.cpp:4722
virtual void output(const COutputInterface::Activity &activity)
void setCurrentState(CState *currentState)
CState * mpCurrentState
const CCopasiMethod::SubType & getSubType() const
#define C_INT32
Definition: copasi.h:90
void processStart(const bool &useInitialValues)
void applyInitialValues()
Definition: CModel.cpp:1236
virtual bool progressItem(const size_t &handle)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
const C_FLOAT64 * mpCurrentTime
CCopasiProblem * mpProblem
Definition: CCopasiTask.h:234
const CVector< C_INT > & getRoots() const
virtual bool proceed()
virtual void addInterface(COutputInterface *pInterface)
#define MCTrajectoryMethod
const C_FLOAT64 & getStepSize() const
bool fle(const C_FLOAT64 &d1, const C_FLOAT64 &d2)
#define MCTrajectoryProblem
bool ble(const C_FLOAT64 &d1, const C_FLOAT64 &d2)
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 setContinueSimultaneousEvents(const bool &continueSimultaneousEvents)
C_FLOAT64 mOutputStartTime
void processRoots(const C_FLOAT64 &time, const bool &equality, const bool &correct, const CVector< C_INT > &roots)
Definition: CModel.cpp:4692
bool(* mpLessOrEqual)(const C_FLOAT64 &, const C_FLOAT64 &)
const Value & getValue() const
virtual void start(const CState *initialState)
virtual bool finishItem(const size_t &handle)
CCopasiParameter * getParameter(const std::string &name)
bool mUpdateModel
Definition: CCopasiTask.h:223
const CTimeSeries & getTimeSeries() const
bool processQueue(const C_FLOAT64 &time, const bool &equality, CProcessQueue::resolveSimultaneousAssignments pResolveSimultaneousAssignments)
Definition: CModel.cpp:4685
static CTrajectoryMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::deterministic)
void load(CReadConfig &configBuffer)
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
bool(* mpLess)(const C_FLOAT64 &, const C_FLOAT64 &)
const C_FLOAT64 & getOutputStartTime() const
bool timeSeriesRequested() const
#define C_FLOAT64
Definition: copasi.h:92
bool bl(const C_FLOAT64 &d1, const C_FLOAT64 &d2)
bool fl(const C_FLOAT64 &d1, const C_FLOAT64 &d2)
static bool isValidMethod(const unsigned int &method, const unsigned int *validMethods)
Definition: CCopasiTask.cpp:89
virtual bool add(CCopasiObject *pObject, const bool &adopt=true)
const C_FLOAT64 & getProcessQueueExecutionTime() const
Definition: CModel.cpp:4702
bool updateInitialValues()
Definition: CModel.cpp:1461
Definition: CModel.h:50
virtual bool restore()
C_INT32 getVariable(const std::string &name, const std::string &type, void *pout, CReadConfig::Mode mode=CReadConfig::NEXT)
Definition: CReadConfig.cpp:81
const CState & getState() const
Definition: CModel.cpp:1771
virtual bool restore()
virtual CCopasiMethod * createMethod(const int &type) const
const CProcessQueue & getProcessQueue() const
Definition: CMathModel.cpp:362
virtual Status step(const double &deltaT)
const bool & getContinueSimultaneousEvents() const
virtual bool setMethodType(const int &type)
CTrajectoryMethod * mpTrajectoryMethod
CModel * getModel() const
#define min(a, b)
Definition: f2c.h:175
CTrajectoryTask(const CCopasiContainer *pParent=NULL)
const bool & getOutputEvent() const