# -*- coding: utf-8 -*-
"""
Chapter 2 — I-O prior 설계
 (1) 산업연관표 부문(30) → 36 GDP 산업 매핑
 (2) 생산유발계수(Leontief 역행렬, 최신연도)로 36x36 연계강도 행렬 W
 (3) 미네소타 prior를 연계강도로 차등화하는 ω 행렬
     · 교차항 사전분산 ∝ λ²/(l²σ_j²) × ω_ij
     · ω = 정규화 W (오프대각 평균 1로 → 총수축량은 미네소타와 동일, 배분만 차등), 하한 0.05
산출: ch2/data/(io_linkage_W.csv, io_prior_omega.csv, io_mapping.csv, 입력데이터 사본),
      ch2/figures/fig1_io_prior_design.png
"""
import os, glob, shutil
import numpy as np, pandas as pd
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm, TwoSlopeNorm

plt.rcParams["font.family"] = "Malgun Gothic"
plt.rcParams["axes.unicode_minus"] = False
plt.rcParams["savefig.dpi"] = 130

HERE = os.path.dirname(os.path.abspath(__file__))         # repo/ch2/code
CH2 = os.path.dirname(HERE)                                # repo/ch2
REPO = os.path.dirname(CH2)                                # repo
DATA = os.path.join(CH2, "data"); os.makedirs(DATA, exist_ok=True)
FIG = os.path.join(CH2, "figures"); os.makedirs(FIG, exist_ok=True)

# --- 입력 데이터를 ch2/data 로 복사(자기완결) ---
shutil.copy(os.path.join(REPO, "ch1", "data", "dlog_rates_sa.csv"),
            os.path.join(DATA, "dlog_rates_sa.csv"))
# 생산유발계수 파일 탐지(대각평균>1 = Leontief 역행렬)
io_src = None
for f in glob.glob(os.path.join(REPO, "data", "raw", "io_*_long.csv")):
    df = pd.read_csv(f, encoding="utf-8-sig")
    if "행코드_i" in df and df[df["행코드_i"] == df["열코드_j"]]["값"].mean() > 1:
        io_src = f; break
shutil.copy(io_src, os.path.join(DATA, "io_production_inducement_long.csv"))

# --- 36 GDP 산업 순서 ---
rate = pd.read_csv(os.path.join(DATA, "dlog_rates_sa.csv"), encoding="utf-8-sig", index_col=0)
inds = [c for c in rate.columns if c != "총량GDP"]
assert len(inds) == 36

# --- 부문 매핑 (GDP 산업 → I-O 코드) ---
MAP = {
 "농림어업": ["A"], "광업": ["B"],
 "음식료품 제조업": ["C01"], "섬유 및 가죽제품 제조업": ["C02"],
 "목재, 종이, 인쇄 및 복제업": ["C03"], "코크스 및 석유정제품 제조업": ["C04"],
 "화학물질 및 화학제품 제조업": ["C05"], "비금속광물제품 제조업": ["C06"],
 "1차금속 제조업": ["C07"], "금속가공제품 제조업": ["C08"],
 "컴퓨터, 전자 및 광학기기 제조업": ["C09"], "전기장비 제조업": ["C10"],
 "기계 및 장비 제조업": ["C11"], "운송장비 제조업": ["C12"],
 "기타 제조업 및 산업용 장비 수리업": ["C13", "C14"],
 "전기업": ["D"], "가스, 증기 및 공기조절 공급업": ["D"],
 "수도, 하수 및 폐기물 처리, 원료 재생업": ["E"],
 "주거용 건물 건설업": ["F"], "비주거용 건물 건설업": ["F"], "토목건설업": ["F"], "건축보수업": ["F"],
 "도소매업": ["G"], "숙박 및 음식점업": ["I"], "운수업": ["H"],
 "금융 및 보험업": ["K"], "부동산업": ["L"],
 "통신업": ["J"], "출판, 방송, 영상, 정보서비스": ["J"],
 "전문, 과학 및 기술관련 서비스업": ["M"], "사업지원 서비스업": ["N"],
 "공공행정, 국방 및 사회보장": ["O"], "교육서비스업": ["P"],
 "의료, 보건업 및 사회복지서비스업": ["Q"],
 "예술, 스포츠 및 여가관련 서비스업": ["R"], "기타 서비스업": ["S"],
}
missing = [i for i in inds if i not in MAP]
assert not missing, f"매핑 누락: {missing}"

# --- 생산유발계수 행렬 M (최신연도) ---
io = pd.read_csv(os.path.join(DATA, "io_production_inducement_long.csv"), encoding="utf-8-sig")
yr = io["연도"].max()
io = io[io["연도"] == yr]
M = io.pivot_table(index="행코드_i", columns="열코드_j", values="값", aggfunc="mean")

