# SPDX-License-Identifier: MIT
"""Build the figures for YS-EN-0005 from results.json and write them into the
note between <!-- fig:name --> markers.

Usage: python build_figures.py   (after analyse.py)
"""
from __future__ import annotations

import json
import sys
from pathlib import Path

HERE = Path(__file__).resolve().parent
sys.path.insert(0, str(HERE.parent))
import svgchart as sc  # noqa: E402

PAGE = HERE.parent.parent / "ys-en-0005-what-a-40-year-mtbf-promises.html"
R = json.loads((HERE / "results.json").read_text())
MODELS = {m["model"]: m for m in R["models"]}
X24_DATASHEET_AFR = 0.0035        # Exos X24 data sheet, "Reliability Rating @ Full 24x7 Operation (AFR)"
MIN_BIN_DRIVE_YEARS = 3000        # don't plot age bins thinner than this


def pct(v, d=1):
    return f"{100 * v:.{d}f}%"


def name(m):
    """Readable model label: maker, model number, capacity."""
    mm = MODELS.get(m)
    tb = f'{mm["capacity_tb"]:g} TB' if mm else ""
    short = m.replace("HGST ", "HGST ").replace("TOSHIBA ", "Toshiba ").replace("WDC ", "WDC ")
    if short.startswith("ST"):
        short = "Seagate " + short
    return f"{short} ({tb})" if tb else short


# ---------------------------------------------------------------- figure 1
def fig_age():
    h = R["hazard_all"]
    e = h["edges"]
    keep = [i for i, dy in enumerate(h["drive_years"]) if dy >= MIN_BIN_DRIVE_YEARS]
    xs = [(e[i] + e[i + 1]) / 2 for i in keep]
    afr = [100 * h["afr"][i] for i in keep]
    lo = [100 * h["lo"][i] for i in keep]
    hi = [100 * h["hi"][i] for i in keep]
    mtbf_rate = 100 * R["mtbf"]["afr_rate"]
    fr = sc.Frame(720, 340, "Annualized failure rate of Backblaze hard drives by drive age, all models pooled")
    fr.x_scale(0, 10)
    fr.y_scale(0, 3.5)
    fr.grid_y([0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5], fmt=lambda v: f"{v:g}%")
    fr.grid_x([0, 2, 4, 6, 8, 10], fmt=lambda v: f"{v:g} yr", lines=False)
    fr.baseline()
    # 95% band
    band = " L".join(f"{fr.x(x):.1f},{fr.y(v):.1f}" for x, v in zip(xs, hi))
    band += " L" + " L".join(f"{fr.x(x):.1f},{fr.y(v):.1f}" for x, v in zip(reversed(xs), reversed(lo)))
    fr.out.append(f'<path class="area s1" d="M{band} Z"/>')
    fr.hline(mtbf_rate, cls="axis")
    fr.note(fr.x(0.15), fr.y(mtbf_rate) - 6, f"a 40-year MTBF: {mtbf_rate:.1f}% a year")
    fr.hline(100 * X24_DATASHEET_AFR, cls="axis")
    fr.note(fr.x(0.15), fr.y(100 * X24_DATASHEET_AFR) - 6, "Exos X24 datasheet: 0.35%")
    fr.titles("Drive age", "Failures per 100 drive-years")
    fr.cursor()
    series = [{"slot": 1, "name": "All drives", "values": afr, "label": False}]
    sc.draw_lines(fr, xs, series)
    sc.column_hits(fr, xs, series, lambda i, x: f"Age {e[keep[i]]:.2f}–{e[keep[i] + 1]:.2f} yr",
                   lambda v: f"{v:.2f}%")
    rows = [[f"{e[i]:.2f}–{e[i + 1]:.2f}", f"{h['drive_years'][i]:,.0f}", f"{h['failures'][i]:,}",
             f"{100 * h['afr'][i]:.2f}%", f"{100 * h['lo'][i]:.2f}–{100 * h['hi'][i]:.2f}%"] for i in keep]
    d = R["data"]
    return sc.figure(
        code="FIG. 1 — FAILURE RATE BY AGE",
        fid="fig-age",
        plot=fr.svg(),
        caption=(f"<b>Figure 1.</b> Failures per 100 drive-years in each quarter-year of drive age, for "
                 f"{d['drives']:,} hard drives of {d['models']} models (Backblaze, April 2013 to June 2026). "
                 f"The shaded band is the 95% Poisson interval. Ages with fewer than {MIN_BIN_DRIVE_YEARS:,} "
                 "drive-years are left out. Every model is pooled here, so the shape mixes ageing with the "
                 "changing mix of models at each age."),
        data_html=sc.table("Failure rate by age, all models", ["Age, yr", "Drive-years", "Failures", "Rate", "95% interval"], rows),
    )


