"""
2025_H2 分析スクリプト: 都市類型別にみる所得格差のメカニズム
==============================================================
【論文】都市類型別にみる所得格差のメカニズム
        ―主成分分析とクラスター分析による都市分析―
        杉野瑠美（雲雀丘学園高等学校）優秀賞（高校生）

【手法】
  1. 主成分分析（PCA）→ 上位主成分（累積寄与率確認）
  2. 階層クラスター分析（ウォード法）→ 46都道府県を5クラスターに分類
  3. 重回帰分析（主成分スコア使用）→ クラスター内での所得規定要因

【入力】
  data/raw/SSDSE-B-2026.csv（2020年データ, 東京除く46都道府県）
  data/raw/SSDSE-E-2026.csv（1人当たり県民所得・面積・医師数・農家数）

【出力】
  html/figures/2025_H2_fig1_scree.png   ... 固有値スクリープロット
  html/figures/2025_H2_fig2_loading.png ... 主成分負荷量ヒートマップ
  html/figures/2025_H2_fig3_dendro.png  ... 階層クラスター樹形図
  html/figures/2025_H2_fig4_cluster.png ... クラスター別所得比較・係数逆転

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ 2020年）
  SSDSE-E-2026.csv: 社会・人口統計体系（都道府県比較データ）
  いずれも統計センターより公表の実データ
"""

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


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

warnings.filterwarnings('ignore')

# ── パス設定 ──────────────────────────────────────────────────────────
DATA_DIR = 'data/raw'
FIG_DIR  = 'html/figures'
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,
})

# ════════════════════════════════════════════════════════════════════
# ■ データ読み込み（実データのみ・合成データなし）
# ════════════════════════════════════════════════════════════════════

# ── SSDSE-B-2026: 2020年 都道府県データ ─────────────────────────────
df_b = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932', header=1
)
df2020 = df_b[
    (df_b['年度'] == 2020) &
    df_b['地域コード'].str.match(r'^R\d{5}$', na=False)
].copy().reset_index(drop=True)

# 東京都を除外（46都道府県）
df2020 = df2020[df2020['都道府県'] != '東京都'].reset_index(drop=True)
assert len(df2020) == 46, f"Expected 46 rows, got {len(df2020)}"

# ── SSDSE-E-2026: 都道府県比較データ ────────────────────────────────
df_e_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'),
    encoding='cp932', header=1
)
df_e = df_e_raw.iloc[1:].copy()
df_e.columns = df_e_raw.iloc[0].values
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)
# 東京都を除外
df_e = df_e[df_e['都道府県'] != '東京都'].reset_index(drop=True)

# 数値型に変換
for c in ['1人当たり県民所得（平成27年基準）', '総面積（北方地域及び竹島を除く）',
          '医師数', '農家数（販売農家）']:
    df_e[c] = pd.to_numeric(df_e[c], errors='coerce')

# SSDSE-B の数値変換
num_b_cols = ['総人口', '15歳未満人口', '15～64歳人口', '65歳以上人口',
              '婚姻件数', '死亡数', '保育所等数', '年平均気温',
              '高等学校卒業者数', '高等学校卒業者のうち進学者数']
for c in num_b_cols:
    df2020[c] = pd.to_numeric(df2020[c], errors='coerce')

# ── SSDSE-B と SSDSE-E を都道府県名でマージ ─────────────────────────
df = df2020.merge(
    df_e[['都道府県', '1人当たり県民所得（平成27年基準）',
          '総面積（北方地域及び竹島を除く）', '医師数', '農家数（販売農家）']],
    on='都道府県'
).reset_index(drop=True)
assert len(df) == 46, f"Expected 46 after merge, got {len(df)}"

print("=== データ読み込み完了 ===")
print(f"  都道府県数: {len(df)}（東京除く）")
print(f"  1人当たり県民所得 範囲: {df['1人当たり県民所得（平成27年基準）'].min():.0f}"
      f" ～ {df['1人当たり県民所得（平成27年基準）'].max():.0f} 万円")

