stack/src/pfs/pipe.py

"""PFS data pipeline — Zotero library → pfs tables → DuckDB.

Discovers CMS Physician Fee Schedule data files stored in Zotero
by tag, reads Excel/CSV with polars, normalizes columns using
``pfs.files`` column maps, and loads into DuckDB (or any lake context).

File discovery uses three Zotero tags:

- ``file:rvu`` — Addendum B (RVU files), PPRRVU quarterly releases
- ``file:gpci`` — Addendum D/E (GAFs/GPCIs), GPCI quarterly files
- ``file:pe`` — PUF files (Labor, Supply, Equipment), Work Time, PEHR

Each file type has a reader that handles CMS header quirks
(multi-row titles, merged cells, varying column names across years).

Usage::

    from pfs.pipe import load_all
    import duckdb

    con = duckdb.connect(path("db.aco"))  # path from conf
    stats = load_all(con)
    for table, info in stats.items():
        print(f"{table}: {info['rows']} rows from {info['files']} files")
"""

from __future__ import annotations

import logging
import re
import sqlite3
from pathlib import Path
from typing import Any

import duckdb
import polars as pl

from pfs.calcs.payment import LIMITING_CHARGE_FACTOR

log = logging.getLogger(__name__)

# ── Column maps from pfs.files (import at use time to avoid cycles) ──

_RVU_COLUMNS = {
    "HCPCS": "hcpcs",
    "CPT1/\nHCPCS": "hcpcs",
    "CPT1/HCPCS": "hcpcs",
    "MOD": "mod",
    "Mod": "mod",
    "DESCRIPTION": "description",
    "STATUS CODE": "status_code",
    "STATUS\nCODE": "status_code",
    "STATUS": "status_code",
    "NOT USED FOR MEDICARE PAYMENT": "not_used_for_medicare",
    "NOT USED\nFOR\nMEDICARE \nPAYMENT": "not_used_for_medicare",
    "WORK RVU": "work_rvu",
    "WORK\nRVU": "work_rvu",
    "\nWork\nRVUs2": "work_rvu",
    "WORK\nRVUs2": "work_rvu",
    "NON-FACILITY PE RVU": "non_fac_pe_rvu",
    "NON-FAC PE RVU": "non_fac_pe_rvu",
    "NON-FAC\nPE RVU": "non_fac_pe_rvu",
    "\nNon-\nFacility\nPE\nRVUs2": "non_fac_pe_rvu",
    "FACILITY PE RVU": "fac_pe_rvu",
    "FAC PE RVU": "fac_pe_rvu",
    "FACILITY\nPE RVU": "fac_pe_rvu",
    "\nFacility\nPE\nRVUs2": "fac_pe_rvu",
    "MALPRACTICE RVU": "mp_rvu",
    "MP RVU": "mp_rvu",
    "MP\nRVU": "mp_rvu",
    "Mal-\nPractice\nRVUs2": "mp_rvu",
    "NON-FACILITY TOTAL": "non_fac_total",
    "NON-FACILITY\nTOTAL": "non_fac_total",
    "\nTotal Non-Facility RVUs2": "non_fac_total",
    "FACILITY TOTAL": "fac_total",
    "FACILITY\nTOTAL": "fac_total",
    "\nTotal Facility RVUs2": "fac_total",
    "PCTC IND": "pctc_ind",
    "PCTC\nIND": "pctc_ind",
    "GLOB DAYS": "glob_days",
    "GLOB\nDAYS": "glob_days",
    "Global": "glob_days",
    "GLOBAL": "glob_days",
    "PRE OP": "pre_op",
    "PRE\nOP": "pre_op",
    "INTRA OP": "intra_op",
    "INTRA\nOP": "intra_op",
    "POST OP": "post_op",
    "POST\nOP": "post_op",
    "MULT PROC": "mult_proc",
    "MULT\nPROC": "mult_proc",
    "BILAT SURG": "bilat_surg",
    "BILAT\nSURG": "bilat_surg",
    "ASST SURG": "asst_surg",
    "ASST\nSURG": "asst_surg",
    "CO SURG": "co_surg",
    "CO-\nSURG": "co_surg",
    "TEAM SURG": "team_surg",
    "TEAM\nSURG": "team_surg",
    "ENDO BASE": "endo_base",
    "ENDO\nBASE": "endo_base",
    "CONV FACTOR": "conv_factor",
    "CONV\nFACTOR": "conv_factor",
    "PHYSICIAN SUPERVISION OF DIAGNOSTIC PROCEDURES": "physician_supervision",
    "PHYSICIAN\nSUPERVISION OF\nDIAGNOSTIC\nPROCEDURES": "physician_supervision",
    "CALCULATION FLAG": "calc_flag",
    "CALCULATION\nFLAG": "calc_flag",
    "DIAGNOSTIC IMAGING FAMILY INDICATOR": "diag_imaging_family",
    "DIAGNOSTIC\nIMAGING\nFAMILY\nINDICATOR": "diag_imaging_family",
    "NON-FACILITY PE USED FOR OPPS PAYMENT AMOUNT": "non_fac_pe_used_for_opps",
    "FACILITY PE USED FOR OPPS PAYMENT AMOUNT": "fac_pe_used_for_opps",
    "MP USED FOR OPPS PAYMENT AMOUNT": "mp_used_for_opps",
    "PRIC IND": "pctc_ind",
    "PRIC\nIND": "pctc_ind",
}

_GPCI_COLUMNS = {
    "MAC": "mac",
    "Medicare Administrative Contractor (MAC)": "mac",
    "CARRIER": "mac",
    "Carrier": "mac",
    "LOCALITY": "locality",
    "Locality Number": "locality",
    "LOCALITY NAME": "locality_name",
    "Locality Name": "locality_name",
    "STATE": "state",
    "State": "state",
}

_LABOR_COLUMNS = {
    "hcpcs": "hcpcs",
    "HCPCS": "hcpcs",
    "labor_code": "activity_code",
    "source": "_source",
    "description": "activity_description",
    "rate_per_minute": "rate_per_minute",
    "nf_pre_svc": "_nf_pre_svc",
    "nf_svc": "_nf_svc",
    "nf_post_svc": "_nf_post_svc",
    "f_pre_svc": "_f_pre_svc",
    "f_svc": "_f_svc",
    "f_post_svc": "_f_post_svc",
    "global_period": "_global_period",
}

_SUPPLY_COLUMNS = {
    "hcpcs": "hcpcs",
    "HCPCS": "hcpcs",
    "cms_code": "supply_code",
    "description": "supply_description",
    "unit": "_unit",
    "price": "unit_price",
    "nf_quantity": "nf_quantity",
    "f_quantity": "f_quantity",
    "source": "_source",
    "category": "_category",
    "global_period": "_global_period",
}

_EQUIPMENT_COLUMNS = {
    "hcpcs": "hcpcs",
    "HCPCS": "hcpcs",
    "cms_code": "equipment_code",
    "description": "equipment_description",
    "useful_life": "useful_life",
    "price": "unit_price",
    "utilization_rate": "_utilization_rate",
    "minutes_per_year": "minutes_per_year",
    "nf_time": "nf_minutes",
    "f_time": "f_minutes",
    "source": "_source",
    "category": "_category",
    "global_period": "_global_period",
}

_WORK_TIME_COLUMNS = {
    "cpt_code": "hcpcs",
    "CPT_CODE": "hcpcs",
    "Cpt_code": "hcpcs",
    "HCPCS": "hcpcs",
    "Pre_Evaluation_Time": "pre_eval_time",
    "Pre_Positioning_time": "pre_positioning_time",
    "Pre_Service_Scrub_Dress_Wait_tim": "pre_scrub_dress_wait",
    "Median_Intra_Service_Time": "intra_service_time",
    "Immediate_post_Service_time": "post_immediate_time",
    "RUC_Recommended_Adjusted_Time": "total_time",
}


# ── Zotero discovery ─────────────────────────────────────────


