"""
教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================
論文タイトル：都市性と医師密度・年齢調整死亡率の関連性の解析
著者：大分工業高等専門学校

【データ】
  SSDSE-B-2026.csv（2022年度 死亡数・人口）
  SSDSE-E-2026.csv（医師数, 総面積）
  ※死亡率は粗死亡率（死亡数/人口×10万）を使用
    （論文では年齢調整死亡率を使用 — 教育目的の近似値）

【分析概要】
  Step1. Pearson相関分析（医師数（10万対）vs 粗死亡率）
  Step2. Mann-Whitney U検定（人口密度上位10県 vs 下位10県）
  Step3. 散布図 + 線形近似
  Step4. 都道府県別ランキング
=================================================================
"""

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


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats
from matplotlib.patches import Patch
import warnings
warnings.filterwarnings('ignore')
import os

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

BASE_DIR = os.path.join(_script_dir, '..')
FIG_DIR  = os.path.join(BASE_DIR, 'html', 'figures')
DATA_DIR = os.path.join(BASE_DIR, 'data', 'raw')

# ================================================================
# ■ データ読み込み
# ================================================================
# SSDSE-B: 死亡数・人口（2022年）
ssdse_b = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'), encoding='cp932', header=1)
df_b = ssdse_b[
    (ssdse_b['年度'] == 2022) &
    ssdse_b['地域コード'].str.match(r'^R\d{5}$', na=False)
].copy().reset_index(drop=True)

# SSDSE-E: 医師数（2022年）, 総面積
ssdse_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'), encoding='cp932', header=1)
ssdse_e = ssdse_e_raw.iloc[1:].copy()
ssdse_e.columns = ssdse_e_raw.iloc[0].values
ssdse_e = ssdse_e[ssdse_e['都道府県'] != '全国'].reset_index(drop=True)

# ================================================================
# ■ 変数の構築
# ================================================================
df = pd.DataFrame()
df['都道府県'] = df_b['都道府県']

pop_total = pd.to_numeric(df_b['総人口'], errors='coerce')
pop_male  = pd.to_numeric(df_b['総人口（男）'], errors='coerce')
pop_female = pd.to_numeric(df_b['総人口（女）'], errors='coerce')
death_total = pd.to_numeric(df_b['死亡数'], errors='coerce')
death_male  = pd.to_numeric(df_b['死亡数（男）'], errors='coerce')
death_female = pd.to_numeric(df_b['死亡数（女）'], errors='coerce')

# 粗死亡率（人口10万あたり）
df['死亡率_総'] = death_total / pop_total * 100000
df['死亡率_男'] = death_male  / pop_male  * 100000
df['死亡率_女'] = death_female / pop_female * 100000

# 医師数（10万対）= SSDSE-E 医師数 / SSDSE-B 総人口 × 10万
doc_map = dict(zip(
    ssdse_e['都道府県'],
    pd.to_numeric(ssdse_e['医師数'], errors='coerce')
))
df['医師数_総数'] = df['都道府県'].map(doc_map)
df['医師数_10万対'] = df['医師数_総数'] / pop_total * 100000

# 人口密度（人/km²）
area_map = dict(zip(
    ssdse_e['都道府県'],
    pd.to_numeric(ssdse_e['総面積（北方地域及び竹島を除く）'], errors='coerce') / 100
))
df['面積_km2'] = df['都道府県'].map(area_map)
df['人口密度']  = pop_total / df['面積_km2']

df = df.dropna().reset_index(drop=True)
N = len(df)

# 都市/地方分類：人口密度上位10 vs 下位10
density_rank = df['人口密度'].rank(ascending=False)
df['都市分類'] = '中間'
df.loc[density_rank <= 10, '都市分類'] = '都市部（上位10県）'
df.loc[density_rank >= N - 9, '都市分類'] = '地方（下位10県）'