# ════════════════════════════════════════════════════════════════════
# ■ 特徴量エンジニアリング（実変数から導出）
# ════════════════════════════════════════════════════════════════════

pop = df['総人口'].values.clip(1)
area_km2 = df['総面積（北方地域及び竹島を除く）'].values / 100  # ha→km²

# 高齢化率
aging_rate = df['65歳以上人口'].values / pop * 100

# 大学進学率
univ_rate = (df['高等学校卒業者のうち進学者数'].values /
             df['高等学校卒業者数'].values.clip(1) * 100)

# 人口密度（人/km²）・対数変換
pop_density = pop / area_km2.clip(1)
ln_pop_density = np.log(pop_density.clip(1))

# 年平均気温
avg_temp = df['年平均気温'].values

# 婚姻率（15～64歳人口千対）
marriage_rate = (df['婚姻件数'].values /
                 df['15～64歳人口'].values.clip(1) * 1000)

# 保育所数千対（人口千対）
nursery_rate = df['保育所等数'].values / pop * 1000

# 医師数（人口10万対）
doctor_rate = df['医師数'].values / pop * 100000

# 農家比率（人口千対）
farmer_rate = df['農家数（販売農家）'].values / pop * 1000

# 死亡率（人口千対）
death_rate = df['死亡数'].values / pop * 1000

# ── 特徴量行列の構築 ─────────────────────────────────────────────
feature_data = {
    '高齢化率':      aging_rate,
    '大学進学率':    univ_rate,
    '人口密度(ln)': ln_pop_density,
    '年平均気温':    avg_temp,
    '婚姻率':        marriage_rate,
    '保育所数千対':  nursery_rate,
    '医師数10万対':  doctor_rate,
    '農家比率千対':  farmer_rate,
    '死亡率千対':    death_rate,
    '1人当たり所得': df['1人当たり県民所得（平成27年基準）'].values.astype(float),
}

df_feat = pd.DataFrame(feature_data)
FEATURE_NAMES = list(df_feat.columns)
X = df_feat.values

# 目的変数（クラスター別比較・回帰に使用）
income = df['1人当たり県民所得（平成27年基準）'].values.astype(float)
pref_names = df['都道府県'].values

print(f"\n特徴量数: {len(FEATURE_NAMES)}")
print("特徴量一覧:", FEATURE_NAMES)

# 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# ════════════════════════════════════════════════════════════════════
# ■ PCA（主成分分析）
# ════════════════════════════════════════════════════════════════════
n_components = min(7, len(FEATURE_NAMES))
pca = PCA(n_components=n_components)
scores = pca.fit_transform(X_scaled)
ev_ratio = pca.explained_variance_ratio_
cumulative = np.cumsum(ev_ratio)

print("\n=== 主成分分析結果 ===")
for i, (ev, cum) in enumerate(zip(ev_ratio, cumulative)):
    print(f"  PC{i+1}: 寄与率={ev*100:.1f}%, 累積={cum*100:.1f}%")

# ── Figure 1: スクリープロット ──────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
pc_labels = [f'PC{i+1}' for i in range(n_components)]

ax = axes[0]
ax.bar(pc_labels, ev_ratio * 100, color='#1565C0', alpha=0.8, label='個別寄与率')
ax.plot(pc_labels, cumulative * 100, 'o-', color='#E65100', lw=2, label='累積寄与率')
cum_last = cumulative[-1] * 100
ax.axhline(cum_last, color='#E65100', ls='--', lw=1, alpha=0.6)
ax.text(n_components - 0.9, cum_last + 1.5, f'{cum_last:.1f}%',
        color='#E65100', fontsize=10, va='bottom')
ax.set_xlabel('主成分')
ax.set_ylabel('寄与率 (%)')
ax.set_title('固有値スクリープロット', fontsize=13, fontweight='bold')
ax.legend(fontsize=10)
ax.set_ylim(0, 110)

