"""
教育用再現コード: 2023年度 統計データ分析コンペティション 統計活用奨励賞（高校生）
=================================================================
論文タイトル：独自指標作成による地方創生の方法論と兵庫県活性化の提案

【分析概要】
  データ：SSDSE-B-2026.csv（社会・人口統計体系 2022年度 47都道府県）

  Step1. 複合指標（創生指数）の構築
         6つの構成指標を z スコア標準化 → サブ指数化 → 総合指数化
           経済活力: 有効求人倍率（F3103/F3102）、消費支出（L3221）
           人口定着: 純移動率（(A5101-A5102)/A1101×100）、合計特殊出生率（A4103）
           子育て環境: 保育所密度（J2503/A1101×10000）、定員余裕率（J2505/J2506）
           教育・人材: 大学進学率（E4602/E4601×100）
           医療・福祉: 病院密度（I510120/A1101×10000）
           観光・交流: 宿泊者 per capita（G7101/A1101）

  Step2. 47都道府県の創生指数ランキング（兵庫県ハイライト）

  Step3. PCA バイプロット（PC1 vs PC2、都道府県ラベル付き、兵庫県強調）

  Step4. 6構成指標の相関ヒートマップ

  Step5. 兵庫県 vs 全国平均 vs Top5 のレーダーチャート

  Key findings（実データから得られた傾向）:
    - 東京・神奈川・愛知など大都市圏が経済活力・観光で高得点
    - 兵庫県は全国中位〜上位圏であり、人口定着と子育て環境に改善余地
    - PC1 は「都市型発展」軸、PC2 は「自然増・定着」軸として解釈可能
    - 経済活力と観光・交流の間に強い正の相関が見られる

【データサイエンス学習ポイント】
  1. 複合指標（コンポジット・インデックス）の構築手順と標準化の意義
  2. 主成分分析（PCA）によるバイプロットと次元削減
  3. レーダーチャートによる多次元比較の可視化
  4. 実公的統計データ（SSDSE-B）の加工・活用方法

【データ】SSDSE-B-2026.csv（実公的統計データ）
         合成データ・np.random は一切使用しない
=================================================================
"""

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


import os
import warnings

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

warnings.filterwarnings('ignore')

# ──────────────────────────────────────────────────────────────
# 共通設定
# ──────────────────────────────────────────────────────────────
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

DATA_DIR = 'data/raw'
FIG_DIR  = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

HYOGO_PREF = '兵庫県'
COLOR_HYOGO = '#C62828'    # 兵庫県強調色（赤）
COLOR_NAT   = '#1565C0'    # 全国平均色（青）
COLOR_TOP5  = '#2E7D32'    # Top5色（緑）
COLOR_DEF   = '#90CAF9'    # デフォルト棒グラフ色

# ================================================================
# ■ Step 0. データ読み込み・前処理
# ================================================================

print("=" * 65)
print("■ Step 0. データ読み込み・前処理")
print("=" * 65)

# SSDSE-B は row0=コード行、row1=ラベル行 → header=0 でコードを列名にする
df_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932',
    header=0,
)
# 列名は変数コード（例: A1101）、値の行0 がラベル行
# 実データは row 1 以降
df_data = df_raw.iloc[1:].copy()
df_data = df_data.rename(columns={
    'SSDSE-B-2026': 'YEAR',
    'Code':         'CODE',
    'Prefecture':   'PREF',
})

# 都道府県レベルの行のみ（CODE が R\d{5} パターン）
mask_pref = df_data['CODE'].astype(str).str.match(r'^R\d{5}$')
df_pref = df_data[mask_pref].copy()

# 2022年度に絞る
df = df_pref[df_pref['YEAR'] == '2022'].copy().reset_index(drop=True)
print(f"2022年度 都道府県数: {len(df)}")

