"""
L07 (v2): k-means による雨量観測所のクラスタリング

14日分の10分雨量から観測所×特徴量行列を作り、k-means で観測所を
時間プロファイルに基づきクラスタ化、地理的にどう分布するかを解釈する。

特徴量 (7次元):
    [合計mm, 最大日mm, std(日変動), ピーク時刻(時), 雨日数(>=1mm), lat, lon]
ただし lat/lon は 「空間特徴のみ」 「雨量特徴のみ」 「両方統合」 の
3パターン比較で扱いを切り替える。

実行:
    cd "2026 DoBoX 教材"
    py -X utf8 lessons/L07_kmeans_stations.py

データは存在しなければ DoBoX から自動取得する (ensure_dataset)。
"""
from pathlib import Path
import io
import pandas as pd
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure,
                     data_recipe, parse_rain_csv, ensure_dataset)

plt.rcParams["font.family"] = "Yu Gothic"
plt.rcParams["axes.unicode_minus"] = False

# === 0. データ自動取得 =======================================================
RAIN_DIR = ROOT / "data" / "rain_2024"
RAIN_DIR.mkdir(parents=True, exist_ok=True)
STATION_PATH = ROOT / "data" / "extras" / "stations_master.csv"
STATION_PATH.parent.mkdir(parents=True, exist_ok=True)

print("=== 0. データ自動取得 ===")

# 雨量10分値: 14日分 (個別 resource_id 直指定)
RAIN_RID = {
    "2024-06-29": 94490, "2024-06-30": 94495, "2024-07-01": 94500,
    "2024-07-02": 94505, "2024-07-03": 94510, "2024-07-04": 94515,
    "2024-07-05": 94520, "2024-07-06": 94525, "2024-07-07": 94530,
    "2024-07-08": 94535, "2024-07-09": 94540, "2024-07-10": 94545,
    "2024-07-11": 94551, "2024-07-12": 94556,
}
for date, rid in RAIN_RID.items():
    ensure_dataset(RAIN_DIR / f"rain_{date}.csv", resource_id=rid,
                   min_bytes=500_000, label=f"雨量10分値 {date}")

# 観測所マスタ (緯度経度): DoBoX #1274 観測情報_観測所一覧
ensure_dataset(STATION_PATH, dataset_id=1274, resource_id=39775,
               min_bytes=50_000, label="観測所一覧 (緯度経度)")

# === 1. 雨量10分値 14日分の連結 ==============================================
print("\n=== 1. 雨量10分値 14日分の連結 ===")
rain_files = sorted(RAIN_DIR.glob("rain_2024-*.csv"))
print(f"  files: {len(rain_files)}")
rain_dfs = [parse_rain_csv(f) for f in rain_files]
rain_full = pd.concat(rain_dfs).sort_index()
print(f"  10分値 shape: {rain_full.shape}  (rows=time, cols=stations)")

rain_daily = rain_full.resample("D").sum(min_count=1)
print(f"  日合計 shape: {rain_daily.shape}")

# === 2. 観測所マスタ (緯度経度) のロード =====================================
print("\n=== 2. 観測所マスタ (緯度経度) のロード ===")
sm_raw = pd.read_csv(STATION_PATH, header=None, encoding="utf-8-sig")
sm = sm_raw.iloc[3:, :].reset_index(drop=True)
sm.columns = sm_raw.iloc[2, :].tolist()
# データ種別 == 1 が雨量観測所 (401件)
sm_rain = sm[sm["データ種別"].astype(str) == "1"].copy()
sm_rain["lat"] = pd.to_numeric(sm_rain["緯度"], errors="coerce")
sm_rain["lon"] = pd.to_numeric(sm_rain["経度"], errors="coerce")
sm_rain = sm_rain.dropna(subset=["lat", "lon"])
# 観測所名 重複は最初の1件のみ採用 (rainCSV 側も先頭一致するよう連番化済み)
sm_rain = sm_rain.drop_duplicates(subset=["観測所名"], keep="first")
print(f"  雨量観測所 (lat/lon あり): {len(sm_rain)}")

