# -*- coding: utf-8 -*-
"""
ICM 2018 D — Out of Gas and Driving on E (for electric, not empty)
真实赛题：COMAP 国际交叉学科建模竞赛（ICM）D 题——电动汽车充电网络架构规划。
  T1 美国当前/增长的特斯拉充电网络；若全美私人乘用车全部电动化，需多少充电站、城乡郊如何分配？
  T2 选一国（韩/爱/乌）做深析：2A 瞬时全电动化的最优数量/布局/分布；2B 从零到全电动的演进提案；
     2C 10%/30%/50%/100% 电动化时间表。
  T3 五国（澳/中/印尼/沙特/新）地理/人口密度/财富各异，方案是否适用？建立分类系统判定增长模型。
  T4 共享出行/自动驾驶/换电/飞行汽车/Hyperloop 等技术对分析的影响。
  T5 给多国领导人的一页非技术备忘（关键要素清单）。

真源模型（纯标准库，确定性；SEED 仅作版本锚定，模型无随机噪声）：
  1. 车辆基数：vehicles = pop_M * veh_per_1000 / 1000。
  2. 年充电能量需求：E_annual(TWh) = vehicles*1e6 * VKT * E_PER_KM / 1e9。
  3. 公共快充（超级充电站）泊位数（需求侧）：
        super_stalls_demand = E_daily * FAST_FRAC / (P_SUPER * HOURS * UTIL_SUPER)
     覆盖侧地板（农村/稀疏区保障可达性）：
        super_stalls_cover = area_km2 * REDUND / (pi * R_COVER^2)
     实际 super_stalls = max(需求, 覆盖)。公共充电站数 stations = super_stalls / STALLS_PER_STATION。
  4. 目的地充电点（家/单位慢充）：dest_points = vehicles*1e6 * (HOME_FRAC + 0.10)（家庭+单位）。
  5. 城乡郊分配：按 zone_frac 切分 super_stalls；农村再加覆盖地板，形成覆盖驱动微超建。
  6. 深析国（韩国）：同上 + 采用率逻辑斯蒂曲线 S(t)=1/(1+exp(-k(t-t0))) 求 10/30/50/99% 时点；
     投资估算 capex = stations * STATION_COST_K。
  7. 分类（T3）：density = pop/area；wealth = gdp_cap；按 (密度,财富) 归类四种增长原型，
     对澳/中/印尼/沙特/新分别推荐适用模型。

输出：assets/problems/data/icm2018d.csv（9 行国家面板 + 表头）
作者说明：国家面板为方法论演示用的确定性合成数据（种子 2018），数值贴近公开量级但非官方实测；
         真实参赛应以 IEA / BloombergNEF / 各国统计年鉴校准；本模型聚焦需求→覆盖→分区→
         演进→分类的可复现链路。
"""
import os
import math
import csv

SEED = 2018
HERE = os.path.dirname(os.path.abspath(__file__))

# ---------------------------------------------------------------------------
# 0. 模型常数（均作为假设显式声明，便于敏感性分析）
# ---------------------------------------------------------------------------
VKT = 15000.0          # 单车年行驶里程 km/veh/yr
E_PER_KM = 0.20        # BEV 能耗 kWh/km
P_DEST = 11.0          # 目的地充电桩功率 kW（家/单位慢充）
P_SUPER = 150.0        # 超级充电桩功率 kW
UTIL_DEST = 0.32       # 目的地桩日利用率
UTIL_SUPER = 0.20      # 超级桩日利用率
HOURS = 24.0
STALLS_PER_STATION = 8 # 每座公共快充站泊位数
HOME_FRAC = 0.90       # 家庭充电桩覆盖的电动乘用车比例
FAST_FRAC = 0.25       # 需公共快充的出行能量占比（长途/途中）
R_COVER = 40.0         # 农村覆盖半径 km（最坏情形可达性）
REDUND = 1.5           # 覆盖冗余系数
ZONE_FRAC = {"urban": 0.45, "suburban": 0.35, "rural": 0.20}
STATION_COST_K = 250.0 # 单座公共快充站投资 千美元（8 泊位硬件+电网接入）

DEEP_CODE = "KR"       # 深析国：韩国

# 采用率曲线（韩国演进时间表）
K_ADOPT = 0.40         # 逻辑斯蒂增长率 /yr
T0_ADOPT = 9.0         # 达到 50% 的年限（自 2023 起）

# 五国分类面板（T3）
TASK3 = ["AU", "CN", "ID", "SA", "SG"]

