#!/usr/bin/env python3
"""I-S3 Cyber Blast Radius (CBR) model — core definitions and reference calculations.

Definitions (all capacities in MW unless noted):

  C_total        = N * c / 1000                      total DER capacity
  CBR_raw(unit)  = capacity reachable from one compromised control unit
                   (device / vendor cloud / aggregator domain / region)
  CBR_eff        = reach * min(CBR_raw, A_max) * (1 - rho * (1 - e))
                   reach : remote-control reach fraction
                   A_max : maximum remote authority per domain (MW)
                   e     : command envelope (max ΔP per command as fraction of rating)
                   rho   : local enforcement reliability (fraction of devices that
                           actually reject out-of-envelope commands)
  x              = s * CBR_eff / P_ref               normalised synchronous disturbance
                   s     : synchronisation factor (fraction acting inside the critical window)
                   P_ref : reference reserve (FCR / primary reserve of the synchronous area)
  Severity(x)    = 1 - exp(-(x / x0)^k)              P(system-level failure | compromise)
  SystemicRisk   = p * Severity(x)                   annual probability of a systemic event
  Exposure       = p * CBR_eff * T_rec               expected MWh/yr under adversary control
  Concentration  = HHI of controllable capacity across control units (0..1)

The same functions are implemented in simulator/cbr_core.js; results/parity_cases.json is
generated here and checked by simulator/parity_test.js.

Run:  python3 model/cbr_model.py
Writes results/*.csv and results/*.json.
"""

from __future__ import annotations

import csv
import json
import math
from pathlib import Path

ROOT = Path(__file__).resolve().parents[1]
DATA = ROOT / "data"
RESULTS = ROOT / "results"


# ------------------------------------------------------------------ assumptions

def load_assumptions(path: Path = DATA / "assumptions.csv") -> dict[str, dict]:
    out: dict[str, dict] = {}
    with path.open(encoding="utf-8") as f:
        for row in csv.DictReader(f):
            out[row["key"]] = {
                "low": float(row["low"]),
                "base": float(row["base"]),
                "high": float(row["high"]),
                "unit": row["unit"],
                "label_ja": row["label_ja"],
                "label_en": row["label_en"],
                "source_id": row["source_id"],
            }
    return out


def base_params(assump: dict[str, dict] | None = None, level: str = "base") -> dict[str, float]:
    assump = assump or load_assumptions()
    return {k: v[level] for k, v in assump.items()}


# ------------------------------------------------------------------ core functions

def total_capacity_mw(n_devices: float, device_kw: float) -> float:
    return n_devices * device_kw / 1000.0


def cbr_raw_mw(c_total_mw: float, n_units: float) -> float:
    """Capacity behind one control unit when the fleet is split into n_units equal units."""
    if n_units <= 0:
        raise ValueError("n_units must be > 0")
    return c_total_mw / n_units


def cbr_effective_mw(
    cbr_raw: float,
    reach: float,
    authority_cap_mw: float,
    envelope_frac: float,
    local_enforce: float,
    local_autonomy: bool = True,
) -> float:
    """Capacity an attacker can actually move after one compromise of that unit."""
    capped = min(cbr_raw, authority_cap_mw) if local_autonomy else cbr_raw
    rho = local_enforce if local_autonomy else 0.0
    e = envelope_frac if local_autonomy else 1.0
    return reach * capped * (1.0 - rho * (1.0 - e))


def sync_factor_from_ramp(ramp_frac_per_min: float, window_s: float, base_sync: float = 1.0) -> float:
    """Fraction of the blast radius that can change inside a critical window given a ramp limit."""
    frac = min(1.0, ramp_frac_per_min * window_s / 60.0)
    return base_sync * frac


def normalised_disturbance(cbr_eff: float, sync: float, fcr_mw: float, frr_mw: float = 0.0) -> float:
    """x = max(fast, slow): a fast (synchronised) change is measured against the primary reserve,
    a ramp-limited change against primary + secondary reserve."""
    fast = sync * cbr_eff / fcr_mw if fcr_mw > 0 else float("inf")
    slow = cbr_eff / (fcr_mw + frr_mw) if (fcr_mw + frr_mw) > 0 else float("inf")
    return max(fast, slow)


def severity(x: float, x0: float, k: float) -> float:
    """P(system-level failure | normalised synchronous disturbance x)."""
    if x <= 0:
        return 0.0
    return 1.0 - math.exp(-((x / x0) ** k))


def systemic_risk(p_compromise: float, sev: float) -> float:
    return p_compromise * sev


