MIMO tangential-Schur/Pick and J-inner benchmark

Tutorial goal

Measure finite MIMO tangential Pick, constant Schur recovery, and Potapov/J-inner diagnostic costs.

Note

New to the terminology? See the lattice DSP concept map and the causality/data-use guide for how online, offline, block, and MIMO examples should be read.

Context

The tangential-Schur layer is a mathematical bridge between MIMO interpolation, J-inner/lossless systems, and the matrix-lattice examples. This benchmark measures the finite definite baseline: constructing the right-tangential Pick matrix, checking its Hermitian spectrum, recovering a compatible constant Schur matrix, and evaluating elementary Potapov J-inner products on a boundary grid.

The benchmark also includes a diagonal-MIMO sanity case. When tangential directions are coordinate vectors and the Schur map is diagonal, the full MIMO Pick matrix decomposes into independent scalar Pick blocks. The benchmark reports both numerical block/eigenvalue agreement and the cost difference between treating the problem as dense MIMO data versus scalar subproblems.

Key idea and equations

For right tangential data

\[S(z_i) U_i = V_i,\]

the Pick blocks are

\[P_{ij} = \frac{U_i^H U_j - V_i^H V_j}{1 - \overline{z_i} z_j}.\]

A finite definite Schur-class certificate requires \(P \succeq 0\). The J-inner diagnostic checks elementary Potapov products on the unit circle:

\[\Theta(e^{j\omega})^H J \Theta(e^{j\omega}) \approx J.\]

In the diagonal case, coordinate directions should make the dense MIMO Pick matrix equal a block diagonal assembly of scalar Pick matrices.

How to read the result

Use the runtime-breakdown figure to see where cost goes, the residual figure to check interpolation and J-inner accuracy, and the diagonal comparison to verify reduction to independent scalar Schur problems.

Run command

python benchmarks/tangential_schur_mimo_benchmark.py --dims 2 3 4 --points 4 8 --diagonal-points-per-channel 4 8 --boundary-grid 256 --repeats 2 --output docs/benchmarks/generated/_artifacts/tangential_schur_mimo/tangential-schur-mimo.json --csv-output docs/benchmarks/generated/_artifacts/tangential_schur_mimo/tangential-schur-mimo.csv

Run status

Return code: 0

Visual and data readout

When the benchmark gallery is built with results, this page embeds PNG summaries generated from the same JSON/CSV artifacts. The raw data stay available below as downloads so exact numbers remain reproducible without making the public page read like console output.

Figures

tangential schur diagonal mimo vs scalar blocks

tangential_schur_diagonal_mimo_vs_scalar_blocks.png

tangential schur diagonal scalar block errors

tangential_schur_diagonal_scalar_block_errors.png

tangential schur mimo error summary

tangential_schur_mimo_error_summary.png

tangential schur mimo quality summary

tangential_schur_mimo_quality_summary.png

tangential schur mimo residuals

tangential_schur_mimo_residuals.png

tangential schur mimo runtime breakdown

tangential_schur_mimo_runtime_breakdown.png

tangential schur mimo runtime summary

tangential_schur_mimo_runtime_summary.png

tangential schur mimo speedup summary

tangential_schur_mimo_speedup_summary.png

Generated data files

