ENH: Add some useful SUNDIALS user functions

MANIFEST.in

@@ -2,7 +2,7 @@ include common.py setup_build.py
 recursive-include scikits *.pyx *.pxd *.pyf *.f
 
 include CONTRIBUTING.md README.md LICENSE.txt
-include tox.ini
+include tox.ini pytest.ini
 include pyproject.toml
 
 prune ci_support

docs/examples/ode/freefall.py

@@ -11,6 +11,7 @@
 from __future__ import print_function
 import numpy as np
 from scikits.odes import ode
+from scikits.odes.sundials import onroot_stop
 
 #data
 g  = 9.81    # gravitational constant
@@ -31,7 +32,7 @@
 # On the other hand experiments 1 and 3 don't use the 'onroot',
 # experiments 2 and 4 do and compute until the time t_end is reached
 # (function onroot_va()). Experiment 5 stops after the first interruption
-# (function onroot_vb()) occurs, whereas experiment 6 stops after the
+# (function onroot_stop()) occurs, whereas experiment 6 stops after the
 # first interruption  at time t>28 (s) (function onroot_vc()).
 # Otherwise all experiments are the same.
 
@@ -66,12 +67,6 @@ def onroot_va(t, y, solver):
 
     return 0
 
-def onroot_vb(t, y, solver):
-    """
-    onroot function to stop solver when root is found
-    """
-    return 1
-
 def onroot_vc(t, y, solver):
     """
     onroot function to reset the solver back at the start, but keep the current
@@ -157,10 +152,10 @@ def print_results(experiment_no, result, require_no_roots=False):
 print_results(4, solver.solve(tspan, y0))
 
 #
-# 5. Solve the problem with onroot function onroot_vb, which behaves similarly
+# 5. Solve the problem with onroot function onroot_stop, which behaves similarly
 # to the default, which is to compute until a root is found.
 #
-solver = ode('cvode', rhs_fn, nr_rootfns=1, rootfn=root_fn, onroot=onroot_vb, old_api=False)
+solver = ode('cvode', rhs_fn, nr_rootfns=1, rootfn=root_fn, onroot=onroot_stop, old_api=False)
 print_results(5, solver.solve(tspan, y0))
 
 #

docs/examples/ode/freefall_tstop.py

@@ -11,6 +11,7 @@
 from __future__ import print_function
 import numpy as np
 from scikits.odes import ode
+from scikits.odes.sundials import ontstop_stop
 
 #data
 g  = 9.81    # gravitational constant
@@ -32,7 +33,7 @@
 # On the other hand experiments 1 and 3 don't use the 'ontstop',
 # experiments 2 and 4 do and compute until the time t_end is reached
 # (function ontstop_va()). Experiment 5 stops after the first interruption
-# (function ontstop_vb()) occurs, whereas experiment 6 stops after the
+# (function ontstop_stop()) occurs, whereas experiment 6 stops after the
 # first interruption  at time t>28 (s) (function ontstop_vc()).
 # Otherwise all experiments are the same.
 
@@ -54,12 +55,6 @@ def ontstop_va(t, y, solver):
 
     return 0
 
-def ontstop_vb(t, y, solver):
-    """
-    ontstop function to stop solver when tstop is found
-    """
-    return 1
-
 def ontstop_vc(t, y, solver):
     """
     ontstop function to reset the solver back at the start, but keep the current
@@ -155,11 +150,11 @@ def print_results(experiment_no, result, require_no_tstop=False):
 print_results(4, solver.solve(tspan, y0))
 
 #
-# 5. Solve the problem with ontstop function ontstop_vb, which behaves similarly
+# 5. Solve the problem with ontstop function ontstop_stop, which behaves similarly
 # to the default, which is to compute until a root is found.
 #
 n = 0
-solver = ode('cvode', rhs_fn, tstop=Y1, ontstop=ontstop_vb, old_api=False)
+solver = ode('cvode', rhs_fn, tstop=Y1, ontstop=ontstop_stop, old_api=False)
 print_results(5, solver.solve(tspan, y0))
 
 #

docs/solvers.rst

@@ -37,3 +37,17 @@ A comparison of different methods is given in following image. In this BDF, RK23
 .. image:: ../ipython_examples/PerformanceTests.png
 
 You can generate above graph via the `Performance notebook <https://github.com/bmcage/odes/blob/master/ipython_examples/Performance%20tests.ipynb>`_.
+
+Solver Specific Options
+#######################
+
+The high level interfaces allow the option of passing solver specific options to
+the solvers. These options are covered in more detail in the `API docs <https://bmcage.github.io/odes>`_, but some features specific to ``odes`` are mentioned below.
+
+SUNDIALS
+========
+
+There are a number of SUNDIALS specific utilities in :py:mod:`scikits.odes.sundials`.
+Firstly there are :py:func:`scikits.odes.sundials.ontstop_stop`, :py:func:`scikits.odes.sundials.ontstop_continue`, :py:func:`scikits.odes.sundials.onroot_stop` and :py:func:`scikits.odes.sundials.onroot_continue`, which can be used with the `ontstop` or `onroot` options to either stop or continue evaluation when tstop or a root is encountered.
+Secondly, there are functions which can be passed to the `err_handler` option to either stop all messages from SUNDIALS being printed (:py:func:`scikits.odes.sundials.drop_all_error_handler`), or to pass them to Python's logging machinery (:py:func:`scikits.odes.sundials.log_error_handler`).
+Finally, the module contains the exceptions which can be caught in user code raised when using the `validate_flags` option.

pytest.ini

