Merge branch 'develop' into pathlib

include/amici/serialization.h

@@ -62,8 +62,10 @@ void serialize(Archive &ar, amici::Solver &s, const unsigned int /*version*/) {
     ar &s.rtol_fsa_;
     ar &s.quad_atol_;
     ar &s.quad_rtol_;
+    ar &s.ss_tol_factor_;
     ar &s.ss_atol_;
     ar &s.ss_rtol_;
+    ar &s.ss_tol_sensi_factor_;
     ar &s.ss_atol_sensi_;
     ar &s.ss_rtol_sensi_;
     ar &s.maxsteps_;

include/amici/solver.h

@@ -44,7 +44,7 @@ namespace amici {
  *
  * NOTE: Any changes in data members here must be propagated to copy ctor,
  * equality operator, serialization functions in serialization.h, and
- * amici::hdf5::readSolverSettingsFromHDF5 in hdf5.cpp.
+ * amici::hdf5::(read/write)SolverSettings(From/To)HDF5 in hdf5.cpp.
  */
 class Solver {
   public:
@@ -384,6 +384,26 @@ class Solver {
      */
     void setAbsoluteToleranceQuadratures(double atol);
 
+    /**
+     * @brief returns the steady state simulation tolerance factor.
+     *
+     * Steady state simulation tolerances are the product of the simulation
+     * tolerances and this factor, unless manually set with
+     * `set(Absolute/Relative)ToleranceSteadyState()`.
+     * @return steady state simulation tolerance factor
+     */
+    double getSteadyStateToleranceFactor() const;
+
+    /**
+     * @brief set the steady state simulation tolerance factor.
+     *
+     * Steady state simulation tolerances are the product of the simulation
+     * tolerances and this factor, unless manually set with
+     * `set(Absolute/Relative)ToleranceSteadyState()`.
+     * @param factor tolerance factor (non-negative number)
+     */
+    void setSteadyStateToleranceFactor(double factor);
+
     /**
      * @brief returns the relative tolerance for the steady state problem
      * @return relative tolerance
@@ -408,6 +428,26 @@ class Solver {
      */
     void setAbsoluteToleranceSteadyState(double atol);
 
+    /**
+     * @brief returns the steady state sensitivity simulation tolerance factor.
+     *
+     * Steady state sensitivity simulation tolerances are the product of the
+     * sensitivity simulation tolerances and this factor, unless manually set
+     * with `set(Absolute/Relative)ToleranceSteadyStateSensi()`.
+     * @return steady state simulation tolerance factor
+     */
+    double getSteadyStateSensiToleranceFactor() const;
+
+    /**
+     * @brief set the steady state sensitivity simulation tolerance factor.
+     *
+     * Steady state sensitivity simulation tolerances are the product of the
+     * sensitivity simulation tolerances and this factor, unless manually set
+     * with `set(Absolute/Relative)ToleranceSteadyStateSensi()`.
+     * @param factor tolerance factor (non-negative number)
+     */
+    void setSteadyStateSensiToleranceFactor(double factor);
+
     /**
      * @brief returns the relative tolerance for the sensitivities of the
      * steady state problem
@@ -855,31 +895,31 @@ class Solver {
     std::vector<int> const& getLastOrder() const {
         return order_;
     }
-    
+
     /**
      * @brief Returns how convergence checks for steadystate computation are performed.
      * @return boolean flag indicating newton step (true) or the right hand side (false)
      */
     bool getNewtonStepSteadyStateCheck() const {
         return newton_step_steadystate_conv_;
     }
-    
+
     /**
      * @brief Returns how convergence checks for steadystate computation are performed.
      * @return boolean flag indicating state and sensitivity equations (true) or only state variables (false).
      */
     bool getSensiSteadyStateCheck() const {
         return check_sensi_steadystate_conv_;
     }
-    
+
     /**
      * @brief Sets how convergence checks for steadystate computation are performed.
      * @param flag boolean flag to pick newton step (true) or the right hand side (false, default)
      */
     void setNewtonStepSteadyStateCheck(bool flag) {
         newton_step_steadystate_conv_ = flag;
     }
-    
+
     /**
      * @brief Sets for which variables convergence checks for steadystate computation are performed.
      * @param flag boolean flag to pick state and sensitivity equations (true, default) or only state variables (false).
@@ -1753,23 +1793,29 @@ class Solver {
     /** relative tolerances for backward quadratures */
     realtype quad_rtol_ {1e-8};
 
+    /** steady state simulation tolerance factor */
+    realtype ss_tol_factor_ {1e2};
+
     /** absolute tolerances for steadystate computation */
     realtype ss_atol_ {NAN};
 
     /** relative tolerances for steadystate computation */
     realtype ss_rtol_ {NAN};
 
-    /** absolute tolerances for steadystate computation */
+    /** steady state sensitivity simulation tolerance factor */
+    realtype ss_tol_sensi_factor_ {1e2};
+
+    /** absolute tolerances for steadystate sensitivity computation */
     realtype ss_atol_sensi_ {NAN};
 
-    /** relative tolerances for steadystate computation */
+    /** relative tolerances for steadystate sensitivity computation */
     realtype ss_rtol_sensi_ {NAN};
 
     RDataReporting rdata_mode_ {RDataReporting::full};
-    
+
     /** whether newton step should be used for convergence steps */
     bool newton_step_steadystate_conv_ {false};
-    
+
     /** whether sensitivities should be checked for convergence to steadystate */
     bool check_sensi_steadystate_conv_ {true};
 

src/hdf5.cpp

@@ -657,6 +657,10 @@ void writeSolverSettingsToHDF5(Solver const& solver,
     H5LTset_attribute_double(file.getId(), hdf5Location.c_str(),
                              "quad_rtol", &dbuffer, 1);
 
+    dbuffer = solver.getSteadyStateToleranceFactor();
+    H5LTset_attribute_double(file.getId(), hdf5Location.c_str(),
+                             "ss_tol_factor", &dbuffer, 1);
+
     dbuffer = solver.getAbsoluteToleranceSteadyState();
     H5LTset_attribute_double(file.getId(), hdf5Location.c_str(),
                              "ss_atol", &dbuffer, 1);
@@ -665,6 +669,10 @@ void writeSolverSettingsToHDF5(Solver const& solver,
     H5LTset_attribute_double(file.getId(), hdf5Location.c_str(),
                              "ss_rtol", &dbuffer, 1);
 
+    dbuffer = solver.getSteadyStateSensiToleranceFactor();
+    H5LTset_attribute_double(file.getId(), hdf5Location.c_str(),
+                             "ss_tol_sensi_factor", &dbuffer, 1);
+
     dbuffer = solver.getAbsoluteToleranceSteadyStateSensi();
     H5LTset_attribute_double(file.getId(), hdf5Location.c_str(),
                              "ss_atol_sensi", &dbuffer, 1);
@@ -773,7 +781,6 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver,
                     getDoubleScalarAttribute(file, datasetPath, "rtol_fsa"));
     }
 
-
     if(attributeExists(file, datasetPath, "atolB")) {
         solver.setAbsoluteToleranceB(
                     getDoubleScalarAttribute(file, datasetPath, "atolB"));
@@ -794,6 +801,11 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver,
                     getDoubleScalarAttribute(file, datasetPath, "quad_rtol"));
     }
 
+    if(attributeExists(file, datasetPath, "ss_tol_factor")) {
+        solver.setSteadyStateToleranceFactor(
+                    getDoubleScalarAttribute(file, datasetPath, "ss_tol_factor"));
+    }
+
     if(attributeExists(file, datasetPath, "ss_atol")) {
         solver.setAbsoluteToleranceSteadyState(
                     getDoubleScalarAttribute(file, datasetPath, "ss_atol"));
@@ -804,6 +816,11 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver,
                     getDoubleScalarAttribute(file, datasetPath, "ss_rtol"));
     }
 
+    if(attributeExists(file, datasetPath, "ss_tol_sensi_factor")) {
+        solver.setSteadyStateSensiToleranceFactor(
+                    getDoubleScalarAttribute(file, datasetPath, "ss_tol_sensi_factor"));
+    }
+
     if(attributeExists(file, datasetPath, "ss_atol_sensi")) {
         solver.setAbsoluteToleranceSteadyStateSensi(
                     getDoubleScalarAttribute(file, datasetPath,

src/solver.cpp

@@ -29,8 +29,10 @@ Solver::Solver(const Solver &other)
       linsol_(other.linsol_), atol_(other.atol_), rtol_(other.rtol_),
       atol_fsa_(other.atol_fsa_), rtol_fsa_(other.rtol_fsa_),
       atolB_(other.atolB_), rtolB_(other.rtolB_), quad_atol_(other.quad_atol_),
-      quad_rtol_(other.quad_rtol_), ss_atol_(other.ss_atol_),
-      ss_rtol_(other.ss_rtol_), ss_atol_sensi_(other.ss_atol_sensi_),
+      quad_rtol_(other.quad_rtol_), ss_tol_factor_(other.ss_tol_factor_),
+      ss_atol_(other.ss_atol_), ss_rtol_(other.ss_rtol_),
+      ss_tol_sensi_factor_(other.ss_tol_sensi_factor_),
+      ss_atol_sensi_(other.ss_atol_sensi_),
       ss_rtol_sensi_(other.ss_rtol_sensi_), rdata_mode_(other.rdata_mode_),
       newton_step_steadystate_conv_(other.newton_step_steadystate_conv_),
       check_sensi_steadystate_conv_(other.check_sensi_steadystate_conv_),
@@ -781,8 +783,19 @@ void Solver::setAbsoluteToleranceQuadratures(const double atol) {
             applyTolerancesASA(iMem);
 }
 
+double Solver::getSteadyStateToleranceFactor() const {
+    return static_cast<double>(ss_tol_factor_);
+}
+
+void Solver::setSteadyStateToleranceFactor(const double ss_tol_factor) {
+    if (ss_tol_factor < 0)
+        throw AmiException("ss_tol_factor must be a non-negative number");
+
+    ss_tol_factor_ = static_cast<realtype>(ss_tol_factor);
+}
+
 double Solver::getRelativeToleranceSteadyState() const {
-    return static_cast<double>(isNaN(ss_rtol_) ? rtol_ : ss_rtol_);
+    return static_cast<double>(isNaN(ss_rtol_) ? rtol_ * ss_tol_factor_ : ss_rtol_);
 }
 
 void Solver::setRelativeToleranceSteadyState(const double rtol) {
@@ -793,7 +806,7 @@ void Solver::setRelativeToleranceSteadyState(const double rtol) {
 }
 
 double Solver::getAbsoluteToleranceSteadyState() const {
-    return static_cast<double>(isNaN(ss_atol_) ? atol_ : ss_atol_);
+    return static_cast<double>(isNaN(ss_atol_) ? atol_ * ss_tol_factor_ : ss_atol_);
 }
 
 void Solver::setAbsoluteToleranceSteadyState(const double atol) {
@@ -803,8 +816,19 @@ void Solver::setAbsoluteToleranceSteadyState(const double atol) {
     ss_atol_ = static_cast<realtype>(atol);
 }
 
+double Solver::getSteadyStateSensiToleranceFactor() const {
+    return static_cast<double>(ss_tol_sensi_factor_);
+}
+
+void Solver::setSteadyStateSensiToleranceFactor(const double ss_tol_sensi_factor) {
+    if (ss_tol_sensi_factor < 0)
+        throw AmiException("ss_tol_sensi_factor must be a non-negative number");
+
+    ss_tol_sensi_factor_ = static_cast<realtype>(ss_tol_sensi_factor);
+}
+
 double Solver::getRelativeToleranceSteadyStateSensi() const {
-    return static_cast<double>(isNaN(ss_rtol_sensi_) ? rtol_ : ss_rtol_sensi_);
+    return static_cast<double>(isNaN(ss_rtol_sensi_) ? rtol_ * ss_tol_sensi_factor_ : ss_rtol_sensi_);
 }
 
 void Solver::setRelativeToleranceSteadyStateSensi(const double rtol) {
@@ -815,7 +839,7 @@ void Solver::setRelativeToleranceSteadyStateSensi(const double rtol) {
 }
 
 double Solver::getAbsoluteToleranceSteadyStateSensi() const {
-    return static_cast<double>(isNaN(ss_atol_sensi_) ? atol_ : ss_atol_sensi_);
+    return static_cast<double>(isNaN(ss_atol_sensi_) ? atol_ * ss_tol_sensi_factor_ : ss_atol_sensi_);
 }
 
 void Solver::setAbsoluteToleranceSteadyStateSensi(const double atol) {

src/steadystateproblem.cpp

@@ -445,7 +445,7 @@ realtype SteadystateProblem::getWrmsNorm(const AmiVector &x,
     N_VAddConst(ewt.getNVector(), atol, ewt.getNVector());
     /* ewt = 1/ewt (ewt = 1/(rtol*x+atol)) */
     N_VInv(ewt.getNVector(), ewt.getNVector());
-    /* wrms = sqrt(sum((xdot/ewt)**2)) */
+    /* wrms = sqrt(sum((xdot/ewt)**2)/n) where n = size of state vector */
     return N_VWrmsNorm(const_cast<N_Vector>(xdot.getNVector()),
                        ewt.getNVector());
 }

tests/cpp/unittests/testMisc.cpp

@@ -448,6 +448,59 @@ testSolverGetterSetters(CVodeSolver solver,
     ASSERT_EQ(solver.getAbsoluteToleranceSteadyState(), tol);
 }
 
+TEST_F(SolverTest, SteadyStateToleranceFactor)
+{
+    CVodeSolver s;
+    // test with unset steadystate tolerances
+    ASSERT_DOUBLE_EQ(
+        s.getRelativeToleranceSteadyState(),
+        s.getSteadyStateToleranceFactor() * s.getRelativeTolerance());
+    ASSERT_DOUBLE_EQ(
+        s.getAbsoluteToleranceSteadyState(),
+        s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance());
+    ASSERT_DOUBLE_EQ(
+        s.getRelativeToleranceSteadyStateSensi(),
+        s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance());
+    ASSERT_DOUBLE_EQ(
+        s.getAbsoluteToleranceSteadyState(),
+        s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance());
+
+    // test with changed steadystate tolerance factor
+    s.setSteadyStateToleranceFactor(5);
+    ASSERT_DOUBLE_EQ(
+        s.getRelativeToleranceSteadyState(),
+        s.getSteadyStateToleranceFactor() * s.getRelativeTolerance());
+    ASSERT_DOUBLE_EQ(
+        s.getAbsoluteToleranceSteadyState(),
+        s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance());
+    s.setSteadyStateSensiToleranceFactor(5);
+    ASSERT_DOUBLE_EQ(
+        s.getRelativeToleranceSteadyStateSensi(),
+        s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance());
+    ASSERT_DOUBLE_EQ(
+        s.getAbsoluteToleranceSteadyState(),
+        s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance());
+
+
+    // test with steadystate tolerance override tolerance factor
+    s.setRelativeToleranceSteadyState(2);
+    ASSERT_NE(s.getRelativeToleranceSteadyState(),
+        s.getSteadyStateToleranceFactor() * s.getRelativeTolerance());
+    ASSERT_EQ(s.getRelativeToleranceSteadyState(), 2);
+    s.setAbsoluteToleranceSteadyState(3);
+    ASSERT_NE(s.getAbsoluteToleranceSteadyState(),
+              s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance());
+    ASSERT_EQ(s.getAbsoluteToleranceSteadyState(), 3);
+    s.setRelativeToleranceSteadyStateSensi(4);
+    ASSERT_NE(s.getRelativeToleranceSteadyStateSensi(),
+              s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance());
+    ASSERT_EQ(s.getRelativeToleranceSteadyStateSensi(), 4);
+    s.setAbsoluteToleranceSteadyStateSensi(5);
+    ASSERT_NE(s.getAbsoluteToleranceSteadyStateSensi(),
+              s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance());
+    ASSERT_EQ(s.getAbsoluteToleranceSteadyStateSensi(), 5);
+}
+
 class AmiVectorTest : public ::testing::Test {
   protected:
     std::vector<double> vec1{ 1, 2, 4, 3 };