def recovery_exposure_mwh(p_compromise: float, cbr_eff: float, t_recovery_h: float) -> float:
    return p_compromise * cbr_eff * t_recovery_h


def hhi(shares: list[float]) -> float:
    s = sum(shares)
    if s <= 0:
        return 0.0
    return sum((v / s) ** 2 for v in shares)


def zipf_shares(n: int, alpha: float = 1.0) -> list[float]:
    w = [1.0 / (i ** alpha) for i in range(1, n + 1)]
    t = sum(w)
    return [v / t for v in w]


def vendor_shares(n_vendors: float, top_share: float) -> list[float]:
    """Shares with a given largest share; remainder spread geometrically over the others."""
    n = max(1, int(round(n_vendors)))
    if n == 1:
        return [1.0]
    top = min(max(top_share, 1.0 / n), 1.0)
    rest = 1.0 - top
    others = [rest * v for v in zipf_shares(n - 1, 0.8)]
    # no other vendor may exceed the largest share: clip and redistribute (water-filling)
    for _ in range(n):
        over = sum(v - top for v in others if v > top)
        if over <= 1e-12:
            break
        below = [i for i, v in enumerate(others) if v < top]
        others = [min(v, top) for v in others]
        if not below:
            break
        w = sum(others[i] for i in below)
        for i in below:
            others[i] += over * (others[i] / w if w > 0 else 1.0 / len(below))
    return [top] + others


def evaluate_architecture(
    p: dict[str, float],
    *,
    n_domains: float,
    local_autonomy: bool,
    envelope_frac: float | None = None,
    authority_cap_mw: float | None = None,
    p_compromise: float | None = None,
    ramp_limit: float | None = None,
    window_s: float = 30.0,
) -> dict[str, float]:
    """Evaluate one architecture. p = parameter dict (base values); keyword args override."""
    c_total = total_capacity_mw(p["n_devices"], p["device_kw"])
    e = p["envelope_frac"] if envelope_frac is None else envelope_frac
    a_max = p["authority_cap_mw"] if authority_cap_mw is None else authority_cap_mw
    pc = p["p_compromise"] if p_compromise is None else p_compromise
    ramp = p["ramp_limit_frac_per_min"] if ramp_limit is None else ramp_limit

    raw_domain = cbr_raw_mw(c_total, n_domains)
    eff_domain = cbr_effective_mw(raw_domain, p["remote_reach"], a_max, e, p["local_enforce"], local_autonomy)

    shares = vendor_shares(p["n_vendors"], p["vendor_top_share"])
    raw_vendor = c_total * shares[0]
    # Vendor channel (vendor cloud / OTA / firmware line): bypasses per-domain credentials, so the
    # authority cap does not bind. The device-side envelope survives only to the extent it is enforced
    # outside the compromised firmware (vendor_hw_enforce).
    rho_vendor = p["local_enforce"] * p.get("vendor_hw_enforce", 0.0) if local_autonomy else 0.0
    eff_vendor = cbr_effective_mw(raw_vendor, p["remote_reach"], float("inf"), e, rho_vendor, local_autonomy)
    p_vendor = pc * p.get("p_vendor_rel", 1.0)

    fcr = p["grid_demand_mw"] * p["fcr_frac"]
    frr = p["grid_demand_mw"] * p.get("frr_frac", 0.0)
    s = sync_factor_from_ramp(ramp if local_autonomy else 100.0, window_s, p["sync_factor"])
    x_domain = normalised_disturbance(eff_domain, s, fcr, frr)
    # The ramp limit is firmware (slew limit on the current reference); it does not survive a vendor
    # firmware compromise even when the floor clamp is enforced by an independent monitor → step change.
    s_vendor = p["sync_factor"]
    x_vendor = normalised_disturbance(eff_vendor, s_vendor, fcr, frr)
    sev_domain = severity(x_domain, p["severity_x0"], p["severity_k"])
    sev_vendor = severity(x_vendor, p["severity_x0"], p["severity_k"])
    # all vendors, each an independent common-cause channel
    vendor_channel_risk = 1.0
    vendor_channel_mw = 0.0
    for sh in shares:
        eff_i = cbr_effective_mw(c_total * sh, p["remote_reach"], float("inf"), e, rho_vendor, local_autonomy)
        sev_i = severity(normalised_disturbance(eff_i, s_vendor, fcr, frr), p["severity_x0"], p["severity_k"])
        vendor_channel_risk *= (1.0 - p_vendor * sev_i)
        vendor_channel_mw += p_vendor * eff_i
    vendor_channel_risk = 1.0 - vendor_channel_risk
    # control-plane channel: independent domain compromises plus a common-cause share β that hits all
    # domains at once (shared identity provider / cloud / CA). Under β the per-domain cap still binds per
    # domain, but all domains move together.
    beta = p.get("domain_beta", 0.0) if n_domains > 1 else 0.0
    eff_all = eff_domain * n_domains
    sev_all = severity(normalised_disturbance(eff_all, s, fcr, frr), p["severity_x0"], p["severity_k"])
    plane_risk = (1.0 - beta) * fleet_prob(pc * sev_domain, n_domains) + beta * pc * sev_all

    # concentration: capacity controllable per domain (equal split) combined with vendor shares
    conc_domain = hhi([1.0] * max(1, int(round(n_domains))))
    conc_vendor = hhi(shares)

    return {
        "c_total_mw": c_total,
        "cbr_raw_domain_mw": raw_domain,
        "cbr_eff_domain_mw": eff_domain,
        "cbr_raw_vendor_mw": raw_vendor,
        "cbr_eff_vendor_mw": eff_vendor,
        "fcr_mw": fcr,
        "sync_factor": s,
        "x_domain": x_domain,
        "x_vendor": x_vendor,
        "severity_domain": sev_domain,
        "severity_vendor": sev_vendor,
        "systemic_risk_domain": systemic_risk(pc, sev_domain),
        "systemic_risk_vendor": systemic_risk(p_vendor, sev_vendor),
        "exposure_domain_mwh": recovery_exposure_mwh(pc, eff_domain, p["t_recovery_h"]),
        "exposure_vendor_mwh": recovery_exposure_mwh(p_vendor, eff_vendor, p["t_recovery_h"]),
        # linear metric: expected capacity under hostile control per year (does not saturate; captures local harm)
        "expected_affected_mw_domain": pc * eff_domain,
        "expected_affected_mw_vendor": p_vendor * eff_vendor,
        "p_vendor": p_vendor,
        "affected_fraction_domain": eff_domain / c_total if c_total else 0.0,
        "concentration_domain": conc_domain,
        "concentration_vendor": conc_vendor,
        "p_compromise": pc,
        # fleet level: n_domains independent attack surfaces, each compromised with probability pc per year.
        # More domains = more (smaller) events. Linear expected MW at fleet level is independent of the split.
        "fleet_p_any_compromise": fleet_prob(pc, n_domains),
        "fleet_systemic_risk": plane_risk,
        "fleet_expected_affected_mw": n_domains * pc * eff_domain,
        "fleet_exposure_mwh": n_domains * pc * eff_domain * p["t_recovery_h"],
        # vendor channel over all vendors, and the total of both channels
        "fleet_systemic_risk_vendor_channel": vendor_channel_risk,
        "fleet_systemic_risk_total": 1.0 - (1.0 - plane_risk) * (1.0 - vendor_channel_risk),
        "cbr_eff_all_domains_mw": eff_all,
        "fleet_expected_affected_mw_total": n_domains * pc * eff_domain + vendor_channel_mw,
    }