Source code

  1"""Benchmark finite MIMO tangential-Schur/Pick and J-inner diagnostics.
  2
  3The benchmark focuses on the pure Python/NumPy tangential-Schur layer that
  4connects MIMO interpolation data with Pick certificates and J-inner Potapov
  5factors.  It measures three public operations:
  6
  7* building and diagonalizing the right-tangential Pick matrix;
  8* recovering a compatible constant Schur matrix and checking interpolation;
  9* building/evaluating elementary J-inner Potapov products on the unit circle.
 10
 11A separate diagonal sanity block compares a full diagonal-MIMO Pick matrix with
 12independent scalar Pick blocks.  That comparison is a benchmark counterpart to
 13the diagonal-MIMO-equals-SISO examples: the full matrix formulation should agree
 14numerically with the independent scalar decomposition while exposing the extra
 15cost of treating a diagonal problem as dense MIMO data.
 16"""
 17
 18from __future__ import annotations
 19
 20import argparse
 21import csv
 22import importlib.util
 23import json
 24import math
 25import platform
 26import statistics
 27import sys
 28import time
 29from pathlib import Path
 30from types import SimpleNamespace
 31from collections.abc import Callable
 32
 33import numpy as np
 34
 35
 36def _load_lattice_api():
 37    """Import public lattice_dsp API, with a source-tree fallback for docs smoke runs.
 38
 39    A source checkout without the compiled extension cannot import the top-level
 40    package because ``lattice_dsp._core`` is missing.  The tangential-Schur module
 41    itself is pure Python, so benchmark/docs generation can still run the finite
 42    Pick/J-inner diagnostics directly from ``lattice_dsp/tangential_schur.py``.
 43    Installed wheels use the public top-level API.
 44    """
 45
 46    def from_source_tree():
 47        repo_root = Path(__file__).resolve().parents[1]
 48        module_path = repo_root / "lattice_dsp" / "tangential_schur.py"
 49        spec = importlib.util.spec_from_file_location(
 50            "_lattice_dsp_tangential_schur_bench", module_path
 51        )
 52        if spec is None or spec.loader is None:
 53            raise ModuleNotFoundError("could not load lattice_dsp/tangential_schur.py")
 54        module = importlib.util.module_from_spec(spec)
 55        sys.modules[spec.name] = module
 56        spec.loader.exec_module(module)
 57        return SimpleNamespace(
 58            RightTangentialSchurData=module.RightTangentialSchurData,
 59            right_tangential_pick_matrix=module.right_tangential_pick_matrix,
 60            pick_matrix_eigenvalues=module.pick_matrix_eigenvalues,
 61            is_tangential_schur_solvable=module.is_tangential_schur_solvable,
 62            constant_schur_solution=module.constant_schur_solution,
 63            max_tangential_residual=module.max_tangential_residual,
 64            potapov_product_from_rank_one_data=module.potapov_product_from_rank_one_data,
 65            j_unitarity_residual=module.j_unitarity_residual,
 66        )
 67
 68    try:
 69        import lattice_dsp as ld  # type: ignore
 70    except ModuleNotFoundError as exc:
 71        if exc.name != "lattice_dsp._core":
 72            raise
 73        return from_source_tree()
 74    # During local docs generation another older lattice_dsp build may be
 75    # importable on sys.path.  Use the source-tree implementation when the
 76    # newly added tangential API is absent.
 77    if not hasattr(ld, "RightTangentialSchurData"):
 78        return from_source_tree()
 79    return ld
 80
 81
 82ld = _load_lattice_api()
 83
 84
 85def median_time(fn: Callable[[], object], repeats: int) -> tuple[float, object]:
 86    values: list[float] = []
 87    result: object = None
 88    for _ in range(repeats):
 89        t0 = time.perf_counter()
 90        result = fn()
 91        values.append(time.perf_counter() - t0)
 92    return statistics.median(values), result
 93
 94
 95def random_points(rng: np.random.Generator, n_points: int, radius: float) -> np.ndarray:
 96    angles = rng.uniform(0.0, 2.0 * np.pi, size=n_points)
 97    # Spread points in the disk but keep them away from the boundary by default.
 98    radii = radius * np.sqrt(rng.uniform(0.05, 1.0, size=n_points))
 99    return radii * np.exp(1j * angles)
