Tooluniverse Computational Biophysics

"""
Burn resuscitation fluid calculator.

Supports three resuscitation formulas:

1. Parkland (Baxter) formula — adults and older children:
    Total 24h = 4 * weight_kg * %TBSA   [mL of Lactated Ringer's]
    First 8h  : half the total (from TIME OF BURN, not admission)
    Next 16h  : other half

2. Modified Brooke formula — alternative for adults/pediatrics:
    Total 24h = 2 * weight_kg * %TBSA   [mL of LR]
    First 8h  : half the total
    Next 16h  : other half
    For pediatric patients: add maintenance fluids (4-2-1 rule) on top.

3. Galveston formula — pediatric (BSA-based):
    Total 24h = 5000 * BSA_m2 * (%TBSA / 100) + 2000 * BSA_m2
    (resuscitation component + maintenance component)
    First 8h  : half the total
    Next 16h  : other half

Pediatric maintenance fluids (4-2-1 rule, added to Parkland/Brooke for children):
    First 10 kg  : 4 mL/kg/h
    Next 10 kg   : 2 mL/kg/h
    Each kg > 20 : 1 mL/kg/h

Holliday-Segar maintenance fluid (daily method, shown for reference):
    <= 10 kg : 100 mL/kg/day
    10-20 kg : 1000 + 50 mL/kg/day above 10 kg
    > 20 kg  : 1500 + 20 mL/kg/day above 20 kg

Second 24h (not auto-calculated here):
    Colloid: 0.3-0.5 mL/kg/%TBSA  (start at ~18-24 h)
    Free water to maintain urine output 0.5-1 mL/kg/h

Usage:
    # Adult (80 kg, 40% TBSA) — Parkland (default for >= 30 kg):
    python burn_fluids.py --weight_kg 80 --tbsa_pct 40

    # Pediatric, modified Brooke + maintenance (25 kg, 45% TBSA):
    python burn_fluids.py --weight_kg 25 --age_years 7 --tbsa_pct 45 --formula brooke

    # Pediatric, Galveston (25 kg, 45% TBSA, BSA 0.95 m2):
    python burn_fluids.py --weight_kg 25 --age_years 7 --tbsa_pct 45 --bsa_m2 0.95

    # Adult, modified Brooke:
    python burn_fluids.py --weight_kg 80 --tbsa_pct 40 --formula brooke
"""

import argparse
import sys


def maintenance_4_2_1(weight_kg: float) -> dict:
    """
    Maintenance fluid rate by the 4-2-1 rule (hourly method).
    Returns hourly rate (mL/h) and daily equivalent (mL/day).

    First 10 kg  : 4 mL/kg/h
    Next 10 kg   : 2 mL/kg/h
    Each kg > 20 : 1 mL/kg/h
    """
    if weight_kg <= 10:
        hourly = 4.0 * weight_kg
    elif weight_kg <= 20:
        hourly = 40.0 + 2.0 * (weight_kg - 10.0)
    else:
        hourly = 60.0 + 1.0 * (weight_kg - 20.0)
    return {"hourly_mL_per_h": hourly, "daily_mL": hourly * 24.0}


def holliday_segar(weight_kg: float) -> dict:
    """
    Maintenance fluid by Holliday-Segar method (daily method).
    Returns daily (mL/day) and hourly (mL/h) rates.
    """
    if weight_kg <= 10:
        daily = 100.0 * weight_kg
    elif weight_kg <= 20:
        daily = 1000.0 + 50.0 * (weight_kg - 10.0)
    else:
        daily = 1500.0 + 20.0 * (weight_kg - 20.0)
    return {"daily_mL": daily, "hourly_mL_per_h": daily / 24.0}


def parkland(weight_kg: float, tbsa_pct: float) -> dict:
    """
    Parkland (Baxter) formula for burn resuscitation.
    Total 24h = 4 * weight_kg * %TBSA mL of LR.
    """
    total_24h = 4.0 * weight_kg * tbsa_pct
    first_8h_total = total_24h / 2.0
    next_16h_total = total_24h / 2.0
    return {
        "formula": "Parkland",
        "coefficient": 4,
        "total_24h_mL": total_24h,
        "first_8h_mL": first_8h_total,
        "next_16h_mL": next_16h_total,
        "rate_first_8h_mL_per_h": first_8h_total / 8.0,
        "rate_next_16h_mL_per_h": next_16h_total / 16.0,
    }


def modified_brooke(weight_kg: float, tbsa_pct: float) -> dict:
    """
    Modified Brooke formula for burn resuscitation.
    Total 24h = 2 * weight_kg * %TBSA mL of LR.
    """
    total_24h = 2.0 * weight_kg * tbsa_pct
    first_8h_total = total_24h / 2.0
    next_16h_total = total_24h / 2.0
    return {
        "formula": "Modified Brooke",
        "coefficient": 2,
        "total_24h_mL": total_24h,
        "first_8h_mL": first_8h_total,
        "next_16h_mL": next_16h_total,
        "rate_first_8h_mL_per_h": first_8h_total / 8.0,
        "rate_next_16h_mL_per_h": next_16h_total / 16.0,
    }


def galveston(bsa_m2: float, tbsa_pct: float) -> dict:
    """
    Galveston formula for pediatric burn resuscitation.
    Total 24h = 5000 * BSA_m2 * (%TBSA/100) + 2000 * BSA_m2 mL of LR.
    """
    resuscitation_component = 5000.0 * bsa_m2 * (tbsa_pct / 100.0)
    maintenance_component = 2000.0 * bsa_m2
    total_24h = resuscitation_component + maintenance_component
    first_8h_total = total_24h / 2.0
    next_16h_total = total_24h / 2.0
    return {
        "formula": "Galveston",
        "resuscitation_component_mL": resuscitation_component,
        "maintenance_component_mL": maintenance_component,
        "total_24h_mL": total_24h,
        "first_8h_mL": first_8h_total,
        "next_16h_mL": next_16h_total,
        "rate_first_8h_mL_per_h": first_8h_total / 8.0,
        "rate_next_16h_mL_per_h": next_16h_total / 16.0,
    }


