"""レビュー（round1〜4）で VALID / PARTIALLY VALID とした論点の感度計算。
dwpt_model.py の関数を使い、results/ に CSV を書く。
実行: python3 model/sensitivities.py
"""
from __future__ import annotations

import csv
import json
import pathlib
import sys

sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
import dwpt_model as M  # noqa: E402

RESULTS = M.RESULTS
low, base, high = M.load_assumptions("low"), M.load_assumptions("base"), M.load_assumptions("high")


def write(name: str, rows):
    M.write_csv(RESULTS / name, rows)
    print(f"saved {name} ({len(rows)} rows)")


def flow_tateyama(a, year, pen=1.0):
    return M.target_flow_per_day(a["aadt_total"], a["heavy_vehicle_fraction"], a[f"ev_fraction_heavy_{year}"], pen, a["lane_fraction"])


# -----------------------------------------------------------------------------
# (A) MCS 側の前提を推進側の主張どおり不利にしたら代替案はいくらになるか（R1-14/16/17/18/19/20, R3-06）
# -----------------------------------------------------------------------------
def mcs_variants():
    rows = []
    a = base
    grid_per_stall = a["mcs_site_grid_capex_yen"] / a["mcs_stalls_per_site"]
    variants = [
        ("BASE（稼働率 15 %・同時率 0.6・系統 1.6 億/6 口・効率 93 %）", dict()),
        ("同時率 1.0（実量制でピーク全額。R1-16）", dict(coincidence=1.0)),
        ("特別高圧連系 8 億円/6 口（R1-17）", dict(grid_capex_yen_per_stall=8e8 / 6)),
        ("充電効率 88 %（電池側損失・冷却込み。R1-18）", dict(eff=0.88)),
        ("C レート 1.2 → 平均出力比 0.55（R1-14）", dict(avg_power_fraction=0.55)),
        ("稼働率 5 %（導入初期。R1-15/R3-06）", dict(utilisation=0.05)),
        ("待ち行列で時間費用 2 倍（R1-20）", dict(time_cost_yen_per_h=a["driver_vehicle_hour_yen"] * 2)),
        ("上の不利条件をすべて同時に（稼働率は 15 %）", dict(coincidence=1.0, grid_capex_yen_per_stall=8e8 / 6, eff=0.88, avg_power_fraction=0.55, time_cost_yen_per_h=a["driver_vehicle_hour_yen"] * 2)),
        ("上の不利条件をすべて同時に（稼働率 5 %）", dict(coincidence=1.0, grid_capex_yen_per_stall=8e8 / 6, eff=0.88, avg_power_fraction=0.55, time_cost_yen_per_h=a["driver_vehicle_hour_yen"] * 2, utilisation=0.05)),
    ]
    for name, kw in variants:
        args = dict(charger_capex_yen=a["mcs_charger_capex_yen"], grid_capex_yen_per_stall=grid_per_stall, rate=a["discount_rate"], lifetime_years=a["mcs_lifetime_years"],
                    om_fraction=a["mcs_om_fraction"], utilisation=a["mcs_utilization"], power_kw=a["mcs_power_kw"], avg_power_fraction=0.7, electricity_yen_per_kwh=a["electricity_yen_per_kwh"],
                    demand_charge_yen_per_kw_month=a["demand_charge_yen_per_kw_month"], eff=0.93, time_cost_yen_per_h=a["driver_vehicle_hour_yen"], logistics_overlap=a["logistics_overlap"])
        args.update(kw)
        m = M.mcs_lcos(**args)
        # 車載 MCS インレット・冷却配線 30 万円/10 年を車両側費用として加える（R1-19）
        veh = M.annualised_cost(300_000, a["discount_rate"], 10.0, 0.0) / (a["truck_kwh_per_km"] * a["truck_km_per_year"] * 0.5)
        bat = M.alternative_lcos(a, 2040)["battery"].total_yen_per_kwh
        total = m.total_yen_per_kwh + veh + bat
        be = M.break_even_flow(a, total)
        rows.append({"variant": name, "mcs_yen_per_kwh": round(m.total_yen_per_kwh, 1), "mcs_vehicle_inlet_yen_per_kwh": round(veh, 1), "battery_2040_yen_per_kwh": round(bat, 1),
                     "alternative_total_yen_per_kwh": round(total, 1), "break_even_flow_base_capex": round(be), "tateyama_2040_pen100_flow": round(flow_tateyama(a, 2040), 1)})
    return rows


