From 6d800ec1a6f233dc1e9d06c97992b8efa7a83d93 Mon Sep 17 00:00:00 2001 From: Jonas Breuling Date: Thu, 20 Feb 2025 13:30:30 +0100 Subject: [PATCH 1/5] Fix github workflow. --- .github/workflows/main.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 7f58478..e24bf92 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -40,9 +40,9 @@ jobs: run: pytest --cov app ${{ env.CODECOV_ATS_TESTS }} - name: Upload coverage data - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: - name: covdata + name: covdata-${{ matrix.os }} path: coverage.xml - name: Upload coverage reports to Codecov From a609bfd3ae68051f9d56b3b65c9d2cccfc7016f1 Mon Sep 17 00:00:00 2001 From: Jonas Breuling Date: Thu, 20 Feb 2025 13:34:30 +0100 Subject: [PATCH 2/5] Add pattern to coverage upload. --- .github/workflows/main.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index e24bf92..6edffa9 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -42,10 +42,13 @@ jobs: - name: Upload coverage data uses: actions/upload-artifact@v4 with: + # see https://github.com/actions/upload-artifact/issues/480 name: covdata-${{ matrix.os }} path: coverage.xml - name: Upload coverage reports to Codecov uses: codecov/codecov-action@v4.0.1 with: + # see https://github.com/actions/upload-artifact/issues/480 + pattern: covdata-* token: ${{ secrets.CODECOV_TOKEN }} From c467b681352e56dd61a361becd87270cfaa46c53 Mon Sep 17 00:00:00 2001 From: Jonas Breuling Date: Thu, 20 Feb 2025 13:37:07 +0100 Subject: [PATCH 3/5] Remove pattern again. --- .github/workflows/main.yml | 2 -- 1 file changed, 2 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 6edffa9..9185b28 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -49,6 +49,4 @@ jobs: - name: Upload coverage reports to Codecov uses: codecov/codecov-action@v4.0.1 with: - # see https://github.com/actions/upload-artifact/issues/480 - pattern: covdata-* token: ${{ secrets.CODECOV_TOKEN }} From 7a8d3b05689d91f4538d8652f132324b9724cf61 Mon Sep 17 00:00:00 2001 From: Jonas Breuling Date: Tue, 25 Feb 2025 07:11:51 +0100 Subject: [PATCH 4/5] Fixed bug in consistent initial conditions and cleanup. --- scipy_dae/integrate/_dae/common.py | 55 +++++++++++++++--------------- 1 file changed, 28 insertions(+), 27 deletions(-) diff --git a/scipy_dae/integrate/_dae/common.py b/scipy_dae/integrate/_dae/common.py index 78d637a..593d92b 100644 --- a/scipy_dae/integrate/_dae/common.py +++ b/scipy_dae/integrate/_dae/common.py @@ -160,18 +160,28 @@ def select_initial_step(t0, y0, yp0, t_bound, rtol, atol, max_step): References ---------- - .. [1] TODO: Find a reference. + .. [1] L. F. Shampine, "Starting an ODE solver", November 1977. """ - min_step = 0.0 + safety = 0.8 + min_step = 16 * EPS * abs(t0) threshold = atol / rtol hspan = abs(t_bound - t0) - # compute an initial step size h using yp = y'(t0) + # compute scaling wt = np.maximum(np.abs(y0), threshold) - rh = 1.25 * np.linalg.norm(yp0 / wt, np.inf) / np.sqrt(rtol) + + # error + e = np.linalg.norm(yp0 / wt, np.inf) + + # reciprocal step size + rh = e / np.sqrt(rtol) / safety + + # compute an initial step size h_abs = min(max_step, hspan) if h_abs * rh > 1: h_abs = 1 / rh + + # ensure h_abs >= min_step h_abs = max(h_abs, min_step) return h_abs @@ -220,8 +230,8 @@ def fun_composite(t, z): if not (isinstance(rtol, float) and rtol > 0): raise ValueError("Relative tolerance must be a positive scalar.") - if rtol < 100 * np.finfo(float).eps: - rtol = 100 * np.finfo(float).eps + if rtol < 100 * EPS: + rtol = 100 * EPS print(f"Relative tolerance increased to {rtol}") if np.any(np.array(atol) <= 0): @@ -235,11 +245,6 @@ def fun_composite(t, z): Jy, Jyp = jac(t0, y0, yp0) scale_f = atol + np.abs(f) * rtol - # z0 = np.concatenate([y0, yp0]) - # scale_z = atol + np.abs(z0) * rtol - # dz_norm_old = None - # rate_z = None - # tol = max(10 * EPS / rtol, min(0.03, rtol ** 0.5)) for _ in range(newton_maxiter): for _ in range(chord_iter): @@ -247,17 +252,6 @@ def fun_composite(t, z): y0 += dy yp0 += dyp - # dz = np.concatenate([dy, dyp]) - # with np.errstate(divide='ignore'): - # dz_norm = norm(dz / scale_z) - # if dz_norm_old is not None: - # rate_z = dz_norm / dz_norm_old - - # if (dz_norm == 0 or (rate_z is not None and rate_z / (1 - rate_z) * dz_norm < safety * tol)): - # return y0, yp0, f - - # dz_norm_old = dz_norm - f = fun(t0, y0, yp0, *args) error = norm(f / scale_f) if error < safety: @@ -271,7 +265,7 @@ def fun_composite(t, z): def qrank(A): """Compute QR-decomposition with column pivoting of A and estimate the rank.""" Q, R, p = qr(A, pivoting=True) - tol = max(A.shape) * np.finfo(float).eps * abs(R[0, 0]) + tol = max(A.shape) * EPS * abs(R[0, 0]) rank = np.sum(abs(np.diag(R)) > tol) return rank, Q, R, p @@ -353,10 +347,17 @@ def solve_underdetermined_system(f, Jy, Jyp, free_y, free_yp): # [S21, S22] [w1] = d2 # [w2] # using column pivoting QR-decomposition - w_ = np.zeros(RS.shape[1]) - w_[:rankS] = solve_triangular(RS[:rankS, :rankS], (QS.T @ d2[:rankS])) - w = np.zeros_like(w_) - w[pS] = w_ + # [RS11, RS12] [v1] = [c1] + # [ 0, 0] [v2] [c2] + # with v2 = 0 this gives + # RS11 @ v1 = c1 + c = QS.T @ d2 + v = np.zeros(RS.shape[1]) + v[:rankS] = solve_triangular(RS[:rankS, :rankS], c[:rankS]) + + # apply permutation + w = np.zeros_like(v) + w[pS] = v # set w2' = 0 and solve the remaining system # [R11] w1' = d1 - [S11, S12] [w1] From 25cbb48d515dd99ead0eedc214385ef61ce924a0 Mon Sep 17 00:00:00 2001 From: Jonas Breuling Date: Tue, 25 Feb 2025 07:18:44 +0100 Subject: [PATCH 5/5] Cleanup. --- scipy_dae/integrate/_dae/common.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/scipy_dae/integrate/_dae/common.py b/scipy_dae/integrate/_dae/common.py index 593d92b..967c439 100644 --- a/scipy_dae/integrate/_dae/common.py +++ b/scipy_dae/integrate/_dae/common.py @@ -126,9 +126,6 @@ def __call__(self, t): return ys, yps -# TODO: Compare this with -# - ddassl.f by Petzold -# - epsode.f by Bryne and Hindmarsh def select_initial_step(t0, y0, yp0, t_bound, rtol, atol, max_step): """Empirically select a good initial step. @@ -194,7 +191,7 @@ def consistent_initial_conditions(fun, t0, y0, yp0, jac=None, fixed_y0=None, References ---------- - .. [1] L. F. Shampine, "Solving 0 = F(t, y(t), y′(t)) in Matlab", Journal + .. [1] L. F. Shampine, "Solving 0 = F(t, y(t), y'(t)) in Matlab", Journal of Numerical Mathematics, vol. 10, no. 4, 2002, pp. 291-310. """ n = len(y0)