# ---------------------------------------------------------------- figure 2
SHOWN = [("ST4000DM000", 1, "Seagate 4 TB"), ("HGST HUH721212ALN604", 3, "HGST 12 TB"),
         ("HGST HMS5C4040BLE640", 2, "HGST 4 TB")]


def fig_models():
    hb = R["hazard_by_model"]
    xs = [k + 0.5 for k in range(10)]
    series = []
    for m, slot, short in SHOWN:
        h = hb[m]
        vals = [100 * a if dy >= 300 else None for a, dy in zip(h["afr"], h["drive_years"])][: len(xs)]
        vals += [None] * (len(xs) - len(vals))
        fit = MODELS[m]
        series.append({"slot": slot, "name": f"{short} {m.split()[-1]} (β {fit['beta']:.2f})",
                       "end": short, "values": vals})
    last = max(i for s in series for i, v in enumerate(s["values"]) if v is not None)
    xs = xs[: last + 1]
    for s in series:
        s["values"] = s["values"][: last + 1]
    mtbf_rate = 100 * R["mtbf"]["afr_rate"]
    fr = sc.Frame(720, 340, "Failure rate by age for three drive models: one roughly constant, one rising, one falling")
    fr.x_scale(0, 10)
    fr.y_scale(0, 6)
    fr.grid_y([0, 1, 2, 3, 4, 5, 6], fmt=lambda v: f"{v:g}%")
    fr.grid_x([0, 2, 4, 6, 8, 10], fmt=lambda v: f"{v:g} yr", lines=False)
    fr.baseline()
    fr.hline(mtbf_rate)
    fr.note(fr.x(10), fr.y(mtbf_rate) + 16, f"40-year MTBF: {mtbf_rate:.1f}%", anchor="end")
    fr.titles("Drive age (one-year bins)", "Failures per 100 drive-years")
    fr.cursor()
    sc.draw_lines(fr, xs, series)
    sc.column_hits(fr, xs, series, lambda i, x: f"Age {x - 0.5:.0f}–{x + 0.5:.0f} yr", lambda v: f"{v:.2f}%")
    rows = []
    for m, slot, short in SHOWN:
        h = hb[m]
        for k, (a, f, dy) in enumerate(zip(h["afr"], h["failures"], h["drive_years"])):
            if dy >= 300:
                rows.append([f"{short} {m.split()[-1]}", f"{k}–{k + 1}", f"{dy:,.0f}", f"{f:,}", f"{100 * a:.2f}%"])
    return sc.figure(
        code="FIG. 2 — THREE MODELS, THREE SHAPES",
        fid="fig-models",
        plot=fr.svg(),
        legend_html=sc.legend(series),
        caption=("<b>Figure 2.</b> Failures per 100 drive-years by year of age for three of the largest models. "
                 f"Over the ages observed, the Seagate 4 TB ({MODELS['ST4000DM000']['drives']:,} drives) and the HGST 12 TB "
                 f"({MODELS['HGST HUH721212ALN604']['drives']:,}) both averaged about one failure per 40 drive-years: an MTBF of "
                 f"{MODELS['ST4000DM000']['mtbf_yr_observed']:.1f} and {MODELS['HGST HUH721212ALN604']['mtbf_yr_observed']:.1f} years. "
                 "The Seagate stays near that rate all its life. The HGST runs at about 0.5% for three years, then "
                 "fails eight to ten times as often. The HGST 4 TB is worst in its first year and improves after that. "
                 "Years with fewer than 300 drive-years are left out."),
        data_html=sc.table("Failure rate by year of age, three models", ["Model", "Age, yr", "Drive-years", "Failures", "Rate"], rows),
    )