# 使用変数コード一覧
USE_COLS = [
    'A1101', 'A5101', 'A5102', 'A4103',
    'F3102', 'F3103', 'L3221',
    'J2503', 'J2505', 'J2506',
    'E4601', 'E4602',
    'I510120',
    'G7101',
]
for col in USE_COLS:
    df[col] = pd.to_numeric(df[col], errors='coerce')

# ────────────────────────────────────────────────
# 構成指標の計算
# ────────────────────────────────────────────────

# 経済活力
df['ind_求人倍率'] = df['F3103'] / df['F3102']            # 有効求人倍率
df['ind_消費支出'] = df['L3221']                           # 消費支出（円/月）

# 人口定着
df['ind_純移動率'] = (df['A5101'] - df['A5102']) / df['A1101'] * 100  # 純移動率（%）
df['ind_TFR']     = df['A4103']                            # 合計特殊出生率

# 子育て環境
df['ind_保育密度']   = df['J2503'] / df['A1101'] * 10000  # 保育所数（人口万対）
df['ind_定員余裕率'] = df['J2505'] / df['J2506']           # 定員 / 在所（余裕率）

# 教育・人材
df['ind_大学進学率'] = df['E4602'] / df['E4601'] * 100    # 高校→大学進学率（%）

# 医療・福祉
df['ind_病院密度'] = df['I510120'] / df['A1101'] * 10000  # 病院数（人口万対）

# 観光・交流
df['ind_宿泊者pc'] = df['G7101'] / df['A1101']            # 宿泊者数 per capita

# ────────────────────────────────────────────────
# z スコア標準化
# ────────────────────────────────────────────────

RAW_INDS = [
    'ind_求人倍率', 'ind_消費支出',
    'ind_純移動率', 'ind_TFR',
    'ind_保育密度', 'ind_定員余裕率',
    'ind_大学進学率',
    'ind_病院密度',
    'ind_宿泊者pc',
]

for col in RAW_INDS:
    mean = df[col].mean()
    std  = df[col].std(ddof=1)
    df[f'z_{col}'] = (df[col] - mean) / std

print("各指標の記述統計（元スケール）:")
print(df[RAW_INDS].describe().round(3).to_string())

# ────────────────────────────────────────────────
# サブ指数（6分野）→ 総合指数
# ────────────────────────────────────────────────

df['sub_経済活力'] = df[['z_ind_求人倍率', 'z_ind_消費支出']].mean(axis=1)
df['sub_人口定着'] = df[['z_ind_純移動率', 'z_ind_TFR']].mean(axis=1)
df['sub_子育て環境'] = df[['z_ind_保育密度', 'z_ind_定員余裕率']].mean(axis=1)
df['sub_教育人材'] = df['z_ind_大学進学率']
df['sub_医療福祉'] = df['z_ind_病院密度']
df['sub_観光交流'] = df['z_ind_宿泊者pc']

SUB_INDICES = [
    'sub_経済活力', 'sub_人口定着', 'sub_子育て環境',
    'sub_教育人材', 'sub_医療福祉', 'sub_観光交流',
]
SUB_LABELS = {
    'sub_経済活力':   '経済活力',
    'sub_人口定着':   '人口定着',
    'sub_子育て環境': '子育て環境',
    'sub_教育人材':   '教育・人材',
    'sub_医療福祉':   '医療・福祉',
    'sub_観光交流':   '観光・交流',
}

# 総合創生指数（6サブ指数の単純平均）
df['創生指数'] = df[SUB_INDICES].mean(axis=1)

# 兵庫県の確認
hyogo_row = df[df['PREF'].str.contains(HYOGO_PREF)].iloc[0]
hyogo_rank = int(df['創生指数'].rank(ascending=False)[hyogo_row.name])
print(f"\n兵庫県 創生指数: {hyogo_row['創生指数']:.4f}  全国順位: {hyogo_rank}位")

