SOTA-vs-SOTA: analytic sparse Jacobians on both sides, Assimulo replaced by CasADi

[1-Bart-1]

Two commits. The first gives the Julia time simulations an analytic sparse Jacobian; the second does the same for Python by replacing Assimulo and every hand-written Jacobian with CasADi. Together they make docs/julia_vs_python.md a comparison of two ecosystems rather than of compiled Julia against interpreted NumPy callbacks.

1. Julia: jac=true, sparse=true

Every example built its ODEProblem without jac, so prob.f.jac was nothing and FBDF rebuilt a dense Jacobian by automatic differentiation at every step. Verified directly:

plain    jac=false  proto=nothing         colorvec=nothing
sparse   jac=false  proto=SparseMatrixCSC colorvec=nothing
jac+sp   jac=true   proto=SparseMatrixCSC colorvec=nothing

ModelingToolkit now generates and compiles the analytic, block-tridiagonal Jacobian ahead of time — which is what paper.md:100 already claimed was happening. Scope: time simulations only; the steady-state solves are untouched.

Median ±IQR, BenchmarkTools, 12 s per configuration:

segments	states	AD dense (before)	jac	jac+sparse (after)	speedup
5	36	13.0 ±0.4	12.7 ±0.7	9.1 ±1.0	1.4x
10	66	31.9 ±1.1	27.9 ±0.2	17.1 ±1.9	1.9x
20	126	485.0 ±6.1	438.4 ±6.1	177.9 ±31.6	2.7x
40	246	4332.9 ±22.5	3765.9 ±154.5	603.1 ±10.5	7.2x

The analytic Jacobian alone buys little — the dense factorization then dominates. It is the two together that pay. Cost: jac=true adds 0.4 s (5 segments) to 4.1 s (40) of build time. autodiff= is dead once a Jacobian is supplied (FBDF() and FBDF(autodiff=AutoForwardDiff()) measure the same), so it is dropped along with the now-unused ADTypes imports.

Trajectories unchanged — every example run before and after, compared on final tether shape and whole-trajectory sum:

example	‖shape‖ before → after	nf before → after
Tether_05	38.4647359 → 38.4647372	2110 → 2154
Tether_06	103.5404041 → 103.5404041	2161 → 1225
Tether_06c	103.5404045 → 103.5404045	2579 → 1304
Tether_07	103.6466433 → 103.6466433	1954 → 1164
Tether_08	155.0078612 → 155.0078612	12156 → 5126
Tether_09	100.3370285 → 100.3370285	1 → 1
Tether_10	155.0078612 → 155.0078612	11208 → 5119

Tether_06b, Tether_06_video and Tether_08_video have no test coverage and were run by hand; all three succeed.

2. Python: Assimulo and the hand-written Jacobians are gone

Tether_01–05 had no Jacobian at all and let IDA finite-difference it. Tether_06, 06c, 07 and 08 carried 86–150 lines of pen-and-paper calculus each. Each model is now one CasADi expression graph; ca.jacobian does the rest, and CVODES integrates it with a sparse Newton solve.

file	before	after
Tether_01.py	72	63
Tether_02.py	84	68
Tether_03.py	103	90
Tether_03b.py	123	138
Tether_04.py	140	108
Tether_05.py	192	162
Tether_06.py	314	195
Tether_06c.py	343	225
Tether_07.py	367	215
Tether_08.py	487	343

All ten reproduce the Julia trajectories within the tolerances test/test_tether_*.jl apply — run locally against freshly generated Julia CSVs:

Tether_01   pos max|d|=1.822e-04 ok   vel max|d|=2.182e-05 ok
Tether_02   pos max|d|=6.362e-05 ok   vel max|d|=4.561e-04 ok
Tether_03   pos max|d|=1.648e-04 ok   vel max|d|=1.190e-03 ok
Tether_03b  pos max|d|=3.043e-03 ok   vel max|d|=2.006e-02 ok   (9 crossings)
Tether_04   pos max|d|=1.205e-04 ok   vel max|d|=5.677e-04 ok
Tether_05   pos max|d|=4.739e-05 ok   vel max|d|=6.615e-04 ok
Tether_06   pos max|d|=2.927e-04 ok   vel max|d|=3.254e-04 ok
Tether_06c  pos max|d|=2.929e-04 ok   vel max|d|=4.315e-03 ok
Tether_07   pos max|d|=1.957e-04 ok   vel max|d|=2.609e-04 ok
Tether_08   pos max|d|=1.099e-04 ok   vel max|d|=1.741e-04 ok