# -----------------------------------------------------------------------------
# (B) 電池側の前提を推進側の主張どおりにしたら（R1-11/12/13, R3-23, R4-07）
# -----------------------------------------------------------------------------
def battery_variants():
    rows = []
    a = base
    for name, year, dk, price, util in [
        ("BASE 2026 ΔE 200 kWh", 2026, 200, None, 1.0),
        ("BASE 2040 ΔE 200 kWh", 2040, 200, None, 1.0),
        ("電池 HIGH 価格 2026（28,800 円/kWh）", 2026, 200, high["battery_pack_yen_per_kwh_2026"], 1.0),
        ("電池 HIGH 価格 2040（16,000 円/kWh）", 2040, 200, high["battery_pack_yen_per_kwh_2040"], 1.0),
        ("顧客負担価格 40,000 円/kWh（R1-12）", 2026, 200, high["battery_customer_yen_per_kwh_2026"], 1.0),
        ("ΔE 500 kWh（推進側の比較軸。R1-13）2026", 2026, 500, None, 1.0),
        ("ΔE 500 kWh 2040", 2040, 500, None, 1.0),
        ("ΔE 実使用率 1/3（最悪日寸法決め。R1-11）2026", 2026, 200, None, 0.33),
        ("ΔE 実使用率 1/3・HIGH 価格・ΔE 500（全部不利）2026", 2026, 500, high["battery_pack_yen_per_kwh_2026"], 0.33),
        ("ΔE 実使用率 1/3・HIGH 価格・ΔE 500（全部不利）2040", 2040, 500, high["battery_pack_yen_per_kwh_2040"], 0.33),
    ]:
        aa = dict(a, delta_kwh_utilisation=util)
        alt = M.alternative_lcos(aa, year, delta_kwh=dk, battery_price=price)
        be = M.break_even_flow(aa, alt["total_yen_per_kwh"])
        rows.append({"variant": name, "year": year, "delta_kwh": dk, "battery_price_yen_per_kwh": price or a[f"battery_pack_yen_per_kwh_{year}"], "delta_utilisation": util,
                     "battery_yen_per_kwh": round(alt["battery"].total_yen_per_kwh, 1), "mcs_yen_per_kwh": round(alt["mcs"].total_yen_per_kwh, 1),
                     "alternative_total_yen_per_kwh": round(alt["total_yen_per_kwh"], 1), "break_even_flow_base_capex": round(be)})
    return rows