def _zotero_pfs_files(
    zotero_db: str | Path,
    zotero_storage: str | Path,
) -> list[dict]:
    """Query Zotero for all module:pfs attachments (xlsx/xls/csv).

    Returns list of dicts with keys: item_key, year, title, path,
    filename, ext.  Files are classified by filename pattern later.
    """
    db = sqlite3.connect(str(zotero_db))
    # All items tagged module:pfs
    rows = db.execute(
        """
        SELECT DISTINCT i.itemID, i.key
        FROM items i
        JOIN itemTags it ON i.itemID = it.itemID
        JOIN tags t ON it.tagID = t.tagID
        WHERE t.name = 'module:pfs'
        """,
    ).fetchall()

    results = []
    for item_id, item_key in rows:
        year_rows = db.execute(
            """
            SELECT t.name FROM itemTags it
            JOIN tags t ON it.tagID = t.tagID
            WHERE it.itemID = ? AND t.name LIKE 'year:%'
            """,
            (item_id,),
        ).fetchall()
        year = int(year_rows[0][0].split(":")[1]) if year_rows else 0

        # Release quarter tag (q1..q4) — used by _load_rvu_files to
        # prefer the latest quarter when multiple releases exist for
        # the same calendar year.
        release_rows = db.execute(
            """
            SELECT t.name FROM itemTags it
            JOIN tags t ON it.tagID = t.tagID
            WHERE it.itemID = ? AND t.name LIKE 'release:%'
            """,
            (item_id,),
        ).fetchall()
        release = release_rows[0][0].split(":", 1)[1] if release_rows else ""

        title_row = db.execute(
            """
            SELECT value FROM itemDataValues idv
            JOIN itemData id ON idv.valueID = id.valueID
            JOIN fields f ON id.fieldID = f.fieldID
            WHERE id.itemID = ? AND f.fieldName = 'title'
            """,
            (item_id,),
        ).fetchone()
        title = title_row[0] if title_row else ""

        attachments = db.execute(
            """
            SELECT ia.itemID, ia.path FROM itemAttachments ia
            WHERE ia.parentItemID = ?
            AND (ia.path LIKE '%.xlsx' OR ia.path LIKE '%.xls'
                 OR ia.path LIKE '%.csv' OR ia.path LIKE '%.txt')
            """,
            (item_id,),
        ).fetchall()

        for att_id, att_path in attachments:
            filename = att_path.replace("storage:", "")
            # Skip macOS resource fork files
            if filename.startswith("._"):
                continue
            storage_key = db.execute(
                "SELECT key FROM items WHERE itemID = ?", (att_id,)
            ).fetchone()
            if not storage_key:
                continue
            full_path = Path(zotero_storage) / storage_key[0] / filename
            if not full_path.exists():
                continue
            ext = full_path.suffix.lower()
            results.append(
                {
                    "item_key": item_key,
                    "year": year,
                    "release": release,  # "q1"/"q2"/"q3"/"q4" or ""
                    "title": title,
                    "path": str(full_path),
                    "filename": filename,
                    "ext": ext,
                }
            )

    db.close()
    return results


# ── Excel/CSV readers with CMS header detection ─────────────


def _find_header_row(path: str, *, min_cols: int = 5) -> int:
    """Find the row index containing column headers in an Excel file.

    Scans for the first row where >= min_cols cells are non-null
    strings that look like column names (contain alpha chars, not
    just copyright/title boilerplate).
    """
    import openpyxl

    wb = openpyxl.load_workbook(path, read_only=True)
    ws = wb.active
    for i, row in enumerate(ws.iter_rows(values_only=True)):
        non_null = [c for c in row if c is not None]
        str_cells = [c for c in non_null if isinstance(c, str) and len(c.strip()) < 60]
        if len(str_cells) >= min_cols:
            wb.close()
            return i
    wb.close()
    return 0


def _read_excel(path: str, header_row: int | None = None) -> pl.DataFrame:
    """Read an Excel file with auto-detected or explicit header row."""
    if header_row is None:
        header_row = _find_header_row(path)
    df = pl.read_excel(path, read_options={"header_row": header_row})
    return df


def _read_pprrvu(path: str) -> pl.DataFrame:
    """Read a PPRRVU quarterly file with its multi-row stacked header.

    PPRRVU files have ~5 title rows, then a 4-row stacked header
    (rows 5-9 typically).  The last header row contains the actual
    short column names (``HCPCS, MOD, DESCRIPTION, CODE, PAYMENT,
    RVU, PE RVU, INDICATOR, PE RVU, INDICATOR, RVU, TOTAL, TOTAL,
    IND, DAYS, OP, OP, OP, PROC, SURG, SURG, SURG, SURG, BASE,
    FACTOR, PROCEDURES, FLAG, INDICATOR, AMOUNT, AMOUNT, AMOUNT``).
    Several short names repeat across unrelated columns (``RVU``,
    ``PE RVU``, ``INDICATOR``, ``TOTAL``, ``SURG``, ``OP``, ``AMOUNT``)
    so a name-based rename can't disambiguate them — we map
    **positionally** to our schema.

    Two historical gotchas this function now guards against:

    1. ``pl.read_excel`` silently ignores ``read_options['header_row']``
       unless the ``calamine`` engine is active — with the default
       xlsx2csv engine it returns a 1-column DataFrame of the title
       row.  We pin ``engine='calamine'``.
    2. The positional rename previously required exactly 32 columns
       (``_pric_ind`` was a phantom entry between ``team_surg`` and
       ``endo_base`` that does not exist in any real PPRRVU file).
       The list now has the correct 31 entries and the rename gate
       uses ``==`` on the true count.
    """
    import openpyxl

    wb = openpyxl.load_workbook(path, read_only=True)
    ws = wb.active
    hcpcs_row = 0
    for i, row in enumerate(ws.iter_rows(values_only=True)):
        cells = [str(c).strip() if c is not None else "" for c in row]
        if cells and cells[0] == "HCPCS":
            hcpcs_row = i
            break
    wb.close()

    # calamine respects read_options['header_row']; xlsx2csv does not.
    df = pl.read_excel(
        path,
        engine="calamine",
        read_options={"header_row": hcpcs_row},
    )

    # CMS PPRRVU schema changed in CY2026 — a "PRIC IND" (pricing
    # indicator) column was inserted at position 23 between TEAM SURG
    # and ENDO BASE, pushing the later columns (ENDO BASE, CONV FACTOR,
    # PHYSICIAN SUPERVISION, ...) out by one.  Detect the layout by
    # column count and pick the matching positional map.  Columns
    # starting with "_" are stripped after renaming.
    _EXPECTED_31 = [
        "hcpcs",  # 0  HCPCS
        "mod",  # 1  MOD
        "description",  # 2  DESCRIPTION
        "status_code",  # 3  STATUS CODE
        "not_used_for_medicare",  # 4  NOT USED FOR MEDICARE PAYMENT
        "work_rvu",  # 5  WORK RVU
        "non_fac_pe_rvu",  # 6  NON-FAC PE RVU
        "_nf_na_indicator",  # 7  NON-FAC NA INDICATOR
        "fac_pe_rvu",  # 8  FACILITY PE RVU
        "_f_na_indicator",  # 9  FACILITY NA INDICATOR
        "mp_rvu",  # 10 MP RVU
        "non_fac_total",  # 11 NON-FACILITY TOTAL
        "fac_total",  # 12 FACILITY TOTAL
        "pctc_ind",  # 13 PCTC IND
        "glob_days",  # 14 GLOB DAYS
        "pre_op",  # 15 PRE OP
        "intra_op",  # 16 INTRA OP
        "post_op",  # 17 POST OP
        "mult_proc",  # 18 MULT PROC
        "bilat_surg",  # 19 BILAT SURG
        "asst_surg",  # 20 ASST SURG
        "co_surg",  # 21 CO-SURG
        "team_surg",  # 22 TEAM SURG
        "endo_base",  # 23 ENDO BASE
        "conv_factor",  # 24 CONV FACTOR
        "physician_supervision",  # 25 PHYSICIAN SUPERVISION
        "calc_flag",  # 26 CALCULATION FLAG
        "diag_imaging_family",  # 27 DIAGNOSTIC IMAGING FAMILY
        "non_fac_pe_used_for_opps",  # 28 NON-FACILITY PE USED FOR OPPS
        "fac_pe_used_for_opps",  # 29 FACILITY PE USED FOR OPPS
        "mp_used_for_opps",  # 30 MP USED FOR OPPS
    ]

    # CY2026+ layout: "PRIC IND" slot at position 23, shifts everything
    # downstream by one.  Mark internal ("_pric_ind") because nothing
    # downstream of ingest uses it yet.
    _EXPECTED_32 = _EXPECTED_31[:23] + ["_pric_ind"] + _EXPECTED_31[23:]

    if len(df.columns) >= len(_EXPECTED_32):
        expected_order = _EXPECTED_32
    elif len(df.columns) >= len(_EXPECTED_31):
        expected_order = _EXPECTED_31
    else:
        log.warning(
            "PPRRVU %s has %d columns, expected 31 or 32 — positional rename skipped",
            path,
            len(df.columns),
        )
        return df

    rename = {
        old: new for old, new in zip(df.columns[: len(expected_order)], expected_order)
    }
    df = df.rename(rename)
    # Keep only the mapped columns — trailing extras are noise.
    df = df.select([c for c in expected_order if not c.startswith("_")])

    return df