print("=" * 60)
print("■ データ概要（2022年, N =", N, "都道府県）")
print("=" * 60)
print(df[['医師数_10万対', '死亡率_男', '死亡率_女', '人口密度']].describe().round(2))
print(f"\n  都市部（上位10県）: {df[df['都市分類']=='都市部（上位10県）']['都道府県'].tolist()}")
print(f"  地方（下位10県）: {df[df['都市分類']=='地方（下位10県）']['都道府県'].tolist()}")

# ================================================================
# ■ Step1. Pearson相関分析
# ================================================================
print("\n" + "=" * 60)
print("■ Step1. Pearson相関分析（医師数10万対 vs 粗死亡率）")
print("=" * 60)
r_m, p_m = stats.pearsonr(df['医師数_10万対'], df['死亡率_男'])
r_f, p_f = stats.pearsonr(df['医師数_10万対'], df['死亡率_女'])
r_t, p_t = stats.pearsonr(df['医師数_10万対'], df['死亡率_総'])
print(f"  総死亡率 vs 医師数: r={r_t:.4f}, p={p_t:.4f}")
print(f"  男性死亡率 vs 医師数: r={r_m:.4f}, p={p_m:.4f}")
print(f"  女性死亡率 vs 医師数: r={r_f:.4f}, p={p_f:.4f}")
print(f"  ※論文参考値: r≈−0.40（年齢調整死亡率を使用）")
print(f"  ※本スクリプトは粗死亡率を使用（高齢化の影響あり）")

# ================================================================
# ■ Step2. Mann-Whitney U検定（都市部 vs 地方）
# ================================================================
print("\n" + "=" * 60)
print("■ Step2. Mann-Whitney U検定（都市部 vs 地方）")
print("=" * 60)
urban = df[df['都市分類'] == '都市部（上位10県）']['死亡率_男'].values
rural = df[df['都市分類'] == '地方（下位10県）']['死亡率_男'].values
u_stat, p_mw = stats.mannwhitneyu(urban, rural, alternative='two-sided')
print(f"  男性死亡率: U={u_stat:.1f}, p={p_mw:.4f}")
print(f"  → {'有意差あり (p<0.05)' if p_mw < 0.05 else '統計的有意差なし (p≥0.05)'}")

urban_f = df[df['都市分類'] == '都市部（上位10県）']['死亡率_女'].values
rural_f = df[df['都市分類'] == '地方（下位10県）']['死亡率_女'].values
u_f, p_f2 = stats.mannwhitneyu(urban_f, rural_f, alternative='two-sided')
print(f"  女性死亡率: U={u_f:.1f}, p={p_f2:.4f}")

# ================================================================
# ■ 図の生成（4枚）
# ================================================================
color_map = {'都市部（上位10県）': '#C62828', '中間': '#9E9E9E', '地方（下位10県）': '#2E7D32'}

# ──────────────────────────────────────────────────────────
# 図1: 医師数 vs 粗死亡率の散布図（男女別）
# ──────────────────────────────────────────────────────────
fig1, axes1 = plt.subplots(1, 2, figsize=(13, 5))
fig1.suptitle('医師数（人口10万対）と粗死亡率の関係\n（2022年 SSDSE-B/E）',
              fontsize=13, fontweight='bold')

for ax, (death_col, gender) in zip(axes1,
        [('死亡率_男', '男性'), ('死亡率_女', '女性')]):
    for cls_name, color_c in color_map.items():
        mask = df['都市分類'] == cls_name
        ax.scatter(df.loc[mask, '医師数_10万対'], df.loc[mask, death_col],
                   color=color_c, s=60, alpha=0.75, edgecolors='white', linewidth=0.5,
                   label=cls_name, zorder=3)
    z = np.polyfit(df['医師数_10万対'], df[death_col], 1)
    xs = np.linspace(df['医師数_10万対'].min() - 10, df['医師数_10万対'].max() + 10, 100)
    ax.plot(xs, np.poly1d(z)(xs), color='#FF8F00', linewidth=2.0, linestyle='--',
            label='線形近似', alpha=0.8, zorder=4)
    r, p = stats.pearsonr(df['医師数_10万対'], df[death_col])
    sig = '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    ax.set_xlabel('医師数（人口10万あたり）', fontsize=11)
    ax.set_ylabel('粗死亡率（人口10万あたり）', fontsize=11)
    ax.set_title(f'{gender}死亡率  r={r:.3f} {sig}', fontsize=11, fontweight='bold')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_H5_5_fig1_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2025_H5_5_fig1_scatter.png")

