#!/usr/bin/env python3
"""
Utility engine for the Assessment Value and Business Case Builder skill.

This is a lean, two-stream model built from textbook personnel-selection utility
analysis. It estimates the economic value of moving from a current selection or
promotion process to a stronger, assessment-informed one. Value flows from the
gain in predictive validity, which is the one parameter with strong meta-analytic
grounding (Sackett et al. 2022).

Two value streams, computed per segment (a segment is one decision population,
for example external hiring or internal promotion):

  1. Performance value (primary). Brogden-Cronbach-Gleser utility in incremental
     form:  N * delta_r * SD_y * Z_x * T
       N        = decisions per year
       delta_r  = R_proposed - R_current  (the GAIN in validity, floored at 0)
       SD_y     = 0.40 * salary           (Schmidt-Hunter; conservative)
       Z_x      = phi(c)/SR               (mean standardised score of those selected)
       T        = expected tenure (years)
     This is the value one year's decisions generate over the tenure of those
     selected. Figures are UNDISCOUNTED (no NPV), by design.

  2. Retention value (secondary; only when a turnover rate is supplied). Pure
     cost avoidance of voluntary departures among the assessment-selected people:
       expected_departures = N * (1 - (1 - VTR)^T)
       retention_value     = expected_departures * reduction_fraction * cost_per_departure
     This is separate from and additive to the performance stream (which holds
     tenure constant and counts performance, not replacement cost), so there is
     no double counting.

Scenarios (Conservative / Expected / Strong) express the genuine uncertainty in
each stream: the proposed battery's validity (credibility band) for performance,
and the modeled reduction in voluntary turnover for retention. SD_y is held flat
at the conservative 0.40 of salary across all scenarios, and the selection ratio
is a fixed input, not a scenario lever.

All financial figures are USD. Currency conversion (if inputs were non-USD)
happens before the engine; the engine assumes USD.

Usage:
    python3 roi_engine.py inputs.json results.json

See references/roi_model.md, references/validity_coefficients.md, and
references/utility_analysis_math.md for the full specification.
"""

import json
import math
import sys

import numpy as np


# ---- Validity coefficients (Sackett et al. 2022, Simplified worksheet) ----

VALIDITY = {
    "Structured interviews": 0.42,
    "Job knowledge tests": 0.40,
    "Empirically keyed biodata": 0.38,
    "Work sample tests": 0.33,
    "Cognitive ability tests": 0.31,
    "Integrity tests": 0.31,
    "Personality-based EI": 0.30,
    "Assessment centers": 0.29,
    "SJT (knowledge)": 0.26,
    "SJT (behavioral tendency)": 0.26,
    "Conscientiousness (contextualized)": 0.25,
    "Interests": 0.24,
    "Emotional stability (contextualized)": 0.23,
    "Ability-based EI": 0.22,
    "Rationally keyed biodata": 0.22,
    "Extraversion (contextualized)": 0.21,
    "Conscientiousness (overall)": 0.19,
    "Unstructured interviews": 0.19,
    "Agreeableness (contextualized)": 0.19,
    "Openness (contextualized)": 0.12,
    "Extraversion (overall)": 0.10,
    "Agreeableness (overall)": 0.10,
    "Emotional stability (overall)": 0.09,
    "Job experience (years)": 0.07,
    "Openness (overall)": 0.05,
}

SHARED_INTERCORRELATION = 0.2
VALIDITY_CAP = 0.60

