# 조치원 상권 분석용 독립 모듈!

import os

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from sklearn.neighbors import BallTree
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.metrics import roc_auc_score, classification_report

from common import *


# 4단계 : 조치원읍 골목상권 분석
# 입지계수(LQ), 골목 덩어리 찾기, 업종별 소멸률, 골목 vs 대로, 핵심 골목 2곳 비교

def 조치원분석(now, prev):
    # 조치원읍 골목상권을 자세히 분석 (보고서 4장, 활성화 방안의 근거 대부분이 여기서 나옴)
    #   입지계수 : 조치원에 많은/부족한 업종
    #   DBSCAN  : 골목 단위 덩어리 찾고 덩어리별 지표 계산
    #   업종별 소멸률, 골목 vs 대로, 핵심 골목 2곳 비교
    # 반환 파라미터: T(조치원 올해), Tp(조치원 작년), alley(골목 덩어리별 지표 표)
    log("\n" + "=" * 70)
    log(f"5. 대상 지역 정밀 분석: {TARGET}")
    log("=" * 70)

    T = now[now["행정동명"] == TARGET].copy()     # 조치원 올해
    Tp = prev[prev["행정동명"] == TARGET].copy()  # 조치원 작년
    new_city = now[now["권역"] == "동(신도심)"]
    log(f"{TARGET} 업소 {len(T):,}개 (세종 전체의 {len(T)/len(now):.1%}), "
        f"1년 소멸률 {Tp['소멸'].mean():.1%}, 체인 비율 {T['체인'].mean():.1%} "
        f"(신도심 체인 비율 {new_city['체인'].mean():.1%})")

    # '철물, 주점, 세탁, 식료품은 많고 교육, 의원은 부족'이 나옴 -> 방안에 작성 예정
    # 입지계수(LQ) = 조치원에서 그 업종 비율 / 세종 전체에서 그 업종 비율
    # 1보다 크면 조치원에 많은 업종, 작으면 적은 업종
    조치원비율 = T["상권업종중분류명"].value_counts(normalize=True)
    세종비율 = now["상권업종중분류명"].value_counts(normalize=True)
    lq = pd.DataFrame({"LQ": (조치원비율 / 세종비율).dropna()})
    lq["조치원업소수"] = T["상권업종중분류명"].value_counts()
    lq["세종업소수"] = now["상권업종중분류명"].value_counts()
    lq.index.name = "상권업종중분류명"
    lq = lq[lq["세종업소수"] >= 50]   # 업소 수 너무 적은 업종은 빼기
    lq = lq.sort_values("LQ", ascending=False)
    log("\n[입지계수 상위(조치원 특화 업종)]")
    log(lq.head(12).round(2).to_string())
    log("\n[입지계수 하위(조치원에 상대적으로 부족한 업종)]")
    log(lq.tail(12).round(2).to_string())
    lq.round(3).to_csv(os.path.join(결과폴더, "table_조치원_입지계수.csv"), encoding="utf-8-sig")

    보기 = pd.concat([lq.head(10), lq.tail(10)])
    색 = []
    for v in 보기["LQ"]:
        if v >= 1:
            색.append(COLORS[1])
        else:
            색.append(COLORS[0])
    plt.figure(figsize=(7, 7))
    plt.barh(보기.index[::-1], 보기["LQ"][::-1], color=색[::-1])
    plt.axvline(1, color="gray")
    plt.title(f"{TARGET} 업종별 입지계수(LQ) 상·하위 10\n(1보다 크면 세종 평균보다 많음)")
    그림저장("08_조치원_입지계수.png")

    # DBSCAN으로 골목상권 덩어리 찾기
    # 행정동 단위로는 골목별 차이가 안 보여서, 가게들이 붙어 있는 구역을 좌표로 찾음
    # K-Means와 달리 군집 수를 미리 안 정해도 되고, 띄엄띄엄 있는 가게는 -1(산재)로 빠져서 골목 찾기에 맞음
    # eps = 60m (걸어서 1분 정도), 최소 15개
    # (40, 80, 100m 등도 해봤는데 60m가 제일 적당했음 그냥 이걸로 하셈)
    coords = T[["x", "y"]].values
    db = DBSCAN(eps=60, min_samples=15)
    db.fit(coords)
    T["골목군집"] = db.labels_   # -1은 어느 덩어리에도 안 들어간 업소
    골목수 = len(set(db.labels_)) - (1 if -1 in db.labels_ else 0)
    log(f"\nDBSCAN 결과: 골목상권 덩어리 {골목수}개, 군집 외(산재) 업소 {(db.labels_ == -1).sum()}개")

    # 작년 업소는 군집 번호가 없어서 제일 가까운 올해 업소(60m 안)의 번호를 따라감
    tree = BallTree(coords)
    dist, idx = tree.query(Tp[["x", "y"]].values, k=1)
    Tp["골목군집"] = np.where(dist[:, 0] <= 60, T["골목군집"].values[idx[:, 0]], -1)

    # 골목 덩어리마다 지표 계산
    # 덩어리 하나씩 돌면서 크기, 업종 비율, 소멸률, 신규비율, 청년소비형/생활밀착형 비율을 구해서 표로 만듦
    # 여기서 '행복길 = 음식 54%, 소멸 17%', '새내로 = 안정적이지만 신규 11%'가 나옴 -> 방안에 활용ㄱㄱ
    YOUTH = ["비알코올", "주점", "유원지·오락"]   # 카페, 주점, 오락
    DAILY = ["이용·미용", "식료품 소매", "섬유·의복·신발 소매", "세탁", "종합 소매"]   # 단순 소매점들
    골목목록 = sorted(T.loc[T["골목군집"] >= 0, "골목군집"].unique())
    행들 = []
    for k in 골목목록:
        a = T[T["골목군집"] == k]    # 올해 이 골목 업소
        b = Tp[Tp["골목군집"] == k]  # 작년 이 골목 업소
        행 = {}
        행["골목군집"] = k
        행["업소수"] = len(a)
        행["대표도로"] = ", ".join(a["도로"].value_counts().index[:3])
        행["음식비율"] = (a["상권업종대분류명"] == "음식").mean()
        행["소매비율"] = (a["상권업종대분류명"] == "소매").mean()
        행["수리개인비율"] = (a["상권업종대분류명"] == "수리·개인").mean()
        행["체인비율"] = a["체인"].mean()
        행["신규비율"] = a["신규"].mean()
        행["위도"] = a["위도"].mean()
        행["경도"] = a["경도"].mean()
        행["소멸률"] = b["소멸"].mean()
        행["업종다양성"] = 다양성(a["상권업종중분류명"])
        행["주요업종"] = ", ".join(a["상권업종중분류명"].value_counts().index[:4])
        행["업소수_2025"] = len(b)
        행["순증가율"] = len(a) / len(b) - 1
        행["청년소비형비율"] = a["상권업종중분류명"].isin(YOUTH).mean()
        행["생활밀착형비율"] = a["상권업종중분류명"].isin(DAILY).mean()
        행들.append(행)
    alley = pd.DataFrame(행들).set_index("골목군집")
    alley = alley.sort_values("업소수", ascending=False)
    log("\n[조치원 골목상권 덩어리별 지표]")
    log(alley.drop(columns=["위도", "경도"]).round(3).to_string())
    alley.round(4).to_csv(os.path.join(결과폴더, "table_조치원_골목군집.csv"), encoding="utf-8-sig")

    # 조치원 업종(중분류)별 소멸률
    # 어떤 업종이 잘 망하고 어떤 업종이 오래 가는지 알아보는 코드
    mid = pd.DataFrame()
    mid["업소수_2025"] = Tp.groupby("상권업종중분류명").size()
    mid["소멸률"] = Tp.groupby("상권업종중분류명")["소멸"].mean()
    mid = mid.join(T.groupby("상권업종중분류명").size().rename("업소수_2026"), how="outer")
    mid = mid.join(T.groupby("상권업종중분류명")["신규"].sum().rename("신규수"), how="outer")
    mid = mid.fillna(0)
    mid["순증감"] = mid["업소수_2026"] - mid["업소수_2025"]
    mid = mid[mid["업소수_2025"] >= 20]   # 작년에 20개 이상 있던 업종만
    mid = mid.sort_values("소멸률", ascending=False)
    log("\n[조치원 업종별 1년 소멸률(2025년 업소 20개 이상 업종)]")
    log(mid.round(3).to_string())
    mid.round(4).to_csv(os.path.join(결과폴더, "table_조치원_업종별_소멸률.csv"), encoding="utf-8-sig")

    # 골목(~길로 되어있는거) vs 대로(~로 단, 다솜로 같이 작게 되어 있는 것도 있어서 주의하세용) 비교
    T["도로유형"] = np.where(T["도로"].str.endswith("길"), "길(골목)", "로(대로)")
    Tp["도로유형"] = np.where(Tp["도로"].str.endswith("길"), "길(골목)", "로(대로)")
    road = pd.DataFrame()
    road["업소수"] = T.groupby("도로유형").size()
    road["체인비율"] = T.groupby("도로유형")["체인"].mean()
    road["신규비율"] = T.groupby("도로유형")["신규"].mean()
    road["소멸률"] = Tp.groupby("도로유형")["소멸"].mean()
    road["음식비율"] = (T["상권업종대분류명"] == "음식").groupby(T["도로유형"]).mean()
    log("\n[조치원 도로유형별 비교]")
    log(road.round(3).to_string())

    # 제일 큰 골목 2개 자세히 보기 (어떤 업종이 없어지고 생겼는지)
    # 행복길에서 없어진 업종과 생긴 업종이 거의 똑같다는 걸 여기서 확인함
    core = alley.index[:2].tolist()
    for k in core:
        a = T[T["골목군집"] == k]
        b = Tp[Tp["골목군집"] == k]
        일층 = a.loc[a["1층여부"] >= 0, "1층여부"].mean()
        log(f"\n[핵심 골목 {k}: {alley.loc[k, '대표도로']}] 업소 {len(a)}개, "
            f"1년 소멸 {b['소멸'].sum()}개 / 신규 {a['신규'].sum()}개, 1층 비율 {일층:.0%}")
        없어진업종 = b[b["소멸"] == 1]["상권업종중분류명"].value_counts().head(6)
        생긴업종 = a[a["신규"] == 1]["상권업종중분류명"].value_counts().head(6)
        log("  - 소멸 상위 업종: " + ", ".join(f"{i}({v})" for i, v in 없어진업종.items()))
        log("  - 신규 상위 업종: " + ", ".join(f"{i}({v})" for i, v in 생긴업종.items()))

    # 핵심 골목 2개vs조치원 전체vs동지역 전체 비교표
    nc_prev = prev[~prev["행정동명"].str.endswith(("읍", "면"))]   # 작년 동지역
    cmp = pd.DataFrame(columns=["1년 소멸률", "신규 비율", "체인 비율", "청년소비형 비율"])
    for k in core:
        이름 = f"골목{k}\n({alley.loc[k, '대표도로'].split(',')[0]})"
        cmp.loc[이름] = [alley.loc[k, "소멸률"], alley.loc[k, "신규비율"], alley.loc[k, "체인비율"], alley.loc[k, "청년소비형비율"]]
    cmp.loc[f"{TARGET} 전체"] = [Tp["소멸"].mean(), T["신규"].mean(), T["체인"].mean(),
                               T["상권업종중분류명"].isin(YOUTH).mean()]
    cmp.loc["신도심(동) 전체"] = [nc_prev["소멸"].mean(), new_city["신규"].mean(), new_city["체인"].mean(),
                              new_city["상권업종중분류명"].isin(YOUTH).mean()]
    cmp = cmp.astype(float)
    log("\n[핵심 골목 비교표]")
    log(cmp.round(3).to_string())

    plt.figure(figsize=(12, 3.4))
    for n, col in enumerate(cmp.columns):
        plt.subplot(1, 4, n + 1)
        plt.barh(cmp.index[::-1], cmp[col][::-1], color=[COLORS[0], COLORS[0], COLORS[1], COLORS[1]])
        for i, v in enumerate(cmp[col][::-1]):
            plt.text(v, i, f" {v:.0%}", va="center", fontsize=8)
        plt.title(col)
        plt.xticks([])
        if n > 0:
            plt.yticks([])
    그림저장("10_핵심골목_비교.png")

    # 조치원 업종별 소멸률 그래프
    mm = mid.sort_values("소멸률")
    plt.figure(figsize=(7, 6))
    plt.barh(mm.index, mm["소멸률"], color=COLORS[0])
    plt.axvline(Tp["소멸"].mean(), color="gray", linestyle=":")
    plt.title(f"{TARGET} 업종별 1년 소멸률 (점선 = 조치원 평균)")
    그림저장("11_조치원_업종별_소멸률.png")

    return T, Tp, alley


