COPASI API  4.16.103
CILDMMethod.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) 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include "copasi.h"
16 #include "CILDMMethod.h"
17 #include "CTSSATask.h"
18 #include "CTSSAProblem.h"
19 
20 #include "model/CReaction.h"
21 
24 #include "model/CModel.h"
25 #include "model/CMetab.h"
26 //#include "model/CState.h"
27 //#include "utilities/CMatrix.h"
28 //#include "utilities/CAnnotatedMatrix.h"
29 //#include "report/CCopasiObjectReference.h"
30 
31 #include "lapack/lapackwrap.h" // CLAPACK
32 #include "lapack/blaswrap.h" // BLAS
33 
34 //#define ILDMDEBUG
35 
37  CTSSAMethod(CCopasiMethod::tssILDM, pParent) //,
38  // mpState(NULL),
39  // mY(NULL)
40 {
41  //assert((void *) &mData == (void *) &mData.dim);
42 
43  // addObjectReference("Number of slow variables", mSlow, CCopasiObject::ValueInt);
44  // addMatrixReference("Contribution of Metabolites to Slow Space", mVslow, CCopasiObject::ValueDbl);
45 
46  mData.pMethod = this;
48 }
49 
51  const CCopasiContainer * pParent):
52  CTSSAMethod(src, pParent) //,
53  //mpState(NULL),
54  //mY(NULL)
55 {
56  //assert((void *) &mData == (void *) &mData.dim);
57 
58  mData.pMethod = this;
60 }
61 
63 {
65 }
66 
68 {
69  addObjectReference("Number of slow variables", mSlow, CCopasiObject::ValueInt);
70  addMatrixReference("Contribution of Species to Slow Space", mVslow, CCopasiObject::ValueDbl);
71 
73 
74  assertParameter("Deuflhard Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-6);
75  mReducedModel = true;
76 
77  //mData.dim = mpState->getNumIndependent();
78 
80  emptyVectors();
81 }
82 
83 void CILDMMethod::step(const double & deltaT)
84 {
85  C_INT failed_while = 0;
86 
87  C_INT dim = mData.dim;
88  C_INT fast = 0;
89  C_INT slow = dim - fast;
90 
91  C_INT i, j, k, info_schur = 0;
92  mY_initial.resize(dim);
93  mJacobian_initial.resize(dim, dim);
94  mQ.resize(dim, dim);
95  mR.resize(dim, dim);
96 
97  mQ_desc.resize(dim, dim);
98  mR_desc.resize(dim, dim);
99 
100  mTd.resize(dim, dim);
101  mTdInverse.resize(dim, dim);
102  mQz.resize(dim, dim);
103 
104  mTd_save.resize(dim, dim);
105  mTdInverse_save.resize(dim, dim);
106 
108  // TO REMOVE : mpModel->applyAssignments();
109  mpModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
110 
111  C_INT flag_jacob;
112  flag_jacob = 1; // Set flag_jacob=0 to printing Jacobian
113 
114  // To get the reduced Stoichiometry Matrix;
115 
116  CMatrix<C_FLOAT64> Stoichiom;
117  Stoichiom = mpModel -> getRedStoi();
118 
119  // const CCopasiVector< CReaction > & reacs = copasiModel->getReactions();
120  C_INT reacs_size = (C_INT) mpModel->getRedStoi().size();
121  reacs_size = reacs_size / dim; //TODO what is this?
122 
123  /* the vector mY is the current state of the system*/
124 
125  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
126  //C_FLOAT62 number2conc = 1.;
127 
128  //this is an ugly hack that only makes sense if all metabs are in the same compartment
129  //at the moment is is the only case the algorithm deals with
130 
131  CVector<C_FLOAT64> Xconc; //current state converted to concentrations
132  Xconc.resize(dim);
133 
134  for (i = 0; i < dim; ++i)
135  Xconc[i] = mY[i] * number2conc;
136 
137  for (i = 0; i < dim; i++)
138  mY_initial[i] = mY[i];
139 
140  CVector<C_FLOAT64> Xconc_initial; //current state converted to concentrations
141  Xconc_initial.resize(dim);
142 
143  for (i = 0; i < dim; ++i)
144  Xconc_initial[i] = mY_initial[i] * number2conc;
145 
146  // save initial Jacobian before next time step
147  for (i = 0; i < dim; i++)
148  for (j = 0; j < dim; j++)
149  mJacobian_initial(i, j) = mJacobian(i, j);
150 
151  // Next time step
152  integrationStep(deltaT);
153 
155  // TO REMOVE : mpModel->applyAssignments();
156 
157  // Calculate Jacobian for time step control
158  mpModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
159 
160 #ifdef ILDMDEBUG
161 
162  if (flag_jacob == 0)
163  {
164  std::cout << "Jacobian_next:" << std::endl;
165  std::cout << mJacobian << std::endl;
166  }
167 
168 #endif
169 
170 // to be removed
171  CMatrix<C_FLOAT64> Jacobian_for_test;
172  Jacobian_for_test.resize(dim, dim);
173 
174  for (i = 0; i < dim; i++)
175  for (j = 0; j < dim; j++)
176  Jacobian_for_test(i, j) = mJacobian_initial(i, j);
177 
178  // save initial Jacobian before next time step
179  for (i = 0; i < dim; i++)
180  for (j = 0; j < dim; j++)
181  mJacobian_initial(i, j) = mJacobian(i, j);
182 
183  schur(info_schur);
184 
185 #ifdef ILDMDEBUG
186  std::cout << "Eigenvalues of next Jacobian" << std::endl;
187 
188  for (i = 0; i < dim; i++)
189  std::cout << mR(i, i) << std::endl;
190 
191 #endif
192 
193  for (i = 0; i < dim; i++)
194  for (j = 0; j < dim; j++)
195  mJacobian_initial(i, j) = Jacobian_for_test(i, j);
196 
197  //end to be removed
198 
199  for (i = 0; i < dim; i++)
200  for (j = 0; j < dim; j++)
201  {
202  mTd_save(i, j) = 0;
203  mTdInverse_save(i, j) = 0;
204  }
205 
206  for (i = 0; i < dim; i++)
207  for (j = 0; j < dim; j++)
208  {
209  mTd(i, j) = 0;
210  mTdInverse(i, j) = 0;
211  }
212 
213  /** Schur Decomposition of Jacobian (reordered).
214  Output: mQ - transformation matrix mR - block upper triangular matrix (with ordered eigenvalues) */
215 
216  C_INT failed = 0;
217 // C_INT info_schur = 0;
218 
219  schur(info_schur); // TO DO : move the test to the TSSAMethod
220 
221  if (info_schur)
222  {
224  MCTSSAMethod + 5, mTime - deltaT);
225 
226  goto integration;
227  }
228 
229  C_INT flag_schur;
230  flag_schur = 0;
231 
232 #ifdef ILDMDEBUG
233 
234  std::cout << "Eigenvalues of initial Jacobian" << std::endl;
235 
236  for (i = 0; i < dim; i++)
237  std::cout << mR(i, i) << std::endl;
238 
239  if (flag_schur == 1)
240  {
241  std::cout << "Schur Decomposition" << std::endl;
242  std::cout << "mR - block upper triangular matrix :" << std::endl;
243  std::cout << mR << std::endl;
244  }
245 
246 #endif
247 
248  /* If complex eigenvalues */
249 
250  //BUG 873
251  if (mR(dim - 1, dim - 1) == mR(dim - 2 , dim - 2))
252  if (dim == 2)
253  {
254  slow = dim;
255  goto integration;
256  }
257 
258  // No reduction if the smallest eigenvalue is positive
259 
260  if (mR(dim - 1, dim - 1) >= 0)
261  {
262  slow = dim;
263  fast = 0;
264 
266  MCTSSAMethod + 6, mTime - deltaT);
267 
268  failed = 1;
269  goto integration;
270  }
271 
272  // C_INT number, k;
273  /** Classical ILDM iterations. The number of slow variables is decreased until the Deuflhard criterium holds */
274  /* do START slow iterations */
275 
276  while (slow > 1)
277  {
278 
279  fast = fast + 1;
280  slow = dim - fast;
281 
282  if (fast < dim - 1)
283  if (mR(slow, slow) == mR(slow - 1 , slow - 1))
284  fast = fast + 1;
285 
286  slow = dim - fast;
287 
288  for (i = 0; i < dim; i++)
289  for (j = 0; j < dim; j++)
290  {
291  mTd_save(i, j) = mTd(i, j);
292  mTdInverse_save(i, j) = mTdInverse(i, j);
293  }
294 
295  C_INT info = 0;
296  failed_while = 0;
297 
298  if (slow == 0)
299  {
300  for (i = 0; i < dim; i++)
301  for (j = 0; j < dim; j++)
302  {
303  mTdInverse(i, j) = mQ(j, i);
304  mTd(i, j) = mQ(i, j);
305  mQz(i, j) = mR(i, j);
306  }
307  }
308  else
309  {
310  /** Solution of Sylvester equation for given slow, mQ,mR
311  Output: mTd, mTdinverse and mQz (mQz is used later for newton iterations) */
312  sylvester(slow, info);
313 
314  if (info)
315  {
317  MCTSSAMethod + 7, slow, mTime - deltaT);
318 
319  failed_while = 1;
320  goto integration;
321  }
322  }
323 
324  /* Check real parts of eigenvalues of Jacobian */
325 
326  for (i = slow ; i < dim; i++)
327  if (mR(i , i) >= 0)
328  {
329  failed_while = 1;
330  goto integration;
331  }
332 
333  if (fast > 0)
334  mEPS = 1 / fabs(mR(slow , slow));
335 
336  mCfast.resize(fast);
337 
338  /** Deuflhard Iteration: Prove Deuflhard criteria, find consistent initial value for DAE
339  output: info - if Deuflhard is satisfied for this slow;
340  transformation matrices mTd and mTdinverse */
341 
342  info = 0;
343 
344  C_INT help;
345  help = 0;
346 
347  deuflhard(slow, info);
348  help = help + 1;
349 
350  failed_while = 0;
351 
352  if (info)
353  {
354  failed_while = 1;
355  goto integration;
356  }
357  }
358 
359  /** end of iterations to find the number of slow modes */
360 
361 integration:
362 
363  if ((failed == 1) || (failed_while == 1))
364  {
365  if (slow < dim)
366  {
367  fast = fast - 1;
368  slow = dim - fast;
369 
370  if ((fast >= 1) && (mR(slow - 1, slow - 1) == mR(slow , slow)))
371  fast = fast - 1;
372 
373  slow = dim - fast;
374 
375  for (i = 0; i < dim; i++)
376  for (j = 0; j < dim; j++)
377  {
378  mTd(i, j) = mTd_save(i, j);
379  mTdInverse(i, j) = mTdInverse_save(i, j);
380  }
381  }
382  }
383 
384  mSlow = slow;
385 
386  if (slow == dim)
388  MCTSSAMethod + 8, mTime);
389 
390  // test for orthogonality of the slow space
391 
392  C_INT flag_orthog = 1;
393 
394  if (flag_orthog == 0)
395  {
396  CMatrix<C_FLOAT64> orthog_prove;
397  orthog_prove.resize(dim, dim);
398 
399  for (i = 0; i < dim; i++)
400  for (j = 0; j < dim; j++)
401  orthog_prove(i, j) = orthog(i, j);
402  }
403 
404  C_INT flag_develop;
405  flag_develop = 0; //1;
406 
407  if (flag_develop == 0)
408  {
409  C_INT info_schur_desc = 0;
410 
411  /** Schur Decomposition of Jacobian (with another sorting: from fast to slow).
412  Output: mQ_desc - transformation matrix
413  mR_desc - block upper triangular matrix (with ordered eigenvalues) */
414 
415  schur_desc(info_schur_desc);
416 
417  for (i = 0; i < dim; i++)
418  for (j = 0; j < dim; j++)
419  {
420  mTd_save(i, j) = mTd(i, j);
421  mTdInverse_save(i, j) = mTdInverse(i, j);
422  }
423 
424  for (i = 0; i < dim; i++)
425  for (j = 0; j < dim; j++)
426  {
427  mTdInverse(i, j) = mQ_desc(j, i);
428  mTd(i, j) = mQ_desc(i, j);
429  }
430 
431  C_INT slow_desc;
432  slow_desc = fast;
433 
434  mat_anal_fast_space(slow_desc);
435  mat_anal_mod_space(slow_desc);
436 
437  CVector<C_FLOAT64> Reac_slow_space_orth;
438  Reac_slow_space_orth.resize(reacs_size);
439 
440  for (i = 0; i < reacs_size; i ++)
441  {
442  Reac_slow_space_orth[i] = 0;
443 
444  for (j = 0; j < dim; j++)
445  {
446  Reac_slow_space_orth[i] = Reac_slow_space_orth[i] + mVfast_space[j] * Stoichiom(j, i);
447  }
448  }
449 
450  for (i = 0; i < dim; i++)
451  for (j = 0; j < dim; j++)
452  {
453  mTd(i, j) = mTd_save(i, j);
454  mTdInverse(i, j) = mTdInverse_save(i, j);
455  }
456  }
457 
458  // end of test
459 
460  // Post Analysis
461 
462  mat_anal_mod(slow);
463  mat_anal_metab(slow);
464  mat_anal_mod_space(slow);
465  mat_anal_fast_space(slow);
466 
467  /** This part corresponds to the investigation of fast and slow reactions. In development at
468  the moment. To activate the calculations take flag_develop = 0;
469  */
470 
471  if (flag_develop == 0)
472  {
473  CMatrix<C_FLOAT64> Slow_react_contr;
474  Slow_react_contr.resize(dim, reacs_size);
475 
476  CMatrix<C_FLOAT64> Mode_react_contr;
477  Mode_react_contr.resize(dim, reacs_size);
478 
479  CVector<C_FLOAT64> Reac_slow_space;
480  Reac_slow_space.resize(reacs_size);
481 
482  CVector<C_FLOAT64> Reac_fast_space;
483  Reac_fast_space.resize(reacs_size);
484 
485  mReacSlowSpace.resize(reacs_size); //NEW TAB
486 
487  for (i = 0; i < dim; i ++)
488  for (j = 0; j < reacs_size; j++)
489  {
490  Slow_react_contr(i, j) = 0;
491 
492  for (k = 0; k < dim; k++)
493  Slow_react_contr(i, j) = Slow_react_contr(i, j) + mVslow(i, k) * Stoichiom(k, j);
494  }
495 
496  CVector<C_FLOAT64> denom_mode;
497  denom_mode.resize(dim);
498 
499  for (j = 0; j < dim; j++)
500  denom_mode[j] = 0;
501 
502  for (i = 0; i < dim; i++)
503  for (j = 0; j < reacs_size; j++)
504  denom_mode[i] = denom_mode[i] + fabs(Slow_react_contr(i, j));
505 
506  for (i = 0; i < dim; i++)
507  for (j = 0; j < reacs_size; j++)
508  {
509  if (denom_mode[i] == 0)
510  Mode_react_contr(i, j) = 0;
511  else
512  Mode_react_contr(i, j) = (Slow_react_contr(i, j)) / denom_mode[i] * 100;
513  }
514 
515  CVector<C_FLOAT64> denom_reac;
516  denom_reac.resize(reacs_size);
517 
518  for (j = 0; j < reacs_size; j++)
519  denom_reac[j] = 0;
520 
521  for (i = 0; i < reacs_size; i++)
522  for (j = 0; j < dim; j++)
523  denom_reac[i] = denom_reac[i] + fabs(Slow_react_contr(j, i));
524 
525  for (i = 0; i < reacs_size; i++)
526  for (j = 0; j < dim; j++)
527  {
528  if (denom_reac[i] == 0)
529  Slow_react_contr(j, i) = 0;
530  else
531  Slow_react_contr(j, i) = (Slow_react_contr(j , i)) / denom_reac[i] * 100;
532  }
533 
534  for (i = 0; i < reacs_size; i ++)
535  {
536  Reac_slow_space[i] = 0;
537 
538  for (j = 0; j < dim; j++)
539  {
540  Reac_slow_space[i] = Reac_slow_space[i] + mVslow_space[j] * Stoichiom(j, i);
541  }
542  }
543 
544  C_FLOAT64 length;
545  length = 0;
546 
547  for (i = 0; i < reacs_size; i++)
548  length = length + fabs(Reac_slow_space[i]);
549 
550  for (i = 0; i < reacs_size; i++)
551  if (length > 0)
552  mReacSlowSpace[i] = Reac_slow_space[i] / length * 100;
553  else
554  mReacSlowSpace[i] = 0;
555 
556  for (i = 0; i < reacs_size; i ++)
557  {
558  Reac_fast_space[i] = 0;
559 
560  for (j = 0; j < dim; j++)
561  Reac_fast_space[i] = Reac_fast_space[i] + mVfast_space[j] * Stoichiom(j, i);
562  }
563 
564  length = 0;
565 
566  for (i = 0; i < reacs_size; i++)
567  length = length + fabs(Reac_fast_space[i]);
568 
569  for (i = 0; i < reacs_size; i++)
570  if (length > 0)
571  Reac_fast_space[i] = Reac_fast_space[i] / length * 100;
572  else
573  Reac_fast_space[i] = 0;
574 
575 #ifdef ILDMDEBUG
576 
577  std::cout << "**********************************************************" << std::endl;
578  std::cout << "**********************************************************" << std::endl;
579  std::cout << std::endl;
580  std::cout << std::endl;
581 #endif
582 
583  mTMP1.resize(reacs_size, dim);
584  mTMP2.resize(reacs_size, dim);
585  mTMP3.resize(reacs_size, 1);
586 
587  for (i = 0; i < dim; i++)
588  for (j = 0; j < reacs_size; j++)
589  {
590  mTMP1(j, i) = Mode_react_contr(i, j);
591  mTMP2(j, i) = Slow_react_contr(i, j);
592  }
593 
594  for (j = 0; j < reacs_size; j++)
595  mTMP3(j, 0) = Reac_fast_space[j];
596  }
597 
598  // End of reaction analysis
599 
601  // TO REMOVE : mpModel->applyAssignments();
602 
603  // Calculate Jacobian for time step control
604 
605  mpModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
606 
607  // new entry for every entry contains the current data of currently step
608  setVectors(slow);
609 
610  // set the stepcounter
611  mCurrentStep += 1;
612 
613  return;
614 }
615 
616 /** Newton: Looking for consistent initial value for DAE system
617 Output: mCfast, info
618  */
619 
620 void CILDMMethod::newton(C_FLOAT64 *ys, C_INT & slow, C_INT & info)
621 {
622  C_INT i, j, iter, iterations, itermax;
623  C_INT nrhs, ok, fast;
624 
625  C_FLOAT64 tol, err;
626  C_INT dim = mData.dim;
627 
628  fast = dim - slow;
629 
630  CVector<C_INT> ipiv;
631  ipiv.resize(fast);
632 
633  CVector<C_FLOAT64> s_22_array;
634  s_22_array.resize(fast * fast);
635 
636  CVector<C_FLOAT64> gf_newton;
637  gf_newton.resize(fast);
638 
639  CVector<C_FLOAT64> d_yf;
640  d_yf.resize(fast);
641 
642  CVector<C_FLOAT64> y_newton;
643  y_newton.resize(dim);
644 
645  CVector<C_FLOAT64> yf_newton;
646  yf_newton.resize(fast);
647 
648  CVector<C_FLOAT64> x_newton;
649  x_newton.resize(dim);
650 
651  CVector<C_FLOAT64> dxdt_newton;
652  dxdt_newton.resize(dim);
653 
654  CVector<C_FLOAT64> g_newton;
655  g_newton.resize(dim);
656 
657  CMatrix<C_FLOAT64> S_22;
658  S_22.resize(fast, fast);
659 
660  C_FLOAT64 g1, g2 = 0;
661 
662  nrhs = 1;
663  //tol = 1e-9;
664  tol = 1e-9 / mpModel->getNumber2QuantityFactor();
665 
666  err = 10.0 / mpModel->getNumber2QuantityFactor();
667  iter = 0;
668 
669  itermax = 100;
670  iterations = 0;
671 
672  info = 0;
673 
674  for (i = 0; i < fast; i++)
675  for (j = 0; j < fast; j++)
676  S_22(i, j) = mQz(i, j);
677 
678  for (i = 0; i < fast; i++)
679  yf_newton[i] = mCfast[i];
680 
681  for (i = 0; i < fast; i++)
682  for (j = 0; j < fast; j++)
683  s_22_array[j + fast * i] = S_22(j, i);
684 
685  for (i = 0; i < fast; i++)
686  d_yf[i] = 0.;
687 
688  while (err > tol)
689  {
690  iter ++;
691 
692  if (iter > itermax)
693  {
694 
695  info = 1;
696  return;
697  }
698 
699  for (i = 0; i < fast; i++)
700  yf_newton[i] = yf_newton[i] + d_yf[i];
701 
702  /* back transformation */
703 
704  for (i = 0; i < slow; i++)
705  y_newton[i] = ys[i];
706 
707  for (i = slow; i < dim; i++)
708  y_newton[i] = yf_newton[i - slow];
709 
710  for (i = 0; i < dim; i++)
711  {
712  x_newton[i] = 0.0;
713 
714  for (j = 0; j < dim; j++)
715  x_newton[i] = x_newton[i] + mTd(i, j) * y_newton[j];
716  }
717 
718  calculateDerivativesX(x_newton.array(), dxdt_newton.array());
719 
720  for (i = 0; i < dim; i++)
721  {
722  g_newton[i] = 0.;
723 
724  for (j = 0; j < dim; j++)
725  g_newton[i] = g_newton[i] + mTdInverse(i, j) * dxdt_newton[j];
726  }
727 
728  // for (i = 0; i < fast; i++)
729  // gf_newton[i] = -1. * g_newton[i + slow];
730 
731  for (i = 0; i < fast; i++)
732  {
733  gf_newton[i] = -1. * g_newton[i + slow];
734  }
735 
736  /* int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
737  * *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info)
738  *
739  * -- LAPACK driver routine (version 3.0) --
740  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
741  * Courant Institute, Argonne National Lab, and Rice University
742  * March 31, 1993
743  *
744  *
745  * Purpose
746  * =======
747  *
748  * DGESV computes the solution to a real system of linear equations
749  * A * X = B,
750  * where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
751  *
752  * The LU decomposition with partial pivoting and row interchanges is
753  * used to factor A as
754  * A = P * L * U,
755  * where P is a permutation matrix, L is unit lower triangular, and U is
756  * upper triangular. The factored form of A is then used to solve the
757  * system of equations A * X = B.
758  *
759  * Arguments
760  * =========
761  *
762  * N (input) INTEGER
763  * The number of linear equations, i.e., the order of the
764  * matrix A. N >= 0.
765  *
766  * NRHS (input) INTEGER
767  * The number of right hand sides, i.e., the number of columns
768  * of the matrix B. NRHS >= 0.
769  *
770  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
771  * On entry, the N-by-N coefficient matrix A.
772  * On exit, the factors L and U from the factorization
773  * A = P*L*U; the unit diagonal elements of L are not stored.
774  *
775  * LDA (input) INTEGER
776  * The leading dimension of the array A. LDA >= max(1,N).
777  *
778  * IPIV (output) INTEGER array, dimension (N)
779  * The pivot indices that define the permutation matrix P;
780  * row i of the matrix was interchanged with row IPIV(i).
781  *
782  * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
783  * On entry, the N-by-NRHS matrix of right hand side matrix B.
784  * On exit, if INFO = 0, the N-by-NRHS solution matrix X.
785  *
786  * LDB (input) INTEGER
787  * The leading dimension of the array B. LDB >= max(1,N).
788  *
789  *
790  * INFO (output) INTEGER
791  * = 0: successful exit
792  * < 0: if INFO = -i, the i-th argument had an illegal value
793  * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
794  * has been completed, but the factor U is exactly
795  * singular, so the solution could not be computed.
796  */
797 
798  dgesv_(&fast, &nrhs, s_22_array.array(), &fast, ipiv.array(), gf_newton.array(), &fast, &ok);
799 
800  if (ok != 0)
801  {
802 
803  info = 2;
804  break;
805  }
806 
807  for (i = 0; i < fast; i++)
808  d_yf[i] = gf_newton[i];
809 
810  err = -10.;
811 
812  /* for (i = 0; i < fast; i++)
813  {
814  gf_newton[i] = fabs(gf_newton[i]);
815 
816  if (err < gf_newton[i])
817  err = gf_newton[i];
818  }
819 
820  */
821 
822  for (i = 0; i < fast; i++)
823  {
824  gf_newton[i] = fabs(gf_newton[i]);
825 
826  if (err < gf_newton[i])
827  err = gf_newton[i];
828  }
829 
830  iterations = iterations + 1;
831 
832  /* stop criterion of newton method */
833 
834  // C_FLOAT64 g1, g2;
835 
836  // g2 = err;
837 
838  if (iter == 1)
839  g1 = 3.0 * err;
840  else
841  g1 = g2;
842 
843  g2 = err;
844 
845  if (g2 / g1 > 1.0)
846  {
847 
848  info = 1;
849  break;
850  }
851  } /* end while */
852 
853  for (i = 0; i < fast; i++)
854  mCfast[i] = yf_newton[i];
855 
856  return;
857 }
858 
859 void CILDMMethod::start(const CState * initialState)
860 {
861  mReducedModel = true;
862 
863  integrationMethodStart(initialState);
864 
865  mDtol = * getValue("Deuflhard Tolerance").pUDOUBLE;
866 
867  /* ILDM related staff */
868 
869  //mDtol = mpProblem->getDeufelhardTol();
870 
875 
876  //createAnnotationsM();
877  emptyVectors();
878 
879  return;
880 }
881 
882 /**
883  Deuflhard Iteration: Prove Deuflhard criteria, find consistent initial value for DAE
884  output: info - if Deuflhard is satisfied for given slow; transformation matrices
885  mTd and mTdinverse
886  */
887 void CILDMMethod::deuflhard(C_INT & slow, C_INT & info)
888 {
889  C_INT i, j;
890  C_INT dim = mData.dim;
891  C_INT fast = dim - slow;
892  C_INT flag_deufl;
893 
894  flag_deufl = 1; // set flag_deufl = 0 to print the results of calculations
895 
896  /* calculations before relaxing yf to slow manifold */
897 
898  CVector<C_FLOAT64> c_full;
899  c_full.resize(dim);
900 
901  CVector<C_FLOAT64> c_slow;
902  c_slow.resize(slow);
903 
904  /* the vector mY is the current state of the system*/
905 
906  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
907  //C_FLOAT62 number2conc = 1.;
908 
909  //this is an ugly hack that only makes sense if all metabs are in the same compartment
910  //at the moment is is the only case the algorithm deals with
911 
912  CVector<C_FLOAT64> Xconc; //current state converted to concentrations
913  Xconc.resize(dim);
914 
915  for (i = 0; i < dim; ++i)
916  Xconc[i] = mY_initial[i] * number2conc;
917 
918  for (i = 0; i < dim; i++)
919  {
920  c_full[i] = 0.0;
921 
922  for (j = 0; j < dim; j++)
923  c_full[i] = c_full[i] + mTdInverse(i, j) * Xconc[j];
924  }
925 
926  for (j = 0; j < slow; j++)
927  c_slow[j] = c_full[j];
928 
929  for (j = 0; j < fast; j++)
930  mCfast[j] = c_full[j + slow];
931 
932  CVector<C_FLOAT64> g_full;
933  g_full.resize(dim);
934 
935  CVector<C_FLOAT64> g_slow;
936  g_slow.resize(slow);
937 
938  CVector<C_FLOAT64> g_fast;
939  g_fast.resize(fast);
940 
941  CVector<C_FLOAT64> dxdt;
942  dxdt.resize(dim);
943 
945 
946  for (j = 0; j < dim; j++)
947  dxdt[j] = 0.;
948 
949  CVector<C_FLOAT64> x_help;
950  x_help.resize(dim);
951 
952  for (j = 0; j < dim; j++)
953  {
954  x_help[j] = mY_initial[j] * number2conc;
955  }
956 
957  calculateDerivativesX(x_help.array(), dxdt.array());
958 
959  for (i = 0; i < dim; i++)
960  {
961  g_full[i] = 0.0;
962 
963  for (j = 0; j < dim; j++)
964  g_full[i] = g_full[i] + mTdInverse(i, j) * dxdt[j];
965  }
966 
967  for (j = 0; j < slow; j++)
968  g_slow[j] = g_full[j];
969 
970  info = 0;
971 
972  /** NEWTON: Looking for consistent initial value for DAE system
973  Output: mCfast, info */
974  newton(c_slow.array(), slow, info);
975 
976  if (info)
977  {
978  /* TODO */
979 
980  return;
981  }
982 
983  /* calculation of g_relax at point x_relax (after relaxing yf to slow manifold)*/
984 
985  CVector<C_FLOAT64> c_relax;
986  c_relax.resize(dim);
987 
988  CVector<C_FLOAT64> x_relax;
989  x_relax.resize(dim);
990 
991  CVector<C_FLOAT64> dxdt_relax;
992  dxdt_relax.resize(dim);
993 
994  CVector<C_FLOAT64> g_relax;
995  g_relax.resize(dim);
996 
997  for (i = 0; i < slow; i++)
998  c_relax[i] = c_slow[i];
999 
1000  for (i = slow; i < dim; i++)
1001  c_relax[i] = mCfast[i - slow];
1002 
1003  for (i = 0; i < dim; i++)
1004  {
1005  x_relax[i] = 0.0;
1006 
1007  for (j = 0; j < dim; j++)
1008  x_relax[i] = x_relax[i] + mTd(i, j) * c_relax[j];
1009  }
1010 
1011  calculateDerivativesX(x_relax.array(), dxdt_relax.array());
1012 
1013  for (i = 0; i < dim; i++)
1014  {
1015  g_relax[i] = 0.0;
1016 
1017  for (j = 0; j < dim; j++)
1018  g_relax[i] = g_relax[i] + mTdInverse(i, j) * dxdt_relax[j];
1019  }
1020 
1021  CVector<C_FLOAT64> re;
1022  re.resize(slow);
1023 
1024  /* stop criterion for slow reaction modes */
1025 
1026  for (i = 0; i < slow; i++)
1027  {
1028  re[i] = fabs(g_relax[i] - g_slow[i]);
1029  re[i] = re[i] * mEPS;
1030  }
1031 
1032  C_FLOAT64 max = 0.;
1033 
1034  for (i = 0; i < slow; i++)
1035  if (max < re[i])
1036  max = re[i];
1037 
1038  C_FLOAT64 max1;
1039  C_FLOAT64 norm = 0;
1040 
1041  for (i = 0; i < slow; i++)
1042  norm = norm + fabs(g_relax[i] - g_slow[i]);
1043 
1044  max1 = norm * mEPS;
1045 
1046  if (max >= mDtol / mpModel->getNumber2QuantityFactor())
1047  info = 1;
1048  else
1049  info = 0;
1050 
1051  return;
1052 }
1053 
1054 /**
1055  * Empty every vector to be able to fill them with new values for a new calculation.
1056  * Also nullify the step counter.
1057  **/
1059 {
1060  mCurrentStep = 0;
1061  mVec_mVslow.erase(mVec_mVslow.begin(), mVec_mVslow.end());
1062  mVec_TimeScale.erase(mVec_TimeScale.begin(), mVec_TimeScale.end());
1063  mVec_mVslowMetab.erase(mVec_mVslowMetab.begin(), mVec_mVslowMetab.end());
1064  mVec_mVslowSpace.erase(mVec_mVslowSpace.begin(), mVec_mVslowSpace.end());
1065  mVec_SlowModes.erase(mVec_SlowModes.begin(), mVec_SlowModes.end());
1066 
1067  /* temporary tabs */
1068 
1069  mVec_mTMP1.erase(mVec_mTMP1.begin(), mVec_mTMP1.end());
1070  mVec_mTMP2.erase(mVec_mTMP2.begin(), mVec_mTMP2.end());
1071  mVec_mTMP3.erase(mVec_mTMP3.begin(), mVec_mTMP3.end());
1072 }
1073 
1074 /**
1075  *upgrade all vectors with values from actually calculalion for current step
1076  **/
1077 void CILDMMethod::setVectors(int slowMode)
1078 {
1079  mVec_mVslow.push_back(mCurrentStep);
1082 
1083  mVec_TimeScale.push_back(mCurrentStep);
1085  size_t i;
1086 
1087  for (i = 0; i < (size_t) mData.dim; i++)
1088  mVec_TimeScale[mCurrentStep][i] = -1 / mR(i, i);
1089 
1090  mVec_mVslowMetab.push_back(mCurrentStep);
1093 
1094  mVec_mVslowSpace.push_back(mCurrentStep);
1096 
1097  for (i = 0; i < mVslow_space.size(); i++)
1098  {
1100  }
1101 
1102  mVec_mVfastSpace.push_back(mCurrentStep);
1105 
1106  mVec_SlowModes.push_back(mCurrentStep);
1107  mVec_SlowModes[mCurrentStep] = slowMode;
1108 
1109  mCurrentTime.push_back(mCurrentStep);
1111 
1112  // NEW TAB
1113 
1114  mVec_mReacSlowSpace.push_back(mCurrentStep);
1117 
1118  /* temporary tabs */
1119 
1120  size_t reacs_size = mpModel->getReactions().size();
1121 
1122  mVec_mTMP1.push_back(mCurrentStep);
1123  mVec_mTMP1[mCurrentStep].resize(reacs_size, mData.dim);
1125 
1126  mVec_mTMP2.push_back(mCurrentStep);
1127  mVec_mTMP2[mCurrentStep].resize(reacs_size, mData.dim);
1129 
1130  mVec_mTMP3.push_back(mCurrentStep);
1131  mVec_mTMP3[mCurrentStep].resize(reacs_size, 1);
1133 }
1134 /**
1135  * Create the CArraAnnotations for every ILDM-tab in the CQTSSAResultSubWidget.
1136  * Input for each CArraAnnotations is a seperate CMatrix.
1137  **/
1139 {
1140  tableNames.erase(tableNames.begin(), tableNames.end());
1141 
1142  std::string name;
1143 
1144  name = "Contribution of species to modes";
1145  tableNames.push_back(name);
1146 
1148  pTmp1 = new CArrayAnnotation("Contribution of species to modes", this,
1150  pTmp1->setMode(0, pTmp1->STRINGS);
1151  pTmp1->setMode(1, pTmp1->VECTOR);
1152  pTmp1->setDescription(" ");
1153  pTmp1->setDimensionDescription(0, "Contribution to mode (TS - corresponding timescale)");
1154  pTmp1->setDimensionDescription(1, "Species");
1156 
1158 
1159  name = "Modes distribution for species";
1160  tableNames.push_back(name);
1161 
1163  pTmp2 = new CArrayAnnotation("Modes distribution for species", this,
1165  pTmp2->setMode(1, pTmp2->STRINGS);
1166  pTmp2->setMode(0, pTmp2->VECTOR);
1167  pTmp2->setDescription(" ");
1168  pTmp2->setDimensionDescription(0, "Mode distribution for each species");
1169  pTmp2->setDimensionDescription(1, "Modes (TS - corresponding timescale)");
1171 
1173 
1174  name = "Slow space";
1175  tableNames.push_back(name);
1176 
1178  pTmp3 = new CArrayAnnotation("Slow space", this,
1180  pTmp3->setMode(1, pTmp3->STRINGS);
1181  pTmp3->setMode(0, pTmp3->VECTOR);
1182  pTmp3->setDescription(" ");
1183  pTmp3->setDimensionDescription(0, "Species");
1184  pTmp3->setDimensionDescription(1, "Contribution to slow space");
1186 
1188 
1189  name = "Fast space";
1190  tableNames.push_back(name);
1191 
1193  pTmp4 = new CArrayAnnotation("Fast space", this,
1195  pTmp4->setMode(1, pTmp4->STRINGS);
1196  pTmp4->setMode(0, pTmp4->VECTOR);
1197  pTmp4->setDescription(" ");
1198  pTmp4->setDimensionDescription(0, "Species");
1199  pTmp4->setDimensionDescription(1, "Contribution to fast space");
1201 
1203 
1204  // NEW TAB
1205 
1206  name = "Reactions slow space";
1207  tableNames.push_back(name);
1208 
1210  pTmp5 = new CArrayAnnotation("Reactions slow space", this,
1212  pTmp5->setMode(1, pTmp5->STRINGS);
1213  pTmp5->setMode(0, pTmp5->VECTOR);
1214  pTmp5->setDescription(" ");
1215  pTmp5->setDimensionDescription(0, "Reactions");
1216  pTmp5->setDimensionDescription(1, "Contribution to slow space");
1218 
1220 
1221  /* tamporary tabs */
1222 
1223  name = "Reactions contribution to the mode";
1224  tableNames.push_back(name);
1225 
1227  pTMP1 = new CArrayAnnotation("", this,
1229  pTMP1->setMode(0, pTMP1->VECTOR);
1230  pTMP1->setMode(1, pTMP1->STRINGS);
1231  pTMP1->setDescription("Reactions contribution to the mode ");
1232  pTMP1->setDimensionDescription(0, "Reactions");
1233  pTMP1->setDimensionDescription(1, "Modes (TS - corresponding timescale)");
1234  pTMP1PrintAnn = pTMP1;
1235 
1236  mapTableToName[name] = pTMP1PrintAnn;
1237 
1238  name = "Reactions distribution between modes ";
1239  tableNames.push_back(name);
1240 
1242  pTMP2 = new CArrayAnnotation("", this,
1244  pTMP2->setMode(0, pTMP2->VECTOR);
1245  pTMP2->setMode(1, pTMP2->STRINGS);
1246  pTMP2->setDescription("Reactions distribution between modes ");
1247  pTMP2->setDimensionDescription(0, "Reactions");
1248  pTMP2->setDimensionDescription(1, "Modes (TS - corresponding timescale)");
1249  pTMP2PrintAnn = pTMP2;
1250 
1251  mapTableToName[name] = pTMP2PrintAnn;
1252 
1253  name = "Reactions fast space";
1254  tableNames.push_back(name);
1255 
1257  pTMP3 = new CArrayAnnotation("", this,
1259  pTMP3->setMode(0, pTMP3->VECTOR);
1260  pTMP3->setMode(1, pTMP3->STRINGS);
1261  pTMP3->setDescription("Reactions fast space ");
1262  pTMP3->setDimensionDescription(0, "Reactions");
1263  pTMP3->setDimensionDescription(1, "");
1264  pTMP3PrintAnn = pTMP3;
1265 
1266  mapTableToName[name] = pTMP3PrintAnn;
1267 }
1268 /**
1269  * Set the every CArrayAnnotation for the requested step.
1270  * Set also the desription of CArayAnnotation for both dimensions:
1271  * - dimension description could consists of some std::srings
1272  * some strings contain the Time Scale values for requested step
1273  * - dimension description could consists of arrays of CommonNames
1274  **/
1275 //void CILDMMethod::setAnnotationM(int step)
1277 {
1278  if (step == 0) return false;
1279 
1280  if (mVec_mVslow.size() == 0) return false;
1281 
1282  if (step > mVec_mVslow.size()) return false;
1283 
1284  if (step > mVec_SlowModes.size()) return false;
1285 
1286  step -= 1;
1287  double timeScale;
1288  std::string str;
1289  std::stringstream sstr;
1290  sstr.str("");
1291  sstr.clear();
1292  C_INT i;
1293 
1298 
1299  for (i = 0; i < mData.dim; i++)
1300  {
1301  timeScale = mVec_TimeScale[step][i];
1302 
1303  if (i < mVec_SlowModes[step])
1304  sstr << "Slow: ";
1305  else
1306  sstr << "Fast: ";
1307 
1308  sstr << timeScale;
1309  str = sstr.str();
1310  pVslowPrintAnn->setAnnotationString(0, i, str);
1311  sstr.str("");
1312  sstr.clear();
1313  }
1314 
1319 
1320  for (i = 0; i < mData.dim; i++)
1321  {
1322  timeScale = mVec_TimeScale[step][i];
1323 
1324  if (i < mVec_SlowModes[step])
1325  sstr << "Slow: ";
1326  else
1327  sstr << "Fast: ";
1328 
1329  sstr << timeScale;
1330  str = sstr.str();
1332  sstr.str("");
1333  sstr.clear();
1334  }
1335 
1336  sstr << mVec_SlowModes[step];
1337  // if (mVec_SlowModes[step] > 1)
1338  // sstr << " slow modes";
1339  //else
1340  // sstr << " slow mode";
1341  sstr << " slow; ";
1342 
1343  C_INT dim = mData.dim;
1344  sstr << dim - mVec_SlowModes[step];
1345  sstr << " fast";
1346 
1347  str = sstr.str();
1349 
1350  for (i = 0; i < mData.dim; i++)
1352 
1356 
1358 
1359  for (i = 0; i < mData.dim; i++)
1361 
1365 
1366  //sstr.clear();
1367  sstr.str("");
1368  str = sstr.str();
1369 
1371 
1372  for (i = 0; i < (C_INT) mReacSlowSpace.size(); i++)
1374 
1378 
1379  sstr.str("");
1380  str = sstr.str();
1381 
1382  /* temporary tabs */
1383 
1384  size_t reacs_size = mpModel->getReactions().size();
1385 
1386  mTMP1Print.resize(reacs_size, mData.dim);
1388 
1389  pTMP1PrintAnn->resize();
1390 
1391  for (i = 0; i < (C_INT) mVec_mTMP1[step].numCols(); i++)
1392  {
1393  timeScale = mVec_TimeScale[step][i];
1394 
1395  if (i < mVec_SlowModes[step])
1396  sstr << "Slow: ";
1397  else
1398  sstr << "Fast: ";
1399 
1400  sstr << timeScale;
1401  str = sstr.str();
1402  pTMP1PrintAnn->setAnnotationString(1, i, str);
1403  sstr.str("");
1404  sstr.clear();
1405  }
1406 
1408 
1409  mTMP2Print.resize(reacs_size, mData.dim);
1411 
1412  pTMP2PrintAnn->resize();
1413 
1414  for (i = 0; i < (C_INT) mVec_mTMP2[step].numCols(); i++)
1415  {
1416  timeScale = mVec_TimeScale[step][i];
1417 
1418  if (i < mVec_SlowModes[step])
1419  sstr << "Slow: ";
1420  else
1421  sstr << "Fast: ";
1422 
1423  sstr << timeScale;
1424  str = sstr.str();
1425  pTMP2PrintAnn->setAnnotationString(1, i, str);
1426  sstr.str("");
1427  sstr.clear();
1428  }
1429 
1431 
1432  mTMP3Print.resize(reacs_size, 1);
1434  pTMP3PrintAnn->resize();
1436  pTMP3PrintAnn->setAnnotationString(1, 0, str);
1437 
1438  return true;
1439 }
1440 
1441 void CILDMMethod::printResult(std::ostream * ostream) const // temporary tabs are not included in Report !!!
1442 {
1443  std::ostream & os = *ostream;
1444  double timeScale;
1445  C_INT i, j, istep = 0;
1446 
1447  C_INT32 stepNumber;
1448 
1449  this->print(&os);
1450 
1451  assert(getObjectDataModel() != NULL);
1452 
1453  //stepNumber = pProblem->getStepNumber();
1454  stepNumber = mVec_SlowModes.size();
1455 
1456  for (istep = 0; istep < stepNumber; istep++)
1457  {
1458 
1459  os << std::endl;
1460  os << "**************** Time step " << istep + 1 << " ************************** " << std::endl;
1461 
1462  os << std::endl;
1463 
1464  os << "Contribution of species to modes" << std::endl;
1465 
1466  os << "Rows : contribution to mode (TS - corresponding timescale)" << std::endl;
1467  os << "Columns: species ";
1468 
1469  for (j = 0; j < mData.dim; j++)
1470  {
1471  os << mpModel->getMetabolitesX()[j]->getObjectName() << " ";
1472  }
1473 
1474  os << std::endl;
1475 
1476  for (i = 0; i < mData.dim; i++)
1477  {
1478  timeScale = mVec_TimeScale[istep][i];
1479 
1480  if (i < mVec_SlowModes[istep])
1481  os << "Slow (";
1482  else
1483  os << "Fast (";
1484 
1485  os << timeScale << "): ";
1486 
1487  for (j = 0; j < mData.dim; j++)
1488  os << mVec_mVslow[istep][i][j] << " ";
1489 
1490  os << std::endl;
1491  }
1492 
1493  os << std::endl;
1494 
1495  os << "Modes distribution for species" << std::endl;
1496 
1497  os << "Rows : Mode distribution for each species" << std::endl;
1498  os << "Columns: Modes (TS - corresponding timescale) ";
1499  os << std::endl;
1500 
1501  for (i = 0; i < mData.dim; i++)
1502  {
1503  timeScale = mVec_TimeScale[istep][i];
1504 
1505  if (i < mVec_SlowModes[istep])
1506  os << "Slow (";
1507  else
1508  os << "Fast (";
1509 
1510  os << timeScale << ") ";
1511  }
1512 
1513  os << std::endl;
1514 
1515  for (j = 0; j < mData.dim; j++)
1516  {
1517  os << mpModel->getMetabolitesX()[j]->getObjectName() << " ";
1518 
1519  for (i = 0; i < mData.dim; i++)
1520  os << mVec_mVslowMetab[istep][j][i] << " ";
1521 
1522  os << std::endl;
1523  }
1524 
1525  os << std::endl;
1526 
1527  os << "Slow space" << std::endl;
1528 
1529  os << "Rows : Species" << std::endl;
1530  os << "Column: Contribution to slow space ";
1531  os << std::endl;
1532 
1533  os << mVec_SlowModes[istep];
1534  os << " slow; ";
1535 
1536  os << mData.dim - mVec_SlowModes[istep];
1537  os << " fast";
1538  os << std::endl;
1539 
1540  for (j = 0; j < mData.dim; j++)
1541  {
1542  os << mpModel->getMetabolitesX()[j]->getObjectName() << " ";
1543  os << mVec_mVslowSpace[istep][j] << " ";
1544 
1545  os << std::endl;
1546  }
1547 
1548  os << std::endl;
1549  os << "Fast space" << std::endl;
1550 
1551  os << "Rows : Species" << std::endl;
1552  os << "Column: Contribution to fast space ";
1553  os << std::endl;
1554 
1555  os << mVec_SlowModes[istep];
1556  os << " slow; ";
1557 
1558  os << mData.dim - mVec_SlowModes[istep];
1559  os << " fast";
1560  os << std::endl;
1561 
1562  for (j = 0; j < mData.dim; j++)
1563  {
1564  os << mpModel->getMetabolitesX()[j]->getObjectName() << " ";
1565  os << mVec_mVfastSpace[istep][j] << " ";
1566 
1567  os << std::endl;
1568  }
1569 
1570  os << std::endl;
1571  os << "Reactions slow space" << std::endl;
1572 
1573  os << "Rows : Reactions" << std::endl;
1574  os << "Column: Contribution to slow space ";
1575  os << std::endl;
1576 
1577  for (j = 0; j < (C_INT32) mpModel->getReactions().size(); j++)
1578  {
1579  os << mpModel->getReactions()[j]->getObjectName() << " ";
1580  os << mVec_mReacSlowSpace[istep][j] << " ";
1581 
1582  os << std::endl;
1583  }
1584 
1585  os << std::endl;
1586  }
1587 
1588  return;
1589 }
CArrayAnnotation * pTMP3
Definition: CILDMMethod.h:160
void schur(C_INT &info)
void deuflhard(C_INT &slow, C_INT &info)
C_FLOAT64 mEPS
Definition: CTSSAMethod.h:346
#define C_INT
Definition: copasi.h:115
virtual void start(const CState *initialState)
int dgesv_(integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info)
void sylvester(C_INT slow, C_INT &info)
#define MCTSSAMethod
#define pdelete(p)
Definition: copasi.h:215
CMatrix< C_FLOAT64 > mVslowPrint
Definition: CILDMMethod.h:166
void emptyVectors()
std::vector< C_FLOAT64 > mCurrentTime
Definition: CTSSAMethod.h:460
CMatrix< C_FLOAT64 > mVslowSpacePrint
Definition: CILDMMethod.h:167
std::vector< CVector< C_FLOAT64 > > mVec_mVfastSpace
Definition: CILDMMethod.h:117
void newton(C_FLOAT64 *ys, C_INT &slow, C_INT &info)
CMatrix< C_FLOAT64 > mVslow
Definition: CTSSAMethod.h:261
virtual size_t size() const
CMatrix< C_FLOAT64 > mTdInverse
Definition: CTSSAMethod.h:230
CVector< C_FLOAT64 > mY_initial
Definition: CTSSAMethod.h:195
CMatrix< C_FLOAT64 > mTMP3
Definition: CILDMMethod.h:112
std::vector< CMatrix< C_FLOAT64 > > mVec_mVslow
Definition: CILDMMethod.h:114
CState * mpState
Definition: CTSSAMethod.h:175
CArrayAnnotation * pTMP2
Definition: CILDMMethod.h:159
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
CArrayAnnotation * pVslowSpacePrintAnn
Definition: CILDMMethod.h:137
void calculateJacobianX(CMatrix< C_FLOAT64 > &jacobianX, const C_FLOAT64 &derivationFactor, const C_FLOAT64 &resolution)
Definition: CModel.cpp:2082
CModel * mpModel
Definition: CTSSAMethod.h:336
std::map< std::string, CArrayAnnotation * > mapTableToName
Definition: CTSSAMethod.h:97
CVector< C_FLOAT64 > mVslow_space
Definition: CTSSAMethod.h:271
CMatrix< C_FLOAT64 > mVfastSpacePrint
Definition: CILDMMethod.h:168
void initializeIntegrationsParameter()
Definition: CState.h:305
CMatrix< C_FLOAT64 > mTd_save
Definition: CTSSAMethod.h:240
void integrationStep(const double &deltaT)
C_FLOAT64 mDtol
Definition: CTSSAMethod.h:341
CMatrix< C_FLOAT64 > mVslowMetabPrint
Definition: CILDMMethod.h:169
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
void schur_desc(C_INT &info)
std::vector< CMatrix< C_FLOAT64 > > mVec_mTMP2
Definition: CILDMMethod.h:122
void setDescription(const std::string &s)
virtual bool setAnnotationM(size_t step)
#define C_INT32
Definition: copasi.h:90
CArrayAnnotation * pVfastSpacePrintAnn
Definition: CILDMMethod.h:138
double orthog(C_INT &number1, C_INT &number2)
CMatrix< C_FLOAT64 > mQ
Definition: CTSSAMethod.h:215
const CMatrix< C_FLOAT64 > & getRedStoi() const
Definition: CModel.cpp:1154
CArrayAnnotation * pTmp1
Definition: CILDMMethod.h:150
CMatrix< C_FLOAT64 > mTdInverse_save
Definition: CTSSAMethod.h:244
void mat_anal_mod_space(C_INT &slow)
virtual void step(const double &deltaT)
Definition: CILDMMethod.cpp:83
CArrayAnnotation * pTMP3PrintAnn
Definition: CILDMMethod.h:145
CVector< C_FLOAT64 > mVfast_space
Definition: CTSSAMethod.h:276
void setAnnotationString(size_t d, size_t i, const std::string s)
std::vector< CVector< C_FLOAT64 > > mVec_mReacSlowSpace
Definition: CILDMMethod.h:118
CArrayAnnotation * pTMP2PrintAnn
Definition: CILDMMethod.h:144
CMatrix< C_FLOAT64 > mQ_desc
Definition: CTSSAMethod.h:216
CArrayAnnotation * pTMP1
Definition: CILDMMethod.h:158
CMatrix< C_FLOAT64 > mJacobian_initial
Definition: CTSSAMethod.h:210
const C_FLOAT64 & getNumber2QuantityFactor() const
Definition: CModel.cpp:2357
void calculateDerivativesX(C_FLOAT64 *X1, C_FLOAT64 *Y1)
CMatrix< C_FLOAT64 > mR_desc
Definition: CTSSAMethod.h:221
std::vector< CMatrix< C_FLOAT64 > > mVec_mTMP1
Definition: CILDMMethod.h:121
std::vector< std::string > tableNames
Definition: CTSSAMethod.h:98
void mat_anal_fast_space(C_INT &slow)
int mCurrentStep
Definition: CTSSAMethod.h:466
CVector< C_FLOAT64 > mReacSlowSpace
Definition: CILDMMethod.h:106
const Value & getValue() const
bool mReducedModel
Definition: CTSSAMethod.h:291
std::vector< CVector< C_FLOAT64 > > mVec_mVslowSpace
Definition: CILDMMethod.h:116
CArrayAnnotation * pReacSlowSpacePrintAnn
Definition: CILDMMethod.h:139
void setCopasiVector(size_t d, const CCopasiContainer *v)
CMatrix< C_FLOAT64 > mQz
Definition: CTSSAMethod.h:235
CArrayAnnotation * pTmp2
Definition: CILDMMethod.h:151
CILDMMethod(const CCopasiContainer *pParent=NULL)
Definition: CILDMMethod.cpp:36
CArrayAnnotation * pVslowPrintAnn
Definition: CILDMMethod.h:135
void createAnnotationsM()
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
void mat_anal_metab(C_INT &slow)
size_t size() const
Definition: CVector.h:100
void printResult(std::ostream *ostream) const
virtual void initializeParameter()
Definition: CILDMMethod.cpp:67
void setMode(size_t d, Mode m)
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
void setDimensionDescription(size_t d, const std::string &s)
CMatrix< C_FLOAT64 > mR
Definition: CTSSAMethod.h:220
CMatrix< C_FLOAT64 > mTMP1
Definition: CILDMMethod.h:110
std::vector< CVector< C_FLOAT64 > > mVec_TimeScale
Definition: CTSSAMethod.h:461
CMatrix< C_FLOAT64 > mJacobian
Definition: CTSSAMethod.h:205
CArrayAnnotation * pVslowMetabPrintAnn
Definition: CILDMMethod.h:136
CMatrix< C_FLOAT64 > mTMP2
Definition: CILDMMethod.h:111
#define C_FLOAT64
Definition: copasi.h:92
std::vector< CMatrix< C_FLOAT64 > > mVec_mTMP3
Definition: CILDMMethod.h:123
CType * array()
Definition: CVector.h:139
virtual size_t size() const
Definition: CMatrix.h:132
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
CMatrix< C_FLOAT64 > mTd
Definition: CTSSAMethod.h:225
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
CMatrix< C_FLOAT64 > mTMP2Print
Definition: CILDMMethod.h:175
CMatrix< C_FLOAT64 > mReacSlowSpacePrint
Definition: CILDMMethod.h:170
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
CArrayAnnotation * pTMP1PrintAnn
Definition: CILDMMethod.h:143
C_FLOAT64 mTime
Definition: CTSSAMethod.h:200
CCopasiObject * addMatrixReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
CArrayAnnotation * pTmp3
Definition: CILDMMethod.h:152
void integrationMethodStart(const CState *initialState)
std::vector< C_INT > mVec_SlowModes
Definition: CTSSAMethod.h:459
CTSSAMethod * pMethod
Definition: CTSSAMethod.h:169
void mat_anal_mod(C_INT &slow)
CArrayAnnotation * pTmp5
Definition: CILDMMethod.h:154
CMatrix< C_FLOAT64 > mVslow_metab
Definition: CTSSAMethod.h:266
void setVectors(int slowMode)
CArrayAnnotation * pTmp4
Definition: CILDMMethod.h:153
CVector< C_FLOAT64 > mCfast
Definition: CTSSAMethod.h:250
C_FLOAT64 * mY
Definition: CTSSAMethod.h:185
CMatrix< C_FLOAT64 > mTMP3Print
Definition: CILDMMethod.h:176
std::vector< CMatrix< C_FLOAT64 > > mVec_mVslowMetab
Definition: CILDMMethod.h:115
CMatrix< C_FLOAT64 > mTMP1Print
Definition: CILDMMethod.h:174
#define max(a, b)
Definition: f2c.h:176