"""Reproducible analysis for SHWRP-2026-021.

The study separates long-run employment-stock change from contemporary vacancy,
recruitment-duration, and turnover signals. It is descriptive screening rather
than a causal estimate or a forecast of firm-level labor shortages.
"""
from __future__ import annotations

import json
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.font_manager as fm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

ROOT = Path(__file__).resolve().parent
OUT = ROOT / "outputs"
FIG = OUT / "figures"
OUT.mkdir(exist_ok=True)
FIG.mkdir(exist_ok=True)

FONT_PATH = ROOT.parents[1] / "public" / "fonts" / "NotoSansTC-Regular.ttf"
if FONT_PATH.exists():
    fm.fontManager.addfont(str(FONT_PATH))
    plt.rcParams["font.family"] = "Noto Sans TC"
plt.rcParams.update({"axes.unicode_minus": False, "figure.dpi": 170, "savefig.dpi": 190})

LABELS = {
    "agriculture": "農林漁牧",
    "mining": "礦業土石",
    "manufacturing": "製造",
    "electricity_gas": "電力燃氣",
    "water_remediation": "用水污染整治",
    "construction": "營建工程",
    "wholesale_retail": "批發零售",
    "transport_storage": "運輸倉儲",
    "accommodation_food": "住宿餐飲",
    "information_communication": "出版影音資通訊",
    "finance_insurance": "金融保險",
    "real_estate": "不動產",
    "professional_scientific": "專業科學技術",
    "support_services": "支援服務",
    "public_administration": "公共行政國防",
    "education": "教育",
    "health_social": "醫療社會工作",
    "arts_recreation": "藝術休閒",
    "other_services": "其他服務",
}

emp = pd.read_csv(ROOT / "employment-by-industry-2014-2025.csv")
vac = pd.read_csv(ROOT / "vacancy-march-2025.csv")
agr = pd.read_csv(ROOT / "agriculture-succession-demand.csv")

records = []
for col, label in LABELS.items():
    start = float(emp.loc[emp.year == 2014, col].iloc[0])
    end = float(emp.loc[emp.year == 2025, col].iloc[0])
    pct = (end / start - 1) * 100
    cagr = (end / start) ** (1 / 11) - 1
    slope = float(np.polyfit(emp.year, emp[col], 1)[0])
    records.append({
        "industry_code": col,
        "industry": label,
        "employment_2014_thousand": start,
        "employment_2025_thousand": end,
        "change_thousand": end - start,
        "percent_change": pct,
        "cagr_percent": cagr * 100,
        "ols_slope_thousand_per_year": slope,
    })

trend = pd.DataFrame(records).merge(
    vac[["employment_column", "vacancies", "vacancy_rate", "entry_rate", "exit_rate", "recruitment_months"]],
    how="left", left_on="industry_code", right_on="employment_column"
).drop(columns=["employment_column"])
trend["turnover_gap_pp"] = trend["exit_rate"] - trend["entry_rate"]
trend["employment_state"] = np.select(
    [trend.percent_change <= -5, trend.percent_change > 5],
    ["明顯下降", "明顯增加"], default="大致持平"
)
trend["vacancy_pressure"] = (trend.vacancy_rate >= 3.1) | (trend.recruitment_months >= 3.5)

def classify(row: pd.Series) -> str:
    code = row.industry_code
    if code == "agriculture":
        return "人力退場與接班風險（農業專案證據）"
    if code == "public_administration":
        return "資料範圍外／不判定"
    if row.percent_change <= -5:
        return "小基數專業替補風險" if row.recruitment_months >= 3.5 else "需求仍待確認的收縮"
    if row.percent_change > 5 and row.vacancy_pressure:
        if row.turnover_gap_pp >= 0.3:
            return "擴張伴隨留任壓力"
        return "需求擴張快於補實"
    if row.percent_change <= 5 and row.vacancy_pressure:
        return "存量持平但技能／補實壓力"
    if row.turnover_gap_pp >= 0.3:
        return "高流動留任壓力"
    return "本次未見強烈背離"

trend["screening_type"] = trend.apply(classify, axis=1)
trend = trend.sort_values("percent_change")
trend.to_csv(OUT / "industry-screening.csv", index=False, encoding="utf-8-sig")

selected = ["agriculture", "manufacturing", "construction", "transport_storage", "accommodation_food", "health_social"]
indexed = emp[["year"] + selected].copy()
for col in selected:
    indexed[col] = indexed[col] / indexed.loc[indexed.year == 2014, col].iloc[0] * 100
indexed.to_csv(OUT / "employment-index-2014-2025.csv", index=False, encoding="utf-8-sig")

colors = {"brown": "#7a3f22", "orange": "#c97832", "green": "#2f6b4f", "blue": "#376a8a", "gray": "#8c8178", "cream": "#f6efe5", "red": "#a74c43"}