# ---- Method aliases (lay terms -> canonical Sackett method name) ----
# Maps common, plain-language ways people name an assessment type to the canonical
# coefficient. Resolution is case-insensitive. This is a safety net and convenience layer;
# anything that does not resolve raises a clear error rather than being silently dropped,
# so a method is never quietly miscounted. Two of these are deliberate judgment calls,
# flagged in references/validity_coefficients.md and easy to change:
#   - a generic "personality assessment" anchors to contextualized conscientiousness (0.25),
#     the most robust single personality predictor, as a conservative representative value.
#   - a generic "values" or "motivation" assessment anchors to Interests (0.24), the closest
#     construct in the Sackett table (vocational fit), treated as a proxy.
METHOD_ALIASES = {
    "personality": "Conscientiousness (contextualized)",
    "personality assessment": "Conscientiousness (contextualized)",
    "personality test": "Conscientiousness (contextualized)",
    "personality questionnaire": "Conscientiousness (contextualized)",
    "personality inventory": "Conscientiousness (contextualized)",
    "core drivers": "Conscientiousness (contextualized)",
    "conscientiousness": "Conscientiousness (contextualized)",
    "cognitive": "Cognitive ability tests",
    "cognitive ability": "Cognitive ability tests",
    "cognitive ability test": "Cognitive ability tests",
    "cognitive test": "Cognitive ability tests",
    "cognitive assessment": "Cognitive ability tests",
    "reasoning": "Cognitive ability tests",
    "reasoning test": "Cognitive ability tests",
    "core reasoning": "Cognitive ability tests",
    "aptitude": "Cognitive ability tests",
    "aptitude test": "Cognitive ability tests",
    "gma": "Cognitive ability tests",
    "general mental ability": "Cognitive ability tests",
    "structured interview": "Structured interviews",
    "structured interviews": "Structured interviews",
    "unstructured interview": "Unstructured interviews",
    "unstructured interviews": "Unstructured interviews",
    "interview": "Unstructured interviews",
    "work sample": "Work sample tests",
    "work samples": "Work sample tests",
    "work sample test": "Work sample tests",
    "work sample exercise": "Work sample tests",
    "job knowledge": "Job knowledge tests",
    "job knowledge test": "Job knowledge tests",
    "knowledge test": "Job knowledge tests",
    "sjt": "SJT (behavioral tendency)",
    "situational judgment test": "SJT (behavioral tendency)",
    "situational judgement test": "SJT (behavioral tendency)",
    "situational judgment": "SJT (behavioral tendency)",
    "situational judgement": "SJT (behavioral tendency)",
    "integrity": "Integrity tests",
    "integrity test": "Integrity tests",
    "assessment centre": "Assessment centers",
    "assessment center": "Assessment centers",
    "biodata": "Empirically keyed biodata",
    "interests": "Interests",
    "interest inventory": "Interests",
    "values": "Interests",
    "core values": "Interests",
    "values assessment": "Interests",
    "motivation": "Interests",
    "motivation assessment": "Interests",
    "emotional intelligence": "Ability-based EI",
    "ei": "Ability-based EI",
    "eq": "Ability-based EI",
    "job experience": "Job experience (years)",
    "experience": "Job experience (years)",
    "years of experience": "Job experience (years)",
    "resume": "Job experience (years)",
    "resume screen": "Job experience (years)",
    "cv screen": "Job experience (years)",
    "resume review": "Job experience (years)",
}


def resolve_methods(methods):
    """Resolve a list of method names to canonical Sackett keys.

    Returns (resolved, unknown). A name resolves if it is already a canonical key (exact or
    case-insensitive) or matches an alias. Unknown names are returned so the caller can raise
    a clear error rather than silently dropping them.
    """
    canon_lower = {k.lower(): k for k in VALIDITY}
    resolved, unknown = [], []
    for m in (methods or []):
        key = m.strip()
        low = key.lower()
        # Tolerate a leading article (a / an / the) on plain-language names.
        for art in ("a ", "an ", "the "):
            if low.startswith(art):
                low = low[len(art):].strip()
                break
        if key in VALIDITY:
            resolved.append(key)
        elif key.lower() in canon_lower:
            resolved.append(canon_lower[key.lower()])
        elif low in canon_lower:
            resolved.append(canon_lower[low])
        elif key.lower() in METHOD_ALIASES:
            resolved.append(METHOD_ALIASES[key.lower()])
        elif low in METHOD_ALIASES:
            resolved.append(METHOD_ALIASES[low])
        else:
            unknown.append(m)
    return resolved, unknown