def _read_csv_with_header(path: str) -> pl.DataFrame:
    """Read a CSV file, skipping title rows before the actual header."""
    # Read raw to find header
    with open(path, encoding="utf-8-sig") as f:
        import csv

        reader = csv.reader(f)
        for i, row in enumerate(reader):
            non_empty = [c for c in row if c.strip()]
            if len(non_empty) >= 3:
                # Check if this looks like a header (has alpha chars)
                alpha_cells = [c for c in non_empty if any(ch.isalpha() for ch in c)]
                if len(alpha_cells) >= 2:
                    break
        else:
            i = 0

    return pl.read_csv(path, skip_rows=i, encoding="utf8-lossy", infer_schema_length=0)


def _detect_encoding(path: str) -> str:
    """Detect if a file is UTF-16 (LE/BE) or UTF-8."""
    with open(path, "rb") as f:
        bom = f.read(4)
    if bom[:2] in (b"\xff\xfe", b"\xfe\xff"):
        return "utf-16"
    return "utf-8-sig"


def _read_tsv(path: str) -> pl.DataFrame:
    """Read a tab-separated .txt file, skipping leading title rows."""
    enc = _detect_encoding(path)
    with open(path, encoding=enc, errors="replace") as f:
        lines = f.readlines()

    # Strip null bytes that can appear in UTF-16 decoded text
    lines = [line.replace("\x00", "") for line in lines]

    # Find header row — first row with multiple tab-separated alpha tokens
    header_idx = 0
    for i, line in enumerate(lines):
        parts = line.strip().split("\t")
        alpha = [p for p in parts if p.strip() and any(c.isalpha() for c in p)]
        if len(alpha) >= 3:
            header_idx = i
            break

    # Re-join lines from header onward and read with polars
    text = "".join(lines[header_idx:])
    return pl.read_csv(
        text.encode("utf-8"),
        separator="\t",
        encoding="utf8-lossy",
        infer_schema_length=0,
        truncate_ragged_lines=True,
    )


def _normalize_columns(
    df: pl.DataFrame,
    column_map: dict[str, str],
) -> pl.DataFrame:
    """Rename columns using a case-insensitive map, drop unmapped."""
    # Build a normalized lookup: strip whitespace, collapse newlines
    lookup: dict[str, str] = {}
    for raw, target in column_map.items():
        # Normalize the key
        norm = raw.replace("\n", " ").strip()
        norm = re.sub(r"\s+", " ", norm)
        lookup[norm.upper()] = target
        # Also add the raw key
        lookup[raw.upper()] = target

    rename = {}
    keep = []
    for col in df.columns:
        norm = col.replace("\n", " ").strip()
        norm = re.sub(r"\s+", " ", norm)
        target = lookup.get(norm.upper()) or lookup.get(col.upper())
        if target and not target.startswith("_"):
            rename[col] = target
            keep.append(col)
        elif target and target.startswith("_"):
            # Internal columns — keep for now, rename
            rename[col] = target
            keep.append(col)

    if not keep:
        return df

    df = df.select(keep).rename(rename)
    # Drop internal columns
    internal = [c for c in df.columns if c.startswith("_")]
    if internal:
        df = df.drop(internal)
    return df


# ── File type readers ────────────────────────────────────────


def _is_addendum_b(filename: str) -> bool:
    fn = filename.upper()
    return "ADDENDUM" in fn and (
        "B " in fn or "B_" in fn or fn.startswith("ADDENDUM_B")
    )


def _is_pprrvu(filename: str) -> bool:
    return "PPRRVU" in filename.upper()


def _is_gpci_file(filename: str) -> bool:
    fn = filename.upper()
    if "GPCI" in fn and "ADDENDUM B" not in fn:
        return True
    if re.search(r"ADDENDUM[_ ]E\b", fn):
        return True
    return False


def _is_puf_labor(filename: str) -> bool:
    fn = filename.upper()
    return "PUF_LABOR" in fn and "TASK" not in fn and "OLD" not in fn


def _is_puf_supply(filename: str) -> bool:
    return "PUF_SUPPLY" in filename.upper()


def _is_puf_equip(filename: str) -> bool:
    fn = filename.upper()
    return "PUF_EQUIP" in fn or "PUF_EQUIPMENT" in fn


def _is_work_time(filename: str) -> bool:
    fn = filename.upper()
    return "WORK_TIME" in fn or "WORK TIME" in fn or "PHYSTIME" in fn


def _is_carrier_file(filename: str) -> bool:
    """Detect carrier fee schedule files (PFALL*.txt or PF{ST}{YY}*.TXT)."""
    fn = filename.upper()
    if fn.startswith("PFALL") and fn.endswith(".TXT"):
        return True
    # Per-state: PF{2-3 letters}{2-digit}{suffix}.TXT
    # e.g. PFAL26A.TXT, PFNY26R.TXT, PFPRV26.TXT, PFCA262.TXT
    if re.match(r"^PF[A-Z]{2,3}\d{2}", fn) and fn.endswith(".TXT"):
        return True
    return False


def _is_zip_carrier_file(filename: str) -> bool:
    """Detect ZIP code carrier locality files (ZIP5_*.txt)."""
    fn = filename.upper()
    return fn.startswith("ZIP5") and fn.endswith(".TXT") and "LYOUT" not in fn


def _extract_year_from_title(title: str) -> int | None:
    """Extract CY year from Zotero title like 'CY 2026 PFS Final — ...'"""
    m = re.search(r"CY\s*(\d{4})", title)
    return int(m.group(1)) if m else None


def _is_final_rule(title: str) -> bool:
    return "Final" in title or "final" in title


# ── Table loaders ────────────────────────────────────────────


