diff --git a/examples/python_fraction/.gitignore b/examples/python_fraction/.gitignore
new file mode 100644
index 000000000..71bca27c9
--- /dev/null
+++ b/examples/python_fraction/.gitignore
@@ -0,0 +1,9 @@
+# run_experiments.bash writes profiles here by default; they are large
+# (several MB per rank per repetition) and are not worth versioning.
+results/
+
+# Scratch left behind by the cylinders and by LaTeX.
+*.log
+*.aux
+*.out
+__pycache__/
diff --git a/examples/python_fraction/make_scalene_latex_table.py b/examples/python_fraction/make_scalene_latex_table.py
new file mode 100644
index 000000000..1cc80510f
--- /dev/null
+++ b/examples/python_fraction/make_scalene_latex_table.py
@@ -0,0 +1,683 @@
+#!/usr/bin/env python3
+###############################################################################
+# mpi-sppy: MPI-based Stochastic Programming in PYthon
+#
+# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for
+# Sustainable Energy, LLC, The Regents of the University of California, et al.
+# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
+# full copyright and license information.
+###############################################################################
+"""
+make_scalene_latex_table.py
+
+Generate a LaTeX summary table from Scalene MPI rank profiles.
+
+Inputs:
+  scalene_rank_0.json, scalene_rank_1.json, ...
+
+Outputs:
+  A LaTeX file containing:
+    - System information (CPU + memory) collected from standard Unix tools/files
+    - Run parameters extracted from JSON argv (best-effort)
+    - A per-rank table with:
+        Wall (s): elapsed_time_sec from JSON
+        Python (s), Native (s), System (s): wall time times the summed per-line
+          percentages taken straight out of the JSON (see scalene_totals.py)
+        Python (%): Python as a percent of the time Scalene attributed to a line
+
+Totals row:
+  - Job wall time = max wall time across ranks
+  - Sum(Python seconds), Sum(Native seconds), Sum(System seconds) across ranks
+
+Why read the JSON:
+  Earlier versions of this script parsed `scalene view --cli --reduced` and summed the
+  per-line percent columns off the screen. That is fragile: the CLI rounds to whole
+  percents, --reduced hides low-usage lines so the sums undercount, and the output now
+  carries ANSI colour codes that broke the row regex outright. The same numbers are in
+  the JSON at full precision, so that is the default. --from-cli still runs the old path
+  (with the colour codes stripped) if you want to cross-check the two.
+
+Usage:
+  python3 make_scalene_latex_table.py --glob "scalene_rank_*.json" --out scalene_summary.tex
+
+Options:
+  --from-cli      Parse `scalene view --cli` instead of reading the JSON directly
+  --cache-cli     Cache CLI output as <json>.cli.txt (reused if present)
+  --no-view       Report wall time only, no Python/native/system breakdown
+  --reduced       Pass --reduced to scalene view (only affects --from-cli)
+  --columns N     Set terminal width (COLUMNS) for scalene view output (default 200)
+"""
+
+from __future__ import annotations
+
+import argparse
+import glob
+import json
+import os
+import platform
+import re
+import subprocess
+import sys
+from dataclasses import dataclass
+from typing import Any, Dict, List, Optional, Tuple
+
+from scalene_totals import consistency_error, totals_from_json
+
+
+@dataclass
+class RankRow:
+    filename: str
+    rank: Optional[int]
+    wall_time: Optional[float]
+    python_time: Optional[float]
+    native_time: Optional[float]
+    system_time: Optional[float]
+    python_pct: Optional[float] = None
+    accounted_pct: Optional[float] = None
+
+
+# ---------------------------
+# JSON + misc helpers
+# ---------------------------
+
+def _load_json(path: str) -> Dict[str, Any]:
+    with open(path, "r", encoding="utf-8") as f:
+        return json.load(f)
+
+
+def _coerce_float(x: Any) -> Optional[float]:
+    try:
+        if x is None:
+            return None
+        return float(x)
+    except Exception:
+        return None
+
+
+def _parse_rank_from_filename(fname: str) -> Optional[int]:
+    m = re.search(r"rank[_\-]?(\d+)", os.path.basename(fname))
+    if m:
+        return int(m.group(1))
+    return None
+
+
+def _latex_escape(s: str) -> str:
+    replacements = {
+        "\\": r"\textbackslash{}",
+        "&": r"\&",
+        "%": r"\%",
+        "$": r"\$",
+        "#": r"\#",
+        "_": r"\_",
+        "{": r"\{",
+        "}": r"\}",
+        "~": r"\textasciitilde{}",
+        "^": r"\textasciicircum{}",
+    }
+    return "".join(replacements.get(ch, ch) for ch in s)
+
+
+def _get(d: Dict[str, Any], path: List[str]) -> Any:
+    cur: Any = d
+    for k in path:
+        if not isinstance(cur, dict) or k not in cur:
+            return None
+        cur = cur[k]
+    return cur
+
+
+def _first_present(d: Dict[str, Any], candidate_paths: List[List[str]]) -> Any:
+    for p in candidate_paths:
+        v = _get(d, p)
+        if v is not None:
+            return v
+    return None
+
+
+# ---------------------------
+# CLI arg extraction (from JSON)
+# ---------------------------
+
+def _extract_argv(run: Dict[str, Any]) -> Optional[List[str]]:
+    candidates = [
+        ["argv"],
+        ["commandline"],
+        ["command_line"],
+        ["cmdline"],
+        ["cmd_line"],
+        ["command"],
+        ["args"],
+        ["metadata", "argv"],
+        ["metadata", "commandline"],
+        ["metadata", "command_line"],
+        ["meta", "argv"],
+        ["meta", "commandline"],
+        ["meta", "command_line"],
+        ["header", "argv"],
+        ["header", "commandline"],
+        ["header", "command_line"],
+    ]
+    val = _first_present(run, candidates)
+    if isinstance(val, list) and all(isinstance(x, (str, int, float)) for x in val):
+        return [str(x) for x in val]
+    if isinstance(val, str):
+        return val.split()
+    return None
+
+
+def _extract_params_from_argv(argv: List[str]) -> Dict[str, str]:
+    params: Dict[str, str] = {}
+
+    def take_value(i: int) -> Optional[str]:
+        return argv[i + 1] if i + 1 < len(argv) else None
+
+    target = next((tok for tok in argv if isinstance(tok, str) and tok.endswith(".py")), None)
+    if target:
+        params["target"] = target
+
+    i = 0
+    while i < len(argv):
+        tok = argv[i]
+        if tok in (
+            "--module-name",
+            "--num-scens",
+            "--solver-name",
+            "--max-iterations",
+            "--max-solver-threads",
+            "--default-rho",
+            "--rel-gap",
+            "--outfile",
+        ):
+            v = take_value(i)
+            if v is not None:
+                params[tok.lstrip("-")] = v
+                i += 2
+                continue
+        if tok in ("--lagrangian", "--xhatshuffle"):
+            params[tok.lstrip("-")] = "true"
+            i += 1
+            continue
+        i += 1
+
+    return params
+
+
+# ---------------------------
+# System info collection (Unix assumptions)
+# ---------------------------
+
+def _run_cmd(cmd: List[str]) -> Optional[str]:
+    try:
+        p = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, text=True, check=False)
+        out = p.stdout.strip()
+        return out if out else None
+    except Exception:
+        return None
+
+
+def _read_first_existing(paths: List[str]) -> Optional[str]:
+    for p in paths:
+        try:
+            with open(p, "r", encoding="utf-8") as f:
+                return f.read()
+        except Exception:
+            continue
+    return None
+
+
+def _parse_meminfo_kib(meminfo_text: str) -> Dict[str, int]:
+    """
+    Return selected values from /proc/meminfo in KiB
+    """
+    out: Dict[str, int] = {}
+    for line in meminfo_text.splitlines():
+        m = re.match(r"^(\w+):\s+(\d+)\s+kB\s*$", line)
+        if m:
+            out[m.group(1)] = int(m.group(2))
+    return out
+
+
+def _format_bytes(n: Optional[int]) -> str:
+    if n is None:
+        return "unknown"
+    # n is bytes
+    units = ["B", "KiB", "MiB", "GiB", "TiB"]
+    x = float(n)
+    i = 0
+    while x >= 1024.0 and i < len(units) - 1:
+        x /= 1024.0
+        i += 1
+    if i == 0:
+        return f"{int(x)} {units[i]}"
+    return f"{x:.2f} {units[i]}"
+
+
+def _collect_system_info() -> Dict[str, str]:
+    """
+    Best-effort system inventory on Unix:
+      - OS/kernel
+      - CPU count (logical + physical if available)
+      - CPU model
+      - CPU max MHz (or base)
+      - Total memory and (approx) available memory at report time
+
+    Works on Linux best; degrades gracefully on macOS/others.
+    """
+    info: Dict[str, str] = {}
+
+    # OS / kernel
+    info["os"] = f"{platform.system()} {platform.release()} ({platform.machine()})"
+
+    # CPU counts
+    logical = os.cpu_count()
+    info["cpu_logical"] = str(logical) if logical is not None else "unknown"
+
+    # lscpu (Linux)
+    lscpu = _run_cmd(["lscpu"])
+    cpu_model = None
+    cpu_mhz = None
+    cpu_max_mhz = None
+    cpu_sockets = None
+    cores_per_socket = None
+    threads_per_core = None
+
+    if lscpu:
+        for line in lscpu.splitlines():
+            if ":" not in line:
+                continue
+            k, v = [x.strip() for x in line.split(":", 1)]
+            kl = k.lower()
+            if kl == "model name":
+                cpu_model = v
+            elif kl in ("cpu mhz",):
+                cpu_mhz = v
+            elif kl in ("cpu max mhz",):
+                cpu_max_mhz = v
+            elif kl == "socket(s)":
+                cpu_sockets = v
+            elif kl == "core(s) per socket":
+                cores_per_socket = v
+            elif kl == "thread(s) per core":
+                threads_per_core = v
+
+    # sysctl (macOS / BSD)
+    if cpu_model is None:
+        cpu_model = _run_cmd(["sysctl", "-n", "machdep.cpu.brand_string"]) or None
+
+    # Physical cores (best-effort)
+    physical_cores = None
+    if cpu_sockets and cores_per_socket:
+        try:
+            physical_cores = int(cpu_sockets) * int(cores_per_socket)
+        except Exception:
+            physical_cores = None
+    if physical_cores is None:
+        # macOS
+        pc = _run_cmd(["sysctl", "-n", "hw.physicalcpu"])
+        if pc and pc.isdigit():
+            physical_cores = int(pc)
+
+    if physical_cores is not None:
+        info["cpu_physical_cores"] = str(physical_cores)
+
+    if cpu_model:
+        info["cpu_model"] = cpu_model
+
+    # Frequency
+    # Prefer max MHz if available
+    freq = cpu_max_mhz or cpu_mhz
+    if freq:
+        info["cpu_mhz"] = freq
+
+    if threads_per_core:
+        info["threads_per_core"] = threads_per_core
+    if cpu_sockets:
+        info["cpu_sockets"] = cpu_sockets
+    if cores_per_socket:
+        info["cores_per_socket"] = cores_per_socket
+
+    # Memory: Linux /proc/meminfo or macOS sysctl/vm_stat
+    meminfo = _read_first_existing(["/proc/meminfo"])
+    if meminfo:
+        m = _parse_meminfo_kib(meminfo)
+        mem_total_bytes = m.get("MemTotal", 0) * 1024 if "MemTotal" in m else None
+        mem_avail_bytes = m.get("MemAvailable", 0) * 1024 if "MemAvailable" in m else None
+        info["mem_total"] = _format_bytes(mem_total_bytes)
+        info["mem_available"] = _format_bytes(mem_avail_bytes)
+    else:
+        # macOS total
+        mt = _run_cmd(["sysctl", "-n", "hw.memsize"])
+        if mt and mt.isdigit():
+            info["mem_total"] = _format_bytes(int(mt))
+        # macOS available is trickier; best-effort via vm_stat
+        vm = _run_cmd(["vm_stat"])
+        if vm:
+            # Parse page size and free/inactive/speculative, etc.
+            page_size = 4096
+            mps = re.search(r"page size of (\d+) bytes", vm)
+            if mps:
+                page_size = int(mps.group(1))
+            counts = {}
+            for line in vm.splitlines():
+                mm = re.match(r"^([^:]+):\s+(\d+)\.", line.strip())
+                if mm:
+                    counts[mm.group(1).strip()] = int(mm.group(2))
+            # rough estimate: free + inactive + speculative
+            avail_pages = (
+                counts.get("Pages free", 0)
+                + counts.get("Pages inactive", 0)
+                + counts.get("Pages speculative", 0)
+            )
+            info["mem_available"] = _format_bytes(avail_pages * page_size)
+
+    return info
+
+
+# ---------------------------
+# Scalene view parsing
+# ---------------------------
+
+def _run_scalene_view_cli(json_path: str, reduced: bool, columns: int) -> str:
+    cmd = [sys.executable, "-m", "scalene", "view", "--cli"]
+    if reduced:
+        cmd.append("--reduced")
+    cmd.append(json_path)
+
+    env = dict(os.environ)
+    env["COLUMNS"] = str(columns)
+
+    p = subprocess.run(
+        cmd,
+        stdout=subprocess.PIPE,
+        stderr=subprocess.STDOUT,
+        text=True,
+        env=env,
+        check=False,
+    )
+    return p.stdout
+
+
+_ANSI_RE = re.compile(r"\x1b\[[0-9;]*[A-Za-z]")
+
+
+def _strip_ansi(s: str) -> str:
+    return _ANSI_RE.sub("", s)
+
+
+def _parse_cli_percent_totals(cli_text: str) -> Tuple[Optional[float], Optional[float], Optional[float]]:
+    """
+    Parse totals by summing per-line percent columns from Scalene `view --cli` output.
+
+    We match table rows like:
+       15 │    2% │   24% │   2%  │ ...
+
+    Scalene colourizes this table even when its output is a pipe, so the escape
+    sequences have to come out before the row pattern will match anything.
+
+    Returns: (python_percent, native_percent, system_percent)
+    """
+    py_pct = 0.0
+    nat_pct = 0.0
+    sys_pct = 0.0
+    saw_any = False
+
+    row_re = re.compile(
+        r"^\s*\d+\s*│\s*([0-9]+(?:\.[0-9]+)?)?\s*%?\s*│\s*([0-9]+(?:\.[0-9]+)?)?\s*%?\s*│\s*([0-9]+(?:\.[0-9]+)?)?\s*%?\s*│"
+    )
+
+    for line in cli_text.splitlines():
+        m = row_re.match(_strip_ansi(line))
+        if not m:
+            continue
+        saw_any = True
+        a, b, c = m.group(1), m.group(2), m.group(3)
+        py_pct += float(a) if a else 0.0
+        nat_pct += float(b) if b else 0.0
+        sys_pct += float(c) if c else 0.0
+
+    if not saw_any:
+        return None, None, None
+
+    return py_pct, nat_pct, sys_pct
+
+
+def _fmt(x: Optional[float], digits: int = 2, na: str = r"\textemdash") -> str:
+    if x is None:
+        return na
+    return f"{x:.{digits}f}"
+
+
+# ---------------------------
+# Main
+# ---------------------------
+
+def main() -> None:
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--glob", default="scalene_rank_*.json", help="Glob for Scalene JSON files")
+    ap.add_argument("--out", default="scalene_summary.tex", help="Output LaTeX filename")
+    ap.add_argument("--caption", default="Scalene timing summary by MPI rank", help="Table caption")
+    ap.add_argument("--label", default="tab:scalene-summary", help="LaTeX label")
+    ap.add_argument("--no-totals", action="store_true", help="Do not add totals row")
+    ap.add_argument("--no-view", action="store_true", help="Report wall time only, no breakdown")
+    ap.add_argument("--from-cli", action="store_true",
+                    help="Parse `scalene view --cli` instead of reading the JSON directly")
+    ap.add_argument("--cache-cli", action="store_true", help="Cache CLI output to <json>.cli.txt and reuse")
+    ap.add_argument("--reduced", action="store_true", help="Pass --reduced to scalene view --cli (--from-cli only)")
+    ap.add_argument("--columns", type=int, default=200, help="Set COLUMNS for scalene view output (default 200)")
+    args = ap.parse_args()
+
+    files = sorted(glob.glob(args.glob))
+    if not files:
+        raise SystemExit(f"No files matched glob: {args.glob}")
+
+    # System info
+    sysinfo = _collect_system_info()
+
+    # Load JSON profiles
+    runs = [(f, _load_json(f)) for f in files]
+
+    # Use the first JSON file as the "run metadata" source
+    first_json = runs[0][1]
+    argv = _extract_argv(first_json)
+    params: Dict[str, str] = _extract_params_from_argv(argv) if argv else {}
+
+    rows: List[RankRow] = []
+    for f, j in runs:
+        wall = _coerce_float(j.get("elapsed_time_sec"))
+
+        python_s = native_s = system_s = None
+        python_pct = accounted_pct = None
+
+        if not args.no_view and wall is not None:
+            if args.from_cli:
+                cache_path = f"{f}.cli.txt"
+                if args.cache_cli and os.path.exists(cache_path):
+                    with open(cache_path, "r", encoding="utf-8") as cf:
+                        cli_text = cf.read()
+                else:
+                    cli_text = _run_scalene_view_cli(f, reduced=args.reduced, columns=args.columns)
+                    if args.cache_cli:
+                        with open(cache_path, "w", encoding="utf-8") as cf:
+                            cf.write(cli_text)
+
+                py_pct, nat_pct, sys_pct = _parse_cli_percent_totals(cli_text)
+                if py_pct is not None and nat_pct is not None and sys_pct is not None:
+                    python_s = wall * (py_pct / 100.0)
+                    native_s = wall * (nat_pct / 100.0)
+                    system_s = wall * (sys_pct / 100.0)
+                    accounted_pct = py_pct + nat_pct + sys_pct
+                    if accounted_pct > 0.0:
+                        python_pct = 100.0 * py_pct / accounted_pct
+            else:
+                bad = consistency_error(f)
+                if bad:
+                    raise SystemExit(
+                        "Scalene JSON failed its internal consistency check, so its "
+                        "layout has probably changed:\n  " + bad
+                    )
+                t = totals_from_json(f)
+                python_s, native_s, system_s = t.python_sec, t.native_sec, t.system_sec
+                python_pct = t.python_fraction
+                accounted_pct = t.accounted_pct
+
+        rows.append(
+            RankRow(
+                filename=os.path.basename(f),
+                rank=_parse_rank_from_filename(f),
+                wall_time=wall,
+                python_time=python_s,
+                native_time=native_s,
+                system_time=system_s,
+                python_pct=python_pct,
+                accounted_pct=accounted_pct,
+            )
+        )
+
+    rows.sort(key=lambda r: (999999 if r.rank is None else r.rank, r.filename))
+
+    vals_wall = [r.wall_time for r in rows if r.wall_time is not None]
+    vals_py = [r.python_time for r in rows if r.python_time is not None]
+    vals_nat = [r.native_time for r in rows if r.native_time is not None]
+    vals_sys = [r.system_time for r in rows if r.system_time is not None]
+
+    job_wall = max(vals_wall) if vals_wall else None
+    sum_py = sum(vals_py) if vals_py else None
+    sum_nat = sum(vals_nat) if vals_nat else None
+    sum_sys = sum(vals_sys) if vals_sys else None
+
+    any_time_breakdown = bool(vals_py or vals_nat or vals_sys)
+
+    # Build LaTeX
+    lines: List[str] = []
+    lines.append("% Auto-generated by make_scalene_latex_table.py")
+    lines.append("")
+
+    # System info block
+    lines.append(r"\noindent\textbf{System information:}\\")
+    lines.append(r"\begin{itemize}")
+    if "os" in sysinfo:
+        lines.append(rf"  \item \texttt{{os}}={{{_latex_escape(sysinfo['os'])}}}")
+    if "cpu_model" in sysinfo:
+        lines.append(rf"  \item \texttt{{cpu\_model}}={{{_latex_escape(sysinfo['cpu_model'])}}}")
+    if "cpu_logical" in sysinfo:
+        lines.append(rf"  \item \texttt{{cpu\_logical}}={{{_latex_escape(sysinfo['cpu_logical'])}}}")
+    if "cpu_physical_cores" in sysinfo:
+        lines.append(rf"  \item \texttt{{cpu\_physical\_cores}}={{{_latex_escape(sysinfo['cpu_physical_cores'])}}}")
+    if "cpu_mhz" in sysinfo:
+        lines.append(rf"  \item \texttt{{cpu\_mhz}}={{{_latex_escape(sysinfo['cpu_mhz'])}}}")
+    if "cpu_sockets" in sysinfo:
+        lines.append(rf"  \item \texttt{{cpu\_sockets}}={{{_latex_escape(sysinfo['cpu_sockets'])}}}")
+    if "cores_per_socket" in sysinfo:
+        lines.append(rf"  \item \texttt{{cores\_per\_socket}}={{{_latex_escape(sysinfo['cores_per_socket'])}}}")
+    if "threads_per_core" in sysinfo:
+        lines.append(rf"  \item \texttt{{threads\_per\_core}}={{{_latex_escape(sysinfo['threads_per_core'])}}}")
+    if "mem_total" in sysinfo:
+        lines.append(rf"  \item \texttt{{mem\_total}}={{{_latex_escape(sysinfo['mem_total'])}}}")
+    if "mem_available" in sysinfo:
+        lines.append(rf"  \item \texttt{{mem\_available}}={{{_latex_escape(sysinfo['mem_available'])}}}")
+    lines.append(r"\end{itemize}")
+    lines.append("")
+
+    # Run parameters comment header
+    lines.append("% Run parameters extracted from JSON (best-effort):")
+    if argv:
+        lines.append(f"% argv: {_latex_escape(' '.join(argv))}")
+    else:
+        lines.append("% argv: (not found in JSON)")
+    if params:
+        lines.append("% Parsed parameters:")
+        for k in sorted(params.keys()):
+            lines.append(f"%   {k}: {_latex_escape(params[k])}")
+    lines.append("")
+
+    # Run parameters block
+    lines.append(r"\noindent\textbf{Run parameters (from Scalene JSON):}\\")
+    if params:
+        show_keys = [
+            "target",
+            "module-name",
+            "num-scens",
+            "solver-name",
+            "max-iterations",
+            "max-solver-threads",
+            "default-rho",
+            "rel-gap",
+            "lagrangian",
+            "xhatshuffle",
+        ]
+        parts = []
+        for k in show_keys:
+            if k in params:
+                parts.append(rf"\texttt{{{_latex_escape(k)}}}={{{_latex_escape(params[k])}}}")
+        lines.append(r"\begin{itemize}")
+        for p in parts:
+            lines.append(rf"  \item {p}")
+        lines.append(r"\end{itemize}")
+    else:
+        lines.append(r"\emph{(Command line not found in JSON.)}\\")
+    lines.append("")
+
+    if args.no_view:
+        lines.append(
+            r"\noindent\emph{Note: --no-view was given, so only wall time is reported.}"
+        )
+        lines.append("")
+    elif not any_time_breakdown:
+        lines.append(
+            r"\noindent\emph{Note: No per-line time percentages were found in the profiles.}"
+        )
+        lines.append("")
+
+    vals_acct = [r.accounted_pct for r in rows if r.accounted_pct is not None]
+    if vals_acct:
+        lines.append(
+            r"\noindent\emph{Scalene attributed "
+            rf"{min(vals_acct):.1f}--{max(vals_acct):.1f}\% "
+            r"of wall time to a source line; the Python (\%) column is Python as a "
+            r"percent of that attributed time.}"
+        )
+        lines.append("")
+
+    # Table
+    lines.append(r"\begin{table}[ht]")
+    lines.append(r"\centering")
+    lines.append(r"\begin{tabular}{r l r r r r r}")
+    lines.append(r"\hline")
+    lines.append(r"Rank & File & Wall (s) & Python (s) & Native (s) & System (s) & Python (\%) \\")
+    lines.append(r"\hline")
+
+    for r in rows:
+        rank_str = "" if r.rank is None else str(r.rank)
+        lines.append(
+            rf"{rank_str} & {_latex_escape(r.filename)} & {_fmt(r.wall_time)} & {_fmt(r.python_time)} & {_fmt(r.native_time)} & {_fmt(r.system_time)} & {_fmt(r.python_pct, 1)} \\"
+        )
+
+    if not args.no_totals:
+        # The job-level Python percent is computed from the summed seconds, so
+        # that ranks are weighted by how long they actually ran.
+        job_py_pct = None
+        if sum_py is not None and sum_nat is not None and sum_sys is not None:
+            denom = sum_py + sum_nat + sum_sys
+            if denom > 0.0:
+                job_py_pct = 100.0 * sum_py / denom
+        lines.append(r"\hline")
+        lines.append(
+            rf"\textbf{{Job wall (max)}} &  & \textbf{{{_fmt(job_wall)}}} & \textbf{{{_fmt(sum_py)}}} & \textbf{{{_fmt(sum_nat)}}} & \textbf{{{_fmt(sum_sys)}}} & \textbf{{{_fmt(job_py_pct, 1)}}} \\"
+        )
+
+    lines.append(r"\hline")
+    lines.append(r"\end{tabular}")
+    lines.append(rf"\caption{{{_latex_escape(args.caption)}}}")
+    lines.append(rf"\label{{{_latex_escape(args.label)}}}")
+    lines.append(r"\end{table}")
+    lines.append("")
+
+    with open(args.out, "w", encoding="utf-8") as f:
+        f.write("\n".join(lines))
+
+    print(f"Wrote LaTeX to: {args.out}")
+    print(f"Read {len(files)} JSON files matched by: {args.glob}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/python_fraction/make_tables.bash b/examples/python_fraction/make_tables.bash
new file mode 100755
index 000000000..7f9ff1e30
--- /dev/null
+++ b/examples/python_fraction/make_tables.bash
@@ -0,0 +1,48 @@
+#!/bin/bash
+###############################################################################
+# mpi-sppy: MPI-based Stochastic Programming in PYthon
+#
+# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for
+# Sustainable Energy, LLC, The Regents of the University of California, et al.
+# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
+# full copyright and license information.
+###############################################################################
+
+# Regenerate the LaTeX tables in python_fraction.tex from the profiles that
+# run_experiments.bash wrote. Cheap to re-run; does not re-run any experiment.
+#
+# Usage:
+#   ./make_tables.bash [results_dir]
+
+set -euo pipefail
+
+HERE="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+RESULTS="${1:-${HERE}/results}"
+
+# Rank order comes from generic_cylinders: the hub is rank 0 and the spokes
+# follow in the order build_spoke_list appends them, which for --lagrangian
+# --xhatshuffle is lagrangian then xhatshuffle.
+LABELS="PH hub,lagrangian,xhatshuffle"
+
+CASES=(farmer3 farmer60 farmer240 farmer240_bun10
+       sslp_15_45_10 sslp_15_45_15 sslp_15_45_15_bun3 sslp_5_25_50)
+
+# Primary results: the persistent interface.
+python3 "${HERE}/summarize_reps.py" \
+  --results "${RESULTS}" \
+  --solvers gurobi_persistent \
+  --cases "${CASES[@]}" \
+  --rank-labels "${LABELS}" \
+  --out "${HERE}/scalene_summary_persistent.tex"
+
+# Secondary: the file-based interface, for the contrast noted in the writeup.
+if [[ -d "${RESULTS}/gurobi" ]]; then
+  python3 "${HERE}/summarize_reps.py" \
+    --results "${RESULTS}" \
+    --solvers gurobi \
+    --cases "${CASES[@]}" \
+    --rank-labels "${LABELS}" \
+    --out "${HERE}/scalene_summary_file_interface.tex"
+fi
+
+echo "Tables written to ${HERE}"
diff --git a/examples/python_fraction/python_fraction.pdf b/examples/python_fraction/python_fraction.pdf
new file mode 100644
index 000000000..4cb369d0e
Binary files /dev/null and b/examples/python_fraction/python_fraction.pdf differ
diff --git a/examples/python_fraction/python_fraction.tex b/examples/python_fraction/python_fraction.tex
new file mode 100644
index 000000000..e12400c4f
--- /dev/null
+++ b/examples/python_fraction/python_fraction.tex
@@ -0,0 +1,249 @@
+\documentclass{article}
+
+% Language setting
+% Replace `english' with e.g. `spanish' to change the document language
+\usepackage[english]{babel}
+
+% Set page size and margins
+% Replace `letterpaper' with`a4paper' for UK/EU standard size
+\usepackage[letterpaper,top=2cm,bottom=2cm,left=3cm,right=3cm,marginparwidth=1.75cm]{geometry}
+
+% Useful packages
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{amsthm}
+\usepackage{amsfonts}
+\usepackage{graphicx}
+\usepackage{algorithm2e}
+\usepackage{enumitem}
+\usepackage{comment}
+\usepackage{natbib}
+
+\usepackage[colorlinks=true, allcolors=blue]{hyperref}
+
+\newcommand{\sdag}[1]{{#1}^{\dag}}
+
+\title{Fraction of Time MPI-SPPY spends in Python}
+\author{ David L Woodruff\\
+  Graduate School of Management\\
+  \\
+  University of California Davis\\
+  Davis CA 95616 USA}
+\date{\today}
+
+\newtheorem{theorem}{Theorem}
+\newtheorem{lemma}{Lemma}
+
+\begin{document}
+\maketitle
+
+When considering applications in practice, or when comparing to other
+packages, a question arises concerning the fraction of time that
+mpi-sppy spends ``in Python'' as opposed to compiled code written in
+other languages such as C and Fortran.  Since solvers, numpy, and MPI
+are all in the latter category, {\em a priori} one expects that the
+fraction spent in Python will be small.
+
+We have used the tool called {\em scalene}, which is designed to
+attribute time to Python vs native (it can estimate time spent in
+compiled code called from Python, which it calls ``native''). It works
+by sampling.
+
+The short answer is: for hard problems, not much.
+The longer answer is that there is no single number: the fraction
+depends almost entirely on how much work the solver does per
+subproblem. When the subproblems are even somewhat substantial, as in the
+\texttt{sslp\_15\_45} instances, about 9\% of the time is spent in
+Python, which matches the {\em a priori} expectation. For harder subproblems,
+the fraction can be smaller. When the
+subproblems are small enough that the solver returns almost
+immediately, as in farmer, the majority of the time is spent in Python
+--- 76\% in the most extreme case measured here. The crossover is not
+subtle and it is not a sampling artifact.
+
+The practical consequences are:
+
+\begin{itemize}
+  \item If your subproblems keep the solver busy, you cannot get much
+    speed-up just by using another language. Maybe another language
+    will make it easier for you to do your work.
+  \item If your subproblems are small, the per-solve Python cost is
+    what you are paying for, and the fix is to give the solver more
+    work per call rather than to rewrite anything. Bundling is the
+    obvious lever, and it is examined below.
+\end{itemize}
+
+\section{Method}
+
+Every case runs the same three cylinders --- a PH hub, a lagrangian
+spoke and an xhatshuffle spoke, one MPI rank each --- so that only the
+model and the instance change between cases. The convergence tolerance
+is set to zero so that runs end on the iteration limit rather than on a
+gap, which keeps the amount of work per case predictable.
+
+Because scalene samples, every case is run three times and the tables
+report the mean over repetitions with the observed min--max range. The
+spread turns out to be small for runs of a minute or so (typically a
+few tenths of a percentage point) and much larger for the deliberately
+short \texttt{farmer3} case, which is why that case is kept: it shows
+what the numbers look like when there is not enough run time to sample
+properly.
+
+The numbers are read out of scalene's JSON profile rather than scraped
+from the output of \texttt{scalene view --cli}. Summing the per-line
+percentages in the JSON gives the same quantities at full precision,
+and avoids three problems with scraping the terminal output: it rounds
+each line to a whole percent, its \texttt{--reduced} form omits
+low-usage lines so the sums undercount, and it is colourized, which
+silently broke the row-matching in the earlier version of this code.
+Scalene attributes 91--99\% of wall time to some source line; the
+percentages below are shares of that attributed time, so the
+unattributed remainder is divided out rather than being silently
+credited to one of the buckets.
+
+All runs are on a single machine with homogeneous cores, which avoids a
+difficulty that affected an earlier version of these experiments: on a
+machine with a mixture of fast and slow cores, ranks landing on slow
+cores showed inflated Python fractions.
+
+\noindent\textbf{System and software:}\\
+\begin{itemize}
+  \item \texttt{os}={Linux 7.0.0-28-generic (x86\_64)}
+  \item \texttt{cpu\_model}={AMD Ryzen 7 7840HS}
+  \item \texttt{cpu\_physical\_cores}={8}, \texttt{cpu\_logical}={16},
+        \texttt{threads\_per\_core}={2}
+  \item \texttt{cpu\_max\_mhz}={5137.9}
+  \item \texttt{mem\_total}={14.31 GiB}
+  \item \texttt{python}={3.11.13}, \texttt{scalene}={2.0.1},
+        \texttt{pyomo}={6.9.5.dev0}, \texttt{mpi4py}={4.1.0},
+        \texttt{gurobipy}={13.0.2}
+  \item solver interface: \texttt{gurobi\_persistent},
+        \texttt{--max-solver-threads 2}
+\end{itemize}
+
+\section{Results}
+
+Table~\ref{tab:python-fraction-summary} gives the split by case and
+Table~\ref{tab:python-fraction-by-rank} breaks the Python percentage
+out by cylinder. The cases are:
+
+\begin{itemize}
+  \item \texttt{farmer3} --- 3 scenarios, 100 iterations. Deliberately
+    short: about two seconds unprofiled, five under scalene.
+  \item \texttt{farmer60}, \texttt{farmer240} --- 60 and 240
+    scenarios, 400 and 200 iterations. Many trivially small LP
+    subproblems.
+  \item \texttt{farmer240\_bun10} --- \texttt{farmer240} with ten
+    scenarios gathered into each proper bundle, so 24 subproblems
+    instead of 240.
+  \item \texttt{sslp\_15\_45\_10}, \texttt{sslp\_15\_45\_15} --- 10 and
+    15 scenarios, 50 iterations. MIP subproblems that keep the solver
+    busy.
+  \item \texttt{sslp\_15\_45\_15\_bun3} ---
+    \texttt{sslp\_15\_45\_15} with three scenarios per bundle, so 5
+    subproblems instead of 15.
+  \item \texttt{sslp\_5\_25\_50} --- 50 scenarios, 150 iterations. Many
+    small MIP subproblems.
+\end{itemize}
+
+\input{scalene_summary_persistent}
+
+The ordering is by subproblem difficulty, not by instance size. The two
+\texttt{sslp\_15\_45} cases sit at 8.7\% and 9.6\% Python; every case
+whose subproblems are small sits far above that, from 58\% for
+\texttt{sslp\_5\_25\_50} up to 76\% for \texttt{farmer240}. Note that
+\texttt{sslp\_5\_25\_50} has more scenarios than
+\texttt{sslp\_15\_45\_15} and yet spends six times the fraction of its
+time in Python: what matters is the size of each subproblem, not how
+many there are.
+
+Looking at where the samples land makes the mechanism explicit. In both
+regimes essentially all of the time is attributed to one line,
+\texttt{spopt.py:337}, which is the Pyomo \texttt{solve} call. For
+\texttt{sslp\_15\_45\_15} that line is 6.6\% Python and 88.2\% native
+--- the solver is doing the work. For \texttt{farmer60} the same line
+is 48.2\% Python and 5.8\% native, and a further 9.4\% of the run is
+Python time in \texttt{spopt.py:267}, where the proximal objective is
+handed to the solver on each iteration. That objective work scales with
+the size of the model rather than with the difficulty of the solve,
+which is exactly why it is invisible for sslp and dominant for farmer.
+
+The per-cylinder table shows the PH hub and the lagrangian spoke
+agreeing closely within every case, as they should, since they solve
+the same subproblems with different objectives. The xhatshuffle spoke
+is consistently the most Python-heavy of the three, and the gap widens
+as the hub's own Python share falls: 6.8\% against 13.0\% for
+\texttt{sslp\_15\_45\_10}, and 5.1\% against 25.3\% for
+\texttt{sslp\_15\_45\_15\_bun3}. This is expected rather than
+anomalous. The xhatshuffle spoke is looking for a feasible incumbent
+rather than solving the subproblem to optimality, so much of its work
+is Python bookkeeping --- shuffling the scenario order, fixing
+candidate nonanticipative values, and checking the result --- with
+comparatively little time left inside the solver. It is a small
+fraction of total work, since it is one rank of three, but it is worth
+knowing that the cylinders are not interchangeable for this
+measurement.
+
+\subsection{Bundles}
+
+Bundling gives the solver more work per call, so it is the natural
+response to a high Python fraction. It helps, but not reliably, and it
+is worth being clear about why.
+
+For farmer, bundling ten scenarios per subproblem takes the Python
+share from 76.1\% down to 62.6\% and cuts wall time by a factor of
+four, from 82.9 to 20.9 seconds. Fewer, larger solves means fewer
+round trips through Pyomo, and the wall-clock gain is large. The Python
+share nonetheless stays high, because a bundle of ten farmer scenarios
+is still a trivial LP.
+
+For \texttt{sslp\_15\_45\_15}, bundling three scenarios per subproblem
+moves the Python share the other way, from 9.6\% up to 12.5\%, while
+also cutting wall time from 104.8 to 73.7 seconds. Both numbers moved
+in the direction that makes sense once the effect is separated into its
+two parts: bundling reduced the total solver work more than it reduced
+the total Python work, so Python's {\em share} rose even though the run
+got faster.
+
+So bundling should be judged on wall time, where it won in both cases,
+and not on the Python fraction, which it can move either way. A high
+Python fraction is a symptom worth investigating, but driving it down
+is not itself the objective.
+
+\subsection{Caveats}
+
+Scalene's instrumentation is not free, and its cost falls mostly on
+Python, so these Python percentages are upper bounds. The overhead
+column in Table~\ref{tab:python-fraction-summary} is the ratio of
+profiled to unprofiled wall time for the same case, and it lines up
+with the Python column: the solver-bound sslp cases run at
+0.99--1.11$\times$, essentially unaffected, while the Python-heavy
+farmer cases run at 1.42--1.55$\times$. So the true Python fractions
+for the farmer cases are somewhat lower than the table reports. The
+exception is \texttt{farmer3} at 2.36$\times$, which is not a
+Python-share effect at all: the profiler's fixed startup cost is simply
+large next to a two-second run, which is another reason not to trust
+that row for anything but its variability.
+
+The gap between the two regimes is far too large to be an artifact of
+this overhead --- 9\% against 76\% will not be closed by a factor of
+1.5 --- but the individual percentages should be read as approximate.
+
+Two further limits are worth stating. These are Pyomo models and the
+Pyomo time is included in the Python total, so a good deal of what is
+labelled Python here is Pyomo rather than mpi-sppy. And all of this is
+one machine, one solver, and three repetitions per case.
+
+\subsection{A note on non-persistent solvers}
+
+All of the above uses \texttt{gurobi\_persistent}. If you use a
+non-persistent interface instead --- \texttt{--solver-name gurobi} is
+Pyomo's file-based interface, which writes an LP file and parses a
+solution file on every solve --- then Pyomo does a great deal of extra
+Python work per solve, which matters most when the subproblems are
+small: \texttt{sslp\_5\_25\_50} goes from 58.1\% to 75.4\% Python and
+from 74 to 124 seconds, and even the solver-bound
+\texttt{sslp\_15\_45\_10}, whose wall time barely moves, goes from
+8.7\% to 31.8\% Python. Use a persistent interface if you have one.
+
+\end{document}
diff --git a/examples/python_fraction/readme.rst b/examples/python_fraction/readme.rst
new file mode 100644
index 000000000..8f250dc53
--- /dev/null
+++ b/examples/python_fraction/readme.rst
@@ -0,0 +1,72 @@
+Fraction of time in python
+==========================
+
+See ``python_fraction.tex`` for the writeup and the results.
+
+Files
+-----
+
+``run_experiments.bash``
+    Runs every case under scalene, repeating each case so that the writeup can
+    show run-to-run spread instead of a single sample. Repetitions matter here:
+    scalene works by sampling, so one run of one case is weak evidence.
+
+``scalene_wrapper.bash``
+    Rank-aware scalene launcher; ``mpiexec`` runs this rather than python
+    directly, so that each rank can name its own output file. Not normally run
+    by hand.
+
+``summarize_reps.py``
+    Aggregates the repetitions into the LaTeX tables.
+
+``make_tables.bash``
+    Regenerates the tables from profiles already on disk. No experiments rerun.
+
+``make_scalene_latex_table.py``
+    Detail table for a single run (one directory of per-rank profiles).
+
+``scalene_totals.py``
+    Pulls the Python/native/system split out of a scalene JSON profile. Shared
+    by the two table generators.
+
+Running
+-------
+
+Profile the default case list with the persistent solver interface, three
+repetitions each::
+
+    $ SOLVER=gurobi_persistent ./run_experiments.bash
+
+Then, to get the unprofiled wall times that the overhead column needs::
+
+    $ SOLVER=gurobi_persistent PROFILE=0 ./run_experiments.bash
+
+And regenerate the tables::
+
+    $ ./make_tables.bash
+
+Individual cases can be named on the command line, e.g.
+``./run_experiments.bash farmer60 sslp_5_25_50``.
+
+Notes
+-----
+
+Prefer a persistent solver interface. ``--solver-name gurobi`` is Pyomo's
+file-based interface, which writes an LP file and parses a solution file in
+Python on every solve; on small subproblems that Python work dominates
+everything else, and the measured Python fraction then says more about Pyomo's
+file writer than about mpi-sppy.
+
+The numbers come from the scalene JSON rather than from scraping
+``scalene view --cli``. Reading the JSON avoids three problems with scraping:
+the CLI rounds each line to a whole percent, ``--reduced`` hides low-usage
+lines so the sums undercount, and the output is colourized, which silently
+broke the original row-matching regex. ``make_scalene_latex_table.py --from-cli``
+still parses the CLI if you want to compare the two paths.
+
+Scalene occasionally dies during startup with a ``KeyError`` from inside
+``importlib``, before any mpi-sppy code runs; ``run_experiments.bash`` retries a
+repetition that comes back with fewer profiles than ranks.
+
+This code suite is probably fragile because scalene seems to do major updates
+that change the output format. It was last run against scalene 2.0.1.
diff --git a/examples/python_fraction/run_experiments.bash b/examples/python_fraction/run_experiments.bash
new file mode 100755
index 000000000..4c06f2959
--- /dev/null
+++ b/examples/python_fraction/run_experiments.bash
@@ -0,0 +1,188 @@
+#!/bin/bash
+###############################################################################
+# mpi-sppy: MPI-based Stochastic Programming in PYthon
+#
+# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for
+# Sustainable Energy, LLC, The Regents of the University of California, et al.
+# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
+# full copyright and license information.
+###############################################################################
+
+# Run the "fraction of time in Python" experiments under scalene.
+#
+# Every case runs the same three cylinders (PH hub, lagrangian, xhatshuffle) so
+# that only the model and the instance size change. Each case is repeated REPS
+# times so that the report can show run-to-run spread rather than a single
+# sample; scalene works by sampling, so a single run says very little.
+#
+# Usage:
+#   ./run_experiments.bash                 # all cases, REPS reps each
+#   REPS=1 ./run_experiments.bash farmer3  # just one case, one rep
+#
+# Environment:
+#   SOLVER   solver name (default gurobi)
+#   REPS     repetitions per case (default 3)
+#   NP       number of MPI ranks, i.e. cylinders (default 3)
+#   THREADS  --max-solver-threads (default 2)
+#   RESULTS  output directory (default ./results)
+#   TRIES    attempts per rep before giving up (default 3); see the retry note below
+#   PROFILE  1 (default) to run under scalene; 0 to run the identical cases with
+#            no profiler and record only wall time, in <outdir>/wall.txt. The
+#            unprofiled times are what make it possible to state how much of the
+#            measured Python time is scalene's own instrumentation overhead.
+#
+# The solver name is part of the output path, because which Pyomo solver
+# interface is used turns out to dominate the answer: "gurobi" is the
+# file-based interface, which writes an LP file and parses a solution file in
+# Python, while "gurobi_persistent" keeps the model in the solver through its C
+# API. Run the sweep once per interface and compare.
+#
+# Profiles land in $RESULTS/<solver>/<case>/rep<n>/scalene_rank_<rank>.json
+# Then run ./make_tables.bash to regenerate the LaTeX.
+
+set -euo pipefail
+
+HERE="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+REPO="$(cd "${HERE}/../.." && pwd)"
+EXAMPLES="${REPO}/examples"
+DRIVER="${REPO}/mpisppy/generic_cylinders.py"
+
+SOLVER="${SOLVER:-gurobi}"
+REPS="${REPS:-3}"
+NP="${NP:-3}"
+THREADS="${THREADS:-2}"
+RESULTS="${RESULTS:-${HERE}/results}"
+TRIES="${TRIES:-3}"
+PROFILE="${PROFILE:-1}"
+
+# Cylinders and PH settings shared by every case. rel-gap is 0 so that runs end
+# on max-iterations (or exact convergence) instead of on a gap tolerance, which
+# keeps the amount of work per case predictable.
+COMMON=(--solver-name "${SOLVER}"
+        --max-solver-threads "${THREADS}"
+        --default-rho 1
+        --lagrangian
+        --xhatshuffle
+        --rel-gap 0.0)
+
+# Case definitions. Iteration counts were calibrated so that the "long" cases
+# each take roughly a minute of wall time; farmer3 is deliberately left short
+# because the short-run numbers are themselves a finding.
+#
+# The bundled cases exist because subproblem size turns out to drive the answer.
+# farmer240 and farmer240_bun10 are a controlled pair: same instance, same
+# iteration count, differing only in whether ten scenarios are gathered into one
+# subproblem. sslp_15_45_15_bun3 is the same comparison on a MIP.
+case_args() {
+  case "$1" in
+    farmer3)
+      echo "--module-name ${EXAMPLES}/farmer/farmer --num-scens 3 --max-iterations 100"
+      ;;
+    farmer60)
+      echo "--module-name ${EXAMPLES}/farmer/farmer --num-scens 60 --max-iterations 400"
+      ;;
+    farmer240)
+      echo "--module-name ${EXAMPLES}/farmer/farmer --num-scens 240 --max-iterations 200"
+      ;;
+    farmer240_bun10)
+      echo "--module-name ${EXAMPLES}/farmer/farmer --num-scens 240 --scenarios-per-bundle 10 --max-iterations 200"
+      ;;
+    sslp_15_45_10)
+      echo "--module-name ${EXAMPLES}/sslp/sslp --sslp-data-path ${EXAMPLES}/sslp/data --instance-name sslp_15_45_10 --max-iterations 50"
+      ;;
+    sslp_15_45_15)
+      echo "--module-name ${EXAMPLES}/sslp/sslp --sslp-data-path ${EXAMPLES}/sslp/data --instance-name sslp_15_45_15 --max-iterations 50"
+      ;;
+    sslp_15_45_15_bun3)
+      echo "--module-name ${EXAMPLES}/sslp/sslp --sslp-data-path ${EXAMPLES}/sslp/data --instance-name sslp_15_45_15 --scenarios-per-bundle 3 --max-iterations 50"
+      ;;
+    sslp_5_25_50)
+      echo "--module-name ${EXAMPLES}/sslp/sslp --sslp-data-path ${EXAMPLES}/sslp/data --instance-name sslp_5_25_50 --max-iterations 150"
+      ;;
+    *)
+      echo "unknown case: $1" >&2
+      return 1
+      ;;
+  esac
+}
+
+ALL_CASES=(farmer3 farmer60 farmer240 farmer240_bun10
+           sslp_15_45_10 sslp_15_45_15 sslp_15_45_15_bun3 sslp_5_25_50)
+CASES=("${@:-}")
+if [[ -z "${CASES[0]}" ]]; then
+  CASES=("${ALL_CASES[@]}")
+fi
+
+echo "solver=${SOLVER} ranks=${NP} threads=${THREADS} reps=${REPS}"
+echo "results -> ${RESULTS}"
+
+# Scalene instruments the import machinery, and occasionally a rank dies during
+# startup with a KeyError out of importlib's lock handling before any mpi-sppy
+# code runs. That is a profiler startup race, not a property of the case, so a
+# rep that comes back with fewer than NP profiles is simply run again. Failures
+# are counted and reported at the end so a retried-away problem is still
+# visible; a rep that never succeeds is left out of the results rather than
+# silently reported as if it had worked.
+failures=0
+retries=0
+
+# Unprofiled runs go somewhere else so that they cannot overwrite the profiles.
+SUBDIR="${SOLVER}"
+if [[ "${PROFILE}" == "0" ]]; then
+  SUBDIR="${SOLVER}-unprofiled"
+fi
+
+for c in "${CASES[@]}"; do
+  read -r -a cargs <<< "$(case_args "$c")"
+  for rep in $(seq 1 "${REPS}"); do
+    outdir="${RESULTS}/${SUBDIR}/${c}/rep${rep}"
+    echo "^^^ ${c} rep ${rep} ^^^"
+    for try in $(seq 1 "${TRIES}"); do
+      rm -rf "${outdir}"
+      mkdir -p "${outdir}"
+      if [[ "${PROFILE}" == "0" ]]; then
+        # Same case, no profiler: record wall time only.
+        /usr/bin/time -f "%e" -o "${outdir}/wall.txt" \
+          mpiexec --oversubscribe -np "${NP}" \
+          python3 -m mpi4py "${DRIVER}" "${cargs[@]}" "${COMMON[@]}" \
+          > "${outdir}/run.log" 2>&1 || true
+        if grep -q 'Cylinder finalization complete' "${outdir}/run.log"; then
+          n_found="${NP}"
+          echo "  wall $(cat "${outdir}/wall.txt")s (unprofiled)"
+          break
+        fi
+        n_found=0
+      else
+        OUTDIR="${outdir}" mpiexec --oversubscribe -np "${NP}" \
+          "${HERE}/scalene_wrapper.bash" \
+          "${DRIVER}" "${cargs[@]}" "${COMMON[@]}" \
+          > "${outdir}/run.log" 2>&1 || true
+        n_found=$(find "${outdir}" -name 'scalene_rank_*.json' | wc -l)
+        if [[ "${n_found}" -eq "${NP}" ]]; then
+          break
+        fi
+      fi
+      if [[ "${PROFILE}" == "0" ]]; then
+        echo "  attempt ${try}/${TRIES}: run did not complete" \
+             "(see ${outdir}/run.log)" >&2
+      else
+        echo "  attempt ${try}/${TRIES}: expected ${NP} profiles, found ${n_found}" \
+             "(see ${outdir}/run.log)" >&2
+      fi
+      retries=$((retries + 1))
+    done
+    if [[ "${n_found}" -ne "${NP}" ]]; then
+      echo "  GIVING UP on ${c} rep ${rep} after ${TRIES} attempts" >&2
+      failures=$((failures + 1))
+      continue
+    fi
+    if [[ "${PROFILE}" != "0" ]]; then
+      tail -1 "${outdir}/run.log"
+    fi
+  done
+done
+
+echo "done: ${retries} retried attempt(s), ${failures} rep(s) abandoned"
+if [[ "${failures}" -ne 0 ]]; then
+  exit 1
+fi
diff --git a/examples/python_fraction/scalene_summary_file_interface.tex b/examples/python_fraction/scalene_summary_file_interface.tex
new file mode 100644
index 000000000..4f589362a
--- /dev/null
+++ b/examples/python_fraction/scalene_summary_file_interface.tex
@@ -0,0 +1,35 @@
+% Auto-generated by summarize_reps.py -- do not edit by hand
+
+\begin{table}[ht]
+\centering
+\begin{tabular}{l r r r r r}
+\hline
+Case & Reps & Wall (s) & Python (\%) & Native (\%) & System (\%) \\
+\hline
+farmer3 & 3 & 5.8 & 19.8 (18.5--21.5) & 68.1 & 12.1 \\
+farmer60 & 3 & 72.1 & 73.9 (73.1--74.3) & 16.7 & 9.4 \\
+sslp\_15\_45\_10 & 3 & 72.2 & 31.8 (31.7--32.0) & 66.7 & 1.5 \\
+sslp\_15\_45\_15 & 3 & 115.3 & 33.4 (33.3--33.5) & 65.2 & 1.4 \\
+sslp\_5\_25\_50 & 3 & 124.3 & 75.4 (75.3--75.5) & 20.7 & 3.9 \\
+\hline
+\end{tabular}
+\caption{Time split by case, averaged over repetitions, with the min--max range over repetitions shown for the Python percentage. Percentages are of the time Scalene attributed to a source line (95.1--99.0\% of wall time here); wall time is the maximum over ranks.}
+\label{tab:python-fraction-summary}
+\end{table}
+
+\begin{table}[ht]
+\centering
+\begin{tabular}{l r r r}
+\hline
+Case & PH hub & lagrangian & xhatshuffle \\
+\hline
+farmer3 & 18.3 (17.8--18.8) & 21.6 (18.6--25.2) & 19.5 (18.6--20.3) \\
+farmer60 & 71.4 (70.2--72.3) & 76.5 (75.9--77.0) & 73.7 (73.4--74.1) \\
+sslp\_15\_45\_10 & 33.0 (32.9--33.2) & 32.6 (32.3--32.8) & 29.8 (29.5--30.0) \\
+sslp\_15\_45\_15 & 27.0 (26.9--27.0) & 36.3 (36.1--36.4) & 37.1 (36.9--37.3) \\
+sslp\_5\_25\_50 & 72.2 (72.1--72.4) & 73.3 (73.1--73.5) & 80.7 (80.5--81.1) \\
+\hline
+\end{tabular}
+\caption{Python percentage by cylinder: mean over repetitions with the min--max range in parentheses.}
+\label{tab:python-fraction-by-rank}
+\end{table}
diff --git a/examples/python_fraction/scalene_summary_persistent.tex b/examples/python_fraction/scalene_summary_persistent.tex
new file mode 100644
index 000000000..13e03b9c5
--- /dev/null
+++ b/examples/python_fraction/scalene_summary_persistent.tex
@@ -0,0 +1,41 @@
+% Auto-generated by summarize_reps.py -- do not edit by hand
+
+\begin{table}[ht]
+\centering
+\begin{tabular}{l r r r r r r}
+\hline
+Case & Reps & Wall (s) & Python (\%) & Native (\%) & System (\%) & Overhead \\
+\hline
+farmer3 & 3 & 5.3 & 16.2 (14.4--17.6) & 70.8 & 13.0 & 2.36$\times$ \\
+farmer60 & 3 & 43.4 & 71.3 (70.3--72.4) & 16.2 & 12.4 & 1.55$\times$ \\
+farmer240 & 3 & 82.9 & 76.1 (75.2--77.0) & 11.8 & 12.2 & 1.42$\times$ \\
+farmer240\_bun10 & 3 & 20.9 & 62.6 (61.9--63.2) & 28.6 & 8.8 & 1.48$\times$ \\
+sslp\_15\_45\_10 & 3 & 74.0 & 8.7 (8.5--9.0) & 89.3 & 2.0 & 1.09$\times$ \\
+sslp\_15\_45\_15 & 3 & 104.8 & 9.6 (9.4--9.9) & 88.5 & 1.9 & 1.11$\times$ \\
+sslp\_15\_45\_15\_bun3 & 3 & 73.7 & 12.5 (12.3--12.7) & 86.1 & 1.4 & 0.99$\times$ \\
+sslp\_5\_25\_50 & 3 & 74.0 & 58.1 (58.0--58.3) & 37.1 & 4.7 & 1.22$\times$ \\
+\hline
+\end{tabular}
+\caption{Time split by case, averaged over repetitions, with the min--max range over repetitions shown for the Python percentage. Percentages are of the time Scalene attributed to a source line (91.1--99.2\% of wall time here); wall time is the maximum over ranks. The overhead column is profiled wall time divided by unprofiled wall time for the same case; because scalene's instrumentation cost lands mostly on Python, it inflates the Python column.}
+\label{tab:python-fraction-summary}
+\end{table}
+
+\begin{table}[ht]
+\centering
+\begin{tabular}{l r r r}
+\hline
+Case & PH hub & lagrangian & xhatshuffle \\
+\hline
+farmer3 & 14.1 (12.7--15.5) & 18.6 (17.5--20.5) & 15.9 (12.6--18.4) \\
+farmer60 & 72.0 (71.3--72.9) & 72.6 (70.4--73.9) & 69.3 (68.6--70.4) \\
+farmer240 & 75.9 (75.5--76.4) & 77.7 (76.6--78.6) & 74.6 (73.6--76.1) \\
+farmer240\_bun10 & 62.3 (61.1--63.4) & 64.8 (64.1--65.7) & 60.8 (60.2--61.5) \\
+sslp\_15\_45\_10 & 6.8 (6.5--7.3) & 6.3 (6.1--6.5) & 13.0 (12.8--13.2) \\
+sslp\_15\_45\_15 & 7.2 (6.7--7.5) & 6.7 (6.3--7.0) & 14.9 (14.6--15.3) \\
+sslp\_15\_45\_15\_bun3 & 5.1 (5.0--5.3) & 7.1 (6.9--7.2) & 25.3 (25.0--25.6) \\
+sslp\_5\_25\_50 & 48.6 (48.4--48.8) & 48.8 (48.1--49.8) & 78.1 (77.2--79.2) \\
+\hline
+\end{tabular}
+\caption{Python percentage by cylinder: mean over repetitions with the min--max range in parentheses.}
+\label{tab:python-fraction-by-rank}
+\end{table}
diff --git a/examples/python_fraction/scalene_totals.py b/examples/python_fraction/scalene_totals.py
new file mode 100644
index 000000000..68f98a58b
--- /dev/null
+++ b/examples/python_fraction/scalene_totals.py
@@ -0,0 +1,149 @@
+###############################################################################
+# mpi-sppy: MPI-based Stochastic Programming in PYthon
+#
+# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for
+# Sustainable Energy, LLC, The Regents of the University of California, et al.
+# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
+# full copyright and license information.
+###############################################################################
+"""
+scalene_totals.py
+
+Extract Python / native / system time totals from a Scalene JSON profile.
+
+Scalene's JSON gives, for every source line it attributed samples to, the
+percentage of the run's wall time spent on that line in Python, in native
+(compiled) code, and in the system.  Summing those per-line percentages over
+every line of every file yields the whole-run totals, which is exactly what
+this module does.
+
+Reading the JSON is preferred over parsing `scalene view --cli` output:
+
+  * the JSON carries full precision, whereas the CLI rounds each line to a
+    whole percent,
+  * `scalene view --reduced` omits low-usage lines, so summing its rows
+    undercounts,
+  * the CLI emits ANSI colour codes and rearranges its table between Scalene
+    releases, which is what made the original version of this code fragile.
+
+The result is self-checking: the sum of the per-line percentages is compared
+against the sum of Scalene's own per-file ``percent_cpu_time`` values.  The two
+agree exactly for most profiles and differ by at most a couple of tenths of a
+percentage point for the rest, because a few samples are attributed to a file
+without landing on one of the lines the file reports.  ``consistency_error``
+therefore allows a small absolute slack; it exists to catch a change in
+Scalene's JSON layout, which would show up as a large disagreement, not to
+police that last tenth of a point.
+
+Note that the per-file ``functions`` lists are *not* a usable substitute for the
+per-line data: summing them overshoots ``percent_cpu_time`` by as much as 30
+percentage points on these profiles, evidently because nested and wrapped
+functions get counted more than once.
+
+The three buckets together should account for close to 100% of wall time.
+Whatever is missing from 100% is samples Scalene could not attribute at all; it
+is reported as ``accounted_pct`` so a reader can judge how much is unexplained.
+"""
+
+from __future__ import annotations
+
+import json
+from dataclasses import dataclass
+from typing import Any, Dict, List, Optional
+
+
+@dataclass
+class Totals:
+    """Whole-run totals for one Scalene profile (one MPI rank)."""
+
+    wall_sec: Optional[float]
+    python_pct: float  # percent of wall time, Python (interpreted) code
+    native_pct: float  # percent of wall time, compiled code (solver, numpy, MPI)
+    system_pct: float  # percent of wall time, system/kernel
+    argv: Optional[List[str]]
+
+    @property
+    def accounted_pct(self) -> float:
+        """Percent of wall time Scalene attributed to some source line."""
+        return self.python_pct + self.native_pct + self.system_pct
+
+    @property
+    def python_fraction(self) -> Optional[float]:
+        """Python as a percent of attributed time; None if nothing attributed.
+
+        This is the headline number: it divides out the unattributed residual
+        so that the three buckets sum to 100%.
+        """
+        if self.accounted_pct <= 0.0:
+            return None
+        return 100.0 * self.python_pct / self.accounted_pct
+
+    def seconds(self, pct: float) -> Optional[float]:
+        if self.wall_sec is None:
+            return None
+        return self.wall_sec * pct / 100.0
+
+    @property
+    def python_sec(self) -> Optional[float]:
+        return self.seconds(self.python_pct)
+
+    @property
+    def native_sec(self) -> Optional[float]:
+        return self.seconds(self.native_pct)
+
+    @property
+    def system_sec(self) -> Optional[float]:
+        return self.seconds(self.system_pct)
+
+
+def totals_from_json_obj(j: Dict[str, Any]) -> Totals:
+    py = nat = sys_ = 0.0
+    for finfo in (j.get("files") or {}).values():
+        for line in finfo.get("lines") or ():
+            py += line.get("n_cpu_percent_python", 0.0) or 0.0
+            nat += line.get("n_cpu_percent_c", 0.0) or 0.0
+            sys_ += line.get("n_sys_percent", 0.0) or 0.0
+
+    wall = j.get("elapsed_time_sec")
+    try:
+        wall = float(wall) if wall is not None else None
+    except (TypeError, ValueError):
+        wall = None
+
+    argv = j.get("args")
+    if isinstance(argv, str):
+        argv = argv.split()
+    elif not isinstance(argv, list):
+        argv = None
+    else:
+        argv = [str(x) for x in argv]
+
+    return Totals(wall_sec=wall, python_pct=py, native_pct=nat,
+                  system_pct=sys_, argv=argv)
+
+
+def totals_from_json(path: str) -> Totals:
+    with open(path, "r", encoding="utf-8") as f:
+        return totals_from_json_obj(json.load(f))
+
+
+def consistency_error(path: str, tol_pct_points: float = 1.0) -> Optional[str]:
+    """Cross-check the per-line sum against Scalene's own per-file totals.
+
+    Returns None when they agree to within ``tol_pct_points`` percentage points,
+    otherwise a message describing the mismatch.  The tolerance is deliberately
+    loose: a real layout change misses by tens of points, while normal profiles
+    agree exactly or to within a few tenths (see the module docstring).
+    """
+    with open(path, "r", encoding="utf-8") as f:
+        j = json.load(f)
+    t = totals_from_json_obj(j)
+    per_file = sum(
+        (finfo.get("percent_cpu_time") or 0.0)
+        for finfo in (j.get("files") or {}).values()
+    )
+    if abs(per_file - t.accounted_pct) > tol_pct_points:
+        return (f"{path}: per-line sum {t.accounted_pct:.6f}% disagrees with sum of "
+                f"per-file percent_cpu_time {per_file:.6f}% by more than "
+                f"{tol_pct_points} percentage points")
+    return None
diff --git a/examples/python_fraction/scalene_wrapper.bash b/examples/python_fraction/scalene_wrapper.bash
new file mode 100755
index 000000000..d0b521a54
--- /dev/null
+++ b/examples/python_fraction/scalene_wrapper.bash
@@ -0,0 +1,40 @@
+#!/bin/bash
+###############################################################################
+# mpi-sppy: MPI-based Stochastic Programming in PYthon
+#
+# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for
+# Sustainable Energy, LLC, The Regents of the University of California, et al.
+# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
+# full copyright and license information.
+###############################################################################
+
+# Rank-aware scalene launcher, meant to be the program that mpiexec runs.
+# Each rank profiles itself and writes $OUTDIR/scalene_rank_<rank>.json
+#
+# Required environment:
+#   OUTDIR  directory to write the per-rank profile into (must already exist)
+# Arguments:
+#   the script and arguments to profile (use absolute paths; cwd is $OUTDIR)
+#
+# Not normally run by hand; see run_experiments.bash.
+
+set -euo pipefail
+
+if [[ -z "${OUTDIR:-}" ]]; then
+  echo "scalene_wrapper.bash: OUTDIR must be set" >&2
+  exit 1
+fi
+
+# Determine MPI rank from common env vars (OpenMPI / MPICH / Slurm)
+RANK="${OMPI_COMM_WORLD_RANK:-${PMI_RANK:-${SLURM_PROCID:-}}}"
+if [[ -z "${RANK}" ]]; then
+  echo "Could not determine MPI rank from environment" \
+       "(OMPI_COMM_WORLD_RANK / PMI_RANK / SLURM_PROCID)." >&2
+  exit 1
+fi
+
+cd "${OUTDIR}"
+
+exec python3 -m scalene run \
+  --outfile "scalene_rank_${RANK}.json" \
+  "$@"
diff --git a/examples/python_fraction/summarize_reps.py b/examples/python_fraction/summarize_reps.py
new file mode 100644
index 000000000..34e1b8710
--- /dev/null
+++ b/examples/python_fraction/summarize_reps.py
@@ -0,0 +1,319 @@
+#!/usr/bin/env python3
+###############################################################################
+# mpi-sppy: MPI-based Stochastic Programming in PYthon
+#
+# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for
+# Sustainable Energy, LLC, The Regents of the University of California, et al.
+# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
+# full copyright and license information.
+###############################################################################
+"""
+summarize_reps.py
+
+Aggregate the repeated runs produced by run_experiments.bash into LaTeX tables.
+
+Scalene estimates the Python/native/system split by sampling, so one run of one
+case is not evidence of much. This script reads every repetition of every case
+and reports the spread across repetitions, which is what makes it possible to
+say whether a difference between cases is real or just sampling noise.
+
+Expected layout (as written by run_experiments.bash):
+
+    <results>/<solver>/<case>/rep<n>/scalene_rank_<rank>.json
+
+and, for runs made with PROFILE=0,
+
+    <results>/<solver>-unprofiled/<case>/rep<n>/wall.txt
+
+Two tables are produced:
+
+  Summary   one row per case: wall time and the job-level Python/native/system
+            split, averaged over repetitions, with the observed min--max range
+            of the Python percentage. When unprofiled wall times are available,
+            a column reports scalene's overhead, since the profiler's own cost
+            falls mostly on Python and therefore inflates the Python share.
+
+  Per-rank  one row per case and one column per rank: the Python percentage for
+            that cylinder, averaged over repetitions with its range. This is the
+            table that shows whether a particular cylinder is an outlier.
+
+The job-level percentage for one repetition is computed from summed seconds
+across ranks, not by averaging per-rank percentages, so that ranks are weighted
+by how long they actually ran.
+
+Usage:
+  python3 summarize_reps.py --results results --out scalene_summary.tex
+"""
+
+from __future__ import annotations
+
+import argparse
+import glob
+import os
+import re
+from dataclasses import dataclass
+from typing import Dict, List, Optional, Sequence
+
+from scalene_totals import consistency_error, totals_from_json
+
+
+@dataclass
+class RepResult:
+    """One repetition of one case, aggregated over its ranks."""
+
+    rep: str
+    wall_sec: float               # max over ranks: the job's wall time
+    python_pct: float             # job-level, percent of attributed time
+    native_pct: float
+    system_pct: float
+    accounted_pct: float          # mean over ranks, for reporting
+    per_rank_python_pct: Dict[int, float]
+
+
+def _rank_of(path: str) -> Optional[int]:
+    m = re.search(r"rank[_\-]?(\d+)", os.path.basename(path))
+    return int(m.group(1)) if m else None
+
+
+def _mean(xs: Sequence[float]) -> float:
+    return sum(xs) / len(xs)
+
+
+def _latex_escape(s: str) -> str:
+    replacements = {
+        "\\": r"\textbackslash{}",
+        "&": r"\&",
+        "%": r"\%",
+        "$": r"\$",
+        "#": r"\#",
+        "_": r"\_",
+        "{": r"\{",
+        "}": r"\}",
+        "~": r"\textasciitilde{}",
+        "^": r"\textasciicircum{}",
+    }
+    return "".join(replacements.get(ch, ch) for ch in s)
+
+
+def load_rep(rep_dir: str) -> Optional[RepResult]:
+    files = sorted(glob.glob(os.path.join(rep_dir, "scalene_rank_*.json")))
+    if not files:
+        return None
+
+    walls: List[float] = []
+    accounted: List[float] = []
+    sum_py = sum_nat = sum_sys = 0.0
+    per_rank: Dict[int, float] = {}
+
+    for f in files:
+        bad = consistency_error(f)
+        if bad:
+            raise SystemExit(
+                "Scalene JSON failed its internal consistency check, so its "
+                "layout has probably changed:\n  " + bad
+            )
+        t = totals_from_json(f)
+        if t.wall_sec is None:
+            continue
+        walls.append(t.wall_sec)
+        accounted.append(t.accounted_pct)
+        sum_py += t.python_sec or 0.0
+        sum_nat += t.native_sec or 0.0
+        sum_sys += t.system_sec or 0.0
+        r = _rank_of(f)
+        if r is not None and t.python_fraction is not None:
+            per_rank[r] = t.python_fraction
+
+    if not walls:
+        return None
+
+    denom = sum_py + sum_nat + sum_sys
+    if denom <= 0.0:
+        return None
+
+    return RepResult(
+        rep=os.path.basename(rep_dir),
+        wall_sec=max(walls),
+        python_pct=100.0 * sum_py / denom,
+        native_pct=100.0 * sum_nat / denom,
+        system_pct=100.0 * sum_sys / denom,
+        accounted_pct=_mean(accounted),
+        per_rank_python_pct=per_rank,
+    )
+
+
+def load_case(case_dir: str) -> List[RepResult]:
+    reps = []
+    for rep_dir in sorted(glob.glob(os.path.join(case_dir, "rep*"))):
+        if not os.path.isdir(rep_dir):
+            continue
+        r = load_rep(rep_dir)
+        if r is not None:
+            reps.append(r)
+    return reps
+
+
+def load_unprofiled_walls(case_dir: str) -> List[float]:
+    """Wall times from PROFILE=0 runs of one case, if any were made."""
+    walls = []
+    for wf in sorted(glob.glob(os.path.join(case_dir, "rep*", "wall.txt"))):
+        try:
+            with open(wf, "r", encoding="utf-8") as f:
+                walls.append(float(f.read().strip()))
+        except (OSError, ValueError):
+            continue
+    return walls
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--results", default="results", help="Results root directory")
+    ap.add_argument("--out", default="scalene_summary.tex", help="Output LaTeX filename")
+    ap.add_argument("--solvers", nargs="*", default=None,
+                    help="Solver subdirectories to report, in order "
+                         "(default: all, alphabetical)")
+    ap.add_argument("--cases", nargs="*", default=None,
+                    help="Cases in the order to report (default: all, alphabetical)")
+    ap.add_argument("--rank-labels", default="",
+                    help="Comma-separated cylinder names for ranks 0,1,2,... "
+                         "(e.g. 'PH hub,lagrangian,xhatshuffle')")
+    args = ap.parse_args()
+
+    if args.solvers:
+        solvers = args.solvers
+    else:
+        solvers = sorted(
+            d for d in os.listdir(args.results)
+            if os.path.isdir(os.path.join(args.results, d))
+            and not d.endswith("-unprofiled")
+        )
+
+    labels = [s.strip() for s in args.rank_labels.split(",") if s.strip()]
+
+    def rank_label(r: int) -> str:
+        return labels[r] if r < len(labels) else f"rank {r}"
+
+    # loaded holds (solver, case, reps, unprofiled_walls)
+    loaded = []
+    for sv in solvers:
+        sv_dir = os.path.join(args.results, sv)
+        if args.cases:
+            cases = args.cases
+        else:
+            cases = sorted(
+                d for d in os.listdir(sv_dir)
+                if os.path.isdir(os.path.join(sv_dir, d))
+            )
+        for c in cases:
+            case_dir = os.path.join(sv_dir, c)
+            if not os.path.isdir(case_dir):
+                continue
+            reps = load_case(case_dir)
+            if not reps:
+                print(f"warning: no usable repetitions for {sv}/{c}")
+                continue
+            bare = load_unprofiled_walls(
+                os.path.join(args.results, f"{sv}-unprofiled", c)
+            )
+            loaded.append((sv, c, reps, bare))
+
+    if not loaded:
+        raise SystemExit(f"No usable results under {args.results}")
+
+    multi_solver = len({sv for sv, _, _, _ in loaded}) > 1
+    any_bare = any(bare for _, _, _, bare in loaded)
+
+    out: List[str] = []
+    out.append("% Auto-generated by summarize_reps.py -- do not edit by hand")
+    out.append("")
+
+    # ---- Summary table: one row per (solver, case) ----
+    acct_all = [r.accounted_pct for _, _, reps, _ in loaded for r in reps]
+    lead_cols = "l l" if multi_solver else "l"
+    lead_head = "Solver & Case" if multi_solver else "Case"
+
+    def lead_cells(sv: str, c: str) -> str:
+        return (f"{_latex_escape(sv)} & {_latex_escape(c)}" if multi_solver
+                else _latex_escape(c))
+
+    out.append(r"\begin{table}[ht]")
+    out.append(r"\centering")
+    out.append(rf"\begin{{tabular}}{{{lead_cols} r r r r r{' r' if any_bare else ''}}}")
+    out.append(r"\hline")
+    out.append(lead_head + r" & Reps & Wall (s) & Python (\%) & Native (\%) & System (\%)"
+               + (r" & Overhead" if any_bare else "") + r" \\")
+    out.append(r"\hline")
+    for sv, c, reps, bare in loaded:
+        pys = [r.python_pct for r in reps]
+        walls = [r.wall_sec for r in reps]
+        row = (
+            f"{lead_cells(sv, c)} & {len(reps)} & {_mean(walls):.1f} & "
+            f"{_mean(pys):.1f} ({min(pys):.1f}--{max(pys):.1f}) & "
+            f"{_mean([r.native_pct for r in reps]):.1f} & "
+            f"{_mean([r.system_pct for r in reps]):.1f}"
+        )
+        if any_bare:
+            row += (rf" & {_mean(walls) / _mean(bare):.2f}$\times$" if bare
+                    else r" & \textemdash")
+        out.append(row + r" \\")
+    out.append(r"\hline")
+    out.append(r"\end{tabular}")
+    caption = (
+        r"\caption{Time split by case, averaged over repetitions, with the "
+        r"min--max range over repetitions shown for the Python percentage. "
+        r"Percentages are of the time Scalene attributed to a source line "
+        rf"({min(acct_all):.1f}--{max(acct_all):.1f}\% of wall time here); "
+        r"wall time is the maximum over ranks."
+    )
+    if any_bare:
+        caption += (
+            r" The overhead column is profiled wall time divided by unprofiled "
+            r"wall time for the same case; because scalene's instrumentation cost "
+            r"lands mostly on Python, it inflates the Python column."
+        )
+    out.append(caption + r"}")
+    out.append(r"\label{tab:python-fraction-summary}")
+    out.append(r"\end{table}")
+    out.append("")
+
+    # ---- Per-rank table: one row per (solver, case) ----
+    all_ranks = sorted({r for _, _, reps, _ in loaded
+                        for rep in reps for r in rep.per_rank_python_pct})
+    out.append(r"\begin{table}[ht]")
+    out.append(r"\centering")
+    out.append(rf"\begin{{tabular}}{{{lead_cols}" + " r" * len(all_ranks) + r"}")
+    out.append(r"\hline")
+    out.append(lead_head + " & "
+               + " & ".join(_latex_escape(rank_label(r)) for r in all_ranks) + r" \\")
+    out.append(r"\hline")
+    for sv, c, reps, _bare in loaded:
+        cells = []
+        for r in all_ranks:
+            vals = [rep.per_rank_python_pct[r] for rep in reps if r in rep.per_rank_python_pct]
+            cells.append(f"{_mean(vals):.1f} ({min(vals):.1f}--{max(vals):.1f})" if vals
+                         else r"\textemdash")
+        out.append(f"{lead_cells(sv, c)} & " + " & ".join(cells) + r" \\")
+    out.append(r"\hline")
+    out.append(r"\end{tabular}")
+    out.append(
+        r"\caption{Python percentage by cylinder: mean over repetitions with the "
+        r"min--max range in parentheses.}"
+    )
+    out.append(r"\label{tab:python-fraction-by-rank}")
+    out.append(r"\end{table}")
+    out.append("")
+
+    with open(args.out, "w", encoding="utf-8") as f:
+        f.write("\n".join(out))
+
+    print(f"Wrote LaTeX to: {args.out}")
+    for sv, c, reps, bare in loaded:
+        pys = [r.python_pct for r in reps]
+        walls = [r.wall_sec for r in reps]
+        extra = f", overhead {_mean(walls) / _mean(bare):.2f}x" if bare else ""
+        print(f"  {sv}/{c}: {len(reps)} reps, wall {_mean(walls):.1f}s, "
+              f"Python {_mean(pys):.1f}% ({min(pys):.1f}-{max(pys):.1f}){extra}")
+
+
+if __name__ == "__main__":
+    main()