More tests and minor cleanup.

Co-authored-by: Copilot <copilot@github.com>

Author: Jonas Breuling

solve_dae/integrate/_dae/base.py

@@ -8,7 +8,7 @@
 )
 from scipy.optimize._numdiff import approx_derivative
 from scipy.linalg import lu_factor, lu_solve
-from scipy.sparse import csc_matrix, issparse, eye
+from scipy.sparse import csc_matrix, issparse
 from scipy.sparse.linalg import splu
 from .common import select_initial_step
 
@@ -359,8 +359,7 @@ def step(self):
             otherwise.
         """
         if self.status != 'running':
-            raise RuntimeError("Attempt to step on a failed or finished "
-                               "solver.")
+            raise RuntimeError("Attempt to step on a failed or finished solver.")
 
         # if self.n == 0 or self.t == self.t_bound:
         #     # Handle corner cases of empty solver or no integration.
@@ -390,15 +389,14 @@ def dense_output(self):
             Local interpolant over the last successful step.
         """
         if self.t_old is None:
-            raise RuntimeError("Dense output is available after a successful "
-                               "step was made.")
+            raise RuntimeError("Dense output is available after a successful step was made.")
         return self._dense_output_impl()
 
-    def _step_impl(self):
-        raise NotImplementedError
+    # def _step_impl(self):
+    #     raise NotImplementedError
 
-    def _dense_output_impl(self):
-        raise NotImplementedError
+    # def _dense_output_impl(self):
+    #     raise NotImplementedError
 
 
 class DAEDenseOutput:
@@ -438,5 +436,5 @@ def __call__(self, t):
             raise ValueError("`t` must be a float or a 1-D array.")
         return self._call_impl(t)
 
-    def _call_impl(self, t):
-        raise NotImplementedError
+    # def _call_impl(self, t):
+    #     raise NotImplementedError

solve_dae/integrate/_dae/tests/test_base.py

@@ -5,30 +5,66 @@
 from solve_dae.integrate import solve_dae
 
 
-# initial conditions
 y0 = [1]
 yp0 = [0.5 * y0[0]]
-t_span = (0, 1)
+t_span = (0, 1.2)
 
 
 def F(t, y, yp):
     return yp - 0.5 * y
 
 def jac_dense(t, y, yp):
-    return np.eye(len(y)), -0.5 * np.eye(len(yp))
+    return -0.5 * np.eye(len(y)), np.eye(len(yp))
 
 def jac_sparse(t, y, yp):
-    return csc_matrix(np.eye(len(y))), csc_matrix(-0.5 * np.eye(len(yp)))
+    return csc_matrix(-0.5 * np.eye(len(y))), csc_matrix(np.eye(len(yp)))
 
-def jac_wrong_shape_dense(t, y, yp):
-    return np.eye(len(y) + 1), -0.5 * np.eye(len(yp) + 1)
+def jac_wrong_shape_Jyp_dense(t, y, yp):
+    return -0.5 * np.eye(len(y)), np.eye(len(yp) + 1), 
 
-def jac_wrong_shape_sparse(t, y, yp):
-    return csc_matrix(np.eye(len(y) + 1)), csc_matrix(-0.5 * np.eye(len(yp) + 1))
+def jac_wrong_shape_Jy_dense(t, y, yp):
+    return -0.5 * np.eye(len(y) + 1), np.eye(len(yp))
+
+def jac_wrong_shape_Jyp_sparse(t, y, yp):
+    return csc_matrix(-0.5 * np.eye(len(y))), csc_matrix(np.eye(len(yp) + 1))
+
+def jac_wrong_shape_Jy_sparse(t, y, yp):
+    return csc_matrix(-0.5 * np.eye(len(y) + 1)), csc_matrix(np.eye(len(yp)))
+
+jac_wrong_shape_Jyp_dense_constant = (
+    np.eye(2), 
+    -0.5 * np.eye(1),
+)
+
+jac_wrong_shape_Jy_dense_constant = (
+    np.eye(1), 
+    -0.5 * np.eye(2),
+)
+
+jac_wrong_shape_Jyp_sparse_constant = (
+    csc_matrix(np.eye(2)), 
+    csc_matrix(-0.5 * np.eye(1)),
+)
+
+jac_wrong_shape_Jy_sparse_constant = (
+    csc_matrix(np.eye(1)), 
+    csc_matrix(-0.5 * np.eye(2)),
+)
 
 
 parameters_method = ["BDF", "Radau"]
-parameters_jac = [jac_dense, jac_sparse, jac_wrong_shape_dense, jac_wrong_shape_sparse]
+parameters_jac_correct_shape = [jac_dense, jac_sparse]
+parameters_jac_wrong_shape = [
+    jac_wrong_shape_Jyp_dense, 
+    jac_wrong_shape_Jy_dense, 
+    jac_wrong_shape_Jyp_sparse, 
+    jac_wrong_shape_Jy_sparse,
+    jac_wrong_shape_Jyp_dense_constant,
+    jac_wrong_shape_Jy_dense_constant,
+    jac_wrong_shape_Jyp_sparse_constant,
+    jac_wrong_shape_Jy_sparse_constant,
+]
+parameters_jac = parameters_jac_correct_shape + parameters_jac_wrong_shape
 
 
 parameters = product(
@@ -39,7 +75,7 @@ def jac_wrong_shape_sparse(t, y, yp):
 
 @pytest.mark.parametrize("method, jac", parameters)
 def test_jacobian_shape(method, jac):
-    if jac in [jac_wrong_shape_dense, jac_wrong_shape_sparse]:
+    if not callable(jac) or jac in parameters_jac_wrong_shape:
         with pytest.raises(ValueError) as excinfo:
             solve_dae(F, t_span, y0, yp0, method=method, jac=jac)
 
@@ -61,12 +97,34 @@ def test_jacobian_shape(method, jac):
         solve_dae(F, t_span, y0, yp0, method=method, jac=jac)
 
 
-# @pytest.mark.parametrize("method,", parameters_method)
-# def test_step(method):
-#     sol = solve_dae(F, t_span, y0, yp0, method=method)
-#     pass
+@pytest.mark.parametrize("method", parameters_method)
+def test_small_max_step(method):
+    solve_dae(F, t_span, y0, yp0, method=method, max_step=1e-2)
+
+
+def F(t, y, yp):
+    return yp - 1 / (1 - t)
+
+y0 = [0]
+yp0 = [1]
+t_span = (0, 1.2)
+
+
+@pytest.mark.filterwarnings("ignore:divide by zero encountered in scalar divide")
+@pytest.mark.parametrize("method", parameters_method)
+def test_overflow(method):
+    sol = solve_dae(F, t_span, y0, yp0, method=method)
+    assert sol.status == -1
+    assert sol.success == False
+    assert sol.message == "Required step size is less than spacing between numbers."
 
 
 # if __name__ == "__main__":
 #     for params in parameters:
 #         test_jacobian_shape(*params)
+
+#     for params in parameters_method:
+#         test_small_max_step(params)
+
+#     for params in parameters_method:
+#         test_overflow(params)