# 세종시 전체 분석

import os

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.decomposition import PCA

from common import *


# 탐색적 데이터 분석
# 업종 분포, 행정동별 지표, 행정동, 업종 구성비, 원도심/신도심

def 탐색분석(now, prev, CATS):
    # 세종 전체의 업종 구조를 숫자와 그림으로 정리
    log("\n" + "=" * 70)
    log("2. 탐색적 데이터 분석(EDA)")
    log("=" * 70)

    # 업종 대분류별 업소 수
    # 세종 전체에서 어떤 업종이 많은지 먼저 확인
    업종수 = now["상권업종대분류명"].value_counts()
    표 = pd.DataFrame({"업소수": 업종수})
    표["비율"] = (표["업소수"] / 표["업소수"].sum()).round(3)
    표.index.name = "상권업종대분류명"
    log("\n[업종 대분류 분포]")
    log(표.to_string())

    plt.figure(figsize=(7, 4))
    plt.barh(업종수.index[::-1], 업종수.values[::-1], color=COLORS[0])
    plt.title("세종시 업종 대분류별 업소 수 (2026.6)")
    plt.xlabel("업소 수")
    그림저장("01_업종대분류_분포.png")

    # 행정동별 지표
    # 동마다 상권 크기(업소수), 활력(신규비율), 위험(소멸률), 성격(체인, 1층, 다양성)을 한 표로 만듦
    dong = pd.DataFrame()
    dong["업소수"] = now.groupby("행정동명").size()
    dong["권역"] = now.groupby("행정동명")["권역"].first()
    dong["신규비율"] = now.groupby("행정동명")["신규"].mean()
    dong["체인비율"] = now.groupby("행정동명")["체인"].mean()
    dong = dong.sort_values("업소수", ascending=False)
    dong["소멸률"] = prev.groupby("행정동명")["소멸"].mean()
    dong["순증가율"] = now.groupby("행정동명").size() / prev.groupby("행정동명").size() - 1
    층있음 = now[now["1층여부"] >= 0]   # 층정보 있는 업소만
    dong["1층비율"] = 층있음.groupby("행정동명")["1층여부"].mean()
    dong["업종다양성"] = now.groupby("행정동명")["상권업종중분류명"].apply(다양성)
    log("\n[행정동별 주요 지표]")
    log(dong.round(3).to_string())
    dong.round(4).to_csv(os.path.join(결과폴더, "table_행정동별_지표.csv"), encoding="utf-8-sig")

    # 행정동별 업소 수 그래프 (원도심은 주황, 신도심은 파랑)
    색 = []
    for g in dong["권역"]:
        if g == "읍면(원도심)":
            색.append(COLORS[1])
        else:
            색.append(COLORS[0])
    plt.figure(figsize=(7, 6))
    plt.barh(dong.index[::-1], dong["업소수"][::-1], color=색[::-1])
    plt.title("행정동별 업소 수 (주황: 읍면 원도심, 파랑: 동 신도심)")
    그림저장("02_행정동별_업소수.png")

    # 행정동별 히트맵
    # 동마다 업종 비율을 색으로 비교
    comp = pd.crosstab(now["행정동명"], now["상권업종대분류명"], normalize="index")
    comp = comp[CATS]
    comp = comp.loc[dong.index]
    plt.figure(figsize=(9, 8))
    plt.imshow(comp.values, cmap="Blues", aspect="auto")
    plt.xticks(range(len(CATS)), CATS, rotation=45, ha="right")
    plt.yticks(range(len(comp)), comp.index)
    for i in range(comp.shape[0]):
        for j in range(comp.shape[1]):
            글자색 = "black"
            if comp.values[i, j] > 0.25:   # 진한 칸은 글자 안 보여서 흰색으로
                글자색 = "white"
            plt.text(j, i, round(comp.values[i, j] * 100), ha="center", va="center", fontsize=7, color=글자색)
    plt.title("행정동별 업종 대분류 구성비(%)")
    plt.colorbar()
    그림저장("03_행정동_업종구성_히트맵.png")
    comp.round(4).to_csv(os.path.join(결과폴더, "table_행정동별_업종구성비.csv"), encoding="utf-8-sig")

    # 원도심 vs 신도심 업종 구성
    # 을 두 묶음으로 합쳐서 차이를 숫자로 확인
    grp = pd.crosstab(now["상권업종대분류명"], now["권역"], normalize="columns").loc[CATS]
    log("\n[권역별 업종 대분류 구성비]")
    log(grp.round(3).to_string())

    # 행정동별 소멸률, 신규비율 산점도
    # x = 얼마나 없어지나, y = 얼마나 새로 생기나
    plt.figure(figsize=(7, 5))
    신도심 = dong[dong["권역"] == "동(신도심)"]
    원도심 = dong[dong["권역"] == "읍면(원도심)"]
    plt.scatter(신도심["소멸률"], 신도심["신규비율"], s=np.sqrt(신도심["업소수"]) * 6, color=COLORS[0], label="동(신도심)")
    plt.scatter(원도심["소멸률"], 원도심["신규비율"], s=np.sqrt(원도심["업소수"]) * 6, color=COLORS[1], label="읍면(원도심)")
    for 이름 in dong.index:
        if dong.loc[이름, "업소수"] > 400:   # 큰 동만 이름 표시
            plt.text(dong.loc[이름, "소멸률"], dong.loc[이름, "신규비율"], 이름, fontsize=8)
    plt.xlabel("1년 소멸률 (2025.6 업소 중 사라진 비율)")
    plt.ylabel("신규 비율 (2026.6 업소 중 새로 생긴 비율)")
    plt.title("행정동별 소멸률과 신규 비율 (원 크기 = 업소 수)")
    plt.legend()
    그림저장("04_행정동_소멸률_신규비율.png")

    return dong, comp