def _load_rvu_files(
    files: list[dict],
    con: duckdb.DuckDBPyConnection,
) -> dict[str, Any]:
    """Load Addendum B and PPRRVU files into pfs.rvu.

    PPRRVU files are preferred over Addendum B because they include
    ~19 additional columns that Addendum B omits — notably
    ``conv_factor`` (the CMS-published conversion factor for the
    year), pre/intra/post-op indicators, surgical modifiers, and
    physician supervision codes.  See :pincite:`VVBEVYLC` for the
    CF history the PPRRVU column captures.  Addendum B is used as
    a fallback only when no PPRRVU file exists for a given year.
    """
    rvu_files = []
    for f in files:
        fn = f["filename"]
        if _is_addendum_b(fn) or _is_pprrvu(fn):
            # Prefer xlsx over csv/txt
            if f["ext"] in (".xlsx", ".xls"):
                rvu_files.append(f)
            elif f["ext"] == ".csv" and not any(
                g["year"] == f["year"]
                and g["ext"] in (".xlsx", ".xls")
                and (
                    (_is_addendum_b(g["filename"]) and _is_addendum_b(fn))
                    or (_is_pprrvu(g["filename"]) and _is_pprrvu(fn))
                )
                for g in files
            ):
                rvu_files.append(f)

    # Deduplicate: per year, pick the PPRRVU release whose embedded
    # conv_factor matches pfs.rules.RULES[year].conversion_factor —
    # the dominant CF per the AMA history table :pincite:`VVBEVYLC`.
    # For years with a mid-year CAA update (2015, 2024), RULES stores
    # the post-update CF as dominant, so this selects the post-update
    # quarterly release.  For years without a split, every quarter has
    # the same CF and we pick the earliest matching release.
    # Addendum B is used only as a last-resort fallback when no
    # PPRRVU file exists for the year.
    from pfs.rules import RULES

    _RELEASE_RANK = {"q1": 1, "q2": 2, "q3": 3, "q4": 4}

    def _release_sort_key(f: dict) -> tuple[int, int, str]:
        # Sort PPRRVU files so the earliest quarter is scanned first.
        # Addendum B files get release_rank = 99 so they're last;
        # within that tier, Final rule beats Proposed rule.
        is_pprrvu_file = _is_pprrvu(f["filename"])
        if not is_pprrvu_file:
            final_priority = 0 if _is_final_rule(f.get("title", "")) else 1
            return (99, final_priority, f["filename"])
        rank = _RELEASE_RANK.get(f.get("release", ""), 50)
        return (rank, 0, f["filename"])

    def _scan_cf(f: dict) -> float | None:
        """Read the embedded conv_factor from a PPRRVU file."""
        try:
            if _is_pprrvu(f["filename"]) and f["ext"] in (".xlsx", ".xls"):
                df = _read_pprrvu(f["path"])
                if "conv_factor" in df.columns and df.height > 0:
                    val = (
                        df.filter(pl.col("conv_factor").is_not_null())
                        .get_column("conv_factor")
                        .cast(pl.Float64, strict=False)
                        .head(1)
                    )
                    if len(val) > 0:
                        return float(val[0])
        except Exception as exc:
            log.debug("CF scan failed for %s: %s", f["filename"], exc)
        return None

    # Group files by year then pick the best release per year.
    files_by_year: dict[int, list[dict]] = {}
    for f in rvu_files:
        files_by_year.setdefault(f["year"], []).append(f)

    by_year: dict[int, dict] = {}
    for year, candidates in files_by_year.items():
        target_cf = RULES[year].conversion_factor if year in RULES else None
        ordered = sorted(candidates, key=_release_sort_key)
        pprrvu_candidates = [f for f in ordered if _is_pprrvu(f["filename"])]

        chosen = None

        # When multiple PPRRVU releases exist, pick the earliest whose
        # embedded CF matches RULES to within 1¢.  This catches the
        # mid-year split years (2015, 2024) cleanly.  Skip the scan
        # entirely when there is only one PPRRVU candidate — it will
        # be selected by the fallback below, avoiding a redundant read.
        if len(pprrvu_candidates) > 1 and target_cf is not None:
            for f in pprrvu_candidates:
                cf = _scan_cf(f)
                if cf is None:
                    continue
                if abs(cf - target_cf) < 0.005:
                    chosen = f
                    log.info(
                        "RVU %d: matched release=%s (%s) CF=%s to RULES dominant CF=%s",
                        year,
                        f.get("release", "?"),
                        f["filename"],
                        cf,
                        target_cf,
                    )
                    break
            if chosen is None:
                log.warning(
                    "RVU %d: no PPRRVU matched RULES CF %s; "
                    "falling back to earliest PPRRVU",
                    year,
                    target_cf,
                )

        # Fallback: earliest PPRRVU by rank, then Addendum B (Final
        # preferred over Proposed via _release_sort_key).
        if chosen is None and pprrvu_candidates:
            chosen = pprrvu_candidates[0]
        if chosen is None:
            chosen = ordered[0]
            log.warning(
                "RVU %d: no PPRRVU available, using %s (Addendum B)",
                year,
                chosen["filename"],
            )

        by_year[year] = chosen

    total_rows = 0
    file_count = 0
    for year in sorted(by_year):
        f = by_year[year]
        path = f["path"]
        log.info("Loading RVU year=%d from %s", year, f["filename"])

        try:
            is_pprrvu = _is_pprrvu(f["filename"])
            if is_pprrvu and f["ext"] in (".xlsx", ".xls"):
                df = _read_pprrvu(path)
            elif f["ext"] == ".csv":
                df = _read_csv_with_header(path)
            else:
                df = _read_excel(path)

            if "hcpcs" not in df.columns:
                df = _normalize_columns(df, _RVU_COLUMNS)
            if "hcpcs" not in df.columns:
                log.warning("Skipping %s — no hcpcs column", f["filename"])
                continue

            # Cast numeric columns
            for col in [
                "work_rvu",
                "non_fac_pe_rvu",
                "fac_pe_rvu",
                "mp_rvu",
                "non_fac_total",
                "fac_total",
                "pre_op",
                "intra_op",
                "post_op",
                "conv_factor",
            ]:
                if col in df.columns:
                    df = df.with_columns(pl.col(col).cast(pl.Float64, strict=False))

            # Cast string columns
            for col in df.columns:
                if col not in [
                    "work_rvu",
                    "non_fac_pe_rvu",
                    "fac_pe_rvu",
                    "mp_rvu",
                    "non_fac_total",
                    "fac_total",
                    "pre_op",
                    "intra_op",
                    "post_op",
                    "conv_factor",
                    "year",
                ]:
                    df = df.with_columns(pl.col(col).cast(pl.Utf8, strict=False))

            df = df.with_columns(pl.lit(year).alias("year"))

            # Filter out empty/null hcpcs rows
            df = df.filter(pl.col("hcpcs").is_not_null() & (pl.col("hcpcs") != ""))

            _insert_into(con, "pfs", "rvu", df)
            total_rows += len(df)
            file_count += 1

        except Exception as exc:
            log.error("Failed to load RVU %s: %s", f["filename"], exc)

    return {"rows": total_rows, "files": file_count}


def _is_stacked_header_gpci(path: str, ext: str) -> bool:
    """Detect if a GPCI file has a 2-row stacked header (2017 format).

    2017 files have years in one row and GPCI sub-headers in the next:
      Row 1: Carrier | Locality Number | Locality Name | 2016 | | | 2017 | ...
      Row 2: | | | PW GPCI | PE GPCI | MP GPCI | PW GPCI | PE GPCI | ...

    Later files (2020+) have year-prefixed columns in ONE row:
      Row 2: MAC | State | Locality | Name | 2019 PW GPCI | 2019 PE GPCI | ...

    This function only returns True for the stacked (2-row) format.
    """
    if ext == ".txt":
        enc = _detect_encoding(path)
        with open(path, encoding=enc, errors="replace") as f:
            lines = [l.replace("\x00", "") for l in f.readlines()]
        for i, line in enumerate(lines):
            parts = [p.strip() for p in line.strip().split("\t")]
            # A year row has bare year numbers as values
            bare_years = [p for p in parts if re.match(r"^20\d{2}$", p)]
            if len(bare_years) >= 1:
                # Check next row has GPCI sub-headers
                if i + 1 < len(lines):
                    next_parts = lines[i + 1].strip().split("\t")
                    if any("GPCI" in p for p in next_parts):
                        return True
        return False
    elif ext == ".xlsx":
        import openpyxl

        wb = openpyxl.load_workbook(path, read_only=True)
        ws = wb.active
        rows = []
        for i, row in enumerate(ws.iter_rows(values_only=True)):
            if i > 5:
                break
            rows.append([str(c).strip() if c else "" for c in row])
        wb.close()

        for i, cells in enumerate(rows):
            bare_years = [c for c in cells if re.match(r"^20\d{2}$", c)]
            if len(bare_years) >= 1 and i + 1 < len(rows):
                if any("GPCI" in c for c in rows[i + 1]):
                    return True
    return False


def _parse_stacked_gpci(
    year_parts: list[str],
    sub_parts: list[str],
    data_lines: list[list[str]],
    target_year: int,
) -> pl.DataFrame | None:
    """Core parser for stacked-header GPCI data.

    Handles both multi-year (2017: 3 years in one file) and single-year
    (2018/2019: sub-header omits identifier columns) formats.
    """
    # Find where years start in year_parts
    first_year_idx = None
    for idx, val in enumerate(year_parts):
        if re.match(r"^20\d{2}$", val):
            first_year_idx = idx
            break
    if first_year_idx is None:
        return None

    n_id_cols = first_year_idx

    # Map data column index → year.
    # Extend beyond year_parts length to cover data columns —
    # single-year files have fewer header cells than data cells.
    max_data_cols = max((len(row) for row in data_lines), default=len(year_parts))
    col_year: dict[int, int] = {}
    current_year = 0
    for idx in range(max(len(year_parts), max_data_cols)):
        if idx < len(year_parts) and re.match(r"^20\d{2}$", year_parts[idx]):
            current_year = int(year_parts[idx])
        if current_year and idx >= n_id_cols:
            col_year[idx] = current_year

    # Align sub_parts with data columns.
    # Multi-year: sub_parts has leading blanks matching id cols
    # Single-year: sub_parts omits id cols entirely
    sub_offset = 0
    if len(sub_parts) < len(year_parts):
        sub_offset = n_id_cols

    def _sub_label(data_idx: int) -> str:
        si = data_idx - sub_offset
        if 0 <= si < len(sub_parts):
            return sub_parts[si].upper()
        return ""

    rows = []
    for parts in data_lines:
        if len(parts) < n_id_cols + 3:
            continue
        carrier = parts[0].strip().strip('"')
        locality = parts[1].strip().strip('"') if n_id_cols > 1 else ""
        loc_name = parts[2].strip().strip('"') if n_id_cols > 2 else ""
        if not locality:
            continue

        work = pe = mp = None
        for idx, yr in sorted(col_year.items()):
            if yr == target_year and idx < len(parts):
                label = _sub_label(idx)
                val = parts[idx].strip().strip('"')
                if "PW GPCI" in label or "WORK" in label:
                    work = val
                elif "PE GPCI" in label:
                    pe = val
                elif "MP GPCI" in label:
                    mp = val

        if work is not None:
            rows.append(
                {
                    "mac": carrier,
                    "locality": locality,
                    "locality_name": loc_name,
                    "work_gpci": work,
                    "pe_gpci": pe,
                    "mp_gpci": mp,
                }
            )

    return pl.DataFrame(rows) if rows else None