def fleet_prob(p_unit: float, n_units: float) -> float:
    """P(at least one of n independent units has the event in a year)."""
    n = max(1.0, float(n_units))
    if p_unit <= 0:
        return 0.0
    if p_unit >= 1:
        return 1.0
    return 1.0 - (1.0 - p_unit) ** n


def evaluate_a_plus(p: dict[str, float], *, p_compromise: float | None = None) -> dict[str, float]:
    """Centralized plane with realistic server-side controls: aggregate authority cap per critical window
    (independent safety monitor, staged rollout, command quotas). An API-credential compromise is bounded by
    the cap; a full-backend compromise (probability p × p_backend_rel) bypasses it. Devices have no envelope."""
    pc = p["p_compromise"] if p_compromise is None else p_compromise
    capped = evaluate_architecture(p, n_domains=1, local_autonomy=True, envelope_frac=1.0, authority_cap_mw=p["server_cap_mw"], ramp_limit=100.0, p_compromise=pc)
    full = evaluate_architecture(p, n_domains=1, local_autonomy=False, envelope_frac=1.0, authority_cap_mw=1e9, p_compromise=pc * p["p_backend_rel"])
    out = dict(full)
    out["cbr_eff_domain_mw"] = capped["cbr_eff_domain_mw"]
    out["x_domain"] = capped["x_domain"]
    out["severity_domain"] = capped["severity_domain"]
    out["systemic_risk_domain"] = capped["systemic_risk_domain"]
    plane = 1.0 - (1.0 - capped["fleet_systemic_risk"]) * (1.0 - full["fleet_systemic_risk"])
    out["fleet_systemic_risk"] = plane
    out["fleet_systemic_risk_total"] = 1.0 - (1.0 - plane) * (1.0 - full["fleet_systemic_risk_vendor_channel"])
    out["fleet_expected_affected_mw"] = capped["fleet_expected_affected_mw"] + full["fleet_expected_affected_mw"]
    out["fleet_expected_affected_mw_total"] = out["fleet_expected_affected_mw"] + (full["fleet_expected_affected_mw_total"] - full["fleet_expected_affected_mw"])
    out["p_compromise"] = pc
    return out