# Figure 1: indexed employment paths.
fig, ax = plt.subplots(figsize=(10.8, 6.2))
line_colors = [colors["brown"], colors["gray"], colors["orange"], colors["blue"], colors["red"], colors["green"]]
for col, color in zip(selected, line_colors):
    ax.plot(indexed.year, indexed[col], marker="o", lw=2.1, ms=3.8, label=LABELS[col], color=color)
ax.axhline(100, color="#bdb3aa", lw=1, ls="--")
ax.set(title="圖1  2014–2025年六類產業就業人數指數", xlabel="年", ylabel="2014年＝100")
ax.legend(ncol=3, frameon=False, loc="upper left")
ax.grid(axis="y", color="#e5ddd5", lw=.7)
fig.tight_layout(); fig.savefig(FIG / "figure-1-employment-index.png", bbox_inches="tight"); plt.close(fig)

# Figure 2: endpoint changes across all broad industries.
plot = trend.sort_values("percent_change")
fig, ax = plt.subplots(figsize=(10.8, 7.2))
bar_colors = [colors["red"] if v <= -5 else colors["green"] if v > 5 else colors["orange"] for v in plot.percent_change]
ax.barh(plot.industry, plot.percent_change, color=bar_colors)
ax.axvline(-5, color="#8c8178", ls="--", lw=1); ax.axvline(5, color="#8c8178", ls="--", lw=1)
ax.axvline(0, color="#33271f", lw=.8)
for y, v in enumerate(plot.percent_change): ax.text(v + (0.7 if v >= 0 else -0.7), y, f"{v:+.1f}%", ha="left" if v >= 0 else "right", va="center", fontsize=8)
ax.set(title="圖2  19類產業就業人數端點變化", xlabel="2014至2025年變化（%）", ylabel="")
ax.grid(axis="x", color="#e5ddd5", lw=.7)
fig.tight_layout(); fig.savefig(FIG / "figure-2-industry-change.png", bbox_inches="tight"); plt.close(fig)

# Figure 3: vacancy rate and recruitment duration.
joined = vac.merge(trend[["industry_code", "percent_change"]], left_on="employment_column", right_on="industry_code")
fig, ax = plt.subplots(figsize=(10.8, 6.6))
sc = ax.scatter(joined.vacancy_rate, joined.recruitment_months, s=np.sqrt(joined.vacancies) * 5.2,
                c=joined.percent_change, cmap="RdYlGn", vmin=-10, vmax=30, alpha=.8, edgecolor="white", linewidth=.8)
ax.axvline(3.1, color="#8c8178", ls="--", lw=1); ax.axhline(3.5, color="#8c8178", ls="--", lw=1)
for _, r in joined.iterrows():
    if r.vacancy_rate >= 4 or r.recruitment_months >= 4.1 or r.vacancies >= 16000:
        ax.annotate(LABELS[r.employment_column], (r.vacancy_rate, r.recruitment_months), xytext=(4, 4), textcoords="offset points", fontsize=8)
ax.set(title="圖3  職缺率、招募時間與長期就業變化", xlabel="2025年3月底職缺率（%）", ylabel="全時職缺平均招募時間（月）")
fig.colorbar(sc, ax=ax, label="2014–2025就業變化（%）")
ax.grid(color="#e5ddd5", lw=.7)
fig.tight_layout(); fig.savefig(FIG / "figure-3-vacancy-recruitment.png", bbox_inches="tight"); plt.close(fig)

# Figure 4: employment change vs combined demand-pressure index.
joined["demand_pressure_index"] = ((joined.vacancy_rate - 3.1) / joined.vacancy_rate.std(ddof=0) +
                                     (joined.recruitment_months - 3.5) / joined.recruitment_months.std(ddof=0)) / 2
fig, ax = plt.subplots(figsize=(10.8, 6.5))
ax.scatter(joined.percent_change, joined.demand_pressure_index, s=75, color=colors["brown"], alpha=.8)
ax.axvline(-5, color="#8c8178", ls="--"); ax.axvline(5, color="#8c8178", ls="--"); ax.axhline(0, color="#8c8178", ls="--")
for _, r in joined.iterrows():
    ax.annotate(LABELS[r.employment_column], (r.percent_change, r.demand_pressure_index), xytext=(3, 3), textcoords="offset points", fontsize=7.5)
ax.text(-24, 2.0, "人力下降＋補實壓力", color=colors["red"], fontsize=10)
ax.text(12, 2.0, "需求擴張＋補實壓力", color=colors["green"], fontsize=10)
ax.set(title="圖4  非農產業需求—勞動供給背離篩選", xlabel="2014–2025就業變化（%）", ylabel="需求壓力指數（職缺率與招募時間標準化平均）")
ax.grid(color="#e5ddd5", lw=.7)
fig.tight_layout(); fig.savefig(FIG / "figure-4-divergence-matrix.png", bbox_inches="tight"); plt.close(fig)

