from __future__ import annotations

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


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

TIMELINE_PATH = SOURCE / "v7_total_fitness_profile_timeline_v3.csv"
FOOD_PATH = SOURCE / "v7_food_pdf_garmin_overlap_v1.csv"
YEARLY_PATH = SOURCE / "v7_pillar_total_adaptation_ecology_yearly_v1.csv"
MONTHLY_PATH = SOURCE / "v7_pillar_total_adaptation_ecology_monthly_v1.csv"
RUNS_PATH = SOURCE / "v7_pillar_conserved_mechanics_runs_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 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 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 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 build_yearly_overview(yearly_rows: list[dict[str, str]]) -> list[dict[str, object]]:
    rows: list[dict[str, object]] = []
    for row in yearly_rows:
        year = int(row["year"])
        if year < 2020:
            continue
        rows.append(
            {
                "year": year,
                "days_present": int(row["days_present"]),
                "training_days": int(row["training_days"]),
                "running_count": int(row["running_count"]),
                "indoor_conditioning_count": int(row["indoor_conditioning_count"]),
                "walking_count": int(row["walking_count"]),
                "strength_count": int(row["strength_count"]),
                "total_structured_training_hours": round(to_float(row["total_structured_training_hours"]) or 0.0, 3),
                "running_duration_hours": round(to_float(row["running_duration_hours"]) or 0.0, 3),
                "indoor_conditioning_duration_hours": round(to_float(row["indoor_conditioning_duration_hours"]) or 0.0, 3),
                "walking_duration_hours": round(to_float(row["walking_duration_hours"]) or 0.0, 3),
                "strength_duration_hours": round(to_float(row["strength_duration_hours"]) or 0.0, 3),
                "running_distance_miles": round(to_float(row["running_distance_miles"]) or 0.0, 3),
                "walking_distance_miles": round(to_float(row["walking_distance_miles"]) or 0.0, 3),
                "locomotion_distance_miles": round(to_float(row["locomotion_distance_miles"]) or 0.0, 3),
                "running_share_of_structured_hours": round((to_float(row["running_share_of_structured_hours"]) or 0.0) * 100.0, 2),
                "indoor_share_of_structured_hours": round((to_float(row["indoor_share_of_structured_hours"]) or 0.0) * 100.0, 2),
                "strength_share_of_structured_hours": round((to_float(row["strength_share_of_structured_hours"]) or 0.0) * 100.0, 2),
                "high_resolution_run_dynamics_sessions": int(row["high_resolution_run_dynamics_sessions"]),
                "qc_pass_run_dynamics_sessions": int(row["qc_pass_run_dynamics_sessions"]),
            }
        )
    return rows


def build_yearly_recovery_summary(timeline_rows: list[dict[str, str]]) -> list[dict[str, object]]:
    by_year: dict[int, list[dict[str, str]]] = defaultdict(list)
    for row in timeline_rows:
        if row["calendar_date"] < "2020-04-25":
            continue
        year = int(row["calendar_date"][:4])
        by_year[year].append(row)

    summaries = []
    metrics = [
        ("resting_heart_rate", "resting_heart_rate_mean"),
        ("sleep_score", "sleep_score_mean"),
        ("health_HRV", "health_hrv_mean"),
        ("running_hours_28d", "running_hours_28d_mean"),
        ("total_activity_hours_28d_rebuilt", "total_activity_hours_28d_mean"),
    ]
    for year in sorted(by_year):
        summary: dict[str, object] = {"year": year}
        rows = by_year[year]
        for source_key, target_key in metrics:
            values = [to_float(row[source_key]) for row in rows]
            values = [value for value in values if value is not None]
            summary[f"{target_key}_n"] = len(values)
            summary[target_key] = round(mean(values), 3) if values else None
        summaries.append(summary)
    return summaries


def build_phase_model(
    monthly_rows: list[dict[str, str]],
    timeline_rows: list[dict[str, str]],
    run_rows: list[dict[str, str]],
) -> list[dict[str, object]]:
    phases = [
        ("P1", "Running Launch", "2020-05", "2020-11"),
        ("P2", "Mixed Reconstruction", "2021-05", "2023-03"),
        ("P3", "Running Re-expansion", "2023-04", "2024-11"),
        ("P4", "Indoor Pivot And High-Resolution Adaptation", "2025-01", "2025-12"),
        ("P5", "Running Reconsolidation", "2026-01", "2026-06"),
    ]
    rows = []
    for phase_id, label, start_ym, end_ym in phases:
        month_subset = [
            row for row in monthly_rows
            if start_ym <= row["year_month"] <= end_ym and row["year_month"] != "2006-01"
        ]
        timeline_subset = [
            row for row in timeline_rows
            if start_ym <= row["calendar_date"][:7] <= end_ym
        ]
        run_subset = [
            row for row in run_rows
            if start_ym <= row["year_month"] <= end_ym
        ]
        qc_subset = [row for row in run_subset if row["run_dynamics_qc_pass"] in {"1", "True", "true"}]

        running_hours = sum(to_float(row["running_duration_hours"]) or 0.0 for row in month_subset)
        indoor_hours = sum(to_float(row["indoor_conditioning_duration_hours"]) or 0.0 for row in month_subset)
        walking_hours = sum(to_float(row["walking_duration_hours"]) or 0.0 for row in month_subset)
        strength_hours = sum(to_float(row["strength_duration_hours"]) or 0.0 for row in month_subset)
        modality_hours = {
            "running": running_hours,
            "indoor": indoor_hours,
            "walking": walking_hours,
            "strength": strength_hours,
        }
        dominant_modality = max(modality_hours, key=modality_hours.get) if modality_hours else "none"
        treadmill_runs = sum(1 for row in run_subset if row["rebuild_modality"] == "running_treadmill")
        outdoor_runs = sum(1 for row in run_subset if row["rebuild_modality"] == "running_outdoor")

        running_28d_vals = [to_float(row["running_hours_28d"]) for row in timeline_subset if to_float(row["running_hours_28d"]) is not None]
        total_28d_vals = [to_float(row["total_activity_hours_28d_rebuilt"]) for row in timeline_subset if to_float(row["total_activity_hours_28d_rebuilt"]) is not None]
        resting_vals = [to_float(row["resting_heart_rate"]) for row in timeline_subset if to_float(row["resting_heart_rate"]) is not None]

        rows.append(
            {
                "phase_id": phase_id,
                "phase_label": label,
                "start_year_month": start_ym,
                "end_year_month": end_ym,
                "months_in_phase": len(month_subset),
                "dominant_modality": dominant_modality,
                "mean_structured_training_hours_per_month": round(mean([to_float(row["total_structured_training_hours"]) or 0.0 for row in month_subset]), 3),
                "mean_running_share_pct": round(mean([(to_float(row["running_share_of_structured_hours"]) or 0.0) * 100.0 for row in month_subset]), 3),
                "mean_indoor_share_pct": round(mean([(to_float(row["indoor_share_of_structured_hours"]) or 0.0) * 100.0 for row in month_subset]), 3),
                "running_rows": len(run_subset),
                "qc_highres_runs": len(qc_subset),
                "treadmill_run_share_pct": round((treadmill_runs / len(run_subset)) * 100.0, 3) if run_subset else None,
                "outdoor_run_share_pct": round((outdoor_runs / len(run_subset)) * 100.0, 3) if run_subset else None,
                "mean_running_hours_28d": round(mean(running_28d_vals), 3) if running_28d_vals else None,
                "mean_total_activity_hours_28d": round(mean(total_28d_vals), 3) if total_28d_vals else None,
                "mean_resting_hr": round(mean(resting_vals), 3) if resting_vals else None,
            }
        )
    return rows