# ------------------------------------------------------------------ reference tables

ARCHITECTURES = [
    # name, n_domains, local_autonomy, envelope, authority cap, p multiplier, description
    ("A0_centralized_baseline", 1, False, 1.0, 1e9, 1.0, "単一の全国制御面。認証強化なし、権限上限・エンベロープなし"),
    ("A_centralized_certified", 1, False, 1.0, 1e9, 1 / 3, "単一制御面＋機器認証強化（侵害確率1/3）。権限上限・エンベロープなし"),
    ("A_centralized_certified_10x", 1, False, 1.0, 1e9, 1 / 10, "単一制御面＋機器認証強化（侵害確率1/10）"),
    ("A+_centralized_server_cap", 1, "aplus", 1.0, None, 1 / 3, "単一制御面＋認証(p/3)＋サーバー側総量上限1,500 MW／窓・段階配信・独立セーフティモニタ（API侵害では上限が効く、バックエンド全体侵害(p×0.3)では効かない）"),
    ("H_hierarchical_10", 10, False, 1.0, 1e9, 1.0, "地域階層10ドメイン（送配電エリア相当）。ローカル拒否なし"),
    ("F_federated_100", 100, False, 1.0, 1e9, 1.0, "連邦型100ドメイン（配電系統単位）。ローカル拒否なし"),
    ("C_cellular_10000", 10000, False, 1.0, 1e9, 1.0, "セル型10,000ドメイン。ローカル拒否なし（分割だけで守る）"),
    ("B4_single_plane_envelope_no_cap", 1, True, None, 1e9, 1.0, "単一制御面＋機器側エンベロープ（ρ=0.95, e=0.3）。権限上限なし"),
    ("B1_distributed_envelope_no_cap", 100, True, None, 1e9, 1.0, "100ドメイン＋機器側エンベロープ。権限上限なし"),
    ("B_distributed_autonomous", 100, True, None, None, 1.0, "100ドメイン＋機器側エンベロープ＋ドメイン権限上限500 MW"),
    ("B3_distributed_certified", 100, True, None, None, 1 / 3, "B＋機器認証強化（侵害確率1/3）"),
    ("B5_distributed_hw_envelope", 100, True, None, None, 1.0, "B＋エンベロープをハードウェア強制（ベンダーFW侵害後も残る）", {"vendor_hw_enforce": 1.0}),
    ("B6_distributed_hw_vendor_cap12", 100, True, None, None, 1.0, "B5＋首位ベンダーのシェア上限12%", {"vendor_hw_enforce": 1.0, "vendor_top_share": 0.12, "n_vendors": 12}),
    ("B7_distributed_hw_vendor_cap5", 100, True, None, None, 1.0, "B5＋首位ベンダーのシェア上限5%", {"vendor_hw_enforce": 1.0, "vendor_top_share": 0.05, "n_vendors": 25}),
]


def architecture_table(p: dict[str, float]) -> list[dict]:
    rows = []
    for row in ARCHITECTURES:
        name, nd, auto, env, cap, pmul, desc = row[:7]
        overrides = row[7] if len(row) > 7 else {}
        if auto == "aplus":
            r = evaluate_a_plus({**p, **overrides}, p_compromise=p["p_compromise"] * pmul)
            rows.append({"architecture": name, "description": desc, "n_domains": nd, "local_autonomy": False, **r})
            continue
        r = evaluate_architecture(
            {**p, **overrides},
            n_domains=nd,
            local_autonomy=auto,
            envelope_frac=env,
            authority_cap_mw=cap,
            p_compromise=p["p_compromise"] * pmul,
        )
        rows.append({"architecture": name, "description": desc, "n_domains": nd, "local_autonomy": auto, **r})
    return rows


