# -*- coding: utf-8 -*-
"""
ICM 2023 D — Prioritizing the UN Sustainable Development Goals (SDGs)
真实赛题：构建 17 个可持续发展目标之间的关系网络，据此设定推进联合国工作最有效的优先级；
评估每个优先级的有效性（未来 10 年合理可达成的目标）；考察"某 SDG 已实现"后的网络结构与优先级变化；
讨论技术进步 / 全球大流行 / 气候变化 / 地区战争与难民流动等危机对网络与优先级的影响；
说明该网络方法如何帮助其他公司 / 组织设定目标优先级。

真源模型（纯标准库，确定性，SEED 固定）：
  1. SDG 交互矩阵 M[i][j] ∈ {-3,-2,-1,0,1,2,3}，依据 ICSU (2017) SDG 交互七分型
     （Indivisible +3 / Enabling +2 / Reinforcing +1 / Consistent 0 / Constraining -1 /
      Counteracting -2 / Cancelling -3）及公开文献中的已知强关联手工标定，对称矩阵。
  2. 网络指标：度中心性、加权正/负协同度、净协同、特征向量中心性（幂迭代）、
     介数中心性（Brandes 算法，无权图）。
  3. 优先级综合得分：Priority = 0.50·eig_norm + 0.30·posdeg_norm + 0.20·netsyn_norm。
  4. What-if（某 SDG 已实现）：删点后重算 16 节点网络与排名，识别"损失最大的枢纽"。
  5. 四类危机冲击：对交互矩阵施加定向扰动后重算排名，量化优先级漂移。
  6. 10 年协同扩散有效性：以正协同矩阵为传播图，投资优先 SDG 后 10 年 cascade 受益总量。
  7. 企业 ESG 路线图应用：在 ESG 相关 8 个 SDG 子图上重排优先级。

输出：assets/problems/data/icm2023d.csv（17 行优先级排名表 + 表头）
作者说明：交互矩阵为"基于 ICSU 七分型与公开文献的综合性标定"，用于方法演示与可复现分析，
         真实建模应以各国实证交互数据校准；本模型聚焦方法论与网络结构。
"""
import os
import math

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

# ---------------------------------------------------------------------------
# 1. SDG 定义
# ---------------------------------------------------------------------------
SDGS = [
    (1,  "无贫困",                 "No Poverty"),
    (2,  "零饥饿",                 "Zero Hunger"),
    (3,  "良好健康与福祉",         "Good Health and Well-being"),
    (4,  "优质教育",               "Quality Education"),
    (5,  "性别平等",               "Gender Equality"),
    (6,  "清洁饮水与卫生",         "Clean Water and Sanitation"),
    (7,  "经济适用的清洁能源",     "Affordable and Clean Energy"),
    (8,  "体面工作与经济增长",     "Decent Work and Economic Growth"),
    (9,  "产业、创新和基础设施",   "Industry, Innovation and Infrastructure"),
    (10, "减少不平等",             "Reduced Inequalities"),
    (11, "可持续城市和社区",       "Sustainable Cities and Communities"),
    (12, "负责任消费和生产",       "Responsible Consumption and Production"),
    (13, "气候行动",               "Climate Action"),
    (14, "水下生物",               "Life Below Water"),
    (15, "陆地生物",               "Life on Land"),
    (16, "和平、正义与强大机构",   "Peace, Justice and Strong Institutions"),
    (17, "促进目标实现的伙伴关系", "Partnerships for the Goals"),
]
N = len(SDGS)
NAMES = [s[1] for s in SDGS]
ABB = [s[2] for s in SDGS]

