"""Generate SAR3D/GPR references from notebook cells 1/5 and Horns coherent BP.

Requires the original, untracked Multi_depthnorm_clean.ipynb; a missing notebook
is an error. NumPy/SciPy execute its function/class ASTs unchanged. Only parameter
values are replaced; plotting imports and the interactive top-level runs are
excluded. No production SAR implementation or hand-transcribed oracle is used.

Run from the repository root:
    .venv/bin/python data_acq_and_processing/tests/generate_backprojection_reference.py
"""
from __future__ import annotations

import argparse
import ast
import hashlib
import json
from pathlib import Path
import tempfile
from types import SimpleNamespace

import numpy as np
import scipy

ROOT = Path(__file__).resolve().parents[2]
REFERENCE = Path(__file__).with_name("sar3d_reference.npz")


def reference_functions(source_path: Path, parameters: dict, cells: tuple[int, ...] = (1, 5)) -> tuple[dict, dict]:
    contents = source_path.read_text(encoding="utf-8")
    sources = [(str(i), "".join(json.loads(contents)["cells"][i]["source"])) for i in cells] if cells else [("script", contents)]
    namespace = {"__name__": __name__}
    hashes = {}
    for index, source in sources:
        hashes[str(index)] = hashlib.sha256(source.encode()).hexdigest()
        tree = ast.parse(source, filename=f"{source_path.name}:cell{index}")
        definitions = []
        for node in tree.body:
            if isinstance(node, (ast.FunctionDef, ast.ClassDef, ast.ImportFrom)):
                definitions.append(node)
            elif isinstance(node, ast.Import):
                if not any(alias.name.startswith(("matplotlib", "plotly")) for alias in node.names):
                    definitions.append(node)
            elif isinstance(node, ast.Assign) and len(node.targets) == 1 and isinstance(node.targets[0], ast.Name):
                try:
                    namespace[node.targets[0].id] = ast.literal_eval(node.value)
                except (ValueError, TypeError):
                    pass  # Dataset loading and the interactive run are supplied below.
        namespace.update(parameters)
        exec(compile(ast.Module(definitions, type_ignores=[]), f"{source_path.name}:cell{index}",
                     "exec", dont_inherit=True), namespace)
    return namespace, hashes


def scene(f_start=1e9, f_stop=4e9, points=201) -> dict:
    """Deterministic two-echo input with motion, vertical baselines and four pairs."""
    request = dict(
        tx=[[0, -.4, 0, 0], [1, .45, .2, .05]],
        rx=[[0, -.2, .1, 0], [1, .3, .35, .02]], f_start=f_start, f_stop=f_stop,
        min_depth=1., max_depth=3., speed=.4, pitch=13., ground=.3,
        voxel=[.4, .4, .4], radius=2, power=1.5, compensation=True, normalize=True,
        equalize=True, subtract_mean=False, mean_count=1, motion="int_minus",
        bg_gate=True, clean=False,
    )
    # Match the float32 settings/geometry accepted by the native fixture.
    for key in ("tx", "rx"):
        request[key] = np.asarray(request[key], dtype=np.float32).astype(float).tolist()
    frequency = np.linspace(f_start, f_stop, points, dtype=np.float32)
    members = []
    for frame, offset in enumerate((0, 180_000_000, 510_000_000)):
        traces = []
        start = 10_000_000_000 + offset
        for tx in range(2):
            for rx in range(2):
                samples = ((1 + .1 * rx) * np.exp(-2j * np.pi * frequency.astype(float)
                    * (4.8 + .08 * tx - .03 * frame) / 3e8)
                    + .35 * np.exp(-2j * np.pi * frequency.astype(float) * (5.6 - .05 * rx) / 3e8)).astype(np.complex64)
                traces.append(dict(input=rx, output=tx, start=start, end=start + 165_000_000,
                    frequency=frequency.tolist(), s21=np.column_stack((samples.real, samples.imag)).tolist()))
        members.append(dict(id=frame + 1, timestamp=start + 165_000_000, traces=traces))
    request.update(previous=members[:-1], latest=members[-1])
    return request