def build_window_context(run_rows: list[dict[str, str]], timeline_rows: list[dict[str, str]]) -> list[dict[str, object]]:
    qc_rows = [row for row in run_rows if row["run_dynamics_qc_pass"] in {"1", "True", "true"}]
    windows = [
        ("early_qc_window", qc_rows[:30]),
        ("late_qc_window", qc_rows[-30:]),
    ]
    results = []
    for label, subset in windows:
        start_date = subset[0]["calendar_date"]
        end_date = subset[-1]["calendar_date"]
        timeline_subset = [
            row for row in timeline_rows
            if start_date <= row["calendar_date"] <= end_date
        ]
        treadmill_runs = sum(1 for row in subset if row["rebuild_modality"] == "running_treadmill")
        outdoor_runs = sum(1 for row in subset if row["rebuild_modality"] == "running_outdoor")
        running_28d_vals = [to_float(row["running_hours_28d"]) for row in timeline_subset if to_float(row["running_hours_28d"]) is not None]
        total_28d_vals = [to_float(row["total_activity_hours_28d_rebuilt"]) for row in timeline_subset if to_float(row["total_activity_hours_28d_rebuilt"]) is not None]
        resting_vals = [to_float(row["resting_heart_rate"]) for row in timeline_subset if to_float(row["resting_heart_rate"]) is not None]
        sleep_vals = [to_float(row["sleep_score"]) for row in timeline_subset if to_float(row["sleep_score"]) is not None]
        hrv_vals = [to_float(row["health_HRV"]) for row in timeline_subset if to_float(row["health_HRV"]) is not None]
        running_mean = mean(running_28d_vals) if running_28d_vals else None
        total_mean = mean(total_28d_vals) if total_28d_vals else None
        results.append(
            {
                "window_label": label,
                "start_date": start_date,
                "end_date": end_date,
                "qc_run_count": len(subset),
                "treadmill_runs": treadmill_runs,
                "outdoor_runs": outdoor_runs,
                "treadmill_share_pct": round((treadmill_runs / len(subset)) * 100.0, 3),
                "mean_running_hours_28d": round(running_mean, 3) if running_mean is not None else None,
                "mean_total_activity_hours_28d": round(total_mean, 3) if total_mean is not None else None,
                "running_share_of_28d_activity_pct": round((running_mean / total_mean) * 100.0, 3) if running_mean is not None and total_mean not in {None, 0} else None,
                "mean_resting_hr": round(mean(resting_vals), 3) if resting_vals else None,
                "mean_sleep_score": round(mean(sleep_vals), 3) if sleep_vals else None,
                "mean_health_hrv": round(mean(hrv_vals), 3) if hrv_vals else None,
            }
        )
    return results


def build_highres_summary(run_rows: list[dict[str, str]]) -> tuple[list[dict[str, object]], list[dict[str, object]], list[dict[str, object]]]:
    qc_rows = [row for row in run_rows if row["run_dynamics_qc_pass"] in {"1", "True", "true"}]
    early = qc_rows[:30]
    late = qc_rows[-30:]
    metrics = [
        ("pace_min_per_mile_norm", "pace_min_per_mile"),
        ("avg_hr_bpm", "avg_hr_bpm"),
        ("avg_power_w", "avg_power_w"),
        ("cadence_spm_est", "cadence_spm"),
        ("stride_length_m_est", "stride_length_m"),
        ("vertical_oscillation_cm", "vertical_oscillation_cm"),
        ("ground_contact_time_ms", "ground_contact_time_ms"),
        ("vertical_ratio_pct", "vertical_ratio_pct"),
    ]

    early_late_rows = []
    for source_key, label in metrics:
        early_values = [to_float(row[source_key]) for row in early]
        late_values = [to_float(row[source_key]) for row in late]
        early_values = [value for value in early_values if value is not None]
        late_values = [value for value in late_values if value is not None]
        early_mean = mean(early_values)
        late_mean = mean(late_values)
        early_late_rows.append(
            {
                "metric": label,
                "early_mean": round(early_mean, 5),
                "late_mean": round(late_mean, 5),
                "pct_change": round(pct_change(early_mean, late_mean), 5),
            }
        )

    correlation_specs = [
        ("pace_min_per_mile_norm", "cadence_spm_est", "pace_vs_cadence"),
        ("pace_min_per_mile_norm", "stride_length_m_est", "pace_vs_stride"),
        ("pace_min_per_mile_norm", "avg_power_w", "pace_vs_power"),
        ("pace_min_per_mile_norm", "ground_contact_time_ms", "pace_vs_gct"),
        ("pace_min_per_mile_norm", "vertical_oscillation_cm", "pace_vs_vertical_oscillation"),
        ("pace_min_per_mile_norm", "vertical_ratio_pct", "pace_vs_vertical_ratio"),
        ("pace_min_per_mile_norm", "avg_hr_bpm", "pace_vs_avg_hr"),
    ]
    correlation_rows = []
    for x_key, y_key, label in correlation_specs:
        r_value, n_points = correlation(
            [to_float(row[x_key]) for row in qc_rows],
            [to_float(row[y_key]) for row in qc_rows],
        )
        correlation_rows.append(
            {
                "relationship": label,
                "correlation_r": round(r_value, 5) if r_value is not None else None,
                "n_points": n_points,
            }
        )

    monthly_groups: dict[str, list[dict[str, str]]] = defaultdict(list)
    for row in qc_rows:
        monthly_groups[row["year_month"]].append(row)

    monthly_rows = []
    for year_month in sorted(monthly_groups):
        rows = monthly_groups[year_month]
        if len(rows) < 3:
            continue
        monthly_rows.append(
            {
                "year_month": year_month,
                "qc_run_count": len(rows),
                "pace_min_per_mile": round(mean([to_float(row["pace_min_per_mile_norm"]) for row in rows if to_float(row["pace_min_per_mile_norm"]) is not None]), 5),
                "cadence_spm": round(mean([to_float(row["cadence_spm_est"]) for row in rows if to_float(row["cadence_spm_est"]) is not None]), 5),
                "stride_length_m": round(mean([to_float(row["stride_length_m_est"]) for row in rows if to_float(row["stride_length_m_est"]) is not None]), 5),
                "avg_power_w": round(mean([to_float(row["avg_power_w"]) for row in rows if to_float(row["avg_power_w"]) is not None]), 5),
                "ground_contact_time_ms": round(mean([to_float(row["ground_contact_time_ms"]) for row in rows if to_float(row["ground_contact_time_ms"]) is not None]), 5),
                "vertical_oscillation_cm": round(mean([to_float(row["vertical_oscillation_cm"]) for row in rows if to_float(row["vertical_oscillation_cm"]) is not None]), 5),
                "vertical_ratio_pct": round(mean([to_float(row["vertical_ratio_pct"]) for row in rows if to_float(row["vertical_ratio_pct"]) is not None]), 5),
            }
        )
    return early_late_rows, correlation_rows, monthly_rows


def build_speed_decomposition(run_rows: list[dict[str, str]]) -> list[dict[str, object]]:
    qc_rows = [row for row in run_rows if row["run_dynamics_qc_pass"] in {"1", "True", "true"}]
    early = qc_rows[:30]
    late = qc_rows[-30:]

    def metric_mean(rows: list[dict[str, str]], key: str) -> float:
        values = [to_float(row[key]) for row in rows]
        return mean([value for value in values if value is not None])

    early_cadence = metric_mean(early, "cadence_spm_est")
    late_cadence = metric_mean(late, "cadence_spm_est")
    early_stride = metric_mean(early, "stride_length_m_est")
    late_stride = metric_mean(late, "stride_length_m_est")
    early_speed = metric_mean(early, "speed_mps_norm")
    late_speed = metric_mean(late, "speed_mps_norm")

    cadence_log_change = math.log(late_cadence / early_cadence)
    stride_log_change = math.log(late_stride / early_stride)
    total_log_change = cadence_log_change + stride_log_change

    return [
        {
            "window_pair": "first30_vs_last30_qc_runs",
            "early_speed_mps": round(early_speed, 5),
            "late_speed_mps": round(late_speed, 5),
            "speed_gain_pct": round(((late_speed / early_speed) - 1.0) * 100.0, 5),
            "early_cadence_spm": round(early_cadence, 5),
            "late_cadence_spm": round(late_cadence, 5),
            "cadence_pct_change": round(pct_change(early_cadence, late_cadence), 5),
            "early_stride_length_m": round(early_stride, 5),
            "late_stride_length_m": round(late_stride, 5),
            "stride_pct_change": round(pct_change(early_stride, late_stride), 5),
            "cadence_log_share_pct": round((cadence_log_change / total_log_change) * 100.0, 5),
            "stride_log_share_pct": round((stride_log_change / total_log_change) * 100.0, 5),
        }
    ]


def build_monthly_variability(monthly_highres_rows: list[dict[str, object]]) -> list[dict[str, object]]:
    metrics = [
        ("pace_min_per_mile", "adaptation"),
        ("cadence_spm", "adaptation"),
        ("stride_length_m", "adaptation"),
        ("avg_power_w", "adaptation"),
        ("ground_contact_time_ms", "adaptation"),
        ("vertical_oscillation_cm", "conserved_candidate"),
        ("vertical_ratio_pct", "conserved_candidate"),
    ]
    rows = []
    for metric, category in metrics:
        values = [float(row[metric]) for row in monthly_highres_rows if row.get(metric) is not None]
        if not values:
            continue
        metric_mean = mean(values)
        monthly_cv_pct = (sample_std(values) / metric_mean) * 100.0 if metric_mean != 0 else math.nan
        rows.append(
            {
                "metric": metric,
                "metric_category": category,
                "monthly_n": len(values),
                "monthly_mean": round(metric_mean, 5),
                "monthly_cv_pct": round(monthly_cv_pct, 5),
                "monthly_first_last_pct_change": round(pct_change(values[0], values[-1]), 5),
            }
        )

    sorted_rows = sorted(rows, key=lambda row: abs(float(row["monthly_cv_pct"])))
    for idx, row in enumerate(sorted_rows, start=1):
        row["conservation_rank_low_cv"] = idx
    return sorted_rows


