diff --git a/doc/index.rst b/doc/index.rst index 761ffd45..e6c2ee80 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -141,9 +141,12 @@ The following flags exist: * - ``disable_stiffness_check`` - :python:`False` - Set to True to disable stiffness check. + * - ``disable_analytic_solver`` + - :python:`False` + - Set to True to return numerical solver recommendations, and no propagators, even for ODEs that are analytically tractable. * - ``disable_singularity_detection`` - :python:`False` - - Set to True to disable detection of conditions under which numerical singularities (division by zero) could occur. + - Set to True to disable detection of conditions under which numerical singularities (division by zero) could occur in the generated analytic solver. This can be useful for analytic solvers containing a large amount of conditions, which could take a long time to compute. If True, at most one analytic solver will be returned, in which numerical singularities could occur. * - ``use_alternative_expM`` - :python:`False` - If :python:`False`, use the sympy function ``sympy.exp`` to compute the matrix exponential. If :python:`True`, use an alternative function (see :py:func:`odetoolbox.sympy_helpers.expMt` for details). This can be useful as calls to ``sympy.exp`` can sometimes take a very large amount of time. @@ -372,6 +375,10 @@ The following global options are defined. Note that all are typically formatted - :python:`["oo", "zoo", "nan", "NaN", "__h"]` - list of strings - For each forbidden name: emit an error if a variable or parameter by this name occurs in the input. + * - ``use_alternative_expM`` + - :python:`False` + - boolean + - If :python:`False`, use the sympy function ``sympy.exp`` to compute the matrix exponential. If :python:`True`, use an alternative function (see :py:func:`odetoolbox.sympy_helpers.expMt` for details). This can be useful as calls to ``sympy.exp`` can sometimes take a very large amount of time. Output diff --git a/odetoolbox/__init__.py b/odetoolbox/__init__.py index 59ebd51f..ec68ef10 100644 --- a/odetoolbox/__init__.py +++ b/odetoolbox/__init__.py @@ -36,8 +36,8 @@ import pygsl.odeiv as odeiv PYGSL_AVAILABLE = True except ImportError as ie: - logging.warning("PyGSL is not available. The stiffness test will be skipped.") - logging.warning("Error when importing: " + str(ie)) + logging.getLogger(__name__).warning("PyGSL is not available. The stiffness test will be skipped.") + logging.getLogger(__name__).warning("Error when importing: " + str(ie)) PYGSL_AVAILABLE = False if PYGSL_AVAILABLE: @@ -45,6 +45,7 @@ try: import graphviz + logging.getLogger("graphviz").setLevel(logging.ERROR) PLOT_DEPENDENCY_GRAPH = True except ImportError: PLOT_DEPENDENCY_GRAPH = False @@ -61,7 +62,7 @@ def _find_analytically_solvable_equations(shape_sys, shapes, parameters=None): Perform dependency analysis and plot dependency graph. """ - logging.info("Finding analytically solvable equations...") + logging.getLogger(__name__).debug("Finding analytically solvable equations...") dependency_edges = shape_sys.get_dependency_edges() if PLOT_DEPENDENCY_GRAPH: @@ -94,7 +95,7 @@ def _read_global_config(indict): r""" Process global configuration options. """ - logging.info("Processing global options...") + logging.getLogger(__name__).debug("Processing global options...") if "options" in indict.keys(): for key, value in indict["options"].items(): assert key in Config.config.keys(), "Unknown key specified in global options dictionary: \"" + str(key) + "\"" @@ -108,12 +109,13 @@ def _from_json_to_shapes(indict, parameters=None) -> Tuple[List[Shape], Dict[sym :param indict: ODE-toolbox input dictionary. """ - logging.info("Processing input shapes...") + logging.getLogger(__name__).debug("Processing input...") # first run for grabbing all the variable names. Coefficients might be incorrect. all_variable_symbols = [] all_parameter_symbols = set() all_variable_symbols_ = set() + for shape_json in indict["dynamics"]: shape = Shape.from_json(shape_json, parameters=parameters) all_variable_symbols.extend(shape.get_state_variables()) @@ -122,7 +124,6 @@ def _from_json_to_shapes(indict, parameters=None) -> Tuple[List[Shape], Dict[sym all_parameter_symbols -= all_variable_symbols_ del all_variable_symbols_ assert all([_is_sympy_type(sym) for sym in all_variable_symbols]) - logging.info("All known variables: " + str(all_variable_symbols) + ", all parameters used in ODEs: " + str(all_parameter_symbols)) # validate input for forbidden names for var in set(all_variable_symbols) | all_parameter_symbols: @@ -137,13 +138,13 @@ def _from_json_to_shapes(indict, parameters=None) -> Tuple[List[Shape], Dict[sym assert isinstance(param, SympyExpr) if not param in parameters.keys(): # this parameter was used in an ODE, but not explicitly numerically specified - logging.info("No numerical value specified for parameter \"" + str(param) + "\"") # INFO level because this is OK! + logging.getLogger(__name__).info("No numerical value specified for parameter \"" + str(param) + "\"") # INFO level because this is OK! parameters[param] = None # second run with the now-known list of variable symbols shapes = [] for shape_json in indict["dynamics"]: - shape = Shape.from_json(shape_json, all_variable_symbols=all_variable_symbols, parameters=parameters, _debug=True) + shape = Shape.from_json(shape_json, all_variable_symbols=all_variable_symbols, parameters=parameters) shapes.append(shape) return shapes, parameters @@ -220,15 +221,16 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so _init_logging(log_level) - logging.info("Analysing input:") - logging.info(json.dumps(indict, indent=4, sort_keys=True)) + logging.getLogger(__name__).info("Analysing input:") + logging.getLogger(__name__).info(json.dumps(indict, indent=4, sort_keys=True)) if "dynamics" not in indict: - logging.info("Warning: empty input (no dynamical equations found); returning empty output") + logging.getLogger(__name__).info("Warning: empty input (no dynamical equations found); returning empty output") return [], SystemOfShapes.from_shapes([]), [] _read_global_config(indict) + # copy parameters from the input and make sure keys are of type sympy.Symbol parameters = None if "parameters" in indict.keys(): @@ -250,17 +252,17 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so for shape in shapes: if not shape.is_homogeneous() and shape.order > 1: - logging.error("For symbol " + str(shape.symbol) + ": higher-order inhomogeneous ODEs are not supported") + logging.getLogger(__name__).error("For symbol " + str(shape.symbol) + ": higher-order inhomogeneous ODEs are not supported") sys.exit(1) shape_sys = SystemOfShapes.from_shapes(shapes, parameters=parameters) _, node_is_analytically_solvable = _find_analytically_solvable_equations(shape_sys, shapes, parameters=parameters) - logging.debug("System of equations:") - logging.debug("x = " + str(shape_sys.x_)) - logging.debug("A = " + repr(shape_sys.A_)) - logging.debug("b = " + str(shape_sys.b_)) - logging.debug("c = " + str(shape_sys.c_)) + logging.getLogger(__name__).info("System of equations (with dx/dt = Ax + b + c):") + logging.getLogger(__name__).info("x = " + str(shape_sys.x_)) + logging.getLogger(__name__).info("A = " + repr(shape_sys.A_)) + logging.getLogger(__name__).info("b = " + str(shape_sys.b_)) + logging.getLogger(__name__).info("c = " + str(shape_sys.c_)) # # generate analytical solutions (propagators) where possible @@ -274,7 +276,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so analytic_solver_json = None if analytic_syms: - logging.info("Generating propagators for the following symbols: " + ", ".join([str(k) for k in analytic_syms])) + logging.getLogger(__name__).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(disable_singularity_detection=disable_singularity_detection, use_alternative_expM=use_alternative_expM) analytic_solver_json["solver"] = "analytical" @@ -286,7 +288,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so 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])) + logging.getLogger(__name__).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 @@ -294,7 +296,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so if not PYGSL_AVAILABLE: raise Exception("Stiffness test requested, but PyGSL not available") - logging.info("Performing stiffness test...") + logging.getLogger(__name__).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"]) @@ -313,7 +315,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so solver_type = tester.check_stiffness() if not solver_type is None: solver_json["solver"] += "-" + solver_type - logging.info(solver_type + " scheme") + logging.getLogger(__name__).info(solver_type + " scheme") solvers_json.append(solver_json) @@ -385,10 +387,10 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so if preserve_expressions and sym in preserve_expressions: if "analytic" in solver_json["solver"]: - logging.warning("Not preserving expression for variable \"" + sym + "\" as it is solved by propagator solver") + logging.getLogger(__name__).warning("Not preserving expression for variable \"" + sym + "\" as it is solved by propagator solver") continue - logging.info("Preserving expression for variable \"" + sym + "\"") + logging.getLogger(__name__).info("Preserving expression for variable \"" + sym + "\"") var_def_str = _find_variable_definition(indict, sym, order=1) assert var_def_str is not None solver_json["update_expressions"][sym] = var_def_str.replace("'", Config().differential_order_symbol) @@ -417,8 +419,8 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so for sym, expr in cond_solver["propagators"].items(): cond_solver["propagators"][sym] = str(expr) - logging.info("In ode-toolbox: returning outdict = ") - logging.info(json.dumps(solvers_json, indent=4, sort_keys=True)) + logging.getLogger(__name__).info("Final output result:") + logging.getLogger(__name__).info(json.dumps(solvers_json, indent=4, sort_keys=True)) return solvers_json, shape_sys, shapes @@ -429,9 +431,10 @@ def _init_logging(log_level: Union[str, int] = logging.WARNING): :param log_level: Sets the logging threshold. Logging messages which are less severe than ``log_level`` will be ignored. Log levels can be provided as an integer or string, for example "INFO" (more messages) or "WARN" (fewer messages). For a list of valid logging levels, see https://docs.python.org/3/library/logging.html#logging-levels """ - fmt = '%(levelname)s:%(message)s' + fmt = "[ODE-toolbox] %(levelname)s:%(message)s" logging.basicConfig(format=fmt) - logging.getLogger().setLevel(log_level) + logger = logging.getLogger(__name__) + logger.setLevel(log_level) def analysis(indict, disable_stiffness_check: bool = False, disable_analytic_solver: bool = False, disable_singularity_detection: bool = False, use_alternative_expM: bool = False, preserve_expressions: Union[bool, Iterable[str]] = False, log_level: Union[str, int] = logging.WARNING) -> List[Dict]: diff --git a/odetoolbox/config.py b/odetoolbox/config.py index 09444a6c..43f3f6a5 100644 --- a/odetoolbox/config.py +++ b/odetoolbox/config.py @@ -36,7 +36,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_alternative_expM": False } def __getitem__(self, key): diff --git a/odetoolbox/dependency_graph_plotter.py b/odetoolbox/dependency_graph_plotter.py index 0256d1c0..546218af 100644 --- a/odetoolbox/dependency_graph_plotter.py +++ b/odetoolbox/dependency_graph_plotter.py @@ -75,5 +75,5 @@ def plot_graph(cls, shapes, dependency_edges, node_is_lin, fn=None): dot.edge(str(e[0]), str(e[1])) if not fn is None: - logging.info("Saving dependency graph plot to " + fn) + logging.getLogger(__name__).debug("Saving dependency graph plot to " + fn) dot.render(fn) diff --git a/odetoolbox/mixed_integrator.py b/odetoolbox/mixed_integrator.py index c5f53492..7fc71b35 100644 --- a/odetoolbox/mixed_integrator.py +++ b/odetoolbox/mixed_integrator.py @@ -41,8 +41,8 @@ import pygsl.odeiv as odeiv PYGSL_AVAILABLE = True except ImportError as ie: - logging.warning("PyGSL is not available. The stiffness test will be skipped.") - logging.warning("Error when importing: " + str(ie)) + logging.getLogger(__name__).warning("PyGSL is not available. The stiffness test will be skipped.") + logging.getLogger(__name__).warning("Error when importing: " + str(ie)) PYGSL_AVAILABLE = False @@ -263,7 +263,7 @@ def integrate_ode(self, initial_values=None, h_min_lower_bound=5E-9, raise_error if h_min < h_min_lower_bound: estr = "Integration step below %.e (s=%.f). Please check your ODE." % (h_min_lower_bound, h_min) - logging.warning(estr) + logging.getLogger(__name__).warning(estr) if raise_errors: raise Exception(estr) @@ -337,7 +337,7 @@ def integrate_ode(self, initial_values=None, h_min_lower_bound=5E-9, raise_error if self._debug_plot_dir: self.integrator_debug_plot(t_log, h_log, y_log, dir=self._debug_plot_dir) - logging.info("For integrator = " + str(self.numeric_integrator) + ": h_min = " + str(h_min) + ", h_avg = " + str(h_avg) + ", runtime = " + str(runtime)) + logging.getLogger(__name__).info("For integrator = " + str(self.numeric_integrator) + ": h_min = " + str(h_min) + ", h_avg = " + str(h_avg) + ", runtime = " + str(runtime)) sym_list = self._system_of_shapes.x_ @@ -460,10 +460,10 @@ def step(self, t, y, params): # return [ float(self._update_expr[str(sym)].evalf(subs=self._locals)) for sym in self._system_of_shapes.x_ ] # non-wrapped version _ret = [self._update_expr_wrapped[str(sym)](*y) for sym in self._system_of_shapes.x_] except Exception as e: - logging.error("E==>", type(e).__name__ + ": " + str(e)) - logging.error(" Local parameters at time of failure:") + logging.getLogger(__name__).error("E==>", type(e).__name__ + ": " + str(e)) + logging.getLogger(__name__).error(" Local parameters at time of failure:") for k, v in self._locals.items(): - logging.error(" ", k, "=", v) + logging.getLogger(__name__).error(" ", k, "=", v) raise return _ret diff --git a/odetoolbox/shapes.py b/odetoolbox/shapes.py index b22d16e3..849e8273 100644 --- a/odetoolbox/shapes.py +++ b/odetoolbox/shapes.py @@ -178,8 +178,6 @@ def __init__(self, symbol: sympy.Symbol, order, initial_values, derivative_facto self.upper_bound = _custom_simplify_expr(self.upper_bound) - logging.debug("Created Shape with symbol " + str(self.symbol) + ", derivative_factors = " + str(self.derivative_factors) + ", inhom_term = " + str(self.inhom_term) + ", nonlin_term = " + str(self.nonlin_term)) - def __str__(self): s = "Shape \"" + str(self.symbol) + "\" of order " + str(self.order) return s @@ -276,7 +274,7 @@ def _parse_defining_expression(cls, s: str) -> Tuple[str, int, str]: return symbol, order, rhs @classmethod - def from_json(cls, indict, all_variable_symbols=None, parameters=None, _debug=False): + def from_json(cls, indict, all_variable_symbols=None, parameters=None): r""" Create a :python:`Shape` instance from an input dictionary. @@ -358,7 +356,6 @@ def reconstitute_expr(self) -> SympyExpr: assert sym.is_real for sym in derivative_symbol.free_symbols: assert sym.is_real - logging.debug("Shape " + str(self.symbol) + ": reconstituting expression " + str(expr)) for sym in expr.free_symbols: assert sym.is_real return expr @@ -384,8 +381,6 @@ def split_lin_inhom_nonlin(expr, x, parameters=None): if parameters is None: parameters = {} - logging.debug("Splitting expression " + str(expr) + " (symbols " + str(x) + ")") - lin_factors = sympy.zeros(len(x), 1) inhom_term = sympy.Float(0) nonlin_term = sympy.Float(0) @@ -410,10 +405,6 @@ def split_lin_inhom_nonlin(expr, x, parameters=None): if not is_lin: nonlin_term += term - logging.debug("\tlinear factors: " + str(lin_factors)) - logging.debug("\tinhomogeneous term: " + str(inhom_term)) - logging.debug("\tnonlinear term: " + str(nonlin_term)) - for sym in nonlin_term.free_symbols: assert sym.is_real @@ -423,7 +414,7 @@ def split_lin_inhom_nonlin(expr, x, parameters=None): return lin_factors, inhom_term, nonlin_term @classmethod - def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_variable_symbols=None, debug=False) -> Shape: + def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_variable_symbols=None) -> Shape: r""" Create a Shape object given a function of time. @@ -452,8 +443,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari # ``derivatives`` is a list of all derivatives of `shape` up to the order we are checking, starting at 0. derivatives = [definition, sympy.diff(definition, sympy.Symbol(Config().input_time_symbol, real=True))] - logging.debug("Processing function-of-time shape \"" + symbol + "\" with defining expression = \"" + str(definition) + "\"") - # # to avoid a division by zero below, we have to find a `t` so that the shape function is not zero at this `t`. # @@ -464,8 +453,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari t_val = t_ break - logging.debug("Found t: " + str(t_val)) - if t_val is None: # @@ -481,8 +468,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari order = 1 - logging.debug("\tFinding ODE of order 1...") - derivative_factors = [(1 / derivatives[0] * derivatives[1]).subs(sympy.Symbol(Config().input_time_symbol, real=True), t_val)] diff_rhs_lhs = derivatives[1] - derivative_factors[0] * derivatives[0] found_ode = _is_zero(diff_rhs_lhs) @@ -494,8 +479,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari while not found_ode and order < max_order: order += 1 - logging.debug("\tFinding ODE of order " + str(order) + "...") - # Add the next higher derivative to the list derivatives.append(sympy.diff(derivatives[-1], sympy.Symbol(Config().input_time_symbol, real=True))) @@ -535,7 +518,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari # diff_rhs_lhs = 0 - logging.debug("\tchecking whether shape definition is satisfied...") for k in range(order): diff_rhs_lhs -= derivative_factors[k] * derivatives[k] diff_rhs_lhs += derivatives[order] @@ -547,7 +529,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari if not found_ode: raise Exception("Shape does not satisfy any ODE of order <= " + str(max_order)) - logging.debug("Shape satisfies ODE of order = " + str(order)) # # calculate the initial values of the found ODE @@ -558,7 +539,7 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari return cls(sympy.Symbol(symbol, real=True), order, initial_values, derivative_factors) @classmethod - def from_ode(cls, symbol: str, definition: str, initial_values: dict, all_variable_symbols=None, lower_bound=None, upper_bound=None, parameters=None, debug=False, **kwargs) -> Shape: + def from_ode(cls, symbol: str, definition: str, initial_values: dict, all_variable_symbols=None, lower_bound=None, upper_bound=None, parameters=None, **kwargs) -> Shape: r""" Create a :python:`Shape` object given an ODE and initial values. @@ -583,8 +564,6 @@ def from_ode(cls, symbol: str, definition: str, initial_values: dict, all_variab for sym in all_variable_symbols: assert sym.is_real - logging.debug("\nProcessing differential-equation form shape " + str(symbol) + " with defining expression = \"" + str(definition) + "\"") - if all_variable_symbols is None: all_variable_symbols = [] @@ -625,6 +604,5 @@ def from_ode(cls, symbol: str, definition: str, initial_values: dict, all_variab sym = sympy.Symbol(symbol, real=True) shape = cls(sym, order, initial_values, local_derivative_factors, inhom_term, nonlin_term, lower_bound, upper_bound) - logging.debug("\tReturning shape: " + str(shape)) return shape diff --git a/odetoolbox/singularity_detection.py b/odetoolbox/singularity_detection.py index e6b7880c..519bd8c8 100644 --- a/odetoolbox/singularity_detection.py +++ b/odetoolbox/singularity_detection.py @@ -168,7 +168,7 @@ def find_inhomogeneous_singularities(A: sympy.Matrix, b: sympy.Matrix) -> Set[Sy conditions a set with equations, where the left-hand side of each equation is the variable that is to be substituted, and the right-hand side is the expression to put in its place """ - logging.debug("Checking for singularities due to inhomogeneous terms in the system of ODEs...") + logging.getLogger(__name__).debug("Checking for singularities (divisions by zero) in the inhomogeneous part of the update equations...") conditions = set() for row in range(A.shape[0]): @@ -182,6 +182,13 @@ def find_inhomogeneous_singularities(A: sympy.Matrix, b: sympy.Matrix) -> Set[Sy conditions = conditions.union(SingularityDetection._find_singularity_conditions_in_expression(particular_solution)) conditions = SingularityDetection._filter_valid_conditions(conditions, A) # filters out the invalid conditions (invalid means those for which A is not defined) + if conditions: + # if there is one or more condition under which the solution goes to infinity... + + logging.getLogger(__name__).warning("Under certain conditions, one or more inhomogeneous term(s) in the system contain a division by zero.") + logging.getLogger(__name__).warning("List of all conditions that result in a division by zero:") + for cond_set in conditions: + logging.getLogger(__name__).warning("\t" + r" ∧ ".join([str(eq.lhs) + " = " + str(eq.rhs) for eq in cond_set])) return conditions @@ -203,7 +210,7 @@ def find_propagator_singularities(P: sympy.Matrix, A: sympy.Matrix) -> Set[Symme conditions a set with equations, where the left-hand side of each equation is the variable that is to be substituted, and the right-hand side is the expression to put in its place """ - logging.debug("Checking for singularities in the propagator matrix...") + logging.getLogger(__name__).debug("Checking for singularities (divisions by zero) in the propagator matrix...") try: conditions = SingularityDetection._generate_singularity_conditions(P) conditions = SingularityDetection._filter_valid_conditions(conditions, A) # filters out the invalid conditions (invalid means those for which A is not defined) diff --git a/odetoolbox/stiffness.py b/odetoolbox/stiffness.py index 764c50b4..a8da0c15 100644 --- a/odetoolbox/stiffness.py +++ b/odetoolbox/stiffness.py @@ -36,8 +36,8 @@ import pygsl.odeiv as odeiv PYGSL_AVAILABLE = True except ImportError as ie: - logging.warning("PyGSL is not available. The stiffness test will be skipped.") - logging.warning("Error when importing: " + str(ie)) + logging.getLogger(__name__).warning("PyGSL is not available. The stiffness test will be skipped.") + logging.getLogger(__name__).warning("Error when importing: " + str(ie)) PYGSL_AVAILABLE = False @@ -110,10 +110,10 @@ def check_stiffness(self, raise_errors=False): step_min_exp, step_average_exp, runtime_exp = self._evaluate_integrator(odeiv.step_rk4, raise_errors=raise_errors) step_min_imp, step_average_imp, runtime_imp = self._evaluate_integrator(odeiv.step_bsimp, raise_errors=raise_errors) except ParametersIncompleteException: - logging.warning("Stiffness test not possible because numerical values were not specified for all parameters.") + logging.getLogger(__name__).warning("Stiffness test not possible because numerical values were not specified for all parameters.") return None - # logging.info("runtime (imp:exp): %f:%f" % (runtime_imp, runtime_exp)) + # logging.getLogger(__name__).info("runtime (imp:exp): %f:%f" % (runtime_imp, runtime_exp)) return self._draw_decision(step_min_imp, step_min_exp, step_average_imp, step_average_exp) @@ -141,7 +141,7 @@ def _evaluate_integrator(self, integrator, h_min_lower_bound=1E-12, raise_errors # initialise and run mixed integrator # - logging.info("Simulating for " + str(self.sim_time) + " with max_step_size = " + str(self.max_step_size)) + logging.getLogger(__name__).info("Simulating for " + str(self.sim_time) + " with max_step_size = " + str(self.max_step_size)) mixed_integrator = MixedIntegrator(integrator, self.system_of_shapes, @@ -157,7 +157,7 @@ def _evaluate_integrator(self, integrator, h_min_lower_bound=1E-12, raise_errors alias_spikes=self.alias_spikes) h_min, h_avg, runtime = (lambda x: x[:3])(mixed_integrator.integrate_ode(h_min_lower_bound=h_min_lower_bound, raise_errors=raise_errors, debug=debug)) - logging.info("For integrator = " + str(integrator) + ": h_min = " + str(h_min) + ", h_avg = " + str(h_avg) + ", runtime = " + str(runtime)) + logging.getLogger(__name__).info("For integrator = " + str(integrator) + ": h_min = " + str(h_min) + ", h_avg = " + str(h_avg) + ", runtime = " + str(runtime)) return h_min, h_avg, runtime diff --git a/odetoolbox/sympy_helpers.py b/odetoolbox/sympy_helpers.py index dd4e1a5d..928b84cc 100644 --- a/odetoolbox/sympy_helpers.py +++ b/odetoolbox/sympy_helpers.py @@ -149,7 +149,7 @@ def _custom_simplify_expr(expr: str): try: # skip simplification for long expressions if len(str(expr)) > Config().expression_simplification_threshold: - logging.warning("Length of expression \"" + str(expr) + "\" exceeds sympy simplification threshold") + logging.getLogger(__name__).warning("Length of expression \"" + str(expr) + "\" exceeds sympy simplification threshold") _simplify_expr = compile(Config().simplify_expression, filename="", mode="eval") expr_simplified = eval(_simplify_expr) diff --git a/odetoolbox/system_of_shapes.py b/odetoolbox/system_of_shapes.py index 1e086f67..903d27f6 100644 --- a/odetoolbox/system_of_shapes.py +++ b/odetoolbox/system_of_shapes.py @@ -47,6 +47,7 @@ def get_block_diagonal_blocks(A): assert A.shape[0] == A.shape[1], "matrix A should be square" A_connectivity_undirected = (A != 0) | (A.T != 0) # make symmetric (undirected) connectivity graph from the system matrix + A_connectivity_undirected = np.triu(A_connectivity_undirected) graph_components = scipy.sparse.csgraph.connected_components(A_connectivity_undirected)[1] @@ -89,7 +90,6 @@ 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)) assert x.shape[0] == A.shape[0] == A.shape[1] == b.shape[0] == c.shape[0] self.x_ = x self.A_ = A @@ -200,7 +200,7 @@ def get_sub_system(self, symbols): return SystemOfShapes(x_sub, A_sub, b_sub, c_sub, shapes_sub) - def _generate_propagator_matrix(self, A, use_alternative_expM: bool = False): + def _generate_propagator_matrix(self, A, use_alternative_expM: bool = False) -> sympy.Matrix: r"""Generate the propagator matrix by matrix exponentiation.""" if use_alternative_expM: @@ -212,6 +212,12 @@ def _generate_propagator_matrix(self, A, use_alternative_expM: bool = False): # optimized: compute propagators separately for each block diagonal element of ``A`` logging.debug("Computing propagator matrix (block-diagonal optimisation)...") blocks = get_block_diagonal_blocks(np.array(A)) + propagators = [sympy.simplify(expM(sympy.Matrix(block) * sympy.Symbol(Config().output_timestep_symbol, real=True))) for block in blocks] + if Config().use_alternative_expM: + expM = expMt + else: + expM = sympy.exp + propagators = [sympy.simplify(expM(sympy.Matrix(block) * sympy.Symbol(Config().output_timestep_symbol, real=True))) for block in blocks] P = sympy.Matrix(scipy.linalg.block_diag(*propagators)) except GetBlockDiagonalException: @@ -258,7 +264,6 @@ def generate_propagator_solver(self, disable_singularity_detection: bool = False r""" Generate the propagator matrix and symbolic expressions for propagator-based updates; return as JSON. """ - P = self._generate_propagator_matrix(self.A_, use_alternative_expM=use_alternative_expM) # @@ -277,7 +282,6 @@ def generate_propagator_solver(self, disable_singularity_detection: bool = False logging.info("List of all conditions that result in a division by zero:") for cond in conditions: logging.info("\t" + str(cond.lhs) + " = " + str(cond.rhs)) - logging.info("Alternate solvers will be generated for each of these conditions (and combinations thereof).") # generate solver for the base case (with singularity conditions that are not met) default_solver = self.generate_solver_dict_based_on_propagator_matrix_(P) @@ -295,6 +299,7 @@ def generate_propagator_solver(self, disable_singularity_detection: bool = False num_conditions = len(conditions) condition_permutations = list(itertools.product([False, True], repeat=num_conditions)) + logging.info("Alternate solvers will be generated for each of these conditions (and combinations thereof), which amounts to " + str(len(condition_permutations)) + " solvers that will be generated.") for condition_permutation in condition_permutations: # each ``condition_permutation[i]`` is True/False corresponding to condition i @@ -328,7 +333,7 @@ def generate_propagator_solver(self, disable_singularity_detection: bool = False conditional_c = conditional_c.subs(eq.lhs, eq.rhs) conditional_dynamics = SystemOfShapes(self.x_, conditional_A, conditional_b, conditional_c, self.shapes_) - solver_dict_conditional = conditional_dynamics.generate_propagator_solver(disable_singularity_detection=True) + solver_dict_conditional = conditional_dynamics.generate_propagator_solver(disable_singularity_detection=True, use_alternative_expM=use_alternative_expM) solver_dict["conditions"][condition_str] = {"propagators": solver_dict_conditional["propagators"], "update_expressions": solver_dict_conditional["update_expressions"]} @@ -336,8 +341,15 @@ def generate_propagator_solver(self, disable_singularity_detection: bool = False return solver_dict + if conditions: + # if there is one or more condition under which the solution goes to infinity... + + logging.getLogger(__name__).warning("Under certain conditions, the propagator matrix is singular (contains infinities).") + logging.getLogger(__name__).warning("List of all conditions that result in a division by zero:") + for cond_set in conditions: + logging.getLogger(__name__).warning("\t" + r" ∧ ".join([str(eq.lhs) + " = " + str(eq.rhs) for eq in cond_set])) except SingularityDetectionException: - logging.warning("Could not check the propagator matrix for singularities.") + logging.getLogger(__name__).warning("Could not check the propagator matrix for singularities.") return self.generate_solver_dict_based_on_propagator_matrix_(P) @@ -389,7 +401,6 @@ def generate_solver_dict_based_on_propagator_matrix_(self, P: sympy.Matrix): 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.debug("update_expr[" + str(self.x_[row]) + "] = " + str(update_expr[str(self.x_[row])])) 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}