# Point reflectors (X, height Y, Z, amplitude) at heights of 0, 0.6 and 0.25 m.
REFLECTORS = [(.3, 0., 2.2, 1.), (-.8, .6, 2.6, .5), (1.3, .25, 1.9, .7)]


def point_scene(targets: list[tuple[float, float, float, float]], noise: float, seed: int, **settings) -> dict:
    """Settings of one static frame of point reflectors (X, height Y, Z, amplitude) in noise, for CLEAN.

    Three Tx and three Rx at different heights resolve Y in a single frame. The voxel steps
    keep the notebook's proportions, 2.5 : 2 : 2.
    """
    request = dict(
        tx=[[0, -.6, 0., 0.], [1, 0., .35, .02], [2, .6, .1, 0.]],
        rx=[[0, -.3, .2, 0.], [1, .3, .45, .03], [2, .9, 0., .01]],
        min_depth=1., max_depth=3.5, speed=0., pitch=13., ground=1., voxel=[.0625, .05, .05], radius=0,
        clean=True, previous=[],
    )
    request.update(settings)
    for key in ("tx", "rx"):
        request[key] = np.asarray(request[key], dtype=np.float32).astype(float).tolist()
    frequency = np.linspace(1e9, 4e9, 201, dtype=np.float32)
    cos_p, sin_p = np.cos(np.radians(request["pitch"])), np.sin(np.radians(request["pitch"]))
    random = np.random.default_rng(seed)
    traces = []
    for tx in request["tx"]:
        for rx in request["rx"]:
            samples = noise * (random.standard_normal(frequency.size) + 1j * random.standard_normal(frequency.size))
            for x, height, z, amplitude in targets:
                # Scene to radar axes, as sar_scene_to_radar_yz.
                y0 = height - request["ground"]
                point = np.array([x, y0 * cos_p + z * sin_p, -y0 * sin_p + z * cos_p])
                path = np.linalg.norm(point - tx[1:]) + np.linalg.norm(point - rx[1:])
                samples += amplitude * np.exp(-2j * np.pi * frequency.astype(float) * path / 3e8)
            samples = samples.astype(np.complex64)
            traces.append(dict(input=int(rx[0]), output=int(tx[0]), start=10_000_000_000, end=10_165_000_000,
                               frequency=frequency.tolist(), s21=np.column_stack((samples.real, samples.imag)).tolist()))
    request["latest"] = dict(id=1, timestamp=10_165_000_000, traces=traces)
    return request


def write_dataset(directory: Path, request: dict) -> list[Path]:
    """Supply the notebook's original config_profile/meta/NPY input contract."""
    profile = {"gpr": {}, "run": {"combos": []}}
    for kind, position in (("tx", "output_pos"), ("rx", "input_pos")):
        profile["gpr"][f"{kind}_geometry"] = [
            {position: int(index), "x_m": x, "y_m": y, "z_m": z}
            for index, x, y, z in request[kind]
        ]
    profile["run"]["combos"] = [
        {"input": trace["input"], "output": trace["output"]} for trace in request["latest"]["traces"]
    ]
    (directory / "config_profile.json").write_text(json.dumps(profile), encoding="utf-8")
    snapshots = []
    members = [*request["previous"], request["latest"]]
    for index, member in enumerate(members):
        background = previous_mean({**request, "previous": members[:index]})
        snapshot = directory / f'{index:04d}_id{member["id"]}_ns{member["timestamp"]}'
        snapshot.mkdir()
        metadata = {"traces": []}
        for trace in member["traces"]:
            stem = f'i{trace["input"]}_o{trace["output"]}'
            values = np.asarray(trace["s21"])
            samples = values[:, 0] + 1j * values[:, 1] - background.get((trace["output"], trace["input"]), 0.)
            np.save(snapshot / f"{stem}_s21.npy", samples)
            np.save(snapshot / f"{stem}_freq.npy", np.asarray(trace["frequency"], dtype=np.float32))
            metadata["traces"].append(dict(input=trace["input"], output=trace["output"],
                capture_start_ns=trace["start"], capture_end_ns=trace["end"]))
        (snapshot / "meta.json").write_text(json.dumps(metadata), encoding="utf-8")
        snapshots.append(snapshot)
    return snapshots