# Default baseline if a segment does not name its current methods: an unstructured
# interview (Sackett 2022 r = 0.19). Assuming a stronger baseline would only lower
# the estimated value, so this default is the value-maximising end; the methodology
# states it and invites the user to specify a more rigorous current process.
DEFAULT_CURRENT_METHODS = ["Unstructured interviews"]

# ---- Scenario conventions ----

# Validity credibility band applied to the PROPOSED battery's combined R. This is
# the genuine source of uncertainty in the performance stream. Current-process
# validity is held at its point estimate (the user is telling us what they do now).
VALIDITY_LOW_MULT = 0.85
VALIDITY_HIGH_MULT = 1.10

# Modeled relative reduction in voluntary turnover among assessment-selected people.
# This is the softest parameter in the model and is what the retention scenarios flex.
# Deliberately conservative. Grounded in the general meta-analytic finding that better
# person-job and person-organisation fit reduces voluntary turnover (see references).
TURNOVER_REDUCTION = {"low": 0.05, "mid": 0.10, "high": 0.15}

# Dollar value of one SD of performance, as a share of salary. Held FLAT across all
# scenarios at the conservative Schmidt-Hunter lower bound.
SDY_MULT = 0.40

# Defaults for optional inputs.
DEFAULT_SR = 0.20
DEFAULT_TENURE = 3.0
DEFAULT_COST_PER_TURNOVER_MULT = 0.50  # share of salary if no cost figure supplied

SCENARIOS = ("low", "mid", "high")
SCENARIO_LABELS = {"low": "Conservative", "mid": "Expected", "high": "Strong"}
VALIDITY_MULT = {"low": VALIDITY_LOW_MULT, "mid": 1.0, "high": VALIDITY_HIGH_MULT}


# ---- Normal-distribution helpers (for an exact Z_x from any selection ratio) ----

def _norm_pdf(z):
    return math.exp(-0.5 * z * z) / math.sqrt(2.0 * math.pi)


