"""Descriptive figures and tables for WP36 (Phase 4a).

Inputs (read-only): data/processed/listings.parquet, data/processed/price_events.parquet,
data/interim/sync_log.parquet, data/processed/sources_registry.csv,
data/processed/sample_construction.csv.

Outputs: figures/fig01_coverage, fig03_listing_age_and_cuts, fig04_price_paths (PDF + PNG,
copied to paper/figures/), paper/tables/tab01_sources.tex, tab02_sample_construction.tex,
tab03_descriptives.tex, results/descriptives/*.csv|json (every number quoted in captions is
in results/descriptives/summary.json).

Conventions: day index t = UTC calendar days since 2026-08-12 (day 0, partial: crawl started
12:30 UTC); day 81 (2026-11-01) is partial (snapshot 06:49 UTC). "Exit" = observed
disappearance confirmed by >= 7 days of absence before the snapshot — not a sale.
"""
from __future__ import annotations

import hashlib
import json
import math
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl

from . import config as C
from . import plotstyle as PS

OUT = C.RESULTS / "descriptives"
DAY0 = math.floor(C.WINDOW_START_TS / 86400)
T = int(math.floor(C.SNAPSHOT_TS / 86400) - DAY0) + 1          # number of day indices (82)
DAY_LABELS = {0: "2026-08-12", 81: "2026-11-01"}
WITHIN_DAYS = 30                                                 # exit-within-30-days horizon
CLASS_SHORT = {"direct_dealer": "Dealer websites", "full_portal": "Full-inventory portals", "banner": "Dealer groups",
               "sampled_portal": "Sampled portals", "private_marketplace": "Private sellers (Kijiji)", "small_portal": "Small portals"}


# ----------------------------------------------------------------------------- helpers
def day_of(ts: pl.Expr) -> pl.Expr:
    return ((ts / 86400).floor() - DAY0).cast(pl.Int32)


def fmt_int(x) -> str:
    """Integer with thin-space thousands separators (LaTeX)."""
    return f"{int(round(float(x))):,}".replace(",", "\\,")


def fmt1(x) -> str:
    return "--" if x is None or (isinstance(x, float) and np.isnan(x)) else f"{float(x):.1f}"


def fmt_pct(x) -> str:
    return "--" if x is None or (isinstance(x, float) and np.isnan(x)) else f"{100 * float(x):.1f}"


