From bbc729c5949842d109348dafdee4c893b79e0146 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@helmholtz-muenchen.de>
Date: Thu, 17 Oct 2019 00:42:35 +0200
Subject: [PATCH 01/19] Update documentation and FAQ for CBLAS requirement

---
 INSTALL.md           |  4 ++++
 documentation/FAQ.md | 14 ++++++++++++++
 2 files changed, 18 insertions(+)

diff --git a/INSTALL.md b/INSTALL.md
index e8ddcee092..ed49840522 100755
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -308,6 +308,10 @@ implementations such as Accelerate, Intel MKL, cblas, openblas, atlas.
 The matlab interface uses the MATLAB MKL, which requires no separate
 installation.
 
+On Ubuntu, this requirement can be satisfied with
+
+    apt install libatlas-base-dev
+
 #### C++ compiler
 
 All AMICI installations require a C++11-compatible C++ compiler.
diff --git a/documentation/FAQ.md b/documentation/FAQ.md
index cbdaf30353..ea79fe25cf 100644
--- a/documentation/FAQ.md
+++ b/documentation/FAQ.md
@@ -42,3 +42,17 @@ __Q__: The simulation/sensitivities I get are incorrect.
 
 __A__: There are some known issues, especially with adjoint sensitivities, events and DAEs. If your particular problem is not featured in the [issues](https://github.com/ICB-DCM/AMICI/issues) list, please add it!
 
+---
+
+__Q__: I am trying to install the AMICI Python package, but installation fails
+with something like
+
+    amici/src/cblas.cpp:16:13: fatal error: cblas.h: No such file or directory
+    #include <cblas.h>
+             ^~~~~~~~~
+    compilation terminated.
+    error: command 'x86_64-linux-gnu-gcc' failed with exit status 1
+
+__A__: You will have to install a CBLAS-compatible BLAS library and/or set
+`BLAS_CFLAGS` as described in the [installation guide](INSTALL.md).
+

From 09e87d24b0157ff3ab5e9be26ef77d931af475b8 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Tue, 22 Oct 2019 17:12:46 +0200
Subject: [PATCH 02/19] Simplify symbolic sensitivity expressions during Python
 SBML import (#838)

* Simplify symbolic sensitivity expressions during Python SBML import

Fixes Infs in the Jacobian when using Hill-functions with states of 0.0.

Mark all sympy symbols as reals, which sometimes allows for additional simplifications.

* Simplify all equations, but only if asked to do so

* doxygen cannot handle attribute typehints
---
 documentation/development.md |  6 ++++--
 python/amici/ode_export.py   | 34 +++++++++++++++++++++---------
 python/amici/sbml_import.py  | 40 ++++++++++++++++++++----------------
 3 files changed, 50 insertions(+), 30 deletions(-)

diff --git a/documentation/development.md b/documentation/development.md
index f8b951bc91..b4463bdc85 100644
--- a/documentation/development.md
+++ b/documentation/development.md
@@ -98,8 +98,10 @@ described below:
   case for all existing code, any new contributions should do so. 
 
 * We use Python [type hints](https://docs.python.org/3/library/typing.html)
-  for all functions. In Python code type hints should be used instead of
-  doxygen `@type`. (All legacy `@type` attributes are to be removed.)
+  for all functions (but not for class attributes, since they are not supported
+  by the current Python doxygen filter). In Python code type hints should be
+  used instead of doxygen `@type`. (All legacy `@type` attributes are to be
+  removed.)
     
   For function docstrings, follow this format:
   
diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py
index cfcb1fd97d..8d63c5ad68 100644
--- a/python/amici/ode_export.py
+++ b/python/amici/ode_export.py
@@ -20,6 +20,8 @@
 
 from string import Template
 import sympy.printing.ccode as ccode
+from sympy.matrices.immutable import ImmutableDenseMatrix
+from sympy.matrices.dense import MutableDenseMatrix
 
 from . import (
     amiciSwigPath, amiciSrcPath, amiciModulePath, __version__, __commit__
@@ -712,12 +714,15 @@ class ODEModel:
         derivative from a partial derivative call to enforce a partial
         derivative in the next recursion. prevents infinite recursion
 
+        _simplify: if set to true, symbolic derivative expressions will be
+        simplified. NOTE: This does currently not work with PySB symbols.
     """
 
-    def __init__(self):
+    def __init__(self, simplify: bool = False):
         """Create a new ODEModel instance.
 
         Arguments:
+            simplify: see ODEModel._simplify
 
         Raises:
 
@@ -804,6 +809,7 @@ def __init__(self):
         }
 
         self._lock_total_derivative = False
+        self._simplify = simplify
 
     def import_from_sbml_importer(self, si):
         """Imports a model specification from a amici.SBMLImporter instance.
@@ -825,7 +831,8 @@ def import_from_sbml_importer(self, si):
         self._eqs['dxdotdw'] = si.stoichiometricMatrix
         self._eqs['w'] = si.fluxVector
         self._syms['w'] = sp.Matrix(
-            [sp.Symbol(f'flux_r{idx}') for idx in range(len(si.fluxVector))]
+            [sp.Symbol(f'flux_r{idx}', real=True)
+             for idx in range(len(si.fluxVector))]
         )
         self._eqs['dxdotdx'] = sp.zeros(si.stoichiometricMatrix.shape[0])
         if len(si.stoichiometricMatrix):
@@ -1160,12 +1167,12 @@ def _generateSymbol(self, name):
                 for comp in getattr(self, component)
             ])
             self._strippedsyms[name] = sp.Matrix([
-                sp.Symbol(comp.get_name())
+                sp.Symbol(comp.get_name(), real=True)
                 for comp in getattr(self, component)
             ])
             if name == 'y':
                 self._syms['my'] = sp.Matrix([
-                    sp.Symbol(f'm{strip_pysb(comp.get_id())}')
+                    sp.Symbol(f'm{strip_pysb(comp.get_id())}', real=True)
                     for comp in getattr(self, component)
                 ])
             return
@@ -1178,7 +1185,7 @@ def _generateSymbol(self, name):
             return
         elif name == 'dtcldp':
             self._syms[name] = sp.Matrix([
-                sp.Symbol(f's{strip_pysb(tcl.get_id())}')
+                sp.Symbol(f's{strip_pysb(tcl.get_id())}', real=True)
                 for tcl in self._conservationlaws
             ])
             return
@@ -1193,7 +1200,7 @@ def _generateSymbol(self, name):
             length = len(self.eq(name))
 
         self._syms[name] = sp.Matrix([
-            sp.Symbol(f'{name}{i}') for i in range(length)
+            sp.Symbol(f'{name}{i}', real=True) for i in range(length)
         ])
 
     def generateBasicVariables(self):
@@ -1382,6 +1389,10 @@ def _compute_equation(self, name):
 
             for index, formula in enumerate(self.eq('x0_fixedParameters')):
                 if formula == 0 or formula == 0.0:
+                    # sp.simplify returns ImmutableDenseMatrix, if we need to
+                    # change them, they need to be made mutable
+                    if isinstance(self._eqs[name], ImmutableDenseMatrix):
+                        self._eqs[name] = MutableDenseMatrix(self._eqs[name])
                     self._eqs[name][index, :] = \
                         sp.zeros(1, self._eqs[name].shape[1])
 
@@ -1429,6 +1440,9 @@ def _compute_equation(self, name):
             if not self._lock_total_derivative:
                 self._eqs[name] = self._eqs[name].transpose()
 
+        if self._simplify:
+            self._eqs[name] = sp.simplify(self._eqs[name])
+
     def symNames(self):
         """Returns a list of names of generated symbolic variables
 
@@ -1556,7 +1570,7 @@ def _total_derivative(self, name, eq, chainvars, var,
             if dydx_name is None:
                 dydx_name = f'd{eq}d{chainvar}'
             if dxdz_name is None:
-                dxdz_name =  f'd{chainvar}d{var}'
+                dxdz_name = f'd{chainvar}d{var}'
 
             dydx = self.sym_or_eq(name, dydx_name)
             dxdz = self.sym_or_eq(name, dxdz_name)
@@ -1799,7 +1813,7 @@ class ODEExporter:
         modelSwigPath: path to the generated swig files @type str
 
         allow_reinit_fixpar_initcond: indicates whether reinitialization of
-        initial states depending on fixedParmeters is allowed for this model
+        initial states depending on fixedParameters is allowed for this model
         @type bool
     """
 
@@ -2628,7 +2642,7 @@ def strip_pysb(symbol):
     # this ensures that the pysb type specific __repr__ is used when converting
     # to string
     if pysb and isinstance(symbol, pysb.Component):
-        return sp.Symbol(symbol.name)
+        return sp.Symbol(symbol.name, real=True)
     else:
         # in this case we will use sympy specific transform anyways
         return symbol
@@ -2834,7 +2848,7 @@ def csc_matrix(matrix, name, base_index=0):
         for row in range(0, matrix.rows):
             if not (matrix[row, col] == 0):
                 symbolName = f'{name}{symbol_name_idx}'
-                sparseMatrix[row, col] = sp.Symbol(symbolName)
+                sparseMatrix[row, col] = sp.Symbol(symbolName, real=True)
                 symbolList.append(symbolName)
                 sparseList.append(matrix[row, col])
                 symbolRowVals.append(row)
diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index 2e2a74f926..9d9ad024a9 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -244,7 +244,7 @@ def sbml2amici(self,
         self.reset_symbols()
         self.processSBML(constantParameters)
         self.processObservables(observables, sigmas, noise_distributions)
-        ode_model = ODEModel()
+        ode_model = ODEModel(simplify=True)
         ode_model.import_from_sbml_importer(self)
         exporter = ODEExporter(
             ode_model,
@@ -342,10 +342,10 @@ def _gather_locals(self):
 
         """
         for s in self.sbml.getListOfSpecies():
-            self.local_symbols[s.getId()] = sp.Symbol(s.getId())
+            self.local_symbols[s.getId()] = sp.Symbol(s.getId(), real=True)
 
         for p in self.sbml.getListOfParameters():
-            self.local_symbols[p.getId()] = sp.Symbol(p.getId())
+            self.local_symbols[p.getId()] = sp.Symbol(p.getId(), real=True)
 
     def processSpecies(self):
         """Get species information from SBML model.
@@ -365,12 +365,12 @@ def processSpecies(self):
         }
 
         self.symbols['species']['identifier'] = sp.Matrix(
-            [sp.Symbol(spec.getId()) for spec in species]
+            [sp.Symbol(spec.getId(), real=True) for spec in species]
         )
         self.symbols['species']['name'] = [spec.getName() for spec in species]
 
         self.speciesCompartment = sp.Matrix(
-            [sp.Symbol(spec.getCompartment()) for spec in species]
+            [sp.Symbol(spec.getCompartment(), real=True) for spec in species]
         )
 
         self.constantSpecies = [species_element.getId()
@@ -448,7 +448,8 @@ def get_species_initial(index, conc):
         self.symbols['species']['value'] = species_initial
 
         if self.sbml.isSetConversionFactor():
-            conversion_factor = sp.Symbol(self.sbml.getConversionFactor())
+            conversion_factor = sp.Symbol(self.sbml.getConversionFactor(),
+                                          real=True)
         else:
             conversion_factor = 1.0
 
@@ -516,7 +517,7 @@ def processParameters(self, constantParameters: List[str] = None):
             settings = loop_settings[partype]
 
             self.symbols[partype]['identifier'] = sp.Matrix(
-                [sp.Symbol(par.getId()) for par in settings['var']]
+                [sp.Symbol(par.getId(), real=True) for par in settings['var']]
             )
             self.symbols[partype]['name'] = [
                 par.getName() for par in settings['var']
@@ -547,7 +548,7 @@ def processCompartments(self):
         """
         compartments = self.sbml.getListOfCompartments()
         self.compartmentSymbols = sp.Matrix(
-            [sp.Symbol(comp.getId()) for comp in compartments]
+            [sp.Symbol(comp.getId(), real=True) for comp in compartments]
         )
         self.compartmentVolume = sp.Matrix(
             [sp.sympify(comp.getVolume()) if comp.isSetVolume()
@@ -620,7 +621,7 @@ def getElementStoichiometry(element):
                     if symMath is None:
                         symMath = sp.sympify(element.getStoichiometry())
                 elif element.getId() in rulevars:
-                    return sp.Symbol(element.getId())
+                    return sp.Symbol(element.getId(), real=True)
                 else:
                     # dont put the symbol if it wont get replaced by a
                     # rule
@@ -712,7 +713,8 @@ def processRules(self):
         volumevars = self.compartmentVolume.free_symbols
         compartmentvars = self.compartmentSymbols.free_symbols
         parametervars = sp.Matrix([
-            sp.Symbol(par.getId()) for par in self.sbml.getListOfParameters()
+            sp.Symbol(par.getId(), real=True)
+            for par in self.sbml.getListOfParameters()
         ])
         stoichvars = self.stoichiometricMatrix.free_symbols
 
@@ -810,8 +812,8 @@ def processTime(self):
         Raises:
 
         """
-        sbmlTimeSymbol = sp.Symbol('time')
-        amiciTimeSymbol = sp.Symbol('t')
+        sbmlTimeSymbol = sp.Symbol('time', real=True)
+        amiciTimeSymbol = sp.Symbol('t', real=True)
 
         self.replaceInAllExpressions(sbmlTimeSymbol, amiciTimeSymbol)
 
@@ -903,11 +905,13 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
                 f'x{index}' for index in range(len(speciesSyms))
             ]
             observableSyms = sp.Matrix(
-                [sp.symbols(f'y{index}', real=True) for index in range(len(speciesSyms))]
+                [sp.symbols(f'y{index}', real=True)
+                 for index in range(len(speciesSyms))]
             )
 
         sigmaYSyms = sp.Matrix(
-            [sp.symbols(f'sigma{symbol}', real=True) for symbol in observableSyms]
+            [sp.symbols(f'sigma{symbol}', real=True)
+             for symbol in observableSyms]
         )
         sigmaYValues = sp.Matrix(
             [1.0] * len(observableSyms)
@@ -942,7 +946,7 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
         llhYValues = sp.Matrix(llhYValues)
 
         llhYSyms = sp.Matrix(
-            [sp.Symbol(f'J{symbol}') for symbol in observableSyms]
+            [sp.Symbol(f'J{symbol}', real=True) for symbol in observableSyms]
         )
 
         # set symbols
@@ -1002,8 +1006,8 @@ def cleanReservedSymbols(self):
         """
         reservedSymbols = ['k','p','y','w']
         for str in reservedSymbols:
-            old_symbol = sp.Symbol(str)
-            new_symbol = sp.Symbol('amici_' + str)
+            old_symbol = sp.Symbol(str, real=True)
+            new_symbol = sp.Symbol('amici_' + str, real=True)
             self.replaceInAllExpressions(old_symbol, new_symbol)
             for symbol in self.symbols.keys():
                 if 'identifier' in self.symbols[symbol].keys():
@@ -1022,7 +1026,7 @@ def replaceSpecialConstants(self):
 
         """
         constants = [
-            (sp.Symbol('avogadro'), sp.Symbol('6.02214179e23')),
+            (sp.Symbol('avogadro', real=True), sp.Symbol('6.02214179e23')),
         ]
         for constant, value in constants:
             # do not replace if any symbol is shadowing default definition

From d4a889eee59e5a659727db159473b257fefec0e3 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@helmholtz-muenchen.de>
Date: Tue, 22 Oct 2019 18:55:08 +0200
Subject: [PATCH 03/19] Revert "Simplify symbolic sensitivity expressions
 during Python SBML import (#838)"

This reverts commit 09e87d24b0157ff3ab5e9be26ef77d931af475b8.
---
 documentation/development.md |  6 ++----
 python/amici/ode_export.py   | 34 +++++++++---------------------
 python/amici/sbml_import.py  | 40 ++++++++++++++++--------------------
 3 files changed, 30 insertions(+), 50 deletions(-)

diff --git a/documentation/development.md b/documentation/development.md
index b4463bdc85..f8b951bc91 100644
--- a/documentation/development.md
+++ b/documentation/development.md
@@ -98,10 +98,8 @@ described below:
   case for all existing code, any new contributions should do so. 
 
 * We use Python [type hints](https://docs.python.org/3/library/typing.html)
-  for all functions (but not for class attributes, since they are not supported
-  by the current Python doxygen filter). In Python code type hints should be
-  used instead of doxygen `@type`. (All legacy `@type` attributes are to be
-  removed.)
+  for all functions. In Python code type hints should be used instead of
+  doxygen `@type`. (All legacy `@type` attributes are to be removed.)
     
   For function docstrings, follow this format:
   
diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py
index 8d63c5ad68..cfcb1fd97d 100644
--- a/python/amici/ode_export.py
+++ b/python/amici/ode_export.py
@@ -20,8 +20,6 @@
 
 from string import Template
 import sympy.printing.ccode as ccode
-from sympy.matrices.immutable import ImmutableDenseMatrix
-from sympy.matrices.dense import MutableDenseMatrix
 
 from . import (
     amiciSwigPath, amiciSrcPath, amiciModulePath, __version__, __commit__
@@ -714,15 +712,12 @@ class ODEModel:
         derivative from a partial derivative call to enforce a partial
         derivative in the next recursion. prevents infinite recursion
 
-        _simplify: if set to true, symbolic derivative expressions will be
-        simplified. NOTE: This does currently not work with PySB symbols.
     """
 
-    def __init__(self, simplify: bool = False):
+    def __init__(self):
         """Create a new ODEModel instance.
 
         Arguments:
-            simplify: see ODEModel._simplify
 
         Raises:
 
@@ -809,7 +804,6 @@ def __init__(self, simplify: bool = False):
         }
 
         self._lock_total_derivative = False
-        self._simplify = simplify
 
     def import_from_sbml_importer(self, si):
         """Imports a model specification from a amici.SBMLImporter instance.
@@ -831,8 +825,7 @@ def import_from_sbml_importer(self, si):
         self._eqs['dxdotdw'] = si.stoichiometricMatrix
         self._eqs['w'] = si.fluxVector
         self._syms['w'] = sp.Matrix(
-            [sp.Symbol(f'flux_r{idx}', real=True)
-             for idx in range(len(si.fluxVector))]
+            [sp.Symbol(f'flux_r{idx}') for idx in range(len(si.fluxVector))]
         )
         self._eqs['dxdotdx'] = sp.zeros(si.stoichiometricMatrix.shape[0])
         if len(si.stoichiometricMatrix):
@@ -1167,12 +1160,12 @@ def _generateSymbol(self, name):
                 for comp in getattr(self, component)
             ])
             self._strippedsyms[name] = sp.Matrix([
-                sp.Symbol(comp.get_name(), real=True)
+                sp.Symbol(comp.get_name())
                 for comp in getattr(self, component)
             ])
             if name == 'y':
                 self._syms['my'] = sp.Matrix([
-                    sp.Symbol(f'm{strip_pysb(comp.get_id())}', real=True)
+                    sp.Symbol(f'm{strip_pysb(comp.get_id())}')
                     for comp in getattr(self, component)
                 ])
             return
@@ -1185,7 +1178,7 @@ def _generateSymbol(self, name):
             return
         elif name == 'dtcldp':
             self._syms[name] = sp.Matrix([
-                sp.Symbol(f's{strip_pysb(tcl.get_id())}', real=True)
+                sp.Symbol(f's{strip_pysb(tcl.get_id())}')
                 for tcl in self._conservationlaws
             ])
             return
@@ -1200,7 +1193,7 @@ def _generateSymbol(self, name):
             length = len(self.eq(name))
 
         self._syms[name] = sp.Matrix([
-            sp.Symbol(f'{name}{i}', real=True) for i in range(length)
+            sp.Symbol(f'{name}{i}') for i in range(length)
         ])
 
     def generateBasicVariables(self):
@@ -1389,10 +1382,6 @@ def _compute_equation(self, name):
 
             for index, formula in enumerate(self.eq('x0_fixedParameters')):
                 if formula == 0 or formula == 0.0:
-                    # sp.simplify returns ImmutableDenseMatrix, if we need to
-                    # change them, they need to be made mutable
-                    if isinstance(self._eqs[name], ImmutableDenseMatrix):
-                        self._eqs[name] = MutableDenseMatrix(self._eqs[name])
                     self._eqs[name][index, :] = \
                         sp.zeros(1, self._eqs[name].shape[1])
 
@@ -1440,9 +1429,6 @@ def _compute_equation(self, name):
             if not self._lock_total_derivative:
                 self._eqs[name] = self._eqs[name].transpose()
 
-        if self._simplify:
-            self._eqs[name] = sp.simplify(self._eqs[name])
-
     def symNames(self):
         """Returns a list of names of generated symbolic variables
 
@@ -1570,7 +1556,7 @@ def _total_derivative(self, name, eq, chainvars, var,
             if dydx_name is None:
                 dydx_name = f'd{eq}d{chainvar}'
             if dxdz_name is None:
-                dxdz_name = f'd{chainvar}d{var}'
+                dxdz_name =  f'd{chainvar}d{var}'
 
             dydx = self.sym_or_eq(name, dydx_name)
             dxdz = self.sym_or_eq(name, dxdz_name)
@@ -1813,7 +1799,7 @@ class ODEExporter:
         modelSwigPath: path to the generated swig files @type str
 
         allow_reinit_fixpar_initcond: indicates whether reinitialization of
-        initial states depending on fixedParameters is allowed for this model
+        initial states depending on fixedParmeters is allowed for this model
         @type bool
     """
 
@@ -2642,7 +2628,7 @@ def strip_pysb(symbol):
     # this ensures that the pysb type specific __repr__ is used when converting
     # to string
     if pysb and isinstance(symbol, pysb.Component):
-        return sp.Symbol(symbol.name, real=True)
+        return sp.Symbol(symbol.name)
     else:
         # in this case we will use sympy specific transform anyways
         return symbol
@@ -2848,7 +2834,7 @@ def csc_matrix(matrix, name, base_index=0):
         for row in range(0, matrix.rows):
             if not (matrix[row, col] == 0):
                 symbolName = f'{name}{symbol_name_idx}'
-                sparseMatrix[row, col] = sp.Symbol(symbolName, real=True)
+                sparseMatrix[row, col] = sp.Symbol(symbolName)
                 symbolList.append(symbolName)
                 sparseList.append(matrix[row, col])
                 symbolRowVals.append(row)
diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index 9d9ad024a9..2e2a74f926 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -244,7 +244,7 @@ def sbml2amici(self,
         self.reset_symbols()
         self.processSBML(constantParameters)
         self.processObservables(observables, sigmas, noise_distributions)
-        ode_model = ODEModel(simplify=True)
+        ode_model = ODEModel()
         ode_model.import_from_sbml_importer(self)
         exporter = ODEExporter(
             ode_model,
@@ -342,10 +342,10 @@ def _gather_locals(self):
 
         """
         for s in self.sbml.getListOfSpecies():
-            self.local_symbols[s.getId()] = sp.Symbol(s.getId(), real=True)
+            self.local_symbols[s.getId()] = sp.Symbol(s.getId())
 
         for p in self.sbml.getListOfParameters():
-            self.local_symbols[p.getId()] = sp.Symbol(p.getId(), real=True)
+            self.local_symbols[p.getId()] = sp.Symbol(p.getId())
 
     def processSpecies(self):
         """Get species information from SBML model.
@@ -365,12 +365,12 @@ def processSpecies(self):
         }
 
         self.symbols['species']['identifier'] = sp.Matrix(
-            [sp.Symbol(spec.getId(), real=True) for spec in species]
+            [sp.Symbol(spec.getId()) for spec in species]
         )
         self.symbols['species']['name'] = [spec.getName() for spec in species]
 
         self.speciesCompartment = sp.Matrix(
-            [sp.Symbol(spec.getCompartment(), real=True) for spec in species]
+            [sp.Symbol(spec.getCompartment()) for spec in species]
         )
 
         self.constantSpecies = [species_element.getId()
@@ -448,8 +448,7 @@ def get_species_initial(index, conc):
         self.symbols['species']['value'] = species_initial
 
         if self.sbml.isSetConversionFactor():
-            conversion_factor = sp.Symbol(self.sbml.getConversionFactor(),
-                                          real=True)
+            conversion_factor = sp.Symbol(self.sbml.getConversionFactor())
         else:
             conversion_factor = 1.0
 
@@ -517,7 +516,7 @@ def processParameters(self, constantParameters: List[str] = None):
             settings = loop_settings[partype]
 
             self.symbols[partype]['identifier'] = sp.Matrix(
-                [sp.Symbol(par.getId(), real=True) for par in settings['var']]
+                [sp.Symbol(par.getId()) for par in settings['var']]
             )
             self.symbols[partype]['name'] = [
                 par.getName() for par in settings['var']
@@ -548,7 +547,7 @@ def processCompartments(self):
         """
         compartments = self.sbml.getListOfCompartments()
         self.compartmentSymbols = sp.Matrix(
-            [sp.Symbol(comp.getId(), real=True) for comp in compartments]
+            [sp.Symbol(comp.getId()) for comp in compartments]
         )
         self.compartmentVolume = sp.Matrix(
             [sp.sympify(comp.getVolume()) if comp.isSetVolume()
@@ -621,7 +620,7 @@ def getElementStoichiometry(element):
                     if symMath is None:
                         symMath = sp.sympify(element.getStoichiometry())
                 elif element.getId() in rulevars:
-                    return sp.Symbol(element.getId(), real=True)
+                    return sp.Symbol(element.getId())
                 else:
                     # dont put the symbol if it wont get replaced by a
                     # rule
@@ -713,8 +712,7 @@ def processRules(self):
         volumevars = self.compartmentVolume.free_symbols
         compartmentvars = self.compartmentSymbols.free_symbols
         parametervars = sp.Matrix([
-            sp.Symbol(par.getId(), real=True)
-            for par in self.sbml.getListOfParameters()
+            sp.Symbol(par.getId()) for par in self.sbml.getListOfParameters()
         ])
         stoichvars = self.stoichiometricMatrix.free_symbols
 
@@ -812,8 +810,8 @@ def processTime(self):
         Raises:
 
         """
-        sbmlTimeSymbol = sp.Symbol('time', real=True)
-        amiciTimeSymbol = sp.Symbol('t', real=True)
+        sbmlTimeSymbol = sp.Symbol('time')
+        amiciTimeSymbol = sp.Symbol('t')
 
         self.replaceInAllExpressions(sbmlTimeSymbol, amiciTimeSymbol)
 
@@ -905,13 +903,11 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
                 f'x{index}' for index in range(len(speciesSyms))
             ]
             observableSyms = sp.Matrix(
-                [sp.symbols(f'y{index}', real=True)
-                 for index in range(len(speciesSyms))]
+                [sp.symbols(f'y{index}', real=True) for index in range(len(speciesSyms))]
             )
 
         sigmaYSyms = sp.Matrix(
-            [sp.symbols(f'sigma{symbol}', real=True)
-             for symbol in observableSyms]
+            [sp.symbols(f'sigma{symbol}', real=True) for symbol in observableSyms]
         )
         sigmaYValues = sp.Matrix(
             [1.0] * len(observableSyms)
@@ -946,7 +942,7 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
         llhYValues = sp.Matrix(llhYValues)
 
         llhYSyms = sp.Matrix(
-            [sp.Symbol(f'J{symbol}', real=True) for symbol in observableSyms]
+            [sp.Symbol(f'J{symbol}') for symbol in observableSyms]
         )
 
         # set symbols
@@ -1006,8 +1002,8 @@ def cleanReservedSymbols(self):
         """
         reservedSymbols = ['k','p','y','w']
         for str in reservedSymbols:
-            old_symbol = sp.Symbol(str, real=True)
-            new_symbol = sp.Symbol('amici_' + str, real=True)
+            old_symbol = sp.Symbol(str)
+            new_symbol = sp.Symbol('amici_' + str)
             self.replaceInAllExpressions(old_symbol, new_symbol)
             for symbol in self.symbols.keys():
                 if 'identifier' in self.symbols[symbol].keys():
@@ -1026,7 +1022,7 @@ def replaceSpecialConstants(self):
 
         """
         constants = [
-            (sp.Symbol('avogadro', real=True), sp.Symbol('6.02214179e23')),
+            (sp.Symbol('avogadro'), sp.Symbol('6.02214179e23')),
         ]
         for constant, value in constants:
             # do not replace if any symbol is shadowing default definition

From 5a2f57dc599a8582a446561ee059ba85bbc511bd Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@helmholtz-muenchen.de>
Date: Thu, 24 Oct 2019 16:01:33 +0200
Subject: [PATCH 04/19] Performance: Avoid unnecessary repeated function calls
 for SUNMatrixWrapper dimensions

---
 src/sundials_matrix_wrapper.cpp | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp
index dd107c5efa..efd17ba9da 100644
--- a/src/sundials_matrix_wrapper.cpp
+++ b/src/sundials_matrix_wrapper.cpp
@@ -170,28 +170,31 @@ void SUNMatrixWrapper::multiply(gsl::span<realtype> c, gsl::span<const realtype>
     if (!matrix)
         return;
 
-    if (static_cast<sunindextype>(c.size()) != rows())
+    sunindextype nrows = rows();
+    sunindextype ncols = columns();
+
+    if (static_cast<sunindextype>(c.size()) != nrows)
         throw std::invalid_argument("Dimension mismatch between number of rows "
-                                    "in A (" + std::to_string(rows()) + ") and "
+                                    "in A (" + std::to_string(nrows) + ") and "
                                     "elements in c (" + std::to_string(c.size())
                                     + ")");
 
-    if (static_cast<sunindextype>(b.size()) != columns())
+    if (static_cast<sunindextype>(b.size()) != ncols)
         throw std::invalid_argument("Dimension mismatch between number of cols "
-                                    "in A (" + std::to_string(columns())
+                                    "in A (" + std::to_string(ncols)
                                     + ") and elements in b ("
                                     + std::to_string(b.size()) + ")");
 
     switch (SUNMatGetID(matrix)) {
     case SUNMATRIX_DENSE:
-        amici_dgemv(BLASLayout::colMajor, BLASTranspose::noTrans, rows(),
-                    columns(), 1.0, data(), rows(), b.data(), 1, 1.0, c.data(), 1);
+        amici_dgemv(BLASLayout::colMajor, BLASTranspose::noTrans, nrows,
+                    ncols, 1.0, data(), nrows, b.data(), 1, 1.0, c.data(), 1);
         break;
     case SUNMATRIX_SPARSE:
 
         switch (sparsetype()) {
         case CSC_MAT:
-            for (sunindextype i = 0; i < columns(); ++i) {
+            for (sunindextype i = 0; i < ncols; ++i) {
                 for (sunindextype k = indexptrs_ptr[i]; k < indexptrs_ptr[i + 1];
                      ++k) {
                     c[indexvals_ptr[k]] += data_ptr[k] * b[i];
@@ -199,7 +202,7 @@ void SUNMatrixWrapper::multiply(gsl::span<realtype> c, gsl::span<const realtype>
             }
             break;
         case CSR_MAT:
-            for (sunindextype i = 0; i < rows(); ++i) {
+            for (sunindextype i = 0; i < nrows; ++i) {
                 for (sunindextype k = indexptrs_ptr[i]; k < indexptrs_ptr[i + 1];
                      ++k) {
                     c[i] += data_ptr[k] * b[indexvals_ptr[k]];

From 40580cd22f1f65848593f87d449a4fd8bbc752f4 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Fri, 25 Oct 2019 12:11:12 +0200
Subject: [PATCH 05/19] Use SWIG>=4.0 on travis to include PyDoc in sdist /
 pypi package (#841)

... and print version of detected SWIG

Closes #834
---
 .travis.yml                | 4 ++--
 python/amici/setuptools.py | 2 ++
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index a11d417e50..4754ec533e 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -21,7 +21,6 @@ matrix:
             - libatlas-base-dev
             - lcov
             - libboost-serialization-dev
-            - swig3.0
             - g++-5
             - libc6-dbg
       env:
@@ -95,7 +94,8 @@ matrix:
 
 install:
     - export BASE_DIR=`pwd`
-    - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then mkdir -p ~/bin/ && ln -s /usr/bin/swig3.0 ~/bin/swig && export PATH=~/bin/:$PATH; fi # Python distutils only looks for `swig` and does not find `swig3.0`
+    # Build swig4.0 (not yet available with apt) to include pydoc in source distribution for pypi
+    - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then scripts/downloadAndBuildSwig.sh && export PATH=${BASE_DIR}/ThirdParty/swig-4.0.1/install/bin:${PATH}; fi
     - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then export PYTHON_EXECUTABLE=$(which python3); fi # cmake wont be able to find python3 on its own ...
     - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CI_MODE" == "test" ]]; then pip3 install --upgrade pip==9.0.3 setuptools wheel pkgconfig scipy; fi
     - if [[ "$TRAVIS_OS_NAME" != "linux" ]] && [[ "$CI_MODE" == "test" ]]; then pip3 install --user --upgrade pip==9.0.3 setuptools wheel pkgconfig scipy; fi
diff --git a/python/amici/setuptools.py b/python/amici/setuptools.py
index cb13caff39..0d66522c01 100644
--- a/python/amici/setuptools.py
+++ b/python/amici/setuptools.py
@@ -207,6 +207,8 @@ def generateSwigInterfaceFiles():
     swig_exe = find_swig()
     swig_version = get_swig_version(swig_exe)
 
+    print(f"Found SWIG version {swig_version}")
+
     # Swig AMICI interface without HDF5 dependency
     swig_cmd = [swig_exe,
                 '-c++',

From 862cb6eca383a92404e149fee1e26395f83f84aa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Paszkowski?=
 <45224812+paszkow@users.noreply.github.com>
Date: Thu, 31 Oct 2019 14:10:55 +0100
Subject: [PATCH 06/19] Newton solver extensions (#848)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The changes that allow performing Newton tests from the paper:

G. T. Lines, Ł. Paszkowski, L. Schmiester, D. Weindl, P. Stapor, and J. Hasenauer. Efficient computation of steady states in large-scale ODE models of biochemical reaction networks. accepted for Proceedings of the 8th IFAC Conference on Foundations of Systems Biology in Engineering (FOSBE), Valencia, Spain, October 2019.

Closes #842


Was:
* Remove positivity assurance of found states.

It is no longer required to converge. In its form zero steadystate
solutions were never found.

* Add a lower bound for the damping factor.

Raise AmiException when a damping factor becomes too small.

* Support turning on/off a damping factor.

* Add damping factor mode to be read from hdf5 file.

* Make lower bound of damping factor user adjustable

* Fix typo.

* Revert "Remove positivity assurance of found states."

This reverts commit 7e6986a2c5fea4a1451a2e5bd0e9bd2476970e8c.

* Recheck convergence after enforcing positivity of found states.

* Unify default value for lower bound of a damping factor

* Fix solver's copy constructor.

* Fix solver's comparison operator.

* Fix documentation.

* Add solver serialization.

Serialize damping factor related variables.


Co-authored-by: Łukasz Paszkowski <l.paszkowski@ssc.net.pl>
---
 include/amici/defines.h       |  6 ++++++
 include/amici/newton_solver.h | 10 +++++++++-
 include/amici/serialization.h |  2 ++
 include/amici/solver.h        | 30 ++++++++++++++++++++++++++++++
 src/hdf5.cpp                  | 11 +++++++++++
 src/newton_solver.cpp         |  5 ++++-
 src/solver.cpp                | 16 ++++++++++++++++
 src/steadystateproblem.cpp    | 28 ++++++++++++++++++++--------
 8 files changed, 98 insertions(+), 10 deletions(-)

diff --git a/include/amici/defines.h b/include/amici/defines.h
index 7482578645..fb60acff55 100644
--- a/include/amici/defines.h
+++ b/include/amici/defines.h
@@ -148,6 +148,12 @@ enum class NewtonStatus {
     newt_sim_newt=3,
 };
 
+/** Damping factor flag for the Newton method */
+enum class NewtonDampingFactorMode {
+    off = 0,
+    on = 1
+};
+
 /**
  * @brief msgIdAndTxtFp
  * @param identifier string with error message identifier
diff --git a/include/amici/newton_solver.h b/include/amici/newton_solver.h
index c0c1232b31..5482e0a208 100644
--- a/include/amici/newton_solver.h
+++ b/include/amici/newton_solver.h
@@ -46,6 +46,8 @@ class NewtonSolver {
      * @param maxsteps maximum number of allowed Newton steps for steady state computation
      * @param atol absolute tolerance
      * @param rtol relative tolerance
+     * @param dampingFactorMode to switch on/off a use of a damping factor
+     * @param dampingFactorLowerBound a lower bound for the damping factor
      * @return solver NewtonSolver according to the specified linsolType
      */
     static std::unique_ptr<NewtonSolver> getSolver(realtype *t, AmiVector *x,
@@ -54,7 +56,9 @@ class NewtonSolver {
                                                    ReturnData *rdata,
                                                    int maxlinsteps,
                                                    int maxsteps,
-                                                   double atol, double rtol);
+                                                   double atol, double rtol,
+                                                   NewtonDampingFactorMode dampingFactorMode,
+                                                   double dampingFactorLowerBound);
 
     /**
      * Computes the solution of one Newton iteration
@@ -103,6 +107,10 @@ class NewtonSolver {
     double atol = 1e-16;
     /** relative tolerance */
     double rtol = 1e-8;
+    /** damping factor flag */
+    NewtonDampingFactorMode dampingFactorMode = NewtonDampingFactorMode::on;
+    /** damping factor lower bound */
+    double dampingFactorLowerBound = 1e-8;
 
   protected:
     /** time variable */
diff --git a/include/amici/serialization.h b/include/amici/serialization.h
index 15dee2c22d..6fb4779e8d 100644
--- a/include/amici/serialization.h
+++ b/include/amici/serialization.h
@@ -57,6 +57,8 @@ void serialize(Archive &ar, amici::Solver &u, const unsigned int version) {
     ar &u.requires_preequilibration;
     ar &u.newton_maxsteps;
     ar &u.newton_maxlinsteps;
+    ar &u.newton_damping_factor_mode;
+    ar &u.newton_damping_factor_lower_bound;
     ar &u.ism;
     ar &u.sensi_meth;
     ar &u.linsol;
diff --git a/include/amici/solver.h b/include/amici/solver.h
index 1e9c6dfb41..09855309d7 100644
--- a/include/amici/solver.h
+++ b/include/amici/solver.h
@@ -219,6 +219,30 @@ class Solver {
      */
     void setNewtonMaxLinearSteps(int newton_maxlinsteps);
 
+    /**
+     * @brief Get a state of the damping factor used in the Newton solver
+     * @return
+     */
+    NewtonDampingFactorMode getNewtonDampingFactorMode() const;
+
+    /**
+     * @brief Turn on/off a damping factor in the Newton method
+     * @param dampingFactorMode
+     */
+    void setNewtonDampingFactorMode(NewtonDampingFactorMode dampingFactorMode);
+
+    /**
+     * @brief Get a lower bound of the damping factor used in the Newton solver
+     * @return
+     */
+    double getNewtonDampingFactorLowerBound() const;
+
+    /**
+     * @brief Set a lower bound of the damping factor in the Newton solver
+     * @param dampingFactorLowerBound
+     */
+    void setNewtonDampingFactorLowerBound(double dampingFactorLowerBound);
+
     /**
      * @brief Get sensitvity order
      * @return sensitivity order
@@ -1381,6 +1405,12 @@ class Solver {
      * computation */
     long int newton_maxlinsteps = 0;
 
+    /** Damping factor state used int the Newton method */
+    NewtonDampingFactorMode newton_damping_factor_mode = NewtonDampingFactorMode::on;
+
+    /** Lower bound of the damping factor. */
+    double newton_damping_factor_lower_bound = 1e-8;
+
     /** Enable model preequilibration */
     bool requires_preequilibration = false;
 
diff --git a/src/hdf5.cpp b/src/hdf5.cpp
index 510703a0fe..003b1460d0 100644
--- a/src/hdf5.cpp
+++ b/src/hdf5.cpp
@@ -628,6 +628,17 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver,
                     getIntScalarAttribute(file, datasetPath, "newton_preeq"));
     }
 
+    if(attributeExists(file, datasetPath, "newton_damping_factor_mode")) {
+        solver.setNewtonDampingFactorMode(
+                    static_cast<NewtonDampingFactorMode>(
+                        getIntScalarAttribute(file, datasetPath, "newton_damping_factor_mode")));
+    }
+
+    if(attributeExists(file, datasetPath, "newton_damping_factor_lower_bound")) {
+        solver.setNewtonDampingFactorLowerBound(
+                    getDoubleScalarAttribute(file, datasetPath, "newton_damping_factor_lower_bound"));
+    }
+
     if(attributeExists(file, datasetPath, "newton_maxlinsteps")) {
         solver.setNewtonMaxLinearSteps(
                     getIntScalarAttribute(file, datasetPath,
diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp
index 26d1247baa..ca86694b44 100644
--- a/src/newton_solver.cpp
+++ b/src/newton_solver.cpp
@@ -28,7 +28,8 @@ NewtonSolver::NewtonSolver(realtype *t, AmiVector *x, Model *model,
 
 std::unique_ptr<NewtonSolver> NewtonSolver::getSolver(
         realtype *t, AmiVector *x, LinearSolver linsolType, Model *model,
-        ReturnData *rdata, int maxlinsteps, int maxsteps, double atol, double rtol) {
+        ReturnData *rdata, int maxlinsteps, int maxsteps, double atol, double rtol,
+        NewtonDampingFactorMode dampingFactorMode, double dampingFactorLowerBound) {
 
     std::unique_ptr<NewtonSolver> solver;
 
@@ -76,6 +77,8 @@ std::unique_ptr<NewtonSolver> NewtonSolver::getSolver(
     solver->rtol = rtol;
     solver->maxlinsteps = maxlinsteps;
     solver->maxsteps = maxsteps;
+    solver->dampingFactorMode = dampingFactorMode;
+    solver->dampingFactorLowerBound = dampingFactorLowerBound;
 
     return solver;
 }
diff --git a/src/solver.cpp b/src/solver.cpp
index cad80533a6..2f38c0d0e2 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -20,6 +20,8 @@ Solver::Solver(const Solver &other)
       sensi_meth(other.sensi_meth), stldet(other.stldet),
       ordering(other.ordering), newton_maxsteps(other.newton_maxsteps),
       newton_maxlinsteps(other.newton_maxlinsteps),
+      newton_damping_factor_mode(other.newton_damping_factor_mode),
+      newton_damping_factor_lower_bound(other.newton_damping_factor_lower_bound),
       requires_preequilibration(other.requires_preequilibration), linsol(other.linsol),
       atol(other.atol), rtol(other.rtol), atol_fsa(other.atol_fsa),
       rtol_fsa(other.rtol_fsa), atolB(other.atolB), rtolB(other.rtolB),
@@ -378,6 +380,8 @@ bool operator==(const Solver &a, const Solver &b) {
            (a.ordering == b.ordering) &&
            (a.newton_maxsteps == b.newton_maxsteps) &&
            (a.newton_maxlinsteps == b.newton_maxlinsteps) &&
+           (a.newton_damping_factor_mode == b.newton_damping_factor_mode) &&
+           (a.newton_damping_factor_lower_bound == b.newton_damping_factor_lower_bound) &&
            (a.requires_preequilibration == b.requires_preequilibration) && (a.ism == b.ism) &&
            (a.linsol == b.linsol) && (a.atol == b.atol) && (a.rtol == b.rtol) &&
            (a.maxsteps == b.maxsteps) && (a.maxstepsB == b.maxstepsB) &&
@@ -493,6 +497,18 @@ void Solver::setNewtonMaxLinearSteps(const int newton_maxlinsteps) {
     this->newton_maxlinsteps = newton_maxlinsteps;
 }
 
+NewtonDampingFactorMode Solver::getNewtonDampingFactorMode() const { return newton_damping_factor_mode; }
+
+void Solver::setNewtonDampingFactorMode(NewtonDampingFactorMode dampingFactorMode) {
+  this->newton_damping_factor_mode = dampingFactorMode;
+}
+
+double Solver::getNewtonDampingFactorLowerBound() const { return newton_damping_factor_lower_bound; }
+
+void Solver::setNewtonDampingFactorLowerBound(double dampingFactorLowerBound) {
+  this->newton_damping_factor_lower_bound = dampingFactorLowerBound;
+}
+
 SensitivityOrder Solver::getSensitivityOrder() const { return sensi; }
 
 void Solver::setSensitivityOrder(const SensitivityOrder sensi) {
diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp
index 8ea73c352a..3b3890d77d 100644
--- a/src/steadystateproblem.cpp
+++ b/src/steadystateproblem.cpp
@@ -16,7 +16,7 @@
 #include <cvodes/cvodes.h>
 
 namespace amici {
-    
+
 SteadystateProblem::SteadystateProblem(const Solver *solver,
                                        const AmiVector &x0):
     t(solver->gett()), delta(solver->nx()), ewt(solver->nx()),
@@ -48,7 +48,9 @@ void SteadystateProblem::workSteadyStateProblem(ReturnData *rdata,
     auto newtonSolver = NewtonSolver::getSolver(
         &t, &x, solver->getLinearSolver(), model, rdata,
         solver->getNewtonMaxLinearSteps(), solver->getNewtonMaxSteps(),
-        solver->getAbsoluteTolerance(), solver->getRelativeTolerance());
+        solver->getAbsoluteTolerance(), solver->getRelativeTolerance(),
+        solver->getNewtonDampingFactorMode(),
+        solver->getNewtonDampingFactorLowerBound());
 
     auto newton_status = NewtonStatus::failed;
     try {
@@ -207,20 +209,30 @@ void SteadystateProblem::applyNewtonsMethod(ReturnData *rdata, Model *model,
             converged = wrms < RCONST(1.0);
 
             if (converged) {
-                /* Ensure positivity of the found state */
+                /* Ensure positivity of the found state and recheck if
+                   the convergence still holds */
+                bool recheck_convergence = false;
                 for (ix = 0; ix < model->nx_solver; ix++) {
                     if (x[ix] < 0.0) {
                         x[ix] = 0.0;
-                        converged = FALSE;
+                        recheck_convergence = true;
                     }
                 }
-            } else {
+                if (recheck_convergence) {
+                  model->fxdot(t, x, dx, xdot);
+                  wrms = getWrmsNorm(x_newton, xdot, newtonSolver->atol, newtonSolver->rtol);
+                  converged = wrms < RCONST(1.0);
+                }
+            } else if (newtonSolver->dampingFactorMode==NewtonDampingFactorMode::on) {
                 /* increase dampening factor (superfluous, if converged) */
                 gamma = fmin(1.0, 2.0 * gamma);
             }
-        } else {
-            /* Reduce dampening factor */
+        } else if (newtonSolver->dampingFactorMode==NewtonDampingFactorMode::on) {
+            /* Reduce dampening factor and raise an error when becomes too small */
             gamma = gamma / 4.0;
+            if (gamma < newtonSolver->dampingFactorLowerBound)
+              throw AmiException("Newton solver failed: a damping factor reached its lower bound");
+
             /* No new linear solve, only try new dampening */
             compNewStep = FALSE;
         }
@@ -325,7 +337,7 @@ std::unique_ptr<Solver> SteadystateProblem::createSteadystateSimSolver(
 
     return newton_solver;
 }
-    
+
 void SteadystateProblem::writeSolution(realtype *t, AmiVector &x,
                                        AmiVectorArray &sx) const {
     *t = this->t;

From 2cba50c71c02c66c6a7aa0aca93dccb0c0c243de Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Fri, 1 Nov 2019 07:23:39 +0100
Subject: [PATCH 07/19] Setup performance test on github actions (#853)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Setup performance test on github actions (Closes #846)

* Setup github action to import and run model from Fröhlich et al., Cell Systems, 2018 (forward simulation, forward+adjoint sensitivities)
* Don't require use of ModelPtr.get to call ExpData(Model)
* Fix import in generated model Python package
* Setup AMICI standalone scripts as setuptools entrypoints
* Add reference values for wall times and check we are within the limits
---
 .github/workflows/test-large-model.yml | 75 +++++++++++++++++++++
 python/amici/__init__.py               |  2 +
 python/amici/__init__.template.py      |  2 +-
 python/amici/petab_import.py           | 86 ++++++++++++++++++++++++
 python/bin/amici_import_petab.py       | 92 --------------------------
 python/sdist/setup.py                  |  8 ++-
 tests/performance/check_time.sh        | 32 +++++++++
 tests/performance/reference.yml        |  8 +++
 tests/performance/test.py              | 42 ++++++++++++
 9 files changed, 253 insertions(+), 94 deletions(-)
 create mode 100644 .github/workflows/test-large-model.yml
 delete mode 100755 python/bin/amici_import_petab.py
 create mode 100755 tests/performance/check_time.sh
 create mode 100644 tests/performance/reference.yml
 create mode 100755 tests/performance/test.py

diff --git a/.github/workflows/test-large-model.yml b/.github/workflows/test-large-model.yml
new file mode 100644
index 0000000000..36d51edc63
--- /dev/null
+++ b/.github/workflows/test-large-model.yml
@@ -0,0 +1,75 @@
+name: Performance test
+on:
+  push:
+    branches:
+      - develop
+      - master
+      - feature_846_large_model
+  pull_request:
+    branches:
+      - master
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v1
+      with:
+        fetch-depth: 20
+
+    # install dependencies
+    - name: apt
+      run: |
+        sudo apt-get install -y swig3.0 libatlas-base-dev \
+          && sudo ln -s /usr/bin/swig3.0 /usr/bin/swig
+    - name: pip
+      run: |
+        pip3 install --upgrade --user wheel \
+          && pip3 install --upgrade --user setuptools
+    - run: pip3 install petab shyaml
+    - run: |
+        echo ::add-path::${HOME}/.local/bin/
+        echo ::add-path::${GITHUB_WORKSPACE}/tests/performance/
+    # install AMICI
+    - name: Create AMICI sdist
+      run: |
+        cd python/sdist \
+          && check_time.sh create_sdist /usr/bin/python3 setup.py sdist
+    - name: Install AMICI sdist
+      run: |
+        AMICI_PARALLEL_COMPILE=2 check_time.sh \
+            install_sdist pip3 install -v --user \
+            $(ls -t python/sdist/dist/amici-*.tar.gz | head -1)
+
+    # retrieve test model
+    - run: git clone --depth 1 https://github.com/ICB-DCM/CS_Signalling_ERBB_RAS_AKT
+
+    # import test model
+    - name: Import test model
+      run: |
+        cd CS_Signalling_ERBB_RAS_AKT \
+          && check_time.sh \
+          petab_import amici_import_petab.py -v \
+            -s 'PEtab/CS_Signalling_ERBB_RAS_AKT_petab.xml' \
+            -c 'PEtab/conditions_petab.tsv' \
+            -m 'PEtab/measurements_petab.tsv' \
+            -p 'PEtab/parameters_petab.tsv' --no-compile
+
+    # install model package
+    - name: Install test model
+      run: |
+        check_time.sh install_model \
+          python3 CS_Signalling_ERBB_RAS_AKT/CS_Signalling_ERBB_RAS_AKT_petab/setup.py install --user
+
+    # run simulations
+    - name: forward_simulation
+      run: |
+        check_time.sh forward_simulation tests/performance/test.py forward_simulation
+    - name: forward_sensitivities
+      run: |
+        check_time.sh forward_sensitivities tests/performance/test.py forward_sensitivities
+    - name: adjoint_sensitivities
+      run: |
+        check_time.sh adjoint_sensitivities tests/performance/test.py adjoint_sensitivities
diff --git a/python/amici/__init__.py b/python/amici/__init__.py
index c26c9ff281..4ae68672ef 100644
--- a/python/amici/__init__.py
+++ b/python/amici/__init__.py
@@ -157,6 +157,8 @@ def ExpData(*args):
         # the constructor signature may have changed and you are glad this
         # wrapper did not break.
         return amici.ExpData(args[0].get(), *args[:1])
+    elif isinstance(args[0], ModelPtr):
+        return amici.ExpData(args[0].get())
     else:
         return amici.ExpData(*args)
 
diff --git a/python/amici/__init__.template.py b/python/amici/__init__.template.py
index b68a5bae90..5c42765ec7 100644
--- a/python/amici/__init__.template.py
+++ b/python/amici/__init__.template.py
@@ -11,6 +11,6 @@
                        'model version or regenerate the model with the AMICI '
                        'currently in your path.')
 
-from TPL_MODELNAME.TPL_MODELNAME import *
+from TPL_MODELNAME._TPL_MODELNAME import *
 
 __version__ = 'TPL_PACKAGE_VERSION'
diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py
index 8a0bb8f5ba..7ecbbbec85 100644
--- a/python/amici/petab_import.py
+++ b/python/amici/petab_import.py
@@ -10,6 +10,7 @@
 import logging
 from typing import List, Dict
 import pandas as pd
+import argparse
 
 import petab
 from colorama import Fore
@@ -286,3 +287,88 @@ def show_model_info(sbml_model: 'libsbml.Model'):
                + str(len(sbml_model.getListOfParameters())))
     logger.log(logging.INFO,
                f'Reactions: {len(sbml_model.getListOfReactions())}')
+
+
+def parse_cli_args():
+    """Parse command line arguments"""
+
+    parser = argparse.ArgumentParser(
+        description='Import PEtab-format model into AMICI.')
+
+    # General options:
+    parser.add_argument('-v', '--verbose', dest='verbose', action='store_true',
+                        help='More verbose output')
+    parser.add_argument('-o', '--output-dir', dest='model_output_dir',
+                        help='Name of the model directory to create')
+    parser.add_argument('--no-compile', action='store_false',
+                        dest='compile',
+                        help='Only generate model code, do not compile')
+
+    # Call with set of files
+    parser.add_argument('-s', '--sbml', dest='sbml_file_name',
+                        help='SBML model filename')
+    parser.add_argument('-m', '--measurements', dest='measurement_file_name',
+                        help='Measurement table')
+    parser.add_argument('-c', '--conditions', dest='condition_file_name',
+                        help='Conditions table')
+    parser.add_argument('-p', '--parameters', dest='parameter_file_name',
+                        help='Parameter table')
+
+    # or with model name, following default naming
+    parser.add_argument('-n', '--model-name', dest='model_name',
+                        help='Model name where all files are in the working '
+                        'directory and follow PEtab naming convention. '
+                        'Specifying -[smcp] will override defaults')
+
+    args = parser.parse_args()
+
+    if args.model_name:
+        if not args.sbml_file_name:
+            args.sbml_file_name = petab.get_default_sbml_file_name(
+                args.model_name)
+        if not args.measurement_file_name:
+            args.measurement_file_name = \
+                petab.get_default_measurement_file_name(args.model_name)
+        if not args.condition_file_name:
+            args.condition_file_name = petab.get_default_condition_file_name(
+                args.model_name)
+        if not args.parameter_file_name:
+            args.parameter_file_name = petab.get_default_parameter_file_name(
+                args.model_name)
+
+    if not args.model_name and \
+            (not args.sbml_file_name
+             or not args.condition_file_name
+             or not args.measurement_file_name):
+        parser.error('When not specifying a model name, sbml, '
+                     'condition and measurement file must be specified')
+
+    return args
+
+
+def main():
+    """
+    Command line interface to import a model in the PEtab
+    (https://github.com/ICB-DCM/PEtab/) format into sAMICI
+    """
+    args = parse_cli_args()
+
+    # First check for valid PEtab
+    pp = petab.Problem.from_files(
+        sbml_file=args.sbml_file_name,
+        condition_file=args.condition_file_name,
+        measurement_file=args.measurement_file_name,
+        parameter_file=args.parameter_file_name)
+    petab.lint_problem(pp)
+
+    import_model(model_name=args.model_name,
+                 sbml_file=args.sbml_file_name,
+                 condition_file=args.condition_file_name,
+                 measurement_file=args.measurement_file_name,
+                 model_output_dir=args.model_output_dir,
+                 compile=args.compile,
+                 verbose=True)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/python/bin/amici_import_petab.py b/python/bin/amici_import_petab.py
deleted file mode 100755
index 2aa7174cb5..0000000000
--- a/python/bin/amici_import_petab.py
+++ /dev/null
@@ -1,92 +0,0 @@
-#!/usr/bin/env python3
-
-"""
-Command line interface to import a model in the PEtab
-(https://github.com/ICB-DCM/PEtab/) format into sAMICI
-"""
-
-import petab
-import argparse
-
-from amici.petab_import import import_model
-
-
-def parse_cli_args():
-    """Parse command line arguments"""
-
-    parser = argparse.ArgumentParser(
-        description='Import PEtab-format model into AMICI.')
-
-    # General options:
-    parser.add_argument('-v', '--verbose', dest='verbose', action='store_true',
-                        help='More verbose output')
-    parser.add_argument('-o', '--output-dir', dest='model_output_dir',
-                        help='Name of the model directory to create')
-    parser.add_argument('--no-compile', action='store_false',
-                        dest='compile',
-                        help='Only generate model code, do not compile')
-
-    # Call with set of files
-    parser.add_argument('-s', '--sbml', dest='sbml_file_name',
-                        help='SBML model filename')
-    parser.add_argument('-m', '--measurements', dest='measurement_file_name',
-                        help='Measurement table')
-    parser.add_argument('-c', '--conditions', dest='condition_file_name',
-                        help='Conditions table')
-    parser.add_argument('-p', '--parameters', dest='parameter_file_name',
-                        help='Parameter table')
-
-    # or with model name, following default naming
-    parser.add_argument('-n', '--model-name', dest='model_name',
-                        help='Model name where all files are in the working '
-                        'directory and follow PEtab naming convention. '
-                        'Specifying -[smcp] will override defaults')
-
-    args = parser.parse_args()
-
-    if args.model_name:
-        if not args.sbml_file_name:
-            args.sbml_file_name = petab.get_default_sbml_file_name(
-                args.model_name)
-        if not args.measurement_file_name:
-            args.measurement_file_name = \
-                petab.get_default_measurement_file_name(args.model_name)
-        if not args.condition_file_name:
-            args.condition_file_name = petab.get_default_condition_file_name(
-                args.model_name)
-        if not args.parameter_file_name:
-            args.parameter_file_name = petab.get_default_parameter_file_name(
-                args.model_name)
-
-    if not args.model_name and \
-            (not args.sbml_file_name
-             or not args.condition_file_name
-             or not args.measurement_file_name):
-        parser.error('When not specifying a model name, sbml, '
-                     'condition and measurement file must be specified')
-
-    return args
-
-
-def main():
-    args = parse_cli_args()
-
-    # First check for valid PEtab
-    pp = petab.Problem.from_files(
-        sbml_file=args.sbml_file_name,
-        condition_file=args.condition_file_name,
-        measurement_file=args.measurement_file_name,
-        parameter_file=args.parameter_file_name)
-    petab.lint_problem(pp)
-
-    import_model(model_name=args.model_name,
-                 sbml_file=args.sbml_file_name,
-                 condition_file=args.condition_file_name,
-                 measurement_file=args.measurement_file_name,
-                 model_output_dir=args.model_output_dir,
-                 compile=args.compile,
-                 verbose=True)
-
-
-if __name__ == '__main__':
-    main()
diff --git a/python/sdist/setup.py b/python/sdist/setup.py
index 8373809c3f..20462a32dd 100755
--- a/python/sdist/setup.py
+++ b/python/sdist/setup.py
@@ -177,7 +177,13 @@ def main():
                     ],
         packages=find_packages(),
         package_dir={'amici': 'amici'},
-        scripts=['bin/amici_import_petab.py'],
+        entry_points={
+            'console_scripts': [
+                'amici_import_petab = amici.petab_import:main',
+                # for backwards compatibility
+                'amici_import_petab.py = amici.petab_import:main'
+            ]
+        },
         install_requires=['sympy',
                           'python-libsbml',
                           'h5py',
diff --git a/tests/performance/check_time.sh b/tests/performance/check_time.sh
new file mode 100755
index 0000000000..62b8eb3cf5
--- /dev/null
+++ b/tests/performance/check_time.sh
@@ -0,0 +1,32 @@
+#!/bin/bash
+# Execute command and compare wall time to a reference
+set -e
+
+SCRIPT_PATH=$(dirname "${BASH_SOURCE}")
+SCRIPT_PATH=$(cd "${SCRIPT_PATH}" && pwd)
+
+# YAML file with reference wall times
+REFERENCE_FILE="${SCRIPT_PATH}"/reference.yml
+
+# Reference key
+REF="$1"
+# Command to time
+CMD="${@:2}"
+# Logfile
+LOG=$(tempfile)
+
+# Run and time
+/usr/bin/time -f %e ${CMD} 2>&1 | tee "$LOG"
+RET=${PIPESTATUS[0]}
+test "$RET" != 0 && exit "$RET"
+TIME=$(tail -n1 "$LOG")
+
+# Read reference
+REF_TIME=$(shyaml get-value $REF < $REFERENCE_FILE)
+
+if [[ $TIME > $REF_TIME ]]; then
+  echo "TOO LONG: $TIME > $REF_TIME"
+  exit 1
+else
+  echo "OKAY: $TIME < $REF_TIME"
+fi
diff --git a/tests/performance/reference.yml b/tests/performance/reference.yml
new file mode 100644
index 0000000000..ca97ecf69a
--- /dev/null
+++ b/tests/performance/reference.yml
@@ -0,0 +1,8 @@
+# Reference wall times (seconds) with some buffer
+create_sdist: 25
+install_sdist: 120
+petab_import: 1560 # 26'
+install_model: 400
+forward_simulation: 2
+forward_sensitivities: 6
+adjoint_sensitivities: 120 # TODO reduce when #849 fixed
diff --git a/tests/performance/test.py b/tests/performance/test.py
new file mode 100755
index 0000000000..8ca43c4733
--- /dev/null
+++ b/tests/performance/test.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+
+import amici
+import sys
+
+import CS_Signalling_ERBB_RAS_AKT_petab as model_module
+
+
+def main():
+    arg = sys.argv[1]
+
+    model = model_module.getModel()
+    solver = model.getSolver()
+    # TODO
+    edata = amici.ExpData(model)
+    edata.setTimepoints([1e8])
+    edata.setObservedData([1.0])
+
+    if arg == 'forward_simulation':
+        solver.setSensitivityMethod(amici.SensitivityMethod_none)
+        solver.setSensitivityOrder(amici.SensitivityOrder_none)
+    elif arg == 'forward_sensitivities':
+        solver.setSensitivityMethod(amici.SensitivityMethod_forward)
+        solver.setSensitivityOrder(amici.SensitivityOrder_first)
+    elif arg == 'adjoint_sensitivities':
+        solver.setSensitivityMethod(amici.SensitivityMethod_adjoint)
+        solver.setSensitivityOrder(amici.SensitivityOrder_first)
+    else:
+        print("Unknown argument:", arg)
+        sys.exit(1)
+    rdata = amici.runAmiciSimulation(model, solver, edata)
+
+    diagnostics = ['numsteps', 'numstepsB', 'numrhsevals', 'numrhsevalsB',
+                   'numerrtestfails', 'numerrtestfailsB',
+                   'numnonlinsolvconvfails', 'numnonlinsolvconvfailsB']
+    for d in diagnostics:
+        print(d, rdata[d])
+    assert rdata['status'] == amici.AMICI_SUCCESS
+
+
+if __name__ == '__main__':
+    main()

From 4755ca9ec3fb775ff44d00ceb4b36aef2bbca555 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Sun, 10 Nov 2019 16:07:32 +0100
Subject: [PATCH 08/19] Simplify symbolic sensitivity expressions during Python
 SBML import (#845)

* Simplify symbolic sensitivity expressions during Python SBML import

Fixes Infs in the Jacobian when using Hill-functions with states of 0.0.

Mark all sympy symbols as reals, which sometimes allows for additional simplifications.

* Simplify all equations, but only if asked to do so

* doxygen cannot handle attribute typehints


Note that sympy substitutions only work for if IDs *and* predicates agree
---
 .../workflows/sbml-semantic-test-suite.yml    |  2 +-
 documentation/development.md                  |  6 ++-
 python/amici/ode_export.py                    | 34 +++++++++----
 python/amici/sbml_import.py                   | 51 ++++++++++++-------
 tests/performance/check_time.sh               |  2 +-
 tests/testSBMLSuite.py                        |  4 +-
 6 files changed, 65 insertions(+), 34 deletions(-)

diff --git a/.github/workflows/sbml-semantic-test-suite.yml b/.github/workflows/sbml-semantic-test-suite.yml
index ad7d27e48d..e00bc354ec 100644
--- a/.github/workflows/sbml-semantic-test-suite.yml
+++ b/.github/workflows/sbml-semantic-test-suite.yml
@@ -4,7 +4,7 @@ on:
     branches:
       - develop
       - master
-      - test_action
+      - fix_simplify
   pull_request:
     branches:
       - master
diff --git a/documentation/development.md b/documentation/development.md
index f8b951bc91..b4463bdc85 100644
--- a/documentation/development.md
+++ b/documentation/development.md
@@ -98,8 +98,10 @@ described below:
   case for all existing code, any new contributions should do so. 
 
 * We use Python [type hints](https://docs.python.org/3/library/typing.html)
-  for all functions. In Python code type hints should be used instead of
-  doxygen `@type`. (All legacy `@type` attributes are to be removed.)
+  for all functions (but not for class attributes, since they are not supported
+  by the current Python doxygen filter). In Python code type hints should be
+  used instead of doxygen `@type`. (All legacy `@type` attributes are to be
+  removed.)
     
   For function docstrings, follow this format:
   
diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py
index cfcb1fd97d..8d63c5ad68 100644
--- a/python/amici/ode_export.py
+++ b/python/amici/ode_export.py
@@ -20,6 +20,8 @@
 
 from string import Template
 import sympy.printing.ccode as ccode
+from sympy.matrices.immutable import ImmutableDenseMatrix
+from sympy.matrices.dense import MutableDenseMatrix
 
 from . import (
     amiciSwigPath, amiciSrcPath, amiciModulePath, __version__, __commit__
@@ -712,12 +714,15 @@ class ODEModel:
         derivative from a partial derivative call to enforce a partial
         derivative in the next recursion. prevents infinite recursion
 
+        _simplify: if set to true, symbolic derivative expressions will be
+        simplified. NOTE: This does currently not work with PySB symbols.
     """
 
-    def __init__(self):
+    def __init__(self, simplify: bool = False):
         """Create a new ODEModel instance.
 
         Arguments:
+            simplify: see ODEModel._simplify
 
         Raises:
 
@@ -804,6 +809,7 @@ def __init__(self):
         }
 
         self._lock_total_derivative = False
+        self._simplify = simplify
 
     def import_from_sbml_importer(self, si):
         """Imports a model specification from a amici.SBMLImporter instance.
@@ -825,7 +831,8 @@ def import_from_sbml_importer(self, si):
         self._eqs['dxdotdw'] = si.stoichiometricMatrix
         self._eqs['w'] = si.fluxVector
         self._syms['w'] = sp.Matrix(
-            [sp.Symbol(f'flux_r{idx}') for idx in range(len(si.fluxVector))]
+            [sp.Symbol(f'flux_r{idx}', real=True)
+             for idx in range(len(si.fluxVector))]
         )
         self._eqs['dxdotdx'] = sp.zeros(si.stoichiometricMatrix.shape[0])
         if len(si.stoichiometricMatrix):
@@ -1160,12 +1167,12 @@ def _generateSymbol(self, name):
                 for comp in getattr(self, component)
             ])
             self._strippedsyms[name] = sp.Matrix([
-                sp.Symbol(comp.get_name())
+                sp.Symbol(comp.get_name(), real=True)
                 for comp in getattr(self, component)
             ])
             if name == 'y':
                 self._syms['my'] = sp.Matrix([
-                    sp.Symbol(f'm{strip_pysb(comp.get_id())}')
+                    sp.Symbol(f'm{strip_pysb(comp.get_id())}', real=True)
                     for comp in getattr(self, component)
                 ])
             return
@@ -1178,7 +1185,7 @@ def _generateSymbol(self, name):
             return
         elif name == 'dtcldp':
             self._syms[name] = sp.Matrix([
-                sp.Symbol(f's{strip_pysb(tcl.get_id())}')
+                sp.Symbol(f's{strip_pysb(tcl.get_id())}', real=True)
                 for tcl in self._conservationlaws
             ])
             return
@@ -1193,7 +1200,7 @@ def _generateSymbol(self, name):
             length = len(self.eq(name))
 
         self._syms[name] = sp.Matrix([
-            sp.Symbol(f'{name}{i}') for i in range(length)
+            sp.Symbol(f'{name}{i}', real=True) for i in range(length)
         ])
 
     def generateBasicVariables(self):
@@ -1382,6 +1389,10 @@ def _compute_equation(self, name):
 
             for index, formula in enumerate(self.eq('x0_fixedParameters')):
                 if formula == 0 or formula == 0.0:
+                    # sp.simplify returns ImmutableDenseMatrix, if we need to
+                    # change them, they need to be made mutable
+                    if isinstance(self._eqs[name], ImmutableDenseMatrix):
+                        self._eqs[name] = MutableDenseMatrix(self._eqs[name])
                     self._eqs[name][index, :] = \
                         sp.zeros(1, self._eqs[name].shape[1])
 
@@ -1429,6 +1440,9 @@ def _compute_equation(self, name):
             if not self._lock_total_derivative:
                 self._eqs[name] = self._eqs[name].transpose()
 
+        if self._simplify:
+            self._eqs[name] = sp.simplify(self._eqs[name])
+
     def symNames(self):
         """Returns a list of names of generated symbolic variables
 
@@ -1556,7 +1570,7 @@ def _total_derivative(self, name, eq, chainvars, var,
             if dydx_name is None:
                 dydx_name = f'd{eq}d{chainvar}'
             if dxdz_name is None:
-                dxdz_name =  f'd{chainvar}d{var}'
+                dxdz_name = f'd{chainvar}d{var}'
 
             dydx = self.sym_or_eq(name, dydx_name)
             dxdz = self.sym_or_eq(name, dxdz_name)
@@ -1799,7 +1813,7 @@ class ODEExporter:
         modelSwigPath: path to the generated swig files @type str
 
         allow_reinit_fixpar_initcond: indicates whether reinitialization of
-        initial states depending on fixedParmeters is allowed for this model
+        initial states depending on fixedParameters is allowed for this model
         @type bool
     """
 
@@ -2628,7 +2642,7 @@ def strip_pysb(symbol):
     # this ensures that the pysb type specific __repr__ is used when converting
     # to string
     if pysb and isinstance(symbol, pysb.Component):
-        return sp.Symbol(symbol.name)
+        return sp.Symbol(symbol.name, real=True)
     else:
         # in this case we will use sympy specific transform anyways
         return symbol
@@ -2834,7 +2848,7 @@ def csc_matrix(matrix, name, base_index=0):
         for row in range(0, matrix.rows):
             if not (matrix[row, col] == 0):
                 symbolName = f'{name}{symbol_name_idx}'
-                sparseMatrix[row, col] = sp.Symbol(symbolName)
+                sparseMatrix[row, col] = sp.Symbol(symbolName, real=True)
                 symbolList.append(symbolName)
                 sparseList.append(matrix[row, col])
                 symbolRowVals.append(row)
diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index 2e2a74f926..e563b77059 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -244,7 +244,7 @@ def sbml2amici(self,
         self.reset_symbols()
         self.processSBML(constantParameters)
         self.processObservables(observables, sigmas, noise_distributions)
-        ode_model = ODEModel()
+        ode_model = ODEModel(simplify=True)
         ode_model.import_from_sbml_importer(self)
         exporter = ODEExporter(
             ode_model,
@@ -342,10 +342,21 @@ def _gather_locals(self):
 
         """
         for s in self.sbml.getListOfSpecies():
-            self.local_symbols[s.getId()] = sp.Symbol(s.getId())
+            self.local_symbols[s.getId()] = sp.Symbol(s.getId(), real=True)
 
         for p in self.sbml.getListOfParameters():
-            self.local_symbols[p.getId()] = sp.Symbol(p.getId())
+            self.local_symbols[p.getId()] = sp.Symbol(p.getId(), real=True)
+
+        for c in self.sbml.getListOfCompartments():
+            self.local_symbols[c.getId()] = sp.Symbol(c.getId(), real=True)
+
+        for r in self.sbml.getListOfRules():
+            self.local_symbols[r.getVariable()] = sp.Symbol(r.getVariable(),
+                                                            real=True)
+
+        # SBML time symbol + constants
+        self.local_symbols['time'] = sp.Symbol('time', real=True)
+        self.local_symbols['avogadro'] = sp.Symbol('avogadro', real=True)
 
     def processSpecies(self):
         """Get species information from SBML model.
@@ -365,12 +376,12 @@ def processSpecies(self):
         }
 
         self.symbols['species']['identifier'] = sp.Matrix(
-            [sp.Symbol(spec.getId()) for spec in species]
+            [sp.Symbol(spec.getId(), real=True) for spec in species]
         )
         self.symbols['species']['name'] = [spec.getName() for spec in species]
 
         self.speciesCompartment = sp.Matrix(
-            [sp.Symbol(spec.getCompartment()) for spec in species]
+            [sp.Symbol(spec.getCompartment(), real=True) for spec in species]
         )
 
         self.constantSpecies = [species_element.getId()
@@ -448,7 +459,8 @@ def get_species_initial(index, conc):
         self.symbols['species']['value'] = species_initial
 
         if self.sbml.isSetConversionFactor():
-            conversion_factor = sp.Symbol(self.sbml.getConversionFactor())
+            conversion_factor = sp.Symbol(self.sbml.getConversionFactor(),
+                                          real=True)
         else:
             conversion_factor = 1.0
 
@@ -516,7 +528,7 @@ def processParameters(self, constantParameters: List[str] = None):
             settings = loop_settings[partype]
 
             self.symbols[partype]['identifier'] = sp.Matrix(
-                [sp.Symbol(par.getId()) for par in settings['var']]
+                [sp.Symbol(par.getId(), real=True) for par in settings['var']]
             )
             self.symbols[partype]['name'] = [
                 par.getName() for par in settings['var']
@@ -547,7 +559,7 @@ def processCompartments(self):
         """
         compartments = self.sbml.getListOfCompartments()
         self.compartmentSymbols = sp.Matrix(
-            [sp.Symbol(comp.getId()) for comp in compartments]
+            [sp.Symbol(comp.getId(), real=True) for comp in compartments]
         )
         self.compartmentVolume = sp.Matrix(
             [sp.sympify(comp.getVolume()) if comp.isSetVolume()
@@ -620,7 +632,7 @@ def getElementStoichiometry(element):
                     if symMath is None:
                         symMath = sp.sympify(element.getStoichiometry())
                 elif element.getId() in rulevars:
-                    return sp.Symbol(element.getId())
+                    return sp.Symbol(element.getId(), real=True)
                 else:
                     # dont put the symbol if it wont get replaced by a
                     # rule
@@ -712,7 +724,8 @@ def processRules(self):
         volumevars = self.compartmentVolume.free_symbols
         compartmentvars = self.compartmentSymbols.free_symbols
         parametervars = sp.Matrix([
-            sp.Symbol(par.getId()) for par in self.sbml.getListOfParameters()
+            sp.Symbol(par.getId(), real=True)
+            for par in self.sbml.getListOfParameters()
         ])
         stoichvars = self.stoichiometricMatrix.free_symbols
 
@@ -810,8 +823,8 @@ def processTime(self):
         Raises:
 
         """
-        sbmlTimeSymbol = sp.Symbol('time')
-        amiciTimeSymbol = sp.Symbol('t')
+        sbmlTimeSymbol = sp.Symbol('time', real=True)
+        amiciTimeSymbol = sp.Symbol('t', real=True)
 
         self.replaceInAllExpressions(sbmlTimeSymbol, amiciTimeSymbol)
 
@@ -903,11 +916,13 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
                 f'x{index}' for index in range(len(speciesSyms))
             ]
             observableSyms = sp.Matrix(
-                [sp.symbols(f'y{index}', real=True) for index in range(len(speciesSyms))]
+                [sp.symbols(f'y{index}', real=True)
+                 for index in range(len(speciesSyms))]
             )
 
         sigmaYSyms = sp.Matrix(
-            [sp.symbols(f'sigma{symbol}', real=True) for symbol in observableSyms]
+            [sp.symbols(f'sigma{symbol}', real=True)
+             for symbol in observableSyms]
         )
         sigmaYValues = sp.Matrix(
             [1.0] * len(observableSyms)
@@ -942,7 +957,7 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
         llhYValues = sp.Matrix(llhYValues)
 
         llhYSyms = sp.Matrix(
-            [sp.Symbol(f'J{symbol}') for symbol in observableSyms]
+            [sp.Symbol(f'J{symbol}', real=True) for symbol in observableSyms]
         )
 
         # set symbols
@@ -1002,8 +1017,8 @@ def cleanReservedSymbols(self):
         """
         reservedSymbols = ['k','p','y','w']
         for str in reservedSymbols:
-            old_symbol = sp.Symbol(str)
-            new_symbol = sp.Symbol('amici_' + str)
+            old_symbol = sp.Symbol(str, real=True)
+            new_symbol = sp.Symbol('amici_' + str, real=True)
             self.replaceInAllExpressions(old_symbol, new_symbol)
             for symbol in self.symbols.keys():
                 if 'identifier' in self.symbols[symbol].keys():
@@ -1022,7 +1037,7 @@ def replaceSpecialConstants(self):
 
         """
         constants = [
-            (sp.Symbol('avogadro'), sp.Symbol('6.02214179e23')),
+            (sp.Symbol('avogadro', real=True), sp.Symbol('6.02214179e23')),
         ]
         for constant, value in constants:
             # do not replace if any symbol is shadowing default definition
diff --git a/tests/performance/check_time.sh b/tests/performance/check_time.sh
index 62b8eb3cf5..cc40e96f60 100755
--- a/tests/performance/check_time.sh
+++ b/tests/performance/check_time.sh
@@ -24,7 +24,7 @@ TIME=$(tail -n1 "$LOG")
 # Read reference
 REF_TIME=$(shyaml get-value $REF < $REFERENCE_FILE)
 
-if [[ $TIME > $REF_TIME ]]; then
+if (( TIME > REF_TIME )); then
   echo "TOO LONG: $TIME > $REF_TIME"
   exit 1
 else
diff --git a/tests/testSBMLSuite.py b/tests/testSBMLSuite.py
index 0b8c07c52b..0a1816c8e9 100755
--- a/tests/testSBMLSuite.py
+++ b/tests/testSBMLSuite.py
@@ -118,7 +118,7 @@ def concentrations_to_amounts(amount_species, wrapper, model, simulated_x):
                 )
             })
             volume = volume.subs({
-                sp.Symbol(name): value
+                sp.Symbol(name, real=True): value
                 for name, value in zip(
                     model.getParameterIds(),
                     model.getParameters()
@@ -127,7 +127,7 @@ def concentrations_to_amounts(amount_species, wrapper, model, simulated_x):
 
             # required for 525-527, 530 as k is renamed to amici_k
             volume = volume.subs({
-                sp.Symbol(name): value
+                sp.Symbol(name, real=True): value
                 for name, value in zip(
                     model.getParameterNames(),
                     model.getParameters()

From 5dfc4abddf5cb50611e9881fcab074d8859ccc6b Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Tue, 19 Nov 2019 23:11:20 +0100
Subject: [PATCH 09/19] Fix time comparison for performance test (#864)

* Fix bash float comparison

* Raise import time limit
---
 tests/performance/check_time.sh | 2 +-
 tests/performance/reference.yml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/performance/check_time.sh b/tests/performance/check_time.sh
index cc40e96f60..d8f440e1f2 100755
--- a/tests/performance/check_time.sh
+++ b/tests/performance/check_time.sh
@@ -24,7 +24,7 @@ TIME=$(tail -n1 "$LOG")
 # Read reference
 REF_TIME=$(shyaml get-value $REF < $REFERENCE_FILE)
 
-if (( TIME > REF_TIME )); then
+if (( $(echo "$TIME > $REF_TIME" | bc) )); then
   echo "TOO LONG: $TIME > $REF_TIME"
   exit 1
 else
diff --git a/tests/performance/reference.yml b/tests/performance/reference.yml
index ca97ecf69a..afc02704e7 100644
--- a/tests/performance/reference.yml
+++ b/tests/performance/reference.yml
@@ -1,7 +1,7 @@
 # Reference wall times (seconds) with some buffer
 create_sdist: 25
 install_sdist: 120
-petab_import: 1560 # 26'
+petab_import: 4500 # TODO reduce to 1560 after #862  ^= 26'
 install_model: 400
 forward_simulation: 2
 forward_sensitivities: 6

From d4a50b3c33709ae948dd15c2dd4e52e527248a56 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Wed, 20 Nov 2019 09:12:01 +0100
Subject: [PATCH 10/19] =?UTF-8?q?Add=20application=20context=20to=20handle?=
 =?UTF-8?q?=20custom=20output=20functions=20and=20other=20'=E2=80=A6=20(#8?=
 =?UTF-8?q?40)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* Add application context to handle custom output functions for other 'global' settings

Add AmiciApplication class as context for handling so far global settings.
This allows for example setting custom logging functions for concurrent AMICI
runs, e.g. in multi-thread applications (Closes #576).

Use AmiciApplication transparently when calling amici::runAmiciSimulation(s)

* Update swig ignore

* Update Matlab code and cleanup

* Regenerate Matlab example models
---
 include/amici/amici.h                         | 101 +++++++--
 include/amici/defines.h                       |   9 +-
 include/amici/misc.h                          |  17 +-
 include/amici/model.h                         |   6 +
 include/amici/model_dae.h                     |   1 -
 include/amici/model_ode.h                     |   1 -
 include/amici/solver.h                        |  17 +-
 include/amici/solver_cvodes.h                 |   2 -
 include/amici/solver_idas.h                   |   1 -
 models/model_dirac/CMakeLists.txt             |   7 +
 models/model_dirac/main.cpp                   |   3 +-
 models/model_dirac/model_dirac.h              |   4 +-
 models/model_events/CMakeLists.txt            |   7 +
 models/model_events/main.cpp                  |   3 +-
 models/model_events/model_events.h            |   4 +-
 models/model_jakstat_adjoint/CMakeLists.txt   |   7 +
 models/model_jakstat_adjoint/main.cpp         |   3 +-
 .../model_jakstat_adjoint.h                   |   4 +-
 .../model_jakstat_adjoint_o2/CMakeLists.txt   |   7 +
 models/model_jakstat_adjoint_o2/main.cpp      |   3 +-
 .../model_jakstat_adjoint_o2.h                |   4 +-
 models/model_nested_events/CMakeLists.txt     |   7 +
 models/model_nested_events/main.cpp           |   3 +-
 .../model_nested_events/model_nested_events.h |   4 +-
 models/model_neuron/CMakeLists.txt            |   7 +
 models/model_neuron/main.cpp                  |   3 +-
 models/model_neuron/model_neuron.h            |   4 +-
 models/model_neuron_o2/CMakeLists.txt         |   7 +
 models/model_neuron_o2/main.cpp               |   3 +-
 models/model_neuron_o2/model_neuron_o2.h      |   4 +-
 models/model_robertson/CMakeLists.txt         |   7 +
 models/model_robertson/main.cpp               |   3 +-
 models/model_robertson/model_robertson.h      |   4 +-
 models/model_steadystate/CMakeLists.txt       |   7 +
 models/model_steadystate/main.cpp             |   3 +-
 models/model_steadystate/model_steadystate.h  |   4 +-
 src/amici.cpp                                 | 199 +++++++++++++-----
 src/forwardproblem.cpp                        |   6 +-
 src/interface_matlab.cpp                      |  26 ++-
 src/misc.cpp                                  |  40 ++--
 src/model.cpp                                 |  56 ++---
 src/solver.cpp                                |  33 +--
 src/solver_cvodes.cpp                         |  16 +-
 src/solver_idas.cpp                           |  13 +-
 swig/amici.i                                  |   4 +-
 45 files changed, 465 insertions(+), 209 deletions(-)

diff --git a/include/amici/amici.h b/include/amici/amici.h
index a0848e692a..89b4bbca90 100644
--- a/include/amici/amici.h
+++ b/include/amici/amici.h
@@ -13,32 +13,105 @@
 namespace amici {
 
 /*!
- * @brief Prints a specified error message associated to the specified
+ * @brief Prints a specified error message associated with the specified
  * identifier
  *
- * @param identifier error identifier
- * @param format string with error message printf-style format
- * @param ... arguments to be formatted
+ * @param id error identifier
+ * @param message error message
  */
 void
-printErrMsgIdAndTxt(const char* identifier, const char* format, ...);
+printErrMsgIdAndTxt(std::string const& id, std::string const& message);
 
 /*!
- * @brief Prints a specified warning message associated to the specified
+ * @brief Prints a specified warning message associated with the specified
  * identifier
  *
- * @param identifier warning identifier
- * @param format string with error message printf-style format
- * @param ... arguments to be formatted
+ * @param id warning identifier
+ * @param message warning message
  */
 void
-printWarnMsgIdAndTxt(const char* identifier, const char* format, ...);
+printWarnMsgIdAndTxt(std::string const& id, std::string const& message);
 
-/** errMsgIdAndTxt is a function pointer for printErrMsgIdAndTxt  */
-extern msgIdAndTxtFp errMsgIdAndTxt;
+/**
+ * @brief Main class for making calls to AMICI.
+ *
+ * This class is used to provide separate AMICI contexts, for example, for use
+ * in multi-threaded applications where different threads want to use AMICI with
+ * different settings, such custom logging functions.
+ *
+ * NOTE: For this moment, the context object needs to be set manually to any
+ * Model and Solver object. If not set, they will use the default output
+ * channel.
+ */
+class AmiciApplication
+{
+  public:
+    AmiciApplication() = default;
+
+    /**
+     * Core integration routine. Initializes the solver and runs the forward
+     * and backward problem.
+     *
+     * @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,
+                                                   bool rethrow = false);
+
+    /**
+     * Same as runAmiciSimulation, but for multiple ExpData instances.
+     *
+     * @param solver Solver instance
+     * @param edatas experimental data objects
+     * @param model model specification object
+     * @param failfast flag to allow early termination
+     * @param num_threads number of threads for parallel execution
+     * @return vector of pointers to return data objects
+     */
+    std::vector<std::unique_ptr<ReturnData>> runAmiciSimulations(
+      Solver const& solver,
+      const std::vector<ExpData*>& edatas,
+      Model const& model,
+      bool failfast,
+      int num_threads);
+
+    /** Function to process warnings */
+    outputFunctionType warning = printWarnMsgIdAndTxt;
+
+    /** Function to process errors */
+    outputFunctionType error = printErrMsgIdAndTxt;
+
+    /**
+     * @brief printf interface to warning()
+     * @param identifier warning identifier
+     * @param format string with warning message printf-style format
+     * @param ... arguments to be formatted
+     */
+    void warningF(const char* identifier, const char* format, ...);
+
+    /**
+     * @brief printf interface to error()
+     * @param identifier warning identifier
+     * @param format string with error message printf-style format
+     * @param ... arguments to be formatted
+     */
+    void errorF(const char* identifier, const char* format, ...);
 
-/** warnMsgIdAndTxt is a function pointer for printWarnMsgIdAndTxt  */
-extern msgIdAndTxtFp warnMsgIdAndTxt;
+    /**
+     * @brief Checks the values in an array for NaNs and Infs
+     *
+     * @param array array
+     * @param fun name of calling function
+     * @return AMICI_RECOVERABLE_ERROR if a NaN/Inf value was found,
+     * AMICI_SUCCESS otherwise
+     */
+    int checkFinite(gsl::span<const realtype> array, const char* fun);
+};
 
 /*!
  * runAmiciSimulation is the core integration routine. It initializes the solver
diff --git a/include/amici/defines.h b/include/amici/defines.h
index fb60acff55..1c90baa92c 100644
--- a/include/amici/defines.h
+++ b/include/amici/defines.h
@@ -2,6 +2,7 @@
 #define AMICI_DEFINES_H
 
 #include <cmath>
+#include <functional>
 
 namespace amici {
 
@@ -155,12 +156,10 @@ enum class NewtonDampingFactorMode {
 };
 
 /**
- * @brief msgIdAndTxtFp
- * @param identifier string with error message identifier
- * @param format string with error message printf-style format
- * @param ... arguments to be formatted
+ * Type for function to process warnings or error messages.
  */
-using msgIdAndTxtFp = void (*)(const char *, const char *, ...);
+using outputFunctionType = std::function<void(std::string const& identifier,
+                                              std::string const& message)>;
 
 // clang-format on
 
diff --git a/include/amici/misc.h b/include/amici/misc.h
index 6277996699..6284a1d8f9 100644
--- a/include/amici/misc.h
+++ b/include/amici/misc.h
@@ -25,14 +25,6 @@ namespace amici {
  gsl::span<realtype> slice(std::vector<realtype> &data, int index,
                            unsigned size);
 
-/**
- * @brief Checks the values in an array for NaNs and Infs
- *
- * @param array array
- * @param fun name of calling function
- * @return AMICI_RECOVERABLE_ERROR if a NaN/Inf value was found, AMICI_SUCCESS otherwise
- */
-int checkFinite(gsl::span<const realtype> array, const char* fun);
 
 
 /**
@@ -92,6 +84,15 @@ std::string backtraceString(int maxFrames);
  */
 std::string regexErrorToString(std::regex_constants::error_type err_type);
 
+/**
+ * @brief Format printf-style arguments to std::string
+ * @param fmt Format string
+ * @param ap Argument list pointer
+ * @return Formatted String
+ */
+std::string printfToString(const char *fmt, va_list ap);
+
+
 } // namespace amici
 
 #ifndef __cpp_lib_make_unique
diff --git a/include/amici/model.h b/include/amici/model.h
index 351c46f9b0..b08e5492cd 100644
--- a/include/amici/model.h
+++ b/include/amici/model.h
@@ -15,6 +15,9 @@ namespace amici {
 class ExpData;
 class Model;
 class Solver;
+class AmiciApplication;
+
+extern AmiciApplication defaultContext;
 
 } // namespace amici
 
@@ -1154,6 +1157,9 @@ class Model : public AbstractModel {
      * nx_solver, row-major) */
     AmiVectorArray dxdotdp;
 
+    /** AMICI context */
+    AmiciApplication *app = &defaultContext;
+
   protected:
     /**
      * @brief Writes part of a slice to a buffer according to indices specified
diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h
index 3218519900..bb8ed38cd2 100644
--- a/include/amici/model_dae.h
+++ b/include/amici/model_dae.h
@@ -13,7 +13,6 @@
 #include <vector>
 
 namespace amici {
-extern msgIdAndTxtFp warnMsgIdAndTxt;
 
 class ExpData;
 class IDASolver;
diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h
index 41aa2be377..5e223af243 100644
--- a/include/amici/model_ode.h
+++ b/include/amici/model_ode.h
@@ -14,7 +14,6 @@
 #include <vector>
 
 namespace amici {
-extern msgIdAndTxtFp warnMsgIdAndTxt;
 
 class CVodeSolver;
 
diff --git a/include/amici/solver.h b/include/amici/solver.h
index 09855309d7..8ab0b75c17 100644
--- a/include/amici/solver.h
+++ b/include/amici/solver.h
@@ -5,6 +5,7 @@
 #include "amici/sundials_linsol_wrapper.h"
 #include "amici/symbolic_functions.h"
 #include "amici/vector.h"
+#include "amici/amici.h"
 
 #include <functional>
 #include <memory>
@@ -16,6 +17,9 @@ class ForwardProblem;
 class BackwardProblem;
 class Model;
 class Solver;
+class AmiciApplication;
+
+extern AmiciApplication defaultContext;
 } // namespace amici
 
 // for serialization friend in Solver
@@ -44,6 +48,12 @@ class Solver {
   public:
     Solver() = default;
 
+    /**
+     * @brief Constructor
+     * @param app AMICI application context
+     */
+    Solver(AmiciApplication *app);
+
     /**
      * @brief Solver copy constructor
      * @param other
@@ -702,6 +712,9 @@ class Solver {
      */
     friend bool operator==(const Solver &a, const Solver &b);
 
+    /** AMICI context */
+    AmiciApplication *app = &defaultContext;
+
   protected:
     /**
      * @brief Sets a timepoint at which the simulation will be stopped
@@ -1490,13 +1503,13 @@ bool operator==(const Solver &a, const Solver &b);
 
 /**
  * @brief Extracts diagnosis information from solver memory block and
- * writes them into the return data object for the backward problem
+ * passes them to the specified output function
  *
  * @param error_code error identifier
  * @param module name of the module in which the error occured
  * @param function name of the function in which the error occured
  * @param msg error message
- * @param eh_data unused input
+ * @param eh_data amici::Solver as void*
  */
 void wrapErrHandlerFn(int error_code, const char *module,
                       const char *function, char *msg, void *eh_data);
diff --git a/include/amici/solver_cvodes.h b/include/amici/solver_cvodes.h
index ce2278c3d8..72e012aca1 100644
--- a/include/amici/solver_cvodes.h
+++ b/include/amici/solver_cvodes.h
@@ -26,8 +26,6 @@ namespace amici {
 
 class CVodeSolver : public Solver {
   public:
-    CVodeSolver() = default;
-
     ~CVodeSolver() override = default;
 
     /**
diff --git a/include/amici/solver_idas.h b/include/amici/solver_idas.h
index 21648ddd88..5a9bf100d8 100644
--- a/include/amici/solver_idas.h
+++ b/include/amici/solver_idas.h
@@ -24,7 +24,6 @@ namespace amici {
 
 class IDASolver : public Solver {
   public:
-    IDASolver() = default;
     ~IDASolver() override = default;
 
     /**
diff --git a/models/model_dirac/CMakeLists.txt b/models/model_dirac/CMakeLists.txt
index 2fbbab79a1..f73a171551 100644
--- a/models/model_dirac/CMakeLists.txt
+++ b/models/model_dirac/CMakeLists.txt
@@ -52,6 +52,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_dirac/main.cpp b/models/model_dirac/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_dirac/main.cpp
+++ b/models/model_dirac/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_dirac/model_dirac.h b/models/model_dirac/model_dirac.h
index 839b28bb78..389dce41c7 100644
--- a/models/model_dirac/model_dirac.h
+++ b/models/model_dirac/model_dirac.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_dirac_h
 #define _amici_model_dirac_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -61,7 +61,7 @@ class Model_model_dirac : public amici::Model_ODE {
 
     virtual amici::Model* clone() const override { return new Model_model_dirac(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         J_model_dirac(J, t, x, p, k, h, w, dwdx);
diff --git a/models/model_events/CMakeLists.txt b/models/model_events/CMakeLists.txt
index c60998f9e6..150dd43b02 100644
--- a/models/model_events/CMakeLists.txt
+++ b/models/model_events/CMakeLists.txt
@@ -66,6 +66,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_events/main.cpp b/models/model_events/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_events/main.cpp
+++ b/models/model_events/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_events/model_events.h b/models/model_events/model_events.h
index 1ea53269cf..f3d592eb5e 100644
--- a/models/model_events/model_events.h
+++ b/models/model_events/model_events.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_events_h
 #define _amici_model_events_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -75,7 +75,7 @@ class Model_model_events : public amici::Model_ODE {
 
     virtual amici::Model* clone() const override { return new Model_model_events(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         J_model_events(J, t, x, p, k, h, w, dwdx);
diff --git a/models/model_jakstat_adjoint/CMakeLists.txt b/models/model_jakstat_adjoint/CMakeLists.txt
index 27f5fd7c62..76484be708 100644
--- a/models/model_jakstat_adjoint/CMakeLists.txt
+++ b/models/model_jakstat_adjoint/CMakeLists.txt
@@ -55,6 +55,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_jakstat_adjoint/main.cpp b/models/model_jakstat_adjoint/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_jakstat_adjoint/main.cpp
+++ b/models/model_jakstat_adjoint/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint.h b/models/model_jakstat_adjoint/model_jakstat_adjoint.h
index e512ee9394..3cf5fdec25 100644
--- a/models/model_jakstat_adjoint/model_jakstat_adjoint.h
+++ b/models/model_jakstat_adjoint/model_jakstat_adjoint.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_jakstat_adjoint_h
 #define _amici_model_jakstat_adjoint_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -64,7 +64,7 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE {
 
     virtual amici::Model* clone() const override { return new Model_model_jakstat_adjoint(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         J_model_jakstat_adjoint(J, t, x, p, k, h, w, dwdx);
diff --git a/models/model_jakstat_adjoint_o2/CMakeLists.txt b/models/model_jakstat_adjoint_o2/CMakeLists.txt
index 3fad280043..92516553d0 100644
--- a/models/model_jakstat_adjoint_o2/CMakeLists.txt
+++ b/models/model_jakstat_adjoint_o2/CMakeLists.txt
@@ -55,6 +55,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_jakstat_adjoint_o2/main.cpp b/models/model_jakstat_adjoint_o2/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_jakstat_adjoint_o2/main.cpp
+++ b/models/model_jakstat_adjoint_o2/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h
index 52f25fccd9..ef9f6a1638 100644
--- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h
+++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_jakstat_adjoint_o2_h
 #define _amici_model_jakstat_adjoint_o2_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -64,7 +64,7 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE {
 
     virtual amici::Model* clone() const override { return new Model_model_jakstat_adjoint_o2(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         J_model_jakstat_adjoint_o2(J, t, x, p, k, h, w, dwdx);
diff --git a/models/model_nested_events/CMakeLists.txt b/models/model_nested_events/CMakeLists.txt
index b5fa8de311..120e5a08ed 100644
--- a/models/model_nested_events/CMakeLists.txt
+++ b/models/model_nested_events/CMakeLists.txt
@@ -55,6 +55,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_nested_events/main.cpp b/models/model_nested_events/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_nested_events/main.cpp
+++ b/models/model_nested_events/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_nested_events/model_nested_events.h b/models/model_nested_events/model_nested_events.h
index 76a3c9b80f..5ca7091d6f 100644
--- a/models/model_nested_events/model_nested_events.h
+++ b/models/model_nested_events/model_nested_events.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_nested_events_h
 #define _amici_model_nested_events_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -64,7 +64,7 @@ class Model_model_nested_events : public amici::Model_ODE {
 
     virtual amici::Model* clone() const override { return new Model_model_nested_events(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         J_model_nested_events(J, t, x, p, k, h, w, dwdx);
diff --git a/models/model_neuron/CMakeLists.txt b/models/model_neuron/CMakeLists.txt
index b7af8a3d7d..0c63349ca8 100644
--- a/models/model_neuron/CMakeLists.txt
+++ b/models/model_neuron/CMakeLists.txt
@@ -69,6 +69,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_neuron/main.cpp b/models/model_neuron/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_neuron/main.cpp
+++ b/models/model_neuron/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_neuron/model_neuron.h b/models/model_neuron/model_neuron.h
index ddf2f6e014..fa33f518ab 100644
--- a/models/model_neuron/model_neuron.h
+++ b/models/model_neuron/model_neuron.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_neuron_h
 #define _amici_model_neuron_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5dfc4abddf5cb50611e9881fcab074d8859ccc6b */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -78,7 +78,7 @@ class Model_model_neuron : public amici::Model_ODE {
 
     virtual amici::Model* clone() const override { return new Model_model_neuron(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5dfc4abddf5cb50611e9881fcab074d8859ccc6b"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         J_model_neuron(J, t, x, p, k, h, w, dwdx);
diff --git a/models/model_neuron_o2/CMakeLists.txt b/models/model_neuron_o2/CMakeLists.txt
index 7ad5944d56..7c2121cb0f 100644
--- a/models/model_neuron_o2/CMakeLists.txt
+++ b/models/model_neuron_o2/CMakeLists.txt
@@ -71,6 +71,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_neuron_o2/main.cpp b/models/model_neuron_o2/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_neuron_o2/main.cpp
+++ b/models/model_neuron_o2/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_neuron_o2/model_neuron_o2.h b/models/model_neuron_o2/model_neuron_o2.h
index ef16e543d1..a4fb5d03d9 100644
--- a/models/model_neuron_o2/model_neuron_o2.h
+++ b/models/model_neuron_o2/model_neuron_o2.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_neuron_o2_h
 #define _amici_model_neuron_o2_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5dfc4abddf5cb50611e9881fcab074d8859ccc6b */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -80,7 +80,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE {
 
     virtual amici::Model* clone() const override { return new Model_model_neuron_o2(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5dfc4abddf5cb50611e9881fcab074d8859ccc6b"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         J_model_neuron_o2(J, t, x, p, k, h, w, dwdx);
diff --git a/models/model_robertson/CMakeLists.txt b/models/model_robertson/CMakeLists.txt
index 072011db85..7af14aba5b 100644
--- a/models/model_robertson/CMakeLists.txt
+++ b/models/model_robertson/CMakeLists.txt
@@ -53,6 +53,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_robertson/main.cpp b/models/model_robertson/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_robertson/main.cpp
+++ b/models/model_robertson/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_robertson/model_robertson.h b/models/model_robertson/model_robertson.h
index 556cea6e75..c526ad14ad 100644
--- a/models/model_robertson/model_robertson.h
+++ b/models/model_robertson/model_robertson.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_robertson_h
 #define _amici_model_robertson_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -62,7 +62,7 @@ class Model_model_robertson : public amici::Model_DAE {
 
     virtual amici::Model* clone() const override { return new Model_model_robertson(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override {
         J_model_robertson(J, t, x, p, k, h, cj, dx, w, dwdx);
diff --git a/models/model_steadystate/CMakeLists.txt b/models/model_steadystate/CMakeLists.txt
index 0537b4b5b1..0b14584260 100644
--- a/models/model_steadystate/CMakeLists.txt
+++ b/models/model_steadystate/CMakeLists.txt
@@ -52,6 +52,13 @@ ${MODEL_DIR}/wrapfunctions.cpp
 add_library(${PROJECT_NAME} ${SRC_LIST_LIB})
 add_library(model ALIAS ${PROJECT_NAME})
 
+# This option can be helpful when using the Intel compiler and compilation of
+# wrapfunctions.cpp fails due to insufficient memory.
+option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON)
+if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS)
+    set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
 target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/models/model_steadystate/main.cpp b/models/model_steadystate/main.cpp
index 6d6caf54a8..07ff62a322 100644
--- a/models/model_steadystate/main.cpp
+++ b/models/model_steadystate/main.cpp
@@ -10,8 +10,7 @@
 #include "wrapfunctions.h"  /* model-provided functions */
 
 /* This is a scaffold for a stand-alone AMICI simulation executable
- * demonstrating
- * use of the AMICI C++ API.
+ * demonstrating the use of the AMICI C++ API.
  *
  * This program reads AMICI options from an HDF5 file, prints some results
  * and writes additional results to an HDF5 file. The name of the HDF5 file
diff --git a/models/model_steadystate/model_steadystate.h b/models/model_steadystate/model_steadystate.h
index 06a1c2773b..4489e3a664 100644
--- a/models/model_steadystate/model_steadystate.h
+++ b/models/model_steadystate/model_steadystate.h
@@ -1,6 +1,6 @@
 #ifndef _amici_model_steadystate_h
 #define _amici_model_steadystate_h
-/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */
+/* Generated by amiwrap (R2017b) 5dfc4abddf5cb50611e9881fcab074d8859ccc6b */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -61,7 +61,7 @@ class Model_model_steadystate : public amici::Model_ODE {
 
     virtual amici::Model* clone() const override { return new Model_model_steadystate(*this); };
 
-    const  std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; };
+    const  std::string getAmiciCommit() const override { return "5dfc4abddf5cb50611e9881fcab074d8859ccc6b"; };
 
     virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         J_model_steadystate(J, t, x, p, k, h, w, dwdx);
diff --git a/src/amici.cpp b/src/amici.cpp
index d71a03c9eb..15e8a98852 100644
--- a/src/amici.cpp
+++ b/src/amici.cpp
@@ -9,15 +9,16 @@
 #include "amici/forwardproblem.h"
 #include "amici/misc.h"
 
+#include <cvodes/cvodes.h>           //return codes
 #include <sundials/sundials_types.h> //realtype
-#include <cvodes/cvodes.h> //return codes
 
-#include <type_traits>
 #include <cassert>
+#include <cstdarg>
 #include <cstdlib>
 #include <cstring>
-#include <cstdarg>
+#include <iostream>
 #include <memory>
+#include <type_traits>
 
 // ensure definitions are in sync
 static_assert(AMICI_SUCCESS == CV_SUCCESS, "AMICI_SUCCESS != CV_SUCCESS");
@@ -32,53 +33,120 @@ static_assert(AMICI_ONE_STEP == CV_ONE_STEP, "AMICI_ONE_STEP != CV_ONE_STEP");
 static_assert(std::is_same<amici::realtype, realtype>::value,
               "Definition of realtype does not match");
 
-
 namespace amici {
 
-msgIdAndTxtFp errMsgIdAndTxt = &printErrMsgIdAndTxt;
-msgIdAndTxtFp warnMsgIdAndTxt = &printWarnMsgIdAndTxt;
+/** AMICI default application context, kept around for convenience for using
+  * amici::runAmiciSimulation or instantiating Solver and Model without special
+  * needs.
+  */
+AmiciApplication defaultContext = AmiciApplication();
+
+std::unique_ptr<ReturnData>
+runAmiciSimulation(Solver& solver,
+                   const ExpData* edata,
+                   Model& model,
+                   bool rethrow)
+{
+    return defaultContext.runAmiciSimulation(solver, edata, model, rethrow);
+}
+
+void
+printErrMsgIdAndTxt(std::string const& id, std::string const& message)
+{
+    std::cerr << "[Error] ";
+    if (!id.empty()) {
+        std::cerr << id << ": ";
+    }
+    std::cerr << message << std::endl;
+}
+
+void
+printWarnMsgIdAndTxt(std::string const& id, std::string const& message)
+{
+    std::cerr << "[Warning] ";
+    if (!id.empty()) {
+        std::cerr << id << ": ";
+    }
+    std::cerr << message << std::endl;
+}
 
+std::vector<std::unique_ptr<ReturnData>>
+runAmiciSimulations(const Solver& solver,
+                    const std::vector<ExpData*>& edatas,
+                    const Model& model,
+                    const bool failfast,
+#if defined(_OPENMP)
+                    int num_threads
+#else
+                    int /* num_threads */
+#endif
+)
+{
+#if defined(_OPENMP)
+    return defaultContext.runAmiciSimulations(
+      solver, edatas, model, failfast, num_threads);
+#else
+    return defaultContext.runAmiciSimulations(solver, edatas, model, failfast, 1);
+#endif
+}
 
-std::unique_ptr<ReturnData> runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model, bool rethrow) {
+std::unique_ptr<ReturnData>
+AmiciApplication::runAmiciSimulation(Solver& solver,
+                                     const ExpData* edata,
+                                     Model& model,
+                                     bool rethrow)
+{
     std::unique_ptr<ReturnData> rdata;
 
     /* Applies condition-specific model settings and restores them when going
      * out of scope */
     ConditionContext conditionContext(&model, edata);
 
-    try{
+    try {
         rdata = std::unique_ptr<ReturnData>(new ReturnData(solver, model));
 
         if (model.nx_solver <= 0) {
             return rdata;
         }
 
-        auto fwd = std::unique_ptr<ForwardProblem>(new ForwardProblem(rdata.get(),edata,&model,&solver));
+        auto fwd = std::unique_ptr<ForwardProblem>(
+          new ForwardProblem(rdata.get(), edata, &model, &solver));
         fwd->workForwardProblem();
 
-        auto bwd = std::unique_ptr<BackwardProblem>(new BackwardProblem(fwd.get()));
+        auto bwd =
+          std::unique_ptr<BackwardProblem>(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());
+        if (rethrow)
+            throw;
+        warningF("AMICI:mex:simulation",
+                 "AMICI forward simulation failed at t = %f:\n%s\n",
+                 ex.time,
+                 ex.what());
     } catch (amici::IntegrationFailureB const& ex) {
         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"
-                    " (see message above):\n%s\n",
-                    ex.time, ex.what());
+        if (rethrow)
+            throw;
+        warningF(
+          "AMICI:mex:simulation",
+          "AMICI backward simulation failed when trying to solve until t = %f"
+          " (see message above):\n%s\n",
+          ex.time,
+          ex.what());
     } 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());
+        if (rethrow)
+            throw;
+        warningF("AMICI:mex:simulation",
+                 "AMICI simulation failed:\n%s\nError occured in:\n%s",
+                 ex.what(),
+                 ex.getBacktrace());
     }
 
     rdata->applyChainRuleFactorToSimulationResults(&model);
@@ -86,50 +154,28 @@ std::unique_ptr<ReturnData> runAmiciSimulation(Solver &solver, const ExpData *ed
     return rdata;
 }
 
-void printErrMsgIdAndTxt(const char *identifier, const char *format, ...) {
-    if(identifier && *identifier != '\0')
-        fprintf(stderr, "[Error] %s: ", identifier);
-    else
-        fprintf(stderr, "[Error] ");
-    va_list argptr;
-    va_start(argptr,format);
-    vfprintf(stderr, format, argptr);
-    va_end(argptr);
-    fprintf(stderr, "\n");
-}
-
-void printWarnMsgIdAndTxt(const char *identifier, const char *format, ...) {
-    if(identifier && *identifier != '\0')
-        printf("[Warning] %s: ", identifier);
-    else
-        printf("[Warning] ");
-    va_list argptr;
-    va_start(argptr,format);
-    vprintf(format, argptr);
-    va_end(argptr);
-    printf("\n");
-}
-
-std::vector<std::unique_ptr<ReturnData> > runAmiciSimulations(const Solver &solver,
-                                                              const std::vector<ExpData*> &edatas,
-                                                              const Model &model,
-                                                              const bool failfast,
+std::vector<std::unique_ptr<ReturnData>>
+AmiciApplication::runAmiciSimulations(const Solver& solver,
+                                      const std::vector<ExpData*>& edatas,
+                                      const Model& model,
+                                      bool failfast,
 #if defined(_OPENMP)
-                                                              int num_threads
+                                      int num_threads
 #else
-                                                              int /* num_threads */
+                                      int /* num_threads */
 #endif
+
 )
 {
-    std::vector<std::unique_ptr<ReturnData> > results(edatas.size());
+    std::vector<std::unique_ptr<ReturnData>> results(edatas.size());
     // is set to true if one simulation fails and we should skip the rest.
     // shared across threads.
     bool skipThrough = false;
 
 #if defined(_OPENMP)
-    #pragma omp parallel for num_threads(num_threads)
+#pragma omp parallel for num_threads(num_threads)
 #endif
-    for(int i = 0; i < (int)edatas.size(); ++i) {
+    for (int i = 0; i < (int)edatas.size(); ++i) {
         auto mySolver = std::unique_ptr<Solver>(solver.clone());
         auto myModel = std::unique_ptr<Model>(model.clone());
 
@@ -138,7 +184,7 @@ std::vector<std::unique_ptr<ReturnData> > runAmiciSimulations(const Solver &solv
         if (skipThrough) {
             ConditionContext conditionContext(myModel.get(), edatas[i]);
             results[i] =
-                std::unique_ptr<ReturnData>(new ReturnData(solver, model));
+              std::unique_ptr<ReturnData>(new ReturnData(solver, model));
         } else {
             results[i] = runAmiciSimulation(*mySolver, edatas[i], *myModel);
         }
@@ -149,4 +195,49 @@ std::vector<std::unique_ptr<ReturnData> > runAmiciSimulations(const Solver &solv
     return results;
 }
 
+void
+AmiciApplication::warningF(const char* identifier, const char* format, ...)
+{
+    va_list argptr;
+    va_start(argptr, format);
+    auto str = printfToString(format, argptr);
+    va_end(argptr);
+    warning(identifier, str);
+}
+
+void
+AmiciApplication::errorF(const char* identifier, const char* format, ...)
+{
+    va_list argptr;
+    va_start(argptr, format);
+    auto str = printfToString(format, argptr);
+    va_end(argptr);
+    error(identifier, str);
+}
+
+int
+AmiciApplication::checkFinite(gsl::span<const realtype> array, const char* fun)
+{
+
+    for (int idx = 0; idx < (int)array.size(); idx++) {
+        if (isNaN(array[idx])) {
+            warningF("AMICI:NaN",
+                     "AMICI encountered a NaN value at index %i of %i in %s!",
+                     idx,
+                     (int)array.size(),
+                     fun);
+            return AMICI_RECOVERABLE_ERROR;
+        }
+        if (isInf(array[idx])) {
+            warningF("AMICI:Inf",
+                     "AMICI encountered an Inf value at index %i of %i in %s!",
+                     idx,
+                     (int)array.size(),
+                     fun);
+            return AMICI_RECOVERABLE_ERROR;
+        }
+    }
+    return AMICI_SUCCESS;
+}
+
 } // namespace amici
diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp
index 52112e5670..5b8ff882f3 100644
--- a/src/forwardproblem.cpp
+++ b/src/forwardproblem.cpp
@@ -15,9 +15,6 @@
 
 namespace amici {
 
-extern msgIdAndTxtFp warnMsgIdAndTxt;
-
-
 ForwardProblem::ForwardProblem(ReturnData *rdata, const ExpData *edata,
                                Model *model, Solver *solver)
     : model(model),
@@ -310,7 +307,8 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
         discs[iroot] = t;
         ++iroot;
     } else {
-        warnMsgIdAndTxt("AMICI:mex:TOO_MUCH_EVENT",
+        solver->app->warning(
+                    "AMICI:mex:TOO_MUCH_EVENT",
                         "Event was recorded but not reported as the number of "
                         "occured events exceeded (nmaxevents)*(number of "
                         "events in model definition)!");
diff --git a/src/interface_matlab.cpp b/src/interface_matlab.cpp
index 5f785d37ee..800bbc061c 100644
--- a/src/interface_matlab.cpp
+++ b/src/interface_matlab.cpp
@@ -462,17 +462,31 @@ void setModelData(const mxArray *prhs[], int nrhs, Model &model)
  */
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
     // use matlab error reporting
-    amici::warnMsgIdAndTxt = &mexWarnMsgIdAndTxt;
-    amici::errMsgIdAndTxt = &mexErrMsgIdAndTxt;
+    amici::AmiciApplication amiciApp;
+    amiciApp.warning = [](
+            std::string const& identifier,
+            std::string const& message){
+        mexWarnMsgIdAndTxt(identifier.c_str(), message.c_str());
+    };
+    amiciApp.error = [](
+            std::string const& identifier,
+            std::string const& message){
+        mexErrMsgIdAndTxt(identifier.c_str(), message.c_str());
+    };
 
     if (nlhs != 1) {
-        amici::errMsgIdAndTxt("AMICI:mex:setup","Incorrect number of output arguments (must be 1)!");
+        amiciApp.errorF("AMICI:mex:setup",
+                        "Incorrect number of output arguments (must be 1)!");
     } else if(nrhs < amici::RHS_NUMARGS_REQUIRED) {
-        amici::errMsgIdAndTxt("AMICI:mex:setup", "Incorrect number of input arguments (must be at least 7)!");
+        amiciApp.errorF("AMICI:mex:setup",
+                        "Incorrect number of input arguments (must be at least 7)!");
     };
 
     auto model = getModel();
+    model->app = &amiciApp;
+
     auto solver = model->getSolver();
+    solver->app = &amiciApp;
     setModelData(prhs, nrhs, *model);
     setSolverOptions(prhs, nrhs, *solver);
 
@@ -481,11 +495,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
         try {
             edata = amici::expDataFromMatlabCall(prhs, *model);
         } catch (amici::AmiException const& ex) {
-            amici::errMsgIdAndTxt("AMICI:mex:setup","Failed to read experimental data:\n%s",ex.what());
+            amiciApp.errorF("AMICI:mex:setup","Failed to read experimental data:\n%s",ex.what());
         }
     } else if (solver->getSensitivityOrder() >= amici::SensitivityOrder::first &&
                solver->getSensitivityMethod() == amici::SensitivityMethod::adjoint) {
-        amici::errMsgIdAndTxt("AMICI:mex:setup","No data provided!");
+        amiciApp.errorF("AMICI:mex:setup","No data provided!");
     }
 
     /* ensures that plhs[0] is available */
diff --git a/src/misc.cpp b/src/misc.cpp
index 3876c7cbfd..731559c96b 100644
--- a/src/misc.cpp
+++ b/src/misc.cpp
@@ -5,6 +5,8 @@
 #include <cstdio>
 #include <cstring>
 #include <sstream>
+#include <cstdarg>
+
 #if defined(_WIN32)
 #define PLATFORM_WINDOWS // Windows
 #elif defined(_WIN64)
@@ -86,27 +88,6 @@ void scaleParameters(gsl::span<const realtype> bufferUnscaled,
 
 }
 
-int checkFinite(gsl::span<const realtype> array, const char *fun)
-{
-    for (int idx = 0; idx < (int) array.size(); idx++) {
-        if (isNaN(array[idx])) {
-            warnMsgIdAndTxt(
-                "AMICI:NaN",
-                "AMICI encountered a NaN value at index %i of %i in %s!", idx,
-                (int) array.size(), fun);
-            return AMICI_RECOVERABLE_ERROR;
-        }
-        if (isInf(array[idx])) {
-            warnMsgIdAndTxt(
-                "AMICI:Inf",
-                "AMICI encountered an Inf value at index %i of %i in %s!", idx,
-                (int) array.size(), fun);
-            return AMICI_RECOVERABLE_ERROR;
-        }
-    }
-    return AMICI_SUCCESS;
-}
-
 std::string backtraceString(const int maxFrames)
 {
     std::ostringstream trace_buf;
@@ -186,4 +167,21 @@ std::string regexErrorToString(std::regex_constants::error_type err_type)
     }
 }
 
+std::string printfToString(const char *fmt, va_list ap) {
+    // Get size of string
+    va_list ap_count;
+    va_copy(ap_count, ap);
+    auto size = vsnprintf(nullptr, 0, fmt, ap_count);
+    va_end(ap_count);
+    ++size;
+
+    // actual formatting
+    auto buf = new char[size];
+    size = vsnprintf(buf, size, fmt, ap);
+    std::string str(buf, size);
+    delete[] buf;
+
+    return str;
+}
+
 } // namespace amici
diff --git a/src/model.cpp b/src/model.cpp
index 9a35c36f83..24c14ee728 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -1002,7 +1002,7 @@ void Model::addStateEventUpdate(AmiVector &x, const int ie, const realtype t,
             fixedParameters.data(), h.data(), ie, xdot.data(), xdot_old.data());
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(deltax, "deltax");
+        app->checkFinite(deltax, "deltax");
     }
 
     // update
@@ -1027,7 +1027,7 @@ void Model::addStateSensitivityEventUpdate(AmiVectorArray &sx, const int ie,
                  xdot.data(), xdot_old.data(), sx.data(ip), &stau.at(ip));
 
         if (alwaysCheckFinite) {
-            amici::checkFinite(deltasx, "deltasx");
+            app->checkFinite(deltasx, "deltasx");
         }
 
         amici_daxpy(nx_solver, 1.0, deltasx.data(), 1, sx.data(ip), 1);
@@ -1047,7 +1047,7 @@ void Model::addAdjointStateEventUpdate(AmiVector &xB, const int ie,
              xB.data());
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(deltaxB, "deltaxB");
+        app->checkFinite(deltaxB, "deltaxB");
     }
 
     // apply update
@@ -1072,7 +1072,7 @@ void Model::addAdjointQuadratureEventUpdate(
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(deltaqB, "deltaqB");
+        app->checkFinite(deltaqB, "deltaqB");
     }
 }
 
@@ -1089,12 +1089,12 @@ void Model::updateHeavisideB(const int *rootsfound) {
 }
 
 int Model::checkFinite(gsl::span<const realtype> array, const char *fun) const {
-    auto result = amici::checkFinite(array, fun);
+    auto result = app->checkFinite(array, fun);
 
     if (result != AMICI_SUCCESS) {
-        amici::checkFinite(fixedParameters, "k");
-        amici::checkFinite(unscaledParameters, "p");
-        amici::checkFinite(w, "w");
+        app->checkFinite(fixedParameters, "k");
+        app->checkFinite(unscaledParameters, "p");
+        app->checkFinite(w, "w");
     }
 
     return result;
@@ -1240,7 +1240,7 @@ void Model::fy(const realtype t, const AmiVector &x) {
        h.data(), w.data());
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(gsl::make_span(y.data(), ny), "y");
+        app->checkFinite(gsl::make_span(y.data(), ny), "y");
     }
 }
 
@@ -1265,7 +1265,7 @@ void Model::fdydp(const realtype t, const AmiVector &x) {
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dydp, "dydp");
+        app->checkFinite(dydp, "dydp");
     }
 }
 
@@ -1281,7 +1281,7 @@ void Model::fdydx(const realtype t, const AmiVector &x) {
           fixedParameters.data(), h.data(), w.data(), dwdx.data());
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dydx, "dydx");
+        app->checkFinite(dydx, "dydx");
     }
 }
 
@@ -1339,7 +1339,7 @@ void Model::fdsigmaydp(const int it, const ExpData *edata) {
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dsigmaydp, "dsigmaydp");
+        app->checkFinite(dsigmaydp, "dsigmaydp");
     }
 }
 
@@ -1375,7 +1375,7 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) {
                    edata.getObservedDataPtr(it));
 
             if (alwaysCheckFinite) {
-                amici::checkFinite(gsl::make_span(dJydy[iyt].get()), "dJydy");
+                app->checkFinite(gsl::make_span(dJydy[iyt].get()), "dJydy");
             }
         }
     } else {
@@ -1389,7 +1389,7 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) {
         }
         if (alwaysCheckFinite) {
             // get dJydy slice (ny) for current timepoint and observable
-            amici::checkFinite(dJydy_matlab, "dJydy");
+            app->checkFinite(dJydy_matlab, "dJydy");
         }
     }
 }
@@ -1410,7 +1410,7 @@ void Model::fdJydsigma(const int it, const AmiVector &x, const ExpData &edata) {
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dJydsigma, "dJydsigma");
+        app->checkFinite(dJydsigma, "dJydsigma");
     }
 }
 
@@ -1486,7 +1486,7 @@ void Model::fdJydx(const int it, const AmiVector &x, const ExpData &edata) {
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dJydx, "dJydx");
+        app->checkFinite(dJydx, "dJydx");
     }
 }
 
@@ -1508,7 +1508,7 @@ void Model::fdzdp(const int ie, const realtype t, const AmiVector &x) {
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dzdp, "dzdp");
+        app->checkFinite(dzdp, "dzdp");
     }
 }
 
@@ -1520,7 +1520,7 @@ void Model::fdzdx(const int ie, const realtype t, const AmiVector &x) {
           fixedParameters.data(), h.data());
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dzdx, "dzdx");
+        app->checkFinite(dzdx, "dzdx");
     }
 }
 
@@ -1542,7 +1542,7 @@ void Model::fdrzdp(const int ie, const realtype t, const AmiVector &x) {
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(drzdp, "drzdp");
+        app->checkFinite(drzdp, "drzdp");
     }
 }
 
@@ -1554,7 +1554,7 @@ void Model::fdrzdx(const int ie, const realtype t, const AmiVector &x) {
            fixedParameters.data(), h.data());
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(drzdx, "drzdx");
+        app->checkFinite(drzdx, "drzdx");
     }
 }
 
@@ -1613,7 +1613,7 @@ void Model::fdsigmazdp(const int ie, const int nroots, const realtype t,
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dsigmazdp, "dsigmazdp");
+        app->checkFinite(dsigmazdp, "dsigmazdp");
     }
 }
 
@@ -1634,7 +1634,7 @@ void Model::fdJzdz(const int ie, const int nroots, const realtype t,
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dJzdz, "dJzdz");
+        app->checkFinite(dJzdz, "dJzdz");
     }
 }
 
@@ -1656,7 +1656,7 @@ void Model::fdJzdsigma(const int ie, const int nroots, const realtype t,
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dJzdsigma, "dJzdsigma");
+        app->checkFinite(dJzdsigma, "dJzdsigma");
     }
 }
 
@@ -1759,7 +1759,7 @@ void Model::fdJrzdz(const int ie, const int nroots, const realtype t,
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dJrzdz, "dJrzdz");
+        app->checkFinite(dJrzdz, "dJrzdz");
     }
 }
 
@@ -1780,7 +1780,7 @@ void Model::fdJrzdsigma(const int ie, const int nroots, const realtype t,
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dJrzdsigma, "dJrzdsigma");
+        app->checkFinite(dJrzdsigma, "dJrzdsigma");
     }
 }
 
@@ -1790,7 +1790,7 @@ void Model::fw(const realtype t, const realtype *x) {
        h.data(), total_cl.data());
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(w, "w");
+        app->checkFinite(w, "w");
     }
 }
 
@@ -1819,7 +1819,7 @@ void Model::fdwdp(const realtype t, const realtype *x) {
     }
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(dwdp, "dwdp");
+        app->checkFinite(dwdp, "dwdp");
     }
 }
 
@@ -1832,7 +1832,7 @@ void Model::fdwdx(const realtype t, const realtype *x) {
     fdwdx_rowvals(dwdx.indexptrs());
 
     if (alwaysCheckFinite) {
-        amici::checkFinite(gsl::make_span(dwdx.get()), "dwdx");
+        app->checkFinite(gsl::make_span(dwdx.get()), "dwdx");
     }
 }
 
diff --git a/src/solver.cpp b/src/solver.cpp
index 2f38c0d0e2..f86a896fa9 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -12,7 +12,10 @@
 
 namespace amici {
 
-extern msgIdAndTxtFp warnMsgIdAndTxt;
+Solver::Solver(AmiciApplication *app) : app(app)
+{
+
+}
 
 Solver::Solver(const Solver &other)
     : ism(other.ism), lmm(other.lmm), iter(other.iter),
@@ -1016,38 +1019,44 @@ const AmiVector &Solver::getAdjointQuadrature(const int which,
 realtype Solver::gett() const { return t; }
 
 void wrapErrHandlerFn(int error_code, const char *module,
-                      const char *function, char *msg, void * /*eh_data*/) {
-    char buffer[250];
-    char buffid[250];
-    sprintf(buffer, "AMICI ERROR: in module %s in function %s : %s ", module,
+                      const char *function, char *msg, void * eh_data) {
+#define BUF_SIZE 250
+    char buffer[BUF_SIZE];
+    char buffid[BUF_SIZE];
+    snprintf(buffer, BUF_SIZE, "AMICI ERROR: in module %s in function %s : %s ", module,
             function, msg);
     switch (error_code) {
     case 99:
-        sprintf(buffid, "AMICI:mex:%s:%s:WARNING", module, function);
+        snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:WARNING", module, function);
         break;
 
     case -1:
-        sprintf(buffid, "AMICI:mex:%s:%s:TOO_MUCH_WORK", module, function);
+        snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:TOO_MUCH_WORK", module, function);
         break;
 
     case -2:
-        sprintf(buffid, "AMICI:mex:%s:%s:TOO_MUCH_ACC", module, function);
+        snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:TOO_MUCH_ACC", module, function);
         break;
 
     case -3:
-        sprintf(buffid, "AMICI:mex:%s:%s:ERR_FAILURE", module, function);
+        snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:ERR_FAILURE", module, function);
         break;
 
     case -4:
-        sprintf(buffid, "AMICI:mex:%s:%s:CONV_FAILURE", module, function);
+        snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:CONV_FAILURE", module, function);
         break;
 
     default:
-        sprintf(buffid, "AMICI:mex:%s:%s:OTHER", module, function);
+        snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:OTHER", module, function);
         break;
     }
 
-    warnMsgIdAndTxt(buffid, buffer);
+
+    if(!eh_data) {
+        throw std::runtime_error("eh_data unset");
+    }
+    auto solver = static_cast<Solver const*>(eh_data);
+    solver->app->warning(buffid, buffer);
 }
 
 } // namespace amici
diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp
index 56441f7828..34fbfe9ef1 100644
--- a/src/solver_cvodes.cpp
+++ b/src/solver_cvodes.cpp
@@ -307,7 +307,9 @@ void CVodeSolver::setNonLinearSolverB(const int which) const {
 
 void CVodeSolver::setErrHandlerFn() const {
     int status =
-        CVodeSetErrHandlerFn(solverMemory.get(), wrapErrHandlerFn, nullptr);
+        CVodeSetErrHandlerFn(solverMemory.get(), wrapErrHandlerFn,
+                             reinterpret_cast<void*>(
+                                 const_cast<CVodeSolver*>(this)));
     if (status != CV_SUCCESS)
         throw CvodeException(status, "CVodeSetErrHandlerFn");
 }
@@ -342,9 +344,9 @@ void CVodeSolver::setStabLimDetB(const int which, const int stldet) const {
         throw CvodeException(status, "CVodeSetStabLimDetB");
 }
 
-void CVodeSolver::setId(const Model *model) const {}
+void CVodeSolver::setId(const Model */*model*/) const {}
 
-void CVodeSolver::setSuppressAlg(const bool flag) const {}
+void CVodeSolver::setSuppressAlg(const bool /*flag*/) const {}
 
 void CVodeSolver::resetState(void *ami_mem, const_N_Vector y0) const {
 
@@ -570,7 +572,7 @@ void CVodeSolver::allocateSolverB(int *which) const {
         solverMemoryB.resize(*which + 1);
     solverMemoryB.at(*which) =
         std::unique_ptr<void, std::function<void(void *)>>(
-            getAdjBmem(solverMemory.get(), *which), [](void *ptr) {});
+            getAdjBmem(solverMemory.get(), *which), [](void */*ptr*/) {});
     if (status != CV_SUCCESS)
         throw CvodeException(status, "CVodeCreateB");
 }
@@ -686,9 +688,9 @@ void *CVodeSolver::getAdjBmem(void *ami_mem, int which) const {
     return CVodeGetAdjCVodeBmem(ami_mem, which);
 }
 
-void CVodeSolver::calcIC(const realtype tout1) const {};
+void CVodeSolver::calcIC(const realtype /*tout1*/) const {};
 
-void CVodeSolver::calcICB(const int which, const realtype tout1) const {};
+void CVodeSolver::calcICB(const int /*which*/, const realtype /*tout1*/) const {};
 
 void CVodeSolver::setStopTime(const realtype tstop) const {
     int status = CVodeSetStopTime(solverMemory.get(), tstop);
@@ -903,7 +905,7 @@ int froot(realtype t, N_Vector x, realtype *root,
 int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data) {
     auto model = static_cast<Model_ODE *>(user_data);
 
-    if (t > 1e200 && !amici::checkFinite(gsl::make_span(x), "fxdot")) {
+    if (t > 1e200 && !model->checkFinite(gsl::make_span(x), "fxdot")) {
         /* when t is large (typically ~1e300), CVODES may pass all NaN x
            to fxdot from which we typically cannot recover. To save time
            on normal execution, we do not always want to check finiteness
diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp
index 670b22805b..a8b9b51e4f 100644
--- a/src/solver_idas.cpp
+++ b/src/solver_idas.cpp
@@ -244,7 +244,9 @@ void IDASolver::getRootInfo(int *rootsfound) const {
 
 void IDASolver::setErrHandlerFn() const {
     int status =
-        IDASetErrHandlerFn(solverMemory.get(), wrapErrHandlerFn, nullptr);
+        IDASetErrHandlerFn(solverMemory.get(), wrapErrHandlerFn,
+                           reinterpret_cast<void*>(
+                               const_cast<IDASolver*>(this)));
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDASetErrHandlerFn");
 }
@@ -267,9 +269,10 @@ void IDASolver::setMaxNumSteps(const long int mxsteps) const {
         throw IDAException(status, "IDASetMaxNumSteps");
 }
 
-void IDASolver::setStabLimDet(const int stldet) const {}
+void IDASolver::setStabLimDet(const int /*stldet*/) const {}
 
-void IDASolver::setStabLimDetB(const int which, const int stldet) const {}
+void IDASolver::setStabLimDetB(const int /*which*/,
+                               const int /*stldet*/) const {}
 
 void IDASolver::setId(const Model *model) const {
 
@@ -511,7 +514,7 @@ void IDASolver::allocateSolverB(int *which) const {
         solverMemoryB.resize(*which + 1);
     solverMemoryB.at(*which) =
         std::unique_ptr<void, std::function<void(void *)>>(
-            getAdjBmem(solverMemory.get(), *which), [](void *ptr) {});
+            getAdjBmem(solverMemory.get(), *which), [](void */*ptr*/) {});
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDACreateB");
 }
@@ -938,7 +941,7 @@ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot,
                      void *user_data) {
     auto model = static_cast<Model_DAE *>(user_data);
 
-    if (t > 1e200 && !amici::checkFinite(gsl::make_span(x), "fxdot")) {
+    if (t > 1e200 && !model->app->checkFinite(gsl::make_span(x), "fxdot")) {
         /* when t is large (typically ~1e300), CVODES may pass all NaN x
            to fxdot from which we typically cannot recover. To save time
            on normal execution, we do not always want to check finiteness
diff --git a/swig/amici.i b/swig/amici.i
index f57eeb8e50..359520677c 100644
--- a/swig/amici.i
+++ b/swig/amici.i
@@ -87,8 +87,8 @@ wrap_unique_ptr(ExpDataPtr, amici::ExpData)
 
 // Add necessary symbols to generated header
 // Ignore due to https://github.com/swig/swig/issues/1643
-%ignore printErrMsgIdAndTxt;
-%ignore printWarnMsgIdAndTxt;
+%ignore amici::AmiciApplication::warningF;
+%ignore amici::AmiciApplication::errorF;
 %{
 #include "amici/amici.h"
 using namespace amici;

From 78e297ffba30cfff237816c02bace3ae6be26c02 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@helmholtz-muenchen.de>
Date: Thu, 28 Nov 2019 11:43:18 +0100
Subject: [PATCH 11/19] Fix choice of likelihood formula; failed if observable
 names were not equal to observable IDs

---
 python/amici/sbml_import.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index e563b77059..67c4fa7aaf 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -910,11 +910,13 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
             observableSyms = sp.Matrix([
                 sp.symbols(obs, real=True) for obs in observables.keys()
             ])
+            observable_ids = observables.keys()
         else:
             observableValues = speciesSyms
-            observableNames = [
+            observable_ids = [
                 f'x{index}' for index in range(len(speciesSyms))
             ]
+            observableNames = observable_ids[:]
             observableSyms = sp.Matrix(
                 [sp.symbols(f'y{index}', real=True)
                  for index in range(len(speciesSyms))]
@@ -942,7 +944,7 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
 
         # set cost functions
         llhYStrings = []
-        for y_name in observableNames:
+        for y_name in observable_ids:
             llhYStrings.append(noise_distribution_to_cost_function(
                 noise_distributions.get(y_name, 'normal')))
 

From 8b83309e1faa2351b1efe1108ec461401f48ba0d Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@helmholtz-muenchen.de>
Date: Thu, 28 Nov 2019 14:54:35 +0100
Subject: [PATCH 12/19] Fix(sbml-import) Compartment IDs in right-hand side of
 Rules are not replaced and lead to undefined identifiers in c++ files

---
 python/amici/sbml_import.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index 67c4fa7aaf..e9d7c55d25 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -577,10 +577,7 @@ def processCompartments(self):
                     locals=self.local_symbols
                 )
 
-        for comp, vol in zip(self.compartmentSymbols, self.compartmentVolume):
-            self.replaceInAllExpressions(
-               comp, vol
-            )
+
 
     def processReactions(self):
         """Get reactions from SBML model.
@@ -790,6 +787,10 @@ def processRules(self):
                 sp.sympify(variable, locals=self.local_symbols),
                 assignments[variable]
             )
+        for comp, vol in zip(self.compartmentSymbols, self.compartmentVolume):
+            self.replaceInAllExpressions(
+               comp, vol
+            )
 
     def processVolumeConversion(self):
         """Convert equations from amount to volume.

From 31b476a05a9599ac5102aa807921fbdcc3e293e8 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Fri, 29 Nov 2019 01:39:54 +0100
Subject: [PATCH 13/19] Fix choice of likelihood formula; failed if observable
 names were not equal to observable IDs (#868)

---
 python/amici/sbml_import.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index e563b77059..67c4fa7aaf 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -910,11 +910,13 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
             observableSyms = sp.Matrix([
                 sp.symbols(obs, real=True) for obs in observables.keys()
             ])
+            observable_ids = observables.keys()
         else:
             observableValues = speciesSyms
-            observableNames = [
+            observable_ids = [
                 f'x{index}' for index in range(len(speciesSyms))
             ]
+            observableNames = observable_ids[:]
             observableSyms = sp.Matrix(
                 [sp.symbols(f'y{index}', real=True)
                  for index in range(len(speciesSyms))]
@@ -942,7 +944,7 @@ def processObservables(self, observables: Dict[str, Dict[str, str]],
 
         # set cost functions
         llhYStrings = []
-        for y_name in observableNames:
+        for y_name in observable_ids:
             llhYStrings.append(noise_distribution_to_cost_function(
                 noise_distributions.get(y_name, 'normal')))
 

From 737e739f452a89b392e89aadc21fcab8bc4479a7 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@helmholtz-muenchen.de>
Date: Fri, 29 Nov 2019 11:12:06 +0100
Subject: [PATCH 14/19] Bump version number to v0.10.14

---
 version.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/version.txt b/version.txt
index f25c43cdc2..c70613aa09 100644
--- a/version.txt
+++ b/version.txt
@@ -1 +1 @@
-0.10.13
+0.10.14

From 723cbc854bf124a7ca139a0b0edc54c16180dc98 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Mon, 2 Dec 2019 15:21:27 +0100
Subject: [PATCH 15/19] Misc (#870)

* Fix invalid logging level

* Update publication list (Closes #867)

* Update reference list

* doc(install) Update installation guide for installing from GitHub using pip
---
 INSTALL.md                   | 14 +++++++
 documentation/amici_refs.bib | 72 ++++++++++++++++++++++++++++--------
 documentation/references.md  | 17 +++++++--
 python/amici/petab_import.py |  2 +-
 4 files changed, 85 insertions(+), 20 deletions(-)

diff --git a/INSTALL.md b/INSTALL.md
index ed49840522..dd802ccf61 100755
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -52,6 +52,20 @@ and custom installations.
 For Python-AMICI usage see 
 [https://github.com/ICB-DCM/AMICI/blob/master/documentation/PYTHON.md](https://github.com/ICB-DCM/AMICI/blob/master/documentation/PYTHON.md).
 
+
+### Installation of development versions
+
+To install development versions which have not been released to pypi yet,
+you can install AMICI with pip directly from GitHub using:
+
+    pip3 install -e https://github.com/icb-dcm/amici/archive/develop.zip#egg=amici\&subdirectory=python/sdist
+
+Replace `develop` by the branch or commit you want to install.
+
+Note that this will probably not work on Windows which does not support
+symlinks by default
+(https://stackoverflow.com/questions/5917249/git-symlinks-in-windows/49913019#49913019).
+
 ### Light installation 
 
 In case you only want to use the AMICI Python package for generating model code
diff --git a/documentation/amici_refs.bib b/documentation/amici_refs.bib
index 8366aecbe1..0f3c404f39 100644
--- a/documentation/amici_refs.bib
+++ b/documentation/amici_refs.bib
@@ -145,21 +145,6 @@ @article{FroehlichKal2017
 	Year = {2017},
 	Bdsk-Url-1 = {https://doi.org/10.1371/journal.pcbi.1005331}}
 
-@article{FroehlichKes2017,
-	Author = {Fr\"ohlich, Fabian and Kessler, Thomas and Weindl, Daniel and Shadrin, Alexey and Schmiester, Leonard and Hache, Hendrik and Muradyan, Artur and Schuette, Moritz and Lim, Ji-Hyun and Heinig, Matthias and Theis, Fabian and Lehrach, Hans and Wierling, Christoph and Lange, Bodo and Hasenauer, Jan},
-	Date-Added = {2019-03-19 23:04:09 +0100},
-	Date-Modified = {2019-03-19 23:04:09 +0100},
-	Doi = {10.1101/174094},
-	Eprint = {http://www.biorxiv.org/content/early/2017/08/09/174094.full.pdf},
-	Journal = {bioRxiv},
-	Keywords = {drug response; cancer; adjoint sensitivity; ccle},
-	Publisher = {Cold Spring Harbor Labs Journals},
-	Title = {Efficient parameterization of large-scale mechanistic models enables drug response prediction for cancer cell lines},
-	Url = {http://www.biorxiv.org/content/early/2017/08/09/174094},
-	Year = {2017},
-	Bdsk-Url-1 = {http://www.biorxiv.org/content/early/2017/08/09/174094},
-	Bdsk-Url-2 = {https://doi.org/10.1101/174094}}
-
 @article{FroehlichRei2018,
 	Abstract = {Single-cell time-lapse studies have advanced the quantitative understanding of cellular pathways and their inherent cell-to-cell variability. However, parameters retrieved from individual experiments are model dependent and their estimation is limited, if based on solely one kind of experiment. Hence, methods to integrate data collected under different conditions are expected to improve model validation and information content. Here we present a multi-experiment nonlinear mixed effect modeling approach for mechanistic pathway models, which allows the integration of multiple single-cell perturbation experiments. We apply this approach to the translation of green fluorescent protein after transfection using a massively parallel read-out of micropatterned single-cell arrays. We demonstrate that the integration of data from perturbation experiments allows the robust reconstruction of cell-to-cell variability, i.e., parameter densities, while each individual experiment provides insufficient information. Indeed, we show that the integration of the datasets on the population level also improves the estimates for individual cells by breaking symmetries, although each of them is only measured in one experiment. Moreover, we confirmed that the suggested approach is robust with respect to batch effects across experimental replicates and can provide mechanistic insights into the nature of batch effects. We anticipate that the proposed multi-experiment nonlinear mixed effect modeling approach will serve as a basis for the analysis of cellular heterogeneity in single-cell dynamics.},
 	Author = {Fr{\"o}hlich, Fabian and Reiser, Anita and Fink, Laura and Wosch{\'e}e, Daniel and Ligon, Thomas and Theis, Fabian Joachim and R{\"a}dler, Joachim Oskar and Hasenauer, Jan},
@@ -478,4 +463,61 @@ @Article{SchmiesterSch2019
   url      = {https://doi.org/10.1093/bioinformatics/btz581},
 }
 
+@Article{SchmiesterWei2019,
+  author       = {Schmiester, Leonard and Weindl, Daniel and Hasenauer, Jan},
+  title        = {Statistical inference of mechanistic models from qualitative data using an efficient optimal scaling approach},
+  journal      = {bioRxiv},
+  year         = {2019},
+  abstract     = {Quantitative dynamical models facilitate the understanding of biological processes and the prediction of their dynamics. These models usually comprise unknown parameters, which have to be inferred from experimental data. For quantitative experimental data, there are several methods and software tools available. However, for qualitative data the available approaches are limited and computationally demanding.Here, we consider the optimal scaling method which has been developed in statistics for categorical data and has been applied to dynamical systems. This approach turns qualitative variables into quantitative ones, accounting for constraints on their relation. We derive a reduced formulation for the optimization problem defining the optimal scaling. The reduced formulation possesses the same optimal points as the established formulation but requires less degrees of freedom. Parameter estimation for dynamical models of cellular pathways revealed that the reduced formulation improves the robustness and convergence of optimizers. This resulted in substantially reduced computation times.We implemented the proposed approach in the open-source Python Parameter EStimation TOolbox (pyPESTO) to facilitate reuse and extension. The proposed approach enables efficient parameterization of quantitative dynamical models using qualitative data.},
+  doi          = {10.1101/848648},
+  elocation-id = {848648},
+  eprint       = {https://www.biorxiv.org/content/early/2019/11/20/848648.full.pdf},
+  publisher    = {Cold Spring Harbor Laboratory},
+  url          = {https://www.biorxiv.org/content/early/2019/11/20/848648},
+}
+
+@Article{PedretscherKal2019,
+  author   = {B. Pedretscher and B. Kaltenbacher and O. Pfeiler},
+  title    = {Parameter identification and uncertainty quantification in stochastic state space models and its application to texture analysis},
+  journal  = {Applied Numerical Mathematics},
+  year     = {2019},
+  volume   = {146},
+  pages    = {38 - 54},
+  issn     = {0168-9274},
+  abstract = {In this paper, a computational framework, which enables efficient and robust parameter identification, as well as uncertainty quantification in state space models based on Itô stochastic processes, is presented. For optimization, a Maximum Likelihood approach based on the system's corresponding Fokker-Planck equation is followed. Gradient information is included by means of an adjoint approach, which is based on the Lagrangian of the optimization problem. To quantify the uncertainty of the Maximum-A-Posteriori estimates of the model parameters, a Bayesian inference approach based on Markov Chain Monte Carlo simulations, as well as profile likelihoods are implemented and compared in terms of runtime and accuracy. The framework is applied to experimental electron backscatter diffraction data of a fatigued metal film, where the aim is to develop a model, which consistently and physically meaningfully captures the metal's microstructural changes that are caused by external loading.},
+  doi      = {https://doi.org/10.1016/j.apnum.2019.06.020},
+  keywords = {Stochastic state space model, Parameter identification, Uncertainty quantification, Profile likelihood, Adjoint approach, Fokker-Planck equation, Thermo-mechanical fatigue, Texture analysis},
+  url      = {http://www.sciencedirect.com/science/article/pii/S0168927419301722},
+}
+
+@Article{FroehlichKes2018,
+  author    = {Fröhlich, Fabian and Kessler, Thomas and Weindl, Daniel and Shadrin, Alexey and Schmiester, Leonard and Hache, Hendrik and Muradyan, Artur and Schütte, Moritz and Lim, Ji-Hyun and Heinig, Matthias and Theis, Fabian J. and Lehrach, Hans and Wierling, Christoph and Lange, Bodo and Hasenauer, Jan},
+  title     = {Efficient Parameter Estimation Enables the Prediction of Drug Response Using a Mechanistic Pan-Cancer Pathway Model},
+  journal   = {Cell Systems},
+  year      = {2018},
+  volume    = {7},
+  number    = {6},
+  pages     = {567--579.e6},
+  issn      = {2405-4712},
+  abstract  = {Mechanistic models are essential to deepen the understanding of complex diseases at the molecular level. Nowadays, high-throughput molecular and phenotypic characterizations are possible, but the integration of such data with prior knowledge on signaling pathways is limited by the availability of scalable computational methods. Here, we present a computational framework for the parameterization of large-scale mechanistic models and its application to the prediction of drug response of cancer cell lines from exome and transcriptome sequencing data. This framework is over 104 times faster than state-of-the-art methods, which enables modeling at previously infeasible scales. By applying the framework to a model describing major cancer-associated pathways (>1,200 species and >2,600 reactions), we could predict the effect of drug combinations from single drug data. This is the first integration of high-throughput datasets using large-scale mechanistic models. We anticipate this to be the starting point for development of more comprehensive models allowing a deeper mechanistic insight.
+Mechanistic models are essential to deepen the understanding of complex diseases at the molecular level. Nowadays, high-throughput molecular and phenotypic characterizations are possible, but the integration of such data with prior knowledge on signaling pathways is limited by the availability of scalable computational methods. Here, we present a computational framework for the parameterization of large-scale mechanistic models and its application to the prediction of drug response of cancer cell lines from exome and transcriptome sequencing data. This framework is over 104 times faster than state-of-the-art methods, which enables modeling at previously infeasible scales. By applying the framework to a model describing major cancer-associated pathways (>1,200 species and >2,600 reactions), we could predict the effect of drug combinations from single drug data. This is the first integration of high-throughput datasets using large-scale mechanistic models. We anticipate this to be the starting point for development of more comprehensive models allowing a deeper mechanistic insight.},
+  comment   = {doi: 10.1016/j.cels.2018.10.013},
+  doi       = {10.1016/j.cels.2018.10.013},
+  publisher = {Elsevier},
+  url       = {https://doi.org/10.1016/j.cels.2018.10.013},
+}
+
+@Article{StaporSch2019,
+  author       = {Stapor, Paul and Schmiester, Leonard and Wierling, Christoph and Lange, Bodo and Weindl, Daniel and Hasenauer, Jan},
+  title        = {Mini-batch optimization enables training of ODE models on large-scale datasets},
+  journal      = {bioRxiv},
+  year         = {2019},
+  abstract     = {Quantitative dynamical models are widely used to study cellular signal processing. A critical step in modeling is the estimation of unknown model parameters from experimental data. As model sizes and datasets are steadily growing, established parameter optimization approaches for mechanistic models become computationally extremely challenging. However, mini-batch optimization methods, as employed in deep learning, have better scaling properties. In this work, we adapt, apply, and benchmark mini-batch optimization for ordinary differential equation (ODE) models thereby establishing a direct link between dynamic modeling and machine learning. On our main application example, a large-scale model of cancer signaling, we benchmark mini-batch optimization against established methods, achieving better optimization results and reducing computation by more than an order of magnitude. We expect that our work will serve as a first step towards mini-batch optimization tailored to ODE models and enable modeling of even larger and more complex systems than what is currently possible.},
+  doi          = {10.1101/859884},
+  elocation-id = {859884},
+  eprint       = {https://www.biorxiv.org/content/early/2019/11/30/859884.full.pdf},
+  publisher    = {Cold Spring Harbor Laboratory},
+  url          = {https://www.biorxiv.org/content/early/2019/11/30/859884},
+}
+
 @Comment{jabref-meta: databaseType:bibtex;}
diff --git a/documentation/references.md b/documentation/references.md
index 759461b437..d4739c7040 100644
--- a/documentation/references.md
+++ b/documentation/references.md
@@ -1,6 +1,6 @@
 # References
 
-List of publications using AMICI. Total number is 32.
+List of publications using AMICI. Total number is 35.
 
 If you applied AMICI in your work and your publication is missing, please let us know via a new Github issue.
 
@@ -24,12 +24,21 @@ If you applied AMICI in your work and your publication is missing, please let us
 <div id="ref-NousiainenInt2019">
 <p>Nousiainen, Kari, Jukka Intosalmi, and Harri Lähdesmäki. 2019. “A Mathematical Model for Enhancer Activation Kinetics During Cell Differentiation.” In <em>Algorithms for Computational Biology</em>, edited by Ian Holmes, Carlos Martín-Vide, and Miguel A. Vega-Rodríguez, 191–202. Cham: Springer International Publishing.</p>
 </div>
+<div id="ref-PedretscherKal2019">
+<p>Pedretscher, B., B. Kaltenbacher, and O. Pfeiler. 2019. “Parameter Identification and Uncertainty Quantification in Stochastic State Space Models and Its Application to Texture Analysis.” <em>Applied Numerical Mathematics</em> 146: 38–54. <a href="https://doi.org/https://doi.org/10.1016/j.apnum.2019.06.020" class="uri">https://doi.org/https://doi.org/10.1016/j.apnum.2019.06.020</a>.</p>
+</div>
 <div id="ref-PittBan2019">
 <p>Pitt, Jake Alan, and Julio R Banga. 2019. “Parameter Estimation in Models of Biological Oscillators: An Automated Regularised Estimation Approach.” <em>BMC Bioinformatics</em> 20 (1): 82. <a href="https://doi.org/10.1186/s12859-019-2630-y" class="uri">https://doi.org/10.1186/s12859-019-2630-y</a>.</p>
 </div>
 <div id="ref-SchmiesterSch2019">
 <p>Schmiester, Leonard, Yannik Schälte, Fabian Fröhlich, Jan Hasenauer, and Daniel Weindl. 2019. “Efficient parameterization of large-scale dynamic models based on relative measurements.” <em>Bioinformatics</em>, July. <a href="https://doi.org/10.1093/bioinformatics/btz581" class="uri">https://doi.org/10.1093/bioinformatics/btz581</a>.</p>
 </div>
+<div id="ref-SchmiesterWei2019">
+<p>Schmiester, Leonard, Daniel Weindl, and Jan Hasenauer. 2019. “Statistical Inference of Mechanistic Models from Qualitative Data Using an Efficient Optimal Scaling Approach.” <em>bioRxiv</em>. Cold Spring Harbor Laboratory. <a href="https://doi.org/10.1101/848648" class="uri">https://doi.org/10.1101/848648</a>.</p>
+</div>
+<div id="ref-StaporSch2019">
+<p>Stapor, Paul, Leonard Schmiester, Christoph Wierling, Bodo Lange, Daniel Weindl, and Jan Hasenauer. 2019. “Mini-Batch Optimization Enables Training of Ode Models on Large-Scale Datasets.” <em>bioRxiv</em>. Cold Spring Harbor Laboratory. <a href="https://doi.org/10.1101/859884" class="uri">https://doi.org/10.1101/859884</a>.</p>
+</div>
 <div id="ref-VillaverdeRai2019">
 <p>Villaverde, Alejandro F., Elba Raimúndez, Jan Hasenauer, and Julio R. Banga. 2019. “A Comparison of Methods for Quantifying Prediction Uncertainty in Systems Biology.” <em>Accepted for Publication in Proc. Of the Foundations of Syst. Biol. In Engin. (FOSBE)</em>.</p>
 </div>
@@ -45,6 +54,9 @@ If you applied AMICI in your work and your publication is missing, please let us
 <div id="ref-BastCal2018">
 <p>Bast, Lisa, Filippo Calzolari, Michael Strasser, Jan Hasenauer, Fabian J. Theis, Jovica Ninkovic, and Carsten Marr. 2018. “Subtle Changes in Clonal Dynamics Underlie the Age-Related Decline in Neurogenesis.” <em>Cell Reports</em>.</p>
 </div>
+<div id="ref-FroehlichKes2018">
+<p>Fröhlich, Fabian, Thomas Kessler, Daniel Weindl, Alexey Shadrin, Leonard Schmiester, Hendrik Hache, Artur Muradyan, et al. 2018. “Efficient Parameter Estimation Enables the Prediction of Drug Response Using a Mechanistic Pan-Cancer Pathway Model.” <em>Cell Systems</em> 7 (6). Elsevier: 567–579.e6. <a href="https://doi.org/10.1016/j.cels.2018.10.013" class="uri">https://doi.org/10.1016/j.cels.2018.10.013</a>.</p>
+</div>
 <div id="ref-FroehlichRei2018">
 <p>Fröhlich, Fabian, Anita Reiser, Laura Fink, Daniel Woschée, Thomas Ligon, Fabian Joachim Theis, Joachim Oskar Rädler, and Jan Hasenauer. 2018. “Multi-Experiment Nonlinear Mixed Effect Modeling of Single-Cell Translation Kinetics After Transfection.” <em>Npj Systems Biology and Applications</em> 5 (1): 1. <a href="https://doi.org/10.1038/s41540-018-0079-7" class="uri">https://doi.org/10.1038/s41540-018-0079-7</a>.</p>
 </div>
@@ -72,9 +84,6 @@ If you applied AMICI in your work and your publication is missing, please let us
 <div id="ref-BallnusHug2017">
 <p>Ballnus, B., S. Hug, K. Hatz, L. Görlitz, J. Hasenauer, and F. J. Theis. 2017. “Comprehensive Benchmarking of Markov Chain Monte Carlo Methods for Dynamical Systems.” <em>BMC Syst. Biol.</em> 11 (63). <a href="https://doi.org/10.1186/s12918-017-0433-1" class="uri">https://doi.org/10.1186/s12918-017-0433-1</a>.</p>
 </div>
-<div id="ref-FroehlichKes2017">
-<p>Fröhlich, Fabian, Thomas Kessler, Daniel Weindl, Alexey Shadrin, Leonard Schmiester, Hendrik Hache, Artur Muradyan, et al. 2017. “Efficient Parameterization of Large-Scale Mechanistic Models Enables Drug Response Prediction for Cancer Cell Lines.” <em>bioRxiv</em>. Cold Spring Harbor Labs Journals. <a href="https://doi.org/10.1101/174094" class="uri">https://doi.org/10.1101/174094</a>.</p>
-</div>
 <div id="ref-FroehlichKal2017">
 <p>Fröhlich, F., B. Kaltenbacher, F. J. Theis, and J. Hasenauer. 2017. “Scalable Parameter Estimation for Genome-Scale Biochemical Reaction Networks.” <em>PLoS Comput. Biol.</em> 13 (1): e1005331. <a href="https://doi.org/10.1371/journal.pcbi.1005331" class="uri">https://doi.org/10.1371/journal.pcbi.1005331</a>.</p>
 </div>
diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py
index 7ecbbbec85..a5c4c51a22 100644
--- a/python/amici/petab_import.py
+++ b/python/amici/petab_import.py
@@ -82,7 +82,7 @@ def get_fixed_parameters(condition_df: pd.DataFrame,
         # check global parameters
         if not sbml_model.getParameter(fixed_parameter) \
                 and not sbml_model.getSpecies(fixed_parameter):
-            logger.log(logging.warning,
+            logger.log(logging.WARN,
                        f"{Fore.YELLOW}Parameter or species '{fixed_parameter}'"
                        " provided in condition table but not present in"
                        " model.")

From 70a515a98fc8c0c421b6b21f3af8b5143a05c153 Mon Sep 17 00:00:00 2001
From: Jan Hasenauer <jan.hasenauer@helmholtz-muenchen.de>
Date: Mon, 2 Dec 2019 19:30:21 +0100
Subject: [PATCH 16/19] Update README.md (#856)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* Update README.md

Co-Authored-By: Yannik Schälte <31767307+yannikschaelte@users.noreply.github.com>
---
 README.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/README.md b/README.md
index efac331da8..3270e6608d 100644
--- a/README.md
+++ b/README.md
@@ -128,6 +128,8 @@ However, the following features are unlikely to be supported:
 - initial assignments for parameters
 - `delay()` due to missing SUNDIALS solver support
 
+In addition to SBML, we also plan to implement support for the [Simulation Experiment Description Markup Language (SED-ML)](https://sed-ml.org/).
+
 
 ## Current build status
 

From 7d73309aca68e275b607e853060fe69cfd72eb34 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Mon, 2 Dec 2019 21:51:24 +0100
Subject: [PATCH 17/19] Speed up sympy simplification (#871)

* Speed up sympy simplification
* Add test case
* Use powsimp instead of simplify
* Adapt walltime reference for model import (Closes #862)
---
 .github/workflows/test-large-model.yml |  2 +-
 python/amici/ode_export.py             | 11 +++++++----
 python/amici/pysb_import.py            |  2 +-
 python/amici/sbml_import.py            |  2 +-
 tests/performance/reference.yml        |  2 +-
 tests/testMisc.py                      | 20 ++++++++++++++++++++
 6 files changed, 31 insertions(+), 8 deletions(-)

diff --git a/.github/workflows/test-large-model.yml b/.github/workflows/test-large-model.yml
index 36d51edc63..07445a43df 100644
--- a/.github/workflows/test-large-model.yml
+++ b/.github/workflows/test-large-model.yml
@@ -4,7 +4,7 @@ on:
     branches:
       - develop
       - master
-      - feature_846_large_model
+      - fix_862_powsimp
   pull_request:
     branches:
       - master
diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py
index 8d63c5ad68..552a43466f 100644
--- a/python/amici/ode_export.py
+++ b/python/amici/ode_export.py
@@ -18,6 +18,7 @@
     pysb = None
 
 
+from typing import Callable, Optional
 from string import Template
 import sympy.printing.ccode as ccode
 from sympy.matrices.immutable import ImmutableDenseMatrix
@@ -714,11 +715,13 @@ class ODEModel:
         derivative from a partial derivative call to enforce a partial
         derivative in the next recursion. prevents infinite recursion
 
-        _simplify: if set to true, symbolic derivative expressions will be
-        simplified. NOTE: This does currently not work with PySB symbols.
+        _simplify: If not None, this function will be used to simplify symbolic
+         derivative expressions. Receives sympy expressions as only argument.
+         To apply multiple simplifications, wrap them in a lambda expression.
+         NOTE: This does currently not work with PySB symbols.
     """
 
-    def __init__(self, simplify: bool = False):
+    def __init__(self, simplify: Optional[Callable] = sp.powsimp):
         """Create a new ODEModel instance.
 
         Arguments:
@@ -1441,7 +1444,7 @@ def _compute_equation(self, name):
                 self._eqs[name] = self._eqs[name].transpose()
 
         if self._simplify:
-            self._eqs[name] = sp.simplify(self._eqs[name])
+            self._eqs[name] = self._simplify(self._eqs[name])
 
     def symNames(self):
         """Returns a list of names of generated symbolic variables
diff --git a/python/amici/pysb_import.py b/python/amici/pysb_import.py
index f7a0e6bd59..579f5c76d0 100644
--- a/python/amici/pysb_import.py
+++ b/python/amici/pysb_import.py
@@ -115,7 +115,7 @@ def ODEModel_from_pysb_importer(model, constants=None,
 
     """
 
-    ODE = ODEModel()
+    ODE = ODEModel(simplify=None)
 
     if not pysb_available:
         raise ImportError(
diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index 67c4fa7aaf..7121e72aaa 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -244,7 +244,7 @@ def sbml2amici(self,
         self.reset_symbols()
         self.processSBML(constantParameters)
         self.processObservables(observables, sigmas, noise_distributions)
-        ode_model = ODEModel(simplify=True)
+        ode_model = ODEModel(simplify=sp.powsimp)
         ode_model.import_from_sbml_importer(self)
         exporter = ODEExporter(
             ode_model,
diff --git a/tests/performance/reference.yml b/tests/performance/reference.yml
index afc02704e7..b392b8bb6f 100644
--- a/tests/performance/reference.yml
+++ b/tests/performance/reference.yml
@@ -1,7 +1,7 @@
 # Reference wall times (seconds) with some buffer
 create_sdist: 25
 install_sdist: 120
-petab_import: 4500 # TODO reduce to 1560 after #862  ^= 26'
+petab_import: 1680 # ^= 28'
 install_model: 400
 forward_simulation: 2
 forward_sensitivities: 6
diff --git a/tests/testMisc.py b/tests/testMisc.py
index 877afceec4..76c7eecc53 100755
--- a/tests/testMisc.py
+++ b/tests/testMisc.py
@@ -153,6 +153,26 @@ def test_constant_species_to_parameters(self):
         assert len(list(r.getListOfProducts())) == 0
         assert len(list(r.getListOfModifiers())) == 0
 
+    def test_hill_function_dwdx(self):
+        """Kinetic laws with Hill functions, may lead to NaNs in the Jacobian
+        if involved states are zero if not properly arranged symbolically.
+        Test that what we are applying the right sympy simplification."""
+
+        w = sp.sympify('Pow(x1, p1) / (Pow(x1, p1) + a)')
+        dwdx = w.diff(sp.Symbol('x1'))
+
+        # Verify that without simplification we fail
+        with self.assertRaises(ZeroDivisionError):
+            with sp.evaluate(False):
+                res = dwdx.subs({'x1': 0.0})
+            _ = str(res)
+
+        # Test that powsimp does the job
+        dwdx = sp.powsimp(dwdx)
+        with sp.evaluate(False):
+            res = dwdx.subs({'x1': 0.0})
+        _ = str(res)
+
 
 if __name__ == '__main__':
     suite = unittest.TestSuite()

From c5d979f56843fca86e13b2afbd68aa80caef5ac7 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <dweindl@users.noreply.github.com>
Date: Mon, 2 Dec 2019 22:21:57 +0100
Subject: [PATCH 18/19] Update .github/workflows/sbml-semantic-test-suite.yml
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-Authored-By: Fabian Fröhlich <fabian@schaluck.com>
---
 .github/workflows/sbml-semantic-test-suite.yml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/.github/workflows/sbml-semantic-test-suite.yml b/.github/workflows/sbml-semantic-test-suite.yml
index e00bc354ec..1cbdf6c34c 100644
--- a/.github/workflows/sbml-semantic-test-suite.yml
+++ b/.github/workflows/sbml-semantic-test-suite.yml
@@ -4,7 +4,6 @@ on:
     branches:
       - develop
       - master
-      - fix_simplify
   pull_request:
     branches:
       - master

From c673dfabba3b0f77398d7dbed932b834f41b9d95 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@helmholtz-muenchen.de>
Date: Mon, 2 Dec 2019 23:57:54 +0100
Subject: [PATCH 19/19] Enable travis for release-* branches

---
 .travis.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.travis.yml b/.travis.yml
index 4754ec533e..36af7b7c0e 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,6 +2,8 @@ branches:
   only:
     - master
     - develop
+    # Release branches, otherwise they cannot be merged into protected master
+    - /^release.*/
     # Version tags
     - /^v\d+\.\d+(\.\d+)?(-\S*)?$/