def _norm_ppf(p):
    """Inverse normal CDF via Acklam's rational approximation. Good to ~1e-9."""
    if p <= 0.0:
        return -math.inf
    if p >= 1.0:
        return math.inf
    a = [-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
         1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00]
    b = [-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
         6.680131188771972e+01, -1.328068155288572e+01]
    c = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
         -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00]
    d = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
         3.754408661907416e+00]
    plow, phigh = 0.02425, 1.0 - 0.02425
    if p < plow:
        q = math.sqrt(-2.0 * math.log(p))
        return (((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) / \
               ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0)
    if p > phigh:
        q = math.sqrt(-2.0 * math.log(1.0 - p))
        return -(((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) / \
                 ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0)
    q = p - 0.5
    r = q * q
    return (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5]) * q / \
           (((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0)


def zx_from_sr(sr):
    """Mean standardised predictor score of those selected at selection ratio sr.

    Z_x = phi(c) / sr, where c is the cutoff such that P(Z > c) = sr.
    """
    if not sr or sr <= 0.0 or sr >= 1.0:
        sr = DEFAULT_SR
    c = _norm_ppf(1.0 - sr)
    return _norm_pdf(c) / sr


# ---- Validity combination (multiple-R, uniform 0.2 intercorrelation, cap 0.60) ----

def combine_validity(methods):
    """Combine assessment method names into a single multiple-R. Returns (R, [r...])."""
    rs = [VALIDITY[m] for m in (methods or []) if m in VALIDITY]
    if not rs:
        return None, []
    if len(rs) == 1:
        return round(min(rs[0], VALIDITY_CAP), 4), rs
    r = np.array(rs, dtype=float)
    n = len(rs)
    Rxx = np.full((n, n), SHARED_INTERCORRELATION, dtype=float)
    np.fill_diagonal(Rxx, 1.0)
    try:
        Rxx_inv = np.linalg.inv(Rxx)
        R_squared = float(r.T @ Rxx_inv @ r)
        R = R_squared ** 0.5 if R_squared > 0 else max(rs)
    except np.linalg.LinAlgError:
        R = max(rs)
    return round(min(R, VALIDITY_CAP), 4), rs


def _cap(x):
    return round(min(x, VALIDITY_CAP), 4)


# ---- Lay translations (relative to the segment's OWN current process) ----

def _lifted_top_quartile_rate(r, zx):
    """Share of selected who are genuine top-quartile performers (base rate 0.25)."""
    base = 0.25
    return max(0.0, min(base + r * zx * base, 1.0))


# ---- Per-segment computation ----

def compute_segment(seg):
    name = seg.get("name") or "This population"
    decision_type = (seg.get("decision_type") or "hiring").lower()
    n = float(seg.get("n") or 0)
    salary = float(seg.get("salary") or 0)
    sr = seg.get("selection_ratio")
    sr = float(sr) if sr else DEFAULT_SR
    tenure = seg.get("tenure")
    tenure = float(tenure) if tenure else DEFAULT_TENURE

    current_methods_in = seg.get("current_methods") or DEFAULT_CURRENT_METHODS
    current_is_default = not seg.get("current_methods")
    proposed_methods_in = seg.get("proposed_methods") or []

    current_methods, cur_unknown = resolve_methods(current_methods_in)
    proposed_methods, prop_unknown = resolve_methods(proposed_methods_in)
    if cur_unknown:
        raise ValueError(
            f"Segment '{name}': unrecognised current method(s): {cur_unknown}. "
            f"Map them to a name or alias in references/validity_coefficients.md.")
    if prop_unknown:
        raise ValueError(
            f"Segment '{name}': unrecognised proposed method(s): {prop_unknown}. "
            f"Map them to a name or alias in references/validity_coefficients.md.")

    R_current, _ = combine_validity(current_methods)
    if R_current is None:
        R_current = VALIDITY["Unstructured interviews"]
        current_is_default = True
    R_proposed, _ = combine_validity(proposed_methods)
    if R_proposed is None:
        raise ValueError(f"Segment '{name}' has no recognised proposed_methods.")

    zx = zx_from_sr(sr)
    sd_y = SDY_MULT * salary

    # Retention inputs (optional)
    vtr = seg.get("turnover_rate")
    retention_included = vtr is not None
    vtr = float(vtr) if retention_included else None
    if retention_included:
        cpt = seg.get("cost_per_turnover")
        cost_per_departure = float(cpt) if cpt else DEFAULT_COST_PER_TURNOVER_MULT * salary
        cost_is_default = not cpt
        expected_departures = n * (1.0 - (1.0 - vtr) ** tenure)
    else:
        cost_per_departure = None
        cost_is_default = None
        expected_departures = None

    perf = {}
    ret = {}
    total = {}
    for s in SCENARIOS:
        R_prop_s = _cap(R_proposed * VALIDITY_MULT[s])
        delta_r = max(0.0, R_prop_s - R_current)
        perf_s = n * delta_r * sd_y * zx * tenure
        perf[s] = round(perf_s, 2)
        if retention_included:
            ret_s = expected_departures * TURNOVER_REDUCTION[s] * cost_per_departure
            ret[s] = round(ret_s, 2)
        else:
            ret[s] = None
        total[s] = round(perf[s] + (ret[s] or 0.0), 2)

    # Lay translations at the point (Expected) estimates, vs the current process.
    lifted_prop = _lifted_top_quartile_rate(R_proposed, zx)
    lifted_cur = _lifted_top_quartile_rate(R_current, zx)
    top_mult = (lifted_prop / lifted_cur) if lifted_cur > 0 else None
    dq_lift = ((R_proposed - R_current) / R_current) if R_current > 0 else None

    return {
        "name": name,
        "decision_type": decision_type,
        "n": n,
        "salary": salary,
        "selection_ratio": sr,
        "tenure": tenure,
        "current_methods": current_methods,
        "current_is_default": current_is_default,
        "proposed_methods": proposed_methods,
        "R_current": round(R_current, 4),
        "R_proposed": round(R_proposed, 4),
        "delta_r": round(max(0.0, R_proposed - R_current), 4),
        "Z_x": round(zx, 4),
        "SD_y": round(sd_y, 2),
        "retention_included": retention_included,
        "turnover_rate": vtr,
        "cost_per_departure": round(cost_per_departure, 2) if cost_per_departure else None,
        "cost_is_default": cost_is_default,
        "expected_departures_over_tenure": round(expected_departures, 3) if expected_departures is not None else None,
        "performance": perf,
        "retention": ret,
        "total": total,
        "lay": {
            "top_performer_multiplier": round(top_mult, 2) if top_mult else None,
            "decision_quality_lift": round(dq_lift, 4) if dq_lift is not None else None,
            "hit_rate_proposed": round(lifted_prop, 4),
            "hit_rate_current": round(lifted_cur, 4),
        },
    }


def compute(inputs):
    budget = inputs.get("budget")
    segments_in = inputs.get("segments") or []
    if not segments_in:
        raise ValueError("inputs must include at least one segment.")

    segments = [compute_segment(s) for s in segments_in]
    any_retention = any(s["retention_included"] for s in segments)

    # Aggregate per scenario
    scenarios = {}
    for s in SCENARIOS:
        perf_total = round(sum(seg["performance"][s] for seg in segments), 2)
        ret_total = round(sum((seg["retention"][s] or 0.0) for seg in segments), 2) if any_retention else None
        value_total = round(perf_total + (ret_total or 0.0), 2)
        bcr = round(value_total / budget, 4) if budget else None
        net = round(value_total - budget, 2) if budget else None
        scenarios[s] = {
            "label": SCENARIO_LABELS[s],
            "performance_value": perf_total,
            "retention_value": ret_total,
            "value_total": value_total,
            "benefit_cost_ratio": bcr,
            "net_value": net,
        }

    # Decision-type label for the narrative
    dtypes = {seg["decision_type"] for seg in segments}
    if dtypes == {"hiring"}:
        decision_label = "hiring decisions"
    elif dtypes == {"promotion"}:
        decision_label = "promotion and succession decisions"
    else:
        decision_label = "hiring and promotion decisions"

    return {
        "currency": "$",
        "currency_note": inputs.get("currency_note"),
        "client_name": inputs.get("client_name"),
        "org_size": inputs.get("org_size"),
        "budget": budget,
        "decision_label": decision_label,
        "any_retention": any_retention,
        "segments": segments,
        "scenarios": scenarios,
        "model": {
            "sdy_mult": SDY_MULT,
            "default_sr": DEFAULT_SR,
            "default_tenure": DEFAULT_TENURE,
            "default_cost_per_turnover_mult": DEFAULT_COST_PER_TURNOVER_MULT,
            "validity_band": {"low": VALIDITY_LOW_MULT, "high": VALIDITY_HIGH_MULT, "cap": VALIDITY_CAP},
            "turnover_reduction": TURNOVER_REDUCTION,
            "discounted": False,
        },
    }


def main():
    if len(sys.argv) < 3:
        print("Usage: roi_engine.py inputs.json results.json", file=sys.stderr)
        sys.exit(1)
    with open(sys.argv[1]) as f:
        inputs = json.load(f)
    results = compute(inputs)
    with open(sys.argv[2], "w") as f:
        json.dump(results, f, indent=2)
    print(f"Wrote {sys.argv[2]}")
    cur = results["currency"]
    b = results["budget"]
    for s in SCENARIOS:
        sc = results["scenarios"][s]
        line = (f"  {sc['label']:13s} perf {cur}{sc['performance_value']:,.0f}"
                f"  total {cur}{sc['value_total']:,.0f}")
        if sc["retention_value"] is not None:
            line += f"  (retention {cur}{sc['retention_value']:,.0f})"
        if sc["benefit_cost_ratio"] is not None:
            line += f"  BCR {sc['benefit_cost_ratio']:.2f}x"
        print(line)
    for seg in results["segments"]:
        print(f"    segment '{seg['name']}': R_current {seg['R_current']}, "
              f"R_proposed {seg['R_proposed']}, delta_r {seg['delta_r']}, Z_x {seg['Z_x']}")


if __name__ == "__main__":
    main()
