#!/usr/bin/env python3
"""製品認証方式（A）と境界防護方式（B）の費用対効果モデル。

入力: data/assumptions.csv（low/base/high と種別 kind）
出力: results/*.csv, results/summary.json, figures/fig07_cost_comparison.svg

設計方針
- 架空の数字を実績値として使わない。kind が assumption / to-confirm の値は本文でも「仮定」と明示する。
- 二つの方式は同じ脅威を扱っていない。費用は「年あたり・1 サイトあたり」と「1 kW あたり」で並べ、
  防御の到達可否（coverage）は確率ではなく設計上の到達可否（0 / 1 / 条件付き）で別表にする。
- 再現: python3 model/cost_model.py  （numpy 不要、標準ライブラリのみ）
"""

from __future__ import annotations

import csv
import json
import math
from pathlib import Path

HERE = Path(__file__).resolve().parent
ROOT = HERE.parent
DATA = ROOT / "data" / "assumptions.csv"
RESULTS = ROOT / "results"
FIGURES = ROOT / "figures"
RESULTS.mkdir(exist_ok=True)
FIGURES.mkdir(exist_ok=True)


def load_assumptions() -> dict[str, dict]:
    out: dict[str, dict] = {}
    with DATA.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"],
                "kind": row["kind"],
                "note": row["note"],
                "source": row["source"],
            }
    return out


A = load_assumptions()


def p(key: str, case: str = "base") -> float:
    return A[key][case]


# ---------------------------------------------------------------- 年間費用（1 サイト）

SITE_CLASSES = [
    ("lv_res", "低圧・住宅（5 kW）", "site_kw_lv_res", "lv"),
    ("lv_biz", "低圧・事業用（49.5 kW）", "site_kw_lv_biz", "lv"),
    ("hv", "高圧（500 kW）", "site_kw_hv", "hv"),
    ("ehv", "特別高圧（5,000 kW）", "site_kw_ehv", "hv"),
]


def annual_cost_A(kw: float, case: str = "base", override: dict | None = None) -> dict:
    """製品認証方式の 1 サイト年間費用（事業者側に落ちる分）。

    (1) 認証費用の価格転嫁: (申請手数料 + 社内工数) / 出荷台数 を有効期間で年割り
    (2) 交換遅延の逸失電力: 故障率 × 追加停止日数 × kW × 設備利用率 × 24h × 単価
    """
    o = override or {}
    g = lambda k: o.get(k, p(k, case))
    per_unit = (g("cert_fee_model") + g("cert_internal_model")) / g("units_per_model")
    cert_year = per_unit / g("cert_validity_years")
    delay_loss = g("pcs_fail_rate_year") * g("replace_delay_days_A") * kw * g("capacity_factor") * 24 * g("value_kwh")
    return {"cert_passthrough_year": cert_year, "replace_delay_loss_year": delay_loss, "total_year": cert_year + delay_loss}


def annual_cost_B(kw: float, tier: str, case: str = "base", override: dict | None = None) -> dict:
    """境界防護方式の 1 サイト年間費用。機器費を更新周期で年割り + 運用費 + （短い）交換遅延。"""
    o = override or {}
    g = lambda k: o.get(k, p(k, case))
    capex = g(f"boundary_capex_{tier}")
    opex = g(f"boundary_opex_{tier}_year")
    cap_year = capex / g("boundary_life_years")
    delay_loss = g("pcs_fail_rate_year") * g("replace_delay_days_B") * kw * g("capacity_factor") * 24 * g("value_kwh")
    return {"boundary_capex_year": cap_year, "boundary_opex_year": opex, "replace_delay_loss_year": delay_loss, "total_year": cap_year + opex + delay_loss}