def _find_stacked_header(
    lines_or_rows: list,
    is_text: bool = True,
) -> tuple[int | None, int | None, list[str], list[str]]:
    """Find year-header and sub-header rows in stacked GPCI data."""
    year_row_idx = None
    sub_header_idx = None
    for i, item in enumerate(lines_or_rows):
        if is_text:
            parts = [p.strip() for p in item.strip().split("\t")]
        else:
            parts = [str(c).strip() if c else "" for c in item]
        years_found = [p for p in parts if re.match(r"^20\d{2}$", p)]
        if len(years_found) >= 1 and year_row_idx is None:
            year_row_idx = i
            continue
        if year_row_idx is not None and any("GPCI" in p for p in parts):
            sub_header_idx = i
            break

    if year_row_idx is None or sub_header_idx is None:
        return None, None, [], []

    if is_text:
        yp = [p.strip() for p in lines_or_rows[year_row_idx].strip().split("\t")]
        sp = [p.strip() for p in lines_or_rows[sub_header_idx].strip().split("\t")]
    else:
        yp = [str(c).strip() if c else "" for c in lines_or_rows[year_row_idx]]
        sp = [str(c).strip() if c else "" for c in lines_or_rows[sub_header_idx]]
    return year_row_idx, sub_header_idx, yp, sp


def _read_multi_year_gpci(
    path: str,
    ext: str,
    target_year: int,
) -> pl.DataFrame | None:
    """Read a stacked-header GPCI file and extract columns for target_year."""
    if ext == ".txt":
        enc = _detect_encoding(path)
        with open(path, encoding=enc, errors="replace") as f:
            lines = [l.replace("\x00", "") for l in f.readlines()]

        _, sub_idx, year_parts, sub_parts = _find_stacked_header(lines, is_text=True)
        if sub_idx is None:
            return None

        data_lines = [
            line.strip().split("\t")
            for line in lines[sub_idx + 1 :]
            if len(line.strip().split("\t")) >= 4
        ]
        return _parse_stacked_gpci(year_parts, sub_parts, data_lines, target_year)

    elif ext == ".xlsx":
        import openpyxl

        wb = openpyxl.load_workbook(path, read_only=True)
        ws = wb.active
        all_rows = list(ws.iter_rows(values_only=True))
        wb.close()

        _, sub_idx, year_parts, sub_parts = _find_stacked_header(
            all_rows, is_text=False
        )
        if sub_idx is None:
            return None

        data_lines = [
            [str(c).strip() if c else "" for c in row]
            for row in all_rows[sub_idx + 1 :]
            if len(row) >= 4
        ]
        return _parse_stacked_gpci(year_parts, sub_parts, data_lines, target_year)

    return None


def _load_gpci_files(
    files: list[dict],
    con: duckdb.DuckDBPyConnection,
) -> dict[str, Any]:
    """Load GPCI/Addendum E files into pfs.gpci."""
    gpci_files = [
        f
        for f in files
        if _is_gpci_file(f["filename"])
        and f["ext"] in (".xlsx", ".xls", ".csv", ".txt")
    ]

    # Deduplicate by year — prefer Final xlsx > Final txt > Proposed
    by_year: dict[int, dict] = {}
    for f in gpci_files:
        year = f["year"]
        is_final = _is_final_rule(f["title"])
        fn = f["filename"].upper()
        is_xlsx = f["ext"] == ".xlsx"
        is_gpci_quarterly = fn.startswith("GPCI")
        priority = (is_gpci_quarterly, is_final, is_xlsx)

        if year not in by_year or priority > by_year[year]["_priority"]:
            f["_priority"] = priority
            by_year[year] = f

    # Check for stacked-header multi-year files that can fill gaps
    multi_year_files = []
    for f in gpci_files:
        if _is_stacked_header_gpci(f["path"], f["ext"]):
            multi_year_files.append(f)

    total_rows = 0
    file_count = 0
    loaded_years: set[int] = set()

    for year in sorted(by_year):
        f = by_year[year]
        path = f["path"]
        log.info("Loading GPCI year=%d from %s", year, f["filename"])

        try:
            # Check if this is a stacked-header multi-year file (2017)
            if _is_stacked_header_gpci(path, f["ext"]):
                df = _read_multi_year_gpci(path, f["ext"], year)
                if df is None or len(df) == 0:
                    log.warning(
                        "Stacked GPCI %s — no data for year %d",
                        f["filename"],
                        year,
                    )
                    continue
            elif f["ext"] == ".txt":
                df = _read_tsv(path)
            elif f["ext"] == ".csv":
                df = _read_csv_with_header(path)
            else:
                df = _read_excel(path)

            # Normalize column names for non-stacked files
            if "work_gpci" not in df.columns:
                rename_map: dict[str, str] = {}
                yr_str = str(year)
                for col in df.columns:
                    # Collapse all whitespace (handles "Work      GPCI ")
                    norm = col.replace("\n", " ").strip()
                    norm = re.sub(r"\s+", " ", norm).upper()

                    if any(k in norm for k in ["MAC", "CONTRACTOR"]):
                        rename_map[col] = "mac"
                    elif norm in ("CARRIER",):
                        rename_map[col] = "mac"
                    elif norm in ("STATE", "ST"):
                        rename_map[col] = "state"
                    elif "LOCALITY" in norm and "NAME" in norm:
                        rename_map[col] = "locality_name"
                    elif "LOCALITY" in norm:
                        rename_map[col] = "locality"
                    elif "WITHOUT" in norm:
                        continue  # skip without-floor variants
                    # Year-prefixed: "2020 PW GPCI ..." — only match target year
                    elif "PW GPCI" in norm or "WORK GPCI" in norm:
                        yr_match = re.match(r"(\d{4})\s", norm)
                        if yr_match and yr_match.group(1) != yr_str:
                            continue
                        if "work_gpci" not in rename_map.values():
                            rename_map[col] = "work_gpci"
                    elif "PE GPCI" in norm:
                        yr_match = re.match(r"(\d{4})\s", norm)
                        if yr_match and yr_match.group(1) != yr_str:
                            continue
                        if "pe_gpci" not in rename_map.values():
                            rename_map[col] = "pe_gpci"
                    elif "MP GPCI" in norm:
                        yr_match = re.match(r"(\d{4})\s", norm)
                        if yr_match and yr_match.group(1) != yr_str:
                            continue
                        if "mp_gpci" not in rename_map.values():
                            rename_map[col] = "mp_gpci"

                if "locality" not in rename_map.values():
                    log.warning(
                        "Skipping GPCI %s — no locality column",
                        f["filename"],
                    )
                    continue

                keep = list(rename_map.keys())
                df = df.select(keep).rename(rename_map)

            for col in ["work_gpci", "pe_gpci", "mp_gpci"]:
                if col in df.columns:
                    df = df.with_columns(pl.col(col).cast(pl.Float64, strict=False))

            for col in ["mac", "state", "locality", "locality_name"]:
                if col in df.columns:
                    df = df.with_columns(pl.col(col).cast(pl.Utf8, strict=False))

            df = df.with_columns(pl.lit(year).alias("year"))
            df = df.filter(
                pl.col("locality").is_not_null() & (pl.col("locality") != "")
            )

            _insert_into(con, "pfs", "gpci", df)
            total_rows += len(df)
            file_count += 1
            loaded_years.add(year)

        except Exception as exc:
            log.error("Failed to load GPCI %s: %s", f["filename"], exc)

    # Fill gaps from multi-year files
    for mf in multi_year_files:
        if mf["ext"] not in (".xlsx", ".txt"):
            continue
        # Detect what years this file contains
        if mf["ext"] == ".txt":
            enc = _detect_encoding(mf["path"])
            with open(mf["path"], encoding=enc, errors="replace") as fh:
                for line in (l.replace("\x00", "") for l in fh):
                    parts = [p.strip() for p in line.strip().split("\t")]
                    avail_years = [int(p) for p in parts if re.match(r"^20\d{2}$", p)]
                    if len(avail_years) >= 1:
                        break
                else:
                    avail_years = []
        else:
            import openpyxl

            wb = openpyxl.load_workbook(mf["path"], read_only=True)
            ws = wb.active
            avail_years = []
            for i, row in enumerate(ws.iter_rows(values_only=True)):
                if i > 5:
                    break
                cells = [str(c).strip() if c else "" for c in row]
                found = [int(c) for c in cells if re.match(r"^20\d{2}$", c)]
                if len(found) >= 1:
                    avail_years = found
                    break
            wb.close()

        for gap_year in avail_years:
            if gap_year in loaded_years:
                continue
            log.info(
                "Filling GPCI gap year=%d from multi-year %s",
                gap_year,
                mf["filename"],
            )
            try:
                df = _read_multi_year_gpci(mf["path"], mf["ext"], gap_year)
                if df is None or len(df) == 0:
                    continue

                for col in ["work_gpci", "pe_gpci", "mp_gpci"]:
                    if col in df.columns:
                        df = df.with_columns(pl.col(col).cast(pl.Float64, strict=False))
                for col in ["mac", "state", "locality", "locality_name"]:
                    if col in df.columns:
                        df = df.with_columns(pl.col(col).cast(pl.Utf8, strict=False))

                df = df.with_columns(pl.lit(gap_year).alias("year"))
                df = df.filter(
                    pl.col("locality").is_not_null() & (pl.col("locality") != "")
                )

                _insert_into(con, "pfs", "gpci", df)
                total_rows += len(df)
                file_count += 1
                loaded_years.add(gap_year)
            except Exception as exc:
                log.error("Failed to fill GPCI gap year=%d: %s", gap_year, exc)

    return {"rows": total_rows, "files": file_count}


