Merge pull request #54 from bjodah/softer-bounds

Softer bounds checking

pyodesys/core.py

@@ -308,7 +308,7 @@ def predefined(self, y0, xout, params=(), **kwargs):
                                           force_predefined=True, **kwargs)
         return yout, info
 
-    def integrate(self, x, y0, params=(), **kwargs):
+    def integrate(self, x, y0, params=(), atol=1e-8, rtol=1e-8, **kwargs):
         """ Integrate the system of ordinary differential equations.
 
         Solves the initial value problem (IVP).
@@ -372,14 +372,19 @@ def integrate(self, x, y0, params=(), **kwargs):
         if hasattr(self, 'ny'):
             if _y.shape[-1] != self.ny:
                 raise ValueError("Incorrect shape of intern_y0")
-        if isinstance(kwargs.get('atol', None), dict):
-            kwargs['atol'] = [kwargs['atol'][k] for k in self.names]
+        if isinstance(atol, dict):
+            kwargs['atol'] = [atol[k] for k in self.names]
+        else:
+            kwargs['atol'] = atol
+        kwargs['rtol'] = rtol
+
         integrator = kwargs.pop('integrator', None)
         if integrator is None:
             integrator = os.environ.get('PYODESYS_INTEGRATOR', 'scipy')
 
         args = tuple(map(np.atleast_2d, (_x, _y, _p)))
 
+        self._current_integration_kwargs = kwargs
         if isinstance(integrator, str):
             nfo = getattr(self, '_integrate_' + integrator)(*args, **kwargs)
         else:
@@ -741,7 +746,7 @@ def integrate_chained(odes, kw, x, y0, params=(), **kwargs):
     """ Auto-switching between formulations of ODE system.
 
     In case one has a formulation of a system of ODEs which is preferential in
-    the beginning of the intergration this function allows the user to run the
+    the beginning of the integration, this function allows the user to run the
     integration with this system where it takes a user-specified maximum number
     of steps before switching to another formulation (unless final value of the
     independent variables has been reached). Number of systems used i returned

pyodesys/native/sources/_cvode_wrapper.pyx