100
101
102def constant_contraction(rng: np.random.Generator, dim: int, scale: float) -> np.ndarray:
103    raw = rng.normal(size=(dim, dim)) + 1j * rng.normal(size=(dim, dim))
104    sigma_max = np.linalg.svd(raw, compute_uv=False)[0]
105    return scale * raw / sigma_max
106
107
108def full_mimo_data(
109    dim: int, n_points: int, multiplicity: int, *, seed: int, radius: float, scale: float
110):
111    rng = np.random.default_rng(seed)
112    s0 = constant_contraction(rng, dim, scale)
113    points = random_points(rng, n_points, radius)
114    directions = rng.normal(size=(n_points, dim, multiplicity)) + 1j * rng.normal(
115        size=(n_points, dim, multiplicity)
116    )
117    values = np.einsum("oi,nir->nor", s0, directions)
118    data = ld.RightTangentialSchurData(points, directions, values)
119    return data, s0
120
121
122def run_full_mimo_case(
123    *,
124    dim: int,
125    n_points: int,
126    multiplicity: int,
127    repeats: int,
128    boundary_grid: int,
129    seed: int,
130    radius: float,
131    scale: float,
132) -> dict[str, object]:
133    data, s0 = full_mimo_data(
134        dim,
135        n_points,
136        multiplicity,
137        seed=seed,
138        radius=radius,
139        scale=scale,
140    )
141
142    pick_build_s, pick = median_time(lambda: ld.right_tangential_pick_matrix(data), repeats)
143    pick_eig_s, eig = median_time(lambda: ld.pick_matrix_eigenvalues(pick), repeats)
144    psd_check_s, solvable = median_time(lambda: ld.is_tangential_schur_solvable(data), repeats)
145    constant_solution_s, recovered = median_time(lambda: ld.constant_schur_solution(data), repeats)
146    residual = float(ld.max_tangential_residual(data, recovered))
147    constant_error = float(np.linalg.norm(recovered - s0) / max(1.0, np.linalg.norm(s0)))
148    sigma = float(np.linalg.svd(recovered, compute_uv=False)[0])
149
150    potapov_build_s, product = median_time(
151        lambda: ld.potapov_product_from_rank_one_data(data), repeats
152    )
153    omega = np.linspace(0.0, 2.0 * np.pi, boundary_grid, endpoint=False)
154    z_grid = np.exp(1j * omega)
155    potapov_eval_s, theta = median_time(lambda: product.evaluate(z_grid), repeats)
156    j_check_s, j_residuals = median_time(lambda: ld.j_unitarity_residual(theta, product.j), repeats)
157    max_j = float(np.max(j_residuals))
158
159    total_conditions = int(data.total_conditions)
160    total_time_s = float(
161        pick_build_s
162        + pick_eig_s
163        + psd_check_s
164        + constant_solution_s
165        + potapov_build_s
166        + potapov_eval_s
167        + j_check_s
168    )
169    return {
170        "case": "full_mimo_tangential_schur",
171        "dim": dim,
172        "points": n_points,
173        "multiplicity": multiplicity,
174        "total_conditions": total_conditions,
175        "pick_size": int(pick.shape[0]),
176        "j_dimension": int(product.dimension),
177        "boundary_grid": boundary_grid,
178        "pick_build_s": float(pick_build_s),
179        "pick_eig_s": float(pick_eig_s),
180        "psd_check_s": float(psd_check_s),
181        "constant_solution_s": float(constant_solution_s),
182        "potapov_build_s": float(potapov_build_s),
183        "potapov_eval_s": float(potapov_eval_s),
184        "j_check_s": float(j_check_s),
185        "time_s": total_time_s,
186        "min_pick_eigenvalue": float(np.min(eig)),
187        "max_pick_eigenvalue": float(np.max(eig)),
188        "solvable": bool(solvable),
189        "constant_solution_sigma_max": sigma,
190        "max_tangential_residual": residual,
191        "constant_solution_relative_error": constant_error,
192        "j_inner_residual": max_j,
193    }
194
195
196def diagonal_mimo_coordinate_data(
197    dim: int, points_per_channel: int, *, seed: int, radius: float, scale: float
198):
199    rng = np.random.default_rng(seed)
200    gains = scale * np.exp(1j * rng.uniform(0.0, 2.0 * np.pi, size=dim))
201    points_blocks = [random_points(rng, points_per_channel, radius) for _ in range(dim)]
202    points = np.concatenate(points_blocks)
203    directions = np.zeros((dim * points_per_channel, dim), dtype=np.complex128)
204    values = np.zeros_like(directions)
205    for ch in range(dim):
206        start = ch * points_per_channel
207        stop = start + points_per_channel
208        directions[start:stop, ch] = 1.0
209        values[start:stop, ch] = gains[ch]
210    data = ld.RightTangentialSchurData(points, directions, values)
211    return data, points_blocks, gains
212
213
214def scalar_pick_matrix(points: np.ndarray, gain: complex) -> np.ndarray:
215    z = np.asarray(points, dtype=np.complex128)
216    pick = np.empty((z.size, z.size), dtype=np.complex128)
217    for i, zi in enumerate(z):
218        for j, zj in enumerate(z):
219            pick[i, j] = (1.0 - np.conj(gain) * gain) / (1.0 - np.conj(zi) * zj)
220    return 0.5 * (pick + pick.conj().T)
221
222
223def block_diag(mats: list[np.ndarray]) -> np.ndarray:
224    size = sum(mat.shape[0] for mat in mats)
225    out = np.zeros((size, size), dtype=np.complex128)
226    offset = 0
227    for mat in mats:
228        n = mat.shape[0]
229        out[offset : offset + n, offset : offset + n] = mat
230        offset += n
231    return out
232
233
234def run_diagonal_scalar_block_case(
235    *,
236    dim: int,
237    points_per_channel: int,
238    repeats: int,
239    seed: int,
240    radius: float,
241    scale: float,
242) -> dict[str, object]:
243    data, point_blocks, gains = diagonal_mimo_coordinate_data(
244        dim,
245        points_per_channel,
246        seed=seed,
247        radius=radius,
248        scale=scale,
249    )
250
251    full_pick_s, full_pick = median_time(lambda: ld.right_tangential_pick_matrix(data), repeats)
252    full_eig_s, full_eig = median_time(lambda: ld.pick_matrix_eigenvalues(full_pick), repeats)
253
254    def scalar_blocks():
255        blocks = [scalar_pick_matrix(points, gains[ch]) for ch, points in enumerate(point_blocks)]
256        eig = np.concatenate([np.linalg.eigvalsh(block) for block in blocks])
257        return blocks, eig
258
259    scalar_blocks_s, (blocks, scalar_eig) = median_time(scalar_blocks, repeats)
260    expected = block_diag(blocks)
261    block_error = float(np.linalg.norm(full_pick - expected) / max(1.0, np.linalg.norm(expected)))
262    eig_error = float(
263        np.linalg.norm(np.sort(full_eig) - np.sort(scalar_eig))
264        / max(1.0, np.linalg.norm(np.sort(scalar_eig)))
265    )
266    speedup = scalar_blocks_s / full_pick_s if full_pick_s > 0.0 else float("inf")
267    return {
268        "case": "diagonal_mimo_vs_scalar_blocks",
269        "dim": dim,
270        "points": int(dim * points_per_channel),
271        "points_per_channel": points_per_channel,
272        "multiplicity": 1,
273        "total_conditions": int(data.total_conditions),
274        "pick_size": int(full_pick.shape[0]),
275        "full_mimo_pick_s": float(full_pick_s),
276        "scalar_blocks_pick_s": float(scalar_blocks_s),
277        "pick_eig_s": float(full_eig_s),
278        "time_s": float(full_pick_s + full_eig_s),
279        "scalar_block_time_s": float(scalar_blocks_s),
280        "speedup_scalar_blocks_vs_full_mimo": float(speedup),
281        "min_pick_eigenvalue": float(np.min(full_eig)),
282        "diagonal_block_relative_error": block_error,
283        "diagonal_block_eigenvalue_relative_error": eig_error,
284    }
285
286
287def write_csv(rows: list[dict[str, object]], path: Path) -> None:
288    keys: list[str] = []
289    for row in rows:
290        for key in row:
291            if key not in keys:
292                keys.append(key)
293    path.parent.mkdir(parents=True, exist_ok=True)
294    with path.open("w", newline="", encoding="utf-8") as f:
295        writer = csv.DictWriter(f, fieldnames=keys)
296        writer.writeheader()
297        for row in rows:
298            writer.writerow(row)
299
300
301def finite_or_nan(value: object) -> float:
302    try:
303        x = float(value)
304    except (TypeError, ValueError):
305        return float("nan")
306    return x if math.isfinite(x) else float("nan")
307
308
309def plot_benchmark(rows: list[dict[str, object]], artifact_dir: Path) -> list[Path]:
310    try:
311        import matplotlib.pyplot as plt
312        import numpy as np
313    except Exception:
314        return []
315
316    written: list[Path] = []
317    artifact_dir.mkdir(parents=True, exist_ok=True)
318
319    full_rows = [row for row in rows if row.get("case") == "full_mimo_tangential_schur"]
320    if full_rows:
321        labels = [f"d={r['dim']}\nN={r['points']} r={r['multiplicity']}" for r in full_rows]
322        x = np.arange(len(full_rows))
323        timing_keys = [
324            "pick_build_s",
325            "pick_eig_s",
326            "constant_solution_s",
327            "potapov_build_s",
328            "potapov_eval_s",
329        ]
330        fig, ax = plt.subplots(figsize=(max(8.0, 0.65 * len(full_rows) + 3.0), 4.8))
331        bottom = np.zeros(len(full_rows))
332        for key in timing_keys:
333            values = np.array([finite_or_nan(row.get(key)) for row in full_rows], dtype=float)
334            values = np.nan_to_num(values, nan=0.0)
335            ax.bar(x, values, bottom=bottom, label=key.replace("_", " "))
336            bottom += values
337        ax.set_title("MIMO tangential-Schur runtime breakdown")
338        ax.set_ylabel("seconds")
339        ax.set_xticks(x)
340        ax.set_xticklabels(labels, rotation=30, ha="right")
341        ax.set_yscale("log")
342        ax.grid(axis="y", alpha=0.25)
343        ax.legend(fontsize="small")
344        fig.tight_layout()
345        out = artifact_dir / "tangential_schur_mimo_runtime_breakdown.png"
346        fig.savefig(out, dpi=150)
347        plt.close(fig)
348        written.append(out)
349
350        fig, ax = plt.subplots(figsize=(max(8.0, 0.65 * len(full_rows) + 3.0), 4.8))
351        metrics = [
352            "max_tangential_residual",
353            "j_inner_residual",
354            "constant_solution_relative_error",
355        ]
356        width = 0.24
357        for idx, key in enumerate(metrics):
358            values = np.array([finite_or_nan(row.get(key)) for row in full_rows], dtype=float)
359            ax.bar(x + (idx - 1) * width, values, width=width, label=key.replace("_", " "))
360        ax.set_title("MIMO tangential-Schur numerical residuals")
361        ax.set_ylabel("residual")
362        ax.set_xticks(x)
363        ax.set_xticklabels(labels, rotation=30, ha="right")
364        ax.set_yscale("log")
365        ax.grid(axis="y", alpha=0.25)
366        ax.legend(fontsize="small")
367        fig.tight_layout()
368        out = artifact_dir / "tangential_schur_mimo_residuals.png"
369        fig.savefig(out, dpi=150)
370        plt.close(fig)
371        written.append(out)
372
373    diag_rows = [row for row in rows if row.get("case") == "diagonal_mimo_vs_scalar_blocks"]
374    if diag_rows:
375        labels = [f"d={r['dim']}\npts/ch={r['points_per_channel']}" for r in diag_rows]
376        x = np.arange(len(diag_rows))
377        full = np.array(
378            [finite_or_nan(row.get("full_mimo_pick_s")) for row in diag_rows], dtype=float
379        )
380        scalar = np.array(
381            [finite_or_nan(row.get("scalar_blocks_pick_s")) for row in diag_rows], dtype=float
382        )
383        width = 0.35
384        fig, ax = plt.subplots(figsize=(max(7.0, 0.65 * len(diag_rows) + 3.0), 4.6))
385        ax.bar(x - width / 2, full, width=width, label="full dense MIMO Pick")
386        ax.bar(x + width / 2, scalar, width=width, label="independent scalar blocks")
387        ax.set_title("Diagonal MIMO Pick matrix versus scalar-block decomposition")
388        ax.set_ylabel("seconds")
389        ax.set_xticks(x)
390        ax.set_xticklabels(labels, rotation=30, ha="right")
391        ax.set_yscale("log")
392        ax.grid(axis="y", alpha=0.25)
393        ax.legend(fontsize="small")
394        fig.tight_layout()
395        out = artifact_dir / "tangential_schur_diagonal_mimo_vs_scalar_blocks.png"
396        fig.savefig(out, dpi=150)
397        plt.close(fig)
398        written.append(out)
399
400        fig, ax = plt.subplots(figsize=(max(7.0, 0.65 * len(diag_rows) + 3.0), 4.6))
401        errors = np.array(
402            [finite_or_nan(row.get("diagonal_block_relative_error")) for row in diag_rows],
403            dtype=float,
404        )
405        eig_errors = np.array(
406            [
407                finite_or_nan(row.get("diagonal_block_eigenvalue_relative_error"))
408                for row in diag_rows
409            ],
410            dtype=float,
411        )
412        # Exact zeros are possible in the diagonal sanity check.  Plot a tiny
413        # floor so the log-scale figure remains visible without changing the
414        # downloadable numeric data.
415        errors_plot = np.maximum(np.nan_to_num(errors, nan=0.0), 1e-18)
416        eig_errors_plot = np.maximum(np.nan_to_num(eig_errors, nan=0.0), 1e-18)
417        ax.bar(x - width / 2, errors_plot, width=width, label="Pick block error")
418        ax.bar(x + width / 2, eig_errors_plot, width=width, label="eigenvalue error")
419        ax.set_title("Diagonal-MIMO reduction-to-scalar numerical check")
420        ax.set_ylabel("relative error")
421        ax.set_xticks(x)
422        ax.set_xticklabels(labels, rotation=30, ha="right")
423        ax.set_yscale("log")
424        ax.grid(axis="y", alpha=0.25)
425        ax.legend(fontsize="small")
426        fig.tight_layout()
427        out = artifact_dir / "tangential_schur_diagonal_scalar_block_errors.png"
428        fig.savefig(out, dpi=150)
429        plt.close(fig)
430        written.append(out)
431
432    return written
433
434
435def parse_args() -> argparse.Namespace:
436    parser = argparse.ArgumentParser(description=__doc__)
437    parser.add_argument("--dims", type=int, nargs="+", default=[2, 3, 4])
438    parser.add_argument("--points", type=int, nargs="+", default=[4, 8])
439    parser.add_argument("--multiplicity", type=int, default=1)
440    parser.add_argument("--diagonal-points-per-channel", type=int, nargs="+", default=[4, 8])
441    parser.add_argument("--boundary-grid", type=int, default=256)
442    parser.add_argument("--repeats", type=int, default=3)
443    parser.add_argument("--radius", type=float, default=0.65)
444    parser.add_argument("--scale", type=float, default=0.55)
445    parser.add_argument("--seed", type=int, default=1901)
446    parser.add_argument("--output", type=Path, default=Path("reports/tangential-schur-mimo.json"))
447    parser.add_argument(
448        "--csv-output", type=Path, default=Path("reports/tangential-schur-mimo.csv")
449    )
450    parser.add_argument("--no-plots", action="store_true")
451    return parser.parse_args()
452
453
454def main() -> None:
455    args = parse_args()
456    if any(dim <= 0 for dim in args.dims):
457        raise SystemExit("--dims must contain positive integers")
458    if any(n <= 0 for n in args.points):
459        raise SystemExit("--points must contain positive integers")
460    if args.multiplicity <= 0:
461        raise SystemExit("--multiplicity must be positive")
462    if args.repeats <= 0:
463        raise SystemExit("--repeats must be positive")
464    if args.boundary_grid <= 0:
465        raise SystemExit("--boundary-grid must be positive")
466    if not (0.0 < args.radius < 1.0):
467        raise SystemExit("--radius must lie in (0, 1)")
468    if not (0.0 < args.scale < 1.0):
469        raise SystemExit("--scale must lie in (0, 1)")
470
471    rows: list[dict[str, object]] = []
472    for dim in args.dims:
473        for n_points in args.points:
474            rows.append(
475                run_full_mimo_case(
476                    dim=dim,
477                    n_points=n_points,
478                    multiplicity=args.multiplicity,
479                    repeats=args.repeats,
480                    boundary_grid=args.boundary_grid,
481                    seed=args.seed + 1000 * dim + 17 * n_points,
482                    radius=args.radius,
483                    scale=args.scale,
484                )
485            )
486        for points_per_channel in args.diagonal_points_per_channel:
487            rows.append(
488                run_diagonal_scalar_block_case(
489                    dim=dim,
490                    points_per_channel=points_per_channel,
491                    repeats=args.repeats,
492                    seed=args.seed + 3000 * dim + 19 * points_per_channel,
493                    radius=args.radius,
494                    scale=args.scale,
495                )
496            )
497
498    payload = {
499        "python": platform.python_version(),
500        "platform": platform.platform(),
501        "description": "Finite MIMO tangential-Schur/Pick and J-inner diagnostic benchmark.",
502        "dims": args.dims,
503        "points": args.points,
504        "multiplicity": args.multiplicity,
505        "diagonal_points_per_channel": args.diagonal_points_per_channel,
506        "boundary_grid": args.boundary_grid,
507        "repeats": args.repeats,
508        "radius": args.radius,
509        "scale": args.scale,
510        "seed": args.seed,
511        "results": rows,
512    }
513    args.output.parent.mkdir(parents=True, exist_ok=True)
514    args.output.write_text(json.dumps(payload, indent=2), encoding="utf-8")
515    write_csv(rows, args.csv_output)
516
517    written_plots: list[Path] = []
518    if not args.no_plots:
519        written_plots = plot_benchmark(rows, args.output.parent)
520
521    print(json.dumps({k: v for k, v in payload.items() if k != "results"}, indent=2))
522    print()
523    print(
524        f"{'case':>34} {'dim':>4} {'points':>6} {'cond':>6} {'total_s':>10} "
525        f"{'pick_s':>10} {'pot_eval_s':>11} {'j_resid':>10} {'interp_resid':>13}"
526    )
527    print("-" * 124)
528    for row in rows:
529        print(
530            f"{str(row['case']):>34} {int(row['dim']):4d} {int(row['points']):6d} "
531            f"{int(row['total_conditions']):6d} {finite_or_nan(row.get('time_s')):10.4e} "
532            f"{finite_or_nan(row.get('pick_build_s', row.get('full_mimo_pick_s'))):10.4e} "
533            f"{finite_or_nan(row.get('potapov_eval_s')):11.4e} "
534            f"{finite_or_nan(row.get('j_inner_residual')):10.2e} "
535            f"{finite_or_nan(row.get('max_tangential_residual', row.get('diagonal_block_relative_error'))):13.2e}"
536        )
537    print(f"\nWrote {args.output}")
538    print(f"Wrote {args.csv_output}")
539    for path in written_plots:
540        print(f"Wrote {path}")
541
542
543if __name__ == "__main__":
544    main()