def main():
    parser = argparse.ArgumentParser(
        description="Burn resuscitation fluid calculator (Parkland / Modified Brooke / Galveston).",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=__doc__,
    )
    parser.add_argument("--weight_kg", type=float, required=True, help="Patient weight in kg.")
    parser.add_argument(
        "--tbsa_pct",
        type=float,
        required=True,
        help="Total body surface area burned (%%TBSA), e.g. 45 for 45%%.",
    )
    parser.add_argument(
        "--age_years", type=float, default=None, help="Patient age in years (triggers pediatric mode if < 14)."
    )
    parser.add_argument(
        "--bsa_m2",
        type=float,
        default=None,
        help="Body surface area in m2 (required for Galveston formula).",
    )
    parser.add_argument(
        "--formula",
        choices=["parkland", "brooke", "galveston"],
        default=None,
        help="Resuscitation formula: parkland (4 mL/kg/%%TBSA), brooke (2 mL/kg/%%TBSA), "
        "galveston (BSA-based, pediatric). Default: galveston if pediatric + BSA given, "
        "else parkland for adults.",
    )
    args = parser.parse_args()

    weight = args.weight_kg
    tbsa = args.tbsa_pct

    # Validate inputs
    if weight <= 0:
        print("Error: weight_kg must be positive.")
        sys.exit(1)
    if not (0 < tbsa <= 100):
        print(f"Error: tbsa_pct must be between 0 and 100 (got {tbsa}).")
        sys.exit(1)
    if args.bsa_m2 is not None and args.bsa_m2 <= 0:
        print(f"Error: bsa_m2 must be positive (got {args.bsa_m2}).")
        sys.exit(1)

    # Determine if pediatric
    is_pediatric = (args.age_years is not None and args.age_years < 14) or weight < 30

    # Determine formula
    if args.formula is not None:
        formula = args.formula
    elif is_pediatric and args.bsa_m2 is not None:
        formula = "galveston"
    elif is_pediatric:
        formula = "brooke"
    else:
        formula = "parkland"

    if formula == "galveston" and args.bsa_m2 is None:
        print("Error: --bsa_m2 is required for the Galveston formula.")
        print("  Provide patient BSA in m2, or use --formula parkland/brooke.")
        sys.exit(1)

    print("=" * 62)
    print("  Burn Resuscitation Fluid Calculator")
    print("=" * 62)
    age_str = f"{args.age_years} y" if args.age_years is not None else "not specified"
    pop_str = "Pediatric" if is_pediatric else "Adult"
    print(f"  Weight       : {weight} kg")
    print(f"  Age          : {age_str}")
    print(f"  %TBSA burned : {tbsa}%")
    if args.bsa_m2 is not None:
        print(f"  BSA          : {args.bsa_m2} m2")
    print(f"  Population   : {pop_str}")

    # ---- Compute resuscitation ----
    if formula == "galveston":
        r = galveston(args.bsa_m2, tbsa)
    elif formula == "brooke":
        r = modified_brooke(weight, tbsa)
    else:
        r = parkland(weight, tbsa)

    # ---- Display formula details ----
    print()
    if formula == "galveston":
        print("  Galveston Formula (LR):")
        print(f"    Resuscitation component : {r['resuscitation_component_mL']:.1f} mL")
        print(f"      = 5000 x {args.bsa_m2} m2 x {tbsa / 100:.2f}")
        print(f"    Maintenance component   : {r['maintenance_component_mL']:.1f} mL")
        print(f"      = 2000 x {args.bsa_m2} m2")
    elif formula == "brooke":
        print("  Modified Brooke Formula (LR):")
        print(f"    Total = 2 x {weight} kg x {tbsa}% = {r['total_24h_mL']:.1f} mL")
    else:
        print("  Parkland Formula (LR):")
        print(f"    Total = 4 x {weight} kg x {tbsa}% = {r['total_24h_mL']:.1f} mL")

    # ---- Maintenance for pediatric (Parkland/Brooke only) ----
    maint_hourly = 0.0
    add_maintenance = is_pediatric and formula != "galveston"
    if add_maintenance:
        m = maintenance_4_2_1(weight)
        maint_hourly = m["hourly_mL_per_h"]

    # Combined rates
    resus_first_8h = r["rate_first_8h_mL_per_h"]
    resus_next_16h = r["rate_next_16h_mL_per_h"]
    total_first_8h_rate = resus_first_8h + maint_hourly
    total_next_16h_rate = resus_next_16h + maint_hourly

    print()
    print("  +---------------------------------------------------------+")
    print(f"  |  Total resuscitation 24h : {r['total_24h_mL']:>8.1f} mL                |")
    print(f"  |  First 8h  (half)        : {r['first_8h_mL']:>8.1f} mL                |")
    print(f"  |    resuscitation rate     : {resus_first_8h:>8.1f} mL/h              |")
    if add_maintenance:
        print(f"  |    + maintenance (4-2-1) : {maint_hourly:>8.1f} mL/h              |")
        print(f"  |    = TOTAL first 8h rate : {total_first_8h_rate:>8.1f} mL/h              |")
    print(f"  |  Next 16h  (half)        : {r['next_16h_mL']:>8.1f} mL                |")
    print(f"  |    resuscitation rate     : {resus_next_16h:>8.1f} mL/h              |")
    if add_maintenance:
        print(f"  |    + maintenance (4-2-1) : {maint_hourly:>8.1f} mL/h              |")
        print(f"  |    = TOTAL next 16h rate : {total_next_16h_rate:>8.1f} mL/h              |")
    print("  +---------------------------------------------------------+")

    if add_maintenance:
        print()
        _print_421_breakdown(weight)

    print()
    print("  IMPORTANT: The 8h clock starts from TIME OF BURN,")
    print("  not time of admission. Adjust if fluid was already given.")

    # ---- Holliday-Segar maintenance (reference) ----
    hs = holliday_segar(weight)
    print()
    print("  Holliday-Segar Maintenance (reference):")
    if weight <= 10:
        rule = f"100 x {weight} kg"
    elif weight <= 20:
        rule = f"1000 + 50 x {weight - 10:.0f} kg (above 10 kg)"
    else:
        rule = f"1500 + 20 x {weight - 20:.0f} kg (above 20 kg)"
    print(f"    Formula : {rule}")
    print(f"    Daily   : {hs['daily_mL']:.1f} mL/day")
    print(f"    Hourly  : {hs['hourly_mL_per_h']:.2f} mL/h")

    if formula == "galveston":
        print()
        print("  Note: Galveston formula already includes maintenance.")
        print("  Holliday-Segar shown above for reference only.")

    # ---- Urine output target ----
    uo_min = 0.5 * weight
    uo_max = 1.0 * weight
    print()
    print(f"  Urine output target  : {uo_min:.1f}-{uo_max:.1f} mL/h  (0.5-1 mL/kg/h)")
    print(f"  Fluid type           : Lactated Ringer's (isotonic)")
    print()
    print("  Titrate infusion rate to urine output; above values are")
    print("  initial estimates only. Reassess hourly.")
    print("=" * 62)


def _print_421_breakdown(weight_kg: float) -> None:
    """Print the 4-2-1 maintenance calculation breakdown."""
    print("  4-2-1 Maintenance Breakdown:")
    if weight_kg <= 10:
        print(f"    4 x {weight_kg:.0f} kg = {4.0 * weight_kg:.0f} mL/h")
    elif weight_kg <= 20:
        extra = weight_kg - 10.0
        print(f"    4 x 10 kg                = 40 mL/h")
        print(f"    2 x {extra:.0f} kg (next 10 kg)  = {2.0 * extra:.0f} mL/h")
        print(f"    Total                    = {40.0 + 2.0 * extra:.0f} mL/h")
    else:
        extra = weight_kg - 20.0
        print(f"    4 x 10 kg                = 40 mL/h")
        print(f"    2 x 10 kg (next 10 kg)   = 20 mL/h")
        print(f"    1 x {extra:.0f} kg (above 20 kg) = {1.0 * extra:.0f} mL/h")
        print(f"    Total                    = {60.0 + 1.0 * extra:.0f} mL/h")


if __name__ == "__main__":
    main()

"""
Environmental risk assessment calculator for contaminated site exposure.

Computes hazard quotient (HQ) for human exposure to soil contaminants
via food ingestion pathway (soil → soil solution → plant → human).

Calculation steps:
  1. Soil concentration: C_soil = total_mass / soil_mass
  2. Soil-water partitioning: C_sw = C_soil * rho_b / (Kd * rho_b + theta_w)
     where Kd = foc * Koc
  3. Plant uptake: C_plant = C_sw * TSCF * PUF
  4. Daily intake: DI = C_plant * intake_rate_g * bioavailability
  5. Average daily dose: ADD = total_DI / body_weight
  6. Hazard quotient: HQ = ADD / RfD

Key parameter notes:
  - Koc (organic carbon partition coefficient) is compound-specific.
    Common values: PFHxS ~36, PFOS ~200, PFOA ~70 L/kg.
  - TSCF (transpiration stream concentration factor): dimensionless,
    typically 0.1-10 for different compounds.
  - PUF (plant uptake factor): fraction of root-zone contaminant
    transferred to edible parts.

Usage:
    # Basic: single food pathway
    python env_risk_assessment.py \\
        --total_mass_ug 1e9 --area_m2 250000 --depth_m 0.6 \\
        --bulk_density 1500 --theta_w 0.35 --foc 0.03 --Koc 36 \\
        --food "fruit:300:0.5:0.1:5" \\
        --body_weight 80 --RfD 0.02

    # Multiple food pathways
    python env_risk_assessment.py \\
        --total_mass_ug 1e9 --area_m2 250000 --depth_m 0.6 \\
        --bulk_density 1500 --theta_w 0.35 --foc 0.03 --Koc 36 \\
        --food "fruit:300:0.5:0.1:5" --food "legume:50:0.3:0.2:5" \\
        --body_weight 80 --RfD 0.02

    Food format: "name:intake_g_per_day:bioavailability:plant_uptake_factor:TSCF"
"""

import argparse
import sys


def compute_soil_concentration(total_mass_ug: float, area_m2: float,
                                depth_m: float, bulk_density_kg_m3: float) -> float:
    """Soil concentration in ug/kg. C_soil = total_mass / (area * depth * rho_b)."""
    soil_mass_kg = area_m2 * depth_m * bulk_density_kg_m3
    return total_mass_ug / soil_mass_kg


def compute_soil_solution(C_soil_ug_kg: float, bulk_density: float,
                           theta_w: float, Kd: float) -> float:
    """
    Soil solution (pore water) concentration in ug/L.
    C_sw = C_soil * rho_b / (Kd * rho_b + theta_w)
    Note: rho_b in kg/m3, theta_w in L/L (= m3/m3), Kd in L/kg.
    Result: ug/m3 internally, converted to ug/L by dividing by 1000.
    """
    # C_soil [ug/kg] * rho_b [kg/m3] = ug/m3 in bulk soil
    # Kd*rho_b [L/kg * kg/m3 = L/m3], theta_w [L/L = L/(0.001 m3)]
    # Denominator: Kd*rho_b + theta_w  (in consistent units)
    # For the partition: C_sw(ug/L) = C_soil(ug/kg) * rho_b(kg/L) / (Kd(L/kg)*rho_b(kg/L) + theta_w(L/L))
    # rho_b in kg/L = rho_b_kg_m3 / 1000
    rho_b_kg_L = bulk_density / 1000.0
    C_sw = C_soil_ug_kg * rho_b_kg_L / (Kd * rho_b_kg_L + theta_w)
    return C_sw


def compute_daily_intake(C_sw_ug_L: float, TSCF: float, PUF: float,
                          intake_g_day: float, bioavailability: float) -> float:
    """
    Daily intake in ug/day from one food pathway.
    DI = C_sw * TSCF * PUF * intake_g * bioavailability
    """
    return C_sw_ug_L * TSCF * PUF * intake_g_day * bioavailability


