Merge pull request #568 from ICB-DCM/develop

0.9.2 Release

include/amici/defines.h

@@ -19,6 +19,7 @@ constexpr double pi = 3.14159265358979323846;
 
 /* Return codes */
 #define AMICI_RECOVERABLE_ERROR          1
+#define AMICI_UNRECOVERABLE_ERROR      -10
 #define AMICI_TOO_MUCH_WORK             -1
 #define AMICI_TOO_MUCH_ACC              -2
 #define AMICI_ERR_FAILURE               -3

matlab/@amimodel/compileAndLinkModel.m

@@ -183,7 +183,7 @@ function compileAndLinkModel(modelname, modelSourceFolder, coptim, debug, funs,
         'interface_matlab', 'misc', ...
         'solver', 'solver_cvodes', 'solver_idas', ...
         'model', 'model_ode', 'model_dae', 'returndata_matlab', ...
-        'forwardproblem', 'steadystateproblem', 'backwardproblem', 'newton_solver',
+        'forwardproblem', 'steadystateproblem', 'backwardproblem', 'newton_solver', ...
         'abstract_model', 'sundials_matrix_wrapper'
     };
     % to be safe, recompile everything if headers have changed. otherwise

matlab/examples/example_steadystate/example_steadystate.m

@@ -184,10 +184,11 @@
         solss = simulate_model_steadystate(tt,log10(p),k,[],options);
         ssxdot(it,:) = solss.diagnosis.xdot;
     end
-    
+
     % Compute steady state wihtout integration before
     sol = simulate_model_steadystate(inf,log10(p),k,[],options);
 
+
     % Test recapturing in the case of Newton solver failing
     options.newton_maxsteps = 4;
     options.maxsteps = 300;
@@ -218,7 +219,7 @@
 
         subplot(1,3,3);
         hold on;
-        bar(sol_newton_fail.diagnosis.newton_numsteps);
+        bar(sol_newton_fail.diagnosis.newton_numsteps([1, 3]));
         legend boxoff;
         title('Number of Newton steps');
         xlabel('Solver run');

models/model_steadystate/model_steadystate.h

@@ -1,6 +1,6 @@
 #ifndef _amici_model_steadystate_h
 #define _amici_model_steadystate_h
-/* Generated by amiwrap (R2017b) 0ba9d6e4dbd5671e7135302f68426a2b8aa5489e */
+/* Generated by amiwrap (R2017b) eb9573b9ed3211a78dec051b0322af263ee3a520 */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -64,7 +64,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 "0ba9d6e4dbd5671e7135302f68426a2b8aa5489e"; };
+    const  std::string getAmiciCommit() const override { return "eb9573b9ed3211a78dec051b0322af263ee3a520"; };
 
     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);

python/amici/ode_export.py

@@ -895,8 +895,13 @@ def import_from_sbml_importer(self, si):
             [sp.Symbol(f'flux_r{idx}') for idx in range(len(si.fluxVector))]
         )
         self._eqs['dxdotdx'] = sp.zeros(si.stoichiometricMatrix.shape[0])
-        symbols['species']['dt'] = \
-            si.stoichiometricMatrix * self.sym('w')
+        if len(si.stoichiometricMatrix):
+            symbols['species']['dt'] = \
+                si.stoichiometricMatrix * self.sym('w')
+        else:
+            symbols['species']['dt'] = sp.zeros(
+                *symbols['species']['identifier'].shape
+            )
 
         for symbol in [s for s in symbols if s != 'my']:
             # transform dict of lists into a list of dicts
@@ -1395,21 +1400,31 @@ def _compute_equation(self, name):
             # if x0_fixedParameters>0 else 0
             # sx0_fixedParameters = sx+deltasx =
             # dx0_fixedParametersdx*sx+dx0_fixedParametersdp
-            self._eqs[name] = \
-                self.eq('x0_fixedParameters').jacobian(self.sym('p'))
+            if len(self.sym('p')):
+                self._eqs[name] = \
+                    self.eq('x0_fixedParameters').jacobian(self.sym('p'))
+            else:
+                self._eqs[name] = sp.zeros(
+                    len(self.eq('x0_fixedParameters')),
+                    len(self.sym('p'))
+                )
 
