"""
2022_H5_16_shorei.py
====================
教育用ハンズオン分析スクリプト
論文: ごみ排出量とリサイクル率の社会経済的要因：都道府県別パネル分析
受賞: 2022年度 統計データ分析コンペティション 審査員奨励賞（高校生の部）

手法: 相関分析・OLS回帰・時系列分析・散布図
データ: SSDSE-B-2026.csv（都道府県別社会経済統計）

生成図:
  Fig1: ごみ排出量とリサイクル率の時系列推移（全国平均, 2軸グラフ）
  Fig2: 消費支出 vs ごみ排出量 散布図（2022年, 都道府県ラベル付き）
  Fig3: OLS回帰係数プロット（ごみ排出量の決定要因）
  Fig4: リサイクル率の都道府県別ランキング（2022年, 地域色分け）

注意: 実データ（SSDSE-B-2026）のみ使用。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/2022_H5_16_shorei.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 matplotlib import rcParams
import statsmodels.api as sm
from scipy import stats

warnings.filterwarnings('ignore')

# ── パス設定 ─────────────────────────────────────────────────────────────────
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'

os.makedirs(FIG_DIR, exist_ok=True)

# ── フォント設定 ───────────────────────────────────────────────────────────────
rcParams['font.family'] = ['Hiragino Sans', 'Hiragino Kaku Gothic ProN',
                            'AppleGothic', 'Noto Sans CJK JP', 'sans-serif']
rcParams['axes.unicode_minus'] = False
DPI = 150

# ── データ読み込み ─────────────────────────────────────────────────────────────
df_raw = pd.read_csv(DATA_B, encoding='cp932', header=None)
col_names = df_raw.iloc[1].tolist()
df = df_raw.iloc[2:].reset_index(drop=True)
df.columns = col_names

# 数値変換
numeric_cols = [
    '年度', '総人口', '65歳以上人口',
    'ごみ総排出量（総量）', '1人1日当たりの排出量', 'ごみのリサイクル率',
    '旅館営業施設数（ホテルを含む）', '消費支出（二人以上の世帯）',
]
for c in numeric_cols:
    df[c] = pd.to_numeric(df[c], errors='coerce')

df = df.dropna(subset=['1人1日当たりの排出量', 'ごみのリサイクル率'])

# ── 派生変数の作成 ─────────────────────────────────────────────────────────────
# 高齢化率（%）
df['高齢化率'] = df['65歳以上人口'] / df['総人口'] * 100

# 消費支出の対数（円 → log(円)）
df['消費支出_log'] = np.log(df['消費支出（二人以上の世帯）'])

# 旅館密度（1万人あたり施設数）
df['旅館密度'] = df['旅館営業施設数（ホテルを含む）'] / (df['総人口'] / 10000)

# 変数名を分析用に短縮
df['ごみ排出量_g'] = df['1人1日当たりの排出量']       # g/人/日
df['リサイクル率']  = df['ごみのリサイクル率']          # %

# ── 地域ブロック定義 ───────────────────────────────────────────────────────────
REGION_MAP = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部',
    '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国',
    '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国',
    '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
REGION_COLORS = {
    '北海道・東北': '#1565C0',
    '関東':         '#C62828',
    '中部':         '#F9A825',
    '近畿':         '#2E7D32',
    '中国・四国':   '#6A1B9A',
    '九州・沖縄':   '#00838F',
}
df['地域'] = df['都道府県'].map(REGION_MAP).fillna('その他')

print("データ読み込み完了")
print(f"  年度: {sorted(df['年度'].dropna().astype(int).unique())}")
print(f"  都道府県数（2022年）: {len(df[df['年度']==2022])}")

# ═══════════════════════════════════════════════════════════════════════════════
# Fig1: ごみ排出量とリサイクル率の時系列推移（全国平均, 2軸グラフ）
# ═══════════════════════════════════════════════════════════════════════════════
print("\n[Fig1] 時系列推移（twinx）...")

ts = (df.groupby('年度')[['ごみ排出量_g', 'リサイクル率']]
        .mean()
        .reset_index()
        .sort_values('年度'))

fig, ax1 = plt.subplots(figsize=(9, 5), dpi=DPI)
ax2 = ax1.twinx()

color1 = '#1565C0'
color2 = '#C62828'

l1, = ax1.plot(ts['年度'], ts['ごみ排出量_g'],
               color=color1, linewidth=2.5, marker='o', markersize=7,
               label='1人1日当たりごみ排出量（g）')
ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('1人1日当たりごみ排出量（g/人/日）', color=color1, fontsize=11)
ax1.tick_params(axis='y', labelcolor=color1)
ax1.set_ylim(820, 1020)

l2, = ax2.plot(ts['年度'], ts['リサイクル率'],
               color=color2, linewidth=2.5, marker='s', markersize=7,
               linestyle='--', label='ごみのリサイクル率（%）')
ax2.set_ylabel('ごみのリサイクル率（%）', color=color2, fontsize=11)
ax2.tick_params(axis='y', labelcolor=color2)
ax2.set_ylim(16, 24)

ax1.set_title('ごみ排出量とリサイクル率の時系列推移（全国47都道府県平均）\n2012〜2023年度', fontsize=13, fontweight='bold', pad=14)
ax1.set_xticks(ts['年度'])
ax1.set_xticklabels([str(int(y)) for y in ts['年度']], rotation=45, ha='right')
ax1.grid(axis='y', linestyle='--', alpha=0.4)

lines = [l1, l2]
labels = [l.get_label() for l in lines]
ax1.legend(lines, labels, loc='upper right', fontsize=10,
           framealpha=0.9, edgecolor='#ccc')

# 注釈
ax1.annotate('ごみ排出量は\n一貫して減少傾向',
             xy=(2017, ts.loc[ts['年度']==2017, 'ごみ排出量_g'].values[0]),
             xytext=(2014, 875),
             fontsize=9, color=color1,
             arrowprops=dict(arrowstyle='->', color=color1, lw=1.2))

plt.tight_layout()
out1 = os.path.join(FIG_DIR, '2022_H5_16_fig1_timeseries.png')
plt.savefig(out1, dpi=DPI, bbox_inches='tight')
plt.close()
print(f"  保存: {out1}")

# ═══════════════════════════════════════════════════════════════════════════════
# Fig2: 消費支出 vs ごみ排出量 散布図（2022年, 都道府県ラベル付き）
# ═══════════════════════════════════════════════════════════════════════════════
print("\n[Fig2] 散布図（消費支出 vs ごみ排出量）...")

df22 = df[df['年度'] == 2022].dropna(subset=['消費支出_log', 'ごみ排出量_g', '地域']).copy()

fig, ax = plt.subplots(figsize=(11, 8), dpi=DPI)

for region, grp in df22.groupby('地域'):
    color = REGION_COLORS.get(region, '#999')
    ax.scatter(grp['消費支出_log'], grp['ごみ排出量_g'],
               color=color, s=80, alpha=0.85, zorder=3, label=region,
               edgecolors='white', linewidths=0.6)

# 都道府県ラベル（省略形）
LABEL_MAP = {
    '北海道': '北海道', '東京都': '東京', '大阪府': '大阪', '愛知県': '愛知',
    '神奈川県': '神奈川', '沖縄県': '沖縄', '高知県': '高知', '秋田県': '秋田',
    '富山県': '富山', '山形県': '山形', '福井県': '福井', '徳島県': '徳島',
}
for _, row in df22.iterrows():
    pref = row['都道府県']
    short = LABEL_MAP.get(pref, pref.replace('県', '').replace('府', '').replace('都', ''))
    ax.annotate(short,
                xy=(row['消費支出_log'], row['ごみ排出量_g']),
                xytext=(3, 3), textcoords='offset points',
                fontsize=7.5, color='#333', alpha=0.9)

# 回帰直線
x_vals = df22['消費支出_log'].dropna()
y_vals = df22['ごみ排出量_g'].dropna()
idx = df22[['消費支出_log', 'ごみ排出量_g']].dropna().index
xv = df22.loc[idx, '消費支出_log']
yv = df22.loc[idx, 'ごみ排出量_g']
slope, intercept, r_val, p_val, _ = stats.linregress(xv, yv)
x_line = np.linspace(xv.min(), xv.max(), 100)
y_line = slope * x_line + intercept
ax.plot(x_line, y_line, '--', color='#555', linewidth=1.5, alpha=0.7,
        label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})')

ax.set_xlabel('log（消費支出：二人以上の世帯, 円）', fontsize=12)
ax.set_ylabel('1人1日当たりごみ排出量（g/人/日）', fontsize=12)
ax.set_title('消費支出とごみ排出量の関係（2022年度・47都道府県）', fontsize=13, fontweight='bold', pad=12)
ax.legend(loc='upper left', fontsize=9, framealpha=0.9, edgecolor='#ccc')
ax.grid(linestyle='--', alpha=0.35)

plt.tight_layout()
out2 = os.path.join(FIG_DIR, '2022_H5_16_fig2_scatter.png')
plt.savefig(out2, dpi=DPI, bbox_inches='tight')
plt.close()
print(f"  保存: {out2}")
print(f"  相関: r={r_val:.3f}, p={p_val:.4f}")

# ═══════════════════════════════════════════════════════════════════════════════
# Fig3: OLS回帰係数プロット（ごみ排出量の決定要因）
# ═══════════════════════════════════════════════════════════════════════════════
print("\n[Fig3] OLS回帰係数プロット...")

# 2022年データで回帰分析
reg_df = df[df['年度'] == 2022].dropna(
    subset=['ごみ排出量_g', '高齢化率', '消費支出_log', '旅館密度', 'リサイクル率']
).copy()

# 説明変数の標準化（係数の比較のため）
VARS = {
    '高齢化率（%）': '高齢化率',
    'log消費支出（円）': '消費支出_log',
    '旅館密度（1万人対）': '旅館密度',
}

X_data = {}
for label, col in VARS.items():
    X_data[label] = (reg_df[col] - reg_df[col].mean()) / reg_df[col].std()

X_mat = pd.DataFrame(X_data)
y_vec = reg_df['ごみ排出量_g']

X_sm = sm.add_constant(X_mat)
model = sm.OLS(y_vec, X_sm).fit()

print(model.summary())

# 係数・信頼区間・p値
var_labels = list(VARS.keys())
coefs   = model.params[1:].values
conf    = model.conf_int().iloc[1:].values   # [lower, upper]
pvals   = model.pvalues[1:].values

fig, ax = plt.subplots(figsize=(9, 5), dpi=DPI)

colors = ['#C62828' if p < 0.05 else '#9E9E9E' for p in pvals]
y_pos  = range(len(var_labels))

for i, (label, coef, ci, p, col) in enumerate(
        zip(var_labels, coefs, conf, pvals, colors)):
    ax.barh(i, coef, color=col, alpha=0.8, height=0.55, zorder=3)
    ax.errorbar(coef, i,
                xerr=[[coef - ci[0]], [ci[1] - coef]],
                fmt='none', color='#333', capsize=5, linewidth=1.5, zorder=4)
    sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))
    ax.text(coef + (2 if coef >= 0 else -2), i,
            f'{sig}  p={p:.3f}', va='center', fontsize=9,
            color='#333', ha='left' if coef >= 0 else 'right')

ax.axvline(0, color='#555', linewidth=1.2, linestyle='--')
ax.set_yticks(list(y_pos))
ax.set_yticklabels(var_labels, fontsize=11)
ax.set_xlabel('標準化回帰係数（95%CI）', fontsize=11)
ax.set_title(f'ごみ排出量の決定要因（OLS回帰, 2022年度, N={len(reg_df)}都道府県）\n'
             f'R²={model.rsquared:.3f}', fontsize=12, fontweight='bold', pad=12)
ax.grid(axis='x', linestyle='--', alpha=0.35)

sig_patch  = mpatches.Patch(color='#C62828', alpha=0.8, label='有意（p<0.05）')
ns_patch   = mpatches.Patch(color='#9E9E9E', alpha=0.8, label='非有意（p≥0.05）')
ax.legend(handles=[sig_patch, ns_patch], loc='lower right', fontsize=10)

plt.tight_layout()
out3 = os.path.join(FIG_DIR, '2022_H5_16_fig3_coef.png')
plt.savefig(out3, dpi=DPI, bbox_inches='tight')
plt.close()
print(f"  保存: {out3}")

# ═══════════════════════════════════════════════════════════════════════════════
# Fig4: リサイクル率の都道府県別ランキング（2022年, 地域色分け）
# ═══════════════════════════════════════════════════════════════════════════════
print("\n[Fig4] リサイクル率都道府県別ランキング...")

df22_sorted = df22.dropna(subset=['リサイクル率']).sort_values('リサイクル率', ascending=True).reset_index(drop=True)

fig, ax = plt.subplots(figsize=(10, 13), dpi=DPI)

bar_colors = [REGION_COLORS.get(r, '#999') for r in df22_sorted['地域']]
bars = ax.barh(range(len(df22_sorted)), df22_sorted['リサイクル率'],
               color=bar_colors, alpha=0.85, height=0.7, zorder=3)

# 都道府県ラベル
for i, (_, row) in enumerate(df22_sorted.iterrows()):
    pref = row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax.text(row['リサイクル率'] + 0.3, i, f"{pref}  {row['リサイクル率']:.1f}%",
            va='center', fontsize=8.5, color='#333')

ax.set_yticks(range(len(df22_sorted)))
ax.set_yticklabels(
    [row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
     for _, row in df22_sorted.iterrows()],
    fontsize=8
)
ax.set_xlabel('ごみのリサイクル率（%）', fontsize=11)
ax.set_title('都道府県別ごみのリサイクル率ランキング（2022年度）', fontsize=13, fontweight='bold', pad=12)
ax.axvline(df22_sorted['リサイクル率'].mean(), color='#F9A825', linewidth=2,
           linestyle='--', label=f"全国平均 {df22_sorted['リサイクル率'].mean():.1f}%")
ax.grid(axis='x', linestyle='--', alpha=0.35)
ax.set_xlim(0, df22_sorted['リサイクル率'].max() * 1.18)

# 凡例（地域）
handles = [mpatches.Patch(color=col, label=region, alpha=0.85)
           for region, col in REGION_COLORS.items()]
handles.append(mpatches.Patch(color='#F9A825', label='全国平均', alpha=0.9))
ax.legend(handles=handles, loc='lower right', fontsize=8.5,
          framealpha=0.9, edgecolor='#ccc')

plt.tight_layout()
out4 = os.path.join(FIG_DIR, '2022_H5_16_fig4_ranking.png')
plt.savefig(out4, dpi=DPI, bbox_inches='tight')
plt.close()
print(f"  保存: {out4}")

# ── サマリー出力 ───────────────────────────────────────────────────────────────
print("\n===== 分析サマリー（2022年度） =====")
print(f"  ごみ排出量（平均）: {df22['ごみ排出量_g'].mean():.1f} g/人/日")
print(f"  リサイクル率（平均）: {df22['リサイクル率'].mean():.1f} %")
print(f"  リサイクル率 最高: {df22.loc[df22['リサイクル率'].idxmax(), '都道府県']} "
      f"({df22['リサイクル率'].max():.1f}%)")
print(f"  リサイクル率 最低: {df22.loc[df22['リサイクル率'].idxmin(), '都道府県']} "
      f"({df22['リサイクル率'].min():.1f}%)")

# 相関分析
corr_vars = ['高齢化率', '消費支出_log', '旅館密度']
print("\n  ごみ排出量との相関係数:")
for v in corr_vars:
    sub = df22.dropna(subset=[v, 'ごみ排出量_g'])
    r, p = stats.pearsonr(sub[v], sub['ごみ排出量_g'])
    print(f"    {v}: r={r:.3f}, p={p:.4f}")

print("\n全図の生成が完了しました。")
print(f"保存先: {os.path.abspath(FIG_DIR)}")