def _read_tabular(path: str, ext: str) -> pl.DataFrame:
    """Read a tabular file by extension: xlsx/xls → Excel, csv → CSV, txt → TSV."""
    if ext == ".txt":
        return _read_tsv(path)
    elif ext == ".csv":
        return _read_csv_with_header(path)
    else:
        return _read_excel(path)


def _load_pe_files(
    files: list[dict],
    con: duckdb.DuckDBPyConnection,
) -> dict[str, dict[str, Any]]:
    """Load PUF Labor/Supply/Equipment and Work Time files."""
    stats: dict[str, dict[str, Any]] = {
        "clinical_labor": {"rows": 0, "files": 0},
        "medical_supply": {"rows": 0, "files": 0},
        "medical_equipment": {"rows": 0, "files": 0},
        "physician_work_time": {"rows": 0, "files": 0},
    }

    # Deduplicate per (type, year) — prefer Final over Proposed
    typed_files: dict[str, dict[int, dict]] = {
        "labor": {},
        "supply": {},
        "equipment": {},
        "work_time": {},
    }

    for f in files:
        fn = f["filename"]
        year = f["year"]
        is_final = _is_final_rule(f["title"])
        ext_ok = f["ext"] in (".xlsx", ".xls", ".csv", ".txt")
        if not ext_ok:
            continue

        for ftype, checker in [
            ("labor", _is_puf_labor),
            ("supply", _is_puf_supply),
            ("equipment", _is_puf_equip),
            ("work_time", _is_work_time),
        ]:
            if checker(fn):
                # Priority: Final xlsx > Final txt > Proposed xlsx > ...
                is_xlsx = f["ext"] == ".xlsx"
                priority = (is_final, is_xlsx)
                existing = typed_files[ftype].get(year)
                if existing is None or priority > existing.get(
                    "_priority", (False, False)
                ):
                    f["_final"] = is_final
                    f["_priority"] = priority
                    typed_files[ftype][year] = f
                break

    # Load labor
    for year in sorted(typed_files["labor"]):
        f = typed_files["labor"][year]
        log.info("Loading Labor year=%d from %s", year, f["filename"])
        try:
            df = _read_tabular(f["path"], f["ext"])
            df = _normalize_columns(df, _LABOR_COLUMNS)
            if "hcpcs" not in df.columns:
                log.warning("Skipping Labor %s — no hcpcs", f["filename"])
                continue

            # Compute nf_minutes and f_minutes from service period columns
            nf_min_cols = [
                c for c in df.columns if c.startswith("nf_") and c != "nf_minutes"
            ]
            f_min_cols = [
                c for c in df.columns if c.startswith("f_") and c != "f_minutes"
            ]
            # Cast numeric
            for col in ["rate_per_minute"] + nf_min_cols + f_min_cols:
                if col in df.columns:
                    df = df.with_columns(pl.col(col).cast(pl.Float64, strict=False))

            if nf_min_cols and "nf_minutes" not in df.columns:
                df = df.with_columns(
                    pl.sum_horizontal(*[pl.col(c) for c in nf_min_cols]).alias(
                        "nf_minutes"
                    )
                )
            if f_min_cols and "f_minutes" not in df.columns:
                df = df.with_columns(
                    pl.sum_horizontal(*[pl.col(c) for c in f_min_cols]).alias(
                        "f_minutes"
                    )
                )

            df = df.with_columns(pl.lit(year).alias("year"))
            # Keep only columns matching the table model
            keep = [
                c
                for c in [
                    "hcpcs",
                    "activity_code",
                    "activity_description",
                    "nf_minutes",
                    "f_minutes",
                    "rate_per_minute",
                    "staff_type",
                    "year",
                ]
                if c in df.columns
            ]
            df = df.select(keep)
            df = df.filter(pl.col("hcpcs").is_not_null() & (pl.col("hcpcs") != ""))

            _insert_into(con, "pfs", "clinical_labor", df)
            stats["clinical_labor"]["rows"] += len(df)
            stats["clinical_labor"]["files"] += 1
        except Exception as exc:
            log.error("Failed to load Labor %s: %s", f["filename"], exc)

    # Load supply
    for year in sorted(typed_files["supply"]):
        f = typed_files["supply"][year]
        log.info("Loading Supply year=%d from %s", year, f["filename"])
        try:
            df = _read_tabular(f["path"], f["ext"])
            df = _normalize_columns(df, _SUPPLY_COLUMNS)
            if "hcpcs" not in df.columns:
                log.warning("Skipping Supply %s — no hcpcs", f["filename"])
                continue

            for col in ["unit_price", "nf_quantity", "f_quantity"]:
                if col in df.columns:
                    df = df.with_columns(pl.col(col).cast(pl.Float64, strict=False))

            df = df.with_columns(pl.lit(year).alias("year"))
            keep = [
                c
                for c in [
                    "hcpcs",
                    "supply_code",
                    "supply_description",
                    "nf_quantity",
                    "f_quantity",
                    "unit_price",
                    "year",
                ]
                if c in df.columns
            ]
            df = df.select(keep)
            df = df.filter(pl.col("hcpcs").is_not_null() & (pl.col("hcpcs") != ""))

            _insert_into(con, "pfs", "medical_supply", df)
            stats["medical_supply"]["rows"] += len(df)
            stats["medical_supply"]["files"] += 1
        except Exception as exc:
            log.error("Failed to load Supply %s: %s", f["filename"], exc)

    # Load equipment
    for year in sorted(typed_files["equipment"]):
        f = typed_files["equipment"][year]
        log.info("Loading Equipment year=%d from %s", year, f["filename"])
        try:
            df = _read_tabular(f["path"], f["ext"])
            df = _normalize_columns(df, _EQUIPMENT_COLUMNS)
            if "hcpcs" not in df.columns:
                log.warning("Skipping Equipment %s — no hcpcs", f["filename"])
                continue

            for col in [
                "unit_price",
                "useful_life",
                "minutes_per_year",
                "nf_minutes",
                "f_minutes",
                "maintenance_factor",
            ]:
                if col in df.columns:
                    df = df.with_columns(pl.col(col).cast(pl.Float64, strict=False))

            df = df.with_columns(pl.lit(year).alias("year"))
            keep = [
                c
                for c in [
                    "hcpcs",
                    "equipment_code",
                    "equipment_description",
                    "nf_minutes",
                    "f_minutes",
                    "unit_price",
                    "useful_life",
                    "minutes_per_year",
                    "maintenance_factor",
                    "year",
                ]
                if c in df.columns
            ]
            df = df.select(keep)
            df = df.filter(pl.col("hcpcs").is_not_null() & (pl.col("hcpcs") != ""))

            _insert_into(con, "pfs", "medical_equipment", df)
            stats["medical_equipment"]["rows"] += len(df)
            stats["medical_equipment"]["files"] += 1
        except Exception as exc:
            log.error("Failed to load Equipment %s: %s", f["filename"], exc)

    # Load work time
    for year in sorted(typed_files["work_time"]):
        f = typed_files["work_time"][year]
        log.info("Loading WorkTime year=%d from %s", year, f["filename"])
        try:
            df = _read_tabular(f["path"], f["ext"])
            df = _normalize_columns(df, _WORK_TIME_COLUMNS)
            if "hcpcs" not in df.columns:
                log.warning("Skipping WorkTime %s — no hcpcs", f["filename"])
                continue

            # hcpcs might be int (like 100 instead of "00100") — pad
            if df["hcpcs"].dtype in (pl.Int64, pl.Float64, pl.Int32):
                df = df.with_columns(
                    pl.col("hcpcs")
                    .cast(pl.Int64, strict=False)
                    .cast(pl.Utf8)
                    .str.zfill(5)
                    .alias("hcpcs")
                )

            for col in [
                "pre_eval_time",
                "pre_positioning_time",
                "pre_scrub_dress_wait",
                "intra_service_time",
                "post_immediate_time",
                "total_time",
            ]:
                if col in df.columns:
                    df = df.with_columns(pl.col(col).cast(pl.Float64, strict=False))

            # Compute total_time if not present
            time_cols = [
                c
                for c in [
                    "pre_eval_time",
                    "pre_positioning_time",
                    "pre_scrub_dress_wait",
                    "intra_service_time",
                    "post_immediate_time",
                ]
                if c in df.columns
            ]
            if time_cols and "total_time" not in df.columns:
                df = df.with_columns(
                    pl.sum_horizontal(*[pl.col(c) for c in time_cols]).alias(
                        "total_time"
                    )
                )

            df = df.with_columns(pl.lit(year).alias("year"))
            keep = [
                c
                for c in [
                    "hcpcs",
                    "pre_eval_time",
                    "pre_positioning_time",
                    "pre_scrub_dress_wait",
                    "intra_service_time",
                    "post_immediate_time",
                    "total_time",
                    "year",
                ]
                if c in df.columns
            ]
            df = df.select(keep)
            df = df.filter(pl.col("hcpcs").is_not_null() & (pl.col("hcpcs") != ""))

            _insert_into(con, "pfs", "physician_work_time", df)
            stats["physician_work_time"]["rows"] += len(df)
            stats["physician_work_time"]["files"] += 1
        except Exception as exc:
            log.error("Failed to load WorkTime %s: %s", f["filename"], exc)

    return stats