# 交互矩阵（i<j 的 (i,j,v) 三元组；对称；其余为 0 = Consistent）
# 七分型：+3 不可分割 / +2 赋能 / +1 强化 / 0 一致 / -1 约束 / -2 抵消 / -3 互斥
_LINKS = [
    # SDG1 无贫困
    (1,2,3),(1,3,1),(1,4,1),(1,5,1),(1,8,2),(1,10,2),(1,12,1),(1,16,1),(1,17,1),
    # SDG2 零饥饿
    (2,3,1),(2,5,1),(2,6,2),(2,7,1),(2,12,1),(2,13,-1),(2,14,1),(2,15,-1),
    # SDG3 健康
    (3,4,1),(3,5,1),(3,6,2),(3,9,1),(3,11,1),(3,12,1),(3,13,-1),
    # SDG4 教育
    (4,5,2),(4,7,1),(4,8,2),(4,9,1),(4,10,1),(4,16,1),
    # SDG5 性别平等
    (5,6,1),(5,8,1),(5,10,2),(5,16,1),
    # SDG6 水
    (6,7,1),(6,11,1),(6,12,1),(6,13,-1),(6,14,1),(6,15,1),
    # SDG7 能源
    (7,8,1),(7,9,2),(7,11,1),(7,12,1),(7,13,1),(7,14,1),
    # SDG8 工作/增长
    (8,9,2),(8,10,1),(8,11,1),(8,12,1),(8,13,-1),(8,14,1),(8,17,1),
    # SDG9 产业/基建
    (9,11,2),(9,12,1),(9,13,-1),(9,15,-1),(9,17,1),
    # SDG10 不平等
    (10,11,1),(10,16,1),
    # SDG11 城市
    (11,12,1),(11,13,-1),(11,15,-1),(11,17,1),
    # SDG12 消费
    (12,13,1),(12,14,1),(12,15,1),(12,17,1),
    # SDG13 气候
    (13,14,1),(13,15,1),(13,16,1),(13,17,1),
    # SDG14 海洋
    (14,15,1),(14,17,1),
    # SDG15 陆地
    (15,17,1),
    # SDG16 和平/机构
    (16,17,2),
]

def _build_matrix():
    M = [[0]*N for _ in range(N)]
    for i,j,v in _LINKS:
        a,b = i-1, j-1
        M[a][b] = v
        M[b][a] = v
    return M

# 优先级权重
W_EIG = 0.50
W_POS = 0.30
W_NET = 0.20

# ---------------------------------------------------------------------------
# 2. 网络度量
# ---------------------------------------------------------------------------
def _eigenvector(P, iters=2000, tol=1e-12):
    """对非负加权矩阵 P 做幂迭代求特征向量中心性（自动适配矩阵维度）。"""
    n = len(P)
    x = [1.0/n]*n
    for _ in range(iters):
        y = [0.0]*n
        for i in range(n):
            for j in range(n):
                y[i] += P[i][j]*x[j]
        s = math.sqrt(sum(v*v for v in y)) or 1.0
        y = [v/s for v in y]
        d = math.sqrt(sum((y[i]-x[i])**2 for i in range(n)))
        x = y
        if d < tol:
            break
    return x

def _betweenness(adj_sets):
    """Brandes 算法（无权无向图），返回每节点介数中心性。"""
    C = [0.0]*N
    for s in range(N):
        S = []
        P = [[] for _ in range(N)]
        sigma = [0.0]*N; sigma[s] = 1.0
        d = [-1]*N; d[s] = 0
        Q = [s]
        while Q:
            v = Q.pop(0); S.append(v)
            for w in adj_sets[v]:
                if d[w] < 0:
                    d[w] = d[v] + 1; Q.append(w)
                if d[w] == d[v] + 1:
                    sigma[w] += sigma[v]; P[w].append(v)
        delta = [0.0]*N
        while S:
            w = S.pop()
            for v in P[w]:
                delta[v] += (sigma[v]/sigma[w])*(1.0+delta[w])
            if w != s:
                C[w] += delta[w]
    return [c/2.0 for c in C]

def _minmax(vals):
    lo, hi = min(vals), max(vals)
    if hi - lo < 1e-12:
        return [0.5]*len(vals)
    return [(v-lo)/(hi-lo) for v in vals]