# -----------------------------------------------------------------------------
# (C) 乗用車も受電する場合（R1-32, R2-25）
# -----------------------------------------------------------------------------
def passenger_scenario():
    rows = []
    a = base
    for year in (2035, 2040):
        for fit_name, fit in (("LOW", low["passenger_dwpt_fit_rate"]), ("BASE", base["passenger_dwpt_fit_rate"]), ("HIGH", high["passenger_dwpt_fit_rate"])):
            bev_pass = a["passenger_bev_fraction_2040"] * (0.6 if year == 2035 else 1.0)
            pass_flow = a["aadt_total"] * (1 - a["heavy_vehicle_fraction"]) * bev_pass * fit * a["lane_fraction"]
            heavy_flow = flow_tateyama(a, year, 1.0)
            e_h = M.energy_per_pass_kwh(a["peak_power_kw"], 1000.0, a["speed_kmh"], a["coupling_duty"])
            e_p = M.energy_per_pass_kwh(a["passenger_receive_kw"], 1000.0, 100.0, a["coupling_duty"])
            annual = M.annual_energy_kwh(heavy_flow, e_h) + M.annual_energy_kwh(pass_flow, e_p)
            ann = M.annualised_cost(a["capex_dwpt_yen_per_lane_km"], a["discount_rate"], a["lifetime_years"], a["om_fraction"])
            infra = ann / annual
            # 乗用車の受電器は J2954 静止無線装備車の増分 10 万円/10 年、年間 25,000 km 回廊走行
            recv_p = M.annualised_cost(100_000, a["discount_rate"], 10.0, 0.0) / (e_p * a["dwpt_km_per_vehicle_year"])
            elec = a["electricity_yen_per_kwh"] / a["grid_to_battery_eff"]
            demand = M.dwpt_lcos_from(a, heavy_flow + pass_flow * e_p / e_h).demand_charge_yen_per_kwh
            # 乗用車側の代替：自宅・目的地の普通充電 25 円/kWh（基本料金込み）、電池削減便益は無視
            rows.append({"year": year, "passenger_fit_case": fit_name, "passenger_dwpt_vehicles_per_day": round(pass_flow), "heavy_dwpt_vehicles_per_day": round(heavy_flow, 1),
                         "passenger_share_of_kwh": round(M.annual_energy_kwh(pass_flow, e_p) / annual, 2), "infra_yen_per_kwh": round(infra, 1),
                         "total_yen_per_kwh_passenger_side": round(infra + recv_p + elec + demand, 1), "passenger_alternative_home_charging_yen_per_kwh": 25.0,
                         "heavy_only_infra_yen_per_kwh": round(M.dwpt_lcos_from(a, heavy_flow).infra_yen_per_kwh, 1)})
    return rows


# -----------------------------------------------------------------------------
# (D) 時間軸：2027 年着工で 20 年、年ごとの交通量で割る（R3-04/R3-15）
# -----------------------------------------------------------------------------
def time_axis():
    rows = []
    a = base
    years = {2026: a["ev_fraction_heavy_2026"], 2030: a["ev_fraction_heavy_2030"], 2035: a["ev_fraction_heavy_2035"], 2040: a["ev_fraction_heavy_2040"]}
    def ev(y):
        ks = sorted(years)
        if y <= ks[0]: return years[ks[0]]
        if y >= ks[-1]: return years[ks[-1]] * (1 + 0.05) ** (y - ks[-1])  # 2040 以降 年 5 % 増
        for k0, k1 in zip(ks, ks[1:]):
            if k0 <= y <= k1:
                return years[k0] * (years[k1] / years[k0]) ** ((y - k0) / (k1 - k0))
    e1 = M.energy_per_pass_kwh(a["peak_power_kw"], 1000.0, a["speed_kmh"], a["coupling_duty"])
    for build_year in (2027, 2032, 2037):
        for pen in (0.5, 1.0):
            for life, r_, om in ((20, 0.05, 0.02), (30, 0.03, 0.01)):
                pv_cost = a["capex_dwpt_yen_per_lane_km"]
                pv_kwh = 0.0
                pv_alt_value = 0.0
                for t in range(life):
                    y = build_year + t
                    flow = M.target_flow_per_day(a["aadt_total"], a["heavy_vehicle_fraction"], min(ev(y), 1.0), pen, a["lane_fraction"])
                    kwh = M.annual_energy_kwh(flow, e1)
                    disc = (1 + r_) ** (-(t + 0.5))
                    pv_cost += a["capex_dwpt_yen_per_lane_km"] * om * disc
                    pv_kwh += kwh * disc
                    alt_year = 2040 if y >= 2040 else (2035 if y >= 2035 else (2030 if y >= 2030 else 2026))
                    alt = M.alternative_lcos(a, alt_year)["total_yen_per_kwh"]
                    pv_alt_value += kwh * (alt - a["electricity_yen_per_kwh"] / a["grid_to_battery_eff"]) * disc  # 代替案の非電力費を便益とみなす
                lev = pv_cost / pv_kwh if pv_kwh else float("inf")
                rows.append({"build_year": build_year, "life_years": life, "discount_rate": r_, "dwpt_penetration": pen, "levelised_infra_yen_per_kwh": round(lev, 1),
                             "benefit_cost_ratio_vs_alternative": round(pv_alt_value / pv_cost, 3), "pv_kwh_per_lane_km_mwh": round(pv_kwh / 1000, 1)})
    return rows