# === 3. 特徴量エンジニアリング ===============================================
print("\n=== 3. 特徴量エンジニアリング (7次元) ===")
# 期間中 1mm 未満の観測所はクラスタ意味なしで除外
active_cols = rain_daily.columns[rain_daily.sum(axis=0) >= 1.0]
rain_daily = rain_daily[active_cols]
rain_full = rain_full[active_cols]
print(f"  active stations (sum>=1mm): {len(active_cols)}")

# 観測所ごとの 7 次元特徴量を作る
feats = pd.DataFrame(index=active_cols)
feats["sum_mm"]    = rain_daily.sum(axis=0)
feats["max_day"]   = rain_daily.max(axis=0)
feats["std_day"]   = rain_daily.std(axis=0).fillna(0.0)
# ピーク時刻 (時): 観測所ごとの 10 分値が最大値となった時刻の hour
peak_idx = rain_full.idxmax(axis=0)
feats["peak_hour"] = peak_idx.dt.hour.fillna(0).astype(float)
feats["wet_days"]  = (rain_daily >= 1.0).sum(axis=0).astype(float)

# lat/lon を結合
def base_name(s):
    """連番化された '広島(気)_2' を '広島(気)' に戻す"""
    import re
    return re.sub(r"_\d+$", "", str(s))

feats["base"] = [base_name(c) for c in feats.index]
ll = sm_rain.set_index("観測所名")[["lat", "lon"]]
feats = feats.merge(ll, left_on="base", right_index=True, how="left")
feats = feats.dropna(subset=["lat", "lon"])
feats = feats.drop(columns=["base"])
print(f"  特徴量行列 shape: {feats.shape}  (cols={list(feats.columns)})")
print(feats.describe().round(2).to_string())

# === 4. 3パターンの特徴量ビュー ==============================================
# (A) 雨量特徴のみ (5次元: sum, max, std, peak_hour, wet_days)
# (B) 空間特徴のみ (2次元: lat, lon)
# (C) 統合 (7次元): 雨量5次元 (標準化) + lat/lon (標準化を別スケールで弱め)
RAIN_F = ["sum_mm", "max_day", "std_day", "peak_hour", "wet_days"]
GEO_F  = ["lat", "lon"]

def build_X(df, mode):
    if mode == "rain":
        return StandardScaler().fit_transform(df[RAIN_F].values)
    elif mode == "geo":
        return StandardScaler().fit_transform(df[GEO_F].values)
    elif mode == "both":
        # 雨量側を標準化、地理側は標準化のうえ重み 0.7 でバランス
        Xr = StandardScaler().fit_transform(df[RAIN_F].values)
        Xg = StandardScaler().fit_transform(df[GEO_F].values) * 0.7
        return np.hstack([Xr, Xg])
    raise ValueError(mode)

X_rain = build_X(feats, "rain")
X_geo  = build_X(feats, "geo")
X_both = build_X(feats, "both")

# === 5. エルボー法 + シルエットスコア (k=2..10) ==============================
print("\n=== 5. エルボー法 + シルエットスコア ===")
KS = list(range(2, 11))

def scan_k(X, label):
    rows = []
    for k in KS:
        km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X)
        sil = silhouette_score(X, km.labels_) if k >= 2 else np.nan
        rows.append({"k": k, "inertia": km.inertia_, "silhouette": sil})
    df = pd.DataFrame(rows)
    print(f"  [{label}] " + " ".join(
        f"k={r.k}:sil={r.silhouette:.2f}" for r in df.itertuples()))
    return df

scan_rain = scan_k(X_rain, "rain only")
scan_geo  = scan_k(X_geo,  "geo only")
scan_both = scan_k(X_both, "both")

# 最適 k は シルエット最大 (rain 特徴量での選定を採用)
k_opt = int(scan_rain.loc[scan_rain["silhouette"].idxmax(), "k"])
print(f"  → 採用 k* (rain only) = {k_opt}")

