"""
2023年度 統計データ分析コンペティション
総務大臣賞（高校生の部）
論文テーマ：生活の形態と女性の社会進出

データ:
  SSDSE-D-2023.csv ... 社会生活基本調査（時間使用データ）都道府県別・男女別
  SSDSE-B-2026.csv ... 社会・人口統計体系（人口・労働・保育・消費）

分析手法:
  1. 相関分析  — 女性の時間使用 vs 合計特殊出生率 等のヒートマップ
  2. 重回帰分析 — TFR ～ 女性仕事時間 + 保育所密度 + 消費支出
  3. クラスター分析 — Wardリンケージで47都道府県を時間使用パターンでグループ化
  4. グループ別 TFR 箱ひげ図 — クラスター間差異の可視化
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-D-2023.csv
#       → data/raw/SSDSE-D-2023.csv に配置
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-D（社会・人口統計体系 都道府県の指標） の CSV をダウンロード）
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2023_H1_daijin.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.font_manager as fm
from scipy import stats
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from scipy.spatial.distance import pdist
from sklearn.preprocessing import StandardScaler
import statsmodels.api as sm

warnings.filterwarnings('ignore')

# ── パス設定 ────────────────────────────────────────────────────
_dir   = os.path.dirname(os.path.abspath(__file__))
DATA_D = os.path.join(_dir, '..', 'data', 'raw', 'SSDSE-D-2023.csv')
DATA_B = os.path.join(_dir, '..', 'data', 'raw', 'SSDSE-B-2026.csv')
FIG_DIR = os.path.join(_dir, '..', 'html', 'figures')
os.makedirs(FIG_DIR, exist_ok=True)

# ── フォント設定 ─────────────────────────────────────────────────
for fname in ['Hiragino Sans', 'Hiragino Kaku Gothic ProN', 'Yu Gothic',
              'MS Gothic', 'Noto Sans CJK JP']:
    if fname in [f.name for f in fm.fontManager.ttflist]:
        plt.rcParams['font.family'] = fname
        break
plt.rcParams['axes.unicode_minus'] = False

# ── 色パレット ────────────────────────────────────────────────────
C1, C2, C3, C4 = '#1565C0', '#E65100', '#2E7D32', '#6A1B9A'
CLUSTER_COLORS = ['#1565C0', '#E65100', '#2E7D32', '#6A1B9A']

# ════════════════════════════════════════════════════════════════
# 1. データ読み込み
# ════════════════════════════════════════════════════════════════

# ── SSDSE-D（社会生活基本調査 時間使用データ）────────────────────
df_d = pd.read_csv(DATA_D, encoding='cp932', header=1)
df_d_f = df_d[(df_d['男女の別'] == '2_女') & (df_d['地域コード'] != 'R00000')].copy()
# 時間使用列（分/日）
TIME_COLS = ['仕事', '家事', '育児', '睡眠', '介護・看護', '買い物', '休養・くつろぎ']
for c in TIME_COLS:
    df_d_f[c] = pd.to_numeric(df_d_f[c], errors='coerce')

print(f"SSDSE-D 女性データ: {len(df_d_f)} 都道府県")
print("時間使用列（分/日）先頭5行:")
print(df_d_f[['都道府県'] + TIME_COLS].head())

# ── SSDSE-B（社会・人口統計体系 2022年度）────────────────────────
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)) &
    (df_b['年度'] == 2022)
].copy()

# 必要列を数値化
NUM_COLS_B = ['合計特殊出生率', '15～64歳人口（女）', '総人口（女）',
              '保育所等数', '消費支出（二人以上の世帯）']
for c in NUM_COLS_B:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

print(f"\nSSDSE-B 2022年度: {len(df_b)} 都道府県")

# ── 派生変数の計算 ───────────────────────────────────────────────
# 15-64歳女性比率（労働年齢女性比率）%
df_b['女性労働年齢比率'] = df_b['15～64歳人口（女）'] / df_b['総人口（女）'] * 100

# 保育所密度（人口10万人あたり保育所数）
df_b['保育所密度'] = df_b['保育所等数'] / df_b['総人口（女）'] * 100000

# 消費支出（万円）
df_b['消費支出_万円'] = df_b['消費支出（二人以上の世帯）'] / 10000

# ── データ結合 ───────────────────────────────────────────────────
D_USE = ['地域コード', '都道府県'] + TIME_COLS
B_USE = ['地域コード', '合計特殊出生率', '女性労働年齢比率', '保育所密度', '消費支出_万円']
df = pd.merge(df_d_f[D_USE], df_b[B_USE], on='地域コード').dropna()

print(f"\n結合後データ: {len(df)} 都道府県（欠損除外後）")
print(df[['都道府県', '仕事', '家事', '育児', '合計特殊出生率', '保育所密度']].head())


# ════════════════════════════════════════════════════════════════
# 2. 分析
# ════════════════════════════════════════════════════════════════

# ── 2-1. 相関行列 ─────────────────────────────────────────────
CORR_COLS = ['仕事', '家事', '育児', '睡眠', '介護・看護', '買い物',
             '合計特殊出生率', '保育所密度', '消費支出_万円']
CORR_LABELS = ['仕事時間', '家事時間', '育児時間', '睡眠時間', '介護時間', '買物時間',
               '合計特殊出生率', '保育所密度', '消費支出']
corr_df = df[CORR_COLS].corr()

# ── 2-2. 重回帰: TFR ～ 仕事時間 + 保育所密度 + 消費支出 ──────
X_cols = ['仕事', '保育所密度', '消費支出_万円']
X_data = sm.add_constant(df[X_cols])
y_data = df['合計特殊出生率']
reg_model = sm.OLS(y_data, X_data).fit()

print("\n重回帰分析結果:")
print(reg_model.summary())

# ── 2-3. クラスター分析（Ward法）─────────────────────────────
cluster_vars = ['仕事', '家事', '育児', '睡眠', '介護・看護', '買い物', '休養・くつろぎ']
scaler = StandardScaler()
X_scaled = scaler.fit_transform(df[cluster_vars])
Z = linkage(X_scaled, method='ward')

# 4クラスター
n_clusters = 4
df['cluster'] = fcluster(Z, n_clusters, criterion='maxclust')

print("\nクラスター割り当て（グループ別都道府県数）:")
for g in sorted(df['cluster'].unique()):
    prefs = df[df['cluster'] == g]['都道府県'].tolist()
    print(f"  グループ {g} ({len(prefs)}都道府県): {', '.join(prefs)}")

print("\nグループ別 合計特殊出生率:")
print(df.groupby('cluster')['合計特殊出生率'].agg(['mean', 'std', 'count']).round(3))


# ════════════════════════════════════════════════════════════════
# 3. 図の作成
# ════════════════════════════════════════════════════════════════

def save_fig(fig, name):
    path = os.path.join(FIG_DIR, name)
    fig.savefig(path, dpi=150, bbox_inches='tight')
    plt.close(fig)
    print(f"  保存: {path}")

# ── 図1: 相関ヒートマップ ─────────────────────────────────────
fig, ax = plt.subplots(figsize=(9, 7))

n = len(CORR_LABELS)
mat = corr_df.values

im = ax.imshow(mat, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='相関係数 r')

ax.set_xticks(range(n))
ax.set_yticks(range(n))
ax.set_xticklabels(CORR_LABELS, rotation=45, ha='right', fontsize=10)
ax.set_yticklabels(CORR_LABELS, fontsize=10)

for i in range(n):
    for j in range(n):
        val = mat[i, j]
        color = 'white' if abs(val) > 0.5 else 'black'
        ax.text(j, i, f'{val:.2f}', ha='center', va='center',
                fontsize=8.5, color=color, fontweight='bold')

ax.set_title('図1：女性時間使用・TFR・保育所密度の相関ヒートマップ\n（都道府県別、2022〜2023年）',
             fontsize=13, fontweight='bold', pad=14)
fig.tight_layout()
save_fig(fig, '2023_H1_fig1_heatmap.png')

# ── 図2: 散布図（女性仕事時間 vs TFR）+ 回帰直線 ─────────────
fig, ax = plt.subplots(figsize=(8, 6))

x_val = df['仕事'].values
y_val = df['合計特殊出生率'].values
r, p = stats.pearsonr(x_val, y_val)

# 散布図（クラスター色）
for g in sorted(df['cluster'].unique()):
    mask = df['cluster'] == g
    ax.scatter(df.loc[mask, '仕事'], df.loc[mask, '合計特殊出生率'],
               s=70, color=CLUSTER_COLORS[g - 1], alpha=0.85, zorder=3,
               label=f'グループ {g}')

# 回帰直線
xfit = np.linspace(x_val.min(), x_val.max(), 100)
slope, intercept, *_ = stats.linregress(x_val, y_val)
ax.plot(xfit, slope * xfit + intercept, color='#C62828', linewidth=2,
        linestyle='--', label=f'回帰直線（r={r:.3f}, p={p:.3f}）')

# 都道府県ラベル（外れ値のみ表示）
mean_x, std_x = x_val.mean(), x_val.std()
mean_y, std_y = y_val.mean(), y_val.std()
for _, row in df.iterrows():
    if (abs(row['仕事'] - mean_x) > 1.3 * std_x or
            abs(row['合計特殊出生率'] - mean_y) > 1.3 * std_y):
        ax.annotate(row['都道府県'],
                    (row['仕事'], row['合計特殊出生率']),
                    textcoords='offset points', xytext=(4, 4),
                    fontsize=8.5, color='#333')

ax.set_xlabel('女性仕事時間（分/日）', fontsize=12)
ax.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax.set_title('図2：女性仕事時間と合計特殊出生率の関係\n（都道府県別、N=47）',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=10, loc='upper right')
ax.grid(True, linestyle='--', alpha=0.4)

sig_str = '有意（p<0.05）' if p < 0.05 else '非有意（p≥0.05）'
ax.text(0.02, 0.05, f'r = {r:.3f}  {sig_str}',
        transform=ax.transAxes, fontsize=10,
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#FFF9C4', alpha=0.9))

fig.tight_layout()
save_fig(fig, '2023_H1_fig2_scatter.png')

# ── 図3: Ward デンドログラム ──────────────────────────────────
fig, ax = plt.subplots(figsize=(14, 5.5))

pref_labels = df['都道府県'].tolist()
dend = dendrogram(Z, labels=pref_labels, ax=ax,
                  color_threshold=Z[-(n_clusters - 1), 2],
                  leaf_rotation=90, leaf_font_size=9.5)

ax.axhline(y=Z[-(n_clusters - 1), 2], color='#C62828',
           linestyle='--', linewidth=1.4, alpha=0.8, label='4クラスター分割点')
ax.set_title('図3：女性時間使用パターンによる47都道府県のWardクラスタリング\n（標準化後7変数）',
             fontsize=13, fontweight='bold')
ax.set_xlabel('都道府県', fontsize=11)
ax.set_ylabel('Ward距離', fontsize=11)
ax.legend(fontsize=10)
ax.grid(True, axis='y', linestyle='--', alpha=0.3)

fig.tight_layout()
save_fig(fig, '2023_H1_fig3_dendrogram.png')

# ── 図4: グループ別TFR箱ひげ図 ───────────────────────────────
fig, ax = plt.subplots(figsize=(8, 6))

groups = sorted(df['cluster'].unique())
data_by_group = [df[df['cluster'] == g]['合計特殊出生率'].values for g in groups]
labels_g = [f'グループ {g}\n(N={len(d)})' for g, d in zip(groups, data_by_group)]

bp = ax.boxplot(data_by_group, patch_artist=True, notch=False,
                medianprops=dict(color='white', linewidth=2.5),
                whiskerprops=dict(linewidth=1.4),
                capprops=dict(linewidth=1.4),
                flierprops=dict(marker='o', markersize=5, alpha=0.6))

for patch, color in zip(bp['boxes'], CLUSTER_COLORS):
    patch.set_facecolor(color)
    patch.set_alpha(0.8)

# グループ別平均の時間使用特徴をラベルとして付加
group_char = {
    1: '（仕事時間\n多め）',
    2: '（家事・育児\n多め）',
    3: '（バランス型）',
    4: '（睡眠・休養\n多め）'
}
# 実データに基づく特徴ラベル（グループ平均で再計算）
g_means = df.groupby('cluster')[cluster_vars].mean()
for i, g in enumerate(groups):
    row = g_means.loc[g]
    top_var = row.idxmax()
    var_map = {'仕事': '仕事時間↑', '家事': '家事↑', '育児': '育児↑',
               '睡眠': '睡眠↑', '介護・看護': '介護↑',
               '買い物': '買物↑', '休養・くつろぎ': '休養↑'}
    char_str = var_map.get(top_var, top_var)
    ax.text(i + 1, data_by_group[i].min() - 0.04,
            f'({char_str})', ha='center', va='top', fontsize=8.5, color='#555')

ax.set_xticklabels(labels_g, fontsize=11)
ax.set_xlabel('クラスターグループ（女性時間使用パターン）', fontsize=12)
ax.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax.set_title('図4：女性生活時間パターンのクラスター別 合計特殊出生率の分布\n（都道府県別、N=47）',
             fontsize=13, fontweight='bold')
ax.grid(True, axis='y', linestyle='--', alpha=0.4)

# Kruskal-Wallis 検定
stat_kw, p_kw = stats.kruskal(*data_by_group)
sig_str_kw = f'Kruskal-Wallis H={stat_kw:.2f}, p={p_kw:.3f}'
ax.text(0.98, 0.97, sig_str_kw, transform=ax.transAxes,
        ha='right', va='top', fontsize=9,
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#E3F2FD', alpha=0.9))

fig.tight_layout()
save_fig(fig, '2023_H1_fig4_boxplot.png')

print("\n全図保存完了。")

# ════════════════════════════════════════════════════════════════
# 4. 回帰結果サマリー出力
# ════════════════════════════════════════════════════════════════
print("\n" + "="*60)
print("重回帰分析サマリー（TFR ～ 女性仕事時間 + 保育所密度 + 消費支出）")
print("="*60)
print(f"  R² = {reg_model.rsquared:.4f}")
print(f"  調整済R² = {reg_model.rsquared_adj:.4f}")
print(f"  F統計量 = {reg_model.fvalue:.3f}, p = {reg_model.f_pvalue:.4f}")
print("\n  係数:")
for name, coef, pval in zip(reg_model.params.index, reg_model.params, reg_model.pvalues):
    sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else ''
    print(f"    {name:20s}: {coef:+.6f}  p={pval:.4f} {sig}")

print("\n" + "="*60)
print("相関係数（女性仕事時間 vs 合計特殊出生率）")
print("="*60)
r_work_tfr, p_work_tfr = stats.pearsonr(df['仕事'], df['合計特殊出生率'])
print(f"  r = {r_work_tfr:.4f}, p = {p_work_tfr:.4f}")

r_house_tfr, p_house_tfr = stats.pearsonr(df['家事'], df['合計特殊出生率'])
print(f"\n相関係数（家事時間 vs 合計特殊出生率）")
print(f"  r = {r_house_tfr:.4f}, p = {p_house_tfr:.4f}")