def cbr_decomposition(p: dict[str, float]) -> list[dict]:
    """CBR_device / vendor / cloud / aggregator / region for the base fleet, with and without autonomy."""
    c_total = total_capacity_mw(p["n_devices"], p["device_kw"])
    shares = vendor_shares(p["n_vendors"], p["vendor_top_share"])
    units = [
        ("CBR_device", p["device_kw"] / 1000.0, "1台"),
        ("CBR_aggregator", cbr_raw_mw(c_total, p["n_domains"]), f"{int(p['n_domains'])}ドメイン均等"),
        ("CBR_region", c_total / 10.0, "10地域均等（一般送配電事業者エリア相当）"),
        ("CBR_vendor", c_total * shares[0], f"最大ベンダー{shares[0]:.0%}"),
        ("CBR_cloud", c_total, "単一の全国制御面／共通OTA基盤"),
    ]
    rows = []
    for name, raw, note in units:
        eff_auto = cbr_effective_mw(raw, p["remote_reach"], p["authority_cap_mw"], p["envelope_frac"], p["local_enforce"], True)
        eff_none = cbr_effective_mw(raw, p["remote_reach"], p["authority_cap_mw"], p["envelope_frac"], p["local_enforce"], False)
        rows.append({"unit": name, "note": note, "cbr_raw_mw": raw, "cbr_eff_no_autonomy_mw": eff_none, "cbr_eff_with_autonomy_mw": eff_auto})
    return rows


def break_even_grid(p: dict[str, float]) -> list[dict]:
    """Does systemic risk fall when p drops by factor X but blast radius rises by factor Y?"""
    fcr = p["grid_demand_mw"] * p["fcr_frac"]
    frr = p["grid_demand_mw"] * p.get("frr_frac", 0.0)
    base_cbr = 10.0  # MW, the reference distributed blast radius
    rows = []
    for p_factor in (1, 2, 3, 5, 10, 30, 100, 1000):
        for cbr_mult in (1, 3, 10, 30, 100, 300, 1000):
            x_b = normalised_disturbance(base_cbr, p["sync_factor"], fcr, frr)
            x_a = normalised_disturbance(base_cbr * cbr_mult, p["sync_factor"], fcr, frr)
            r_b = systemic_risk(p["p_compromise"], severity(x_b, p["severity_x0"], p["severity_k"]))
            r_a = systemic_risk(p["p_compromise"] / p_factor, severity(x_a, p["severity_x0"], p["severity_k"]))
            lin_b = p["p_compromise"] * base_cbr
            lin_a = p["p_compromise"] / p_factor * base_cbr * cbr_mult
            rows.append({
                "p_reduction_factor": p_factor,
                "blast_radius_multiplier": cbr_mult,
                "cbr_a_mw": base_cbr * cbr_mult,
                "x_a": x_a,
                "risk_centralized": r_a,
                "risk_distributed": r_b,
                "risk_ratio_A_over_B": (r_a / r_b) if r_b > 0 else float("inf"),
                "centralized_safer_systemic": r_a < r_b,
                "expected_mw_centralized": lin_a,
                "expected_mw_distributed": lin_b,
                "centralized_safer_linear": lin_a < lin_b,
            })
    return rows


def required_p_reduction(p: dict[str, float], cbr_a_mw: float, cbr_b_mw: float) -> float:
    """Factor by which p must fall for architecture A (bigger CBR) to match B's systemic risk."""
    fcr = p["grid_demand_mw"] * p["fcr_frac"]
    frr = p["grid_demand_mw"] * p.get("frr_frac", 0.0)
    sa = severity(normalised_disturbance(cbr_a_mw, p["sync_factor"], fcr, frr), p["severity_x0"], p["severity_k"])
    sb = severity(normalised_disturbance(cbr_b_mw, p["sync_factor"], fcr, frr), p["severity_x0"], p["severity_k"])
    return sa / sb if sb > 0 else float("inf")


def extreme_example(p: dict[str, float]) -> dict:
    fcr = p["grid_demand_mw"] * p["fcr_frac"]
    frr = p["grid_demand_mw"] * p.get("frr_frac", 0.0)
    a = {"p": 0.001, "cbr_mw": 5000.0}
    b = {"p": 0.01, "cbr_mw": 10.0}
    for d in (a, b):
        d["x"] = normalised_disturbance(d["cbr_mw"], p["sync_factor"], fcr, frr)
        d["severity"] = severity(d["x"], p["severity_x0"], p["severity_k"])
        d["systemic_risk_per_year"] = systemic_risk(d["p"], d["severity"])
        d["exposure_mwh_per_year"] = recovery_exposure_mwh(d["p"], d["cbr_mw"], p["t_recovery_h"])
        d["expected_affected_mw_per_year"] = d["p"] * d["cbr_mw"]
    return {
        "fcr_mw": fcr,
        "architecture_A": a,
        "architecture_B": b,
        "risk_ratio_A_over_B": a["systemic_risk_per_year"] / b["systemic_risk_per_year"] if b["systemic_risk_per_year"] > 0 else float("inf"),
        "linear_ratio_A_over_B": a["expected_affected_mw_per_year"] / b["expected_affected_mw_per_year"],
        "p_reduction_needed_for_A_to_match_B": required_p_reduction(p, a["cbr_mw"], b["cbr_mw"]),
        "note": "数字は例示。p・CBR・系統パラメータは data/assumptions.csv と simulator で変更できる。",
    }


