From c417b7d97c977145583ade738f407bcfc39f2b8a Mon Sep 17 00:00:00 2001 From: Thomas Duvernay Date: Mon, 10 Oct 2022 18:08:03 +1100 Subject: [PATCH] TST: Update segment selection and add a test for LSODA interpolant. --- scipy/integrate/_ivp/common.py | 22 ++++++---- scipy/integrate/_ivp/tests/test_ivp.py | 59 ++++++++++++++++++++------ 2 files changed, 60 insertions(+), 21 deletions(-) diff --git a/scipy/integrate/_ivp/common.py b/scipy/integrate/_ivp/common.py index 7f6d2331100a..271e73ea9905 100644 --- a/scipy/integrate/_ivp/common.py +++ b/scipy/integrate/_ivp/common.py @@ -144,13 +144,21 @@ class OdeSolution: interpolants : list of DenseOutput with n_segments elements Local interpolants. An i-th interpolant is assumed to be defined between ``ts[i]`` and ``ts[i + 1]``. + alt_segment : boolean + Requests the alternative interpolant segment selection scheme. At each + solver integration point, two interpolant segments are available. The + default (False) and alternative (True) behaviours select the segment + for which the requested time corresponded to ``t`` and ``t_old``, + respectively. This functionality is only relevant for testing the + interpolants' accuracy: different integrators use different + construction strategies. Attributes ---------- t_min, t_max : float Time range of the interpolation. """ - def __init__(self, ts, interpolants): + def __init__(self, ts, interpolants, alt_segment=False): ts = np.asarray(ts) d = np.diff(ts) # The first case covers integration on zero segment. @@ -169,20 +177,19 @@ def __init__(self, ts, interpolants): self.t_min = ts[0] self.t_max = ts[-1] self.ascending = True + self.side = "right" if alt_segment else "left" self.ts_sorted = ts else: self.t_min = ts[-1] self.t_max = ts[0] self.ascending = False + self.side = "left" if alt_segment else "right" self.ts_sorted = ts[::-1] def _call_single(self, t): # Here we preserve a certain symmetry that when t is in self.ts, # then we prioritize a segment with a lower index. - if self.ascending: - ind = np.searchsorted(self.ts_sorted, t, side='left') - else: - ind = np.searchsorted(self.ts_sorted, t, side='right') + ind = np.searchsorted(self.ts_sorted, t, side=self.side) segment = min(max(ind - 1, 0), self.n_segments - 1) if not self.ascending: @@ -215,10 +222,7 @@ def __call__(self, t): t_sorted = t[order] # See comment in self._call_single. - if self.ascending: - segments = np.searchsorted(self.ts_sorted, t_sorted, side='left') - else: - segments = np.searchsorted(self.ts_sorted, t_sorted, side='right') + segments = np.searchsorted(self.ts_sorted, t_sorted, side=self.side) segments -= 1 segments[segments < 0] = 0 segments[segments > self.n_segments - 1] = self.n_segments - 1 diff --git a/scipy/integrate/_ivp/tests/test_ivp.py b/scipy/integrate/_ivp/tests/test_ivp.py index 2fb294823af5..8f0b4ec89319 100644 --- a/scipy/integrate/_ivp/tests/test_ivp.py +++ b/scipy/integrate/_ivp/tests/test_ivp.py @@ -1,5 +1,5 @@ from itertools import product -from numpy.testing import (assert_, assert_allclose, +from numpy.testing import (assert_, assert_allclose, assert_array_less, assert_equal, assert_no_warnings, suppress_warnings) import pytest from pytest import raises as assert_raises @@ -135,6 +135,18 @@ def sol_complex(t): return y.reshape((1, -1)) +def fun_dense_output_LSODA(t, y): + return y * (t - 2) + + +def jac_dense_output_LSODA(t, y): + return t - 2 + + +def sol_dense_output_LSODA(t): + return np.exp(t ** 2 / 2 - 2 * t + np.log(0.05) - 6) + + def compute_error(y, y_true, rtol, atol): e = (y - y_true) / (atol + rtol * np.abs(y_true)) return np.linalg.norm(e, axis=0) / np.sqrt(e.shape[0]) @@ -201,11 +213,7 @@ def test_integration(): e = compute_error(yc, yc_true, rtol, atol) assert_(np.all(e < 5)) - # LSODA for some reasons doesn't pass the polynomial through the - # previous points exactly after the order change. It might be some - # bug in LSOSA implementation or maybe we missing something. - if method != 'LSODA': - assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15) + assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15) def test_integration_complex(): @@ -538,9 +546,7 @@ def test_max_step(): e = compute_error(yc, yc_true, rtol, atol) assert_(np.all(e < 5)) - # See comment in test_integration. - if method is not LSODA: - assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15) + assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15) assert_raises(ValueError, method, fun_rational, t_span[0], y0, t_span[1], max_step=-1) @@ -584,9 +590,7 @@ def test_first_step(): e = compute_error(yc, yc_true, rtol, atol) assert_(np.all(e < 5)) - # See comment in test_integration. - if method is not LSODA: - assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15) + assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15) assert_raises(ValueError, method, fun_rational, t_span[0], y0, t_span[1], first_step=-1) @@ -711,6 +715,37 @@ def early_event(t, y): assert res.t_events[0][0] == 7 +def test_dense_output_LSODA(): + rtol = 1e-3 + atol = 1e-6 + y0 = [0.05] + t_span = [-2, 2] + first_step = 1e-3 + res = solve_ivp(fun_dense_output_LSODA, t_span, y0, method="LSODA", + dense_output=True, first_step=first_step, max_step=1, + rtol=rtol, atol=atol, jac=jac_dense_output_LSODA) + + assert_equal(res.t[0], t_span[0]) + assert_equal(res.t[-1], t_span[-1]) + assert_allclose(first_step, np.abs(res.t[1] - t_span[0])) + assert res.t_events is None + assert res.success + assert_equal(res.status, 0) + + y_true = sol_dense_output_LSODA(res.t) + e = compute_error(res.y, y_true, rtol, atol) + assert_array_less(e, 5) + + tc = np.linspace(*t_span) + yc_true = sol_dense_output_LSODA(tc) + yc = res.sol(tc) + + e = compute_error(yc, yc_true, rtol, atol) + assert_array_less(e, 5) + + assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15) + + def test_no_integration(): for method in ['RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA']: sol = solve_ivp(lambda t, y: -y, [4, 4], [2, 3],