From ebd3d10eac543fb4780edebfc17647b26ba22216 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Fri, 9 Apr 2021 16:31:09 -0400 Subject: [PATCH 1/5] construct condition specific plist in parameter mapping (#1487) * construct condition specific plist in parameter mapping * account for edata plist in gradient check * fix gradient check * fix gradient check --- python/amici/gradient_check.py | 6 ++++++ python/amici/parameter_mapping.py | 10 ++++++++++ 2 files changed, 16 insertions(+) diff --git a/python/amici/gradient_check.py b/python/amici/gradient_check.py index badf98afef..4490614331 100644 --- a/python/amici/gradient_check.py +++ b/python/amici/gradient_check.py @@ -67,6 +67,8 @@ 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) @@ -74,6 +76,8 @@ def check_finite_difference(x0: Sequence[float], 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, diff --git a/python/amici/parameter_mapping.py b/python/amici/parameter_mapping.py index 9f91c6b1a2..80baf2b4f1 100644 --- a/python/amici/parameter_mapping.py +++ b/python/amici/parameter_mapping.py @@ -236,12 +236,22 @@ def _get_par(model_par, value): scales = [petab_to_amici_scale(scale_map_sim_var[par_id]) for par_id in amici_model.getParameterIds()] + # plist + plist = [ + amici_model.getParameterIds().index(name) + for name, val in parameter_mapping.map_sim_var.items() + if not isinstance(val, numbers.Number) + ] + if parameters: edata.parameters = parameters if scales: edata.pscale = amici.parameterScalingFromIntVector(scales) + if plist: + edata.plist = plist + ########################################################################## # fixed parameters preequilibration if map_preeq_fix: From 3d61ac66d2a455f567d8982eb15dc5308ca3442a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Sat, 10 Apr 2021 16:26:59 -0400 Subject: [PATCH 2/5] Mapping plist (#1488) * construct condition specific plist in parameter mapping * account for edata plist in gradient check * fix gradient check * fix gradient check * fix plist ordering --- python/amici/parameter_mapping.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/python/amici/parameter_mapping.py b/python/amici/parameter_mapping.py index 80baf2b4f1..e7dcc0d2ec 100644 --- a/python/amici/parameter_mapping.py +++ b/python/amici/parameter_mapping.py @@ -238,9 +238,8 @@ def _get_par(model_par, value): # plist plist = [ - amici_model.getParameterIds().index(name) - for name, val in parameter_mapping.map_sim_var.items() - if not isinstance(val, numbers.Number) + ip for ip, par_id in enumerate(amici_model.getParameterIds()) + if isinstance(parameter_mapping.map_sim_var[par_id], str) ] if parameters: From 84675433d518f2e8c069e87a926c92857bc97816 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 13 Apr 2021 17:15:43 +0200 Subject: [PATCH 3/5] Fix serialization (#1490) Previously the archive was truncated under certain circumstances. The archive must be destroyed *before* returning. --- include/amici/serialization.h | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/include/amici/serialization.h b/include/amici/serialization.h index 7906a2cc48..0e2e6e0902 100644 --- 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> 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 serializeToStdVec(T const& data) { std::vector>> 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; } From 0a6a4ca9a0dcf38ba7276ac50486f21fb71d6fa5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Tue, 13 Apr 2021 12:19:42 -0400 Subject: [PATCH 4/5] 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 * fix doc * fix doc * fix doc? Co-authored-by: Daniel Weindl --- include/amici/model.h | 79 +++++++++++++++++++----- include/amici/rdata.h | 16 ++++- python/amici/gradient_check.py | 6 +- python/amici/numpy.py | 6 +- python/tests/test_pregenerated_models.py | 39 ++++++++++-- src/rdata.cpp | 62 +++++++++++++++---- 6 files changed, 165 insertions(+), 43 deletions(-) diff --git a/include/amici/model.h b/include/amici/model.h index 4d998e3c2c..1e47808443 100644 --- 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,7 +1344,7 @@ 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 @@ -1313,7 +1352,7 @@ class Model : public AbstractModel, public ModelDimensions { 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 index c0c42b78a5..ec8fe146b1 100644 --- a/include/amici/rdata.h +++ b/include/amici/rdata.h @@ -53,6 +53,8 @@ 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 ts, ModelDimensions const& model_dimensions, @@ -60,7 +62,8 @@ class ReturnData: public ModelDimensions { int newton_maxsteps, std::vector 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 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 index 4490614331..6055e2bdd2 100644 --- 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 index adace3b18f..5d8fd51d0e 100644 --- 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 index e339e03a5f..6bd63768b2 100755 --- 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') diff --git a/src/rdata.cpp b/src/rdata.cpp index 6e78a4b83d..50400d5af2 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -11,6 +11,7 @@ #include "amici/symbolic_functions.h" #include +#include namespace amici { @@ -21,8 +22,9 @@ ReturnData::ReturnData(Solver const &solver, const Model &model) solver.getNewtonMaxSteps(), model.getParameterScale(), model.o2mode, solver.getSensitivityOrder(), solver.getSensitivityMethod(), - solver.getReturnDataReportingMode(), model.hasQuadraticLLH()) { -} + solver.getReturnDataReportingMode(), model.hasQuadraticLLH(), + model.getAddSigmaResiduals(), + model.getMinimumSigmaResiduals()) {} ReturnData::ReturnData(std::vector ts, ModelDimensions const& model_dimensions, @@ -31,12 +33,14 @@ ReturnData::ReturnData(std::vector ts, std::vector pscale, SecondOrderMode o2mode, SensitivityOrder sensi, SensitivityMethod sensi_meth, RDataReporting rdrm, - bool quadratic_llh) + bool quadratic_llh, bool sigma_res, + realtype sigma_offset) : ModelDimensions(model_dimensions), ts(std::move(ts)), nx(nx_rdata), nxtrue(nxtrue_rdata), nplist(nplist), nmaxevent(nmaxevent), nt(nt), newton_maxsteps(newton_maxsteps), pscale(std::move(pscale)), o2mode(o2mode), sensi(sensi), - sensi_meth(sensi_meth), rdata_reporting(rdrm), x_solver_(nx_solver), + sensi_meth(sensi_meth), rdata_reporting(rdrm), sigma_res(sigma_res), + sigma_offset(sigma_offset), x_solver_(nx_solver), sx_solver_(nx_solver, nplist), x_rdata_(nx), sx_rdata_(nx, nplist), nroots_(ne) { switch (rdata_reporting) { @@ -72,7 +76,8 @@ void ReturnData::initializeResidualReporting(bool enable_res) { y.resize(nt * ny, 0.0); sigmay.resize(nt * ny, 0.0); if (enable_res) - res.resize(nt * nytrue, 0.0); + res.resize((sigma_res ? 2 : 1) * nt * nytrue, 0.0); + if ((sensi_meth == SensitivityMethod::forward && sensi >= SensitivityOrder::first) || sensi >= SensitivityOrder::second) { @@ -80,7 +85,7 @@ void ReturnData::initializeResidualReporting(bool enable_res) { sy.resize(nt * ny * nplist, 0.0); ssigmay.resize(nt * ny * nplist, 0.0); if (enable_res) - sres.resize(nt * nytrue * nplist, 0.0); + sres.resize((sigma_res ? 2 : 1) * nt * nytrue * nplist, 0.0); } } @@ -293,7 +298,7 @@ void ReturnData::getDataOutput(int it, Model &model, ExpData const *edata) { if (!isNaN(llh)) model.addObservableObjective(llh, it, x_solver_, *edata); fres(it, model, *edata); - fchi2(it); + fchi2(it, *edata); } if (sensi >= SensitivityOrder::first && nplist > 0) { @@ -664,9 +669,10 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model &model) { if (!sres.empty()) - for (int iyt = 0; iyt < nytrue * nt; ++iyt) + for (int ires = 0; ires < static_cast(res.size()); ++ires) for (int ip = 0; ip < nplist; ++ip) - sres.at((iyt * nplist + ip)) *= pcoefficient.at(ip); + sres.at((ires * nplist + ip)) *= pcoefficient.at(ip); + if(!FIM.empty()) for (int ip = 0; ip < nplist; ++ip) @@ -775,9 +781,8 @@ static realtype fres(realtype y, realtype my, realtype sigma_y) { return (y - my) / sigma_y; } -static realtype fsres(realtype y, realtype sy, realtype my, - realtype sigma_y, realtype ssigma_y) { - return (sy - ssigma_y * fres(y, my, sigma_y)) / sigma_y; +static realtype fres_error(realtype sigma_y, realtype sigma_offset) { + return sqrt(2*log(sigma_y) + sigma_offset); } void ReturnData::fres(const int it, Model &model, const ExpData &edata) { @@ -797,19 +802,35 @@ void ReturnData::fres(const int it, Model &model, const ExpData &edata) { continue; res.at(iyt) = amici::fres(y_it.at(iy), observedData[iy], sigmay_it.at(iy)); + if (sigma_res) + res.at(iyt + nt * nytrue) = fres_error(sigmay_it.at(iy), + sigma_offset); } } -void ReturnData::fchi2(const int it) { +void ReturnData::fchi2(const int it, const ExpData &edata) { if (res.empty() || isNaN(chi2)) return; for (int iy = 0; iy < nytrue; ++iy) { int iyt_true = iy + it * nytrue; chi2 += pow(res.at(iyt_true), 2); + if (sigma_res && edata.isSetObservedData(it, iy)) + chi2 += pow(res.at(iyt_true + nt * nytrue), 2) - sigma_offset; } } +static realtype fsres(realtype y, realtype sy, realtype my, + realtype sigma_y, realtype ssigma_y) { + return (sy - ssigma_y * fres(y, my, sigma_y)) / sigma_y; +} + +static realtype fsres_error(realtype sigma_y, realtype ssigma_y, + realtype sigma_offset) { + return ssigma_y / ( fres_error(sigma_y, sigma_offset) * sigma_y); +} + + void ReturnData::fsres(const int it, Model &model, const ExpData &edata) { if (sres.empty()) return; @@ -833,6 +854,14 @@ void ReturnData::fsres(const int it, Model &model, const ExpData &edata) { sres.at(idx) = amici::fsres(y_it.at(iy), sy_it.at(iy + ny * ip), observedData[iy], sigmay_it.at(iy), ssigmay_it.at(iy + ny * ip)); + if (sigma_res) { + int idx_res = + (iy + it * edata.nytrue() + edata.nytrue() * edata.nt()) * + nplist + ip; + sres.at(idx_res) = amici::fsres_error(sigmay_it.at(iy), + ssigmay_it.at(iy + ny * ip), + sigma_offset); + } } } } @@ -917,11 +946,18 @@ void ReturnData::fFIM(int it, Model &model, const ExpData &edata) { auto dy_i = sy_it.at(iy + ny * ip); auto ds_i = ssigmay_it.at(iy + ny * ip); auto sr_i = amici::fsres(y, dy_i, m, s, ds_i); + realtype sre_i; + if (sigma_res) + sre_i = amici::fsres_error(s, ds_i, sigma_offset); for (int jp = 0; jp < nplist; ++jp) { auto dy_j = sy_it.at(iy + ny * jp); auto ds_j = ssigmay_it.at(iy + ny * jp); auto sr_j = amici::fsres(y, dy_j, m, s, ds_j); FIM.at(ip + nplist * jp) += sr_i*sr_j; + if (sigma_res) { + auto sre_j = amici::fsres_error(s, ds_j, sigma_offset); + FIM.at(ip + nplist * jp) += sre_i * sre_j; + } /*+ ds_i*ds_j*(2*pow(r/pow(s,2.0), 2.0) - 1/pow(s,2.0));*/ } } From 669fd588823b5c9dd6a9cad0702b6af3eb3a9c1e Mon Sep 17 00:00:00 2001 From: FFroehlich Date: Tue, 13 Apr 2021 12:35:33 -0400 Subject: [PATCH 5/5] version bump and add changelog --- CHANGELOG.md | 9 +++++++++ version.txt | 2 +- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f3e981e425..3a0980a33e 100644 --- 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/version.txt b/version.txt index 4aa0906934..db07b54af2 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.11.15 +0.11.16