"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部]
=================================================================
論文タイトル：味覚から眺める地域別「ふるさとの味」
             -関西が薄味，関東が濃い味は本当なのか-

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別データ，2012〜2023年）
  指標：食料費割合（食料費 / 消費支出）を「味の濃淡」の代理変数として使用

  Step1. 6地域別 食料費割合の箱ひげ図（Kruskal-Wallis検定）
  Step2. 関東 vs 関西の複数指標比較棒グラフ
  Step3. 食料費割合の時系列（関東・関西・全国平均）
  Step4. 住宅地価格 vs 食料費割合の散布図（地域カラー）

【地域区分（SSDSE-Bの地域コード）】
  北海道   : R01000
  東北     : R02000 - R07000（青森〜福島）
  関東     : R08000 - R14000（茨城〜神奈川）
  中部     : R15000 - R23000（新潟〜愛知）
  関西（近畿）: R24000 - R30000（三重〜和歌山）
  中国・四国: R31000 - R39000（鳥取〜高知）
  九州・沖縄: R40000 - R47000（福岡〜沖縄）

【分析の焦点：関東 vs 関西】
  関東: 東京・神奈川・埼玉・千葉・茨城・栃木・群馬（R08〜R14）
  関西: 大阪・京都・兵庫・奈良・滋賀・和歌山・三重（R24〜R30）
  ※ 論文の問い：「関西が薄味，関東が濃い味」は食費データに現れるか？

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計数理研究所より公表の実データ

【データサイエンス学習ポイント】
  1. 食料費割合（食料費/消費支出）を「味の濃淡代理変数」として解釈
  2. Kruskal-Wallis検定（ノンパラメトリック地域差検定）
  3. 時系列データの地域別比較
  4. 住宅地価格と食費行動の相関分析