def reference(notebook: Path, request: dict, directory: Path) -> tuple[dict, dict]:
    snapshots = write_dataset(directory, request)
    center = snapshots[-1]
    parameters = dict(
        SAR_Y_RANGE_M=(request.get("y_min", -.5), request.get("y_max", 2.)),
        BG_GATE_ENABLE=request["bg_gate"], SAR_CLEAN_ENABLE=request["clean"], SAR_MARK_ENABLE=True,
        INPUT_IDX=[int(row[0]) for row in request["rx"]],
        OUTPUT_IDX=[int(row[0]) for row in request["tx"]],
        F_START=request["f_start"], F_STOP=request["f_stop"],
        MIN_DEPTH=request["min_depth"], MAX_DEPTH=request["max_depth"],
        PAIR_NORM_DEPTH_MIN=request["min_depth"], PAIR_NORM_DEPTH_MAX=request["max_depth"],
        SPEED_M_S=request["speed"], LOOK_ANGLE_DEG=request["pitch"],
        BP_GROUND_DEPTH_M=request["ground"], COMP_RANGE_POWER=request["power"],
        SAR_VOXEL_M=tuple(request["voxel"]), SAR_FRAME_RADIUS=request["radius"],
        SAR_USE_COMPENSATION=request["compensation"], SAR_NORMALIZE_EACH_FRAME=request["normalize"],
        PAIR_DEPTH_EQUALIZE=request["equalize"],
        SAR_SCAN_ROOT=directory, DATA_PATH=center, BG_PATH=directory, v=3e8,
    )
    ns, hashes = reference_functions(notebook, parameters)
    tx, rx, order, tx_order, pair_count = ns["load_antenna_config_profile"](
        directory / "config_profile.json", ns["INPUT_IDX"], ns["OUTPUT_IDX"])
    ns["N_tx"], ns["N_rx"] = len(tx), len(rx)
    for kind, positions in (("tx", tx), ("rx", rx)):
        for axis, values in zip("xyz", ns["positions_to_xyz"](positions).T):
            ns[f"{axis}_{kind}"] = values
    antenna_x = np.concatenate((ns["x_tx"], ns["x_rx"]))
    ns["x_ant"] = antenna_x
    ns["x_min"], ns["x_max"] = antenna_x.min() - 2, antenna_x.max() + 2
    ns["x_grid_bp"] = np.linspace(ns["x_min"], ns["x_max"], ns["NX_BP"])
    sweep_time, _, _ = ns["estimate_sweep_time_from_meta"](directory, n_pairs=pair_count)
    config = ns["MotionConfig"](
        speed_m_s=ns["SPEED_M_S"], look_angle_deg=ns["LOOK_ANGLE_DEG"],
        tx_sweep_time_s=sweep_time, frame_acq_time_s=sweep_time * pair_count,
        pair_order_phys=order, tx_order_full=tuple(tx_order),
        n_freq_hint=len(request["latest"]["traces"][0]["frequency"]),
        apply_pair_motion_correction=ns["APPLY_PAIR_MOTION_CORRECTION"],
        apply_freq_phase_correction=ns["APPLY_FREQ_PHASE_CORRECTION"],
        pair_motion_correction_mode=ns["PAIR_MOTION_CORRECTION_MODE"],
        motion_phase_sign=ns["MOTION_PHASE_SIGN"],
    )
    timestamps, _ = ns["compute_pair_timestamps_for_snapshot"](center, config)
    grid = ns["sar_build_full_region_grid"](config)
    for index, axis in ((0, "x"), (2, "z")):
        size, step = request.get(f"{axis}_size", 0.), request["voxel"][index]
        if size:
            center_axis = .5 * (grid[axis][0] + grid[axis][-1])
            grid[axis] = np.linspace(center_axis - size/2, center_axis + size/2, max(round(size/step)+1, 2))
    grid["XX"], grid["ZZ"] = np.meshgrid(grid["x"], grid["z"])
    ns["x_min"], ns["x_max"] = grid["x"][0], grid["x"][-1]
    result = ns["sar_compute_volume"](center, None, config, timestamps, grid=grid, verbose=False)
    # What the processor publishes: the frame volume, or the CLEAN display when CLEAN found
    # objects, and the peak it marks; the background gain stays to the GUI.
    ns["sar_prepare_display"](result, None, timestamps)
    objects = result["clean_objects"] or []
    volume, grid = result["display_volume"] if objects else result["volume"], result["grid"]
    iy, iz, ix = result["display_max_idx"]
    expected = {axis: grid[axis] for axis in "xyz"}
    expected.update(front=volume.max(axis=1), side=volume.max(axis=2), top=volume.max(axis=0),
                    peak=np.array([result["level"], grid["x"][ix], grid["y"][iy], grid["z"][iz]]))
    if not objects:
        # Width of the response through the marked maximum; objects are marked instead.
        from scipy.signal import peak_widths
        widths = []
        for axis, peak_index, profile in (("x", ix, volume[iy, iz, :]), ("y", iy, volume[:, iz, ix]),
                                          ("z", iz, volume[iy, :, ix])):
            profile = np.asarray(profile, dtype=np.float64)
            width, *_ = peak_widths(profile, [peak_index], prominence_data=(np.array([profile[peak_index]]),
                                        np.array([0], dtype=np.intp), np.array([len(profile)-1], dtype=np.intp)))
            widths.append(max(float(width[0]), 1.) * (grid[axis][1]-grid[axis][0]))
        expected["target_width"] = np.asarray(widths)
    if request["bg_gate"]:
        expected["background"] = np.array([result["bg_ratio"], ns["BG_GATE_3D_SHOW"], ns["BG_GATE_3D_FULL"]])
    if result["clean_info"] is not None:
        expected["clean"] = np.array([result["clean_info"]["steps"], result["clean_info"]["residual_ratio"]])
    if objects:
        kinds = {"low": 0, "ground": 1, "up": 2, "edge": 3}
        expected["objects"] = np.array([[item["x"], item["y"], item["z"], item["level_bg"]] for item in objects])
        expected["marks"] = np.array([[kinds[group["kind"]], *group["box"]] for group in result["mark_groups"]])
    return expected, hashes