# ── ZIP code carrier locality loader ─────────────────────────

# Fixed-width column spec for ZIP5 files (from ZIP5lyout.txt)
_ZIP5_COLUMNS = [
    ("state", 0, 2),
    ("zip_code", 2, 7),
    ("carrier", 7, 12),
    ("locality", 12, 14),
    ("rural_indicator", 14, 15),
    # 15-17: lab CB locality (skip)
    # 17: rural indicator 2 (skip)
    # 18-19: filler
    ("plus_four_flag", 20, 21),
    # 21: filler
    ("part_b_indicator", 22, 23),
    # 23-74: filler
    ("year_quarter", 75, 80),
]


def _load_zip_carrier_files(
    files: list[dict],
    con: duckdb.DuckDBPyConnection,
) -> dict[str, Any]:
    """Load ZIP5 carrier locality files into pfs.zip_carrier_locality."""
    zip_files = [
        f for f in files if _is_zip_carrier_file(f["filename"]) and f["ext"] == ".txt"
    ]

    # Deduplicate by year — prefer latest version
    by_year: dict[int, dict] = {}
    for f in zip_files:
        year = f["year"]
        existing = by_year.get(year)
        if existing is None or f["filename"] > existing["filename"]:
            by_year[year] = f

    total_rows = 0
    file_count = 0

    for year in sorted(by_year):
        f = by_year[year]
        path = f["path"]
        log.info("Loading ZipCarrier year=%d from %s", year, f["filename"])

        try:
            with open(path, encoding="utf-8-sig") as fh:
                lines = fh.readlines()

            rows = []
            for line in lines:
                if len(line.rstrip()) < 14:
                    continue
                row = {}
                for col_name, start, end in _ZIP5_COLUMNS:
                    row[col_name] = line[start:end].strip()
                rows.append(row)

            if not rows:
                log.warning("Skipping %s — no data rows", f["filename"])
                continue

            df = pl.DataFrame(rows)

            # Parse year/quarter from the YYYYQ field
            if "year_quarter" in df.columns:
                df = df.with_columns(
                    pl.col("year_quarter")
                    .str.slice(0, 4)
                    .cast(pl.Int64, strict=False)
                    .alias("year"),
                    pl.col("year_quarter")
                    .str.slice(4, 1)
                    .cast(pl.Int64, strict=False)
                    .alias("quarter"),
                )
                df = df.drop("year_quarter")
            else:
                df = df.with_columns(
                    pl.lit(year).alias("year"),
                    pl.lit(None, dtype=pl.Int64).alias("quarter"),
                )

            # Filter empty zip codes
            df = df.filter(
                pl.col("zip_code").is_not_null() & (pl.col("zip_code") != "")
            )

            _insert_into(con, "pfs", "zip_carrier_locality", df)
            total_rows += len(df)
            file_count += 1
            log.info("  ZipCarrier year=%d: %d rows", year, len(df))

        except Exception as exc:
            log.error("Failed to load ZipCarrier %s: %s", f["filename"], exc)

    return {"rows": total_rows, "files": file_count}


# ── Carrier file loader ──────────────────────────────────────

# Positional columns for the 16-column format (2011+)
_CARRIER_POS_16 = [
    "year",
    "mac",
    "locality",
    "hcpcs",
    "mod",
    "non_fac_fee",
    "fac_fee",
    "_filler",
    "_pctc",
    "_status",
    "_mult_surg",
    "_therapy_nf",
    "_therapy_f",
    "_opps_ind",
    "_opps_nf",
    "_opps_f",
]

# Positional columns for the 8-column format (2002-2009)
_CARRIER_POS_8 = [
    "year",
    "mac",
    "locality",
    "hcpcs",
    "mod",
    "non_fac_fee",
    "fac_fee",
    "_opps_ind",
]

# Limiting-charge factor (0.95 × 1.15 = 1.0925) is defined in
# pfs.calcs.payment.LIMITING_CHARGE_FACTOR — single source of truth.


