"""Frequency, severity and supplier common-shock model for Japan's 2026 breach notices.
Run: python3 japan-breach-risk-model.py   (Python 3.9 or later, standard library only, fixed seed)
Reads the two CSVs that sit beside this file and writes japan-breach-risk-model-results.json there.
Source: https://secopsnomad.si/blog/japan-breaches-how-timeline-risk-model/
Records are counted in each notice's own unit (records, accounts or people); nothing here is money.
"""
import csv
import datetime as dt
import json
import math
import random
import statistics as st
from pathlib import Path

HERE = Path(__file__).resolve().parent
NOTICES = HERE / "japan-breach-notices-2026-08-10-to-2026-10-25.csv"
YEAR = HERE / "japan-breach-disclosures-100k-2025-09-to-2026-10.csv"
OUT = HERE / "japan-breach-risk-model-results.json"
SEED = 20261025
START, END = dt.date(2026, 8, 10), dt.date(2026, 10, 25)
DAYS = (END - START).days + 1  # 77 days, eleven Monday weeks
WEEKS = DAYS // 7
THRESHOLD = 100_000          # tail threshold, the same line as the 12-month list
CAP = 122_650_000            # illustrative cap; record/account counts can exceed the population
LDA_YEARS = 40_000
SHOCK_YEARS = 200_000
BOOTSTRAP = 2_000


def read(path):
    with open(path, encoding="utf-8-sig", newline="") as f:  # the CSVs carry a BOM so Excel reads them as UTF-8
        return list(csv.DictReader(f))


# ---------- distributions (stdlib has no binomial or Poisson sampler in 3.11) ----------
def binomial(rng, n, p):
    if n == 0 or p <= 0:
        return 0
    if p > 0.5:
        return n - binomial(rng, n, 1 - p)
    u, k = rng.random(), 0
    pk = math.exp(n * math.log1p(-p))
    cdf, ratio = pk, p / (1 - p)
    while u > cdf and k < n:
        pk *= (n - k) / (k + 1) * ratio
        k += 1
        cdf += pk
    return k


def poisson(rng, lam):
    if lam > 500:  # split so exp(-lam) cannot underflow
        half = lam / 2
        return poisson(rng, half) + poisson(rng, half)
    u, k = rng.random(), 0
    pk = math.exp(-lam)
    cdf = pk
    while u > cdf:
        k += 1
        pk *= lam / k
        cdf += pk
    return k


def chi2_sf(x, dof):
    """Upper tail of the chi-squared distribution: regularised upper incomplete gamma Q(dof/2, x/2)."""
    a, z = dof / 2, x / 2
    if z < a + 1:  # series for P, then Q = 1 - P
        term = total = 1 / a
        for i in range(1, 500):
            term *= z / (a + i)
            total += term
        return 1 - total * math.exp(-z + a * math.log(z) - math.lgamma(a))
    b, c, d = z + 1 - a, 1e300, 1 / (z + 1 - a)  # continued fraction for Q
    h = d
    for i in range(1, 500):
        an = -i * (i - a)
        b += 2
        d = 1 / (an * d + b)
        c = b + an / c
        h *= d * c
    return h * math.exp(-z + a * math.log(z) - math.lgamma(a))


def pct(sorted_values, q):
    """Nearest-rank percentile of an already sorted list."""
    return sorted_values[min(len(sorted_values) - 1, max(0, math.ceil(q * len(sorted_values)) - 1))]


# ---------- 1. frequency: weekly counts, dispersion, trend ----------
def attack(r):
    return r["category"] == "attack_personal_data"


def weekly(rows, keep):
    counts = [0] * WEEKS
    for r in rows:
        i = (dt.date.fromisoformat(r["notice_date"]) - START).days // 7
        if 0 <= i < WEEKS and keep(r):
            counts[i] += 1
    return counts