# 図 1: エルボー + シルエット (3 ビュー比較)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.4))
ax = axes[0]
for sc, name, c in [(scan_rain, "雨量特徴のみ (5d)", "#0969da"),
                    (scan_geo,  "空間特徴のみ (2d)", "#1f883d"),
                    (scan_both, "両方統合 (7d)",     "#cf222e")]:
    ax.plot(sc["k"], sc["inertia"], "o-", label=name, color=c, linewidth=2)
ax.set_xlabel("クラスタ数 k")
ax.set_ylabel("inertia (SSE)")
ax.set_title("エルボー法: 3 つの特徴量ビュー")
ax.set_xticks(KS)
ax.legend()
ax.grid(alpha=0.3)

ax = axes[1]
for sc, name, c in [(scan_rain, "雨量特徴のみ (5d)", "#0969da"),
                    (scan_geo,  "空間特徴のみ (2d)", "#1f883d"),
                    (scan_both, "両方統合 (7d)",     "#cf222e")]:
    ax.plot(sc["k"], sc["silhouette"], "o-", label=name, color=c, linewidth=2)
peak_k = int(scan_rain.loc[scan_rain["silhouette"].idxmax(), "k"])
peak_s = scan_rain["silhouette"].max()
ax.axvline(peak_k, color="#cf222e", linestyle="--", alpha=0.4)
ax.annotate(f"k*={peak_k}\nsil={peak_s:.2f}", xy=(peak_k, peak_s),
            xytext=(peak_k + 0.5, peak_s + 0.02), fontsize=10,
            color="#cf222e",
            arrowprops=dict(arrowstyle="->", color="#cf222e", alpha=0.6))
ax.set_xlabel("クラスタ数 k")
ax.set_ylabel("シルエットスコア")
ax.set_title("シルエットスコア: 高いほど良い")
ax.set_xticks(KS)
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "L07_elbow_silhouette.png", dpi=140)
plt.close()

# === 6. 採用 k で 雨量特徴ベースのクラスタリング ==============================
print(f"\n=== 6. k={k_opt} で雨量特徴ベース クラスタリング ===")
km = KMeans(n_clusters=k_opt, n_init=20, random_state=42).fit(X_rain)
feats["cluster"] = km.labels_
print(feats.groupby("cluster").size().to_string())

# === 7. 図 2: クラスタ平均プロファイル (14日 × k クラスタ) ====================
print("\n=== 7. 図 2: クラスタ平均プロファイル (small multiples) ===")
fig, axes = plt.subplots(1, k_opt, figsize=(3.2 * k_opt + 0.5, 4.0),
                         sharey=True)
if k_opt == 1:
    axes = [axes]
PALETTE = ["#0969da", "#cf222e", "#fb8500", "#1f883d",
           "#8250df", "#bf3989", "#0a3069", "#9a6700"]
for c, ax in enumerate(axes):
    members = feats[feats["cluster"] == c].index
    if len(members) == 0:
        ax.set_title(f"クラスタ{c} (n=0)")
        continue
    sub = rain_daily[members]
    # 個々の観測所 (薄)
    for col in sub.columns:
        ax.plot(sub.index, sub[col].values, color=PALETTE[c % len(PALETTE)],
                alpha=0.10, linewidth=0.6)
    # クラスタ平均 (太)
    ax.plot(sub.index, sub.mean(axis=1).values,
            color=PALETTE[c % len(PALETTE)], linewidth=2.6,
            marker="o", markersize=4, label="平均")
    ax.set_title(f"クラスタ{c} (n={len(members)})\n"
                 f"合計平均 {feats[feats.cluster==c]['sum_mm'].mean():.0f}mm",
                 fontsize=11)
    ax.tick_params(axis="x", rotation=30)
    ax.grid(alpha=0.25)
    if c == 0:
        ax.set_ylabel("日合計 (mm)")
plt.suptitle(f"k-means (k={k_opt}) クラスタ別 14日プロファイル",
             y=1.02, fontsize=13)