def severity_calibration(p: dict[str, float]) -> list[dict]:
    """Observed / literature points used to sanity-check the severity curve (see research/grid_sensitivity_calibration.md)."""
    pts = [
        ("FR-2023 two nuclear units trip", 2660, 3000, "no UFLS, nadir 49.88 Hz (SPE/DNV 2025)", 0.05),
        ("ENTSO-E reference incident", 3000, 3000, "design basis; FCR sized to contain", 0.2),
        ("Dabrowski et al. 4.5 GW", 4500, 3000, "frequency below 49 Hz → UFLS (Forescout 2025 citing)", 0.5),
        ("SPE/DNV critical 10 GW", 10000, 3000, "estimated cascading-outage threshold (SPE 2025)", 0.95),
        ("Hokkaido 2018 (3,090 MW system)", 1160, 150, "N-3 + line faults → UFLS 1,300 MW → blackout (OCCTO)", 1.0),
    ]
    rows = []
    for name, dp, pref, note, expected in pts:
        x = dp / pref
        rows.append({"case": name, "delta_p_mw": dp, "p_ref_mw": pref, "x": x, "model_severity": severity(x, p["severity_x0"], p["severity_k"]), "expected_order": expected, "note": note})
    return rows


def vendor_share_limits(p: dict[str, float]) -> list[dict]:
    """Largest vendor share that keeps a vendor-channel compromise within x<1 / x<1.5 (section 6 definition of x)."""
    c_total = total_capacity_mw(p["n_devices"], p["device_kw"])
    fcr = p["grid_demand_mw"] * p["fcr_frac"]
    frr = p["grid_demand_mw"] * p.get("frr_frac", 0.0)
    s_ramp = sync_factor_from_ramp(p["ramp_limit_frac_per_min"], 30.0, p["sync_factor"])
    cases = [
        ("software-only: floor clamp and ramp both lost with firmware", 0.0, p["sync_factor"]),
        ("independent-monitor floor clamp survives, ramp lost (firmware)", p["local_enforce"], p["sync_factor"]),
        ("hypothetical: floor clamp and ramp both survive", p["local_enforce"], s_ramp),
    ]
    rows = []
    for name, rho, s in cases:
        factor = p["remote_reach"] * (1.0 - rho * (1.0 - p["envelope_frac"]))
        out = {"case": name, "fcr_mw": fcr, "frr_mw": frr, "c_total_mw": c_total}
        for thr in (1.0, 1.5):
            # x = max(s*E/FCR, E/(FCR+FRR)) < thr  →  E < thr*min(FCR/s, FCR+FRR)
            e_max = thr * min(fcr / s if s > 0 else float("inf"), fcr + frr)
            out[f"max_vendor_share_for_x_below_{str(thr).replace('.', '_')}"] = e_max / (factor * c_total)
            out[f"max_vendor_reach_mw_for_x_below_{str(thr).replace('.', '_')}"] = e_max
        rows.append(out)
    return rows


