Bug 3260: Use the full model when required.

copasi · Jan 10, 2025 · 665ee37 · 665ee37
1 parent fecdf1e
commit 665ee37
Show file tree

Hide file tree

Showing 9 changed files with 38 additions and 46 deletions.
diff --git a/copasi/bindings/swig/CModel.i b/copasi/bindings/swig/CModel.i
@@ -1,4 +1,4 @@
-// Copyright (C) 2019 - 2024 by Pedro Mendes, Rector and Visitors of the 
+// Copyright (C) 2019 - 2025 by Pedro Mendes, Rector and Visitors of the 
 // University of Virginia, University of Heidelberg, and University 
 // of Connecticut School of Medicine. 
 // All rights reserved. 
@@ -113,7 +113,6 @@ typedef std::map< std::string, double > StringDoubleMap;
 %ignore CModel::processQueue;
 %ignore CModel::processRoots;
 %ignore CModel::calculateDerivatives;
-%ignore CModel::calculateDerivativesX;
 %ignore CModel::getListOfInitialRefreshes;
 %ignore CModel::getListOfSimulatedRefreshes;
 %ignore CModel::getListOfNonSimulatedRefreshes;

diff --git a/copasi/math/CMathContainer.cpp b/copasi/math/CMathContainer.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2019 - 2024 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2019 - 2025 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -800,7 +800,8 @@ CVector< C_FLOAT64 > CMathContainer::initializeAtolVector(const C_FLOAT64 & atol
   C_FLOAT64 * pAtol = Atol.array();
   C_FLOAT64 * pAtolEnd = pAtol + Atol.size();
   const C_FLOAT64 * pInitialValue = mInitialState.array() + mSize.nFixed;
-  const CMathObject * pObject = getMathObject(getState(reduced).array());
+  const CMathObject * pObject = getMathObject(pInitialValue);
+  const C_FLOAT64 Quantity2NumberFactor = * (C_FLOAT64 *) mpQuantity2NumberFactor->getValuePointer();
 
   for (; pAtol != pAtolEnd; ++pAtol, ++pObject, ++pInitialValue)
     {
@@ -812,12 +813,7 @@ CVector< C_FLOAT64 > CMathContainer::initializeAtolVector(const C_FLOAT64 & atol
         {
           case CMath::EntityType::Species:
           {
-            const CMetab * pMetab = static_cast< const CMetab * >(pObject->getDataObject()->getObjectParent());
-            std::map< const CDataObject *, CMathObject * >::const_iterator itFound
-              = mDataObject2MathObject.find(pMetab->getCompartment()->getInitialValueReference());
-
-            C_FLOAT64 Limit = fabs(* (C_FLOAT64 *) itFound->second->getValuePointer()) *
-                              * (C_FLOAT64 *) mpQuantity2NumberFactor->getValuePointer();
+            C_FLOAT64 Limit = fabs(*pObject->getCompartmentValue()) * Quantity2NumberFactor;
 
             if (InitialValue != 0.0)
               *pAtol *= std::min(Limit, InitialValue);
@@ -3143,6 +3139,8 @@ void CMathContainer::calculateJacobian(CMatrix< C_FLOAT64 > & jacobian,
                                        const bool & reduced,
                                        const bool & includeTime)
 {
+  CVector< C_FLOAT64 > State = getState(false);
+
   size_t Rows = getState(reduced).size() - mSize.nFixedEventTargets - 1;
   size_t Columns = getState(reduced).size() - mSize.nFixedEventTargets - (includeTime ? 0 : 1);
   jacobian.resize(Rows, Columns);
@@ -3212,6 +3210,7 @@ void CMathContainer::calculateJacobian(CMatrix< C_FLOAT64 > & jacobian,
     }
 
   updateSimulatedValues(reduced);
+  setState(State);
 }
 
 void CMathContainer::calculateJacobianDependencies(CMatrix< C_INT32 > & jacobianDependencies,

diff --git a/copasi/steadystate/CNewtonMethod.cpp b/copasi/steadystate/CNewtonMethod.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2019 - 2024 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2019 - 2025 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -318,7 +318,6 @@ CNewtonMethod::NewtonResultCode CNewtonMethod::doIntegration(bool forward)
           break;
         }
 
-      calculateDerivativesX();
       C_FLOAT64 value = targetFunction();
 
       if (isSteadyState(value))
@@ -437,7 +436,7 @@ CNewtonMethod::NewtonResultCode CNewtonMethod::doNewtonStep(C_FLOAT64 & currentV
 
   // DebugFile << "Iteration: " << k << std::endl;
 
-  calculateJacobian(currentValue, true);
+  calculateJacobian(true);
 
   if (CLeastSquareSolution::solve(*mpJacobian, mdxdt, mH) != mpJacobian->numCols())
     {
@@ -480,7 +479,7 @@ CNewtonMethod::NewtonResultCode CNewtonMethod::doNewtonStep(C_FLOAT64 & currentV
           (*pHit) *= 0.5;
         }
 
-      calculateDerivativesX();
+      mpContainer->updateSimulatedValues(true);
       newValue = targetFunction();
 
       // mpParentTask->output(COutputInterface::DURING);
@@ -491,7 +490,7 @@ CNewtonMethod::NewtonResultCode CNewtonMethod::doNewtonStep(C_FLOAT64 & currentV
       //discard the step
       memcpy(mpX, mXold.array(), mDimension * sizeof(C_FLOAT64));
 
-      calculateDerivativesX();
+      mpContainer->updateSimulatedValues(true);
       currentValue = targetFunction();
 
       if (mKeepProtocol)
@@ -547,7 +546,7 @@ CNewtonMethod::NewtonResultCode CNewtonMethod::processNewton()
                                       k,
                                       &mIterationLimit);
 
-  calculateDerivativesX();
+  mpContainer->updateSimulatedValues(true);
   C_FLOAT64 targetValue = targetFunction();
 
   {
@@ -640,11 +639,6 @@ CNewtonMethod::NewtonResultCode CNewtonMethod::processNewton()
   return result;
 }
 
-void CNewtonMethod::calculateDerivativesX()
-{
-  mpContainer->updateSimulatedValues(true);
-}
-
 bool CNewtonMethod::containsNaN() const
 {
   return !mpContainer->isStateValid();
@@ -665,11 +659,11 @@ C_FLOAT64 CNewtonMethod::targetFunction()
 {
   if (mTargetCriterion != eTargetCriterion::Rate)
     {
-      calculateJacobian(std::max(mTargetRate, mTargetDistance), true);
+      calculateJacobian(false);
     }
 
   // Assure that all values are updated.
-  mpContainer->updateSimulatedValues(true);
+  mpContainer->updateSimulatedValues(false);
   mpContainer->applyUpdateSequence(mUpdateConcentrations);
 
   mTargetRate = targetFunctionRate();
@@ -754,6 +748,7 @@ C_FLOAT64 CNewtonMethod::targetFunctionDistance()
 
   CMatrix< C_FLOAT64 > JacobianWithTime;
   mpContainer->calculateJacobian(JacobianWithTime, *mpDerivationFactor, true, !mpContainer->isAutonomous());
+  mpContainer->updateSimulatedValues(false);
 
   CLeastSquareSolution::ResultInfo Info = CLeastSquareSolution::solve(JacobianWithTime, mdxdt, mAtol, mCompartmentVolumes, mpContainer->getQuantity2NumberFactor(), Distance);
 

diff --git a/copasi/steadystate/CNewtonMethod.h b/copasi/steadystate/CNewtonMethod.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2019 - 2022 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2019 - 2025 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -211,8 +211,6 @@ class CNewtonMethod : public CSteadyStateMethod
    */
   C_FLOAT64 targetFunctionDistance();
 
-  void calculateDerivativesX();
-
   bool containsNaN() const;
 
   void cleanup();

diff --git a/copasi/steadystate/CSteadyStateMethod.cpp b/copasi/steadystate/CSteadyStateMethod.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2019 - 2024 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2019 - 2025 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -214,8 +214,6 @@ bool CSteadyStateMethod::isEquilibrium(const C_FLOAT64 & resolution) const
 
 bool CSteadyStateMethod::allPositive()
 {
-  mpContainer->updateSimulatedValues(true);
-
   const C_FLOAT64 * pValue = mContainerState.array();
   const C_FLOAT64 * pValueEnd = pValue + mContainerState.size();
   pValue += mpContainer->getCountFixedEventTargets() + 1; // + 1 for time
@@ -294,9 +292,13 @@ void CSteadyStateMethod::doJacobian(CMatrix< C_FLOAT64 > & jacobian,
                                     CMatrix< C_FLOAT64 > & jacobianX)
 {
   mpContainer->setState(mSteadyState);
-
   mpContainer->calculateJacobian(jacobian, *mpDerivationFactor, false);
+
+  mpContainer->setState(mSteadyState);
   mpContainer->calculateJacobian(jacobianX, *mpDerivationFactor, true);
+
+  mpContainer->setState(mSteadyState);
+  mpContainer->updateSimulatedValues(false);
 }
 
 C_FLOAT64 CSteadyStateMethod::getStabilityResolution()
@@ -306,7 +308,7 @@ C_FLOAT64 CSteadyStateMethod::getStabilityResolution()
   return *pTmp;
 }
 
-void CSteadyStateMethod::calculateJacobian(const C_FLOAT64 & oldMaxRate, const bool & reduced)
+void CSteadyStateMethod::calculateJacobian(const bool & reduced)
 {
   if (reduced)
     {

diff --git a/copasi/steadystate/CSteadyStateMethod.h b/copasi/steadystate/CSteadyStateMethod.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2019 - 2022 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2019 - 2025 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -229,7 +229,7 @@ class CSteadyStateMethod : public CCopasiMethod
    */
   bool allPositive();
 
-  void calculateJacobian(const C_FLOAT64 & oldMaxRate, const bool & reduced);
+  void calculateJacobian(const bool & reduced);
 };
 
 #include "CNewtonMethod.h"

diff --git a/copasi/steadystate/CSteadyStateTask.cpp b/copasi/steadystate/CSteadyStateTask.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2019 - 2022 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2019 - 2025 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -347,7 +347,7 @@ bool CSteadyStateTask::process(const bool & useInitialValues)
 #endif
 
   mpContainer->setState(mSteadyState);
-  mpContainer->updateSimulatedValues(true);
+  mpContainer->updateSimulatedValues(false);
   mpContainer->updateTransientDataValues();
   mpContainer->pushAllTransientValues();
 
@@ -419,7 +419,7 @@ std::ostream &operator<<(std::ostream &os, const CSteadyStateTask &A)
   os << std::endl;
 
   A.mpContainer->setState(A.mSteadyState);
-  A.mpContainer->updateSimulatedValues(true);
+  A.mpContainer->updateSimulatedValues(false);
   A.mpContainer->updateTransientDataValues();
   A.mpContainer->pushAllTransientValues();
 

diff --git a/copasi/tssanalysis/CILDMModifiedMethod.cpp b/copasi/tssanalysis/CILDMModifiedMethod.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2019 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2019 - 2025 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -229,7 +229,7 @@ void CILDMModifiedMethod::step(const double & deltaT)
 
   /* If complex eigenvalues */
   //BUG 873
-  if (mR(dim - 1, dim - 1) == mR(dim - 2 , dim - 2))
+  if (mR(dim - 1, dim - 1) == mR(dim - 2, dim - 2))
     if (dim == 2)
       {
         slow = dim;
@@ -288,7 +288,7 @@ void CILDMModifiedMethod::step(const double & deltaT)
           fast = fast - 1;
           slow = dim - fast;
 
-          if ((fast >= 1) && (mR(slow - 1, slow - 1) == mR(slow , slow)))
+          if ((fast >= 1) && (mR(slow - 1, slow - 1) == mR(slow, slow)))
             fast = fast - 1;
 
           slow = dim - fast;
@@ -399,7 +399,7 @@ void CILDMModifiedMethod::step(const double & deltaT)
                   re.resize(dim);
 
                   C_FLOAT64 eps;
-                  eps = 1 / fabs(mR(dim - k - 1 , dim - k - 1));
+                  eps = 1 / fabs(mR(dim - k - 1, dim - k - 1));
 
                   // stop criterion for slow reaction modes
 
@@ -524,7 +524,7 @@ void CILDMModifiedMethod::deuflhard_metab(C_INT & slow, C_INT & info)
   index_temp.resize(dim);
 
   C_FLOAT64 eps;
-  eps = 1 / fabs(mR(dim - fast , dim - fast));
+  eps = 1 / fabs(mR(dim - fast, dim - fast));
   //eps = fabs(mR(dim - fast - 1, dim - fast - 1)) / fabs(mR(dim - fast , dim - fast));
 
   mat_anal_fast_space(slow);
@@ -563,7 +563,6 @@ void CILDMModifiedMethod::deuflhard_metab(C_INT & slow, C_INT & info)
       x_help[j] = mY_initial[j] * mNumber2Concentration;
     }
 
-  // Model.calculateDerivativesX(dxdt.array());
   calculateDerivatives(x_help.array(), dxdt.array(), true);
 
   info_newton = 0;

diff --git a/copasi/utilities/CLeastSquareSolution.h b/copasi/utilities/CLeastSquareSolution.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2020 - 2022 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2020 - 2025 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -14,9 +14,9 @@ class CLeastSquareSolution
 public:
   struct ResultInfo
   {
-    size_t rank;
-    C_FLOAT64 relativeError;
-    C_FLOAT64 absoluteError;
+    size_t rank = 0;
+    C_FLOAT64 relativeError = 0;
+    C_FLOAT64 absoluteError = 0;
   };
 
   /**