"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部]
=================================================================
論文タイトル：生活時間と幸福感：SSDSE-D「社会生活基本調査」を用いた年代別分析
受賞：審査員奨励賞（大学生・一般の部）

【分析概要】
  データ：SSDSE-D-2023.csv（社会生活基本調査 都道府県別生活時間, 2023年公開）
  対象：全47都道府県（+全国計）× 男女別（総数・男・女）
  単位：1週間の平均時間（分）

  分析の流れ
  1. 生活時間配分の都道府県比較（仕事・家事・余暇の三角形）
  2. 性別による生活時間差の可視化（男女比較棒グラフ）
  3. 都道府県別「仕事時間」と「趣味・娯楽時間」の相関・散布図
  4. 主成分分析 (PCA) による生活時間パターンの可視化

【主要時間活動変数】
  睡眠, 身の回りの用事, 食事, 通勤・通学, 仕事, 家事, 介護・看護, 育児,
  買い物, 移動(通勤・通学を除く), テレビ・ラジオ・新聞・雑誌, 休養・くつろぎ,
  学習・自己啓発・訓練(学業以外), 趣味・娯楽, スポーツ, 交際・付き合い

【データ出典】
  SSDSE-D-2023.csv: 社会・人口統計体系（都道府県別生活時間）
  総務省統計局 社会生活基本調査（2021年実施）をもとに作成
