"""
教育用再現コード: 2024年 統計データ分析コンペティション 審査員奨励賞（大学生・一般）
=================================================================================
論文タイトル：地方創生を推進する多角的複合指標の提案
              ―SVMに基づく主観的でない変数選択と重み付け手法の検討―
受賞区分  ：審査員奨励賞 ［大学生・一般の部］
著者      ：衣川凌太（神戸大学国際人間科学部）

【分析概要】
  データ：SSDSE-B-2026.csv（47都道府県, 2022年度）
  目的   ：SVM（サポートベクターマシン）による客観的な変数選択と重み付けを通じて
            地方創生度の複合指標を構築し、都道府県間比較を行う

  Step1. SSDSE-Bから実データを読み込み、地方創生関連指標を算出
  Step2. LinearSVRで変数の重要度（重み）を算出
          目的変数：合計特殊出生率（地域活力・持続可能性の代理変数）
  Step3. 複合指標（Composite Index）の構築とランキング作成
  Step4. PCAとの比較（バイプロット）

【使用変数】
  目的変数 (y)：合計特殊出生率 (A4103)
  説明変数 (X)：
    転入率         = A5101 / A1101 × 100
    転出率         = A5102 / A1101 × 100
    保育所密度     = J2503 / A1101 × 10000
    保育定員密度   = J2505 / A1101 × 10000
    高校→大学進学率 = E4602 / E4601 × 100
    消費支出       = L3221（二人以上世帯, 円）
    住宅地標準価格 = C5401（円/m²）
    有効求人倍率   = F3103 / F3102
    高齢化率       = A1303 / A1101 × 100
    一般病院密度   = I510120 / A1101 × 10000
    宿泊者数per人  = G7101 / A1101

【図の出力】
  html/figures/2024_U5_2_fig1_svm_weights.png   ... SVM重み上位変数棒グラフ
  html/figures/2024_U5_2_fig2_index_rank.png    ... 複合指標ランキング
  html/figures/2024_U5_2_fig3_pca.png           ... PCAバイプロット
  html/figures/2024_U5_2_fig4_scatter.png       ... SVM指標 vs PCA指標散布図
=================================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2024_U5_2_shorei.py
# ============================================================


import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

from sklearn.svm import LinearSVR
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from scipy import stats as scipy_stats

# ── パス設定 ─────────────────────────────────────────────────────────────────
BASE_DIR = os.path.join(_script_dir, '..')
FIG_DIR  = os.path.join(BASE_DIR, 'html', 'figures')
DATA_PATH = os.path.join(BASE_DIR, 'data', 'raw', 'SSDSE-B-2026.csv')
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family':        'Hiragino Sans',
    'axes.unicode_minus': False,
    'figure.dpi':         150,
    'axes.spines.top':    False,
    'axes.spines.right':  False,
})

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step1. 実データ読み込み（SSDSE-B-2026, 2022年度, 47都道府県）
# ═══════════════════════════════════════════════════════════════════════════════
print("=" * 60)
print("■ Step1. SSDSE-B-2026 実データ読み込み（2022年度）")
print("=" * 60)

df_b = pd.read_csv(DATA_PATH, encoding='cp932', header=1)
# 都道府県行のみ（市区町村を除外）
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)]
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
print(f"  読み込み完了: {len(df)} 都道府県")

# 都道府県名を短縮（「県」「都」「府」「道」を保持しつつ表示用に整理）
PREFS = df['都道府県'].tolist()

# ── 変数算出 ──────────────────────────────────────────────────────────────────
pop   = df['総人口']                                      # A1101
tfr   = df['合計特殊出生率']                             # A4103  ← 目的変数
inflow    = df['転入者数（日本人移動者）']               # A5101
outflow   = df['転出者数（日本人移動者）']               # A5102
nursery_n = df['保育所等数']                             # J2503
nursery_c = df['保育所等定員数']                         # J2505
univ_grads= df['高等学校卒業者のうち進学者数']          # E4602
hs_grads  = df['高等学校卒業者数']                       # E4601
consumption= df['消費支出（二人以上の世帯）']           # L3221
land_price = df['標準価格（平均価格）（住宅地）']       # C5401
job_open   = df['月間有効求人数（一般）']                # F3103
job_seek   = df['月間有効求職者数（一般）']              # F3102
elderly    = df['65歳以上人口']                          # A1303
hospitals  = df['一般病院数']                            # I510120
tourists   = df['延べ宿泊者数']                          # G7101

# 派生指標算出
in_rate      = inflow    / pop * 100         # 転入率 (%)
out_rate     = outflow   / pop * 100         # 転出率 (%)
nursery_dens = nursery_n / pop * 10000       # 保育所密度 (10000人当たり)
nursery_cap  = nursery_c / pop * 10000       # 保育定員密度 (10000人当たり)
univ_rate    = univ_grads / hs_grads * 100   # 高校→大学進学率 (%)
job_ratio    = job_open   / job_seek         # 有効求人倍率
aging_rate   = elderly    / pop * 100        # 高齢化率 (%)
hosp_dens    = hospitals  / pop * 10000      # 一般病院密度 (10000人当たり)
tourist_pc   = tourists   / pop             # 宿泊者数per capita

# 説明変数・変数名を整理
VAR_NAMES = [
    '転入率(%)',
    '転出率(%)',
    '保育所密度',
    '保育定員密度',
    '大学進学率(%)',
    '消費支出(万円)',
    '住宅地標準価格',
    '有効求人倍率',
    '高齢化率(%)',
    '病院密度',
    '宿泊者数/人',
]

X_raw = np.column_stack([
    in_rate.values,
    out_rate.values,
    nursery_dens.values,
    nursery_cap.values,
    univ_rate.values,
    consumption.values / 10000,   # 万円単位にスケール調整
    land_price.values,
    job_ratio.values,
    aging_rate.values,
    hosp_dens.values,
    tourist_pc.values,
])

y = tfr.values   # 目的変数：合計特殊出生率

print(f"  説明変数 ({len(VAR_NAMES)}個): {VAR_NAMES}")
print(f"  目的変数: 合計特殊出生率（TFR）, 範囲 {y.min():.2f} – {y.max():.2f}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step2. LinearSVRによる変数重み算出
#   線形カーネルSVRの係数の絶対値を「変数重要度」とみなす
# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 60)
print("■ Step2. LinearSVR変数重み算出")
print("=" * 60)

scaler   = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)

svr = LinearSVR(C=1.0, max_iter=10000)
svr.fit(X_scaled, y)

svm_weights = np.abs(svr.coef_)   # shape: (n_features,)
weight_df = pd.DataFrame({
    '変数':          VAR_NAMES,
    'SVM重み（絶対値）': svm_weights,
}).sort_values('SVM重み（絶対値）', ascending=False).reset_index(drop=True)

print("【SVM変数重み（降順）】")
print(weight_df.round(4).to_string(index=False))

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step3. 複合指標の構築
#   SVM重みで加重平均 → 複合指標（地方創生度スコア）を計算
# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 60)
print("■ Step3. 複合指標（SVM重み付き）の構築")
print("=" * 60)

weights_norm   = svm_weights / svm_weights.sum()
composite_index = X_scaled @ weights_norm   # (47,)

rank_df = pd.DataFrame({
    '都道府県':     PREFS,
    '複合指標（SVM）': composite_index,
}).sort_values('複合指標（SVM）', ascending=False).reset_index(drop=True)
rank_df['順位'] = range(1, len(PREFS) + 1)

print("【複合指標ランキング（上位10）】")
print(rank_df.head(10).to_string(index=False))

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step4. PCAとの比較
# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 60)
print("■ Step4. PCA（主成分分析）との比較")
print("=" * 60)

pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)
pca_index = X_pca[:, 0]   # 第1主成分をPCA指標として使用

print(f"  PCA説明分散比: PC1={pca.explained_variance_ratio_[0]:.3f}, "
      f"PC2={pca.explained_variance_ratio_[1]:.3f}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ 図の生成（4枚）
# ═══════════════════════════════════════════════════════════════════════════════

# ────────────────────────────────────────────────────────────────────────────
# 図1：SVM重み棒グラフ（全変数、重み順）
# ────────────────────────────────────────────────────────────────────────────
print("\n図1: SVM重み棒グラフを作成中...")

n_vars = len(VAR_NAMES)
colors_bar = ['#1565C0' if i < 3 else '#90CAF9' for i in range(n_vars)]

fig1, ax1 = plt.subplots(figsize=(10, 6))
# weight_df は降順 → 横棒グラフは上から重み大の順になるよう逆転
bars = ax1.barh(
    weight_df['変数'][::-1],
    weight_df['SVM重み（絶対値）'][::-1],
    color=colors_bar[::-1],
    edgecolor='white',
    alpha=0.88,
)
mean_w = weight_df['SVM重み（絶対値）'].mean()
ax1.axvline(mean_w, color='gray', linestyle='--', linewidth=1.2,
            label=f'平均 = {mean_w:.3f}')
ax1.set_xlabel('LinearSVR係数の絶対値（変数重要度）', fontsize=12)
ax1.set_title(
    'SVM（線形カーネル）による変数重要度\n地方創生関連変数のランキング（2022年度, 47都道府県）',
    fontsize=13, fontweight='bold',
)
ax1.legend(fontsize=10)
ax1.grid(axis='x', alpha=0.3)
for bar, val in zip(bars, weight_df['SVM重み（絶対値）'][::-1]):
    ax1.text(val + mean_w * 0.02, bar.get_y() + bar.get_height() / 2,
             f'{val:.3f}', va='center', fontsize=9, color='#333')

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2024_U5_2_fig1_svm_weights.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  → 2024_U5_2_fig1_svm_weights.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図2：複合指標ランキング（上位・下位 15都道府県）
# ────────────────────────────────────────────────────────────────────────────
print("図2: 複合指標ランキングを作成中...")

top15 = rank_df.head(15).copy()
bot15 = rank_df.tail(15).iloc[::-1].copy()

# 上位15の最高スコア都道府県を強調
top1_name = rank_df.iloc[0]['都道府県']
last_name  = rank_df.iloc[-1]['都道府県']

fig2, axes2 = plt.subplots(1, 2, figsize=(14, 7))
fig2.suptitle('地方創生度 複合指標ランキング（LinearSVR重み付き, 2022年度）',
              fontsize=13, fontweight='bold')

for ax, sub_df, title, clr_hi, clr_lo in zip(
    axes2,
    [top15, bot15],
    ['上位15都道府県', '下位15都道府県'],
    ['#1565C0', '#C62828'],
    ['#90CAF9',  '#EF9A9A'],
):
    bar_colors = [
        clr_hi if row['都道府県'] == top1_name or row['都道府県'] == last_name
        else clr_lo
        for _, row in sub_df.iterrows()
    ]
    ax.barh(sub_df['都道府県'], sub_df['複合指標（SVM）'],
            color=bar_colors, edgecolor='white', alpha=0.88)
    ax.axvline(0, color='black', linewidth=0.8)
    ax.set_xlabel('複合指標スコア（標準化加重和）', fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.grid(axis='x', alpha=0.3)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2024_U5_2_fig2_index_rank.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  → 2024_U5_2_fig2_index_rank.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図3：PCAバイプロット（PC1 vs PC2, ローディング矢印付き）
# ────────────────────────────────────────────────────────────────────────────
print("図3: PCAバイプロットを作成中...")

# 都道府県ラベル表示: 複合指標上位10・下位10
top10_prefs = set(rank_df.head(10)['都道府県'].tolist())
bot10_prefs = set(rank_df.tail(10)['都道府県'].tolist())
label_prefs = top10_prefs | bot10_prefs

fig3, ax3 = plt.subplots(figsize=(11, 8))
sc = ax3.scatter(
    X_pca[:, 0], X_pca[:, 1],
    c=composite_index, cmap='RdYlBu',
    s=65, alpha=0.82, zorder=3,
)
plt.colorbar(sc, ax=ax3, label='SVM複合指標スコア')

for i, pref in enumerate(PREFS):
    if pref in label_prefs:
        ax3.annotate(pref, (X_pca[i, 0], X_pca[i, 1]),
                     textcoords='offset points', xytext=(6, 3),
                     fontsize=7.5, color='#333')

# ローディングベクトル（矢印）
loadings = pca.components_.T   # shape: (n_features, 2)
scale = 2.2
for j, var in enumerate(VAR_NAMES):
    ax3.annotate(
        '', xy=(loadings[j, 0] * scale, loadings[j, 1] * scale), xytext=(0, 0),
        arrowprops=dict(arrowstyle='->', color='#E65100', lw=1.5),
    )
    ax3.text(
        loadings[j, 0] * scale * 1.12,
        loadings[j, 1] * scale * 1.12,
        var, color='#E65100', fontsize=8,
    )

ax3.axhline(0, color='gray', linewidth=0.5)
ax3.axvline(0, color='gray', linewidth=0.5)
ax3.set_xlabel(f'PC1（説明率 {pca.explained_variance_ratio_[0]*100:.1f}%）', fontsize=12)
ax3.set_ylabel(f'PC2（説明率 {pca.explained_variance_ratio_[1]*100:.1f}%）', fontsize=12)
ax3.set_title(
    'PCAバイプロット（都道府県 × 地方創生変数, 2022年度）\n色：SVM複合指標スコア',
    fontsize=13, fontweight='bold',
)
ax3.grid(True, alpha=0.2)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2024_U5_2_fig3_pca.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  → 2024_U5_2_fig3_pca.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図4：SVM複合指標 vs PCA第1主成分 散布図
# ────────────────────────────────────────────────────────────────────────────
print("図4: SVM指標 vs PCA指標散布図を作成中...")

# 両指標を標準化（方向を揃えるため）
svm_std = (composite_index - composite_index.mean()) / composite_index.std()
pca_std = (pca_index - pca_index.mean()) / pca_index.std()

# PCAとSVM指標の方向が逆の場合は符号を揃える（正の相関になるよう）
if scipy_stats.pearsonr(svm_std, pca_std)[0] < 0:
    pca_std = -pca_std

r_val, p_val = scipy_stats.pearsonr(svm_std, pca_std)

# 上位・下位5都道府県をラベル表示
top5_prefs  = set(rank_df.head(5)['都道府県'].tolist())
bot5_prefs  = set(rank_df.tail(5)['都道府県'].tolist())
label5_prefs = top5_prefs | bot5_prefs

scatter_c = [
    '#E53935' if p in top5_prefs else
    '#1E88E5' if p in bot5_prefs else
    '#BDBDBD'
    for p in PREFS
]

fig4, ax4 = plt.subplots(figsize=(9, 7))
ax4.scatter(svm_std, pca_std, c=scatter_c, s=70, alpha=0.85, zorder=3)

for i, pref in enumerate(PREFS):
    if pref in label5_prefs:
        ax4.annotate(pref, (svm_std[i], pca_std[i]),
                     textcoords='offset points', xytext=(7, 3),
                     fontsize=9.5, fontweight='bold',
                     color='#C62828' if pref in top5_prefs else '#0D47A1')

# 回帰直線
x_fit = np.linspace(svm_std.min(), svm_std.max(), 200)
coef  = np.polyfit(svm_std, pca_std, 1)
ax4.plot(x_fit, np.polyval(coef, x_fit),
         color='#444', linewidth=1.8, linestyle='--')

ax4.axhline(0, color='gray', linewidth=0.5)
ax4.axvline(0, color='gray', linewidth=0.5)
ax4.set_xlabel('SVM複合指標（標準化）', fontsize=12)
ax4.set_ylabel('PCA第1主成分スコア（標準化）', fontsize=12)

p_str = '<0.001' if p_val < 0.001 else f'={p_val:.3f}'
ax4.set_title(
    f'SVM複合指標 vs PCA第1主成分（2022年度, 47都道府県）\nr = {r_val:.3f} (p{p_str})',
    fontsize=13, fontweight='bold',
)
ax4.grid(True, alpha=0.2)
ax4.text(
    0.05, 0.95, f'相関係数 r = {r_val:.3f}',
    transform=ax4.transAxes, fontsize=11, va='top',
    bbox=dict(boxstyle='round', facecolor='#E3F2FD', alpha=0.8),
)

# 凡例
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#E53935', label='SVM複合指標 上位5'),
    Patch(facecolor='#1E88E5', label='SVM複合指標 下位5'),
    Patch(facecolor='#BDBDBD', label='その他'),
]
ax4.legend(handles=legend_elements, fontsize=9, loc='lower right')

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2024_U5_2_fig4_scatter.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  → 2024_U5_2_fig4_scatter.png 保存完了")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ まとめ
# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 60)
print("全図の生成完了（4枚）")
print("=" * 60)
print("\n【主要知見】")
print(f"  SVM変数重み上位3変数: {weight_df['変数'].head(3).tolist()}")
print(f"  複合指標1位:  {rank_df.iloc[0]['都道府県']}（スコア = {rank_df.iloc[0]['複合指標（SVM）']:.3f}）")
print(f"  複合指標最下位: {rank_df.iloc[-1]['都道府県']}（スコア = {rank_df.iloc[-1]['複合指標（SVM）']:.3f}）")
print(f"  SVM指標 vs PCA指標 相関: r = {r_val:.3f}")
print(f"  PCA説明率: PC1={pca.explained_variance_ratio_[0]*100:.1f}%, "
      f"PC2={pca.explained_variance_ratio_[1]*100:.1f}%")