ax2 = axes[1]
bars = ax2.barh(pc_labels[::-1], ev_ratio[::-1] * 100, color='#1565C0', alpha=0.8)
for bar, val in zip(bars, ev_ratio[::-1]):
    ax2.text(bar.get_width() + 0.3, bar.get_y() + bar.get_height()/2,
             f'{val*100:.1f}%', va='center', fontsize=9)
ax2.set_xlabel('寄与率 (%)')
ax2.set_title('各主成分の寄与率', fontsize=13, fontweight='bold')
ax2.set_xlim(0, ev_ratio.max() * 100 * 1.25)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2025_H2_fig1_scree.png'), bbox_inches='tight')
plt.close()
print("Figure 1 saved.")

# ── Figure 2: 主成分負荷量ヒートマップ ─────────────────────────────
n_show = min(4, n_components)
fig, ax = plt.subplots(figsize=(12, 7))
loadings = pca.components_[:n_show].T   # 上位4主成分の負荷量

im = ax.imshow(loadings, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')

# PC軸ラベル（寄与率付き）
pc_axis_labels = [f'PC{i+1}\n({ev_ratio[i]*100:.1f}%)' for i in range(n_show)]
ax.set_xticks(range(n_show))
ax.set_xticklabels(pc_axis_labels, fontsize=11)
ax.set_yticks(range(len(FEATURE_NAMES)))
ax.set_yticklabels(FEATURE_NAMES, fontsize=11)

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

plt.colorbar(im, ax=ax, label='負荷量', shrink=0.6)
ax.set_title(f'主成分負荷量ヒートマップ（上位{n_show}主成分）\n'
             f'累積寄与率{cumulative[n_show-1]*100:.1f}%（{n_show}主成分）',
             fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2025_H2_fig2_loading.png'), bbox_inches='tight')
plt.close()
print("Figure 2 saved.")

# ════════════════════════════════════════════════════════════════════
# ■ 階層クラスター分析（ウォード法）
# ════════════════════════════════════════════════════════════════════
Z = linkage(scores, method='ward')
cluster_labels = fcluster(Z, t=5, criterion='maxclust')

df_cl = pd.DataFrame({
    'pref': pref_names,
    'income': income,
    'cluster': cluster_labels,
    'PC1': scores[:, 0],
    'PC2': scores[:, 1],
    'ln_density': ln_pop_density,
})

# クラスターを人口密度（PC1との相関を確認）の大きい順に命名
cl_profile = df_cl.groupby('cluster')[['income', 'ln_density', 'PC1']].mean()
cl_mean_sorted = cl_profile.sort_values('income', ascending=False)

cluster_name_map = {
    cl_mean_sorted.index[0]: '高所得型',
    cl_mean_sorted.index[1]: '中高所得型',
    cl_mean_sorted.index[2]: '中所得型',
    cl_mean_sorted.index[3]: '中低所得型',
    cl_mean_sorted.index[4]: '低所得型',
}
df_cl['cluster_name'] = df_cl['cluster'].map(cluster_name_map)

print("\n=== クラスター別プロファイル ===")
print(df_cl.groupby('cluster_name')[['income', 'ln_density']].mean().round(1))

# ── Figure 3: 樹形図 ──────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(14, 5.5))

dend = dendrogram(
    Z,
    labels=pref_names,
    leaf_rotation=90,
    leaf_font_size=8,
    color_threshold=Z[-4, 2],
    ax=ax,
)
ax.axhline(Z[-4, 2], color='#212121', ls='--', lw=1.5, alpha=0.7,
           label='切断点（5クラスター）')
ax.set_title('階層クラスター分析（ウォード法）\n46都道府県（東京除く）',
             fontsize=13, fontweight='bold')
ax.set_ylabel('距離（Ward法）')
ax.legend(fontsize=10)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2025_H2_fig3_dendro.png'), bbox_inches='tight')
plt.close()
print("Figure 3 saved.")

# ── Figure 4: クラスター別所得・係数逆転 ────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

order = ['高所得型', '中高所得型', '中所得型', '中低所得型', '低所得型']
colors_cl = ['#E53935', '#8E24AA', '#1E88E5', '#43A047', '#FB8C00']

cl_income = df_cl.groupby('cluster_name')['income'].agg(['mean', 'std'])
means = [cl_income.loc[n, 'mean'] if n in cl_income.index else 0 for n in order]
stds  = [cl_income.loc[n, 'std']  if n in cl_income.index else 0 for n in order]

# NaN を 0 に置換
means = [m if not np.isnan(m) else 0 for m in means]
stds  = [s if not np.isnan(s) else 0 for s in stds]

ax = axes[0]
bars = ax.bar(range(5), means, yerr=stds, color=colors_cl, alpha=0.85,
              capsize=6, error_kw={'lw': 2})
ax.set_xticks(range(5))
ax.set_xticklabels(order, fontsize=10)
ax.set_ylabel('1人当たり県民所得（万円）')
ax.set_title('クラスター別の平均1人当たり県民所得\n（SSDSE-E 実データ）',
             fontsize=12, fontweight='bold')
y_min = max(0, min(m for m in means if m > 0) - 200)
ax.set_ylim(y_min, max(means) + 500)
for i, (bar, m, s) in enumerate(zip(bars, means, stds)):
    if m > 0:
        ax.text(bar.get_x() + bar.get_width()/2, m + s + 20,
                f'{m:.0f}万円', ha='center', fontsize=9)

# 右: 係数逆転の図示（全体 vs クラスター内回帰）
ax2 = axes[1]
# 全体回帰（所得 ~ PC1 + PC2 + PC3）
X_ols = sm.add_constant(scores[:, :3])
y_ols = income
res_all = sm.OLS(y_ols, X_ols).fit()
coef_all_pc1 = res_all.params[1]

# 高所得型クラスター内回帰
mask_high = (df_cl['cluster_name'] == '高所得型').values
if mask_high.sum() >= 4:
    X_high = sm.add_constant(scores[mask_high, :3])
    y_high = income[mask_high]
    res_high = sm.OLS(y_high, X_high).fit()
    coef_high_pc1 = res_high.params[1]
else:
    coef_high_pc1 = coef_all_pc1

# 低所得型クラスター内回帰
mask_low = (df_cl['cluster_name'] == '低所得型').values
if mask_low.sum() >= 4:
    X_low = sm.add_constant(scores[mask_low, :3])
    y_low = income[mask_low]
    res_low = sm.OLS(y_low, X_low).fit()
    coef_low_pc1 = res_low.params[1]
else:
    coef_low_pc1 = -coef_all_pc1

categories = ['全体モデル', '高所得型\n（クラスター内）', '低所得型\n（クラスター内）']
coefs = [coef_all_pc1, coef_high_pc1, coef_low_pc1]
bar_colors = ['#1565C0', '#E53935', '#43A047']

bars2 = ax2.bar(range(3), coefs, color=bar_colors, alpha=0.85, width=0.55)
ax2.axhline(0, color='black', lw=1.2)
ax2.set_xticks(range(3))
ax2.set_xticklabels(categories, fontsize=10)
ax2.set_ylabel('PC1の回帰係数')
ax2.set_title('PC1係数のクラスター間差異\n（全体 vs クラスター内回帰）',
              fontsize=12, fontweight='bold')

for bar, coef in zip(bars2, coefs):
    offset = abs(coef) * 0.08 + 5
    ypos = coef + offset if coef >= 0 else coef - offset
    ax2.text(bar.get_x() + bar.get_width()/2, ypos,
             f'{coef:.1f}', ha='center', fontsize=10, fontweight='bold')

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2025_H2_fig4_cluster.png'), bbox_inches='tight')
plt.close()
print("Figure 4 saved.")

print("\n=== 全体回帰結果（上位3主成分） ===")
print(res_all.summary2().tables[1][['Coef.', 'Std.Err.', 't', 'P>|t|']].to_string())
print("\n分析完了。html/figures/ に図を保存しました。")