def build_window_change(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 = [
        ("treadmill_share_pct", "treadmill_share_pct"),
        ("running_share_of_28d_activity_pct", "running_share_of_28d_activity_pct"),
        ("mean_running_hours_28d", "running_hours_28d"),
        ("mean_total_activity_hours_28d", "total_activity_hours_28d"),
        ("mean_resting_hr", "resting_hr"),
        ("mean_sleep_score", "sleep_score"),
    ]
    rows = []
    for source_key, label in metrics:
        early_value = early[source_key]
        late_value = late[source_key]
        if early_value is None or late_value is None:
            continue
        rows.append(
            {
                "metric": label,
                "early_value": round(float(early_value), 5),
                "late_value": round(float(late_value), 5),
                "pct_change": round(pct_change(float(early_value), float(late_value)), 5),
            }
        )
    return rows


def build_timeline_correlations(timeline_rows: list[dict[str, str]]) -> list[dict[str, object]]:
    filtered = [row for row in timeline_rows if row["calendar_date"] >= "2020-04-25"]
    specs = [
        ("running_hours_28d", "resting_heart_rate", "running_28d_vs_resting_hr"),
        ("running_hours_28d", "sleep_score", "running_28d_vs_sleep_score"),
        ("running_hours_28d", "health_HRV", "running_28d_vs_hrv"),
        ("total_activity_hours_28d_rebuilt", "resting_heart_rate", "total_activity_28d_vs_resting_hr"),
    ]
    rows = []
    for x_key, y_key, label in specs:
        r_value, n_points = correlation(
            [to_float(row[x_key]) for row in filtered],
            [to_float(row[y_key]) for row in filtered],
        )
        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_food_summary(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_row = {
        "days_in_overlap": len(food_rows),
        "days_with_calories": sum(1 for row in food_rows if to_float(row["food_total_calories"]) is not None),
        "days_with_net_carbs": sum(1 for row in food_rows if to_float(row["food_total_net_carbs_g"]) is not None),
        "days_with_food_item_rows": len(food_item_days),
        "days_with_sleep_score": sum(1 for row in food_rows if to_float(row["sleep_score"]) is not None),
        "days_with_hrv": sum(1 for row in food_rows if to_float(row["health_HRV"]) is not None),
        "mean_food_total_calories_all_overlap_days": round(mean([to_float(row["food_total_calories"]) for row in food_rows if to_float(row["food_total_calories"]) is not None]), 3),
        "mean_food_total_net_carbs_g_all_overlap_days": 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]), 3),
        "mean_food_total_calories_food_item_days_only": round(mean([to_float(row["food_total_calories"]) for row in food_item_days if to_float(row["food_total_calories"]) is not None]), 3),
        "mean_food_total_net_carbs_g_food_item_days_only": 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]), 3),
        "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]), 3),
        "mean_total_activity_hours_7d": round(mean([to_float(row["total_activity_hours_7d_rebuilt"]) for row in food_rows if to_float(row["total_activity_hours_7d_rebuilt"]) is not None]), 3),
    }

    corr_specs = [
        ("food_total_calories", "running_hours_7d", "food_calories_vs_running_hours_7d"),
        ("food_total_net_carbs_g", "running_hours_7d", "net_carbs_vs_running_hours_7d"),
        ("food_total_calories", "resting_heart_rate", "food_calories_vs_resting_hr"),
        ("food_total_net_carbs_g", "sleep_score", "net_carbs_vs_sleep_score"),
    ]
    corr_rows = []
    for x_key, y_key, label in corr_specs:
        r_value, n_points = correlation(
            [to_float(row[x_key]) for row in food_rows],
            [to_float(row[y_key]) for row in food_rows],
        )
        corr_rows.append(
            {
                "relationship": label,
                "correlation_r": round(r_value, 5) if r_value is not None else None,
                "n_points": n_points,
            }
        )
    return [summary_row], corr_rows


def build_scope_claims() -> list[dict[str, object]]:
    return [
        {
            "category": "answered",
            "question": "Did workload, modality mix, and performance change across the six-year record?",
            "answer": "Yes. The account shows major year-to-year changes in structured training volume, modality share, running volume, and later high-resolution running performance.",
        },
        {
            "category": "answered",
            "question": "Did all measured running mechanics change in parallel with performance?",
            "answer": "No. Multiple variables changed substantially, but some mechanics remained comparatively more conserved than pace, cadence, power, and GCT.",
        },
        {
            "category": "supported_not_proven",
            "question": "What adaptation pathway does the dataset most strongly suggest?",
            "answer": "The strongest within-subject signal is turnover-led adaptation, with cadence tracking pace more strongly than stride length in the high-resolution window.",
        },
        {
            "category": "not_answered",
            "question": "How does this system compare directly to a non-compromised or normative comparison system?",
            "answer": "This study cannot answer that directly because it has no matched external comparator cohort or literature-normalized control framework embedded in the package.",
        },
        {
            "category": "not_answered",
            "question": "What exact biomechanical mechanism caused the adaptation pattern?",
            "answer": "This study cannot prove mechanism. It identifies a longitudinal pattern and the variables that changed together, but not the exact tissue-level or control-level cause.",
        },
    ]


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, dash: str | None = None) -> str:
    dash_attr = f' stroke-dasharray="{dash}"' if dash else ""
    return f'<line x1="{x1:.1f}" y1="{y1:.1f}" x2="{x2:.1f}" y2="{y2:.1f}" stroke="{stroke}" stroke-width="{width:.1f}"{dash_attr} />'


def svg_circle(cx: float, cy: float, r: float, fill: str) -> str:
    return f'<circle cx="{cx:.1f}" cy="{cy:.1f}" r="{r:.1f}" fill="{fill}" />'


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_yearly_ecology(yearly_rows: list[dict[str, object]], path: Path) -> None:
    width, height = 960, 520
    left, top, chart_w, chart_h = 80, 70, 760, 340
    years = [row["year"] for row in yearly_rows]
    totals = [row["total_structured_training_hours"] for row in yearly_rows]
    max_total = max(totals) if totals else 1.0
    colors = {
        "running_duration_hours": "#1768ac",
        "indoor_conditioning_duration_hours": "#f26419",
        "walking_duration_hours": "#86bbd8",
        "strength_duration_hours": "#33658a",
    }
    labels = {
        "running_duration_hours": "Running",
        "indoor_conditioning_duration_hours": "Indoor",
        "walking_duration_hours": "Walking",
        "strength_duration_hours": "Strength",
    }

    body = []
    body.append(svg_text(width / 2, 32, "Figure 1. Six-Year Structured Training Ecology", 22, "bold", "middle"))
    body.append(svg_text(width / 2, 52, "Stacked annual hours by primary modality", 13, "normal", "middle", "#555"))
    body.append(svg_line(left, top + chart_h, left + chart_w, top + chart_h, "#444", 1.5))
    body.append(svg_line(left, top, left, top + chart_h, "#444", 1.5))

    bar_w = chart_w / max(len(years), 1) * 0.62
    gap = chart_w / max(len(years), 1)
    keys = list(colors.keys())
    for idx, row in enumerate(yearly_rows):
        x = left + idx * gap + (gap - bar_w) / 2
        cumulative = 0.0
        for key in keys:
            value = float(row[key])
            h = 0 if max_total == 0 else (value / max_total) * chart_h
            y = top + chart_h - cumulative - h
            body.append(svg_rect(x, y, bar_w, h, colors[key]))
            cumulative += h
        body.append(svg_text(x + bar_w / 2, top + chart_h + 22, str(row["year"]), 11, "normal", "middle"))
        body.append(svg_text(x + bar_w / 2, top + chart_h - cumulative - 6, f'{row["total_structured_training_hours"]:.0f}', 10, "normal", "middle", "#333"))

    for tick in range(5):
        value = max_total * tick / 4
        y = top + chart_h - (chart_h * tick / 4)
        body.append(svg_line(left - 4, y, left + chart_w, y, "#d9d3c5", 1.0))
        body.append(svg_text(left - 10, y + 4, f"{value:.0f}", 10, "normal", "end", "#555"))

    legend_x = left + chart_w + 20
    legend_y = top + 20
    for idx, key in enumerate(keys):
        body.append(svg_rect(legend_x, legend_y + idx * 24, 16, 16, colors[key]))
        body.append(svg_text(legend_x + 24, legend_y + idx * 24 + 13, labels[key], 11))

    write_text(path, svg_canvas(width, height, "".join(body)))