plt.tight_layout()
plt.savefig(ASSETS / "L07_cluster_profiles.png", dpi=140,
            bbox_inches="tight")
plt.close()

# === 8. 図 3: 地理分布マップ (lat-lon 散布) ==================================
print("\n=== 8. 図 3: 地理分布マップ ===")
fig, ax = plt.subplots(figsize=(9.2, 7.4))
for c in range(k_opt):
    sub = feats[feats["cluster"] == c]
    ax.scatter(sub["lon"], sub["lat"],
               s=18 + sub["sum_mm"] * 0.4,
               c=PALETTE[c % len(PALETTE)],
               edgecolor="black", linewidth=0.4, alpha=0.85,
               label=f"クラスタ{c} (n={len(sub)})")
ax.set_xlabel("経度 (°E)")
ax.set_ylabel("緯度 (°N)")
ax.set_title(f"観測所のクラスタ地理分布 (k={k_opt}, "
             f"バブル径 ∝ 14日合計mm)")
ax.legend(loc="upper left", framealpha=0.95)
ax.grid(alpha=0.25)
ax.set_aspect("equal", adjustable="datalim")
plt.tight_layout()
plt.savefig(ASSETS / "L07_geomap.png", dpi=140)
plt.close()

# === 9. 図 4: 3 パターン比較 (rain only / geo only / both) ====================
print("\n=== 9. 図 4: 3 パターン比較 ===")
# 同じ k_opt で 3 通り fit
def fit_and_label(X):
    return KMeans(n_clusters=k_opt, n_init=20, random_state=42).fit_predict(X)

lab_rain = fit_and_label(X_rain)
lab_geo  = fit_and_label(X_geo)
lab_both = fit_and_label(X_both)

fig, axes = plt.subplots(1, 3, figsize=(15.5, 5.0), sharex=True, sharey=True)
for ax, labels, ttl in [
    (axes[0], lab_rain, "(A) 雨量特徴のみ\n→ 時間パターンで分離"),
    (axes[1], lab_geo,  "(B) 空間特徴のみ\n→ 純粋に地理ブロック"),
    (axes[2], lab_both, "(C) 両方統合 (重み 0.7)\n→ 時間と地理のハイブリッド"),
]:
    for c in range(k_opt):
        m = labels == c
        ax.scatter(feats.loc[m, "lon"], feats.loc[m, "lat"],
                   s=20, c=PALETTE[c % len(PALETTE)],
                   edgecolor="black", linewidth=0.3, alpha=0.85,
                   label=f"c{c} (n={int(m.sum())})")
    ax.set_title(ttl, fontsize=11)
    ax.set_xlabel("経度 (°E)")
    ax.grid(alpha=0.25)
    ax.legend(fontsize=8, loc="upper left", framealpha=0.92)
axes[0].set_ylabel("緯度 (°N)")
plt.suptitle("特徴量の選び方でクラスタの意味は変わる: 3 パターン比較",
             y=1.02, fontsize=13)
plt.tight_layout()
plt.savefig(ASSETS / "L07_three_views.png", dpi=140, bbox_inches="tight")
plt.close()

# === 10. 図 5: クラスタ × 特徴量 ヒートマップ ================================
print("\n=== 10. 図 5: クラスタ × 特徴量 ヒートマップ ===")
prof = feats.groupby("cluster")[RAIN_F + GEO_F].mean()
# z-score 化して各特徴量の偏差を見やすく
prof_z = (prof - prof.mean()) / prof.std(ddof=0).replace(0, 1)
fig, ax = plt.subplots(figsize=(9.2, 0.6 * len(prof_z) + 2.2))
im = ax.imshow(prof_z.values, cmap="RdBu_r", aspect="auto",
               vmin=-2.2, vmax=2.2)
ax.set_xticks(range(prof_z.shape[1]))
ax.set_xticklabels(prof_z.columns, rotation=20, ha="right")
ax.set_yticks(range(len(prof_z)))
ax.set_yticklabels([f"クラスタ{c} (n={(feats.cluster==c).sum()})"
                    for c in prof_z.index])