# 머신러닝 - 랜덤포레스트로 폐업(소멸) 위험 보기
# 작년(2025.6) 업소가 1년 뒤에 없어졌는지(1) 남아있는지(0) 예측
# 매출 같은 정보가 없어서 정확도는 낮을 것 같음

def 소멸예측(prev, Tp, alley):
    # 랜덤포레스트로 1년 뒤 소멸 여부를 예측하고, 어떤 조건이 소멸과 관련 있는지 봄
    # 출력: 그림 09, table_조치원_골목군집.csv
    log("\n" + "=" * 70)
    log("6. 머신러닝 ②: 업소 소멸 위험 분류 (Random Forest)")
    log("=" * 70)

    P_ = prev.copy()

    # 모델 입력 변수 만들기
    # 작년 업소 기준이라 전처리에서 now 에 만든 것들을 prev 에도 똑같이 만들어 줌
    # 1층 여부 (작년 데이터도 같음!)
    일층 = []
    for 층 in P_["층정보"]:
        if pd.isna(층):
            일층.append(-1)
        elif str(층).strip() == "1":
            일층.append(1)
        else:
            일층.append(0)
    P_["1층여부"] = 일층

    이름개수2 = P_["상호명"].value_counts()
    P_["체인"] = ((P_["지점명"].notna()) | (P_["상호명"].map(이름개수2) >= 3)).astype(int)
    P_["원도심"] = P_["행정동명"].str.endswith(("읍", "면")).astype(int)
    P_["골목길"] = P_["도로"].astype(str).str.endswith("길").astype(int)

    # 주변 100m 안에 업소가 몇 개 있는지 (자기 자신은 빼려고 -1)
    tree = BallTree(P_[["x", "y"]].values)
    P_["주변100m_업소수"] = tree.query_radius(P_[["x", "y"]].values, r=100, count_only=True) - 1

    # 주변 100m 안에 같은 업종(중분류)이 몇 개 있는지 => 경쟁 보려고 
    same = np.zeros(len(P_), dtype=int)
    for 업종, 그룹 in P_.groupby("상권업종중분류명"):
        t = BallTree(그룹[["x", "y"]].values)
        개수 = t.query_radius(그룹[["x", "y"]].values, r=100, count_only=True) - 1
        same[P_.index.get_indexer(그룹.index)] = 개수
    P_["주변100m_동종업소수"] = same
    P_["동종비중"] = P_["주변100m_동종업소수"] / (P_["주변100m_업소수"] + 1)

    # 업종은 글자라서 0, 1로 바꿈
    업종더미 = pd.get_dummies(P_["상권업종대분류명"], prefix="업종").astype(int)
    num_cols = ["1층여부", "체인", "원도심", "골목길", "주변100m_업소수", "주변100m_동종업소수", "동종비중"]
    Xm = pd.concat([P_[num_cols], 업종더미], axis=1)
    y = P_["소멸"].values

    # 학습이랑 평가
    # 소멸이 13.7%뿐이라 그냥 정확도는 의미가 없어서 ROC-AUC로 평가
    # 로지스틱 회귀도 비교해 봤는데 AUC 비슷했음
    # 학습 75%, 검증 25%로 나누기
    X_tr, X_te, y_tr, y_te = train_test_split(Xm, y, test_size=0.25, stratify=y, random_state=SEED)

    # 소멸이 13%밖에 안 돼서 class_weight="balanced" 넣음 힘들다....
    rf = RandomForestClassifier(n_estimators=400, max_depth=10, min_samples_leaf=20,
                                class_weight="balanced", random_state=SEED, n_jobs=-1)
    rf.fit(X_tr, y_tr)
    proba = rf.predict_proba(X_te)[:, 1]   # 소멸할 확률
    auc = roc_auc_score(y_te, proba)
    log(f"학습 {len(X_tr):,}건 / 검증 {len(X_te):,}건, 소멸 비율 {y.mean():.1%}")
    log(f"검증 ROC-AUC = {auc:.3f} (0.5=무작위)")
    예측 = (proba >= 0.5).astype(int)
    log(classification_report(y_te, 예측, target_names=["유지", "소멸"], digits=3))

    # 해석
    # 순열 중요도로 중요한 변수를 찾고, 그 변수별 실제 소멸률을 직접 계산해서 맞는지 다시 확인
    # 변수 중요도 (순열 중요도: 변수 하나를 섞었을 때 성능이 얼마나 떨어지는지)
    pi = permutation_importance(rf, X_te, y_te, scoring="roc_auc", n_repeats=10, random_state=SEED, n_jobs=-1)
    imp = pd.Series(pi.importances_mean, index=Xm.columns).sort_values(ascending=False)
    log("[순열 중요도 상위 10]")
    log(imp.head(10).round(4).to_string())

    top = imp.head(10)
    plt.figure(figsize=(7, 4.5))
    plt.barh(top.index[::-1], top.values[::-1], color=COLORS[0])
    plt.title(f"소멸 위험 모델 변수 중요도 (순열 중요도, AUC={auc:.2f})")
    그림저장("09_소멸위험_변수중요도.png")

    # 중요하다고 나온 변수들이 실제로도 소멸률 차이가 있는지 직접 확인
    P_["주변밀도구간"] = pd.cut(P_["주변100m_업소수"], [-1, 5, 20, 50, 1000], labels=["0~5", "6~20", "21~50", "51+"])
    P_["동종비중구간"] = pd.cut(P_["동종비중"], [-0.01, 0.05, 0.15, 1], labels=["5%이하", "5~15%", "15%초과"])
    log("\n[조건별 실제 1년 소멸률]")
    for c in ["상권업종대분류명", "체인", "1층여부", "원도심", "골목길", "주변밀도구간", "동종비중구간"]:
        t = P_.groupby(c, observed=True)["소멸"].agg(["mean", "size"])
        글 = []
        for k, v in t.iterrows():
            글.append(f"{k}={v['mean']:.1%}(n={int(v['size'])})")
        log(f"- {c}: " + ", ".join(글))

    # 완료!
    P_["예측위험"] = rf.predict_proba(Xm)[:, 1]
    Tp["예측위험"] = P_.loc[Tp.index, "예측위험"]
    alley["평균예측위험"] = Tp[Tp["골목군집"] >= 0].groupby("골목군집")["예측위험"].mean()
    log("\n[조치원 골목군집별 평균 예측 위험도]")
    log(alley[["대표도로", "업소수", "소멸률", "평균예측위험"]].round(3).to_string())
    alley.round(4).to_csv(os.path.join(결과폴더, "table_조치원_골목군집.csv"), encoding="utf-8-sig")