def write_tex(path: Path, lines: list[str]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text("\n".join(lines) + "\n", encoding="utf-8")


def load() -> tuple[pl.DataFrame, pl.DataFrame, pl.DataFrame, pd.DataFrame, pd.DataFrame]:
    L = pl.read_parquet(C.PROCESSED / "listings.parquet")
    E = pl.read_parquet(C.PROCESSED / "price_events.parquet")
    S = pl.read_parquet(C.INTERIM / "sync_log.parquet")
    reg = pd.read_csv(C.PROCESSED / "sources_registry.csv")
    sc = pd.read_csv(C.PROCESSED / "sample_construction.csv")
    L = L.with_columns(day_of(pl.col("entry_ts")).alias("e_day"),
                       pl.min_horizontal(day_of(pl.col("end_ts")), pl.lit(T - 1)).alias("x_day"))
    return L, E, S, reg, sc


# ----------------------------------------------------------------------------- fig01
def daily_coverage(L: pl.DataFrame, S: pl.DataFrame) -> tuple[pd.DataFrame, dict]:
    """Daily active stock / entries / confirmed exits in S2L, and crawl cycles per day."""
    U = L.filter(pl.col("s2l_dynamic_with_stock"))
    e = U["e_day"].to_numpy(); x = U["x_day"].to_numpy(); ex = U["exit_confirmed"].to_numpy()
    days = np.arange(T)
    entries = np.bincount(e, minlength=T)[:T]
    exits = np.bincount(x[ex], minlength=T)[:T]
    ends = np.bincount(x, minlength=T)[:T]
    # observed-active on day t: entered on or before t and not yet ended before t (ends are inclusive days)
    active = np.cumsum(entries) - np.concatenate([[0], np.cumsum(ends)[:-1]])
    # crawl cycles per day = sync rows / distinct sources logged that day (one row per source per cycle)
    sd = S.with_columns(day_of(pl.col("ts")).alias("day"))
    agg = (sd.group_by("day").agg(pl.len().alias("rows"), pl.col("source").n_unique().alias("n_sources"),
                                  (pl.col("ok") == 0).sum().alias("failed_rows"))
           .with_columns((pl.col("rows") / pl.col("n_sources")).alias("cycles")).sort("day").to_pandas())
    cyc = np.zeros(T); cyc[agg["day"].to_numpy()] = agg["cycles"].to_numpy()
    nsrc = np.zeros(T, dtype=int); nsrc[agg["day"].to_numpy()] = agg["n_sources"].to_numpy()
    fail = np.zeros(T); fail[agg["day"].to_numpy()] = (agg["failed_rows"] / agg["rows"]).to_numpy()
    # largest gap between consecutive crawls of the same source (hours)
    gaps = (sd.sort(["source", "ts"]).with_columns(pl.col("ts").diff().over("source").alias("dt"))
            .filter(pl.col("dt").is_not_null()).select((pl.col("dt") / 3600).alias("h")))
    # connector start-up stock: for portal sources first crawled after day 0, the first crawl day's entries are
    # initial stock, not an entry flow (not flagged left_censored because the cut-off is 2026-08-13 00:00 UTC)
    first_day = sd.group_by("source").agg(pl.col("day").min().alias("first")).filter(pl.col("first") > 0)
    st = (U.join(first_day, on="source", how="inner").filter(pl.col("e_day") == pl.col("first"))
          .group_by(["source", "e_day"]).agg(pl.len().alias("n")).filter(pl.col("n") >= 1000).sort("e_day"))
    startup_src = {str(int(r["e_day"])): r["source"] for r in st.iter_rows(named=True)}
    startup_n = {str(int(r["e_day"])): int(entries[int(r["e_day"])]) for r in st.iter_rows(named=True)}
    startup_listings = {r["source"]: int(r["n"]) for r in st.iter_rows(named=True)}
    cov = pd.DataFrame({"day": days, "active_s2l": active, "entries_s2l": entries, "exits_confirmed_s2l": exits,
                        "crawl_cycles": cyc, "sources_logged": nsrc, "failed_share": fail})
    cov["date_utc"] = pd.to_datetime((DAY0 + days) * 86400, unit="s").strftime("%Y-%m-%d")
    full = cov[(cov.day >= 1) & (cov.day <= T - 2)]
    info = {"T_days": int(T), "day0_initial_stock_s2l": int(entries[0]),
            "connector_startup_entries": startup_n, "connector_startup_source": startup_src, "connector_startup_listings_s2l": startup_listings,
            "connector_startup_listings_s2_total": int(L.filter(pl.col("s2_dynamic")).join(first_day, on="source", how="inner")
                                                       .filter((pl.col("e_day") == pl.col("first")) & pl.col("source").is_in(list(startup_listings))).height),
            "entries_per_full_day_mean_excl_startup": float(full[~full.day.isin([int(k) for k in startup_n])].entries_s2l.mean()), "entries_day0_left_censored": int(U["left_censored"].sum()),
            "active_s2l_min": int(full.active_s2l.min()), "active_s2l_max": int(full.active_s2l.max()),
            "active_s2l_mean_full_days": float(full.active_s2l.mean()), "active_s2l_at_last_full_day": int(active[T - 2]),
            "entries_per_full_day_mean": float(full.entries_s2l.mean()), "entries_per_full_day_median": float(full.entries_s2l.median()),
            "exits_per_full_day_mean_confirmable": float(cov[(cov.day >= 1) & (cov.day <= T - 1 - C.EXIT_CONFIRM_DAYS)].exits_confirmed_s2l.mean()),
            "crawl_cycles_per_full_day_mean": float(full.crawl_cycles.mean()), "crawl_cycles_per_full_day_min": float(full.crawl_cycles.min()),
            "crawl_cycles_per_full_day_max": float(full.crawl_cycles.max()), "crawl_cycles_total": int(round(cyc.sum())),
            "days_with_fewer_than_4_cycles": [int(d) for d in full.day[full.crawl_cycles < 4]],
            "max_gap_hours_same_source": float(gaps["h"].max()), "median_gap_hours_same_source": float(gaps["h"].median()),
            "sources_logged_min": int(full.sources_logged.min()), "sources_logged_max": int(full.sources_logged.max()),
            "failed_source_cycles_share": float((S["ok"] == 0).mean())}
    return cov, info


def region_stock(L: pl.DataFrame) -> pd.DataFrame:
    """Active canonical auto listings by region at the snapshot (all sources)."""
    A = L.filter((pl.col("active") == 1) & pl.col("is_canonical") & (pl.col("kind") == "auto"))
    r = (A.group_by(pl.col("region").fill_null("Region missing")).agg(pl.len().alias("active_canonical"))
         .sort("active_canonical", descending=True).to_pandas())
    r["share"] = r.active_canonical / r.active_canonical.sum()
    return r


def fig01(cov: pd.DataFrame, reg: pd.DataFrame, info: dict) -> None:
    fig = plt.figure(figsize=(6.3, 4.6))
    gs = fig.add_gridspec(3, 2, width_ratios=[1.35, 1.0])
    ax_a1 = fig.add_subplot(gs[0, 0]); ax_a2 = fig.add_subplot(gs[1, 0], sharex=ax_a1); ax_b = fig.add_subplot(gs[2, 0], sharex=ax_a1)
    ax_c = fig.add_subplot(gs[:, 1])
    d = cov.day.to_numpy()
    c_act, c_ent, c_ex = PS.CAT[0], PS.CAT[2], PS.CAT[1]
    # (a1) active stock
    ax_a1.plot(d, cov.active_s2l / 1000, color=c_act, lw=1.4)
    ax_a1.set_ylabel("Active listings\n(thousands)")
    ax_a1.set_ylim(0, None)
    ax_a1.set_title("(a) Dynamic-reliable universe (S2L): stock and flows", loc="left")
    # (a2) entries / exits (day 0 = initial left-censored stock, not an entry flow)
    m = d >= 1
    ent = cov.entries_s2l.to_numpy().astype(float)
    spikes = {int(k): int(v) for k, v in info["connector_startup_entries"].items()}
    ycap = float(np.quantile(ent[m & ~np.isin(d, list(spikes))], 1.0)) * 1.5
    ax_a2.plot(d[m], np.minimum(ent[m], ycap), color=c_ent, lw=1.1, label="Entries")
    mx = m & (d <= T - 1 - C.EXIT_CONFIRM_DAYS)
    ax_a2.plot(d[mx], cov.exits_confirmed_s2l[mx], color=c_ex, lw=1.1, label="Confirmed exits")
    ax_a2.plot(list(spikes), [ycap] * len(spikes), ls="none", marker="^", ms=4, color=c_ent, clip_on=False)
    note = "\u25b2 truncated: connector\n   start-up stock\n" + "\n".join(
        f"   day {k}: {v:,} ({info['connector_startup_source'][str(k)]})".replace(",", "\u2009") for k, v in spikes.items())
    ax_a2.text(0.14, 0.97, note, transform=ax_a2.transAxes, fontsize=6, color=PS.TEXT2, ha="left", va="top", linespacing=1.1)
    ax_a2.axvspan(T - 1 - C.EXIT_CONFIRM_DAYS, T - 1, color=PS.GRID, alpha=0.6, lw=0, zorder=0)
    ax_a2.set_ylabel("Listings per day")
    ax_a2.set_ylim(0, ycap)
    ax_a2.legend(loc="upper right", bbox_to_anchor=(1.0, 1.02), ncol=1, handlelength=1.4, labelspacing=0.3)
    # (b) crawl cycles per day
    ax_b.bar(d, cov.crawl_cycles, width=0.8, color=PS.CAT[0], alpha=0.85)
    ax_b.set_ylabel("Crawl cycles\nper day")
    ax_b.set_xlabel("Day index (0 = 2026-08-12 UTC; days 0 and 81 are partial)")
    ax_b.set_title("(b) Crawl cycles per day, all sources", loc="left")
    ax_b.set_ylim(0, max(12, cov.crawl_cycles.max() * 1.1))
    ax_b.set_xlim(-1, T)
    for ax in (ax_a1, ax_a2):
        plt.setp(ax.get_xticklabels(), visible=False)
    # (c) regions
    rr = reg.sort_values("active_canonical", ascending=True)
    cols = [PS.color_for(PS.REGIONS, k) for k in rr.region]
    ax_c.barh(rr.region, rr.active_canonical / 1000, color=cols, height=0.72)
    ax_c.set_xlabel("Active canonical listings (thousands)")
    ax_c.set_title("(c) Stock by region, 2026-11-01", loc="left")
    ax_c.tick_params(axis="y", labelsize=7)
    ax_c.grid(axis="y", visible=False)
    for y, v in enumerate(rr.active_canonical):
        ax_c.text(v / 1000 + 0.08, y, f"{v:,}".replace(",", " "), va="center", fontsize=6, color=PS.TEXT2)
    ax_c.set_xlim(0, rr.active_canonical.max() / 1000 * 1.22)
    PS.save(fig, "fig01_coverage")


# ----------------------------------------------------------------------------- fig03
def duration_and_cuts(L: pl.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame, dict]:
    S2 = L.filter(pl.col("s2_dynamic"))
    done = S2.filter(pl.col("exit_confirmed"))["duration_days"].to_numpy()
    done = np.sort(done)
    ecdf = pd.DataFrame({"duration_days": done, "ecdf": np.arange(1, done.size + 1) / done.size})
    fc = -100 * S2.filter(pl.col("n_cuts") >= 1)["first_cut_pct"].to_numpy()
    cuts = pd.DataFrame({"first_cut_pct": np.sort(fc)})
    bc = (S2.group_by("source_class").agg(pl.len().alias("n"), (pl.col("n_cuts") >= 1).sum().alias("n_with_cut"),
                                          pl.col("exit_confirmed").sum().alias("n_exit_confirmed"))
          .with_columns((pl.col("n_with_cut") / pl.col("n")).alias("share_with_cut")).to_pandas())
    bc["order"] = bc.source_class.map(PS.CLASS_ORDER.index); bc = bc.sort_values("order").drop(columns="order")
    S2L = L.filter(pl.col("s2l_dynamic_with_stock"))
    q = lambda a, p: float(np.quantile(a, p))
    info = {"S2_n": S2.height, "S2_exit_confirmed": int(S2["exit_confirmed"].sum()),
            "S2_right_censored_share": float(1 - S2["exit_confirmed"].mean()),
            "S2L_n": S2L.height, "S2L_left_censored_n": int(S2L["left_censored"].sum()), "S2L_left_censored_share": float(S2L["left_censored"].mean()),
            "S2_completed_duration_median_days": q(done, .5), "S2_completed_duration_p25_days": q(done, .25), "S2_completed_duration_p75_days": q(done, .75),
            "S2_completed_duration_mean_days": float(done.mean()),
            "S2_completed_share_within_7d": float((done <= 7).mean()), "S2_completed_share_within_30d": float((done <= 30).mean()),
            "S2_share_with_cut": float((S2["n_cuts"] >= 1).mean()), "S2_n_with_cut": int((S2["n_cuts"] >= 1).sum()),
            "S2_share_with_2plus_cuts": float((S2["n_cuts"] >= 2).mean()),
            "S2_first_cut_pct_median": q(fc, .5), "S2_first_cut_pct_mean": float(fc.mean()), "S2_first_cut_pct_p25": q(fc, .25), "S2_first_cut_pct_p75": q(fc, .75),
            "S2_first_cut_share_le_5pct": float((fc <= 5).mean()), "S2_first_cut_share_gt_10pct": float((fc > 10).mean()),
            "S2_share_with_cut_by_class": {r.source_class: float(r.share_with_cut) for r in bc.itertuples()}}
    return ecdf, cuts, bc, info


def fig03(ecdf: pd.DataFrame, cuts: pd.DataFrame, bc: pd.DataFrame, info: dict) -> None:
    fig, axes = plt.subplots(1, 3, figsize=(6.3, 2.3), gridspec_kw={"width_ratios": [1.15, 1.15, 1.0]})
    a, b, c = axes
    # (a) ECDF of completed durations
    x = ecdf.duration_days.to_numpy(); y = ecdf.ecdf.to_numpy()
    a.step(x, y, where="post", color=PS.CAT[0], lw=1.2)
    med = info["S2_completed_duration_median_days"]
    a.axvline(med, color=PS.NEUTRAL, lw=0.7, ls=(0, (3, 2)))
    a.text(med + 1.5, 0.08, f"median {med:.1f} d", fontsize=7, color=PS.TEXT2)
    a.set_xlim(0, 75); a.set_ylim(0, 1)
    a.set_xlabel("Completed duration (days)"); a.set_ylabel("Cumulative share of spells")
    a.set_title("(a) Completed spells, S2", loc="left")
    # (b) first markdown size
    v = cuts.first_cut_pct.to_numpy()
    bins = np.arange(1, 25.5, 0.5)
    b.hist(np.clip(v, 1, 25.4), bins=bins, color=PS.CAT[1], alpha=0.9, lw=0, density=False, weights=np.full(v.size, 100 / v.size))
    medc = info["S2_first_cut_pct_median"]
    b.axvline(medc, color=PS.NEUTRAL, lw=0.7, ls=(0, (3, 2)))
    b.text(medc + 0.6, b.get_ylim()[1] * 0.97, f"median {medc:.1f} %", fontsize=7, color=PS.TEXT2, va="top")
    b.set_xlim(1, 25.5)
    b.set_xlabel("First markdown (% of entry price; \u2265 25 pooled)"); b.set_ylabel("Share of first markdowns (%)")
    b.set_title("(b) First markdown size, S2", loc="left")
    # (c) share with >=1 cut by source class
    labels = [CLASS_SHORT[k] for k in bc.source_class]
    cols = [PS.color_for(PS.CLASS_ORDER, k) for k in bc.source_class]
    c.bar(range(len(bc)), 100 * bc.share_with_cut, color=cols, width=0.65)
    c.set_xticks(range(len(bc)), [l.replace(" ", "\n", 1) for l in labels], fontsize=7)
    for i, (s, n) in enumerate(zip(bc.share_with_cut, bc.n)):
        c.text(i, 100 * s + 0.6, f"{100 * s:.1f} %", ha="center", fontsize=7, color=PS.TEXT2)
    c.set_ylabel("Listings with ≥1 markdown (%)")
    c.set_ylim(0, max(100 * bc.share_with_cut) * 1.25)
    c.grid(axis="x", visible=False)
    c.set_title("(c) Markdown incidence", loc="left"); c.set_xlabel("Source class (S2)")
    PS.save(fig, "fig03_listing_age_and_cuts")


# ----------------------------------------------------------------------------- fig04
def select_paths(L: pl.DataFrame) -> tuple[pl.DataFrame, dict]:
    """Deterministic selection: S2 listings with >= 3 sustained cuts and 30 <= duration <= 70 days,
    2 per body type among {SUV, Sedan, Pickup} (the three largest body types), chosen by the
    lowest SHA-256 of the anonymised listing_id (hex, ascending)."""
    bodies = ["SUV", "Sedan", "Pickup"]
    cand = (L.filter(pl.col("s2_dynamic") & (pl.col("n_cuts") >= 3) & (pl.col("duration_days") >= 30) & (pl.col("duration_days") <= 70)
                     & pl.col("body_filled").is_in(bodies))
            .with_columns(pl.col("listing_id").map_elements(lambda s: hashlib.sha256(s.encode()).hexdigest(), return_dtype=pl.String).alias("h")))
    picks = (cand.sort(["body_filled", "h"]).group_by("body_filled", maintain_order=True).head(2)
             .with_columns(pl.col("body_filled").replace_strict({b: i for i, b in enumerate(bodies)}, return_dtype=pl.Int64).alias("bo"))
             .sort(["bo", "h"]).drop("bo"))
    picks = picks.with_columns(pl.Series("label", [f"Listing {chr(65 + i)}" for i in range(picks.height)]))
    rule = {"rule": "S2 listings with n_cuts >= 3 and 30 <= duration_days <= 70, body in {SUV, Sedan, Pickup}; "
                    "within each body type the two listings with the lowest SHA-256(listing_id) hex digest; labelled A-F in body order.",
            "n_candidates": cand.height, "n_candidates_by_body": {b: int(cand.filter(pl.col("body_filled") == b).height) for b in bodies},
            "selected": [{"label": r["label"], "listing_id": r["listing_id"], "sha256_prefix": r["h"][:12], "make": r["make"], "body": r["body_filled"],
                          "year": int(r["year"]), "n_cuts": int(r["n_cuts"]), "duration_days": round(float(r["duration_days"]), 1),
                          "exit_confirmed": bool(r["exit_confirmed"]), "entry_price": float(r["entry_price"]), "last_price": float(r["last_price"]),
                          "total_change_pct": round(100 * float(r["total_change_pct"]), 1)}
                         for r in picks.iter_rows(named=True)]}
    return picks, rule


def fig04(picks: pl.DataFrame, E: pl.DataFrame) -> pd.DataFrame:
    ev = (E.filter(pl.col("listing_id").is_in(picks["listing_id"].to_list())).join(
        picks.select(["listing_id", "label", "entry_ts", "end_ts", "exit_confirmed", "make", "body_filled", "year"]), on="listing_id")
        .with_columns(((pl.col("ts") - pl.col("entry_ts")) / 86400).alias("days"))
        .sort(["label", "ts"]))
    fig, axes = plt.subplots(2, 3, figsize=(6.3, 3.7), sharex=True)
    rows = []
    for ax, (lab, g) in zip(axes.flat, ev.group_by("label", maintain_order=True)):
        lab = lab[0] if isinstance(lab, tuple) else lab
        d = g["days"].to_numpy(); p = g["price"].to_numpy() / 1000
        d0 = np.concatenate([[0.0], d[1:]]) if d.size else d        # entry level starts at day 0
        dur = float((g["end_ts"][0] - g["entry_ts"][0]) / 86400)
        body = g["body_filled"][0]; col = PS.color_for(PS.BODY_ORDER, body)
        ax.step(np.concatenate([d0, [dur]]), np.concatenate([p, [p[-1]]]), where="post", color=col, lw=1.3)
        ax.plot(d[1:], p[1:], ls="none", marker="o", ms=2.6, color=col, mfc=PS.SURFACE, mew=0.9)
        if bool(g["exit_confirmed"][0]):
            ax.plot([dur], [p[-1]], ls="none", marker="x", ms=4.5, color=PS.TEXT2, mew=0.9)
        ax.set_title(f"{lab} — {g['make'][0]} {body} {int(g['year'][0])}", loc="left", fontsize=8)
        ax.set_xlim(0, 72)
        lo, hi = p.min(), p.max(); pad = max(0.6, 0.12 * (hi - lo))
        ax.set_ylim(lo - pad, hi + pad)
        for r in g.iter_rows(named=True):
            rows.append({"label": lab, "days_since_entry": round(r["days"], 3), "price_cad": r["price"], "kind": r["kind"]})
        rows.append({"label": lab, "days_since_entry": round(dur, 3), "price_cad": g["price"][-1], "kind": "exit" if bool(g["exit_confirmed"][0]) else "censored"})
    for ax in axes[1]:
        ax.set_xlabel("Days since entry")
    fig.supylabel("Asking price (CAD, thousands)", fontsize=9)
    from matplotlib.lines import Line2D
    h = [Line2D([], [], color=PS.TEXT2, lw=1.2, label="Sustained price level"),
         Line2D([], [], ls="none", marker="o", ms=2.6, color=PS.TEXT2, mfc=PS.SURFACE, label="Price change"),
         Line2D([], [], ls="none", marker="x", ms=4.5, color=PS.TEXT2, label="Confirmed exit")]
    fig.legend(handles=h, loc="outside lower center", ncol=3, fontsize=7, handlelength=1.6, columnspacing=1.6)
    PS.save(fig, "fig04_price_paths")
    return pd.DataFrame(rows)


# ----------------------------------------------------------------------------- tables
def tab01(L: pl.DataFrame, reg: pd.DataFrame) -> pd.DataFrame:
    by = (L.group_by("source_class").agg(pl.col("source").n_unique().alias("n_sources_with_listings"), pl.col("dealer_id").n_unique().alias("n_dealers"),
                                         pl.len().alias("listings"), pl.col("active").sum().alias("active_at_snapshot"),
                                         pl.col("s1_hedonic").any().alias("in_s1"), pl.col("s2_dynamic").any().alias("in_s2"),
                                         pl.col("source").filter(pl.col("s2_dynamic")).n_unique().alias("n_sources_in_s2"),
                                         pl.col("s1_hedonic").sum().alias("n_s1"), pl.col("s2_dynamic").sum().alias("n_s2"))
          .to_pandas().set_index("source_class"))
    rg = reg.groupby("class").agg(n_sources=("source", "nunique"), found_per_cycle=("mean_found", "sum"), mean_found_per_source=("mean_found", "mean"))
    t = by.join(rg).reindex(PS.CLASS_ORDER)
    t["label"] = [PS.CLASS_LABEL[k] for k in t.index]
    t.index.name = "source_class"
    rows = [r"\begin{tabular}{lrrrrrcc}", r"\toprule",
            r"Source class & Sources & Dealers & Listings & Active & Found/cycle & S1 & S2 \\", r"\midrule"]
    for k, r in t.iterrows():
        s1 = "Yes" if r.in_s1 else "No"
        s2 = ("Yes" if r.n_sources_in_s2 == r.n_sources else f"Yes ({int(r.n_sources_in_s2)} of {int(r.n_sources)})") if r.in_s2 else "No"
        rows.append(f"{CLASS_SHORT[k]} & {fmt_int(r.n_sources)} & {fmt_int(r.n_dealers)} & {fmt_int(r.listings)} & {fmt_int(r.active_at_snapshot)} & {fmt_int(r.found_per_cycle)} & {s1} & {s2} \\\\")
    rows += [r"\midrule",
             f"All sources & {fmt_int(t.n_sources.sum())} & {fmt_int(L['dealer_id'].n_unique())} & {fmt_int(t.listings.sum())} & {fmt_int(t.active_at_snapshot.sum())} & {fmt_int(t.found_per_cycle.sum())} & & \\\\",
             r"\bottomrule", r"\end{tabular}"]
    write_tex(C.TABLES / "tab01_sources.tex", rows)
    return t.reset_index()


def tab02(sc: pd.DataFrame) -> pd.DataFrame:
    sc = sc.copy()
    sc["dropped"] = (-sc.listings_remaining.diff()).fillna(0).astype(int)
    rows = [r"\begin{tabular}{lrr}", r"\toprule", r"Step & Remaining & Dropped \\", r"\midrule"]
    marks = {"− make or model missing": None, "− uid identity collision": "S1 (all sources)",
             "− sampled portal (autotrader/cargurus) or private marketplace (kijiji)": "S1 hedonic",
             "− oscillating price path": "S2L", "− left-censored day-1 stock": "S2 dynamic"}
    for _, r in sc.iterrows():
        step = r.step.replace("−", "$-$").replace("[", "[").replace("$", r"\$").replace(r"\$-\$", "$-$")
        tag = marks.get(r.step)
        step = step if tag is None else f"{step} \\hfill $\\Rightarrow$ {tag}"
        rows.append(f"{step} & {fmt_int(r.listings_remaining)} & {fmt_int(r.dropped) if r.dropped else '--'} \\\\")
    rows += [r"\bottomrule", r"\end{tabular}"]
    write_tex(C.TABLES / "tab02_sample_construction.tex", rows)
    return sc


def tab03(L: pl.DataFrame) -> tuple[dict, pd.DataFrame]:
    S1 = L.filter(pl.col("s1_hedonic")); S2 = L.filter(pl.col("s2_dynamic"))

    def block(D: pl.DataFrame, dyn_completed: pl.DataFrame) -> dict:
        n = D.height
        d = {"N": n, "entry_price_mean": D["entry_price"].mean(), "entry_price_median": D["entry_price"].median(),
             "mileage_mean_km": D["mileage_km"].mean(), "mileage_median_km": D["mileage_km"].median(),
             "age_mean_years": D["age_years_entry"].mean(), "age_median_years": D["age_years_entry"].median(),
             "dealers": D["dealer_id"].n_unique(), "active_share": D["active"].mean(),
             "powertrain_share": {k: D.filter(pl.col("powertrain") == k).height / n for k in PS.POWERTRAIN_ORDER},
             "body_share": {k: D.filter(pl.col("body_filled") == k).height / n for k in PS.BODY_ORDER},
             "body_missing_share": D["body_filled"].null_count() / n,
             "duration_median_completed": dyn_completed["duration_days"].median(), "n_completed": dyn_completed.height,
             "share_with_cut": (D["n_cuts"] >= 1).mean(),
             "first_cut_pct_median": -100 * D.filter(pl.col("n_cuts") >= 1)["first_cut_pct"].median()}
        return d
    s1 = block(S1, S1.filter(pl.col("exit_confirmed") & ~pl.col("left_censored")))
    s2 = block(S2, S2.filter(pl.col("exit_confirmed")))
    top_bodies = sorted(PS.BODY_ORDER, key=lambda k: -s1["body_share"][k])[:5]
    # by region (S1), with exit-within-30d among S2 entrants observed long enough to confirm an exit at day 30
    horizon = (WITHIN_DAYS + C.EXIT_CONFIRM_DAYS) * 86400
    S2e = S2.filter(pl.col("entry_ts") <= C.SNAPSHOT_TS - horizon).with_columns(
        (pl.col("exit_confirmed") & (pl.col("duration_days") <= WITHIN_DAYS)).alias("exit30"))
    rs1 = S1.group_by(pl.col("region").fill_null("Region missing")).agg(pl.len().alias("N"), pl.col("entry_price").median().alias("median_price"),
                                                                        pl.col("mileage_km").median().alias("median_km"), pl.col("age_years_entry").median().alias("median_age"))
    rs2 = S2e.group_by(pl.col("region").fill_null("Region missing")).agg(pl.len().alias("n_s2_entrants_37d"), pl.col("exit30").mean().alias("share_exit_30d"))
    R = rs1.join(rs2, on="region", how="left").to_pandas().set_index("region")
    order = [r for r in PS.REGIONS if r in R.index] + [r for r in R.index if r not in PS.REGIONS]
    R = R.reindex(order).reset_index()
    tot = {"N": S1.height, "median_price": S1["entry_price"].median(), "median_km": S1["mileage_km"].median(), "median_age": S1["age_years_entry"].median(),
           "n_s2_entrants_37d": S2e.height, "share_exit_30d": S2e["exit30"].mean()}

    rows = [r"\begin{tabular}{lrr}", r"\toprule", r" & S1 hedonic & S2 dynamic \\", r"\midrule",
            f"Listings (N) & {fmt_int(s1['N'])} & {fmt_int(s2['N'])} \\\\",
            f"Dealers & {fmt_int(s1['dealers'])} & {fmt_int(s2['dealers'])} \\\\",
            f"Entry asking price, mean (\\$) & {fmt_int(s1['entry_price_mean'])} & {fmt_int(s2['entry_price_mean'])} \\\\",
            f"Entry asking price, median (\\$) & {fmt_int(s1['entry_price_median'])} & {fmt_int(s2['entry_price_median'])} \\\\",
            f"Mileage, mean (km) & {fmt_int(s1['mileage_mean_km'])} & {fmt_int(s2['mileage_mean_km'])} \\\\",
            f"Vehicle age at entry, mean (years) & {fmt1(s1['age_mean_years'])} & {fmt1(s2['age_mean_years'])} \\\\",
            r"\addlinespace", r"\multicolumn{3}{l}{\textit{Powertrain (\% of listings)}} \\"]
    for k in PS.POWERTRAIN_ORDER:
        rows.append(f"\\quad {k if k != 'unknown' else 'Unknown'} & {fmt_pct(s1['powertrain_share'][k])} & {fmt_pct(s2['powertrain_share'][k])} \\\\")
    rows += [r"\addlinespace", r"\multicolumn{3}{l}{\textit{Body type, top 5 (\% of listings)}} \\"]
    for k in top_bodies:
        rows.append(f"\\quad {k} & {fmt_pct(s1['body_share'][k])} & {fmt_pct(s2['body_share'][k])} \\\\")
    rows += [r"\addlinespace", r"\multicolumn{3}{l}{\textit{Listing dynamics}} \\",
             f"\\quad Completed spells (confirmed exits) & {fmt_int(s1['n_completed'])} & {fmt_int(s2['n_completed'])} \\\\",
             f"\\quad Duration of completed spells, median (days) & {fmt1(s1['duration_median_completed'])} & {fmt1(s2['duration_median_completed'])} \\\\",
             f"\\quad Listings with $\\geq 1$ sustained markdown (\\%) & {fmt_pct(s1['share_with_cut'])} & {fmt_pct(s2['share_with_cut'])} \\\\",
             f"\\quad First markdown, median (\\% of entry price) & {fmt1(s1['first_cut_pct_median'])} & {fmt1(s2['first_cut_pct_median'])} \\\\",
             r"\bottomrule", r"\end{tabular}", "",
             r"\begin{tabular}{lrrrrr}", r"\toprule",
             r"Region & N (S1) & Median price (\$) & Median km & Median age (y) & Exit $\leq$ 30 d (\%) \\", r"\midrule"]
    for _, r in R.iterrows():
        rows.append(f"{r.region} & {fmt_int(r.N)} & {fmt_int(r.median_price)} & {fmt_int(r.median_km)} & {fmt1(r.median_age)} & {fmt_pct(r.share_exit_30d)} \\\\")
    rows += [r"\midrule",
             f"All regions & {fmt_int(tot['N'])} & {fmt_int(tot['median_price'])} & {fmt_int(tot['median_km'])} & {fmt1(tot['median_age'])} & {fmt_pct(tot['share_exit_30d'])} \\\\",
             r"\bottomrule", r"\end{tabular}"]
    write_tex(C.TABLES / "tab03_descriptives.tex", rows)
    info = {"S1": s1, "S2": s2, "top_bodies_S1": top_bodies, "region_total": tot,
            "exit30_definition": f"S2 listings with entry_ts <= snapshot - {WITHIN_DAYS + C.EXIT_CONFIRM_DAYS} days (so that an exit at day {WITHIN_DAYS} is confirmable); "
                                 f"exit within {WITHIN_DAYS} d = exit_confirmed and duration_days <= {WITHIN_DAYS}",
            "S1_duration_note": "S1 completed spells exclude left-censored day-1 stock (duration unknown) but include all S1 sources"}
    return info, R


# ----------------------------------------------------------------------------- main
def run() -> dict:
    PS.setup()
    OUT.mkdir(parents=True, exist_ok=True)
    L, E, S, reg, sc = load()
    cov, cinfo = daily_coverage(L, S)
    rstock = region_stock(L)
    fig01(cov, rstock, cinfo)
    cov.to_csv(OUT / "fig01_daily_coverage.csv", index=False)
    rstock.to_csv(OUT / "fig01_region_stock.csv", index=False)

    ecdf, cuts, bc, dinfo = duration_and_cuts(L)
    fig03(ecdf, cuts, bc, dinfo)
    ecdf.iloc[::25].to_csv(OUT / "fig03_duration_ecdf.csv", index=False)
    bc.to_csv(OUT / "fig03_cut_share_by_class.csv", index=False)

    picks, rule = select_paths(L)
    paths = fig04(picks, E)
    (OUT / "fig04_selection.json").write_text(json.dumps(rule, indent=2, ensure_ascii=False), encoding="utf-8")
    paths.to_csv(OUT / "fig04_price_paths.csv", index=False)

    t1 = tab01(L, reg); t1.to_csv(OUT / "tab01_sources.csv", index=False)
    t2 = tab02(sc); t2.to_csv(OUT / "tab02_sample_construction.csv", index=False)
    t3info, R = tab03(L); R.to_csv(OUT / "tab03_by_region.csv", index=False)

    A = L.filter((pl.col("active") == 1) & pl.col("is_canonical") & (pl.col("kind") == "auto"))
    summary = {
        "snapshot_date": C.SNAPSHOT_DATE, "window_start_utc": "2026-08-12 12:30 UTC", "window_end_utc": "2026-11-01 06:49 UTC",
        "window_days": round((C.SNAPSHOT_TS - C.WINDOW_START_TS) / 86400, 1),
        "S0_n": L.height, "S1_n": t3info["S1"]["N"], "S2_n": t3info["S2"]["N"], "S2L_n": dinfo["S2L_n"],
        "active_canonical_auto_at_snapshot": A.height, "active_canonical_region_missing": int(A["region"].null_count()),
        "active_canonical_top_regions": rstock.head(5).to_dict(orient="records"),
        "sources_total": int(reg["source"].nunique()), "dealers_total": int(L["dealer_id"].n_unique()),
        "coverage": cinfo, "dynamics": dinfo, "fig04": rule, "tab01": t1.to_dict(orient="records"), "tab03": t3info,
        "tab03_by_region": R.to_dict(orient="records"),
        "notes": ["exit = observed disappearance from the source confirmed by >= 7 days of absence before the snapshot; it is NOT a sale",
                  "day index 0 = 2026-08-12 UTC (partial day, crawl started 12:30 UTC); day 81 = 2026-11-01 (partial, snapshot 06:49 UTC)",
                  "crawl cycles per day = sync_log rows / distinct sources logged that day (one row per source per cycle)",
                  "S2 contains no left-censored listing by construction; the left-censored share refers to S2L"]}
    (OUT / "summary.json").write_text(json.dumps(summary, indent=2, ensure_ascii=False, default=float), encoding="utf-8")
    return summary