From 58c95a901fe06dec7f4a162a6b5f10a8d62cf41a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 30 Jan 2019 14:34:39 +0100 Subject: [PATCH 1/8] fix(ci) Fix premature exit of scripts/runNotebook.sh --- scripts/runNotebook.sh | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/scripts/runNotebook.sh b/scripts/runNotebook.sh index fc94c6f993..90d3952c60 100755 --- a/scripts/runNotebook.sh +++ b/scripts/runNotebook.sh @@ -12,9 +12,11 @@ runNotebook () { tempfile=$(mktemp) jupyter nbconvert --debug --stdout --execute --ExecutePreprocessor.timeout=300 --to markdown $@ &> $tempfile ret=$? - if [[ $ret != 0 ]]; then cat $tempfile; fi + if [[ $ret != 0 ]]; then + cat $tempfile + exit $ret + fi rm $tempfile - exit $ret set -e } From fa616fdb808547062251465663461fc533124221 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 31 Jan 2019 08:10:35 +0100 Subject: [PATCH 2/8] Fix(deploy) Update pyenv shims to find twine --- scripts/deployPyPi.sh | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/scripts/deployPyPi.sh b/scripts/deployPyPi.sh index e43c07f38b..804fdb5b25 100755 --- a/scripts/deployPyPi.sh +++ b/scripts/deployPyPi.sh @@ -8,6 +8,10 @@ SCRIPT_PATH=$(dirname $BASH_SOURCE) AMICI_PATH=$(cd $SCRIPT_PATH/.. && pwd) pip3 install twine +# in case we are running with pyenv, we need to update pyenv shims after installing packages with binaries +if [[ -z "${PYENV_VERSION}" ]]; then + pyenv rehash +fi # authentication via env variables set in travis twine upload $(ls -t ${AMICI_PATH}/build/python/amici-*.tar.gz | head -1) From 15bf57066bbaf151f657ace2a64fb3c6d428efd0 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 22 Jan 2019 15:43:41 +0100 Subject: [PATCH 3/8] Redirect C/C++ output in stdout is redirected (e.g. in ipython notebooks) (Closes #456) --- python/amici/__init__.py | 29 +++++++++++++++++++++-------- python/amici/setup.template.py | 7 ++++--- python/sdist/setup.py | 9 +++++++-- 3 files changed, 32 insertions(+), 13 deletions(-) mode change 100644 => 100755 python/sdist/setup.py diff --git a/python/amici/__init__.py b/python/amici/__init__.py index 79f6d9e240..9f9676af8d 100644 --- a/python/amici/__init__.py +++ b/python/amici/__init__.py @@ -27,6 +27,17 @@ import os import re +import sys +from contextlib import suppress + +# redirect C/C++ stdout to python stdout if python stdout is redirected, +# e.g. in ipython notebook +sys_pipes = suppress +if sys.stdout != sys.__stdout__: + try: + from wurlitzer import sys_pipes + except ModuleNotFoundError: + pass hdf5_enabled = False has_clibs = False @@ -103,12 +114,13 @@ def runAmiciSimulation(model, solver, edata=None): ReturnData object with simulation results Raises: - + """ if edata and edata.__class__.__name__ == 'ExpDataPtr': edata = edata.get() - rdata = amici.runAmiciSimulation(solver.get(), edata, model.get()) + with sys_pipes(): + rdata = amici.runAmiciSimulation(solver.get(), edata, model.get()) return rdataToNumPyArrays(rdata) @@ -134,7 +146,7 @@ def ExpData(*args): return amici.ExpData(args[0].get(), *args[:1]) else: return amici.ExpData(*args) - + def runAmiciSimulations(model, solver, edata_list, num_threads=1): """ Convenience wrapper for loops of amici.runAmiciSimulation @@ -152,9 +164,10 @@ def runAmiciSimulations(model, solver, edata_list, num_threads=1): Raises: """ - edata_ptr_vector = amici.ExpDataPtrVector(edata_list) - rdata_ptr_list = amici.runAmiciSimulations(solver.get(), - edata_ptr_vector, - model.get(), - num_threads) + with sys_pipes(): + edata_ptr_vector = amici.ExpDataPtrVector(edata_list) + rdata_ptr_list = amici.runAmiciSimulations(solver.get(), + edata_ptr_vector, + model.get(), + num_threads) return [numpy.rdataToNumPyArrays(r) for r in rdata_ptr_list] diff --git a/python/amici/setup.template.py b/python/amici/setup.template.py index 2ae5beeb0a..5ddf14d00b 100644 --- a/python/amici/setup.template.py +++ b/python/amici/setup.template.py @@ -48,7 +48,7 @@ def getAmiciLibs(): libraries.extend(['hdf5_hl_cpp', 'hdf5_hl', 'hdf5_cpp', 'hdf5']) sources = ['swig/TPL_MODELNAME.i', *getModelSources()] - + # Remove the "-Wstrict-prototypes" compiler option, which isn't valid for # C++ to fix warnings. @@ -67,7 +67,7 @@ def getAmiciLibs(): model_module = Extension('TPL_MODELNAME._TPL_MODELNAME', sources=sources, include_dirs=[os.getcwd(), - os.path.join(amici_path, 'include'), + os.path.join(amici_path, 'include'), os.path.join(amici_path, 'ThirdParty/sundials/include'), os.path.join(amici_path, 'ThirdParty/SuiteSparse/include'), *h5pkgcfg['include_dirs'], @@ -98,7 +98,8 @@ def getAmiciLibs(): packages=find_packages(), # TODO: should specify amici version with which the model was generated install_requires=['amici'], - python_requires='>=3', + extras_require={'wurlitzer': ['wurlitzer']}, + python_requires='>=3.6', package_data={ }, zip_safe = False, diff --git a/python/sdist/setup.py b/python/sdist/setup.py old mode 100644 new mode 100755 index 3151b7fe3a..d8dbf91e93 --- a/python/sdist/setup.py +++ b/python/sdist/setup.py @@ -1,6 +1,6 @@ """Setuptools file for creating AMICI module -This file is based on setuptools alone and does not require CMake. +This file is based on setuptools alone and does not require CMake. All sources are compiled anew. This file expects to be run from within its directory. @@ -272,8 +272,13 @@ def main(): ], packages=find_packages(), package_dir={'amici': 'amici'}, - install_requires=['sympy', 'python-libsbml', 'h5py', 'pandas', 'setuptools>=40.6.3'], + install_requires=['sympy', + 'python-libsbml', + 'h5py', + 'pandas', + 'setuptools>=40.6.3'], python_requires='>=3.6', + extras_require={'wurlitzer': ['wurlitzer']}, package_data={ 'amici': ['amici/include/amici/*', 'src/*template*', From ea38f43d1868df5b985aa86dcd874743d381c110 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Sat, 2 Feb 2019 11:54:56 +0100 Subject: [PATCH 4/8] feature(core) Add option to rethrow AMICI exception (Closes #552) - Introduce ConditionContext for setting and restoring condition-specific data - Don't invalidate ReturnData::llh and ~chi2 on backwards integration problems - Doxygen to header --- include/amici/amici.h | 3 +- include/amici/edata.h | 46 ++++++++- include/amici/rdata.h | 52 +++++++--- src/amici.cpp | 49 ++++------ src/edata.cpp | 133 +++++++++++++++++--------- src/forwardproblem.cpp | 1 + src/rdata.cpp | 61 ++++-------- tests/cpputest/steadystate/tests1.cpp | 22 ++++- 8 files changed, 228 insertions(+), 139 deletions(-) diff --git a/include/amici/amici.h b/include/amici/amici.h index 034f966c06..a94605cb4f 100644 --- a/include/amici/amici.h +++ b/include/amici/amici.h @@ -27,9 +27,10 @@ extern msgIdAndTxtFp warnMsgIdAndTxt; * @param solver Solver instance * @param edata pointer to experimental data object * @param model model specification object + * @param rethrow rethrow integration exceptions? * @return rdata pointer to return data object */ -std::unique_ptr runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model); +std::unique_ptr runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model, bool rethrow=false); /*! * runAmiciSimulations does the same as runAmiciSimulation, but for multiple ExpData instances. diff --git a/include/amici/edata.h b/include/amici/edata.h index 56cb1e44d1..484f418df0 100644 --- a/include/amici/edata.h +++ b/include/amici/edata.h @@ -38,7 +38,7 @@ class ExpData { */ ExpData(int nytrue, int nztrue, int nmaxevent, std::vector const& ts); - + /** * constructor that initializes timepoints and fixed parameters from vectors * @@ -351,18 +351,18 @@ class ExpData { realtype t_presim = 0; protected: - + /** * resizes observedData, observedDataStdDev, observedEvents and * observedEventsStdDev */ void applyDimensions(); - + /** * resizes observedData and observedDataStdDev */ void applyDataDimension(); - + /** * resizes observedEvents and observedEventsStdDev */ @@ -425,6 +425,44 @@ void checkSigmaPositivity(std::vector const& sigmaVector, const char * */ void checkSigmaPositivity(const realtype sigma, const char *sigmaName); + +/** + * @brief The ConditionContext class applies condition-specific amici::Model + * settings and restores them when going out of scope + */ +class ConditionContext { +public: + /** + * @brief Apply condition-specific settings from edata to model while + * keeping a backup of the original values. + * @param model + * @param edata + */ + ConditionContext(Model *model, const ExpData *edata = nullptr); + + ~ConditionContext(); + + /** + * @brief Apply condition-specific settings from edata to the + * constructor-supplied model, not changing the settings which were + * backed-up in the constructor call. + * @param edata + */ + void applyCondition(const ExpData *edata); + + /** + * @brief Restore original settings on constructor-supplied amici::Model. + * Will be called during destruction. Explicit call is generally not + * necessary. + */ + void restore(); + +private: + Model *model = nullptr; + std::vector originalFixedParameters; + std::vector originalTimepoints; +}; + } // namespace amici #endif /* AMICI_EDATA_H */ diff --git a/include/amici/rdata.h b/include/amici/rdata.h index 5f58be2914..9d77ebeac8 100644 --- a/include/amici/rdata.h +++ b/include/amici/rdata.h @@ -19,14 +19,16 @@ void serialize(Archive &ar, amici::ReturnData &u, const unsigned int version); namespace amici { -/** @brief class that stores all data which is later returned by the mex - * function +/** @brief Stores all data to be returned by amici::runAmiciSimulation. * * NOTE: multidimensional arrays are stored in row-major order * (FORTRAN-style) */ class ReturnData { public: + /** + * @brief default constructor + */ ReturnData(); /** @@ -59,6 +61,13 @@ class ReturnData { int nt, int newton_maxsteps, std::vector pscale, SecondOrderMode o2mode, SensitivityOrder sensi, SensitivityMethod sensi_meth); + /** + * @brief constructor that uses information from model and solver to + * appropriately initialize fields + * @param solver solver + * @param model pointer to model specification object + * bool + */ ReturnData(Solver const& solver, const Model *model); ~ReturnData() = default; @@ -68,10 +77,31 @@ class ReturnData { */ void initializeObjectiveFunction(); + /** + * @brief Set likelihood, state variables, outputs and respective + * sensitivities to NaN (typically after integration failure) + * @param t time of integration failure + */ void invalidate(const realtype t); + + /** + * @brief Set likelihood and chi2 to NaN + * (typically after integration failure) + */ void invalidateLLH(); - void + /** + * @brief Set likelihood sensitivities to NaN + * (typically after integration failure) + */ + void invalidateSLLH(); + + /** + * @brief applies the chain rule to account for parameter transformation + * in the sensitivities of simulation results + * @param model Model from which the ReturnData was obtained + */ + void applyChainRuleFactorToSimulationResults(const Model *model); /** timepoints (dimension: nt) */ @@ -130,14 +160,14 @@ class ReturnData { /** parameter derivative of observable standard deviation (dimension: nt x * nplist x ny, row-major) */ std::vector ssigmay; - + /** observable (dimension: nt*ny, row-major) */ std::vector res; /** parameter derivative of residual (dimension: nt*ny x nplist, * row-major) */ std::vector sres; - + /** fisher information matrix (dimension: nplist x nplist, * row-major) */ std::vector FIM; @@ -184,28 +214,28 @@ class ReturnData { /** number of linear steps by Newton step for steady state problem. this will only be filled for iterative solvers (length = newton_maxsteps * 2) */ std::vector newton_numlinsteps; - + /** time at which steadystate was reached in the simulation based approach */ realtype t_steadystate = NAN; - + /** weighted root-mean-square of the rhs when steadystate was reached*/ realtype wrms_steadystate = NAN; - + /** weighted root-mean-square of the rhs when steadystate was reached*/ realtype wrms_sensi_steadystate = NAN; - + /** initial state (dimension: nx) */ std::vector x0; - + /** preequilibration steady state found by Newton solver (dimension: nx) */ std::vector x_ss; /** initial sensitivities (dimension: nplist x nx, row-major) */ std::vector sx0; - + /** preequilibration sensitivities found by Newton solver (dimension: nplist x nx, row-major) */ std::vector sx_ss; diff --git a/src/amici.cpp b/src/amici.cpp index 57803f8655..031a3422e3 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -37,50 +37,38 @@ namespace amici { msgIdAndTxtFp errMsgIdAndTxt = &printErrMsgIdAndTxt; /** warnMsgIdAndTxt is a function pointer for printWarnMsgIdAndTxt */ msgIdAndTxtFp warnMsgIdAndTxt = &printWarnMsgIdAndTxt; - -std::unique_ptr runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model) { - + +std::unique_ptr runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model, bool rethrow) { std::unique_ptr rdata; - - auto originalFixedParameters = model.getFixedParameters(); // to restore after simulation - auto originalTimepoints = model.getTimepoints(); - if(edata) { - if(!edata->fixedParameters.empty()) { - // fixed parameter in model are superseded by those provided in edata - if(edata->fixedParameters.size() != (unsigned) model.nk()) - throw AmiException("Number of fixed parameters (%d) in model does not match ExpData (%zd).", - model.nk(), edata->fixedParameters.size()); - model.setFixedParameters(edata->fixedParameters); - } - if(edata->nt()) { - // fixed parameter in model are superseded by those provided in edata - model.setTimepoints(edata->getTimepoints()); - } - } - + + /* Applies condition-specific model settings and restores them when going + * out of scope */ + ConditionContext conditionContext(&model, edata); + try{ rdata = std::unique_ptr(new ReturnData(solver,&model)); + if (model.nx_solver <= 0) { - model.setFixedParameters(originalFixedParameters); - model.setTimepoints(originalTimepoints); return rdata; } - + auto fwd = std::unique_ptr(new ForwardProblem(rdata.get(),edata,&model,&solver)); fwd->workForwardProblem(); auto bwd = std::unique_ptr(new BackwardProblem(fwd.get())); bwd->workBackwardProblem(); - + rdata->status = AMICI_SUCCESS; } catch (amici::IntegrationFailure const& ex) { rdata->invalidate(ex.time); rdata->status = ex.error_code; + if(rethrow) throw; amici::warnMsgIdAndTxt("AMICI:mex:simulation","AMICI forward simulation failed at t = %f:\n%s\n",ex.time,ex.what()); } catch (amici::IntegrationFailureB const& ex) { - rdata->invalidateLLH(); + rdata->invalidateSLLH(); rdata->status = ex.error_code; + if(rethrow) throw; amici::warnMsgIdAndTxt( "AMICI:mex:simulation", "AMICI backward simulation failed when trying to solve until t = %f" @@ -89,19 +77,14 @@ std::unique_ptr runAmiciSimulation(Solver &solver, const ExpData *ed } catch (amici::AmiException const& ex) { rdata->invalidate(model.t0()); rdata->status = AMICI_ERROR; + if(rethrow) throw; amici::warnMsgIdAndTxt("AMICI:mex:simulation","AMICI simulation failed:\n%s\nError occured in:\n%s",ex.what(),ex.getBacktrace()); - } catch (std::exception const& ex) { - model.setFixedParameters(originalFixedParameters); - model.setTimepoints(originalTimepoints); - throw; } catch (...) { - model.setFixedParameters(originalFixedParameters); - model.setTimepoints(originalTimepoints); throw std::runtime_error("Unknown internal error occured!"); } - model.setFixedParameters(originalFixedParameters); - model.setTimepoints(originalTimepoints); + rdata->applyChainRuleFactorToSimulationResults(&model); + return rdata; } diff --git a/src/edata.cpp b/src/edata.cpp index 649d0b7b46..b997e485b3 100644 --- a/src/edata.cpp +++ b/src/edata.cpp @@ -25,7 +25,7 @@ ExpData::ExpData(int nytrue, int nztrue, int nmaxevent, { applyDimensions(); } - + ExpData::ExpData(int nytrue, int nztrue, int nmaxevent, std::vector const& ts_, std::vector const& fixedParameters_ @@ -53,28 +53,28 @@ ExpData::ExpData(int nytrue, int nztrue, int nmaxevent, ExpData::ExpData(Model const& model) : ExpData(model.nytrue, model.nztrue, model.nMaxEvent(), model.getTimepoints(), model.getFixedParameters()) {} - + ExpData::ExpData(ReturnData const& rdata, realtype sigma_y, realtype sigma_z) : ExpData(rdata, std::vector(rdata.nytrue*rdata.nt, sigma_y), std::vector(rdata.nztrue*rdata.nmaxevent, sigma_z)) {} - + ExpData::ExpData(ReturnData const& rdata, std::vector sigma_y, std::vector sigma_z) : ExpData(rdata.nytrue, rdata.nztrue, rdata.nmaxevent, rdata.ts) { if (sigma_y.size() != (unsigned) nytrue_ && sigma_y.size() != (unsigned) nytrue_*nt()) throw AmiException("Dimension of sigma_y must be %d or %d, was %d", nytrue_, nytrue_*nt(), sigma_y.size()); - + if (sigma_z.size() != (unsigned) nztrue_ && sigma_z.size() != (unsigned) nztrue_*nmaxevent_) throw AmiException("Dimension of sigma_z must be %d or %d, was %d", nztrue_, nztrue_*nmaxevent_, sigma_z.size()); - + std::random_device rd{}; std::mt19937 gen{rd()}; - + realtype sigma; - + checkSigmaPositivity(sigma_y, "sigma_y"); checkSigmaPositivity(sigma_z, "sigma_z"); - + for (int iy = 0; iy < nytrue_; ++iy) { for (int it = 0; it < nt(); ++it) { sigma = sigma_y.size() == (unsigned) nytrue_ ? sigma_y.at(iy) : sigma_y.at(iy + nytrue_ * it); @@ -83,7 +83,7 @@ ExpData::ExpData(ReturnData const& rdata, std::vector sigma_y, observedDataStdDev.at(iy + nytrue_ * it) = sigma; } } - + for (int iz = 0; iz < nztrue_; ++iz) { for (int ie = 0; ie < nmaxevent_; ++ie) { sigma = sigma_z.size() == (unsigned) nztrue_ ? sigma_z.at(iz) : sigma_z.at(iz + nztrue_ * ie); @@ -93,14 +93,14 @@ ExpData::ExpData(ReturnData const& rdata, std::vector sigma_y, } } } - + void ExpData::setTimepoints(const std::vector &ts) { if (!std::is_sorted(ts.begin(), ts.end())) throw AmiException("Encountered non-monotonic timepoints, please order timepoints such that they are monotonically increasing!"); this->ts = ts; applyDataDimension(); } - + std::vector const& ExpData::getTimepoints() const { return ts; } @@ -108,77 +108,77 @@ std::vector const& ExpData::getTimepoints() const { int ExpData::nt() const { return ts.size(); } - + realtype ExpData::getTimepoint(int it) const { return ts.at(it); } - + void ExpData::setObservedData(const std::vector &observedData) { checkDataDimension(observedData, "observedData"); - + if (observedData.size() == (unsigned) nt()*nytrue_) this->observedData = observedData; else if (observedData.empty()) this->observedData.clear(); } - + void ExpData::setObservedData(const std::vector &observedData, int iy) { if (observedData.size() != (unsigned) nt()) throw AmiException("Input observedData did not match dimensions nt (%i), was %i", nt(), observedData.size()); - + for (int it = 0; it < nt(); ++it) this->observedData.at(iy + it*nytrue_) = observedData.at(it); } - + bool ExpData::isSetObservedData(int it, int iy) const { return !observedData.empty() && !isNaN(observedData.at(it * nytrue_ + iy)); } - + std::vector const& ExpData::getObservedData() const { return observedData; } - + const realtype *ExpData::getObservedDataPtr(int it) const { if (!observedData.empty()) return &observedData.at(it*nytrue_); - + return nullptr; } void ExpData::setObservedDataStdDev(const std::vector &observedDataStdDev) { checkDataDimension(observedDataStdDev, "observedDataStdDev"); checkSigmaPositivity(observedDataStdDev, "observedDataStdDev"); - + if (observedDataStdDev.size() == (unsigned) nt()*nytrue_) this->observedDataStdDev = observedDataStdDev; else if (observedDataStdDev.empty()) this->observedDataStdDev.clear(); } - + void ExpData::setObservedDataStdDev(const realtype stdDev) { checkSigmaPositivity(stdDev, "stdDev"); std::fill(observedDataStdDev.begin() ,observedDataStdDev.end(), stdDev); } - + void ExpData::setObservedDataStdDev(const std::vector &observedDataStdDev, int iy) { if (observedDataStdDev.size() != (unsigned) nt()) throw AmiException("Input observedDataStdDev did not match dimensions nt (%i), was %i", nt(), observedDataStdDev.size()); checkSigmaPositivity(observedDataStdDev, "observedDataStdDev"); - + for (int it = 0; it < nt(); ++it) this->observedDataStdDev.at(iy + it*nytrue_) = observedDataStdDev.at(it); } - + void ExpData::setObservedDataStdDev(const realtype stdDev, int iy) { checkSigmaPositivity(stdDev, "stdDev"); for (int it = 0; it < nt(); ++it) observedDataStdDev.at(iy + it*nytrue_) = stdDev; } - + bool ExpData::isSetObservedDataStdDev(int it, int iy) const { return !observedDataStdDev.empty() && !isNaN(observedDataStdDev.at(it * nytrue_ + iy)); } - + std::vector const& ExpData::getObservedDataStdDev() const { return observedDataStdDev; } @@ -186,28 +186,28 @@ std::vector const& ExpData::getObservedDataStdDev() const { const realtype *ExpData::getObservedDataStdDevPtr(int it) const { if (!observedDataStdDev.empty()) return &observedDataStdDev.at(it*nytrue_); - + return nullptr; } - + void ExpData::setObservedEvents(const std::vector &observedEvents) { checkEventsDimension(observedEvents, "observedEvents"); - + if (observedEvents.size() == (unsigned) nmaxevent_*nztrue_) this->observedEvents = observedEvents; else if (observedEvents.empty()) this->observedEvents.clear(); } - + void ExpData::setObservedEvents(const std::vector &observedEvents, int iz) { if (observedEvents.size() != (unsigned) nmaxevent_) { throw AmiException("Input observedEvents did not match dimensions nmaxevent (%i), was %i", nmaxevent_, observedEvents.size()); } - + for (int ie = 0; ie < nmaxevent_; ++ie) this->observedEvents.at(iz + ie*nztrue_) = observedEvents.at(ie); } - + bool ExpData::isSetObservedEvents(int ie, int iz) const { return !observedEvents.size() && !isNaN(observedEvents.at(ie * nztrue_ + iz)); } @@ -219,20 +219,20 @@ std::vector const& ExpData::getObservedEvents() const { const realtype *ExpData::getObservedEventsPtr(int ie) const { if (!observedEvents.empty()) return &observedEvents.at(ie*nztrue_); - + return nullptr; } - + void ExpData::setObservedEventsStdDev(const std::vector &observedEventsStdDev) { checkEventsDimension(observedEventsStdDev, "observedEventsStdDev"); checkSigmaPositivity(observedEventsStdDev, "observedEventsStdDev"); - + if (observedEventsStdDev.size() == (unsigned) nmaxevent_*nztrue_) this->observedEventsStdDev = observedEventsStdDev; else if (observedEventsStdDev.empty()) this->observedEventsStdDev.clear(); } - + void ExpData::setObservedEventsStdDev(const realtype stdDev) { checkSigmaPositivity(stdDev, "stdDev"); std::fill(observedEventsStdDev.begin() ,observedEventsStdDev.end(), stdDev); @@ -242,22 +242,22 @@ void ExpData::setObservedEventsStdDev(const std::vector &observedEvent if (observedEventsStdDev.size() != (unsigned) nmaxevent_) throw AmiException("Input observedEventsStdDev did not match dimensions nmaxevent (%i), was %i", nmaxevent_, observedEventsStdDev.size()); checkSigmaPositivity(observedEventsStdDev, "observedEventsStdDev"); - + for (int ie = 0; ie < nmaxevent_; ++ie) this->observedEventsStdDev.at(iz + ie*nztrue_) = observedEventsStdDev.at(ie); } void ExpData::setObservedEventsStdDev(const realtype stdDev, int iz) { checkSigmaPositivity(stdDev, "stdDev"); - + for (int ie = 0; ie < nmaxevent_; ++ie) observedEventsStdDev.at(iz + ie*nztrue_) = stdDev; } - + bool ExpData::isSetObservedEventsStdDev(int ie, int iz) const { if (!observedEventsStdDev.empty()) // avoid out of bound memory access return !isNaN(observedEventsStdDev.at(ie * nztrue_ + iz)); - + return false; } @@ -268,10 +268,10 @@ std::vector const& ExpData::getObservedEventsStdDev() const { const realtype *ExpData::getObservedEventsStdDevPtr(int ie) const { if (!observedEventsStdDev.empty()) return &observedEventsStdDev.at(ie*nztrue_); - + return nullptr; } - + void ExpData::applyDimensions() { applyDataDimension(); applyEventDimension(); @@ -281,17 +281,17 @@ void ExpData::applyDataDimension() { observedData.resize(nt()*nytrue_, getNaN()); observedDataStdDev.resize(nt()*nytrue_, getNaN()); } - + void ExpData::applyEventDimension() { observedEvents.resize(nmaxevent_*nztrue_, getNaN()); observedEventsStdDev.resize(nmaxevent_*nztrue_, getNaN()); } - + void ExpData::checkDataDimension(std::vector input, const char *fieldname) const { if (input.size() != (unsigned) nt()*nytrue_ && !input.empty()) throw AmiException("Input %s did not match dimensions nt (%i) x nytrue (%i), was %i", fieldname, nt(), nytrue_, input.size()); } - + void ExpData::checkEventsDimension(std::vector input, const char *fieldname) const { if (input.size() != (unsigned) nmaxevent_*nztrue_ && !input.empty()) throw AmiException("Input %s did not match dimensions nt (%i) x nytrue (%i), was %i", fieldname, nmaxevent_, nztrue_, input.size()); @@ -321,5 +321,44 @@ int ExpData::nmaxevent() const { return nmaxevent_; } - + +ConditionContext::ConditionContext(Model *model, const ExpData *edata) + : model(model), + originalFixedParameters(model->getFixedParameters()), + originalTimepoints(model->getTimepoints()) +{ + applyCondition(edata); +} + +ConditionContext::~ConditionContext() +{ + restore(); +} + +void ConditionContext::applyCondition(const ExpData *edata) +{ + if(!edata) + return; + + if(!edata->fixedParameters.empty()) { + // fixed parameter in model are superseded by those provided in edata + if(edata->fixedParameters.size() != (unsigned) model->nk()) + throw AmiException("Number of fixed parameters (%d) in model does not match ExpData (%zd).", + model->nk(), edata->fixedParameters.size()); + model->setFixedParameters(edata->fixedParameters); + } + + if(edata->nt()) { + // fixed parameter in model are superseded by those provided in edata + model->setTimepoints(edata->getTimepoints()); + } +} + +void ConditionContext::restore() +{ + model->setFixedParameters(originalFixedParameters); + model->setTimepoints(originalTimepoints); +} + + } // namespace amici diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index ffe4bb41a7..91467c0e98 100644 --- a/src/forwardproblem.cpp +++ b/src/forwardproblem.cpp @@ -170,6 +170,7 @@ void ForwardProblem::workForwardProblem() { // set likelihood if (!edata) { rdata->invalidateLLH(); + rdata->invalidateSLLH(); } storeJacobianAndDerivativeInReturnData(); diff --git a/src/rdata.cpp b/src/rdata.cpp index d2b3df39db..0bd0c1abbe 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -11,9 +11,6 @@ namespace amici { ReturnData::ReturnData() - /** - * @brief default constructor - */ : np(0), nk(0), nx(0), nx_solver(0), nxtrue(0), ny(0), nytrue(0), nz(0), nztrue(0), ne(0), nJ(0), nplist(0), nmaxevent(0), nt(0), newton_maxsteps(0), pscale(std::vector(0, ParameterScaling::none)), o2mode(SecondOrderMode::none), @@ -27,13 +24,6 @@ ReturnData::ReturnData(Solver const& solver, const Model *model) solver.getNewtonMaxSteps(), model->getParameterScale(), model->o2mode, solver.getSensitivityOrder(), static_cast(solver.getSensitivityMethod())) { - /** - * @brief constructor that uses information from model and solver to - * appropriately initialize fields - * @param solver solver - * @param model pointer to model specification object - * bool - */ } @@ -73,7 +63,7 @@ ReturnData::ReturnData( numerrtestfails.resize(nt, 0); numnonlinsolvconvfails.resize(nt, 0); order.resize(nt, 0); - + if (sensi_meth == SensitivityMethod::adjoint && sensi >= SensitivityOrder::first) { numstepsB.resize(nt, 0); numrhsevalsB.resize(nt, 0); @@ -83,9 +73,9 @@ ReturnData::ReturnData( } x0.resize(nx, getNaN()); - + x_ss.resize(nx, getNaN()); - + res.resize(nt * nytrue, 0.0); llh = getNaN(); @@ -94,18 +84,18 @@ ReturnData::ReturnData( sllh.resize(nplist, getNaN()); sx0.resize(nx * nplist, getNaN()); sx_ss.resize(nx * nplist, getNaN()); - + if (sensi_meth == SensitivityMethod::forward || sensi >= SensitivityOrder::second){ // for second order we can fill in from the augmented states sx.resize(nt * nx * nplist, 0.0); sy.resize(nt * ny * nplist, 0.0); sz.resize(nmaxevent * nz * nplist, 0.0); srz.resize(nmaxevent * nz * nplist, 0.0); - + FIM.resize(nplist * nplist, 0.0); sres.resize(nt * nytrue * nplist, 0.0); } - + ssigmay.resize(nt * ny * nplist, 0.0); ssigmaz.resize(nmaxevent * nz * nplist, 0.0); if (sensi >= SensitivityOrder::second) { @@ -114,23 +104,19 @@ ReturnData::ReturnData( s2rz.resize(nmaxevent * nztrue * nplist * nplist, 0.0); } } - + } void ReturnData::invalidate(const realtype t) { - /** - * @brief routine to set likelihood, state variables, outputs and respective sensitivities to NaN - * (typically after integration failure) - * @param t time of integration failure - */ invalidateLLH(); - + invalidateSLLH(); + // find it corresponding to datapoint after integration failure int it_start; for (it_start = 0; it_start < nt; it_start++) if(ts.at(it_start)>t) break; - + for (int it = it_start; it < nt; it++){ for (int ix = 0; ix < nx; ix++) x.at(ix + nx * it) = getNaN(); @@ -155,26 +141,19 @@ void ReturnData::invalidate(const realtype t) { } } } - -void ReturnData::invalidateLLH() { - /** - * @brief routine to set likelihood and respective sensitivities to NaN - * (typically after integration failure) - */ +void ReturnData::invalidateLLH() +{ llh = getNaN(); chi2 = getNaN(); - std::fill(sllh.begin(),sllh.end(),getNaN()); - std::fill(s2llh.begin(),s2llh.end(),getNaN()); } -void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { - /** - * @brief applies the chain rule to account for parameter transformation - * in the sensitivities of simulation results - * @param model Model from which the ReturnData was obtained - */ +void ReturnData::invalidateSLLH() { + std::fill(sllh.begin(), sllh.end(), getNaN()); + std::fill(s2llh.begin(), s2llh.end(), getNaN()); +} +void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { // chain-rule factor: multiplier for am_p std::vector coefficient(nplist, 1.0); std::vector pcoefficient(nplist, 1.0); @@ -183,7 +162,7 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { unscaleParameters(unscaledParameters, model->getParameterScale(), unscaledParameters); std::vector augcoefficient(np, 1.0); - + if (sensi == SensitivityOrder::second && o2mode == SecondOrderMode::full) { for (int ip = 0; ip < np; ++ip) { switch (pscale[ip]) { @@ -240,12 +219,12 @@ void ReturnData::applyChainRuleFactorToSimulationResults(const Model *model) { for (int ip = 0; ip < nplist; ++ip) sllh.at(ip) *= pcoefficient.at(ip); - + if(!sres.empty()) for (int iyt = 0; iyt < nytrue*nt; ++iyt) for (int ip = 0; ip < nplist; ++ip) sres.at((iyt * nplist + ip)) *= pcoefficient.at(ip); - + if(!FIM.empty()) for (int ip = 0; ip < nplist; ++ip) for (int jp = 0; jp < nplist; ++jp) diff --git a/tests/cpputest/steadystate/tests1.cpp b/tests/cpputest/steadystate/tests1.cpp index 31d7bc2471..8ccbb927ac 100644 --- a/tests/cpputest/steadystate/tests1.cpp +++ b/tests/cpputest/steadystate/tests1.cpp @@ -68,10 +68,10 @@ TEST(groupSteadystate, testCloneModel) { TEST(groupSteadystate, testExpDataFromReturnData) { auto model = getModel(); auto solver = model->getSolver(); - + amici::hdf5::readModelDataFromHDF5(NEW_OPTION_FILE, *model, "/model_steadystate/nosensi/options"); amici::hdf5::readSolverSettingsFromHDF5(NEW_OPTION_FILE, *solver, "/model_steadystate/nosensi/options"); - + auto rdata = runAmiciSimulation(*solver, nullptr, *model); auto edata = amici::ExpData(*rdata, 0.1, 0.1); runAmiciSimulation(*solver, &edata, *model); @@ -88,6 +88,24 @@ TEST(groupSteadystate, testReuseSolver) { runAmiciSimulation(*solver, nullptr, *model); } +TEST(groupSteadystate, testRethrow) { + auto model = getModel(); + auto solver = model->getSolver(); + + amici::hdf5::readModelDataFromHDF5(NEW_OPTION_FILE, *model, "/model_steadystate/nosensi/options"); + amici::hdf5::readSolverSettingsFromHDF5(NEW_OPTION_FILE, *solver, "/model_steadystate/nosensi/options"); + + // p = NaN will raise amici::IntegrationFailure + auto p = model->getParameters(); + std::fill(p.begin(), p.end(), std::nan("")); + model->setParameters(p); + + // must not throw + runAmiciSimulation(*solver, nullptr, *model); + + // must throw + CHECK_THROWS(amici::IntegrationFailure, runAmiciSimulation(*solver, nullptr, *model, true)); +} TEST(groupSteadystate, testSimulation) { amici::simulateVerifyWrite("/model_steadystate/nosensi/"); From fed09f78c7b346b25387d928cb4749277cc7e5de Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Sat, 2 Feb 2019 12:00:15 +0100 Subject: [PATCH 5/8] fix(python) Fix doc and rename sys_pipes to something more meaningful --- python/amici/__init__.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/python/amici/__init__.py b/python/amici/__init__.py index 9f9676af8d..d45153c006 100644 --- a/python/amici/__init__.py +++ b/python/amici/__init__.py @@ -23,6 +23,8 @@ hdf5_enabled: boolean indicating if amici was compiled with hdf5 support has_clibs: boolean indicating if this is the full package with swig interface or the raw package without + capture_cstdout: context to redirect C/C++ stdout to python stdout if + python stdout was redirected (doing nothing if not redirected). """ import os @@ -32,10 +34,10 @@ # redirect C/C++ stdout to python stdout if python stdout is redirected, # e.g. in ipython notebook -sys_pipes = suppress +capture_cstdout = suppress if sys.stdout != sys.__stdout__: try: - from wurlitzer import sys_pipes + from wurlitzer import sys_pipes as capture_cstdout except ModuleNotFoundError: pass @@ -119,7 +121,7 @@ def runAmiciSimulation(model, solver, edata=None): if edata and edata.__class__.__name__ == 'ExpDataPtr': edata = edata.get() - with sys_pipes(): + with capture_cstdout(): rdata = amici.runAmiciSimulation(solver.get(), edata, model.get()) return rdataToNumPyArrays(rdata) @@ -164,7 +166,7 @@ def runAmiciSimulations(model, solver, edata_list, num_threads=1): Raises: """ - with sys_pipes(): + with capture_cstdout(): edata_ptr_vector = amici.ExpDataPtrVector(edata_list) rdata_ptr_list = amici.runAmiciSimulations(solver.get(), edata_ptr_vector, From c99c8a5dc7de249d6588dcb77976f943fec66a3f Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Sun, 3 Feb 2019 19:45:36 +0100 Subject: [PATCH 6/8] feature(python) Check for matching AMICI versions when importing model (Closes #556). Set exact AMICI version as model package requirement. --- python/amici/__init__.template.py | 13 +++++++++++++ python/amici/ode_export.py | 3 ++- python/amici/setup.template.py | 5 ++--- 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/python/amici/__init__.template.py b/python/amici/__init__.template.py index 40a385209f..b68a5bae90 100644 --- a/python/amici/__init__.template.py +++ b/python/amici/__init__.template.py @@ -1,3 +1,16 @@ """AMICI-generated module for model TPL_MODELNAME""" +import amici + +# Ensure we are binary-compatible, see #556 +if 'TPL_AMICI_VERSION' != amici.__version__: + raise RuntimeError('Cannot use model TPL_MODELNAME, generated with AMICI ' + 'version TPL_AMICI_VERSION, together with AMICI version' + f' {amici.__version__} which is present in your ' + 'PYTHONPATH. Install the AMICI package matching the ' + 'model version or regenerate the model with the AMICI ' + 'currently in your path.') + from TPL_MODELNAME.TPL_MODELNAME import * + +__version__ = 'TPL_PACKAGE_VERSION' diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index 3318ef272b..b840f8346d 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -2350,7 +2350,8 @@ def _writeModuleSetup(self): """ templateData = {'MODELNAME': self.modelName, - 'VERSION': '0.1.0'} + 'AMICI_VERSION': __version__, + 'PACKAGE_VERSION': '0.1.0'} applyTemplate(os.path.join(amiciModulePath, 'setup.template.py'), os.path.join(self.modelPath, 'setup.py'), templateData) diff --git a/python/amici/setup.template.py b/python/amici/setup.template.py index 5ddf14d00b..94c5e4d9ec 100644 --- a/python/amici/setup.template.py +++ b/python/amici/setup.template.py @@ -88,7 +88,7 @@ def getAmiciLibs(): # Install setup( name='TPL_MODELNAME', - version='TPL_VERSION', + version='TPL_PACKAGE_VERSION', description='AMICI-generated module for model TPL_MODELNAME', url='https://github.com/ICB-DCM/AMICI', author='model-author-todo', @@ -96,8 +96,7 @@ def getAmiciLibs(): #license = 'BSD', ext_modules=[model_module], packages=find_packages(), - # TODO: should specify amici version with which the model was generated - install_requires=['amici'], + install_requires=['amici==TPL_AMICI_VERSION'], extras_require={'wurlitzer': ['wurlitzer']}, python_requires='>=3.6', package_data={ From 0e44dac6dd4f129f434baa80c6348cf14fc13749 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Wed, 6 Feb 2019 18:10:37 -0500 Subject: [PATCH 7/8] fix(python) fix symbolic computations for adjoint (#583) --- python/amici/ode_export.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index b840f8346d..969a40527b 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -866,7 +866,6 @@ def __init__(self): 'xBdot': { 'x': 'JB', 'y': 'xB', - 'sign': -1, }, } @@ -1430,7 +1429,7 @@ def _compute_equation(self, name): sp.zeros(1, self._eqs[name].shape[1]) elif name == 'JB': - self._eqs[name] = self.eq('J').transpose() + self._eqs[name] = -self.eq('J').transpose() elif name == 'JDiag': self._eqs[name] = getSymbolicDiagonal(self.eq('J')) From ba5cab59b93f50d19975bab8f6ff292404838faa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Wed, 6 Feb 2019 17:07:35 -0500 Subject: [PATCH 8/8] version bump --- version.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/version.txt b/version.txt index 2003b639c4..965065db5b 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.9.2 +0.9.3