def cost_table() -> list[dict]:
    rows = []
    for cid, label, kwkey, tier in SITE_CLASSES:
        for case in ("low", "base", "high"):
            kw = p(kwkey, "base")
            a = annual_cost_A(kw, case)
            b = annual_cost_B(kw, tier, case)
            rows.append({
                "site_class": cid, "label": label, "case": case, "kw": kw,
                "A_cert_passthrough_year": round(a["cert_passthrough_year"]),
                "A_replace_delay_loss_year": round(a["replace_delay_loss_year"]),
                "A_total_year": round(a["total_year"]),
                "A_total_year_per_kw": round(a["total_year"] / kw, 1),
                "B_boundary_capex_year": round(b["boundary_capex_year"]),
                "B_boundary_opex_year": round(b["boundary_opex_year"]),
                "B_replace_delay_loss_year": round(b["replace_delay_loss_year"]),
                "B_total_year": round(b["total_year"]),
                "B_total_year_per_kw": round(b["total_year"] / kw, 1),
                "B_minus_A_year": round(b["total_year"] - a["total_year"]),
            })
    return rows


# ---------------------------------------------------------------- 損益分岐

def break_even_capex(kw: float, tier: str) -> float:
    """B の年間費用が A と等しくなる境界装置 capex（base の他条件）。"""
    a = annual_cost_A(kw)["total_year"]
    b0 = annual_cost_B(kw, tier, override={f"boundary_capex_{tier}": 0.0})["total_year"]
    return max(0.0, (a - b0) * p("boundary_life_years"))


def break_even_table() -> list[dict]:
    rows = []
    for cid, label, kwkey, tier in SITE_CLASSES:
        kw = p(kwkey, "base")
        rows.append({
            "site_class": cid, "label": label, "kw": kw,
            "A_total_year": round(annual_cost_A(kw)["total_year"]),
            "break_even_boundary_capex": round(break_even_capex(kw, tier)),
            "base_boundary_capex": p(f"boundary_capex_{tier}"),
            "note": "境界装置 capex がこの値以下なら、年間費用で B が A を下回る（他条件は base）",
        })
    # 交換遅延日数を振ったときの損益分岐（高圧）
    for days in (0, 7, 14, 30, 45, 60, 90, 120):
        kw = p("site_kw_hv")
        a = annual_cost_A(kw, override={"replace_delay_days_A": days})["total_year"]
        b0 = annual_cost_B(kw, "hv", override={"boundary_capex_hv": 0.0})["total_year"]
        rows.append({
            "site_class": "hv_sweep_delay", "label": f"高圧・交換遅延 {days} 日", "kw": kw,
            "A_total_year": round(a), "break_even_boundary_capex": round(max(0.0, (a - b0) * p("boundary_life_years"))),
            "base_boundary_capex": p("boundary_capex_hv"), "note": "replace_delay_days_A を振った感度",
        })
    return rows


# ---------------------------------------------------------------- 感度（トルネード）

SENS_KEYS = [
    "cert_fee_model", "cert_internal_model", "units_per_model", "cert_validity_years",
    "pcs_fail_rate_year", "replace_delay_days_A", "replace_delay_days_B", "capacity_factor", "value_kwh",
    "boundary_life_years",
]


def tornado(site: str) -> list[dict]:
    cid, label, kwkey, tier = next(s for s in SITE_CLASSES if s[0] == site)
    kw = p(kwkey)
    keys = SENS_KEYS + [f"boundary_capex_{tier}", f"boundary_opex_{tier}_year"]
    base_diff = annual_cost_B(kw, tier)["total_year"] - annual_cost_A(kw)["total_year"]
    rows = []
    for k in keys:
        vals = {}
        for case in ("low", "high"):
            o = {k: p(k, case)}
            vals[case] = annual_cost_B(kw, tier, override=o)["total_year"] - annual_cost_A(kw, override=o)["total_year"]
        rows.append({
            "site_class": cid, "param": k, "kind": A[k]["kind"], "base_diff_B_minus_A": round(base_diff),
            "diff_at_low": round(vals["low"]), "diff_at_high": round(vals["high"]),
            "swing": round(abs(vals["high"] - vals["low"])),
        })
    rows.sort(key=lambda r: -r["swing"])
    return rows


# ---------------------------------------------------------------- 到達可否（coverage）