def _load_carrier_files(
    files: list[dict],
    con: duckdb.DuckDBPyConnection,
) -> dict[str, Any]:
    """Load carrier fee schedule files into pfs.carrier_locality."""
    carrier_files = [
        f for f in files if _is_carrier_file(f["filename"]) and f["ext"] == ".txt"
    ]

    # Group by year
    by_year: dict[int, list[dict]] = {}
    for f in carrier_files:
        by_year.setdefault(f["year"], []).append(f)

    # For 2026, prefer NonQP over QP (standard payment amounts)
    for year, year_files in by_year.items():
        titles = {f["title"] for f in year_files}
        if any("Nonqp" in t or "NONQP" in t.upper() for t in titles):
            # Keep only files from the NonQP item
            by_year[year] = [
                f
                for f in year_files
                if "Nonqp" in f["title"] or "NONQP" in f["title"].upper()
            ]

    total_rows = 0
    file_count = 0

    for year in sorted(by_year):
        year_files = by_year[year]
        log.info("Loading Carrier year=%d from %d file(s)", year, len(year_files))

        year_frames: list[pl.DataFrame] = []
        for f in year_files:
            try:
                df = pl.read_csv(
                    f["path"],
                    has_header=False,
                    infer_schema_length=0,
                    truncate_ragged_lines=True,
                )
            except Exception as exc:
                log.warning("Skipping carrier file %s: %s", f["filename"], exc)
                continue

            n_cols = len(df.columns)
            if n_cols >= 16:
                pos = _CARRIER_POS_16
            elif n_cols >= 8:
                pos = _CARRIER_POS_8
            else:
                log.warning("Skipping %s — only %d columns", f["filename"], n_cols)
                continue

            rename = {old: new for old, new in zip(df.columns[: len(pos)], pos)}
            df = df.select(list(rename.keys())).rename(rename)

            # Drop internal columns
            internal = [c for c in df.columns if c.startswith("_")]
            if internal:
                df = df.drop(internal)

            # Strip quotes and whitespace from string fields
            for col in ["mac", "locality", "hcpcs", "mod"]:
                if col in df.columns:
                    df = df.with_columns(pl.col(col).str.strip_chars('" ').alias(col))

            # Cast fee columns to float
            for col in ["non_fac_fee", "fac_fee"]:
                if col in df.columns:
                    df = df.with_columns(
                        pl.col(col)
                        .str.strip_chars('" ')
                        .cast(pl.Float64, strict=False)
                        .alias(col)
                    )

            # Compute limiting charges (PAR fee × 0.95 × 1.15 = 1.0925).
            # Constant lives in pfs.calcs.payment for a single source of truth.
            if "non_fac_fee" in df.columns:
                df = df.with_columns(
                    (pl.col("non_fac_fee") * LIMITING_CHARGE_FACTOR)
                    .round(2)
                    .alias("non_fac_limiting_charge")
                )
            if "fac_fee" in df.columns:
                df = df.with_columns(
                    (pl.col("fac_fee") * LIMITING_CHARGE_FACTOR)
                    .round(2)
                    .alias("fac_limiting_charge")
                )

            # Override year from metadata (more reliable than file content)
            df = df.with_columns(pl.lit(year).alias("year"))

            # Filter empty HCPCS rows
            df = df.filter(pl.col("hcpcs").is_not_null() & (pl.col("hcpcs") != ""))

            year_frames.append(df)

        if not year_frames:
            continue

        combined = pl.concat(year_frames)
        _insert_into(con, "pfs", "carrier_locality", combined)
        total_rows += len(combined)
        file_count += 1
        log.info(
            "  Carrier year=%d: %d rows from %d file(s)",
            year,
            len(combined),
            len(year_frames),
        )

    return {"rows": total_rows, "files": file_count}


# ── DuckDB insertion ─────────────────────────────────────────


def _insert_into(
    con: duckdb.DuckDBPyConnection,
    schema: str,
    table: str,
    df: pl.DataFrame,
) -> None:
    """INSERT INTO schema.table from a polars DataFrame."""
    qualified = f'"{schema}"."{table}"'
    # Ensure schema exists
    con.execute(f'CREATE SCHEMA IF NOT EXISTS "{schema}"')
    # Create table if not exists from the DataFrame
    try:
        con.execute(f"SELECT 1 FROM {qualified} LIMIT 0")
    except duckdb.CatalogException:
        con.execute(f"CREATE TABLE {qualified} AS SELECT * FROM df WHERE 1=0")
    # Insert — match columns by name
    cols = df.columns
    # Only insert columns that exist in target
    target_cols = [
        r[0]
        for r in con.execute(
            f"SELECT column_name FROM information_schema.columns "
            f"WHERE table_schema = '{schema}' AND table_name = '{table}'"
        ).fetchall()
    ]
    shared = [c for c in cols if c in target_cols]
    if not shared:
        # Table doesn't have overlapping cols — create fresh
        con.execute(f"DROP TABLE IF EXISTS {qualified}")
        con.execute(f"CREATE TABLE {qualified} AS SELECT * FROM df")
        return

    missing_in_target = [c for c in cols if c not in target_cols]
    if missing_in_target:
        for mc in missing_in_target:
            dtype = df[mc].dtype
            sql_type = "VARCHAR"
            if dtype in (pl.Float64, pl.Float32):
                sql_type = "DOUBLE"
            elif dtype in (pl.Int64, pl.Int32, pl.Int16, pl.Int8):
                sql_type = "BIGINT"
            try:
                con.execute(f'ALTER TABLE {qualified} ADD COLUMN "{mc}" {sql_type}')
            except Exception:  # pragma: no cover
                pass  # Column may already exist

    col_list = ", ".join(f'"{c}"' for c in cols)
    con.execute(f"INSERT INTO {qualified} ({col_list}) SELECT {col_list} FROM df")


# ── Main entry point ─────────────────────────────────────────


def load_all(
    con: duckdb.DuckDBPyConnection,
    *,
    zotero_db: str | Path | None = None,
    zotero_storage: str | Path | None = None,
    years: tuple[int, int] | None = None,
) -> dict[str, dict[str, Any]]:
    """Load all PFS data from Zotero into DuckDB tables.

    Parameters
    ----------
    con : duckdb.DuckDBPyConnection
        DuckDB connection to load into.
    zotero_db : str | Path, optional
        Path to Zotero SQLite database.  Defaults to ``db.zotero``
        from ``stack.toml``.
    zotero_storage : str | Path, optional
        Path to Zotero storage directory.  Defaults to
        ``storage.zotero`` from ``stack.toml``.
    years : tuple[int, int] | None
        Optional (start, end) year range filter (inclusive).

    Returns
    -------
    dict[str, dict]
        Per-table stats: ``{"pfs.rvu": {"rows": N, "files": N}, ...}``
    """
    if zotero_db is None or zotero_storage is None:
        from conf import path

        if zotero_db is None:
            zotero_db = path("db.zotero")
        if zotero_storage is None:
            zotero_storage = path("storage.zotero")
    log.info("Discovering PFS files from Zotero...")

    all_files = _zotero_pfs_files(zotero_db, zotero_storage)

    # Filter by year range
    if years:
        lo, hi = years
        all_files = [f for f in all_files if lo <= f["year"] <= hi]

    # Classify files by filename pattern
    rvu_files = [
        f
        for f in all_files
        if _is_addendum_b(f["filename"]) or _is_pprrvu(f["filename"])
    ]
    gpci_files = [f for f in all_files if _is_gpci_file(f["filename"])]
    pe_files = [
        f
        for f in all_files
        if (
            _is_puf_labor(f["filename"])
            or _is_puf_supply(f["filename"])
            or _is_puf_equip(f["filename"])
            or _is_work_time(f["filename"])
        )
    ]
    carrier_files = [f for f in all_files if _is_carrier_file(f["filename"])]
    zip_carrier_files = [f for f in all_files if _is_zip_carrier_file(f["filename"])]

    log.info(
        "Found %d RVU, %d GPCI, %d PE, %d Carrier, %d ZipCarrier "
        "attachment(s) from %d total",
        len(rvu_files),
        len(gpci_files),
        len(pe_files),
        len(carrier_files),
        len(zip_carrier_files),
        len(all_files),
    )

    con.execute('CREATE SCHEMA IF NOT EXISTS "pfs"')

    # Drop existing tables for clean load
    for tbl in [
        "rvu",
        "gpci",
        "clinical_labor",
        "medical_supply",
        "medical_equipment",
        "physician_work_time",
        "carrier_locality",
        "zip_carrier_locality",
    ]:
        con.execute(f'DROP TABLE IF EXISTS "pfs"."{tbl}"')

    stats: dict[str, dict[str, Any]] = {}

    # Load RVU files
    stats["pfs.rvu"] = _load_rvu_files(rvu_files, con)

    # Load GPCI files
    stats["pfs.gpci"] = _load_gpci_files(gpci_files, con)

    # Load PE files (labor, supply, equipment, work time)
    pe_stats = _load_pe_files(pe_files, con)
    for table_name, table_stats in pe_stats.items():
        stats[f"pfs.{table_name}"] = table_stats

    # Load Carrier files
    stats["pfs.carrier_locality"] = _load_carrier_files(carrier_files, con)

    # Load ZIP code carrier locality files
    stats["pfs.zip_carrier_locality"] = _load_zip_carrier_files(zip_carrier_files, con)

    return stats