def _priority_from_matrix(M):
    """给定交互矩阵，返回每 SDG 的完整度量 + 排名。"""
    pos = [[max(M[i][j],0) for j in range(N)] for i in range(N)]
    neg = [[max(-M[i][j],0) for j in range(N)] for i in range(N)]
    posdeg = [sum(pos[i]) for i in range(N)]
    negdeg = [sum(neg[i]) for i in range(N)]
    net    = [posdeg[i]-negdeg[i] for i in range(N)]
    nlinks = [sum(1 for j in range(N) if M[i][j]!=0) for i in range(N)]
    eig = _eigenvector(pos)
    adj = [set(j for j in range(N) if M[i][j]!=0) for i in range(N)]
    bet = _betweenness(adj)
    en = _minmax(eig); pn = _minmax(posdeg); nn = _minmax(net)
    prio = [W_EIG*en[i] + W_POS*pn[i] + W_NET*nn[i] for i in range(N)]
    order = sorted(range(N), key=lambda i: prio[i], reverse=True)
    rank = [0]*N
    for r, i in enumerate(order):
        rank[i] = r + 1
    return {
        "posdeg": posdeg, "negdeg": negdeg, "net": net, "nlinks": nlinks,
        "eig": eig, "bet": bet, "prio": prio, "rank": rank, "order": order,
    }

# ---------------------------------------------------------------------------
# 3. 10 年协同扩散有效性
# ---------------------------------------------------------------------------
def _diffusion_effectiveness(M, years=10, alpha=0.5):
    """以正协同（行归一）为传播图，投资 SDG p 后 10 年 cascade 受益总量。"""
    row = [[0.0]*N for _ in range(N)]
    for i in range(N):
        s = sum(max(M[i][j],0) for j in range(N))
        if s > 0:
            for j in range(N):
                row[i][j] = max(M[i][j],0)/s
    eff = []
    for p in range(N):
        b = [0.0]*N; b[p] = 1.0
        total = sum(b)
        for _ in range(years):
            nb = b[:]
            for j in range(N):
                inc = alpha*(1.0-b[j])*sum(row[i][j]*b[i] for i in range(N))
                nb[j] = min(1.0, b[j]+inc)
            b = nb
            total += sum(b)
        eff.append(total)
    en = _minmax(eff)
    order = sorted(range(N), key=lambda i: eff[i], reverse=True)
    rank = [0]*N
    for r,i in enumerate(order): rank[i] = r+1
    return {"eff": eff, "norm": en, "rank": rank, "order": order}

def _diffusion_series(M, p, years=10, alpha=0.5):
    """投资 SDG p 后，每年全网累计受益总量（用于 cascade 曲线）。"""
    row = [[0.0]*N for _ in range(N)]
    for i in range(N):
        s = sum(max(M[i][j],0) for j in range(N))
        if s > 0:
            for j in range(N):
                row[i][j] = max(M[i][j],0)/s
    b = [0.0]*N; b[p] = 1.0
    series = [sum(b)]
    for _ in range(years):
        nb = b[:]
        for j in range(N):
            inc = alpha*(1.0-b[j])*sum(row[i][j]*b[i] for i in range(N))
            nb[j] = min(1.0, b[j]+inc)
        b = nb
        series.append(sum(b))
    return series

def _all_diffusion_series(M, years=10):
    return {p: _diffusion_series(M, p, years) for p in range(N)}

# ---------------------------------------------------------------------------
# 4. 危机冲击
# ---------------------------------------------------------------------------
_SHOCKS = {
    "tech":   [(9,2),(7,1),(4,1),(17,1)],                       # 技术进步：强化产业/能源/教育
    "pand":   [(3,-2),(1,-1),(8,-1),(4,-1),(10,-1)],            # 大流行：重创健康，连累贫困/就业/教育
    "climate":[(2,-1),(6,-1),(11,-1),(13,-1),(14,-1),(15,-1)],  # 气候变化：多地/资源承压
    "war":    [(16,-2),(1,-1),(2,-1),(10,-1),(17,-1)],          # 战争/难民：冲击机构/贫困/平等
}

