Merge pull request #1491 from AMICI-dev/release_0_11_16

Release 0.11.16
AMICI-dev · Apr 13, 2021 · 57a7b11 · 57a7b11
2 parents 4e671a5 + 669fd58
commit 57a7b11
Show file tree

Hide file tree

Showing 10 changed files with 201 additions and 54 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,15 @@
 
 ## v0.X Series
 
+### v0.11.16 (2021-04-13)
+
+Fixes:
+* Fixed serialization bug  (#1490)
+
+New:
+* Construction of condition specific plist for parameter mappings (#1487, #1488)
+* **Add support for error residuals** (#1489)
+
 ### v0.11.15 (2021-03-31)
 
 Fixes:

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/include/amici/serialization.h b/include/amici/serialization.h
@@ -299,12 +299,13 @@ T deserializeFromChar(const char *buffer, int size) {
     T data;
 
     try {
+        // archive must be destroyed BEFORE returning
         ba::binary_iarchive iar(s);
         iar >> data;
-        return data;
     } catch(ba::archive_exception const& e) {
         throw AmiException("Deserialization from char failed: %s", e.what());
     }
+    return data;
 }
 
 /**
@@ -325,14 +326,14 @@ std::string serializeToString(T const& data) {
     bio::stream<bio::back_insert_device<std::string>> os(inserter);
 
     try {
+        // archive must be destroyed BEFORE returning
         ba::binary_oarchive oar(os);
-
         oar << data;
-
-        return serialized;
     } catch(ba::archive_exception const& e) {
         throw AmiException("Serialization to string failed: %s", e.what());
     }
+
+    return serialized;
 }
 
 /**
@@ -354,14 +355,14 @@ std::vector<char> serializeToStdVec(T const& data) {
             std::vector<char>>> os(buffer);
 
     try{
+        // archive must be destroyed BEFORE returning
         ba::binary_oarchive oar(os);
-
         oar << data;
-
-        return buffer;
     } catch(ba::archive_exception const& e) {
         throw AmiException("Serialization to std::vector failed: %s", e.what());
     }
+
+    return buffer;
 }
 
 /**
@@ -382,15 +383,15 @@ T deserializeFromString(std::string const& serialized) {
     T deserialized;
 
     try{
+        // archive must be destroyed BEFORE returning
         ba::binary_iarchive iar(os);
-
         iar >> deserialized;
-
-        return deserialized;
     } catch(ba::archive_exception const& e) {
         throw AmiException("Deserialization from std::string failed: %s",
                            e.what());
     }
+
+    return deserialized;
 }
 
 

diff --git a/python/amici/gradient_check.py b/python/amici/gradient_check.py
@@ -67,13 +67,17 @@ def check_finite_difference(x0: Sequence[float],
     og_sensitivity_order = solver.getSensitivityOrder()
     og_parameters = model.getParameters()
     og_plist = model.getParameterList()
+    if edata:
+        og_eplist = edata.plist
 
     # sensitivity
     p = copy.deepcopy(x0)
     plist = [ip]
 
     model.setParameters(p)
     model.setParameterList(plist)
+    if edata:
+        edata.plist = plist
 
     # simulation with gradient
     if int(og_sensitivity_order) < int(SensitivityOrder.first):
@@ -122,6 +126,8 @@ def check_finite_difference(x0: Sequence[float],
     solver.setSensitivityOrder(og_sensitivity_order)
     model.setParameters(og_parameters)
     model.setParameterList(og_plist)
+    if edata:
+        edata.plist = og_eplist
 
 
 def check_derivatives(model: Model,
@@ -190,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):