# PCA 用データ（9変数 z スコア）
Z_COLS = [f'z_{c}' for c in RAW_INDS]
Z_LABELS = [
    '求人倍率', '消費支出',
    '純移動率', 'TFR',
    '保育密度', '定員余裕率',
    '大学進学率',
    '病院密度',
    '宿泊者pc',
]

df_z = df[Z_COLS].copy()

# ================================================================
# ■ 図1: 47都道府県の創生指数ランキング（横棒グラフ）
# ================================================================

print("\n" + "=" * 65)
print("■ 図1: 47都道府県の創生指数ランキング")
print("=" * 65)

df_rank = df[['PREF', '創生指数']].copy()
df_rank = df_rank.sort_values('創生指数', ascending=True).reset_index(drop=True)

fig1, ax1 = plt.subplots(figsize=(10, 14))

bar_colors = [
    COLOR_HYOGO if '兵庫' in str(row['PREF']) else COLOR_DEF
    for _, row in df_rank.iterrows()
]
bars = ax1.barh(
    df_rank['PREF'],
    df_rank['創生指数'],
    color=bar_colors,
    edgecolor='white',
    linewidth=0.5,
)

# 兵庫県のバーに値ラベル
for bar, (_, row) in zip(bars, df_rank.iterrows()):
    if '兵庫' in str(row['PREF']):
        ax1.text(
            row['創生指数'] + 0.02, bar.get_y() + bar.get_height() / 2,
            f"  {row['創生指数']:.3f}（全国{hyogo_rank}位）",
            va='center', fontsize=9, color=COLOR_HYOGO, fontweight='bold',
        )

ax1.axvline(0, color='black', linewidth=0.8, linestyle='--', alpha=0.5)
ax1.set_xlabel('創生指数（標準化済み、6分野平均）', fontsize=11)
ax1.set_title(
    '地方創生複合指標：47都道府県の創生指数ランキング（2022年度）\n'
    '（出典：SSDSE-B-2026 実データ）',
    fontsize=12, fontweight='bold', pad=12,
)
ax1.grid(axis='x', alpha=0.25)
ax1.tick_params(axis='y', labelsize=9)

legend_patches = [
    mpatches.Patch(color=COLOR_HYOGO, label='兵庫県'),
    mpatches.Patch(color=COLOR_DEF,   label='その他都道府県'),
]
ax1.legend(handles=legend_patches, fontsize=9, loc='lower right')

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2023_H4_fig1_ranking.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  -> 2023_H4_fig1_ranking.png 保存完了")

# ================================================================
# ■ 図2: レーダーチャート（兵庫県 vs 全国平均 vs Top5）
# ================================================================

print("\n■ 図2: レーダーチャート（兵庫県 vs 全国平均 vs Top5）")

# Top5 平均（創生指数上位5都道府県）
top5_prefs = df.nlargest(5, '創生指数')['PREF'].tolist()
print(f"  Top5: {top5_prefs}")

hyogo_vals = [hyogo_row[s] for s in SUB_INDICES]
nat_vals   = [df[s].mean() for s in SUB_INDICES]
top5_vals  = [df[df['PREF'].isin(top5_prefs)][s].mean() for s in SUB_INDICES]

categories = [SUB_LABELS[s] for s in SUB_INDICES]
N = len(categories)

import matplotlib.pyplot as plt  # already imported but re-reference for clarity

angles = np.linspace(0, 2 * np.pi, N, endpoint=False).tolist()
angles += angles[:1]

def close_loop(vals):
    return vals + vals[:1]

hyogo_v = close_loop(hyogo_vals)
nat_v   = close_loop(nat_vals)
top5_v  = close_loop(top5_vals)

fig2, ax2 = plt.subplots(figsize=(8, 8), subplot_kw=dict(polar=True))

ax2.plot(angles, hyogo_v, color=COLOR_HYOGO, linewidth=2.5, label='兵庫県')
ax2.fill(angles, hyogo_v, color=COLOR_HYOGO, alpha=0.18)

