From 47208b65ced3de07599be26b46dca7c72883c79d Mon Sep 17 00:00:00 2001 From: KoyoNoguchi Date: Wed, 14 Oct 2026 19:06:18 +0900 Subject: [PATCH] =?UTF-8?q?=E7=B2=BE=E5=BA=A6=E8=A9=95=E4=BE=A1=E5=AE=9F?= =?UTF-8?q?=E9=A8=93=E3=81=AE=E5=88=86=E6=9E=90=E3=82=92=E8=A4=87=E6=95=B0?= =?UTF-8?q?=E8=A2=AB=E9=A8=93=E8=80=85=E5=AF=BE=E5=BF=9C?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- config/experiment.yaml | 2 +- src/accuracy_report.py | 516 +++++++++++++++++++++++++++-- src/accuracy_report.py.orig | 475 ++++++++++++++++++++++++++ src/accuracy_session.py | 21 +- tests/test_accuracy_report.py | 119 ++++++- tests/test_accuracy_report.py.orig | 139 ++++++++ 6 files changed, 1226 insertions(+), 46 deletions(-) create mode 100644 src/accuracy_report.py.orig create mode 100644 tests/test_accuracy_report.py.orig diff --git a/config/experiment.yaml b/config/experiment.yaml index 3c6a1d8..cc1aa13 100644 --- a/config/experiment.yaml +++ b/config/experiment.yaml @@ -1,5 +1,5 @@ experiment: - experiment_id: phase1_test_001 + experiment_id: test location: laboratory operator: student notes: "Phase 2 pedestrian detection and tracking test" diff --git a/src/accuracy_report.py b/src/accuracy_report.py index ec5b323..850ae12 100644 --- a/src/accuracy_report.py +++ b/src/accuracy_report.py @@ -14,6 +14,26 @@ Threshold: ``theta_th = |bias| + k * SD`` of the frame errors for conditions wit ``|true_yaw| <= --pool-max-angle`` (pitch 0), per distance. ``--basis worst`` (default) takes the worst single condition, ``pooled`` pools all of them. k=3 and the pooling range are provisional choices to be agreed on with the supervisor. + +Multiple subjects +----------------- +Give every subject's marker file(s); the subject name is taken from, in this order, a +``LABEL=path`` prefix on ``--trials``, the ``subject`` column written by ``accuracy_session plan +--subject``, or the marker file name:: + + python -m src.accuracy_report --trials A=data/analysis/accuracy/trials_a.csv ^ + B=data/analysis/accuracy/trials_b.csv C=data/analysis/accuracy/trials_c.csv ^ + --face-csv "data/logs/face/*_face_frames.csv" + +The thresholds, regression and ``accuracy_conditions.csv`` pool the frames of all subjects (so the +SD contains the between-subject bias differences, which is how individual bias is averaged out). +With two or more subjects the report additionally contains per-subject thresholds +(``accuracy_subject_thresholds.csv`` / ``accuracy_subject_conditions.csv``), a leave-one-subject-out +check (``accuracy_loso.csv``), the +/- direction asymmetry of the gain (``accuracy_asymmetry*.csv``), +a between-subject comparison of the gain (``accuracy_subject_differences.csv``) and per-subject plots. +The face logs of different subjects must not overlap in time (they are joined by timestamp). +``--exclude-subjects``, ``--subjects``, ``--alias OLD=NEW``, ``--flip-subject NAME`` and +``--exclude-seq 12 S01:7`` (seq of one subject only) adjust what is analysed. """ from __future__ import annotations @@ -23,23 +43,33 @@ import bisect import csv import glob import math +import random import statistics import sys from collections import Counter from dataclasses import dataclass from datetime import datetime from pathlib import Path -from typing import Sequence +from typing import Mapping, Sequence NAN = float("nan") TRIAL_COLUMNS = ["session_id", "seq", "condition_id", "distance_m", "true_yaw_deg", "pitch_deg", "repeat", "subject", "status", "track_id", "baseline_n", "baseline_mean", "baseline_sd", "hold_n", "hold_mean", "hold_sd", "delta_deg", "error_deg", "detection_rate"] -CONDITION_COLUMNS = ["distance_m", "pitch_deg", "true_yaw_deg", "n_trials", "n_frames", "bias_deg", "sd_trials_deg", - "sd_within_deg", "sd_total_deg", "rmse_frame_deg", "p95_abs_error_deg", "mean_detection_rate", - "exceed_rate"] -THRESHOLD_COLUMNS = ["distance_m", "n_frames", "bias_deg", "sd_deg", "theta_pooled_deg", "theta_worst_deg", - "theta_th_deg", "slope", "intercept_deg", "r2", "min_detectable_deg"] +CONDITION_BASE_COLUMNS = ["distance_m", "pitch_deg", "true_yaw_deg", "n_trials", "n_frames", "bias_deg", + "sd_trials_deg", "sd_within_deg", "sd_total_deg", "rmse_frame_deg", "p95_abs_error_deg", + "mean_detection_rate", "exceed_rate"] +CONDITION_COLUMNS = CONDITION_BASE_COLUMNS + ["n_subjects", "sd_subjects_deg"] +THRESHOLD_BASE_COLUMNS = ["distance_m", "n_frames", "bias_deg", "sd_deg", "theta_pooled_deg", "theta_worst_deg", + "theta_th_deg", "slope", "intercept_deg", "r2", "min_detectable_deg"] +THRESHOLD_COLUMNS = THRESHOLD_BASE_COLUMNS + ["n_subjects"] +SUBJECT_CONDITION_COLUMNS = ["subject"] + CONDITION_BASE_COLUMNS +SUBJECT_THRESHOLD_COLUMNS = ["subject"] + THRESHOLD_BASE_COLUMNS +LOSO_COLUMNS = ["held_out", "distance_m", "n_train_subjects", "theta_train_deg", "false_alarm_rate", + "min_detectable_deg"] +ASYMMETRY_CELL_COLUMNS = ["subject", "distance_m", "abs_angle_deg", "gain_pos", "gain_neg", "gain_diff"] +ASYMMETRY_COLUMNS = ["subject", "n_cells", "gain_pos", "gain_neg", "gain_diff", "p_signflip"] +DIFFERENCE_COLUMNS = ["distance_m", "n_subjects", "kruskal_h", "p_permutation"] @dataclass(frozen=True) @@ -57,6 +87,7 @@ class Options: basis: str = "worst" detect_ratio: float = 0.95 pitch_ref: float = 0.0 + flip_subjects: tuple[str, ...] = () @dataclass @@ -113,6 +144,14 @@ def expand_paths(patterns: Sequence[str]) -> list[Path]: return paths +def split_label(token: str) -> tuple[str, str]: + """``LABEL=path`` -> (LABEL, path); anything else -> ("", token).""" + label, sep, rest = token.partition("=") + if sep and label and rest and not any(c in label for c in "\\/:*?[]"): + return label, rest + return "", token + + def load_face_series(paths: Sequence[Path], default_fps: float = 30.0) -> FaceSeries: rows: list[tuple[float, int, bool, float | None, float | None, float]] = [] for path in paths: @@ -139,14 +178,38 @@ def load_face_series(paths: Sequence[Path], default_fps: float = 30.0) -> FaceSe [r[4] for r in rows], [r[5] for r in rows]) -def load_trials(paths: Sequence[Path], exclude_seqs: Sequence[int] = ()) -> list[dict[str, str]]: - """Last log row per (session, seq); only status ``ok`` survives.""" - last: dict[tuple[str, int], dict[str, str]] = {} +def _parse_exclusions(tokens: Sequence[int | str]) -> tuple[set[int], set[tuple[str, int]]]: + """``12`` excludes seq 12 of every subject, ``S01:12`` only that of subject S01.""" + general: set[int] = set() + specific: set[tuple[str, int]] = set() + for token in tokens: + text = str(token) + if ":" in text: + name, _, number = text.rpartition(":") + specific.add((name, int(number))) + else: + general.add(int(text)) + return general, specific + + +def load_trials(paths: Sequence[Path], exclude_seqs: Sequence[int | str] = (), + labels: Mapping[Path, str] | None = None, + aliases: Mapping[str, str] | None = None) -> list[dict[str, str]]: + """Last log row per (subject, session, seq); only status ``ok`` survives. + + The subject is the ``LABEL`` given for the file, else the ``subject`` column, else the file name. + """ + labels, aliases = labels or {}, aliases or {} + last: dict[tuple[str, str, int], dict[str, str]] = {} for path in paths: with path.open(encoding="utf-8-sig", newline="") as stream: for row in csv.DictReader(stream): - last[(row.get("session_id", ""), int(row["seq"]))] = row - return [row for (_, seq), row in sorted(last.items()) if row["status"] == "ok" and seq not in set(exclude_seqs)] + subject = labels.get(path) or (row.get("subject") or "").strip() or path.stem + row["subject"] = aliases.get(subject, subject) + last[(row["subject"], row.get("session_id", ""), int(row["seq"]))] = row + general, specific = _parse_exclusions(exclude_seqs) + return [row for (subject, _, seq), row in sorted(last.items()) + if row["status"] == "ok" and seq not in general and (subject, seq) not in specific] def evaluate_trial(trial: dict[str, str], series: FaceSeries, opt: Options) -> tuple[dict[str, object], list[float]]: @@ -191,7 +254,8 @@ def evaluate_trial(trial: dict[str, str], series: FaceSeries, opt: Options) -> t if len(hold) < opt.min_frames: record["status"] = "hold_insufficient" return record, [] - sign = -1.0 if opt.flip_sign else 1.0 + flipped = opt.flip_sign or str(record["subject"]) in opt.flip_subjects + sign = -1.0 if flipped else 1.0 base_mean = float(record["baseline_mean"]) # type: ignore[arg-type] deltas = [_wrap(sign * (y - base_mean)) for y in hold_yaw] fps = opt.fps or series.fps[hold[0]] @@ -209,6 +273,17 @@ def _errors(record: dict[str, object], deltas: Sequence[float]) -> list[float]: return [d - float(record["true_yaw_deg"]) for d in deltas] # type: ignore[arg-type] +def _subject(record: dict[str, object]) -> str: + return str(record.get("subject") or "") + + +def split_by_subject(usable: Usable) -> dict[str, Usable]: + groups: dict[str, Usable] = {} + for item in usable: + groups.setdefault(_subject(item[0]), []).append(item) + return groups + + def summarize_conditions(usable: Usable) -> list[dict[str, object]]: groups: dict[tuple[float, float, float], Usable] = {} for record, deltas in usable: @@ -218,6 +293,9 @@ def summarize_conditions(usable: Usable) -> list[dict[str, object]]: for (distance, pitch, true), items in sorted(groups.items()): frame_errors = [e for record, deltas in items for e in _errors(record, deltas)] within = [float(r["hold_sd"]) ** 2 for r, _ in items if not math.isnan(float(r["hold_sd"]))] # type: ignore[arg-type] + by_subject: dict[str, list[float]] = {} + for r, _ in items: + by_subject.setdefault(_subject(r), []).append(float(r["error_deg"])) # type: ignore[arg-type] rows.append({ "distance_m": distance, "pitch_deg": pitch, "true_yaw_deg": true, "n_trials": len(items), "n_frames": len(frame_errors), "bias_deg": _mean([float(r["error_deg"]) for r, _ in items]), # type: ignore[arg-type] @@ -226,7 +304,8 @@ def summarize_conditions(usable: Usable) -> list[dict[str, object]]: "rmse_frame_deg": math.sqrt(_mean([e * e for e in frame_errors])), "p95_abs_error_deg": percentile(sorted(abs(e) for e in frame_errors), 95), "mean_detection_rate": _mean([float(r["detection_rate"]) for r, _ in items]), # type: ignore[arg-type] - "exceed_rate": NAN, + "exceed_rate": NAN, "n_subjects": len(by_subject), + "sd_subjects_deg": _sd([_mean(v) for v in by_subject.values()]), }) return rows @@ -264,6 +343,7 @@ def compute_thresholds(usable: Usable, opt: Options) -> dict[float, dict[str, fl "theta_pooled_deg": pooled, "theta_worst_deg": worst, "theta_th_deg": worst if opt.basis == "worst" else pooled, "slope": slope, "intercept_deg": intercept, "r2": r2, "min_detectable_deg": NAN, + "n_subjects": len({_subject(r) for r, _ in pooled_sel}), } return result @@ -288,6 +368,162 @@ def add_exceed_rates(conditions: list[dict[str, object]], usable: Usable, thresh break +# --------------------------------------------------------------------------------------------- +# Multi-subject analysis +# --------------------------------------------------------------------------------------------- +def subject_analysis(usable: Usable, opt: Options) -> dict[str, dict[str, object]]: + """Conditions and thresholds computed from each subject alone.""" + result: dict[str, dict[str, object]] = {} + for subject, items in sorted(split_by_subject(usable).items()): + conditions = summarize_conditions(items) + thresholds = compute_thresholds(items, opt) + add_exceed_rates(conditions, items, thresholds, opt) + result[subject] = {"conditions": conditions, "thresholds": thresholds} + return result + + +def leave_one_subject_out(usable: Usable, opt: Options) -> list[dict[str, object]]: + """Threshold made from the other subjects, applied to the held-out one (false alarm at 0 deg, min angle).""" + groups = split_by_subject(usable) + if len(groups) < 2: + return [] + rows: list[dict[str, object]] = [] + for held, test in sorted(groups.items()): + train = [item for subject, items in groups.items() if subject != held for item in items] + trained = compute_thresholds(train, opt) + carried = {d: {"theta_th_deg": t["theta_th_deg"], "min_detectable_deg": NAN} for d, t in trained.items()} + conditions = summarize_conditions(test) + add_exceed_rates(conditions, test, carried, opt) + for distance, info in sorted(carried.items()): + zero = [c for c in conditions if c["distance_m"] == distance and c["pitch_deg"] == opt.pitch_ref + and float(c["true_yaw_deg"]) == 0.0] # type: ignore[arg-type] + if not any(c["distance_m"] == distance for c in conditions): + continue + rows.append({"held_out": held, "distance_m": distance, + "n_train_subjects": len({_subject(r) for r, _ in train}), + "theta_train_deg": info["theta_th_deg"], + "false_alarm_rate": zero[0]["exceed_rate"] if zero else NAN, + "min_detectable_deg": info["min_detectable_deg"]}) + return rows + + +def gain_records(usable: Usable, opt: Options) -> list[dict[str, object]]: + """Per-trial gain = (estimated change toward the set direction) / |set angle|; 1 is ideal.""" + rows: list[dict[str, object]] = [] + for record, _ in usable: + true = float(record["true_yaw_deg"]) # type: ignore[arg-type] + if float(record["pitch_deg"]) != opt.pitch_ref or true == 0.0: # type: ignore[arg-type] + continue + rows.append({"subject": _subject(record), "distance_m": float(record["distance_m"]), # type: ignore[arg-type] + "abs_angle_deg": abs(true), "side": "pos" if true > 0 else "neg", + "gain": float(record["delta_deg"]) * math.copysign(1.0, true) / abs(true)}) # type: ignore[arg-type] + return rows + + +def asymmetry_cells(gains: Sequence[dict[str, object]]) -> list[dict[str, object]]: + cells: dict[tuple[str, float, float], dict[str, list[float]]] = {} + for g in gains: + key = (str(g["subject"]), float(g["distance_m"]), float(g["abs_angle_deg"])) # type: ignore[arg-type] + cells.setdefault(key, {"pos": [], "neg": []})[str(g["side"])].append(float(g["gain"])) # type: ignore[arg-type] + rows: list[dict[str, object]] = [] + for (subject, distance, angle), sides in sorted(cells.items()): + if sides["pos"] and sides["neg"]: + pos, neg = _mean(sides["pos"]), _mean(sides["neg"]) + rows.append({"subject": subject, "distance_m": distance, "abs_angle_deg": angle, + "gain_pos": pos, "gain_neg": neg, "gain_diff": pos - neg}) + return rows + + +def _signflip_p(diffs: Sequence[float], *, seed: int = 0, exact_max: int = 16, samples: int = 20000) -> float: + """Two-sided paired sign-flip permutation test of mean(diffs) = 0 (exact for small n).""" + n = len(diffs) + if n < 2 or all(d == 0 for d in diffs): + return NAN + observed, eps = abs(sum(diffs)), 1e-12 + if n <= exact_max: + hits = sum(abs(sum(d if (mask >> i) & 1 else -d for i, d in enumerate(diffs))) >= observed - eps + for mask in range(1 << n)) + return hits / (1 << n) + rng = random.Random(seed) + hits = sum(abs(sum(d if rng.random() < 0.5 else -d for d in diffs)) >= observed - eps for _ in range(samples)) + return (hits + 1) / (samples + 1) + + +def asymmetry_summary(cells: Sequence[dict[str, object]]) -> list[dict[str, object]]: + """Per-subject and overall +/- gain difference. The p value treats cells as independent (reference only).""" + rows: list[dict[str, object]] = [] + subjects = sorted({str(c["subject"]) for c in cells}) + for name, selected in [(s, [c for c in cells if c["subject"] == s]) for s in subjects] + [("ALL", list(cells))]: + if not selected: + continue + diffs = [float(c["gain_diff"]) for c in selected] # type: ignore[arg-type] + rows.append({"subject": name, "n_cells": len(selected), + "gain_pos": _mean([float(c["gain_pos"]) for c in selected]), # type: ignore[arg-type] + "gain_neg": _mean([float(c["gain_neg"]) for c in selected]), # type: ignore[arg-type] + "gain_diff": _mean(diffs), "p_signflip": _signflip_p(diffs)}) + return rows + + +def _ranks(values: Sequence[float]) -> list[float]: + order = sorted(range(len(values)), key=lambda i: values[i]) + ranks = [0.0] * len(values) + i = 0 + while i < len(order): + j = i + while j + 1 < len(order) and values[order[j + 1]] == values[order[i]]: + j += 1 + for k in range(i, j + 1): + ranks[order[k]] = (i + j) / 2.0 + 1.0 + i = j + 1 + return ranks + + +def _kruskal_h(rank_groups: Sequence[Sequence[float]]) -> float: + total = sum(len(g) for g in rank_groups) + return 12.0 / (total * (total + 1)) * sum(sum(g) ** 2 / len(g) for g in rank_groups) - 3.0 * (total + 1) + + +def kruskal_permutation(groups: Sequence[Sequence[float]], *, seed: int = 0, samples: int = 5000) -> tuple[float, float]: + """Kruskal-Wallis H and its permutation p value (no scipy needed, valid with ties).""" + flat = [v for g in groups for v in g] + ranks = _ranks(flat) + sizes = [len(g) for g in groups] + + def split(values: Sequence[float]) -> list[list[float]]: + out, start = [], 0 + for size in sizes: + out.append(list(values[start:start + size])) + start += size + return out + + observed = _kruskal_h(split(ranks)) + rng, hits = random.Random(seed), 0 + shuffled = list(ranks) + for _ in range(samples): + rng.shuffle(shuffled) + hits += _kruskal_h(split(shuffled)) >= observed - 1e-12 + return observed, (hits + 1) / (samples + 1) + + +def subject_differences(gains: Sequence[dict[str, object]], min_abs_angle: float = 10.0) -> list[dict[str, object]]: + """Per distance: does the gain (|set angle| >= 10 deg) differ between subjects?""" + rows: list[dict[str, object]] = [] + for distance in sorted({float(g["distance_m"]) for g in gains}): # type: ignore[arg-type] + by_subject: dict[str, list[float]] = {} + for g in gains: + if float(g["distance_m"]) == distance and float(g["abs_angle_deg"]) >= min_abs_angle: # type: ignore[arg-type] + by_subject.setdefault(str(g["subject"]), []).append(float(g["gain"])) # type: ignore[arg-type] + groups = [v for v in by_subject.values() if len(v) >= 2] + if len(groups) < 2: + continue + h, p = kruskal_permutation(groups) + rows.append({"distance_m": distance, "n_subjects": len(groups), "kruskal_h": h, "p_permutation": p}) + return rows + + +# --------------------------------------------------------------------------------------------- +# Output +# --------------------------------------------------------------------------------------------- def _fmt(value: object, digits: int = 2) -> str: if isinstance(value, float): return "-" if math.isnan(value) else f"{value:.{digits}f}" @@ -303,18 +539,24 @@ def _write_csv(path: Path, columns: Sequence[str], rows: Sequence[dict[str, obje def build_markdown(records: list[dict[str, object]], conditions: list[dict[str, object]], - thresholds: dict[float, dict[str, float]], opt: Options) -> str: + thresholds: dict[float, dict[str, float]], opt: Options, + multi: dict[str, object] | None = None) -> str: lines = ["# 振り向き角度推定の精度実験 結果", ""] status = Counter(str(r["status"]) for r in records) lines += [f"- 記録された有効試行: {len(records)}、解析に使用: {status['ok']}", "- 除外: " + ( ", ".join(f"{k}={v}" for k, v in status.items() if k != "ok") or "なし")] - bad = [str(r["seq"]) for r in records if r["status"] != "ok"] + bad = [str(r["seq"]) if not multi else f"{_subject(r)}:{r['seq']}" for r in records if r["status"] != "ok"] if bad: lines.append(f"- 除外された seq: {', '.join(bad)}") - lines += [f"- 解析設定: signal={opt.signal}, flip_sign={opt.flip_sign}, settle={opt.settle:g}s, " + if multi: + counts: dict[str, Counter[str]] = multi["counts"] # type: ignore[assignment] + lines.append("- 被験者: " + ", ".join(f"{s} ({c['ok']}試行)" for s, c in sorted(counts.items()))) + flip = f", flip_subjects={','.join(opt.flip_subjects)}" if opt.flip_subjects else "" + lines += [f"- 解析設定: signal={opt.signal}, flip_sign={opt.flip_sign}{flip}, settle={opt.settle:g}s, " f"baseline_trim={opt.baseline_trim:g}s, k={opt.k:g}, pool_max_angle={opt.pool_max_angle:g}°, " - f"basis={opt.basis}", "", "誤差の定義: error = (保持区間の平均yaw − 基準yaw) − 真値角(度)。", ""] - lines += ["## しきい値 (距離別・pitch 基準条件)", "", + f"basis={opt.basis}", "", "誤差の定義: error = (保持区間の平均yaw − 基準yaw) − 真値角(度)。" + + ("閾値・回帰は全被験者のフレームをプールして算出 (SDには被験者間のバイアス差が含まれる)。" if multi else ""), ""] + lines += ["## しきい値 (距離別・pitch 基準条件" + (", 全被験者プール)" if multi else ")"), "", "| 距離[m] | frames | bias | SD | θ_pooled | θ_worst | **θ_th** | 傾き | 切片 | R² | 判定可能な最小角 |", "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"] warn = [] @@ -328,12 +570,17 @@ def build_markdown(records: list[dict[str, object]], conditions: list[dict[str, if thresholds: lines += ["", f"全距離で共通に使う保守的な値: θ_th = {max(t['theta_th_deg'] for t in thresholds.values()):.2f}°"] lines += [f"- {w}" for w in warn] - lines += ["", "## 条件別の結果", "", - "| 距離[m] | pitch | 真値[°] | 試行 | frames | bias | SD(試行間) | SD(試行内) | SD(全体) | RMSE | 95%誤差 | 検出率 | θ超過率 |", - "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"] + if multi: + lines += _multi_markdown(multi) + extra_head = " 被験者 | SD(被験者間) |" if multi else "" + extra_rule = "---:|---:|" if multi else "" + lines += ["", "## 条件別の結果" + (" (全被験者プール)" if multi else ""), "", + "| 距離[m] | pitch | 真値[°] | 試行 |" + extra_head + " frames | bias | SD(試行間) | SD(試行内) | SD(全体) | RMSE | 95%誤差 | 検出率 | θ超過率 |", + "|---:|---:|---:|---:|" + extra_rule + "---:|---:|---:|---:|---:|---:|---:|---:|---:|"] for c in conditions: + extra = f" {c['n_subjects']} | {_fmt(c['sd_subjects_deg'])} |" if multi else "" lines.append(f"| {float(c['distance_m']):g} | {float(c['pitch_deg']):+g} | {float(c['true_yaw_deg']):+g} | " # type: ignore[arg-type] - f"{c['n_trials']} | {c['n_frames']} | {_fmt(c['bias_deg'])} | {_fmt(c['sd_trials_deg'])} | " + f"{c['n_trials']} |{extra} {c['n_frames']} | {_fmt(c['bias_deg'])} | {_fmt(c['sd_trials_deg'])} | " f"{_fmt(c['sd_within_deg'])} | {_fmt(c['sd_total_deg'])} | {_fmt(c['rmse_frame_deg'])} | " f"{_fmt(c['p95_abs_error_deg'])} | {_fmt(c['mean_detection_rate'], 3)} | {_fmt(c['exceed_rate'], 3)} |") lines += ["", "注: θ超過率は保持区間のフレームのうち |Δyaw| ≥ θ_th の割合 (真値0°では誤検出率、それ以外では検出率)。", @@ -341,6 +588,49 @@ def build_markdown(records: list[dict[str, object]], conditions: list[dict[str, return "\n".join(lines) + "\n" +def _multi_markdown(multi: dict[str, object]) -> list[str]: + lines: list[str] = ["", "## 被験者別しきい値 (各被験者単独で算出)", "", + "| 被験者 | 距離[m] | frames | bias | SD | θ_th | 傾き | R² | 判定可能な最小角 |", + "|---|---:|---:|---:|---:|---:|---:|---:|---:|"] + per_subject: dict[str, dict[str, object]] = multi["per_subject"] # type: ignore[assignment] + by_distance: dict[float, list[float]] = {} + warn: list[str] = [] + for subject, info in per_subject.items(): + for distance, t in sorted(info["thresholds"].items()): # type: ignore[attr-defined] + lines.append(f"| {subject} | {distance:g} | {int(t['n_frames'])} | {_fmt(t['bias_deg'])} | {_fmt(t['sd_deg'])} | " + f"**{_fmt(t['theta_th_deg'])}** | {_fmt(t['slope'], 3)} | {_fmt(t['r2'], 3)} | " + f"{_fmt(t['min_detectable_deg'])} |") + by_distance.setdefault(distance, []).append(t["theta_th_deg"]) + if not math.isnan(t["slope"]) and t["slope"] < 0: + warn.append(f"{subject} / {distance:g} m: 傾きが負です。符号定義が逆の可能性があります (--flip-subject {subject})。") + if by_distance: + lines += ["", "距離別の被験者間レンジ(θ_th): " + ", ".join( + f"{d:g} m: {min(v):.1f}〜{max(v):.1f}°" for d, v in sorted(by_distance.items())), + f"全被験者・全距離で最も保守的な値: θ_th = {max(max(v) for v in by_distance.values()):.2f}°"] + lines += [f"- {w}" for w in warn] + loso: list[dict[str, object]] = multi["loso"] # type: ignore[assignment] + if loso: + lines += ["", "## 被験者除外 (leave-one-subject-out) 検証", "", + "他の被験者だけで求めたθ_thを除外した被験者に適用した場合の、真値0°の誤検出率と判定可能な最小角。", "", + "| 除外した被験者 | 距離[m] | 学習θ_th | 誤検出率(0°) | 判定可能な最小角 |", "|---|---:|---:|---:|---:|"] + lines += [f"| {r['held_out']} | {float(r['distance_m']):g} | {_fmt(r['theta_train_deg'])} | " # type: ignore[arg-type] + f"{_fmt(r['false_alarm_rate'], 3)} | {_fmt(r['min_detectable_deg'])} |" for r in loso] + asym: list[dict[str, object]] = multi["asymmetry"] # type: ignore[assignment] + if asym: + lines += ["", "## 正負方向の非対称 (ゲイン = 設定方向への推定変化量 / 設定角, 1が理想)", "", + "| 被験者 | セル数 | ゲイン(正) | ゲイン(負) | 差(正−負) | p(符号反転検定) |", "|---|---:|---:|---:|---:|---:|"] + lines += [f"| {r['subject']} | {r['n_cells']} | {_fmt(r['gain_pos'], 3)} | {_fmt(r['gain_neg'], 3)} | " + f"{_fmt(r['gain_diff'], 3)} | {_fmt(r['p_signflip'], 3)} |" for r in asym] + lines += ["", "セル = (被験者×距離×|設定角|)。pはセルを独立とみなした参考値で、被験者数が少ない場合は記述統計として扱う。"] + diffs: list[dict[str, object]] = multi["differences"] # type: ignore[assignment] + if diffs: + lines += ["", "## 被験者間差 (|設定角|≥10°のゲイン, Kruskal-Wallis, 並べ替え検定)", "", + "| 距離[m] | 被験者数 | H | p |", "|---:|---:|---:|---:|"] + lines += [f"| {float(r['distance_m']):g} | {r['n_subjects']} | {_fmt(r['kruskal_h'])} | {_fmt(r['p_permutation'], 3)} |" # type: ignore[arg-type] + for r in diffs] + return lines + + def make_plots(usable: Usable, thresholds: dict[float, dict[str, float]], out_dir: Path, opt: Options) -> list[Path]: try: import matplotlib @@ -413,9 +703,120 @@ def make_plots(usable: Usable, thresholds: dict[float, dict[str, float]], out_di return paths +def make_subject_plots(usable: Usable, thresholds: dict[float, dict[str, float]], cells: Sequence[dict[str, object]], + out_dir: Path, opt: Options) -> list[Path]: + """Per-subject versions of the figures (only meaningful with two or more subjects).""" + try: + import matplotlib + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt + except ImportError: + return [] + ref = [(r, d) for r, d in usable if float(r["pitch_deg"]) == opt.pitch_ref] # type: ignore[arg-type] + subjects = sorted({_subject(r) for r, _ in ref}) + distances = sorted({float(r["distance_m"]) for r, _ in ref}) # type: ignore[arg-type] + if len(subjects) < 2 or not distances: + return [] + out_dir.mkdir(parents=True, exist_ok=True) + axis_label_size, tick_label_size = 14, 13 + cmap = plt.get_cmap("tab10") + color = {s: cmap(i % 10) for i, s in enumerate(subjects)} + paths: list[Path] = [] + n = len(distances) + + fig, axes = plt.subplots(1, n, figsize=(5.2 * n, 5), sharex=True, sharey=True, squeeze=False) + lim = max(abs(float(r["true_yaw_deg"])) for r, _ in ref) * 1.15 + 1 # type: ignore[arg-type] + for ax, distance in zip(axes[0], distances): + ax.plot([-lim, lim], [-lim, lim], "k--", linewidth=1) + theta = thresholds.get(distance, {}).get("theta_th_deg") + if theta is not None: + ax.axhspan(-theta, theta, color="red", alpha=0.08) + for i, subject in enumerate(subjects): + jitter = (i - (len(subjects) - 1) / 2) * 0.8 + pts = [(float(r["true_yaw_deg"]) + jitter, float(r["delta_deg"])) for r, _ in ref # type: ignore[arg-type] + if _subject(r) == subject and float(r["distance_m"]) == distance] # type: ignore[arg-type] + ax.scatter([p[0] for p in pts], [p[1] for p in pts], s=18, alpha=0.8, color=color[subject], label=subject) + ax.set_title(f"{distance:g} m") + ax.set_xlabel("true yaw [deg]", fontsize=axis_label_size) + ax.tick_params(axis="both", labelsize=tick_label_size) + ax.grid(alpha=0.3) + axes[0][0].set_ylabel("estimated yaw change [deg]", fontsize=axis_label_size) + axes[0][-1].legend() + fig.tight_layout() + paths.append(out_dir / "accuracy_scatter_by_subject.png") + fig.savefig(paths[-1]) + plt.close(fig) + + fig, axes = plt.subplots(n, 1, figsize=(10, 3.4 * n), sharex=True, squeeze=False) + angles = sorted({float(r["true_yaw_deg"]) for r, _ in ref}) # type: ignore[arg-type] + step = min((b - a for a, b in zip(angles, angles[1:])), default=5.0) + width = 0.8 * step / len(subjects) + for ax, distance in zip(axes[:, 0], distances): + for i, subject in enumerate(subjects): + offset = (i - (len(subjects) - 1) / 2) * width + data, positions = [], [] + for angle in angles: + values = [float(r["error_deg"]) for r, _ in ref if _subject(r) == subject # type: ignore[arg-type] + and float(r["distance_m"]) == distance and float(r["true_yaw_deg"]) == angle] # type: ignore[arg-type] + if values: + data.append(values) + positions.append(angle + offset) + if data: + box = ax.boxplot(data, positions=positions, widths=width * 0.9, patch_artist=True, showfliers=False, + manage_ticks=False) + for patch in box["boxes"]: + patch.set_facecolor(color[subject]) + patch.set_alpha(0.6) + theta = thresholds.get(distance, {}).get("theta_th_deg") + ax.axhline(0, color="k", linewidth=0.8) + if theta is not None: + ax.axhline(theta, color="r", linestyle=":") + ax.axhline(-theta, color="r", linestyle=":") + ax.set_ylabel(f"error [deg] ({distance:g} m)", fontsize=axis_label_size) + ax.tick_params(axis="both", labelsize=tick_label_size) + ax.grid(alpha=0.3) + axes[-1, 0].set_xticks(angles) + axes[-1, 0].set_xlabel("true yaw [deg]", fontsize=axis_label_size) + handles = [plt.Line2D([0], [0], color=color[s], lw=6, alpha=0.6) for s in subjects] + axes[0, 0].legend(handles, subjects, ncol=len(subjects)) + fig.tight_layout() + paths.append(out_dir / "accuracy_error_by_subject.png") + fig.savefig(paths[-1]) + plt.close(fig) + + if cells: + fig, axes = plt.subplots(1, n, figsize=(5.2 * n, 4.4), sharey=True, squeeze=False) + for ax, distance in zip(axes[0], distances): + for subject in subjects: + rows = sorted((c for c in cells if c["subject"] == subject and c["distance_m"] == distance), + key=lambda c: float(c["abs_angle_deg"])) # type: ignore[arg-type] + ax.plot([c["abs_angle_deg"] for c in rows], [c["gain_pos"] for c in rows], "-o", color=color[subject], + markersize=4, label=f"{subject} (+)") + ax.plot([c["abs_angle_deg"] for c in rows], [c["gain_neg"] for c in rows], "--s", color=color[subject], + markersize=4, label=f"{subject} (-)") + ax.axhline(1.0, color="k", linewidth=0.8) + ax.set_title(f"{distance:g} m") + ax.set_xlabel("|true yaw| [deg]", fontsize=axis_label_size) + ax.tick_params(axis="both", labelsize=tick_label_size) + ax.grid(alpha=0.3) + axes[0][0].set_ylabel("gain (estimated / true)", fontsize=axis_label_size) + axes[0][-1].legend(fontsize=10) + fig.tight_layout() + paths.append(out_dir / "accuracy_gain_by_side.png") + fig.savefig(paths[-1]) + plt.close(fig) + return paths + + def run_report(trial_paths: Sequence[Path], face_paths: Sequence[Path], out_dir: Path, opt: Options, - exclude_seqs: Sequence[int] = (), plots: bool = True) -> dict[str, object]: - trials = load_trials(trial_paths, exclude_seqs) + exclude_seqs: Sequence[int | str] = (), plots: bool = True, + labels: Mapping[Path, str] | None = None, include_subjects: Sequence[str] = (), + exclude_subjects: Sequence[str] = (), aliases: Mapping[str, str] | None = None) -> dict[str, object]: + trials = load_trials(trial_paths, exclude_seqs, labels, aliases) + if include_subjects: + trials = [t for t in trials if t["subject"] in set(include_subjects)] + if exclude_subjects: + trials = [t for t in trials if t["subject"] not in set(exclude_subjects)] series = load_face_series(face_paths) evaluated = [evaluate_trial(t, series, opt) for t in trials] records = [rec for rec, _ in evaluated] @@ -423,25 +824,54 @@ def run_report(trial_paths: Sequence[Path], face_paths: Sequence[Path], out_dir: conditions = summarize_conditions(usable) thresholds = compute_thresholds(usable, opt) add_exceed_rates(conditions, usable, thresholds, opt) + subjects = sorted({_subject(r) for r in records}) + multi: dict[str, object] | None = None + cells: list[dict[str, object]] = [] + if len(split_by_subject(usable)) >= 2: + counts: dict[str, Counter[str]] = {} + for r in records: + counts.setdefault(_subject(r), Counter())[str(r["status"])] += 1 + gains = gain_records(usable, opt) + cells = asymmetry_cells(gains) + multi = {"counts": counts, "per_subject": subject_analysis(usable, opt), + "loso": leave_one_subject_out(usable, opt), "cells": cells, "asymmetry": asymmetry_summary(cells), + "differences": subject_differences(gains)} out_dir.mkdir(parents=True, exist_ok=True) _write_csv(out_dir / "accuracy_trials.csv", TRIAL_COLUMNS, records) _write_csv(out_dir / "accuracy_conditions.csv", CONDITION_COLUMNS, conditions) _write_csv(out_dir / "accuracy_thresholds.csv", THRESHOLD_COLUMNS, [{"distance_m": d, **t} for d, t in sorted(thresholds.items())]) - markdown = build_markdown(records, conditions, thresholds, opt) + if multi: + per: dict[str, dict[str, object]] = multi["per_subject"] # type: ignore[assignment] + _write_csv(out_dir / "accuracy_subject_conditions.csv", SUBJECT_CONDITION_COLUMNS, + [{"subject": s, **c} for s, info in per.items() for c in info["conditions"]]) # type: ignore[attr-defined] + _write_csv(out_dir / "accuracy_subject_thresholds.csv", SUBJECT_THRESHOLD_COLUMNS, + [{"subject": s, "distance_m": d, **t} for s, info in per.items() + for d, t in sorted(info["thresholds"].items())]) # type: ignore[attr-defined] + _write_csv(out_dir / "accuracy_loso.csv", LOSO_COLUMNS, multi["loso"]) # type: ignore[arg-type] + _write_csv(out_dir / "accuracy_asymmetry_cells.csv", ASYMMETRY_CELL_COLUMNS, cells) + _write_csv(out_dir / "accuracy_asymmetry.csv", ASYMMETRY_COLUMNS, multi["asymmetry"]) # type: ignore[arg-type] + _write_csv(out_dir / "accuracy_subject_differences.csv", DIFFERENCE_COLUMNS, multi["differences"]) # type: ignore[arg-type] + markdown = build_markdown(records, conditions, thresholds, opt, multi) (out_dir / "accuracy_report.md").write_text(markdown, encoding="utf-8") - figures = make_plots(usable, thresholds, out_dir / "plots", opt) if plots else [] + figures: list[Path] = [] + if plots: + figures = make_plots(usable, thresholds, out_dir / "plots", opt) + if multi: + figures += make_subject_plots(usable, thresholds, cells, out_dir / "plots", opt) return {"records": records, "conditions": conditions, "thresholds": thresholds, "markdown": markdown, - "figures": figures} + "figures": figures, "subjects": subjects, "multi": multi} def main(argv: Sequence[str] | None = None, out=print) -> int: # type: ignore[no-untyped-def] parser = argparse.ArgumentParser(description="Accuracy report from trial markers and face_frames logs") - parser.add_argument("--trials", nargs="+", required=True, help="trials CSV file(s) from accuracy_session") + parser.add_argument("--trials", nargs="+", required=True, + help="trials CSV file(s) from accuracy_session; prefix LABEL= to name the subject (A=trials_a.csv)") parser.add_argument("--face-csv", nargs="+", required=True, help="*_face_frames.csv file(s); globs allowed") parser.add_argument("--out", type=Path, default=Path("data/analysis/accuracy/report")) parser.add_argument("--signal", choices=("yaw", "raw_yaw"), default="yaw") - parser.add_argument("--flip-sign", action="store_true") + parser.add_argument("--flip-sign", action="store_true", help="flip the yaw sign of every subject") + parser.add_argument("--flip-subject", nargs="*", default=[], help="flip the yaw sign of these subjects only") parser.add_argument("--settle", type=float, default=3.0, help="seconds after the cue ignored (turn in progress)") parser.add_argument("--baseline-trim", type=float, default=0.5) parser.add_argument("--min-frames", type=int, default=10) @@ -452,17 +882,37 @@ def main(argv: Sequence[str] | None = None, out=print) -> int: # type: ignore[n parser.add_argument("--pool-max-angle", type=float, default=15.0) parser.add_argument("--basis", choices=("worst", "pooled"), default="worst") parser.add_argument("--detect-ratio", type=float, default=0.95) - parser.add_argument("--exclude-seq", type=int, nargs="*", default=[]) + parser.add_argument("--exclude-seq", nargs="*", default=[], + help="seq numbers to drop; SUBJECT:seq drops it for that subject only") + parser.add_argument("--subjects", nargs="*", default=[], help="analyse only these subjects") + parser.add_argument("--exclude-subjects", nargs="*", default=[], help="drop these subjects") + parser.add_argument("--alias", nargs="*", default=[], help="merge subject names, e.g. S01b=S01") parser.add_argument("--no-plots", action="store_true") args = parser.parse_args(argv) - trial_paths, face_paths = expand_paths(args.trials), expand_paths(args.face_csv) + labels: dict[Path, str] = {} + trial_paths: list[Path] = [] + for token in args.trials: + label, pattern = split_label(token) + found = expand_paths([pattern]) + trial_paths.extend(found) + if label: + labels.update({path: label for path in found}) + face_paths = expand_paths(args.face_csv) missing = [str(p) for p in trial_paths + face_paths if not p.is_file()] if missing or not trial_paths or not face_paths: print(f"file not found: {', '.join(missing) or 'no input files'}", file=sys.stderr) return 2 + try: + aliases = dict(item.split("=", 1) for item in args.alias) + _parse_exclusions(args.exclude_seq) + except ValueError: + print("invalid --alias (use OLD=NEW) or --exclude-seq (use N or SUBJECT:N)", file=sys.stderr) + return 2 opt = Options(args.signal, args.flip_sign, args.settle, args.baseline_trim, args.min_frames, args.face_every, - args.track_id, args.fps, args.k, args.pool_max_angle, args.basis, args.detect_ratio) - result = run_report(trial_paths, face_paths, args.out, opt, args.exclude_seq, not args.no_plots) + args.track_id, args.fps, args.k, args.pool_max_angle, args.basis, args.detect_ratio, + flip_subjects=tuple(args.flip_subject)) + result = run_report(trial_paths, face_paths, args.out, opt, args.exclude_seq, not args.no_plots, labels, + args.subjects, args.exclude_subjects, aliases) if not result["records"]: print("no usable trials (status ok) in the trials file", file=sys.stderr) return 1 diff --git a/src/accuracy_report.py.orig b/src/accuracy_report.py.orig new file mode 100644 index 0000000..ec5b323 --- /dev/null +++ b/src/accuracy_report.py.orig @@ -0,0 +1,475 @@ +"""Accuracy-experiment report: join trial markers with ``*_face_frames.csv`` by timestamp. + +Usage (after ``src.accuracy_session`` recorded the markers):: + + python -m src.accuracy_report --trials data/analysis/accuracy/trials_.csv ^ + --face-csv "data/logs/face/*_face_frames.csv" + +Per trial the baseline yaw (subject facing the camera) is averaged, the hold-phase yaw is +compared against it, and ``error = (hold - baseline) - true_yaw``. Outputs (``--out``): +``accuracy_trials.csv``, ``accuracy_conditions.csv``, ``accuracy_thresholds.csv``, +``accuracy_report.md`` and ``plots/*.png``. + +Threshold: ``theta_th = |bias| + k * SD`` of the frame errors for conditions with +``|true_yaw| <= --pool-max-angle`` (pitch 0), per distance. ``--basis worst`` (default) takes the +worst single condition, ``pooled`` pools all of them. k=3 and the pooling range are provisional +choices to be agreed on with the supervisor. +""" + +from __future__ import annotations + +import argparse +import bisect +import csv +import glob +import math +import statistics +import sys +from collections import Counter +from dataclasses import dataclass +from datetime import datetime +from pathlib import Path +from typing import Sequence + +NAN = float("nan") +TRIAL_COLUMNS = ["session_id", "seq", "condition_id", "distance_m", "true_yaw_deg", "pitch_deg", "repeat", "subject", + "status", "track_id", "baseline_n", "baseline_mean", "baseline_sd", "hold_n", "hold_mean", "hold_sd", + "delta_deg", "error_deg", "detection_rate"] +CONDITION_COLUMNS = ["distance_m", "pitch_deg", "true_yaw_deg", "n_trials", "n_frames", "bias_deg", "sd_trials_deg", + "sd_within_deg", "sd_total_deg", "rmse_frame_deg", "p95_abs_error_deg", "mean_detection_rate", + "exceed_rate"] +THRESHOLD_COLUMNS = ["distance_m", "n_frames", "bias_deg", "sd_deg", "theta_pooled_deg", "theta_worst_deg", + "theta_th_deg", "slope", "intercept_deg", "r2", "min_detectable_deg"] + + +@dataclass(frozen=True) +class Options: + signal: str = "yaw" + flip_sign: bool = False + settle: float = 1.5 + baseline_trim: float = 0.5 + min_frames: int = 10 + face_every: int = 1 + track_id: int | None = None + fps: float | None = None + k: float = 3.0 + pool_max_angle: float = 15.0 + basis: str = "worst" + detect_ratio: float = 0.95 + pitch_ref: float = 0.0 + + +@dataclass +class FaceSeries: + times: list[float] + track_ids: list[int] + pose: list[bool] + yaw: list[float | None] + raw_yaw: list[float | None] + fps: list[float] + + +def _ts(value: object) -> float | None: + try: + return datetime.fromisoformat(str(value).strip()).timestamp() + except ValueError: + return None + + +def _float(value: object) -> float | None: + try: + number = float(str(value)) + except ValueError: + return None + return number if math.isfinite(number) else None + + +def _mean(values: Sequence[float]) -> float: + return statistics.fmean(values) if values else NAN + + +def _sd(values: Sequence[float]) -> float: + return statistics.stdev(values) if len(values) >= 2 else NAN + + +def _wrap(angle: float) -> float: + return (angle + 180.0) % 360.0 - 180.0 + + +def percentile(sorted_values: Sequence[float], q: float) -> float: + if not sorted_values: + return NAN + position = (len(sorted_values) - 1) * q / 100.0 + lower, upper = math.floor(position), math.ceil(position) + return float(sorted_values[lower] + (sorted_values[upper] - sorted_values[lower]) * (position - lower)) + + +def expand_paths(patterns: Sequence[str]) -> list[Path]: + """Expand globs ourselves because Windows shells do not.""" + paths: list[Path] = [] + for pattern in patterns: + matches = sorted(glob.glob(pattern)) if any(c in pattern for c in "*?[") else [pattern] + paths.extend(Path(m) for m in matches) + return paths + + +def load_face_series(paths: Sequence[Path], default_fps: float = 30.0) -> FaceSeries: + rows: list[tuple[float, int, bool, float | None, float | None, float]] = [] + for path in paths: + file_rows: list[tuple[float, int, bool, float | None, float | None]] = [] + frames: list[tuple[float, float]] = [] + with path.open(encoding="utf-8-sig", newline="") as stream: + for row in csv.DictReader(stream): + t, track = _ts(row.get("timestamp")), _float(row.get("track_id")) + if t is None or track is None: + continue + file_rows.append((t, int(track), str(row.get("pose_estimated", "")).lower() == "true", + _float(row.get("yaw")), _float(row.get("raw_yaw")))) + frame = _float(row.get("frame_number")) + if frame is not None: + frames.append((t, frame)) + fps = default_fps + if len(frames) >= 2: + dt, df = frames[-1][0] - frames[0][0], frames[-1][1] - frames[0][1] + if dt > 0 and df > 0: + fps = df / dt + rows.extend((*r, fps) for r in file_rows) + rows.sort(key=lambda r: r[0]) + return FaceSeries([r[0] for r in rows], [r[1] for r in rows], [r[2] for r in rows], [r[3] for r in rows], + [r[4] for r in rows], [r[5] for r in rows]) + + +def load_trials(paths: Sequence[Path], exclude_seqs: Sequence[int] = ()) -> list[dict[str, str]]: + """Last log row per (session, seq); only status ``ok`` survives.""" + last: dict[tuple[str, int], dict[str, str]] = {} + for path in paths: + with path.open(encoding="utf-8-sig", newline="") as stream: + for row in csv.DictReader(stream): + last[(row.get("session_id", ""), int(row["seq"]))] = row + return [row for (_, seq), row in sorted(last.items()) if row["status"] == "ok" and seq not in set(exclude_seqs)] + + +def evaluate_trial(trial: dict[str, str], series: FaceSeries, opt: Options) -> tuple[dict[str, object], list[float]]: + """Return the trial record and the frame-level deltas (degrees, signed, hold phase).""" + record: dict[str, object] = {column: NAN for column in TRIAL_COLUMNS} + record.update({ + "session_id": trial.get("session_id", ""), "seq": int(trial["seq"]), "condition_id": trial["condition_id"], + "distance_m": float(trial["distance_m"]), "true_yaw_deg": float(trial["true_yaw_deg"]), + "pitch_deg": float(trial["pitch_deg"]), "repeat": trial.get("repeat", ""), "subject": trial.get("subject", ""), + "status": "ok", "track_id": "", + }) + t_base, t_cue, t_end = (_ts(trial.get(key)) for key in ("t_baseline_start", "t_turn_cue", "t_hold_end")) + if t_base is None or t_cue is None or t_end is None: + record["status"] = "bad_timestamps" + return record, [] + b0, b1, h0, h1 = t_base + opt.baseline_trim, t_cue, t_cue + opt.settle, t_end + values = series.raw_yaw if opt.signal == "raw_yaw" else series.yaw + + def span(a: float, b: float) -> tuple[int, int]: + return bisect.bisect_left(series.times, a), bisect.bisect_right(series.times, b) + + lo, hi = span(b0, h1) + counts = Counter(series.track_ids[i] for i in range(lo, hi) if series.pose[i] and values[i] is not None) + track = opt.track_id if opt.track_id is not None else (counts.most_common(1)[0][0] if counts else None) + if track is None or counts.get(track, 0) == 0: + record["status"] = "no_track" + return record, [] + record["track_id"] = track + + def pick(a: float, b: float) -> list[int]: + s, e = span(a, b) + return [i for i in range(s, e) if series.track_ids[i] == track and series.pose[i] and values[i] is not None] + + base, hold = pick(b0, b1), pick(h0, h1) + base_yaw = [float(values[i]) for i in base] # type: ignore[arg-type] + hold_yaw = [float(values[i]) for i in hold] # type: ignore[arg-type] + record.update({"baseline_n": len(base), "hold_n": len(hold), "baseline_mean": _mean(base_yaw), + "baseline_sd": _sd(base_yaw), "hold_mean": _mean(hold_yaw), "hold_sd": _sd(hold_yaw)}) + if len(base) < opt.min_frames: + record["status"] = "baseline_insufficient" + return record, [] + if len(hold) < opt.min_frames: + record["status"] = "hold_insufficient" + return record, [] + sign = -1.0 if opt.flip_sign else 1.0 + base_mean = float(record["baseline_mean"]) # type: ignore[arg-type] + deltas = [_wrap(sign * (y - base_mean)) for y in hold_yaw] + fps = opt.fps or series.fps[hold[0]] + expected = max(1.0, (h1 - h0) * fps / max(1, opt.face_every)) + record["detection_rate"] = min(1.0, len(hold) / expected) + record["delta_deg"] = _mean(deltas) + record["error_deg"] = float(record["delta_deg"]) - float(record["true_yaw_deg"]) # type: ignore[arg-type] + return record, deltas + + +Usable = list[tuple[dict[str, object], list[float]]] + + +def _errors(record: dict[str, object], deltas: Sequence[float]) -> list[float]: + return [d - float(record["true_yaw_deg"]) for d in deltas] # type: ignore[arg-type] + + +def summarize_conditions(usable: Usable) -> list[dict[str, object]]: + groups: dict[tuple[float, float, float], Usable] = {} + for record, deltas in usable: + key = (float(record["distance_m"]), float(record["pitch_deg"]), float(record["true_yaw_deg"])) # type: ignore[arg-type] + groups.setdefault(key, []).append((record, deltas)) + rows: list[dict[str, object]] = [] + for (distance, pitch, true), items in sorted(groups.items()): + frame_errors = [e for record, deltas in items for e in _errors(record, deltas)] + within = [float(r["hold_sd"]) ** 2 for r, _ in items if not math.isnan(float(r["hold_sd"]))] # type: ignore[arg-type] + rows.append({ + "distance_m": distance, "pitch_deg": pitch, "true_yaw_deg": true, "n_trials": len(items), + "n_frames": len(frame_errors), "bias_deg": _mean([float(r["error_deg"]) for r, _ in items]), # type: ignore[arg-type] + "sd_trials_deg": _sd([float(r["error_deg"]) for r, _ in items]), # type: ignore[arg-type] + "sd_within_deg": math.sqrt(_mean(within)) if within else NAN, "sd_total_deg": _sd(frame_errors), + "rmse_frame_deg": math.sqrt(_mean([e * e for e in frame_errors])), + "p95_abs_error_deg": percentile(sorted(abs(e) for e in frame_errors), 95), + "mean_detection_rate": _mean([float(r["detection_rate"]) for r, _ in items]), # type: ignore[arg-type] + "exceed_rate": NAN, + }) + return rows + + +def _line_fit(points: Sequence[tuple[float, float]]) -> tuple[float, float, float]: + xs, ys = [p[0] for p in points], [p[1] for p in points] + if len(points) < 2 or len(set(xs)) < 2: + return NAN, NAN, NAN + mx, my = statistics.fmean(xs), statistics.fmean(ys) + sxx = sum((x - mx) ** 2 for x in xs) + sxy = sum((x - mx) * (y - my) for x, y in points) + syy = sum((y - my) ** 2 for y in ys) + slope = sxy / sxx + return slope, my - slope * mx, (sxy * sxy / (sxx * syy) if syy > 0 else NAN) + + +def compute_thresholds(usable: Usable, opt: Options) -> dict[float, dict[str, float]]: + result: dict[float, dict[str, float]] = {} + for distance in sorted({float(r["distance_m"]) for r, _ in usable}): # type: ignore[arg-type] + ref = [(r, d) for r, d in usable if float(r["distance_m"]) == distance # type: ignore[arg-type] + and float(r["pitch_deg"]) == opt.pitch_ref] # type: ignore[arg-type] + pooled_sel = [(r, d) for r, d in ref if abs(float(r["true_yaw_deg"])) <= opt.pool_max_angle] # type: ignore[arg-type] + all_errors = [e for r, d in pooled_sel for e in _errors(r, d)] + if len(all_errors) < 2: + continue + pooled = abs(_mean(all_errors)) + opt.k * _sd(all_errors) + by_condition: dict[float, list[float]] = {} + for r, d in pooled_sel: + by_condition.setdefault(float(r["true_yaw_deg"]), []).extend(_errors(r, d)) # type: ignore[arg-type] + worst_values = [abs(_mean(e)) + opt.k * _sd(e) for e in by_condition.values() if len(e) >= 2] + worst = max(worst_values) if worst_values else pooled + slope, intercept, r2 = _line_fit([(float(r["true_yaw_deg"]), float(r["delta_deg"])) for r, _ in ref]) # type: ignore[arg-type] + result[distance] = { + "n_frames": float(len(all_errors)), "bias_deg": _mean(all_errors), "sd_deg": _sd(all_errors), + "theta_pooled_deg": pooled, "theta_worst_deg": worst, + "theta_th_deg": worst if opt.basis == "worst" else pooled, + "slope": slope, "intercept_deg": intercept, "r2": r2, "min_detectable_deg": NAN, + } + return result + + +def add_exceed_rates(conditions: list[dict[str, object]], usable: Usable, thresholds: dict[float, dict[str, float]], + opt: Options) -> None: + """exceed_rate = share of hold frames with |delta| >= theta_th (false-positive rate at 0 deg).""" + frames: dict[tuple[float, float, float], list[float]] = {} + for record, deltas in usable: + key = (float(record["distance_m"]), float(record["pitch_deg"]), float(record["true_yaw_deg"])) # type: ignore[arg-type] + frames.setdefault(key, []).extend(deltas) + for row in conditions: + theta = thresholds.get(float(row["distance_m"]), {}).get("theta_th_deg") # type: ignore[arg-type] + data = frames[(float(row["distance_m"]), float(row["pitch_deg"]), float(row["true_yaw_deg"]))] # type: ignore[arg-type] + row["exceed_rate"] = sum(abs(d) >= theta for d in data) / len(data) if theta is not None and data else NAN + for distance, info in thresholds.items(): + ref = [c for c in conditions if c["distance_m"] == distance and c["pitch_deg"] == opt.pitch_ref + and float(c["true_yaw_deg"]) != 0.0 and not math.isnan(float(c["exceed_rate"]))] # type: ignore[arg-type] + for angle in sorted({abs(float(c["true_yaw_deg"])) for c in ref}): # type: ignore[arg-type] + if all(float(c["exceed_rate"]) >= opt.detect_ratio for c in ref if abs(float(c["true_yaw_deg"])) >= angle): # type: ignore[arg-type] + info["min_detectable_deg"] = angle + break + + +def _fmt(value: object, digits: int = 2) -> str: + if isinstance(value, float): + return "-" if math.isnan(value) else f"{value:.{digits}f}" + return str(value) + + +def _write_csv(path: Path, columns: Sequence[str], rows: Sequence[dict[str, object]]) -> None: + with path.open("w", encoding="utf-8-sig", newline="") as stream: + writer = csv.DictWriter(stream, fieldnames=list(columns), extrasaction="ignore") + writer.writeheader() + for row in rows: + writer.writerow({k: (_fmt(v, 4) if isinstance(v, float) else v) for k, v in row.items()}) + + +def build_markdown(records: list[dict[str, object]], conditions: list[dict[str, object]], + thresholds: dict[float, dict[str, float]], opt: Options) -> str: + lines = ["# 振り向き角度推定の精度実験 結果", ""] + status = Counter(str(r["status"]) for r in records) + lines += [f"- 記録された有効試行: {len(records)}、解析に使用: {status['ok']}", "- 除外: " + ( + ", ".join(f"{k}={v}" for k, v in status.items() if k != "ok") or "なし")] + bad = [str(r["seq"]) for r in records if r["status"] != "ok"] + if bad: + lines.append(f"- 除外された seq: {', '.join(bad)}") + lines += [f"- 解析設定: signal={opt.signal}, flip_sign={opt.flip_sign}, settle={opt.settle:g}s, " + f"baseline_trim={opt.baseline_trim:g}s, k={opt.k:g}, pool_max_angle={opt.pool_max_angle:g}°, " + f"basis={opt.basis}", "", "誤差の定義: error = (保持区間の平均yaw − 基準yaw) − 真値角(度)。", ""] + lines += ["## しきい値 (距離別・pitch 基準条件)", "", + "| 距離[m] | frames | bias | SD | θ_pooled | θ_worst | **θ_th** | 傾き | 切片 | R² | 判定可能な最小角 |", + "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"] + warn = [] + for distance, t in sorted(thresholds.items()): + lines.append(f"| {distance:g} | {int(t['n_frames'])} | {_fmt(t['bias_deg'])} | {_fmt(t['sd_deg'])} | " + f"{_fmt(t['theta_pooled_deg'])} | {_fmt(t['theta_worst_deg'])} | **{_fmt(t['theta_th_deg'])}** | " + f"{_fmt(t['slope'], 3)} | {_fmt(t['intercept_deg'])} | {_fmt(t['r2'], 3)} | " + f"{_fmt(t['min_detectable_deg'])} |") + if not math.isnan(t["slope"]) and t["slope"] < 0: + warn.append(f"{distance:g} m: 傾きが負です。角度の符号定義が逆の可能性があります (--flip-sign)。") + if thresholds: + lines += ["", f"全距離で共通に使う保守的な値: θ_th = {max(t['theta_th_deg'] for t in thresholds.values()):.2f}°"] + lines += [f"- {w}" for w in warn] + lines += ["", "## 条件別の結果", "", + "| 距離[m] | pitch | 真値[°] | 試行 | frames | bias | SD(試行間) | SD(試行内) | SD(全体) | RMSE | 95%誤差 | 検出率 | θ超過率 |", + "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"] + for c in conditions: + lines.append(f"| {float(c['distance_m']):g} | {float(c['pitch_deg']):+g} | {float(c['true_yaw_deg']):+g} | " # type: ignore[arg-type] + f"{c['n_trials']} | {c['n_frames']} | {_fmt(c['bias_deg'])} | {_fmt(c['sd_trials_deg'])} | " + f"{_fmt(c['sd_within_deg'])} | {_fmt(c['sd_total_deg'])} | {_fmt(c['rmse_frame_deg'])} | " + f"{_fmt(c['p95_abs_error_deg'])} | {_fmt(c['mean_detection_rate'], 3)} | {_fmt(c['exceed_rate'], 3)} |") + lines += ["", "注: θ超過率は保持区間のフレームのうち |Δyaw| ≥ θ_th の割合 (真値0°では誤検出率、それ以外では検出率)。", + "検出率は姿勢推定に成功したフレーム数 / 期待フレーム数 (フレームレートはログから推定)。"] + return "\n".join(lines) + "\n" + + +def make_plots(usable: Usable, thresholds: dict[float, dict[str, float]], out_dir: Path, opt: Options) -> list[Path]: + try: + import matplotlib + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt + except ImportError: + return [] + out_dir.mkdir(parents=True, exist_ok=True) + axis_label_size = 14 + tick_label_size = 13 + ref = [(r, d) for r, d in usable if float(r["pitch_deg"]) == opt.pitch_ref] # type: ignore[arg-type] + distances = sorted({float(r["distance_m"]) for r, _ in ref}) # type: ignore[arg-type] + if not distances: + return [] + paths: list[Path] = [] + fig, ax = plt.subplots(figsize=(6, 6)) + for distance in distances: + pts = [(float(r["true_yaw_deg"]), float(r["delta_deg"])) for r, _ in ref if float(r["distance_m"]) == distance] # type: ignore[arg-type] + ax.scatter([p[0] for p in pts], [p[1] for p in pts], s=18, alpha=0.7, label=f"{distance:g} m") + lim = max(abs(float(r["true_yaw_deg"])) for r, _ in ref) * 1.15 + 1 # type: ignore[arg-type] + ax.plot([-lim, lim], [-lim, lim], "k--", linewidth=1, label="y = x") + ax.set(title="Estimated vs true") + ax.set_xlabel("true yaw [deg]", fontsize=axis_label_size) + ax.set_ylabel("estimated yaw change [deg]", fontsize=axis_label_size) + ax.tick_params(axis="both", labelsize=tick_label_size) + ax.legend() + ax.grid(alpha=0.3) + fig.tight_layout() + paths.append(out_dir / "accuracy_scatter.png") + fig.savefig(paths[-1]) + plt.close(fig) + fig, axes = plt.subplots(len(distances), 1, figsize=(8, 3.2 * len(distances)), squeeze=False) + for ax, distance in zip(axes[:, 0], distances): + by_angle: dict[float, list[float]] = {} + for r, d in ref: + if float(r["distance_m"]) == distance: # type: ignore[arg-type] + by_angle.setdefault(float(r["true_yaw_deg"]), []).extend(_errors(r, d)) # type: ignore[arg-type] + angles = sorted(by_angle) + ax.boxplot([by_angle[a] for a in angles], positions=list(range(len(angles))), showfliers=False) + ax.set_xticks(list(range(len(angles)))) + ax.set_xticklabels([f"{a:+g}" for a in angles]) + theta = thresholds.get(distance, {}).get("theta_th_deg") + ax.axhline(0, color="k", linewidth=0.8) + if theta is not None: + ax.axhline(theta, color="r", linestyle=":") + ax.axhline(-theta, color="r", linestyle=":") + ax.set(title=f"Frame error by true yaw ({distance:g} m)") + ax.set_xlabel("true yaw [deg]", fontsize=axis_label_size) + ax.set_ylabel("error [deg]", fontsize=axis_label_size) + ax.tick_params(axis="both", labelsize=tick_label_size) + fig.tight_layout() + paths.append(out_dir / "accuracy_error_by_condition.png") + fig.savefig(paths[-1]) + plt.close(fig) + fig, ax = plt.subplots(figsize=(7, 4)) + for distance in distances: + errors = [e for r, d in ref if float(r["distance_m"]) == distance # type: ignore[arg-type] + and abs(float(r["true_yaw_deg"])) <= opt.pool_max_angle for e in _errors(r, d)] # type: ignore[arg-type] + if errors: + ax.hist(errors, bins=40, histtype="step", density=True, label=f"{distance:g} m") + ax.set(title=f"Frame error distribution (|true| <= {opt.pool_max_angle:g} deg)") + ax.set_xlabel("error [deg]", fontsize=axis_label_size) + ax.set_ylabel("density", fontsize=axis_label_size) + ax.tick_params(axis="both", labelsize=tick_label_size) + ax.legend() + fig.tight_layout() + paths.append(out_dir / "accuracy_error_hist.png") + fig.savefig(paths[-1]) + plt.close(fig) + return paths + + +def run_report(trial_paths: Sequence[Path], face_paths: Sequence[Path], out_dir: Path, opt: Options, + exclude_seqs: Sequence[int] = (), plots: bool = True) -> dict[str, object]: + trials = load_trials(trial_paths, exclude_seqs) + series = load_face_series(face_paths) + evaluated = [evaluate_trial(t, series, opt) for t in trials] + records = [rec for rec, _ in evaluated] + usable: Usable = [(rec, deltas) for rec, deltas in evaluated if rec["status"] == "ok"] + conditions = summarize_conditions(usable) + thresholds = compute_thresholds(usable, opt) + add_exceed_rates(conditions, usable, thresholds, opt) + out_dir.mkdir(parents=True, exist_ok=True) + _write_csv(out_dir / "accuracy_trials.csv", TRIAL_COLUMNS, records) + _write_csv(out_dir / "accuracy_conditions.csv", CONDITION_COLUMNS, conditions) + _write_csv(out_dir / "accuracy_thresholds.csv", THRESHOLD_COLUMNS, + [{"distance_m": d, **t} for d, t in sorted(thresholds.items())]) + markdown = build_markdown(records, conditions, thresholds, opt) + (out_dir / "accuracy_report.md").write_text(markdown, encoding="utf-8") + figures = make_plots(usable, thresholds, out_dir / "plots", opt) if plots else [] + return {"records": records, "conditions": conditions, "thresholds": thresholds, "markdown": markdown, + "figures": figures} + + +def main(argv: Sequence[str] | None = None, out=print) -> int: # type: ignore[no-untyped-def] + parser = argparse.ArgumentParser(description="Accuracy report from trial markers and face_frames logs") + parser.add_argument("--trials", nargs="+", required=True, help="trials CSV file(s) from accuracy_session") + parser.add_argument("--face-csv", nargs="+", required=True, help="*_face_frames.csv file(s); globs allowed") + parser.add_argument("--out", type=Path, default=Path("data/analysis/accuracy/report")) + parser.add_argument("--signal", choices=("yaw", "raw_yaw"), default="yaw") + parser.add_argument("--flip-sign", action="store_true") + parser.add_argument("--settle", type=float, default=3.0, help="seconds after the cue ignored (turn in progress)") + parser.add_argument("--baseline-trim", type=float, default=0.5) + parser.add_argument("--min-frames", type=int, default=10) + parser.add_argument("--face-every", type=int, default=1, help="face.process_every_n_frames of the run") + parser.add_argument("--track-id", type=int, default=None) + parser.add_argument("--fps", type=float, default=None, help="override the frame rate used for detection rate") + parser.add_argument("--k", type=float, default=3.0) + parser.add_argument("--pool-max-angle", type=float, default=15.0) + parser.add_argument("--basis", choices=("worst", "pooled"), default="worst") + parser.add_argument("--detect-ratio", type=float, default=0.95) + parser.add_argument("--exclude-seq", type=int, nargs="*", default=[]) + parser.add_argument("--no-plots", action="store_true") + args = parser.parse_args(argv) + trial_paths, face_paths = expand_paths(args.trials), expand_paths(args.face_csv) + missing = [str(p) for p in trial_paths + face_paths if not p.is_file()] + if missing or not trial_paths or not face_paths: + print(f"file not found: {', '.join(missing) or 'no input files'}", file=sys.stderr) + return 2 + opt = Options(args.signal, args.flip_sign, args.settle, args.baseline_trim, args.min_frames, args.face_every, + args.track_id, args.fps, args.k, args.pool_max_angle, args.basis, args.detect_ratio) + result = run_report(trial_paths, face_paths, args.out, opt, args.exclude_seq, not args.no_plots) + if not result["records"]: + print("no usable trials (status ok) in the trials file", file=sys.stderr) + return 1 + out(str(result["markdown"])) + out(f"-> {args.out}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/accuracy_session.py b/src/accuracy_session.py index 17ac898..48ee873 100644 --- a/src/accuracy_session.py +++ b/src/accuracy_session.py @@ -44,6 +44,11 @@ INTERVAL_AUDIO = { "block": Path("assets/audio/block.wav"), "auto": Path("assets/audio/auto_interval.wav"), } +TRIAL_AUDIO = { + "trial_start": Path("assets/audio/1.wav"), + "turn": Path("assets/audio/2.wav"), + "trial_end": Path("assets/audio/3.wav"), +} # Angle announcement is composed from separate sign and number WAV files. # Prepare these files under assets/audio/ using exactly these names. ANGLE_SIGN_AUDIO = { @@ -184,6 +189,14 @@ class IntervalAudioPlayer: self._sounds[kind] = mixer.Sound(str(path)) except Exception as exc: LOGGER.error("Failed to load interval audio %s: %s", path, exc) + for kind, path in TRIAL_AUDIO.items(): + if not path.is_file(): + LOGGER.warning("Trial audio file does not exist: %s", path) + continue + try: + self._sounds[kind] = mixer.Sound(str(path)) + except Exception as exc: + LOGGER.error("Failed to load trial audio %s: %s", path, exc) for sign, path in ANGLE_SIGN_AUDIO.items(): if not path.is_file(): LOGGER.warning("Angle sign audio file does not exist (%s): %s", sign, path) @@ -347,7 +360,7 @@ class SessionRunner: continue if first_trial_move_pending: self.interval_audio.play("block") - self._countdown(20.0, label="initial move") + self._countdown(10.0, label="initial move") first_trial_move_pending = False elif block_change and self.block_interval > 0: self.interval_audio.play("block") @@ -374,17 +387,17 @@ class SessionRunner: def _execute(self, row: PlanRow) -> tuple[str, str, str]: t_base = self.clock() - self.beep("start") + self.interval_audio.play("trial_start") self.out(f" baseline {self.baseline_sec:g} s: face the camera and stay still") self.sleep(self.baseline_sec) t_cue = self.clock() - self.beep("turn") + self.interval_audio.play("turn") self.out(f" TURN to {_sign(row.true_yaw_deg)} deg (head only) — you have {self.turn_grace_sec:g} s to turn") self.sleep(self.turn_grace_sec) self.out(f" HOLD {self.hold_sec:g} s") self.sleep(self.hold_sec) t_end = self.clock() - self.beep("end") + #self.interval_audio.play("trial_end") self.out(" return to the front") return t_base, t_cue, t_end diff --git a/tests/test_accuracy_report.py b/tests/test_accuracy_report.py index fc5887c..fcbe75e 100644 --- a/tests/test_accuracy_report.py +++ b/tests/test_accuracy_report.py @@ -7,7 +7,8 @@ from pathlib import Path import pytest from src.accuracy_report import ( - Options, evaluate_trial, load_face_series, load_trials, main, percentile, run_report, + Options, evaluate_trial, kruskal_permutation, load_face_series, load_trials, main, percentile, run_report, + split_label, _signflip_p, ) from src.accuracy_session import SessionRunner, build_plan @@ -16,8 +17,8 @@ FPS = 30.0 class FakeTime: - def __init__(self) -> None: - self.now = datetime(2026, 10, 4, 12, 0, 0, tzinfo=TZ) + def __init__(self, hour: int = 12) -> None: + self.now = datetime(2026, 10, 4, hour, 0, 0, tzinfo=TZ) def sleep(self, seconds: float) -> None: self.now += timedelta(seconds=seconds) @@ -26,10 +27,15 @@ class FakeTime: return self.now.isoformat(timespec="milliseconds") -def _simulate(tmp_path: Path, scale: float = 0.95, offset: float = -19.0, noise: float = 0.5, face_dropout: float = 0.0): - """Run a fake session and write a matching face_frames CSV (30 fps, yaw ramps after each cue).""" +def _simulate(tmp_path: Path, scale: float = 0.95, offset: float = -19.0, noise: float = 0.5, face_dropout: float = 0.0, + hour: int = 12, seed: int = 0): + """Run a fake session and write a matching face_frames CSV (30 fps, yaw ramps after each cue). + + ``hour`` shifts the wall clock so that several simulated subjects never overlap in time. + """ + tmp_path.mkdir(parents=True, exist_ok=True) plan = build_plan([1, 3], [0, 5, -5, 10, -10], 3, seed=3) - fake, trials = FakeTime(), tmp_path / "trials.csv" + fake, trials = FakeTime(hour), tmp_path / "trials.csv" start = fake.now def operator_pause(_prompt: str) -> str: @@ -39,10 +45,10 @@ def _simulate(tmp_path: Path, scale: float = 0.95, offset: float = -19.0, noise: SessionRunner(plan, trials, "S", input_fn=operator_pause, out=lambda _m: None, sleep=fake.sleep, clock=fake.clock, beep=lambda _k: None).run() end = fake.now + timedelta(seconds=2) - rng = random.Random(0) + rng = random.Random(seed) marks = [(datetime.fromisoformat(r["t_turn_cue"]), float(r["true_yaw_deg"]), datetime.fromisoformat(r["t_hold_end"])) for r in csv.DictReader(trials.open(encoding="utf-8-sig"))] - path = tmp_path / "20261004_120000_exp_face_frames.csv" + path = tmp_path / f"20261004_{hour:02d}0000_exp_face_frames.csv" with path.open("w", encoding="utf-8-sig", newline="") as stream: writer = csv.writer(stream) writer.writerow(["frame_number", "timestamp", "track_id", "pose_estimated", "yaw", "raw_yaw"]) @@ -83,6 +89,9 @@ def test_full_pipeline_recovers_scale_bias_and_threshold(tmp_path: Path) -> None for name in ("accuracy_trials.csv", "accuracy_conditions.csv", "accuracy_thresholds.csv", "accuracy_report.md"): assert (tmp_path / "out" / name).is_file() assert "θ_th" in result["markdown"] + # one subject: no multi-subject sections or files + assert result["multi"] is None and "被験者別" not in result["markdown"] + assert not (tmp_path / "out" / "accuracy_loso.csv").exists() def test_flip_sign_is_flagged_by_negative_slope(tmp_path: Path) -> None: @@ -137,3 +146,97 @@ def test_cli_end_to_end_with_glob_and_plots(tmp_path: Path, capsys: pytest.Captu assert "しきい値" in capsys.readouterr().out assert (out / "plots" / "accuracy_scatter.png").is_file() # matplotlib is in requirements.txt assert main(["--trials", str(trials), "--face-csv", str(tmp_path / "missing.csv")]) == 2 + + +# --------------------------------------------------------------------------------------------- +# multiple subjects +# --------------------------------------------------------------------------------------------- +def _two_subjects(tmp_path: Path, scale_b: float = 0.80): + ta, fa = _simulate(tmp_path / "a", scale=0.95, offset=-19.0, hour=12) + tb, fb = _simulate(tmp_path / "b", scale=scale_b, offset=-5.0, hour=14, seed=1) + return ta, fa, tb, fb + + +def test_split_label() -> None: + assert split_label("A=trials_a.csv") == ("A", "trials_a.csv") + assert split_label(r"C:\data\x=1.csv") == ("", r"C:\data\x=1.csv") + assert split_label("trials.csv") == ("", "trials.csv") + + +def test_subject_labels_aliases_and_exclusions(tmp_path: Path) -> None: + header = "session_id,seq,block,condition_id,distance_m,true_yaw_deg,pitch_deg,repeat,subject,status,t_baseline_start,t_turn_cue,t_hold_end\n" + a, b = tmp_path / "a.csv", tmp_path / "b.csv" + a.write_text(header + "S,1,1,c,1,10,0,1,S01,ok,a,b,c\nS,2,1,c,1,10,0,2,S01,ok,a,b,c\n", encoding="utf-8-sig") + b.write_text(header + "T,1,1,c,1,10,0,1,,ok,a,b,c\nT,2,1,c,1,10,0,2,,ok,a,b,c\n", encoding="utf-8-sig") + assert {r["subject"] for r in load_trials([a, b])} == {"S01", "b"} # column, then file name + rows = load_trials([a, b], labels={b: "S02"}, aliases={"S01": "A"}) + assert {r["subject"] for r in rows} == {"A", "S02"} and len(rows) == 4 + kept = load_trials([a, b], exclude_seqs=["S02:2", 1], labels={b: "S02"}) + assert [(r["subject"], r["seq"]) for r in kept] == [("S01", "2")] + only = load_trials([a, b], exclude_seqs=["S02:2"], labels={b: "S02"}) + assert [(r["subject"], r["seq"]) for r in only] == [("S01", "1"), ("S01", "2"), ("S02", "1")] + + +def test_multi_subject_pools_frames_and_reports_each_subject(tmp_path: Path) -> None: + ta, fa, tb, fb = _two_subjects(tmp_path) + result = run_report([ta, tb], [fa, fb], tmp_path / "out", Options(), plots=False, labels={ta: "A", tb: "B"}) + assert result["subjects"] == ["A", "B"] and len(result["records"]) == 60 + assert all(r["status"] == "ok" for r in result["records"]) + multi = result["multi"] + per = multi["per_subject"] + assert abs(per["A"]["thresholds"][1.0]["slope"] - 0.95) < 0.05 + assert abs(per["B"]["thresholds"][1.0]["slope"] - 0.80) < 0.05 + pooled = result["thresholds"][1.0] + assert pooled["n_subjects"] == 2 and 0.80 < pooled["slope"] < 0.95 + cond = {(c["distance_m"], c["true_yaw_deg"]): c for c in result["conditions"]} + assert cond[(1.0, 10.0)]["n_subjects"] == 2 and cond[(1.0, 10.0)]["n_trials"] == 6 + assert cond[(1.0, 10.0)]["sd_subjects_deg"] > 0.5 # +10 deg reads 9.5 (A) vs 8.0 (B) + # leave-one-subject-out: one row per held-out subject and distance + loso = multi["loso"] + assert {(r["held_out"], r["distance_m"]) for r in loso} == {("A", 1.0), ("A", 3.0), ("B", 1.0), ("B", 3.0)} + assert all(r["n_train_subjects"] == 1 for r in loso) + # the +/- gain is symmetric in this simulation; the per-subject gains differ + asym = {r["subject"]: r for r in multi["asymmetry"]} + assert set(asym) == {"A", "B", "ALL"} and abs(asym["ALL"]["gain_diff"]) < 0.05 + assert asym["A"]["gain_pos"] > asym["B"]["gain_pos"] + 0.05 + assert multi["differences"] and all(r["p_permutation"] < 0.05 for r in multi["differences"]) + for name in ("accuracy_subject_conditions.csv", "accuracy_subject_thresholds.csv", "accuracy_loso.csv", + "accuracy_asymmetry.csv", "accuracy_asymmetry_cells.csv", "accuracy_subject_differences.csv"): + assert (tmp_path / "out" / name).is_file(), name + header = (tmp_path / "out" / "accuracy_conditions.csv").read_text(encoding="utf-8-sig").splitlines()[0] + assert header.endswith("exceed_rate,n_subjects,sd_subjects_deg") + for text in ("被験者別しきい値", "leave-one-subject-out", "正負方向の非対称", "A (30試行)", "B (30試行)"): + assert text in result["markdown"] + + +def test_flip_only_the_named_subject(tmp_path: Path) -> None: + ta, fa, tb, fb = _two_subjects(tmp_path, scale_b=-0.95) + labels = {ta: "A", tb: "B"} + raw = run_report([ta, tb], [fa, fb], tmp_path / "o1", Options(), plots=False, labels=labels) + assert raw["multi"]["per_subject"]["B"]["thresholds"][1.0]["slope"] < 0 + assert "--flip-subject B" in raw["markdown"] + fixed = run_report([ta, tb], [fa, fb], tmp_path / "o2", Options(flip_subjects=("B",)), plots=False, labels=labels) + per = fixed["multi"]["per_subject"] + assert per["A"]["thresholds"][1.0]["slope"] > 0.9 and per["B"]["thresholds"][1.0]["slope"] > 0.9 + assert "--flip-subject B" not in fixed["markdown"] + only_a = run_report([ta, tb], [fa, fb], tmp_path / "o3", Options(), plots=False, labels=labels, + exclude_subjects=["B"]) + assert only_a["subjects"] == ["A"] and only_a["multi"] is None + + +def test_cli_multi_subject_with_labels_and_subject_plots(tmp_path: Path, capsys: pytest.CaptureFixture[str]) -> None: + ta, fa, tb, fb = _two_subjects(tmp_path) + out = tmp_path / "report" + assert main(["--trials", f"A={ta}", f"B={tb}", "--face-csv", str(fa), str(fb), "--out", str(out)]) == 0 + assert "被験者別しきい値" in capsys.readouterr().out + for name in ("accuracy_scatter.png", "accuracy_scatter_by_subject.png", "accuracy_error_by_subject.png", + "accuracy_gain_by_side.png"): + assert (out / "plots" / name).is_file(), name + assert main(["--trials", f"A={ta}", "--face-csv", str(fa), "--alias", "broken"]) == 2 + + +def test_permutation_helpers() -> None: + assert _signflip_p([1.0, 1.0, 1.0, 1.0]) == pytest.approx(2 / 16) + assert math.isnan(_signflip_p([0.0, 0.0, 0.0])) + h, p = kruskal_permutation([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) + assert h == pytest.approx(3.857, abs=0.01) and 0.05 < p < 0.2 diff --git a/tests/test_accuracy_report.py.orig b/tests/test_accuracy_report.py.orig new file mode 100644 index 0000000..fc5887c --- /dev/null +++ b/tests/test_accuracy_report.py.orig @@ -0,0 +1,139 @@ +import csv +import math +import random +from datetime import datetime, timedelta, timezone +from pathlib import Path + +import pytest + +from src.accuracy_report import ( + Options, evaluate_trial, load_face_series, load_trials, main, percentile, run_report, +) +from src.accuracy_session import SessionRunner, build_plan + +TZ = timezone(timedelta(hours=9)) +FPS = 30.0 + + +class FakeTime: + def __init__(self) -> None: + self.now = datetime(2026, 10, 4, 12, 0, 0, tzinfo=TZ) + + def sleep(self, seconds: float) -> None: + self.now += timedelta(seconds=seconds) + + def clock(self) -> str: + return self.now.isoformat(timespec="milliseconds") + + +def _simulate(tmp_path: Path, scale: float = 0.95, offset: float = -19.0, noise: float = 0.5, face_dropout: float = 0.0): + """Run a fake session and write a matching face_frames CSV (30 fps, yaw ramps after each cue).""" + plan = build_plan([1, 3], [0, 5, -5, 10, -10], 3, seed=3) + fake, trials = FakeTime(), tmp_path / "trials.csv" + start = fake.now + + def operator_pause(_prompt: str) -> str: + fake.sleep(4.0) # time the subject needs to return and the operator needs to press Enter + return "" + + SessionRunner(plan, trials, "S", input_fn=operator_pause, out=lambda _m: None, sleep=fake.sleep, clock=fake.clock, + beep=lambda _k: None).run() + end = fake.now + timedelta(seconds=2) + rng = random.Random(0) + marks = [(datetime.fromisoformat(r["t_turn_cue"]), float(r["true_yaw_deg"]), datetime.fromisoformat(r["t_hold_end"])) + for r in csv.DictReader(trials.open(encoding="utf-8-sig"))] + path = tmp_path / "20261004_120000_exp_face_frames.csv" + with path.open("w", encoding="utf-8-sig", newline="") as stream: + writer = csv.writer(stream) + writer.writerow(["frame_number", "timestamp", "track_id", "pose_estimated", "yaw", "raw_yaw"]) + t, frame = start - timedelta(seconds=1), 0 + while t < end: + frame += 1 + true = 0.0 + for cue, angle, stop in marks: + if cue <= t <= stop: + true = angle * min(1.0, (t - cue).total_seconds() / 0.8) + yaw = offset + scale * true + rng.gauss(0, noise) + ok = rng.random() >= face_dropout + writer.writerow([frame, t.isoformat(timespec="milliseconds"), 1, ok, f"{yaw:.4f}" if ok else "", f"{yaw:.4f}" if ok else ""]) + t += timedelta(seconds=1 / FPS) + return trials, path + + +def test_percentile() -> None: + assert percentile([0.0, 10.0], 50) == 5.0 and math.isnan(percentile([], 50)) + + +def test_full_pipeline_recovers_scale_bias_and_threshold(tmp_path: Path) -> None: + trials, face = _simulate(tmp_path) + result = run_report([trials], [face], tmp_path / "out", Options(), plots=False) + records = result["records"] + assert len(records) == 30 and all(r["status"] == "ok" for r in records) + # a constant -19 deg offset must cancel; 0.95 scale makes +10 read about 9.5 -> bias about -0.5 + cond = {(c["distance_m"], c["true_yaw_deg"]): c for c in result["conditions"]} + assert abs(cond[(1.0, 0.0)]["bias_deg"]) < 0.3 + assert abs(cond[(3.0, 10.0)]["bias_deg"] + 0.5) < 0.3 + th = result["thresholds"][1.0] + assert 0.95 - 0.05 < th["slope"] < 0.95 + 0.05 and th["r2"] > 0.99 + assert 1.0 < th["theta_th_deg"] < 3.5 + assert th["theta_worst_deg"] >= th["theta_pooled_deg"] - 1e-9 + assert th["min_detectable_deg"] == 5.0 # 5 deg turns exceed theta_th in >= 95 % of frames + assert cond[(1.0, 0.0)]["exceed_rate"] < 0.05 # false-positive rate at 0 deg + assert cond[(1.0, 10.0)]["mean_detection_rate"] > 0.95 + for name in ("accuracy_trials.csv", "accuracy_conditions.csv", "accuracy_thresholds.csv", "accuracy_report.md"): + assert (tmp_path / "out" / name).is_file() + assert "θ_th" in result["markdown"] + + +def test_flip_sign_is_flagged_by_negative_slope(tmp_path: Path) -> None: + trials, face = _simulate(tmp_path, scale=-1.0) + result = run_report([trials], [face], tmp_path / "out", Options(), plots=False) + assert result["thresholds"][1.0]["slope"] < 0 and "--flip-sign" in result["markdown"] + fixed = run_report([trials], [face], tmp_path / "out2", Options(flip_sign=True), plots=False) + assert fixed["thresholds"][1.0]["slope"] > 0.9 + + +def test_dropout_lowers_detection_rate_and_missing_data_is_excluded(tmp_path: Path) -> None: + trials, face = _simulate(tmp_path, face_dropout=0.3) + result = run_report([trials], [face], tmp_path / "out", Options(), plots=False) + rates = [r["detection_rate"] for r in result["records"] if r["status"] == "ok"] + assert 0.55 < sum(rates) / len(rates) < 0.85 + # a face log that ends early leaves later trials without any track + series = load_face_series([face]) + cut = series.times[len(series.times) // 2] + short = tmp_path / "short_face_frames.csv" + with face.open(encoding="utf-8-sig", newline="") as src, short.open("w", encoding="utf-8-sig", newline="") as dst: + reader, writer = csv.reader(src), csv.writer(dst) + for index, row in enumerate(reader): + if index == 0 or datetime.fromisoformat(row[1]).timestamp() <= cut: + writer.writerow(row) + partial = run_report([trials], [short], tmp_path / "out3", Options(), plots=False) + statuses = {r["status"] for r in partial["records"]} + assert "no_track" in statuses and "ok" in statuses + + +def test_last_row_wins_and_exclude_seq(tmp_path: Path) -> None: + path = tmp_path / "t.csv" + header = "session_id,seq,block,condition_id,distance_m,true_yaw_deg,pitch_deg,repeat,subject,status,t_baseline_start,t_turn_cue,t_hold_end\n" + path.write_text(header + "S,1,1,c,1,10,0,1,,ok,a,b,c\nS,1,1,c,1,10,0,1,,discarded,,,\nS,2,1,c,1,10,0,2,,ok,a,b,c\n" + "S,3,1,c,1,10,0,3,,ok,a,b,c\n", encoding="utf-8-sig") + assert [r["seq"] for r in load_trials([path])] == ["2", "3"] + assert [r["seq"] for r in load_trials([path], exclude_seqs=[3])] == ["2"] + + +def test_bad_timestamps_do_not_crash(tmp_path: Path) -> None: + trials, face = _simulate(tmp_path) + series = load_face_series([face]) + record, deltas = evaluate_trial({"seq": "1", "condition_id": "x", "distance_m": "1", "true_yaw_deg": "5", + "pitch_deg": "0", "t_baseline_start": "", "t_turn_cue": "", "t_hold_end": ""}, + series, Options()) + assert record["status"] == "bad_timestamps" and deltas == [] + + +def test_cli_end_to_end_with_glob_and_plots(tmp_path: Path, capsys: pytest.CaptureFixture[str]) -> None: + trials, face = _simulate(tmp_path) + out = tmp_path / "report" + assert main(["--trials", str(trials), "--face-csv", str(tmp_path / "*_face_frames.csv"), "--out", str(out)]) == 0 + assert "しきい値" in capsys.readouterr().out + assert (out / "plots" / "accuracy_scatter.png").is_file() # matplotlib is in requirements.txt + assert main(["--trials", str(trials), "--face-csv", str(tmp_path / "missing.csv")]) == 2