# 値は確率ではない。設計上、その経路の不正指令が機器に「届くか」。
#   0 = 評価・遮断の対象外（届く）、1 = 設計上届かない、"cond" = 条件付き（本文で条件を書く）
COVERAGE = [
    ("T1", "機器一台の既定パスワード・露出した管理画面", "1(要件あり)", "1(露出しない)", "★1 は既定パスワード廃止・更新機能を要求。境界は管理画面を外に出さない"),
    ("T2", "メーカーのクラウド制御面の侵害（正規 API からの一斉指令）", "0(対象外)", "cond(許可リスト・書込遮断)", "IPA 回答: ラベルは設置・運用を示さない。境界側は書込系ファンクションを遮断できる構成なら届かない。ρ で割り引く"),
    ("T3", "更新配信サーバー・署名鍵の侵害（正規署名の不正更新）", "0(対象外)", "cond(段階適用・更新経路制限)", "鍵が正しければ機器側検証は通る。境界側は配信先・時期を制限できるが内容は検証できない"),
    ("T4", "正規管理者アカウントの乗っ取り（保守・遠隔監視）", "0(対象外)", "cond(最小権限・操作制限)", "認証は機能の有無を見る。境界側で操作種別と到達台数を制限する"),
    ("T5", "認証後に公開された未知の脆弱性の悪用", "cond(更新機能)", "1(到達しない)", "★1 の更新機能は公開後の修正適用を速める。境界内にあれば公開前後とも外から届かない"),
    ("T6", "第三国の基盤・踏み台を経由した遠隔攻撃", "0(対象外)", "cond(経路制御)", "製造国の選別は経路を限定しない。境界側は通信先を列挙する"),
    ("T7", "内部関係者・物理アクセス", "0(対象外)", "0〜cond(権限分離・監査)", "どちらも単独では防げない。境界側は操作範囲を狭める"),
    ("T8", "境界装置自身の脆弱性", "n/a", "0(新しい攻撃面)", "B にだけ存在する。型式の多様化と段階更新で共通原因を抑える（下の blast_radius 表）"),
]


def coverage_table() -> list[dict]:
    return [{"id": i, "path": pth, "A_product_certification": a, "B_boundary_protection": b, "note": n} for i, pth, a, b, n in COVERAGE]


# ---------------------------------------------------------------- 同時に動く容量（Blast radius）

def blast_radius() -> list[dict]:
    rows = []
    for case in ("low", "base", "high"):
        fleet = p("fleet_mw", case)
        sv = p("share_vendor_cloud", case)
        rho = p("rho_boundary", case)
        sb = p("share_boundary_vendor", case)
        pb = p("p_boundary_vuln_year", case)
        a_cloud = fleet * sv  # クラウド侵害 1 件で届く容量（A）
        b_cloud = fleet * sv * (1 - rho)  # 境界が設計どおり働かないサイト分だけ届く（B）
        b_boundary_cm = fleet * sb  # 境界装置型式の脆弱性が一斉悪用されたときの上限（B 固有）
        rows.append({
            "case": case, "fleet_mw": fleet,
            "A_cloud_compromise_mw": round(a_cloud), "B_cloud_compromise_mw": round(b_cloud),
            "B_boundary_common_mode_mw": round(b_boundary_cm), "B_boundary_common_mode_expected_mw_year": round(b_boundary_cm * pb, 1),
            "note": "A: 制御面に属する容量がそのまま動く。B: (1-ρ) のサイトだけ届く。B 固有: 境界装置型式の共通原因（年率 p で期待値化）",
        })
    # 境界装置の型式数で共通原因を分散
    for k in (1, 2, 3, 5, 8):
        fleet = p("fleet_mw")
        rows.append({
            "case": f"boundary_vendors_{k}", "fleet_mw": fleet,
            "A_cloud_compromise_mw": round(fleet * p("share_vendor_cloud")), "B_cloud_compromise_mw": round(fleet * p("share_vendor_cloud") * (1 - p("rho_boundary"))),
            "B_boundary_common_mode_mw": round(fleet / k), "B_boundary_common_mode_expected_mw_year": round(fleet / k * p("p_boundary_vuln_year"), 1),
            "note": f"境界装置を {k} 型式に分散（均等）したときの共通原因の上限",
        })
    return rows