# ---------------------------------------------------------------------------
# 1. 国家面板（确定性合成，种子 2018；量级贴近公开数据但非官方实测）
#    (code, name_cn, name_en, region, pop_M, area_km2, gdp_cap_k, veh_per_1000, ev_share, urban_frac)
# ---------------------------------------------------------------------------
COUNTRIES = [
    ("US", "美国",       "United States", "美洲",   331.0, 9_834_000, 63.0, 850, 0.020, 0.83),
    ("KR", "韩国",       "South Korea",    "亚洲",    51.7,   100_000, 32.0, 470, 0.030, 0.82),
    ("IE", "爱尔兰",     "Ireland",        "欧洲",     5.0,    70_000, 84.0, 520, 0.025, 0.64),
    ("UY", "乌拉圭",     "Uruguay",        "美洲",     3.5,   176_000, 17.0, 600, 0.015, 0.77),
    ("AU", "澳大利亚",   "Australia",      "大洋洲",  25.7, 7_692_000, 52.0, 740, 0.020, 0.86),
    ("CN", "中国",       "China",          "亚洲",  1412.0, 9_597_000, 12.0, 220, 0.040, 0.65),
    ("ID", "印度尼西亚", "Indonesia",      "亚洲",   273.0, 1_911_000,  4.8,  90, 0.010, 0.57),
    ("SA", "沙特",       "Saudi Arabia",   "中东",    35.0, 2_250_000, 23.0, 330, 0.010, 0.84),
    ("SG", "新加坡",     "Singapore",      "亚洲",     5.9,       730, 73.0, 110, 0.050, 1.00),
]


def _clip(x, lo, hi):
    return max(lo, min(hi, x))


def _compute(c):
    code, name_cn, name_en, region, pop_M, area_km2, gdp_cap_k, veh_per_1000, ev_share, urban_frac = c
    vehicles_M = pop_M * veh_per_1000 / 1000.0
    vehicles = vehicles_M * 1e6
    E_annual_TWh = vehicles * VKT * E_PER_KM / 1e9
    E_daily_kWh = E_annual_TWh * 1e9 / 365.0
    E_fast_daily = E_daily_kWh * FAST_FRAC
    super_demand = E_fast_daily / (P_SUPER * HOURS * UTIL_SUPER)
    cell_area = math.pi * R_COVER ** 2
    super_cover = area_km2 * REDUND / cell_area
    super_stalls = max(super_demand, super_cover)
    stations = super_stalls / STALLS_PER_STATION
    dest_points = vehicles * (HOME_FRAC + 0.10)
    # 城乡郊分配
    base_u = super_stalls * ZONE_FRAC["urban"]
    base_s = super_stalls * ZONE_FRAC["suburban"]
    base_r = super_stalls * ZONE_FRAC["rural"]
    rural_floor = super_cover * 0.30
    rural_extra = max(0.0, rural_floor - base_r)
    rural_stalls = base_r + rural_extra
    total_stalls = base_u + base_s + rural_stalls
    capex_k = stations * STATION_COST_K  # 千美元
    density = pop_M * 1e6 / area_km2    # 人/km2
    return dict(
        code=code, name_cn=name_cn, name_en=name_en, region=region,
        pop_M=pop_M, area_km2=area_km2, gdp_cap_k=gdp_cap_k,
        veh_per_1000=veh_per_1000, ev_share=ev_share, urban_frac=urban_frac,
        vehicles_M=vehicles_M, vehicles=vehicles,
        E_annual_TWh=E_annual_TWh, E_daily_kWh=E_daily_kWh,
        super_demand=super_demand, super_cover=super_cover,
        super_stalls=super_stalls, stations=stations,
        dest_points=dest_points,
        zone={"urban": base_u, "suburban": base_s, "rural": rural_stalls},
        total_stalls=total_stalls, capex_k=capex_k, density=density,
    )


def _timeline():
    """韩国演进时间表：逻辑斯蒂 S(t)=1/(1+exp(-k(t-t0))) 求达到 p 的年限。"""
    out = {}
    for p in (0.10, 0.30, 0.50, 0.99):
        t = T0_ADOPT - math.log(1.0 / p - 1.0) / K_ADOPT
        out["p%02d" % int(p * 100)] = round(2023 + t, 1)
    return out


def _classify(row):
    """(密度,财富) 二维原型归类。"""
    d, w = row["density"], row["gdp_cap_k"]
    if d >= 200 and w >= 40:
        return "高密度-高财富：城市优先（city-first）"
    if d >= 200 and w < 40:
        return "高密度-低财富：规模化公共充电（mass-rapid）"
    if d < 200 and w >= 40:
        return "低密度-高财富：走廊+平衡（corridor-balanced）"
    return "低密度-低财富：渐进式（gradual）"


