"""
2020年度 統計数理賞（大学生・一般の部）
「地域産業構造と経済成長：主成分分析とクラスタリングによる類型化」
再現教育用スクリプト

使用データ: SSDSE-B-2026.csv（2022年断面データ、47都道府県）
手法: PCA → Ward法クラスタリング（k=4）→ クラスター別産業特性分析
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2020_U3_suri.py
# ============================================================


import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ─────────────────────────────────────────────
# データ読込・前処理
# ─────────────────────────────────────────────
df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}', na=False)].copy()
df_b['年度'] = df_b['年度'].astype(int)

# 2022年断面データ（47都道府県）
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
prefs = df['都道府県'].tolist()
print(f"分析対象: {len(df)} 都道府県（2022年度）")

# ─────────────────────────────────────────────
# 産業構造指標の構築（LQ類似指標: 全国平均=1.0に正規化）
# ─────────────────────────────────────────────
# 指標1: 高齢化率（農村・第一次産業地域の代理）
lq_agri = df['65歳以上人口'] / df['総人口']
lq_agri = lq_agri / lq_agri.mean()

# 指標2: 建設・着工密度（第二次産業・製造業系の代理）
lq_mfg = df['着工建築物数'] / df['総人口']
lq_mfg = lq_mfg / lq_mfg.mean()

# 指標3: 消費・サービス支出水準（第三次産業代理）
lq_svc = df['消費支出（二人以上の世帯）']
lq_svc = lq_svc / lq_svc.mean()

# 指標4: 医療密度（医療・福祉特化）
lq_med = (df['一般病院数'] + df['一般診療所数']) / df['総人口']
lq_med = lq_med / lq_med.mean()

# 指標5: 教育費支出（知識産業・高学歴地域代理）
lq_edu = df['教育費（二人以上の世帯）']
lq_edu = lq_edu / lq_edu.mean()

# 指標6: 観光施設密度（旅館・ホテル数/人口: 観光特化）
lq_tour = df['旅館営業施設数（ホテルを含む）'] / df['総人口']
lq_tour = lq_tour / lq_tour.mean()

# 指標行列の作成
feat_names = [
    '農村型高齢化指数',
    '建設・製造活発度',
    'サービス消費水準',
    '医療特化度',
    '教育・知識集積度',
    '観光特化度',
]
X_raw = pd.DataFrame({
    feat_names[0]: lq_agri.values,
    feat_names[1]: lq_mfg.values,
    feat_names[2]: lq_svc.values,
    feat_names[3]: lq_med.values,
    feat_names[4]: lq_edu.values,
    feat_names[5]: lq_tour.values,
}, index=prefs)

print("\nLQ類似指標（全国平均=1.0）の記述統計:")
print(X_raw.describe().round(3))

# ─────────────────────────────────────────────
# 標準化 → PCA
# ─────────────────────────────────────────────
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)

pca = PCA(n_components=3)
X_pca = pca.fit_transform(X_scaled)

explained = pca.explained_variance_ratio_
cumulative = np.cumsum(explained)
loadings = pca.components_  # shape: (3, 6)

print(f"\nPCA 累積寄与率:")
for i, (e, c) in enumerate(zip(explained, cumulative)):
    print(f"  PC{i+1}: {e*100:.1f}%  累積 {c*100:.1f}%")

print(f"\n主成分負荷量 (PC1〜PC3):")
load_df = pd.DataFrame(loadings.T, index=feat_names,
                        columns=['PC1', 'PC2', 'PC3'])
print(load_df.round(3))

# ─────────────────────────────────────────────
# Ward法クラスタリング
# ─────────────────────────────────────────────
Z = linkage(X_scaled, method='ward')

k = 4
labels = fcluster(Z, k, criterion='maxclust')
df_result = X_raw.copy()
df_result['PC1'] = X_pca[:, 0]
df_result['PC2'] = X_pca[:, 1]
df_result['PC3'] = X_pca[:, 2]
df_result['Cluster'] = labels
df_result['都道府県'] = prefs

print(f"\nWard法 k={k} クラスター割当:")
for c in sorted(df_result['Cluster'].unique()):
    prefs_in = df_result[df_result['Cluster'] == c]['都道府県'].tolist()
    print(f"  クラスター{c} ({len(prefs_in)}府県): {', '.join(prefs_in)}")

# ─────────────────────────────────────────────
# 地域マップ
# ─────────────────────────────────────────────
region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東',
    '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄'
}
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12'
}

df_result['地方'] = df_result['都道府県'].map(region_map)

# ─────────────────────────────────────────────
# 図1: PCAバイプロット（PC1 vs PC2）
# ─────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(11, 9))

regions_ordered = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
for reg in regions_ordered:
    mask = df_result['地方'] == reg
    ax.scatter(
        df_result.loc[mask, 'PC1'],
        df_result.loc[mask, 'PC2'],
        c=region_colors[reg], label=reg,
        s=70, alpha=0.85, zorder=3
    )

# 都道府県ラベル（全47）
for _, row in df_result.iterrows():
    pref_short = row['都道府県'].replace('都', '').replace('道', '').replace('府', '').replace('県', '')
    ax.annotate(pref_short, (row['PC1'], row['PC2']),
                fontsize=6.5, ha='center', va='bottom',
                xytext=(0, 4), textcoords='offset points', color='#333333')

# 変数ベクトル（biplot arrows）
scale = 2.5
for i, name in enumerate(feat_names):
    vx = loadings[0, i] * scale
    vy = loadings[1, i] * scale
    ax.annotate('', xy=(vx, vy), xytext=(0, 0),
                arrowprops=dict(arrowstyle='->', color='#CC0000', lw=1.6))
    offset_x = 0.08 if vx >= 0 else -0.08
    offset_y = 0.08 if vy >= 0 else -0.08
    ax.text(vx + offset_x, vy + offset_y, name,
            fontsize=8, color='#CC0000', fontweight='bold',
            ha='center', va='center')

ax.axhline(0, color='gray', lw=0.5, ls='--', alpha=0.5)
ax.axvline(0, color='gray', lw=0.5, ls='--', alpha=0.5)
ax.set_xlabel(f'PC1（寄与率 {explained[0]*100:.1f}%）', fontsize=12)
ax.set_ylabel(f'PC2（寄与率 {explained[1]*100:.1f}%）', fontsize=12)
ax.set_title('図1: PCAバイプロット（47都道府県 × 産業構造指標）', fontsize=13, pad=14)
ax.legend(loc='lower right', fontsize=9, framealpha=0.85)
ax.grid(True, alpha=0.2)

plt.tight_layout()
fig_path1 = os.path.join(FIG_DIR, '2020_U3_fig1.png')
plt.savefig(fig_path1, dpi=150, bbox_inches='tight')
plt.close()
print(f"\n図1 保存: {fig_path1}")

# ─────────────────────────────────────────────
# 図2: 主成分負荷量ヒートマップ（PC1〜PC3）
# ─────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(8, 5))

im = ax.imshow(loadings.T, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax, label='負荷量', shrink=0.8)

ax.set_xticks(range(3))
ax.set_xticklabels(['PC1', 'PC2', 'PC3'], fontsize=12)
ax.set_yticks(range(len(feat_names)))
ax.set_yticklabels(feat_names, fontsize=10)

for i in range(3):
    for j in range(len(feat_names)):
        val = loadings[i, j]
        text_color = 'white' if abs(val) > 0.5 else 'black'
        ax.text(i, j, f'{val:.3f}', ha='center', va='center',
                fontsize=9, color=text_color, fontweight='bold')

ax.set_title('図2: 主成分負荷量ヒートマップ（PC1〜PC3）', fontsize=13, pad=12)
plt.tight_layout()
fig_path2 = os.path.join(FIG_DIR, '2020_U3_fig2.png')
plt.savefig(fig_path2, dpi=150, bbox_inches='tight')
plt.close()
print(f"図2 保存: {fig_path2}")

# ─────────────────────────────────────────────
# 図3: Ward法デンドログラム
# ─────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(14, 6))

# カット高さの取得
heights = Z[:, 2]
sorted_h = np.sort(heights)
cut_height = (sorted_h[-k] + sorted_h[-k+1]) / 2.0

dn = dendrogram(
    Z,
    labels=prefs,
    ax=ax,
    leaf_rotation=90,
    leaf_font_size=7.5,
    color_threshold=cut_height,
)

ax.axhline(y=cut_height, color='red', lw=1.8, ls='--', label=f'k={k}カット線')
ax.set_title('図3: Ward法デンドログラム（47都道府県 産業構造クラスタリング）',
             fontsize=12, pad=10)
ax.set_ylabel('ユークリッド距離（Ward連結）', fontsize=10)
ax.set_xlabel('都道府県', fontsize=10)
ax.legend(fontsize=10)
ax.grid(axis='y', alpha=0.2)

plt.tight_layout()
fig_path3 = os.path.join(FIG_DIR, '2020_U3_fig3.png')
plt.savefig(fig_path3, dpi=150, bbox_inches='tight')
plt.close()
print(f"図3 保存: {fig_path3}")

# ─────────────────────────────────────────────
# 図4: クラスター別レーダーチャート（2×2配置）
# ─────────────────────────────────────────────
# クラスター別平均（LQ値）
cluster_means = {}
for c in sorted(df_result['Cluster'].unique()):
    mask = df_result['Cluster'] == c
    cluster_means[c] = X_raw[mask].mean()

# クラスター名の解釈（値から自動設定）
cluster_labels = {}
for c, means in cluster_means.items():
    top_feat = means.idxmax()
    cluster_labels[c] = f'クラスター{c}\n（{top_feat}型）'

N_feat = len(feat_names)
angles = np.linspace(0, 2 * np.pi, N_feat, endpoint=False).tolist()
angles += angles[:1]  # 閉じる

feat_short = ['農村型\n高齢化', '建設・\n製造', 'サービス\n消費', '医療\n特化', '教育・\n知識', '観光\n特化']

cluster_colors_radar = ['#4e9af1', '#e05c5c', '#5cb85c', '#f0a500']

fig, axes = plt.subplots(2, 2, figsize=(11, 9),
                          subplot_kw=dict(projection='polar'))
axes = axes.flatten()

for idx, c in enumerate(sorted(cluster_means.keys())):
    ax = axes[idx]
    vals = cluster_means[c].values.tolist()
    vals += vals[:1]

    ax.plot(angles, vals, color=cluster_colors_radar[idx], lw=2, marker='o', markersize=5)
    ax.fill(angles, vals, color=cluster_colors_radar[idx], alpha=0.22)

    # 全国平均（LQ=1.0）
    ax.plot(angles, [1.0] * (N_feat + 1), color='gray', lw=1, ls='--', alpha=0.6)

    ax.set_xticks(angles[:-1])
    ax.set_xticklabels(feat_short, fontsize=8.5)
    ax.set_ylim(0, 2.0)
    ax.set_yticks([0.5, 1.0, 1.5, 2.0])
    ax.set_yticklabels(['0.5', '1.0', '1.5', '2.0'], fontsize=7, color='gray')

    n_pref = (df_result['Cluster'] == c).sum()
    ax.set_title(
        f'クラスター {c}  ({n_pref}府県)\n{cluster_means[c].idxmax()} 型',
        fontsize=9, pad=12, fontweight='bold', color=cluster_colors_radar[idx]
    )

    # 代表都道府県をタイトル下に表示
    pref_list = df_result[df_result['Cluster'] == c]['都道府県'].tolist()
    pref_short_str = '・'.join(
        [p.replace('都', '').replace('道', '').replace('府', '').replace('県', '')
         for p in pref_list[:6]]
    )
    if len(pref_list) > 6:
        pref_short_str += '…'
    ax.text(0, -0.18, pref_short_str, transform=ax.transAxes,
            fontsize=6.5, ha='center', va='top', color='#555555')

fig.suptitle('図4: クラスター別産業特性レーダーチャート（全国平均=1.0）',
             fontsize=12, y=1.01, fontweight='bold')
plt.tight_layout()
fig_path4 = os.path.join(FIG_DIR, '2020_U3_fig4.png')
plt.savefig(fig_path4, dpi=150, bbox_inches='tight')
plt.close()
print(f"図4 保存: {fig_path4}")

# ─────────────────────────────────────────────
# 重回帰分析: クラスター番号 → 消費水準（経済成長代理）
# ─────────────────────────────────────────────
# クラスターダミーを使った重回帰（経済成長の代理: 消費支出水準）
y_reg = df['消費支出（二人以上の世帯）'].values.astype(float)
# PC1, PC2, PC3 を説明変数
X_reg = sm.add_constant(X_pca[:, :2])
ols_result = sm.OLS(y_reg, X_reg).fit()
print("\n重回帰分析結果（目的変数: 消費支出水準）:")
print(ols_result.summary())

# ─────────────────────────────────────────────
# 統計サマリー出力
# ─────────────────────────────────────────────
print("\n" + "=" * 60)
print("【分析サマリー】")
print(f"分析年度: 2022年度 / 対象: 47都道府県")
print(f"PCA 第1主成分 寄与率: {explained[0]*100:.1f}%")
print(f"PCA 第2主成分 寄与率: {explained[1]*100:.1f}%")
print(f"PCA 第3主成分 寄与率: {explained[2]*100:.1f}%")
print(f"PC1〜PC3 累積寄与率: {cumulative[2]*100:.1f}%")
print()
print("PC1 上位負荷量変数:")
pc1_load = pd.Series(loadings[0], index=feat_names).abs().sort_values(ascending=False)
for feat, val in pc1_load.items():
    sign = '+' if loadings[0, feat_names.index(feat)] >= 0 else '-'
    print(f"  {sign}{feat}: {val:.3f}")
print()
print("PC2 上位負荷量変数:")
pc2_load = pd.Series(loadings[1], index=feat_names).abs().sort_values(ascending=False)
for feat, val in pc2_load.items():
    sign = '+' if loadings[1, feat_names.index(feat)] >= 0 else '-'
    print(f"  {sign}{feat}: {val:.3f}")
print()
print(f"重回帰 R² = {ols_result.rsquared:.3f}, 調整済みR² = {ols_result.rsquared_adj:.3f}")
print(f"F統計量 p値 = {ols_result.f_pvalue:.4f}")
print()
print("クラスター別都道府県数:")
for c in sorted(df_result['Cluster'].unique()):
    n = (df_result['Cluster'] == c).sum()
    top = cluster_means[c].idxmax()
    print(f"  クラスター{c} ({n}府県, {top}型)")
print("=" * 60)