def policy_kpis(p: dict[str, float]) -> list[dict]:
    """Candidate system-level KPIs with the value the base case implies for architectures A and B."""
    a = evaluate_architecture(p, n_domains=1, local_autonomy=False, envelope_frac=1.0, authority_cap_mw=1e9)
    b = evaluate_architecture(p, n_domains=100, local_autonomy=True)
    return [
        {"kpi": "Cyber Blast Radius (max, MW)", "A_centralized": a["cbr_eff_domain_mw"], "B_distributed": b["cbr_eff_domain_mw"], "definition": "単一侵害で同時操作できる最大容量"},
        {"kpi": "Maximum Remote Authority (MW per credential/domain)", "A_centralized": a["cbr_raw_domain_mw"], "B_distributed": min(b["cbr_raw_domain_mw"], p["authority_cap_mw"]), "definition": "1つの認証情報・制御ドメインに許す遠隔権限の上限"},
        {"kpi": "Failure Domain Capacity (MW)", "A_centralized": a["cbr_raw_domain_mw"], "B_distributed": b["cbr_raw_domain_mw"], "definition": "独立した故障ドメイン1つの容量"},
        {"kpi": "Blast radius / FCR (x)", "A_centralized": a["x_domain"], "B_distributed": b["x_domain"], "definition": "一次調整力に対する比。1.5超でUFLS級"},
        {"kpi": "Vendor concentration (HHI)", "A_centralized": a["concentration_vendor"], "B_distributed": b["concentration_vendor"], "definition": "設置容量シェアのHHI。1/HHI＝有効ベンダー数"},
        {"kpi": "Firmware-line remote reach (MW, largest line)", "A_centralized": a["cbr_eff_vendor_mw"], "B_distributed": b["cbr_eff_vendor_mw"], "definition": "同一鍵・同一OTAで届く最大容量（ベンダー経路の CBR）。ソフトのみの B では A と同じ"},
        {"kpi": "Autonomous survival (runs on last valid schedule without upper link)", "A_centralized": 1.0, "B_distributed": 1.0, "definition": "通信断時に最後の有効スケジュールで運転する機能。日本の出力制御 PCS は A・B 共通に持つ（差分ではない）"},
        {"kpi": "Local command rejection (floor clamp enforced on device)", "A_centralized": 0.0, "B_distributed": p["local_enforce"], "definition": "枠外の遠隔命令を機器側で拒否する割合 ρ。監査可能な属性（遠隔停止クラスの無効化・独立監視MCU・OTA分離）に分解して申告する"},
        {"kpi": "Expected affected capacity (MW/yr)", "A_centralized": a["expected_affected_mw_domain"], "B_distributed": b["expected_affected_mw_domain"], "definition": "p×CBR（飽和しない線形指標。局所被害も拾う）"},
        {"kpi": "Recovery exposure (MWh/yr)", "A_centralized": a["exposure_domain_mwh"], "B_distributed": b["exposure_domain_mwh"], "definition": "p×CBR×復旧時間"},
        {"kpi": "Systemic risk per compromised unit (events/yr)", "A_centralized": a["systemic_risk_domain"], "B_distributed": b["systemic_risk_domain"], "definition": "p×Severity(x)：1つの制御単位が侵害されたときの広域障害確率"},
        {"kpi": "Fleet systemic risk, control-plane channel (events/yr)", "A_centralized": a["fleet_systemic_risk"], "B_distributed": b["fleet_systemic_risk"], "definition": "1−(1−p·S)^n：ドメイン数ぶんの攻撃面を数えた制御面経路の広域障害確率"},
        {"kpi": "Fleet systemic risk, vendor channel (events/yr)", "A_centralized": a["fleet_systemic_risk_vendor_channel"], "B_distributed": b["fleet_systemic_risk_vendor_channel"], "definition": "ベンダークラウド／OTA／ファームウェア経路。権限上限を迂回し、ソフトのみのエンベロープは消える"},
        {"kpi": "Fleet systemic risk, total (events/yr)", "A_centralized": a["fleet_systemic_risk_total"], "B_distributed": b["fleet_systemic_risk_total"], "definition": "両経路の合算"},
        {"kpi": "Fleet P(any compromise)/yr", "A_centralized": a["fleet_p_any_compromise"], "B_distributed": b["fleet_p_any_compromise"], "definition": "1−(1−p)^n：年にどこかが侵害される確率。分散側は高い"},
        {"kpi": "Fleet expected affected capacity (MW/yr)", "A_centralized": a["fleet_expected_affected_mw"], "B_distributed": b["fleet_expected_affected_mw"], "definition": "n·p·CBR_eff：分割では変わらず、エンベロープ・到達率でしか下がらない"},
    ]


def parity_cases(p: dict[str, float]) -> list[dict]:
    """Cases the JS simulator must reproduce (checked by simulator/parity_test.js)."""
    cases = []
    grid = [
        dict(n_domains=1, local_autonomy=False, envelope_frac=1.0, authority_cap_mw=1e9),
        dict(n_domains=100, local_autonomy=True),
        dict(n_domains=10000, local_autonomy=True, envelope_frac=0.1, authority_cap_mw=50.0),
        dict(n_domains=1, local_autonomy=True, envelope_frac=0.3, authority_cap_mw=500.0, p_compromise=0.1),
        dict(n_domains=10, local_autonomy=False, p_compromise=0.001),
    ]
    for i, kw in enumerate(grid):
        r = evaluate_architecture(p, **kw)
        cases.append({"id": f"case{i+1}", "inputs": {**p, **kw}, "expected": r})
    hw = {**p, "vendor_hw_enforce": 1.0, "vendor_top_share": 0.12, "n_vendors": 12}
    cases.append({"id": "case_hw_cap12", "inputs": {**hw, "n_domains": 100, "local_autonomy": True}, "expected": evaluate_architecture(hw, n_domains=100, local_autonomy=True)})
    cases.append({"id": "case_aplus", "inputs": {**p, "aplus": True, "p_compromise": p["p_compromise"] / 3}, "expected": evaluate_a_plus(p, p_compromise=p["p_compromise"] / 3)})
    return cases


