# -*- coding: utf-8 -*-
"""
신인 걸그룹 '첫 음악방송 1위' 결정요인 실증분석 (2018~2024 데뷔 전수)
+ A2O LITE 사전확률 봉인

실행: python -X utf8 idol_model.py
입력: ../01_로우데이터/girlgroups_2018_2024_raw.csv
출력: outputs/ (코딩된 자료, 표 CSV, 그림 PNG, key_results.json, T7_봉인판단.csv)
작성: 조봉수 / AI(Claude) 협업 — 원자료·코드 공개
"""
import json, os, warnings
import numpy as np, pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm
from scipy import stats
import statsmodels.api as sm

warnings.filterwarnings("ignore")
HERE = os.path.dirname(os.path.abspath(__file__))
RAW = os.path.join(HERE, "..", "01_로우데이터", "girlgroups_2018_2024_raw.csv")
OUT = os.path.join(HERE, "outputs"); os.makedirs(OUT, exist_ok=True)
CUTOFF = pd.Timestamp("2026-10-15")       # 관측 기준일
WINDOW = 365                              # 주 결과변수: 데뷔 후 365일 이내 첫 1위

# ── 한글 글꼴
for p in ["/usr/share/fonts/opentype/noto/NotoSansCJK-Regular.ttc", "C:/Windows/Fonts/malgun.ttf"]:
    if os.path.exists(p):
        fm.fontManager.addfont(p); plt.rcParams["font.family"] = fm.FontProperties(fname=p).get_name(); break
plt.rcParams["axes.unicode_minus"] = False
FIRE, INK, GREY = "#C8862A", "#1F2430", "#9AA0A8"

# ── 1. 코딩 규칙 (논문 3장과 동일)
NOT_NEW = {  # 신인 그룹이 아닌 것: 기존 그룹 유닛·프로젝트·개명·해외 선데뷔·밴드
    "Fanatics-Flavor": "기존 팀의 데뷔 전 유닛", "Gugudan SeMiNa": "기존 그룹 유닛", "Oh!GG": "기존 그룹 유닛",
    "Pristin V": "기존 그룹 유닛", "WJMK": "기존 그룹 멤버 프로젝트", "Rolling Quartz": "아이돌이 아닌 밴드·실제 데뷔 2020-12",
    "Blackswan": "기존 그룹(라니아) 개명", "Red Velvet - Irene & Seulgi": "기존 그룹 유닛",
    "Refund Sisters": "예능 프로젝트(기성 가수)", "WJSN Chocome": "기존 그룹 유닛", "WJSN The Black": "기존 그룹 유닛",
    "Acid Angel from Asia": "트리플에스 데뷔 전 유닛", "Apink Chobom": "기존 그룹 유닛",
    "Got the Beat": "기성 가수 프로젝트", "Mamamoo+": "기존 그룹 유닛", "BB Girls": "기존 그룹(브레이브걸스) 개명",
    "NiziU": "일본 선데뷔(2020)", "Aria (tripleS sub-unit)": "기존 그룹 유닛",
}
VETERAN = {"Viviz", "Loossemble", "Artms"}           # 새 그룹이지만 멤버 전원이 기존 아이돌 출신
BIG4_KEYS = ["SM Entertainment", "JYP", "YG Entertainment", "HYBE"]  # 4대 기획사(하이브 계열 레이블 포함)
YG_AFFIL = {"Meovv"}                                  # 더블랙레이블(YG 관계사) — 민감도 분석용

def code(df):
    d = df.copy()
    d["debut"] = pd.to_datetime(d["debut_date"]); d["win"] = pd.to_datetime(d["first_win_date"])
    d["new_group"] = (~d["group"].isin(NOT_NEW)).astype(int)
    d["exclude_reason"] = d["group"].map(NOT_NEW).fillna("")
    d["veteran"] = d["group"].isin(VETERAN).astype(int)
    d["big4"] = d["agency"].apply(lambda a: int(any(k in str(a) for k in BIG4_KEYS)))
    d["big4_plus"] = ((d["big4"] == 1) | d["group"].isin(YG_AFFIL)).astype(int)
    d["survival"] = (d["survival_show"].str.lower() == "yes").astype(int)
    d["members"] = pd.to_numeric(d["members_at_debut"], errors="coerce")
    d["days_to_win"] = (d["win"] - d["debut"]).dt.days
    # 개명 그룹(브브걸)의 이전 이름 1위는 이 그룹 1위로 보지 않음
    d.loc[d["group"] == "BB Girls", ["days_to_win"]] = np.nan
    d["won_ever"] = d["days_to_win"].notna().astype(int)
    d["win_12m"] = ((d["days_to_win"] >= 0) & (d["days_to_win"] <= WINDOW)).astype(int)
    d["follow_days"] = (CUTOFF - d["debut"]).dt.days
    d["era"] = np.where(d["year"] <= 2021, "2018-21", "2022-24")
    return d

