COPASI API  4.16.103
CLsodaMethod.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 #include "copasi.h"
16 
17 #include "CLsodaMethod.h"
18 #include "CTrajectoryProblem.h"
19 
22 #include "model/CModel.h"
23 #include "model/CState.h"
24 
26  const CCopasiContainer * pParent):
27  CTrajectoryMethod(subType, pParent),
28  mMethodState(),
29  mY(NULL),
30  mRootMask(),
31  mTargetTime(0.0),
32  mRootCounter(0),
33  mPeekAheadMode(false)
34 {
35  assert((void *) &mData == (void *) &mData.dim);
36 
37  mData.pMethod = this;
39 }
40 
42  const CCopasiContainer * pParent):
43  CTrajectoryMethod(src, pParent),
44  mMethodState(),
45  mY(NULL),
46  mRootMask(src.mRootMask),
47  mTargetTime(src.mTargetTime),
48  mRootCounter(src.mRootCounter),
49  mPeekAheadMode(src.mPeekAheadMode)
50 {
51  assert((void *) &mData == (void *) &mData.dim);
52 
53  mData.pMethod = this;
55 }
56 
58 {}
59 
61 {
62  CCopasiParameter *pParm;
63 
65  assertParameter("Integrate Reduced Model", CCopasiParameter::BOOL, (bool) false)->getValue().pBOOL;
67  assertParameter("Relative Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-6)->getValue().pUDOUBLE;
69  assertParameter("Absolute Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-12)->getValue().pUDOUBLE;
71  assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) 10000)->getValue().pUINT;
72 
73  // Check whether we have a method with the old parameter names
74  if ((pParm = getParameter("LSODA.RelativeTolerance")) != NULL)
75  {
77  removeParameter("LSODA.RelativeTolerance");
78 
79  if ((pParm = getParameter("LSODA.AbsoluteTolerance")) != NULL)
80  {
82  removeParameter("LSODA.AbsoluteTolerance");
83  }
84 
85  if ((pParm = getParameter("LSODA.AdamsMaxOrder")) != NULL)
86  {
87  removeParameter("LSODA.AdamsMaxOrder");
88  }
89 
90  if ((pParm = getParameter("LSODA.BDFMaxOrder")) != NULL)
91  {
92  removeParameter("LSODA.BDFMaxOrder");
93  }
94 
95  if ((pParm = getParameter("LSODA.MaxStepsInternal")) != NULL)
96  {
97  *mpMaxInternalSteps = *pParm->getValue().pUINT;
98  removeParameter("LSODA.MaxStepsInternal");
99  }
100  }
101 
102  // Check whether we have a method with "Use Default Absolute Tolerance"
103  if ((pParm = getParameter("Use Default Absolute Tolerance")) != NULL)
104  {
105  C_FLOAT64 NewValue;
106 
107  if (*pParm->getValue().pBOOL)
108  {
109  // The default
110  NewValue = 1.e-12;
111  }
112  else
113  {
114  C_FLOAT64 OldValue = *mpAbsoluteTolerance;
115  CCopasiDataModel* pDataModel = getObjectDataModel();
116  assert(pDataModel != NULL);
117  CModel * pModel = pDataModel->getModel();
118 
119  if (pModel == NULL)
120  // The default
121  NewValue = 1.e-12;
122  else
123  {
124  const CCopasiVectorNS< CCompartment > & Compartment = pModel->getCompartments();
125  size_t i, imax;
127 
128  for (i = 0, imax = Compartment.size(); i < imax; i++)
129  if (Compartment[i]->getValue() < Volume)
130  Volume = Compartment[i]->getValue();
131 
133  // The default
134  NewValue = 1.e-12;
135  else
136  // Invert the scaling as best as we can
137  NewValue = OldValue / (Volume * pModel->getQuantity2NumberFactor());
138  }
139  }
140 
141  *mpAbsoluteTolerance = NewValue;
142  removeParameter("Use Default Absolute Tolerance");
143  }
144 
145  // These parameters are no longer supported.
146  removeParameter("Adams Max Order");
147  removeParameter("BDF Max Order");
148 }
149 
151 {
153  return true;
154 }
155 
156 // virtual
158 {
159  if (!mNoODE && mLsodaStatus != 1)
160  {
161  // Compare the independent state variables
162  // This an be done directly by comparing mMethodState and *mpCurrentState
163  C_FLOAT64 *pMethod = mY;
164  C_FLOAT64 *pMethodEnd = pMethod + mData.dim;
166 
167  for (; pMethod != pMethodEnd; ++pMethod, ++pCurrent)
168  {
169  // We may need to use the absolute and relative tolerance
170  if (*pMethod != *pCurrent)
171  {
172  mLsodaStatus = 1;
173 
174  break;
175  }
176  }
177 
178  // Bug 2106 It is not sufficient to check for the state variables we also need
179  // to check for fixed values which are event targets.
180  pMethod = mMethodState.beginFixed();
181  pMethodEnd = mMethodState.endFixed();
182  pCurrent = mpCurrentState->beginFixed();
183 
184  for (; pMethod != pMethodEnd; ++pMethod, ++pCurrent)
185  {
186  // We may need to use the absolute and relative tolerance
187  if (*pMethod != *pCurrent)
188  {
189  mLsodaStatus = 1;
190 
191  break;
192  }
193  }
194 
195  if (mLsodaStatus != 1)
196  {
197  // Compare the rates of the independent state variables
198  // we need to call evalF for mMethodState and *mpCurrentState and compare the returned rates.
199  CVector< C_FLOAT64 > MethodRate(mData.dim);
200  CVector< C_FLOAT64 > CurrentRate(mData.dim);
201 
202  evalF(&mTime, mY, MethodRate.array());
203 
205  mTime = mMethodState.getTime();
206 
207  evalF(&mTime, mY, CurrentRate.array());
208 
209  pMethod = MethodRate.array();
210  pMethodEnd = pMethod + mData.dim;
211  pCurrent = CurrentRate.array();
212 
213  for (; pMethod != pMethodEnd; ++pMethod, ++pCurrent)
214  {
215  // We may need to use the absolute and relative tolerance
216  if (*pMethod != *pCurrent)
217  {
218  mLsodaStatus = 1;
219  break;
220  }
221  }
222  }
223  }
224 
227 
228  mPeekAheadMode = false;
229  destroyRootMask();
230 }
231 
233 {
234  if (mData.dim == 0 && mNumRoots == 0) //just do nothing if there are no variables
235  {
236  mTime = mTime + deltaT;
239 
240  return NORMAL;
241  }
242 
243  C_FLOAT64 EndTime = mTime + deltaT;
244 
245  if (mTargetTime != EndTime)
246  {
247  // We have a new end time and reset the root counter.
248  mTargetTime = EndTime;
249  mRootCounter = 0;
250  }
251  else
252  {
253  // We are called with the same end time which means a root has previously been
254  // found. We increase the root counter and check whether the limit is reached.
255  mRootCounter++;
256 
258  {
259  return FAILURE;
260  }
261  }
262 
263  C_INT ITOL = 2; // mRtol scalar, mAtol vector
264  C_INT one = 1;
265  C_INT DSize = (C_INT) mDWork.size();
266  C_INT ISize = (C_INT) mIWork.size();
267 
268  // The return status of the integrator.
269  Status Status = NORMAL;
270 
271  if (mRoots.size() > 0)
272  {
273  mLSODAR(&EvalF, // 1. evaluate F
274  &mData.dim, // 2. number of variables
275  mY, // 3. the array of current concentrations
276  &mTime, // 4. the current time
277  &EndTime, // 5. the final time
278  &ITOL, // 6. error control
279  &mRtol, // 7. relative tolerance array
280  mAtol.array(), // 8. absolute tolerance array
281  &mState, // 9. output by overshoot & interpolation
282  &mLsodaStatus, // 10. the state control variable
283  &one, // 11. further options (one)
284  mDWork.array(), // 12. the double work array
285  &DSize, // 13. the double work array size
286  mIWork.array(), // 14. the int work array
287  &ISize, // 15. the int work array size
288  &EvalJ, // 16. evaluate J (not given)
289  &mJType, // 17. type of j evaluation 2 internal full matrix
290  &EvalR, // 18. evaluate constraint functions
291  &mNumRoots, // 19. number of constraint functions g(i)
292  mRoots.array()); // 20. integer array of length NG for output of root information
293 
294  // There exist situations where LSODAR reports status = 3, which are actually status = -33
295  // Obviously the trivial case is where LSODAR did not advance at all, i.e, the start time
296  // equals the current time. It may also happen that a very small steps has been taken in
297  // we reset short before we reach the internal step limit.
298  if (mLsodaStatus == 3 &&
299  (mRootCounter > 0.99 * *mpMaxInternalSteps ||
301  {
302  mLsodaStatus = -33;
303  mRootCounter = 0;
304  }
305 
306  switch (mLsodaStatus)
307  {
308  case -33:
309 
310  switch (mRootMasking)
311  {
312  case NONE:
313  case DISCRETE:
314  // Reset the integrator to the state before the failed integration.
315 
318  mLsodaStatus = 1;
319 
320  // Create a mask which hides all roots being constant and zero.
321  createRootMask();
322  break;
323 
324  case ALL:
325  break;
326  }
327 
328  break;
329 
330  case 3:
331 
332  // If mLsodaStatus == 3 we have found a root. This needs to be indicated to
333  // the caller as it is not sufficient to rely on the fact that T < TOUT
334  if (mLsodaStatus == 3)
335  {
336  // It is sufficient to switch to 2. Eventual state changes due to events
337  // are indicated via the method stateChanged()
338  mLsodaStatus = 2;
339  Status = ROOT;
340  }
341 
342  // To detect simultaneous roots we have to peek ahead, i.e., continue
343  // integration until the state changes are larger that the relative
344  if (!mPeekAheadMode)
345  {
346  Status = peekAhead();
347  }
348 
349  // The break statement is intentionally missing since we
350  // have to continue to check the root masking state.
351  default:
352 
353  switch (mRootMasking)
354  {
355  case NONE:
356  case DISCRETE:
357  break;
358 
359  case ALL:
360  {
361  const bool * pDiscrete = mDiscreteRoots.array();
362  bool * pMask = mRootMask.array();
363  bool * pMaskEnd = pMask + mNumRoots;
364  bool Destroy = true;
365 
366  for (; pMask != pMaskEnd; ++pMask, ++pDiscrete)
367  {
368  if (*pMask)
369  {
370  if (*pDiscrete)
371  {
372  Destroy = false;
373  }
374  else
375  {
376  *pMask = false;
377  }
378  }
379  }
380 
381  if (Destroy)
382  {
383  destroyRootMask();
384  }
385  else
386  {
388  }
389 
390  // We have to restart the integrator
391  mLsodaStatus = 1;
392  }
393 
394  break;
395  }
396 
397  break;
398  }
399  }
400  else
401  {
402  mLSODA(&EvalF, // 1. evaluate F
403  &mData.dim, // 2. number of variables
404  mY, // 3. the array of current concentrations
405  &mTime, // 4. the current time
406  &EndTime, // 5. the final time
407  &ITOL, // 6. error control
408  &mRtol, // 7. relative tolerance array
409  mAtol.array(), // 8. absolute tolerance array
410  &mState, // 9. output by overshoot & interpolation
411  &mLsodaStatus, // 10. the state control variable
412  &one, // 11. further options (one)
413  mDWork.array(), // 12. the double work array
414  &DSize, // 13. the double work array size
415  mIWork.array(), // 14. the int work array
416  &ISize, // 15. the int work array size
417  EvalJ, // 16. evaluate J (not given)
418  &mJType); // 17. the type of jacobian calculate (2)
419  }
420 
421  // Why did we ignore this error?
422  // if (mLsodaStatus == -1) mLsodaStatus = 2;
423 
425 
426  if (!mPeekAheadMode)
427  {
429  }
430 
431  if ((mLsodaStatus <= 0))
432  {
433  Status = FAILURE;
434  mPeekAheadMode = false;
436  }
437 
438  if (!mMethodState.isValid())
439  {
440  Status = FAILURE;
441  mPeekAheadMode = false;
443  }
444 
445  return Status;
446 }
447 
448 void CLsodaMethod::start(const CState * initialState)
449 {
450  /* Retrieve the model to calculate */
452 
453  /* Reset lsoda */
454  mLsodaStatus = 1;
455  mState = 1;
456  mJType = 2;
457  mErrorMsg.str("");
458 
459  /* Release previous state and make the initialState the current */
460  mMethodState = *initialState;
462  mTargetTime = mTime;
463  mRootCounter = 0;
464  mPeekAheadMode = false;
465 
468  destroyRootMask();
469 
470  if (*mpReducedModel)
472  else
474 
475  // When we have roots we need to add an artificial ODE dDummy/dt = 1
476  if (mData.dim == 0 && mNumRoots != 0)
477  {
478  mData.dim = 1;
479  mNoODE = true;
480  mAtol.resize(1);
482  mDummy = 0;
483  mY = &mDummy;
484  }
485  else
486  {
487  mNoODE = false;
490  }
491 
493 
494  /* Configure lsoda(r) */
496 
497  mDWork.resize(22 + mData.dim * std::max<C_INT>(16, mData.dim + 9) + 3 * mNumRoots);
498  mDWork[4] = mDWork[5] = mDWork[6] = mDWork[7] = mDWork[8] = mDWork[9] = 0.0;
499  mIWork.resize(20 + mData.dim);
500  mIWork[4] = mIWork[6] = mIWork[9] = 0;
501 
503  mIWork[7] = 12;
504  mIWork[8] = 5;
505 
506  if (mNumRoots > 0)
507  {
510 
511  CMathTrigger::CRootFinder * const* ppRootFinder = mpModel->getRootFinders().array();
512  CMathTrigger::CRootFinder * const* ppRootFinderEnd = ppRootFinder + mNumRoots;
513  bool * pDiscrete = mDiscreteRoots.array();
514 
515  for (; ppRootFinder != ppRootFinderEnd; ++ppRootFinder, ++pDiscrete)
516  {
517  *pDiscrete = (*ppRootFinder)->isDiscrete();
518  }
519  }
520  else
521  {
523  }
524 
525  return;
526 }
527 
528 void CLsodaMethod::EvalF(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot)
529 {static_cast<Data *>((void *) n)->pMethod->evalF(t, y, ydot);}
530 
531 void CLsodaMethod::evalF(const C_FLOAT64 * t, const C_FLOAT64 * /* y */, C_FLOAT64 * ydot)
532 {
533  // If we have no ODEs add a constant one.
534  if (mNoODE)
535  {
536  *ydot = 1.0;
537  return;
538  }
539 
540  mMethodState.setTime(*t);
541 
544 
545  if (*mpReducedModel)
547  else
549 
550  return;
551 }
552 
553 void CLsodaMethod::EvalR(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y,
554  const C_INT * nr, C_FLOAT64 * r)
555 {static_cast<Data *>((void *) n)->pMethod->evalR(t, y, nr, r);}
556 
557 void CLsodaMethod::evalR(const C_FLOAT64 * t, const C_FLOAT64 * /* y */,
558  const C_INT * nr, C_FLOAT64 * r)
559 {
560  assert(*nr == (C_INT) mRoots.size());
561 
562  mMethodState.setTime(*t);
563 
565 
566  if (*mpReducedModel)
567  {
569  }
570 
571  CVectorCore< C_FLOAT64 > RootValues(*nr, r);
572 
573  mpModel->evaluateRoots(RootValues, true);
574 
575  if (mRootMasking != NONE)
576  {
577  maskRoots(RootValues);
578  }
579 };
580 
581 // static
582 void CLsodaMethod::EvalJ(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y,
583  const C_INT * ml, const C_INT * mu, C_FLOAT64 * pd, const C_INT * nRowPD)
584 {static_cast<Data *>((void *) n)->pMethod->evalJ(t, y, ml, mu, pd, nRowPD);}
585 
586 // virtual
587 void CLsodaMethod::evalJ(const C_FLOAT64 * t, const C_FLOAT64 * y,
588  const C_INT * ml, const C_INT * mu, C_FLOAT64 * pd, const C_INT * nRowPD)
589 {
590  // TODO Implement me.
591 }
592 
594 {
595  const bool *pMask = mRootMask.array();
596  const bool *pMaskEnd = pMask + mRootMask.size();
597  C_FLOAT64 * pRoot = rootValues.array();
598 
599  for (; pMask != pMaskEnd; ++pMask, ++pRoot)
600  {
601  if (*pMask)
602  {
603  *pRoot = 1.0;
604  }
605  }
606 }
607 
609 {
610  size_t NumRoots = mRoots.size();
611  mRootMask.resize(NumRoots);
612  CVector< C_FLOAT64 > RootValues;
613  RootValues.resize(NumRoots);
614  CVector< C_FLOAT64 > RootDerivatives;
615  RootDerivatives.resize(NumRoots);
616 
618 
619  if (*mpReducedModel)
620  {
622  }
623 
624  mpModel->evaluateRoots(RootValues, true);
625  mpModel->calculateRootDerivatives(RootDerivatives);
626 
627  bool *pMask = mRootMask.array();
628  bool *pMaskEnd = pMask + mRootMask.size();
629  C_FLOAT64 * pRootValue = RootValues.array();
630  C_FLOAT64 * pRootDerivative = RootDerivatives.array();
631 
632  for (; pMask != pMaskEnd; ++pMask, ++pRootValue, ++pRootDerivative)
633  {
634  *pMask = (fabs(*pRootDerivative) < *mpAbsoluteTolerance ||
635  fabs(*pRootValue) < 1e3 * std::numeric_limits< C_FLOAT64 >::min()) ? true : false;
636  }
637 
638  mRootMasking = ALL;
639 }
640 
642 {
643  mRootMask.resize(0);
644  mRootMasking = NONE;
645 }
646 
648 {
649  // Save the current state
651 
652  CState StartState = mMethodState;
653 
654  CState ResetState = mMethodState;
655  mLSODAR.saveState();
656  CVector< C_FLOAT64 > ResetDWork = mDWork;
657  CVector< C_INT > ResetIWork = mIWork;
658 
659  mPeekAheadMode = true;
660  Status PeekAheadStatus = ROOT;
661 
662  CVector< C_INT > CombinedRoots = mRoots;
663 
664  C_FLOAT64 MaxPeekAheadTime = std::max(mTargetTime, mTime * (1.0 + 2.0 * *mpRelativeTolerance));
665 
666  while (mPeekAheadMode)
667  {
668  switch (step(MaxPeekAheadTime - mTime))
669  {
670  case ROOT:
671  {
672  // Check whether the new state is within the tolerances
673  C_FLOAT64 * pOld = StartState.beginIndependent();
674  C_FLOAT64 * pOldEnd = pOld + mData.dim;
676  C_FLOAT64 * pAtol = mAtol.array();
677 
678  for (; pOld != pOldEnd; ++pOld, ++pNew, ++pAtol)
679  {
680  if ((2.0 * fabs(*pNew - *pOld) > (fabs(*pNew) + fabs(*pOld)) * *mpRelativeTolerance) &&
681  fabs(*pNew) > *pAtol &&
682  fabs(*pOld) > *pAtol)
683  {
684  break;
685  }
686  }
687 
688  if (pOld != pOldEnd ||
689  (mTime > StartState.getTime() * (1.0 + *mpRelativeTolerance) &&
690  fabs(mTime) > *mpAbsoluteTolerance &&
691  fabs(StartState.getTime()) > *mpAbsoluteTolerance))
692  {
693  mPeekAheadMode = false;
694  }
695  else
696  {
697  ResetState = mMethodState;
698  mLSODAR.saveState();
699  ResetDWork = mDWork;
700  ResetIWork = mIWork;
701 
702  // Combine all the roots
703  C_INT * pRoot = mRoots.array();
704  C_INT * pRootEnd = pRoot + mRoots.size();
705  C_INT * pCombinedRoot = CombinedRoots.array();
706 
707  for (; pRoot != pRootEnd; ++pRoot, ++pCombinedRoot)
708  {
709  if (*pRoot > 0)
710  {
711  *pCombinedRoot = 1;
712  }
713  }
714  }
715  }
716 
717  break;
718 
719  case FAILURE:
720  mPeekAheadMode = false;
721  PeekAheadStatus = FAILURE;
722  break;
723 
724  case NORMAL:
725 
726  if (mRootMasking != ALL)
727  {
728  mPeekAheadMode = false;
729  }
730 
731  break;
732  }
733  }
734 
735  // Reset the integrator to the saved state
736  mMethodState = ResetState;
739  mDWork = ResetDWork;
740  mIWork = ResetIWork;
741 
742  mPeekAheadMode = false;
743 
744  mRoots = CombinedRoots;
745 
746  return PeekAheadStatus;
747 }
CCopasiDataModel * getObjectDataModel()
#define C_INT
Definition: copasi.h:115
C_INT mLsodaStatus
Definition: CLsodaMethod.h:107
static void EvalR(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, const C_INT *nr, C_FLOAT64 *r)
CVector< C_FLOAT64 > mAtol
Definition: CLsodaMethod.h:118
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
C_FLOAT64 mDummy
Definition: CLsodaMethod.h:170
virtual size_t size() const
const CVector< CMathTrigger::CRootFinder * > & getRootFinders() const
Definition: CModel.cpp:4717
C_FLOAT64 * mpAbsoluteTolerance
Definition: CLsodaMethod.h:63
CVector< bool > mRootMask
Definition: CLsodaMethod.h:175
void evaluateRoots(CVectorCore< C_FLOAT64 > &rootValues, const bool &ignoreDiscrete)
Definition: CModel.cpp:4679
void setOstream(std::ostream &os)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
virtual void start(const CState *initialState)
C_FLOAT64 mRtol
Definition: CLsodaMethod.h:113
Definition: CState.h:305
CTrajectoryProblem * mpProblem
CVector< C_FLOAT64 > mDWork
Definition: CLsodaMethod.h:143
virtual void evalR(const C_FLOAT64 *t, const C_FLOAT64 *y, const C_INT *nr, C_FLOAT64 *r)
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
CVector< bool > mDiscreteRoots
Definition: CLsodaMethod.h:180
C_FLOAT64 * mY
Definition: CLsodaMethod.h:86
CModel * mpModel
Definition: CLsodaMethod.h:159
#define C_INT32
Definition: copasi.h:90
void maskRoots(CVectorCore< C_FLOAT64 > &rootValues)
unsigned C_INT32 mRootCounter
Definition: CLsodaMethod.h:199
static void EvalJ(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, const C_INT *ml, const C_INT *mu, C_FLOAT64 *pd, const C_INT *nRowPD)
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
void calculateRootDerivatives(CVector< C_FLOAT64 > &rootDerivatives)
Definition: CModel.cpp:4712
virtual Status step(const double &deltaT)
CState mMethodState
Definition: CLsodaMethod.h:74
C_FLOAT64 mTargetTime
Definition: CLsodaMethod.h:193
std::ostringstream mErrorMsg
Definition: CLsodaMethod.h:123
C_FLOAT64 * endFixed()
Definition: CState.cpp:333
CVector< C_INT > mRoots
bool removeParameter(const std::string &name)
#define MCTrajectoryMethod
C_INT mNumRoots
Definition: CLsodaMethod.h:96
const C_FLOAT64 & getQuantity2NumberFactor() const
Definition: CModel.cpp:2354
CLsodaMethod(const CCopasiMethod::SubType &subType=deterministic, const CCopasiContainer *pParent=NULL)
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
virtual bool elevateChildren()
void setState(const CState &state)
Definition: CModel.cpp:1785
bool * mpReducedModel
Definition: CLsodaMethod.h:52
const Value & getValue() const
CLsodaMethod * pMethod
Definition: CLsodaMethod.h:37
unsigned C_INT32 * pUINT
size_t getNumIndependent() const
Definition: CState.cpp:342
CCopasiParameter * getParameter(const std::string &name)
void destroyRootMask()
bool isValid() const
Definition: CState.cpp:351
C_FLOAT64 mTime
Definition: CLsodaMethod.h:102
unsigned C_INT32 * mpMaxInternalSteps
Definition: CLsodaMethod.h:68
size_t size() const
Definition: CVector.h:100
void calculateDerivativesX(C_FLOAT64 *derivativesX)
Definition: CModel.cpp:1941
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
RootMasking mRootMasking
Definition: CLsodaMethod.h:186
bool mPeekAheadMode
Definition: CLsodaMethod.h:204
CType * array()
Definition: CVector.h:139
Status peekAhead()
void createRootMask()
virtual void stateChanged()
void calculateDerivatives(C_FLOAT64 *derivatives)
Definition: CModel.cpp:1903
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
CLSODAR mLSODAR
Definition: CLsodaMethod.h:133
virtual void evalJ(const C_FLOAT64 *t, const C_FLOAT64 *y, const C_INT *ml, const C_INT *mu, C_FLOAT64 *pd, const C_INT *nRowPD)
CVector< C_FLOAT64 > initializeAtolVector(const C_FLOAT64 &baseTolerance, const bool &reducedModel) const
Definition: CModel.cpp:4368
Definition: CModel.h:50
C_FLOAT64 * beginFixed()
Definition: CState.cpp:332
CVector< C_FLOAT64 > mYdot
Definition: CLsodaMethod.h:91
void initializeParameter()
virtual void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
CVector< C_INT > mIWork
Definition: CLsodaMethod.h:148
size_t getNumRoots() const
Definition: CModel.cpp:4707
CModel * getModel() const
#define min(a, b)
Definition: f2c.h:175
C_FLOAT64 * mpRelativeTolerance
Definition: CLsodaMethod.h:58
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
#define max(a, b)
Definition: f2c.h:176