81 if ((pParm =
getParameter(
"Newton.UseIntegration")) != NULL)
87 if ((pParm =
getParameter(
"Newton.UseBackIntegration")) != NULL)
93 if ((pParm =
getParameter(
"Newton.acceptNegativeConcentrations")) != NULL)
99 if ((pParm =
getParameter(
"Newton.IterationLimit")) != NULL)
143 setValue(
"Use Back Integration",
false);
149 setValue(
"Use Back Integration",
false);
155 setValue(
"Use Back Integration",
false);
161 setValue(
"Use Back Integration",
true);
169 configBuffer.
getVariable(
"SSBackIntegration",
"bool", &Bool);
170 setValue(
"Use Back Integration", Bool);
172 configBuffer.
getVariable(
"NewtonLimit",
"C_INT32", &Int,
176 configBuffer.
getVariable(
"SSResoltion",
"C_FLOAT64", &Dbl);
188 C_FLOAT64 iterationFactor = forward ? 10.0 : 2.0;
198 MaxSteps = (
unsigned C_INT32) ceil(log(maxDuration / minDuration) / log(iterationFactor));
200 std::string tmpstring = forward ?
"forward integrating..." :
"backward integrating...";
216 assert(pTrajectoryProblem);
220 bool stepLimitReached =
false;
223 for (duration = minDuration; fabs(duration) <= fabs(maxDuration); duration *= iterationFactor, Step++)
241 mMethodLog <<
" Integration with duration " << duration <<
" failed (Exception).\n\n";
251 mMethodLog <<
" Integration with duration " << duration <<
" failed (NaN).\n\n";
259 mMethodLog <<
" Integration with duration " << duration
260 <<
" resulted in negative concentrations.\n\n";
273 mMethodLog <<
" Integration with duration " << duration
274 <<
". Criterium matched by " << value <<
".\n\n";
281 mMethodLog <<
" Integration with duration " << duration
282 <<
". Criterium not matched by " << value <<
".\n\n";
303 if (stepLimitReached)
306 mMethodLog <<
" Integration with duration " << duration
307 <<
" reached internal step limit.\n";
385 for (; pH != pHEnd; ++pH)
387 if (fabs(*pH) > 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon())
394 mMethodLog <<
" Newton step failed. Jacobian could not be inverted.\n\n";
400 C_FLOAT64 newValue = currentValue * 1.001;
406 for (i = 0; (i < 32) && !((newValue < currentValue)); i++)
413 for (; pHit != pEnd; ++pHit, ++pXit, ++pXoldIt)
415 *pXit = *pXoldIt - *pHit;
450 mMethodLog <<
" Newton step failed. Negative volume or concentration found.\n\n";
455 currentValue = newValue;
460 mMethodLog <<
" Regular Newton step. New value: " << currentValue <<
"\n";
462 mMethodLog <<
" Newton step with damping. New value: " << currentValue
463 <<
" (" << i - 1 <<
" damping iteration(s))\n";
517 mMethodLog <<
" Success: Target criterium matched by " << targetValue <<
".\n";
519 mMethodLog <<
" Failed: Target criterium not matched after reaching iteration limit. " << targetValue <<
"\n";
580 for (; pIt != pEnd; ++pIt)
618 const CMetab * pMetab = NULL;
626 for (; pDistance != pDistanceEnd; ++pDistance, ++pCurrentState, ++pAtol, ++ppEntity)
629 tmp = fabs(*pDistance) /
std::max(fabs(*pCurrentState), *pAtol);
630 RelativeDistance += tmp * tmp;
632 tmp = fabs(*pDistance);
634 if ((pMetab = dynamic_cast< const CMetab * >(*ppEntity)) != NULL)
639 AbsoluteDistance += tmp * tmp;
643 isnan(RelativeDistance) ? std::numeric_limits< C_FLOAT64 >::infinity() : sqrt(RelativeDistance);
645 isnan(AbsoluteDistance) ? std::numeric_limits< C_FLOAT64 >::infinity() : sqrt(AbsoluteDistance);
649 if (Error < TargetValue)
650 return TargetValue * (1.0 + Error);
668 if (!((*
getValue(
"Use Newton").pBOOL)
669 || (*
getValue(
"Use Integration").pBOOL)
670 || (*
getValue(
"Use Back Integration").pBOOL)))
677 if (*
getValue(
"Maximum duration for forward integration").pUDOUBLE <= 0)
683 if (*
getValue(
"Maximum duration for backward integration").pUDOUBLE <= 0)
708 if (*
getValue(
"Use Integration").pBOOL)
711 if (*
getValue(
"Use Back Integration").pBOOL)
714 if (*
getValue(
"Accept Negative Concentrations").pBOOL)
749 assert(pDataModel != NULL);
760 assert(pTrajectoryProblem);
767 assert(pTrajectoryMethod);
785 if (M == 0 || N == 0 || M != N)
787 return std::numeric_limits< C_FLOAT64 >::infinity();
790 C_INT LDA = std::max< C_INT >(1, M);
797 C_FLOAT64 * mpJTcolumnEnd = mpJTcolumn + M;
801 for (; mpJTcolumn != mpJTcolumnEnd; ++mpJTcolumn)
805 for (; mpJT < mpJTEnd; mpJT += M, ++mpJ)
814 C_FLOAT64 RCOND = 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon();
959 dgelsy_(&M, &N, &NRHS, JT.
array(), &LDA, X.
array(), &LDA, JPVT.
array(), &RCOND, &RANK,
960 WORK.
array(), &LWORK, &INFO);
964 return std::numeric_limits< C_FLOAT64 >::infinity();
967 LWORK = (
C_INT) WORK[0];
970 dgelsy_(&M, &N, &NRHS, JT.
array(), &LDA, X.
array(), &LDA, JPVT.
array(), &RCOND, &RANK,
971 WORK.
array(), &LWORK, &INFO);
975 return std::numeric_limits< C_FLOAT64 >::infinity();
1005 const CMetab * pMetab = NULL;
1013 for (; pAx != pAxEnd; ++pAx, ++pB, ++pCurrentState, ++pAtol, ++ppEntity)
1016 tmp = fabs(*pAx - *pB) /
std::max(fabs(*pCurrentState), *pAtol);
1017 RelativeDistance += tmp * tmp;
1019 tmp = fabs(*pAx - *pB);
1021 if ((pMetab = dynamic_cast< const CMetab * >(*ppEntity)) != NULL)
1026 AbsoluteDistance += tmp * tmp;
1030 isnan(RelativeDistance) ? std::numeric_limits< C_FLOAT64 >::infinity() : sqrt(RelativeDistance);
1032 isnan(AbsoluteDistance) ? std::numeric_limits< C_FLOAT64 >::infinity() : sqrt(AbsoluteDistance);
1034 Error =
std::max(RelativeDistance, AbsoluteDistance);
CCopasiDataModel * getObjectDataModel()
virtual bool setCallBack(CProcessReport *pCallBack)
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
CMatrix< C_FLOAT64 > * mpJacobianX
CCopasiProblem * getProblem()
C_FLOAT64 mMaxDurationBackward
int dgelsy_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *jpvt, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, integer *info)
size_t getNumIndependent() const
virtual bool isValidProblem(const CCopasiProblem *pProblem)
virtual bool elevateChildren()
virtual bool setName(const std::string &name)
void updateSimulatedValues(const bool &updateMoieties)
virtual size_t numRows() const
void initializeParameter()
virtual bool process(const bool &useInitialValues)
std::ostringstream mMethodLog
void resize(size_t size, const bool ©=false)
void setDuration(const C_FLOAT64 &duration)
CVector< C_FLOAT64 > mdxdt
const CCopasiMethod::SubType & getSubType() const
C_FLOAT64 targetFunction(const CVector< C_FLOAT64 > &particleFluxes)
bool isSteadyState(C_FLOAT64 value)
virtual bool progressItem(const size_t &handle)
virtual bool setModel(CModel *pModel)
C_FLOAT64 * endIndependent()
bool removeParameter(const std::string &name)
CVector< C_FLOAT64 > mXold
const C_FLOAT64 & getNumber2QuantityFactor() const
void setState(const CState &state)
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
void calculateJacobianX(const C_FLOAT64 &oldMaxRate)
virtual CSteadyStateMethod::ReturnCode processInternal()
CNewtonMethod::NewtonResultCode processNewton()
const bool & isAutonomous() const
const Value & getValue() const
CCopasiVectorN< CCopasiTask > * getTaskList()
void setStepNumber(const unsigned C_INT32 &stepNumber)
bool setValue(const std::string &name, const CType &value)
virtual CSteadyStateMethod::ReturnCode returnProcess(bool steadyStateFound)
CVector< C_FLOAT64 > mAtol
virtual bool finishItem(const size_t &handle)
CModelEntity ** beginIndependent()
void initializeParameter()
virtual bool isValidProblem(const CCopasiProblem *pProblem)
CCopasiParameter * getParameter(const std::string &name)
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
CTrajectoryTask * mpTrajectory
virtual bool initialize(const CSteadyStateProblem *pProblem)
const CSteadyStateProblem * mpProblem
unsigned C_INT32 mIterationLimit
void calculateDerivativesX(C_FLOAT64 *derivativesX)
const C_FLOAT64 & getValue() const
virtual void load(CReadConfig &configBuffer, CReadConfig::Mode mode=CReadConfig::SEARCH)
CNewtonMethod(const CCopasiContainer *pParent=NULL)
CCopasiMethod * getMethod()
void calculateDerivativesX()
CNewtonMethod::NewtonResultCode doNewtonStep(C_FLOAT64 ¤tValue)
virtual bool initialize(const CSteadyStateProblem *pProblem)
C_FLOAT64 solveJacobianXeqB(CVector< C_FLOAT64 > &X, const CVector< C_FLOAT64 > &B) const
const CStateTemplate & getStateTemplate() const
virtual size_t size() const
C_FLOAT64 * mpSSResolution
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
CVector< C_FLOAT64 > initializeAtolVector(const C_FLOAT64 &baseTolerance, const bool &reducedModel) const
C_INT32 getVariable(const std::string &name, const std::string &type, void *pout, CReadConfig::Mode mode=CReadConfig::NEXT)
CNewtonMethod::NewtonResultCode doIntegration(bool forward)
virtual size_t numCols() const
virtual bool setMethodType(const int &type)
CProcessReport * mpCallBack
CModel * getModel() const
C_FLOAT64 mMaxDurationForward
C_FLOAT64 * beginIndependent()