85 C_INT failed_while = 0;
89 C_INT slow = dim - fast;
91 C_INT i, j, k, info_schur = 0;
117 Stoichiom =
mpModel -> getRedStoi();
121 reacs_size = reacs_size / dim;
134 for (i = 0; i < dim; ++i)
135 Xconc[i] =
mY[i] * number2conc;
137 for (i = 0; i < dim; i++)
141 Xconc_initial.
resize(dim);
143 for (i = 0; i < dim; ++i)
144 Xconc_initial[i] =
mY_initial[i] * number2conc;
147 for (i = 0; i < dim; i++)
148 for (j = 0; j < dim; j++)
164 std::cout <<
"Jacobian_next:" << std::endl;
172 Jacobian_for_test.
resize(dim, dim);
174 for (i = 0; i < dim; i++)
175 for (j = 0; j < dim; j++)
179 for (i = 0; i < dim; i++)
180 for (j = 0; j < dim; j++)
186 std::cout <<
"Eigenvalues of next Jacobian" << std::endl;
188 for (i = 0; i < dim; i++)
189 std::cout <<
mR(i, i) << std::endl;
193 for (i = 0; i < dim; i++)
194 for (j = 0; j < dim; j++)
199 for (i = 0; i < dim; i++)
200 for (j = 0; j < dim; j++)
206 for (i = 0; i < dim; i++)
207 for (j = 0; j < dim; j++)
234 std::cout <<
"Eigenvalues of initial Jacobian" << std::endl;
236 for (i = 0; i < dim; i++)
237 std::cout <<
mR(i, i) << std::endl;
241 std::cout <<
"Schur Decomposition" << std::endl;
242 std::cout <<
"mR - block upper triangular matrix :" << std::endl;
243 std::cout <<
mR << std::endl;
251 if (
mR(dim - 1, dim - 1) ==
mR(dim - 2 , dim - 2))
260 if (
mR(dim - 1, dim - 1) >= 0)
283 if (
mR(slow, slow) ==
mR(slow - 1 , slow - 1))
288 for (i = 0; i < dim; i++)
289 for (j = 0; j < dim; j++)
300 for (i = 0; i < dim; i++)
301 for (j = 0; j < dim; j++)
304 mTd(i, j) =
mQ(i, j);
305 mQz(i, j) =
mR(i, j);
326 for (i = slow ; i < dim; i++)
334 mEPS = 1 / fabs(
mR(slow , slow));
363 if ((failed == 1) || (failed_while == 1))
370 if ((fast >= 1) && (
mR(slow - 1, slow - 1) ==
mR(slow , slow)))
375 for (i = 0; i < dim; i++)
376 for (j = 0; j < dim; j++)
392 C_INT flag_orthog = 1;
394 if (flag_orthog == 0)
397 orthog_prove.
resize(dim, dim);
399 for (i = 0; i < dim; i++)
400 for (j = 0; j < dim; j++)
401 orthog_prove(i, j) =
orthog(i, j);
407 if (flag_develop == 0)
409 C_INT info_schur_desc = 0;
417 for (i = 0; i < dim; i++)
418 for (j = 0; j < dim; j++)
424 for (i = 0; i < dim; i++)
425 for (j = 0; j < dim; j++)
438 Reac_slow_space_orth.
resize(reacs_size);
440 for (i = 0; i < reacs_size; i ++)
442 Reac_slow_space_orth[i] = 0;
444 for (j = 0; j < dim; j++)
446 Reac_slow_space_orth[i] = Reac_slow_space_orth[i] +
mVfast_space[j] * Stoichiom(j, i);
450 for (i = 0; i < dim; i++)
451 for (j = 0; j < dim; j++)
471 if (flag_develop == 0)
474 Slow_react_contr.
resize(dim, reacs_size);
477 Mode_react_contr.
resize(dim, reacs_size);
480 Reac_slow_space.
resize(reacs_size);
483 Reac_fast_space.
resize(reacs_size);
487 for (i = 0; i < dim; i ++)
488 for (j = 0; j < reacs_size; j++)
490 Slow_react_contr(i, j) = 0;
492 for (k = 0; k < dim; k++)
493 Slow_react_contr(i, j) = Slow_react_contr(i, j) +
mVslow(i, k) * Stoichiom(k, j);
499 for (j = 0; j < dim; j++)
502 for (i = 0; i < dim; i++)
503 for (j = 0; j < reacs_size; j++)
504 denom_mode[i] = denom_mode[i] + fabs(Slow_react_contr(i, j));
506 for (i = 0; i < dim; i++)
507 for (j = 0; j < reacs_size; j++)
509 if (denom_mode[i] == 0)
510 Mode_react_contr(i, j) = 0;
512 Mode_react_contr(i, j) = (Slow_react_contr(i, j)) / denom_mode[i] * 100;
516 denom_reac.
resize(reacs_size);
518 for (j = 0; j < reacs_size; j++)
521 for (i = 0; i < reacs_size; i++)
522 for (j = 0; j < dim; j++)
523 denom_reac[i] = denom_reac[i] + fabs(Slow_react_contr(j, i));
525 for (i = 0; i < reacs_size; i++)
526 for (j = 0; j < dim; j++)
528 if (denom_reac[i] == 0)
529 Slow_react_contr(j, i) = 0;
531 Slow_react_contr(j, i) = (Slow_react_contr(j , i)) / denom_reac[i] * 100;
534 for (i = 0; i < reacs_size; i ++)
536 Reac_slow_space[i] = 0;
538 for (j = 0; j < dim; j++)
540 Reac_slow_space[i] = Reac_slow_space[i] +
mVslow_space[j] * Stoichiom(j, i);
547 for (i = 0; i < reacs_size; i++)
548 length = length + fabs(Reac_slow_space[i]);
550 for (i = 0; i < reacs_size; i++)
556 for (i = 0; i < reacs_size; i ++)
558 Reac_fast_space[i] = 0;
560 for (j = 0; j < dim; j++)
561 Reac_fast_space[i] = Reac_fast_space[i] +
mVfast_space[j] * Stoichiom(j, i);
566 for (i = 0; i < reacs_size; i++)
567 length = length + fabs(Reac_fast_space[i]);
569 for (i = 0; i < reacs_size; i++)
571 Reac_fast_space[i] = Reac_fast_space[i] / length * 100;
573 Reac_fast_space[i] = 0;
577 std::cout <<
"**********************************************************" << std::endl;
578 std::cout <<
"**********************************************************" << std::endl;
579 std::cout << std::endl;
580 std::cout << std::endl;
587 for (i = 0; i < dim; i++)
588 for (j = 0; j < reacs_size; j++)
590 mTMP1(j, i) = Mode_react_contr(i, j);
591 mTMP2(j, i) = Slow_react_contr(i, j);
594 for (j = 0; j < reacs_size; j++)
595 mTMP3(j, 0) = Reac_fast_space[j];
622 C_INT i, j, iter, iterations, itermax;
623 C_INT nrhs, ok, fast;
634 s_22_array.
resize(fast * fast);
674 for (i = 0; i < fast; i++)
675 for (j = 0; j < fast; j++)
676 S_22(i, j) =
mQz(i, j);
678 for (i = 0; i < fast; i++)
681 for (i = 0; i < fast; i++)
682 for (j = 0; j < fast; j++)
683 s_22_array[j + fast * i] = S_22(j, i);
685 for (i = 0; i < fast; i++)
699 for (i = 0; i < fast; i++)
700 yf_newton[i] = yf_newton[i] + d_yf[i];
704 for (i = 0; i < slow; i++)
707 for (i = slow; i < dim; i++)
708 y_newton[i] = yf_newton[i - slow];
710 for (i = 0; i < dim; i++)
714 for (j = 0; j < dim; j++)
715 x_newton[i] = x_newton[i] +
mTd(i, j) * y_newton[j];
720 for (i = 0; i < dim; i++)
724 for (j = 0; j < dim; j++)
725 g_newton[i] = g_newton[i] +
mTdInverse(i, j) * dxdt_newton[j];
731 for (i = 0; i < fast; i++)
733 gf_newton[i] = -1. * g_newton[i + slow];
807 for (i = 0; i < fast; i++)
808 d_yf[i] = gf_newton[i];
822 for (i = 0; i < fast; i++)
824 gf_newton[i] = fabs(gf_newton[i]);
826 if (err < gf_newton[i])
830 iterations = iterations + 1;
853 for (i = 0; i < fast; i++)
891 C_INT fast = dim - slow;
915 for (i = 0; i < dim; ++i)
918 for (i = 0; i < dim; i++)
922 for (j = 0; j < dim; j++)
923 c_full[i] = c_full[i] +
mTdInverse(i, j) * Xconc[j];
926 for (j = 0; j < slow; j++)
927 c_slow[j] = c_full[j];
929 for (j = 0; j < fast; j++)
930 mCfast[j] = c_full[j + slow];
946 for (j = 0; j < dim; j++)
952 for (j = 0; j < dim; j++)
959 for (i = 0; i < dim; i++)
963 for (j = 0; j < dim; j++)
964 g_full[i] = g_full[i] +
mTdInverse(i, j) * dxdt[j];
967 for (j = 0; j < slow; j++)
968 g_slow[j] = g_full[j];
997 for (i = 0; i < slow; i++)
998 c_relax[i] = c_slow[i];
1000 for (i = slow; i < dim; i++)
1001 c_relax[i] =
mCfast[i - slow];
1003 for (i = 0; i < dim; i++)
1007 for (j = 0; j < dim; j++)
1008 x_relax[i] = x_relax[i] +
mTd(i, j) * c_relax[j];
1013 for (i = 0; i < dim; i++)
1017 for (j = 0; j < dim; j++)
1018 g_relax[i] = g_relax[i] +
mTdInverse(i, j) * dxdt_relax[j];
1026 for (i = 0; i < slow; i++)
1028 re[i] = fabs(g_relax[i] - g_slow[i]);
1029 re[i] = re[i] *
mEPS;
1034 for (i = 0; i < slow; i++)
1041 for (i = 0; i < slow; i++)
1042 norm = norm + fabs(g_relax[i] - g_slow[i]);
1087 for (i = 0; i < (size_t)
mData.
dim; i++)
1144 name =
"Contribution of species to modes";
1159 name =
"Modes distribution for species";
1174 name =
"Slow space";
1189 name =
"Fast space";
1206 name =
"Reactions slow space";
1223 name =
"Reactions contribution to the mode";
1238 name =
"Reactions distribution between modes ";
1253 name =
"Reactions fast space";
1278 if (step == 0)
return false;
1289 std::stringstream sstr;
1344 sstr << dim - mVec_SlowModes[
step];
1395 if (i < mVec_SlowModes[step])
1418 if (i < mVec_SlowModes[step])
1443 std::ostream & os = *ostream;
1445 C_INT i, j, istep = 0;
1451 assert(getObjectDataModel() != NULL);
1454 stepNumber = mVec_SlowModes.size();
1456 for (istep = 0; istep < stepNumber; istep++)
1460 os <<
"**************** Time step " << istep + 1 <<
" ************************** " << std::endl;
1464 os <<
"Contribution of species to modes" << std::endl;
1466 os <<
"Rows : contribution to mode (TS - corresponding timescale)" << std::endl;
1467 os <<
"Columns: species ";
1469 for (j = 0; j < mData.dim; j++)
1471 os << mpModel->getMetabolitesX()[j]->getObjectName() <<
" ";
1476 for (i = 0; i < mData.dim; i++)
1478 timeScale = mVec_TimeScale[istep][i];
1480 if (i < mVec_SlowModes[istep])
1485 os << timeScale <<
"): ";
1487 for (j = 0; j < mData.dim; j++)
1488 os << mVec_mVslow[istep][i][j] <<
" ";
1495 os <<
"Modes distribution for species" << std::endl;
1497 os <<
"Rows : Mode distribution for each species" << std::endl;
1498 os <<
"Columns: Modes (TS - corresponding timescale) ";
1501 for (i = 0; i < mData.dim; i++)
1503 timeScale = mVec_TimeScale[istep][i];
1505 if (i < mVec_SlowModes[istep])
1510 os << timeScale <<
") ";
1515 for (j = 0; j < mData.dim; j++)
1517 os << mpModel->getMetabolitesX()[j]->getObjectName() <<
" ";
1519 for (i = 0; i < mData.dim; i++)
1520 os << mVec_mVslowMetab[istep][j][i] <<
" ";
1527 os <<
"Slow space" << std::endl;
1529 os <<
"Rows : Species" << std::endl;
1530 os <<
"Column: Contribution to slow space ";
1533 os << mVec_SlowModes[istep];
1536 os << mData.dim - mVec_SlowModes[istep];
1540 for (j = 0; j < mData.dim; j++)
1542 os << mpModel->getMetabolitesX()[j]->getObjectName() <<
" ";
1543 os << mVec_mVslowSpace[istep][j] <<
" ";
1549 os <<
"Fast space" << std::endl;
1551 os <<
"Rows : Species" << std::endl;
1552 os <<
"Column: Contribution to fast space ";
1555 os << mVec_SlowModes[istep];
1558 os << mData.dim - mVec_SlowModes[istep];
1562 for (j = 0; j < mData.dim; j++)
1564 os << mpModel->getMetabolitesX()[j]->getObjectName() <<
" ";
1565 os << mVec_mVfastSpace[istep][j] <<
" ";
1571 os <<
"Reactions slow space" << std::endl;
1573 os <<
"Rows : Reactions" << std::endl;
1574 os <<
"Column: Contribution to slow space ";
1577 for (j = 0; j < (
C_INT32) mpModel->getReactions().size(); j++)
1579 os << mpModel->getReactions()[j]->getObjectName() <<
" ";
1580 os << mVec_mReacSlowSpace[istep][j] <<
" ";
void deuflhard(C_INT &slow, C_INT &info)
virtual void start(const CState *initialState)
int dgesv_(integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info)
void sylvester(C_INT slow, C_INT &info)
CMatrix< C_FLOAT64 > mVslowPrint
std::vector< C_FLOAT64 > mCurrentTime
CMatrix< C_FLOAT64 > mVslowSpacePrint
std::vector< CVector< C_FLOAT64 > > mVec_mVfastSpace
void newton(C_FLOAT64 *ys, C_INT &slow, C_INT &info)
CMatrix< C_FLOAT64 > mVslow
virtual size_t size() const
CMatrix< C_FLOAT64 > mTdInverse
CVector< C_FLOAT64 > mY_initial
CMatrix< C_FLOAT64 > mTMP3
std::vector< CMatrix< C_FLOAT64 > > mVec_mVslow
void updateSimulatedValues(const bool &updateMoieties)
CArrayAnnotation * pVslowSpacePrintAnn
void calculateJacobianX(CMatrix< C_FLOAT64 > &jacobianX, const C_FLOAT64 &derivationFactor, const C_FLOAT64 &resolution)
std::map< std::string, CArrayAnnotation * > mapTableToName
CVector< C_FLOAT64 > mVslow_space
CMatrix< C_FLOAT64 > mVfastSpacePrint
void initializeIntegrationsParameter()
CMatrix< C_FLOAT64 > mTd_save
void integrationStep(const double &deltaT)
CMatrix< C_FLOAT64 > mVslowMetabPrint
void resize(size_t size, const bool ©=false)
void schur_desc(C_INT &info)
std::vector< CMatrix< C_FLOAT64 > > mVec_mTMP2
void setDescription(const std::string &s)
virtual bool setAnnotationM(size_t step)
CArrayAnnotation * pVfastSpacePrintAnn
double orthog(C_INT &number1, C_INT &number2)
const CMatrix< C_FLOAT64 > & getRedStoi() const
CMatrix< C_FLOAT64 > mTdInverse_save
void mat_anal_mod_space(C_INT &slow)
virtual void step(const double &deltaT)
CArrayAnnotation * pTMP3PrintAnn
CVector< C_FLOAT64 > mVfast_space
void setAnnotationString(size_t d, size_t i, const std::string s)
std::vector< CVector< C_FLOAT64 > > mVec_mReacSlowSpace
CArrayAnnotation * pTMP2PrintAnn
CMatrix< C_FLOAT64 > mQ_desc
CMatrix< C_FLOAT64 > mJacobian_initial
const C_FLOAT64 & getNumber2QuantityFactor() const
void calculateDerivativesX(C_FLOAT64 *X1, C_FLOAT64 *Y1)
CMatrix< C_FLOAT64 > mR_desc
std::vector< CMatrix< C_FLOAT64 > > mVec_mTMP1
std::vector< std::string > tableNames
void mat_anal_fast_space(C_INT &slow)
CVector< C_FLOAT64 > mReacSlowSpace
const Value & getValue() const
std::vector< CVector< C_FLOAT64 > > mVec_mVslowSpace
CArrayAnnotation * pReacSlowSpacePrintAnn
void setCopasiVector(size_t d, const CCopasiContainer *v)
CILDMMethod(const CCopasiContainer *pParent=NULL)
CArrayAnnotation * pVslowPrintAnn
void createAnnotationsM()
virtual void resize(size_t rows, size_t cols, const bool ©=false)
void mat_anal_metab(C_INT &slow)
void printResult(std::ostream *ostream) const
virtual void initializeParameter()
void setMode(size_t d, Mode m)
CCopasiVectorNS< CCompartment > & getCompartments()
void setDimensionDescription(size_t d, const std::string &s)
CMatrix< C_FLOAT64 > mTMP1
std::vector< CVector< C_FLOAT64 > > mVec_TimeScale
CMatrix< C_FLOAT64 > mJacobian
CArrayAnnotation * pVslowMetabPrintAnn
CMatrix< C_FLOAT64 > mTMP2
std::vector< CMatrix< C_FLOAT64 > > mVec_mTMP3
virtual size_t size() const
CCopasiVectorNS< CReaction > & getReactions()
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
CMatrix< C_FLOAT64 > mTMP2Print
CMatrix< C_FLOAT64 > mReacSlowSpacePrint
const CCopasiVector< CMetab > & getMetabolitesX() const
CArrayAnnotation * pTMP1PrintAnn
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)
void integrationMethodStart(const CState *initialState)
std::vector< C_INT > mVec_SlowModes
void mat_anal_mod(C_INT &slow)
CMatrix< C_FLOAT64 > mVslow_metab
void setVectors(int slowMode)
CVector< C_FLOAT64 > mCfast
CMatrix< C_FLOAT64 > mTMP3Print
std::vector< CMatrix< C_FLOAT64 > > mVec_mVslowMetab
CMatrix< C_FLOAT64 > mTMP1Print