CondaPkg.toml swaps assimulo for casadi and adds scipy. A side benefit: the whole Python side is now pip-installable and runnable without a Fortran/SUNDIALS build.

Where Python is genuinely behind

Tether_03b is the one place. CasADi does have event detection — a zero entry in the DAE dictionary — but it is experimental, and on that model it aborts with "tout too far back in direction of integration": after restarting at an event it cannot interpolate back to the requested output point. It fails whether the run starts on, above or below the threshold, and stepping interval-by-interval does not help. So Tether_03b.py walks the output grid itself, watches the sign of the indicator, bisects for the crossing and restarts there — about twenty lines against four for Julia's ContinuousCallback. Tether_06c.py, whose events fire once per segment, uses CasADi's event support successfully.

3. The comparison, now like for like

With both sides analytic and sparse, Julia FBDF and SUNDIALS CVODES-BDF land within about 20% of each other:

segments	Julia jac+sparse	CasADi CVODES sparse
5	9.1 ±1.0	12.3 ±0.5
10	17.1 ±1.9	20.7 ±0.2
20	177.9 ±31.6	170.0 ±4.9
40	603.1 ±10.5	617.1 ±13.8

The previously reported 13–30x was compiled Julia against a hand-written Jacobian evaluated in interpreted NumPy — 707 µs per call against 33 µs for the same Jacobian as a CasADi function. It measured the binding, not the language. docs/julia_vs_python.md, README.md, docs/src/index.md, docs/src/examples.md, docs/src/references.md and docs/src/python.md are updated accordingly.

4. Quasi-steady vs dynamic

examples/quasisteady/benchmark_scaling.jl + docs/images/qsm_vs_dynamic.png, both per simulated second:

segments	quasi-steady [ms/sim s]	dynamic [ms/sim s]	ratio
4	0.174 ±0.004	0.97 ±0.11	5.6
8	0.227 ±0.002	2.53 ±0.37	11
16	0.340 ±0.003	8.73 ±1.84	26
32	0.572 ±0.014	27.94 ±6.82	49

Two traps, both hit while writing this and both now recorded in docs/quasisteady.md: comparing a per-step! cost against a per-simulated-second cost inflates the ratio by 50x at dt = 0.02; and benchmarking step! in a loop without moving the kite measures nothing, because step! writes the converged answer back into te.state_vec.

This bears on Uwe's "only 1.5x faster" against Williams' 1000x: that is a whole-script comparison where the dynamic side's one-off symbolic build dominates, while these are solve times only.

Other findings

• Tether_11 is left alone. Its time simulation goes Unstable with any non-finite-difference Jacobian, including plain FBDF() with no jac at all — so this predates the change. Worth its own issue.
• Tether_09 does no integration (nf=1, both ends fixed, starts at steady state), so the change is inert there.
• sparse=true never sets a colorvec, so sparse AD still costs N residual evaluations per Jacobian; only the factorization benefits.
• .zenodo.json listed one creator while CITATION.cff listed two. Both now list three.

Authorship

Bart van de Lint added as third author (paper.md, CITATION.cff, .zenodo.json), affiliation Open Source AWE. No ORCID — add one if you have it.

Verified locally

Full suite in a clean Julia session, including the CondaPkg resolve and every Python
example through run_python:

Test Summary:                            | Pass  Total   Time
Tethers.jl                               |  341    341  4m51.1s

CondaPkg installs CasADi 3.7.2 from conda-forge, not the 3.8.0 used while developing;
the zero event support and both csparse/lapacklu linear solvers are present in 3.7.2,
checked explicitly.

One bug surfaced this way: examples/quasisteady/benchmark_scaling.jl imported Tether
from QuasiSteady while also importing TetherComponents, which defines a component of
the same name. menu3.jl runs every quasi-steady example into the same Main, so that
would have shadowed the component for anything running afterwards. Fixed in the third
commit, which also adds the benchmark to menu3.jl.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Lw7QQPcUvFcYNV4pSBs7Lg

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 92.28%. Comparing base (7738ec9) to head (2fdd814).

Additional details and impacted files

@@           Coverage Diff           @@
##             main      #33   +/-   ##
=======================================
  Coverage   92.28%   92.28%           
=======================================
  Files           4        4           
  Lines         350      350           
=======================================
  Hits          323      323           
  Misses         27       27           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.