Add residuals for quadratic nllhs with parameter dependent sigma (#1489)

* fixes and add tests * update doc * fix doc * fix doc? * fix doc * fixup * fixup and fix tests * Apply suggestions from code review Co-authored-by: Daniel Weindl <[email protected]> * fix doc * fix doc * fix doc? Co-authored-by: Daniel Weindl <[email protected]>
AMICI-dev · Apr 13, 2021 · 0a6a4ca · 0a6a4ca
1 parent 8467543
commit 0a6a4ca
Show file tree

Hide file tree

Showing 6 changed files with 165 additions and 43 deletions.
diff --git a/include/amici/model.h b/include/amici/model.h
@@ -642,6 +642,45 @@ class Model : public AbstractModel, public ModelDimensions {
             throw AmiException("Mismatch in conservation law sensitivity size");
         state_ = state;
     };
+
+    /**
+     * @brief Sets the estimated lower boundary for sigma_y. When :meth:`setAddSigmaResiduals` is
+     * activated, this lower boundary must ensure that log(sigma) + min_sigma > 0.
+     * @param min_sigma lower boundary
+     */
+    void setMinimumSigmaResiduals(double min_sigma) {
+        min_sigma_ = min_sigma;
+    }
+
+    /**
+     * @brief Gets the specified estimated lower boundary for sigma_y.
+     * @return lower boundary
+     */
+    realtype getMinimumSigmaResiduals() const {
+        return min_sigma_;
+    }
+
+    /**
+     * @brief Specifies whether residuals should be added to account for parameter dependent sigma.
+     *
+     * If set to true, additional residuals of the form \f$ \sqrt{\log(\sigma) + C} \f$ will be added.
+     * This enables least-squares optimization for variables with Gaussian noise assumption and parameter
+     * dependent standard deviation sigma. The constant \f$ C \f$ can be set via
+     * :meth:`setMinimumSigmaResiduals`.
+     *
+     * @param sigma_res if true, additional residuals are added
+     */
+    void setAddSigmaResiduals(bool sigma_res) {
+        sigma_res_ = sigma_res;
+    }
+
+    /**
+     * @brief Checks whether residuals should be added to account for parameter dependent sigma.
+     * @return sigma_res
+     */
+    bool getAddSigmaResiduals() const {
+        return sigma_res_;
+    }
 
     /**
      * @brief Get the list of parameters for which sensitivities are computed.
@@ -767,7 +806,7 @@ class Model : public AbstractModel, public ModelDimensions {
     /**
      * @brief Get sensitivity of time-resolved observables.
      *
-     * Total derivative \f$sy = dydx * sx + dydp\f$
+     * Total derivative \f$ sy = dydx * sx + dydp\f$
      * (only for forward sensitivities).
      * @param sy buffer (shape `ny` x `nplist`, row-major)
      * @param t Timpoint
@@ -801,7 +840,7 @@ class Model : public AbstractModel, public ModelDimensions {
                                        const int it, const ExpData *edata);
 
     /**
-     * @brief Add time-resolved measurement negative log-likelihood \f$Jy\f$.
+     * @brief Add time-resolved measurement negative log-likelihood \f$ Jy \f$.
      * @param Jy Buffer (shape 1)
      * @param it Timepoint index
      * @param x State variables
@@ -812,7 +851,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Add sensitivity of time-resolved measurement negative log-likelihood
-     * \f$Jy\f$.
+     * \f$ Jy \f$.
      *
      * @param sllh First-order buffer (shape `nplist`)
      * @param s2llh Second-order buffer (shape `nJ - 1` x `nplist`, row-major)
@@ -829,7 +868,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Add sensitivity of time-resolved measurement negative
-     * log-likelihood \f$Jy\f$.
+     * log-likelihood \f$ Jy \f$.
      *
      * Partial derivative (to be used with adjoint sensitivities).
      *
@@ -847,7 +886,7 @@ class Model : public AbstractModel, public ModelDimensions {
                                                   const ExpData &edata);
 
     /**
-     * @brief Get state sensitivity of the negative loglikelihood \f$Jy\f$,
+     * @brief Get state sensitivity of the negative loglikelihood \f$ Jy \f$,
      * partial derivative (to be used with adjoint sensitivities).
      *
      * @param dJydx Output buffer (shape `nJ` x `nx_solver`, row-major)
@@ -978,7 +1017,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Add sensitivity of time-resolved measurement negative
-     * log-likelihood \f$Jy\f$.
+     * log-likelihood \f$ Jy \f$.
      *
      * Total derivative (to be used with forward sensitivities).
      *
@@ -1001,7 +1040,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Add sensitivity of time-resolved measurement negative
-     * log-likelihood \f$Jy\f$.
+     * log-likelihood \f$ Jy \f$.
      *
      * Partial derivative (to be used with adjoint sensitivities).
      *
@@ -1022,7 +1061,7 @@ class Model : public AbstractModel, public ModelDimensions {
                                              const ExpData &edata);
 
     /**
-     * @brief State sensitivity of the negative loglikelihood \f$Jz\f$.
+     * @brief State sensitivity of the negative loglikelihood \f$ Jz \f$.
      *
      * Partial derivative (to be used with adjoint sensitivities).
      *
@@ -1305,15 +1344,15 @@ class Model : public AbstractModel, public ModelDimensions {
     void fy(realtype t, const AmiVector &x);
 
     /**
-     * @brief Compute partial derivative of observables \f$y\f$ w.r.t. model
+     * @brief Compute partial derivative of observables \f$ y \f$ w.r.t. model
      * parameters `p`.
      * @param t Current timepoint
      * @param x Current state
      */
     void fdydp(realtype t, const AmiVector &x);
 
     /**
-     * @brief Compute partial derivative of observables \f$y\f$ w.r.t. state
+     * @brief Compute partial derivative of observables \f$ y \f$ w.r.t. state
      * variables `x`.
      * @param t Current timepoint
      * @param x Current state
@@ -1336,7 +1375,7 @@ class Model : public AbstractModel, public ModelDimensions {
     void fdsigmaydp(int it, const ExpData *edata);
 
     /**
-     * @brief Compute negative log-likelihood of measurements \f$y\f$.
+     * @brief Compute negative log-likelihood of measurements \f$ y \f$.
      *
      * @param Jy Variable to which llh will be added
      * @param it Timepoint index
@@ -1347,7 +1386,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Compute partial derivative of time-resolved measurement negative
-     * log-likelihood \f$Jy\f$.
+     * log-likelihood \f$ Jy \f$.
      * @param it timepoint index
      * @param x state variables
      * @param edata Pointer to experimental data
@@ -1365,7 +1404,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Compute sensitivity of time-resolved measurement negative
-     * log-likelihood \f$Jy\f$ w.r.t. parameters for the given timepoint.
+     * log-likelihood \f$ Jy \f$ w.r.t. parameters for the given timepoint.
      * @param it timepoint index
      * @param x state variables
      * @param edata pointer to experimental data instance
@@ -1374,7 +1413,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Sensitivity of time-resolved measurement negative log-likelihood
-     * \f$Jy\f$ w.r.t. state variables.
+     * \f$ Jy \f$ w.r.t. state variables.
      * @param it Timepoint index
      * @param x State variables
      * @param edata Pointer to experimental data instance
@@ -1428,7 +1467,7 @@ class Model : public AbstractModel, public ModelDimensions {
     void fdrzdp(int ie, realtype t, const AmiVector &x);
 
     /**
-     * @brief Compute sensitivity of event-resolved measurements \f$rz\f$ w.r.t.
+     * @brief Compute sensitivity of event-resolved measurements \f$ rz \f$ w.r.t.
      * model states `x`.
      * @param ie Event index
      * @param t Current timepoint
@@ -1469,7 +1508,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Compute partial derivative of event measurement negative
-     * log-likelihood \f$Jz\f$.
+     * log-likelihood \f$ Jz \f$.
      * @param ie Event index
      * @param nroots Event index
      * @param t Current timepoint
@@ -1481,7 +1520,7 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /**
      * @brief Compute sensitivity of event measurement negative log-likelihood
-     * \f$Jz\f$ w.r.t. standard deviation sigmaz.
+     * \f$ Jz \f$ w.r.t. standard deviation sigmaz.
      * @param ie Event index
      * @param nroots Event index
      * @param t Current timepoint
@@ -1704,6 +1743,12 @@ class Model : public AbstractModel, public ModelDimensions {
      * checked for finiteness
      */
     bool always_check_finite_ {false};
+
+    /** indicates whether sigma residuals are to be added for every datapoint  */
+    bool sigma_res_ {false};
+
+    /** offset to ensure positivity of sigma residuals, only has an effect when `sigma_res_` is `true`  */
+    realtype min_sigma_ {50.0};
 
   private:
     /** Sparse dwdp implicit temporary storage (shape `ndwdp`) */

diff --git a/include/amici/rdata.h b/include/amici/rdata.h
@@ -53,14 +53,17 @@ class ReturnData: public ModelDimensions {
      * @param rdrm see amici::Solver::rdata_reporting
      * @param quadratic_llh whether model defines a quadratic nllh and
      * computing res, sres and FIM makes sense
+     * @param sigma_res indicates whether additional residuals are to be added for each sigma
+     * @param sigma_offset offset to ensure real-valuedness of sigma residuals
      */
     ReturnData(std::vector<realtype> ts,
                ModelDimensions const& model_dimensions,
                int nplist, int nmaxevent, int nt,
                int newton_maxsteps,
                std::vector<ParameterScaling> pscale, SecondOrderMode o2mode,
                SensitivityOrder sensi, SensitivityMethod sensi_meth,
-               RDataReporting rdrm, bool quadratic_llh);
+               RDataReporting rdrm, bool quadratic_llh, bool sigma_res,
+               realtype sigma_offset);
 
     /**
      * @brief constructor that uses information from model and solver to
@@ -372,9 +375,15 @@ class ReturnData: public ModelDimensions {
     template <class Archive>
     friend void boost::serialization::serialize(Archive &ar, ReturnData &r,
                                                 unsigned int version);
+
+    /** boolean indicating whether residuals for standard deviations have been added */
+    bool sigma_res;
 
   protected:
-
+
+    /** offset for sigma_residuals */
+    realtype sigma_offset;
+
     /** timepoint for model evaluation*/
     realtype t_;
 
@@ -517,8 +526,9 @@ class ReturnData: public ModelDimensions {
     /**
      * @brief Chi-squared function
      * @param it time index
+     * @param edata ExpData instance containing observable data
      */
-    void fchi2(int it);
+    void fchi2(int it, const ExpData &edata);
 
     /**
      * @brief Residual sensitivity function

diff --git a/python/amici/gradient_check.py b/python/amici/gradient_check.py
@@ -196,18 +196,18 @@ def check_derivatives(model: Model,
 
     if 'ssigmay' in rdata.keys() \
             and rdata['ssigmay'] is not None \
-            and rdata['ssigmay'].any():
+            and rdata['ssigmay'].any() and not model.getAddSigmaResiduals():
         leastsquares_applicable = False
 
     if check_least_squares and leastsquares_applicable:
         fields += ['res', 'y']
 
         check_results(rdata, 'FIM',
-                      np.dot(rdata['sres'].transpose(), rdata['sres']),
+                      np.dot(rdata['sres'].T, rdata['sres']),
                       assert_fun,
                       1e-8, 1e-4)
         check_results(rdata, 'sllh',
-                      -np.dot(rdata['res'].transpose(), rdata['sres']),
+                      -np.dot(rdata['res'].T, rdata['sres']),
                       assert_fun,
                       1e-8, 1e-4)
     for ip, pval in enumerate(p):

diff --git a/python/amici/numpy.py b/python/amici/numpy.py
@@ -195,8 +195,10 @@ def __init__(self, rdata: Union[ReturnDataPtr, ReturnData]):
             'sllh': [rdata.nplist],
             's2llh': [rdata.np, rdata.nplist],
 
-            'res': [rdata.nt * rdata.nytrue],
-            'sres': [rdata.nt * rdata.nytrue, rdata.nplist],
+            'res': [rdata.nt * rdata.nytrue *
+                    (2 if rdata.sigma_res else 1)],
+            'sres': [rdata.nt * rdata.nytrue *
+                     (2 if rdata.sigma_res else 1), rdata.nplist],
             'FIM': [rdata.nplist, rdata.nplist],
 
             # diagnosis

diff --git a/python/tests/test_pregenerated_models.py b/python/tests/test_pregenerated_models.py
@@ -3,19 +3,17 @@
 """Run simulations with Matlab-AMICI pre-generated models and verify using
 saved expectations."""
 
-import sys
 import h5py
 import amici
-import unittest
-import importlib
 import os
-import copy
 from amici.gradient_check import check_derivatives, check_results
 import pytest
 
+import numpy as np
+
 
 options_file = os.path.join(os.path.dirname(__file__), '..', '..',
-                                     'tests', 'cpputest', 'testOptions.h5')
+                            'tests', 'cpputest', 'testOptions.h5')
 expected_results_file = os.path.join(os.path.dirname(__file__), '..', '..',
                                      'tests', 'cpputest', 'expectedResults.h5')
 expected_results = h5py.File(expected_results_file, 'r')
@@ -147,6 +145,8 @@ def test_pregenerated_model(sub_test, case):
         with pytest.raises(RuntimeError):
             solver.setSensitivityMethod(amici.SensitivityMethod.adjoint)
 
+    chi2_ref = rdata.chi2
+
     # test likelihood mode
     solver.setReturnDataReportingMode(amici.RDataReporting.likelihood)
     rdata = amici.runAmiciSimulation(model, solver, edata)
@@ -155,6 +155,35 @@ def test_pregenerated_model(sub_test, case):
         fields=['t', 'llh', 'sllh', 's2llh', 'FIM'], **verify_simulation_opts
     )
 
+    # test sigma residuals
+
+    if model_name == 'model_jakstat_adjoint' and \
+            solver.getSensitivityMethod() != amici.SensitivityMethod.adjoint:
+        model.setAddSigmaResiduals(True)
+        solver.setReturnDataReportingMode(amici.RDataReporting.full)
+        rdata = amici.runAmiciSimulation(model, solver, edata)
+        # check whether activation changes chi2
+        assert chi2_ref != rdata.chi2
+
+        if edata and solver.getSensitivityMethod() and \
+                solver.getSensitivityOrder() and len(model.getParameterList()):
+            check_derivatives(model, solver, edata, assert_fun,
+                              **check_derivative_opts)
+
+        chi2_ref = rdata.chi2
+        res_ref = rdata.res
+
+        model.setMinimumSigmaResiduals(100)
+        rdata = amici.runAmiciSimulation(model, solver, edata)
+        # check whether changing the minimum changes res but not chi2
+        assert np.isclose(chi2_ref, rdata.chi2)
+        assert not np.allclose(res_ref, rdata.res)
+
+        model.setMinimumSigmaResiduals(-10)
+        rdata = amici.runAmiciSimulation(model, solver, edata)
+        # check whether having a bad minimum results in nan chi2
+        assert np.isnan(rdata.chi2)
+
     with pytest.raises(RuntimeError):
         model.getParameterByName('thisParameterDoesNotExist')