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

.github/workflows/draft-paper.yml

@@ -0,0 +1,24 @@
+name: Draft paper
+
+on:
+  push:
+    paths:
+      - paper/**
+      - .github/workflows/draft-paper.yml
+  workflow_dispatch:
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    name: Build the JOSS paper
+    steps:
+      - uses: actions/checkout@v4
+      - name: Build the draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          paper-path: paper/paper.md
+      - uses: actions/upload-artifact@v4
+        with:
+          name: paper
+          path: paper/paper.pdf

.gitignore

@@ -29,3 +29,4 @@ bin/kps-image-.so
 bin/kps-image-1.12.so.bak
 paper/paper.bib
 paper/paper.pdf
+paper/paper_tether.pdf

.zenodo.json

@@ -6,6 +6,15 @@
       "name": "Fechner, Uwe",
       "affiliation": "Delft University of Technology",
       "orcid": "0009-0008-2532-9458"
+    },
+    {
+      "name": "Bertozzi, Andrea",
+      "affiliation": "Politecnico di Milano",
+      "orcid": "0009-0005-1214-1737"
+    },
+    {
+      "name": "van de Lint, Bart",
+      "affiliation": "Open Source AWE"
     }
   ],
   "keywords": [

CHANGELOG.md

@@ -1,3 +1,81 @@
+### Tethers v2.1.0 (unreleased)
+
+#### Added
+
+- `examples/python/bench_casadi.py`, deriving the Jacobian of the `Tether_08`
+  model with `ca.jacobian` instead of by hand, checking it against the
+  hand-derived `calc_acc_jac` of `Tether_08.py`, and timing the SUNDIALS
+  integrators CasADi ships; needs `pip install casadi` and is not part of the
+  tutorial
+- `examples/quasisteady/benchmark_scaling.jl` and
+  `docs/images/qsm_vs_dynamic.png`, comparing the cost of the quasi-steady model
+  against the dynamic one over 4 to 32 tether segments, both per simulated
+  second; the quasi-steady solve is 5.6x cheaper at four segments and 49x at
+  thirty-two
+- `.github/workflows/draft-paper.yml`, building the JOSS paper with the
+  journal's own `openjournals/inara` image on every push that touches `paper/`
+  and uploading it as the `paper` artifact, so no local Docker or LaTeX
+  toolchain is needed
+- Bart van de Lint as a third author of the JOSS paper, and Andrea Bertozzi to
+  `.zenodo.json`, which listed only one creator
+
+#### Changed
+
+- every example's time simulation now builds its `ODEProblem` with
+  `jac=true, sparse=true`, so ModelingToolkit generates and compiles the
+  analytic, sparse Jacobian ahead of time instead of the solver rebuilding a
+  dense one by automatic differentiation at every step; the `autodiff=` keyword
+  these solves passed to `FBDF` is now redundant and was dropped, and with it
+  the `ADTypes` import where nothing else used it
+- `examples/Tether_11.jl` keeps its finite-difference Jacobian: with a forward-mode or an analytic one, `FBDF` drives dt below eps at t = 0 and returns `Unstable`. This predates the change - plain `FBDF()` with no Jacobian at all already fails there - so the example is left as it was
+- the steady-state solves are unchanged; only the time simulations were touched
+- the Python examples no longer use Assimulo. Each model is now written once as
+  a CasADi expression graph, CasADi derives its exact Jacobian and sparsity, and
+  SUNDIALS' CVODES integrates it with a sparse Newton solve - the same
+  fixed-leading-coefficient BDF family as the Julia examples' `FBDF`. The
+  hand-written Jacobians of `Tether_06.py`, `Tether_06c.py`, `Tether_07.py` and
+  `Tether_08.py`, 86 to 150 lines of calculus each, are gone; the ten examples
+  lost about 700 lines between them. All ten still reproduce the Julia
+  trajectories within the tolerances `test/test_tether_*.jl` apply
+- `CondaPkg.toml` swaps `assimulo` for `casadi` and adds `scipy`, which
+  `Tether_08.py`'s steady-state solve uses
+- `examples/python/bench_casadi.py` now imports the model from `Tether_08.py`
+  instead of defining its own copy
+- `Tether_03b.py` locates the taut/slack crossings itself, by bisection between
+  output points: CasADi's event detection (a `zero` entry in the DAE dictionary)
+  is experimental and aborts on this model with "tout too far back in direction
+  of integration". `Tether_06c.py`, whose events fire once per segment, uses it
+  successfully
+- `paper/paper.md` describes what the package now does: the Python examples use
+  CasADi and SUNDIALS rather than Assimulo's IDA, both implementations use an
+  analytic sparse Jacobian, and the performance comparison reports them within
+  about 20% of each other instead of the earlier 13-30x, with an explanation of
+  where that number came from. The AI usage disclosure covers the Jacobian and
+  CasADi work
+- `paper/paper.bib` cites CasADi [Andersson2019] instead of Assimulo, which the
+  examples no longer use
+- `paper/build` no longer needs a container. It used the third-party
+  `openbases/openbases-pdf` image, which does not produce the journal's layout
+  anyway; it now builds a readable PDF with pandoc and a local LaTeX
+  installation in a couple of seconds, and points at the CI workflow for the
+  journal's own layout
+- `docs/julia_vs_python.md` and `paper/paper.md` record what CasADi's `jit=True`
+  is worth: it compiles the expression tape to machine code instead of walking
+  it inside `libcasadi`, gains a further 1.7-1.9x with bit-identical
+  trajectories, and costs 0.6 to 7.3 s of compilation. The examples leave it
+  off, since a tutorial example is run once
+- `README.md` and `docs/src/index.md` claimed Julia was "about 16 times faster"
+  and "13 to 30 times faster" than Python, and carried a comparison table whose
+  Python timings predate the CasADi rewrite. The claims are corrected and the
+  table now gives lines of code only, with the timings measured in one place
+  instead of duplicated
+- the benchmarks in `docs/` were labelled as run on a Ryzen 9 7950X, which is
+  the machine the older numbers in `README.md` came from; the ones added here
+  were run on an Intel Core i7-11850H
+- `docs/julia_vs_python.md` now compares like with like: both sides get an
+  analytic sparse Jacobian and a BDF integrator, which closes most of the gap it
+  used to report
+
 ### Tethers v2.0.0 2026-09-11
 #### Added
 - `examples/Project.toml` and `test/Project.toml`, joined to the root package as workspace members on Julia 1.12, so the example and test dependencies no longer bloat the main `Project.toml`

CITATION.cff

@@ -16,6 +16,9 @@ authors:
     given-names: "Andrea"
     email: andrea.bertozzi@polimi.it
     orcid: "https://orcid.org/0009-0005-1214-1737"
+  - family-names: "van de Lint"
+    given-names: "Bart"
+    email: bart@vandelint.net
 doi: 10.5281/zenodo.19220562
 repository-code: "https://github.com/ufechner7/Tethers.jl"
 url: "https://ufechner7.github.io/Tethers.jl/dev/"

CondaPkg.toml

@@ -1,5 +1,6 @@
 [deps]
 ipython = ""
-assimulo = "==3.4.3"
+casadi = ""
 matplotlib = "==3.8.2"
 numpy = ""
+scipy = ""

README.md

@@ -126,7 +126,7 @@ From the Julia prompt execute:
 run_python("Tether_01")
 ```
 
-This will install Python, Matplotlib and Assimulo and execute the script `Tether_01.py`.
+This will install Python, Matplotlib, NumPy, SciPy and CasADi and execute the script `Tether_01.py`.
 
 **Python code:** [Tether_01.py](https://github.com/ufechner7/Tethers.jl/blob/main/examples/python/Tether_01.py)
 
@@ -142,35 +142,42 @@ at the Julia prompt.
 If you compare the Python and the Julia scripts you can see that:
 
 - the Julia script is shorter and easier to read
-- Julia is about 16 times faster when running the simulation  
+- both reach the same speed, but only once each is given an analytic, sparse Jacobian
 
-For a stiff, segmented tether (example 6 and 7) the Julia solvers are 13 to 30 times faster than Python.
+For a stiff, segmented tether the two are within about 20% of each other; see
+[docs/julia_vs_python.md](https://github.com/ufechner7/Tethers.jl/blob/main/docs/julia_vs_python.md)
+for the measurements and for what each ecosystem needs to get there.
 
 Have a look at the [Examples](https://ufechner7.github.io/Tethers.jl/dev/examples/) that teach you how to construct a full tether model step by step.
 
 ## Overall comparison
 
-Execution time for a simulation of 10s duration with logging the state every 20ms.
-Relative and absolute tolerance: $1.0^{-6}$. CPU: Ryzen 9 7950X.
-LOC counts exclude blank lines and comment lines.
-
-| Test-case                   | Lines of code (LOC) Julia | LOC Python | Time Julia [ms] | Time Python [ms] |
-|:----------------------------------|:-------------------:|:----------:|:---------------:|:---:|
-|Falling mass (1)                   |     36              | 49         | 0.17            | 2.6 |
-|Non-linear Spring damper (3)       |     45              | 80         | 0.64            | 20  |
-|ditto with callbacks (3b, 3c)      |     53              | 97         | 1.0             | 25  |
-|swinging tether, 5 segments (5)    |    103              | 153        | 2.7             | 47  |
-|Dyneema tether, reeling out (6)    |    115              | 240        | 2.4             | 37   |
-|ditto with callbacks       (6c)    |    201              | 271        | 4.2             | 56   |
-|Dyneema, reeling out with drag (7) |    169              | 282        | 2.4             | 70   |  
+Lines of code, excluding blank lines and comment lines. Execution times, which depend on
+the Jacobian each side is given, are measured in
+[docs/julia_vs_python.md](https://github.com/ufechner7/Tethers.jl/blob/main/docs/julia_vs_python.md)
+rather than duplicated here.
+
+| Test-case                          | LOC Julia | LOC Python |
+|:-----------------------------------|:---------:|:----------:|
+|Falling mass (1)                    |     36    |     48     |
+|Non-linear Spring damper (3)        |     45    |     68     |
+|ditto with callbacks (3b)           |     53    |    108     |
+|swinging tether, 5 segments (5)     |    102    |    130     |
+|Dyneema tether, reeling out (6)     |    115    |    154     |
+|ditto with callbacks       (6c)     |    201    |    180     |
+|Dyneema, reeling out with drag (7)  |    168    |    173     |
 
 **Tradeoff Julia vs Python:** In Julia, the code is compiled before it is executed, which can cause about 5 to 30 seconds delay when running a simulation the first time, but speeds up the execution a lot afterward. In addition,
 the Julia code is much more compact and better readable due to the use of symbolic differential equations.
 
-In Python, the IDA solver needs an analytic Jacobian to handle a simulation with the very
-stiff Dyneema tether at all, which increases code size and complexity. Even then, the Julia solvers achieve 13 to 30 times the performance. On the other hand, installing the required packages in Python takes less
-than a minute, while installing and compiling the Julia software might take half an hour. The reason is,
-that currently Julia packages are only distributed as source code and have to be compiled locally.
+Both sides need an exact, sparse Jacobian to handle the very stiff Dyneema tether at all.
+Julia gets one from ModelingToolkit with `jac=true, sparse=true`; Python gets one from CasADi,
+which differentiates the model expression and finds its sparsity by itself. With both compiled
+and sparse the two land within about 20% of each other - see
+[docs/julia_vs_python.md](https://github.com/ufechner7/Tethers.jl/blob/main/docs/julia_vs_python.md).
+The remaining trade-off is setup: installing the required packages in Python takes less than a
+minute, while installing and compiling the Julia software might take half an hour, because Julia
+packages are distributed as source code and have to be compiled locally.
 
 See also: [Why Julia?](https://ufechner7.github.io/2022/08/13/why-julia.html) and read the [documentation](https://ufechner7.github.io/Tethers.jl/dev/) or go straight to the [examples](https://ufechner7.github.io/Tethers.jl/dev/examples/).
 

docs/julia_vs_python.md

@@ -1,44 +1,101 @@
-# Julia vs. Python performance — Tether_06
-
-Comparison of the two `Tether_06` implementations, both solving the same DAE (5-segment
-damped tether, taut/slack spring with a smoothed 1% compression floor) with an exact
-analytic Jacobian.
-
-| | Julia (`FBDF`, ModelingToolkit) | Python (`IDA`, Assimulo) |
-|---|---|---|
-| Solve time | ~0.002–0.003s | ~0.041s |
-| Output points | 501 (`saveat`, exact) | 501 (resampled onto exact grid) |
-| Steps taken | not reported by SciML by default | 436 |
-| Jacobian evals | analytic, symbolically generated (compiled once, ~free per call) | 83 (analytic, hand-derived) |
-| Setup/compile cost | `mtkcompile` symbolic pipeline (one-time, amortized across repeat solves — not counted here) | negligible (Python interpreted, no compile step) |
-| Lines of code | 140 (`examples/Tether_06.jl`) | 313 (`examples/python/Tether_06.py`) |
-
-**Julia is about 15–20x faster** than Python for this problem, even with both now using
-exact analytic Jacobians. The Julia code is also much more compact, at less than half the
-line count of the Python version.
-
-## Why the gap remains
-
-Despite both being "Newton + analytic Jacobian" DAE solvers:
-
-- **Compiled vs. interpreted inner loop**: ModelingToolkit generates and compiles native
-  Julia functions for the residual and Jacobian ahead of time; Python's `res`/`jac`
-  methods re-run interpreted NumPy code (with per-call array allocation, Python-level
-  loops, function-call overhead) on every single Newton iteration.
-- **Allocation overhead**: Julia's solve reused ~5,000 small allocations total across the
-  whole 10s run; Python's `res`/`jac` allocate fresh NumPy arrays (segment vectors, outer
-  products, etc.) on every call — 436 steps × several evaluations per step × several
-  small array allocations each adds up in Python's interpreter/GC overhead in a way it
-  doesn't in compiled Julia.
-- **Fixed cost per Python↔C boundary crossing**: Assimulo's IDA is itself compiled
-  Sundials C code, but every residual/Jacobian evaluation has to cross back into the
-  Python interpreter, whereas Julia's whole pipeline (solver + residual + Jacobian) is
-  native compiled code with no language-boundary crossings.
-
-## Context
-
-Both are fast in absolute terms for this problem (2ms vs 41ms is imperceptible either
-way). The practically important result is that the Python side went from **79.2s to
-0.04s** (~2000x) over the course of this optimization work, closing nearly all of the
-gap that existed purely from using a finite-difference Jacobian instead of an analytic
-one.
+# Julia vs. Python performance
+
+Both ecosystems can give a stiff tether model an analytic, symbolically
+generated, sparse Jacobian, and neither is fast without one. This page compares
+them on the `Tether_08` model (ten segments unless stated otherwise, tether
+drag, reel-out, one fixed and one free end point), over 10 s, sampled every 20
+ms, at a relative and absolute tolerance of $10^{-6}$, from a straight-line
+initial condition. CPU: Intel Core i7-11850H; Julia 1.13, ModelingToolkit 11,
+CasADi 3.8.
+
+|            | Julia                          | Python                                           |
+| ---------- | ------------------------------ | ------------------------------------------------ |
+| Model      | ModelingToolkit `System`       | CasADi `SX` expression graph                     |
+| Jacobian   | `ODEProblem(sys, …; jac=true)` | `ca.jacobian(ydot, y)`                           |
+| Sparsity   | `sparse=true`                  | found by CasADi, always                          |
+| Integrator | `FBDF`                         | SUNDIALS CVODES, `linear_multistep_method='bdf'` |
+
+`FBDF` is a fixed-leading-coefficient BDF, which is the variant SUNDIALS
+implements, so CVODES-BDF is the like-for-like match. IDAS solves the same
+formulas for an implicit DAE; it is 3-11% slower than CVODES here, so the
+examples use CVODES.
+
+## The Jacobian dominates
+
+Solve time in ms, median and interquartile range over the samples
+`BenchmarkTools` fits into twelve seconds per configuration:
+
+| segments | states |     Julia AD |   Julia `jac` | Julia `jac+sparse` | CasADi sparse | CasADi sparse + `jit` |
+| -------: | -----: | -----------: | ------------: | -----------------: | ------------: | --------------------: |
+|        5 |     36 |    13.0 ±0.4 |     12.7 ±0.7 |           9.1 ±1.0 |     12.4 ±1.2 |              7.2 ±0.5 |
+|       10 |     66 |    31.9 ±1.1 |     27.9 ±0.2 |          17.1 ±1.9 |     20.4 ±1.5 |             11.9 ±0.7 |
+|       20 |    126 |   485.0 ±6.1 |    438.4 ±6.1 |        177.9 ±31.6 |    164.8 ±1.6 |             90.6 ±2.0 |
+|       40 |    246 | 4332.9 ±22.5 | 3765.9 ±154.5 |        603.1 ±10.5 |    560.5 ±8.7 |            299.5 ±2.0 |
+
+`jit=True` makes CasADi emit C for its expression tape and compile it with gcc,
+instead of walking the tape inside `libcasadi`. It is worth a further 1.7-1.9x
+and the trajectories come out bit-identical, but it needs a C compiler at run
+time and costs 0.6 s of compilation at five segments and 7.3 s at forty, so the
+examples leave it off: they are run once, and the compilation costs more than
+the solve saves. Turn it on for a model that is solved repeatedly:
+
+```python
+sim = ca.integrator('sim', 'cvodes', dae, 0.0, grid, {
+    'jit': True, 'compiler': 'shell',
+    'jit_options': {'flags': ['-O3'], 'compiler': 'gcc'}, ...})
+```
+
+CasADi with a dense linear solver instead costs 13.9, 28.7, 1037.2 and 4551.8
+ms, so the sparse factorization is worth 7.4x at forty segments there. The Julia
+column comes from `ODEProblem(sys, …; jac, sparse)`, the CasADi one from
+`examples/python/bench_casadi.py`.
+
+The Jacobian is block-tridiagonal: 591 of 4356 entries at ten segments, 2391 of
+60516 at forty. The analytic Jacobian alone buys little, because the dense
+factorization then dominates; it is the two together that pay, by 1.4x at five
+segments and 7.2x at forty; the same holds on the CasADi side.
+
+With `jac=true` the `autodiff` keyword no longer affects anything: the solver
+uses the supplied Jacobian and never differentiates. `FBDF()` and
+`FBDF(autodiff=AutoForwardDiff())` measure the same within noise.
+
+Head to head, with both sides analytic and sparse, Julia's `FBDF` and SUNDIALS
+CVODES-BDF are within 20% of each other, Julia ahead on the small models and
+CVODES on the large ones. CVODES beats IDAS at every size here (10.8 against
+13.2 ms at five segments, 601 against 703 at forty), so the implicit-DAE
+formulation costs a little.
+
+Generating the Jacobian is not free. ModelingToolkit needs 0.5 s at five
+segments and 4.1 s at forty, on top of `mtkcompile`; `ca.jacobian` needs 2 ms
+and 16 ms. For a script that solves once, that build cost outweighs the saving;
+it pays back when a compiled model is re-solved.
+
+## Where the old numbers came from
+
+Earlier versions of this page reported Julia as 13 to 30 times faster than
+Python. That compared compiled Julia against a hand-derived Jacobian evaluated
+in interpreted NumPy, which costs 707 µs per call against 33 µs for the same
+Jacobian as a CasADi function — a factor of 21 in the inner loop of every Newton
+iteration. It measured the binding, not the language. With both sides compiled
+and sparse, they are within 20% of each other.
+
+The code size changed with it. The Python examples used to carry the Jacobian of
+every model as hand-written calculus - 86 to 150 lines per file in
+`Tether_06.py`, `Tether_06c.py`, `Tether_07.py` and `Tether_08.py`.
+`ca.jacobian` replaced all of it, and the ten examples lost about 700 lines
+between them. The Julia examples never wrote a Jacobian at all.
+
+## Caveats
+
+- `FBDF` and CVODES are different BDF implementations, with their own step-size
+  and order heuristics; they take different paths through the same problem.
+- The CasADi timings cross into Python once per solve, not once per Newton
+  iteration, because the whole model is compiled into the integrator. A model
+  driven from a Python callback per step pays an interpreter cost that the table
+  above avoids.
+- The tether's taut/slack switch is a hard step. Its exact derivative is zero on
+  either side, which is what ModelingToolkit, CasADi's `if_else` and ForwardDiff
+  all report. Finite differences instead return a large secant slope across the
+  step, and some configurations only converge because of it — see the
+  steady-state solves of `Tether_08`–`Tether_11`, and `Tether_11`'s time
+  simulation.