# ---------------------------------------------------------------- 図 7（SVG を直接書く）

def _fmt_yen(v: float) -> str:
    if v >= 1e6:
        return f"{v/1e6:.1f}百万円"
    if v >= 1e4:
        return f"{v/1e4:.0f}万円"
    return f"{v:,.0f}円"


def fig_cost_svg(rows: list[dict]) -> str:
    base = [r for r in rows if r["case"] == "base"]
    low = {r["site_class"]: r for r in rows if r["case"] == "low"}
    high = {r["site_class"]: r for r in rows if r["case"] == "high"}
    W, H = 760, 420
    ml, mr, mt, mb = 70, 20, 50, 70
    pw, ph = W - ml - mr, H - mt - mb
    # 対数軸（円/kW/年）
    vals = []
    for r in base:
        vals += [r["A_total_year_per_kw"], r["B_total_year_per_kw"]]
    for d in (low, high):
        for r in d.values():
            vals += [r["A_total_year_per_kw"], r["B_total_year_per_kw"]]
    vmin = max(1.0, min(v for v in vals if v > 0) / 2)
    vmax = max(vals) * 2
    lo, hi = math.log10(vmin), math.log10(vmax)

    def y(v: float) -> float:
        v = max(v, vmin)
        return mt + ph - (math.log10(v) - lo) / (hi - lo) * ph

    n = len(base)
    gw = pw / n
    bw = gw * 0.28
    out = [f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {W} {H}" role="img" aria-labelledby="fig-cost-title fig-cost-desc" class="i3wp-chart">',
           '<title id="fig-cost-title">製品認証方式と境界防護方式の 1 kW あたり年間費用（対数軸・仮定値）</title>',
           '<desc id="fig-cost-desc">横軸はサイト規模 4 区分、縦軸は 1 kW あたりの年間費用（円、対数）。青が製品認証方式 A（認証費用の価格転嫁＋交換遅延の逸失電力）、緑が境界防護方式 B（境界装置の年割り＋運用費）。細線は low〜high の幅。</desc>',
           f'<rect x="0" y="0" width="{W}" height="{H}" fill="#ffffff" rx="10"/>']
    # grid
    for e in range(math.ceil(lo), math.floor(hi) + 1):
        yy = y(10 ** e)
        out.append(f'<line x1="{ml}" y1="{yy:.1f}" x2="{W-mr}" y2="{yy:.1f}" stroke="#e5e5ea" stroke-width="1"/>')
        out.append(f'<text x="{ml-8}" y="{yy+4:.1f}" text-anchor="end" font-size="11" fill="#6e6e73">{10**e:,.0f}</text>')
    out.append(f'<text x="18" y="{mt+ph/2:.0f}" text-anchor="middle" font-size="11" fill="#6e6e73" transform="rotate(-90 18 {mt+ph/2:.0f})">円 / kW / 年（対数）</text>')
    for i, r in enumerate(base):
        cx = ml + gw * i + gw / 2
        sc = r["site_class"]
        for j, (key, color, label) in enumerate((("A_total_year_per_kw", "#0071e3", "A"), ("B_total_year_per_kw", "#2f7a4d", "B"))):
            x0 = cx - bw - 4 if j == 0 else cx + 4
            v = r[key]
            top = y(v)
            out.append(f'<rect x="{x0:.1f}" y="{top:.1f}" width="{bw:.1f}" height="{mt+ph-top:.1f}" fill="{color}" opacity="0.85" rx="3"/>')
            vl, vh = low[sc][key], high[sc][key]
            xm = x0 + bw / 2
            out.append(f'<line x1="{xm:.1f}" y1="{y(vl):.1f}" x2="{xm:.1f}" y2="{y(vh):.1f}" stroke="#1d1d1f" stroke-width="1.2"/>')
            out.append(f'<line x1="{xm-4:.1f}" y1="{y(vl):.1f}" x2="{xm+4:.1f}" y2="{y(vl):.1f}" stroke="#1d1d1f" stroke-width="1.2"/>')
            out.append(f'<line x1="{xm-4:.1f}" y1="{y(vh):.1f}" x2="{xm+4:.1f}" y2="{y(vh):.1f}" stroke="#1d1d1f" stroke-width="1.2"/>')
            out.append(f'<text x="{xm:.1f}" y="{top-6:.1f}" text-anchor="middle" font-size="10.5" fill="{color}">{v:,.0f}</text>')
        out.append(f'<text x="{cx:.1f}" y="{mt+ph+20}" text-anchor="middle" font-size="11.5" fill="#1d1d1f">{r["label"]}</text>')
    # legend
    lx = ml
    ly = mt - 28
    out.append(f'<rect x="{lx}" y="{ly}" width="14" height="14" fill="#0071e3" rx="3"/><text x="{lx+20}" y="{ly+12}" font-size="11.5" fill="#1d1d1f">A 製品認証方式（認証費用の転嫁＋同型式交換が重くなる分の逸失電力）</text>')
    out.append(f'<rect x="{lx}" y="{ly+18}" width="14" height="14" fill="#2f7a4d" rx="3"/><text x="{lx+20}" y="{ly+30}" font-size="11.5" fill="#1d1d1f">B 境界防護方式（境界装置の年割り＋運用費＋短い交換遅延）　細線: low〜high</text>')
    out.append(f'<text x="{W/2:.0f}" y="{H-14}" text-anchor="middle" font-size="10.5" fill="#6e6e73">入力はすべて data/assumptions.csv の仮定値。kind=to-confirm の項目（申請手数料・有効期間）は公表値で置き換える。</text>')
    out.append("</svg>")
    return "\n".join(out)


