from __future__ import annotations

import csv
import hashlib
import json
import math
import statistics
from collections import defaultdict
from datetime import date, datetime, timedelta
from pathlib import Path


ROOT = Path(__file__).resolve().parents[1]
SOURCE = ROOT / "source_tables"

RUNS_PATH = SOURCE / "v7_pillar_conserved_mechanics_runs_v1.csv"
TIMELINE_PATH = SOURCE / "v7_total_fitness_profile_timeline_v3.csv"
FOOD_PATH = SOURCE / "v7_food_pdf_garmin_overlap_v1.csv"


def read_csv(path: Path) -> list[dict[str, str]]:
    with path.open("r", newline="", encoding="utf-8-sig") as handle:
        return list(csv.DictReader(handle))


def write_csv(path: Path, rows: list[dict[str, object]], headers: list[str]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.DictWriter(handle, fieldnames=headers)
        writer.writeheader()
        writer.writerows(rows)


def write_text(path: Path, text: str) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(text, encoding="utf-8")


def to_float(value: str | None) -> float | None:
    if value is None:
        return None
    value = value.strip()
    if not value or value.lower() == "nan":
        return None
    return float(value)


def mean(values: list[float]) -> float:
    return statistics.fmean(values)


def sample_std(values: list[float]) -> float:
    if len(values) < 2:
        return 0.0
    return statistics.stdev(values)


def pct_change(first: float, last: float) -> float:
    if first == 0:
        return math.nan
    return ((last - first) / abs(first)) * 100.0


def fmt(value: object, digits: int = 2) -> str:
    if value is None:
        return ""
    if isinstance(value, float):
        if math.isnan(value):
            return ""
        return f"{value:.{digits}f}"
    return str(value)


def format_metric_value(metric: str, value: object) -> str:
    if value is None:
        return ""
    numeric = float(value)
    if metric in {"speed_per_hr_mean", "speed_per_hr"}:
        return f"{numeric:.5f}"
    if metric in {"power_per_hr_mean", "power_per_hr", "speed_mps_mean", "speed_mps_norm", "hr_residual_bpm_mean", "hr_residual_bpm"}:
        return f"{numeric:.5f}"
    return f"{numeric:.2f}"


def markdown_table(rows: list[dict[str, object]], headers: list[str]) -> str:
    if not rows:
        return "_No rows_"
    lines = []
    lines.append("| " + " | ".join(headers) + " |")
    lines.append("|" + "|".join(["---"] * len(headers)) + "|")
    for row in rows:
        values = []
        for header in headers:
            value = row.get(header, "")
            values.append(fmt(value) if isinstance(value, float) else str(value))
        lines.append("| " + " | ".join(values) + " |")
    return "\n".join(lines)


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


def build_source_manifest(paths: list[Path]) -> list[dict[str, object]]:
    rows = []
    for path in paths:
        rows.append(
            {
                "table_name": path.name,
                "relative_path": str(path.relative_to(ROOT)).replace("\\", "/"),
                "bytes": path.stat().st_size,
                "sha256": sha256(path),
            }
        )
    return rows


def correlation(xs: list[float | None], ys: list[float | None]) -> tuple[float | None, int]:
    points = [(x, y) for x, y in zip(xs, ys) if x is not None and y is not None]
    if len(points) < 3:
        return None, len(points)
    xs_clean = [x for x, _ in points]
    ys_clean = [y for _, y in points]
    mx = mean(xs_clean)
    my = mean(ys_clean)
    sx = sum((x - mx) ** 2 for x in xs_clean)
    sy = sum((y - my) ** 2 for y in ys_clean)
    if sx == 0 or sy == 0:
        return None, len(points)
    cov = sum((x - mx) * (y - my) for x, y in points)
    return cov / math.sqrt(sx * sy), len(points)


def solve_ols_hr_model(run_rows: list[dict[str, str]]) -> dict[str, float]:
    xs = [(1.0, float(row["speed_mps_norm"]), float(row["avg_power_w"])) for row in run_rows]
    ys = [float(row["avg_hr_bpm"]) for row in run_rows]
    matrix = [[0.0] * 3 for _ in range(3)]
    vector = [0.0] * 3
    for xi, yi in zip(xs, ys):
        for i in range(3):
            vector[i] += xi[i] * yi
            for j in range(3):
                matrix[i][j] += xi[i] * xi[j]
    augmented = [row[:] + [b] for row, b in zip(matrix, vector)]
    for i in range(3):
        pivot = max(range(i, 3), key=lambda idx: abs(augmented[idx][i]))
        augmented[i], augmented[pivot] = augmented[pivot], augmented[i]
        factor = augmented[i][i]
        for j in range(i, 4):
            augmented[i][j] /= factor
        for r in range(3):
            if r == i:
                continue
            factor = augmented[r][i]
            for j in range(i, 4):
                augmented[r][j] -= factor * augmented[i][j]
    return {
        "intercept": augmented[0][3],
        "speed_coef": augmented[1][3],
        "power_coef": augmented[2][3],
    }


def build_run_level_dataset(run_rows: list[dict[str, str]], timeline_by_date: dict[str, dict[str, str]]) -> list[dict[str, object]]:
    coeffs = solve_ols_hr_model(run_rows)
    rows = []
    for row in run_rows:
        calendar_date = row["calendar_date"]
        same_day = timeline_by_date.get(calendar_date)
        next_day = timeline_by_date.get((date.fromisoformat(calendar_date) + timedelta(days=1)).isoformat())
        speed = float(row["speed_mps_norm"])
        hr = float(row["avg_hr_bpm"])
        power = float(row["avg_power_w"])
        predicted_hr = coeffs["intercept"] + coeffs["speed_coef"] * speed + coeffs["power_coef"] * power
        total_28d = to_float(same_day["total_activity_hours_28d_rebuilt"]) if same_day else None
        run_28d = to_float(same_day["running_hours_28d"]) if same_day else None
        run_share = (run_28d / total_28d) * 100.0 if run_28d is not None and total_28d not in {None, 0} else None
        rows.append(
            {
                "calendar_date": calendar_date,
                "year_month": row["year_month"],
                "rebuild_modality": row["rebuild_modality"],
                "distance_miles_norm": round(float(row["distance_miles_norm"]), 5),
                "speed_mps_norm": round(speed, 5),
                "pace_min_per_mile_norm": round(float(row["pace_min_per_mile_norm"]), 5),
                "avg_hr_bpm": round(hr, 5),
                "avg_power_w": round(power, 5),
                "cadence_spm_est": round(float(row["cadence_spm_est"]), 5),
                "stride_length_m_est": round(float(row["stride_length_m_est"]), 5),
                "vertical_ratio_pct": round(float(row["vertical_ratio_pct"]), 5),
                "speed_per_hr": round(speed / hr, 8),
                "power_per_hr": round(power / hr, 8),
                "predicted_hr_from_speed_power": round(predicted_hr, 5),
                "hr_residual_bpm": round(hr - predicted_hr, 5),
                "running_hours_28d": round(run_28d, 5) if run_28d is not None else None,
                "total_activity_hours_28d_rebuilt": round(total_28d, 5) if total_28d is not None else None,
                "running_share_28d_activity_pct": round(run_share, 5) if run_share is not None else None,
                "next_day_resting_hr": round(to_float(next_day["resting_heart_rate"]), 5) if next_day and next_day["resting_heart_rate"] else None,
                "next_day_sleep_score": round(to_float(next_day["sleep_score"]), 5) if next_day and next_day["sleep_score"] else None,
                "next_day_hrv": round(to_float(next_day["health_HRV"]), 5) if next_day and next_day["health_HRV"] else None,
            }
        )
    return rows


def summarize_window(rows: list[dict[str, object]], label: str) -> dict[str, object]:
    metrics = {
        "speed_mps_mean": "speed_mps_norm",
        "pace_min_per_mile_mean": "pace_min_per_mile_norm",
        "avg_hr_bpm_mean": "avg_hr_bpm",
        "avg_power_w_mean": "avg_power_w",
        "cadence_spm_mean": "cadence_spm_est",
        "stride_length_m_mean": "stride_length_m_est",
        "speed_per_hr_mean": "speed_per_hr",
        "power_per_hr_mean": "power_per_hr",
        "hr_residual_bpm_mean": "hr_residual_bpm",
        "running_hours_28d_mean": "running_hours_28d",
        "total_activity_hours_28d_mean": "total_activity_hours_28d_rebuilt",
        "running_share_28d_activity_pct_mean": "running_share_28d_activity_pct",
        "next_day_resting_hr_mean": "next_day_resting_hr",
        "next_day_sleep_score_mean": "next_day_sleep_score",
        "next_day_hrv_mean": "next_day_hrv",
    }
    output: dict[str, object] = {
        "window_label": label,
        "run_count": len(rows),
        "start_date": rows[0]["calendar_date"],
        "end_date": rows[-1]["calendar_date"],
    }
    for out_key, source_key in metrics.items():
        values = [row[source_key] for row in rows if row[source_key] is not None]
        output[out_key] = round(mean([float(value) for value in values]), 5) if values else None
    return output


def build_window_summary(run_level_rows: list[dict[str, object]]) -> list[dict[str, object]]:
    return [
        summarize_window(run_level_rows[:30], "early_qc_window"),
        summarize_window(run_level_rows[-30:], "late_qc_window"),
    ]


def build_window_change_rows(window_rows: list[dict[str, object]]) -> list[dict[str, object]]:
    early = next(row for row in window_rows if row["window_label"] == "early_qc_window")
    late = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    metrics = [
        "speed_mps_mean",
        "pace_min_per_mile_mean",
        "avg_hr_bpm_mean",
        "avg_power_w_mean",
        "cadence_spm_mean",
        "stride_length_m_mean",
        "speed_per_hr_mean",
        "power_per_hr_mean",
        "hr_residual_bpm_mean",
        "running_share_28d_activity_pct_mean",
        "next_day_resting_hr_mean",
        "next_day_sleep_score_mean",
    ]
    rows = []
    for metric in metrics:
        early_value = early[metric]
        late_value = late[metric]
        if early_value is None or late_value is None:
            continue
        rows.append(
            {
                "metric": metric,
                "early_value": round(float(early_value), 5),
                "late_value": round(float(late_value), 5),
                "change_value": round(float(late_value) - float(early_value), 5),
                "change_units": "delta_bpm" if metric == "hr_residual_bpm_mean" else "pct_change",
                "pct_change": round(pct_change(float(early_value), float(late_value)), 5),
            }
        )
    return rows


def build_monthly_summary(run_level_rows: list[dict[str, object]]) -> list[dict[str, object]]:
    by_month: dict[str, list[dict[str, object]]] = defaultdict(list)
    for row in run_level_rows:
        by_month[str(row["year_month"])].append(row)
    rows = []
    for year_month in sorted(by_month):
        month_rows = by_month[year_month]
        if len(month_rows) < 3:
            continue
        summary = {
            "year_month": year_month,
            "qc_run_count": len(month_rows),
        }
        for metric in [
            "speed_mps_norm",
            "avg_hr_bpm",
            "avg_power_w",
            "speed_per_hr",
            "power_per_hr",
            "hr_residual_bpm",
            "running_share_28d_activity_pct",
            "next_day_resting_hr",
            "next_day_sleep_score",
            "next_day_hrv",
        ]:
            values = [row[metric] for row in month_rows if row[metric] is not None]
            summary[metric] = round(mean([float(value) for value in values]), 5) if values else None
        rows.append(summary)
    return rows


def build_specialization_correlations(run_level_rows: list[dict[str, object]]) -> list[dict[str, object]]:
    run_share = [to_float(str(row["running_share_28d_activity_pct"])) if row["running_share_28d_activity_pct"] is not None else None for row in run_level_rows]
    specs = [
        ("running_share_vs_speed_per_hr", run_share, [row["speed_per_hr"] for row in run_level_rows]),
        ("running_share_vs_power_per_hr", run_share, [row["power_per_hr"] for row in run_level_rows]),
        ("running_share_vs_hr_residual", run_share, [row["hr_residual_bpm"] for row in run_level_rows]),
        ("running_share_vs_nextday_resting_hr", run_share, [row["next_day_resting_hr"] for row in run_level_rows]),
        ("running_share_vs_nextday_sleep_score", run_share, [row["next_day_sleep_score"] for row in run_level_rows]),
        ("running_share_vs_nextday_hrv", run_share, [row["next_day_hrv"] for row in run_level_rows]),
    ]
    rows = []
    for label, xs, ys in specs:
        r_value, n_points = correlation(xs, ys)
        rows.append(
            {
                "relationship": label,
                "correlation_r": round(r_value, 5) if r_value is not None else None,
                "n_points": n_points,
            }
        )
    return rows


def build_hr_model_summary(run_level_rows: list[dict[str, object]]) -> list[dict[str, object]]:
    coeffs = solve_ols_hr_model(
        [
            {
                "speed_mps_norm": row["speed_mps_norm"],
                "avg_power_w": row["avg_power_w"],
                "avg_hr_bpm": row["avg_hr_bpm"],
            }
            for row in run_level_rows
        ]
    )
    early = run_level_rows[:30]
    late = run_level_rows[-30:]
    early_resid = [float(row["hr_residual_bpm"]) for row in early if row["hr_residual_bpm"] is not None]
    late_resid = [float(row["hr_residual_bpm"]) for row in late if row["hr_residual_bpm"] is not None]
    return [
        {
            "model_name": "hr_from_speed_power",
            "intercept": round(coeffs["intercept"], 5),
            "speed_coef": round(coeffs["speed_coef"], 5),
            "power_coef": round(coeffs["power_coef"], 5),
            "early_mean_residual_bpm": round(mean(early_resid), 5),
            "late_mean_residual_bpm": round(mean(late_resid), 5),
            "early_residual_sd_bpm": round(sample_std(early_resid), 5),
            "late_residual_sd_bpm": round(sample_std(late_resid), 5),
        }
    ]


def build_food_context(food_rows: list[dict[str, str]]) -> tuple[list[dict[str, object]], list[dict[str, object]]]:
    food_item_days = [row for row in food_rows if to_float(row["food_item_count"]) is not None and to_float(row["food_item_count"]) > 0]
    summary = [
        {
            "subset": "all_overlap_days",
            "n_days": len(food_rows),
            "mean_food_calories": round(mean([to_float(row["food_total_calories"]) for row in food_rows if to_float(row["food_total_calories"]) is not None]), 5),
            "mean_net_carbs_g": round(mean([to_float(row["food_total_net_carbs_g"]) for row in food_rows if to_float(row["food_total_net_carbs_g"]) is not None]), 5),
            "mean_sleep_score": round(mean([to_float(row["sleep_score"]) for row in food_rows if to_float(row["sleep_score"]) is not None]), 5),
            "mean_hrv": round(mean([to_float(row["health_HRV"]) for row in food_rows if to_float(row["health_HRV"]) is not None]), 5),
            "mean_running_hours_7d": round(mean([to_float(row["running_hours_7d"]) for row in food_rows if to_float(row["running_hours_7d"]) is not None]), 5),
        },
        {
            "subset": "food_item_days_only",
            "n_days": len(food_item_days),
            "mean_food_calories": round(mean([to_float(row["food_total_calories"]) for row in food_item_days if to_float(row["food_total_calories"]) is not None]), 5),
            "mean_net_carbs_g": round(mean([to_float(row["food_total_net_carbs_g"]) for row in food_item_days if to_float(row["food_total_net_carbs_g"]) is not None]), 5),
            "mean_sleep_score": round(mean([to_float(row["sleep_score"]) for row in food_item_days if to_float(row["sleep_score"]) is not None]), 5),
            "mean_hrv": round(mean([to_float(row["health_HRV"]) for row in food_item_days if to_float(row["health_HRV"]) is not None]), 5),
            "mean_running_hours_7d": round(mean([to_float(row["running_hours_7d"]) for row in food_item_days if to_float(row["running_hours_7d"]) is not None]), 5),
        },
    ]
    corr_specs = [
        ("food_calories_vs_running_hours_7d", [to_float(row["food_total_calories"]) for row in food_rows], [to_float(row["running_hours_7d"]) for row in food_rows]),
        ("net_carbs_vs_sleep_score", [to_float(row["food_total_net_carbs_g"]) for row in food_rows], [to_float(row["sleep_score"]) for row in food_rows]),
        ("food_calories_vs_hrv", [to_float(row["food_total_calories"]) for row in food_rows], [to_float(row["health_HRV"]) for row in food_rows]),
    ]
    correlations = []
    for label, xs, ys in corr_specs:
        r_value, n_points = correlation(xs, ys)
        correlations.append(
            {
                "relationship": label,
                "correlation_r": round(r_value, 5) if r_value is not None else None,
                "n_points": n_points,
            }
        )
    return summary, correlations


def build_scope_claims() -> list[dict[str, object]]:
    return [
        {
            "category": "answered",
            "question": "Did adaptive efficiency improve as specialization increased in the later running window?",
            "answer": "Yes, at least at the speed-per-heart-rate level. The late window shows better external speed per unit heart-rate cost than the early window.",
        },
        {
            "category": "answered",
            "question": "Did background recovery markers move in a favorable direction as specialization increased?",
            "answer": "Yes, descriptively for next-day resting heart rate and next-day sleep score. HRV is late-window only and should be treated as contextual rather than symmetric early-late evidence.",
        },
        {
            "category": "supported_not_proven",
            "question": "Did unexplained session HR burden disappear as adaptation advanced?",
            "answer": "No. Session heart rate remained high, and a simple speed-plus-power HR model shows slightly positive late-window unexplained HR burden rather than a full disappearance of burden.",
        },
        {
            "category": "not_answered",
            "question": "Can this study prove causal physiology or exact fatigue mechanism?",
            "answer": "No. The study characterizes a mixed pattern of improved efficiency and persistent session-level burden, but it does not prove mechanism.",
        },
    ]


def svg_canvas(width: int, height: int, body: str) -> str:
    return (
        f'<svg xmlns="http://www.w3.org/2000/svg" width="{width}" height="{height}" '
        f'viewBox="0 0 {width} {height}"><rect width="100%" height="100%" fill="#f7f5ef"/>'
        f"{body}</svg>"
    )


def svg_text(x: float, y: float, text: str, size: int = 12, weight: str = "normal", anchor: str = "start", fill: str = "#222") -> str:
    safe = text.replace("&", "&amp;").replace("<", "&lt;").replace(">", "&gt;")
    return f'<text x="{x:.1f}" y="{y:.1f}" font-family="Arial" font-size="{size}" font-weight="{weight}" text-anchor="{anchor}" fill="{fill}">{safe}</text>'


def svg_rect(x: float, y: float, w: float, h: float, fill: str) -> str:
    return f'<rect x="{x:.1f}" y="{y:.1f}" width="{w:.1f}" height="{h:.1f}" fill="{fill}" />'


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


def svg_polyline(points: list[tuple[float, float]], stroke: str, width: float = 2.0, fill: str = "none") -> str:
    point_str = " ".join(f"{x:.1f},{y:.1f}" for x, y in points)
    return f'<polyline points="{point_str}" fill="{fill}" stroke="{stroke}" stroke-width="{width:.1f}" />'


def build_figure_window_pct_change(change_rows: list[dict[str, object]], path: Path) -> None:
    focus = {
        "speed_mps_mean": "#1768ac",
        "avg_hr_bpm_mean": "#f26419",
        "avg_power_w_mean": "#33658a",
        "speed_per_hr_mean": "#2a9d8f",
        "power_per_hr_mean": "#86bbd8",
        "next_day_resting_hr_mean": "#6d597a",
        "next_day_sleep_score_mean": "#8ab17d",
    }
    selected = [row for row in change_rows if row["metric"] in focus]
    width, height = 960, 520
    left, top, chart_w, chart_h = 90, 80, 760, 320
    max_abs = max(abs(float(row["pct_change"])) for row in selected)
    body = []
    body.append(svg_text(width / 2, 30, "Figure 1. Early-to-Late Change In Output, Cost, And Recovery", 21, "bold", "middle"))
    body.append(svg_text(width / 2, 52, "Efficiency improved while session HR remained elevated", 13, "normal", "middle", "#555"))
    mid_y = top + chart_h / 2
    body.append(svg_line(left, mid_y, left + chart_w, mid_y, "#444", 1.2))
    body.append(svg_line(left, top, left, top + chart_h, "#444", 1.2))
    gap = chart_w / max(len(selected), 1)
    bar_w = gap * 0.62
    for idx, row in enumerate(selected):
        x = left + idx * gap + (gap - bar_w) / 2
        pct = float(row["pct_change"])
        h = (abs(pct) / max_abs) * (chart_h / 2 - 18) if max_abs else 0.0
        y = mid_y - h if pct >= 0 else mid_y
        body.append(svg_rect(x, y, bar_w, h, focus[row["metric"]]))
        body.append(svg_text(x + bar_w / 2, top + chart_h + 26, row["metric"].replace("_mean", "").replace("_", " "), 10, "normal", "middle"))
        body.append(svg_text(x + bar_w / 2, y - 6 if pct >= 0 else y + h + 14, f"{pct:.1f}%", 10, "normal", "middle"))
    write_text(path, svg_canvas(width, height, "".join(body)))


def build_figure_monthly_efficiency(monthly_rows: list[dict[str, object]], path: Path) -> None:
    width, height = 980, 520
    left, top, chart_w, chart_h = 100, 90, 760, 300
    xs = list(range(len(monthly_rows)))
    run_share = [float(row["running_share_28d_activity_pct"]) for row in monthly_rows if row["running_share_28d_activity_pct"] is not None]
    speed_eff = [float(row["speed_per_hr"]) for row in monthly_rows if row["speed_per_hr"] is not None]
    max_run_share = max(run_share) if run_share else 1.0
    min_eff = min(speed_eff) if speed_eff else 0.0
    max_eff = max(speed_eff) if speed_eff else 1.0
    body = []
    body.append(svg_text(width / 2, 30, "Figure 2. Monthly Specialization And Speed Efficiency", 21, "bold", "middle"))
    body.append(svg_text(width / 2, 52, "Speed per heart-rate cost rose as running share increased", 13, "normal", "middle", "#555"))
    body.append(svg_line(left, top + chart_h, left + chart_w, top + chart_h, "#444", 1.2))
    body.append(svg_line(left, top, left, top + chart_h, "#444", 1.2))
    eff_points = []
    share_points = []
    for idx, row in enumerate(monthly_rows):
        x = left + (idx / max(len(monthly_rows) - 1, 1)) * chart_w
        eff = float(row["speed_per_hr"])
        share = float(row["running_share_28d_activity_pct"])
        eff_y = top + chart_h - ((eff - min_eff) / max(max_eff - min_eff, 1e-9)) * chart_h
        share_y = top + chart_h - (share / max_run_share) * chart_h
        eff_points.append((x, eff_y))
        share_points.append((x, share_y))
        body.append(svg_text(x, top + chart_h + 22, row["year_month"], 9, "normal", "middle", "#555"))
    body.append(svg_polyline(eff_points, "#1768ac", 2.5))
    body.append(svg_polyline(share_points, "#f26419", 2.5))
    body.append(svg_text(820, 110, "Speed/HR", 11, "bold", "start", "#1768ac"))
    body.append(svg_text(820, 132, "Running share 28d", 11, "bold", "start", "#f26419"))
    write_text(path, svg_canvas(width, height, "".join(body)))


def build_figure_hr_residual(monthly_rows: list[dict[str, object]], path: Path) -> None:
    width, height = 920, 460
    left, top, chart_w, chart_h = 90, 90, 720, 260
    vals = [float(row["hr_residual_bpm"]) for row in monthly_rows if row["hr_residual_bpm"] is not None]
    min_v, max_v = min(vals), max(vals)
    body = []
    body.append(svg_text(width / 2, 30, "Figure 3. Monthly Session HR Residual", 21, "bold", "middle"))
    body.append(svg_text(width / 2, 52, "Residual burden did not disappear even as efficiency improved", 13, "normal", "middle", "#555"))
    zero_y = top + chart_h - ((0 - min_v) / max(max_v - min_v, 1e-9)) * chart_h
    body.append(svg_line(left, zero_y, left + chart_w, zero_y, "#444", 1.2))
    body.append(svg_line(left, top, left, top + chart_h, "#444", 1.2))
    points = []
    for idx, row in enumerate(monthly_rows):
        x = left + (idx / max(len(monthly_rows) - 1, 1)) * chart_w
        value = float(row["hr_residual_bpm"])
        y = top + chart_h - ((value - min_v) / max(max_v - min_v, 1e-9)) * chart_h
        points.append((x, y))
        body.append(svg_text(x, top + chart_h + 22, row["year_month"], 9, "normal", "middle", "#555"))
    body.append(svg_polyline(points, "#6d597a", 2.5))
    write_text(path, svg_canvas(width, height, "".join(body)))


def build_figure_recovery(window_rows: list[dict[str, object]], path: Path) -> None:
    early = next(row for row in window_rows if row["window_label"] == "early_qc_window")
    late = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    metrics = [
        ("next_day_resting_hr_mean", "Next-day RHR", "#1768ac"),
        ("next_day_sleep_score_mean", "Next-day Sleep", "#f26419"),
        ("running_share_28d_activity_pct_mean", "Run share 28d %", "#2a9d8f"),
    ]
    width, height = 980, 500
    left, top, chart_w, chart_h = 100, 90, 760, 280
    max_value = max(float(row[key]) for row in [early, late] for key, _, _ in metrics if row[key] is not None)
    body = []
    body.append(svg_text(width / 2, 30, "Figure 4. Recovery And Specialization Shift", 21, "bold", "middle"))
    body.append(svg_text(width / 2, 52, "Background recovery improved as running became a larger share of activity", 13, "normal", "middle", "#555"))
    body.append(svg_line(left, top + chart_h, left + chart_w, top + chart_h, "#444", 1.2))
    body.append(svg_line(left, top, left, top + chart_h, "#444", 1.2))
    group_gap = chart_w / len(metrics)
    bar_w = group_gap * 0.22
    for idx, (key, label, color) in enumerate(metrics):
        center = left + idx * group_gap + group_gap / 2
        for offset, row, fill in [(-bar_w * 0.75, early, color), (bar_w * 0.75, late, "#555f6d")]:
            value = float(row[key]) if row[key] is not None else 0.0
            h = (value / max_value) * chart_h if max_value else 0.0
            x = center + offset - bar_w / 2
            y = top + chart_h - h
            body.append(svg_rect(x, y, bar_w, h, fill))
            body.append(svg_text(x + bar_w / 2, y - 6, f"{value:.1f}", 10, "normal", "middle"))
        body.append(svg_text(center, top + chart_h + 24, label, 10, "normal", "middle"))
    body.append(svg_rect(760, 100, 16, 16, "#1768ac"))
    body.append(svg_text(784, 113, "Early", 11))
    body.append(svg_rect(760, 126, 16, 16, "#555f6d"))
    body.append(svg_text(784, 139, "Late", 11))
    write_text(path, svg_canvas(width, height, "".join(body)))


def build_methods_markdown(source_manifest: list[dict[str, object]]) -> str:
    manifest_table = markdown_table(source_manifest, ["table_name", "relative_path", "bytes", "sha256"])
    return f"""# Study 000B Methods

## Design

Primary flagship companion study conducted after `Study 000A` to examine whether later specialization coincided with adaptive efficiency, persistent unexplained HR burden, or both.

## Relationship to Study 000A

`Study 000A` established that the system adapted through selective flexibility and specialization.

`Study 000B` does not reopen that question. It stands alongside the flagship and asks a second primary program question:

- did external output become more efficient relative to cardiovascular cost?
- did some background recovery markers move in a favorable direction as specialization increased?
- did session-level internal burden disappear, or did part of it persist?

## Core question

As running specialization increased in the later high-resolution window, did the system show better efficiency, persistent unexplained HR burden, or both together?

## Packaged source tables

{manifest_table}

## Analytic structure

1. QC-pass high-resolution runs were linked back to the unified daily timeline.
2. Run-level efficiency metrics were derived as `speed_per_hr` and `power_per_hr`.
3. Early and late 30-run windows were compared on output, cost, specialization, and next-day recovery. These windows were matched by run count rather than elapsed time, spanning `93` calendar days early and `43` calendar days late.
4. A simple descriptive HR model (`HR ~ speed + power`) was fit across QC-pass runs to estimate unexplained session HR burden after basic external-output adjustment.
5. Monthly summaries and specialization correlations were generated to characterize broader late-window trends.
6. The food overlap layer was retained only as partial exploratory late-window context.

## Interpretation boundary

This study is descriptive and within-subject. It does not prove physiology mechanism or clinical causation. It asks whether the data support a mixed pattern of better efficiency together with persistent unexplained HR burden.
"""


def build_results_markdown(
    window_rows: list[dict[str, object]],
    change_rows: list[dict[str, object]],
    hr_model_rows: list[dict[str, object]],
    monthly_rows: list[dict[str, object]],
    corr_rows: list[dict[str, object]],
    food_summary: list[dict[str, object]],
    food_corr_rows: list[dict[str, object]],
    scope_rows: list[dict[str, object]],
) -> str:
    window_display = []
    for row in window_rows:
        display = dict(row)
        for key in ["speed_mps_mean", "speed_per_hr_mean", "power_per_hr_mean", "hr_residual_bpm_mean"]:
            if display.get(key) is not None:
                display[key] = format_metric_value(key, display[key])
        window_display.append(display)
    window_table = markdown_table(window_display, list(window_rows[0].keys()))
    change_display = []
    for row in change_rows:
        display = {"metric": row["metric"]}
        display["early_value"] = format_metric_value(str(row["metric"]), row["early_value"])
        display["late_value"] = format_metric_value(str(row["metric"]), row["late_value"])
        if row["change_units"] == "delta_bpm":
            display["change_value"] = f"{float(row['change_value']):.2f}"
            display["change_units"] = "bpm_shift"
        else:
            display["change_value"] = f"{float(row['pct_change']):.2f}"
            display["change_units"] = "pct_change"
        change_display.append(display)
    change_table = markdown_table(change_display, ["metric", "early_value", "late_value", "change_value", "change_units"])
    hr_model_table = markdown_table(hr_model_rows, list(hr_model_rows[0].keys()))
    monthly_display = []
    for row in monthly_rows:
        display = dict(row)
        for key in ["speed_mps_norm", "speed_per_hr", "power_per_hr", "hr_residual_bpm", "running_share_28d_activity_pct"]:
            if display.get(key) is not None:
                display[key] = format_metric_value(key, display[key])
        monthly_display.append(display)
    monthly_table = markdown_table(
        monthly_display,
        [
            "year_month",
            "qc_run_count",
            "speed_mps_norm",
            "avg_hr_bpm",
            "avg_power_w",
            "speed_per_hr",
            "power_per_hr",
            "hr_residual_bpm",
            "running_share_28d_activity_pct",
            "next_day_resting_hr",
            "next_day_sleep_score",
            "next_day_hrv",
        ],
    )
    corr_table = markdown_table(corr_rows, ["relationship", "correlation_r", "n_points"])
    food_summary_table = markdown_table(food_summary, list(food_summary[0].keys()))
    food_corr_table = markdown_table(food_corr_rows, ["relationship", "correlation_r", "n_points"])
    scope_table = markdown_table(scope_rows, ["category", "question", "answer"])
    return f"""# Study 000B Results

## Summary

This flagship companion study supports a mixed answer rather than a simple one:

- later specialization coincided with better speed-per-heart-rate efficiency
- background recovery markers moved in a favorable direction as running became a larger share of recent activity
- session heart rate did not simply normalize away
- a small positive late-window unexplained HR burden remained even after simple speed-plus-power adjustment

## Early-to-late adaptive efficiency window

{window_table}

The main pattern is visible immediately:

- speed increased
- heart rate also increased
- speed per heart-rate cost improved
- power per heart-rate cost improved only slightly
- running specialization rose sharply
- next-day resting heart rate improved and next-day sleep score improved

These two windows are matched by run count, not by elapsed time. The early window spans `93` calendar days, while the late window spans `43`, so specialization and rolling-load contrasts should be read as denser-late versus broader-early context rather than as perfectly matched periods.

## Early-to-late change table

{change_table}

This table shows why the study is not just a "higher HR" story:

- speed improved by more than heart rate did
- `speed_per_hr_mean` improved by double digits
- `power_per_hr_mean` improved only modestly
- `hr_residual_bpm_mean` shifted upward rather than disappearing

That is the core mixed signal: the system became more efficient in one sense, but not cost-free.

## Session HR residual model

{hr_model_table}

The HR model is intentionally simple and descriptive. It uses speed and power to estimate expected session HR. The unexplained HR pattern matters more than the coefficients themselves:

- early-window residual burden was slightly negative on average
- late-window residual burden was slightly positive on average

That means the late window does not look like a pure internal-burden reduction story even though output efficiency improved.

## Monthly efficiency and specialization summary

{monthly_table}

At the monthly level:

- running share of recent activity rose markedly into 2026
- speed-per-HR improved into the later months
- HR residual burden remained mixed rather than collapsing monotonically

This supports a specialization-with-mixed-cost reading better than a simple fitness-improvement reading.

## Specialization correlations

{corr_table}

These correlations are descriptive and likely time-confounded, but they still help position the system:

- higher running share tracked better speed-per-HR efficiency
- higher running share tracked lower next-day resting HR
- higher running share tracked higher available HRV in the late window
- specialization did not meaningfully reduce the simple run-level unexplained HR burden by itself

## Partial late-window food context

{food_summary_table}

{food_corr_table}

The food layer remains partial and late-window only. It is useful as confound context, but not strong enough to anchor the central claim of this study.

## What this study can and cannot answer

{scope_table}
"""


def build_discussion_markdown(window_rows: list[dict[str, object]], change_rows: list[dict[str, object]], corr_rows: list[dict[str, object]]) -> str:
    early = next(row for row in window_rows if row["window_label"] == "early_qc_window")
    late = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    speed_eff = next(row for row in change_rows if row["metric"] == "speed_per_hr_mean")
    power_eff = next(row for row in change_rows if row["metric"] == "power_per_hr_mean")
    next_rhr = next(row for row in change_rows if row["metric"] == "next_day_resting_hr_mean")
    next_sleep = next(row for row in change_rows if row["metric"] == "next_day_sleep_score_mean")
    runshare_speed = next(row for row in corr_rows if row["relationship"] == "running_share_vs_speed_per_hr")
    runshare_rhr = next(row for row in corr_rows if row["relationship"] == "running_share_vs_nextday_resting_hr")
    residual_shift = float(late["hr_residual_bpm_mean"]) - float(early["hr_residual_bpm_mean"])
    return f"""# Study 000B Discussion

## What this flagship companion study adds to Study 000A

`Study 000A` established that later adaptation occurred through selective flexibility under increasing specialization. `Study 000B` asks what that specialization coincided with in terms of internal burden.

The answer is not "cost disappeared" and not "everything got worse." The answer is mixed.

## Best-supported interpretation

The clearest result is that adaptive efficiency improved.

Across the first and last 30 QC-pass runs:

- speed-per-heart-rate efficiency improved by `{speed_eff['pct_change']:.2f}%`
- power-per-heart-rate efficiency improved by only `{power_eff['pct_change']:.2f}%`
- running share of 28-day activity rose from `{early['running_share_28d_activity_pct_mean']:.2f}%` to `{late['running_share_28d_activity_pct_mean']:.2f}%`

That makes it hard to say the later system was merely working harder in a crude sense. It was producing more speed per unit heart-rate cost.

## Why the study still finds persistent burden

The mixed part of the picture is equally important.

The simple speed-plus-power HR model shows that average unexplained HR burden shifted from `{early['hr_residual_bpm_mean']:.2f}` bpm to `{late['hr_residual_bpm_mean']:.2f}` bpm, an upward shift of `{residual_shift:.2f}` bpm. In plain language, late sessions still carried slightly more HR than the simple external-output model would predict.

So the late system looks more efficient, but not frictionless.

## Recovery moved in a more favorable direction

Some background recovery markers moved in a favorable direction alongside specialization:

- next-day resting heart rate changed `{next_rhr['pct_change']:.2f}%`
- next-day sleep score changed `{next_sleep['pct_change']:.2f}%`
- late-window HRV became available, but early-window missingness prevents a symmetric early-late comparison

That means the later system may have been better supported in the background even while harder sessions still carried real acute burden.

## Specialization appears relevant, but not sufficient by itself

Specialization tracked better speed efficiency (`r = {runshare_speed['correlation_r']:.4f}`) and lower next-day resting HR (`r = {runshare_rhr['correlation_r']:.4f}`), but it did not by itself erase the simple run-level unexplained HR burden.

That is the most honest synthesis of the study:

`specialization coincided with better adaptive efficiency and more favorable background recovery markers, while leaving a persistent session-level unexplained HR burden`

## Why this matters for the broader program

This is valuable because it keeps the research program from collapsing into a false choice.

The system did not simply become more costly.
The system did not simply become more efficient.

It appears to have become more efficient in output terms while still carrying some unexplained burden in session terms.

That mixed pattern is likely one of the most important things this dataset can currently say.
"""


def build_limitations_markdown() -> str:
    return """# Study 000B Limitations

1. This is a single-subject within-system study and does not support population inference.
2. The HR residual model is intentionally simple and descriptive, not a validated physiology model.
3. The early and late comparison windows are matched by run count rather than elapsed time, spanning `93` versus `43` calendar days.
4. High-resolution running mechanics and session HR are concentrated in the later device-supported window.
5. Specialization correlations are descriptive and likely partly time-confounded.
6. HRV is available only late, which limits symmetric early-late comparison.
7. Sleep availability is slightly uneven across windows, and recovery markers are not equally mature in both periods.
8. The food layer is partial and late-window only, so it can only provide confound context rather than causal adjudication.
9. The study cannot prove exact physiological mechanism for the unexplained HR burden signal.
"""


def build_abstract_markdown(window_rows: list[dict[str, object]], change_rows: list[dict[str, object]]) -> str:
    early = next(row for row in window_rows if row["window_label"] == "early_qc_window")
    late = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    speed_eff = next(row for row in change_rows if row["metric"] == "speed_per_hr_mean")
    power_eff = next(row for row in change_rows if row["metric"] == "power_per_hr_mean")
    hr_resid = next(row for row in change_rows if row["metric"] == "hr_residual_bpm_mean")
    return f"""# Study 000B Abstract

## Title

Adaptive Efficiency And Internal Cost Under Specialization: A Flagship Companion Study Following Study 000A

## Abstract

This flagship companion study was conducted after `Study 000A` to determine whether later running specialization in the high-resolution window coincided with adaptive efficiency, persistent unexplained HR burden, or both. Using bundled QC-pass running rows and the unified daily timeline, run-level efficiency metrics, next-day recovery measures, and a simple descriptive heart-rate residual model were evaluated.

Across the first and last 30 QC-pass runs, running share of 28-day activity increased from `{early['running_share_28d_activity_pct_mean']:.2f}%` to `{late['running_share_28d_activity_pct_mean']:.2f}%`. Speed-per-heart-rate efficiency improved by `{speed_eff['pct_change']:.2f}%`, while power-per-heart-rate efficiency improved by only `{power_eff['pct_change']:.2f}%`. Next-day resting heart rate and next-day sleep score moved in a favorable direction. However, the simple speed-plus-power HR model showed that mean residual burden shifted upward rather than disappearing, from `{early['hr_residual_bpm_mean']:.2f}` bpm to `{late['hr_residual_bpm_mean']:.2f}` bpm.

These results support a mixed interpretation. Specialization coincided with better adaptive efficiency and more favorable background recovery markers, but it did not eliminate session-level unexplained HR burden. As a result, `Study 000B` strengthens the broader program by showing that later adaptation in this system was not reducible either to pure hidden cost or to pure efficiency gain.
"""


def build_plain_language_summary_markdown(change_rows: list[dict[str, object]]) -> str:
    speed_eff = next(row for row in change_rows if row["metric"] == "speed_per_hr_mean")
    power_eff = next(row for row in change_rows if row["metric"] == "power_per_hr_mean")
    hr_resid = next(row for row in change_rows if row["metric"] == "hr_residual_bpm_mean")
    return f"""# Study 000B Plain-Language Summary

## Why this exists

`Study 000A` showed that the system adapted under specialization. This study asks what that specialization coincided with in terms of internal burden.

## What this study found

The answer is mixed:

- speed per heart-rate cost improved by `{speed_eff['pct_change']:.2f}%`
- power per heart-rate cost improved by only `{power_eff['pct_change']:.2f}%`
- some next-day recovery markers moved in a favorable direction
- but mean unexplained session HR burden still shifted upward by `{hr_resid['change_value']:.2f}` bpm after simple speed-plus-power adjustment

## What that means

The later system looks more efficient, but not free of burden. It appears to have gained output efficiency while still carrying some unexplained session-level HR burden.

## Why that matters

This is valuable because it prevents the research program from oversimplifying the adaptation pattern. The system did not just get fitter, and it did not just get more strained. It appears to have become more efficient in some ways while preserving some burden in others.
"""


def build_submission_guidance_markdown() -> str:
    return """# Study 000B Submission Guidance

## Study type

Primary flagship companion study following the flagship longitudinal synthesis.

## Best use

- as the physiology-facing flagship companion to `Study 000A`
- as a paper or appendix showing that specialization coincided with better efficiency without fully removing unexplained session HR burden
- as a bridge between biomechanical adaptation and recovery-oriented follow-up work

## Safe claims

- specialization coincided with better speed-per-HR efficiency
- some background recovery markers moved in a favorable direction as specialization increased
- session-level HR burden did not disappear entirely
- the internal-cost picture is mixed rather than one-directional

## Avoid

- exact physiology mechanism claims
- clinical overreach
- strong fueling causation claims from the partial overlap layer
"""


def build_appendix_markdown(source_manifest: list[dict[str, object]], scope_rows: list[dict[str, object]]) -> str:
    source_table = markdown_table(source_manifest, ["table_name", "relative_path", "bytes", "sha256"])
    scope_table = markdown_table(scope_rows, ["category", "question", "answer"])
    return f"""# Appendix A: Packaged Source Tables And Scope

## Scope

{scope_table}

## Packaged source tables

{source_table}
"""


def build_audit_markdown(
    source_manifest: list[dict[str, object]],
    window_rows: list[dict[str, object]],
    change_rows: list[dict[str, object]],
    hr_model_rows: list[dict[str, object]],
) -> str:
    early = next(row for row in window_rows if row["window_label"] == "early_qc_window")
    late = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    speed_eff = next(row for row in change_rows if row["metric"] == "speed_per_hr_mean")
    hr_resid = next(row for row in change_rows if row["metric"] == "hr_residual_bpm_mean")
    model_row = hr_model_rows[0]
    return f"""# Study 000B Audit

## Structural audit

- Pass: the package is self-contained and reads only from bundled `source_tables`.
- Pass: the source manifest covers `{len(source_manifest)}` packaged input tables.
- Pass: the package includes manuscript, methods, results, discussion, limitations, figures, outputs, appendix, and audit report.

## Scientific audit

- Pass: the study clearly follows `Study 000A` while standing as a second primary program study focused on internal cost and efficiency.
- Pass: the central mixed finding is data-backed. `speed_per_hr_mean` changed `{speed_eff['pct_change']:.2f}%`.
- Pass: the persistent-burden claim is constrained and explicit. Mean unexplained HR burden shifted from `{early['hr_residual_bpm_mean']:.2f}` to `{late['hr_residual_bpm_mean']:.2f}` bpm.
- Pass: the descriptive HR model is documented rather than hidden, with coefficients recorded in the package.
- Pass: the late-window interpretation is limited by partial HRV, partial food coverage, and unmatched window span.

## Residual cautions

- The HR model is descriptive rather than physiologically validated.
- The early and late windows are not matched in elapsed time.
- Specialization correlations are likely partly time-confounded.
- Food context remains partial and cannot adjudicate the main finding.

## Audit verdict

`Study 000B` is structurally sound and scientifically coherent. It supports a valuable secondary conclusion for the program: specialization coincided with better adaptive efficiency and more favorable recovery markers, but did not fully remove session-level unexplained HR burden. That mixed pattern is well supported and appropriately bounded.
"""


def build_manuscript_markdown() -> str:
    parts = [
        (ROOT / "reports" / "STUDY000B_ABSTRACT.md").read_text(encoding="utf-8").strip(),
        (ROOT / "reports" / "STUDY000B_METHODS.md").read_text(encoding="utf-8").strip(),
        (ROOT / "reports" / "STUDY000B_RESULTS.md").read_text(encoding="utf-8").strip(),
        (ROOT / "reports" / "STUDY000B_DISCUSSION.md").read_text(encoding="utf-8").strip(),
        (ROOT / "reports" / "STUDY000B_LIMITATIONS.md").read_text(encoding="utf-8").strip(),
    ]
    return "\n\n".join(parts) + "\n"


def build_readme_markdown(project_state: dict[str, object]) -> str:
    findings = "\n".join(f"- {item}" for item in project_state["headline_findings"])
    return f"""# STUDY 000B

This folder contains the flagship companion study after `Study 000A`.

## Why this study exists

The flagship study established specialization and selective adaptation. This package asks whether the later specialized system became more efficient, more costly, or both.

## Headline findings

{findings}

## Main outputs

- `manuscript/STUDY000B_MANUSCRIPT.md`
- `reports/STUDY000B_ABSTRACT.md`
- `reports/STUDY000B_METHODS.md`
- `reports/STUDY000B_RESULTS.md`
- `reports/STUDY000B_DISCUSSION.md`
- `reports/STUDY000B_LIMITATIONS.md`
- `reports/STUDY000B_PLAIN_LANGUAGE_SUMMARY.md`
- `reports/STUDY000B_SUBMISSION_GUIDANCE.md`
- `reports/STUDY000B_AUDIT.md`
- `appendices/APPENDIX_A_SCOPE_AND_SOURCES.md`
- `figures/`
- `outputs/`
- `source_tables/`
"""


def build_project_state(change_rows: list[dict[str, object]], window_rows: list[dict[str, object]]) -> dict[str, object]:
    speed_eff = next(row for row in change_rows if row["metric"] == "speed_per_hr_mean")
    power_eff = next(row for row in change_rows if row["metric"] == "power_per_hr_mean")
    hr_resid = next(row for row in change_rows if row["metric"] == "hr_residual_bpm_mean")
    next_rhr = next(row for row in change_rows if row["metric"] == "next_day_resting_hr_mean")
    early = next(row for row in window_rows if row["window_label"] == "early_qc_window")
    late = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    return {
        "study_id": "STUDY-000B",
        "title": "Adaptive Efficiency And Internal Cost Under Specialization",
        "created": datetime.now().strftime("%Y-%m-%d %H:%M:%S"),
        "package_root": ".",
        "source_root": "source_tables",
        "source_manifest_file": "manifest/source_table_manifest.csv",
        "status": "flagship_companion_study_complete",
        "headline_findings": [
            f"Running share of 28-day activity increased from {early['running_share_28d_activity_pct_mean']:.2f}% to {late['running_share_28d_activity_pct_mean']:.2f}%.",
            f"Speed-per-heart-rate efficiency improved by {speed_eff['pct_change']:.2f}%.",
            f"Power-per-heart-rate efficiency improved by {power_eff['pct_change']:.2f}%, much less than speed-per-HR.",
            f"Next-day resting heart rate improved by {next_rhr['pct_change']:.2f}%.",
            f"Mean unexplained HR burden shifted upward by {hr_resid['change_value']:.2f} bpm, indicating persistent session-level burden.",
        ],
    }


def main() -> None:
    run_rows = [row for row in read_csv(RUNS_PATH) if row["run_dynamics_qc_pass"] in {"1", "True", "true"}]
    run_rows.sort(key=lambda row: row["calendar_date"])
    timeline_rows = read_csv(TIMELINE_PATH)
    timeline_by_date = {row["calendar_date"]: row for row in timeline_rows}
    food_rows = read_csv(FOOD_PATH)

    source_manifest = build_source_manifest([RUNS_PATH, TIMELINE_PATH, FOOD_PATH])
    run_level_rows = build_run_level_dataset(run_rows, timeline_by_date)
    window_rows = build_window_summary(run_level_rows)
    change_rows = build_window_change_rows(window_rows)
    monthly_rows = build_monthly_summary(run_level_rows)
    corr_rows = build_specialization_correlations(run_level_rows)
    hr_model_rows = build_hr_model_summary(run_level_rows)
    food_summary, food_corr_rows = build_food_context(food_rows)
    scope_rows = build_scope_claims()
    project_state = build_project_state(change_rows, window_rows)

    write_csv(ROOT / "manifest" / "source_table_manifest.csv", source_manifest, ["table_name", "relative_path", "bytes", "sha256"])
    write_csv(ROOT / "outputs" / "study000b_run_level_dataset.csv", run_level_rows, list(run_level_rows[0].keys()))
    write_csv(ROOT / "outputs" / "study000b_window_summary.csv", window_rows, list(window_rows[0].keys()))
    write_csv(ROOT / "outputs" / "study000b_window_change.csv", change_rows, ["metric", "early_value", "late_value", "change_value", "change_units", "pct_change"])
    write_csv(ROOT / "outputs" / "study000b_monthly_summary.csv", monthly_rows, list(monthly_rows[0].keys()))
    write_csv(ROOT / "outputs" / "study000b_specialization_correlations.csv", corr_rows, ["relationship", "correlation_r", "n_points"])
    write_csv(ROOT / "outputs" / "study000b_hr_model_summary.csv", hr_model_rows, list(hr_model_rows[0].keys()))
    write_csv(ROOT / "outputs" / "study000b_food_context_summary.csv", food_summary, list(food_summary[0].keys()))
    write_csv(ROOT / "outputs" / "study000b_food_context_correlations.csv", food_corr_rows, ["relationship", "correlation_r", "n_points"])
    write_csv(ROOT / "outputs" / "study000b_scope_claims.csv", scope_rows, ["category", "question", "answer"])

    build_figure_window_pct_change(change_rows, ROOT / "figures" / "figure01_output_cost_change.svg")
    build_figure_monthly_efficiency(monthly_rows, ROOT / "figures" / "figure02_monthly_specialization_efficiency.svg")
    build_figure_hr_residual(monthly_rows, ROOT / "figures" / "figure03_monthly_hr_residual.svg")
    build_figure_recovery(window_rows, ROOT / "figures" / "figure04_recovery_shift.svg")

    write_text(ROOT / "reports" / "STUDY000B_METHODS.md", build_methods_markdown(source_manifest))
    write_text(ROOT / "reports" / "STUDY000B_RESULTS.md", build_results_markdown(window_rows, change_rows, hr_model_rows, monthly_rows, corr_rows, food_summary, food_corr_rows, scope_rows))
    write_text(ROOT / "reports" / "STUDY000B_DISCUSSION.md", build_discussion_markdown(window_rows, change_rows, corr_rows))
    write_text(ROOT / "reports" / "STUDY000B_LIMITATIONS.md", build_limitations_markdown())
    write_text(ROOT / "reports" / "STUDY000B_ABSTRACT.md", build_abstract_markdown(window_rows, change_rows))
    write_text(ROOT / "reports" / "STUDY000B_PLAIN_LANGUAGE_SUMMARY.md", build_plain_language_summary_markdown(change_rows))
    write_text(ROOT / "reports" / "STUDY000B_SUBMISSION_GUIDANCE.md", build_submission_guidance_markdown())
    write_text(ROOT / "reports" / "STUDY000B_AUDIT.md", build_audit_markdown(source_manifest, window_rows, change_rows, hr_model_rows))
    write_text(ROOT / "appendices" / "APPENDIX_A_SCOPE_AND_SOURCES.md", build_appendix_markdown(source_manifest, scope_rows))
    write_text(ROOT / "manuscript" / "STUDY000B_MANUSCRIPT.md", build_manuscript_markdown())
    write_text(ROOT / "README.md", build_readme_markdown(project_state))
    write_text(ROOT / "manifest" / "study000b_project_state.json", json.dumps(project_state, indent=2))


if __name__ == "__main__":
    main()
