Merge pull request #243 from ICB-DCM/feature_more_adjoint_options

Feature more adjoint options

Author: paulstapor

@amioption/amioption.m

@@ -13,6 +13,12 @@
         rtol = 1e-8;
         % maximum number of integration steps
         maxsteps = 1e4;
+        % absolute integration tolerace
+        quad_atol = 1e-12;
+        % relative integration tolerace
+        quad_rtol = 1e-8;
+        % maximum number of integration steps
+        maxstepsB = 0;
         % index of parameters for which the sensitivities are computed
         sens_ind = double.empty();
         % index of states for which positivity should be enforced

include/udata.h

@@ -2,6 +2,7 @@
 #define AMICI_UDATA_H
 
 #include "include/amici_defines.h"
+#include "include/symbolic_functions.h" //getNaN
 #include <cmath>
 #include <vector>
 
@@ -319,6 +320,15 @@ class UserData {
     /** maximum number of allowed integration steps */
     int maxsteps = 0;
 
+    /** absolute tolerances for backward quadratures */
+    double quad_atol = 1e-12;
+    
+    /** relative tolerances for backward quadratures */
+    double quad_rtol = 1e-8;
+    
+    /** maximum number of allowed integration steps for backward problem */
+    int maxstepsB = 0;
+    
     /** flag indicating whether sensitivities are supposed to be computed */
     AMICI_sensi_order sensi = AMICI_SENSI_ORDER_NONE;
 

src/amici_interface_matlab.cpp

@@ -124,6 +124,9 @@ UserData *userDataFromMatlabCall(const mxArray *prhs[], int nrhs) {
         readOptionScalarDouble(atol);
         readOptionScalarDouble(rtol);
         readOptionScalarInt(maxsteps, int);
+        readOptionScalarDouble(quad_atol);
+        readOptionScalarDouble(quad_rtol);
+        readOptionScalarInt(maxstepsB, int);
         readOptionScalarInt(lmm, LinearMultistepMethod);
         readOptionScalarInt(iter, NonlinearSolverIteration);
         readOptionScalarInt(interpType, InterpolationType);

src/amici_solver.cpp

@@ -148,18 +148,25 @@ void Solver::setupAMIB(BackwardProblem *bwd, const UserData *udata, Model *model
     AMISetUserDataB(bwd->getwhich(), model);
 
     /* Number of maximal internal steps */
-    AMISetMaxNumStepsB(bwd->getwhich(), 100 * udata->maxsteps);
+    AMISetMaxNumStepsB(bwd->getwhich(), (udata->maxstepsB == 0) ? udata->maxsteps * 100 : udata->maxstepsB);
 
     setLinearSolverB(udata, model, bwd->getwhich());
 
     /* Initialise quadrature calculation */
     qbinit(bwd->getwhich(), bwd->getxQBptr());
-
+    
+    double quad_rtol = isNaN(udata->quad_rtol) ? udata->rtol : udata->quad_rtol;
+    double quad_atol = isNaN(udata->quad_atol) ? udata->atol : udata->quad_atol;
+    
     /* Enable Quadrature Error Control */
-    AMISetQuadErrConB(bwd->getwhich(), TRUE);
+    if (std::isinf(quad_atol) || std::isinf(quad_rtol)) {
+        AMISetQuadErrConB(bwd->getwhich(), FALSE);
+    } else {
+        AMISetQuadErrConB(bwd->getwhich(), TRUE);
+    }
 
-    AMIQuadSStolerancesB(bwd->getwhich(), RCONST(udata->rtol),
-                                  RCONST(udata->atol));
+    AMIQuadSStolerancesB(bwd->getwhich(), RCONST(quad_rtol),
+                         RCONST(quad_atol));
 
     AMISetStabLimDetB(bwd->getwhich(), udata->getStabilityLimitFlag());
 }

src/forwardproblem.cpp

@@ -610,9 +610,12 @@ void ForwardProblem::getDataSensisFSA(int it) {
             }
         }
     }
-    model->fsy(it, rdata);
-    if (edata) {
-        model->fsJy(it, dJydx, rdata);
+    
+    if (model->ny > 0) {
+        model->fsy(it, rdata);
+        if (edata) {
+            model->fsJy(it, dJydx, rdata);
+        }
     }
 }
 

src/steadystateproblem.cpp

@@ -294,4 +294,4 @@ void SteadystateProblem::getNewtonSimulation(const UserData *udata,
     }
 }
 
-} // namespace amici
+} // namespace amici

src/udata.cpp

@@ -37,7 +37,10 @@ UserData::UserData(const UserData &other) : UserData(other.np(), other.nk(), oth
     sensi = other.sensi;
     atol = other.atol;
     rtol = other.rtol;
+    quad_atol = other.quad_atol;
+    quad_rtol = other.quad_rtol;
     maxsteps = other.maxsteps;
+    maxstepsB = other.maxstepsB;
     newton_maxsteps = other.newton_maxsteps;
     newton_maxlinsteps = other.newton_maxlinsteps;
     newton_preeq = other.newton_preeq;
@@ -237,6 +240,9 @@ void UserData::print() const {
     printf("atol: %e\n", atol);
     printf("rtol: %e\n", rtol);
     printf("maxsteps: %d\n", maxsteps);
+    printf("quad_atol: %e\n", quad_atol);
+    printf("quad_rtol: %e\n", quad_rtol);
+    printf("maxstepsB: %d\n", maxstepsB);
     printf("newton_maxsteps: %d\n", newton_maxsteps);
     printf("newton_maxlinsteps: %d\n", newton_maxlinsteps);
     printf("ism: %d\n", ism);

tests/cpputest/unittests/testsSerialization.cpp

@@ -95,8 +95,8 @@ void checkReturnDataEqual(amici::ReturnData const& r, amici::ReturnData const& s
     CHECK_EQUAL(*r.chi2, *s.chi2);
     CHECK_EQUAL(*r.status, *s.status);
 
-    checkEqualArray(r.sllh, s.sllh, r.nplist, 1e-16, 1e-16, "sllh");
-    checkEqualArray(r.s2llh, s.s2llh, r.nplist * r.nplist, 1e-16, 1e-16, "s2llh");
+    checkEqualArray(r.sllh, s.sllh, r.nplist, 1e-5, 1e-5, "sllh");
+    checkEqualArray(r.s2llh, s.s2llh, r.nplist * r.nplist, 1e-5, 1e-5, "s2llh");
 }