# -----------------------------------------------------------------------------
# (E) 外部性：結論が変わる閾値（R3-12, R2-15/16/18）
# -----------------------------------------------------------------------------
def externality_thresholds():
    a = base
    rows = []
    for year, pen in ((2040, 1.0), (2040, 0.5), (2035, 1.0)):
        flow = flow_tateyama(a, year, pen)
        l = M.dwpt_lcos_from(a, flow)
        alt = M.alternative_lcos(a, year)["total_yen_per_kwh"]
        gap = l.total_yen_per_kwh - alt  # 円/kWh：外部便益がこれだけあれば並ぶ
        annual_kwh_300m = l.annual_kwh_per_lane_km * a["section_length_m"] / 1000
        gap_yen_per_year_300m = gap * annual_kwh_300m
        # 電池 200 kWh 削減の軸重 1.2 t 減 → 舗装損傷（4 乗則）：館山道 大型車 1 日 flow 台、25 t 車の軸重 10 t 基準
        # 路面維持費 概算 500 万円/車線 km/年（高速道路の舗装維持管理費の桁）、大型車寄与 90 %
        pav_saving_per_vehicle_km = 5e6 * 0.9 / max(a["aadt_total"] * a["heavy_vehicle_fraction"] * a["lane_fraction"] * 365, 1) * (1 - ((25 - 1.2) / 25) ** 4)
        pav_saving_yen_per_year_300m = pav_saving_per_vehicle_km * flow * 365 * a["section_length_m"] / 1000
        # CO2：DWPT と MCS は同じ系統電力。差は電池製造 CO2（60 kgCO2/kWh × 200 kWh ÷ 10 年）× 炭素価格
        co2_saving_t_per_vehicle_year = 60 * 200 / 1000 / 10
        rows.append({"year": year, "dwpt_penetration": pen, "target_flow_per_day": round(flow, 1), "dwpt_total_yen_per_kwh": round(l.total_yen_per_kwh, 1), "alternative_yen_per_kwh": round(alt, 1),
                     "gap_yen_per_kwh": round(gap, 1), "annual_kwh_300m": round(annual_kwh_300m), "external_benefit_needed_yen_per_year_300m": round(gap_yen_per_year_300m), "external_benefit_needed_yen_per_lane_km_year": round(gap * l.annual_kwh_per_lane_km),
                     "external_benefit_needed_yen_per_dwpt_vehicle_year": round(gap_yen_per_year_300m / max(flow, 1e-9)),
                     "pavement_saving_yen_per_year_300m_4th_power": round(pav_saving_yen_per_year_300m),
                     "carbon_price_needed_yen_per_tCO2_if_only_battery_co2": round(gap_yen_per_year_300m / max(flow * co2_saving_t_per_vehicle_year, 1e-9)) if gap > 0 else 0,
                     "stranding_events_avoided_needed_per_year_at_50man_yen": round(gap_yen_per_year_300m / 500_000, 1) if gap > 0 else 0})
    return rows


