COPASI API  4.16.103
CLyapWolfMethod.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 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) 2006 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <string.h>
16 #include <cmath>
17 
18 #include "copasi.h"
19 
20 #include "CLyapMethod.h"
21 #include "CLyapProblem.h"
22 #include "CLyapTask.h"
23 
25 #include "model/CModel.h"
26 #include "model/CState.h"
27 #include "lapack/blaswrap.h"
28 
31  mpState(NULL)
32 {
33  assert((void *) &mData == (void *) &mData.dim);
34 
35  mData.pMethod = this;
37 }
38 
40  const CCopasiContainer * pParent):
41  CLyapMethod(src, pParent),
42  mpState(NULL)
43 {
44  assert((void *) &mData == (void *) &mData.dim);
45 
46  mData.pMethod = this;
48 }
49 
51 {
53 }
54 
56 {
57  assertParameter("Orthonormalization Interval", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0);
58  assertParameter("Overall time", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1000.0);
59  assertParameter("Relative Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-6);
60  assertParameter("Absolute Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-12);
61  assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) 10000);
62 
63  // Check whether we have an (obsolete) parameter "Use Default Absolute Tolerance"
64  CCopasiParameter *pParm;
65 
66  if ((pParm = getParameter("Use Default Absolute Tolerance")) != NULL)
67  {
68  C_FLOAT64 NewValue;
69 
70  if (*pParm->getValue().pBOOL)
71  {
72  // The default
73  NewValue = 1.e-12;
74  }
75  else
76  {
77  NewValue = *getValue("Absolute Tolerance").pUDOUBLE;
78  }
79 
80  setValue("Absolute Tolerance", NewValue);
81  removeParameter("Use Default Absolute Tolerance");
82  }
83 
84  // These parameters are no longer supported.
85  removeParameter("Adams Max Order");
86  removeParameter("BDF Max Order");
87 }
88 
90 {
92  return true;
93 }
94 
95 double CLyapWolfMethod::step(const double & deltaT)
96 {
97  if (!mData.dim) //just do nothing if there are no variables
98  {
99  mTime = mTime + deltaT;
101 
102  return deltaT;
103  }
104 
105  C_FLOAT64 startTime = mTime;
106  C_FLOAT64 EndTime = mTime + deltaT;
107  C_INT one = 1;
108  C_INT two = 2;
109  C_INT DSize = (C_INT) mDWork.size();
110  C_INT ISize = (C_INT) mIWork.size();
111 
112  mLSODA(&EvalF , // 1. evaluate F
113  &mData.dim , // 2. number of variables
114  mVariables.array(), // 3. the array of current concentrations
115  &mTime , // 4. the current time
116  &EndTime , // 5. the final time
117  &two , // 6. vector absolute error, scalar relative error
118  &mRtol , // 7. relative tolerance array
119  mAtol.array() , // 8. absolute tolerance array
120  &mState , // 9. output by overshoot & interpolation
121  &mLsodaStatus , // 10. the state control variable
122  &one , // 11. futher options (one)
123  mDWork.array() , // 12. the double work array
124  &DSize , // 13. the double work array size
125  mIWork.array() , // 14. the int work array
126  &ISize , // 15. the int work array size
127  NULL , // 16. evaluate J (not given)
128  &mJType); // 17. the type of Jacobian calculate (2)
129 
130  if (mLsodaStatus == -1) mLsodaStatus = 2;
131 
132  if ((mLsodaStatus != 1) && (mLsodaStatus != 2) && (mLsodaStatus != -1))
133  {
135  }
136 
137  return mTime - startTime;
138 }
139 
140 void CLyapWolfMethod::start(/*const CState * initialState*/)
141 {
142  /* Reset lsoda */
143  mLsodaStatus = 1;
144  mState = 1;
145  mJType = 2;
146  mErrorMsg.str("");
148 
149  /* Release previous state and make the initialState the current */
150  pdelete(mpState);
152 
153  //mY = mpState->beginIndependent();
154  mTime = mpState->getTime();
155 
156  mReducedModel = true; /* *getValue("Integrate Reduced Model").pBOOL;*/
157 
158  if (mReducedModel)
160  else
162 
165 
166  //calculate the number of variables for lsoda integration
167  if (mDoDivergence)
168  mData.dim = (C_INT)(mSystemSize * (1 + mNumExp) + 1);
169  else
170  mData.dim = (C_INT)(mSystemSize * (1 + mNumExp));
171 
172  //reserve space for exponents. The vectors in the task are resized by the task because they
173  //need to have a minimum size defined in the task
174  //pTask->mLocalExponents.resize(mNumExp);
177  //pTask->mExponents.resize(mNumExp);
178 
179  //initialize the vector on which lsoda will work
182 
183  //generate base vectors. Just define some arbitrary starting vectors that are not too specific and orthonormalize
184  // first fill the array with 0.1
185  C_FLOAT64 *dbl, *dblEnd = mVariables.array() + mData.dim;
186 
187  for (dbl = mVariables.array() + mSystemSize; dbl != dblEnd; ++dbl)
188  *dbl = 0.01;
189 
190  //now add 1.0
191  if (mNumExp > 0)
192  for (dbl = mVariables.array() + mSystemSize; dbl < dblEnd; dbl += (mSystemSize + 1))
193  * dbl = 1.0;
194 
195  orthonormalize();
196 
197  //reserve space for jacobian
199 
200  size_t i, imax = mNumExp;
201 
202  for (i = 0; i < imax; ++i)
203  {
204  mpTask->mLocalExponents[i] = 0;
205  mSumExponents[i] = 0;
206  mpTask->mExponents[i] = 0;
207  }
208 
210  mSumDivergence = 0;
212 
213  /* Configure lsoda */
214  mRtol = * getValue("Relative Tolerance").pUDOUBLE;
215 
216  C_FLOAT64 * pTolerance = getValue("Absolute Tolerance").pUDOUBLE;
218 
220 
221  for (i = 0; i < mSystemSize; ++i)
222  mAtol[i] = tmpAtol[i];
223 
224  for (i = mSystemSize; (int)i < mData.dim; ++i)
225  mAtol[i] = 1e-20;
226 
227  mDWork.resize(22 + mData.dim * std::max<C_INT>(16, mData.dim + 9));
228  mDWork[4] = mDWork[5] = mDWork[6] = mDWork[7] = mDWork[8] = mDWork[9] = 0.0;
229  mIWork.resize(20 + mData.dim);
230  mIWork[4] = mIWork[6] = mIWork[9] = 0;
231 
232  mIWork[5] = * getValue("Max Internal Steps").pUINT;
233  mIWork[7] = 12;
234  mIWork[8] = 5;
235 
236  return;
237 }
238 
239 void CLyapWolfMethod::EvalF(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot)
240 {static_cast<Data *>((void *) n)->pMethod->evalF(t, y, ydot);}
241 
242 void CLyapWolfMethod::evalF(const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot)
243 {
244  assert(y == mVariables.array());
245 
246  //set time in model
247  mpState->setTime(*t);
248 
249  //copy working array to state
251 
252  //copy state to model
253  CModel * pModel = mpProblem->getModel();
254  pModel->setState(*mpState);
256 
257  //the model
258  if (mReducedModel)
259  pModel->calculateDerivativesX(ydot);
260  else
261  pModel->calculateDerivatives(ydot);
262 
263  //the linearized model
264 
265  //get jacobian
266  pModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
267 
268  //empty dummy entries... to be removed later
269  //C_FLOAT64 *dbl, *dblEnd = ydot + mData.dim;
270  //for (dbl = ydot + mSystemSize; dbl != dblEnd; ++dbl)
271  // *dbl = 0;
272 
273  // char T = 'N';
274  // C_INT M = mSystemSize;
275  // C_INT N = mSystemSize;
276  // C_INT K = 1;
277  // //C_INT LD = mUnscaledElasticities.numCols();
278  //
279  // C_FLOAT64 Alpha = 1.0;
280  // C_FLOAT64 Beta = 0.0;
281  //
282  // dgemm_(&T, &T, &K, &N, &M, &Alpha,
283  // const_cast<C_FLOAT64*>(y) + mSystemSize, &K,
284  // mJacobian.array(), &M,
285  // &Beta,
286  // ydot + mSystemSize, &M);
287 
288  C_FLOAT64 *dbl1;
289  const C_FLOAT64 *dbl2, *dbl3, *dbl1end, *dbl3end;
290 
291  dbl1 = ydot + mSystemSize;
292 
293  size_t i;
294 
295  for (i = 1; i <= mNumExp; ++i)
296  {
297  //dbl1 += mSystemSize;
298  dbl1end = dbl1 + mSystemSize;
299  dbl2 = mJacobian.array();
300 
301  for (; dbl1 != dbl1end; ++dbl1)
302  {
303  *dbl1 = 0.0;
304 
305  dbl3 = y + i * mSystemSize;
306  dbl3end = dbl3 + mSystemSize;
307 
308  for (; dbl3 != dbl3end; ++dbl3, ++dbl2)
309  *dbl1 += *dbl2 * *dbl3;
310  }
311  }
312 
313  if (!mDoDivergence) return;
314 
315  //divergence; trace of jacobian
316  *dbl1 = 0;
317  dbl2 = mJacobian.array();
318 
319  for (i = 0; i < mSystemSize; ++i, dbl2 += (mSystemSize + 1))
320  * dbl1 += *dbl2;
321 
322  return;
323 }
324 
325 //this starts with current state
327 {
328  //get the pointer to the parent task. It's needed for output
329  mpTask = dynamic_cast<CLyapTask*>(getObjectParent());
330  assert(mpTask);
331 
332  //initialize LSODA
333  start();
334 
335  C_FLOAT64 stepSize = *getValue("Orthonormalization Interval").pUDOUBLE;
336  C_FLOAT64 transientTime = mpProblem->getTransientTime() + mTime;
337  C_FLOAT64 endTime = mTime + *getValue("Overall time").pUDOUBLE;
338  C_FLOAT64 startTime = mTime;
339 
340  bool flagProceed = true;
341  C_FLOAT64 handlerFactor = 100.0 / (endTime - startTime);
342 
343  //** do the transient **
344  C_FLOAT64 CompareTime = transientTime - 100.0 * fabs(transientTime) * std::numeric_limits< C_FLOAT64 >::epsilon();
345 
346  if (mTime < CompareTime)
347  {
348  do
349  {
350  step(transientTime - mTime);
351 
352  if (mTime > CompareTime) break;
353 
354  /* Here we will do conditional event processing */
355 
356  /* Currently this is correct since no events are processed. */
357  //CCopasiMessage(CCopasiMessage::EXCEPTION, MCTrajectoryMethod + 12);
358 
359  flagProceed &= mpTask->methodCallback((mTime - startTime) * handlerFactor, true);
360  }
361  while (flagProceed);
362 
363  //mpLyapProblem->getModel()->setState(*mpCurrentState);
364  //mpLyapProblem->getModel()->updateSimulatedValues();
365  }
366  else
367  {
368  }
369 
370  //copy working array to state
373 
374  //copy state to model and do output
377  mpTask->methodCallback((mTime - startTime) * handlerFactor, false);
378  //********
379 
380  orthonormalize();
381 
382  if (mDoDivergence)
383  *(mVariables.array() + mVariables.size() - 1) = 0; //divergence
384 
385  mLsodaStatus = 1; //the state has changed, we need to restart lsoda
386 
387  size_t i;
388 
389  C_FLOAT64 realStepSize;
390 
391  do
392  {
393  realStepSize = step(stepSize);
394 
395  orthonormalize();
396  mLsodaStatus = 1; //the state has changed, we need to restart lsoda
397 
398  //process results of orthonormalisation
399  for (i = 0; i < mNumExp; ++i)
400  {
401  mpTask->mLocalExponents[i] = log(mNorms[i]);
403  mpTask->mLocalExponents[i] = mpTask->mLocalExponents[i] / realStepSize;
404  mpTask->mExponents[i] = mSumExponents[i] / (mTime - transientTime);
405  }
406 
407  //process result of divergence integration
408  if (mDoDivergence)
409  {
411  mpTask->mIntervalDivergence = *(mVariables.array() + mVariables.size() - 1) / realStepSize;
412  *(mVariables.array() + mVariables.size() - 1) = 0;
413  mpTask->mAverageDivergence = mSumDivergence / (mTime - transientTime);
414  }
415 
416  //copy working array to state
419  //copy state to model and do output
422  flagProceed &= mpTask->methodCallback((mTime - startTime) * handlerFactor, false);
423  }
424  while ((mTime < endTime) && flagProceed);
425 
426  return flagProceed;
427 }
428 
430 {
431  if (mNumExp < 1) return;
432 
433  //TODO generalize
434  C_FLOAT64 *dbl, *dblEnd;
435 
436  dbl = mVariables.array() + mSystemSize;
437  dblEnd = dbl + mSystemSize;
438  mNorms[0] = norm(dbl, dblEnd);
439  scalarmult(dbl, dblEnd, 1 / mNorms[0]);
440 
441  size_t i, j;
442 
443  for (i = 1; i < mNumExp; ++i)
444  {
445  dbl += mSystemSize;
446  dblEnd = dbl + mSystemSize;
447 
448  //orthogonalisation
449  for (j = 0; j < i; ++j)
450  {
451  add(dbl, dblEnd,
452  -product(dbl, dblEnd, mVariables.array() + (j + 1)*mSystemSize),
453  mVariables.array() + (j + 1)*mSystemSize);
454  }
455 
456  //normalisation
457  mNorms[i] = norm(dbl, dblEnd);
458  scalarmult(dbl, dblEnd, 1 / mNorms[i]);
459  }
460 }
461 
462 //static
464 {
465  C_FLOAT64 sum = 0;
466 
467  for (; dbl1 != dbl2; ++dbl1)
468  sum += *dbl1 * *dbl1;
469 
470  return sqrt(sum);
471 }
472 
473 //static
475  const C_FLOAT64 & f)
476 {
477  for (; dbl1 != dbl2; ++dbl1)
478  *dbl1 *= f;
479 }
480 
481 //static
483  const C_FLOAT64* dbl2)
484 {
485  C_FLOAT64 sum = 0;
486 
487  for (; dbl1 != dbl1End; ++dbl1, ++dbl2)
488  sum += *dbl1 * *dbl2;
489 
490  return sum;
491 }
492 
493 //static
494 void CLyapWolfMethod::add(C_FLOAT64* dbl1, const C_FLOAT64* dbl1End,
495  const C_FLOAT64 & f, const C_FLOAT64* dbl2)
496 {
497  //calculate v1 = v1 + f * v2
498  for (; dbl1 != dbl1End; ++dbl1, ++dbl2)
499  *dbl1 += f * *dbl2;
500 }
501 
502 //virtual
504 {
505  if (!CLyapMethod::isValidProblem(pProblem)) return false;
506 
507  const CLyapProblem * pLP = dynamic_cast<const CLyapProblem *>(pProblem);
508  assert(pLP);
509 
510  C_FLOAT64 stepSize = *getValue("Orthonormalization Interval").pUDOUBLE;
511  C_FLOAT64 transientTime = pLP->getTransientTime();
512  C_FLOAT64 endTime = *getValue("Overall time").pUDOUBLE;
513 
514  if (transientTime >= endTime)
515  {
516  //
518  return false;
519  }
520 
521  if (stepSize > (endTime - transientTime))
522  {
523  //
525  return false;
526  }
527 
528  return true;
529 }
static C_FLOAT64 norm(const C_FLOAT64 *dbl1, const C_FLOAT64 *dbl2)
virtual bool elevateChildren()
#define C_INT
Definition: copasi.h:115
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
static C_FLOAT64 product(const C_FLOAT64 *dbl1, const C_FLOAT64 *dbl1End, const C_FLOAT64 *dbl2)
#define pdelete(p)
Definition: copasi.h:215
CVector< C_FLOAT64 > mLocalExponents
Definition: CLyapTask.h:58
CLyapWolfMethod * pMethod
CVector< C_FLOAT64 > mAtol
static void add(C_FLOAT64 *dbl1, const C_FLOAT64 *dbl1End, const C_FLOAT64 &f, const C_FLOAT64 *dbl2)
void setOstream(std::ostream &os)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
void calculateJacobianX(CMatrix< C_FLOAT64 > &jacobianX, const C_FLOAT64 &derivationFactor, const C_FLOAT64 &resolution)
Definition: CModel.cpp:2082
Definition: CState.h:305
bool methodCallback(const C_FLOAT64 &percentage, bool onlyProgress)
Definition: CLyapTask.cpp:260
CVector< C_FLOAT64 > mNorms
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
CLyapProblem * mpProblem
Definition: CLyapMethod.h:74
#define C_INT32
Definition: copasi.h:90
C_FLOAT64 mSumDivergence
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
virtual void start()
CVector< C_INT > mIWork
virtual double step(const double &deltaT)
bool removeParameter(const std::string &name)
#define MCLyap
#define MCTrajectoryMethod
CLyapWolfMethod(const CCopasiContainer *pParent=NULL)
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
void setState(const CState &state)
Definition: CModel.cpp:1785
const C_FLOAT64 & getTransientTime() const
C_FLOAT64 mIntervalDivergence
Definition: CLyapTask.h:67
unsigned C_INT32 mNumExp
const Value & getValue() const
static void scalarmult(C_FLOAT64 *dbl1, const C_FLOAT64 *dbl2, const C_FLOAT64 &f)
std::ostringstream mErrorMsg
bool setValue(const std::string &name, const CType &value)
CVector< C_FLOAT64 > mVariables
virtual bool isValidProblem(const CCopasiProblem *pProblem)
unsigned C_INT32 * pUINT
size_t getNumIndependent() const
Definition: CState.cpp:342
CVector< C_FLOAT64 > mExponents
Definition: CLyapTask.h:59
CCopasiParameter * getParameter(const std::string &name)
const unsigned C_INT32 & getExponentNumber() const
CVector< C_FLOAT64 > mDWork
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
CMatrix< C_FLOAT64 > mJacobian
size_t size() const
Definition: CVector.h:100
void calculateDerivativesX(C_FLOAT64 *derivativesX)
Definition: CModel.cpp:1941
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
#define C_FLOAT64
Definition: copasi.h:92
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
CType * array()
Definition: CVector.h:139
bool divergenceRequested() const
CLyapTask * mpTask
void calculateDerivatives(C_FLOAT64 *derivatives)
Definition: CModel.cpp:1903
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
CVector< C_FLOAT64 > initializeAtolVector(const C_FLOAT64 &baseTolerance, const bool &reducedModel) const
Definition: CModel.cpp:4368
Definition: CModel.h:50
C_FLOAT64 mAverageDivergence
Definition: CLyapTask.h:72
const CState & getState() const
Definition: CModel.cpp:1771
virtual bool isValidProblem(const CCopasiProblem *pProblem)
virtual CType * array()
Definition: CMatrix.h:337
CModel * getModel() const
CCopasiContainer * getObjectParent() const
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
CVector< C_FLOAT64 > mSumExponents