# ---------------------------------------------------------------- figure 3
def fig_shapes():
    ms = sorted(R["models"], key=lambda m: m["beta"])
    row_h = 22
    left = 250
    fr = sc.Frame(720, 40 + row_h * len(ms) + 46, "Fitted Weibull shape for each drive model with its 95% interval",
                  margins=(left, 30, 30, 46))
    fr.x_scale(0.5, 5, log=True)
    fr.y_scale(0, 1)
    ticks = [0.5, 0.7, 1, 1.5, 2, 3, 5]
    fr.out.append('<g class="grid">')
    for t in ticks:
        fr.out.append(f'<line x1="{fr.x(t):.1f}" x2="{fr.x(t):.1f}" y1="{fr.mt}" y2="{fr.mt + fr.ph}"/>')
    fr.out.append("</g>")
    for t in ticks:
        fr.out.append(f'<text x="{fr.x(t):.1f}" y="{fr.mt + fr.ph + 18}" text-anchor="middle">{t:g}</text>')
    x1 = fr.x(1)
    fr.out.append(f'<line class="axis" x1="{x1:.1f}" x2="{x1:.1f}" y1="{fr.mt}" y2="{fr.mt + fr.ph}"/>')
    fr.note(x1 - 6, fr.mt - 10, "← failure rate falls with age", anchor="end")
    fr.note(x1 + 6, fr.mt - 10, "constant (exponential) · rises with age →")
    fr.titles("Weibull shape β (log scale)")
    rows = []
    for i, m in enumerate(ms):
        cy = fr.mt + row_h * (i + 0.5) + 4
        slot = 3 if m["beta_lo"] > 1 else (2 if m["beta_hi"] < 1 else 0)
        lab = name(m["model"])
        fr.out.append(f'<text class="ink" x="{fr.ml - 12}" y="{cy + 4:.1f}" text-anchor="end">{sc.esc(lab)}</text>')
        a, b = fr.x(max(0.5, m["beta_lo"])), fr.x(min(5, m["beta_hi"]))
        fr.out.append(f'<line class="s{slot}" x1="{a:.1f}" x2="{b:.1f}" y1="{cy:.1f}" y2="{cy:.1f}" stroke-width="2"/>')
        cx = fr.x(m["beta"])
        fr.out.append(f'<circle class="dot s{slot}" cx="{cx:.1f}" cy="{cy:.1f}" r="5"/>')
        info = (f'{slot}|β {m["beta"]:.2f} ({m["beta_lo"]:.2f}–{m["beta_hi"]:.2f})|'
                f'{m["drives"]:,} drives, {m["failures"]:,} failures, {pct(m["afr"], 2)} a year')
        sc.point_hit(fr, cx, cy, i, lab, info, r=11)
        rows.append([lab, f'{m["drives"]:,}', f'{m["failures"]:,}', f'{m["drive_years"]:,.0f}', pct(m["afr"], 2),
                     f'{m["max_age_yr"]:.1f}', f'{m["beta"]:.2f}', f'{m["beta_lo"]:.2f}–{m["beta_hi"]:.2f}',
                     f'{m["p_exponential"]:.1g}'])
    s = R["model_summary"]
    return sc.figure(
        code="FIG. 3 — IS THE RATE CONSTANT?",
        fid="fig-shapes",
        plot=fr.svg(),
        caption=(f"<b>Figure 3.</b> Weibull shape β fitted to each of the {s['n']} models with at least "
                 f"{R['inputs']['min_drives']:,} drives and {R['inputs']['min_fails']} failures, with 95% profile-likelihood "
                 f"intervals. β = 1 is a constant failure rate, the assumption behind every MTBF. "
                 f"For {s['beta_gt1_sig']} models (coral) the rate rises with age, for {s['beta_lt1_sig']} (violet) it falls, and "
                 f"for {s['beta_ci_covers_1']} (grey) a constant rate cannot be ruled out. A single β is itself a simplification; "
                 "Figure 2 shows shapes it cannot capture."),
        data_html=sc.table("Weibull fits by model (left-truncated, right-censored maximum likelihood)",
                           ["Model", "Drives", "Failures", "Drive-years", "Rate", "Oldest, yr", "β", "95% interval",
                            "p (β = 1)"], rows),
    )