def marks_reference(notebook: Path) -> tuple[dict, np.ndarray]:
    """sar_mark_groups on synthetic CLEAN objects, each of which tests one rule of the boxes.

    Returns the fixture's "marks" request and the expected rows [class, X0, X1, Y0, Y1, Z0, Z1].
    """
    scene = dict(x=[-2.6, 2.9], y=[-.5, 2.], pitch=13., ground=1., min_depth=1., max_depth=3.5)
    objects = [  # [x, y, z, y_low, y_high, level], strongest first
        [0., 0., 2., -.1, .3, 60.],         # low
        [1., .05, 2.2, -.1, .6, 55.],       # centre low, top above SAR_MARK_LOW_TOP_M: ground
        [-1., .3, 2.4, .1, .5, 50.],        # ground
        [-.2, .45, 3., .25, .65, 45.],      # up while its weak member (last) stays out of the centre
        [-1.6, 0., 1.3, -.1, .3, 40.],      # a chain of three: the ends are 1.4 m apart, still one box
        [-.9, .1, 1.3, -.1, .35, 35.],
        [-.2, .05, 1.3, -.15, .3, 30.],
        [-2.5, .2, 2., 0., .4, 28.],        # near the X minimum: edge
        [2.75, .2, 2., 0., .4, 26.],        # near the X maximum: edge
        [2.2, .1, 3.2, -.1, .3, 24.],       # near MAX_DEPTH: edge
        [2.2, 0., .85, -.2, .2, 22.],       # near MIN_DEPTH: edge
        [1.2, 1.7, 2.8, 1.5, 1.9, 20.],     # near the top: edge
        [.2, -.4, 3.1, -.5, -.25, 15.],     # weak member of the "up" box, below 40 % of its top
    ]
    ns, _ = reference_functions(notebook, dict(BP_GROUND_DEPTH_M=scene["ground"], MIN_DEPTH=scene["min_depth"],
                                               MAX_DEPTH=scene["max_depth"]))
    groups = ns["sar_mark_groups"](
        [dict(zip(("x", "y", "z", "y_lo", "y_hi", "level_bg"), row)) for row in objects],
        {axis: np.asarray(scene[axis]) for axis in "xy"}, SimpleNamespace(look_angle_deg=scene["pitch"]))
    kinds = {"low": 0, "ground": 1, "up": 2, "edge": 3}
    return dict(mode="marks", objects=objects, **scene), np.array([[kinds[g["kind"]], *g["box"]] for g in groups])