def poisson_trend_dispersion(y):
    """Fit log(mu) = a + b*week by Newton's method; return the Pearson dispersion of the residuals."""
    a, b = math.log(st.mean(y)), 0.0
    for _ in range(50):
        mu = [math.exp(a + b * t) for t in range(len(y))]
        g0 = sum(yi - m for yi, m in zip(y, mu))
        g1 = sum(t * (yi - m) for t, (yi, m) in enumerate(zip(y, mu)))
        h00 = sum(mu)
        h01 = sum(t * m for t, m in enumerate(mu))
        h11 = sum(t * t * m for t, m in enumerate(mu))
        det = h00 * h11 - h01 * h01
        a += (h11 * g0 - h01 * g1) / det
        b += (h00 * g1 - h01 * g0) / det
    mu = [math.exp(a + b * t) for t in range(len(y))]
    pearson = sum((yi - m) ** 2 / m for yi, m in zip(y, mu))
    return {"weeklyGrowth": round(math.exp(b) - 1, 3), "residualDispersion": round(pearson / (len(y) - 2), 2),
            "pValue": round(chi2_sf(pearson, len(y) - 2), 3)}


def frequency(rows):
    out = {}
    for name, keep in [("attackNotices", attack),
                       ("ownSystem", lambda r: attack(r) and not r["supplier_event_en"]),
                       ("customerNotices", lambda r: attack(r) and r["supplier_event_en"]),
                       ("attackNoticesWithFigure", lambda r: attack(r) and r["figure_value"])]:
        y = weekly(rows, keep)
        m, v = st.mean(y), st.variance(y)
        d = v / m
        out[name] = {"weeks": y, "mean": round(m, 2), "variance": round(v, 2), "dispersion": round(d, 2),
                     "poissonPValue": round(chi2_sf((len(y) - 1) * d, len(y) - 1), 3),
                     "negBinomialSize": round(m * m / (v - m), 2) if v > m else None,
                     "trend": poisson_trend_dispersion(y)}
    return out


# ---------- 2. severity: Hill tail index, bootstrap, lognormal, cross-check ----------
def hill(x_desc, k):
    return k / sum(math.log(x_desc[i] / x_desc[k]) for i in range(k))


def hill_above(values):
    x = sorted(values, reverse=True)
    k = sum(1 for v in x if v > THRESHOLD)
    return hill(x, k) if 1 < k < len(x) else None


def severity(rows, year_rows):
    figs = [r for r in rows if attack(r) and r["figure_value"]]
    x = sorted((int(r["figure_value"]) for r in figs), reverse=True)
    logs = [math.log(v) for v in x]
    k = sum(1 for v in x if v > THRESHOLD)
    alpha = hill(x, k)
    total = sum(x)
    boot_rng = random.Random(SEED + 1)  # its own stream, so the simulations below do not depend on it
    boot = sorted(a for a in (hill_above([boot_rng.choice(x) for _ in x]) for _ in range(BOOTSTRAP)) if a)
    no_idcf = hill_above([int(r["figure_value"]) for r in figs if r["supplier_event_en"] != "IDC Frontier"])
    # Cross-check on the 12-month list: Pareto MLE above the same threshold (all 38 are >= 100,000).
    y = [int(r["figure_value"]) for r in year_rows]
    alpha_year = len(y) / sum(math.log(v / THRESHOLD) for v in y)
    backtest = [{"level": level, "yearExpected": round(len(y) * (THRESHOLD / level) ** alpha, 1),
                 "yearObserved": sum(1 for v in y if v >= level), "expected": round(k * (THRESHOLD / level) ** alpha, 1),
                 "observed": sum(1 for v in x if v >= level)} for level in (1_000_000, 10_000_000)]
    return {"n": len(x), "median": st.median(x), "total": total,
            "top3Share": round(sum(x[:3]) / total, 3), "top1Share": round(x[0] / total, 3),
            "lognormalMu": round(st.mean(logs), 3), "lognormalSigma": round(st.pstdev(logs), 3),
            "threshold": THRESHOLD, "tailCount": k, "hillAlpha": round(alpha, 3),
            # The Hill plot is not flat: the very largest notices alone give a lighter tail (a natural ceiling).
            "hillByK": [[j, x[j], round(hill(x, j), 2)] for j in (15, 20, 30, 40, k)],
            "bootstrap": {"resamples": BOOTSTRAP, "low": round(pct(boot, 0.025), 2), "high": round(pct(boot, 0.975), 2)},
            "hillWithoutIdcFrontierCustomers": round(no_idcf, 3),
            "yearListN": len(y), "yearListAlpha": round(alpha_year, 3), "backtest": backtest}, x