# 머신러닝 - K-Means 로 상권 유형 나누기
# 세종 전체를 250m 격자로 나누고, 격자마다 업종 비율/업소 수/1층 비율/체인 비율을 구해서
# 비슷한 격자끼리 묶기

def 상권유형군집(now, prev, CATS):
    # 250m 격자 단위로 K-Means 군집화해서 상권 유형을 자동으로 나눔
    # 행정동은 너무 커서 한 동 안에 학원가, 먹자골목, 주택가가 섞임 -> 격자로 더 잘게 봄
    log("\n" + "=" * 70)
    log("3. 머신러닝 ①: 격자(250m) 단위 상권 유형 군집화 (K-Means)")
    log("=" * 70)

    # 격자 만들기
    now["gx"] = (now["x"] // GRID_M).astype(int)
    now["gy"] = (now["y"] // GRID_M).astype(int)
    prev["gx"] = (prev["x"] // GRID_M).astype(int)
    prev["gy"] = (prev["y"] // GRID_M).astype(int)

    # 격자별로 값 계산
    g = now.groupby(["gx", "gy"])
    cell = pd.DataFrame()
    cell["업소수"] = g.size()
    cell["체인비율"] = g["체인"].mean()
    cell["신규비율"] = g["신규"].mean()
    cell["위도"] = g["위도"].mean()
    cell["경도"] = g["경도"].mean()
    cell["대표동"] = g["행정동명"].agg(최빈값)
    층있음 = now[now["1층여부"] >= 0]
    cell["1층비율"] = 층있음.groupby(["gx", "gy"])["1층여부"].mean()
    cell["1층비율"] = cell["1층비율"].fillna(cell["1층비율"].median())   # 층정보 없는 격자는 중앙값으로
    업종비율 = pd.crosstab([now["gx"], now["gy"]], now["상권업종대분류명"], normalize="index")
    cell = cell.join(업종비율[CATS])
    cell = cell[cell["업소수"] >= MIN_STORES].copy()   # 업소 적은 격자 빼기
    cell["log업소수"] = np.log1p(cell["업소수"])   # 업소 수 차이가 너무 커서 로그 씌움
    log(f"분석 대상 격자: {len(cell)}개 (업소 {MIN_STORES}개 이상), 포함 업소 {cell['업소수'].sum():,}개")

    # K-Means
    # 변수 단위가 다르니까 표준화한 뒤 군집
    # 군집 수 k는 엘보 + 실루엣으로 봄. 실루엣이 전체적으로 낮아서(0.1대) 경계가 뚜렷하진 않음
    FEATS = CATS + ["log업소수", "1층비율", "체인비율"]
    X = StandardScaler().fit_transform(cell[FEATS])   # 표준화 (단위 맞추기)

    # k 정하기 : k를 3~8까지 해보고 점수 확인
    ks = [3, 4, 5, 6, 7, 8]
    inertias = []
    sils = []
    for k in ks:
        km = KMeans(n_clusters=k, n_init=20, random_state=SEED)
        km.fit(X)
        inertias.append(km.inertia_)
        sils.append(silhouette_score(X, km.labels_))
    log("k별 실루엣 점수: " + ", ".join(f"k={k}:{s:.3f}" for k, s in zip(ks, sils)))

    # k=3은 유형이 너무 적어서 4 이상 중에 실루엣 제일 높은 걸로
    BEST_K = 4
    최고점수 = -1
    for k, s in zip(ks, sils):
        if k >= 4 and s > 최고점수:
            최고점수 = s
            BEST_K = k
    log(f"선택된 군집 수 K = {BEST_K} (k>=4 중 실루엣 최대)")

    plt.figure(figsize=(10, 3.6))
    plt.subplot(1, 2, 1)
    plt.plot(ks, inertias, marker="o")
    plt.title("엘보 곡선")
    plt.xlabel("군집 수 k")
    plt.subplot(1, 2, 2)
    plt.plot(ks, sils, marker="o")
    plt.axvline(BEST_K, color="gray", linestyle=":")
    plt.title("실루엣 점수")
    plt.xlabel("군집 수 k")
    그림저장("05_군집수_선택.png")

    # 최종 K-Means
    km = KMeans(n_clusters=BEST_K, n_init=50, random_state=SEED)
    km.fit(X)
    cell["군집"] = km.labels_

    # 군집별 평균
    prof = cell.groupby("군집")[CATS + ["업소수", "1층비율", "체인비율", "신규비율"]].mean()
    prof["격자수"] = cell.groupby("군집").size()
    전체평균 = cell[CATS].mean()

    # 군집 이름 붙이기
    # 손으로 붙이면 k나 데이터가 바뀔 때마다 고쳐야 하니까 규칙으로 함
    상위25 = cell["업소수"].quantile(0.75)
    중앙값 = cell["업소수"].quantile(0.5)
    이름들 = []
    for c in prof.index:
        r = prof.loc[c]
        if r["업소수"] >= 상위25:
            밀도 = "고밀"
        elif r["업소수"] <= 중앙값:
            밀도 = "저밀"
        else:
            밀도 = "중밀"
        배율 = (r[CATS] / 전체평균).sort_values(ascending=False)
        if 배율.iloc[0] < 1.3:
            이름 = 밀도 + " 업종혼합"
        elif 배율.iloc[1] >= 1.3:
            이름 = f"{밀도} {배율.index[0]}·{배율.index[1]} 특화"
        else:
            이름 = f"{밀도} {배율.index[0]} 특화"
        이름들.append(이름)
    prof["유형명"] = 이름들
    # 같은 이름이 있으면 뒤에 번호 붙이기
    겹침 = prof["유형명"].duplicated(keep=False)
    prof.loc[겹침, "유형명"] = prof.loc[겹침, "유형명"] + " " + prof.loc[겹침].index.astype(str)
    cell["유형명"] = cell["군집"].map(prof["유형명"])
    log("\n[군집(상권 유형) 프로파일]")
    log(prof.round(3).to_string())
    prof.round(4).to_csv(os.path.join(결과폴더, "table_상권유형_군집프로파일.csv"), encoding="utf-8-sig")

    # 동별로 어떤 유형 격자가 몇 개 있는지
    dong_type = pd.crosstab(cell["대표동"], cell["유형명"])
    log("\n[행정동별 상권 유형 격자 수]")
    log(dong_type.to_string())

    # 결과 그림
    # 2차원으로 줄여서 군집 그림 그리기
    pca = PCA(n_components=2, random_state=SEED)
    P = pca.fit_transform(X)
    plt.figure(figsize=(7, 5))
    for c in range(BEST_K):
        m = cell["군집"].values == c
        plt.scatter(P[m, 0], P[m, 1], s=20, color=COLORS[c % 8], label=prof.loc[c, "유형명"])
    plt.title(f"상권 유형 군집 결과 (PCA 2차원, 설명력 {pca.explained_variance_ratio_.sum():.0%})")
    plt.legend(fontsize=8)
    그림저장("06_군집_PCA.png")

    # 군집별 업종 배율 히트맵 (1보다 크면 평균보다 많음)
    lift = prof[CATS] / 전체평균
    plt.figure(figsize=(9, 0.6 * BEST_K + 1.5))
    plt.imshow(lift.values, cmap="Blues", aspect="auto", vmin=0, vmax=3)
    plt.xticks(range(len(CATS)), CATS, rotation=45, ha="right")
    plt.yticks(range(BEST_K), prof["유형명"])
    for i in range(lift.shape[0]):
        for j in range(lift.shape[1]):
            글자색 = "black"
            if lift.values[i, j] > 1.8:
                글자색 = "white"
            plt.text(j, i, round(lift.values[i, j], 1), ha="center", va="center", fontsize=8, color=글자색)
    plt.title("상권 유형별 업종 특화 배율 (전체 평균=1.0)")
    plt.colorbar()
    그림저장("07_상권유형_업종특화배율.png")

    return cell, prof, BEST_K