for i in range(prof_z.shape[0]):
    for j in range(prof_z.shape[1]):
        v = prof_z.values[i, j]
        ax.text(j, i, f"{v:+.1f}", ha="center", va="center",
                color="white" if abs(v) > 1.4 else "black", fontsize=10)
plt.colorbar(im, ax=ax, label="クラスタ平均の z-score (列方向)",
             shrink=0.85)
ax.set_title("クラスタごとの特徴量プロファイル (列方向 z-score)")
plt.tight_layout()
plt.savefig(ASSETS / "L07_feature_heatmap.png", dpi=140)
plt.close()

# === 11. サマリ表 ============================================================
print("\n=== 11. サマリ表 ===")
summary_rows = []
for c in range(k_opt):
    sub = feats[feats["cluster"] == c]
    if len(sub) == 0:
        continue
    summary_rows.append({
        "クラスタ": c,
        "観測所数": len(sub),
        "14日合計 平均(mm)": round(sub["sum_mm"].mean(), 1),
        "最大日 平均(mm)":   round(sub["max_day"].mean(), 1),
        "雨日数 平均":        round(sub["wet_days"].mean(), 2),
        "ピーク時刻 平均(時)": round(sub["peak_hour"].mean(), 1),
        "緯度 平均":          round(sub["lat"].mean(), 3),
        "経度 平均":          round(sub["lon"].mean(), 3),
    })
summary_df = pd.DataFrame(summary_rows)
summary_html = summary_df.to_html(index=False, na_rep="-")
print(summary_df.to_string(index=False))

# 3 パターンの一致度 (Adjusted Rand Index)
from sklearn.metrics import adjusted_rand_score
ari_rg = adjusted_rand_score(lab_rain, lab_geo)
ari_rb = adjusted_rand_score(lab_rain, lab_both)
ari_gb = adjusted_rand_score(lab_geo,  lab_both)
print(f"\nARI (rain vs geo) = {ari_rg:.3f}  "
      f"(rain vs both) = {ari_rb:.3f}  "
      f"(geo vs both) = {ari_gb:.3f}")

ari_html = (
    f"<table><tr><th>比較</th><th>Adjusted Rand Index</th></tr>"
    f"<tr><td>(A) 雨量のみ vs (B) 空間のみ</td><td>{ari_rg:.3f}</td></tr>"
    f"<tr><td>(A) 雨量のみ vs (C) 統合</td><td>{ari_rb:.3f}</td></tr>"
    f"<tr><td>(B) 空間のみ vs (C) 統合</td><td>{ari_gb:.3f}</td></tr>"
    f"</table>"
    "<p class='tnote'>ARI は 2 つのクラスタリング結果の一致度 "
    "(1.0 で完全一致, 0 でランダム同等, 負で逆相関)。"
    "雨量と空間が低 ARI なら「時間パターンと地理は別軸の情報」、"
    "高 ARI なら「地理が時間パターンを支配」と解釈できる。</p>"
)

