"""03_Simulated - Jeddah rate card, simulated monthly P&L, cost per case / pallet and Monte Carlo (version 3, 22 Oct 2026).

Version 3: labour cost is on the HEAD-COUNT basis (sim_labour 6b): payroll-confirmed people are paid for every working day
from their first to their last scan, so the paid days on which they did not scan are a cost too. That cost is not tied to
any task, so it is spread over the service families pro rata, like warehouse overhead. "Worked-time" figures (days with
scans only, as in version 2) are kept next to the head-count ones. Scenarios: Conservative = scan days only, Base =
payroll-confirmed head count, Stretch = every employee id salaried; the Monte Carlo draws the basis 25 / 50 / 25%.

Revenue side now follows Binzagr's own transaction mapping (received 22 Oct 2026): only the activities Binzagr marks
"Revenue Generation" are billed - inbound receipt, damaged-goods receipt, unload / unpick, case pick (per CASE, using
Binzagr's conversion factors), pallet pick, consolidation, re-palletisation, pallet split, returns put-away.
Everything else (put-away, picking scans, loading, shipping, billing echo rows, moves, counts, QC ...) is cost that the
billed activities must recover.

Rate card (Simulated, cost-recovery): the simulated labour cost, posted to activities by task time (sim_labour step 7),
is grouped into five service families (revenue activity + the cost activities that support it); warehouse overhead
(in-house moves, counts, QC, adjustments) is spread pro rata. Each family's cost is split over its billed activities by
number of billed lines, then grossed up for non-labour opex (labour share of opex) and a benchmark 3PL EBIT margin.
By construction the BASELINE P&L earns the benchmark margin - the value of the model is the CHANGE when recoverable
labour is removed. Binzagr's actual P&L and client rates replace this when received (04_Gaps).

Outputs: sim_rate_card, sim_jeddah_pnl_monthly, sim_jeddah_pnl_by_wh_month, sim_cost_per_unit_wh_month, sim_scenarios,
         sim_monte_carlo_draws, sim_monte_carlo_summary
"""
import glob, os
import numpy as np
import pandas as pd
from sim_common import DATA, AUDIT, CLEAN, FX_USD_SAR, write_table, save_stats, load_stats
from parameters import tables, get, pdict
from sim_labour import loaded_cost, TIERS

FAMILIES = {
    'Inbound handling': {'revenue': {'Inbound receipt': (153, 'pallet received', 'events'),
                                     'Damaged-goods receipt': (183, 'damaged pallet received', 'events'),
                                     'Unload / unpick': (391, 'pallet unloaded / unpicked', 'events')},
                         'support': [154, 155, 156, 157, 158, 159, 160, 184, 211, 212, 392, 393, 914]},
    'Outbound handling': {'revenue': {'Case pick': (400, 'case picked', 'cases'),
                                      'Pallet pick': (410, 'pallet picked', 'events')},
                          'support': [301, 302, 305, 306, 317, 318, 321, 322, 323, 330, 331, 340, 341, 343, 344, 345, 346, 347]},
    'Consolidation': {'revenue': {'Consolidation': (350, 'line consolidated', 'events')}, 'support': [351]},
    'Re-palletisation / split': {'revenue': {'Re-palletisation': (204, 'pallet re-palletised', 'events'),
                                             'Pallet split': (206, 'pallet split', 'events')}, 'support': [203, 205]},
    'Returns': {'revenue': {'Returns put-away': (172, 'return put away', 'events')}, 'support': [171]},
}
OVERHEAD = 'Warehouse overhead (in-house / control)'
SERVICES = {svc: (fam, code, unit, basis) for fam, f in FAMILIES.items() for svc, (code, unit, basis) in f['revenue'].items()}
MONTHS_OBSERVED = 8


def tri(rng, lo, mode, hi, n):
    lo, hi = min(lo, hi), max(lo, hi)
    if hi - lo < 1e-12:
        return np.full(n, mode)
    return rng.triangular(lo, min(max(mode, lo), hi), hi, n)


def family_of(code):
    for fam, f in FAMILIES.items():
        if code in f['support'] or code in [c for c, _, _ in f['revenue'].values()]:
            return fam
    return OVERHEAD


def units_frame(tx, keys):
    """Billable units (Actual) per key: events of the revenue code, or cases for case pick."""
    rows = {}
    for svc, (fam, code, unit, basis) in SERVICES.items():
        x = tx[tx.tran_type == code]
        rows[svc] = (x.groupby(keys).qty_cases.sum() if basis == 'cases' else x.groupby(keys).size())
    return pd.DataFrame(rows).fillna(0)