def compute_hazard_quotient(total_DI_ug_day: float, body_weight_kg: float,
                             RfD_ug_kg_day: float) -> float:
    """HQ = (DI / BW) / RfD."""
    ADD = total_DI_ug_day / body_weight_kg
    return ADD / RfD_ug_kg_day


def parse_food(food_str: str) -> dict:
    """Parse food specification: 'name:intake_g:bioavailability:PUF:TSCF'."""
    parts = food_str.split(":")
    if len(parts) != 5:
        raise ValueError(
            f"Food must have 5 colon-separated fields "
            f"(name:intake_g:bioavailability:PUF:TSCF), got: {food_str}"
        )
    return {
        "name": parts[0],
        "intake_g": float(parts[1]),
        "bioavailability": float(parts[2]),
        "PUF": float(parts[3]),
        "TSCF": float(parts[4]),
    }


def main():
    parser = argparse.ArgumentParser(
        description="Environmental risk assessment: hazard quotient from soil contamination.",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=__doc__,
    )
    parser.add_argument("--total_mass_ug", type=float, required=True,
                        help="Total contaminant mass in soil (ug).")
    parser.add_argument("--area_m2", type=float, required=True,
                        help="Contaminated area (m2).")
    parser.add_argument("--depth_m", type=float, required=True,
                        help="Contamination depth (m).")
    parser.add_argument("--bulk_density", type=float, required=True,
                        help="Soil bulk density (kg/m3).")
    parser.add_argument("--theta_w", type=float, required=True,
                        help="Volumetric water content (L water / L soil).")
    parser.add_argument("--foc", type=float, required=True,
                        help="Fraction organic carbon in soil (0-1).")
    parser.add_argument("--Koc", type=float, required=True,
                        help="Organic carbon partition coefficient (L/kg).")
    parser.add_argument("--food", action="append", required=True,
                        help="Food pathway: 'name:intake_g:bioavailability:PUF:TSCF'. "
                             "Can be specified multiple times.")
    parser.add_argument("--body_weight", type=float, required=True,
                        help="Body weight (kg).")
    parser.add_argument("--RfD", type=float, required=True,
                        help="Reference dose (ug/kg body weight/day).")
    args = parser.parse_args()

    # Parse food pathways
    foods = []
    for f_str in args.food:
        try:
            foods.append(parse_food(f_str))
        except ValueError as e:
            print(f"Error: {e}", file=sys.stderr)
            sys.exit(1)

    # Step 1: Soil concentration
    C_soil = compute_soil_concentration(
        args.total_mass_ug, args.area_m2, args.depth_m, args.bulk_density
    )
    soil_mass = args.area_m2 * args.depth_m * args.bulk_density

    # Step 2: Soil-water partition
    Kd = args.foc * args.Koc
    C_sw = compute_soil_solution(C_soil, args.bulk_density, args.theta_w, Kd)

    # Step 3-4: Daily intake per food pathway
    food_results = []
    total_DI = 0.0
    for food in foods:
        DI = compute_daily_intake(C_sw, food["TSCF"], food["PUF"],
                                   food["intake_g"], food["bioavailability"])
        food_results.append({**food, "DI": DI})
        total_DI += DI

    # Step 5-6: Dose and HQ
    ADD = total_DI / args.body_weight
    HQ = compute_hazard_quotient(total_DI, args.body_weight, args.RfD)

    # --- Output ---
    print("=" * 65)
    print("  Environmental Risk Assessment — Hazard Quotient")
    print("=" * 65)

    print()
    print("  Step 1: Soil Concentration")
    print(f"    Total contaminant  : {args.total_mass_ug:.2e} ug")
    print(f"    Area               : {args.area_m2:,.0f} m2")
    print(f"    Depth              : {args.depth_m} m")
    print(f"    Bulk density       : {args.bulk_density} kg/m3")
    print(f"    Soil mass          : {soil_mass:,.0f} kg")
    print(f"    C_soil = {args.total_mass_ug:.2e} / {soil_mass:.2e}")
    print(f"           = {C_soil:.4f} ug/kg")

    print()
    print("  Step 2: Soil Solution Concentration (Kd-corrected)")
    print(f"    foc = {args.foc},  Koc = {args.Koc} L/kg")
    print(f"    Kd = foc * Koc = {Kd:.4f} L/kg")
    print(f"    theta_w = {args.theta_w} L/L")
    rho_b_L = args.bulk_density / 1000.0
    print(f"    rho_b = {rho_b_L:.3f} kg/L")
    denom = Kd * rho_b_L + args.theta_w
    print(f"    C_sw = C_soil * rho_b / (Kd*rho_b + theta_w)")
    print(f"         = {C_soil:.4f} * {rho_b_L:.3f} / ({Kd:.4f}*{rho_b_L:.3f} + {args.theta_w})")
    print(f"         = {C_soil * rho_b_L:.4f} / {denom:.4f}")
    print(f"         = {C_sw:.4f} ug/L")

    print()
    print("  Step 3-4: Daily Intake by Food Pathway")
    for fr in food_results:
        print(f"    {fr['name']}:")
        print(f"      Intake = {fr['intake_g']:.0f} g/day,  BF = {fr['bioavailability']},  "
              f"PUF = {fr['PUF']},  TSCF = {fr['TSCF']}")
        print(f"      DI = {C_sw:.4f} * {fr['TSCF']} * {fr['PUF']} * {fr['intake_g']:.0f} * {fr['bioavailability']}")
        print(f"         = {fr['DI']:.4f} ug/day")
    print(f"    ----------------------------------------")
    print(f"    Total daily intake = {total_DI:.4f} ug/day")

    print()
    print("  Step 5: Average Daily Dose")
    print(f"    ADD = {total_DI:.4f} / {args.body_weight} = {ADD:.4f} ug/kg/day")

    print()
    print("  Step 6: Hazard Quotient")
    print(f"    RfD = {args.RfD} ug/kg/day")
    print(f"    HQ  = ADD / RfD = {ADD:.4f} / {args.RfD}")
    print(f"        = {HQ:.1f}")
    print()

    if HQ > 1:
        print(f"    RESULT: HQ = {HQ:.1f} >> 1 — significant health risk.")
    else:
        print(f"    RESULT: HQ = {HQ:.2f} < 1 — within acceptable risk level.")

    print()
    print("  Verification:")
    HQ_check = (C_sw * sum(f["TSCF"] * f["PUF"] * f["intake_g"] * f["bioavailability"]
                            for f in foods)) / args.body_weight / args.RfD
    print(f"    HQ recomputed: {HQ_check:.1f}  (matches: {abs(HQ_check - HQ) < 0.1})")
    print("=" * 65)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Enzyme kinetics calculator: Km/Vmax, Hill coefficient, and Ki determination.

No external dependencies — uses only math stdlib.

Usage:
  python enzyme_kinetics.py --type km_vmax --substrate "1,2,5,10,20" --velocity "0.5,0.8,1.2,1.5,1.7"
  python enzyme_kinetics.py --type hill --substrate "0.1,0.5,1,2,5,10,50" --velocity "0.02,0.1,0.2,0.35,0.6,0.8,0.95"
  python enzyme_kinetics.py --type ki --substrate "1,2,5,10,20" --velocity_no_inh "0.5,0.8,1.2,1.5,1.7" --velocity_inh "0.3,0.5,0.8,1.0,1.1" --inhibitor_conc 5 --inhibition_type competitive
