"""
教育用再現コード: 2022年 統計データ分析コンペティション 優秀賞 [高校生の部]
=================================================================
論文タイトル：ボランティア活動の決定要因：地域特性・社会資本の相関分析
受賞：優秀賞（高校生の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別データ, 2022年度）

  ※ ボランティア活動率の直接データはSSDSE-Bに含まれないため、
    社会資本・地域コミュニティの代理指標として以下を使用：
    - 保育所等数（地域コミュニティ密度の代理）
    - 婚姻件数（社会関係の活発さの代理）
    - 幼稚園・保育所等利用率（地域参加意識の代理）

  分析の流れ
  1. 相関ヒートマップ：社会資本関連指標の相関構造
  2. 散布図：高齢化率 vs 保育所密度（地域コミュニティ構造）
  3. OLS回帰：婚姻率の決定要因（社会的絆の指標）
  4. 地域別比較：6地域ブロックの社会資本指標箱ひげ図

【代理変数設計】
  社会的絆 → 婚姻率 = 婚姻件数 / 総人口 × 1000
  コミュニティ密度 → 保育所密度 = 保育所等数 / 総人口 × 10000
  高齢化率 = 65歳以上人口 / 総人口 × 100
  教育費（二人以上世帯）
  消費支出_log = log(消費支出（二人以上世帯）)
  大学進学率 = 高校卒業者のうち進学者数 / 高校卒業者数 × 100
  離婚率 = 離婚件数 / 総人口 × 1000

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
=================================================================
"""

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


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.api as sm
import warnings
warnings.filterwarnings('ignore')

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

import os
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)
df_2022 = df_b[df_b['年度'] == 2022].copy()

# 変数生成
df_2022['婚姻率'] = df_2022['婚姻件数'] / df_2022['総人口'] * 1000
df_2022['保育所密度'] = df_2022['保育所等数'] / df_2022['総人口'] * 10000
df_2022['高齢化率'] = df_2022['65歳以上人口'] / df_2022['総人口'] * 100
df_2022['大学進学率'] = df_2022['高等学校卒業者のうち進学者数'] / df_2022['高等学校卒業者数'].replace(0, np.nan) * 100
df_2022['消費支出_log'] = np.log(df_2022['消費支出（二人以上の世帯）'].clip(lower=1))
df_2022['教育費_万円'] = df_2022['教育費（二人以上の世帯）'] / 10000
df_2022['離婚率'] = df_2022['離婚件数'] / df_2022['総人口'] * 1000

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

# Fig1: 相関ヒートマップ
analysis_vars = ['婚姻率', '保育所密度', '高齢化率', '大学進学率', '消費支出_log', '教育費_万円', '離婚率']
df_corr = df_2022[analysis_vars].dropna()
corr = df_corr.corr()
fig, ax = plt.subplots(figsize=(9, 7))
im = ax.imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)
ax.set_xticks(range(len(analysis_vars)))
ax.set_yticks(range(len(analysis_vars)))
ax.set_xticklabels(analysis_vars, rotation=45, ha='right', fontsize=9)
ax.set_yticklabels(analysis_vars, fontsize=9)
for i in range(len(analysis_vars)):
    for j in range(len(analysis_vars)):
        val = corr.iloc[i, j]
        ax.text(j, i, f'{val:.2f}', ha='center', va='center',
                fontsize=7.5, color='white' if abs(val) > 0.5 else 'black')