ax2.plot(angles, nat_v, color=COLOR_NAT, linewidth=2.0, linestyle='--', label='全国平均')
ax2.fill(angles, nat_v, color=COLOR_NAT, alpha=0.10)

top5_label = 'Top5平均\n（' + '、'.join([p[:3] for p in top5_prefs[:3]]) + '…）'
ax2.plot(angles, top5_v, color=COLOR_TOP5, linewidth=2.0, linestyle=':', label=top5_label)
ax2.fill(angles, top5_v, color=COLOR_TOP5, alpha=0.10)

ax2.set_xticks(angles[:-1])
ax2.set_xticklabels(categories, fontsize=11, fontweight='bold')
ax2.set_ylim(
    min(min(hyogo_vals), min(nat_vals), min(top5_vals)) - 0.3,
    max(max(hyogo_vals), max(nat_vals), max(top5_vals)) + 0.3,
)
ax2.set_title(
    '地方創生サブ指数：兵庫県 vs 全国平均 vs Top5（2022年度）\n'
    '（出典：SSDSE-B-2026 実データ）',
    fontsize=11, fontweight='bold', pad=22,
)
ax2.legend(loc='upper right', bbox_to_anchor=(1.32, 1.15), fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2023_H4_fig2_radar.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  -> 2023_H4_fig2_radar.png 保存完了")

# ================================================================
# ■ 図3: PCA バイプロット（PC1 vs PC2）
# ================================================================

print("\n■ 図3: PCA バイプロット")

X_pca = df_z.values
pca   = PCA(n_components=2)
scores = pca.fit_transform(X_pca)

pc1_var = pca.explained_variance_ratio_[0] * 100
pc2_var = pca.explained_variance_ratio_[1] * 100
loadings = pca.components_  # shape (2, n_features)

print(f"  PC1 寄与率: {pc1_var:.1f}%  PC2 寄与率: {pc2_var:.1f}%")
print(f"  累積寄与率: {pc1_var + pc2_var:.1f}%")

fig3, ax3 = plt.subplots(figsize=(11, 9))

# 都道府県スコア散布図
for i, (x, y) in enumerate(scores):
    pref = df.iloc[i]['PREF']
    is_hyogo = '兵庫' in str(pref)
    color  = COLOR_HYOGO if is_hyogo else '#90CAF9'
    zorder = 5 if is_hyogo else 2
    size   = 100 if is_hyogo else 55
    ax3.scatter(x, y, color=color, s=size, alpha=0.85,
                edgecolors='white', linewidth=0.7, zorder=zorder)
    fontsize = 9.5 if is_hyogo else 7
    fw = 'bold' if is_hyogo else 'normal'
    fc = COLOR_HYOGO if is_hyogo else '#333333'
    ax3.annotate(pref, (x, y), fontsize=fontsize, color=fc,
                 fontweight=fw, xytext=(4, 4), textcoords='offset points',
                 zorder=zorder)

# ローディングベクトル
scale = 3.0   # ベクトルを見やすい倍率にスケール
for j, label in enumerate(Z_LABELS):
    lx = loadings[0, j] * scale
    ly = loadings[1, j] * scale
    ax3.annotate(
        '', xy=(lx, ly), xytext=(0, 0),
        arrowprops=dict(arrowstyle='->', color='#E65100', lw=1.8),
    )
    ax3.text(lx * 1.08, ly * 1.08, label,
             fontsize=8.5, color='#E65100', fontweight='bold',
             ha='center', va='center')

ax3.axhline(0, color='grey', linewidth=0.7, linestyle='--', alpha=0.5)
ax3.axvline(0, color='grey', linewidth=0.7, linestyle='--', alpha=0.5)
ax3.set_xlabel(f'PC1（寄与率 {pc1_var:.1f}%）', fontsize=12)
ax3.set_ylabel(f'PC2（寄与率 {pc2_var:.1f}%）', fontsize=12)
ax3.set_title(
    '主成分分析バイプロット：地方創生指標（2022年度 47都道府県）\n'
    '（出典：SSDSE-B-2026 実データ）',
    fontsize=11, fontweight='bold', pad=12,
)
ax3.grid(alpha=0.2)

legend_patches = [
    mpatches.Patch(color=COLOR_HYOGO, label='兵庫県'),
    mpatches.Patch(color='#90CAF9',   label='その他都道府県'),
    mpatches.FancyArrow(0, 0, 1, 0, color='#E65100', width=0.01,
                        label='変数ローディング'),
]
ax3.legend(handles=legend_patches, fontsize=9, loc='lower right')

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2023_H4_fig3_pca_biplot.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  -> 2023_H4_fig3_pca_biplot.png 保存完了")

# ================================================================
# ■ 図4: 6構成指標の相関ヒートマップ
# ================================================================

print("\n■ 図4: 6構成指標の相関ヒートマップ")

# サブ指数の相関行列
corr_data = df[SUB_INDICES].copy()
corr_mat  = corr_data.corr()
labels_disp = [SUB_LABELS[s] for s in SUB_INDICES]

fig4, ax4 = plt.subplots(figsize=(8, 7))
im4 = ax4.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')

ax4.set_xticks(range(len(SUB_INDICES)))
ax4.set_xticklabels(labels_disp, fontsize=10, rotation=30, ha='right')
ax4.set_yticks(range(len(SUB_INDICES)))
ax4.set_yticklabels(labels_disp, fontsize=10)

for i in range(len(SUB_INDICES)):
    for j in range(len(SUB_INDICES)):
        val = corr_mat.values[i, j]
        txt_color = 'white' if abs(val) > 0.55 else 'black'
        ax4.text(j, i, f'{val:.2f}', ha='center', va='center',
                 fontsize=9.5, color=txt_color, fontweight='bold')

plt.colorbar(im4, ax=ax4, fraction=0.046, pad=0.04,
             label='ピアソン相関係数')
ax4.set_title(
    '地方創生6サブ指数の相関行列（2022年度 47都道府県）\n'
    '（出典：SSDSE-B-2026 実データ）',
    fontsize=11, fontweight='bold', pad=12,
)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2023_H4_fig4_corr_heatmap.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  -> 2023_H4_fig4_corr_heatmap.png 保存完了")

# ================================================================
# ■ 完了サマリ
# ================================================================

print("\n" + "=" * 65)
print("完了: 全図の生成完了（4枚）")
print("=" * 65)
print(f"\n保存先: {os.path.abspath(FIG_DIR)}")
print("  2023_H4_fig1_ranking.png        - 47都道府県の創生指数ランキング")
print("  2023_H4_fig2_radar.png          - 兵庫県 vs 全国平均 vs Top5 レーダー")
print("  2023_H4_fig3_pca_biplot.png     - PCA バイプロット（PC1 vs PC2）")
print("  2023_H4_fig4_corr_heatmap.png   - 6サブ指数の相関ヒートマップ")
print()
print("使用データ: SSDSE-B-2026.csv（社会・人口統計体系 2022年度 47都道府県）")
print("合成データ: なし（np.random 未使用）")
print()

# 兵庫県詳細レポート
print("=" * 65)
print("■ 兵庫県 創生指数詳細")
print("=" * 65)
print(f"  総合創生指数:  {hyogo_row['創生指数']:.4f}  （全国{hyogo_rank}位）")
for s in SUB_INDICES:
    nat_avg = df[s].mean()
    pref_rank = int(df[s].rank(ascending=False)[hyogo_row.name])
    print(f"  {SUB_LABELS[s]:<10}: {hyogo_row[s]:+.4f}  全国平均 {nat_avg:+.4f}  （{pref_rank}位）")