"""

import argparse
import math
import sys


# ---------------------------------------------------------------------------
# Km / Vmax via Lineweaver-Burk linear regression (1/v vs 1/[S])
# ---------------------------------------------------------------------------

def _linreg(xs, ys):
    """Simple linear regression: y = slope*x + intercept. Returns (slope, intercept, r2)."""
    n = len(xs)
    sx = sum(xs)
    sy = sum(ys)
    sxx = sum(x * x for x in xs)
    sxy = sum(x * y for x, y in zip(xs, ys))
    denom = n * sxx - sx * sx
    if abs(denom) < 1e-30:
        raise ValueError("Degenerate data — all x values identical.")
    slope = (n * sxy - sx * sy) / denom
    intercept = (sy * sxx - sx * sxy) / denom
    y_mean = sy / n
    ss_tot = sum((y - y_mean) ** 2 for y in ys)
    ss_res = sum((y - (slope * x + intercept)) ** 2 for x, y in zip(xs, ys))
    r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0
    return slope, intercept, r2


def calc_km_vmax(substrate, velocity):
    """Determine Km and Vmax using Lineweaver-Burk and direct nonlinear grid search."""
    if len(substrate) != len(velocity):
        raise ValueError("substrate and velocity must have the same length.")
    if len(substrate) < 3:
        raise ValueError("Need at least 3 data points.")

    # --- Method 1: Lineweaver-Burk ---
    inv_s = [1.0 / s for s in substrate]
    inv_v = [1.0 / v for v in velocity]
    slope, intercept, r2 = _linreg(inv_s, inv_v)

    if intercept <= 0 or slope <= 0:
        lb_vmax = None
        lb_km = None
        lb_note = "Lineweaver-Burk gave non-physical (negative) parameters — data may not follow Michaelis-Menten."
    else:
        lb_vmax = 1.0 / intercept
        lb_km = slope * lb_vmax
        lb_note = None

    # --- Method 2: Nonlinear least-squares via grid + refinement ---
    vmax_est = max(velocity) * 1.2
    km_est = substrate[len(substrate) // 2]
    best_sse = float("inf")
    best_km, best_vmax = km_est, vmax_est

    for vmax_mult in [x * 0.1 for x in range(5, 31)]:
        for km_mult in [x * 0.1 for x in range(1, 51)]:
            vm = max(velocity) * vmax_mult
            km = max(substrate) * km_mult * 0.1
            sse = sum((v - vm * s / (km + s)) ** 2 for s, v in zip(substrate, velocity))
            if sse < best_sse:
                best_sse = sse
                best_km, best_vmax = km, vm

    # Refine with smaller steps around best
    for _ in range(3):
        step_v = best_vmax * 0.05
        step_k = best_km * 0.05
        improved = False
        for dv in [-step_v, 0, step_v]:
            for dk in [-step_k, 0, step_k]:
                vm = best_vmax + dv
                km = best_km + dk
                if vm <= 0 or km <= 0:
                    continue
                sse = sum((v - vm * s / (km + s)) ** 2 for s, v in zip(substrate, velocity))
                if sse < best_sse:
                    best_sse = sse
                    best_km, best_vmax = km, vm
                    improved = True
        if not improved:
            break

    # R-squared for nonlinear fit
    v_mean = sum(velocity) / len(velocity)
    ss_tot = sum((v - v_mean) ** 2 for v in velocity)
    nl_r2 = 1.0 - best_sse / ss_tot if ss_tot > 0 else 0.0

    print("=" * 60)
    print("MICHAELIS-MENTEN KINETICS: Km and Vmax")
    print("=" * 60)
    print(f"\nData points: {len(substrate)}")
    print(f"[S] range: {min(substrate):.4g} - {max(substrate):.4g}")
    print(f"v range:   {min(velocity):.4g} - {max(velocity):.4g}")

    print("\n--- Lineweaver-Burk (1/v vs 1/[S]) ---")
    if lb_vmax is not None:
        print(f"  Vmax = {lb_vmax:.4g}")
        print(f"  Km   = {lb_km:.4g}")
        print(f"  R²   = {r2:.4f}")
    else:
        print(f"  {lb_note}")

    print("\n--- Nonlinear fit (grid search) ---")
    print(f"  Vmax = {best_vmax:.4g}")
    print(f"  Km   = {best_km:.4g}")
    print(f"  R²   = {nl_r2:.4f}")
    print(f"  SSE  = {best_sse:.4g}")

    print("\n--- Predicted vs Observed ---")
    print(f"  {'[S]':>10s}  {'v_obs':>10s}  {'v_pred':>10s}  {'residual':>10s}")
    for s, v in zip(substrate, velocity):
        v_pred = best_vmax * s / (best_km + s)
        print(f"  {s:10.4g}  {v:10.4g}  {v_pred:10.4g}  {v - v_pred:10.4g}")

    print("\n--- Catalytic efficiency ---")
    print(f"  kcat/Km = Vmax/Km = {best_vmax / best_km:.4g} (units depend on [E] normalization)")
    print("\nNote: For publication-quality fits, use scipy.optimize.curve_fit with proper error estimates.")


# ---------------------------------------------------------------------------
# Hill coefficient from cooperative binding data
# ---------------------------------------------------------------------------

def calc_hill(substrate, velocity):
    """Determine Hill coefficient from log-log linearization of binding data."""
    if len(substrate) != len(velocity):
        raise ValueError("substrate and velocity must have the same length.")
    if len(substrate) < 3:
        raise ValueError("Need at least 3 data points.")

    vmax_est = max(velocity) * 1.1

    # Hill linearization: log(v / (Vmax - v)) = nH * log([S]) - nH * log(K0.5)
    # Use data points where 0.1*Vmax < v < 0.9*Vmax for reliable linearization
    log_s = []
    log_y = []
    for s, v in zip(substrate, velocity):
        if 0.1 * vmax_est < v < 0.9 * vmax_est:
            y = v / (vmax_est - v)
            if y > 0:
                log_s.append(math.log10(s))
                log_y.append(math.log10(y))

    if len(log_s) < 2:
        # Try broader range with adjusted Vmax
        vmax_est = max(velocity) * 1.3
        for s, v in zip(substrate, velocity):
            if 0.05 * vmax_est < v < 0.95 * vmax_est:
                y = v / (vmax_est - v)
                if y > 0:
                    log_s.append(math.log10(s))
                    log_y.append(math.log10(y))

    if len(log_s) < 2:
        print("ERROR: Insufficient data points in the 10-90% saturation range for Hill analysis.")
        print("Provide more data points spanning a wider concentration range.")
        sys.exit(1)

    slope, intercept, r2 = _linreg(log_s, log_y)
    nH = slope
    k05 = 10 ** (-intercept / nH) if abs(nH) > 1e-10 else float("inf")

    print("=" * 60)
    print("HILL ANALYSIS: Cooperative Binding")
    print("=" * 60)
    print(f"\nData points used for Hill plot: {len(log_s)} of {len(substrate)}")
    print(f"Estimated Vmax: {vmax_est:.4g}")
    print(f"\n--- Hill Parameters ---")
    print(f"  Hill coefficient (nH) = {nH:.3f}")
    print(f"  K0.5                  = {k05:.4g}")
    print(f"  R² (Hill plot)        = {r2:.4f}")

    print(f"\n--- Interpretation ---")
    if nH > 1.05:
        print(f"  nH = {nH:.2f} > 1 → POSITIVE cooperativity")
        print(f"  Binding at one site increases affinity at other sites.")
    elif nH < 0.95:
        print(f"  nH = {nH:.2f} < 1 → NEGATIVE cooperativity")
        print(f"  Binding at one site decreases affinity at other sites.")
    else:
        print(f"  nH = {nH:.2f} ≈ 1 → NO cooperativity (independent sites)")

    print(f"\n--- Predicted vs Observed ---")
    print(f"  {'[S]':>10s}  {'v_obs':>10s}  {'v_pred':>10s}")
    for s, v in zip(substrate, velocity):
        v_pred = vmax_est * (s ** nH) / (k05 ** nH + s ** nH)
        print(f"  {s:10.4g}  {v:10.4g}  {v_pred:10.4g}")

    print(f"\nNote: nH is an empirical parameter, not the number of binding sites.")
    print(f"True number of sites must be determined by stoichiometry (ITC, AUC).")


# ---------------------------------------------------------------------------
# Ki from inhibition data
# ---------------------------------------------------------------------------

def calc_ki(substrate, velocity_no_inh, velocity_inh, inhibitor_conc, inhibition_type):
    """Determine Ki from paired velocity data with and without inhibitor."""
    n = len(substrate)
    if len(velocity_no_inh) != n or len(velocity_inh) != n:
        raise ValueError("All data arrays must have the same length.")
    if n < 3:
        raise ValueError("Need at least 3 data points.")

    # First, get Km and Vmax from uninhibited data
    inv_s = [1.0 / s for s in substrate]
    inv_v0 = [1.0 / v for v in velocity_no_inh]
    slope0, intercept0, r2_0 = _linreg(inv_s, inv_v0)

    if intercept0 <= 0 or slope0 <= 0:
        print("ERROR: Uninhibited data does not give valid Km/Vmax. Check data.")
        sys.exit(1)

    vmax = 1.0 / intercept0
    km = slope0 * vmax

    # Get apparent parameters from inhibited data
    inv_vi = [1.0 / v for v in velocity_inh]
    slope_i, intercept_i, r2_i = _linreg(inv_s, inv_vi)

    if intercept_i <= 0 or slope_i <= 0:
        print("ERROR: Inhibited data does not give valid apparent parameters. Check data.")
        sys.exit(1)

    vmax_app = 1.0 / intercept_i
    km_app = slope_i * vmax_app

    inh_type = inhibition_type.lower().strip()
    ki = None

    print("=" * 60)
    print(f"ENZYME INHIBITION ANALYSIS: {inhibition_type.title()}")
    print("=" * 60)
    print(f"\nInhibitor concentration: {inhibitor_conc}")
    print(f"\n--- Uninhibited Parameters ---")
    print(f"  Vmax = {vmax:.4g}")
    print(f"  Km   = {km:.4g}")
    print(f"  R²   = {r2_0:.4f}")
    print(f"\n--- Inhibited Apparent Parameters ---")
    print(f"  Vmax_app = {vmax_app:.4g}")
    print(f"  Km_app   = {km_app:.4g}")
    print(f"  R²       = {r2_i:.4f}")

    if inh_type == "competitive":
        # Competitive: Km_app = Km * (1 + [I]/Ki), Vmax unchanged
        # Ki = [I] / (Km_app/Km - 1)
        ratio = km_app / km
        if ratio <= 1.0:
            print(f"\nWARNING: Km_app ({km_app:.4g}) <= Km ({km:.4g}). Not consistent with competitive inhibition.")
        else:
            ki = inhibitor_conc / (ratio - 1.0)
        print(f"\n--- Competitive Inhibition ---")
        print(f"  Km_app / Km = {ratio:.4f}")
        print(f"  Expected: Vmax_app ≈ Vmax (ratio = {vmax_app / vmax:.3f})")

    elif inh_type == "uncompetitive":
        # Uncompetitive: Vmax_app = Vmax / (1 + [I]/Ki), Km_app = Km / (1 + [I]/Ki)
        ratio = vmax / vmax_app
        if ratio <= 1.0:
            print(f"\nWARNING: Vmax_app ({vmax_app:.4g}) >= Vmax ({vmax:.4g}). Not consistent with uncompetitive inhibition.")
        else:
            ki = inhibitor_conc / (ratio - 1.0)
        print(f"\n--- Uncompetitive Inhibition ---")
        print(f"  Vmax / Vmax_app = {ratio:.4f}")
        print(f"  Km / Km_app = {km / km_app:.4f} (should ≈ {ratio:.4f})")

    elif inh_type in ("noncompetitive", "non-competitive"):
        # Pure noncompetitive: Vmax_app = Vmax / (1 + [I]/Ki), Km unchanged
        ratio = vmax / vmax_app
        if ratio <= 1.0:
            print(f"\nWARNING: Vmax_app ({vmax_app:.4g}) >= Vmax ({vmax:.4g}). Not consistent with noncompetitive inhibition.")
        else:
            ki = inhibitor_conc / (ratio - 1.0)
        print(f"\n--- Noncompetitive Inhibition ---")
        print(f"  Vmax / Vmax_app = {ratio:.4f}")
        print(f"  Expected: Km_app ≈ Km (ratio = {km_app / km:.3f})")

    else:
        print(f"\nERROR: Unknown inhibition type '{inhibition_type}'. Use: competitive, uncompetitive, noncompetitive.")
        sys.exit(1)

    if ki is not None:
        print(f"\n  Ki = {ki:.4g}")
        print(f"\n--- Verification ---")
        print(f"  {'[S]':>10s}  {'v_obs':>10s}  {'v_pred':>10s}  {'residual':>10s}")
        for s, v_obs in zip(substrate, velocity_inh):
            if inh_type == "competitive":
                v_pred = vmax * s / (km * (1 + inhibitor_conc / ki) + s)
            elif inh_type == "uncompetitive":
                factor = 1 + inhibitor_conc / ki
                v_pred = (vmax / factor) * s / (km / factor + s)
            else:  # noncompetitive
                v_pred = (vmax / (1 + inhibitor_conc / ki)) * s / (km + s)
            print(f"  {s:10.4g}  {v_obs:10.4g}  {v_pred:10.4g}  {v_obs - v_pred:10.4g}")
    else:
        print(f"\n  Ki could not be determined — check data consistency with {inhibition_type} model.")

    print(f"\nNote: Lineweaver-Burk analysis amplifies error at low [S]. For publication, use nonlinear regression.")


# ---------------------------------------------------------------------------
# CLI
# ---------------------------------------------------------------------------

def parse_floats(s):
    """Parse comma-separated float string."""
    return [float(x.strip()) for x in s.split(",")]


def main():
    parser = argparse.ArgumentParser(description="Enzyme kinetics calculator")
    parser.add_argument("--type", required=True, choices=["km_vmax", "hill", "ki"],
                        help="Calculation type")
    parser.add_argument("--substrate", help="Comma-separated substrate concentrations")
    parser.add_argument("--velocity", help="Comma-separated velocity values")
    parser.add_argument("--velocity_no_inh", help="Velocities without inhibitor (for ki)")
    parser.add_argument("--velocity_inh", help="Velocities with inhibitor (for ki)")
    parser.add_argument("--inhibitor_conc", type=float, help="Inhibitor concentration (for ki)")
    parser.add_argument("--inhibition_type", default="competitive",
                        help="Type of inhibition: competitive, uncompetitive, noncompetitive")

    args = parser.parse_args()

    if args.type == "km_vmax":
        if not args.substrate or not args.velocity:
            parser.error("--substrate and --velocity required for km_vmax")
        calc_km_vmax(parse_floats(args.substrate), parse_floats(args.velocity))

    elif args.type == "hill":
        if not args.substrate or not args.velocity:
            parser.error("--substrate and --velocity required for hill")
        calc_hill(parse_floats(args.substrate), parse_floats(args.velocity))

    elif args.type == "ki":
        if not args.substrate or not args.velocity_no_inh or not args.velocity_inh:
            parser.error("--substrate, --velocity_no_inh, and --velocity_inh required for ki")
        if args.inhibitor_conc is None:
            parser.error("--inhibitor_conc required for ki")
        calc_ki(
            parse_floats(args.substrate),
            parse_floats(args.velocity_no_inh),
            parse_floats(args.velocity_inh),
            args.inhibitor_conc,
            args.inhibition_type,
        )


if __name__ == "__main__":
    main()

"""Epidemiology calculations: R0/herd immunity, NNT, diagnostic tests, Bayesian post-test probability,
and vaccine coverage threshold from field data.