# --- 36x36 연계강도 W  (W[a,b] = j(=b) 가 i(=a) 를 유발하는 강도) ---
n = 36
W = np.zeros((n, n))
for a, ia in enumerate(inds):
    for b, ib in enumerate(inds):
        vals = [M.loc[ri, cj] for ri in MAP[ia] for cj in MAP[ib]
                if ri in M.index and cj in M.columns]
        W[a, b] = np.mean(vals) if vals else 0.0
Wdf = pd.DataFrame(W, index=inds, columns=inds)

# --- ω : 미네소타 차등화 계수 (오프대각 평균 1, 하한 0.05) ---
off = ~np.eye(n, dtype=bool)
omega = W / W[off].mean()          # 오프대각 평균 = 1
omega = np.clip(omega, 0.2, 5.0)   # 극단치 제한(안정)·순위 보존
np.fill_diagonal(omega, 1.0)        # 자기항은 표준 미네소타
Odf = pd.DataFrame(omega, index=inds, columns=inds)

# --- 저장 ---
Wdf.round(4).to_csv(os.path.join(DATA, "io_linkage_W.csv"), encoding="utf-8-sig")
Odf.round(4).to_csv(os.path.join(DATA, "io_prior_omega.csv"), encoding="utf-8-sig")
pd.DataFrame({"GDP산업": list(MAP), "IO코드": [",".join(v) for v in MAP.values()]}).to_csv(
    os.path.join(DATA, "io_mapping.csv"), index=False, encoding="utf-8-sig")

# --- 검증 출력 ---
def short(s): return (s.replace(" 제조업", "").replace(" 및 ", "·").replace("서비스업", "서비스"))[:9]
print(f"생산유발계수 기준연도: {yr}  | 부문 매핑 36/36 완료 (누락 없음)")
Woff = Wdf.where(off)
top = Woff.stack().sort_values(ascending=False).head(8)
print("\n[연계강도 상위 8쌍  (j → i)]")
for (i, j), v in top.items():
    print(f"  {v:6.2f}   {short(j)}  →  {short(i)}")
ce = Wdf.loc["전기장비 제조업", "컴퓨터, 전자 및 광학기기 제조업"]
print(f"\n검증) 컴퓨터·전자 → 전기장비 연계 = {ce:.2f}  (챕터1 예제쌍, 상호 강함 기대)")
print(f"ω 범위: {omega[off].min():.3f} ~ {omega[off].max():.3f}  (오프대각 평균 {omega[off].mean():.3f})")
print(f"ω<0.5(강수축) 쌍 비율: {(omega[off]<0.5).mean()*100:.0f}%   ω>2(느슨) 쌍 비율: {(omega[off]>2).mean()*100:.0f}%")

# --- 그림: (a) 연계강도 W  (b) prior 차등화 ω ---
lab = [short(s) for s in inds]
fig, ax = plt.subplots(1, 2, figsize=(16, 8))
fig.suptitle(f"그림 1. I-O prior 설계 — 산업연관표(생산유발계수 {yr})로 미네소타 수축을 차등화",
             fontsize=14, fontweight="bold")
Wp = W.copy(); Wp[Wp <= 0] = np.nan
im0 = ax[0].imshow(Wp, cmap="YlOrRd", norm=LogNorm(vmin=max(1e-2, np.nanmin(Wp)), vmax=np.nanmax(Wp)), aspect="auto")
ax[0].set_title("(a) 연계강도 W (생산유발, 로그색)\n열 j → 행 i 유발")
im1 = ax[1].imshow(omega, cmap="RdBu_r", norm=TwoSlopeNorm(1.0, omega[off].min(), min(4, omega[off].max())), aspect="auto")
ax[1].set_title("(b) prior 차등화 ω (붉음=느슨/연계강, 청=수축/연계약)\n미네소타 대비 교차항 배분")
for a_ in ax:
    a_.set_xticks(range(n)); a_.set_xticklabels(lab, rotation=90, fontsize=5)
    a_.set_yticks(range(n)); a_.set_yticklabels(lab, fontsize=5)
fig.colorbar(im0, ax=ax[0], shrink=.8, pad=.02); fig.colorbar(im1, ax=ax[1], shrink=.8, pad=.02)
fig.tight_layout(rect=[0, 0, 1, 0.95])
fig.savefig(os.path.join(FIG, "fig1_io_prior_design.png"), bbox_inches="tight")
print("\n저장: ch2/data/(W, ω, mapping, 입력사본), ch2/figures/fig1_io_prior_design.png")