def km(times, events):
    """카플란-마이어 (직접 구현): 반환 (시점, 미달성 비율)"""
    t = np.asarray(times, float); e = np.asarray(events, int)
    order = np.argsort(t); t, e = t[order], e[order]
    xs, ys, s, n = [0.0], [1.0], 1.0, len(t)
    for u in np.unique(t[e == 1]):
        at_risk = (t >= u).sum(); d_ = ((t == u) & (e == 1)).sum()
        s *= 1 - d_ / at_risk; xs.append(u); ys.append(s)
    return np.array(xs), np.array(ys)

def rate_table(d, col):
    g = d.groupby(col).agg(팀수=("group", "size"), 일년내1위=("win_12m", "sum"), 언젠가1위=("won_ever", "sum"))
    g["일년내_비율"] = (g["일년내1위"] / g["팀수"]).round(3); g["언젠가_비율"] = (g["언젠가1위"] / g["팀수"]).round(3)
    return g

def fisher(d, col):
    t = pd.crosstab(d[col], d["win_12m"]).reindex(index=[0, 1], columns=[0, 1], fill_value=0)
    orr, p = stats.fisher_exact(t.values)
    return {"table": t.values.tolist(), "odds_ratio": None if not np.isfinite(orr) else round(float(orr), 2), "p": round(float(p), 4)}

def logit(d, cols, y="win_12m"):
    X = sm.add_constant(d[cols].astype(float)); m = sm.Logit(d[y], X).fit(disp=0)
    ci = m.conf_int()
    rows = [{"변수": c, "계수": round(m.params[c], 3), "표준오차": round(m.bse[c], 3), "p값": round(m.pvalues[c], 4),
             "오즈비": round(float(np.exp(m.params[c])), 2),
             "오즈비_95%하한": round(float(np.exp(ci.loc[c, 0])), 2), "오즈비_95%상한": round(float(np.exp(ci.loc[c, 1])), 2)} for c in m.params.index]
    return m, pd.DataFrame(rows), {"n": int(m.nobs), "pseudo_r2": round(float(m.prsquared), 3), "llr_p": round(float(m.llr_pvalue), 4)}