# ──────────────────────────────────────────────────────────
# 図2: 都市/地方別の箱ひげ図（Mann-Whitney U検定）
# ──────────────────────────────────────────────────────────
fig2, axes2 = plt.subplots(1, 2, figsize=(13, 5))
fig2.suptitle('都市部・地方別の死亡率・医師数分布（Mann-Whitney U検定）',
              fontsize=13, fontweight='bold')

grp_order = ['都市部（上位10県）', '中間', '地方（下位10県）']
box_colors2 = ['#C62828', '#9E9E9E', '#2E7D32']

ax2a = axes2[0]
data_m = [df[df['都市分類'] == g]['死亡率_男'].values for g in grp_order]
bp = ax2a.boxplot(data_m, labels=['都市部\n(上位10)', '中間', '地方\n(下位10)'],
                  patch_artist=True, medianprops=dict(color='white', linewidth=2))
for patch, col in zip(bp['boxes'], box_colors2):
    patch.set_facecolor(col); patch.set_alpha(0.6)
sig_txt = f'Mann-Whitney U={u_stat:.0f}\np={p_mw:.3f} {"n.s." if p_mw >= 0.05 else "*"}'
ax2a.text(0.97, 0.97, sig_txt, transform=ax2a.transAxes, fontsize=9,
          va='top', ha='right', bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.9))
ax2a.set_ylabel('粗死亡率（男性）', fontsize=11)
ax2a.set_title('男性死亡率の分布', fontsize=11, fontweight='bold')
ax2a.grid(axis='y', alpha=0.3)

ax2b = axes2[1]
data_doc = [df[df['都市分類'] == g]['医師数_10万対'].values for g in grp_order]
bp2 = ax2b.boxplot(data_doc, labels=['都市部\n(上位10)', '中間', '地方\n(下位10)'],
                   patch_artist=True, medianprops=dict(color='white', linewidth=2))
for patch, col in zip(bp2['boxes'], box_colors2):
    patch.set_facecolor(col); patch.set_alpha(0.6)
u_doc, p_doc = stats.mannwhitneyu(
    df[df['都市分類'] == '都市部（上位10県）']['医師数_10万対'].values,
    df[df['都市分類'] == '地方（下位10県）']['医師数_10万対'].values,
    alternative='two-sided'
)
sig_doc = f'Mann-Whitney U={u_doc:.0f}\np={p_doc:.3f} {"*" if p_doc < 0.05 else "n.s."}'
ax2b.text(0.97, 0.97, sig_doc, transform=ax2b.transAxes, fontsize=9,
          va='top', ha='right', bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.9))
ax2b.set_ylabel('医師数（人口10万あたり）', fontsize=11)
ax2b.set_title('医師数の分布', fontsize=11, fontweight='bold')
ax2b.grid(axis='y', alpha=0.3)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_H5_5_fig2_box.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2025_H5_5_fig2_box.png")

# ──────────────────────────────────────────────────────────
# 図3: 人口密度と医師数・死亡率の関係
# ──────────────────────────────────────────────────────────
fig3, axes3 = plt.subplots(1, 2, figsize=(13, 5))
fig3.suptitle('人口密度と医師数・粗死亡率の関係（2022年）', fontsize=13, fontweight='bold')

log_density = np.log(df['人口密度'].values)

ax3a = axes3[0]
for cls_name, color_c in color_map.items():
    mask = df['都市分類'] == cls_name
    ax3a.scatter(log_density[mask], df.loc[mask, '医師数_10万対'],
                 color=color_c, s=60, alpha=0.75, edgecolors='white', linewidth=0.5,
                 label=cls_name, zorder=3)