def build_figure_adaptation_scatter(run_rows: list[dict[str, str]], path: Path) -> None:
    qc_rows = [row for row in run_rows if row["run_dynamics_qc_pass"] in {"1", "True", "true"}]
    width, height = 980, 460
    body = []
    body.append(svg_text(width / 2, 28, "Figure 2. Pace Tracking By Cadence And Stride Length", 21, "bold", "middle"))
    body.append(svg_text(width / 2, 48, "QC-pass high-resolution sessions only", 13, "normal", "middle", "#555"))

    def panel(x0: int, y0: int, w: int, h: int, x_key: str, y_key: str, title: str, point_color: str) -> str:
        xs = [to_float(r[x_key]) for r in qc_rows if to_float(r[x_key]) is not None and to_float(r[y_key]) is not None]
        ys = [to_float(r[y_key]) for r in qc_rows if to_float(r[x_key]) is not None and to_float(r[y_key]) is not None]
        if not xs or not ys:
            return ""
        min_x, max_x = min(xs), max(xs)
        min_y, max_y = min(ys), max(ys)
        pad_x = (max_x - min_x) * 0.05 if max_x > min_x else 1.0
        pad_y = (max_y - min_y) * 0.05 if max_y > min_y else 1.0
        min_x -= pad_x
        max_x += pad_x
        min_y -= pad_y
        max_y += pad_y

        parts = [svg_line(x0, y0 + h, x0 + w, y0 + h, "#444", 1.4), svg_line(x0, y0, x0, y0 + h, "#444", 1.4)]
        parts.append(svg_text(x0 + w / 2, y0 - 12, title, 14, "bold", "middle"))
        for x_val, y_val in zip(xs, ys):
            px = x0 + ((x_val - min_x) / (max_x - min_x)) * w
            py = y0 + h - ((y_val - min_y) / (max_y - min_y)) * h
            parts.append(svg_circle(px, py, 3.0, point_color))
        parts.append(svg_text(x0 + w / 2, y0 + h + 28, "Pace (min/mile)", 12, "normal", "middle"))
        parts.append(svg_text(x0 - 36, y0 + h / 2, y_key.replace("_", " "), 12, "normal", "middle"))
        return "".join(parts)

    body.append(panel(70, 90, 360, 280, "pace_min_per_mile_norm", "cadence_spm_est", "Pace vs Cadence", "#1768ac"))
    body.append(panel(540, 90, 360, 280, "pace_min_per_mile_norm", "stride_length_m_est", "Pace vs Stride Length", "#f26419"))
    write_text(path, svg_canvas(width, height, "".join(body)))