def main():
    src, par = tables()
    L = load_stats('labour')
    vol = pd.read_parquet(os.path.join(DATA, 'sim_activity_volumes.parquet'))
    sess = pd.read_parquet(os.path.join(DATA, 'sim_session_attendance.parquet'))
    pm = pd.read_parquet(os.path.join(DATA, 'sim_person_month_labour.parquet'))
    hcm = pd.read_parquet(os.path.join(DATA, 'inf_jeddah_headcount.parquet')).set_index('month')
    pmb = pm[pm.in_headcount_base]
    noscan_m = pmb.groupby('month').labour_cost_without_scan_sar.sum()                 # base cost parameters
    noscan_whm = pmb.groupby(['main_wh', 'month']).labour_cost_without_scan_sar.sum()
    noscan_8m = float(noscan_m.sum())
    vol['family'] = vol.tran_type.map(family_of)
    st = {'families': {f: sorted(int(c) for c in vol[vol.family == f].tran_type.unique()) for f in list(FAMILIES) + [OVERHEAD]}}

    # ---------------- billable units (Actual, all logins incl. system ones) from 02_Cleaned
    tx = pd.concat([pd.read_parquet(p, columns=['tran_type', 'wh_id', 'month', 'qty_cases', 'employee_id_norm'])
                    for p in sorted(glob.glob(os.path.join(CLEAN, 'jeddah_wms_transactions', '*.parquet')))], ignore_index=True)
    rev_codes = [code for (_, code, _, _) in SERVICES.values()]
    tx = tx[tx.tran_type.isin(rev_codes)].copy()
    tx['wh_id'] = tx.wh_id.astype(str)
    st['revenue_rows'] = int(len(tx))
    st['revenue_rows_by_system_logins'] = int(tx.employee_id_norm.isin(['SAPIMPORT', 'ADMIN', 'SHORTPICK']).sum())
    cp = tx[tx.tran_type == 400]
    st['case_pick_rows_with_case_factor_share'] = float(cp.qty_cases.notna().mean())
    u_tot = units_frame(tx.assign(_all='all'), '_all').iloc[0]
    st['units_8m'] = {k: float(v) for k, v in u_tot.items()}

    # ---------------- cost by family (base posting) -> service cost shares
    fam_cost = vol.groupby('family').posted_labour_cost_sar.sum()
    direct = fam_cost.drop(OVERHEAD, errors='ignore')
    overhead = float(fam_cost.get(OVERHEAD, 0))
    fam_total = direct + (overhead + noscan_8m) * direct / direct.sum()          # head-count basis: paid days without a scan spread pro rata
    fam_total_w = direct + overhead * direct / direct.sum()                        # worked time only (days with scans; version 2)
    lines = {svc: float(tx[tx.tran_type == code].shape[0]) for svc, (_, code, _, _) in SERVICES.items()}
    svc_cost, svc_cost_w = {}, {}
    for fam, f in FAMILIES.items():
        n = sum(lines[s] for s in f['revenue'])
        for s in f['revenue']:
            svc_cost[s] = float(fam_total.get(fam, 0)) * (lines[s] / n if n else 0)
            svc_cost_w[s] = float(fam_total_w.get(fam, 0)) * (lines[s] / n if n else 0)
    labour_scan_8m = float(sess.labour_cost_base_sar.sum())
    labour_base_8m = labour_scan_8m + noscan_8m
    svc_share = {s: v / labour_base_8m for s, v in svc_cost.items()}
    st['family_cost_base_8m'] = {k: float(v) for k, v in fam_cost.items()}
    st['family_total_incl_overhead_base_8m'] = {k: float(v) for k, v in fam_total.items()}
    st['overhead_share_of_labour'] = overhead / labour_base_8m
    st['paid_days_without_scan_share_of_labour'] = noscan_8m / labour_base_8m
    st['labour_8m_by_basis_base_costs'] = {'scan_days': labour_scan_8m, 'headcount_base': labour_base_8m}
    st['service_cost_share'] = svc_share

    # ---------------- rate card (low / base / high)
    def labour_total(pv, basis='headcount_base'):
        c = loaded_cost(pv)
        shift = sess.cost_tier.map({t: c[t]['cost_per_worked_shift'] for t in TIERS})
        ot = sess.cost_tier.map({t: c[t]['overtime_cost_per_hour'] for t in TIERS})
        x = pm if basis == 'all_ids' else (pmb if basis == 'headcount_base' else pm.iloc[:0])
        return float((sess.paid_shifts * shift + sess.overtime_h * ot).sum() + (x.paid_days_without_scan * x.cost_tier.map(
            {t: c[t]['cost_per_worked_shift'] for t in TIERS})).sum())
    rc_rows = []
    lab_scen = {}
    for scen in ('low', 'base', 'high'):
        pv = pdict(par, scen)
        lab = labour_total(pv); lab_scen[scen] = lab
        gross = 1 / pv['labour_share_of_warehouse_opex'] / (1 - pv['ebit_margin_3pl'])
        for s, (fam, code, unit, basis) in SERVICES.items():
            u = float(u_tot[s]); cost = lab * svc_share[s]
            rc_rows.append({'scenario': scen, 'service': s, 'labour_cost_per_unit_sar': cost / u if u else np.nan,
                            'worked_time_labour_per_unit_sar': svc_cost_w[s] / labour_scan_8m * labour_total(pv, 'scan_days') / u if u else np.nan,
                            'rate_sar_per_unit': cost * gross / u if u else np.nan})
    rc = pd.DataFrame(rc_rows)
    rcw = rc.pivot(index='service', columns='scenario', values='rate_sar_per_unit').reset_index()
    rcw = rcw.rename(columns={'low': 'rate_low_sar', 'base': 'rate_base_sar', 'high': 'rate_high_sar'})
    rcw = rcw.merge(rc[rc.scenario == 'base'][['service', 'labour_cost_per_unit_sar', 'worked_time_labour_per_unit_sar']], on='service')
    rcw['family'] = rcw.service.map(lambda s: SERVICES[s][0])
    rcw['tran_type'] = rcw.service.map(lambda s: SERVICES[s][1])
    rcw['unit'] = rcw.service.map(lambda s: SERVICES[s][2])
    rcw['units_8_months'] = rcw.service.map(lambda s: float(u_tot[s]))
    rcw['share_of_labour_cost'] = rcw.service.map(svc_share)
    rcw['binzagr_movement_type'] = rcw.tran_type.map(vol.groupby('tran_type').movement_type.first())
    # external checks: receiving (US survey conversion) and case pick (Binzagr's own example rates)
    pvb = pdict(par, 'base'); cb = loaded_cost(pvb)
    ksa_bill_rate = cb['operative']['cost_per_paid_hour'] / pvb['labour_share_of_warehouse_opex'] / (1 - pvb['ebit_margin_3pl'])
    us_equiv_h = get(par, 'us_receiving_fee_per_pallet') / get(par, 'us_receiving_labour_rate')
    lo_c, b_c, hi_c = (get(par, 'client_case_rate_example', s) for s in ('low', 'base', 'high'))
    rcw['external_check_sar'] = np.select([rcw.service.eq('Inbound receipt'), rcw.service.eq('Case pick')], [us_equiv_h * ksa_bill_rate, b_c], np.nan)
    rcw['external_check_note'] = np.select(
        [rcw.service.eq('Inbound receipt'), rcw.service.eq('Case pick')],
        [f'US survey receiving fee / US labour rate = {us_equiv_h:.3f} billing-hours per pallet x KSA billing rate SAR {ksa_bill_rate:.1f}/h (S17)',
         f'Binzagr example: the same case activity is billed SAR {lo_c:.2f}-{hi_c:.2f} depending on the client (S39)'], '')
    rcw['caveat'] = np.select([rcw.service.isin(['Re-palletisation', 'Pallet split', 'Damaged-goods receipt', 'Unload / unpick'])],
                              ['Small volume, mostly manual work the scans see poorly - replace with the client rate card'], '')
    rcw['basis'] = ('Simulated (cost-recovery on head-count-basis labour: posted to activities by task time, grouped by service family, overhead '
                    'and paid days without a scan pro rata, split by billed lines, grossed up for non-labour opex and margin; worked-time labour = '
                    'days with scans only); units = Actual')
    rcw = rcw[['service', 'family', 'tran_type', 'unit', 'units_8_months', 'share_of_labour_cost', 'labour_cost_per_unit_sar',
               'worked_time_labour_per_unit_sar', 'rate_low_sar',
               'rate_base_sar', 'rate_high_sar', 'external_check_sar', 'external_check_note', 'binzagr_movement_type', 'caveat', 'basis']]
    st['ksa_billing_labour_rate_base'] = ksa_bill_rate
    st['labour_total_8m'] = lab_scen
    cpk = rcw.set_index('service').loc['Case pick']
    st['case_pick_rate_base'] = float(cpk.rate_base_sar)
    st['case_pick_rate_vs_client_example'] = {'low': lo_c, 'high': hi_c, 'inside': bool(lo_c <= cpk.rate_base_sar <= hi_c)}

    # ---------------- monthly P&L (base) - Jeddah total and per warehouse
    rate_base = rcw.set_index('service').rate_base_sar
    pv = pdict(par, 'base')
    ls = pv['labour_share_of_warehouse_opex']
    u_m = units_frame(tx, 'month')
    rev_m = (u_m * rate_base).sum(axis=1)
    lab_scan_m = sess.groupby('month').labour_cost_base_sar.sum()
    lab_m = lab_scan_m.add(noscan_m, fill_value=0)
    rec_m = sess.groupby('month').recoverable_cost_base_sar.sum()
    P = pd.DataFrame({'month': rev_m.index, 'handling_revenue_sar': rev_m.values})
    P = P.merge(u_m.add_prefix('units_').reset_index().rename(columns={'index': 'month'}), on='month')
    P['labour_cost_sar'] = P.month.map(lab_m)
    P['labour_cost_scan_days_sar'] = P.month.map(lab_scan_m)
    P['labour_cost_paid_days_without_scan_sar'] = P.month.map(noscan_m).fillna(0)
    P['other_opex_sar'] = P.labour_cost_sar * (1 / ls - 1)
    P['ebit_sar'] = P.handling_revenue_sar - P.labour_cost_sar - P.other_opex_sar
    P['recoverable_labour_sar'] = P.month.map(rec_m)
    P['labour_cost_optimised_sar'] = P.labour_cost_sar - P.recoverable_labour_sar
    P['ebit_optimised_sar'] = P.ebit_sar + P.recoverable_labour_sar
    P['ebit_margin'] = P.ebit_sar / P.handling_revenue_sar
    P['ebit_margin_optimised'] = P.ebit_optimised_sar / P.handling_revenue_sar
    P['labour_pct_of_revenue'] = P.labour_cost_sar / P.handling_revenue_sar
    P['labour_pct_of_revenue_optimised'] = P.labour_cost_optimised_sar / P.handling_revenue_sar
    fte = sess.groupby('month').paid_shifts.sum() / pv['working_days_per_month']
    P['fte_equivalent'] = P.month.map(hcm.fte_headcount_base)                      # head count (base): payroll-confirmed people + other logins' scan days
    P['fte_scan_days'] = P.month.map(fte)
    P['employee_ids_scanning'] = P.month.map(hcm.employee_ids_scanning)
    P['registered_active_people'] = P.month.map(hcm.registered_active_people_log_sites)
    P['fte_recoverable'] = P.month.map(sess.groupby('month').apply(lambda d: (d.recoverable_h_base / d.legal_hours).sum(), include_groups=False)) / pv['working_days_per_month']
    P['basis'] = 'units = Actual; revenue, costs, EBIT = Simulated (base scenario, head-count basis); head counts = Actual / Inferred'
    tot = P[[c for c in P.columns if c.endswith('_sar') or c.startswith('units_')]].sum()
    tot['ebit_margin'] = tot.ebit_sar / tot.handling_revenue_sar
    tot['ebit_margin_optimised'] = tot.ebit_optimised_sar / tot.handling_revenue_sar
    tot['labour_pct_of_revenue'] = tot.labour_cost_sar / tot.handling_revenue_sar
    tot['labour_pct_of_revenue_optimised'] = tot.labour_cost_optimised_sar / tot.handling_revenue_sar
    tot['fte_equivalent_avg'] = P.fte_equivalent.mean(); tot['fte_recoverable_avg'] = P.fte_recoverable.mean()
    tot['fte_scan_days_avg'] = P.fte_scan_days.mean(); tot['employee_ids_scanning_avg'] = P.employee_ids_scanning.mean()
    st['pnl_8m_base'] = {k: float(v) for k, v in tot.items()}
    st['pnl_annualised_base'] = {k: float(v) * 12 / MONTHS_OBSERVED for k, v in tot.items() if k.endswith('_sar')}

    uw = units_frame(tx, ['wh_id', 'month']).reset_index()
    uw['handling_revenue_sar'] = (uw[list(SERVICES)] * rate_base).sum(axis=1)
    lw = sess.groupby(['main_wh', 'month']).agg(labour_cost_sar=('labour_cost_base_sar', 'sum'), recoverable_labour_sar=('recoverable_cost_base_sar', 'sum'),
                                                paid_shifts=('paid_shifts', 'sum')).reset_index().rename(columns={'main_wh': 'wh_id'})
    nsw = noscan_whm.rename('labour_cost_paid_days_without_scan_sar').reset_index().rename(columns={'main_wh': 'wh_id'})
    lw = lw.rename(columns={'labour_cost_sar': 'labour_cost_scan_days_sar'}).merge(nsw, on=['wh_id', 'month'], how='outer').fillna(0)
    lw['labour_cost_sar'] = lw.labour_cost_scan_days_sar + lw.labour_cost_paid_days_without_scan_sar
    W = uw.merge(lw, on=['wh_id', 'month'], how='outer').fillna(0)
    W['labour_pct_of_revenue'] = np.where(W.handling_revenue_sar > 0, W.labour_cost_sar / W.handling_revenue_sar, np.nan)
    W['basis'] = 'units = Actual; money = Simulated (base, head-count basis: paid days without a scan by the person\'s main site)'
    W = W.rename(columns={s: f'units_{s}' for s in SERVICES})

    # ---------------- labour cost per case / per pallet (Binzagr costs this way at facility level) - wh x month, base
    vw = vol.groupby(['wh_id', 'month', 'family']).posted_labour_cost_sar.sum().unstack(fill_value=0)
    fams = [f for f in FAMILIES if f in vw.columns]
    d = vw[fams].sum(axis=1)
    oh = vw[OVERHEAD] if OVERHEAD in vw.columns else 0
    nsx = noscan_whm.copy(); nsx.index.names = ['wh_id', 'month']
    ns = nsx.reindex(vw.index).fillna(0)
    shr = vw[fams].div(d.where(d > 0, np.nan), axis=0).fillna(0)
    fam_full = vw[fams].add(shr.mul(oh + ns, axis=0), fill_value=0)              # head-count basis
    fam_full_w = vw[fams].add(shr.mul(oh, axis=0), fill_value=0)                 # worked time (days with scans)
    cpu = pd.DataFrame(index=vw.index)
    cpu['labour_cost_sar'] = vw.sum(axis=1) + ns
    cpu['labour_cost_worked_time_sar'] = vw.sum(axis=1)
    cpu['inbound_cost_incl_overhead_sar'] = fam_full.get('Inbound handling', 0)
    cpu['outbound_cost_incl_overhead_sar'] = fam_full.get('Outbound handling', 0)
    cpu['inbound_worked_time_cost_sar'] = fam_full_w.get('Inbound handling', 0)
    cpu['outbound_worked_time_cost_sar'] = fam_full_w.get('Outbound handling', 0)

    def agg(codes, how):
        x = tx[tx.tran_type.isin(codes)].groupby(['wh_id', 'month'])
        return (x.qty_cases.sum() if how == 'cases' else x.size()).reindex(cpu.index).fillna(0)

    def ratio(a, b):
        return np.where(b > 0, a / b.where(b > 0, 1), np.nan)
    cpu['cases_shipped'] = agg([400, 410], 'cases')          # every case picked: loose (Case Pick) + in full pallets (Pallet Pick)
    cpu['loose_cases_picked'] = agg([400], 'cases')
    cpu['case_pick_lines'] = agg([400], 'lines')
    cpu['pallets_picked'] = agg([410], 'lines')
    cpu['pallets_received'] = agg([153], 'lines')
    other_in = agg([183], 'lines') + agg([391], 'lines')
    sh_case = ratio(cpu.case_pick_lines, cpu.case_pick_lines + cpu.pallets_picked)   # same split as the rate card: billed lines
    sh_rcpt = ratio(cpu.pallets_received, cpu.pallets_received + other_in)
    cpu['case_pick_cost_sar'] = cpu.outbound_cost_incl_overhead_sar * np.nan_to_num(sh_case)
    cpu['pallet_pick_cost_sar'] = cpu.outbound_cost_incl_overhead_sar * np.nan_to_num(1 - sh_case)
    cpu['receipt_cost_sar'] = cpu.inbound_cost_incl_overhead_sar * np.nan_to_num(sh_rcpt)
    cpu['outbound_labour_per_case_shipped_sar'] = ratio(cpu.outbound_cost_incl_overhead_sar, cpu.cases_shipped)
    cpu['case_pick_labour_per_loose_case_sar'] = ratio(cpu.case_pick_cost_sar, cpu.loose_cases_picked)
    cpu['pallet_pick_labour_per_pallet_sar'] = ratio(cpu.pallet_pick_cost_sar, cpu.pallets_picked)
    cpu['receipt_labour_per_pallet_received_sar'] = ratio(cpu.receipt_cost_sar, cpu.pallets_received)
    cpu['all_labour_per_case_shipped_sar'] = ratio(cpu.labour_cost_sar, cpu.cases_shipped)
    cpu['case_pick_worked_time_labour_per_loose_case_sar'] = ratio(cpu.outbound_worked_time_cost_sar * np.nan_to_num(sh_case), cpu.loose_cases_picked)
    cpu['pallet_pick_worked_time_labour_per_pallet_sar'] = ratio(cpu.outbound_worked_time_cost_sar * np.nan_to_num(1 - sh_case), cpu.pallets_picked)
    cpu['receipt_worked_time_labour_per_pallet_received_sar'] = ratio(cpu.inbound_worked_time_cost_sar * np.nan_to_num(sh_rcpt), cpu.pallets_received)
    cpu = cpu.reset_index()
    cpu['basis'] = ('Simulated (labour cost posted by task time, base; each site-month\'s overhead and paid days without a scan spread over its '
                    'own families; worked-time columns = days with scans only) / Actual (cases, pallets from Binzagr conversion factors)')
    T8 = cpu[['labour_cost_sar', 'cases_shipped', 'loose_cases_picked', 'case_pick_lines', 'pallets_picked', 'pallets_received']].sum()
    rcs = rcw.set_index('service')
    # Jeddah totals use the rate-card allocation (overhead spread over all families together), so they equal the rate card's labour per unit
    st['cost_per_unit_8m'] = {
        'outbound_labour_per_case_shipped_sar': float(fam_total.get('Outbound handling', 0) / T8.cases_shipped),
        'case_pick_labour_per_loose_case_sar': float(rcs.loc['Case pick', 'labour_cost_per_unit_sar']),
        'pallet_pick_labour_per_pallet_sar': float(rcs.loc['Pallet pick', 'labour_cost_per_unit_sar']),
        'receipt_labour_per_pallet_received_sar': float(rcs.loc['Inbound receipt', 'labour_cost_per_unit_sar']),
        'all_labour_per_case_shipped_sar': float(T8.labour_cost_sar / T8.cases_shipped),
        'outbound_worked_time_labour_per_case_shipped_sar': float(fam_total_w.get('Outbound handling', 0) / T8.cases_shipped),
        'case_pick_worked_time_labour_per_loose_case_sar': float(rcs.loc['Case pick', 'worked_time_labour_per_unit_sar']),
        'pallet_pick_worked_time_labour_per_pallet_sar': float(rcs.loc['Pallet pick', 'worked_time_labour_per_unit_sar']),
        'receipt_worked_time_labour_per_pallet_received_sar': float(rcs.loc['Inbound receipt', 'worked_time_labour_per_unit_sar']),
        'cases_shipped': float(T8.cases_shipped), 'loose_cases_picked': float(T8.loose_cases_picked),
        'loose_share_of_cases_shipped': float(T8.loose_cases_picked / T8.cases_shipped),
        'case_pick_lines': float(T8.case_pick_lines), 'pallets_picked': float(T8.pallets_picked), 'pallets_received': float(T8.pallets_received)}

    # ---------------- test of the split by billed lines: Picking (Pick) scans linked to the Case Pick / Pallet Pick row stamped when they end
    pick_cols = ['tran_type', 'employee_id_norm', 'item_number', 'location_id', 'start_ts', 'end_ts', 'dur_s']
    linked = []
    for p in sorted(glob.glob(os.path.join(CLEAN, 'jeddah_wms_transactions', '*.parquet'))):
        z = pd.read_parquet(p, columns=pick_cols, filters=[('tran_type', 'in', [301, 400, 410])])
        pk = z[(z.tran_type == 301) & z.end_ts.notna()].sort_values('end_ts')
        bl = z[z.tran_type.isin([400, 410])][['employee_id_norm', 'item_number', 'location_id', 'start_ts', 'tran_type']]
        bl = bl.rename(columns={'start_ts': 'bill_ts', 'tran_type': 'bill_code'}).dropna(subset=['bill_ts']).sort_values('bill_ts')
        m = pd.merge_asof(pk, bl, left_on='end_ts', right_on='bill_ts', by=['employee_id_norm', 'item_number', 'location_id'],
                          direction='nearest', tolerance=pd.Timedelta('3s'))
        linked.append(m[['bill_code', 'dur_s']])
    lk = pd.concat(linked, ignore_index=True)
    lk['work_s'] = lk.dur_s.clip(lower=0, upper=3600)
    t = lk.groupby('bill_code').work_s.sum()
    st['outbound_split_check'] = {
        'picking_scans': int(len(lk)), 'linked_share': float(lk.bill_code.notna().mean()),
        'case_share_by_pick_time': float(t.get(400, 0) / (t.get(400, 0) + t.get(410, 0))),
        'case_share_by_billed_lines': float(T8.case_pick_lines / (T8.case_pick_lines + T8.pallets_picked)),
        'median_pick_seconds': {int(k): float(v) for k, v in lk.groupby('bill_code').dur_s.median().items()}}

    # ---------------- deterministic scenarios
    npz = np.load(os.path.join(AUDIT, 'labour_variants.npz'), allow_pickle=True)
    paid_shifts, ot, legal = npz['paid_shifts'], npz['overtime_h'], npz['legal_hours']
    ns_days, ns_tier, ns_conf = npz['pm_noscan_days'], npz['pm_tier'], npz['pm_confirmed'].astype(bool)
    BASES = {'scan_days': np.zeros(len(ns_days), bool), 'headcount_base': ns_conf, 'all_ids': np.ones(len(ns_days), bool)}
    FTE_BASIS = {'scan_days': L['headcount']['avg_per_month']['fte_scan_days'], 'headcount_base': L['headcount']['avg_per_month']['fte_headcount_base'],
                 'all_ids': L['headcount']['avg_per_month']['fte_all_ids']}
    observed, elig, role, tier = npz['observed_h'], npz['eligible'], npz['role'], npz['tier']
    roles = np.unique(role[elig])

    def recover(earned, q, capture):
        sp = earned / np.maximum(observed, 0.25)
        rec = np.zeros_like(earned)
        for r in roles:
            m = elig & (role == r)
            if m.sum() < 30:
                continue
            t = np.quantile(sp[m], q)
            need = earned[m] / t
            rec[m] = np.maximum(0, observed[m] - need) * capture
        return rec

    def evaluate(pv, cap, stdname, q, capture, basis='headcount_base'):
        c = loaded_cost(pv)
        shift = np.select([tier == t for t in TIERS], [c[t]['cost_per_worked_shift'] for t in TIERS])
        otc = np.select([tier == t for t in TIERS], [c[t]['overtime_cost_per_hour'] for t in TIERS])
        nsc = np.select([ns_tier == t for t in TIERS], [c[t]['cost_per_worked_shift'] for t in TIERS])
        earned = npz[f'earned_{cap}_{stdname}']
        rec_h = recover(earned, q, capture)
        lab = float((paid_shifts * shift + ot * otc).sum() + (ns_days * nsc)[BASES[basis]].sum())
        rec_sar = float((rec_h * shift / legal).sum())
        rev = lab / pv['labour_share_of_warehouse_opex'] / (1 - pv['ebit_margin_3pl'])
        ebit = rev * pv['ebit_margin_3pl']
        f = 12 / MONTHS_OBSERVED
        return {'labour_cost_annual_sar': lab * f, 'recoverable_labour_annual_sar': rec_sar * f, 'recoverable_pct_of_labour': rec_sar / lab,
                'recoverable_hours_annual': float(rec_h.sum()) * f, 'fte_recoverable': float((rec_h / legal).sum()) / MONTHS_OBSERVED / pv['working_days_per_month'],
                'handling_revenue_annual_sar': rev * f, 'ebit_annual_sar': ebit * f, 'ebit_optimised_annual_sar': (ebit + rec_sar) * f,
                'ebit_margin': ebit / rev, 'ebit_margin_optimised': (ebit + rec_sar) / rev,
                'labour_pct_of_revenue': lab / rev, 'labour_pct_of_revenue_optimised': (lab - rec_sar) / rev, 'fte_paid': FTE_BASIS[basis]}
    scen_def = {'Conservative': ('low', 30, 'lenient_p60', 0.60, 0.5, 'scan_days'), 'Base': ('base', 30, 'median_p50', 0.75, 0.7, 'headcount_base'),
                'Stretch': ('high', 30, 'tight_p40', 0.90, 0.9, 'all_ids')}
    sc_rows = []
    for nm, (cs, cap, sn, q, cp_, bs) in scen_def.items():
        r = evaluate(pdict(par, cs), cap, sn, q, cp_, bs)
        sc_rows.append({'scenario': nm, 'cost_parameters': cs, 'paid_time_basis': bs, 'idle_cap_min': cap, 'standard': sn, 'target_percentile': q,
                        'capture_rate': cp_, **r})
    SC = pd.DataFrame(sc_rows)
    SC['recoverable_labour_annual_usd'] = SC.recoverable_labour_annual_sar / FX_USD_SAR
    SC['basis'] = 'Simulated'

    # ---------------- Monte Carlo
    N = 5000
    rng = np.random.default_rng(20261021)
    pi = par.set_index('parameter')
    cont = [p for p in par.parameter if pi.loc[p, 'group'] in ('labour_expat', 'labour_saudi', 'labour_mix', 'labour_roles', 'cost')
            and not p.startswith('us_') and not p.startswith('ksa_')]
    draws = {p: tri(rng, *pi.loc[p, ['low', 'base', 'high']].astype(float).tolist(), N) for p in cont}
    q_draw = tri(rng, 0.60, 0.75, 0.90, N); cap_draw = tri(rng, 0.5, 0.7, 0.9, N)
    std_choice = rng.choice(['lenient_p60', 'median_p50', 'tight_p40'], size=N, p=[0.25, 0.5, 0.25])
    cap_choice = rng.choice([15, 30, 45], size=N, p=[0.25, 0.5, 0.25])
    basis_choice = rng.choice(['scan_days', 'headcount_base', 'all_ids'], size=N, p=[0.25, 0.5, 0.25])   # drawn last: earlier draws unchanged
    base_pv = pdict(par, 'base')
    mc = []
    for i in range(N):
        pvi = dict(base_pv); pvi.update({p: float(draws[p][i]) for p in cont})
        r = evaluate(pvi, int(cap_choice[i]), std_choice[i], float(q_draw[i]), float(cap_draw[i]), basis_choice[i])
        mc.append({'draw': i, 'target_percentile': q_draw[i], 'capture_rate': cap_draw[i], 'standard': std_choice[i], 'idle_cap_min': cap_choice[i],
                   'paid_time_basis': basis_choice[i],
                   **{p: draws[p][i] for p in cont}, **r})
    MC = pd.DataFrame(mc)
    MC['basis'] = 'Simulated'
    keys = ['labour_cost_annual_sar', 'recoverable_labour_annual_sar', 'recoverable_pct_of_labour', 'recoverable_hours_annual', 'fte_recoverable',
            'handling_revenue_annual_sar', 'ebit_annual_sar', 'ebit_optimised_annual_sar', 'ebit_margin', 'ebit_margin_optimised',
            'labour_pct_of_revenue', 'labour_pct_of_revenue_optimised']
    MS = pd.DataFrame([{'metric': k, 'p10': MC[k].quantile(.1), 'p50': MC[k].quantile(.5), 'p90': MC[k].quantile(.9), 'mean': MC[k].mean()} for k in keys])
    MS['basis'] = 'Simulated (Monte Carlo, 5,000 draws)'
    varying = [p for p in cont + ['target_percentile', 'capture_rate', 'idle_cap_min'] if MC[p].nunique() > 1]
    sens = {p: float(MC[p].rank().corr(MC.recoverable_labour_annual_sar.rank())) for p in varying}
    st['mc_median_recoverable_by_standard'] = {k: float(v) for k, v in MC.groupby('standard').recoverable_labour_annual_sar.median().items()}
    st['mc_sensitivity_spearman'] = dict(sorted(sens.items(), key=lambda kv: -abs(kv[1]))[:10])
    st['mc_labour_cost_median_by_basis'] = {k: float(v) for k, v in MC.groupby('paid_time_basis').labour_cost_annual_sar.median().items()}
    st['mc_summary'] = MS.set_index('metric')[['p10', 'p50', 'p90']].to_dict('index')
    st['scenarios'] = SC.set_index('scenario').drop(columns=['basis']).to_dict('index')

    files = []
    files += write_table(rcw, 'sim_rate_card')
    files += write_table(P, 'sim_jeddah_pnl_monthly')
    files += write_table(W, 'sim_jeddah_pnl_by_wh_month')
    files += write_table(cpu, 'sim_cost_per_unit_wh_month')
    files += write_table(SC, 'sim_scenarios')
    files += write_table(MC, 'sim_monte_carlo_draws')
    files += write_table(MS, 'sim_monte_carlo_summary')
    st['files'] = files
    save_stats('pnl', st)
    print({k: v for k, v in st.items() if k not in ('files',)})


if __name__ == '__main__':
    main()