# -----------------------------------------------------------------------------
# (F) 両案で未計上の費用と方向（R3-22）— 表としてそのまま記事に載せる
# -----------------------------------------------------------------------------
UNCOUNTED = [
    ("DWPT", "受変電・系統引込・制御・通信の区間固定費", "DWPT 不利（300 m では特に）", "length_scaling の with_fixed_cost 列で 0.6 億円/区間として感度計上"),
    ("DWPT", "舗装打換え周期とコイル再敷設の同期費", "DWPT 不利", "未計上。寿命 10 年 HIGH で代理"),
    ("DWPT", "舗装打換え時に敷設すれば土木費の一部が既存予算と重なる", "DWPT 有利", "未計上（R1-07）。CAPEX LOW 側で代理"),
    ("DWPT", "受電器の車両側統合（車検・保安基準対応）", "DWPT 不利", "receiver_cost HIGH 150 万円で代理"),
    ("DWPT", "乗用車の受電による kWh 増", "DWPT 有利", "passenger_scenario.csv で感度計上"),
    ("DWPT", "電池削減による車両価格低下・積載増（円/kWh では不完全）", "DWPT 有利", "fig15 円/車両 km・円/t km で併記。TCO は未計上"),
    ("DWPT", "電欠・立ち往生の削減、電池資源の対外依存低減", "DWPT 有利", "externality_thresholds.csv で必要便益額を提示"),
    ("MCS", "ピーク時の待ち行列", "MCS 不利", "mcs_variants 時間費用 2 倍で感度計上"),
    ("MCS", "特別高圧連系（20 kV 以上）の費用", "MCS 不利", "mcs_variants 8 億円/6 口で感度計上"),
    ("MCS", "車載インレット・冷却配線", "MCS 不利", "mcs_variants 30 万円/台で計上"),
    ("MCS", "高 C レート充電による電池劣化加速", "MCS 不利", "C レート 1.2 で平均出力比 0.55、サイクル寿命 LOW 2,000 で代理"),
    ("MCS", "SA/PA の大型車駐車マス不足（口数制約）", "MCS 不利", "未計上（R1-21）"),
    ("MCS", "蓄電池併設による契約電力削減", "MCS 有利", "mcs_hub で感度計上"),
    ("電池", "残存価値（10 年後の二次利用）", "代替案有利", "未計上（R3-23）。償却は車両寿命 10 年"),
    ("電池", "顧客負担価格がパック卸値より高い", "代替案不利", "battery_variants 40,000 円/kWh で感度計上"),
    ("電池", "最悪日寸法決めによる ΔE 実使用率低下", "代替案不利", "battery_variants 実使用率 1/3 で感度計上"),
    ("共通", "電力の時間帯別単価（昼間再エネ余剰）", "両案とも享受可能。MCS は蓄電池で選択可", "未計上（R1-30）"),
    ("共通", "価格基準年と USD→円の物価補正", "DWPT CAPEX は 2015〜2021 年ドル建て。補正すれば DWPT に 2〜3 割不利", "未計上（R2-08/R3-11）。LOW 側で代理"),
]


# -----------------------------------------------------------------------------
# (G) 館山道交通量での一変数感度（R5-02, R5-06）
# -----------------------------------------------------------------------------
def dream_tech_at_tateyama_tornado():
    """Dream の技術・資金前提＋館山道 BASE 交通量＋搭載率 100% を起点に、DWPT 側の変数を 1 つずつ BASE へ戻す。"""
    d = M.dream_case_assumptions(low, base, high)
    dt = dict(d, aadt_total=base["aadt_total"], heavy_vehicle_fraction=base["heavy_vehicle_fraction"], ev_fraction_heavy_2040=base["ev_fraction_heavy_2040"], lane_fraction=base["lane_fraction"])
    f = flow_tateyama(dt, 2040, 1.0)
    alt_base = M.alternative_lcos(base, 2040)["total_yen_per_kwh"]
    rows = [{"variable_reset_to_base": "起点（Dream 技術・資金前提、館山道 BASE 交通量、搭載率 100 %）", "dwpt_total_yen_per_kwh": round(M.dwpt_lcos_from(dt, f).total_yen_per_kwh, 1), "alt_base_2040_yen_per_kwh": round(alt_base, 1), "target_flow_per_day": round(f, 1)}]
    for name, upd in [
        ("受電出力 150 kW", dict(peak_power_kw=150.0)), ("CAPEX BASE 3.0 億円", dict(capex_dwpt_yen_per_lane_km=base["capex_dwpt_yen_per_lane_km"])),
        ("寿命 20 年・5 %・O&M 2 %", dict(lifetime_years=20.0, discount_rate=0.05, om_fraction=0.02)), ("速度 90 km/h", dict(speed_kmh=90.0)),
        ("実効出力比 0.75", dict(coupling_duty=0.75)), ("効率 0.85", dict(grid_to_battery_eff=0.85)), ("区間固定費 0.6 億円", dict(fixed_capex_per_segment_yen=6e7)),
        ("受電器 80 万円・10 年", dict(receiver_cost_yen=8e5, receiver_lifetime_years=10.0)), ("搭載率 50 %", dict(_pen=0.5)),
    ]:
        pen = upd.pop("_pen", 1.0)
        a = dict(dt, **upd)
        rows.append({"variable_reset_to_base": name, "dwpt_total_yen_per_kwh": round(M.dwpt_lcos_from(a, f * pen).total_yen_per_kwh, 1), "alt_base_2040_yen_per_kwh": round(alt_base, 1), "target_flow_per_day": round(f * pen, 1)})
    return rows


