#!/usr/bin/env python3
"""Measure redundancy and effective dimensionality across CSI subcarriers.

The input container is:
    <float64 timestamp> <uint16 payload-length> <payload>

The 24-byte CSI payload header is:
    <I timestamp_us 6s transmitter_mac b rssi B snr B bandwidth
     B primary_channel B receive_mode B sequence/padding
     H subcarrier_count H tx_chain H rx_chain H chain_info>

For every closed measurement group this script selects tx0/rx0 only. It never
substitutes another chain. Each 256-bin amplitude vector is divided by its
median so that a shared gain change does not dominate the comparison. The
remaining analysis asks whether the 256 bins carry 256 independent dimensions
in these two captures; it does not assign physical antenna or frequency
positions to the array indices.
"""

from __future__ import annotations

import argparse
import csv
import hashlib
import html
import json
import math
import struct
from pathlib import Path

import numpy as np


RECORD = struct.Struct("<dH")
HEADER = struct.Struct("<I6sbBBBBBHHHH")
EXPECTED_SUBCARRIERS = 256
EXPECTED_CHAIN = (0, 0)


def sha256_file(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as handle:
        for chunk in iter(lambda: handle.read(1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()


def rounded(value: float, digits: int = 6) -> float:
    return round(float(value), digits)


def amplitude(iq: np.ndarray) -> np.ndarray:
    return np.hypot(iq[0::2].astype(np.float64), iq[1::2].astype(np.float64))


def parse_capture(path: Path) -> dict:
    data = path.read_bytes()
    offset = 0
    outer_records = 0
    malformed_records = 0
    pending: list[dict] = []
    groups: list[list[dict]] = []

    while offset + RECORD.size <= len(data):
        outer_ts, length = RECORD.unpack_from(data, offset)
        payload_start = offset + RECORD.size
        payload_end = payload_start + length
        if payload_end > len(data) or length < HEADER.size:
            malformed_records += 1
            break

        header = HEADER.unpack_from(data, payload_start)
        (
            payload_ts,
            station,
            rssi,
            snr,
            bandwidth,
            channel,
            receive_mode,
            sequence,
            subcarrier_count,
            tx,
            rx,
            chain_info,
        ) = header

        iq_offset = payload_start + HEADER.size
        available_pairs = max(0, (payload_end - iq_offset) // 4)
        pairs = min(subcarrier_count, available_pairs)
        iq = np.frombuffer(
            data,
            dtype="<i2",
            count=pairs * 2,
            offset=iq_offset,
        )
        pending.append(
            {
                "outer_ts": outer_ts,
                "payload_ts": payload_ts,
                "station": station,
                "rssi": rssi,
                "snr": snr,
                "bandwidth": bandwidth,
                "channel": channel,
                "receive_mode": receive_mode,
                "sequence": sequence,
                "subcarrier_count": subcarrier_count,
                "tx": tx,
                "rx": rx,
                "chain_info": chain_info,
                "pairs": pairs,
                "amplitude": amplitude(iq),
            }
        )
        outer_records += 1

        if chain_info & 0x8000:
            groups.append(pending)
            pending = []
        offset = payload_end

    return {
        "outer_records": outer_records,
        "malformed_records": malformed_records,
        "incomplete_tail_records": len(pending),
        "groups": groups,
    }


def select_fixed_chain(parsed: dict) -> dict:
    rows: list[np.ndarray] = []
    times: list[float] = []
    selected_counts: dict[str, int] = {}
    missing_fixed_chain = 0
    duplicate_fixed_chain = 0
    invalid_fixed_chain = 0

    for group in parsed["groups"]:
        matches = [
            record
            for record in group
            if (int(record["tx"]), int(record["rx"])) == EXPECTED_CHAIN
        ]
        if not matches:
            missing_fixed_chain += 1
            continue
        if len(matches) > 1:
            duplicate_fixed_chain += 1
            continue

        chosen = matches[0]
        values = chosen["amplitude"]
        if (
            int(chosen["subcarrier_count"]) != EXPECTED_SUBCARRIERS
            or int(chosen["pairs"]) != EXPECTED_SUBCARRIERS
            or values.size != EXPECTED_SUBCARRIERS
            or not np.any(values)
        ):
            invalid_fixed_chain += 1
            continue

        rows.append(values)
        times.append(float(chosen["outer_ts"]))
        chain = f"tx{chosen['tx']}/rx{chosen['rx']}"
        selected_counts[chain] = selected_counts.get(chain, 0) + 1

    if len(rows) < 10:
        raise SystemExit("not enough valid tx0/rx0 measurement groups")

    matrix = np.vstack(rows)
    normalized = normalize_rows(matrix)
    times_array = np.asarray(times, dtype=np.float64)
    return {
        "matrix": matrix,
        "normalized": normalized,
        "times_s": times_array - times_array[0],
        "outer_records": parsed["outer_records"],
        "closed_groups": len(parsed["groups"]),
        "rows": int(matrix.shape[0]),
        "subcarriers": int(matrix.shape[1]),
        "duration_s": rounded(times_array[-1] - times_array[0]),
        "malformed_records": parsed["malformed_records"],
        "incomplete_tail_records": parsed["incomplete_tail_records"],
        "missing_fixed_chain": missing_fixed_chain,
        "duplicate_fixed_chain": duplicate_fixed_chain,
        "invalid_fixed_chain": invalid_fixed_chain,
        "selected_chains": selected_counts,
    }


def normalize_rows(matrix: np.ndarray) -> np.ndarray:
    scale = np.median(matrix, axis=1)
    if np.any(scale <= 0):
        raise ValueError("row median must be positive")
    normalized = matrix / scale[:, None]
    if not np.all(np.isfinite(normalized)):
        raise ValueError("normalized matrix contains non-finite values")
    return normalized


def per_subcarrier_metrics(normalized: np.ndarray, raw: np.ndarray) -> dict:
    means = normalized.mean(axis=0)
    stds = normalized.std(axis=0)
    cv = np.divide(stds, means, out=np.zeros_like(stds), where=np.abs(means) > 1e-12)
    return {
        "median_profile": np.median(normalized, axis=0),
        "p05_profile": np.percentile(normalized, 5, axis=0),
        "p95_profile": np.percentile(normalized, 95, axis=0),
        "median_raw_amplitude": np.median(raw, axis=0),
        "zero_fraction": np.mean(raw == 0, axis=0),
        "temporal_cv": cv,
    }


def pearson_correlation(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    centered = matrix - matrix.mean(axis=0, keepdims=True)
    stds = centered.std(axis=0)
    active = stds > 1e-12
    if int(np.sum(active)) < 2:
        raise ValueError("fewer than two active subcarriers")
    standardized = centered[:, active] / stds[active]
    correlation = np.corrcoef(standardized, rowvar=False)
    correlation = np.nan_to_num(correlation, nan=0.0, posinf=0.0, neginf=0.0)
    return correlation, active


def dimension_metrics(normalized: np.ndarray) -> dict:
    correlation, active = pearson_correlation(normalized)
    eigenvalues = np.linalg.eigvalsh(correlation)[::-1]
    eigenvalues = np.clip(eigenvalues, 0.0, None)
    total = float(eigenvalues.sum())
    ratios = eigenvalues / total
    cumulative = np.cumsum(ratios)

    def first_at(threshold: float) -> int:
        return int(np.searchsorted(cumulative, threshold, side="left") + 1)

    off_diagonal = correlation[~np.eye(correlation.shape[0], dtype=bool)]
    abs_off = np.abs(off_diagonal)
    adjacent = np.diag(correlation, 1)
    participation_rank = float(total**2 / np.sum(eigenvalues**2))

    return {
        "correlation": correlation,
        "active_mask": active,
        "eigenvalues": eigenvalues,
        "variance_ratio": ratios,
        "cumulative": cumulative,
        "summary": {
            "active_subcarriers": int(np.sum(active)),
            "median_signed": rounded(np.median(off_diagonal)),
            "median_abs": rounded(np.median(abs_off)),
            "p95_abs": rounded(np.percentile(abs_off, 95)),
            "max_abs": rounded(np.max(abs_off)),
            "adjacent_signed_median": rounded(np.median(adjacent)),
            "adjacent_abs_mean": rounded(np.mean(np.abs(adjacent))),
            "fraction_abs_ge_0_75": rounded(np.mean(abs_off >= 0.75)),
            "fraction_abs_ge_0_9": rounded(np.mean(abs_off >= 0.9)),
        },
        "pca_dimensions": {
            "variance_90": first_at(0.90),
            "variance_95": first_at(0.95),
            "variance_99": first_at(0.99),
            "effective_rank": rounded(participation_rank),
            "first_component_variance": rounded(ratios[0]),
            "top_10_variance": rounded(np.sum(ratios[:10])),
        },
    }


def correlation_between(left: np.ndarray, right: np.ndarray) -> float:
    if left.size != right.size or left.size < 2:
        return 0.0
    value = float(np.corrcoef(left, right)[0, 1])
    return rounded(value) if math.isfinite(value) else 0.0


def cross_session_metrics(primary: dict, secondary: dict) -> dict:
    return {
        "profile_pearson": correlation_between(
            primary["metrics"]["median_profile"],
            secondary["metrics"]["median_profile"],
        ),
        "temporal_cv_pearson": correlation_between(
            primary["metrics"]["temporal_cv"],
            secondary["metrics"]["temporal_cv"],
        ),
        "median_profile_l1_distance": rounded(
            np.mean(
                np.abs(
                    primary["metrics"]["median_profile"]
                    - secondary["metrics"]["median_profile"]
                )
            )
        ),
    }


def write_metrics_csv(
    path: Path,
    primary: dict,
    secondary: dict,
) -> None:
    fields = [
        "subcarrier",
        "median_profile_10m",
        "p05_profile_10m",
        "p95_profile_10m",
        "temporal_cv_10m",
        "median_raw_amplitude_10m",
        "zero_fraction_10m",
        "median_profile_60s",
        "p05_profile_60s",
        "p95_profile_60s",
        "temporal_cv_60s",
        "median_raw_amplitude_60s",
        "zero_fraction_60s",
        "profile_delta_60s_minus_10m",
    ]
    with path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.DictWriter(handle, fieldnames=fields)
        writer.writeheader()
        for index in range(EXPECTED_SUBCARRIERS):
            pm = primary["metrics"]
            sm = secondary["metrics"]
            writer.writerow(
                {
                    "subcarrier": index,
                    "median_profile_10m": rounded(pm["median_profile"][index], 9),
                    "p05_profile_10m": rounded(pm["p05_profile"][index], 9),
                    "p95_profile_10m": rounded(pm["p95_profile"][index], 9),
                    "temporal_cv_10m": rounded(pm["temporal_cv"][index], 9),
                    "median_raw_amplitude_10m": rounded(
                        pm["median_raw_amplitude"][index], 6
                    ),
                    "zero_fraction_10m": rounded(pm["zero_fraction"][index], 9),
                    "median_profile_60s": rounded(sm["median_profile"][index], 9),
                    "p05_profile_60s": rounded(sm["p05_profile"][index], 9),
                    "p95_profile_60s": rounded(sm["p95_profile"][index], 9),
                    "temporal_cv_60s": rounded(sm["temporal_cv"][index], 9),
                    "median_raw_amplitude_60s": rounded(
                        sm["median_raw_amplitude"][index], 6
                    ),
                    "zero_fraction_60s": rounded(sm["zero_fraction"][index], 9),
                    "profile_delta_60s_minus_10m": rounded(
                        sm["median_profile"][index] - pm["median_profile"][index], 9
                    ),
                }
            )


def write_pca_csv(path: Path, primary_dimensions: dict) -> None:
    with path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.writer(handle)
        writer.writerow(
            [
                "component",
                "eigenvalue",
                "variance_ratio",
                "cumulative_variance",
            ]
        )
        for index, (eigenvalue, ratio, cumulative) in enumerate(
            zip(
                primary_dimensions["eigenvalues"],
                primary_dimensions["variance_ratio"],
                primary_dimensions["cumulative"],
            ),
            start=1,
        ):
            writer.writerow(
                [
                    index,
                    rounded(eigenvalue, 9),
                    rounded(ratio, 12),
                    rounded(cumulative, 12),
                ]
            )


def svg_text(
    x: float,
    y: float,
    value: object,
    size: int = 12,
    color: str = "#dbeafe",
    anchor: str = "start",
    weight: int = 400,
) -> str:
    return (
        f'<text x="{x:.2f}" y="{y:.2f}" fill="{color}" font-size="{size}" '
        f'font-family="ui-sans-serif,system-ui,sans-serif" '
        f'text-anchor="{anchor}" font-weight="{weight}">{html.escape(str(value))}</text>'
    )


def svg_line(
    x1: float,
    y1: float,
    x2: float,
    y2: float,
    color: str = "#475569",
    width: float = 1,
) -> str:
    return (
        f'<line x1="{x1:.2f}" y1="{y1:.2f}" x2="{x2:.2f}" y2="{y2:.2f}" '
        f'stroke="{color}" stroke-width="{width}"/>'
    )


def polyline(
    values: np.ndarray,
    x: float,
    y: float,
    width: float,
    height: float,
    y_min: float,
    y_max: float,
    color: str,
    max_points: int = 256,
) -> str:
    array = np.asarray(values, dtype=np.float64)
    if len(array) > max_points:
        index = np.linspace(0, len(array) - 1, max_points).astype(int)
        array = array[index]
    x_values = x + np.linspace(0, width, len(array))
    y_values = y + height * (1.0 - (array - y_min) / (y_max - y_min))
    points = " ".join(f"{xv:.2f},{yv:.2f}" for xv, yv in zip(x_values, y_values))
    return (
        f'<polyline points="{points}" fill="none" stroke="{color}" '
        f'stroke-width="2.2" stroke-linejoin="round"/>'
    )


def correlation_color(value: float) -> str:
    value = max(-1.0, min(1.0, value))
    if value >= 0:
        low = np.array([15, 23, 42], dtype=float)
        high = np.array([251, 191, 36], dtype=float)
        rgb = low + (high - low) * value
    else:
        low = np.array([15, 23, 42], dtype=float)
        high = np.array([14, 165, 233], dtype=float)
        rgb = low + (high - low) * abs(value)
    return "#" + "".join(f"{int(channel):02x}" for channel in rgb)


def build_svg(
    path: Path,
    primary: dict,
    secondary: dict,
) -> None:
    width, height = 1120, 760
    panel_w, panel_h = 470, 270
    left_x, right_x = 70, 610
    top_y, bottom_y = 90, 430
    parts = [
        f'<svg xmlns="http://www.w3.org/2000/svg" width="{width}" height="{height}" '
        f'viewBox="0 0 {width} {height}" role="img" '
        f'aria-label="Subcarrier redundancy and effective-dimensionality analysis">',
        '<rect width="100%" height="100%" fill="#0f172a"/>',
        svg_text(40, 38, "256 CSI subcarriers are not 256 independent dimensions", 24, "#f8fafc", weight=700),
        svg_text(
            40,
            63,
            "GL-MT6000 tx0/rx0 amplitude; rows normalised by their 256-bin median",
            13,
            "#94a3b8",
        ),
    ]

    panels = [
        (left_x, top_y, "Normalised amplitude profile"),
        (right_x, top_y, "Subcarrier correlation (4 x 4 block mean)"),
        (left_x, bottom_y, "PCA variance on the 10-minute capture"),
        (right_x, bottom_y, "Absolute off-diagonal correlation"),
    ]
    for px, py, title in panels:
        parts.append(f'<rect x="{px}" y="{py}" width="{panel_w}" height="{panel_h}" fill="#111c33" rx="8"/>')
        parts.append(svg_text(px + 14, py + 24, title, 15, "#e2e8f0", weight=650))

    # Panel 1: median normalised subcarrier profile in both sessions.
    px, py = left_x, top_y
    plot_x, plot_y, plot_w, plot_h = px + 45, py + 48, panel_w - 65, panel_h - 86
    y_min, y_max = 0.0, max(
        float(np.percentile(primary["metrics"]["median_profile"], 98)),
        float(np.percentile(secondary["metrics"]["median_profile"], 98)),
        1.0,
    ) * 1.12
    for tick in np.linspace(y_min, y_max, 5):
        y = plot_y + plot_h * (1.0 - (tick - y_min) / (y_max - y_min))
        parts.append(svg_line(plot_x, y, plot_x + plot_w, y))
        parts.append(svg_text(plot_x - 8, y + 4, f"{tick:.1f}", 10, "#94a3b8", "end"))
    parts.append(polyline(primary["metrics"]["median_profile"], plot_x, plot_y, plot_w, plot_h, y_min, y_max, "#38bdf8"))
    parts.append(polyline(secondary["metrics"]["median_profile"], plot_x, plot_y, plot_w, plot_h, y_min, y_max, "#fbbf24"))
    parts.append(svg_line(plot_x, plot_y, plot_x, plot_y + plot_h, "#94a3b8"))
    parts.append(svg_line(plot_x, plot_y + plot_h, plot_x + plot_w, plot_y + plot_h, "#94a3b8"))
    parts.append(svg_text(plot_x, py + panel_h - 13, "subcarrier index 0", 10, "#94a3b8"))
    parts.append(svg_text(plot_x + plot_w, py + panel_h - 13, "255", 10, "#94a3b8", "end"))
    parts.append(svg_line(px + 330, py + 27, px + 355, py + 27, "#38bdf8", 3))
    parts.append(svg_text(px + 362, py + 31, "10 min", 10, "#e2e8f0"))
    parts.append(svg_line(px + 405, py + 27, px + 430, py + 27, "#fbbf24", 3))
    parts.append(svg_text(px + 437, py + 31, "60 s", 10, "#e2e8f0"))

    # Panel 2: block-averaged correlation matrix.
    px, py = right_x, top_y
    correlation = primary["dimensions"]["correlation"]
    block = 4
    blocks = correlation.shape[0] // block
    matrix_x, matrix_y, matrix_size = px + 35, py + 42, panel_w - 70
    cell = matrix_size / blocks
    for row in range(blocks):
        for col in range(blocks):
            value = float(
                correlation[
                    row * block : (row + 1) * block,
                    col * block : (col + 1) * block,
                ].mean()
            )
            parts.append(
                f'<rect x="{matrix_x + col * cell:.2f}" y="{matrix_y + row * cell:.2f}" '
                f'width="{cell + 0.2:.2f}" height="{cell + 0.2:.2f}" '
                f'fill="{correlation_color(value)}"/>'
            )
    parts.append(svg_text(matrix_x, py + panel_h - 13, "correlation -1", 10, "#94a3b8"))
    parts.append(svg_text(matrix_x + matrix_size, py + panel_h - 13, "+1", 10, "#94a3b8", "end"))

    # Panel 3: scree/cumulative variance.
    px, py = left_x, bottom_y
    plot_x, plot_y, plot_w, plot_h = px + 50, py + 45, panel_w - 75, panel_h - 83
    components = min(30, len(primary["dimensions"]["cumulative"]))
    cumulative = primary["dimensions"]["cumulative"][:components]
    parts.append(svg_line(plot_x, plot_y, plot_x, plot_y + plot_h, "#94a3b8"))
    parts.append(svg_line(plot_x, plot_y + plot_h, plot_x + plot_w, plot_y + plot_h, "#94a3b8"))
    for threshold in (0.90, 0.95, 0.99):
        y = plot_y + plot_h * (1.0 - threshold)
        parts.append(svg_line(plot_x, y, plot_x + plot_w, y, "#475569"))
        parts.append(svg_text(plot_x - 8, y + 4, f"{threshold:.0%}", 9, "#94a3b8", "end"))
    parts.append(polyline(cumulative, plot_x, plot_y, plot_w, plot_h, 0.0, 1.0, "#34d399"))
    dims = primary["dimensions"]["pca_dimensions"]
    parts.append(
        svg_text(
            plot_x + plot_w,
            py + panel_h - 13,
            f"90%={dims['variance_90']}  95%={dims['variance_95']}  99%={dims['variance_99']}",
            10,
            "#e2e8f0",
            "end",
        )
    )
    parts.append(svg_text(plot_x, py + 38, "first 30 components", 10, "#94a3b8"))

    # Panel 4: histogram of absolute correlations.
    px, py = right_x, bottom_y
    plot_x, plot_y, plot_w, plot_h = px + 45, py + 48, panel_w - 65, panel_h - 87
    corr = primary["dimensions"]["correlation"]
    off = np.abs(corr[~np.eye(corr.shape[0], dtype=bool)])
    counts, edges = np.histogram(off, bins=50, range=(0.0, 1.0))
    max_count = max(int(counts.max()), 1)
    bar_w = plot_w / len(counts)
    for index, count in enumerate(counts):
        bar_h = plot_h * count / max_count
        parts.append(
            f'<rect x="{plot_x + index * bar_w:.2f}" '
            f'y="{plot_y + plot_h - bar_h:.2f}" width="{bar_w * 0.88:.2f}" '
            f'height="{bar_h:.2f}" fill="#a78bfa"/>'
        )
    parts.append(svg_line(plot_x, plot_y + plot_h, plot_x + plot_w, plot_y + plot_h, "#94a3b8"))
    parts.append(svg_text(plot_x, py + panel_h - 13, "|r| = 0", 10, "#94a3b8"))
    parts.append(svg_text(plot_x + plot_w, py + panel_h - 13, "1", 10, "#94a3b8", "end"))
    corr_summary = primary["dimensions"]["summary"]
    parts.append(
        svg_text(
            plot_x + plot_w,
            py + 38,
            f"median |r|={corr_summary['median_abs']:.3f}",
            10,
            "#e2e8f0",
            "end",
        )
    )

    parts.append("</svg>")
    path.write_text("\n".join(parts) + "\n", encoding="utf-8")


def build_capture_summary(capture: dict) -> dict:
    return {
        "outer_records": capture["outer_records"],
        "closed_groups": capture["closed_groups"],
        "rows": capture["rows"],
        "subcarriers": capture["subcarriers"],
        "duration_s": capture["duration_s"],
        "malformed_records": capture["malformed_records"],
        "incomplete_tail_records": capture["incomplete_tail_records"],
        "missing_fixed_chain": capture["missing_fixed_chain"],
        "duplicate_fixed_chain": capture["duplicate_fixed_chain"],
        "invalid_fixed_chain": capture["invalid_fixed_chain"],
        "selected_chains": capture["selected_chains"],
    }


def analysis_for_capture(capture: dict) -> dict:
    metrics = per_subcarrier_metrics(capture["normalized"], capture["matrix"])
    dimensions = dimension_metrics(capture["normalized"])
    return {**capture, "metrics": metrics, "dimensions": dimensions}


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("capture_10m", type=Path)
    parser.add_argument("capture_60s", type=Path)
    parser.add_argument("--out-dir", type=Path, required=True)
    parser.add_argument("--img", type=Path, required=True)
    args = parser.parse_args()

    args.out_dir.mkdir(parents=True, exist_ok=True)
    args.img.parent.mkdir(parents=True, exist_ok=True)

    primary = analysis_for_capture(select_fixed_chain(parse_capture(args.capture_10m)))
    secondary = analysis_for_capture(select_fixed_chain(parse_capture(args.capture_60s)))

    vicinity = cross_session_metrics(primary, secondary)
    summary = {
        "source_10m": {
            "file": args.capture_10m.name,
            "size_bytes": args.capture_10m.stat().st_size,
            "sha256": sha256_file(args.capture_10m),
        },
        "source_60s": {
            "file": args.capture_60s.name,
            "size_bytes": args.capture_60s.stat().st_size,
            "sha256": sha256_file(args.capture_60s),
        },
        "capture_10m": build_capture_summary(primary),
        "capture_60s": build_capture_summary(secondary),
        "subcarrier_correlation": primary["dimensions"]["summary"],
        "pca_dimensions": primary["dimensions"]["pca_dimensions"],
        "cross_session": vicinity,
        "subcarrier_temporal_cv_10m": {
            "min": rounded(np.min(primary["metrics"]["temporal_cv"])),
            "p25": rounded(np.percentile(primary["metrics"]["temporal_cv"], 25)),
            "median": rounded(np.median(primary["metrics"]["temporal_cv"])),
            "p75": rounded(np.percentile(primary["metrics"]["temporal_cv"], 75)),
            "p95": rounded(np.percentile(primary["metrics"]["temporal_cv"], 95)),
            "max": rounded(np.max(primary["metrics"]["temporal_cv"])),
            "fraction_above_1": rounded(
                np.mean(primary["metrics"]["temporal_cv"] > 1.0)
            ),
        },
        "subcarrier_temporal_cv_60s": {
            "min": rounded(np.min(secondary["metrics"]["temporal_cv"])),
            "p25": rounded(np.percentile(secondary["metrics"]["temporal_cv"], 25)),
            "median": rounded(np.median(secondary["metrics"]["temporal_cv"])),
            "p75": rounded(np.percentile(secondary["metrics"]["temporal_cv"], 75)),
            "p95": rounded(np.percentile(secondary["metrics"]["temporal_cv"], 95)),
            "max": rounded(np.max(secondary["metrics"]["temporal_cv"])),
            "fraction_above_1": rounded(
                np.mean(secondary["metrics"]["temporal_cv"] > 1.0)
            ),
        },
        "method": {
            "group_boundary": "chain_info bit 15",
            "chain_selection": "tx0/rx0 only; no fallback chain",
            "row_normalisation": "divide each 256-bin amplitude vector by its median",
            "correlation": "Pearson correlation after centring and scaling each subcarrier over time",
            "pca": "eigendecomposition of the 256 x 256 subcarrier correlation matrix",
            "dimension_threshold": "first component count reaching 90%, 95%, or 99% cumulative variance",
            "effective_rank": "participation ratio: (sum eigenvalues)^2 / sum(eigenvalues^2)",
        },
        "limitations": [
            "The subcarrier array index is not mapped to a physical frequency, antenna position, or 802.11 resource unit.",
            "High amplitude correlation can come from the radio front end, shared gain, multipath, or the data format; correlation alone does not establish a causal source.",
            "The two captures come from one GL-MT6000, one client, one home, one bandwidth, and similar firmware; they do not establish cross-device generality.",
            "There is no noise-only baseline, controlled subject position, or human-event ground truth.",
            "PCA variance depends on row normalisation and standardisation choices; the published matrix and script are required for comparison.",
            "The result supports data-redundancy checks and dimension reporting, not breathing, fall, occupancy, identity, or medical claims.",
        ],
    }

    (args.out_dir / "summary.json").write_text(
        json.dumps(summary, ensure_ascii=False, indent=2) + "\n",
        encoding="utf-8",
    )
    write_metrics_csv(args.out_dir / "subcarrier_metrics.csv", primary, secondary)
    write_pca_csv(args.out_dir / "pca_variance.csv", primary["dimensions"])
    build_svg(args.img, primary, secondary)

    print(f"10-minute rows: {primary['rows']} x {primary['subcarriers']}")
    print(f"60-second rows: {secondary['rows']} x {secondary['subcarriers']}")
    print(
        "PCA dimensions: "
        f"90%={summary['pca_dimensions']['variance_90']}, "
        f"95%={summary['pca_dimensions']['variance_95']}, "
        f"99%={summary['pca_dimensions']['variance_99']}, "
        f"effective rank={summary['pca_dimensions']['effective_rank']}"
    )
    print(
        "Median |correlation|: "
        f"{summary['subcarrier_correlation']['median_abs']}"
    )
    print(f"Wrote {args.out_dir / 'summary.json'}")
    print(f"Wrote {args.out_dir / 'subcarrier_metrics.csv'}")
    print(f"Wrote {args.out_dir / 'pca_variance.csv'}")
    print(f"Wrote {args.img}")


if __name__ == "__main__":
    main()