def _apply_shock(M, shock):
    M2 = [row[:] for row in M]
    for i,v in _SHOCKS[shock]:
        a = i-1
        for j in range(N):
            if M2[a][j] != 0:
                M2[a][j] = max(-3, min(3, M2[a][j]+v))
                M2[j][a] = M2[a][j]
    return M2

# ---------------------------------------------------------------------------
# 5. 企业 ESG 子图
# ---------------------------------------------------------------------------
ESG_SUB = [6,7,8,9,11,12,13,17]  # 水/能源/就业/产业/城市/消费/气候/伙伴

def _subgraph_rank(M, sub):
    n = len(sub)
    idx = {sdg-1: k for k, sdg in enumerate(sub)}
    Ms = [[0]*n for _ in range(n)]
    for a in range(n):
        for b in range(n):
            Ms[a][b] = M[sub[a]-1][sub[b]-1]
    pos = [[max(Ms[i][j],0) for j in range(n)] for i in range(n)]
    posdeg = [sum(pos[i]) for i in range(n)]
    net = [posdeg[i]-sum(max(-Ms[i][j],0) for j in range(n)) for i in range(n)]
    eig = _eigenvector(pos)
    en = _minmax(eig); pn = _minmax(posdeg); nn = _minmax(net)
    prio = [W_EIG*en[i]+W_POS*pn[i]+W_NET*nn[i] for i in range(n)]
    order = sorted(range(n), key=lambda i: prio[i], reverse=True)
    return [sub[i] for i in order], prio

# ---------------------------------------------------------------------------
# 6. 主生成
# ---------------------------------------------------------------------------
def gen_icm2023d(seed=SEED, out_csv="icm2023d.csv"):
    M = _build_matrix()
    base = _priority_from_matrix(M)
    diff = _diffusion_effectiveness(M)

    # What-if：依次"实现"每个 SDG，看其余节点中心性损失最大的"受害者"
    biggest_loser = {}
    biggest_loss = {}
    for k in range(N):
        idxs = [i for i in range(N) if i != k]
        Ms = [[M[i][j] for j in idxs] for i in idxs]
        sub = _priority_from_matrix_sub(Ms, len(idxs))
        # 被删节点 k 的邻居中，删点后特征向量中心性下降最多者
        loss = {}
        for pos, i in enumerate(idxs):
            loss[i] = base["eig"][i] - sub["eig"][pos]
        worst = max(loss, key=lambda i: loss[i])
        biggest_loser[k] = worst
        biggest_loss[k] = loss[worst]

    # 危机冲击排名漂移
    shocks = {}
    for name in _SHOCKS:
        Ms = _apply_shock(M, name)
        r = _priority_from_matrix(Ms)
        # 与前 5 对比，记录排名变化绝对值之和（前 5 节点）
        drift = sum(abs(r["rank"][i]-base["rank"][i]) for i in range(N))
        top5_before = base["order"][:5]
        top5_after = r["order"][:5]
        shocks[name] = {
            "top5_before": top5_before, "top5_after": top5_after,
            "drift": drift, "order": r["order"], "prio": r["prio"],
            "rank": r["rank"],
        }

    # ESG 子图排名
    esg_order, esg_prio = _subgraph_rank(M, ESG_SUB)

    # 汇总权威数字
    D = {
        "N": N,
        "M": M,
        "base": base,
        "order": base["order"],
        "rank": base["rank"],
        "prio": base["prio"],
        "eig": base["eig"],
        "bet": base["bet"],
        "posdeg": base["posdeg"],
        "negdeg": base["negdeg"],
        "net": base["net"],
        "nlinks": base["nlinks"],
        "diff_eff": diff["eff"],
        "diff_order": diff["order"],
        "diff_rank": diff["rank"],
        "diff_series": _all_diffusion_series(M),
        "biggest_loser": biggest_loser,
        "biggest_loss": biggest_loss,
        "shocks": shocks,
        "esg_sub": ESG_SUB,
        "esg_order": esg_order,
        "esg_prio": esg_prio,
        "W_EIG": W_EIG, "W_POS": W_POS, "W_NET": W_NET,
    }

    # 写 CSV：17 行优先级排名表
    posdeg = base["posdeg"]; negdeg = base["negdeg"]; net = base["net"]
    nlinks = base["nlinks"]; eig = base["eig"]; bet = base["bet"]; prio = base["prio"]
    out = os.path.normpath(os.path.join(HERE, "..", "assets", "problems", "data", out_csv))
    os.makedirs(os.path.dirname(out), exist_ok=True)
    lines = ["sdg,abbr_cn,posdeg,negdeg,net,nlinks,eig,betweenness,priority,rank,diff10y_eff,diff10y_rank"]
    for i in range(N):
        lines.append("%d,%s,%.4f,%.4f,%.4f,%d,%.6f,%.6f,%.6f,%d,%.4f,%d"
                     % (i+1, NAMES[i], posdeg[i], negdeg[i], net[i], nlinks[i],
                        eig[i], bet[i], prio[i], base["rank"][i], diff["eff"][i], diff["rank"][i]))
    with open(out, "w", encoding="utf-8-sig") as f:
        f.write("\n".join(lines) + "\n")

    D["csv"] = out
    return D

