Add support for exponential Euler solver

odetoolbox/__init__.py

@@ -23,6 +23,7 @@
 import json
 import logging
 import sys
+import numpy as np
 import sympy
 from sympy.core.expr import Expr as SympyExpr
 
@@ -231,64 +232,74 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
             sys.exit(1)
 
     shape_sys = SystemOfShapes.from_shapes(shapes, parameters=parameters)
-    _, node_is_analytically_solvable = _find_analytically_solvable_equations(shape_sys, shapes, parameters=parameters)
 
+    if Config().use_exponential_euler_solver:
+        solvers_json = []
 
-    #
-    #   generate analytical solutions (propagators) where possible
-    #
+        # compute exponential Euler solver
+        solver_json = shape_sys.generate_exponential_euler_solver()
+        solvers_json.append(solver_json)
+        logging.info("Generating propagators for the following symbols: " + ", ".join([str(k) for k in shape_sys.x_]))
 
-    solvers_json = []
-    if disable_analytic_solver:
-        analytic_syms = []
     else:
-        analytic_syms = [node_sym for node_sym, _node_is_analytically_solvable in node_is_analytically_solvable.items() if _node_is_analytically_solvable]
+        _, node_is_analytically_solvable = _find_analytically_solvable_equations(shape_sys, shapes, parameters=parameters)
 
-    analytic_solver_json = None
-    if analytic_syms:
-        logging.info("Generating propagators for the following symbols: " + ", ".join([str(k) for k in analytic_syms]))
-        sub_sys = shape_sys.get_sub_system(analytic_syms)
-        analytic_solver_json = sub_sys.generate_propagator_solver()
-        analytic_solver_json["solver"] = "analytical"
-        solvers_json.append(analytic_solver_json)
 
+        #
+        #   generate analytical solutions (propagators) where possible
+        #
 
-    #
-    #   generate numerical solvers for the remainder
-    #
-
-    if len(analytic_syms) < len(shape_sys.x_):
-        numeric_syms = list(set(shape_sys.x_) - set(analytic_syms))
-        logging.info("Generating numerical solver for the following symbols: " + ", ".join([str(sym) for sym in numeric_syms]))
-        sub_sys = shape_sys.get_sub_system(numeric_syms)
-        solver_json = sub_sys.generate_numeric_solver(state_variables=shape_sys.x_)
-        solver_json["solver"] = "numeric"   # will be appended to if stiffness testing is used
-        if not disable_stiffness_check:
-            if not PYGSL_AVAILABLE:
-                raise Exception("Stiffness test requested, but PyGSL not available")
-
-            logging.info("Performing stiffness test...")
-            kwargs = {}   # type: Dict[str, Any]
-            if "options" in indict.keys() and "random_seed" in indict["options"].keys():
-                random_seed = int(indict["options"]["random_seed"])
-                assert random_seed >= 0, "Random seed needs to be a non-negative integer"
-                kwargs["random_seed"] = random_seed
-            if "parameters" in indict.keys():
-                kwargs["parameters"] = indict["parameters"]
-            if "stimuli" in indict.keys():
-                kwargs["stimuli"] = indict["stimuli"]
-            for key in ["sim_time", "max_step_size", "integration_accuracy_abs", "integration_accuracy_rel"]:
-                if "options" in indict.keys() and key in Config().keys():
-                    kwargs[key] = float(Config()[key])
-            if not analytic_solver_json is None:
-                kwargs["analytic_solver_dict"] = analytic_solver_json
-            tester = StiffnessTester(sub_sys, shapes, **kwargs)
-            solver_type = tester.check_stiffness()
-            if not solver_type is None:
-                solver_json["solver"] += "-" + solver_type
-                logging.info(solver_type + " scheme")
-
-        solvers_json.append(solver_json)
+        solvers_json = []
+        if disable_analytic_solver:
+            analytic_syms = []
+        else:
+            analytic_syms = [node_sym for node_sym, _node_is_analytically_solvable in node_is_analytically_solvable.items() if _node_is_analytically_solvable]
+
+        analytic_solver_json = None
+        if analytic_syms:
+            logging.info("Generating propagators for the following symbols: " + ", ".join([str(k) for k in analytic_syms]))
+            sub_sys = shape_sys.get_sub_system(analytic_syms)
+            analytic_solver_json = sub_sys.generate_propagator_solver()
+            analytic_solver_json["solver"] = "analytical"
+            solvers_json.append(analytic_solver_json)
+
+
+        #
+        #   generate numerical solvers for the remainder
+        #
+
+        if len(analytic_syms) < len(shape_sys.x_):
+            numeric_syms = list(set(shape_sys.x_) - set(analytic_syms))
+            logging.info("Generating numerical solver for the following symbols: " + ", ".join([str(sym) for sym in numeric_syms]))
+            sub_sys = shape_sys.get_sub_system(numeric_syms)
+            solver_json = sub_sys.generate_numeric_solver(state_variables=shape_sys.x_)
+            solver_json["solver"] = "numeric"   # will be appended to if stiffness testing is used
+            if not disable_stiffness_check:
+                if not PYGSL_AVAILABLE:
+                    raise Exception("Stiffness test requested, but PyGSL not available")
+
+                logging.info("Performing stiffness test...")
+                kwargs = {}   # type: Dict[str, Any]
+                if "options" in indict.keys() and "random_seed" in indict["options"].keys():
+                    random_seed = int(indict["options"]["random_seed"])
+                    assert random_seed >= 0, "Random seed needs to be a non-negative integer"
+                    kwargs["random_seed"] = random_seed
+                if "parameters" in indict.keys():
+                    kwargs["parameters"] = indict["parameters"]
+                if "stimuli" in indict.keys():
+                    kwargs["stimuli"] = indict["stimuli"]
+                for key in ["sim_time", "max_step_size", "integration_accuracy_abs", "integration_accuracy_rel"]:
+                    if "options" in indict.keys() and key in Config().keys():
+                        kwargs[key] = float(Config()[key])
+                if not analytic_solver_json is None:
+                    kwargs["analytic_solver_dict"] = analytic_solver_json
+                tester = StiffnessTester(sub_sys, shapes, **kwargs)
+                solver_type = tester.check_stiffness()
+                if not solver_type is None:
+                    solver_json["solver"] += "-" + solver_type
+                    logging.info(solver_type + " scheme")
+
+            solvers_json.append(solver_json)
 
 
     #