=================================================================
"""

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


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats
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)

# ================================================================
# ■ 実データ読み込み（SSDSE-B-2026：都道府県データ）
# ================================================================
df_b_raw = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b_raw[df_b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

print("=" * 65)
print("■ SSDSE-B-2026 都道府県データ読み込み")
print("=" * 65)
print(f"  行数: {len(df_b)}, 年度範囲: {df_b['年度'].min()}〜{df_b['年度'].max()}")
print(f"  都道府県数（年平均）: {df_b['年度'].value_counts().median():.0f}")

# ──────────────────────────────────────────────────────────────
# 地域コードからプレフィックス番号を取得
# ──────────────────────────────────────────────────────────────
df_b['pref_num'] = df_b['地域コード'].str.extract(r'R(\d{2})').astype(int)

# 地域区分の定義
def assign_region(n):
    if n == 1:
        return '北海道'
    elif 2 <= n <= 7:
        return '東北'
    elif 8 <= n <= 14:
        return '関東'
    elif 15 <= n <= 23:
        return '中部'
    elif 24 <= n <= 30:
        return '関西'
    elif 31 <= n <= 39:
        return '中国・四国'
    elif 40 <= n <= 47:
        return '九州・沖縄'
    else:
        return 'その他'

df_b['地域'] = df_b['pref_num'].apply(assign_region)

# ──────────────────────────────────────────────────────────────
# 主要指標の計算
# ──────────────────────────────────────────────────────────────
# 食料費割合（食料費 / 消費支出）: 論文の中心指標
df_b['食料費割合'] = df_b['食料費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100

# 光熱・水道費割合（調理・生活エネルギー使用の代理変数）
df_b['光熱・水道費割合'] = df_b['光熱・水道費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100

# 教養娯楽費割合（生活スタイルの指標）
df_b['教養娯楽費割合'] = df_b['教養娯楽費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100

print(f"\n  食料費割合の基本統計:")
print(df_b['食料費割合'].describe().round(2))

# 関東・関西フラグ
KANTO_PREFS  = list(range(8, 15))   # R08〜R14（茨城〜神奈川）
KANSAI_PREFS = list(range(24, 31))  # R24〜R30（三重〜和歌山）

df_b['地域2'] = df_b['pref_num'].apply(
    lambda n: '関東' if n in KANTO_PREFS else ('関西' if n in KANSAI_PREFS else 'その他')
)

# ──────────────────────────────────────────────────────────────
# 最新年度のデータ（図2・図4用）
# ──────────────────────────────────────────────────────────────
latest_year = df_b['年度'].max()
df_latest = df_b[df_b['年度'] == latest_year].copy()
print(f"\n  最新年度: {latest_year}年 ({len(df_latest)}都道府県)")

# ──────────────────────────────────────────────────────────────
# Kruskal-Wallis検定（6地域間の食料費割合の差）
# ──────────────────────────────────────────────────────────────
print("\n" + "=" * 65)
print("■ Kruskal-Wallis検定：6地域別 食料費割合の差")
print("=" * 65)

REGION_ORDER = ['北海道', '東北', '関東', '中部', '関西', '中国・四国', '九州・沖縄']

# 全年度データを使って地域別に集計
region_groups = [
    df_b[df_b['地域'] == r]['食料費割合'].dropna().values
    for r in REGION_ORDER
]

kw_stat, kw_p = stats.kruskal(*region_groups)
print(f"  H統計量={kw_stat:.4f}, p値={kw_p:.6f}")
sig_kw = '*** (p<0.001)' if kw_p < 0.001 else '** (p<0.01)' if kw_p < 0.01 else '* (p<0.05)' if kw_p < 0.05 else 'n.s.'
print(f"  判定: {sig_kw}")

# 地域別の記述統計
print(f"\n  {'地域':<12} {'N':>5} {'中央値':>8} {'平均':>8} {'SD':>7}")
print("  " + "-" * 45)
for r, grp in zip(REGION_ORDER, region_groups):
    print(f"  {r:<12} {len(grp):>5} {np.median(grp):>8.2f} {np.mean(grp):>8.2f} {np.std(grp):>7.2f}")

# ──────────────────────────────────────────────────────────────
# 関東 vs 関西 Mann-Whitney U検定
# ──────────────────────────────────────────────────────────────
print("\n" + "=" * 65)
print("■ Mann-Whitney U検定：関東 vs 関西（全年度）")
print("=" * 65)

kanto_vals  = df_b[df_b['地域2'] == '関東']['食料費割合'].dropna().values
kansai_vals = df_b[df_b['地域2'] == '関西']['食料費割合'].dropna().values

u_stat, u_p = stats.mannwhitneyu(kanto_vals, kansai_vals, alternative='two-sided')
print(f"  U統計量={u_stat:.1f}, p値={u_p:.6f}")
sig_u = '*** (p<0.001)' if u_p < 0.001 else '** (p<0.01)' if u_p < 0.01 else '* (p<0.05)' if u_p < 0.05 else 'n.s.'
print(f"  判定: {sig_u}")
print(f"  関東 中央値={np.median(kanto_vals):.2f}%, 平均={np.mean(kanto_vals):.2f}%")
print(f"  関西 中央値={np.median(kansai_vals):.2f}%, 平均={np.mean(kansai_vals):.2f}%")

# 最新年度の関東・関西比較
kanto_latest  = df_latest[df_latest['地域2'] == '関東']
kansai_latest = df_latest[df_latest['地域2'] == '関西']
metrics = {
    '食料費割合（%）':    ('食料費割合', 2),
    '光熱・水道費割合（%）': ('光熱・水道費割合', 2),
    '教養娯楽費割合（%）': ('教養娯楽費割合', 2),
}
print(f"\n  {latest_year}年 関東 vs 関西 比較:")
print(f"  {'指標':<22} {'関東':>10} {'関西':>10}")
print("  " + "-" * 45)
for label, (col, nd) in metrics.items():
    v_k = kanto_latest[col].mean()
    v_ks = kansai_latest[col].mean()
    print(f"  {label:<22} {v_k:>10.{nd}f} {v_ks:>10.{nd}f}")

# ──────────────────────────────────────────────────────────────
# 相関分析：住宅地価格 vs 食料費割合（最新年度）
# ──────────────────────────────────────────────────────────────
print("\n" + "=" * 65)
print("■ 相関分析：住宅地価格 vs 食料費割合（最新年度）")
print("=" * 65)

scatter_df = df_latest[['都道府県', '地域', '食料費割合',
                          '標準価格（平均価格）（住宅地）']].dropna()
r_scatter, p_scatter = stats.pearsonr(
    scatter_df['標準価格（平均価格）（住宅地）'],
    scatter_df['食料費割合']
)
print(f"  r={r_scatter:.4f}, p={p_scatter:.4f}")
sig_s = '*** (p<0.001)' if p_scatter < 0.001 else '** (p<0.01)' if p_scatter < 0.01 else '* (p<0.05)' if p_scatter < 0.05 else 'n.s.'
print(f"  判定: {sig_s}")

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

# 地域別カラーパレット
REGION_COLORS = {
    '北海道':    '#1565C0',
    '東北':      '#2E7D32',
    '関東':      '#C62828',
    '中部':      '#E65100',
    '関西':      '#6A1B9A',
    '中国・四国': '#00695C',
    '九州・沖縄': '#F9A825',
}

# ────────────────────────────────────────────────────────────────
# 図1: 6地域別 食料費割合 箱ひげ図
# ────────────────────────────────────────────────────────────────
fig1, ax1 = plt.subplots(figsize=(11, 6))

bp_data = [grp for grp in region_groups]
bp = ax1.boxplot(
    bp_data,
    labels=REGION_ORDER,
    patch_artist=True,
    medianprops=dict(color='black', linewidth=2),
    flierprops=dict(marker='o', markersize=4, alpha=0.6),
    boxprops=dict(linewidth=1.2),
    whiskerprops=dict(linewidth=1.2),
    capprops=dict(linewidth=1.5),
)

for patch, region in zip(bp['boxes'], REGION_ORDER):
    c = REGION_COLORS[region]
    patch.set_facecolor(c)
    patch.set_alpha(0.65)

for flier, region in zip(bp['fliers'], REGION_ORDER):
    flier.set(markerfacecolor=REGION_COLORS[region], alpha=0.6)

ax1.set_ylabel('食料費割合（食料費 / 消費支出，%）', fontsize=12)
ax1.set_xlabel('地域', fontsize=12)
ax1.set_title(
    f'6地域別 食料費割合の分布（2012〜{latest_year}年）\n'
    f'Kruskal-Wallis: H={kw_stat:.2f}, {sig_kw}',
    fontsize=13, fontweight='bold'
)
ax1.grid(axis='y', alpha=0.3)

# 各地域の中央値をテキスト表示
for i, grp in enumerate(region_groups, start=1):
    med = np.median(grp)
    ax1.text(i, med + 0.1, f'{med:.1f}%', ha='center', va='bottom', fontsize=9,
             fontweight='bold', color=REGION_COLORS[REGION_ORDER[i-1]])

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

# ────────────────────────────────────────────────────────────────
# 図2: 関東 vs 関西 比較棒グラフ（複数指標）
# ────────────────────────────────────────────────────────────────
fig2, axes2 = plt.subplots(1, 3, figsize=(13, 5))
fig2.suptitle(
    f'関東 vs 関西 消費支出内訳の比較（{latest_year}年，各都道府県平均）',
    fontsize=13, fontweight='bold'
)

metric_list = [
    ('食料費割合', '食料費割合\n（食料費 / 消費支出，%）',  '#C62828', '#6A1B9A'),
    ('光熱・水道費割合', '光熱・水道費割合\n（光熱・水道費 / 消費支出，%）', '#E65100', '#9C27B0'),
    ('教養娯楽費割合', '教養娯楽費割合\n（教養娯楽費 / 消費支出，%）', '#1565C0', '#00695C'),
]

for ax, (col, ylabel, c_kanto, c_kansai) in zip(axes2, metric_list):
    v_kanto  = kanto_latest[col].mean()
    v_kansai = kansai_latest[col].mean()
    se_kanto  = kanto_latest[col].sem()
    se_kansai = kansai_latest[col].sem()

    bars = ax.bar(['関東', '関西'],
                  [v_kanto, v_kansai],
                  color=[c_kanto, c_kansai], alpha=0.75,
                  edgecolor='white', linewidth=1.2, width=0.5)
    ax.errorbar(['関東', '関西'], [v_kanto, v_kansai],
                yerr=[1.96*se_kanto, 1.96*se_kansai],
                fmt='none', color='#333', capsize=5, linewidth=1.5)

    # 値のラベル
    for bar, val in zip(bars, [v_kanto, v_kansai]):
        ax.text(bar.get_x() + bar.get_width() / 2,
                bar.get_height() + 0.05,
                f'{val:.2f}%', ha='center', va='bottom',
                fontsize=11, fontweight='bold')

    # Mann-Whitney検定
    col_kanto_series  = kanto_latest[col].dropna().values
    col_kansai_series = kansai_latest[col].dropna().values
    if len(col_kanto_series) >= 2 and len(col_kansai_series) >= 2:
        _, p_tmp = stats.mannwhitneyu(col_kanto_series, col_kansai_series, alternative='two-sided')
        sig_tmp = '***' if p_tmp < 0.001 else '**' if p_tmp < 0.01 else '*' if p_tmp < 0.05 else 'n.s.'
        ax.set_title(f'Mann-Whitney: {sig_tmp}', fontsize=10)

    ax.set_ylabel(ylabel, fontsize=9)
    ax.set_ylim(0, max(v_kanto, v_kansai) * 1.25)
    ax.grid(axis='y', alpha=0.3)
    ax.tick_params(axis='x', labelsize=12)

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

# ────────────────────────────────────────────────────────────────
# 図3: 食料費割合の時系列（関東・関西・全国平均）
# ────────────────────────────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(11, 5))

years = sorted(df_b['年度'].unique())

# 全国平均
national_ts = df_b.groupby('年度')['食料費割合'].mean()

# 関東平均
kanto_ts  = df_b[df_b['地域2'] == '関東'].groupby('年度')['食料費割合'].mean()
kansai_ts = df_b[df_b['地域2'] == '関西'].groupby('年度')['食料費割合'].mean()

# 関東・関西の信頼区間（SE×1.96）
kanto_sem  = df_b[df_b['地域2'] == '関東'].groupby('年度')['食料費割合'].sem()
kansai_sem = df_b[df_b['地域2'] == '関西'].groupby('年度')['食料費割合'].sem()

ax3.fill_between(kanto_ts.index,
                 kanto_ts - 1.96 * kanto_sem,
                 kanto_ts + 1.96 * kanto_sem,
                 color='#C62828', alpha=0.12)
ax3.fill_between(kansai_ts.index,
                 kansai_ts - 1.96 * kansai_sem,
                 kansai_ts + 1.96 * kansai_sem,
                 color='#6A1B9A', alpha=0.12)

ax3.plot(national_ts.index, national_ts.values,
         color='#455A64', linewidth=2.0, linestyle='--',
         marker='o', markersize=5, label='全国平均')
ax3.plot(kanto_ts.index, kanto_ts.values,
         color='#C62828', linewidth=2.5,
         marker='s', markersize=6, label='関東（7都県平均）')
ax3.plot(kansai_ts.index, kansai_ts.values,
         color='#6A1B9A', linewidth=2.5,
         marker='^', markersize=6, label='関西（7府県平均）')

ax3.set_xlabel('年度', fontsize=12)
ax3.set_ylabel('食料費割合（食料費 / 消費支出，%）', fontsize=12)
ax3.set_title(
    '食料費割合の時系列推移：関東・関西・全国平均（SSDSE-B 実データ）',
    fontsize=13, fontweight='bold'
)
ax3.set_xticks(years)
ax3.set_xticklabels([str(y) for y in years], rotation=45, fontsize=9)
ax3.legend(fontsize=11)
ax3.grid(alpha=0.3)

# 直近年の差を注釈
diff_latest = kanto_ts[latest_year] - kansai_ts[latest_year]
arrow_x = latest_year
ax3.annotate(
    f'差={diff_latest:+.2f}pp',
    xy=(arrow_x, kanto_ts[latest_year]),
    xytext=(arrow_x - 1.5, kanto_ts[latest_year] + 0.3),
    fontsize=9, color='#333',
    arrowprops=dict(arrowstyle='->', color='#333', lw=1.2)
)

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

# ────────────────────────────────────────────────────────────────
# 図4: 住宅地価格 vs 食料費割合 散布図（地域カラー）
# ────────────────────────────────────────────────────────────────
fig4, ax4 = plt.subplots(figsize=(10, 7))

for region in REGION_ORDER:
    sub = scatter_df[scatter_df['地域'] == region]
    ax4.scatter(
        sub['標準価格（平均価格）（住宅地）'] / 1000,   # 千円/m² 単位
        sub['食料費割合'],
        color=REGION_COLORS[region],
        s=80, alpha=0.75, zorder=3,
        label=region, edgecolors='white', linewidths=0.5
    )

# 全体回帰直線
x_sc = scatter_df['標準価格（平均価格）（住宅地）'].values / 1000
y_sc = scatter_df['食料費割合'].values
z_sc = np.polyfit(x_sc, y_sc, 1)
xs_sc = np.linspace(x_sc.min(), x_sc.max(), 200)
ax4.plot(xs_sc, np.poly1d(z_sc)(xs_sc),
         color='#333', linewidth=1.8, linestyle='--', zorder=2,
         label=f'回帰直線 r={r_scatter:.3f} ({sig_s})')

# 注目都道府県をラベル表示
highlight_prefs = ['東京都', '大阪府', '京都府', '神奈川県', '北海道', '沖縄県', '秋田県']
for _, row in scatter_df.iterrows():
    if row['都道府県'] in highlight_prefs:
        ax4.annotate(
            row['都道府県'],
            xy=(row['標準価格（平均価格）（住宅地）'] / 1000, row['食料費割合']),
            xytext=(6, 4), textcoords='offset points',
            fontsize=8.5, color='#333',
        )

ax4.set_xlabel('標準価格・住宅地（千円/m²）', fontsize=12)
ax4.set_ylabel('食料費割合（食料費 / 消費支出，%）', fontsize=12)
ax4.set_title(
    f'住宅地価格 vs 食料費割合の関係（{latest_year}年，47都道府県）\n'
    f'相関係数 r={r_scatter:.3f}，{sig_s}',
    fontsize=13, fontweight='bold'
)
ax4.legend(fontsize=9, loc='upper right', framealpha=0.85)
ax4.grid(alpha=0.3)

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

print("\n全図の生成完了（4枚）")
print(f"  2022_U5_4_fig1_region.png      : 6地域別 食料費割合 箱ひげ図")
print(f"  2022_U5_4_fig2_kanto_kansai.png: 関東 vs 関西 複数指標比較棒グラフ")
print(f"  2022_U5_4_fig3_ts.png          : 食料費割合の時系列推移")
print(f"  2022_U5_4_fig4_scatter.png     : 住宅地価格 vs 食料費割合 散布図")

print("\n" + "=" * 65)
print("■ 分析サマリー（論文の問いへの答え）")
print("=" * 65)
print(f"  ・6地域間の食料費割合の差 → Kruskal-Wallis {sig_kw}")
print(f"  ・関東 vs 関西（食料費割合）→ Mann-Whitney {sig_u}")
print(f"  ・関東 食料費割合平均: {np.mean(kanto_vals):.2f}%")
print(f"  ・関西 食料費割合平均: {np.mean(kansai_vals):.2f}%")
diff_mean = np.mean(kanto_vals) - np.mean(kansai_vals)
print(f"  ・差（関東－関西）: {diff_mean:+.2f}pp")
print(f"  ・住宅地価格 vs 食料費割合 相関: r={r_scatter:.3f} {sig_s}")
print("=" * 65)