def previous_mean(request: dict) -> dict:
    samples = {}
    if request["subtract_mean"] and request["mean_count"]:
        for member in request["previous"][-request["mean_count"]:]:
            for trace in member["traces"]:
                values = np.asarray(trace["s21"])
                samples.setdefault((trace["output"], trace["input"]), []).append(values[:, 0] + 1j * values[:, 1])
    return {key: np.mean(values, axis=0) for key, values in samples.items()}


def gpr_reference(notebook: Path, script: Path, request: dict) -> tuple[np.ndarray, dict]:
    """Execute Horns coherent BP unchanged, adding the notebook's equalizer before its scalar norm."""
    parameters = dict(F_START=request["f_start"], F_STOP=request["f_stop"], MIN_DEPTH=request["min_depth"],
                      MAX_DEPTH=request["max_depth"], LOOK_ANGLE_DEG=request["pitch"],
                      BP_GROUND_DEPTH_M=request["ground"], COMP_RANGE_POWER=request["power"], v=3e8)
    ns, hashes = reference_functions(script, parameters, ())
    ns["N_tx"], ns["N_rx"] = len(request["tx"]), len(request["rx"])
    ns["OUTPUT_IDX"], ns["INPUT_IDX"] = [int(row[0]) for row in request["tx"]], [int(row[0]) for row in request["rx"]]
    for kind in ("tx", "rx"):
        for axis, values in zip("xyz", np.asarray(request[kind])[:, 1:].T):
            ns[f"{axis}_{kind}"] = values
    xs = np.concatenate((ns["x_tx"], ns["x_rx"]))
    ns["x_ant"] = xs
    ns["x_min"], ns["x_max"] = xs.min()-2, xs.max()+2
    ns["x_grid_bp"] = np.linspace(ns["x_min"], ns["x_max"], ns["NX_BP"])
    ns["z_grid_bp"] = np.linspace(.1, request["max_depth"], ns["NZ_BP"])
    ns["XX_bp"], ns["ZZ_bp"] = np.meshgrid(ns["x_grid_bp"], ns["z_grid_bp"])
    ns["depth_gate"] = (ns["ZZ_bp"] >= request["min_depth"]) & (ns["ZZ_bp"] <= request["max_depth"])
    member = request["latest"]
    duration = (max(t["end"] for t in member["traces"]) - min(t["start"] for t in member["traces"])) * 1e-9
    config = ns["MotionConfig"](speed_m_s=request["speed"], look_angle_deg=request["pitch"],
        tx_per_segment=len(request["tx"]), tx_order_full=tuple(ns["OUTPUT_IDX"]),
        frame_acq_time_s=duration, n_freq_hint=len(member["traces"][0]["frequency"]),
        intra_phase_focus_only=request["motion"] == "int_focus", motion_phase_sign=-1.,
        pair_order_phys=[(t["output"], t["input"]) for t in member["traces"]])
    clocks, _ = ns["compute_pair_timestamps"](config, ns["INPUT_IDX"], ns["OUTPUT_IDX"])
    mean = previous_mean(request)
    H, T = {}, {}
    for trace in member["traces"]:
        key = (trace["output"], trace["input"])
        values = np.asarray(trace["s21"], dtype=np.float32)
        samples = values[:, 0].astype(float) + 1j * values[:, 1].astype(float)
        samples -= mean.get(key, 0.)
        frequency = np.asarray(trace["frequency"], dtype=np.float32)
        samples, _ = ns["apply_intra_sweep_phase_correction"](samples, frequency, clocks[key], config,
                                                            ns["v"], ns["F_START"], ns["F_STOP"])
        t, _, h, _, _ = ns["compute_ascan_bp"](samples, frequency, ns["F_START"], ns["F_STOP"])
        H[key], T[key] = h, t
    if request["equalize"]:
        eq, notebook_hashes = reference_functions(notebook, {**parameters, "PAIR_DEPTH_EQUALIZE": True})
        hashes.update({f"notebook_cell_{cell}": value for cell, value in notebook_hashes.items()})
        eq.update(x_ant=xs, x_min=ns["x_min"], x_max=ns["x_max"])
        H = eq["equalize_pairs_by_depth"](H, T)
    for key, h in H.items():
        depth = T[key] * ns["v"] / 2
        mask = (depth >= request["min_depth"]) & (depth <= request["max_depth"])
        scale = np.percentile(np.abs(h)[mask] if mask.any() else np.abs(h), 50.)
        H[key] = h / ((scale if np.isfinite(scale) and scale > 1e-15 else 1.) + 1e-15)
    raw, _ = ns["backproject_coherent"](H, T)
    return ns["normalize_bp_map"](raw).astype(np.float32), hashes