def gen_icm2018d(seed=SEED, out_csv="icm2018d.csv"):
    rows = [_compute(c) for c in COUNTRIES]
    by_code = {r["code"]: r for r in rows}

    # 深析国（韩国）
    kr = by_code[DEEP_CODE]
    timeline = _timeline()
    deep = dict(kr)
    deep["timeline"] = timeline
    deep["adopt_k"] = K_ADOPT
    deep["adopt_t0"] = T0_ADOPT

    # 分类（T3 五国）
    classification = []
    for code in TASK3:
        r = by_code[code]
        classification.append(dict(
            code=code, name_cn=r["name_cn"], density=round(r["density"], 1),
            wealth=r["gdp_cap_k"], archetype=_classify(r),
            stations=round(r["stations"], 1), super_stalls=round(r["super_stalls"], 1),
        ))

    # 美国（T1 主对象）速览
    us = by_code["US"]

    # ---- 写 CSV（9 行国家面板 + 表头）----
    OUT_DIR = os.path.normpath(os.path.join(HERE, "..", "assets", "problems", "data"))
    os.makedirs(OUT_DIR, exist_ok=True)
    csv_path = os.path.join(OUT_DIR, out_csv)
    cols = ["code", "name_cn", "region", "pop_M", "area_km2", "gdp_cap_k",
            "veh_per_1000", "ev_share", "urban_frac", "vehicles_M",
            "E_annual_TWh", "super_stalls", "stations", "dest_points",
            "zone_urban", "zone_suburban", "zone_rural", "capex_k", "density"]
    with open(csv_path, "w", newline="", encoding="utf-8-sig") as f:
        w = csv.writer(f)
        w.writerow(cols)
        for r in rows:
            w.writerow([
                r["code"], r["name_cn"], r["region"], round(r["pop_M"], 1),
                r["area_km2"], r["gdp_cap_k"], r["veh_per_1000"], r["ev_share"],
                r["urban_frac"], round(r["vehicles_M"], 3),
                round(r["E_annual_TWh"], 2), round(r["super_stalls"], 1),
                round(r["stations"], 1), round(r["dest_points"], 1),
                round(r["zone"]["urban"], 1), round(r["zone"]["suburban"], 1),
                round(r["zone"]["rural"], 1), round(r["capex_k"], 1),
                round(r["density"], 2),
            ])

    D = dict(
        SEED=SEED, VKT=VKT, E_PER_KM=E_PER_KM, P_DEST=P_DEST, P_SUPER=P_SUPER,
        UTIL_DEST=UTIL_DEST, UTIL_SUPER=UTIL_SUPER, HOURS=HOURS,
        STALLS_PER_STATION=STALLS_PER_STATION, HOME_FRAC=HOME_FRAC,
        FAST_FRAC=FAST_FRAC, R_COVER=R_COVER, REDUND=REDUND,
        ZONE_FRAC=ZONE_FRAC, STATION_COST_K=STATION_COST_K,
        K_ADOPT=K_ADOPT, T0_ADOPT=T0_ADOPT, DEEP_CODE=DEEP_CODE,
        countries=rows, by_code=by_code, deep=deep,
        classification=classification, us=us,
        csv_path=csv_path,
        highlights=dict(
            us_stations=round(us["stations"], 1),
            us_super_stalls=round(us["super_stalls"], 1),
            us_dest_points=round(us["dest_points"], 1),
            us_E_annual_TWh=round(us["E_annual_TWh"], 1),
            kr_stations=round(kr["stations"], 1),
            kr_super_stalls=round(kr["super_stalls"], 1),
            kr_dest_points=round(kr["dest_points"], 1),
            kr_capex_k=round(kr["capex_k"], 1),
            kr_t10=timeline["p10"], kr_t30=timeline["p30"],
            kr_t50=timeline["p50"], kr_t99=timeline["p99"],
        ),
    )
    return D


if __name__ == "__main__":
    D = gen_icm2018d()
    h = D["highlights"]
    print("面板国家数 N =", len(D["countries"]))
    print("美国(EV 全电动化): 公共快充站=%.1f 座, 超充泊位=%.1f, 目的地充电点=%.1f, 年充电能量=%.1f TWh"
          % (h["us_stations"], h["us_super_stalls"], h["us_dest_points"], h["us_E_annual_TWh"]))
    print("韩国(深析): 公共快充站=%.1f 座, 超充泊位=%.1f, 目的地点=%.1f, 投资=%.1f 千美元"
          % (h["kr_stations"], h["kr_super_stalls"], h["kr_dest_points"], h["kr_capex_k"]))
    print("韩国演进: 10%%@%.1f / 30%%@%.1f / 50%%@%.1f / 99%%@%.1f (年)"
          % (h["kr_t10"], h["kr_t30"], h["kr_t50"], h["kr_t99"]))
    print("T3 分类:")
    for c in D["classification"]:
        print("  %s %s: 密度=%.1f 人/km2 财富=%d k$ -> %s"
              % (c["code"], c["name_cn"], c["density"], c["wealth"], c["archetype"]))
    print("CSV ->", D["csv_path"])