# ---------------------------------------------------------------- figure 4
def fig_extrap(age="5"):
    pts = [m for m in R["models"] if age in m["cum_failed"]]
    fr = sc.Frame(720, 400, f"Share of drives failed by age {age}: predicted from the first year's rate versus observed",
                  margins=(58, 30, 30, 46))
    top = 20
    fr.x_scale(0, top)
    fr.y_scale(0, top)
    tk = [0, 5, 10, 15, 20]
    fr.grid_y(tk, fmt=lambda v: f"{v:g}%")
    fr.grid_x(tk, fmt=lambda v: f"{v:g}%")
    fr.baseline()
    fr.out.append(f'<line class="axis" x1="{fr.x(0):.1f}" y1="{fr.y(0):.1f}" x2="{fr.x(top):.1f}" y2="{fr.y(top):.1f}"/>')
    fr.note(fr.x(top) - 4, fr.y(top) + 16, "prediction = outcome", anchor="end")
    fr.titles(f"Predicted from the first year, as a constant rate: share failed by year {age}",
              f"Observed share failed by year {age}")
    rows = []
    for i, m in enumerate(sorted(pts, key=lambda m: m["cum_failed"][age]["first_year_rate_predicts"])):
        c = m["cum_failed"][age]
        p, o = 100 * c["first_year_rate_predicts"], 100 * c["observed_failed"]
        slot = 3 if o > 1.25 * p else (2 if o < 0.8 * p else 1)
        cx, cy = fr.x(min(p, top)), fr.y(min(o, top))
        fr.out.append(f'<circle class="dot s{slot}" cx="{cx:.1f}" cy="{cy:.1f}" r="5"/>')
        lab = name(m["model"])
        sc.point_hit(fr, cx, cy, i, lab, f"{slot}|Predicted {p:.1f}%, observed {o:.1f}%|ratio {o / p:.1f}×", r=10)
        rows.append([lab, f"{100 * m['first_year_afr']:.2f}%", f"{p:.1f}%", f"{o:.1f}%", f"{o / p:.2f}×"])
    under = sum(1 for r in rows if float(r[4][:-1]) > 1.25)
    over = sum(1 for r in rows if float(r[4][:-1]) < 0.8)
    return sc.figure(
        code="FIG. 4 — EXTRAPOLATING THE FIRST YEAR",
        fid="fig-extrap",
        plot=fr.svg(),
        caption=(f"<b>Figure 4.</b> For each model followed to age {age}: the share of drives that failed by then, against "
                 "the share a constant rate measured in the first year of life would predict. That is how a "
                 "short test turns into an MTBF. "
                 f"Of {len(rows)} models, {under} failed more than 1.25× as often as predicted (coral) and {over} "
                 "less than 0.8× as often (violet). Cumulative shares come from the piecewise-constant hazard in "
                 "quarter-year bins, which handles drives that joined or left partway."),
        data_html=sc.table(f"First-year extrapolation versus outcome at year {age}",
                           ["Model", "First-year rate", f"Predicted by yr {age}", f"Observed by yr {age}", "Observed ÷ predicted"],
                           rows),
    )


# ---------------------------------------------------------------- table: the 4077 under real shapes
LIKE = {"0.66": "HGST 4 TB HMS5C4040BLE640", "1": "Seagate 4 TB ST4000DM000", "1.45": "Seagate 12 TB ST12000NM0008",
        "2.18": "HGST 12 TB HUH721212ALN604", "4": "Seagate 3 TB ST3000DM001"}


def table_4077():
    w = R["mtbf"]["weibull_same_mtbf"]
    n = R["inputs"]["fleet"]
    rows = []
    for b, v in w.items():
        y1 = v["afr_year1"]
        rows.append([f"β {float(b):g}", LIKE[b], f'{v["fail_by_10"]:.0f}', f'{v["fail_by_40"]:.0f}',
                     (f"{100 * y1:.1f}%" if y1 >= 0.001 else f"{100 * y1:.3f}%" if y1 >= 1e-5 else "<0.001%"),
                     f'{100 * v["afr_year30"]:.1f}%', f'{v["mrl_at_30"]:.0f} yr'])
    return sc.table(f"{n} cores with a 40-year MTBF, unrepaired, under five Weibull shapes fitted to real drive models",
                    ["Shape", "Fitted to", "Failed by year 10", "Failed by year 40", "Fails in year 1",
                     "Fails in year 30, if still running", "Expected further life at 30"], rows, numeric_from=2)


def main() -> None:
    sc.inject(PAGE, "table4077", table_4077())
    for nm, fn in (("age", fig_age), ("models", fig_models), ("shapes", fig_shapes), ("extrap", fig_extrap)):
        sc.inject(PAGE, nm, fn())
    print("figures written to", PAGE.name)


if __name__ == "__main__":
    main()