104 mReducedModel(false),
129 mpCurrentState(src.mpCurrentState),
130 mpProblem(src.mpProblem),
156 mReducedModel(src.mReducedModel),
285 for (; it != end; ++it)
321 return std::numeric_limits<C_FLOAT64>::quiet_NaN();
337 if ((pParm =
getParameter(
"LSODA.RelativeTolerance")) != NULL)
342 if ((pParm =
getParameter(
"LSODA.AbsoluteTolerance")) != NULL)
348 if ((pParm =
getParameter(
"LSODA.AdamsMaxOrder")) != NULL)
353 if ((pParm =
getParameter(
"LSODA.BDFMaxOrder")) != NULL)
358 if ((pParm =
getParameter(
"LSODA.MaxStepsInternal")) != NULL)
366 if ((pParm =
getParameter(
"Use Default Absolute Tolerance")) != NULL)
380 assert(pDataModel != NULL);
392 for (i = 0, imax = Compartment.
size(); i < imax; i++)
393 if (Compartment[i]->
getValue() < Volume)
394 Volume = Compartment[i]->
getValue();
405 setValue(
"Absolute Tolerance", NewValue);
465 if ((mLsodaStatus <= 0))
470 if (mLsodaStatus == 3)
496 Matrix_anal.
resize(dim, dim);
502 for (j = 0; j < dim; j++)
505 for (i = 0; i < dim; i++)
507 for (j = 0; j < dim; j++)
511 for (i = 0; i < dim; i++)
513 for (j = 0; j < dim; j++)
520 for (i = 0; i < dim; i++)
521 for (j = 0; j < dim; j++)
541 for (j = 0; j < dim; j++)
544 for (i = 0; i < dim; i++)
545 for (j = 0; j < dim; j++)
546 denom[i] = denom[i] + fabs(
mTd(i, j));
548 for (i = 0; i < dim; i++)
549 for (j = 0; j < dim; j++)
553 for (i = 0; i < dim; i++)
554 for (j = 0; j < dim; j++)
574 Matrix_anal.
resize(dim, dim);
576 for (j = 0; j < dim; j++)
580 for (i = 0; i < dim; i++)
583 length = sqrt(length);
586 for (i = 0; i < dim; i++)
587 Matrix_anal(i, j) =
mTdInverse(i, j) / length;
590 if ((slow < dim) & (slow > 0))
594 for (i = 0; i < dim; i++)
596 for (j = 0; j < slow; j++)
597 denom = denom + fabs(Matrix_anal(j, i));
600 for (i = 0; i < dim; i++)
603 for (j = 0; j < dim; j++)
605 for (i = 0; i < slow; i++)
612 for (i = 0; i < dim; i++)
632 Matrix_anal.
resize(dim, dim);
634 for (j = 0; j < dim; j++)
638 for (i = 0; i < dim; i++)
641 length = sqrt(length);
644 for (i = 0; i < dim; i++)
645 Matrix_anal(i, j) =
mTdInverse(i, j) / length;
652 for (i = 0; i < dim; i++)
654 for (j = slow; j < dim; j++)
655 denom = denom + fabs(Matrix_anal(j, i));
658 for (i = 0; i < dim; i++)
661 for (j = 0; j < dim; j++)
663 for (i = slow; i < dim; i++)
670 for (i = 0; i < dim; i++)
682 C_FLOAT64 scalar_product, absolute_value_1;
698 for (i = 0; i < dim; ++i)
704 d_transformed.
resize(dim);
706 for (j = 0; j < dim; j++)
708 for (i = 0; i < slow; i++)
711 for (i = slow; i < dim; i++)
712 d_full[i] =
mQ(j, i) * Xconc[j];
714 for (i = 0; i < dim; i ++)
716 d_transformed[i] = 0.0;
718 for (k = 0; k < dim; k++)
719 d_transformed[i] +=
mQ(i, k) * d_full[k];
722 scalar_product = d_transformed[j];
723 absolute_value_1 = 0.0;
725 for (i = 0; i < dim; i++)
726 absolute_value_1 += d_transformed[i] * d_transformed[i];
728 absolute_value_1 = sqrt(absolute_value_1);
730 if (absolute_value_1 * Xconc[j] > 0.0)
731 scalar_product = scalar_product / absolute_value_1;
733 scalar_product = -2.0;
748 for (k = 0; k < dim; k++)
901 for (i = 0; i < dim; i++)
902 for (j = 0; j < dim; j++)
913 C_INT lwork = 10 * dim;
920 dgees_(&V, &N, NULL, &dim, R.
array(), &dim, &SDIM, eval_r.
array(), eval_i.
array(), Q.
array(), &dim, work.
array(), & lwork, Bwork.
array(), &info);
935 for (i = 0; i < dim; i++)
936 eval_reor[i] = eval_r [i];
938 for (i = 0; i < dim; i++) index[i] = 0;
954 for (i = 0; i < dim; i++)
964 while (changed ==
true)
969 while (count < dim - 1)
974 bool diagorder =
false;
983 diagorder = (index[count + 1] < index[count]);
987 diagorder = (index[count + 1] < index[count]);
991 diagorder = (index[count + 1] > index[count]);
1010 for (i = 0; i < dim; i++)
1011 j_diag[i] = R.
array()[i * dim + i];
1095 dtrexc_(&V, &dim, R.
array(), &dim, Q.
array(), &dim, &first, &second, work1.array(), &info);
1105 if (nid[count] == 0)
1107 if (pid[count] == 0)
1112 index[count] = index[count + 1];
1113 index[count + 1] = tmp;
1124 tmp = index[count - 1];
1125 index[count - 1] = index[count + 1];
1126 index[count + 1] = tmp;
1136 if (pid[count] == 0)
1141 index[count] = index[count + 2];
1142 index[count + 2] = tmp;
1154 index[count - 1] = index[count + 1];
1155 index[count] = index[count + 1];
1156 index[count + 1] = tmp;
1157 index[count + 2] = tmp;
1171 for (i = 0; i < dim; i++)
1172 for (j = 0; j < dim; j++)
1174 mQ(j, i) = Q.
array()[j + dim * i];
1175 mR(j, i) = R.
array()[j + dim * i];
1328 for (i = 0; i < dim; i++)
1329 for (j = 0; j < dim; j++)
1340 C_INT lwork = 10 * dim;
1347 dgees_(&V, &N, NULL, &dim, R.
array(), &dim, &SDIM, eval_r.
array(), eval_i.
array(), Q.
array(), &dim, work.
array(), & lwork, Bwork.
array(), &info);
1362 for (i = 0; i < dim; i++)
1363 eval_reor[i] = eval_r [i];
1365 for (i = 0; i < dim; i++) index[i] = 0;
1381 for (i = 0; i < dim; i++)
1388 bool changed =
true;
1391 while (changed ==
true)
1396 while (count < dim - 1)
1401 if (index[count + 1] < index[count])
1410 for (i = 0; i < dim; i++)
1411 j_diag[i] = R.
array()[i * dim + i];
1495 dtrexc_(&V, &dim, R.
array(), &dim, Q.
array(), &dim, &first, &second, work1.array(), &info);
1505 if (nid[count] == 0)
1507 if (pid[count] == 0)
1512 index[count] = index[count + 1];
1513 index[count + 1] = tmp;
1524 tmp = index[count - 1];
1525 index[count - 1] = index[count + 1];
1526 index[count + 1] = tmp;
1536 if (pid[count] == 0)
1541 index[count] = index[count + 2];
1542 index[count + 2] = tmp;
1554 index[count - 1] = index[count + 1];
1555 index[count] = index[count + 1];
1556 index[count + 1] = tmp;
1557 index[count + 2] = tmp;
1571 for (i = 0; i < dim; i++)
1572 for (j = 0; j < dim; j++)
1595 C_INT fast = dim - slow;
1603 st_slow.
resize(slow * slow);
1606 st_fast.
resize(fast * fast);
1609 st_coup.
resize(slow * fast);
1612 S_coup.
resize(slow, fast);
1614 for (i = 0; i < slow; i++)
1615 for (j = 0; j < slow; j++)
1616 st_slow[j + slow * i] =
mR(j, i);
1618 for (i = 0; i < fast; i++)
1619 for (j = 0; j < fast; j++)
1620 st_fast[j + fast * i] =
mR(j + slow, i + slow);
1622 for (i = 0; i < slow; i++)
1623 for (j = 0; j < fast; j++)
1624 S_coup(i, j) =
mR(i, j + slow);
1626 for (j = 0; j < fast; j++)
1627 for (i = 0; i < slow; i++)
1628 st_coup[i + slow * j] = -S_coup(i, j);
1716 dtrsyl_(&N1, &N2, &isgn, &slow, &fast, st_slow.
array(), &slow, st_fast.
array(), &fast, st_coup.
array(), &slow, &scale, &info);
1730 for (i = 0; i < dim; i++)
1731 for (j = 0; j < dim; j++)
1734 for (j = 0; j < fast; j++)
1735 for (i = 0; i < slow; i++)
1736 Cmat(i, j + slow) = st_coup[i + slow * j];
1738 for (i = 0; i < dim; i++)
1739 for (j = 0; j < dim; j++)
1741 C(i, j) = 1. + Cmat(i, j);
1743 C(i, j) = Cmat(i, j);
1745 for (i = 0; i < dim; i++)
1747 for (j = 0; j < dim; j++)
1751 for (k = 0; k < dim; k++)
1752 mTd(i, j) =
mTd(i, j) +
mQ(i, k) * C(k, j);
1756 for (i = 0; i < dim; i++)
1757 for (j = 0; j < dim; j++)
1759 C(i, j) = 1. - Cmat(i, j);
1761 C(i, j) = - Cmat(i, j);
1764 transp_Qt.
resize(dim, dim);
1766 for (i = 0; i < dim; i++)
1767 for (j = 0; j < dim; j++)
1768 transp_Qt(i, j) =
mQ(j, i);
1770 for (i = 0; i < dim; i++)
1772 for (j = 0; j < dim; j++)
1776 for (k = 0; k < dim; k++)
1787 for (i = 0; i < dim; i++)
1789 for (j = 0; j < dim; j++)
1793 for (k = 0; k < dim; k++)
1798 for (i = 0; i < dim; i++)
1800 for (j = 0; j < dim; j++)
1804 for (k = 0; k < dim; k++)
1805 S(i, j) = S(i, j) +
mTdInverse(i, k) * E(k, j);
1809 C_INT flag_sylvester;
1813 for (i = 0; i < dim; i++)
1814 for (j = 0; j < dim; j++)
1817 for (i = 0; i < fast; i++)
1818 for (j = 0; j < fast; j++)
1819 mQz(i, j) = S(i + slow, j + slow);
1836 for (i = 0, imax = indep; i < imax; i++)
1843 for (i = 0, imax = indep; i < imax; i++)
1854 for (i = 0; i < imax; ++i)
1855 Y1[i] *= number2conc;
1858 for (i = 0, imax = indep; i < imax; i++)
1878 for (i = 0, imax = nmetab; i < imax; i++)
1885 for (i = 0, imax = nmetab; i < imax; i++)
1894 for (i = 0; i < imax; ++i)
1895 Y1[i] *= number2conc;
1898 for (i = 0, imax = nmetab; i < imax; i++)
1918 max_value = eval_r[0];
1920 for (i = 1; i < dim; i++)
1921 if (eval_r[i] > max_value)
1922 max_value = eval_r[i];
1937 for (i = 0; i < dim; i++)
1941 abs_eval_r[i] = (eval_r[i]);
1946 for (i = 0; i < dim; i++)
1950 for (j = 0; j < dim; j++)
1953 if (abs_eval_r[j] < abs_eval_r[max])
1958 abs_eval_r[
max] = factor * max_value;
1962 for (i = 0; i < dim - 1; i++)
1963 if (eval_r[i] == eval_r[i + 1])
1964 index[i + 1] = index[i];
1983 min_value = eval_r[0];
1985 for (i = 1; i < dim; i++)
1986 if (eval_r[i] < min_value)
1987 min_value = eval_r[i];
1989 for (i = 0; i < dim; i++)
1993 abs_eval_r[i] = (eval_r[i]);
1998 for (i = 0; i < dim; i++)
2002 for (j = 0; j < dim; j++)
2005 if (abs_eval_r[j] >= abs_eval_r[min])
2010 abs_eval_r[
min] = factor * min_value;
2014 for (i = 0; i < dim - 1; i++)
2015 if (eval_r[i] == eval_r[i + 1])
2016 index[i + 1] = index[i];
2028 for (k = 0; k < dim; k++)
2031 for (k = 1; k < dim - 1; k++)
2032 if (index[k] == index[k + 1])
2045 for (k = 0; k < dim; k++)
2048 for (k = 1; k < dim; k++)
2049 if (index[k] == index[k - 1])
2115 for (; pAtol != pEnd; ++pAtol, ++ppEntity)
2120 if ((pMetab = dynamic_cast< const CMetab * >(*ppEntity)) != NULL)
CCopasiDataModel * getObjectDataModel()
void update_pid(C_INT *index, C_INT *pid, const C_INT &dim)
void mat_anal_fast_space_thomas(C_INT &slow)
void sylvester(C_INT slow, C_INT &info)
void update_nid(C_INT *index, C_INT *nid, const C_INT &dim)
CCopasiVectorN< CEvent > & getEvents()
void setCurrentState(CState *currentState)
std::vector< C_FLOAT64 > mCurrentTime
int dtrsyl_(char *trana, char *tranb, integer *isgn, integer *m, integer *n, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *c__, integer *ldc, doublereal *scale, integer *info)
const CCopasiVector< CMetab > & getMetabolites() const
const CCopasiVectorN< CModelValue > & getModelValues() const
CVector< C_FLOAT64 > mYdot
CMatrix< C_FLOAT64 > mVslow
const int & getCurrentStep() const
virtual size_t size() const
CMatrix< C_FLOAT64 > mTdInverse
virtual void initializeParameter()
CVector< C_FLOAT64 > mY_initial
void setOstream(std::ostream &os)
void updateSimulatedValues(const bool &updateMoieties)
void createAnnotationsM()
CVector< C_FLOAT64 > mVslow_space
CVector< C_FLOAT64 > mDWork
virtual void predifineAnnotation()
void initializeIntegrationsParameter()
static CTSSAMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::unset)
void integrationStep(const double &deltaT)
CVector< C_FLOAT64 > mAtol
void resize(size_t size, const bool ©=false)
void schur_desc(C_INT &info)
void setVectors(int slowMode)
std::ostringstream mErrorMsg
size_t getNumODEMetabs() const
const CCopasiMethod::SubType & getSubType() const
void map_index(C_FLOAT64 *eval_r, C_INT *index, const C_INT &dim)
double orthog(C_INT &number1, C_INT &number2)
void mat_anal_mod_space(C_INT &slow)
void setTime(const C_FLOAT64 &time)
CVector< C_FLOAT64 > mVfast_space
void map_index_desc(C_FLOAT64 *eval_r, C_INT *index, const C_INT &dim)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
CMatrix< C_FLOAT64 > mQ_desc
CMatrix< C_FLOAT64 > mJacobian_initial
bool removeParameter(const std::string &name)
#define MCTrajectoryMethod
std::vector< CType * >::const_iterator const_iterator
const C_FLOAT64 & getQuantity2NumberFactor() const
const C_FLOAT64 & getNumber2QuantityFactor() const
size_t getNumDependentReactionMetabs() const
void calculateDerivativesX(C_FLOAT64 *X1, C_FLOAT64 *Y1)
CMatrix< C_FLOAT64 > mR_desc
void setState(const CState &state)
void calculateDerivatives(C_FLOAT64 *X1, C_FLOAT64 *Y1)
void mat_anal_fast_space(C_INT &slow)
int dtrexc_(char *compq, integer *n, doublereal *t, integer *ldt, doublereal *q, integer *ldq, integer *ifst, integer *ilst, doublereal *work, integer *info)
const Value & getValue() const
bool setValue(const std::string &name, const CType &value)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
size_t getNumIndependent() const
void setModel(CModel *model)
CModelEntity ** beginIndependent()
CCopasiParameter * getParameter(const std::string &name)
CVector< C_FLOAT64 > getVec_TimeScale(int step)
virtual void resize(size_t rows, size_t cols, const bool ©=false)
size_t getNumIndependentReactionMetabs() const
void mat_anal_metab(C_INT &slow)
void calculateDerivativesX(C_FLOAT64 *derivativesX)
CCopasiVectorNS< CCompartment > & getCompartments()
const C_FLOAT64 & getValue() const
const C_FLOAT64 & getTime() const
std::vector< CVector< C_FLOAT64 > > mVec_TimeScale
CMatrix< C_FLOAT64 > mJacobian
C_FLOAT64 returnCurrentTime(int step)
void setProblem(CTSSAProblem *problem)
size_t getNumAssignmentMetabs() const
const CStateTemplate & getStateTemplate() const
virtual void start(const CState *initialState)
int dgees_(char *jobvs, char *sort, L_fp select, integer *n, doublereal *a, integer *lda, integer *sdim, doublereal *wr, doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work, integer *lwork, logical *bwork, integer *info)
void calculateDerivatives(C_FLOAT64 *derivatives)
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
const CCopasiVector< CMetab > & getMetabolitesX() const
virtual void step(const double &deltaT)
void integrationMethodStart(const CState *initialState)
CModel * getModel() const
void mat_anal_mod(C_INT &slow)
CMatrix< C_FLOAT64 > mVslow_metab
C_FLOAT64 * beginIndependent()
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)