COPASI API  4.16.103
CCSPMethod.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 //#define CSPDEBUG
16 
17 #include <math.h>
18 
19 #include "copasi.h"
20 
21 #include "CCSPMethod.h"
22 #include "CTSSAProblem.h"
23 #include "CTSSATask.h"
26 #include "model/CModel.h"
27 #include "model/CMetab.h"
28 #include "model/CState.h"
29 #include "utilities/CMatrix.h"
32 
33 #include "lapack/lapackwrap.h" // CLAPACK
34 #include "lapack/blaswrap.h" // BLAS
35 
37  CTSSAMethod(CCopasiMethod::tssCSP, pParent) //,
38  //mpState(NULL),
39  //mY(NULL)
40 {
41  assert((void *) &mData == (void *) &mData.dim);
42 
43  mData.pMethod = this;
45 
47  //emptyVectors();
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 
62  //emptyVectors();
63 }
64 
66 {
68 }
69 
71 {
73 
74  assertParameter("Integrate Reduced Model", CCopasiParameter::BOOL, (bool) false);//->getValue().pBOOL;
75  assertParameter("Ratio of Modes Separation", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-2);
76  assertParameter("Maximum Relative Error", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-5);
77  assertParameter("Maximum Absolute Error", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-10);
78  assertParameter("Refinement Iterations Number", CCopasiParameter::UINT, (unsigned C_INT32) 1000);
79 
80  //createAnnotationsM();
81  //emptyVectors();
82 }
83 
84 /* multiply submatrix */
86 {
87  C_INT i, j, k;
88 
89  for (i = 0; i < n1 ; i++)
90  for (j = 0; j < n3; j++)
91  {
92  C(i, j) = 0.;
93 
94  for (k = 0; k < n2; k++)
95  C(i, j) += A(i, k) * B(k, j);
96  }
97 
98  return;
99 }
100 
101 /* subtract submatrix */
103 {
104  C_INT i, j;
105 
106  for (i = 0; i < n1 ; i++)
107  for (j = 0; j < n2; j++)
108  C(i, j) = A(i, j) - B(i, j);
109 
110  return;
111 }
112 
113 /* add submatrix */
115 {
116  C_INT i, j;
117 
118  for (i = 0; i < n1 ; i++)
119  for (j = 0; j < n2; j++)
120  C(i, j) = A(i, j) + B(i, j);
121 
122  return;
123 }
124 
125 /* normalize submatrix */
127 {
128  C_INT i, j;
129  C_FLOAT64 c, d;
130 
131  for (j = 0; j < n1 ; j++)
132  {
133  c = 0.0;
134 
135  for (i = 0; i < n ; i++)
136  {
137  d = fabs(A(i, j));
138 
139  if (d > c) c = d;
140  }
141 
142  for (i = 0; i < n ; i++)
143  {
144  A(i, j) = A(i, j) / c;
145  B(j, i) = B(j, i) * c;
146  }
147  }
148 
149  return;
150 }
151 
152 /* perturbate basis */
154 {
155  C_INT i, j;
156 
157  for (j = 0; j < n ; j++)
158  for (i = 0; i < n ; i++)
159  A(i, j) = A(i, j) * delta;
160 
161  return;
162 }
163 
165 {
166 
167  /* int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
168  * *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info)
169  *
170  * -- LAPACK driver routine (version 3.0) --
171  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
172  * Courant Institute, Argonne National Lab, and Rice University
173  * March 31, 1993
174  *
175  *
176  * Purpose
177  * =======
178  *
179  * DGESV computes the solution to a real system of linear equations
180  * A * X = B,
181  * where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
182  *
183  * The LU decomposition with partial pivoting and row interchanges is
184  * used to factor A as
185  * A = P * L * U,
186  * where P is a permutation matrix, L is unit lower triangular, and U is
187  * upper triangular. The factored form of A is then used to solve the
188  * system of equations A * X = B.
189  *
190  * Arguments
191  * =========
192  *
193  * N (input) INTEGER
194  * The number of linear equations, i.e., the order of the
195  * matrix A. N >= 0.
196  *
197  * NRHS (input) INTEGER
198  * The number of right hand sides, i.e., the number of columns
199  * of the matrix B. NRHS >= 0.
200  *
201  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
202  * On entry, the N-by-N coefficient matrix A.
203  * On exit, the factors L and U from the factorization
204  * A = P*L*U; the unit diagonal elements of L are not stored.
205  *
206  * LDA (input) INTEGER
207  * The leading dimension of the array A. LDA >= max(1,N).
208  *
209  * IPIV (output) INTEGER array, dimension (N)
210  * The pivot indices that define the permutation matrix P;
211  * row i of the matrix was interchanged with row IPIV(i).
212  *
213  * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
214  * On entry, the N-by-NRHS matrix of right hand side matrix B.
215  * On exit, if INFO = 0, the N-by-NRHS solution matrix X.
216  *
217  * LDB (input) INTEGER
218  * The leading dimension of the array B. LDB >= max(1,N).
219  *
220  *
221  * INFO (output) INTEGER
222  * = 0: successful exit
223  * < 0: if INFO = -i, the i-th argument had an illegal value
224  * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
225  * has been completed, but the factor U is exactly
226  * singular, so the solution could not be computed.
227  */
228 
229  C_INT info = 0;
230 
231  C_INT * ipiv;
232  ipiv = new C_INT [n];
233 
234  C_INT N = n;
235 
236  CMatrix<C_FLOAT64> TMP;
237  TMP.resize(N, N);
238 
239  TMP = A;
240 
241  C_INT i, j;
242 
243  for (i = 0; i < N; i++)
244  for (j = 0; j < N; j++)
245  B(i, j) = 0.;
246 
247  for (i = 0; i < N; i++)
248  B(i, i) = 1.;
249 
250  dgesv_(&N, &N, TMP.array(), &N, ipiv, B.array(), &N, &info);
251 
252  if (info != 0)
253  {
254  return;
255  }
256 
257  return;
258 }
259 #if 0
260 /* find the new number of fast according to the time-scale separation ratio */
261 /* TODO : equal,complex eigenvalues are presenting !!! */
263 {
264 
265  C_INT i;
266  C_FLOAT64 tmp;
267 
268  k = 0;
269  i = 0;
270 
271  if (mTsc)
272  for (i = 0; i < n; i++)
273  {
274  tmp = 1. / eigen[i];
275 
276  if (tmp < 0 && fabs(tmp) < mTsc)
277  k++;
278 
279  else
280  {
281  if (tmp > 0) info = 1;
282 
283  break;
284  }
285  }
286  else
287  for (i = 0; i < n - 1; i++)
288  {
289  if (eigen[i] != eigen[i + 1])
290  {
291  tmp = eigen[i + 1 ] / eigen[i];
292 
293 #if 0
294 
295  std::cout << "tsc[" << i << "]/tsc[" << i + 1 << "] " << tmp << std::endl;
296  std::cout << "mEps " << mEps << std::endl;
297 #endif
298 
299  if (tmp > 0 && tmp < mEps)
300  {
301  k++;
302 
303  if (i)
304  if (eigen(i) == eigen(i - 1)) k++;
305  }
306  else
307  {
308  if (tmp < 0) info = 1;
309 
310  break;
311  }
312  }
313  else
314  {
315 
316 #if 0
317  std::cout << "the following time scales are equal: " << std::endl;
318  std::cout << "tsc[" << i << "] = tsc[" << i + 1 << "] " << std::endl;
319 #endif
320  }
321  }
322 
323  return;
324 }
325 #endif
326 
327 /* find the number of candidates to fast according to the time-scale separation ratio */
329 {
330 
331  C_INT i;
332  C_FLOAT64 tmp;
333 
334  k = 0;
335  i = 0;
336 
337  if (eigen[0] == 0.) return;
338 
339 #ifdef CSPDEBUG
340 
341  for (i = 0; i < n ; i++)
342  std::cout << "eigen[" << i << "] " << eigen[i] << std::endl;
343 
344  for (i = 0; i < n ; i++)
345  std::cout << "tsc[" << i << "] -" << 1 / eigen[i] << std::endl;
346 
347 #endif
348 
349  for (i = 0; i < n - 1; i++)
350  {
351 
352  if (eigen[i + 1] == 0.) return;
353 
354  if (eigen[i] != eigen[i + 1])
355  {
356  tmp = eigen[i + 1 ] / eigen[i];
357 
358 #ifdef CSPDEBUG
359 
360  std::cout << "tsc[" << i << "]/tsc[" << i + 1 << "] " << tmp << std::endl;
361  std::cout << "mEps " << mEps << std::endl;
362 #endif
363 
364  if (tmp > 0 && tmp < mEps)
365  {
366  k++;
367 
368 #ifdef CSPDEBUG
369  std::cout << "k " << k << std::endl;
370 #endif
371 
372  if (i)
373  if (eigen[i] == eigen[i - 1]) k++;
374  }
375  else
376  {
377  if (tmp < 0) info = 1;
378 
379  //if (tmp >= mEps) k++;
380 
381  break;
382  }
383  }
384  else
385  {
386 
387 #ifdef CSPDEBUG
388  std::cout << "the following time scales are equal: " << std::endl;
389  std::cout << "tsc[" << i << "] = tsc[" << i + 1 << "] " << std::endl;
390 #endif
391  }
392  }
393 
394 #ifdef CSPDEBUG
395  std::cout << "k " << k << std::endl;
396 #endif
397 
398  return;
399 }
400 
401 void CCSPMethod::cspstep(const double & /* deltaT */, C_INT & N, C_INT & M, CMatrix< C_FLOAT64 > & A, CMatrix< C_FLOAT64 > & B)
402 {
403 
404  C_INT reacs_size = (C_INT) mpModel->getReactions().size();
405  emptyOutputData(N, N, reacs_size);
406 
407 #ifdef CSPDEBUG
408  std::cout << " ********************* New time step **********************" << std::endl;
409  std::cout << "mTime " << mTime << std::endl;
410 
411 #endif
412 
415 
419 
420  g.resize(N);
421  y.resize(N);
422 
423  A0.resize(N, N);
424  B0.resize(N, N);
425  J.resize(N, N);
426 
427  C_INT i, j;
428 
429  for (j = 0; j < N; j++)
430  {
431  y[j] = mY[j];
432  g[j] = mG[j];
433 
434  CModelEntity* tmp;
435 
437 
438  CMetab* metab = dynamic_cast< CMetab * >(tmp);
439 
440  const CCompartment* comp = metab->getCompartment();
441 
442  mYerror[j] = mRerror * y[j] + mAerror * comp->getInitialValue();
443  }
444 
445  J = mJacobian;
446 
447 #ifdef CSPDEBUG
448  std::cout << "current Jacobian " << std::endl;
449  std::cout << J << std::endl;
450 #endif
451 
452  CMatrix<C_FLOAT64> ALA;
454 
455  ALA.resize(N, N);
456  F.resize(N, 1);
457 
458  /* csp iterations */
459  C_INT iter = 0;
460 
462  CMatrix<C_FLOAT64> QSL;
463 
464  QF.resize(N, N);
465  QSL.resize(N, N);
466 
468  mQ.resize(N, N);
469  mR.resize(N, N);
470 
471  mQ = 0.;
472  mR = 0.;
473 
474  mJacobian_initial = J;
475 
476  C_INT info = 0;
477 
478  schur(info);
479 
480  if (info)
481  {
483  MCTSSAMethod + 9, mTime);
484  return;
485  }
486 
487  /* trial basis vectors */
488 
489  /* use the matrix of Schur vectors */
490 
491  A0 = mQ;
492  B0 = 0;
493 
494  A = A0;
495  B = B0;
496 
497  //perturbateA(N, A0, 0.99); // TEST
498 
499  smnorm(N, A0, B0, N);
500  sminverse(N, A0, B0);
501 
502  /* ordered real parts of Eigen values */
503 
504  CVector<C_FLOAT64> eigen;
505  CVector<C_FLOAT64> tsc;
506 
507  eigen.resize(N);
508  tsc.resize(N);
509 
510  for (i = 0; i < N; i++)
511  {
512  eigen[i] = mR(i, i);
513 
514  tsc[i] = 1. / eigen[i];
515  }
516 
517 #ifdef CSPDEBUG
518  std::cout << "CSP iteration: " << std::endl;
519  std::cout << "Eigen values ordered by increasing " << std::endl;
520 
521  for (i = 0; i < N; i++)
522  std::cout << "eigen[" << i << "] " << eigen[i] << std::endl;
523 
524 #endif
525 
526 #ifdef CSPDEBUG
527  std::cout << "time scales : " << std::endl;
528 
529  for (i = 0; i < N; i++)
530  std::cout << fabs(tsc[i]) << std::endl;
531 
532 #endif
533 
534  /* find the number of candidate to fast */
535 
536  info = 0;
537 
538  //findTimeScaleSeparation(N, M, eigen, info);
539 
540  findCandidatesNumber(N, M, eigen, info);
541 
542 #ifdef CSPDEBUG
543  std::cout << " Candidates number " << M << std::endl;
544 #endif
545 
546  if (info)
547  {
548  /* the fastest of slow modes has positive eigen value */
550  MCTSSAMethod + 15, mTime);
551  }
552 
553  if (M == N)
554  {
555 
556  /* this case is not possible : */
557  /* if the ratio of time scale separation nearly 1, the time scale */
558  /* considered as slow */
559 #ifdef CSPDEBUG
560  std::cout << "After time scales separation : " << std::endl;
561  std::cout << "the number of candidates to fast modes is equal to the total number of modes" << std::endl;
562 #endif
563 
564  return;
565  }
566 
567  if (M == 0)
568  {
569 
570  /* After time scales separation : */
571  /* no candidates to fast modes */
572 
574  MCTSSAMethod + 12, mTime);
575 
576  return;
577  }
578 
579 analyseMmodes:
580 
581 #ifdef CSPDEBUG
582  std::cout << " ************************************** Number of candidates to fast " << M << " ************************" << std::endl;
583 #endif
584 
585  iter = 0;
586 
587  A = A0;
588  B = B0;
589 
590  /* */
591  /* ALA = B*J*A */
592 
593  CMatrix<C_FLOAT64> TMP;
594  TMP.resize(N, N);
595 
596 #if 0
597  CMatrix<C_FLOAT64> DBDT;
598 
599  DBDT.resize(N, N);
600 
601  DBDT = 0.;
602 
603  if (mTStep)
604  for (i = 0; i < N; i++)
605  for (j = 0; j < N; j++)
606  {
607  DBDT(i, j) = (B(i, j) - mB(i, j)) / deltaT;
608  }
609 
610 #ifdef CSPDEBUG
611  std::cout << "time derivatives of B " << std::endl;
612  std::cout << DBDT << std::endl;
613 #endif
614 #endif
615 
616  smmult(B, J, TMP, N, N, N);
617 
618 #if 0
619  /* TEST: time derivatives are present */
620 
621  if (mTStep)
622  for (i = 0; i < N; i++)
623  for (j = 0; j < N; j++)
624  {
625  TMP(i, j) += DBDT(i, j);
626  }
627 
628 #endif
629 
630  smmult(TMP, A, ALA, N, N, N);
631 
632 #ifdef CSPDEBUG
633  std::cout << "B*J*A should converge to block-diagonal for an ideal basis:" << std::endl;
634  std::cout << ALA << std::endl;
635 
636  std::cout << "considered time resolution of the solution " << std::endl;
637  std::cout << fabs(tsc[M]) << std::endl; // to check this
638 #endif
639 
640  CMatrix<C_FLOAT64> TAUM;
641  CMatrix<C_FLOAT64> ALAM;
642 
643  ALAM.resize(M, M);
644  TAUM.resize(M, M);
645 
646  TAUM = 0;
647  ALAM = 0;
648 
649  for (i = 0; i < M; i++)
650  for (j = 0; j < M; j++)
651  ALAM(i, j) = ALA(i, j);
652 
653  if (M > 1)
654  sminverse(M, ALAM, TAUM);
655  else
656  TAUM(0, 0) = 1. / ALA(0, 0);
657 
658 #if 1
659  modesAmplitude(N, M, g, B, F);
660 
661 #ifdef CSPDEBUG
662 
663  std::cout << " scaled amplitudes via trial basis : " << std::endl;
664 
665  for (i = 0; i < M; i++)
666  {
667 
668  std::cout << F(i, 0) << std::endl;
669  }
670 
671 #endif
672 
673 #ifdef CSPDEBUG
674  std::cout << " |A(i,m) * F(m,0) * tsc[M - 1]| , mYerror[i] " << std::endl;
675 
676  C_FLOAT64 tmp;
677 
678  for (j = 0; j < M; j++)
679  {
680  std::cout << " m " << j << std::endl;
681 
682  for (i = 0; i < N; i++)
683  {
684  tmp = fabs(A(i, j) * F(j, 0) * tsc[M - 1]);
685  std::cout << A(i, j) << " * " << F(j, 0) << " * " << tsc[M - 1] << " = " << tmp << " " << mYerror[i] << std::endl;
686  }
687  }
688 
689 #endif
690 
691 #endif
692 
693 cspiteration:
694 
695  emptyOutputData(N, M, reacs_size);
696 
697 #ifdef CSPDEBUG
698  std::cout << "*********************************** CSP refinement iteration " << iter << "*******************************" << std::endl;
699 #endif
700 
703 
704  A1.resize(N, N);
705  B1.resize(N, N);
706 
707  basisRefinement(N, M, ALA, TAUM, A, B, A1, B1);
708 
709  /* recompute ALA */
710 
711  ALA = 0;
712  TMP = 0;
713 
714  smmult(B1, J, TMP, N, N, N);
715 
716 #if 0
717  DBDT = 0.;
718 
719  if (mTStep)
720  for (i = 0; i < N; i++)
721  for (j = 0; j < N; j++)
722  {
723  DBDT(i, j) = (B1(i, j) - mB(i, j)) / deltaT;
724  }
725 
726 #ifdef CSPDEBUG
727  std::cout << "time derivatives of B " << std::endl;
728  std::cout << DBDT << std::endl;
729 #endif
730 
731  /* TEST: time derivatives are present */
732 
733  if (mTStep)
734  for (i = 0; i < N; i++)
735  for (j = 0; j < N; j++)
736  {
737  TMP(i, j) += DBDT(i, j);
738  }
739 
740 #endif
741 
742  smmult(TMP, A1, ALA, N, N, N);
743 
744 #ifdef CSPDEBUG
745  std::cout << "B1*J*A1 : " << std::endl;
746  std::cout << ALA << std::endl;
747 #endif
748 
749  C_INT result;
750  result = isBlockDiagonal(N, M, ALA, 1.e-8);
751 #ifdef CSPDEBUG
752  std::cout << " iteration result: " << std::endl;
753  std::cout << " 0 - convergence is achieved, 1 - convergence is not achieved, -1 - iterations crucially disconverged " << std::endl;
754  std::cout << result << std::endl;
755 #endif
756 
757  switch (result)
758  {
759  case 0: // iterations converged
760 
761  if (modesAreExhausted(N, M, tsc[M - 1], tsc[M] , g, A1, B1, F))
762  {
763 
764 #ifdef CSPDEBUG
765 
766  std::cout << " the M modes are exhausted, M = " << M << std::endl;
767 #endif
768 
769  //mAmplitude.resize(M);
770 
771  for (j = 0; j < N; j++)
772  //for (j = 0; j < M; j++)
773  mAmplitude[j] = F(j, 0);
774 
775  CSPradicalPointer(N, M, A1, B1);
776  /**
777  * compute CSP Participation Index:
778  * a measure of participation of the r-th elementary reaction to the balancing act of the i-th mode
779  * It is assumed that forward and reverse reactions are counted as distinct
780  *
781  **/
782 
783  CSPParticipationIndex(N, M, tsc[M], B1); //NEW
784 
785  /**
786  * compute CSP Importance Index :
787  * a measure of relative importance of the contribution of r-th elementary reaction
788  * to the current reaction rate of the i-th species
789  *
790  **/
791 
792  /* fast subspace projection matrix */
793  smmult(A1, B1, QF, N, M, N);
794 
795  /* slow subspace projection matrix */
796  smsubst(mI, QF, QSL, N, N);
797 
798  CSPImportanceIndex(N, tsc[M], QSL);
799 
800 #ifdef CSPDEBUG
801  CSPOutput(N, M, reacs_size);
802 #endif
803 
804  //setVectors(M);
805 
806  A = A1;
807  B = B1;
808 
809  return;
810  }
811  else if (M > 1)
812  {
813 
814  M --;
815 
816  if (tsc[M] == tsc[M - 1]) M --;
817 
818 #ifdef CSPDEBUG
819  std::cout << " iterations converged, but the criterion is not sutisfied for the given absolute an relative errors, the new number of candidates to fast " << M << std::endl;
820 #endif
821 
822  if (M) goto analyseMmodes;
823  else
824  {
825 
826  /* No any fast exhausted modes was found on this time step */
828  MCTSSAMethod + 12, mTime);
829 
830  return;
831  }
832  }
833 
834  break;
835 
836  case 1: // next iteration
837 
838  if (iter < mIter)
839  {
840 
841  modesAmplitude(N, M, g, B1, F);
842 //#ifdef CSPDEBUG
843 #if 0
844 
845  std::cout << "scaled amplitudes via refined basis : " << std::endl;
846 
847  for (i = 0; i < M; i++)
848  {
849 
850  std::cout << F(i, 0) << std::endl;
851  }
852 
853  std::cout << " |A(i,m) * F(m,0) * tsc[M - 1]| , mYerror[i] " << std::endl;
854 
855  C_FLOAT64 tmp;
856 
857  for (j = 0; j < M; j++)
858  {
859  std::cout << " m " << j << std::endl;
860 
861  for (i = 0; i < N; i++)
862  {
863  tmp = fabs(A1(i, j) * F(j, 0) * tsc[M - 1]);
864  std::cout << A1(i, j) << " * " << F(j, 0) << " * " << tsc[M - 1] << " = " << tmp << " " << mYerror[i] << std::endl;
865  }
866  }
867 
868 #endif
869  iter ++;
870  A = A1;
871  B = B1;
872 
873  goto cspiteration;
874  }
875  else if (iter >= mIter)
876  {
877  if (M > 1)
878  {
879  M --;
880 
881  if (tsc[M] == tsc[M - 1]) M --;
882 
883  if (M) goto analyseMmodes;
884  else
885  {
886 
887  /* No any fast exhausted modes was found on this time step */
889  MCTSSAMethod + 12, mTime);
890  return;
891  }
892  }
893  }
894 
895  break;
896 
897  case -1: // iterations crucially disconverged
898 
899  if (M > 1)
900  {
901  M --;
902 
903  if (tsc[M] == tsc[M - 1]) M --;
904 
905  if (M) goto analyseMmodes;
906  else
907  {
908 
909  /* No any fast exhausted modes was found on this time step */
910 
912  MCTSSAMethod + 12, mTime);
913  return;
914  }
915  }
916 
917  break;
918 
919  default:
920  /* No any fast exhausted modes was found on this time step */
921  M = 0;
922 
924  MCTSSAMethod + 12, mTime);
925  break;
926  }
927 
928  /* No any fast exhausted modes was found on this time step */
929 
930  M = 0;
931 
933  MCTSSAMethod + 12, mTime);
934 
935  return;
936 }
937 /* compute the norm C of the off-diagonal blocks */
939 {
940  C_INT i, j, imax, jmax, imaxl, jmaxl;
941  C_FLOAT64 max = -1., maxl = -1.;
942 #ifdef CSPDEBUG
943  std::cout << "blocks of ALA : " << std::endl;
944 
945  std::cout << "upper - left : " << std::endl;
946 
947  for (i = 0; i < M; i++)
948  {
949  for (j = 0 ; j < M; j++)
950  std::cout << ALA(i, j) << " ";
951 
952  std::cout << std::endl;
953  }
954 
955  std::cout << "upper - right : " << std::endl;
956 
957  for (i = 0; i < M; i++)
958  {
959  for (j = M ; j < N; j++)
960  std::cout << ALA(i, j) << " ";
961 
962  std::cout << std::endl;
963  }
964 
965  std::cout << "low - left : " << std::endl;
966 
967  for (i = M; i < N; i++)
968  {
969  for (j = 0 ; j < M; j++)
970  std::cout << ALA(i, j) << " ";
971 
972  std::cout << std::endl;
973  }
974 
975  std::cout << "low - right : " << std::endl;
976 
977  for (i = M; i < N; i++)
978  {
979  for (j = M ; j < N; j++)
980  std::cout << ALA(i, j) << " ";
981 
982  std::cout << std::endl;
983  }
984 
985 #endif
986  /* step #1: upper-right block */
987 
988  for (i = 0; i < M; i++)
989  for (j = M ; j < N; j++)
990  if (fabs(ALA(i, j)) > max)
991  {
992  max = fabs(ALA(i, j));
993  imax = i; jmax = j;
994  }
995 
996 #if 1
997  /* step #2: lower-left block */
998 
999  for (i = M; i < N; i++)
1000  for (j = 0 ; j < M; j++)
1001  if (fabs(ALA(i, j)) > maxl)
1002  {
1003  maxl = fabs(ALA(i, j));
1004  imaxl = i ; jmaxl = j;
1005  }
1006 
1007 #ifdef CSPDEBUG
1008  std::cout << "norm C of the lower-left block of ALA is ALA(" << imaxl << "," << jmaxl << ") = " << maxl << std::endl;
1009  std::cout << "the low-left block : " << maxl << std::endl;
1010 #endif
1011 
1012 #endif
1013 
1014 #ifdef CSPDEBUG
1015  std::cout << "the norm C of the upper-right block of ALA is ALA(" << imax << "," << jmax << ") = " << max << std::endl;
1016  std::cout << "the upper-right block : " << max << std::endl;
1017 #endif
1018 
1019  C_INT result;
1020 
1021  result = 1;
1022 
1023  if (fabs(max) >= std::numeric_limits< C_FLOAT64 >::max() || fabs(maxl) >= std::numeric_limits< C_FLOAT64 >::max() || max < 0 || maxl < 0)
1024  {
1025 #ifdef CSPDEBUG
1026  std::cout << " iterations crucially disconverged " << std::endl;
1027 #endif
1028 
1029  result = -1;
1030  }
1031  else if (max <= SMALL) result = 0;
1032 
1033  return result;
1034 }
1036 {
1037 
1038  C_INT i, m, r;
1039  //const CCopasiVector< CReaction > & reacs = mpModel->getReactions();
1040 
1041  for (m = 0; m < M; m++)
1042  for (i = 0; i < N; i++)
1043  {
1044  mAmplitude[i] = 0.;
1045  mRadicalPointer(i, m) = 0;
1046  }
1047 
1048  for (m = 0; m < M; m++)
1049  for (r = 0; r < R; r++)
1050  mFastReactionPointer(r, m) = 0;
1051 
1052  for (m = 0; m < M; m++)
1053  for (r = 0; r < R; r++)
1054  mFastReactionPointerNormed(r, m) = 0;
1055 
1056  for (i = 0; i < N; i++)
1057  for (r = 0; r < R; r++)
1058  mParticipationIndex(r, i) = 0;
1059 
1060  for (i = 0; i < N; i++)
1061  for (r = 0; r < R; r++)
1062  mParticipationIndexNormedRow(r, i) = 0;
1063 
1064  for (i = 0; i < N; i++)
1065  for (r = 0; r < R; r++)
1067 
1068  for (r = 0; r < R; r++)
1069  mFastParticipationIndex[r] = 0;
1070 
1071  for (r = 0; r < R; r++)
1072  mSlowParticipationIndex[r] = 0;
1073 
1074  for (i = 0; i < N; i++)
1075  for (r = 0; r < R; r++)
1076  mImportanceIndex(r, i) = 0;
1077 
1078  for (i = 0; i < N; i++)
1079  for (r = 0; r < R; r++)
1080  mImportanceIndexNormedRow(r, i) = 0;
1081 
1082  return;
1083 }
1084 
1085 void CCSPMethod::step(const double & deltaT)
1086 {
1087  C_INT N = mData.dim;
1088 
1089  C_INT M = 0;
1090 
1093 
1094  A.resize(N, N);
1095  B.resize(N, N);
1096 
1097  C_INT i, j;
1098 
1099  for (i = 0; i < N; i++)
1100  for (j = 0; j < N; j++)
1101  {
1102  A(i, j) = 0;
1103  B(i, j) = 0;
1104  }
1105 
1107 
1108  for (j = 0; j < N; j++)
1109  mG[j] = 0.;
1110 
1111  if (mReducedModel)
1113  else
1115 
1116  if (mReducedModel)
1117  mpModel->calculateJacobianX(mJacobian, 1e-6, 1e-12);
1118  else
1119  mpModel->calculateJacobian(mJacobian, 1e-6, 1e-12);
1120 
1121  cspstep(deltaT, N, M, A, B);
1122 
1123  mB = B;
1124  mTStep = 1;
1125 
1126  if (M > 0)
1127  setVectors(M);
1128  else
1129  setVectorsToNaN();
1130 
1131  setAnnotationM(mCurrentStep); // need for ploting
1132 
1133  mCurrentStep += 1;
1134 
1135  /* integrate one time step */
1136 
1137  integrationStep(deltaT);
1138 
1139  return;
1140 }
1141 
1142 void CCSPMethod::start(const CState * initialState)
1143 {
1144 
1145  mReducedModel = *getValue("Integrate Reduced Model").pBOOL;
1146 
1147  emptyVectors();
1148 
1149  integrationMethodStart(initialState);
1150 
1151  /* CSP related staff */
1152 
1153  mG.resize(mData.dim);
1155  mEps = * getValue("Ratio of Modes Separation").pUDOUBLE;
1156  mRerror = * getValue("Maximum Relative Error").pUDOUBLE;
1157  mAerror = * getValue("Maximum Absolute Error").pUDOUBLE;
1158  mIter = * getValue("Refinement Iterations Number").pUINT;
1159 
1161 
1162  mI.resize(mData.dim, mData.dim);
1163  mB.resize(mData.dim, mData.dim);
1164 
1165  C_INT i, j;
1166 
1167  for (i = 0; i < mData.dim; i++)
1168  for (j = 0; j < mData.dim; j++)
1169  {
1170  mI(i, j) = 0.;
1171  mB(i, j) = 0.;
1172  }
1173 
1174  for (i = 0; i < mData.dim; i++)
1175  mI(i, i) = 1.;
1176 
1177  mTStep = 0;
1178  mCSPbasis = 0;
1179 
1181 
1182  /* CSP Output */
1183 
1184  size_t reacs_size = mpModel->getReactions().size();
1185 
1188  mParticipationIndex.resize(reacs_size, mData.dim);
1189  mImportanceIndex.resize(reacs_size, mData.dim);
1190  mFastReactionPointer.resize(reacs_size, mData.dim);
1191 
1194 
1195  mFastParticipationIndex.resize(reacs_size);
1196  mSlowParticipationIndex.resize(reacs_size);
1197 
1200 
1201  CCopasiVector<CMetab> metabs;
1202  metabs.resize(mData.dim);
1203 
1204 #if 0
1205  CModelEntity* tmp;
1206 
1207  for (j = 0; j < mData.dim; j++)
1208  {
1209 
1211 
1212  CMetab* metab = dynamic_cast< CMetab * >(tmp);
1213 
1214  metabs[j] = metab;
1215  }
1216 
1217  mImportanceIndexTab.resize(reacs_size, mData.dim);
1220  pImportanceIndexAnn->setCopasiVector(1, &metabs);
1221 
1222  mFastParticipationIndexTab.resize(reacs_size, 1);
1225 
1226  mSlowParticipationIndexTab.resize(reacs_size, 1);
1229 
1230  mVec_mFastParticipationIndex[0].resize(reacs_size, 1);
1231  mVec_mSlowParticipationIndex[0].resize(reacs_size, 1);
1232  mVec_mImportanceIndex[0].resize(reacs_size, mData.dim);
1233 
1237 #endif
1238 
1239  mSetVectors = 0;
1240 
1241  return;
1242 }
1243 
1245 {
1246 
1247  C_INT i, m, r;
1248  const CCopasiVector< CReaction > & reacs = mpModel->getReactions();
1249 
1250  std::cout << "Amplitudes of reaction modes :" << std::endl;
1251 
1252  for (m = 0; m < M; m++)
1253  {
1254  std::cout << "reaction mode " << m << " :" << std::endl;
1255 
1256  for (i = 0; i < N; i++)
1257  {
1258  std::cout << " mode " << i << " : " << mAmplitude[i];
1259 
1260  std::cout << std::endl;
1261  }
1262 
1263  std::cout << std::endl;
1264  std::cout << "Radical Pointer: whenever is not a small number, species k is said to be CSP radical" << std::endl;
1265 
1266  for (i = 0; i < N; i++)
1267  std::cout <<
1269  << " : " << mRadicalPointer(i, m) << std::endl;
1270  }
1271 
1272  std::cout << std::endl;
1273  std::cout << " Fast Reaction Pointer of the m-th reaction mode : whenever is not a small number, " << std::endl;
1274  std::cout << " the r-th reaction is said to be a fast reaction " << std::endl;
1275 
1276  for (m = 0; m < M; m++)
1277  {
1278  std::cout << "reaction mode " << m << " :" << std::endl;
1279 
1280  for (r = 0; r < R; r++)
1281  std::cout << reacs[r]->getObjectName() << " :" << mFastReactionPointer(r, m) << std::endl;
1282  }
1283 
1284  std::cout << std::endl;
1285  std::cout << " Participation Index : is a mesure of participation of the r-th elementary reaction " << std::endl;
1286  std::cout << " to the balancing act of the i-th mode " << std::endl;
1287 
1288  for (i = 0; i < N; i++)
1289  {
1290  std::cout << "reaction mode " << i << " :" << std::endl;
1291 
1292  for (r = 0; r < R; r++)
1293  std::cout << reacs[r]->getObjectName() << " :" << mParticipationIndex(r, i) << std::endl;
1294  }
1295 
1296  std::cout << std::endl;
1297  std::cout << " Importance Index: is a mesure of relative importance of the contribution of r-th elementary " << std::endl;
1298  std::cout << " reaction to the current reaction rate of i-th spiecies " << std::endl;
1299 
1300  for (i = 0; i < N; i++)
1301  {
1302  std::cout <<
1304  << " :" << std::endl;
1305 
1306  for (r = 0; r < R; r++)
1307  std::cout << reacs[r]->getObjectName() << " :" << mImportanceIndex(r, i) << std::endl;
1308  }
1309 
1310  return;
1311 }
1312 
1313 /* compute CSP radical pointer and fast reaction pointer */
1315 {
1316 
1317  C_INT i, j, m, r;
1318  C_INT reacs_size = (C_INT) mpModel->getReactions().size();
1319  //const CCopasiVector< CReaction > & reacs = mpModel->getReactions();
1320 
1321  CMatrix< C_FLOAT64 > redStoi;
1322 
1323  if (mReducedModel)
1324  redStoi = mpModel->getRedStoi();
1325  else
1326  redStoi = mpModel->getStoi();
1327 
1328  CMatrix<C_FLOAT64> A0;
1329  CMatrix<C_FLOAT64> B0;
1330  CMatrix<C_FLOAT64> TMP;
1331 
1332  A0.resize(N, reacs_size);
1333  B0.resize(reacs_size, N);
1334  TMP.resize(N, reacs_size);
1335 
1336  A0 = 0;
1337  B0 = 0;
1338 
1339  C_FLOAT64 tmp;
1340 
1341  for (r = 0; r < reacs_size; r++)
1342  {
1343  tmp = 0;
1344 
1345  for (i = 0; i < N; i++)
1346  {
1347  A0(i, r) = redStoi(i, r);
1348  tmp += A0(i, r) * A0(i, r);
1349  }
1350 
1351  for (i = 0; i < N; i++)
1352  B0(r, i) = A0(i, r) / tmp;
1353  }
1354 
1355  /* m-th fast mode projection matrix */
1356 
1357  CMatrix<C_FLOAT64> QM;
1358  QM.resize(N, N);
1359 
1360  QM = 0.;
1361 
1362  for (m = 0; m < M ; m++)
1363  {
1364 
1365  CMatrix<C_FLOAT64> Qm;
1366  Qm.resize(N, N);
1367 
1368  for (i = 0; i < N; ++i)
1369  for (j = 0; j < N; ++j)
1370  {
1371  Qm(i, j) = A(i, m) * B(m, j);
1372  }
1373 
1374  /**
1375  * some comments on the Qm matrix: Qm(i,i) , i = 0,1,...,N,
1376  * is a measure of projection of i-th unit vector in the m-th mode,
1377  * whenever Qm(i,i) is not a small number, species m is said to be a CSP radical
1378  **/
1379 
1380  C_FLOAT64 tmp = 0.;
1381 
1382  for (i = 0; i < N ; i++)
1383  {
1384  mRadicalPointer(i, m) = Qm(i, i);
1385 
1386  tmp += Qm(i, i);
1387  }
1388 
1389  /* use stoichiometric vectors to build the fast reaction pointer */
1390 
1391  /**
1392  * Pmr is a measure of projection of r-th stoichiometric vector in the m-th mode,
1393  * whenever Pmr is not a small number, the r-th reaction is said to be a fast reaction
1394  **/
1395 
1396  for (r = 0; r < reacs_size; r++)
1397  {
1398  C_FLOAT64 Pmr = 0.;
1399 
1400  for (i = 0; i < N; i++)
1401  {
1402  TMP(i, r) = 0.;
1403 
1404  for (j = 0; j < N; j++)
1405  TMP(i, r) += Qm(i, j) * A0(j, r);
1406  }
1407 
1408  for (j = 0; j < N; j++)
1409  Pmr += B0(r, j) * TMP(j, r);
1410 
1411  mFastReactionPointer(r, m) = Pmr;
1412  }
1413  }
1414 
1415  /* the Fast Reaction Pointer normed by column */
1416 
1417  for (m = 0; m < M ; m++)
1418  {
1419  C_FLOAT64 sum = 0.;
1420 
1421  for (r = 0; r < reacs_size; r++)
1422  sum += fabs(mFastReactionPointer(r, m));
1423 
1424  for (r = 0; r < reacs_size; r++)
1425  mFastReactionPointerNormed(r, m) = 100.*mFastReactionPointer(r, m) / sum;
1426  }
1427 
1428  return;
1429 }
1430 
1431 /**
1432  * compute CSP Participation Index:
1433  * a measure of participation of the r-th elementary reaction to the balancing act of the i-th mode
1434  * It is assumed that forward and reverse reactions are counted as distinct
1435  **/
1437 {
1438 
1439  C_INT i, r, j;
1440  C_INT reacs_size = (C_INT) mpModel->getReactions().size();
1441  const CCopasiVector< CReaction > & reacs = mpModel->getReactions();
1442 
1443  CMatrix< C_FLOAT64 > redStoi;
1444 
1445  if (mReducedModel)
1446  redStoi = mpModel->getRedStoi();
1447  else
1448  redStoi = mpModel->getStoi();
1449 
1450  CVector<C_FLOAT64> flux;
1451  flux.resize(reacs_size);
1452 
1454  P.resize(N, reacs_size);
1455 
1456  CVector<C_FLOAT64> estim;
1457  estim.resize(N);
1458 
1459  CVector<C_FLOAT64> ampl;
1460  ampl.resize(N);
1461 
1462  for (r = 0; r < reacs_size; ++r)
1463  flux[r] = reacs[r]->calculateParticleFlux();
1464 
1465  for (i = 0; i < N; ++i)
1466  {
1467  ampl[i] = 0;
1468 
1469  for (r = 0; r < reacs_size; ++r)
1470  {
1471 
1472  P(i, r) = 0;
1473 
1474  for (j = 0; j < N; ++j)
1475  P(i, r) += B0(i, j) * redStoi(j, r);
1476 
1477  ampl[i] += fabs(P(i, r) * flux[r]);
1478  }
1479 
1480  C_FLOAT64 tmp = 0.0;
1481 
1482  for (j = 0; j < N; ++j)
1483  {
1484  tmp += B0(i, j) * mYerror[j];
1485  }
1486 
1487  estim[i] = fabs(tmp / tauM1);
1488  }
1489 
1490  for (i = 0; i < N; ++i)
1491  {
1492 
1493  for (r = 0; r < reacs_size; ++r)
1494  {
1495  P(i, r) *= flux[r] / (ampl[i] + estim[i]);
1496 
1497  mParticipationIndex(r, i) = P(i, r);
1498  }
1499  }
1500 
1501  for (r = 0; r < reacs_size; ++r)
1502  {
1503  C_FLOAT64 PI = 0.;
1504 
1505  for (i = 0; i < M; ++i)
1506  PI += fabs(mParticipationIndex(r, i));
1507 
1509 
1510  for (i = M; i < N; ++i)
1511  PI += fabs(mParticipationIndex(r, i));
1512 
1514  // TEST: mFastParticipationIndex[r] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
1515 
1516  for (i = M; i < N; ++i)
1517  mSlowParticipationIndex[r] += fabs(mParticipationIndex(r, i));
1518 
1520  }
1521 
1522  /* the Participation Index normed by column */
1523 
1524  for (i = 0; i < N; ++i)
1525  {
1526  C_FLOAT64 sum = 0.;
1527 
1528  for (r = 0; r < reacs_size; ++r)
1529  sum += fabs(mParticipationIndex(r, i));
1530 
1531  for (r = 0; r < reacs_size; ++r)
1532  mParticipationIndexNormedColumn(r, i) = 100.* mParticipationIndex(r, i) / sum;
1533  }
1534 
1535  /* the Participation Index normed by row */
1536 
1537  for (r = 0; r < reacs_size; ++r)
1538  {
1539  C_FLOAT64 sum = 0.;
1540 
1541  for (i = 0; i < N; ++i)
1542  sum += fabs(mParticipationIndex(r, i));
1543 
1544  for (i = 0; i < N; ++i)
1545  mParticipationIndexNormedRow(r, i) = 100.* mParticipationIndex(r, i) / sum;
1546  }
1547 
1548  return;
1549 }
1550 
1551 /**
1552  * compute CSP Importance Index :
1553  * a measure of relative importance of the contribution of r-th elementary reaction
1554  * to the current reaction rate of the i-th species
1555  **/
1556 
1558 {
1559 
1560  C_INT i, r;
1561  C_INT reacs_size = (C_INT) mpModel->getReactions().size();
1562  const CCopasiVector< CReaction > & reacs = mpModel->getReactions();
1563 
1564  CMatrix< C_FLOAT64 > redStoi;
1565 
1566  if (mReducedModel)
1567  redStoi = mpModel->getRedStoi();
1568  else
1569  redStoi = mpModel->getStoi();
1570 
1571  CVector<C_FLOAT64> flux;
1572  flux.resize(reacs_size);
1573 
1574  CMatrix<C_FLOAT64> S0;
1575  S0.resize(N, reacs_size);
1576 
1578  S.resize(N, reacs_size);
1579 
1581  I.resize(N, reacs_size);
1582 
1583  CVector<C_FLOAT64> estim;
1584  estim.resize(N);
1585 
1587  g.resize(N);
1588 
1589  S = redStoi;
1590 
1591  smmult(Q, S, S0, N, N, reacs_size);
1592 
1593  for (r = 0; r < reacs_size; ++r)
1594  flux[r] = reacs[r]->calculateParticleFlux();
1595 
1596  for (i = 0; i < N; ++i)
1597  {
1598  g[i] = 0;
1599 
1600  for (r = 0; r < reacs_size; ++r)
1601  g[i] += fabs(S0(i, r) * flux[r]);
1602 
1603  estim[i] = fabs(mYerror[i] / tauM1);
1604  }
1605 
1606  for (i = 0; i < N; ++i)
1607  {
1608 
1609  for (r = 0; r < reacs_size; ++r)
1610  {
1611  I(i, r) = S0(i, r) * flux[r] / (g[i] + estim[i]);
1612 
1613  mImportanceIndex(r, i) = I(i, r);
1614  }
1615  }
1616 
1617  /* the Importance Index normed by row */
1618 
1619  for (r = 0; r < reacs_size; ++r)
1620  {
1621  C_FLOAT64 sum = 0.;
1622 
1623  for (i = 0; i < N; ++i)
1624  sum += fabs(mImportanceIndex(r, i));
1625 
1626  for (i = 0; i < N; ++i)
1627  mImportanceIndexNormedRow(r, i) = 100.* mImportanceIndex(r, i) / sum;
1628  }
1629 
1630  return;
1631 }
1632 
1633 /* compute amplitudes of fast and slow modes */
1635 {
1636 
1637  C_INT i, j;
1638 
1639  /* evaluate amplitudes */
1640 
1641  for (i = 0; i < N; i++)
1642  {
1643  F(i, 0) = 0.;
1644 
1645  for (j = 0; j < N; j++)
1646 
1647  for (j = 0; j < N; j++)
1648  {
1649  F(i, 0) += B(i, j) * g[j];
1650  }
1651  }
1652 
1653  return;
1654 }
1655 
1656 /* correct for the contribution of the fast time-scales to y */
1658 {
1659  CMatrix<C_FLOAT64> TMP;
1660  CVector<C_FLOAT64> dy;
1661 
1662  TMP.resize(N, M);
1663  dy.resize(N);
1664 
1665  C_INT i, k;
1666 
1667  /* A*TAU*F */
1668 
1669  smmult(A, TAUM, TMP, N, M, M);
1670 
1671  for (i = 0; i < N ; i++) dy[i] = 0.;
1672 
1673  for (i = 0; i < N ; i++)
1674  for (k = 0; k < M; k++)
1675  dy[i] += TMP(i, k) * F(k, 0);
1676 
1677  /* contribution of fast time-scales */
1678 
1679  for (i = 0; i < N; i++)
1680  {
1681  y[i] -= dy[i];
1682  }
1683 
1684  return;
1685 }
1686 
1687 /* Refinement Procedure :
1688  * Lamm, Combustion Science and Technology, 1993.
1689  **/
1691 {
1692 
1693  C_INT i, j, n, m;
1694 
1697 
1698  P.resize(N, N);
1699  Q.resize(N, N);
1700 
1701  P = 0.;
1702  Q = 0.;
1703 
1704  for (m = 0; m < M; m++)
1705  for (j = M; j < N; j++)
1706  {
1707  for (n = 0; n < M; n++)
1708  P(m, j) += TAU(m, n) * ALA(n, j);
1709  }
1710 
1711  for (j = M ; j < N; j++)
1712  for (m = 0; m < M; m++)
1713  {
1714  for (n = 0; n < M; n++)
1715  Q(j, m) += ALA(j, n) * TAU(n, m);
1716  }
1717 
1718  A0 = A;
1719  B0 = B;
1720 
1721  /* step #1 */
1722 
1723  for (m = 0; m < M; m++)
1724  for (i = 0; i < N; i++)
1725  for (j = M ; j < N; j++)
1726  B0(m, i) += P(m, j) * B(j, i);
1727 
1728  for (i = 0; i < N; i++)
1729  for (j = M ; j < N; j++)
1730  for (n = 0; n < M ; n++)
1731  A0(i, j) -= A(i, n) * P(n, j);
1732 
1733 #if 1
1734  /* step #2 */
1735 
1736  A = A0;
1737  B = B0;
1738 
1739  for (i = M; i < N; i++)
1740  for (j = 0; j < N; j++)
1741  for (n = 0; n < M; n++)
1742  B0(i, j) -= Q(i, n) * B(n, j);
1743 
1744  for (i = 0; i < N; i++)
1745  for (m = 0; m < M; m++)
1746  for (j = M; j < N; j++)
1747  A0(i, m) += A(i, j) * Q(j, m);
1748 
1749 #endif
1750 
1751  //smnorm(N, A0, B0, N);
1752 
1753  return;
1754 }
1755 
1756 /* "true" if each of the analyzed M modes is exhausted */
1758 {
1759  C_INT i, j;
1760  C_FLOAT64 tmp;
1761 
1762  bool exhausted = true;
1763 
1764  modesAmplitude(N, M, g, B, F);
1765 
1766 #ifdef CSPDEBUG
1767  std::cout << " |A(i,m) * F(m,0) * tsc[M - 1]| , mYerror[i] " << std::endl;
1768 #endif
1769 
1770  for (j = 0; j < M; j++)
1771  {
1772 
1773  for (i = 0; i < N; i++)
1774  {
1775  tmp = fabs(A(i, j) * F(j, 0) * tauM);
1776 
1777 #ifdef CSPDEBUG
1778  std::cout << A(i, j) << " * " << F(j, 0) << " * " << tauM << " = " << tmp << " " << mYerror[i] << std::endl;
1779 #endif
1780 
1781  if (tmp >= mYerror[i]) exhausted = false;
1782  }
1783  }
1784 
1785  return exhausted;
1786 }
1787 
1788 /**
1789  * Predifine the CArrayAnnotation for plots
1790  **/
1791 
1793 {
1794  mReducedModel = *getValue("Integrate Reduced Model").pBOOL;
1795 
1796  C_INT N;
1797 
1798  if (mReducedModel)
1799  {
1801  }
1802  else
1803  {
1805  }
1806 
1807  CCopasiVector<CMetab> metabs;
1808  metabs.resize(N);
1809 
1810  CModelEntity* tmp;
1811  C_INT j;
1812 
1813  for (j = 0; j < N; j++)
1814  {
1815 
1817 
1818  CMetab* metab = dynamic_cast< CMetab * >(tmp);
1819 
1820  metabs[j] = metab;
1821  }
1822 
1824  N);
1825 
1826  pImportanceIndexAnn->resize(); // fixed
1828  pImportanceIndexAnn->setCopasiVector(1, &metabs);
1829 
1832 
1834 
1837 
1839 }
1840 
1841 /**
1842  * Create the CArraAnnotations for tables in the CQTSSAResultSubWidget.
1843  * Input for each CArraAnnotations is a separate CMatrix.
1844  **/
1846 {
1847  tableNames.erase(tableNames.begin(), tableNames.end());
1848 
1849  std::string name;
1850 
1851  /* this table is not visible more */
1852 #if 0
1853 
1855  pTmp1 = new CArrayAnnotation("Amplitude", this,
1857  pTmp1->setMode(0, pTmp1->STRINGS);
1858  pTmp1->setMode(1, pTmp1->STRINGS);
1859  pTmp1->setDescription(" ");
1860  pTmp1->setDimensionDescription(0, "Fast Reaction Modes");
1861  pTmp1->setDimensionDescription(1, "Amplitudes ");
1862  pAmplitudeAnn = pTmp1;
1863 
1864 #endif
1865 
1866  name = "Radical Pointer";
1867  tableNames.push_back(name);
1868 
1870  pTmp2 = new CArrayAnnotation("Radical Pointer", this,
1872  pTmp2->setMode(0, pTmp2->VECTOR);
1873  pTmp2->setMode(1, pTmp2->STRINGS);
1874  pTmp2->setDescription("Radical Pointer: whenever is not a small number, species k is said to be CSP radical ");
1875  pTmp2->setDimensionDescription(0, "Species");
1876  pTmp2->setDimensionDescription(1, "Fast Time Scales");
1878 
1880 
1881  name = "Fast Reaction Pointer";
1882  tableNames.push_back(name);
1883 
1885  pTmp3 = new CArrayAnnotation("Fast Reaction Pointer", this,
1887  pTmp3->setMode(0, pTmp3->VECTOR);
1888  pTmp3->setMode(1, pTmp3->STRINGS);
1889  pTmp3->setDescription("Fast Reaction Pointer of the m-th reaction mode : whenever is not a small number, the r-th reaction is said to be a fast reaction");
1890  pTmp3->setDimensionDescription(0, "Reactions");
1891  pTmp3->setDimensionDescription(1, "Fast Time Scales");
1893 
1895 
1896  name = "Normed Fast Reaction Pointer";
1897  tableNames.push_back(name);
1898 
1900  pTmp3Normed = new CArrayAnnotation("Normed Fast Reaction Pointer", this,
1902  pTmp3Normed->setMode(0, pTmp3Normed->VECTOR);
1903  pTmp3Normed->setMode(1, pTmp3Normed->STRINGS);
1904  pTmp3Normed->setDescription("Fast Reaction Pointer of the m-th reaction mode : whenever is not a small number, the r-th reaction is said to be a fast reaction");
1905  pTmp3Normed->setDimensionDescription(0, "Reactions");
1906  pTmp3Normed->setDimensionDescription(1, "Fast Time Scales");
1908 
1910 
1911  name = "Participation Index";
1912  tableNames.push_back(name);
1913 
1915  pTmp4 = new CArrayAnnotation("Participation Index", this,
1917  pTmp4->setMode(1, pTmp4->STRINGS);
1918  pTmp4->setMode(0, pTmp4->VECTOR);
1919  pTmp4->setDescription("Participation Index : is a measure of participation of the r-th elementary reaction to the balancing act of the i-th mode");
1920  pTmp4->setDimensionDescription(0, "Reactions");
1921  pTmp4->setDimensionDescription(1, "Time Scales");
1923 
1925 
1926  name = "Normed Participation Index (by column)";
1927  tableNames.push_back(name);
1928 
1930  pTmp4NormedColumn = new CArrayAnnotation("Normed Participation Index (by column)", this,
1932  pTmp4NormedColumn->setMode(1, pTmp4NormedColumn->STRINGS);
1933  pTmp4NormedColumn->setMode(0, pTmp4NormedColumn->VECTOR);
1934  pTmp4NormedColumn->setDescription("Participation Index : is a measure of participation of the r-th elementary reaction to the balancing act of the i-th mode");
1935  pTmp4NormedColumn->setDimensionDescription(0, "Reactions");
1936  pTmp4NormedColumn->setDimensionDescription(1, "Time Scales");
1938 
1940 
1941  name = "Normed Participation Index (by row)";
1942  tableNames.push_back(name);
1943 
1945  pTmp4NormedRow = new CArrayAnnotation("Normed Participation Index (by row)", this,
1947  pTmp4NormedRow->setMode(1, pTmp4NormedRow->STRINGS);
1948  pTmp4NormedRow->setMode(0, pTmp4NormedRow->VECTOR);
1949  pTmp4NormedRow->setDescription("Participation Index : is a measure of participation of the r-th elementary reaction to the balancing act of the i-th mode");
1950  pTmp4NormedRow->setDimensionDescription(0, "Reactions");
1951  pTmp4NormedRow->setDimensionDescription(1, "Time Scales");
1953 
1955 
1956  name = "Fast Participation Index";
1957  tableNames.push_back(name);
1958 
1960  pTmp4Fast = new CArrayAnnotation("Fast Participation Index", this,
1962  pTmp4Fast->setMode(0, pTmp4Fast->VECTOR);
1963  pTmp4Fast->setMode(1, pTmp4Fast->STRINGS);
1964  pTmp4Fast->setDescription(" Fast Participation Index : is a measure of participation of the r-th elementary reaction to the balancing act of fast modes");
1965  pTmp4Fast->setDimensionDescription(0, "Reactions");
1966  pTmp4Fast->setDimensionDescription(1, " ");
1968 
1970 
1971  name = "Slow Participation Index";
1972  tableNames.push_back(name);
1973 
1975  pTmp4Slow = new CArrayAnnotation("Slow Participation Index", this,
1977 
1978  pTmp4Slow->setMode(0, pTmp4Slow->VECTOR);
1979  pTmp4Slow->setMode(1, pTmp4Slow->STRINGS);
1980  pTmp4Slow->setDescription("Slow Participation Index : is a measure of participation of the r-th elementary reaction to the balancing act of slow modes");
1981  pTmp4Slow->setDimensionDescription(0, "Reactions");
1982  pTmp4Slow->setDimensionDescription(1, " ");
1984 
1986 
1987  name = "Importance Index";
1988  tableNames.push_back(name);
1989 
1991  pTmp5 = new CArrayAnnotation("Importance Index", this,
1993  pTmp5->setMode(1, pTmp5->VECTOR);
1994  pTmp5->setMode(0, pTmp5->VECTOR);
1995  pTmp5->setDescription("Importance Index: is a measure of relative importance of the contribution of r-th elementary reaction to the current reaction rate of i-th species");
1996  pTmp5->setDimensionDescription(0, "Reactions");
1997  pTmp5->setDimensionDescription(1, "Species");
1999 
2001 
2002  name = "Normed Importance Index (by row)";
2003  tableNames.push_back(name);
2004 
2006  pTmp5NormedRow = new CArrayAnnotation("Normed Importance Index (by row)", this,
2008  pTmp5NormedRow->setMode(1, pTmp5NormedRow->VECTOR);
2009  pTmp5NormedRow->setMode(0, pTmp5NormedRow->VECTOR);
2010  pTmp5NormedRow->setDescription("Importance Index: is a measure of relative importance of the contribution of r-th elementary reaction to the current reaction rate of i-th species");
2011  pTmp5NormedRow->setDimensionDescription(0, "Reactions");
2012  pTmp5NormedRow->setDimensionDescription(1, "Species");
2014 
2016 }
2017 
2018 /**
2019  * Set the every CArrayAnnotation for the requested step.
2020  * Set also the description of CArayAnnotation for both dimensions:
2021  * - dimension description could consists of some std::strings
2022  * some strings contain the Time Scale values for requested step
2023  * - dimension description could consists of arrays of CommonNames
2024  **/
2026 {
2027  std::string str;
2028  std::stringstream sstr;
2029  sstr.str("");
2030  sstr.clear();
2031  unsigned C_INT32 i;
2032  double timeScale;
2033  unsigned C_INT32 M;
2034 
2035  C_INT N;
2036 
2037  if (mReducedModel)
2038  {
2040  }
2041  else
2042  {
2044  }
2045 
2046  //TEST : if (step == 0) return false;
2047 
2048  //if (mVec_SlowModes.size() == 0) return false;
2049 
2050  //if (step > mVec_SlowModes.size()) return false;
2051 
2052  // TEST : step -= 1;
2053 
2054  if (mVec_SlowModes.size() <= step)
2055  {
2056  return false;
2057  }
2058 
2059  M = mVec_SlowModes[step];
2060 
2061  // fill pAmplitudeAnn
2062 
2063 #if 0
2064  mAmplitudeTab.resize(mVec_mAmplitude[step].numRows(), 1);
2065 
2066  for (i = 0; i < mVec_mAmplitude[step].numRows(); i++)
2067  mAmplitudeTab(i, 0) = mVec_mAmplitude[step][i][0];
2068 
2069  pAmplitudeAnn->resize();
2070 
2071  for (i = 0; i < mVec_mAmplitude[step].numRows(); i++)
2072  {
2073  timeScale = mVec_TimeScale[step][i];
2074 
2075  if (i < M)
2076  sstr << "Fast: ";
2077  else
2078  sstr << "Slow: ";
2079 
2080  sstr << timeScale;
2081  str = sstr.str();
2082  pAmplitudeAnn->setAnnotationString(0, i, str);
2083  sstr.str("");
2084  sstr.clear();
2085  }
2086 
2087  sstr.str("");
2088  str = sstr.str();
2089  pAmplitudeAnn->setAnnotationString(1, 0, str);
2090 
2091 #endif
2092 
2093  // fill pRadicalPointerAnn
2095  mVec_mRadicalPointer[step].numRows());
2098 
2099  CCopasiVector<CMetab> metabs;
2100  //FIXED : metabs.resize(mData.dim);
2101  metabs.resize(N);
2102 
2103  CModelEntity* tmp;
2104  C_INT j;
2105 
2106  // FIXED: for (j = 0; j < mData.dim; j++)
2107  for (j = 0; j < N; j++)
2108  {
2109 
2111 
2112  CMetab* metab = dynamic_cast< CMetab * >(tmp);
2113 
2114  metabs[j] = metab;
2115  }
2116 
2117  pRadicalPointerAnn->setCopasiVector(0, &metabs);
2118 
2119 #if 0
2120  std::cout << "metab.size " << mData.dim << std::endl;
2121  std::cout << "step " << step << std::endl;
2122  std::cout << "mVec_mRadicalPointer[step].numCols() " << mVec_mRadicalPointer[step].numCols() << std::endl;
2123  std::cout << "mVec_mRadicalPointer[step].numRows() " << mVec_mRadicalPointer[step].numRows() << std::endl;
2124  std::cout << "+++++++++++++++++++++++++++++++++++++++++++ " << std::endl;
2125 
2126 #endif
2127 
2128  for (i = 0; i < mVec_mRadicalPointer[step].numCols(); i++)
2129  {
2130  timeScale = mVec_TimeScale[step][i];
2131 
2132  if (i < M)
2133  sstr << "Fast: ";
2134  else
2135  sstr << "Slow: ";
2136 
2137  sstr << timeScale;
2138  str = sstr.str();
2140  sstr.str("");
2141  sstr.clear();
2142  }
2143 
2144  // fill pFastReactionPointerAnn
2146  mVec_mFastReactionPointer[step].numRows());
2149 
2150  for (i = 0; i < mVec_mFastReactionPointer[step].numCols(); i++)
2151  {
2152  timeScale = mVec_TimeScale[step][i];
2153 
2154  if (i < M)
2155  sstr << "Fast: ";
2156  else
2157  sstr << "Slow: ";
2158 
2159  sstr << timeScale;
2160  str = sstr.str();
2162  sstr.str("");
2163  sstr.clear();
2164  }
2165 
2167 
2168  // fill pFastReactionPointerNormedAnn
2170  mVec_mFastReactionPointerNormed[step].numRows());
2173 
2174  for (i = 0; i < mVec_mFastReactionPointerNormed[step].numCols(); i++)
2175  {
2176  timeScale = mVec_TimeScale[step][i];
2177 
2178  if (i < M)
2179  sstr << "Fast: ";
2180  else
2181  sstr << "Slow: ";
2182 
2183  sstr << timeScale;
2184  str = sstr.str();
2186  sstr.str("");
2187  sstr.clear();
2188  }
2189 
2191 
2192  // fill pParticipationIndexAnn
2194  mVec_mParticipationIndex[step].numRows());
2197 
2198  for (i = 0; i < mVec_mParticipationIndex[step].numCols(); i++)
2199  {
2200  timeScale = mVec_TimeScale[step][i];
2201 
2202  if (i < M)
2203  sstr << "Fast: ";
2204  else
2205  sstr << "Slow: ";
2206 
2207  sstr << timeScale;
2208  str = sstr.str();
2210  sstr.str("");
2211  sstr.clear();
2212  }
2213 
2215 
2216  // fill pParticipationIndexNormedColumnAnn
2218  mVec_mParticipationIndexNormedColumn[step].numRows());
2221 
2222  for (i = 0; i < mVec_mParticipationIndexNormedColumn[step].numCols(); i++)
2223  {
2224  timeScale = mVec_TimeScale[step][i];
2225 
2226  if (i < M)
2227  sstr << "Fast: ";
2228  else
2229  sstr << "Slow: ";
2230 
2231  sstr << timeScale;
2232  str = sstr.str();
2234  sstr.str("");
2235  sstr.clear();
2236  }
2237 
2239 
2240  // fill pParticipationIndexNormedRowAnn
2242  mVec_mParticipationIndexNormedRow[step].numRows());
2245 
2246  for (i = 0; i < mVec_mParticipationIndexNormedRow[step].numCols(); i++)
2247  {
2248  timeScale = mVec_TimeScale[step][i];
2249 
2250  if (i < M)
2251  sstr << "Fast: ";
2252  else
2253  sstr << "Slow: ";
2254 
2255  sstr << timeScale;
2256  str = sstr.str();
2258  sstr.str("");
2259  sstr.clear();
2260  }
2261 
2263 
2264  // fill pFastParticipationIndexAnn
2265 
2269 
2270  sstr << " ";
2271  str = sstr.str();
2274 
2275  // fill pSlowParticipationIndexAnn
2276 
2277  //mSlowParticipationIndexTab.resize(mpModel->getReactions().size(),1);
2279 
2280  //pSlowParticipationIndexAnn->resize();
2281 
2282  //sstr << " ";
2283  //str = sstr.str();
2284  //pSlowParticipationIndexAnn->setAnnotationString(1, 0, str);
2285  //pSlowParticipationIndexAnn->setCopasiVector(0, &mpModel->getReactions());
2286 
2287  // fill pmImportanceIndexAnn
2289  mVec_mImportanceIndex[step].numRows());
2291  pImportanceIndexAnn->resize(); // fixed
2293  pImportanceIndexAnn->setCopasiVector(1, &metabs);
2294 
2295 // fill pmImportanceIndexNormedRowAnn
2297  mVec_mImportanceIndexNormedRow[step].numRows());
2302 
2303  return true;
2304 }
2305 
2306 /**
2307  * set vectors to NaN when the reduction was not possible
2308  **/
2310 {
2311 
2312  mVec_TimeScale.push_back(mCurrentStep);
2314  C_INT i, r, m, fast;
2315  C_INT reacs_size = (C_INT) mpModel->getReactions().size();
2316 
2317  fast = mData.dim;
2318 
2319  for (i = 0; i < mData.dim; i++)
2320  mVec_TimeScale[mCurrentStep][i] = -1 / mR(i, i);
2321 
2322  mVec_SlowModes.push_back(mCurrentStep);
2324 
2326 
2327  mVec_mRadicalPointer[mCurrentStep].resize(mData.dim, fast);
2328 
2329  for (m = 0; m < fast; m++)
2330  for (i = 0; i < mData.dim; i++)
2331  mVec_mRadicalPointer[mCurrentStep][i][m] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2332 
2334 
2335  mVec_mFastReactionPointer[mCurrentStep].resize(reacs_size, fast);
2336 
2337  for (r = 0; r < reacs_size; r++)
2338  for (i = 0; i < fast; i++)
2339  mVec_mFastReactionPointer[mCurrentStep][r][i] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2340 
2342 
2343  mVec_mFastReactionPointerNormed[mCurrentStep].resize(reacs_size, fast);
2344 
2345  for (r = 0; r < reacs_size; r++)
2346  for (i = 0; i < fast; i++)
2347  mVec_mFastReactionPointerNormed[mCurrentStep][r][i] = std::numeric_limits<C_FLOAT64>::quiet_NaN();;
2348 
2350  mVec_mParticipationIndex[mCurrentStep].resize(reacs_size, mData.dim);
2351 
2352  for (r = 0; r < reacs_size; r++)
2353  for (i = 0; i < fast; i++)
2354  mVec_mParticipationIndex[mCurrentStep][r][i] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2355 
2357  mVec_mFastParticipationIndex[mCurrentStep].resize(reacs_size, 1);
2358 
2359  for (i = 0; i < reacs_size; i++)
2360  mVec_mFastParticipationIndex[mCurrentStep][i][0] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2361 
2363  mVec_mSlowParticipationIndex[mCurrentStep].resize(reacs_size, 1);
2364 
2365  for (i = 0; i < reacs_size; i++)
2366  mVec_mSlowParticipationIndex[mCurrentStep][i][0] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2367 
2370 
2371  for (r = 0; r < reacs_size; r++)
2372  for (i = 0; i < fast; i++)
2373  mVec_mParticipationIndexNormedColumn[mCurrentStep][r][i] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2374 
2377 
2378  for (r = 0; r < reacs_size; r++)
2379  for (i = 0; i < fast; i++)
2380  mVec_mParticipationIndexNormedRow[mCurrentStep][r][i] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2381 
2383  mVec_mImportanceIndex[mCurrentStep].resize(reacs_size, mData.dim);
2384 
2385  for (r = 0; r < reacs_size; r++)
2386  for (i = 0; i < mData.dim; i++)
2387  mVec_mImportanceIndex[mCurrentStep][r][i] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2388 
2391 
2392  for (r = 0; r < reacs_size; r++)
2393  for (i = 0; i < mData.dim; i++)
2394  mVec_mImportanceIndexNormedRow[mCurrentStep][r][i] = std::numeric_limits<C_FLOAT64>::quiet_NaN();
2395 
2396  mCurrentTime.push_back(mCurrentStep);
2398 }
2399 /**
2400  *upgrade all vectors with values from actually calculation for current step
2401  **/
2403 {
2404 
2405  mVec_TimeScale.push_back(mCurrentStep);
2407  C_INT i, r, m;
2408  C_INT reacs_size = (C_INT) mpModel->getReactions().size();
2409 
2410  for (i = 0; i < mData.dim; i++)
2411  mVec_TimeScale[mCurrentStep][i] = -1 / mR(i, i);
2412 
2413  mVec_SlowModes.push_back(mCurrentStep);
2414  mVec_SlowModes[mCurrentStep] = fast;
2415 
2416 #if 0
2417  mVec_mAmplitude.push_back(mCurrentStep);
2418  //mVec_mAmplitude[mCurrentStep].resize(mAmplitude.size(), 1);
2419  mVec_mAmplitude[mCurrentStep].resize(fast, 1);
2420 
2421  for (i = 0; i < fast; i++)
2423 
2424 #endif
2425 
2427  //mVec_mRadicalPointer[mCurrentStep].resize(mRadicalPointer.numCols(), mRadicalPointer.numRows());
2428  //mVec_mRadicalPointer[mCurrentStep] = mRadicalPointer;
2429 
2430  mVec_mRadicalPointer[mCurrentStep].resize(mData.dim, fast);
2431 
2432  for (m = 0; m < fast; m++)
2433  for (i = 0; i < mData.dim; i++)
2435 
2437  //mVec_mFastReactionPointer[mCurrentStep].resize(mFastReactionPointer.numCols(), mFastReactionPointer.numRows());
2438  //mVec_mFastReactionPointer[mCurrentStep] = mFastReactionPointer;
2439 
2440  mVec_mFastReactionPointer[mCurrentStep].resize(reacs_size, fast);
2441 
2442  for (r = 0; r < reacs_size; r++)
2443  for (i = 0; i < fast; i++)
2445 
2447  //mVec_mFastReactionPointerNormed[mCurrentStep].resize(mFastReactionPointerNormed.numCols(), mFastReactionPointerNormed.numRows());
2448  //mVec_mFastReactionPointerNormed[mCurrentStep] = mFastReactionPointerNormed;
2449 
2450  mVec_mFastReactionPointerNormed[mCurrentStep].resize(reacs_size, fast);
2451 
2452  for (r = 0; r < reacs_size; r++)
2453  for (i = 0; i < fast; i++)
2455 
2459 
2462 
2463  for (i = 0; i < reacs_size; i++)
2465 
2468 
2469  for (i = 0; i < reacs_size; i++)
2471 
2475 
2479 
2483 
2487 
2488  mCurrentTime.push_back(mCurrentStep);
2490 }
2491 
2492 /**
2493  * Empty every vector to be able to fill them with new values for a new calculation.
2494  * Also nullify the step counter.
2495  **/
2497 {
2498 
2499  mCurrentStep = 0;
2500  mVec_TimeScale.erase(mVec_TimeScale.begin(), mVec_TimeScale.end());
2501  mVec_SlowModes.erase(mVec_SlowModes.begin(), mVec_SlowModes.end());
2502 
2503 #if 0
2504  mVec_mAmplitude.erase(mVec_mAmplitude.begin(), mVec_mAmplitude.end());
2505 #endif
2506 
2517 }
2518 
2519 void CCSPMethod::printResult(std::ostream * ostream) const
2520 {
2521  std::ostream & os = *ostream;
2522  C_INT M, i, m, r, istep = 0;
2523 
2524  C_INT32 stepNumber;
2525  //double timeScale;
2526 
2528  CTSSATask* pTask =
2529  dynamic_cast<CTSSATask *>((*(*CCopasiRootContainer::getDatamodelList())[0]->getTaskList())["Time Scale Separation Analysis"]);
2530 
2531  CTSSAProblem* pProblem = dynamic_cast<CTSSAProblem*>(pTask->getProblem());
2532 
2533  stepNumber = (int)mVec_SlowModes.size();
2534 
2535  this->print(&os);
2536 
2537  const CCopasiVector< CReaction > & reacs = mpModel->getReactions();
2538 
2539  os << std::endl;
2540  os << " Radical Pointer: whenever is not a small number, species k is said to be CSP radical" << std::endl;
2541  os << std::endl;
2542 
2543  os << " Fast Reaction Pointer of the m-th reaction mode : whenever is not a small number, " << std::endl;
2544  os << " the r-th reaction is said to be a fast reaction " << std::endl;
2545  os << std::endl;
2546 
2547  os << " Participation Index : is a measure of participation of the r-th elementary reaction " << std::endl;
2548  os << " to the balancing act of the i-th mode " << std::endl;
2549  os << std::endl;
2550 
2551  os << " Importance Index: is a measure of relative importance of the contribution of r-th elementary " << std::endl;
2552  os << " reaction to the current reaction rate of i-th species " << std::endl;
2553  os << std::endl;
2554 
2555  os << " Species : " << std::endl;
2556 
2557  for (i = 0; i < mData.dim; i++)
2558  os << mpModel->getStateTemplate().beginIndependent()[i]->getObjectName() << std::endl;
2559 
2560  os << std::endl;
2561 
2562  os << " Reactions : " << std::endl;
2563 
2564  for (r = 0; r < (C_INT) reacs.size(); r++)
2565  os << reacs[r]->getObjectName() << std::endl;
2566 
2567  os << std::endl;
2568 
2569  os << "%%% Number of fast modes: " << std::endl;
2570 
2571  for (istep = 0; istep < stepNumber; istep++)
2572  {
2573 
2574  M = mVec_SlowModes[istep];
2575 
2576  os << std::endl;
2577  os << "%%% Time step " << istep + 1 << std::endl;
2578  os << std::endl;
2579 
2580  os << M << std::endl;
2581 
2582  os << std::endl;
2583  }
2584 
2585  os << std::endl;
2586  os << "%%% Time scales: " << std::endl;
2587 
2588  for (istep = 0; istep < stepNumber; istep++)
2589  {
2590 
2591  M = mVec_SlowModes[istep];
2592 
2593  os << std::endl;
2594  os << "%%% Time step " << istep + 1 << std::endl;
2595  os << std::endl;
2596 
2597  for (i = 0; i < mData.dim; i++)
2598  {
2599  os << mVec_TimeScale[istep][i] << " ";
2600  }
2601 
2602  os << std::endl;
2603  }
2604 
2605  os << std::endl;
2606  os << "% Radical Pointer " << std::endl;
2607 
2608  CMatrix<C_FLOAT64> RP;
2609  RP.resize(mData.dim, mData.dim);
2610 
2611  for (istep = 0; istep < stepNumber; istep++)
2612  {
2613 
2614  M = mVec_SlowModes[istep];
2615 
2616  os << std::endl;
2617  os << "%%% Time step " << istep + 1 << std::endl;
2618  os << std::endl;
2619 
2620  for (m = 0; m < M; m++)
2621  {
2622  for (i = 0; i < mData.dim; i++)
2623  RP(i, m) = mVec_mRadicalPointer[istep][i][m];
2624  }
2625 
2626  for (m = M; m < mData.dim; m++)
2627  {
2628  for (i = 0; i < mData.dim; i++)
2629  RP(i, m) = 0;
2630  }
2631 
2632  os << std::endl;
2633 
2634  for (i = 0; i < mData.dim; i++)
2635  {
2636  for (m = 0; m < mData.dim; m++)
2637  os << RP(i, m) << " ";
2638 
2639  os << std::endl;
2640  }
2641 
2642  os << std::endl;
2643  }
2644 
2645  os << std::endl;
2646  os << "%%%% Participation Index : " << std::endl;
2647 
2648  for (istep = 0; istep < stepNumber; istep++)
2649  {
2650 
2651  M = mVec_SlowModes[istep];
2652 
2653  os << std::endl;
2654  os << "%%% Time step " << istep + 1 << std::endl;
2655  os << std::endl;
2656 
2657  // os << istep + 1 << " ";
2658 
2659  for (r = 0; r < (C_INT) reacs.size(); r++)
2660  {
2661 
2662  for (i = 0; i < mData.dim; i++)
2663  os << mVec_mParticipationIndex[istep][r][i] << " ";
2664 
2665  os << std::endl;
2666  }
2667 
2668  os << std::endl;
2669  }
2670 
2671  os << std::endl;
2672  os << "%%% Importance Index " << std::endl;
2673 
2674  for (istep = 0; istep < stepNumber; istep++)
2675  {
2676 
2677  M = mVec_SlowModes[istep];
2678 
2679  os << std::endl;
2680  os << "%%% Time step " << istep + 1 << std::endl;
2681  os << std::endl;
2682 
2683  for (r = 0; r < (C_INT) reacs.size(); r++)
2684  {
2685 
2686  for (i = 0; i < mData.dim; i++)
2687  os << mVec_mImportanceIndex[istep][r][i] << " ";
2688 
2689  os << std::endl;
2690  }
2691  }
2692 
2693  return;
2694 }
void CSPParticipationIndex(C_INT &N, C_INT &M, C_FLOAT64 &tauM1, CMatrix< C_FLOAT64 > &B0)
void schur(C_INT &info)
std::vector< CMatrix< C_FLOAT64 > > mVec_mImportanceIndex
Definition: CCSPMethod.h:321
#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)
CMatrix< C_FLOAT64 > mAmplitudeTab
Definition: CCSPMethod.h:366
#define MCTSSAMethod
CCSPMethod(const CCopasiContainer *pParent=NULL)
Definition: CCSPMethod.cpp:36
CMatrix< C_FLOAT64 > mParticipationIndexNormedColumn
Definition: CCSPMethod.h:132
CArrayAnnotation * pParticipationIndexNormedRowAnn
Definition: CCSPMethod.h:338
#define pdelete(p)
Definition: copasi.h:215
CArrayAnnotation * pSlowParticipationIndexAnn
Definition: CCSPMethod.h:341
CCopasiProblem * getProblem()
std::vector< C_FLOAT64 > mCurrentTime
Definition: CTSSAMethod.h:460
CMatrix< C_FLOAT64 > mParticipationIndexTab
Definition: CCSPMethod.h:370
std::vector< CMatrix< C_FLOAT64 > > mVec_mFastParticipationIndex
Definition: CCSPMethod.h:319
const std::string & getObjectName() const
void findCandidatesNumber(C_INT &n, C_INT &k, CVector< C_FLOAT64 > &tsc, C_INT &info)
Definition: CCSPMethod.cpp:328
CArrayAnnotation * pTmp5NormedRow
Definition: CCSPMethod.h:360
C_FLOAT64 mTsc
Definition: CCSPMethod.h:61
virtual size_t size() const
void CSPOutput(C_INT &N, C_INT &M, C_INT &R)
virtual void step(const double &deltaT)
void findTimeScaleSeparation(C_INT &n, C_INT &k, CVector< C_FLOAT64 > &tsc, C_INT &info)
#define TAU
CState * mpState
Definition: CTSSAMethod.h:175
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
virtual size_t numRows() const
Definition: CMatrix.h:138
std::vector< CMatrix< C_FLOAT64 > > mVec_mRadicalPointer
Definition: CCSPMethod.h:313
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
void sminverse(C_INT &n, CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B)
Definition: CCSPMethod.cpp:164
void initializeIntegrationsParameter()
Definition: CState.h:305
void integrationStep(const double &deltaT)
std::vector< CMatrix< C_FLOAT64 > > mVec_mSlowParticipationIndex
Definition: CCSPMethod.h:320
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
bool modesAreExhausted(C_INT &N, C_INT &M, C_FLOAT64 &tauM, C_FLOAT64 &tauM1, CVector< C_FLOAT64 > &g, CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B, CMatrix< C_FLOAT64 > &F)
void calculateJacobian(CMatrix< C_FLOAT64 > &jacobian, const C_FLOAT64 &derivationFactor, const C_FLOAT64 &resolution)
Definition: CModel.cpp:2006
CMatrix< C_FLOAT64 > mFastParticipationIndexTab
Definition: CCSPMethod.h:376
void CSPImportanceIndex(C_INT &N, C_FLOAT64 &tauM1, CMatrix< C_FLOAT64 > &Q)
void setDescription(const std::string &s)
void CSPradicalPointer(C_INT &N, C_INT &M, CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B)
CMatrix< C_FLOAT64 > mImportanceIndexNormedRow
Definition: CCSPMethod.h:142
#define C_INT32
Definition: copasi.h:90
void smadd(CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B, CMatrix< C_FLOAT64 > &C, C_INT &n1, C_INT &n2)
Definition: CCSPMethod.cpp:114
CMatrix< C_FLOAT64 > mFastReactionPointerNormedTab
Definition: CCSPMethod.h:369
Definition: CMetab.h:178
CMatrix< C_FLOAT64 > mQ
Definition: CTSSAMethod.h:215
const CMatrix< C_FLOAT64 > & getRedStoi() const
Definition: CModel.cpp:1154
std::vector< CMatrix< C_FLOAT64 > > mVec_mParticipationIndexNormedRow
Definition: CCSPMethod.h:317
void yCorrection(C_INT &N, C_INT &M, CVector< C_FLOAT64 > &y, CMatrix< C_FLOAT64 > &TAUM, CMatrix< C_FLOAT64 > &F, CMatrix< C_FLOAT64 > &A)
CMatrix< C_FLOAT64 > mImportanceIndex
Definition: CCSPMethod.h:141
CVector< C_FLOAT64 > mFastParticipationIndex
Definition: CCSPMethod.h:133
void setAnnotationString(size_t d, size_t i, const std::string s)
C_INT mCSPbasis
Definition: CCSPMethod.h:97
CMatrix< C_FLOAT64 > mI
Definition: CCSPMethod.h:51
CArrayAnnotation * pTmp3
Definition: CCSPMethod.h:350
CTSSATask * pTask
CArrayAnnotation * pFastReactionPointerNormedAnn
Definition: CCSPMethod.h:336
const C_FLOAT64 & getInitialValue() const
void modesAmplitude(C_INT &N, C_INT &M, CVector< C_FLOAT64 > &g, CMatrix< C_FLOAT64 > &B, CMatrix< C_FLOAT64 > &F)
std::vector< CMatrix< C_FLOAT64 > > mVec_mAmplitude
Definition: CCSPMethod.h:312
CMatrix< C_FLOAT64 > mJacobian_initial
Definition: CTSSAMethod.h:210
CArrayAnnotation * pTmp4NormedColumn
Definition: CCSPMethod.h:353
CVector< C_FLOAT64 > mSlowParticipationIndex
Definition: CCSPMethod.h:134
CMatrix< C_FLOAT64 > mParticipationIndexNormedRowTab
Definition: CCSPMethod.h:371
void emptyOutputData(C_INT &N, C_INT &M, C_INT &R)
CArrayAnnotation * pParticipationIndexAnn
Definition: CCSPMethod.h:337
CMatrix< C_FLOAT64 > mRadicalPointer
Definition: CCSPMethod.h:115
const C_FLOAT64 & getNumber2QuantityFactor() const
Definition: CModel.cpp:2357
CArrayAnnotation * pTmp2
Definition: CCSPMethod.h:349
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
CArrayAnnotation * pTmp1
Definition: CCSPMethod.h:348
CVector< C_FLOAT64 > mG
Definition: CCSPMethod.h:81
std::vector< std::string > tableNames
Definition: CTSSAMethod.h:98
int mCurrentStep
Definition: CTSSAMethod.h:466
void createAnnotationsM()
const Value & getValue() const
bool mReducedModel
Definition: CTSSAMethod.h:291
CMatrix< C_FLOAT64 > mParticipationIndexNormedRow
Definition: CCSPMethod.h:131
static CCopasiVector< CCopasiDataModel > * getDatamodelList()
unsigned C_INT32 * pUINT
std::vector< CMatrix< C_FLOAT64 > > mVec_mParticipationIndex
Definition: CCSPMethod.h:316
CMatrix< C_FLOAT64 > mFastReactionPointerNormed
Definition: CCSPMethod.h:123
CArrayAnnotation * pImportanceIndexNormedRowAnn
Definition: CCSPMethod.h:343
void setCopasiVector(size_t d, const CCopasiContainer *v)
virtual void print(std::ostream *ostream) const
C_FLOAT64 mEps
Definition: CCSPMethod.h:56
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
CArrayAnnotation * pTmp4Fast
Definition: CCSPMethod.h:356
CVector< C_FLOAT64 > mYerror
Definition: CCSPMethod.h:86
#define PI
void smmult(CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B, CMatrix< C_FLOAT64 > &C, C_INT &n1, C_INT &n2, C_INT &n3)
Definition: CCSPMethod.cpp:85
CArrayAnnotation * pTmp4NormedRow
Definition: CCSPMethod.h:354
CMatrix< C_FLOAT64 > mSlowParticipationIndexTab
Definition: CCSPMethod.h:377
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
C_FLOAT64 mRerror
Definition: CCSPMethod.h:66
CArrayAnnotation * pAmplitudeAnn
Definition: CCSPMethod.h:333
virtual void resize(const size_t &newSize)
void perturbateA(C_INT &n, CMatrix< C_FLOAT64 > &A, C_FLOAT64 delta)
Definition: CCSPMethod.cpp:153
size_t size() const
Definition: CVector.h:100
CArrayAnnotation * pTmp4
Definition: CCSPMethod.h:352
void initializeParameter()
Definition: CCSPMethod.cpp:70
CArrayAnnotation * pImportanceIndexAnn
Definition: CCSPMethod.h:342
C_INT mTStep
Definition: CCSPMethod.h:93
void setMode(size_t d, Mode m)
void calculateDerivativesX(C_FLOAT64 *derivativesX)
Definition: CModel.cpp:1941
CMatrix< C_FLOAT64 > mParticipationIndex
Definition: CCSPMethod.h:130
void setDimensionDescription(size_t d, const std::string &s)
CArrayAnnotation * pTmp3Normed
Definition: CCSPMethod.h:351
CMatrix< C_FLOAT64 > mFastReactionPointerTab
Definition: CCSPMethod.h:368
CMatrix< C_FLOAT64 > mR
Definition: CTSSAMethod.h:220
CArrayAnnotation * pTmp5
Definition: CCSPMethod.h:359
void emptyVectors()
std::vector< CVector< C_FLOAT64 > > mVec_TimeScale
Definition: CTSSAMethod.h:461
CMatrix< C_FLOAT64 > mJacobian
Definition: CTSSAMethod.h:205
Header file of class CArrayAnnotation.
CMatrix< C_FLOAT64 > mImportanceIndexTab
Definition: CCSPMethod.h:373
#define C_FLOAT64
Definition: copasi.h:92
CArrayAnnotation * pTmp4Slow
Definition: CCSPMethod.h:357
std::vector< CMatrix< C_FLOAT64 > > mVec_mImportanceIndexNormedRow
Definition: CCSPMethod.h:322
CArrayAnnotation * pParticipationIndexNormedColumnAnn
Definition: CCSPMethod.h:339
CType * array()
Definition: CVector.h:139
CMatrix< C_FLOAT64 > mB
Definition: CCSPMethod.h:91
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
CArrayAnnotation * pFastReactionPointerAnn
Definition: CCSPMethod.h:335
C_INT isBlockDiagonal(C_INT &N, C_INT &M, CMatrix< C_FLOAT64 > &ALA, C_FLOAT64 SMALL)
Definition: CCSPMethod.cpp:938
void smsubst(CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B, CMatrix< C_FLOAT64 > &C, C_INT &n1, C_INT &n2)
Definition: CCSPMethod.cpp:102
CCopasiVectorNS< CReaction > & getReactions()
Definition: CModel.cpp:1039
virtual void predifineAnnotation()
std::vector< CMatrix< C_FLOAT64 > > mVec_mFastReactionPointerNormed
Definition: CCSPMethod.h:315
void calculateDerivatives(C_FLOAT64 *derivatives)
Definition: CModel.cpp:1903
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
virtual bool setAnnotationM(size_t step)
CArrayAnnotation * pRadicalPointerAnn
Definition: CCSPMethod.h:334
C_INT mSetVectors
Definition: CCSPMethod.h:103
CMatrix< C_FLOAT64 > mImportanceIndexNormedRowTab
Definition: CCSPMethod.h:374
const CMatrix< C_FLOAT64 > & getStoi() const
Definition: CModel.cpp:1160
C_FLOAT64 mAerror
Definition: CCSPMethod.h:71
void cspstep(const double &deltaT, C_INT &n, C_INT &m, CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B)
Definition: CCSPMethod.cpp:401
void setVectors(int fast)
virtual size_t numCols() const
Definition: CMatrix.h:144
CMatrix< C_FLOAT64 > mRadicalPointerTab
Definition: CCSPMethod.h:367
C_INT mIter
Definition: CCSPMethod.h:76
std::vector< CMatrix< C_FLOAT64 > > mVec_mParticipationIndexNormedColumn
Definition: CCSPMethod.h:318
CArrayAnnotation * pFastParticipationIndexAnn
Definition: CCSPMethod.h:340
C_FLOAT64 mTime
Definition: CTSSAMethod.h:200
CMatrix< C_FLOAT64 > mFastReactionPointer
Definition: CCSPMethod.h:122
const CCompartment * getCompartment() const
Definition: CMetab.cpp:222
void integrationMethodStart(const CState *initialState)
virtual CType * array()
Definition: CMatrix.h:337
std::vector< C_INT > mVec_SlowModes
Definition: CTSSAMethod.h:459
CTSSAMethod * pMethod
Definition: CTSSAMethod.h:169
void setVectorsToNaN()
void smnorm(C_INT &n, CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B, C_INT &n1)
Definition: CCSPMethod.cpp:126
void printResult(std::ostream *ostream) const
C_FLOAT64 * mY
Definition: CTSSAMethod.h:185
void basisRefinement(C_INT &N, C_INT &M, CMatrix< C_FLOAT64 > &ALA, CMatrix< C_FLOAT64 > &TAU, CMatrix< C_FLOAT64 > &A, CMatrix< C_FLOAT64 > &B, CMatrix< C_FLOAT64 > &A0, CMatrix< C_FLOAT64 > &B0)
std::vector< CMatrix< C_FLOAT64 > > mVec_mFastReactionPointer
Definition: CCSPMethod.h:314
virtual void start(const CState *initialState)
CVector< C_FLOAT64 > mAmplitude
Definition: CCSPMethod.h:109
CMatrix< C_FLOAT64 > mParticipationIndexNormedColumnTab
Definition: CCSPMethod.h:372
#define max(a, b)
Definition: f2c.h:176