-            dx0_fixedParametersdx = \
-                self.eq('x0_fixedParameters').jacobian(self.sym('x'))
+            if len(self.sym('x')):
+                dx0_fixedParametersdx = \
+                    self.eq('x0_fixedParameters').jacobian(self.sym('x'))
+            else:
+                dx0_fixedParametersdx = sp.zeros(
+                    len(self.eq('x0_fixedParameters')),
+                    len(self.sym('x'))
+                )
 
             if dx0_fixedParametersdx.is_zero is not True:
                 for ip in range(self._eqs[name].shape[1]):
                     self._eqs[name][:,ip] += \
                         dx0_fixedParametersdx \
                         * self.sym('sx0') \
 
-            for index, formula in enumerate(
-                    self.eq('x0_fixedParameters')
-            ):
+            for index, formula in enumerate(self.eq('x0_fixedParameters')):
                 if formula == 0 or formula == 0.0:
                     self._eqs[name][index, :] = \
                         sp.zeros(1, self._eqs[name].shape[1])
@@ -1639,7 +1654,11 @@ def _multiplication(self, name, x, y,
         else:
             xx = variables[x]
 
-        self._eqs[name] = sign * xx * variables[y]
+        if xx.is_zero is not True and variables[y].is_zero is not True \
+                and len(xx) and len(variables[y]):
+            self._eqs[name] = sign * xx * variables[y]
+        else:
+            self._eqs[name] = sp.zeros(len(xx), len(variables[y]))
 
     def _equationFromComponent(self, name, component):
         """Generates the formulas of a symbolic variable from the attributes
@@ -1715,11 +1734,6 @@ def _generateName(self, name):
             component = self._variable_prototype[name]
         elif name in self._equation_prototype:
             component = self._equation_prototype[name]
-        elif name == 'x':
-            self._names[name]= [
-                s._name for s in self._states if s.conservation_law is None
-            ]
-            return
         else:
             raise Exception(f'No names for {name}')
 
@@ -2213,15 +2227,15 @@ def _writeModelHeader(self):
             'PARAMETER_NAMES_INITIALIZER_LIST':
                 self._getSymbolNameInitializerList('p'),
             'STATE_NAMES_INITIALIZER_LIST':
-                self._getSymbolNameInitializerList('x'),
+                self._getSymbolNameInitializerList('x_rdata'),
             'FIXED_PARAMETER_NAMES_INITIALIZER_LIST':
                 self._getSymbolNameInitializerList('k'),
             'OBSERVABLE_NAMES_INITIALIZER_LIST':
                 self._getSymbolNameInitializerList('y'),
             'PARAMETER_IDS_INITIALIZER_LIST':
                 self._getSymbolIDInitializerList('p'),
             'STATE_IDS_INITIALIZER_LIST':
-                self._getSymbolIDInitializerList('x'),
+                self._getSymbolIDInitializerList('x_rdata'),
             'FIXED_PARAMETER_IDS_INITIALIZER_LIST':
                 self._getSymbolIDInitializerList('k'),
             'OBSERVABLE_IDS_INITIALIZER_LIST':

python/amici/sbml_import.py

@@ -6,6 +6,8 @@
 import libsbml as sbml
 import re
 import math
+from sympy.logic.boolalg import BooleanTrue as spTrue
+from sympy.logic.boolalg import BooleanFalse as spFalse
 
 from .ode_export import ODEExporter, ODEModel
 
@@ -343,23 +345,39 @@ def getSpeciesInitial(index, conc):
                 index = species_ids.index(
                         initial_assignment.getId()
                     )
-                speciesInitial[index] = sp.sympify(
+                symMath = sp.sympify(
                     sbml.formulaToL3String(initial_assignment.getMath())
                 )
-                _check_unsupported_functions(speciesInitial[index])
+                if symMath is not None:
+                    _check_unsupported_functions(symMath, 'InitialAssignment')
+                    speciesInitial[index] = symMath
 
         # flatten initSpecies
         while any([species in speciesInitial.free_symbols
                    for species in self.symbols['species']['identifier']]):
-            speciesInitial = speciesInitial.subs(
-                self.symbols['species']['identifier'],
-                speciesInitial
-            )
+            identical_assignment = next((
+                symbol == init
+                for symbol, init in zip(
+                    self.symbols['species']['identifier'], speciesInitial
+                )
+            ), None)
+            if identical_assignment is not None:
+                raise SBMLException('Species without initial assignment are '
+                                    'currently not supported (this is error '
+                                    'is likely to be due to the existence of '
+                                    'a species assignment rule, which is '
+                                    'also not supported)!')
+            speciesInitial = speciesInitial.subs([
+                (symbol, init)
+                for symbol, init in zip(
+                    self.symbols['species']['identifier'], speciesInitial
+                )
+            ])
 
         self.symbols['species']['value'] = speciesInitial
 
         if self.sbml.isSetConversionFactor():
-            conversion_factor = self.sbml.getConversionFactor()
+            conversion_factor = sp.Symbol(self.sbml.getConversionFactor())
         else:
             conversion_factor = 1.0
 
@@ -481,7 +499,7 @@ def processCompartments(self):
                     sbml.formulaToL3String(initial_assignment.getMath())
                 )
 
-        for comp, vol in zip(self.compartmentSymbols,self.compartmentVolume):
+        for comp, vol in zip(self.compartmentSymbols, self.compartmentVolume):
             self.replaceInAllExpressions(
                comp, vol
             )
@@ -519,29 +537,34 @@ def getElementFromAssignment(element_id):
             assignment = self.sbml.getInitialAssignment(
                 element_id
             )
-            sym = sp.sympify(assignment.getFormula())
+            sym = sp.sympify(sbml.formulaToL3String(assignment.getMath()))
             # this is an initial assignment so we need to use
             # initial conditions
-            sym = sym.subs(
-                self.symbols['species']['identifier'],
-                self.symbols['species']['value']
-            )
+            if sym is not None:
+                sym = sym.subs(
+                    self.symbols['species']['identifier'],
+                    self.symbols['species']['value']
+                )
             return sym
 
         def getElementStoichiometry(element):
             if element.isSetId():
                 if element.getId() in assignment_ids:
-                    return getElementFromAssignment(element.getId())
+                    symMath = getElementFromAssignment(element.getId())
+                    if symMath is None:
+                        symMath = sp.sympify(element.getStoichiometry())
                 elif element.getId() in rulevars:
-                    return sp.sympify(element.getId())
+                    return sp.Symbol(element.getId())
                 else:
                     # dont put the symbol if it wont get replaced by a
                     # rule
-                    return sp.sympify(element.getStoichiometry())
+                    symMath = sp.sympify(element.getStoichiometry())
             elif element.isSetStoichiometry():
-                return sp.sympify(element.getStoichiometry())
+                symMath = sp.sympify(element.getStoichiometry())
             else:
                 return sp.sympify(1.0)
+            _check_unsupported_functions(symMath, 'Stoichiometry')
+            return symMath
 
         def isConstant(specie):
             return specie in self.constantSpecies or \
@@ -578,7 +601,7 @@ def isConstant(specie):
             except:
                 raise SBMLException(f'Kinetic law "{math}" contains an '
                                     'unsupported expression!')
-            _check_unsupported_functions(symMath)
+            _check_unsupported_functions(symMath, 'KineticLaw')
             for r in reactions:
                 elements = list(r.getListOfReactants()) \
                            + list(r.getListOfProducts())
@@ -629,7 +652,7 @@ def processRules(self):
             variable = sp.sympify(rule.getVariable())
             # avoid incorrect parsing of pow(x, -1) in symengine
             formula = sp.sympify(sbml.formulaToL3String(rule.getMath()))
-            _check_unsupported_functions(formula)
+            _check_unsupported_functions(formula, 'Rule')
 
             if variable in stoichvars:
                 self.stoichiometricMatrix = \
@@ -1112,35 +1135,53 @@ def checkLibSBMLErrors(sbml_doc, show_warnings=False):
         )
 
 
-def _check_unsupported_functions(sym):
+def _check_unsupported_functions(sym, expression_type, full_sym=None):
     """Recursively checks the symbolic expression for unsupported symbolic
     functions
 
         Arguments:
             sym: symbolic expressions @type sympy.Basic
+            expression_type: type of expression
 
         Returns:
 
         Raises:
             raises SBMLException if an unsupported function is encountered
     """
+    if full_sym is None:
+        full_sym = sym
 
     unsupported_functions = [
-        sp.functions.factorial, sp.functions.ceiling,
-        sp.functions.Piecewise,
+        sp.functions.factorial, sp.functions.ceiling, sp.functions.floor,
+        sp.functions.Piecewise, spTrue, spFalse, sp.function.UndefinedFunction
     ]
 
-    for arg in sym._args:
-        unsupp_fun = next(
+    unsupp_fun_type = next(
+        (
+            fun_type
+            for fun_type in unsupported_functions
+            if isinstance(sym.func, fun_type)
+        ),
+        None
+    )
+    if unsupp_fun_type:
+        raise SBMLException(f'Encountered unsupported expression '
+                            f'"{sym.func}" of type '
+                            f'"{unsupp_fun_type}" as part of a '
+                            f'{expression_type}: "{full_sym}"!')
+    for fun in list(sym._args) + [sym]:
+        unsupp_fun_type = next(
             (
-                fun
-                for fun in unsupported_functions
-                if isinstance(arg, fun)
+                fun_type
+                for fun_type in unsupported_functions
+                if isinstance(fun, fun_type)
             ),
             None
         )
-        if unsupp_fun:
-            raise SBMLException(f'Encountered unsupported function type '
-                                f'"{unsupp_fun}" in subexpression "{sym}"')
-
-        _check_unsupported_functions(arg)
+        if unsupp_fun_type:
+            raise SBMLException(f'Encountered unsupported expression '
+                                f'"{fun}" of type '
+                                f'"{unsupp_fun_type}" as part of a '
+                                f'{expression_type}: "{full_sym}"!')
+        if fun is not sym:
+            _check_unsupported_functions(fun, expression_type)

scripts/deployPyPi.sh

@@ -7,5 +7,7 @@ set -e
 SCRIPT_PATH=$(dirname $BASH_SOURCE)
 AMICI_PATH=$(cd $SCRIPT_PATH/.. && pwd)
 
+pip3 install twine
+
 # authentication via env variables set in travis
 twine upload $(ls -t ${AMICI_PATH}/build/python/amici-*.tar.gz | head -1)

src/misc.cpp

@@ -7,21 +7,26 @@
 
 namespace amici {
 
-int checkFinite(const int N,const realtype *array, const char* fun){
-    for(int idx = 0; idx < N; idx++) {
-        if(isNaN(array[idx])) {
-            warnMsgIdAndTxt("AMICI:mex:NaN","AMICI encountered a NaN value at index %i of %i in %s! Trying to recover ... ",idx,N,fun);
-            return(AMICI_RECOVERABLE_ERROR);
+int checkFinite(const int N, const realtype *array, const char *fun) {
+    for (int idx = 0; idx < N; idx++) {
+        if (isNaN(array[idx])) {
+            warnMsgIdAndTxt(
+                "AMICI:NaN",
+                "AMICI encountered a NaN value at index %i of %i in %s!", idx,
+                N, fun);
+            return AMICI_RECOVERABLE_ERROR;
         }
-        if(isInf(array[idx])) {
-            warnMsgIdAndTxt("AMICI:mex:Inf","AMICI encountered an Inf value at index %i of %i in %s! Trying to recover ... ",idx,N,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,
+                N, fun);
+            return AMICI_RECOVERABLE_ERROR;
         }
     }
-    return(AMICI_SUCCESS);
+    return AMICI_SUCCESS;
 }
 
-
 double getUnscaledParameter(double scaledParameter, ParameterScaling scaling)
 {
     switch (scaling) {