Merge pull request #585 from ICB-DCM/develop

Release 0.9.3

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<ReturnData> runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model);
+std::unique_ptr<ReturnData> runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model, bool rethrow=false);
 
 /*!
  * runAmiciSimulations does the same as runAmiciSimulation, but for multiple ExpData instances.

include/amici/edata.h

@@ -38,7 +38,7 @@ class ExpData {
      */
     ExpData(int nytrue, int nztrue, int nmaxevent,
             std::vector<realtype>  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<realtype> 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<realtype> originalFixedParameters;
+    std::vector<realtype> originalTimepoints;
+};
+
 } // namespace amici
 
 #endif /* AMICI_EDATA_H */

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<ParameterScaling> 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<realtype> ssigmay;
-    
+
     /** observable (dimension: nt*ny, row-major) */
     std::vector<realtype> res;
 
     /** parameter derivative of residual (dimension: nt*ny x nplist,
      * row-major) */
     std::vector<realtype> sres;
-    
+
     /** fisher information matrix (dimension: nplist x nplist,
      * row-major) */
     std::vector<realtype> 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<int> 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<realtype> x0;
-    
+
     /** preequilibration steady state found by Newton solver (dimension: nx) */
     std::vector<realtype> x_ss;
 
     /** initial sensitivities (dimension: nplist x nx, row-major) */
     std::vector<realtype> sx0;
-    
+
     /** preequilibration sensitivities found by Newton solver (dimension: nplist x nx, row-major) */
     std::vector<realtype> sx_ss;
 

python/amici/__init__.py

@@ -23,10 +23,23 @@
     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
 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
+capture_cstdout = suppress
+if sys.stdout != sys.__stdout__:
+    try:
+        from wurlitzer import sys_pipes as capture_cstdout
+    except ModuleNotFoundError:
+        pass
 
 hdf5_enabled = False
 has_clibs = False
@@ -103,12 +116,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 capture_cstdout():
+        rdata = amici.runAmiciSimulation(solver.get(), edata, model.get())
     return rdataToNumPyArrays(rdata)
 
 
@@ -134,7 +148,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 +166,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 capture_cstdout():
+        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]

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'

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'))
@@ -2350,7 +2349,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)
 

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'],
@@ -88,17 +88,17 @@ 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',
     author_email='model-author-todo',
     #license = 'BSD',
     ext_modules=[model_module],
     packages=find_packages(),
-    # TODO: should specify amici version with which the model was generated
-    install_requires=['amici'],
-    python_requires='>=3',
+    install_requires=['amici==TPL_AMICI_VERSION'],
+    extras_require={'wurlitzer': ['wurlitzer']},
+    python_requires='>=3.6',
     package_data={
     },
     zip_safe = False,

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*',

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)