plt.colorbar(im, ax=ax, label='Pearson相関係数')
ax.set_title('社会資本指標間の相関ヒートマップ（2022年, N=47）', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H2_fig1_corr.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

# Fig2: 散布図 高齢化率 vs 保育所密度（地域別色分け）
df_sc = df_2022[['都道府県', '高齢化率', '保育所密度', '地域']].dropna()
fig, ax = plt.subplots(figsize=(9, 6))
for reg, grp in df_sc.groupby('地域'):
    ax.scatter(grp['高齢化率'], grp['保育所密度'],
               color=region_colors.get(reg, 'gray'), label=reg, s=65, alpha=0.8)
for _, row in df_sc.iterrows():
    ax.annotate(row['都道府県'][:2], (row['高齢化率'], row['保育所密度']),
                fontsize=7, alpha=0.7, xytext=(3, 3), textcoords='offset points')
slope, intercept, r, p, _ = stats.linregress(df_sc['高齢化率'], df_sc['保育所密度'])
xr = np.linspace(df_sc['高齢化率'].min(), df_sc['高齢化率'].max(), 100)
ax.plot(xr, slope * xr + intercept, 'k--', linewidth=1.5,
        label=f'回帰直線 (r={r:.2f}, p={p:.3f})')
ax.set_xlabel('高齢化率（%）', fontsize=12)
ax.set_ylabel('保育所密度（/人口万対）', fontsize=12)
ax.set_title('高齢化率 vs 保育所密度（2022年, N=47）\n（地域コミュニティ構造の可視化）', fontsize=12, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H2_fig2_scatter.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

# Fig3: OLS回帰 婚姻率の決定要因
xvars = ['高齢化率', '大学進学率', '消費支出_log', '保育所密度']
df_reg = df_2022[['婚姻率'] + xvars].dropna()
X = sm.add_constant(df_reg[xvars])
res = sm.OLS(df_reg['婚姻率'], X).fit()
print(res.summary())
coefs = res.params.drop('const')
ses = res.bse.drop('const')
pvals = res.pvalues.drop('const')
fig, ax = plt.subplots(figsize=(8, 5))
colors_c = ['#e05c5c' if p < 0.05 else '#888888' for p in pvals]
ax.barh(range(len(coefs)), coefs, xerr=1.96 * ses, color=colors_c, alpha=0.8,
        error_kw={'elinewidth': 1.5, 'capsize': 4})
ax.set_yticks(range(len(coefs)))
ax.set_yticklabels(coefs.index, fontsize=10)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('OLS回帰係数', fontsize=12)
ax.set_title(f'婚姻率の決定要因（社会的絆の代理, N=47）\nR²={res.rsquared:.3f}（赤=p<0.05）', fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H2_fig3_ols.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved")

# Fig4: 地域別箱ひげ図（保育所密度）
df_box = df_2022[['地域', '保育所密度', '婚姻率']].dropna()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
bxdata1 = [df_box[df_box['地域'] == reg]['保育所密度'].values for reg in region_order]
bp1 = ax1.boxplot(bxdata1, patch_artist=True,
                   labels=[r[:4] for r in region_order])
for patch, reg in zip(bp1['boxes'], region_order):
    patch.set_facecolor(region_colors.get(reg, 'gray'))
    patch.set_alpha(0.7)
ax1.set_xlabel('地域', fontsize=11)
ax1.set_ylabel('保育所密度（/万人）', fontsize=11)
ax1.set_title('地域別 保育所密度の分布', fontsize=12, fontweight='bold')
ax1.grid(axis='y', alpha=0.3)
bxdata2 = [df_box[df_box['地域'] == reg]['婚姻率'].values for reg in region_order]
bp2 = ax2.boxplot(bxdata2, patch_artist=True,
                   labels=[r[:4] for r in region_order])
for patch, reg in zip(bp2['boxes'], region_order):
    patch.set_facecolor(region_colors.get(reg, 'gray'))
    patch.set_alpha(0.7)
ax2.set_xlabel('地域', fontsize=11)
ax2.set_ylabel('婚姻率（件/千人）', fontsize=11)
ax2.set_title('地域別 婚姻率の分布', fontsize=12, fontweight='bold')
ax2.grid(axis='y', alpha=0.3)
plt.suptitle('地域ブロック別 社会資本指標の分布（2022年）', fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H2_fig4_boxplot.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")