"""The second domain: adaptive step-size control for ODEs. What is being tested is the transfer, not the integrator: everything except the controller is shared, so a difference can only come from the controller. """ from __future__ import annotations import math import os import sys import pytest sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) from odestep.controllers import (ALL, Classical, Elementary, GustafssonPI, MeasuredSafety) from odestep.driver import integrate, validate_against_scipy from stage15_ode import PROBLEMS, TOLS, reference, scaled_err def test_the_driver_reproduces_the_library(): """The instrument before the experiment. The first version of this work carried a hand-typed Dormand-Prince tableau and validated it by measuring convergence orders. That was both an unnecessary risk and, as it turned out, a material one: with the hand-written tableau the measured-safety controller appeared to beat the five-constant PI controller by 8%, and with SciPy's it is 2.6% behind. The integrator is now SciPy's throughout and the check is direct.""" for name, mine, ref, endpoint in validate_against_scipy(): assert abs(mine - ref) / ref < 0.05, (name, mine, ref) assert endpoint < 1e-6, (name, endpoint) def test_every_controller_solves_an_easy_problem(): f = lambda t, y: [-y[0]] for c in ALL: r = integrate(f, 0.0, 10.0, [1.0], c, atol=1e-8, rtol=1e-8) assert r["reached_end"], c.name assert abs(r["y"][0] - math.exp(-10.0)) < 1e-6, c.name def test_the_constant_free_baseline_is_the_one_to_beat(): """The elementary controller is what the error model alone dictates and what every practical code refuses to use. If it were not much worse, there would be nothing for a safety factor - chosen or measured - to buy.""" f = lambda t, y: [y[1], 5.0 * (1 - y[0] * y[0]) * y[1] - y[0]] elem = integrate(f, 0.0, 25.0, [2.0, 0.0], Elementary(), atol=1e-8, rtol=1e-8) meas = integrate(f, 0.0, 25.0, [2.0, 0.0], MeasuredSafety(), atol=1e-8, rtol=1e-8) # margins measured with SciPy's integrator. An earlier version demanded # a factor of two on evaluations and five on rejections; those came from # a hand-written tableau and are not what the library produces. assert elem["nfev"] > 1.2 * meas["nfev"], (elem["nfev"], meas["nfev"]) assert elem["rejected"] > 1.5 * meas["rejected"], (elem["rejected"], meas["rejected"]) def test_where_the_measured_safety_factor_actually_stands(): """The honest position, over the whole grid, with SciPy's integrator. It is NOT ahead of the PI controller - an earlier version of this test asserted that, on numbers produced by a hand-written tableau. What it is: level with a five-constant controller to within a few percent while having none, and clearly ahead of the constant-free elementary rule that every practical code refuses to use. It is behind the three-constant classical controller and the test says so.""" tot = {c.name: 0 for c in ALL} fails = {c.name: 0 for c in ALL} refs = {name: (exact if exact is not None else reference(f, t0, t1, y0)) for name, f, t0, t1, y0, exact in PROBLEMS} for tol in TOLS: for name, f, t0, t1, y0, _exact in PROBLEMS: for c in ALL: r = integrate(f, t0, t1, y0, c, atol=tol, rtol=tol) err = scaled_err(r["y"], refs[name], tol, tol) bad = (not r["reached_end"]) or (math.isfinite(err) and err > 1000.0) tot[c.name] += r["nfev"] fails[c.name] += bad # nobody wins by failing more assert len(set(fails.values())) == 1, fails m = tot[MeasuredSafety().name] pi = tot[GustafssonPI().name] # level with the five-constant PI controller, within five percent assert abs(m - pi) / pi < 0.05, (m, pi) # clearly ahead of the constant-free rule the error model alone dictates assert m < 0.9 * tot[Elementary().name], (m, tot[Elementary().name]) # and honestly behind the three-constant classical one assert m > tot[Classical().name], (m, tot[Classical().name]) def test_the_controller_holds_no_unused_state(): """It was called a credit controller while banking a margin that took no part in any decision - the ledger was vestigial and the name unearned. Nothing may accumulate here that the step factor does not read.""" c = MeasuredSafety() c.reset() fields = set(vars(c)) assert fields == {"predicted", "delta"}, fields