Consolidates the most commonly needed epidemiology formulas into one script with
verified output. No external dependencies — pure Python stdlib.

Calculation types (--type):

  r0_herd     R0, effective reproduction number, and herd immunity threshold.
              Requires: --R0  (and optionally --VE for vaccine-adjusted threshold,
              --coverage to evaluate Re at a given vaccination fraction)

  vaccine_coverage
              Derive VE from field surveillance data using the screening method
              (Farrington 1993), then compute required vaccination coverage Vc.
              Requires: --R0 --PCV --PPV
              PCV = proportion of cases that were vaccinated (0-1)
              PPV = proportion of population vaccinated (0-1)
              Formula: VE = 1 - [PCV*(1-PPV)] / [(1-PCV)*PPV]
              Then:    Vc = (1 - 1/R0) / VE

  nnt         Number needed to treat / number needed to harm.
              Requires: --control_rate  --treatment_rate

  diagnostic  Sensitivity, specificity, PPV, NPV, accuracy, likelihood ratios.
              Requires: --tp --fp --tn --fn  (2×2 contingency counts)

  bayesian    Post-test probability via Bayes' theorem.
              Requires: --prevalence --sensitivity --specificity
              Optional: --test_result {positive,negative}  (default: positive)

Usage:
    python epidemiology.py --type r0_herd --R0 3.5 --VE 0.90
    python epidemiology.py --type r0_herd --R0 14 --VE 0.97 --coverage 0.95
    python epidemiology.py --type vaccine_coverage --R0 3.0 --PCV 0.06 --PPV 0.35
    python epidemiology.py --type nnt --control_rate 0.30 --treatment_rate 0.20
    python epidemiology.py --type diagnostic --tp 90 --fp 10 --tn 880 --fn 20
    python epidemiology.py --type bayesian --prevalence 0.01 --sensitivity 0.95 --specificity 0.90
    python epidemiology.py --type bayesian --prevalence 0.05 --sensitivity 0.80 --specificity 0.95 --test_result negative