@@ -0,0 +1,5 @@
+[pytest]
+log_cli = 1
+log_cli_level = INFO
+log_cli_format = %(asctime)s %(levelname)s %(message)s
+log_cli_date_format = %H:%M:%S

scikits/odes/sundials/__init__.py

@@ -3,6 +3,11 @@
 #
 
 import inspect
+from logging import getLogger
+
+logger = getLogger(__name__)
+DEFAULT_LOG_FORMAT = "SUNDIALS message in %s:%s: %s"
+
 
 class CVODESolveException(Exception):
     """Base class for exceptions raised by ``CVODE.validate_flags``."""
@@ -64,3 +69,86 @@ def _get_num_args(func):
         return numargs
     else:
         return len(inspect.getargspec(func).args)
+
+
+def drop_all_error_handler(error_code, module, func, msg, user_data):
+    """
+    Drop all CVODE/IDA messages, rather than printing them.
+
+    Examples
+    --------
+    >>> scikits.odes.ode('cvode', rhsfuc, err_handler=drop_all_error_handler)
+    """
+    # pylint: disable=unused-argument
+    pass
+
+
+def log_error_handler(error_code, module, func, msg, user_data):
+    """
+    Log all CVODE/IDA messages using the builtin python logging.
+
+    Examples
+    --------
+    >>> scikits.odes.ode('cvode', rhsfuc, err_handler=log_error_handler)
+    """
+    # pylint: disable=unused-argument
+    if error_code > 0:
+        logger.warning(DEFAULT_LOG_FORMAT, module, func, msg)
+    else:
+        logger.error(DEFAULT_LOG_FORMAT, module, func, msg)
+
+
+def onroot_continue(*args):
+    """
+    Always continue after finding root.
+
+    Examples
+    --------
+    >>> scikits.odes.ode(
+    ...     'cvode', rhsfuc, rootfn=rootfn, nr_rootfns=nroots,
+    ...     onroot=onroot_continue
+    ... )
+    """
+    # pylint: disable=unused-argument
+    return 0
+
+
+def onroot_stop(*args):
+    """
+    Always stop after finding root.
+
+    Examples
+    --------
+    >>> scikits.odes.ode(
+    ...     'cvode', rhsfuc, rootfn=rootfn, nr_rootfns=nroots,
+    ...     onroot=onroot_stop
+    ... )
+    """
+    # pylint: disable=unused-argument
+    return 1
+
+
+def ontstop_continue(*args):
+    """
+    Always continue after finding tstop.
+
+    Examples
+    --------
+    >>> scikits.odes.ode(
+    ...     'cvode', rhsfuc, tstop=tstop, ontstop=ontstop_continue
+    ... )
+    """
+    # pylint: disable=unused-argument
+    return 0
+
+
+def ontstop_stop(*args):
+    """
+    Always stop after finding tstop.
+
+    Examples
+    --------
+    >>> scikits.odes.ode('cvode', rhsfuc, tstop=tstop, ontstop=ontstop_stop)
+    """
+    # pylint: disable=unused-argument
+    return 1

scikits/odes/tests/test_dae.py

@@ -11,8 +11,8 @@
 
 from numpy.testing import TestCase, run_module_suite
 from scipy.integrate import ode as Iode
-from scikits.odes import ode,dae
-from scikits.odes.sundials.common_defs import DTYPE
+from .. import ode, dae
+from ..sundials.common_defs import DTYPE
 
 class TestDae(TestCase):
     """

scikits/odes/tests/test_dop.py

@@ -12,8 +12,8 @@
 from numpy.testing import (
     assert_, TestCase, run_module_suite, assert_array_almost_equal,
     assert_raises, assert_allclose, assert_array_equal, assert_equal)
-from scikits.odes import ode
-from scikits.odes.dopri5 import StatusEnumDOP
+from .. import ode
+from ..dopri5 import StatusEnumDOP
 
 
 class SimpleOscillator():

scikits/odes/tests/test_get_info.py

@@ -1,7 +1,13 @@
 from __future__ import print_function
 import numpy as np
 import unittest
-from scikits.odes import ode
+from .. import ode
+from ..sundials import log_error_handler
+
+COMMON_ARGS = {
+    "old_api": False,
+    "err_handler": log_error_handler
+}
 
 
 xs = np.linspace(1, 10, 10)
@@ -20,7 +26,7 @@ def rhs(x, y, ydot):
 
 class GetInfoTest(unittest.TestCase):
     def setUp(self):
-        self.ode = ode('cvode', rhs, old_api=False)
+        self.ode = ode('cvode', rhs, **COMMON_ARGS)
         self.solution = self.ode.solve(xs, np.array([1]))
 
     def test_we_integrated_correctly(self):
@@ -47,7 +53,7 @@ def test_ode_exposes_num_rhs_evals(self):
 
 class GetInfoTestSpils(unittest.TestCase):
     def setUp(self):
-        self.ode = ode('cvode', rhs, linsolver="spgmr", old_api=False)
+        self.ode = ode('cvode', rhs, linsolver="spgmr", **COMMON_ARGS)
         self.solution = self.ode.solve(xs, np.array([1]))
 
     def test_ode_exposes_num_njtimes_evals(self):

scikits/odes/tests/test_odeint.py

@@ -19,8 +19,9 @@
     assert_, TestCase, run_module_suite, assert_array_almost_equal,
     assert_raises, assert_allclose, assert_array_equal, assert_equal)
 
-from scikits.odes.odeint import odeint
-from scikits.odes.sundials.common_defs import DTYPE
+from ..odeint import odeint
+from ..sundials import log_error_handler
+from ..sundials.common_defs import DTYPE
 
 TEST_LAPACK = DTYPE == np.double