# Figure 5: entry and exit rates.
fig, ax = plt.subplots(figsize=(9.8, 6.5))
ax.scatter(vac.entry_rate, vac.exit_rate, s=np.sqrt(vac.vacancies) * 5, color=colors["orange"], alpha=.8, edgecolor="white")
lim = max(vac.entry_rate.max(), vac.exit_rate.max()) + .3
ax.plot([0, lim], [0, lim], color="#8c8178", ls="--", lw=1)
for _, r in vac.iterrows():
    gap = r.exit_rate - r.entry_rate
    if abs(gap) >= .3 or r.vacancies >= 16000:
        ax.annotate(LABELS[r.employment_column], (r.entry_rate, r.exit_rate), xytext=(4, 3), textcoords="offset points", fontsize=8)
ax.set(title="圖5  進入率與退出率：把留任問題與存量下降分開", xlabel="2025年3月進入率（%）", ylabel="2025年3月退出率（%）", xlim=(0, lim), ylim=(0, lim))
ax.grid(color="#e5ddd5", lw=.7)
fig.tight_layout(); fig.savefig(FIG / "figure-5-entry-exit.png", bbox_inches="tight"); plt.close(fig)

# Figure 6: agricultural succession risk.
age = agr[agr.indicator == "manager_age"].copy()
fig, ax1 = plt.subplots(figsize=(10.2, 6.1))
ax2 = ax1.twinx()
ax1.plot(age.year, age.manager_mean_age, color=colors["brown"], marker="o", lw=2.5, label="平均年齡")
ax2.bar(age.year, age.share_age_65_plus, width=2.6, color="#d9c7b5", alpha=.75, label="65歲以上占比")
ax1.set(title="圖6  農業經營管理者高齡化與接班風險", xlabel="農林漁牧業普查年", ylabel="平均年齡（歲）")
ax2.set_ylabel("65歲以上占比（%）")
ax1.set_xticks(age.year)
lines, labels = ax1.get_legend_handles_labels(); bars, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines + bars, labels + labels2, frameon=False, loc="upper left")
ax1.grid(axis="y", color="#e5ddd5", lw=.7)
fig.tight_layout(); fig.savefig(FIG / "figure-6-agriculture-succession.png", bbox_inches="tight"); plt.close(fig)

ag_start = float(emp.loc[emp.year == 2014, "agriculture"].iloc[0])
ag_end = float(emp.loc[emp.year == 2025, "agriculture"].iloc[0])
mining = trend.loc[trend.industry_code == "mining"].iloc[0]
results = {
    "report_number": "SHWRP-2026-021",
    "period": [2014, 2025],
    "industry_count": len(LABELS),
    "vacancy_industry_count": len(vac),
    "employment_total_change_percent": float((emp.total.iloc[-1] / emp.total.iloc[0] - 1) * 100),
    "agriculture": {
        "employment_2014_thousand": ag_start,
        "employment_2025_thousand": ag_end,
        "change_thousand": ag_end - ag_start,
        "percent_change": float((ag_end / ag_start - 1) * 100),
        "manager_mean_age_2020": 64.13,
        "share_age_65_plus_2020": 48.76,
        "farm_sales_change_2015_2020_percent": 7.0,
        "rice_area_change_2015_2020_percent": 7.6,
        "screening_conclusion": "high-confidence labor-exit and succession-risk candidate; demand persistence is supported by separate agricultural indicators, not by the vacancy survey"
    },
    "mining": {
        "employment_change_percent": float(mining.percent_change),
        "employment_2025_thousand": float(mining.employment_2025_thousand),
        "vacancy_rate": float(mining.vacancy_rate),
        "recruitment_months": float(mining.recruitment_months),
        "screening_conclusion": "small-base specialist replacement watchlist; broad demand persistence not established"
    },
    "selected_expansion_pressure": trend[trend.screening_type.str.contains("需求擴張", na=False)][["industry", "percent_change", "vacancy_rate", "recruitment_months", "screening_type"]].to_dict("records"),
    "method_thresholds": {"employment_decline_percent": -5, "employment_growth_percent": 5, "vacancy_rate_percent": 3.1, "recruitment_months": 3.5},
    "limitations": [
        "Vacancy data are a March 2025 cross-section and exclude agriculture and public administration.",
        "Industry classifications changed across the period; the official linked series is used but discontinuity risk remains.",
        "Agricultural output value and sales are nominal and do not alone identify real demand or profitability.",
        "Broad-industry averages can hide occupation, region, contract, wage, and skill heterogeneity."
    ]
}
(ROOT / "analysis-results.json").write_text(json.dumps(results, ensure_ascii=False, indent=2), encoding="utf-8")
print(json.dumps(results, ensure_ascii=False, indent=2))