# ---------- 3. loss-distribution approach on reported records ----------
def lda(freq, x_desc, rng):
    n = len(x_desc)
    lam = n * 365 / DAYS                          # notices with a figure per year
    disp = freq["attackNoticesWithFigure"]["dispersion"]
    r = lam / (disp - 1)                          # negative binomial as a Poisson-gamma mixture
    k = sum(1 for v in x_desc if v > THRESHOLD)
    alpha = hill(x_desc, k)
    body = x_desc[k:]
    mu = st.mean(math.log(v) for v in x_desc)
    sigma = st.pstdev(math.log(v) for v in x_desc)
    models = {
        "empirical": lambda: rng.choice(x_desc),
        "lognormal": lambda: min(CAP, rng.lognormvariate(mu, sigma)),
        "spliced": lambda: (min(CAP, THRESHOLD * rng.random() ** (-1 / alpha))
                            if rng.random() < k / n else rng.choice(body)),
    }
    out = {"noticesPerYear": round(lam, 1), "dispersion": disp, "years": LDA_YEARS, "cap": CAP}
    for name, draw in models.items():
        totals = []
        for _ in range(LDA_YEARS):
            count = poisson(rng, rng.gammavariate(r, lam / r))
            totals.append(sum(draw() for _ in range(count)))
        totals.sort()
        out[name] = {"mean": round(st.mean(totals)), "p50": round(pct(totals, 0.5)), "p90": round(pct(totals, 0.9)),
                     "p995": round(pct(totals, 0.995)), "p999": round(pct(totals, 0.999))}
    # With a tail index below 1 the mean is set by the cap, not the data: E[min(X, C)] for the Pareto part.
    body_mean = st.mean(body)
    out["capSensitivity"] = [{"cap": c, "meanPerYear": round(lam * ((k / n) * (THRESHOLD + THRESHOLD * ((c / THRESHOLD) ** (1 - alpha) - 1) / (1 - alpha))
                                                              + (1 - k / n) * body_mean))}
                             for c in (53_230_000, CAP, 1_000_000_000)]
    return out


# ---------- 4. supplier common shock: same mean, different concentration ----------
SHOCK = {"firms": 1000, "ownProbability": 0.03, "supplierProbability": 0.02, "hitRate": 0.5}
SCENARIOS = [
    ("own-supplier", None),                      # every firm has its own supplier: independent
    ("50-equal", [1 / 50] * 50),
    ("10-equal", [1 / 10] * 10),
    ("concentrated", [0.40, 0.20, 0.10, 0.10, 0.05, 0.05, 0.05, 0.05]),
    ("one-supplier", [1.0]),
]


def shock(rng):
    # Sum firm-level own/supplier impacts; the same firm can contribute twice (not unique firms).
    m, p0, q, h = (SHOCK[k] for k in ("firms", "ownProbability", "supplierProbability", "hitRate"))
    own = [binomial(rng, m, p0) for _ in range(SHOCK_YEARS)]   # common random numbers across scenarios
    out = {"parameters": SHOCK, "years": SHOCK_YEARS, "scenarios": []}
    for name, shares in SCENARIOS:
        sizes = [1] * m if shares is None else [round(s * m) for s in shares]
        hhi = 1 / m if shares is None else sum(s * s for s in shares)
        totals = []
        for year in range(SHOCK_YEARS):
            if shares is None:
                hit = binomial(rng, m, q * h)
            else:
                hit = sum(binomial(rng, n_j, h) for n_j in sizes if rng.random() < q)
            totals.append(own[year] + hit)
        totals.sort()
        var = m * p0 * (1 - p0) + sum(q * n * h * (1 - h) + q * (1 - q) * (n * h) ** 2 for n in sizes)
        mean = st.mean(totals)
        out["scenarios"].append({
            "name": name, "hhi": round(hhi, 3), "mean": round(mean, 1),
            "sd": round(st.pstdev(totals), 1), "sdFormula": round(math.sqrt(var), 1),
            "p50": pct(totals, 0.5), "p995": pct(totals, 0.995), "p999": pct(totals, 0.999),
            "pAtLeast100": round(sum(1 for t in totals if t >= 100) / SHOCK_YEARS, 4),
            "unexpectedMultiple999": round((pct(totals, 0.999) - mean) / mean, 1),
        })
    return out