# ---------------------------------------------------------------- 書き出し

def write_csv(path: Path, rows: list[dict]) -> None:
    with path.open("w", encoding="utf-8", newline="") as f:
        w = csv.DictWriter(f, fieldnames=list(rows[0].keys()))
        w.writeheader()
        w.writerows(rows)


def main() -> None:
    costs = cost_table()
    be = break_even_table()
    cov = coverage_table()
    br = blast_radius()
    torn = tornado("hv") + tornado("lv_res")
    write_csv(RESULTS / "cost_per_site.csv", costs)
    write_csv(RESULTS / "break_even.csv", be)
    write_csv(RESULTS / "coverage_matrix.csv", cov)
    write_csv(RESULTS / "blast_radius.csv", br)
    write_csv(RESULTS / "tornado.csv", torn)
    base = {r["site_class"]: r for r in costs if r["case"] == "base"}
    summary = {
        "version": "1.0",
        "base": {k: {"A_total_year": v["A_total_year"], "B_total_year": v["B_total_year"], "A_per_kw": v["A_total_year_per_kw"], "B_per_kw": v["B_total_year_per_kw"]} for k, v in base.items()},
        "break_even_capex": {r["site_class"]: r["break_even_boundary_capex"] for r in be if not r["site_class"].startswith("hv_sweep")},
        "blast_radius_base": next(r for r in br if r["case"] == "base"),
        "to_confirm": [k for k, v in A.items() if v["kind"] == "to-confirm"],
        "assumptions": [k for k, v in A.items() if v["kind"] == "assumption"],
    }
    (RESULTS / "summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2), encoding="utf-8")
    for k, v in base.items():
        print(f"{v['label']:<22} A {v['A_total_year']:>10,} 円/年 ({v['A_total_year_per_kw']:>8,} 円/kW)   B {v['B_total_year']:>10,} 円/年 ({v['B_total_year_per_kw']:>8,} 円/kW)")
    print("break-even capex:", summary["break_even_capex"])
    print("blast radius (base):", summary["blast_radius_base"])


if __name__ == "__main__":
    main()