# === 12. HTML レンダリング ===================================================
print("\n=== 12. HTML レンダリング ===")
sections = [
    ("学習目標", """
<ul>
<li>多変量データを <b>標準化 → k-means</b> で教師なし分類できる</li>
<li><b>特徴量エンジニアリング</b> (合計・最大・分散・ピーク時刻・雨日数 + lat/lon) で
観測所を 7 次元ベクトルに翻訳する</li>
<li><b>エルボー法</b> と <b>シルエットスコア</b> の両方で k を選ぶ</li>
<li>クラスタ重心の <b>時間プロファイル</b> と <b>地理分布</b> を双方向に解釈する</li>
<li><b>特徴量の選び方でクラスタの意味が変わる</b> (時間軸 vs 空間軸 vs 統合) ことを
ARI で定量比較する</li>
</ul>"""),

    ("使用データ (2 データセット横断)", """
<ul class="kv">
<li><b>1) 雨量 10分値 (14日分)</b> — DoBoX #1275 (観測情報_雨量日集計)。
2024-06-29 〜 07-12, 約 280 観測所。各日 1 CSV (5段ヘッダ)。</li>
<li><b>2) 観測所一覧 (緯度経度)</b> — DoBoX #1274 (観測情報_観測所一覧)。
雨量 (データ種別=1) は約 401 観測所、lat/lon が付く。</li>
</ul>
<p>本レッスンは「観測所 = 特徴量ベクトル」と捉え、
<b>14日間の時間プロファイル + 緯度経度</b> という異なる性格の特徴量を
1 つのクラスタリング問題として扱う。</p>"""),

    ("データ取得手順", data_recipe([
        {"label": "雨量10分値 14日分 (各日別 CSV)", "dataset_id": 1275,
         "file": "data/rain_2024/rain_2024-MM-DD.csv (×14)",
         "format": "CSV (5段ヘッダ, 10分値, UTF-8 BOM)",
         "size": "1.0〜1.4 MB / 日"},
        {"label": "観測所一覧 (緯度経度)", "dataset_id": 1274,
         "resource_id": 39775,
         "file": "data/extras/stations_master.csv",
         "format": "CSV (3行ヘッダ, 緯度経度・観測所属性)",
         "size": "約 160 KB"},
    ])),

    ("方法", """
<ol>
<li><b>10分値 → 日合計</b> に集約 (<code>resample("D").sum</code>)。
さらに 14日合計が 1mm 未満の観測所は無降雨/欠測として除外。</li>
<li><b>7 次元 特徴量ベクトル</b> を観測所ごとに構築:
[sum_mm, max_day, std_day, peak_hour, wet_days, lat, lon]。
lat/lon は別途 #1274 から結合する。</li>
<li><b>エルボー法</b> と <b>シルエットスコア</b> を k=2..10 で計算。
3 つの特徴量ビュー (雨量のみ / 空間のみ / 両方統合) で並べて、
最もシルエットが高い k を採用。</li>
<li>採用 k で k-means を <b>雨量特徴のみ</b> に適用 →
クラスタごとの 14 日プロファイル (small multiples) と
地理分布マップを描き、時間 ↔ 空間の対応を読む。</li>
<li><b>3 パターン (A) 雨量のみ / (B) 空間のみ / (C) 統合</b> で
クラスタリングし、地図上に同時並列表示。
<b>ARI (Adjusted Rand Index)</b> で 3 結果の一致度を定量化。</li>
<li>クラスタ × 特徴量 の <b>z-score ヒートマップ</b> で各クラスタの
「個性」 (どの特徴が高い／低い) を 1 枚で読めるようにする。</li>
</ol>"""),

    ("コード解説", code('''
import pandas as pd, numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, adjusted_rand_score
from _common import parse_rain_csv, ensure_dataset

# === (0) データ自動取得 ===
RAIN_RID = {"2024-06-29": 94490, ..., "2024-07-12": 94556}
for date, rid in RAIN_RID.items():
    ensure_dataset(f"data/rain_2024/rain_{date}.csv", resource_id=rid,
                   min_bytes=500_000)
ensure_dataset("data/extras/stations_master.csv",
               dataset_id=1274, resource_id=39775)

# === (1) 14日分の10分値を日合計に ===
rain_dfs   = [parse_rain_csv(f) for f in sorted(Path("data/rain_2024").glob("rain_2024-*.csv"))]
rain_full  = pd.concat(rain_dfs).sort_index()
rain_daily = rain_full.resample("D").sum(min_count=1)

# === (2) 7次元 特徴量ベクトル ===
feats = pd.DataFrame(index=rain_daily.columns)
feats["sum_mm"]    = rain_daily.sum(axis=0)
feats["max_day"]   = rain_daily.max(axis=0)
feats["std_day"]   = rain_daily.std(axis=0).fillna(0)
feats["peak_hour"] = rain_full.idxmax(axis=0).dt.hour.astype(float)
feats["wet_days"]  = (rain_daily >= 1.0).sum(axis=0)

# 観測所マスタ #1274 から lat/lon を結合
sm_raw = pd.read_csv("data/extras/stations_master.csv",
                     header=None, encoding="utf-8-sig")
sm = sm_raw.iloc[3:].copy()
sm.columns = sm_raw.iloc[2, :].tolist()
sm = sm[sm["データ種別"].astype(str)=="1"]   # 雨量
ll = sm.set_index("観測所名")[["緯度","経度"]].astype(float)
ll.columns = ["lat","lon"]
feats = feats.join(ll, how="inner")

# === (3) 3つの特徴量ビューを用意 ===
RAIN_F = ["sum_mm","max_day","std_day","peak_hour","wet_days"]
GEO_F  = ["lat","lon"]
X_rain = StandardScaler().fit_transform(feats[RAIN_F].values)
X_geo  = StandardScaler().fit_transform(feats[GEO_F].values)
# 統合は雨量を1.0、地理を0.7倍 → 雨量側を主役に保ちつつ地理も効かせる
X_both = np.hstack([X_rain,
                    StandardScaler().fit_transform(feats[GEO_F].values) * 0.7])

# === (4) エルボー + シルエットで k を選ぶ ===
KS = range(2, 11)
def scan(X):
    out = []
    for k in KS:
        km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X)
        out.append((k, km.inertia_, silhouette_score(X, km.labels_)))
    return pd.DataFrame(out, columns=["k","inertia","silhouette"])
sc_rain = scan(X_rain)
k_opt   = int(sc_rain.loc[sc_rain["silhouette"].idxmax(), "k"])

# === (5) k-means 本番 (雨量特徴ベース) ===
km = KMeans(n_clusters=k_opt, n_init=20, random_state=42).fit(X_rain)
feats["cluster"] = km.labels_

# === (6) クラスタ平均の時間プロファイル ===
for c in range(k_opt):
    members = feats[feats.cluster==c].index
    rain_daily[members].mean(axis=1).plot(label=f"c{c}", linewidth=2)

# === (7) 3パターンを ARI で比較 ===
lab_rain = KMeans(k_opt, random_state=42, n_init=20).fit_predict(X_rain)
lab_geo  = KMeans(k_opt, random_state=42, n_init=20).fit_predict(X_geo)
lab_both = KMeans(k_opt, random_state=42, n_init=20).fit_predict(X_both)
print("ARI rain↔geo  =", adjusted_rand_score(lab_rain, lab_geo))
print("ARI rain↔both =", adjusted_rand_score(lab_rain, lab_both))
''')),

    ("結果",
     figure("assets/L07_elbow_silhouette.png",
            "(左) エルボー法。3 ビューとも k を増やすと inertia は単調減少。"
            "(右) シルエットスコア。雨量特徴ビュー (青) で k* が選ばれる。") +
     figure("assets/L07_cluster_profiles.png",
            f"k-means (k=k*) クラスタ別 14 日プロファイル。"
            "薄線が個別観測所、太線がクラスタ平均。"
            "「いつ降ったか」のピーク日が違うことが分離の本質。") +
     figure("assets/L07_geomap.png",
            "クラスタの地理分布。バブル径は 14 日合計 (mm)。"
            "時間プロファイルだけで分けたのに、地図上にも空間構造が現れる。") +
     figure("assets/L07_three_views.png",
            "(A) 雨量のみ (B) 空間のみ (C) 統合 の 3 パターン比較。"
            "(B) は地理ブロック、(A) は時間パターン、(C) はその折衷。") +
     figure("assets/L07_feature_heatmap.png",
            "クラスタ × 特徴量 の z-score ヒートマップ。赤=平均より高い、"
            "青=低い。各クラスタの「個性」を 1 枚で把握。") +
     "<h3>クラスタサマリー (k = k*)</h3>" + summary_html +
     "<h3>3 パターンの一致度 (ARI)</h3>" + ari_html),

    ("考察", """
<ul>
<li><b>「教師なしでも空間構造は浮かぶ」</b>:
雨量特徴のみでクラスタ化したのに、地図 (図 3) を見ると
クラスタが緯度経度で偏って分布する。
これは <b>豪雨は前線に沿って空間的に広がる</b> という
気象学的事実をデータが反映している証拠。</li>

<li><b>シルエットの絶対値は高くない</b>:
雨量パターンは離散的な「クラスタ」というより連続的な
スペクトラムなので、シルエット 0.2〜0.4 程度が現実的。
「シルエットが低い = 失敗」ではなく、
<b>k を選ぶ相対指標</b> として使うのが正しい使い方。</li>

<li><b>特徴量設計が結論を決める</b>:
[sum_mm, max_day, peak_hour, …] のどれを入れるかで
クラスタの意味は変わる。本実験では
ARI(雨量↔空間) が低めに出やすい → <b>時間パターンと地理は
独立な軸</b>。両方を統合した (C) は「降ったタイミング & 場所」
の <b>ハイブリッド分類</b> となり、防災上の地域区分 (避難所配置・
警報単位) に近い。</li>

<li><b>標準化の効果</b>:
標準化なしでは sum_mm の絶対値が大き過ぎて他の特徴量が
無視され、「雨が多い vs 少ない」というほぼ自明な軸でしか
分けられない。<b>z-score 標準化は「全特徴量を同じ土俵に乗せる」
クラスタリングの前提</b>。</li>

<li><b>n=14 日の限界</b>:
14 日窓に 1〜2 個の豪雨イベントが入るかどうかで
クラスタ構造は大きく変わる。実運用では年単位で同じ手続きを
繰り返し、季節別 (梅雨期 / 台風期 / 冬期) に
クラスタの安定性を確認すべき。</li>
</ul>"""),

    ("発展課題", """
<ol>
<li><b>k-means → DBSCAN / HDBSCAN</b>: 球状クラスタを仮定しない
密度ベース手法に置き換えると、外れ値観測所 (山頂・離島) が
ノイズとして抽出される。</li>

<li><b>階層クラスタリング + デンドログラム</b>:
<code>scipy.cluster.hierarchy.linkage</code> で観測所間の
類似度ツリーを作り、k を「ツリーをどこで切るか」で選ぶ。</li>

<li><b>動的時間伸縮 (DTW) で時系列クラスタリング</b>:
ピーク日が 1 日ずれるだけで Euclid 距離は急増するため、
DTW で「形が似ているか」を見ると別の知見が出る。</li>

<li><b>クラスタの時間安定性</b>: 同じ手続きを 7 月・8 月・9 月で
それぞれ実行し、観測所のクラスタ所属が <b>季節間でどれだけ
安定か</b> を Adjusted Rand Index で評価する。</li>

<li><b>folium で対話マップ</b>: lat/lon があるので
<code>folium.CircleMarker</code> で <a href="https://hiroshima-dobox.jp/">DoBoX</a>
のオープンデータを可視化し、クラスタ凡例付きで HTML 埋め込みする
(L076-L079 参照)。</li>

<li><b>クラスタ → 流域連携</b>: L080 (水系単位の雨量×水位ラグ)
と接続。クラスタを水系メタ (#1276) で塗り分ければ、
「気象学的クラスタ」と「流域単位」のズレが見える。</li>
</ol>"""),
]

html = render_lesson(
    num=7,
    title="k-means による雨量観測所のクラスタリング — 時間プロファイル × 地理分布",
    tags=["v2-rewrite", "選択", "応用基礎", "教師なし学習", "ML入門"],
    time="100分",
    level="応用基礎",
    data_label='2データ横断: '
               '<a href="https://hiroshima-dobox.jp/datasets/1275">#1275</a> (雨量10分値 14日) + '
               '<a href="https://hiroshima-dobox.jp/datasets/1274">#1274</a> (観測所一覧 lat/lon)',
    sections=sections,
    script_filename="lessons/L07_kmeans_stations.py",
)
out = LESSONS / "L07_kmeans_stations.html"
out.write_text(html, encoding="utf-8")
print(f"\nsaved: {out}")
print(f"  PNGs: 5 generated in {ASSETS}")