def generate_gpr(notebook: Path, script: Path, output: Path) -> None:
    cases = [dict(equalize=False), dict(equalize=True),
             dict(equalize=True, subtract_mean=True, mean_count=1, speed=-.2, pitch=-7., motion="int_focus", power=2.)]
    arrays = {}
    for i, settings in enumerate(cases):
        request = scene(1.01e9, 4.01e9, 151) if i == 1 else scene()
        request.update(settings, mode="gpr")
        arrays[f"case_{i}_image"], hashes = gpr_reference(notebook, script, request)
        arrays[f"case_{i}_request"] = np.array(json.dumps(request))
    arrays["provenance"] = np.array(json.dumps(dict(script=script.name, source_sha256=hashes,
                                                  notebook=notebook.name, numpy=np.__version__, scipy=scipy.__version__)))
    np.savez_compressed(output, **arrays)
    print(f"Wrote {output}: {output.stat().st_size} bytes, {len(cases)} GPR references")


def generate(notebook: Path, output: Path) -> None:
    # Read up front so a missing notebook cannot produce an incomplete artifact.
    notebook.read_text(encoding="utf-8")
    cases = [
        dict(compensation=True, normalize=True),
        dict(compensation=True, normalize=False, equalize=False, speed=-.2, pitch=-7., ground=.5),
        # Three different voxel steps: each axis keeps its own.
        dict(compensation=False, normalize=True, f_start=1.3e9, f_stop=3.7e9, voxel=[.45, .4, .35]),
        dict(compensation=False, normalize=False, speed=0., pitch=0., radius=0),
        dict(compensation=True, normalize=True, subtract_mean=True, mean_count=1),
        dict(compensation=True, normalize=True, x_size=3., z_size=2., y_min=-.2, y_max=1.4),
        # CLEAN stops below its thresholds; the background gate hides the residual.
        point_scene(REFLECTORS, noise=.02, seed=1),
        # CLEAN stops at SAR_CLEAN_MAX_STEPS; without the gate the residual stays visible.
        point_scene(REFLECTORS, noise=.01, seed=3, bg_gate=False),
        # The stop at the background of each range (SAR_CLEAN_STOP_RANGE_BG) adds two steps to
        # those of the stop at the background of the whole volume.
        point_scene([(.2, 0., 2., 1.), (-.5, .1, 4.4, .08)], noise=.01, seed=4, max_depth=5.),
        # A narrowed frame: CLEAN still searches the whole scene, so the reflector at X = -0.8 m
        # just outside the frame is subtracted instead of its sidelobe becoming an object.
        point_scene(REFLECTORS, noise=.02, seed=19, x_size=1.6, z_size=.8),
    ]
    arrays, configurations = {}, []
    for index, settings in enumerate(cases):
        request = scene()
        request.update(settings)
        # Match the float32 live settings of the native fixture.
        for key in ("min_depth", "max_depth", "speed", "pitch", "ground", "power",
                    "x_size", "z_size", "y_min", "y_max"):
            if key in request:
                request[key] = float(np.float32(request[key]))
        request["voxel"] = [float(np.float32(step)) for step in request["voxel"]]
        if index == 1:
            # A separate, unequal capture window per output instead of interleaving.
            for member in [*request["previous"], request["latest"]]:
                for trace in member["traces"]:
                    if trace["output"] == 0:
                        trace["end"] = trace["start"] + 80_000_000
                    else:
                        trace["start"] += 95_000_000
        with tempfile.TemporaryDirectory(prefix="sar3d_notebook_") as directory:
            expected, hashes = reference(notebook, request, Path(directory))
        if index == 0:
            arrays["request"] = np.array(json.dumps(request, separators=(",", ":")))
        configurations.append(dict(
            settings={key: request[key] for key in settings},
            windows=[[[trace["start"], trace["end"]] for trace in member["traces"]]
                     for member in [*request["previous"], request["latest"]]],
        ))
        arrays.update({f"case_{index}_{key}": value for key, value in expected.items()})
    arrays["cases"] = np.array(json.dumps(configurations, separators=(",", ":")))
    marks_request, arrays["marks"] = marks_reference(notebook)
    arrays["marks_request"] = np.array(json.dumps(marks_request))
    arrays["provenance"] = np.array(json.dumps(dict(
        notebook=notebook.name, cell_source_sha256=hashes, cells=[1, 5], cases=len(cases),
        numpy=np.__version__, scipy=scipy.__version__,
        method="Unchanged notebook AST functions; per-frame NumPy previous-mean subtraction; "
               "max projections of the sar_prepare_display volume, its background, CLEAN objects and marks",
    ), sort_keys=True))
    output.parent.mkdir(parents=True, exist_ok=True)
    np.savez_compressed(output, **arrays)
    print(f"Wrote {output}: {output.stat().st_size} bytes, {len(cases)} notebook references")


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--notebook", type=Path, default=ROOT / "Multi_depthnorm_clean.ipynb")
    parser.add_argument("--output", type=Path, default=REFERENCE)
    parser.add_argument("--gpr-script", type=Path, default=ROOT / "Horns_motion_4x4_ground.py")
    parser.add_argument("--gpr-output", type=Path, default=REFERENCE.with_name("gpr_reference.npz"))
    args = parser.parse_args()
    generate(args.notebook, args.output)
    generate_gpr(args.notebook, args.gpr_script, args.gpr_output)