#!/usr/bin/env python3
"""Score Frey and Osborne's 2013 automation probabilities against BLS employment change.

Reproduces the figures in the Save The 1 Percent video "It Took Twelve Years to Show"
(savethe1percent.com/videos/002): the six-job scoreboard, and employment growth by 2013
risk band on the occupations counted under the same code at every date.

Inputs are downloaded into --workdir:
  - Frey and Osborne (2013) working paper, Appendix ranking of 702 occupations
  - BLS Occupational Employment and Wage Statistics, national, May 2013 and a later year

Only occupations whose Standard Occupational Classification (SOC) code is identical in the
2010 and 2018 SOC are compared. Reclassified occupations drop out; the script reports how many.

The BLS site refuses downloads without a contact address, so set your own first:
    export CONTACT_EMAIL=you@example.org
Needs pandas, openpyxl, xlrd, and pdftotext on PATH:
    uv run --with pandas --with openpyxl --with xlrd fo_vs_oews.py
    uv run --with pandas --with openpyxl --with xlrd fo_vs_oews.py --panel 18,22,25

Copyright (c) 2026 Stephen Angeloni. MIT License: see LICENSE.
"""
import argparse
import os
import re
import subprocess
import urllib.request
import zipfile
from pathlib import Path

import numpy as np
import pandas as pd

FO_URL = "https://www.oxfordmartin.ox.ac.uk/downloads/academic/The_Future_of_Employment.pdf"
OEWS_URL = "https://www.bls.gov/oes/special-requests/oesm{yy}nat.zip"


def fetch(url, dest):
    if not dest.exists():
        contact = os.environ.get("CONTACT_EMAIL")
        if not contact:
            raise SystemExit("Set CONTACT_EMAIL to your email address: the BLS site refuses "
                             "downloads without one.")
        req = urllib.request.Request(url, headers={"User-Agent": f"research-script ({contact})"})
        dest.write_bytes(urllib.request.urlopen(req).read())
    return dest


def frey_osborne(work):
    pdf = fetch(FO_URL, work / "fo2013.pdf")
    text = subprocess.run(["pdftotext", "-layout", str(pdf), "-"],
                          capture_output=True, text=True, check=True).stdout
    rows = []
    for line in text.splitlines():
        m = re.match(r"\s*(\d+)\.\s+([01]?\.\d+|[01])\s+(?:[01]\s+)?(\d{2}-\d{4})\s+(.*\S)", line)
        if m:
            rows.append((int(m[1]), float(m[2]), m[3], m[4]))
    fo = pd.DataFrame(rows, columns=["rank", "p", "code", "title"]).drop_duplicates("rank")
    assert set(fo["rank"]) == set(range(1, 703)), "appendix parse incomplete"
    return fo


def oews(work, yy):
    z = fetch(OEWS_URL.format(yy=yy), work / f"oesm{yy}nat.zip")
    with zipfile.ZipFile(z) as zf:
        name = next(n for n in zf.namelist() if re.search(r"national_M20\d\d_dl\.xlsx?$", n))
        zf.extract(name, work)
    df = pd.read_excel(work / name, dtype=str)
    df.columns = [c.upper() for c in df.columns]
    group = "OCC_GROUP" if "OCC_GROUP" in df else "O_GROUP"
    df = df[df[group] == "detailed"]
    return df.set_index("OCC_CODE")["TOT_EMP"].pipe(pd.to_numeric, errors="coerce")


def panel(fo, work, years):
    e13 = oews(work, "13")
    ends = {y: oews(work, y) for y in years}
    d = fo.assign(e13=fo.code.map(e13), **{f"e{y}": fo.code.map(e) for y, e in ends.items()}).dropna()
    d["band"] = pd.cut(d.p, [-1, 0.3, 0.7, 1.1], labels=["low (<=0.3)", "medium", "high (>0.7)"])
    print(f"Occupations with an identical SOC code in May 2013 and every end year: {len(d)}")
    rows = {}
    for y in years:
        g = d.groupby("band", observed=True)[["e13", f"e{y}"]].sum()
        rows[f"May 20{y}"] = (g[f"e{y}"] / g.e13 - 1) * 100
    t = pd.DataFrame(rows).T
    t["gap, low minus high"] = t.iloc[:, 0] - t.iloc[:, -1]
    print("Employment growth since May 2013, percent")
    print(t.round(1).to_string())
    print(d.groupby("band", observed=True).size().to_string())


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--workdir", default="/tmp/fo_vs_oews")
    ap.add_argument("--end", default="25", help="two-digit OEWS year to compare with May 2013")
    ap.add_argument("--panel", help="comma-separated end years, e.g. 18,22,25: band growth on the "
                                    "occupations present in every year")
    args = ap.parse_args()
    work = Path(args.workdir)
    work.mkdir(parents=True, exist_ok=True)

    fo = frey_osborne(work)
    if args.panel:
        return panel(fo, work, args.panel.split(","))
    e13, eN = oews(work, "13"), oews(work, args.end)
    m = fo.assign(e13=fo.code.map(e13), eN=fo.code.map(eN)).dropna(subset=["e13", "eN"])
    m["logchg"] = np.log(m.eN / m.e13)
    print(f"Matched {len(m)} of {len(fo)} occupations (identical SOC code, both years published)")
    print(f"Spearman rank correlation, probability vs log employment change: "
          f"{m[['p', 'logchg']].corr(method='spearman').iloc[0, 1]:.3f}")

    m["band"] = pd.cut(m.p, [-1, 0.3, 0.7, 1.1], labels=["low (<=0.3)", "medium", "high (>0.7)"])
    bands = m.groupby("band", observed=True).agg(n=("code", "size"), e13=("e13", "sum"),
                                                 eN=("eN", "sum"),
                                                 shrank=("logchg", lambda x: int((x < 0).sum())))
    bands["growth_pct"] = (bands.eN / bands.e13 - 1) * 100
    print(bands.round(1).to_string())

    m["decile"] = pd.qcut(m.p.rank(method="first"), 10, labels=False) + 1
    dec = m.groupby("decile").agg(p_min=("p", "min"), p_max=("p", "max"), n=("code", "size"),
                                  e13=("e13", "sum"), eN=("eN", "sum"),
                                  shrank_pct=("logchg", lambda x: (x < 0).mean() * 100))
    dec["growth_pct"] = (dec.eN / dec.e13 - 1) * 100
    print(dec.round(2).to_string())
    m.to_csv(work / f"fo_vs_oews_20{args.end}.csv", index=False)


if __name__ == "__main__":
    main()