# ------------------------------------------------------------------ io

def write_csv(path: Path, rows: list[dict]) -> None:
    if not rows:
        return
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", encoding="utf-8", newline="") as f:
        w = csv.DictWriter(f, fieldnames=list(rows[0].keys()))
        w.writeheader()
        for r in rows:
            w.writerow({k: (f"{v:.6g}" if isinstance(v, float) else v) for k, v in r.items()})


def main() -> None:
    assump = load_assumptions()
    p = base_params(assump)
    RESULTS.mkdir(parents=True, exist_ok=True)

    write_csv(RESULTS / "cbr_decomposition.csv", cbr_decomposition(p))
    write_csv(RESULTS / "architecture_comparison.csv", architecture_table(p))
    write_csv(RESULTS / "break_even_grid.csv", break_even_grid(p))
    write_csv(RESULTS / "severity_calibration.csv", severity_calibration(p))
    write_csv(RESULTS / "policy_kpis.csv", policy_kpis(p))
    write_csv(RESULTS / "vendor_share_limit.csv", vendor_share_limits(p))
    (RESULTS / "extreme_example.json").write_text(json.dumps(extreme_example(p), ensure_ascii=False, indent=2), encoding="utf-8")
    (RESULTS / "parity_cases.json").write_text(json.dumps(parity_cases(p), ensure_ascii=False, indent=2), encoding="utf-8")

    # sensitivity of required p-reduction to blast radius for the base grid
    rows = []
    for cbr_a in (10, 30, 100, 300, 1000, 3000, 5000, 10000, 35000):
        rows.append({"cbr_a_mw": cbr_a, "cbr_b_mw": 10, "p_reduction_needed": required_p_reduction(p, cbr_a, 10)})
    write_csv(RESULTS / "required_p_reduction.csv", rows)

    # how much the conclusion depends on the severity exponent k (k=1 ≈ linear expected-MW criterion)
    rows = []
    for k in (1.0, 1.5, 2.0, 3.0, 4.0):
        pk = {**p, "severity_k": k}
        for cbr_a in (30, 100, 300, 1000, 3000, 10000):
            rows.append({"severity_k": k, "cbr_a_mw": cbr_a, "cbr_b_mw": 10, "p_reduction_needed": required_p_reduction(pk, cbr_a, 10)})
    write_csv(RESULTS / "required_p_reduction_by_k.csv", rows)

    # envelope / enforcement sweep for a single national control plane
    rows = []
    c_total = total_capacity_mw(p["n_devices"], p["device_kw"])
    for e in (0.05, 0.1, 0.2, 0.3, 0.5, 1.0):
        for rho in (0.8, 0.9, 0.95, 0.99, 0.999):
            eff = cbr_effective_mw(c_total, p["remote_reach"], 1e9, e, rho, True)
            rows.append({"envelope_frac": e, "local_enforce": rho, "cbr_eff_single_plane_mw": eff, "x": normalised_disturbance(eff, p["sync_factor"], p["grid_demand_mw"] * p["fcr_frac"], p["grid_demand_mw"] * p.get("frr_frac", 0.0))})
    write_csv(RESULTS / "envelope_sweep.csv", rows)

    ex = extreme_example(p)
    print(f"C_total = {c_total:,.0f} MW; FCR = {ex['fcr_mw']:,.0f} MW")
    print(f"Extreme example: risk A/B = {ex['risk_ratio_A_over_B']:.2f}; p-reduction needed = {ex['p_reduction_needed_for_A_to_match_B']:.0f}x")
    for r in architecture_table(p):
        print(f"{r['architecture']:<32} CBR dom={r['cbr_eff_domain_mw']:>9,.0f} vend={r['cbr_eff_vendor_mw']:>9,.0f} MW  x_d={r['x_domain']:>6.2f} x_v={r['x_vendor']:>6.2f}  fleet: plane={r['fleet_systemic_risk']:.1e} vendor={r['fleet_systemic_risk_vendor_channel']:.1e} total={r['fleet_systemic_risk_total']:.1e}  E[MW]={r['fleet_expected_affected_mw_total']:.0f}")


if __name__ == "__main__":
    main()