# ---------- 5. what the notices show about suppliers, timing and data ----------
def notices_detail(rows):
    day = dt.date.fromisoformat
    attacks = [r for r in rows if attack(r)]
    customers = [r for r in attacks if r["supplier_event_en"]]  # attack notices from a supplier's customers
    cascades = {}
    for r in customers:
        cascades[r["supplier_event_en"]] = cascades.get(r["supplier_event_en"], 0) + 1
    cascades = sorted(cascades.items(), key=lambda kv: (-kv[1], kv[0]))
    bursts = []  # customers of one supplier publish together
    for name, n in cascades:
        if n < 5:
            continue
        mine = [r for r in customers if r["supplier_event_en"] == name]
        dates = sorted(day(r["notice_date"]) for r in mine)
        busiest = min(set(dates), key=lambda d: (-dates.count(d), d))
        bursts.append({"supplier": name, "notices": n, "busiestDay": busiest.isoformat(), "onBusiestDay": dates.count(busiest),
                       "spanDays": (dates[-1] - dates[0]).days + 1, "withFoundDate": sum(1 for r in mine if r["found"]),
                       "foundDates": sorted({r["found"] for r in mine if r["found"]})})
    def lag(a, b):
        v = sorted((day(r[b]) - day(r[a])).days for r in rows if r[a] and r[b])
        return {"n": len(v), "median": st.median(v), "p75": pct(v, 0.75), "max": v[-1]}
    coded = [r for r in rows if r["data_fields"]]
    fields = lambda r: set(r["data_fields"].split("; "))
    late = sorted(([r["organisation_en"], r["notice_date"], r["figure_date"], (day(r["figure_date"]) - day(r["notice_date"])).days]
                   for r in rows if r["figure_date"]), key=lambda x: -x[3])
    return {"customerNotices": len(customers), "customerShareOfAttacks": round(len(customers) / len(attacks), 3),
            "supplierIncidents": len(cascades), "cascades": [[name, n] for name, n in cascades],
            "largestFourShare": round(sum(n for _, n in cascades[:4]) / len(customers), 3),
            "bursts": bursts,
            "lags": {"foundToNotice": lag("found", "notice_date"), "accessToFound": lag("access_start", "found"),
                     "accessToNotice": lag("access_start", "notice_date")},
            "dataFieldsCoded": len(coded),
            "identityFields": sum(1 for r in coded if {"date of birth", "name"} <= fields(r) and {"phone", "email"} & fields(r)),
            "passwords": {c: sum(1 for r in rows if r["passwords"] == c) for c in ("hashed", "unknown-storage", "not-confirmed", "none")},
            "withSecondSource": sum(1 for r in rows if r["second_source"]),
            "lateFigures": late}


def main():
    rows, year_rows = read(NOTICES), read(YEAR)
    rng = random.Random(SEED)
    by = lambda key, subset: {c: sum(1 for r in subset if r[key] == c) for c in sorted({r[key] for r in subset})}
    attacks = [r for r in rows if attack(r)]
    own = [r for r in attacks if not r["supplier_event_en"]]
    counts = {"notices": len(rows), "byCategory": by("category", rows), "attackByCause": by("cause", attacks),
              "ownSystemByCause": by("cause", own), "ownSystemAttacks": len(own),
              "attackWithFigure": sum(1 for r in attacks if r["figure_value"]),
              "byMonth": {m: sum(r["notice_date"].startswith(m) for r in rows) for m in ("2026-08", "2026-09", "2026-10")}}
    freq = frequency(rows)
    sev, x = severity(rows, year_rows)
    result = {"seed": SEED, "window": [START.isoformat(), END.isoformat()], "counts": counts,
              "frequency": freq, "severity": sev, "lda": lda(freq, x, rng), "commonShock": shock(rng),
              "notices": notices_detail(rows)}
    OUT.write_text(json.dumps(result, ensure_ascii=False, indent=1) + "\n", encoding="utf-8")
    print(json.dumps(result, ensure_ascii=False, indent=1))


if __name__ == "__main__":
    main()