=================================================================
"""

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


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.patches import FancyArrowPatch
import warnings
warnings.filterwarnings('ignore')

from scipy import stats
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

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

import os
FIG_DIR = 'html/figures'
DATA_D  = 'data/raw/SSDSE-D-2023.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ------------------------------------------------------------------ #
# データ読み込み
# ------------------------------------------------------------------ #
df = pd.read_csv(DATA_D, encoding='cp932', header=1)

# 都道府県別（全国計 R00000 を除く）
total = df[(df['男女の別'] == '0_総数') & (df['地域コード'] != 'R00000')].copy()
male  = df[(df['男女の別'] == '1_男')  & (df['地域コード'] != 'R00000')].copy()
female= df[(df['男女の別'] == '2_女')  & (df['地域コード'] != 'R00000')].copy()

# 全国計
nat_total  = df[(df['男女の別'] == '0_総数') & (df['地域コード'] == 'R00000')].iloc[0]
nat_male   = df[(df['男女の別'] == '1_男')   & (df['地域コード'] == 'R00000')].iloc[0]
nat_female = df[(df['男女の別'] == '2_女')   & (df['地域コード'] == 'R00000')].iloc[0]

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

# ------------------------------------------------------------------ #
# Fig1: 都道府県別 生活時間配分 積み上げ横棒グラフ（上位・下位 仕事時間）
# ------------------------------------------------------------------ #
time_cols_main = ['睡眠', '仕事', '家事', '育児', '買い物',
                  'テレビ・ラジオ・新聞・雑誌', '休養・くつろぎ', '趣味・娯楽', 'スポーツ', '交際・付き合い']
colors_main = ['#5C85D6', '#E05C5C', '#F4A261', '#76C442', '#9B59B6',
               '#3D9970', '#F9A825', '#1ABC9C', '#E74C3C', '#95A5A6']

# 仕事時間で並び替え（47都道府県）
df_sorted = total.sort_values('仕事', ascending=True).reset_index(drop=True)

fig, ax = plt.subplots(figsize=(10, 13))
lefts = np.zeros(len(df_sorted))
for col, color in zip(time_cols_main, colors_main):
    vals = df_sorted[col].values
    ax.barh(df_sorted['都道府県'], vals, left=lefts, color=color, label=col, alpha=0.85)
    lefts += vals

ax.set_xlabel('1週間の合計時間（分）', fontsize=12)
ax.set_title('都道府県別 主要生活時間配分\n（仕事時間の短い順）', fontsize=14, fontweight='bold')
ax.legend(loc='lower right', fontsize=8, ncol=2)
ax.grid(axis='x', alpha=0.3)
ax.tick_params(axis='y', labelsize=8)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_13_fig1_time_stack.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved: 生活時間配分積み上げ棒グラフ")

# ------------------------------------------------------------------ #
# Fig2: 性別による生活時間差（全国計）
# ------------------------------------------------------------------ #
gender_cols = ['睡眠', '仕事', '家事', '育児', '通勤・通学',
               '買い物', 'テレビ・ラジオ・新聞・雑誌', '休養・くつろぎ',
               '趣味・娯楽', 'スポーツ', '学習・自己啓発・訓練(学業以外)', '交際・付き合い']
label_map = {
    '睡眠': '睡眠',
    '仕事': '仕事',
    '家事': '家事',
    '育児': '育児',
    '通勤・通学': '通勤・通学',
    '買い物': '買い物',
    'テレビ・ラジオ・新聞・雑誌': 'TV・ラジオ等',
    '休養・くつろぎ': '休養・くつろぎ',
    '趣味・娯楽': '趣味・娯楽',
    'スポーツ': 'スポーツ',
    '学習・自己啓発・訓練(学業以外)': '学習・自己啓発',
    '交際・付き合い': '交際・付き合い',
}
male_vals   = [nat_male[c]   for c in gender_cols]
female_vals = [nat_female[c] for c in gender_cols]
labels      = [label_map[c]  for c in gender_cols]

x = np.arange(len(labels))
width = 0.38

fig, ax = plt.subplots(figsize=(12, 6))
bars_m = ax.bar(x - width/2, male_vals,   width, label='男性', color='#4e9af1', alpha=0.85)
bars_f = ax.bar(x + width/2, female_vals, width, label='女性', color='#f06292', alpha=0.85)

# 差分アノテーション
for i, (m, f) in enumerate(zip(male_vals, female_vals)):
    diff = m - f
    sign = '+' if diff > 0 else ''
    color = '#1565C0' if diff > 0 else '#C62828'
    ax.text(x[i], max(m, f) + 4, f'{sign}{diff}分', ha='center', va='bottom',
            fontsize=7.5, color=color, fontweight='bold')

ax.set_xticks(x)
ax.set_xticklabels(labels, rotation=30, ha='right', fontsize=10)
ax.set_ylabel('1週間の平均時間（分）', fontsize=12)
ax.set_title('生活時間の男女差（全国計）\n（数値は男性 − 女性、青=男性が多い、赤=女性が多い）',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_13_fig2_gender.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved: 性別生活時間差")

# ------------------------------------------------------------------ #
# Fig3: 仕事時間 vs 趣味・娯楽時間 散布図（47都道府県）
# ------------------------------------------------------------------ #
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 左: 仕事時間 vs 趣味・娯楽時間
ax = axes[0]
for reg, grp in total.groupby('地域'):
    ax.scatter(grp['仕事'], grp['趣味・娯楽'],
               color=region_colors.get(reg, 'gray'), label=reg, s=55, alpha=0.85)
for _, row in total.iterrows():
    ax.annotate(row['都道府県'][:2], (row['仕事'], row['趣味・娯楽']),
                fontsize=6.5, alpha=0.7, xytext=(2, 2), textcoords='offset points')
x_vals = total['仕事'].values
y_vals = total['趣味・娯楽'].values
slope, intercept, r, p_val, _ = stats.linregress(x_vals, y_vals)
xr = np.linspace(x_vals.min(), x_vals.max(), 100)
ax.plot(xr, slope * xr + intercept, 'k--', linewidth=1.5,
        label=f'回帰直線\nr={r:.2f}, p={p_val:.3f}')
ax.set_xlabel('仕事時間（分/週）', fontsize=12)
ax.set_ylabel('趣味・娯楽時間（分/週）', fontsize=12)
ax.set_title('仕事時間 vs 趣味・娯楽時間\n（N=47都道府県, 総数）', fontsize=12, fontweight='bold')
ax.legend(fontsize=7.5, ncol=2)
ax.grid(alpha=0.3)

# 右: 家事時間 vs テレビ・ラジオ時間
ax = axes[1]
for reg, grp in total.groupby('地域'):
    ax.scatter(grp['家事'], grp['テレビ・ラジオ・新聞・雑誌'],
               color=region_colors.get(reg, 'gray'), label=reg, s=55, alpha=0.85)
for _, row in total.iterrows():
    ax.annotate(row['都道府県'][:2], (row['家事'], row['テレビ・ラジオ・新聞・雑誌']),
                fontsize=6.5, alpha=0.7, xytext=(2, 2), textcoords='offset points')
x_vals2 = total['家事'].values
y_vals2 = total['テレビ・ラジオ・新聞・雑誌'].values
slope2, intercept2, r2, p_val2, _ = stats.linregress(x_vals2, y_vals2)
xr2 = np.linspace(x_vals2.min(), x_vals2.max(), 100)
ax.plot(xr2, slope2 * xr2 + intercept2, 'k--', linewidth=1.5,
        label=f'回帰直線\nr={r2:.2f}, p={p_val2:.3f}')
ax.set_xlabel('家事時間（分/週）', fontsize=12)
ax.set_ylabel('TV・ラジオ等時間（分/週）', fontsize=12)
ax.set_title('家事時間 vs TV・ラジオ等時間\n（N=47都道府県, 総数）', fontsize=12, fontweight='bold')
ax.legend(fontsize=7.5, ncol=2)
ax.grid(alpha=0.3)

plt.suptitle('生活時間の関係性：都道府県間散布図', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_13_fig3_scatter.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved: 生活時間散布図")

# ------------------------------------------------------------------ #
# Fig4: 主成分分析 (PCA) による都道府県の生活時間パターン
# ------------------------------------------------------------------ #
pca_cols = ['睡眠', '仕事', '家事', '育児', '通勤・通学',
            'テレビ・ラジオ・新聞・雑誌', '休養・くつろぎ', '趣味・娯楽',
            'スポーツ', '買い物', '学習・自己啓発・訓練(学業以外)', '交際・付き合い']
pca_label_map = {
    '睡眠': '睡眠', '仕事': '仕事', '家事': '家事', '育児': '育児',
    '通勤・通学': '通勤', 'テレビ・ラジオ・新聞・雑誌': 'TV等',
    '休養・くつろぎ': '休養', '趣味・娯楽': '趣味', 'スポーツ': 'スポーツ',
    '買い物': '買い物', '学習・自己啓発・訓練(学業以外)': '学習', '交際・付き合い': '交際',
}

df_pca = total[['都道府県', '地域'] + pca_cols].dropna().copy()
X_raw = df_pca[pca_cols].values
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)

pca = PCA(n_components=2, random_state=42)
scores = pca.fit_transform(X_scaled)
loadings = pca.components_  # shape (2, n_features)
explained = pca.explained_variance_ratio_

df_pca['PC1'] = scores[:, 0]
df_pca['PC2'] = scores[:, 1]

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 左: スコアプロット
ax = axes[0]
for reg, grp in df_pca.groupby('地域'):
    ax.scatter(grp['PC1'], grp['PC2'],
               color=region_colors.get(reg, 'gray'), label=reg, s=60, alpha=0.85)
for _, row in df_pca.iterrows():
    ax.annotate(row['都道府県'][:2], (row['PC1'], row['PC2']),
                fontsize=6.5, alpha=0.75, xytext=(2, 2), textcoords='offset points')
ax.axhline(0, color='gray', linewidth=0.7, linestyle='--')
ax.axvline(0, color='gray', linewidth=0.7, linestyle='--')
ax.set_xlabel(f'第1主成分（寄与率 {explained[0]*100:.1f}%）', fontsize=11)
ax.set_ylabel(f'第2主成分（寄与率 {explained[1]*100:.1f}%）', fontsize=11)
ax.set_title('PCAスコアプロット\n（都道府県の生活時間パターン）', fontsize=12, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)

# 右: ローディングプロット（バイプロット補助）
ax = axes[1]
scale = 3.0
for i, col in enumerate(pca_cols):
    lx, ly = loadings[0, i] * scale, loadings[1, i] * scale
    ax.arrow(0, 0, lx, ly, head_width=0.06, head_length=0.05,
             fc='#1565C0', ec='#1565C0', alpha=0.8)
    ax.text(lx * 1.12, ly * 1.12, pca_label_map[col],
            fontsize=9, ha='center', va='center', color='#0D47A1', fontweight='bold')
ax.axhline(0, color='gray', linewidth=0.7, linestyle='--')
ax.axvline(0, color='gray', linewidth=0.7, linestyle='--')
ax.set_xlim(-scale * 0.7, scale * 0.7)
ax.set_ylim(-scale * 0.7, scale * 0.7)
ax.set_xlabel(f'第1主成分（寄与率 {explained[0]*100:.1f}%）', fontsize=11)
ax.set_ylabel(f'第2主成分（寄与率 {explained[1]*100:.1f}%）', fontsize=11)
ax.set_title('PCAローディングプロット\n（各活動が主成分に与える方向）', fontsize=12, fontweight='bold')
ax.grid(alpha=0.3)

plt.suptitle('主成分分析（PCA）による生活時間パターンの可視化', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_13_fig4_pca.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved: PCA分析")

# ------------------------------------------------------------------ #
# 結果サマリー
# ------------------------------------------------------------------ #
print("\n========== 分析結果サマリー ==========")
print(f"対象都道府県数: {len(total)}")
print(f"\n【全国計 生活時間 (分/週)】")
for c in ['睡眠', '仕事', '家事', '趣味・娯楽', 'テレビ・ラジオ・新聞・雑誌']:
    m, f = nat_male[c], nat_female[c]
    print(f"  {c:20s}: 男={m:4.0f}, 女={f:4.0f}, 差={m-f:+.0f}")

print(f"\n【仕事時間 最長上位5県】")
print(total.nlargest(5, '仕事')[['都道府県', '仕事', '趣味・娯楽', '家事']].to_string(index=False))

print(f"\n【仕事時間 最短下位5県】")
print(total.nsmallest(5, '仕事')[['都道府県', '仕事', '趣味・娯楽', '家事']].to_string(index=False))

x_corr = total['仕事'].values
y_corr = total['趣味・娯楽'].values
r_corr, p_corr = stats.pearsonr(x_corr, y_corr)
print(f"\n【仕事時間 vs 趣味・娯楽時間】r={r_corr:.3f}, p={p_corr:.3f}")

print(f"\n【PCA寄与率】PC1={explained[0]*100:.1f}%, PC2={explained[1]*100:.1f}%")
print(f"  累積寄与率: {(explained[0]+explained[1])*100:.1f}%")

print("\nAll figures saved to:", FIG_DIR)
print("Done!")