r3a, _ = stats.pearsonr(log_density, df['医師数_10万対'])
z3a = np.polyfit(log_density, df['医師数_10万対'], 1)
xs3a = np.linspace(log_density.min() - 0.2, log_density.max() + 0.2, 100)
ax3a.plot(xs3a, np.poly1d(z3a)(xs3a), 'k--', linewidth=1.5, alpha=0.6)
ax3a.set_xlabel('log(人口密度)', fontsize=11)
ax3a.set_ylabel('医師数（10万対）', fontsize=11)
ax3a.set_title(f'人口密度 vs 医師数  r={r3a:.3f}', fontsize=11, fontweight='bold')
ax3a.legend(fontsize=8)
ax3a.grid(True, alpha=0.3)

ax3b = axes3[1]
for cls_name, color_c in color_map.items():
    mask = df['都市分類'] == cls_name
    ax3b.scatter(log_density[mask], df.loc[mask, '死亡率_男'],
                 color=color_c, s=60, alpha=0.75, edgecolors='white', linewidth=0.5,
                 label=cls_name, zorder=3)
r3b, _ = stats.pearsonr(log_density, df['死亡率_男'])
z3b = np.polyfit(log_density, df['死亡率_男'], 1)
xs3b = np.linspace(log_density.min() - 0.2, log_density.max() + 0.2, 100)
ax3b.plot(xs3b, np.poly1d(z3b)(xs3b), 'k--', linewidth=1.5, alpha=0.6)
ax3b.set_xlabel('log(人口密度)', fontsize=11)
ax3b.set_ylabel('粗死亡率（男性, 10万対）', fontsize=11)
ax3b.set_title(f'人口密度 vs 男性死亡率  r={r3b:.3f}', fontsize=11, fontweight='bold')
ax3b.legend(fontsize=8)
ax3b.grid(True, alpha=0.3)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_H5_5_fig3_density.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2025_H5_5_fig3_density.png")

# ──────────────────────────────────────────────────────────
# 図4: 都道府県別死亡率のランキング
# ──────────────────────────────────────────────────────────
fig4, axes4 = plt.subplots(1, 2, figsize=(14, 8))
fig4.suptitle('都道府県別粗死亡率（2022年）', fontsize=13, fontweight='bold')

for ax, (col, gender) in zip(axes4, [('死亡率_男', '男性'), ('死亡率_女', '女性')]):
    sorted_idx = np.argsort(df[col].values)
    vals_sorted  = df[col].values[sorted_idx]
    prefs_sorted = df['都道府県'].values[sorted_idx]
    cls_sorted   = df['都市分類'].values[sorted_idx]
    bar_cols4 = ['#C62828' if c == '都市部（上位10県）'
                 else '#2E7D32' if c == '地方（下位10県）'
                 else '#9E9E9E' for c in cls_sorted]
    ax.barh(range(N), vals_sorted, color=bar_cols4, alpha=0.75, edgecolor='white', height=0.7)
    ax.set_yticks(range(N))
    ax.set_yticklabels(prefs_sorted, fontsize=6.5)
    ax.set_xlabel('粗死亡率（人口10万あたり）', fontsize=10)
    ax.set_title(f'{gender}粗死亡率の都道府県ランキング', fontsize=11, fontweight='bold')
    ax.axvline(df[col].mean(), color='black', linestyle='--', linewidth=1.0,
               label=f'平均={df[col].mean():.0f}')
    ax.legend(fontsize=8)
    ax.grid(axis='x', alpha=0.3)

legend_els4 = [Patch(color='#C62828', alpha=0.75, label='都市部（上位10県）'),
               Patch(color='#9E9E9E', alpha=0.75, label='中間'),
               Patch(color='#2E7D32', alpha=0.75, label='地方（下位10県）')]
axes4[0].legend(handles=legend_els4, fontsize=7, loc='lower right')

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_H5_5_fig4_rank.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2025_H5_5_fig4_rank.png")

print("\n全図の生成完了（4枚）")
print(f"\nデータ出典: SSDSE-B-2026.csv（2022年度）, SSDSE-E-2026.csv（統計センター）")