@@ -52,7 +52,7 @@ def integrate_adaptive(cnp.ndarray[cnp.float64_t, ndim=2, mode='c'] y0,
                        long int mxsteps=0,
                        str iter_type='undecided', int linear_solver=0, str method='BDF',
                        bool with_jacobian=True, bool return_on_root=False,
-                       int autorestart=0, bool return_on_error=False,
+                       int autorestart=0, bool return_on_error=False, bool with_jtimes=False,
                        bool record_rhs_xvals=False, bool record_jac_xvals=False,
                        bool record_order=False, bool record_fpe=False,
                        double get_dx_max_factor=-1.0, bool error_outside_bounds=False,
@@ -119,7 +119,7 @@ def integrate_adaptive(cnp.ndarray[cnp.float64_t, ndim=2, mode='c'] y0,
         systems, atol, rtol, lmm_from_name(_lmm), <double *>y0.data,
         <double *>x0.data, <double *>xend.data, mxsteps,
         &_dx0[0], &_dx_min[0], &_dx_max[0], with_jacobian, iter_type_from_name(_iter_t), linear_solver,
-        maxl, eps_lin, nderiv, return_on_root, autorestart, return_on_error
+        maxl, eps_lin, nderiv, return_on_root, autorestart, return_on_error, with_jtimes
     )
 
     xout, yout = [], []
@@ -162,6 +162,7 @@ def integrate_predefined(cnp.ndarray[cnp.float64_t, ndim=2, mode='c'] y0,
                          long int mxsteps=0,
                          str iter_type='undecided', int linear_solver=0, str method='BDF',
                          bool with_jacobian=True, int autorestart=0, bool return_on_error=False,
+                         bool with_jtimes=False,
                          bool record_rhs_xvals=False, bool record_jac_xvals=False,
                          bool record_order=False, bool record_fpe=False,
                          double get_dx_max_factor=0.0, bool error_outside_bounds=False,
@@ -231,7 +232,8 @@ def integrate_predefined(cnp.ndarray[cnp.float64_t, ndim=2, mode='c'] y0,
     result = multi_predefined[OdeSys](
         systems, atol, rtol, lmm_from_name(_lmm), <double *>y0.data, xout.shape[1], <double *>xout.data,
         <double *>yout.data, mxsteps, &_dx0[0], &_dx_min[0], &_dx_max[0], with_jacobian,
-        iter_type_from_name(_iter_t), linear_solver, maxl, eps_lin, nderiv, autorestart, return_on_error)
+        iter_type_from_name(_iter_t), linear_solver, maxl, eps_lin, nderiv, autorestart,
+        return_on_error, with_jtimes)
 
     for idx in range(y0.shape[0]):
         nreached = result[idx].first

pyodesys/native/sources/standalone_template.cpp

@@ -23,6 +23,7 @@ int main(int argc, char *argv[]){
         ("return-on-root", po::value<bool>()->default_value(false), "Return on root")
         ("autorestart", po::value<int>()->default_value(0), "Autorestart (autonomous)")
         ("return-on-error", po::value<bool>()->default_value(false), "Return on error")
+        ("with-jtimes", po::value<bool>()->default_value(false), "With jtimes")
         ("get-dx-max-factor", po::value<realtype>()->default_value(1.0), "get_dx_max multiplicative factor")
         ("error-outside-bounds", po::value<bool>()->default_value(false), "Return recoverable error to solver when outside bounds")
         ("max-invariant-violation", po::value<realtype>()->default_value(0.0), "Limit at which to return recoverable error when supported by integrator.")
@@ -68,6 +69,7 @@ int main(int argc, char *argv[]){
     bool return_on_root(vm["return-on-root"].as<bool>());
     int autorestart(vm["autorestart"].as<int>());
     bool return_on_error(vm["return-on-error"].as<bool>());
+    bool with_jtimes(vm["with-jtimes"].as<bool>());
 
     std::vector<realtype> params;
     const realtype get_dx_max_factor(vm["get-dx-max-factor"].as<realtype>());
@@ -141,14 +143,14 @@ int main(int argc, char *argv[]){
         xy_ri = cvodes_anyode_parallel::multi_adaptive(
         systems, atol, rtol, lmm, &y0[0], &t0[0], &tend[0], mxsteps, &dx0[0],
         &dx_min[0], &dx_max[0], with_jacobian, iter_type, linear_solver, maxl,
-        eps_lin, nderiv, return_on_root, autorestart, return_on_error);
+        eps_lin, nderiv, return_on_root, autorestart, return_on_error, with_jtimes);
 
     } else {
         yout.resize(systems.size()*ny*nout);
         ri_ro = cvodes_anyode_parallel::multi_predefined(
         systems, atol, rtol, lmm, &y0[0], nout, &tout[0], &yout[0], mxsteps, &dx0[0],
         &dx_min[0], &dx_max[0], with_jacobian, iter_type, linear_solver, maxl,
-        eps_lin, nderiv, autorestart, return_on_error);
+        eps_lin, nderiv, autorestart, return_on_error, with_jtimes);
     }
     // Output:
     for (int si=0; si<systems.size(); ++si){

pyodesys/symbolic.py

@@ -228,8 +228,8 @@ def __init__(self, dep_exprs, indep=None, params=None, jac=True, dfdx=True, firs
             kwargs['param_names'] = [p.name for p in self.params]
 
         self.band = kwargs.get('band', None)  # needed by get_j_ty_callback
-        self.lower_bounds = lower_bounds  # needed by get_f_ty_callback
-        self.upper_bounds = upper_bounds  # needed by get_f_ty_callback
+        self.lower_bounds = None if lower_bounds is None else np.array(lower_bounds)  # needed by get_f_ty_callback
+        self.upper_bounds = None if upper_bounds is None else np.array(upper_bounds)  # needed by get_f_ty_callback
 
         super(SymbolicSys, self).__init__(
             self.get_f_ty_callback(),
@@ -555,11 +555,15 @@ def get_f_ty_callback(self):
         if lb is not None or ub is not None:
             def _bounds_wrapper(t, y, p=(), be=None):
                 if lb is not None:
-                    if np.any(y < lb):
+                    if np.any(y < lb - 10*self._current_integration_kwargs['atol']):
                         raise RecoverableError
+                    y = np.array(y)
+                    y[y < lb] = lb[y < lb]
                 if ub is not None:
-                    if np.any(y > ub):
+                    if np.any(y > ub + 10*self._current_integration_kwargs['atol']):
                         raise RecoverableError
+                    y = np.array(y)
+                    y[y > ub] = ub[y > ub]
                 return cb(t, y, p, be)
             return _bounds_wrapper
         else:

setup.py

@@ -82,7 +82,7 @@ def _path_under_setup(*args):
 _author, _author_email = open(_path_under_setup('AUTHORS'), 'rt').readline().split('<')
 
 extras_req = {
-    'integrators': ['scipy>=0.19.1', 'pyodeint>=0.9.0', 'pycvodes>=0.8.0', 'pygslodeiv2>=0.8.0'],
+    'integrators': ['scipy>=0.19.1', 'pyodeint>=0.9.0', 'pycvodes>=0.8.3', 'pygslodeiv2>=0.8.0'],
     'symbolic': ['sym', 'sympy>=1.1.1'],
     'native': ['pycompilation>=0.4.3', 'pycodeexport>=0.1.1', 'appdirs'],
     'docs': ['Sphinx', 'sphinx_rtd_theme', 'numpydoc'],