"""

import argparse
import math
import sys


# ---------------------------------------------------------------------------
# R0 / herd immunity
# ---------------------------------------------------------------------------

def r0_herd(R0: float, VE: float = 1.0, coverage: float | None = None) -> dict:
    """
    Compute herd immunity threshold and effective reproduction number.

    Args:
        R0:       Basic reproduction number (must be > 1).
        VE:       Vaccine efficacy as a fraction [0, 1]. Default 1.0 (perfect vaccine).
        coverage: Optional vaccination fraction [0, 1] to evaluate Re at.

    Returns dict with:
        herd_threshold_perfect: Herd immunity threshold assuming 100% efficacy.
        herd_threshold_ve:      VE-adjusted minimum vaccination fraction.
        Re_at_coverage:         Effective R at the given coverage (if provided).
    """
    if R0 <= 1:
        raise ValueError(f"R0 must be > 1 for epidemic spread (got {R0}).")
    if not 0 < VE <= 1:
        raise ValueError(f"VE must be in (0, 1] (got {VE}).")

    hc_perfect = 1.0 - 1.0 / R0
    hc_ve = hc_perfect / VE  # = (1 - 1/R0) / VE

    result = {
        "R0": R0,
        "VE": VE,
        "herd_threshold_perfect": hc_perfect,
        "herd_threshold_ve_adjusted": hc_ve,
    }

    if coverage is not None:
        if not 0 <= coverage <= 1:
            raise ValueError(f"coverage must be in [0, 1] (got {coverage}).")
        # Re = R0 * (1 - VE * coverage)
        Re = R0 * (1.0 - VE * coverage)
        result["coverage"] = coverage
        result["Re_at_coverage"] = Re
        result["epidemic_suppressed"] = Re < 1.0

    return result


def print_r0_herd(res: dict) -> None:
    print("=" * 60)
    print("  R0 and Herd Immunity Analysis")
    print("=" * 60)
    print(f"  Basic reproduction number (R0)  : {res['R0']}")
    print(f"  Vaccine efficacy (VE)           : {res['VE'] * 100:.1f}%")
    print()
    hc = res["herd_threshold_perfect"]
    hc_ve = res["herd_threshold_ve_adjusted"]
    print("  Herd Immunity Threshold (Hc = 1 - 1/R0):")
    print(f"    Perfect vaccine (VE=100%)     : {hc * 100:.2f}%  ({hc:.4f})")
    print()
    print("  VE-adjusted minimum vaccination coverage:")
    print(f"    Vc = Hc / VE = {hc:.4f} / {res['VE']:.2f}")
    print(f"       = {hc_ve:.4f}  ({hc_ve * 100:.2f}%)")
    if hc_ve > 1.0:
        print("    WARNING: VE too low to achieve herd immunity even at 100% coverage.")

    if "coverage" in res:
        cov = res["coverage"]
        Re = res["Re_at_coverage"]
        print()
        print(f"  Effective R (Re) at {cov * 100:.1f}% vaccination coverage:")
        print(f"    Re = R0 × (1 - VE × coverage)")
        print(f"       = {res['R0']} × (1 - {res['VE']:.2f} × {cov:.2f})")
        print(f"       = {Re:.4f}")
        suppressed = res["epidemic_suppressed"]
        status = "SUPPRESSED (Re < 1)" if suppressed else "NOT suppressed (Re ≥ 1)"
        print(f"    Epidemic status: {status}")

    print()
    print("  Verification:")
    check = 1.0 - 1.0 / res["R0"]
    assert abs(check - res["herd_threshold_perfect"]) < 1e-9, "Herd threshold mismatch"
    print(f"    1 - 1/R0 = {check:.6f}  ✓")
    print("=" * 60)


# ---------------------------------------------------------------------------
# Vaccine coverage threshold from field data (screening method)
# ---------------------------------------------------------------------------


def vaccine_coverage(R0: float, PCV: float, PPV: float) -> dict:
    """
    Derive VE from observational field data using the screening method
    (Farrington 1993), then compute the required vaccination coverage Vc.

    The screening method estimates vaccine effectiveness from two readily
    available surveillance numbers:
        PCV = proportion of cases that were vaccinated  [0, 1]
        PPV = proportion of the population vaccinated   [0, 1]

    Formula:
        VE  = 1 - [PCV * (1 - PPV)] / [(1 - PCV) * PPV]
        Hc  = 1 - 1/R0          (herd immunity threshold, perfect vaccine)
        Vc  = Hc / VE           (minimum coverage with real-world VE)

    Args:
        R0:  Basic reproduction number (must be > 1).
        PCV: Proportion of disease cases that are vaccinated (0, 1).
        PPV: Proportion of total population that is vaccinated (0, 1).

    Returns dict with VE, Hc, Vc, and intermediate values.
    """
    if R0 <= 1:
        raise ValueError(f"R0 must be > 1 for epidemic spread (got {R0}).")
    if not 0 < PCV < 1:
        raise ValueError(f"PCV must be in (0, 1) exclusive (got {PCV}).")
    if not 0 < PPV < 1:
        raise ValueError(f"PPV must be in (0, 1) exclusive (got {PPV}).")

    # Screening method: VE = 1 - [PCV*(1-PPV)] / [(1-PCV)*PPV]
    numerator = PCV * (1.0 - PPV)
    denominator = (1.0 - PCV) * PPV
    VE = 1.0 - numerator / denominator

    if VE <= 0:
        raise ValueError(
            f"Derived VE = {VE:.4f} <= 0 — vaccination appears ineffective or harmful "
            f"with PCV={PCV}, PPV={PPV}."
        )

    Hc = 1.0 - 1.0 / R0
    Vc = Hc / VE

    # Re at the current PPV with derived VE
    Re_current = R0 * (1.0 - VE * PPV)

    return {
        "R0": R0,
        "PCV": PCV,
        "PPV": PPV,
        "VE": VE,
        "herd_threshold_perfect": Hc,
        "Vc": Vc,
        "Re_current": Re_current,
    }


def print_vaccine_coverage(res: dict) -> None:
    print("=" * 60)
    print("  Vaccine Coverage Threshold (Screening Method)")
    print("=" * 60)
    print(f"  Basic reproduction number (R0)       : {res['R0']}")
    print(f"  Proportion of cases vaccinated (PCV)  : {res['PCV']}")
    print(f"  Proportion of population vaccinated (PPV): {res['PPV']}")
    print()

    PCV = res["PCV"]
    PPV = res["PPV"]
    VE = res["VE"]
    Hc = res["herd_threshold_perfect"]
    Vc = res["Vc"]

    print("  Step 1 — Derive VE via screening method (Farrington 1993):")
    print(f"    VE = 1 - [PCV * (1-PPV)] / [(1-PCV) * PPV]")
    print(f"       = 1 - [{PCV} * {1-PPV:.4f}] / [{1-PCV:.4f} * {PPV}]")
    num = PCV * (1 - PPV)
    den = (1 - PCV) * PPV
    print(f"       = 1 - {num:.6f} / {den:.6f}")
    print(f"       = 1 - {num/den:.6f}")
    print(f"       = {VE:.4f}  ({VE * 100:.2f}%)")
    print()

    print("  Step 2 — Herd immunity threshold (perfect vaccine):")
    print(f"    Hc = 1 - 1/R0 = 1 - 1/{res['R0']} = {Hc:.4f}  ({Hc * 100:.2f}%)")
    print()

    print("  Step 3 — Required vaccination coverage:")
    print(f"    Vc = Hc / VE = {Hc:.4f} / {VE:.4f}")
    print(f"       = {Vc:.4f}  ({Vc * 100:.1f}%)")

    if Vc > 1.0:
        print("    WARNING: Vc > 100% — herd immunity unachievable at this VE.")
    print()

    Re = res["Re_current"]
    print(f"  Current situation (PPV = {PPV:.0%}):")
    print(f"    Re = R0 * (1 - VE * PPV) = {res['R0']} * (1 - {VE:.4f} * {PPV})")
    print(f"       = {Re:.4f}")
    status = "suppressed" if Re < 1.0 else "NOT suppressed"
    print(f"    Epidemic: {status}")
    print()

    print("  Verification:")
    Re_at_Vc = res["R0"] * (1.0 - VE * Vc)
    print(f"    Re at Vc = R0*(1 - VE*Vc) = {Re_at_Vc:.6f}  (should be ~1.00)")
    VE_check = 1.0 - (PCV * (1 - PPV)) / ((1 - PCV) * PPV)
    assert abs(VE_check - VE) < 1e-9, "VE mismatch"
    print(f"    VE recomputed: {VE_check:.6f}  matches")
    print("=" * 60)


# ---------------------------------------------------------------------------
# NNT / NNH
# ---------------------------------------------------------------------------

def nnt(control_rate: float, treatment_rate: float) -> dict:
    """
    Compute NNT, ARR, RRR, RR, and odds ratio.

    Args:
        control_rate:   Event rate in the control group [0, 1].
        treatment_rate: Event rate in the treatment group [0, 1].

    Returns dict with NNT/NNH, ARR, RRR, RR, OR.
    """
    for name, val in [("control_rate", control_rate), ("treatment_rate", treatment_rate)]:
        if not 0 <= val <= 1:
            raise ValueError(f"{name} must be in [0, 1] (got {val}).")

    ARR = control_rate - treatment_rate  # Absolute Risk Reduction (positive = benefit)
    RR = treatment_rate / control_rate if control_rate > 0 else float("nan")
    RRR = 1.0 - RR if not math.isnan(RR) else float("nan")

    # Odds ratio
    odds_control = control_rate / (1.0 - control_rate) if control_rate < 1 else float("inf")
    odds_treatment = treatment_rate / (1.0 - treatment_rate) if treatment_rate < 1 else float("inf")
    OR = odds_treatment / odds_control if odds_control > 0 else float("nan")

    if abs(ARR) < 1e-12:
        NNT = float("inf")
        label = "NNT"
    elif ARR > 0:
        NNT = 1.0 / ARR
        label = "NNT"  # treatment reduces risk — number needed to treat
    else:
        NNT = 1.0 / abs(ARR)
        label = "NNH"  # treatment increases risk — number needed to harm

    return {
        "control_rate": control_rate,
        "treatment_rate": treatment_rate,
        "ARR": ARR,
        "RR": RR,
        "RRR": RRR,
        "OR": OR,
        "NNT": NNT,
        "NNT_label": label,
    }


def print_nnt(res: dict) -> None:
    print("=" * 60)
    print("  NNT / NNH Analysis")
    print("=" * 60)
    print(f"  Control group event rate   : {res['control_rate'] * 100:.2f}%")
    print(f"  Treatment group event rate : {res['treatment_rate'] * 100:.2f}%")
    print()
    arr = res["ARR"]
    print(f"  Absolute Risk Reduction (ARR)        : {arr * 100:.4f}%")
    rr = res["RR"]
    if not math.isnan(rr):
        print(f"  Relative Risk (RR)                   : {rr:.4f}")
        rrr = res["RRR"]
        print(f"  Relative Risk Reduction (RRR)        : {rrr * 100:.2f}%")
    or_ = res["OR"]
    if not math.isnan(or_) and not math.isinf(or_):
        print(f"  Odds Ratio (OR)                      : {or_:.4f}")
    print()
    label = res["NNT_label"]
    nnt_val = res["NNT"]
    if math.isinf(nnt_val):
        print(f"  {label}: ∞  (treatment has no effect on event rate)")
    else:
        print(f"  {label} = 1 / |ARR| = 1 / {abs(arr):.6f} = {nnt_val:.2f}")
        print(f"  Rounded up: {math.ceil(nnt_val)} patients")
        if label == "NNT":
            print(f"  Interpretation: Treat {math.ceil(nnt_val)} patients to prevent 1 event.")
        else:
            print(f"  Interpretation: Treat {math.ceil(nnt_val)} patients to cause 1 additional event (harm).")
    print()
    print("  Verification:")
    check = abs(1.0 / arr) if abs(arr) > 1e-12 else float("inf")
    if not math.isinf(check):
        assert abs(check - nnt_val) < 1e-6, "NNT mismatch"
    print(f"    1 / |ARR| = {check:.4f}  ✓")
    print("=" * 60)


# ---------------------------------------------------------------------------
# Diagnostic test metrics
# ---------------------------------------------------------------------------

def diagnostic(tp: int, fp: int, tn: int, fn: int) -> dict:
    """
    Compute diagnostic test performance metrics from a 2×2 table.

    Args:
        tp: True positives  (disease+, test+)
        fp: False positives (disease-, test+)
        tn: True negatives  (disease-, test-)
        fn: False negatives (disease+, test-)

    Returns dict with sensitivity, specificity, PPV, NPV, accuracy, LR+, LR-.
    """
    for name, val in [("tp", tp), ("fp", fp), ("tn", tn), ("fn", fn)]:
        if val < 0:
            raise ValueError(f"{name} must be >= 0 (got {val}).")

    n_disease_pos = tp + fn
    n_disease_neg = fp + tn
    n_total = n_disease_pos + n_disease_neg

    if n_disease_pos == 0:
        raise ValueError("No disease-positive cases (tp + fn = 0). Cannot compute sensitivity.")
    if n_disease_neg == 0:
        raise ValueError("No disease-negative cases (fp + tn = 0). Cannot compute specificity.")

    sensitivity = tp / n_disease_pos  # True Positive Rate
    specificity = tn / n_disease_neg  # True Negative Rate
    prevalence = n_disease_pos / n_total

    PPV = tp / (tp + fp) if (tp + fp) > 0 else float("nan")  # Positive Predictive Value
    NPV = tn / (tn + fn) if (tn + fn) > 0 else float("nan")  # Negative Predictive Value
    accuracy = (tp + tn) / n_total

    LR_pos = sensitivity / (1.0 - specificity) if specificity < 1.0 else float("inf")
    LR_neg = (1.0 - sensitivity) / specificity if specificity > 0.0 else float("nan")
    DOR = LR_pos / LR_neg if LR_neg and LR_neg > 0 else float("nan")

    return {
        "tp": tp, "fp": fp, "tn": tn, "fn": fn,
        "n_disease_pos": n_disease_pos,
        "n_disease_neg": n_disease_neg,
        "n_total": n_total,
        "prevalence": prevalence,
        "sensitivity": sensitivity,
        "specificity": specificity,
        "PPV": PPV,
        "NPV": NPV,
        "accuracy": accuracy,
        "LR_pos": LR_pos,
        "LR_neg": LR_neg,
        "DOR": DOR,
    }


def print_diagnostic(res: dict) -> None:
    print("=" * 60)
    print("  Diagnostic Test Performance")
    print("=" * 60)
    print("  2×2 Contingency Table:")
    print(f"    True  Positives (TP) : {res['tp']:>6}")
    print(f"    False Positives (FP) : {res['fp']:>6}")
    print(f"    True  Negatives (TN) : {res['tn']:>6}")
    print(f"    False Negatives (FN) : {res['fn']:>6}")
    print(f"    Total                : {res['n_total']:>6}")
    print(f"    Disease prevalence   : {res['prevalence'] * 100:.2f}%")
    print()
    print("  Performance Metrics:")
    print(f"    Sensitivity (TPR)  = TP / (TP+FN) = {res['tp']} / {res['n_disease_pos']}")
    print(f"                       = {res['sensitivity']:.4f}  ({res['sensitivity'] * 100:.2f}%)")
    print(f"    Specificity (TNR)  = TN / (TN+FP) = {res['tn']} / {res['n_disease_neg']}")
    print(f"                       = {res['specificity']:.4f}  ({res['specificity'] * 100:.2f}%)")
    ppv = res["PPV"]
    npv = res["NPV"]
    n_pos_test = res["tp"] + res["fp"]
    n_neg_test = res["tn"] + res["fn"]
    if not math.isnan(ppv):
        print(f"    PPV                = TP / (TP+FP) = {res['tp']} / {n_pos_test}")
        print(f"                       = {ppv:.4f}  ({ppv * 100:.2f}%)")
    if not math.isnan(npv):
        print(f"    NPV                = TN / (TN+FN) = {res['tn']} / {n_neg_test}")
        print(f"                       = {npv:.4f}  ({npv * 100:.2f}%)")
    print(f"    Accuracy           = (TP+TN) / N  = {res['tp'] + res['tn']} / {res['n_total']}")
    print(f"                       = {res['accuracy']:.4f}  ({res['accuracy'] * 100:.2f}%)")
    print()
    lrp = res["LR_pos"]
    lrn = res["LR_neg"]
    if not math.isinf(lrp):
        print(f"    Likelihood Ratio+  = Sens / (1-Spec)  = {lrp:.4f}")
    if not math.isnan(lrn) and lrn:
        print(f"    Likelihood Ratio-  = (1-Sens) / Spec  = {lrn:.4f}")
    dor = res["DOR"]
    if not math.isnan(dor) and not math.isinf(dor):
        print(f"    Diagnostic Odds Ratio (DOR)           = {dor:.2f}")
    print()
    print("  Interpretation:")
    se = res["sensitivity"]
    sp = res["specificity"]
    if se >= 0.95:
        print(f"    High sensitivity ({se * 100:.1f}%): few missed cases; good rule-out test.")
    elif se < 0.80:
        print(f"    Low sensitivity ({se * 100:.1f}%): many missed cases (false negatives).")
    if sp >= 0.95:
        print(f"    High specificity ({sp * 100:.1f}%): few false alarms; good rule-in test.")
    elif sp < 0.80:
        print(f"    Low specificity ({sp * 100:.1f}%): many false positives.")
    if not math.isnan(ppv):
        print(f"    PPV {ppv * 100:.1f}%: of those testing positive, {ppv * 100:.1f}% truly have the disease.")
    if not math.isnan(npv):
        print(f"    NPV {npv * 100:.1f}%: of those testing negative, {npv * 100:.1f}% truly do not.")
    print()
    print("  Verification:")
    check_se = res["tp"] / (res["tp"] + res["fn"])
    check_sp = res["tn"] / (res["tn"] + res["fp"])
    assert abs(check_se - se) < 1e-9, "Sensitivity mismatch"
    assert abs(check_sp - sp) < 1e-9, "Specificity mismatch"
    print(f"    Sensitivity recomputed: {check_se:.6f}  ✓")
    print(f"    Specificity recomputed: {check_sp:.6f}  ✓")
    print("=" * 60)


# ---------------------------------------------------------------------------
# Bayesian post-test probability
# ---------------------------------------------------------------------------

def bayesian(prevalence: float, sensitivity: float, specificity: float,
             test_result: str = "positive") -> dict:
    """
    Compute post-test probability via Bayes' theorem.

    Args:
        prevalence:  Pre-test probability (prior) of disease [0, 1].
        sensitivity: P(test+ | disease+) [0, 1].
        specificity: P(test- | disease-) [0, 1].
        test_result: "positive" or "negative".

    Returns dict with pre-test odds, LR, post-test odds, and post-test probability.
    """
    for name, val in [("prevalence", prevalence), ("sensitivity", sensitivity),
                      ("specificity", specificity)]:
        if not 0 <= val <= 1:
            raise ValueError(f"{name} must be in [0, 1] (got {val}).")
    if test_result not in ("positive", "negative"):
        raise ValueError("test_result must be 'positive' or 'negative'.")

    pre_test_odds = prevalence / (1.0 - prevalence) if prevalence < 1 else float("inf")

    if test_result == "positive":
        LR = sensitivity / (1.0 - specificity) if specificity < 1.0 else float("inf")
    else:
        LR = (1.0 - sensitivity) / specificity if specificity > 0.0 else float("nan")

    post_test_odds = pre_test_odds * LR
    post_test_prob = post_test_odds / (1.0 + post_test_odds) if not math.isinf(post_test_odds) else 1.0

    # Alternative: direct formula via Bayes
    if test_result == "positive":
        p_test_given_disease = sensitivity
        p_test_given_no_disease = 1.0 - specificity
    else:
        p_test_given_disease = 1.0 - sensitivity
        p_test_given_no_disease = specificity

    p_test = p_test_given_disease * prevalence + p_test_given_no_disease * (1.0 - prevalence)
    ptp_direct = (p_test_given_disease * prevalence / p_test) if p_test > 0 else float("nan")

    return {
        "prevalence": prevalence,
        "sensitivity": sensitivity,
        "specificity": specificity,
        "test_result": test_result,
        "pre_test_odds": pre_test_odds,
        "LR": LR,
        "post_test_odds": post_test_odds,
        "post_test_prob": post_test_prob,
        "post_test_prob_direct": ptp_direct,
        "p_test": p_test,
    }


def print_bayesian(res: dict) -> None:
    print("=" * 60)
    print("  Bayesian Post-Test Probability")
    print("=" * 60)
    print(f"  Pre-test probability (prevalence) : {res['prevalence'] * 100:.2f}%")
    print(f"  Sensitivity                       : {res['sensitivity'] * 100:.2f}%")
    print(f"  Specificity                       : {res['specificity'] * 100:.2f}%")
    print(f"  Test result                       : {res['test_result'].upper()}")
    print()

    pre_odds = res["pre_test_odds"]
    LR = res["LR"]
    post_odds = res["post_test_odds"]
    post_prob = res["post_test_prob"]

    print("  Step 1 — Convert prevalence to pre-test odds:")
    print(f"    Pre-test odds = {res['prevalence']:.4f} / (1 - {res['prevalence']:.4f})")
    print(f"                  = {pre_odds:.4f}")
    print()

    lr_label = "LR+" if res["test_result"] == "positive" else "LR-"
    if res["test_result"] == "positive":
        lr_formula = f"Sensitivity / (1 - Specificity) = {res['sensitivity']:.4f} / {1 - res['specificity']:.4f}"
    else:
        lr_formula = f"(1 - Sensitivity) / Specificity = {1 - res['sensitivity']:.4f} / {res['specificity']:.4f}"

    print(f"  Step 2 — Likelihood Ratio ({lr_label}):")
    print(f"    {lr_label} = {lr_formula}")
    if not math.isnan(LR) and not math.isinf(LR):
        print(f"       = {LR:.4f}")
    elif math.isinf(LR):
        print(f"       = ∞  (perfect specificity)")
    print()

    print("  Step 3 — Post-test odds = pre-test odds × LR:")
    if not math.isinf(post_odds):
        print(f"    = {pre_odds:.4f} × {LR:.4f} = {post_odds:.4f}")
    else:
        print(f"    = ∞")
    print()

    print("  Step 4 — Convert post-test odds to probability:")
    if not math.isinf(post_odds):
        print(f"    Post-test prob = post-test odds / (1 + post-test odds)")
        print(f"                   = {post_odds:.4f} / {1 + post_odds:.4f}")
        print(f"                   = {post_prob:.4f}  ({post_prob * 100:.2f}%)")
    else:
        print(f"    Post-test prob = 100% (certain)")
    print()

    print("  Interpretation:")
    delta = post_prob - res["prevalence"]
    direction = "increased" if delta > 0 else "decreased"
    print(
        f"    A {res['test_result']} test result {direction} disease probability "
        f"from {res['prevalence'] * 100:.2f}% to {post_prob * 100:.2f}%."
    )
    if res["test_result"] == "positive" and post_prob < 0.5:
        print("    PPV < 50%: most positive tests are false positives at this prevalence.")
    elif res["test_result"] == "negative" and (1.0 - post_prob) > 0.99:
        print("    NPV > 99%: a negative result nearly rules out disease.")
    print()
    print("  Verification (direct Bayes formula):")
    direct = res["post_test_prob_direct"]
    if not math.isnan(direct):
        print(f"    P(disease | test {res['test_result']}) via direct formula = {direct:.6f}")
        assert abs(direct - post_prob) < 1e-6, f"Bayesian mismatch: {direct} vs {post_prob}"
        print("    Matches odds-ratio method  ✓")
    print("=" * 60)


# ---------------------------------------------------------------------------
# CLI
# ---------------------------------------------------------------------------

def _build_parser() -> argparse.ArgumentParser:
    p = argparse.ArgumentParser(
        description="Epidemiology calculations: R0/herd immunity, NNT, diagnostics, Bayes.",
        epilog=(
            "Examples:\n"
            "  python epidemiology.py --type r0_herd --R0 3.5 --VE 0.90\n"
            "  python epidemiology.py --type nnt --control_rate 0.30 --treatment_rate 0.20\n"
            "  python epidemiology.py --type diagnostic --tp 90 --fp 10 --tn 880 --fn 20\n"
            "  python epidemiology.py --type bayesian --prevalence 0.01 "
            "--sensitivity 0.95 --specificity 0.90\n"
        ),
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    p.add_argument(
        "--type",
        required=True,
        choices=["r0_herd", "vaccine_coverage", "nnt", "diagnostic", "bayesian"],
        help="Calculation type.",
    )
    # r0_herd
    p.add_argument("--R0", type=float, help="Basic reproduction number (r0_herd).")
    p.add_argument("--VE", type=float, default=1.0, help="Vaccine efficacy [0,1] (r0_herd).")
    p.add_argument("--coverage", type=float, help="Vaccination coverage fraction to evaluate Re (r0_herd).")
    # vaccine_coverage (screening method)
    p.add_argument("--PCV", type=float, help="Proportion of cases vaccinated (vaccine_coverage).")
    p.add_argument("--PPV", type=float, help="Proportion of population vaccinated (vaccine_coverage).")
    # nnt
    p.add_argument("--control_rate", type=float, help="Event rate in control group (nnt).")
    p.add_argument("--treatment_rate", type=float, help="Event rate in treatment group (nnt).")
    # diagnostic
    p.add_argument("--tp", type=int, help="True positives (diagnostic).")
    p.add_argument("--fp", type=int, help="False positives (diagnostic).")
    p.add_argument("--tn", type=int, help="True negatives (diagnostic).")
    p.add_argument("--fn", type=int, help="False negatives (diagnostic).")
    # bayesian
    p.add_argument("--prevalence", type=float, help="Pre-test probability / disease prevalence (bayesian).")
    p.add_argument("--sensitivity", type=float, help="Test sensitivity (bayesian, diagnostic).")
    p.add_argument("--specificity", type=float, help="Test specificity (bayesian, diagnostic).")
    p.add_argument(
        "--test_result",
        choices=["positive", "negative"],
        default="positive",
        help="Test result to condition on (bayesian, default: positive).",
    )
    return p


def main() -> None:
    parser = _build_parser()
    args = parser.parse_args()

    calc_type = args.type

    try:
        if calc_type == "r0_herd":
            if args.R0 is None:
                parser.error("--type r0_herd requires --R0.")
            res = r0_herd(args.R0, VE=args.VE, coverage=args.coverage)
            print_r0_herd(res)

        elif calc_type == "vaccine_coverage":
            if args.R0 is None:
                parser.error("--type vaccine_coverage requires --R0.")
            if args.PCV is None or args.PPV is None:
                parser.error("--type vaccine_coverage requires --PCV and --PPV.")
            res = vaccine_coverage(args.R0, args.PCV, args.PPV)
            print_vaccine_coverage(res)

        elif calc_type == "nnt":
            if args.control_rate is None or args.treatment_rate is None:
                parser.error("--type nnt requires --control_rate and --treatment_rate.")
            res = nnt(args.control_rate, args.treatment_rate)
            print_nnt(res)

        elif calc_type == "diagnostic":
            for flag in ("--tp", "--fp", "--tn", "--fn"):
                attr = flag.lstrip("-")
                if getattr(args, attr) is None:
                    parser.error(f"--type diagnostic requires {flag}.")
            res = diagnostic(args.tp, args.fp, args.tn, args.fn)
            print_diagnostic(res)

        elif calc_type == "bayesian":
            for flag, attr in [("--prevalence", "prevalence"),
                                ("--sensitivity", "sensitivity"),
                                ("--specificity", "specificity")]:
                if getattr(args, attr) is None:
                    parser.error(f"--type bayesian requires {flag}.")
            res = bayesian(args.prevalence, args.sensitivity, args.specificity, args.test_result)
            print_bayesian(res)

    except ValueError as exc:
        print(f"ERROR: {exc}", file=sys.stderr)
        sys.exit(1)


if __name__ == "__main__":
    main()