Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 1 addition & 13 deletions odetoolbox/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,17 +32,6 @@
from .shapes import MalformedInputException, Shape


try:
import pygsl.odeiv as odeiv
PYGSL_AVAILABLE = True
except ImportError as 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:
from .stiffness import StiffnessTester

try:
logging.getLogger("graphviz").setLevel(logging.ERROR)
import graphviz
Expand Down Expand Up @@ -292,8 +281,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
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")
from .stiffness import StiffnessTester

logging.getLogger(__name__).info("Performing stiffness test...")
kwargs = {} # type: Dict[str, Any]
Expand Down
31 changes: 12 additions & 19 deletions odetoolbox/stiffness.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,6 @@

import logging
import numpy as np
import numpy.random
import sympy

from odetoolbox.sympy_helpers import _sympy_parse_real

Expand All @@ -32,15 +30,6 @@
from .spike_generator import SpikeGenerator


try:
import pygsl.odeiv as odeiv
PYGSL_AVAILABLE = True
except ImportError as 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


class StiffnessTester:

def __init__(self, system_of_shapes, shapes, analytic_solver_dict=None, parameters=None, stimuli=None, random_seed=123, max_step_size=np.inf, integration_accuracy_abs=1E-6, integration_accuracy_rel=1E-6, sim_time=100., alias_spikes=False):
Expand Down Expand Up @@ -95,23 +84,29 @@ def random_seed(self, value):
assert value >= 0
self._random_seed = value

def check_stiffness(self, raise_errors=False):
def check_stiffness(self, raise_errors=False) -> str:
r"""
Perform stiffness testing: use implicit and explicit solvers to simulate the dynamical system, then decide which is the better solver to use.

For details, see https://ode-toolbox.readthedocs.io/en/latest/index.html#numeric-solver-selection-criteria

:return: Either :python:`"implicit"`, :python:`"explicit"` or :python:`"warning"`.
:rtype: str
"""
assert PYGSL_AVAILABLE

try:
import pygsl.odeiv as odeiv
except ImportError as ie:
error_msg = "Stiffness test requested, but PyGSL is not available"
logging.getLogger(__name__).error(error_msg)
raise Exception(error_msg)

try:
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.getLogger(__name__).warning("Stiffness test not possible because numerical values were not specified for all parameters.")
return None
error_msg = "Stiffness test not possible because numerical values were not specified for all parameters."
logging.getLogger(__name__).error(error_msg)
raise Exception(error_msg)

# logging.getLogger(__name__).info("runtime (imp:exp): %f:%f" % (runtime_imp, runtime_exp))

Expand All @@ -131,8 +126,6 @@ def _evaluate_integrator(self, integrator, h_min_lower_bound=1E-12, raise_errors
:return h_avg: Average recommended step size.
:return runtime: Wall clock runtime.
"""
assert PYGSL_AVAILABLE

np.random.seed(self.random_seed)

spike_times = SpikeGenerator.spike_times_from_json(self._stimuli, self.sim_time)
Expand Down Expand Up @@ -161,7 +154,7 @@ def _evaluate_integrator(self, integrator, h_min_lower_bound=1E-12, raise_errors

return h_min, h_avg, runtime

def _draw_decision(self, step_min_imp, step_min_exp, step_average_imp, step_average_exp, machine_precision_dist_ratio=10, avg_step_size_ratio=6):
def _draw_decision(self, step_min_imp, step_min_exp, step_average_imp, step_average_exp, machine_precision_dist_ratio=10, avg_step_size_ratio=6) -> str:
r"""
Decide which is the best integrator to use.

Expand Down
Loading