def tateyama_base_tornado():
    """館山道 2040 年・搭載率 100 %・BASE（477 円）を起点に、1 変数だけ最良値へ動かす。"""
    alt_base = M.alternative_lcos(base, 2040)["total_yen_per_kwh"]
    f0 = flow_tateyama(base, 2040, 1.0)
    rows = [{"variable_set_to_best": "起点（BASE）", "dwpt_total_yen_per_kwh": round(M.dwpt_lcos_from(base, f0).total_yen_per_kwh, 1), "alt_base_2040_yen_per_kwh": round(alt_base, 1), "target_flow_per_day": round(f0, 1)}]
    for name, upd in [
        ("BEV 比率 45 %", dict(ev_fraction_heavy_2040=0.45)), ("受電出力 300 kW", dict(peak_power_kw=300.0)), ("寿命 50 年・3 %", dict(lifetime_years=50.0, discount_rate=0.03)),
        ("寿命 30 年・3 %・O&M 1 %", dict(lifetime_years=30.0, discount_rate=0.03, om_fraction=0.01)), ("CAPEX LOW 1.9 億円", dict(capex_dwpt_yen_per_lane_km=low["capex_dwpt_yen_per_lane_km"])),
        ("交通量 HIGH（16,000 台・大型車 15 %）", dict(aadt_total=16000.0, heavy_vehicle_fraction=0.15)), ("実効出力比 0.9", dict(coupling_duty=0.9)), ("速度 80 km/h", dict(speed_kmh=80.0)),
        ("敷設車線 50 %", dict(lane_fraction=0.5)), ("受電器 30 万円・15 年", dict(receiver_cost_yen=3e5, receiver_lifetime_years=15.0)), ("効率 0.92", dict(grid_to_battery_eff=0.92)),
        ("寿命 200 年（上限確認）", dict(lifetime_years=200.0)),
    ]:
        a = dict(base, **upd)
        f = flow_tateyama(a, 2040, 1.0)
        rows.append({"variable_set_to_best": name, "dwpt_total_yen_per_kwh": round(M.dwpt_lcos_from(a, f).total_yen_per_kwh, 1), "alt_base_2040_yen_per_kwh": round(alt_base, 1), "target_flow_per_day": round(f, 1)})
    rows[1:] = sorted(rows[1:], key=lambda r: r["dwpt_total_yen_per_kwh"])
    return rows


def main():
    write("mcs_variants.csv", mcs_variants())
    write("battery_variants.csv", battery_variants())
    write("passenger_scenario.csv", passenger_scenario())
    write("time_axis_npv.csv", time_axis())
    write("externality_thresholds.csv", externality_thresholds())
    write("dream_at_tateyama_tornado.csv", dream_tech_at_tateyama_tornado())
    write("tateyama_base_tornado.csv", tateyama_base_tornado())
    write("uncounted_costs.csv", [{"side": s, "item": i, "direction": d, "treatment": t} for s, i, d, t in UNCOUNTED])


if __name__ == "__main__":
    main()