odetoolbox/config.py

@@ -35,7 +35,8 @@ class Config:
         "max_step_size": 999.,
         "integration_accuracy_abs": 1E-6,
         "integration_accuracy_rel": 1E-6,
-        "forbidden_names": ["oo", "zoo", "nan", "NaN", "__h"]
+        "forbidden_names": ["oo", "zoo", "nan", "NaN", "__h"],
+        "use_exponential_euler_solver": True
     }
 
     def __getitem__(self, key):

odetoolbox/system_of_shapes.py

@@ -88,14 +88,20 @@ def __init__(self, x: sympy.Matrix, A: sympy.Matrix, b: sympy.Matrix, c: sympy.M
         :param b: Vector containing inhomogeneous part (constant term).
         :param c: Vector containing nonlinear part.
         """
-        logging.debug("Initializing system of shapes with x = " + str(x) + ", A = " + str(A) + ", b = " + str(b) + ", c = " + str(c))
+
+        #logging.debug("Initializing system of shapes with x = " + str(x) + ", A = " + str(A) + ", b = " + str(b) + ", c = " + str(c))
         assert x.shape[0] == A.shape[0] == A.shape[1] == b.shape[0] == c.shape[0]
         self.x_ = x
         self.A_ = A
         self.b_ = b
         self.c_ = c
         self.shapes_ = shapes
 
+        logging.debug("System of equations:")
+        logging.debug("x = " + str(self.x_))
+        logging.debug("A = " + repr(self.A_))
+        logging.debug("b = " + str(self.b_))
+        logging.debug("c = " + str(self.c_))
 
     def get_shape_by_symbol(self, sym: str) -> Optional[Shape]:
         for shape in self.shapes_:
@@ -259,11 +265,11 @@ def generate_propagator_solver(self):
         update_expr = {}    # keys are str(variable symbol), values are str(expressions) that evaluate to the new value of the corresponding key
         for row in range(P_sym.shape[0]):
             # assemble update expression for symbol ``self.x_[row]``
-            if not _is_zero(self.c_[row]):
-                raise PropagatorGenerationException("For symbol " + str(self.x_[row]) + ": nonlinear part should be zero for propagators")
+            # if not _is_zero(self.c_[row]):
+            #     raise PropagatorGenerationException("For symbol " + str(self.x_[row]) + ": nonlinear part should be zero for propagators")
 
-            if not _is_zero(self.b_[row]) and self.shape_order_from_system_matrix(row) > 1:
-                raise PropagatorGenerationException("For symbol " + str(self.x_[row]) + ": higher-order inhomogeneous ODEs are not supported")
+            # if not _is_zero(self.b_[row]) and self.shape_order_from_system_matrix(row) > 1:
+            #     raise PropagatorGenerationException("For symbol " + str(self.x_[row]) + ": higher-order inhomogeneous ODEs are not supported")
 
             update_expr_terms = []
             for col in range(P_sym.shape[1]):
@@ -305,6 +311,7 @@ def generate_propagator_solver(self):
         return solver_dict
 
 
+
     def generate_numeric_solver(self, state_variables=None):
         r"""
         Generate the symbolic expressions for numeric integration state updates; return as JSON.
@@ -319,6 +326,102 @@ def generate_numeric_solver(self, state_variables=None):
 
         return solver_dict
 
+    def generate_update_expr_from_propagator_matrix(self, P, prop_prefix: str):
+        r"""generate symbols for each nonzero entry of the nonlin_P matrix"""
+
+        P_sym = sympy.zeros(*P.shape)   # each entry in the propagator matrix is assigned its own symbol
+        P_expr = {}     # the expression corresponding to each propagator symbol
+        update_expr = {}    # keys are str(variable symbol), values are str(expressions) that evaluate to the new value of the corresponding key
+        for row in range(P_sym.shape[0]):
+            # assemble update expression for symbol ``self.x_[row]``
+            update_expr_terms = []
+            for col in range(P_sym.shape[1]):
+                if not _is_zero(P[row, col]):
+                    sym_str = "__NP__{}__{}".format(str(self.x_[row]), str(self.x_[col]))
+                    P_sym[row, col] = sympy.parsing.sympy_parser.parse_expr(sym_str, global_dict=Shape._sympy_globals)
+                    P_expr[sym_str] = P[row, col]
+                    update_expr_terms.append(sym_str + " * " + str(self.x_[col]))
+
+            # if not _is_zero(self.b_[row]):
+            #     # this is an inhomogeneous ODE
+            #     if _is_zero(self.A_[row, row]):
+            #         # of the form x' = const
+            #         update_expr_terms.append(Config().output_timestep_symbol + " * " + str(self.b_[row]))
+            #     else:
+            #         particular_solution = -self.b_[row] / self.A_[row, row]
+            #         sym_str = "__P__{}__{}".format(str(self.x_[row]), str(self.x_[row]))
+            #         update_expr_terms.append("-" + sym_str + " * " + str(self.x_[row]))    # remove the term (add its inverse) that would have corresponded to a homogeneous solution and that was added in the ``for col...`` loop above
+            #         update_expr_terms.append(sym_str + " * (" + str(self.x_[row]) + " - (" + str(particular_solution) + "))" + " + (" + str(particular_solution) + ")")
+
+            update_expr[str(self.x_[row])] = " + ".join(update_expr_terms)
+            update_expr[str(self.x_[row])] = sympy.parsing.sympy_parser.parse_expr(update_expr[str(self.x_[row])], global_dict=Shape._sympy_globals)
+            if not _is_zero(self.b_[row]):
+                # only simplify in case an inhomogeneous term is present
+                update_expr[str(self.x_[row])] = _custom_simplify_expr(update_expr[str(self.x_[row])])
+            logging.info("update_expr[" + str(self.x_[row]) + "] = " + str(update_expr[str(self.x_[row])]))
+
+        return update_expr
+
+
+
+
+
+    def generate_nonlin_propagators(self, P):
+        #
+        #   generate symbols for each nonzero entry of the propagator matrix
+        #
+
+        P_sym = sympy.zeros(*P.shape)   # each entry in the propagator matrix is assigned its own symbol
+        P_expr = {}     # the expression corresponding to each propagator symbol
+        for row in range(P_sym.shape[0]):
+            for col in range(P_sym.shape[1]):
+                if not _is_zero(P[row, col]):
+                    sym_str = "__NP__{}__{}".format(str(self.x_[row]), str(self.x_[col]))
+                    P_sym[row, col] = sympy.parsing.sympy_parser.parse_expr(sym_str, global_dict=Shape._sympy_globals)
+                    P_expr[sym_str] = P[row, col]
+
+        return P_expr
+
+    def generate_exponential_euler_solver(self):
+        r"""
+        """
+
+        # for the linear part, compute propagators
+        linear_propagator_matrix = self._generate_propagator_matrix(self.A_)
+        linear_solver_json = self.generate_propagator_solver()
+
+        # for the nonlinear part
+        nonlin_P = self.A_.inv() - self.A_.inv() * linear_propagator_matrix
+        nonlin_update_expr = self.generate_update_expr_from_propagator_matrix(nonlin_P, prop_prefix="__NP")
+        nonlin_propagators = self.generate_nonlin_propagators(nonlin_P)
+
+        all_state_symbols = [str(sym) for sym in self.x_]
+        initial_values = {sym: str(self.get_initial_value(sym)) for sym in all_state_symbols}
+
+        combined_update_expr = linear_solver_json["update_expressions"]
+        for sym in nonlin_update_expr.keys():
+            if sym in combined_update_expr.keys():
+                combined_update_expr[sym] += nonlin_update_expr[sym]
+            else:
+                combined_update_expr[sym] = nonlin_update_expr[sym]
+
+        print("linear_solver_json = " + str(linear_solver_json))
+        print("nonlin_update_expr = " + str(nonlin_update_expr))
+
+        solver_dict = {"solver": "analytic",    # actually should be "exponential_euler"; but "analytic" fools NESTML into accepting the update equations and propagators with the existing software infrastructure (no need for a special case)
+                       "update_expressions": combined_update_expr,
+                       "state_variables": all_state_symbols,
+                       "propagators": linear_solver_json["propagators"] | nonlin_propagators,
+                       "initial_values": initial_values }
+        print(solver_dict)
+        # logging.info(json.dumps(solver_dict, indent=4, sort_keys=True))
+        return solver_dict
+
+
+
+
+
+
 
     def reconstitute_expr(self, state_variables=None):
         r"""