def build_figure_pct_change(change_rows: list[dict[str, object]], path: Path) -> None:
    width, height = 960, 520
    left, top, chart_w, chart_h = 90, 80, 760, 320
    focus = [
        "pace_min_per_mile",
        "cadence_spm",
        "stride_length_m",
        "avg_power_w",
        "ground_contact_time_ms",
        "vertical_oscillation_cm",
        "vertical_ratio_pct",
    ]
    selected = [row for row in change_rows if row["metric"] in focus]
    max_abs = max(abs(float(row["pct_change"])) for row in selected)
    body = []
    body.append(svg_text(width / 2, 32, "Figure 3. Early-to-Late Percent Change In High-Resolution Metrics", 21, "bold", "middle"))
    body.append(svg_text(width / 2, 52, "Adaptation variables changed more than conserved candidates", 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
    colors = {
        "pace_min_per_mile": "#f26419",
        "cadence_spm": "#1768ac",
        "stride_length_m": "#86bbd8",
        "avg_power_w": "#33658a",
        "ground_contact_time_ms": "#6d597a",
        "vertical_oscillation_cm": "#2a9d8f",
        "vertical_ratio_pct": "#8ab17d",
    }
    for idx, row in enumerate(selected):
        x = left + idx * gap + (gap - bar_w) / 2
        pct = float(row["pct_change"])
        h = 0 if max_abs == 0 else (abs(pct) / max_abs) * (chart_h / 2 - 18)
        y = mid_y - h if pct >= 0 else mid_y
        body.append(svg_rect(x, y, bar_w, h, colors[row["metric"]]))
        body.append(svg_text(x + bar_w / 2, top + chart_h + 24, row["metric"].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_window_context(window_rows: list[dict[str, object]], path: Path) -> None:
    width, height = 980, 520
    left, top, chart_w, chart_h = 100, 90, 760, 300
    metrics = [
        ("treadmill_share_pct", "Treadmill Share %", "#1768ac"),
        ("running_share_of_28d_activity_pct", "Running Share Of 28d Activity %", "#f26419"),
        ("mean_resting_hr", "Mean Resting HR", "#2a9d8f"),
    ]
    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")
    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. High-Resolution Window Context Shift", 21, "bold", "middle"))
    body.append(svg_text(width / 2, 50, "Specialization increased while treadmill dominance remained high", 13, "normal", "middle", "#555"))
    body.append(svg_line(left, top + chart_h, left + chart_w, top + chart_h, "#444", 1.3))
    body.append(svg_line(left, top, left, top + chart_h, "#444", 1.3))

    group_gap = chart_w / len(metrics)
    bar_w = group_gap * 0.22
    for idx, (key, label, color) in enumerate(metrics):
        x_center = left + idx * group_gap + group_gap / 2
        for offset, row, name, fill in [(-bar_w * 0.7, early, "Early", color), (bar_w * 0.7, late, "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 = 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(x_center, top + chart_h + 24, label, 10, "normal", "middle"))

    legend_x = left + chart_w - 120
    legend_y = top + 15
    body.append(svg_rect(legend_x, legend_y, 16, 16, "#1768ac"))
    body.append(svg_text(legend_x + 24, legend_y + 13, "Early", 11))
    body.append(svg_rect(legend_x, legend_y + 24, 16, 16, "#555f6d"))
    body.append(svg_text(legend_x + 24, legend_y + 37, "Late", 11))

    write_text(path, svg_canvas(width, height, "".join(body)))


def build_figure_speed_decomposition(speed_rows: list[dict[str, object]], path: Path) -> None:
    row = speed_rows[0]
    width, height = 920, 420
    left, top, chart_w, chart_h = 110, 90, 640, 180
    cadence_share = float(row["cadence_log_share_pct"])
    stride_share = float(row["stride_log_share_pct"])
    body = []
    body.append(svg_text(width / 2, 30, "Figure 5. Decomposition Of Late-Window Speed Gain", 21, "bold", "middle"))
    body.append(svg_text(width / 2, 50, "Cadence contributed more than stride length to the observed speed gain", 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))

    full_w = chart_w * 0.78
    bar_x = left + 30
    bar_y = top + 52
    bar_h = 54
    cadence_w = full_w * (cadence_share / 100.0)
    stride_w = full_w * (stride_share / 100.0)

    body.append(svg_rect(bar_x, bar_y, cadence_w, bar_h, "#1768ac"))
    body.append(svg_rect(bar_x + cadence_w, bar_y, stride_w, bar_h, "#f26419"))
    body.append(svg_text(bar_x + cadence_w / 2, bar_y + 32, f"Cadence {cadence_share:.1f}%", 14, "bold", "middle", "#ffffff"))
    body.append(svg_text(bar_x + cadence_w + stride_w / 2, bar_y + 32, f"Stride {stride_share:.1f}%", 14, "bold", "middle", "#ffffff"))

    body.append(svg_text(bar_x, bar_y + 92, f"Observed speed gain: {float(row['speed_gain_pct']):.1f}%", 13, "bold"))
    body.append(svg_text(bar_x, bar_y + 114, f"Cadence: {float(row['early_cadence_spm']):.1f} to {float(row['late_cadence_spm']):.1f} spm", 12))
    body.append(svg_text(bar_x, bar_y + 136, f"Stride length: {float(row['early_stride_length_m']):.3f} to {float(row['late_stride_length_m']):.3f} m", 12))
    body.append(svg_text(bar_x, bar_y + 158, "Contribution shares use log-change decomposition of cadence and stride length across the first and last 30 QC-pass runs.", 11, "normal", "start", "#555"))

    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 000A Methods

## Design

Single-subject longitudinal synthesis study built from the Stage 5 Garmin main package and designed to answer a broader six-year question than the earlier foundation studies.

## Core question

Across six years of self-tracked data, what changed, what remained comparatively anchored, and what adaptation pathway does the dataset support in the setting of altered biomechanics and a nonstandard movement system?

## Packaged source tables

{manifest_table}

## Analytic structure

1. Annual ecology was summarized from the rebuilt master to describe total structured training, modality mix, and locomotion context.
2. A descriptive phase model was constructed from sustained modality shifts and later running-context transitions.
3. High-resolution running mechanics were analyzed only inside the run-dynamics-supported window and restricted to QC-pass sessions for the main adaptation-pathway screen.
4. Early and late QC-pass windows were linked back to the daily master timeline to quantify specialization, treadmill dominance, and recovery context around the performance shift.
5. Recovery context was summarized from the unified timeline using resting heart rate, sleep score, HRV coverage, and rolling load variables where available.
6. Food overlap was treated as partial late-window context only, not as complete six-year fueling truth.
7. Cadence-stride speed decomposition was treated as descriptive rather than causal, because pace and speed are kinematically related to cadence and stride length.

## High-resolution mechanics window

The high-resolution mechanics layer covers the later device-supported running window and uses the already-built `run_dynamics_qc_pass` indicator from the master package.

The synthesis uses:

- annual ecology for the six-year trajectory
- a descriptive phase model for the major training eras
- QC-pass high-resolution runs for direct mechanics and adaptation-pathway analysis
- first 30 versus last 30 QC-pass runs for a simple early-to-late contrast
- timeline-linked window context to measure specialization and treadmill dominance around the mechanics shift
- correlation screens to identify which variables tracked pace most strongly

## Interpretation boundary

This study is not a clinical or physician-guided paper. It interprets a six-year longitudinal pattern in a person with known altered biomechanics, but it does not claim diagnosis, treatment effect, or tissue-level mechanism. It is strongest as a human-performance, biomechanics, and wearable-interpretation synthesis.
"""


def build_results_markdown(
    yearly_overview: list[dict[str, object]],
    phase_rows: list[dict[str, object]],
    window_rows: list[dict[str, object]],
    yearly_recovery: list[dict[str, object]],
    highres_changes: list[dict[str, object]],
    correlation_rows: list[dict[str, object]],
    speed_rows: list[dict[str, object]],
    variability_rows: list[dict[str, object]],
    window_change_rows: list[dict[str, object]],
    timeline_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:
    yearly_table = markdown_table(
        yearly_overview,
        [
            "year",
            "training_days",
            "running_count",
            "indoor_conditioning_count",
            "walking_count",
            "strength_count",
            "total_structured_training_hours",
            "running_duration_hours",
            "indoor_conditioning_duration_hours",
            "running_share_of_structured_hours",
            "indoor_share_of_structured_hours",
        ],
    )
    phase_table = markdown_table(
        phase_rows,
        [
            "phase_id",
            "phase_label",
            "start_year_month",
            "end_year_month",
            "dominant_modality",
            "mean_structured_training_hours_per_month",
            "mean_running_share_pct",
            "mean_indoor_share_pct",
            "treadmill_run_share_pct",
            "qc_highres_runs",
        ],
    )
    window_table = markdown_table(
        window_rows,
        [
            "window_label",
            "start_date",
            "end_date",
            "qc_run_count",
            "treadmill_share_pct",
            "running_share_of_28d_activity_pct",
            "mean_running_hours_28d",
            "mean_total_activity_hours_28d",
            "mean_resting_hr",
            "mean_sleep_score",
            "mean_health_hrv",
        ],
    )
    recovery_table = markdown_table(
        yearly_recovery,
        [
            "year",
            "resting_heart_rate_mean_n",
            "resting_heart_rate_mean",
            "sleep_score_mean_n",
            "sleep_score_mean",
            "health_hrv_mean_n",
            "health_hrv_mean",
            "running_hours_28d_mean",
            "total_activity_hours_28d_mean",
        ],
    )
    change_table = markdown_table(
        highres_changes,
        ["metric", "early_mean", "late_mean", "pct_change"],
    )
    corr_table = markdown_table(correlation_rows, ["relationship", "correlation_r", "n_points"])
    speed_table = markdown_table(speed_rows, list(speed_rows[0].keys()))
    variability_table = markdown_table(
        variability_rows,
        [
            "metric",
            "metric_category",
            "monthly_n",
            "monthly_mean",
            "monthly_cv_pct",
            "monthly_first_last_pct_change",
            "conservation_rank_low_cv",
        ],
    )
    window_change_table = markdown_table(window_change_rows, ["metric", "early_value", "late_value", "pct_change"])
    timeline_corr_table = markdown_table(timeline_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 000A Results

## Summary

This six-year synthesis supports a broader answer than the earlier foundation studies:

- the training ecosystem changed dramatically across years
- the adaptation pathway shifted by era, specialization, and running environment, not just by speed alone
- the later high-resolution running window shows turnover-led adaptation with modest stride expansion
- that later mechanics window is overwhelmingly treadmill-dominant
- some mechanics remained comparatively more conserved than pace, cadence, power, and ground contact time
- recovery context improved in some ways even while late-window running heart rate remained high

## Six-year training ecology

{yearly_table}

The largest whole-system change was not simply "more running." The dataset moved through distinct training eras:

- 2020 was overwhelmingly running-dominant
- 2021 through 2024 showed broader modality mixing and indoor-conditioning prominence
- 2025 became a high-volume indoor-conditioning year with meaningful running preserved
- 2026 shifted back toward a strongly running-dominant structured profile

That means the longitudinal adaptation story is ecological, not single-modality.

## Descriptive phase model

{phase_table}

The six-year record can be organized into a clearer phase model:

- a running-launch phase
- a mixed reconstruction phase
- a running re-expansion phase
- an indoor pivot and high-resolution adaptation phase
- a running reconsolidation phase

This phase model makes the adaptation story more concrete: the system did not just improve. It reorganized repeatedly.

## Six-year recovery context

{recovery_table}

Within the data that are available, later years show lower resting heart rate and better sleep-score coverage than the earlier years. HRV only becomes available late, so it supports later-window interpretation rather than six-year symmetric comparison.

## High-resolution early-to-late mechanics contrast

{change_table}

This is where the adaptation pathway becomes clearer. In the later device-supported window:

- pace improved substantially
- cadence increased strongly
- stride length increased, but less than cadence
- power rose
- ground contact time fell
- vertical oscillation and vertical ratio changed less than the main adaptation variables

The strongest reading is not that one metric stayed "perfectly fixed." It is that the mechanics changed unequally, with some features remaining more anchored than others.

## Cadence-stride decomposition of speed gain

{speed_table}

This table makes the turnover-led claim more concrete. Across the first and last 30 QC-pass runs:

- observed speed increased by more than one third
- cadence contributed the larger share of the cadence-stride speed gain
- stride length still contributed meaningfully, but less than cadence

That is stronger than saying cadence merely "correlated better." Descriptively, it shows that the later speed gain was assembled more through turnover than through stride expansion.

## Monthly variability and conservation ranking

{variability_table}

The conservation signal is not only about early-versus-late contrast. It also appears in month-to-month variability:

- `vertical_ratio_pct` had the lowest monthly coefficient of variation
- `vertical_oscillation_cm` was the lowest-variability direct mechanics candidate
- pace, power, and ground contact time were materially more volatile than the conserved candidates

That makes the "comparatively anchored" interpretation more defensible than a single before-and-after contrast alone.

## High-resolution window context

{window_table}

## High-resolution window change summary

{window_change_table}

This table sharpens the interpretation of the mechanics window:

- the QC-pass high-resolution dataset is treadmill-dominant overall
- the early 30-run window is already mostly treadmill-based
- the late 30-run window is entirely treadmill-based
- running-specific 28-day load increased while total 28-day activity fell
- running became a much larger share of the total recent activity profile

That means the later adaptation signal was not just "more fitness." It was also a specialization shift.

## Adaptation-pathway correlations

{corr_table}

The pace relationship pattern is important:

- pace tracked cadence very strongly
- pace also tracked power and vertical oscillation strongly
- pace tracked stride length strongly, but less strongly than cadence
- pace had only a weak relationship with vertical ratio

That supports a turnover-led adaptation interpretation rather than a simple "stride got longer and everything else followed" story.

## Timeline load and recovery correlations

{timeline_corr_table}

These relationships are exploratory rather than definitive, but they still help define the system:

- later running-specific load shows a meaningful positive relationship with late-window HRV coverage
- the resting-HR relationship is weaker across the full six-year timeline than the later-window within-run story might imply
- total activity load and running-specific load are not interchangeable in this dataset
- high-resolution running improvement happened alongside greater specialization, not alongside higher total activity

## Partial fueling overlap

{food_summary_table}

{food_corr_table}

The food overlap layer is valuable but limited:

- it covers only the late overlap window
- only a minority of overlap days contain parsed food-item rows
- it should be treated as partial context, not full six-year fueling truth
- it is strong enough to justify later branch studies, but not strong enough to anchor the whole flagship paper

## What this flagship study can and cannot answer

{scope_table}
"""


def build_discussion_markdown(
    yearly_overview: list[dict[str, object]],
    phase_rows: list[dict[str, object]],
    window_rows: list[dict[str, object]],
    highres_changes: list[dict[str, object]],
    correlation_rows: list[dict[str, object]],
    speed_rows: list[dict[str, object]],
    variability_rows: list[dict[str, object]],
) -> str:
    running_2020 = next(row for row in yearly_overview if row["year"] == 2020)
    indoor_2025 = next(row for row in yearly_overview if row["year"] == 2025)
    running_2026 = next(row for row in yearly_overview if row["year"] == 2026)
    pace_row = next(row for row in highres_changes if row["metric"] == "pace_min_per_mile")
    cadence_row = next(row for row in highres_changes if row["metric"] == "cadence_spm")
    stride_row = next(row for row in highres_changes if row["metric"] == "stride_length_m")
    vo_row = next(row for row in highres_changes if row["metric"] == "vertical_oscillation_cm")
    vr_row = next(row for row in highres_changes if row["metric"] == "vertical_ratio_pct")
    phase_2025 = next(row for row in phase_rows if row["phase_id"] == "P4")
    late_window = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    early_window = next(row for row in window_rows if row["window_label"] == "early_qc_window")
    cadence_corr = next(row for row in correlation_rows if row["relationship"] == "pace_vs_cadence")
    stride_corr = next(row for row in correlation_rows if row["relationship"] == "pace_vs_stride")
    speed_row = speed_rows[0]
    vo_var = next(row for row in variability_rows if row["metric"] == "vertical_oscillation_cm")
    vr_var = next(row for row in variability_rows if row["metric"] == "vertical_ratio_pct")

    return f"""# Study 000A Discussion

## What the six-year synthesis now says

This dataset does more than show that "things changed." It shows how they changed together.

At the whole-system level, the account moved through distinct training ecologies. {running_2020['year']} was nearly all running. {indoor_2025['year']} became the most indoor-conditioning-heavy year, with only {indoor_2025['running_share_of_structured_hours']:.1f}% of structured hours coming from running. By {running_2026['year']}, the system shifted back toward running dominance, with {running_2026['running_share_of_structured_hours']:.1f}% of structured hours coming from running.

That matters because the adaptation pattern cannot be reduced to "a faster runner over time." The training system itself reconfigured.

The phase model makes that concrete. The most intense restructuring happened during `{phase_2025['phase_label']}`, where indoor share averaged {phase_2025['mean_indoor_share_pct']:.1f}% and treadmill run share averaged {phase_2025['treadmill_run_share_pct']:.1f}% across the running rows in that phase.

## What changed inside the running system

Inside the later high-resolution window:

- pace changed {pace_row['pct_change']:.2f}%
- cadence changed {cadence_row['pct_change']:.2f}%
- stride length changed {stride_row['pct_change']:.2f}%
- vertical oscillation changed {vo_row['pct_change']:.2f}%
- vertical ratio changed {vr_row['pct_change']:.2f}%

This is the core insight. The system adapted, but not uniformly.

It also adapted inside a specific environment. The high-resolution window was not a balanced outdoor-running sample. Overall, QC-pass runs were treadmill-dominant, and the late comparison window was {late_window['treadmill_share_pct']:.1f}% treadmill. That is not a flaw, but it is a concrete fact about where the strongest adaptation signal emerged.

## Best-supported adaptation pathway

The strongest within-subject interpretation is descriptively turnover-led adaptation.

Pace tracked cadence at `r = {cadence_corr['correlation_r']:.4f}` and stride length at `r = {stride_corr['correlation_r']:.4f}`. Both matter, but cadence tracked pace more strongly than stride length. That suggests the system's most available lever of speed change was turnover rather than large stride expansion.

The decomposition analysis makes that more concrete. Across the first and last 30 QC-pass runs, observed speed increased by {speed_row['speed_gain_pct']:.2f}%. Of the cadence-stride speed gain, {speed_row['cadence_log_share_pct']:.1f}% came from cadence change and {speed_row['stride_log_share_pct']:.1f}% from stride-length change. Because speed is kinematically related to cadence and stride length, this is not proof of mechanism. It is a direct descriptive account of how the later speed gain was assembled inside the measured running system.

The window-context results deepen that conclusion. Recent running-specific load rose from {early_window['mean_running_hours_28d']:.2f} to {late_window['mean_running_hours_28d']:.2f} hours per 28 days, while total recent activity fell from {early_window['mean_total_activity_hours_28d']:.2f} to {late_window['mean_total_activity_hours_28d']:.2f}. Running therefore went from {early_window['running_share_of_28d_activity_pct']:.1f}% to {late_window['running_share_of_28d_activity_pct']:.1f}% of total recent activity. That is a specialization shift, not just more exercise.

For a study framed in the setting of altered biomechanics, that is a meaningful result. It suggests that the adaptation strategy may have depended more on the variables the system could flex most readily, while other movement features remained comparatively anchored.

## What this may mean biomechanically

The study does not prove exact mechanism, but it supports a plausible interpretation:

`adaptation occurred through selective flexibility under increasing specialization rather than across-the-board mechanical remodeling`

In plain language, the system changed a lot, but not everywhere equally. Some levers moved substantially. Others stayed comparatively more stable.

The monthly variability ranking strengthens that reading. `vertical_ratio_pct` had the lowest monthly coefficient of variation at {vr_var['monthly_cv_pct']:.2f}%, and `vertical_oscillation_cm` was the lowest-variability direct mechanics candidate at {vo_var['monthly_cv_pct']:.2f}%. So the conserved-mechanics signal shows up both in early-versus-late contrast and in month-to-month stability.

That is the main scientific value of the flagship study.

## What it still cannot prove

This paper still cannot tell us exactly how a non-compromised comparison system would have adapted under the same conditions, because there is no matched comparator inside the dataset. It also cannot prove why the conserved features remained more stable. Those are the next-study questions, not failures of the present study.
"""


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

1. This is still a single-subject longitudinal study and should not be interpreted as population inference.
2. The study is built from wearable-derived and rebuilt master tables rather than laboratory motion-capture or force-plate data.
3. High-resolution running mechanics are concentrated in the later device-supported window, not evenly across all six years.
4. The study can characterize within-subject adaptation under altered biomechanics, but it cannot directly quantify difference from a non-compromised comparison system without an external comparator.
5. Food overlap is partial and late-window only, so it is not strong enough to anchor six-year causal inference about fueling.
6. Recovery signals are uneven across years, especially HRV, which appears only in the later window.
7. The cadence-stride speed decomposition is descriptive, not causal, because cadence and stride length are mathematically related to speed and pace.
8. The study supports an adaptation-pathway interpretation, but it does not prove exact biomechanical or physiological mechanism.
9. The study should be read as a flagship synthesis of what the data can answer now, with branch studies needed for deeper mechanism, fueling, and comparison questions.
"""


def build_abstract_markdown(
    yearly_overview: list[dict[str, object]],
    phase_rows: list[dict[str, object]],
    window_rows: list[dict[str, object]],
    highres_changes: list[dict[str, object]],
    correlation_rows: list[dict[str, object]],
    speed_rows: list[dict[str, object]],
) -> str:
    pace_row = next(row for row in highres_changes if row["metric"] == "pace_min_per_mile")
    cadence_row = next(row for row in highres_changes if row["metric"] == "cadence_spm")
    stride_row = next(row for row in highres_changes if row["metric"] == "stride_length_m")
    vr_row = next(row for row in highres_changes if row["metric"] == "vertical_ratio_pct")
    vo_row = next(row for row in highres_changes if row["metric"] == "vertical_oscillation_cm")
    phase_2025 = next(row for row in phase_rows if row["phase_id"] == "P4")
    late_window = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    cadence_corr = next(row for row in correlation_rows if row["relationship"] == "pace_vs_cadence")
    stride_corr = next(row for row in correlation_rows if row["relationship"] == "pace_vs_stride")
    speed_row = speed_rows[0]
    return f"""# Study 000A Abstract

## Title

Six-Year Longitudinal Adaptation Under Altered Biomechanics: Workload, Performance, And Conserved Mechanics In A Nonstandard Movement System

## Abstract

This single-subject flagship synthesis study evaluated a broader six-year question than the earlier foundation studies: what changed, what remained comparatively anchored, and what adaptation pathway the data support in the setting of altered biomechanics and a nonstandard movement system. The study was built from the Stage 5 Garmin main package and packaged source tables covering annual ecology, daily master timeline context, late-window food overlap, and QC-pass high-resolution running mechanics.

The six-year record showed large ecological change, not just faster running. Structured training shifted from heavily running-dominant in 2020, through indoor-conditioning-heavy years, to a strongly running-dominant profile again in 2026. The strongest restructuring occurred in `{phase_2025['phase_label']}`, where indoor share averaged `{phase_2025['mean_indoor_share_pct']:.1f}%`. In the high-resolution mechanics window, adaptation was substantial. Comparing the first and last 30 QC-pass sessions, pace changed by `{pace_row['pct_change']:.2f}%`, cadence by `{cadence_row['pct_change']:.2f}%`, and stride length by `{stride_row['pct_change']:.2f}%`. At the same time, selected mechanics changed less, including `vertical_oscillation_cm` (`{vo_row['pct_change']:.2f}%`) and `vertical_ratio_pct` (`{vr_row['pct_change']:.2f}%`).

The strongest within-subject adaptation-pathway signal was descriptively turnover-led change. Pace tracked cadence at `r = {cadence_corr['correlation_r']:.4f}`, compared with `r = {stride_corr['correlation_r']:.4f}` for stride length, and cadence accounted for `{speed_row['cadence_log_share_pct']:.1f}%` of the cadence-stride speed gain. The later window was also `{late_window['treadmill_share_pct']:.1f}%` treadmill and showed markedly greater running specialization than the early window. This suggests that the system adapted through selective flexibility under specialization rather than uniform mechanical remodeling.

The study does not prove exact mechanism or direct difference from a non-compromised comparator system. It does, however, support a scientifically meaningful flagship conclusion: over six years, this altered system adapted substantially through phased ecological reconfiguration, treadmill-dominant specialization, and turnover-led running change, while selected movement features remained comparatively conserved.
"""


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

## Study type

Flagship longitudinal single-subject biomechanics, human-performance, and wearable-interpretation synthesis.

## Strongest framing

- six-year longitudinal adaptation under altered biomechanics
- ecological and mechanical adaptation in a nonstandard movement system
- selective flexibility and comparatively conserved movement features under changing load

## Best-fit outlets

- biomechanics or gait-research readers willing to engage a strong single-subject synthesis
- wearable-measurement and sports-technology audiences
- human-performance and self-research venues
- conference abstracts in biomechanics, locomotion, or measurement science

## What this paper can claim

- the six-year system changed substantially
- the changes were ecological as well as mechanical
- later running improvement occurred inside a treadmill-dominant specialization phase
- the strongest within-subject adaptation pathway was turnover-led rather than pure stride-expansion-led
- some movement features remained comparatively more conserved than the main adaptation variables

## What this paper should not claim

- direct clinical mechanism
- treatment implication
- population generalizability
- exact difference from a normative comparison system

## Best next branch studies

1. Physiological-cost branch study
   - Does the adaptation pathway carry a distinct exertion or recovery burden?

2. Fueling branch study
   - What does the partial intake layer suggest about performance under limited fueling support?

3. Comparator branch study
   - How does this within-subject pattern compare with literature-normalized or matched comparison profiles?
"""


def build_plain_language_summary_markdown(
    yearly_overview: list[dict[str, object]],
    phase_rows: list[dict[str, object]],
    window_rows: list[dict[str, object]],
    highres_changes: list[dict[str, object]],
    speed_rows: list[dict[str, object]],
) -> str:
    pace_row = next(row for row in highres_changes if row["metric"] == "pace_min_per_mile")
    cadence_row = next(row for row in highres_changes if row["metric"] == "cadence_spm")
    stride_row = next(row for row in highres_changes if row["metric"] == "stride_length_m")
    phase_2025 = next(row for row in phase_rows if row["phase_id"] == "P4")
    late_window = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    speed_row = speed_rows[0]
    return f"""# Study 000A Plain-Language Summary

## What this study asked

Over roughly six years of data, what changed, what stayed comparatively more stable, and what does that say about how this movement system adapted?

## What the study found

The answer is not just "running got better." The whole training system changed.

- the activity mix changed across years
- some years were more indoor-conditioning-heavy
- later years became much more running-focused again
- the biggest indoor pivot happened in `{phase_2025['phase_label']}`
- in the high-resolution running window, pace improved by `{pace_row['pct_change']:.2f}%`
- cadence increased by `{cadence_row['pct_change']:.2f}%`
- stride length increased by `{stride_row['pct_change']:.2f}%`
- about `{speed_row['cadence_log_share_pct']:.1f}%` of the cadence-stride speed gain came from cadence rather than stride length
- the late high-resolution window was `{late_window['treadmill_share_pct']:.1f}%` treadmill

At the same time, some movement features changed less than those bigger performance variables.

## Why that matters

The main insight is that adaptation did not happen evenly. The system appears to have changed through the levers it could move most effectively, while some movement features remained more anchored. The data also suggest that later improvement happened during a more specialized, treadmill-dominant running phase rather than during a simple rise in all activity.

The best-supported pattern is:

`phased adaptation plus comparatively conserved mechanics`

## What this does not prove

- it does not prove exact biomechanical mechanism
- it does not prove a clinical conclusion
- it does not prove how a fully normal comparison runner would look under the same conditions

## Best use of this study

This is the main synthesis study for the six-year record. It establishes the broad pattern. Later branch studies should explain cost, fueling, and comparator questions in more detail.
"""


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

## What this flagship study can answer now

{scope_table}

## Packaged source tables

{source_table}
"""


def build_manuscript_markdown() -> str:
    parts = [
        (ROOT / "reports" / "STUDY000A_ABSTRACT.md").read_text(encoding="utf-8").strip(),
        (ROOT / "reports" / "STUDY000A_METHODS.md").read_text(encoding="utf-8").strip(),
        (ROOT / "reports" / "STUDY000A_RESULTS.md").read_text(encoding="utf-8").strip(),
        (ROOT / "reports" / "STUDY000A_DISCUSSION.md").read_text(encoding="utf-8").strip(),
        (ROOT / "reports" / "STUDY000A_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 000A

This folder contains the flagship six-year longitudinal synthesis study built from the current Garmin Stage 5 master package.

## Why this study exists

Earlier studies established narrower foundation questions. This package is the broader synthesis study that asks what changed across the whole six-year system, what stayed comparatively anchored, and what adaptation pathway the data support.

## Headline findings

{findings}

## Main outputs

- `manuscript/STUDY000A_MANUSCRIPT.md`
- `reports/STUDY000A_ABSTRACT.md`
- `reports/STUDY000A_METHODS.md`
- `reports/STUDY000A_RESULTS.md`
- `reports/STUDY000A_DISCUSSION.md`
- `reports/STUDY000A_LIMITATIONS.md`
- `reports/STUDY000A_SUBMISSION_GUIDANCE.md`
- `reports/STUDY000A_PLAIN_LANGUAGE_SUMMARY.md`
- `appendices/APPENDIX_A_SCOPE_AND_SOURCES.md`
- `figures/`
- `outputs/`
- `source_tables/`
"""


def build_audit_markdown(
    source_manifest: list[dict[str, object]],
    window_rows: list[dict[str, object]],
    highres_changes: list[dict[str, object]],
    speed_rows: list[dict[str, object]],
    variability_rows: list[dict[str, object]],
) -> str:
    late_window = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    speed_row = speed_rows[0]
    vr_var = next(row for row in variability_rows if row["metric"] == "vertical_ratio_pct")
    vo_var = next(row for row in variability_rows if row["metric"] == "vertical_oscillation_cm")
    pace_row = next(row for row in highres_changes if row["metric"] == "pace_min_per_mile")
    cadence_row = next(row for row in highres_changes if row["metric"] == "cadence_spm")
    stride_row = next(row for row in highres_changes if row["metric"] == "stride_length_m")

    return f"""# Study 000A Audit

## Structural audit

- Pass: the package is self-contained and reads only from bundled `source_tables`.
- Pass: the packaged source manifest covers `{len(source_manifest)}` study input tables.
- Pass: the package includes manuscript, methods, results, discussion, limitations, figures, outputs, appendix, and source manifest.
- Pass: the regenerated package includes five figures and expanded flagship outputs for phase model, window context, speed decomposition, and monthly variability ranking.

## Scientific audit

- Pass: the flagship claim is broader and more concrete than before. Pace changed `{pace_row['pct_change']:.2f}%`, cadence `{cadence_row['pct_change']:.2f}%`, and stride length `{stride_row['pct_change']:.2f}%` in the QC-pass high-resolution window.
- Pass: the conserved-mechanics claim is supported in two ways: early-late contrast and monthly variability ranking.
- Pass: `vertical_ratio_pct` has the lowest monthly coefficient of variation at `{vr_var['monthly_cv_pct']:.2f}%`, and `vertical_oscillation_cm` is the lowest-variability direct mechanics candidate at `{vo_var['monthly_cv_pct']:.2f}%`.
- Pass: the turnover-led interpretation is more concrete because cadence accounts for `{speed_row['cadence_log_share_pct']:.1f}%` of the cadence-stride speed gain, while stride accounts for `{speed_row['stride_log_share_pct']:.1f}%`.
- Pass: treadmill specialization is explicit rather than hidden. The late comparison window is `{late_window['treadmill_share_pct']:.1f}%` treadmill.

## Residual cautions

- The cadence-stride speed decomposition is descriptive rather than causal.
- High-resolution mechanics remain concentrated in the later device-supported window.
- The study still has no external comparator and therefore cannot directly quantify difference from a non-compromised system.
- The food layer remains partial late-window context only.

## Audit verdict

`Study 000A` is structurally sound and scientifically stronger after revision. The flagship now concludes more concretely than "things changed": it supports phased ecological reconfiguration, treadmill-dominant specialization, uneven mechanical adaptation, and descriptively turnover-led speed gain in a nonstandard movement system. Remaining gaps are mechanism and comparator gaps, not package-integrity failures.
"""


def build_project_state(
    source_manifest: list[dict[str, object]],
    yearly_overview: list[dict[str, object]],
    phase_rows: list[dict[str, object]],
    window_rows: list[dict[str, object]],
    highres_changes: list[dict[str, object]],
    correlation_rows: list[dict[str, object]],
    speed_rows: list[dict[str, object]],
) -> dict[str, object]:
    pace_row = next(row for row in highres_changes if row["metric"] == "pace_min_per_mile")
    cadence_row = next(row for row in highres_changes if row["metric"] == "cadence_spm")
    stride_row = next(row for row in highres_changes if row["metric"] == "stride_length_m")
    running_2026 = next(row for row in yearly_overview if row["year"] == 2026)
    phase_2025 = next(row for row in phase_rows if row["phase_id"] == "P4")
    late_window = next(row for row in window_rows if row["window_label"] == "late_qc_window")
    cadence_corr = next(row for row in correlation_rows if row["relationship"] == "pace_vs_cadence")
    stride_corr = next(row for row in correlation_rows if row["relationship"] == "pace_vs_stride")
    speed_row = speed_rows[0]
    return {
        "study_id": "STUDY-000A",
        "title": "Six-Year Longitudinal Adaptation Under Altered Biomechanics",
        "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_synthesis_complete",
        "source_manifest": source_manifest,
        "headline_findings": [
            "The six-year record shows ecological reconfiguration, not just faster running.",
            f"The most intense restructuring occurred in {phase_2025['phase_label']}, where indoor share averaged {phase_2025['mean_indoor_share_pct']:.1f}%.",
            f"Pace changed {pace_row['pct_change']:.2f}%, cadence changed {cadence_row['pct_change']:.2f}%, and stride length changed {stride_row['pct_change']:.2f}% in the QC-pass high-resolution window.",
            f"Pace tracked cadence more strongly (r={cadence_corr['correlation_r']:.4f}) than stride length (r={stride_corr['correlation_r']:.4f}); cadence accounted for {speed_row['cadence_log_share_pct']:.1f}% of the cadence-stride speed gain.",
            f"The late high-resolution window was {late_window['treadmill_share_pct']:.1f}% treadmill, showing that the strongest adaptation signal emerged in a specialized treadmill-dominant context.",
            f"By 2026, running accounted for {running_2026['running_share_of_structured_hours']:.1f}% of structured hours, showing a major shift back toward running dominance.",
            "Some movement features remained comparatively more conserved than the main adaptation variables.",
        ],
        "branch_study_targets": [
            "Physiological cost under altered mechanics",
            "Fueling, recovery, and performance overlap",
            "Comparator study versus normative or matched external profiles",
        ],
    }


def main() -> None:
    ROOT.mkdir(parents=True, exist_ok=True)

    timeline_rows = read_csv(TIMELINE_PATH)
    food_rows = read_csv(FOOD_PATH)
    yearly_rows = read_csv(YEARLY_PATH)
    monthly_rows = read_csv(MONTHLY_PATH)
    run_rows = read_csv(RUNS_PATH)

    source_manifest = build_source_manifest([TIMELINE_PATH, FOOD_PATH, YEARLY_PATH, MONTHLY_PATH, RUNS_PATH])
    yearly_overview = build_yearly_overview(yearly_rows)
    phase_rows = build_phase_model(monthly_rows, timeline_rows, run_rows)
    window_rows = build_window_context(run_rows, timeline_rows)
    yearly_recovery = build_yearly_recovery_summary(timeline_rows)
    highres_changes, correlation_rows, monthly_highres_rows = build_highres_summary(run_rows)
    speed_rows = build_speed_decomposition(run_rows)
    variability_rows = build_monthly_variability(monthly_highres_rows)
    window_change_rows = build_window_change(window_rows)
    timeline_corr_rows = build_timeline_correlations(timeline_rows)
    food_summary, food_corr_rows = build_food_summary(food_rows)
    scope_rows = build_scope_claims()
    project_state = build_project_state(source_manifest, yearly_overview, phase_rows, window_rows, highres_changes, correlation_rows, speed_rows)

    write_csv(
        ROOT / "manifest" / "source_table_manifest.csv",
        source_manifest,
        ["table_name", "relative_path", "bytes", "sha256"],
    )
    write_csv(
        ROOT / "outputs" / "study000a_yearly_overview.csv",
        yearly_overview,
        list(yearly_overview[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_yearly_recovery_summary.csv",
        yearly_recovery,
        list(yearly_recovery[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_phase_model.csv",
        phase_rows,
        list(phase_rows[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_window_context.csv",
        window_rows,
        list(window_rows[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_highres_early_late_changes.csv",
        highres_changes,
        ["metric", "early_mean", "late_mean", "pct_change"],
    )
    write_csv(
        ROOT / "outputs" / "study000a_highres_correlations.csv",
        correlation_rows,
        ["relationship", "correlation_r", "n_points"],
    )
    write_csv(
        ROOT / "outputs" / "study000a_speed_gain_decomposition.csv",
        speed_rows,
        list(speed_rows[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_monthly_highres_summary.csv",
        monthly_highres_rows,
        list(monthly_highres_rows[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_monthly_variability_rank.csv",
        variability_rows,
        list(variability_rows[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_window_change_summary.csv",
        window_change_rows,
        list(window_change_rows[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_timeline_correlations.csv",
        timeline_corr_rows,
        ["relationship", "correlation_r", "n_points"],
    )
    write_csv(
        ROOT / "outputs" / "study000a_food_overlap_summary.csv",
        food_summary,
        list(food_summary[0].keys()),
    )
    write_csv(
        ROOT / "outputs" / "study000a_food_overlap_correlations.csv",
        food_corr_rows,
        ["relationship", "correlation_r", "n_points"],
    )
    write_csv(
        ROOT / "outputs" / "study000a_scope_claims.csv",
        scope_rows,
        ["category", "question", "answer"],
    )

    build_figure_yearly_ecology(yearly_overview, ROOT / "figures" / "figure01_yearly_training_ecology.svg")
    build_figure_adaptation_scatter(run_rows, ROOT / "figures" / "figure02_adaptation_pathways_scatter.svg")
    build_figure_pct_change(highres_changes, ROOT / "figures" / "figure03_highres_pct_change.svg")
    build_figure_window_context(window_rows, ROOT / "figures" / "figure04_window_context_shift.svg")
    build_figure_speed_decomposition(speed_rows, ROOT / "figures" / "figure05_speed_gain_decomposition.svg")

    write_text(ROOT / "reports" / "STUDY000A_METHODS.md", build_methods_markdown(source_manifest))
    write_text(
        ROOT / "reports" / "STUDY000A_RESULTS.md",
        build_results_markdown(
            yearly_overview,
            phase_rows,
            window_rows,
            yearly_recovery,
            highres_changes,
            correlation_rows,
            speed_rows,
            variability_rows,
            window_change_rows,
            timeline_corr_rows,
            food_summary,
            food_corr_rows,
            scope_rows,
        ),
    )
    write_text(
        ROOT / "reports" / "STUDY000A_DISCUSSION.md",
        build_discussion_markdown(yearly_overview, phase_rows, window_rows, highres_changes, correlation_rows, speed_rows, variability_rows),
    )
    write_text(ROOT / "reports" / "STUDY000A_LIMITATIONS.md", build_limitations_markdown())
    write_text(
        ROOT / "reports" / "STUDY000A_ABSTRACT.md",
        build_abstract_markdown(yearly_overview, phase_rows, window_rows, highres_changes, correlation_rows, speed_rows),
    )
    write_text(ROOT / "reports" / "STUDY000A_SUBMISSION_GUIDANCE.md", build_submission_guidance_markdown())
    write_text(
        ROOT / "reports" / "STUDY000A_PLAIN_LANGUAGE_SUMMARY.md",
        build_plain_language_summary_markdown(yearly_overview, phase_rows, window_rows, highres_changes, speed_rows),
    )
    write_text(
        ROOT / "appendices" / "APPENDIX_A_SCOPE_AND_SOURCES.md",
        build_appendix_markdown(scope_rows, source_manifest),
    )
    write_text(
        ROOT / "reports" / "STUDY000A_AUDIT.md",
        build_audit_markdown(source_manifest, window_rows, highres_changes, speed_rows, variability_rows),
    )
    write_text(ROOT / "manuscript" / "STUDY000A_MANUSCRIPT.md", build_manuscript_markdown())
    write_text(ROOT / "README.md", build_readme_markdown(project_state))
    write_text(ROOT / "manifest" / "study000a_project_state.json", json.dumps(project_state, indent=2))


if __name__ == "__main__":
    main()
