Faster solver: Brent search, parallel figure tabs, optional Numba (results unchanged)

.github/workflows/ci.yml

@@ -20,7 +20,7 @@ jobs:
           cache: pip
 
       - name: Install
-        run: pip install -e ".[dev]"
+        run: pip install -e ".[dev,fast]"
 
       - name: Lint
         run: |

CLAUDE.md

@@ -19,7 +19,7 @@ Read these before substantial changes:
 
 ## Environment
 - Conda env: `outflow` (Python 3.12)
-- Install: `conda activate outflow; pip install -e ".[dev]"`
+- Install: `conda activate outflow; pip install -e ".[dev,fast]"` (`fast` = Numba, optional)
 - Run the app: `streamlit run app/streamlit_app.py`
 - Tests: `pytest` (fast); `pytest -m slow` for the full paper-reproduction set
 - Lint: `ruff check . ; ruff format --check .`
@@ -36,7 +36,17 @@ Read these before substantial changes:
   output, now stored as package data in `src/outflow/data/paper/`, not `legacy/outflow.py` (which has a temperature
   bug and auto-enables cooling).
 - Defaults follow the paper code exactly (values, rounded constants, tolerances, numerical
-  method). Changing any of them is a decision for the user.
+  method). Changing any of them is a decision for the user. One such decision so far: the last
+  stage of the ρ₀ search uses Brent's method (`docs/decisions.md`, "Speed").
+- A speed-up that changes results is a numerics change (the user's decision). To show a
+  speed-up changes nothing, use `tests/unit/test_integrator.py` (bit-identical RK step) and
+  compare all 14 paper runs before and after.
+- `src/outflow/_compiled.py` is a Numba copy of the trial integration. Any change to
+  `physics.py`, `equations.py`, `integrator.py` or `boundary._sonic_trial_python` must be made
+  there too, keeping every exponent opaque (`x ** p[TWO]`, never `x**2`) and sums compensated;
+  `pytest tests/unit/test_compiled.py` (and `-m slow`) must stay bit-identical.
+- The app script keeps its `if __name__ == "__main__": main()` guard: the figure tabs' worker
+  processes re-import it.
 - Cooling is off by default and is not exposed in the UI.
 
 ## Conventions

README.md

@@ -37,10 +37,12 @@ Python 3.12. With conda:
 ```powershell
 conda create -n outflow python=3.12
 conda activate outflow
-pip install -e ".[dev]"
+pip install -e ".[dev,fast]"
 ```
 
-Any other Python 3.12 environment works the same way with `pip install -e ".[dev]"`.
+Any other Python 3.12 environment works the same way. The `fast` extra installs Numba, which
+compiles the slowest part of the solver: same results to the last digit, about 30 times faster
+(a run takes about 0.1 s instead of 1.4 s). Without it, everything works in plain Python.
 
 ## Quickstart
 
@@ -57,7 +59,7 @@ profile.to_csv("wind.csv")  # r, v, rho, P, Pc, T, vA with a header of the input
 cr_driven = WindParameters.from_astro_units(thermal_fraction=0.1, cr_fraction=0.8)
 ```
 
-The run takes a few seconds. Main entry points:
+The run takes about a second (0.1 s with the `fast` extra). Main entry points:
 
 - `WindParameters` – the inputs (cgs), with `paper_default()`, `from_astro_units()` and `replace()`.
 - `solve_wind(params, outer_radius=None)` – returns a `WindProfile`: arrays `r, v, rho, P, Pc,
@@ -78,7 +80,8 @@ streamlit run app/streamlit_app.py
 
 Tabs: **Explore** (change any input, overlay runs, download CSV), **Fig. 6–9** (recompute the
 paper's run sets and compare them with the published curves) and **About**. Every run is
-computed when you press its button; a figure tab takes 20–40 s.
+computed when you press its button. With the `fast` extra a figure tab takes about a second;
+without it, the tab solves its runs in parallel and takes about 5 s on a multi-core computer.
 
 ## How the results compare with the paper
 
@@ -87,6 +90,11 @@ published runs are recomputed and compared between 0.2 r_sb and the outer radius
 relative differences are about 3e-5 in velocity and 2e-4 in pressure (the tests allow 1e-3). The
 thermal-only case matches the analytic Chevalier & Clegg wind to 3e-6.
 
+One step differs from the paper code, for speed: the last stage of the search for the central
+density uses Brent's method instead of bisection, to the same precision (central densities agree
+to 2e-12). `outflow.boundary.find_central_density(params, method="bisection")` gives the paper
+code's search.
+
 **A note on Ė.** The published Figs. 6–9 were made with an energy injection rate
 Ė = 4.2e41 erg/s (the value in the paper's code), although the paper's Table 1 lists
 1e42 erg/s. The defaults here use 4.2e41 erg/s so that they reproduce the figures.

app/streamlit_app.py

@@ -11,8 +11,11 @@
 """
 
 import math
+import os
 import tempfile
 import traceback
+from concurrent.futures import ProcessPoolExecutor, as_completed
+from concurrent.futures.process import BrokenProcessPool
 from dataclasses import dataclass
 from pathlib import Path
 
@@ -28,6 +31,7 @@
     paper_runs,
     solve_wind,
 )
+from outflow.boundary import compiled_trials_available
 from outflow.solver import DEFAULT_OUTER_RADIUS_IN_SB
 
 TAB_NAMES = (
@@ -60,10 +64,18 @@ def solve_safely(build_params, outer_radius_in_sb: float) -> RunOutcome:
     ``build_params`` is a function returning :class:`WindParameters`, so that invalid inputs
     are reported the same way as solver failures.
     """
-    try:
+
+    def solve():
         params = build_params()
-        outer_radius = outer_radius_in_sb * params.starburst_radius
-        return RunOutcome(profile=solve_wind(params, outer_radius=outer_radius))
+        return solve_wind(params, outer_radius=outer_radius_in_sb * params.starburst_radius)
+
+    return outcome_of(solve)
+
+
+def outcome_of(compute) -> RunOutcome:
+    """Call ``compute`` (which returns a profile); turn every failure into a plain message."""
+    try:
+        return RunOutcome(profile=compute())
     except NoSolutionError as error:
         return RunOutcome(error=error.message, details=error.details)
     except ValueError as error:
@@ -469,16 +481,56 @@ class FigureRunResult:
 
 
 def run_figure(spec: FigureSpec):
-    """Recompute every run of one figure, showing "run k of n" while it works."""
+    """Recompute every run of one figure, showing "k of n finished".
+
+    With the compiled search (Numba installed) a run takes a fraction of a second, less than
+    starting a worker process, so the runs are solved here one after another. Without it they
+    are solved in parallel processes.
+    """
     runs = paper_runs(spec.number)
-    progress = st.progress(0.0, text=f"Run 1 of {len(runs)}: {runs[0].label}…")
-    results = []
-    for index, run in enumerate(runs):
-        progress.progress(index / len(runs), text=f"Run {index + 1} of {len(runs)}: {run.label}…")
-        outcome = solve_safely(lambda run=run: run.params, run.outer_radius_in_sb)
-        results.append(FigureRunResult(run, outcome, load_paper_profile(run)))
+    progress = st.progress(0.0, text=f"Solving the {len(runs)} runs…")
+
+    def show_progress(finished: int):
+        progress.progress(finished / len(runs), text=f"{finished} of {len(runs)} runs finished…")
+
+    if compiled_trials_available():
+        outcomes = solve_runs_one_by_one(runs, show_progress)
+    else:
+        outcomes = solve_runs_in_parallel(runs, show_progress)
     progress.empty()
-    st.session_state[f"fig{spec.number}_results"] = results
+    st.session_state[f"fig{spec.number}_results"] = [
+        FigureRunResult(run, outcome, load_paper_profile(run))
+        for run, outcome in zip(runs, outcomes, strict=True)
+    ]
+
+
+def solve_runs_in_parallel(runs: list[PaperRun], on_finished) -> list[RunOutcome]:
+    """Solve the runs in separate processes (one per run, up to the number of CPU cores).
+
+    The runs are independent, so this changes only the waiting time, not the results. If
+    worker processes cannot be started (some hosting platforms forbid them), the runs are
+    solved one after another instead. ``on_finished(k)`` is called as each run finishes.
+    """
+    workers = min(len(runs), os.cpu_count() or 1)
+    try:
+        with ProcessPoolExecutor(max_workers=workers) as pool:
+            futures = [
+                pool.submit(solve_wind, run.params, outer_radius=run.outer_radius) for run in runs
+            ]
+            for finished, _ in enumerate(as_completed(futures), start=1):
+                on_finished(finished)
+            return [outcome_of(future.result) for future in futures]
+    except (BrokenProcessPool, OSError, NotImplementedError):
+        return solve_runs_one_by_one(runs, on_finished)
+
+
+def solve_runs_one_by_one(runs: list[PaperRun], on_finished) -> list[RunOutcome]:
+    """Solve the runs in this process, one after another."""
+    outcomes = []
+    for finished, run in enumerate(runs, start=1):
+        outcomes.append(solve_safely(lambda run=run: run.params, run.outer_radius_in_sb))
+        on_finished(finished)
+    return outcomes
 
 
 def figure_series(results: list[FigureRunResult], show_published: bool) -> list[plots.Series]:
@@ -526,7 +578,8 @@ def figure_tab(spec: FigureSpec):
         f"Run the {len(runs)} cases of Fig. {spec.number}",
         type="primary",
         key=f"fig{spec.number}_run",
-        help="Each case takes a few seconds; the whole figure about half a minute.",
+        help="About a second with Numba installed; otherwise the cases are solved in parallel "
+        "and take a few seconds on a multi-core computer.",
     )
     if pressed:
         run_figure(spec)
@@ -636,4 +689,8 @@ def main():
         about_tab()
 
 
-main()
+# Streamlit runs this script as "__main__". The worker processes of the figure tabs load it
+# under another name (Windows and macOS start workers by re-importing the main script), and
+# this guard stops them from drawing the page or starting workers of their own.
+if __name__ == "__main__":
+    main()

docs/architecture.md

@@ -37,7 +37,8 @@ paper.load_paper_profile (published runs, package data CSV)
 | `physics.py` | single formulas of the paper (density, Alfvén speed, c\*², luminosity, radiative force and heating, NFW gravity) |
 | `equations.py` | `WindEquations`: d(v, P, P_C)/dr, inner and outer equations |
 | `integrator.py` | Cash–Karp step, error measure and step-size rule of the paper code |
-| `boundary.py` | central values (eqs. 44–45) and the search for ρ₀ |
+| `boundary.py` | central values (eqs. 44–45) and the search for ρ₀ (last stage: Brent's method by default, or the paper code's bisection) |
+| `_compiled.py` | optional Numba copy of the trial integration (physics, equations, Cash–Karp step, `_sonic_trial_python`), bit-identical, ~30× faster; loaded on first use only if Numba is installed |
 | `solver.py` | `solve_wind` and the outward integration, stall detection |
 | `profile.py` | `WindProfile`: arrays, derived quantities, interpolation, comparison, CSV read/write |
 | `paper.py` | the 14 published runs (`paper_runs`) and their saved profiles (`load_paper_profile`) |
@@ -61,6 +62,13 @@ sit on top.
 Runs are computed fresh on each button press (no caching). Explore keeps up to 8 runs of the
 session in `st.session_state` for overlay. Cooling is not offered in the app.
 
+A figure tab solves its runs one after another in the app's own process when the compiled
+search is available (a run then takes about 0.1 s, less than starting a worker process).
+Without Numba it solves them in parallel (`solve_runs_in_parallel`: a `ProcessPoolExecutor`
+with one worker per run, up to the number of CPU cores), and falls back to one after another
+if worker processes cannot be started. Workers re-import the app script, so its page code runs
+only under `if __name__ == "__main__"`.
+
 ## Other folders
 
 | Path | Role |

docs/decisions.md

@@ -32,7 +32,9 @@ minor decisions made without asking; the user can overturn them.
   thermal-only limit. (user, default accepted)
 - No command-line tool. The scriptable interface is the Python API; the UI is the only app.
   This overrides the "CLI/UI parity" app convention. (user)
-- Speed improvements are out of scope. (Claude's call; never raised as a goal)
+- Speed improvements are out of scope. (Claude's call; never raised as a goal) **Reversed
+  2026-10-15:** the user asked for speed-ups 1, 2 and 4 of the proposal (see "Speed"); Numba
+  (compiling the hot loop) is left for later.
 
 ## Physics and numerics (follow the paper code)
 
@@ -60,6 +62,8 @@ minor decisions made without asking; the user can overturn them.
     (dv/dr = q/3ρ₀, dv_A/dr = B₀/(r_sb √(4πρ₀))).
   - ρ₀ found as in `CRline.CRline()`: start at 0.1 (natural units), multiply by 10 until the
     trial solution turns over or reaches r_sb, bisect for the peak, then bisect for r_s = r_sb.
+    **Changed 2026-10-15 (user):** the last stage uses Brent's method by default (see "Speed");
+    the paper code's bisection stays available as `method="bisection"`.
   - Source-term ("inner") equations are used while v² ≤ c*² (subsonic), as the paper code does;
     density and Alfvén speed switch at r_sb (ρ: r < r_sb; v_A: r ≤ r_sb). Documented in
     `docs/physics.md` as the paper code's convention.
@@ -90,6 +94,55 @@ minor decisions made without asking; the user can overturn them.
 - The temperature is computed from each radius's own density (the `outflow.py` bug is not
   carried over).
 
+## Speed (2026-10-15, user: "implement 1, 2, 4, leave numba for later")
+
+Baseline: the 14 paper runs took 95 s in total; 96 % of a run is the ρ₀ search (81 trial
+integrations for the paper default), and half of each trial was the Runge–Kutta step's Python
+overhead. Built on branch `feature/faster-solver`:
+
+1. **Faster Runge–Kutta step, same results.** `cash_karp_step` computes its weighted sums with
+   `sum(map(mul, weights, increments))`. The products and their order are unchanged, so every
+   result is bit-for-bit the same as before (all 14 runs checked; a unit test pins it against
+   the plain loop). Writing the sums as plain `a + b + c` would be faster still (~1.6× vs
+   1.36×) but not identical: since Python 3.12, `sum()` adds floats with compensated
+   summation. (Claude's call: keep results identical.) 95 s → 70 s.
+2. **Figure tabs solve their runs in parallel**, one process per run (up to the number of CPU
+   cores), with a fallback to one after another if processes cannot be started. Results are
+   unchanged. Fig. 7 tab: 20–40 s before, 5.5 s now (measured in the real app on a 20-core
+   machine). The app script is guarded by `if __name__ == "__main__"`, because worker
+   processes re-import it.
+3. **Brent's method for the last stage of the ρ₀ search** (finding r_s = r_sb), stopping at the
+   same absolute precision as the bisection (1e-12 of the search range). The first two stages
+   (×10 scan, peak bisection) are unchanged: the peak is where the trial switches from wind to
+   breeze, a jump that suits bisection only. Measured over the 14 runs: ρ₀ differs from the
+   bisection result by ≤ 1.8e-12; profiles move by ≤ 2.1e-4 relative (the transonic solution
+   is that sensitive to ρ₀ at this precision); the match to the published runs stays within
+   1.7e-4 (was 2.1e-4). 70 s → 43.5 s. The paper code's search remains available as
+   `find_central_density(..., method="bisection")`.
+
+Overall: 95 s → 43.5 s for the 14 runs (2.2×); a default Explore run about 1.4 s.
+
+4. **Numba-compiled trial integration** (user, 2026-10-15: "try numba now"; branch
+   `experiment/numba`). `src/outflow/_compiled.py` is a line-by-line compiled copy of the trial
+   integration of the ρ₀ search (physics, equations without cooling, Cash–Karp step,
+   `_sonic_trial`). The Python modules stay the readable reference. Decisions (Claude's calls):
+   - **Bit-identical, not just close.** Python's `x**2` calls the C library's `pow`, which
+     differs from the exactly rounded `x*x` in ~0.06 % of cases (`x**3`: ~26 %), and a
+     compiler replaces constant powers by multiplications. So every exponent is read from the
+     parameter array at run time (the compiler cannot fold it), and sums repeat Python 3.12's
+     compensated `sum()`. Result: all 14 runs bit-for-bit identical; tests compare compiled and
+     Python trials on 45 inputs × 2 first-step tolerances, and all 14 full runs (slow).
+   - **Optional:** extra `fast` (`pip install -e ".[fast]"`), in `environment.yml` and in CI;
+     without Numba the Python code runs. Loaded on first use, not at import.
+   - **Only the trials are compiled** (96 % of the time). The outward integration stays in
+     Python (~0.1 s per run); cooled winds use the Python code.
+   - **App:** with the compiled search, figure tabs solve their runs one after another in the
+     app's process (worker start-up would take longer than the runs); the parallel pool is
+     kept for installs without Numba.
+   Measured: 14 runs 43 s → 1.4 s; default solve 1.37 s → 0.10 s; Fig. 7 tab 1.5 s (3.6 s on
+   the first press after the app starts, which loads Numba). Compiling takes ~5 s once after
+   installing; the compiled code is cached on disk (~0.2 s to load).
+
 ## Verification
 
 - "Same results" = differences not visible by eye. Test threshold: every reference run matches
@@ -158,7 +211,8 @@ From the builder check on these documents. All minor; the user can overturn any.
   made in this session are kept in the page's memory so they can be overlaid, and can be
   cleared.
 - **App, progress:** a progress indicator shows "run k of n" while a figure tab computes
-  (about 20–40 s locally).
+  (about 20–40 s locally). Since 2026-10-15 the runs are solved in parallel and it shows
+  "k of n runs finished" (about 5 s locally).
 - **App, axes:** as in the published figures. x is r/r_sb on a linear axis (0–5 for
   Figs. 6–8, 0–30 for Fig. 9); v in km/s, n in cm⁻³ (log), T in K (log), P in erg cm⁻³ (log).
   The Explore tab has a log-r toggle.

docs/known-issues.md

@@ -69,6 +69,29 @@ The central-density search stops after 40 increases by ×10 (`NoSolutionError`).
 integration stops after 2 000 000 steps: inside 1.1 r_sb this raises `NoSolutionError`, beyond
 it the profile is returned as stalled. Neither limit is reached by the published runs.
 
+## Profiles are sensitive to ρ₀ at the 1e-4 level
+
+Near the sonic point the solution is very sensitive to the central density: changing ρ₀ by
+about 1e-12 (the search's own precision) moves the profiles by up to 2e-4 relative. This is why
+Brent's method and the paper code's bisection, which agree on ρ₀ to 2e-12, give profiles that
+differ by up to 2e-4, and why the match to the published runs cannot be much better than
+that. Differences of this size are far below what is visible in a plot.
+
+## The compiled search is a second copy of the physics
+
+`_compiled.py` repeats the arithmetic of `physics.py`, `equations.py`, `integrator.py` and the
+trial integration. A change to any of those must be made in both places;
+`tests/unit/test_compiled.py` fails if the two stop giving identical results for the tested
+inputs. A change that never alters a result on those inputs (for example `v*v` instead of
+`v**2` in one place) can go unnoticed, so edit both copies the same way. The first run after
+installing compiles the code (about 5 s); later runs load it from Numba's on-disk cache.
+
+## Figure tabs without Numba on machines with few cores
+
+Without Numba the figure tabs solve their runs in parallel, one process per CPU core; on a
+machine or host with one or two cores (e.g. free cloud hosting) the gain is smaller. Where
+worker processes cannot be started at all, the runs are solved one after another.
+
 ## Repeated radii in the published data
 
 The published profiles repeat some radii (where the step size fell below round-off near r_sb,

docs/open-items.md

@@ -61,6 +61,9 @@ here; a limit or trap we live with goes in `known-issues.md`. Each item in one p
 
 ## Later
 
+- Compile the outward integration too (now ~0.1 s per run in Python; only worth it if many
+  runs are needed, e.g. parameter scans).
+
 - Port Figs. 1–5 scripts and the §3.5 line-emission work (Figs. 10–11).
 - Validate the experimental cooling option against `Silich.py` (Fig. 2).
 - Possible erratum on Ė in Table 1 (user's decision).