COPASI API  4.16.103
CILDMModifiedMethod.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 #include "copasi.h"
12 
13 #include "CILDMModifiedMethod.h"
14 #include "CTSSAProblem.h"
15 #include "CTSSATask.h"
16 
19 #include "model/CModel.h"
20 #include "model/CMetab.h"
21 //#include "model/CState.h"
22 //#include "utilities/CMatrix.h"
23 //#include "utilities/CAnnotatedMatrix.h"
24 //#include "report/CCopasiObjectReference.h"
25 
26 #include "lapack/lapackwrap.h" // CLAPACK
27 #include "lapack/blaswrap.h" // BLAS
28 
30  CTSSAMethod(CCopasiMethod::tssILDMModified, pParent) //,
31  // mpState(NULL),
32  // mY(NULL)
33 {
34  // assert((void *) &mData == (void *) &mData.dim);
35 
36  // addObjectReference("Number of slow variables", mSlow, CCopasiObject::ValueInt);
37  // addMatrixReference("Contribution of Metabolites to Slow Space", mVslow, CCopasiObject::ValueDbl);
38 
39  mData.pMethod = this;
41 }
42 
44  const CCopasiContainer * pParent):
45  CTSSAMethod(src, pParent) //,
46  //mpState(NULL),
47  //mY(NULL)
48 {
49  // assert((void *) &mData == (void *) &mData.dim);
50 
51  mData.pMethod = this;
53 }
54 
56 {
58 }
59 
61 {
62  addObjectReference("Number of slow variables", mSlow, CCopasiObject::ValueInt);
63  addMatrixReference("Contribution of Species to Slow Space", mVslow, CCopasiObject::ValueDbl);
64 
66 
67  assertParameter("Deuflhard Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-6);
68  mReducedModel = true;
69 
70  //mData.dim = mpState->getNumIndependent();
71 
73  emptyVectors();
74 }
75 
76 void CILDMModifiedMethod::start(const CState * initialState)
77 {
78  mReducedModel = true;
79 
80  integrationMethodStart(initialState);
81 
82  /* ILDM related staff */
83 
84  mDtol = * getValue("Deuflhard Tolerance").pUDOUBLE;
85 
90 
91  emptyVectors();
92 
93  return;
94 }
95 
96 void CILDMModifiedMethod::step(const double & deltaT)
97 {
98 
99  C_INT dim = mData.dim;
100  C_INT fast = 0;
101  C_INT slow = dim - fast;
102 
103  C_INT slow2, fast2;
104 
105  slow2 = dim;
106  fast2 = dim - slow2;
107 
108  C_INT i, j;
109 
110  mY_initial.resize(dim);
111  mJacobian_initial.resize(dim, dim);
112  mQ.resize(dim, dim);
113  mR.resize(dim, dim);
114  mTd.resize(dim, dim);
115  mTdInverse.resize(dim, dim);
116  mQz.resize(dim, dim);
117 
118  mTd_save.resize(dim, dim);
119 
120  mTdInverse_save.resize(dim, dim);
121 
123  // TO REMOVE : mpModel->applyAssignments();
124  mpModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
125 
126  C_INT flag_jacob;
127  flag_jacob = 1; // Set flag_jacob=0 to print Jacobian
128 
129  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
130  //C_FLOAT64 number2conc = 1.;
131 
132  //this is an ugly hack that only makes sense if all metabs are in the same compartment
133  //at the moment is is the only case the algorithm deals with
134 
135  CVector<C_FLOAT64> Xconc; //current state converted to concentrations
136  Xconc.resize(dim);
137 
138  for (i = 0; i < dim; ++i)
139  Xconc[i] = mY[i] * number2conc;
140 
141  for (i = 0; i < dim; i++)
142  mY_initial[i] = mY[i];
143 
144  CVector<C_FLOAT64> Xconc_initial; //current state converted to concentrations
145  Xconc_initial.resize(dim);
146 
147  for (i = 0; i < dim; ++i)
148  Xconc_initial[i] = mY_initial[i] * number2conc;
149 
150  // save initial Jacobian before next time step
151  for (i = 0; i < dim; i++)
152  for (j = 0; j < dim; j++)
153  mJacobian_initial(i, j) = mJacobian(i, j);
154 
155  // Next time step
156  integrationStep(deltaT);
157 
159  // TO REMOVE : mpModel->applyAssignments();
160 
161  // Calculate Jacobian for time step control
162  mpModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
163 
164  //CMatrix<C_FLOAT64> mTd_save;
165  for (i = 0; i < dim; i++)
166  for (j = 0; j < dim; j++)
167  {
168  mTd_save(i, j) = 0;
169  mTdInverse_save(i, j) = 0;
170  }
171 
172  for (i = 0; i < dim; i++)
173  for (j = 0; j < dim; j++)
174  {
175  mTd(i, j) = 0;
176  mTdInverse(i, j) = 0;
177  }
178 
179  /** Schur Decomposition of Jacobian (reordered).
180  Output: mQ - transformation matrix mR - block upper triangular matrix (with ordered eigenvalues) */
181 
182  C_INT failed = 0;
183  C_INT info_schur = 0;
184 
185  C_INT number, k;
186  C_INT failed_while = 0;
187 
188  C_INT flag_deufl;
189  flag_deufl = 1;
190 
191  C_FLOAT64 max = 0.;
192 
193  //C_FLOAT64 max = 0;
195  CVector<C_FLOAT64> dxdt_relax;
196  CVector<C_FLOAT64> x_relax;
197  CVector<C_FLOAT64> x_help;
198  CVector<C_FLOAT64> dxdt;
199 
200  CVector<C_FLOAT64> x_zero;
201 
202  CVector<C_FLOAT64> dxdt_zero;
203 
204  CVector<C_FLOAT64> dxdt_real;
205 
206  CVector<C_FLOAT64> help;
207 
208  CVector<C_INT> index;
209  CVector<C_INT> index_temp;
210  CMatrix<C_FLOAT64> orthog_prove;
211  orthog_prove.resize(dim, dim);
212 
213  C_INT info;
214 
215  CVector<C_INT> index_metab;
216  index_metab.resize(dim);
217 
218  /** Schur transformation of Jacobian */
219  schur(info_schur);
220 
221  if (info_schur)
222  {
224  MCTSSAMethod + 9, mTime - deltaT);
225 
226  goto integration;
227  }
228 
229  for (i = 0; i < dim; i++)
230  for (j = 0; j < dim; j++)
231  mTdInverse(i, j) = mQ(j, i);
232 
233  C_INT flag_schur;
234 
235  flag_schur = 1; // set flag_schur = 0 to print Schur decomposition of jacobian (matrices mR and transformation mQ)
236 
237  mY_cons.resize(dim); // consistent initial vector for DAE
238 
239  /* If complex eigenvalues */
240  //BUG 873
241  if (mR(dim - 1, dim - 1) == mR(dim - 2 , dim - 2))
242  if (dim == 2)
243  {
244  slow = dim;
245  goto integration;
246  }
247 
248  // If positive eigenvalues
249 
250  if (mR(dim - 1, dim - 1) >= 0)
251  {
252  slow = dim;
253  fast = 0;
255  MCTSSAMethod + 10, mTime - deltaT);
256 
257  failed = 1;
258  goto integration;
259  }
260 
261  // Iterations to determine the number of slow metabolites
262 
263  while ((slow2 > 1))
264  {
265  slow2 = slow2 - 1;
266  fast2 = dim - slow2;
267 
268  if (mR(dim - fast2, dim - fast2) >= 0)
269  {
270  failed = 1;
271  goto integration;
272  }
273 
274  deuflhard_metab(slow2, info);
275 
276  if (info)
277  {
278  failed_while = 1;
279  goto integration;
280  }
281  }
282 
283  //end of iterations to determine the number of slow metabolites
284 
285  /** end of the block %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
286  /** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
287 
288 integration:
289 
290  slow = slow2;
291  fast = fast2;
292 
293  if ((failed == 1) || (failed_while == 1))
294  {
295  if (slow < dim)
296  {
297  fast = fast - 1;
298  slow = dim - fast;
299 
300  if ((fast >= 1) && (mR(slow - 1, slow - 1) == mR(slow , slow)))
301  fast = fast - 1;
302 
303  slow = dim - fast;
304  }
305  }
306 
307  mSlow = slow;
308 
309  if (slow == dim)
311  MCTSSAMethod + 11, mTime);
312 
313  for (i = 0; i < dim; i++)
314  for (j = 0; j < dim; j++)
315  // mTdInverse(i,j) = mQ(i,j);
316  mTd(i, j) = mQ(i, j);
317 
318  // Flag for print Tabs
319 
320  mat_anal_mod(slow);
321  mat_anal_metab(slow);
322  mat_anal_mod_space(slow);
323  mat_anal_fast_space(slow);
324 
325  // This block proves which metabolite could be considered as QSS. In development
326  /** Begin of the block %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
327 
328  C_INT flag_dev;
329  flag_dev = 1;
330 
331  if (flag_dev == 0)
332  {
333  C_INT temp;
334  temp = dim - 1;
335 
336  // Mode Analysis to find the dominant metabolites in the modes
337  mat_anal_mod(temp);
338 
339  //index_metab = -1, if there is no dominant metabolites in corresponding mode,
340  //and is equal
341  // to the number of dominant metabolite in opposite case
342  // Dominant metabolite - its contribution to the mode is larger as 70%
343 
344  for (j = 0; j < dim; j++)
345  index_metab[j] = -1;
346 
347  for (i = 0; i < dim ; i ++)
348  for (j = 0; j < dim; j++)
349  if (mVslow(dim - i - 1, j) > 70)
350  index_metab[i] = j;
351 
352  C_FLOAT64 y_cons;
353 
354  info = 0;
355  k = 0;
356  number = index_metab[k];
357 
358  if (number > - 1)
359  newton_for_timestep(number, y_cons, info);
360 
361  while (k < dim - 1)
362  {
363  if (number > -1)
364  {
365  dxdt.resize(dim);
366 
367  for (j = 0; j < dim; j++)
368  dxdt[j] = 0.;
369 
370  //CVector<C_FLOAT64> x_help;
371  x_help.resize(dim);
372 
373  for (j = 0; j < dim; j++)
374  {
375  x_help[j] = mY_initial[j] * number2conc;
376  }
377 
378  calculateDerivativesX(x_help.array(), dxdt.array());
379  info = 0;
380 
381  //NEWTON: Looking for consistent initial value for DAE system
382  //Output: y_cons, info
383 
384  newton_for_timestep(number, y_cons, info);
385 
386  if (info)
387  {
388  // TODO info: newton iteration stop
389  }
390 
391  if (info == 0)
392  {
393  // CVector<C_FLOAT64> x_relax;
394  x_relax.resize(dim);
395 
396  for (i = 0; i < dim; i ++)
397  if (i == number)
398  x_relax[i] = y_cons;
399  else
400  x_relax[i] = x_help[i];
401 
402  //CVector<C_FLOAT64> dxdt_relax;
403  dxdt_relax.resize(dim);
404 
405  calculateDerivativesX(x_relax.array(), dxdt_relax.array());
406 
407  //CVector<C_FLOAT64> re;
408  re.resize(dim);
409 
410  C_FLOAT64 eps;
411  eps = 1 / fabs(mR(dim - k - 1 , dim - k - 1));
412 
413  // stop criterion for slow reaction modes
414 
415  for (i = 0; i < dim; i++)
416  {
417  if (i == number)
418  re[i] = 0;
419  else
420  {
421  re[i] = fabs(dxdt_relax[i] - dxdt[i]);
422  re[i] = re[i] * eps;
423  }
424  }
425 
426  //C_FLOAT64 max = 0.;
427  for (i = 0; i < dim; i++)
428  if (max < re[i])
429  max = re[i];
430 
431  if (max >= mDtol)
432  info = 1;
433  else
434  info = 0;
435  }
436  }
437 
438  k = k + 1;
439  number = index_metab[k];
440  max = 0;
441  }
442  }
443 
444  /** end of the of block %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
445  /** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
446 
448  // TO REMOVE : mpModel->applyAssignments();
449 
450  // Calculate Jacobian for time step control
451 
452  mpModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
453 
454  // new entry for every entry contains the current data of currently step
455  setVectors(slow);
456 
457  // set the stepcounter
458  mCurrentStep += 1;
459 
460  return;
461 }
462 
463 /**
464  EVALSORT for vector sorting
465  */
466 
467 void CILDMModifiedMethod::evalsort(C_FLOAT64 *reval, C_INT *index, const C_INT & dim)
468 {
469  C_INT i, j, min, tmp_ind;
470  C_FLOAT64 tmp1r;
471 
472  for (i = 0; i < dim - 1; i++)
473  {
474  min = i;
475 
476  for (j = i + 1; j < dim; j++)
477  {
478  if (reval[j] <= reval[min])
479  min = j;
480  }
481 
482  tmp1r = reval[min];
483  tmp_ind = index[min];
484  reval[min] = reval[i];
485  index[min] = index[i];
486  reval[i] = tmp1r;
487  index[i] = tmp_ind;
488  }
489 
490  return;
491 }
492 
493 /**
494  Deuflhard Iteration: Prove Deuflhard criteria, find consistent initial value for DAE
495  output: info - if Deuflhard is satisfied
496  */
498 {
499  C_INT i, j, info_newton;
500  C_INT dim = mData.dim;
501  C_INT fast = dim - slow;
502 
503  C_INT flag_deufl;
504 
505  flag_deufl = 1; // set flag_deufl=0 to print temporaly steps for this function
506 
507  C_FLOAT64 max = 0;
509  CVector<C_FLOAT64> dxdt_relax;
510  CVector<C_FLOAT64> x_relax;
511  CVector<C_FLOAT64> x_help;
512  CVector<C_FLOAT64> dxdt;
513 
514  CVector<C_FLOAT64> help;
515  help.resize(dim);
516  CVector<C_INT> index;
517  index.resize(dim);
518  CVector<C_INT> index_temp;
519  index_temp.resize(dim);
520 
521  C_FLOAT64 eps;
522  eps = 1 / fabs(mR(dim - fast , dim - fast));
523  //eps = fabs(mR(dim - fast - 1, dim - fast - 1)) / fabs(mR(dim - fast , dim - fast));
524 
525  mat_anal_fast_space(slow);
526 
527  for (i = 0; i < dim; i++)
528  {
529  index[i] = i;
530  index_temp[i] = i;
531  }
532 
533  for (i = 0; i < dim; i++)
534  {
535  help[i] = mVfast_space[i];
536  }
537 
538  evalsort(help.array(), index.array(), dim);
539 
540  for (i = 0; i < dim; i++)
541  index_temp[i] = index[i];
542 
543  for (i = 0; i < dim; i++)
544  {
545  index[i] = index_temp[dim - i - 1];
546  }
547 
548  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
549  //C_FLOAT64 number2conc = 1.;
550 
551  dxdt.resize(dim);
552 
553  for (j = 0; j < dim; j++)
554  dxdt[j] = 0.;
555 
556  //CVector<C_FLOAT64> x_help;
557  x_help.resize(dim);
558 
559  for (j = 0; j < dim; j++)
560  {
561  x_help[j] = mY_initial[j] * number2conc;
562  }
563 
564  // mpModel->calculateDerivativesX(dxdt.array());
565  calculateDerivativesX(x_help.array(), dxdt.array());
566 
567  info_newton = 0;
568 
569  newton_new(index.array(), slow, info_newton);
570 
571  if (info_newton)
572  {
573  // TODO
574  info = 1;
575  return;
576  }
577 
578  x_relax.resize(dim);
579 
580  for (i = 0; i < dim; i++)
581  x_relax[i] = mY_cons[i];
582 
583  //CVector<C_FLOAT64> dxdt_relax;
584  dxdt_relax.resize(dim);
585 
586  calculateDerivativesX(x_relax.array(), dxdt_relax.array());
587 
588  //CVector<C_FLOAT64> re;
589  re.resize(dim);
590 
591  // stop criterion for slow reaction modes
592 
593  for (i = 0; i < dim; i++)
594  {
595  re[i] = fabs(dxdt_relax[i] - dxdt[i]);
596  re[i] = re[i] * eps;
597 
598  for (j = 0; j < fast; j ++)
599  if (i == index[j])
600  re[i] = 0;
601  }
602 
603  for (i = 0; i < dim; i++)
604  if (max < re[i])
605  max = re[i];
606 
607  if (max >= mDtol)
608  info = 1;
609  else
610  info = 0;
611 
612  return;
613 }
614 
615 void CILDMModifiedMethod::newton_new(C_INT *index_metab, C_INT & slow, C_INT & info)
616 {
617  C_INT i, j, k, m, iter, iterations, itermax;
618  C_INT nrhs, ok, fast, flag_newton;
619 
620  flag_newton = 1; // set flag_newton=1 to print temporaly steps of newton iterations
621 
622  C_FLOAT64 tol, err;
623  C_INT dim = mData.dim;
624 
625  fast = dim - slow;
626 
627  CVector<C_INT> ipiv;
628  ipiv.resize(fast);
629 
630  CVector<C_FLOAT64> s_22_array;
631  s_22_array.resize(fast * fast);
632 
633  CVector<C_FLOAT64> gf_newton;
634  gf_newton.resize(fast);
635 
636  CVector<C_FLOAT64> d_yf;
637  d_yf.resize(dim);
638 
639  CVector<C_FLOAT64> y_newton;
640  y_newton.resize(dim);
641 
642  CVector<C_FLOAT64> yf_newton;
643  yf_newton.resize(fast);
644 
645  CVector<C_FLOAT64> dydt_newton;
646  dydt_newton.resize(dim);
647 
648  CVector<C_FLOAT64> g_newton;
649  g_newton.resize(dim);
650 
651  CMatrix<C_FLOAT64> Jac_fast;
652  Jac_fast.resize(fast, fast);
653 
654  C_FLOAT64 g1, g2 = 0.0;
655  mY_cons.resize(dim);
656 
657  nrhs = 1;
658  tol = 1e-6;
659  err = 10.0;
660  iter = 0;
661 
662  itermax = 150;
663  iterations = 0;
664 
665  info = 0;
666 
667  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
668  //C_FLOAT64 number2conc = 1.;
669 
670  for (i = 0; i < fast; i++)
671  {
672  m = index_metab[i];
673 
674  for (j = 0; j < fast; j++)
675  {
676  k = index_metab[ j];
677 
678  if ((m > -1) & (k > -1))
679  Jac_fast(i, j) = mJacobian_initial(m, k);
680  else
681  {
682  info = 3;
683  return;
684  }
685  }
686  }
687 
688  for (i = 0; i < dim; i++)
689  y_newton[i] = mY_initial[i] * number2conc;
690 
691  for (i = 0; i < fast; i++)
692  for (j = 0; j < fast; j++)
693  s_22_array[j + fast * i] = Jac_fast(j, i);
694 
695  for (i = 0; i < dim; i++)
696  d_yf[i] = 0.;
697 
698  while (err > tol)
699  {
700  iter ++;
701 
702  if (iter > itermax)
703  {
704  info = 1;
705  return;
706  }
707 
708  for (i = 0; i < dim; i++)
709  y_newton[i] = y_newton[i] + d_yf[i];
710 
711  calculateDerivativesX(y_newton.array(), dydt_newton.array());
712 
713  for (i = 0; i < fast; i++)
714  {
715  j = index_metab[i];
716  gf_newton[i] = -1. * dydt_newton[j];
717  }
718 
719  /* int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
720  * *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info)
721  *
722  * -- LAPACK driver routine (version 3.0) --
723  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
724  * Courant Institute, Argonne National Lab, and Rice University
725  * March 31, 1993
726  *
727  *
728  * Purpose
729  * =======
730  *
731  * DGESV computes the solution to a real system of linear equations
732  * A * X = B,
733  * where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
734  *
735  * The LU decomposition with partial pivoting and row interchanges is
736  * used to factor A as
737  * A = P * L * U,
738  * where P is a permutation matrix, L is unit lower triangular, and U is
739  * upper triangular. The factored form of A is then used to solve the
740  * system of equations A * X = B.
741  *
742  * Arguments
743  * =========
744  *
745  * N (input) INTEGER
746  * The number of linear equations, i.e., the order of the
747  * matrix A. N >= 0.
748  *
749  * NRHS (input) INTEGER
750  * The number of right hand sides, i.e., the number of columns
751  * of the matrix B. NRHS >= 0.
752  *
753  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
754  * On entry, the N-by-N coefficient matrix A.
755  * On exit, the factors L and U from the factorization
756  * A = P*L*U; the unit diagonal elements of L are not stored.
757  *
758  * LDA (input) INTEGER
759  * The leading dimension of the array A. LDA >= max(1,N).
760  *
761  * IPIV (output) INTEGER array, dimension (N)
762  * The pivot indices that define the permutation matrix P;
763  * row i of the matrix was interchanged with row IPIV(i).
764  *
765  * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
766  * On entry, the N-by-NRHS matrix of right hand side matrix B.
767  * On exit, if INFO = 0, the N-by-NRHS solution matrix X.
768  *
769  * LDB (input) INTEGER
770  * The leading dimension of the array B. LDB >= max(1,N).
771  *
772  *
773  * INFO (output) INTEGER
774  * = 0: successful exit
775  * < 0: if INFO = -i, the i-th argument had an illegal value
776  * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
777  * has been completed, but the factor U is exactly
778  * singular, so the solution could not be computed.
779  */
780 
781  dgesv_(&fast, &nrhs, s_22_array.array(), &fast, ipiv.array(), gf_newton.array(), &fast, &ok);
782 
783  if (ok != 0)
784  {
785  info = 2;
786  return;
787  }
788 
789  for (i = 0; i < fast; i++)
790  d_yf[i] = 0;
791 
792  for (j = 0; j < fast; j++)
793  {
794  k = index_metab[j];
795 
796  for (i = 0; i < dim; i ++)
797  {
798  if (i == k)
799  d_yf[k] = gf_newton[j];
800  }
801  }
802 
803  err = -10.;
804 
805  for (i = 0; i < fast; i++)
806  {
807  gf_newton[i] = fabs(gf_newton[i]);
808 
809  if (err < gf_newton[i])
810  err = gf_newton[i];
811  }
812 
813  iterations = iterations + 1;
814 
815  // stop criterion of newton method
816 
817  // C_FLOAT64 g1, g2 = 0.0;
818 
819  // g2 = err;
820 
821  if (iter == 1)
822  g1 = 3.0 * err;
823  else
824  g1 = g2;
825 
826  g2 = err;
827 
828  if (g2 / g1 > 1.0)
829  {
830  info = 1;
831  return;
832  }
833  } /* end while */
834 
835  for (i = 0; i < dim; i++)
836  mY_cons[i] = y_newton[i];
837 
838  info = 0;
839 
840  return;
841 }
842 
843 /**
844 NEWTON for "postprove": Prove of "fast" varibles
845 Output: y_consistent, info
846  */
847 
848 void CILDMModifiedMethod::newton_for_timestep(C_INT metabolite_number, C_FLOAT64 & y_consistent, C_INT & info)
849 {
850  C_INT i, iter, itermax, flag_newton;
851  iter = 0;
852  itermax = 150;
853 
854  flag_newton = 1; // set flag_newton=0 to print the iteration steps
855 
856  C_FLOAT64 tol, err;
857  tol = 1e-6;
858  err = 10.0;
859 
860  C_INT dim = mData.dim;
861 
862  C_FLOAT64 d_y, deriv;
863  d_y = 0;
864  deriv = mJacobian_initial(metabolite_number, metabolite_number);
865 
866  if (deriv == 0)
867  {
868  return;
869  }
870 
871  info = 0;
872 
873  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
874  //C_FLOAT62number2conc = 1.;
875 
876  //this is an ugly hack that only makes sense if all metabs are in the same compartment
877  //at the moment is is the only case the algorithm deals with
878 
879  CVector<C_FLOAT64> y_newton; //current state converted to concentrations
880  y_newton.resize(dim);
881 
882  for (i = 0; i < dim; ++i)
883  y_newton[i] = mY_initial[i] * number2conc;
884 
885  CVector<C_FLOAT64> dydt;
886  dydt.resize(dim);
887 
888  while (err > tol)
889  {
890  iter ++;
891 
892  if (iter > itermax)
893  {
894  info = 1;
895  break;
896  }
897 
898  y_newton[metabolite_number] = y_newton[metabolite_number] + d_y;
899 
900  calculateDerivativesX(y_newton.array(), dydt.array());
901 
902  d_y = - 1 / deriv * dydt[metabolite_number];
903 
904  if (err > fabs(d_y))
905  err = fabs(d_y);
906  }
907 
908  y_consistent = y_newton[metabolite_number];
909 
910  return;
911 }
912 
913 /**
914  * Empty every vector to be able to fill them with new values for a new calculation.
915  * Also nullify the step counter.
916  **/
918 {
919  mCurrentStep = 0;
920  mVec_mVslow.erase(mVec_mVslow.begin(), mVec_mVslow.end());
921  mVec_TimeScale.erase(mVec_TimeScale.begin(), mVec_TimeScale.end());
922  mVec_mVslowMetab.erase(mVec_mVslowMetab.begin(), mVec_mVslowMetab.end());
923  mVec_mVslowSpace.erase(mVec_mVslowSpace.begin(), mVec_mVslowSpace.end());
924  mVec_SlowModes.erase(mVec_SlowModes.begin(), mVec_SlowModes.end());
925 }
926 
927 /**
928  *upgrade all vectors with values from actually calculalion for current step
929  **/
931 {
932  mVec_mVslow.push_back(mCurrentStep);
935 
936  mVec_TimeScale.push_back(mCurrentStep);
938  int i;
939 
940  for (i = 0; i < mData.dim; i++)
941  mVec_TimeScale[mCurrentStep][i] = -1 / mR(i, i);
942 
943  mVec_mVslowMetab.push_back(mCurrentStep);
946 
947  mVec_mVslowSpace.push_back(mCurrentStep);
950 
951  mVec_mVfastSpace.push_back(mCurrentStep);
954 
955  mVec_SlowModes.push_back(mCurrentStep);
956  mVec_SlowModes[mCurrentStep] = slowMode;
957 
958  mCurrentTime.push_back(mCurrentStep);
960 }
961 
962 /**
963  * Create the CArraAnnotations for every ILDM-tab in the CQTSSAResultSubWidget.
964  * Input for each CArraAnnotations is a seperate CMatrix.
965  **/
967 {
968 
969  tableNames.erase(tableNames.begin(), tableNames.end());
970 
971  std::string name;
972 
973  name = "Contribution of species to modes";
974  tableNames.push_back(name);
975 
977  pTmp1 = new CArrayAnnotation("Contribution of species to modes", this,
979  pTmp1->setMode(0, pTmp1->STRINGS);
980  pTmp1->setMode(1, pTmp1->VECTOR);
981  pTmp1->setDescription(" ");
982  pTmp1->setDimensionDescription(0, "Contribution to mode (TS - corresponding timescale)");
983  pTmp1->setDimensionDescription(1, "Species");
985 
987 
988  name = "Modes distribution for species";
989  tableNames.push_back(name);
990 
992  pTmp2 = new CArrayAnnotation("Modes distribution for species", this,
994  pTmp2->setMode(1, pTmp2->STRINGS);
995  pTmp2->setMode(0, pTmp2->VECTOR);
996  pTmp2->setDescription(" ");
997  pTmp2->setDimensionDescription(0, "Mode distribution for each metabolite");
998  pTmp2->setDimensionDescription(1, "modes (TS - corresponding timescale)");
1000 
1002 
1003  name = "Slow space";
1004  tableNames.push_back(name);
1005 
1007  pTmp3 = new CArrayAnnotation("Slow space", this,
1009  pTmp3->setMode(1, pTmp3->STRINGS);
1010  pTmp3->setMode(0, pTmp3->VECTOR);
1011  pTmp3->setDescription(" ");
1012  pTmp3->setDimensionDescription(0, "Species");
1013  pTmp3->setDimensionDescription(1, "Contribution to slow space");
1015 
1017 
1018  name = "Fast space";
1019  tableNames.push_back(name);
1020 
1022  pTmp4 = new CArrayAnnotation("Fast space", this,
1024  pTmp4->setMode(1, pTmp4->STRINGS);
1025  pTmp4->setMode(0, pTmp4->VECTOR);
1026  pTmp4->setDescription(" ");
1027  pTmp4->setDimensionDescription(0, "Species");
1028  pTmp4->setDimensionDescription(1, "Contribution to fast space");
1030 
1032 }
1033 /**
1034  * Set the every CArrayAnnotation for the requested step.
1035  * Set also the desription of CArayAnnotation for both dimensions:
1036  * - dimension description could consists of some std::srings
1037  * some strings contain the Time Scale values for requested step
1038  * - dimension description could consists of arrays of CommonNames
1039  **/
1040 //void CILDMModifiedMethod::setAnnotationM(int step)
1042 {
1043  if (step == 0) return false;
1044 
1045  if (mVec_mVslow.size() == 0) return false;
1046 
1047  //if (step > mVec_mVslow.size()) return false;
1048 
1049  if (step > mVec_SlowModes.size()) return false;
1050 
1051  step -= 1;
1052  double timeScale;
1053  std::string str;
1054  std::stringstream sstr;
1055  sstr.str("");
1056  sstr.clear();
1057  int i;
1058 
1063 
1064  for (i = 0; i < mData.dim; i++)
1065  {
1066  timeScale = mVec_TimeScale[step][i];
1067 
1068  if (i < mVec_SlowModes[step])
1069  sstr << "Slow: ";
1070  else
1071  sstr << "Fast: ";
1072 
1073  sstr << timeScale;
1074  str = sstr.str();
1075  pVslowPrintAnn->setAnnotationString(0, i, str);
1076  sstr.str("");
1077  sstr.clear();
1078  }
1079 
1084 
1085  for (i = 0; i < mData.dim; i++)
1086  {
1087  timeScale = mVec_TimeScale[step][i];
1088 
1089  if (i < mVec_SlowModes[step])
1090  sstr << "Slow: ";
1091  else
1092  sstr << "Fast: ";
1093 
1094  sstr << timeScale;
1095  str = sstr.str();
1097  sstr.str("");
1098  sstr.clear();
1099  }
1100 
1101  sstr << mVec_SlowModes[step];
1102  // if (mVec_SlowModes[step] > 1)
1103  // sstr << " slow modes";
1104  //else
1105  // sstr << " slow mode";
1106  sstr << " slow; ";
1107 
1108  C_INT dim = mData.dim;
1109  sstr << dim - mVec_SlowModes[step];
1110  sstr << " fast";
1111 
1112  str = sstr.str();
1114 
1115  for (i = 0; i < mData.dim; i++)
1117 
1121 
1123 
1124  for (i = 0; i < mData.dim; i++)
1126 
1130 
1131  return true;
1132 }
1133 
1134 void CILDMModifiedMethod::printResult(std::ostream * ostream) const
1135 {
1136  std::ostream & os = *ostream;
1137  double timeScale;
1138  C_INT i, j, istep = 0;
1139 
1140  this->print(&os);
1141 
1142  C_INT32 stepNumber;
1143 
1144  assert(getObjectDataModel() != NULL);
1145 
1146  //stepNumber = pProblem->getStepNumber();
1147  stepNumber = mVec_SlowModes.size();
1148 
1149  for (istep = 0; istep < stepNumber; istep++)
1150  {
1151 
1152  os << std::endl;
1153  os << "**************** Time step " << istep + 1 << " ************************** " << std::endl;
1154 
1155  os << std::endl;
1156 
1157  os << "Contribution of species to modes" << std::endl;
1158 
1159  os << "Rows : contribution to mode (TS - corresponding timescale)" << std::endl;
1160  os << "Columns: species ";
1161 
1162  for (j = 0; j < mData.dim; j++)
1163  {
1164  os << mpModel->getMetabolitesX()[j]->getObjectName() << " ";
1165  }
1166 
1167  os << std::endl;
1168 
1169  for (i = 0; i < mData.dim; i++)
1170  {
1171  timeScale = mVec_TimeScale[istep][i];
1172 
1173  if (i < mVec_SlowModes[istep])
1174  os << "Slow (";
1175  else
1176  os << "Fast (";
1177 
1178  os << timeScale << "): ";
1179 
1180  for (j = 0; j < mData.dim; j++)
1181  os << mVec_mVslow[istep][i][j] << " ";
1182 
1183  os << std::endl;
1184  }
1185 
1186  os << std::endl;
1187 
1188  os << "Modes distribution for species" << std::endl;
1189 
1190  os << "Rows : Mode distribution for each species" << std::endl;
1191  os << "Columns: Modes (TS - corresponding timescale) ";
1192  os << std::endl;
1193 
1194  for (i = 0; i < mData.dim; i++)
1195  {
1196  timeScale = mVec_TimeScale[istep][i];
1197 
1198  if (i < mVec_SlowModes[istep])
1199  os << "Slow (";
1200  else
1201  os << "Fast (";
1202 
1203  os << timeScale << ") ";
1204  }
1205 
1206  os << std::endl;
1207 
1208  for (j = 0; j < mData.dim; j++)
1209  {
1210  os << mpModel->getMetabolitesX()[j]->getObjectName() << " ";
1211 
1212  for (i = 0; i < mData.dim; i++)
1213  os << mVec_mVslowMetab[istep][j][i] << " ";
1214 
1215  os << std::endl;
1216  }
1217 
1218  os << std::endl;
1219 
1220  os << "Slow space" << std::endl;
1221 
1222  os << "Rows : Species" << std::endl;
1223  os << "Column: Contribution to slow space ";
1224  os << std::endl;
1225 
1226  os << mVec_SlowModes[istep];
1227  os << " slow; ";
1228 
1229  os << mData.dim - mVec_SlowModes[istep];
1230  os << " fast";
1231  os << std::endl;
1232 
1233  for (j = 0; j < mData.dim; j++)
1234  {
1235  os << mpModel->getMetabolitesX()[j]->getObjectName() << " ";
1236  os << mVec_mVslowSpace[istep][j] << " ";
1237 
1238  os << std::endl;
1239  }
1240 
1241  os << std::endl;
1242  os << "Fast space" << std::endl;
1243 
1244  os << "Rows : Species" << std::endl;
1245  os << "Column: Contribution to fast space ";
1246  os << std::endl;
1247 
1248  os << mVec_SlowModes[istep];
1249  os << " slow; ";
1250 
1251  os << mData.dim - mVec_SlowModes[istep];
1252  os << " fast";
1253  os << std::endl;
1254 
1255  for (j = 0; j < mData.dim; j++)
1256  {
1257  os << mpModel->getMetabolitesX()[j]->getObjectName() << " ";
1258  os << mVec_mVfastSpace[istep][j] << " ";
1259 
1260  os << std::endl;
1261  }
1262 
1263  os << std::endl;
1264  }
1265 
1266  return;
1267 }
void schur(C_INT &info)
CCopasiDataModel * getObjectDataModel()
std::vector< CMatrix< C_FLOAT64 > > mVec_mVslowMetab
#define C_INT
Definition: copasi.h:115
int dgesv_(integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info)
#define MCTSSAMethod
CArrayAnnotation * pTmp4
#define pdelete(p)
Definition: copasi.h:215
CArrayAnnotation * pTmp2
std::vector< C_FLOAT64 > mCurrentTime
Definition: CTSSAMethod.h:460
CMatrix< C_FLOAT64 > mVslowPrint
const std::string & getObjectName() const
CMatrix< C_FLOAT64 > mVslow
Definition: CTSSAMethod.h:261
CMatrix< C_FLOAT64 > mTdInverse
Definition: CTSSAMethod.h:230
CVector< C_FLOAT64 > mY_initial
Definition: CTSSAMethod.h:195
CState * mpState
Definition: CTSSAMethod.h:175
void evalsort(C_FLOAT64 *reval, C_INT *index, const C_INT &dim_x)
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
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
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
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
void setDescription(const std::string &s)
CArrayAnnotation * pVslowPrintAnn
void deuflhard_metab(C_INT &slow, C_INT &info)
std::vector< CMatrix< C_FLOAT64 > > mVec_mVslow
#define C_INT32
Definition: copasi.h:90
virtual void step(const double &deltaT)
virtual void initializeParameter()
CMatrix< C_FLOAT64 > mQ
Definition: CTSSAMethod.h:215
CMatrix< C_FLOAT64 > mTdInverse_save
Definition: CTSSAMethod.h:244
void mat_anal_mod_space(C_INT &slow)
CMatrix< C_FLOAT64 > mVfastSpacePrint
CVector< C_FLOAT64 > mVfast_space
Definition: CTSSAMethod.h:276
void setAnnotationString(size_t d, size_t i, const std::string s)
CArrayAnnotation * pVslowMetabPrintAnn
CMatrix< C_FLOAT64 > mJacobian_initial
Definition: CTSSAMethod.h:210
CArrayAnnotation * pTmp1
CArrayAnnotation * pVslowSpacePrintAnn
virtual bool setAnnotationM(size_t step)
const C_FLOAT64 & getNumber2QuantityFactor() const
Definition: CModel.cpp:2357
void calculateDerivativesX(C_FLOAT64 *X1, C_FLOAT64 *Y1)
std::vector< std::string > tableNames
Definition: CTSSAMethod.h:98
void mat_anal_fast_space(C_INT &slow)
int mCurrentStep
Definition: CTSSAMethod.h:466
void setVectors(int slowMode)
const Value & getValue() const
bool mReducedModel
Definition: CTSSAMethod.h:291
CVector< C_FLOAT64 > mY_cons
Definition: CTSSAMethod.h:256
void setCopasiVector(size_t d, const CCopasiContainer *v)
CMatrix< C_FLOAT64 > mQz
Definition: CTSSAMethod.h:235
CMatrix< C_FLOAT64 > mVslowMetabPrint
virtual void print(std::ostream *ostream) const
void newton_new(C_INT *index_metab, C_INT &slow, C_INT &info)
virtual void start(const CState *initialState)
void newton_for_timestep(C_INT metabolite_number, C_FLOAT64 &y_consistent, C_INT &info)
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
void mat_anal_metab(C_INT &slow)
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
std::vector< CVector< C_FLOAT64 > > mVec_TimeScale
Definition: CTSSAMethod.h:461
CMatrix< C_FLOAT64 > mJacobian
Definition: CTSSAMethod.h:205
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
std::vector< CVector< C_FLOAT64 > > mVec_mVslowSpace
CMatrix< C_FLOAT64 > mTd
Definition: CTSSAMethod.h:225
void printResult(std::ostream *ostream) const
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
CILDMModifiedMethod(const CCopasiContainer *pParent=NULL)
CMatrix< C_FLOAT64 > mVslowSpacePrint
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
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)
std::vector< CVector< C_FLOAT64 > > mVec_mVfastSpace
void integrationMethodStart(const CState *initialState)
#define min(a, b)
Definition: f2c.h:175
std::vector< C_INT > mVec_SlowModes
Definition: CTSSAMethod.h:459
CTSSAMethod * pMethod
Definition: CTSSAMethod.h:169
void mat_anal_mod(C_INT &slow)
CArrayAnnotation * pTmp3
CArrayAnnotation * pVfastSpacePrintAnn
CMatrix< C_FLOAT64 > mVslow_metab
Definition: CTSSAMethod.h:266
C_FLOAT64 * mY
Definition: CTSSAMethod.h:185
#define max(a, b)
Definition: f2c.h:176