def main():
    raw = pd.read_csv(RAW); d = code(raw)
    d.drop(columns=["debut", "win"]).to_csv(os.path.join(OUT, "T1_코딩자료_전체96.csv"), index=False, encoding="utf-8-sig")
    m = d[d["new_group"] == 1].copy()
    K = {"관측기준일": str(CUTOFF.date()), "전체_행수": len(d), "제외_행수": int((d.new_group == 0).sum()), "분석표본": len(m)}

    # 기초 통계
    K["일년내1위_팀수"] = int(m.win_12m.sum()); K["일년내1위_비율"] = round(m.win_12m.mean(), 3)
    K["언젠가1위_팀수"] = int(m.won_ever.sum()); K["언젠가1위_비율"] = round(m.won_ever.mean(), 3)
    K["1위까지_중앙값일수(1위팀)"] = float(m.loc[m.won_ever == 1, "days_to_win"].median())
    K["미검증_무1위_팀수(2018-20)"] = int(((m.win_status == "unverified") & (m.won_ever == 0)).sum())

    t2 = pd.concat({"4대기획사": rate_table(m, "big4"), "서바이벌출신": rate_table(m, "survival"),
                    "기성멤버재결성": rate_table(m, "veteran"), "시기": rate_table(m, "era")})
    t2.to_csv(os.path.join(OUT, "T2_집단별_1위비율.csv"), encoding="utf-8-sig")
    K["fisher_big4"] = fisher(m, "big4"); K["fisher_survival"] = fisher(m, "survival")
    big = m[m.big4 == 1]
    K["4대기획사_팀"] = big[["group", "agency", "days_to_win"]].to_dict("records")

    # 로지스틱 회귀 — 주모형 (비4대 기획사만: 4대는 7팀 중 6팀 1위로 준분리)
    nb = m[m.big4 == 0].copy(); nb["members_c"] = nb["members"] - nb["members"].mean()
    nb["late"] = (nb["year"] >= 2022).astype(int)
    nb = nb.dropna(subset=["members"])
    mod, coef, fit = logit(nb, ["survival", "members_c", "late", "veteran"])
    coef.to_csv(os.path.join(OUT, "T3_로지스틱_비4대기획사.csv"), index=False, encoding="utf-8-sig")
    K["로지스틱_비4대"] = fit; K["로지스틱_비4대_계수"] = coef.to_dict("records")
    K["비4대_인원평균"] = round(float(m.loc[m.big4 == 0, "members"].mean()), 2)

    # 민감도: (a) 전체 표본 + big4 포함(준분리 주의) (b) 미검증 연도 제외 (2021~24만) (c) 기성멤버 재결성 제외
    sens = {}
    full = m.dropna(subset=["members"]).copy(); full["members_c"] = full["members"] - full["members"].mean()
    full["late"] = (full["year"] >= 2022).astype(int)
    try:
        _, c_a, f_a = logit(full, ["big4", "survival", "members_c", "late"]); sens["a_전체_big4포함"] = {"fit": f_a, "coef": c_a.to_dict("records")}
    except Exception as e:
        sens["a_전체_big4포함"] = {"error": repr(e)}
    v = nb[nb.year >= 2021]
    try:
        _, c_b, f_b = logit(v, ["survival", "members_c", "veteran"]); sens["b_2021-24만"] = {"fit": f_b, "coef": c_b.to_dict("records")}
    except Exception as e:
        sens["b_2021-24만"] = {"error": repr(e)}
    sens["c_비4대_무서바이벌_비기성_1위율"] = round(float(nb[(nb.survival == 0) & (nb.veteran == 0)].win_12m.mean()), 3)
    sens["d_big4_plus(Meovv포함)_1위율"] = round(float(m[m.big4_plus == 1].win_12m.mean()), 3)
    K["민감도"] = sens

    # ── A2O LITE 사전확률 (데뷔 시점 정보만 사용)
    a2o = {"big4": 0, "survival": 0, "members": 9, "late": 1, "veteran": 0}
    xa = pd.DataFrame([{"const": 1.0, "survival": 0.0, "members_c": 9 - nb["members"].mean(), "late": 1.0, "veteran": 0.0}])
    p_model = float(mod.predict(xa[["const", "survival", "members_c", "late", "veteran"]])[0])
    # 예측구간: 계수 공분산으로 선형예측 표준오차 → 로짓 척도 95%
    Xv = xa[["const", "survival", "members_c", "late", "veteran"]].values[0]
    se = float(np.sqrt(Xv @ mod.cov_params().values @ Xv)); lp = float(Xv @ mod.params.values)
    lo, hi = 1 / (1 + np.exp(-(lp - 1.96 * se))), 1 / (1 + np.exp(-(lp + 1.96 * se)))
    ref = nb[(nb.survival == 0) & (nb.veteran == 0) & (nb.year >= 2022)]
    K["A2O_LITE"] = {"입력": a2o, "모형확률_1년내1위": round(p_model, 3), "95%구간": [round(lo, 3), round(hi, 3)],
                     "비교집단_단순비율(비4대·비서바이벌·비기성·2022-24)": round(float(ref.win_12m.mean()), 3), "비교집단_팀수": len(ref),
                     "주의": "이수만·유영진 등 제작자 이력은 표본에 변수로 없음 → 모형이 반영하지 못함"}

    # 봉인 판단
    seal = pd.DataFrame([
        {"봉인ID": "A2O-1", "봉인일": "2026-10-16", "대상": "A2O LITE(데뷔 2026-10-13)",
         "명제": "데뷔 후 365일 이내(~2027-10-13) 6개 지상파·케이블 주간 음악방송 중 하나에서 첫 1위",
         "봉수_확률": round(p_model, 3), "근거": "비4대 기획사 신인 걸그룹 로지스틱 모형(2018~2024 전수)",
         "채점일": "2027-10-14", "중간점검": "2026-11-10(데뷔 4주차 지표로 갱신 확률 별도 봉인)", "결과": "", "브라이어점수": ""},
        {"봉인ID": "A2O-2", "봉인일": "2026-10-16", "대상": "A2O LITE",
         "명제": "데뷔 후 30일 이내(~2026-11-12) 첫 1위", "봉수_확률": round(float((nb.days_to_win <= 30).mean()), 3),
         "근거": "비4대 기획사 신인 중 30일 내 1위 비율(단순 기저율)", "채점일": "2026-11-13", "중간점검": "", "결과": "", "브라이어점수": ""},
    ])
    seal.to_csv(os.path.join(OUT, "T7_봉인판단.csv"), index=False, encoding="utf-8-sig")
    K["봉인"] = seal.to_dict("records")

    # ── 그림 1: 카플란-마이어 (4대 vs 비4대)
    fig, ax = plt.subplots(figsize=(8, 4.8), dpi=150)
    for flag, lab, col in [(1, "4대 기획사(7팀)", FIRE), (0, f"그 외({(m.big4==0).sum()}팀)", INK)]:
        s = m[m.big4 == flag]; t = s["days_to_win"].fillna(s["follow_days"]); e = s["won_ever"]
        x, y = km(t, e); x = np.append(x, 730); y = np.append(y, y[-1])
        ax.step(x, 1 - y, where="post", color=col, lw=2.2, label=lab)
    ax.axvline(365, color=GREY, ls="--", lw=1); ax.text(370, 0.03, "1년", color=GREY)
    ax.set_xlim(0, 730); ax.set_ylim(0, 1); ax.set_xlabel("데뷔 후 일수"); ax.set_ylabel("첫 1위를 한 팀의 누적 비율")
    ax.set_title("그림 1. 데뷔 후 첫 음악방송 1위까지 걸린 시간 (카플란-마이어)"); ax.legend(frameon=False); ax.grid(alpha=.25)
    fig.tight_layout(); fig.savefig(os.path.join(OUT, "fig1_km.png")); plt.close(fig)

    # 그림 2: 집단별 1년 내 1위 비율
    groups = [("전체", m), ("4대 기획사", m[m.big4 == 1]), ("그 외 기획사", m[m.big4 == 0]),
              ("서바이벌 출신", m[m.survival == 1]), ("서바이벌 아님", m[m.survival == 0]),
              ("2018~21 데뷔", m[m.era == "2018-21"]), ("2022~24 데뷔", m[m.era == "2022-24"])]
    fig, ax = plt.subplots(figsize=(8, 4.6), dpi=150)
    names = [g[0] for g in groups]; vals = [g[1].win_12m.mean() for g in groups]; ns = [len(g[1]) for g in groups]
    bars = ax.barh(names[::-1], vals[::-1], color=[FIRE if "4대" in n and "그 외" not in n else INK for n in names[::-1]])
    for b, v_, n_ in zip(bars, vals[::-1], ns[::-1]):
        ax.text(v_ + .01, b.get_y() + b.get_height() / 2, f"{v_*100:.0f}%  (n={n_})", va="center", fontsize=9)
    ax.set_xlim(0, 1.05); ax.set_xlabel("데뷔 1년 안에 첫 1위를 한 비율"); ax.set_title("그림 2. 집단별 '1년 내 첫 1위' 비율"); ax.grid(axis="x", alpha=.25)
    fig.tight_layout(); fig.savefig(os.path.join(OUT, "fig2_rates.png")); plt.close(fig)

    # 그림 3: A2O LITE 사전확률 위치
    fig, ax = plt.subplots(figsize=(8, 2.6), dpi=150)
    ax.axhline(0, color=GREY, lw=1)
    for lab, val, col, yy in [("비4대 신인 전체", nb.win_12m.mean(), INK, .5), ("4대 기획사 신인", m[m.big4 == 1].win_12m.mean(), FIRE, .5),
                              ("A2O LITE 봉인확률", p_model, "#B03A2E", -.5)]:
        ax.plot(val, 0, "o", ms=11, color=col); ax.annotate(f"{lab}\n{val*100:.0f}%", (val, 0), (val, yy), ha="center", fontsize=9,
                                                           textcoords="data", arrowprops=dict(arrowstyle="-", color=col))
    ax.plot([lo, hi], [0, 0], color="#B03A2E", lw=5, alpha=.35)
    ax.set_xlim(0, 1); ax.set_ylim(-1, 1); ax.set_yticks([]); ax.set_xlabel("데뷔 1년 안에 첫 1위를 할 확률")
    ax.set_title("그림 3. A2O LITE 봉인 확률의 위치 (막대 = 95% 구간)"); [ax.spines[s].set_visible(False) for s in ["left", "right", "top"]]
    fig.tight_layout(); fig.savefig(os.path.join(OUT, "fig3_a2o.png")); plt.close(fig)

    with open(os.path.join(OUT, "key_results.json"), "w", encoding="utf-8") as f:
        json.dump(K, f, ensure_ascii=False, indent=2, default=lambda o: o.item() if hasattr(o, "item") else str(o))
    print(json.dumps({k: K[k] for k in ["분석표본", "일년내1위_비율", "언젠가1위_비율", "fisher_big4", "fisher_survival", "로지스틱_비4대", "A2O_LITE"]}, ensure_ascii=False, indent=1))

if __name__ == "__main__":
    main()