def _priority_from_matrix_sub(Ms, n):
    """子图（n 节点）的优先级度量（供 what-if 删点用）。"""
    pos = [[max(Ms[i][j],0) for j in range(n)] for i in range(n)]
    neg = [[max(-Ms[i][j],0) for j in range(n)] for i in range(n)]
    posdeg = [sum(pos[i]) for i in range(n)]
    negdeg = [sum(neg[i]) for i in range(n)]
    net = [posdeg[i]-negdeg[i] for i in range(n)]
    nlinks = [sum(1 for j in range(n) if Ms[i][j]!=0) for i in range(n)]
    # 局部 eig（仅 n 维幂迭代）
    x = [1.0/n]*n
    for _ in range(2000):
        y = [sum(pos[i][j]*x[j] for j in range(n)) for i in range(n)]
        s = math.sqrt(sum(v*v for v in y)) or 1.0
        y = [v/s for v in y]
        if math.sqrt(sum((y[i]-x[i])**2 for i in range(n))) < 1e-12:
            x = y; break
        x = y
    en = _minmax(x); pn = _minmax(posdeg); nn = _minmax(net)
    prio = [W_EIG*en[i]+W_POS*pn[i]+W_NET*nn[i] for i in range(n)]
    order = sorted(range(n), key=lambda i: prio[i], reverse=True)
    rank = [0]*n
    for r,i in enumerate(order): rank[i] = r+1
    return {"posdeg":posdeg,"negdeg":negdeg,"net":net,"nlinks":nlinks,
            "eig":x,"prio":prio,"rank":rank,"order":order}

if __name__ == "__main__":
    d = gen_icm2023d()
    print("ICM2023D SDG 优先级网络 生成完成 ->", d["csv"])
    print("前 5 优先级 SDG:", [d["order"][i]+1 for i in range(5)])
    print("  ", [NAMES[i] for i in d["order"][:5]])
    print("后 3 优先级 SDG:", [d["order"][i]+1 for i in range(N-3,N)])
    print("特征向量中心性最高:", NAMES[d["eig"].index(max(d["eig"]))])
    print("介数中心性最高:", NAMES[d["bet"].index(max(d["bet"]))])
    print("10年扩散有效性最高(最适合优先投资):", NAMES[d["diff_order"][0]], "=SDG", d["diff_order"][0]+1)
    print("企业 ESG 子图优先级:", [NAMES[i-1] for i in d["esg_order"]])
    for name in d["shocks"]:
        s = d["shocks"][name]
        print("  危机[%s] 排名漂移=%d  新前5=%s" % (name, s["drift"], [i+1 for i in s["top5_after"]]))