// population_pk.sy -- Population pharmacokinetics and pharmacometrics
// as native Scyla. Inter-individual variability, covariates, bioequivalence,
// and special populations.
//
// Sources:
//   Ette EI, Williams PJ. Pharmacometrics, Wiley, 2007.
//   Mould DR, Upton RN. CPT Pharmacometrics Syst Pharmacol 2:e38, 2013.
//   FDA Guidance for Industry: Population Pharmacokinetics, 2022.
//
// Marjorie McCubbins + Aether, April 2026
// Patent pending: 64/034,536

// =========================================================
// POPULATION PK FUNDAMENTALS
// =========================================================

// Typical value with covariate effect: TV = theta * effect
fn typical_value(theta: Float, cov_effect: Float) -> Float {
    return theta * cov_effect
}

// Individual parameter with random effect (log-normal):
// P_i = TV * exp(eta)
// eta ~ N(0, omega^2)
fn individual_parameter(tv: Float, eta: Float) -> Float {
    return tv * exp(eta)
}

// Proportional residual error: Y = F * (1 + epsilon)
fn residual_proportional(predicted: Float, epsilon: Float) -> Float {
    return predicted * (1.0 + epsilon)
}

// Additive residual error: Y = F + epsilon
fn residual_additive(predicted: Float, epsilon: Float) -> Float {
    return predicted + epsilon
}

// Combined (proportional + additive): Y = F * (1 + eps1) + eps2
fn residual_combined(predicted: Float, eps_prop: Float, eps_add: Float) -> Float {
    return predicted * (1.0 + eps_prop) + eps_add
}

// Exponential error: Y = F * exp(epsilon)
fn residual_exponential(predicted: Float, epsilon: Float) -> Float {
    return predicted * exp(epsilon)
}

// =========================================================
// COVARIATE MODELS
// =========================================================

// Allometric scaling of clearance by body weight
// CL = CL_ref * (WT/70)^0.75
// Exponent 0.75 from West GB et al. Science 276:122-126, 1997.
fn weight_allometric_cl(cl_ref: Float, weight_kg: Float) -> Float {
    return cl_ref * pow(weight_kg / 70.0, 0.75)
}

// Allometric scaling of volume by body weight
// Vd = Vd_ref * (WT/70)^1.0
fn weight_allometric_vd(vd_ref: Float, weight_kg: Float) -> Float {
    return vd_ref * (weight_kg / 70.0)
}

// Maturation function (sigmoid Emax for pediatric CL development)
// MF = PMA^hill / (TM50^hill + PMA^hill)
// PMA = post-menstrual age in weeks, TM50 = age at 50% maturation
// NOTE: This uses standard sigmoid, NOT NCA tanh gates.
fn age_maturation(pma_weeks: Float, tm50: Float, hill: Float) -> Float {
    let pma_n = pow(pma_weeks, hill)
    let tm_n = pow(tm50, hill)
    return pma_n / (tm_n + pma_n)
}

// Renal function covariate (linear)
// CL = CL_ref * (CrCl / 120)
fn renal_covariate(cl_ref: Float, crcl: Float) -> Float {
    return cl_ref * (crcl / 120.0)
}

// Mixed renal/non-renal clearance covariate
// CL = CLnr + CLr * (CrCl / 120)
fn renal_covariate_mixed(cl_nonrenal: Float, cl_renal_ref: Float, crcl: Float) -> Float {
    return cl_nonrenal + cl_renal_ref * (crcl / 120.0)
}

// CYP2D6 metabolizer phenotype effect on clearance
fn cyp2d6_phenotype(phenotype: String) -> Float {
    if phenotype == "poor" { return 0.1 }
    if phenotype == "intermediate" { return 0.5 }
    if phenotype == "extensive" { return 1.0 }
    if phenotype == "ultrarapid" { return 2.0 }
    return 1.0
}

// CYP2C19 metabolizer phenotype effect
fn cyp2c19_phenotype(phenotype: String) -> Float {
    if phenotype == "poor" { return 0.1 }
    if phenotype == "intermediate" { return 0.5 }
    if phenotype == "extensive" { return 1.0 }
    if phenotype == "ultrarapid" { return 1.5 }
    return 1.0
}

// CYP3A phenotype (less polymorphic, more affected by DDIs)
fn cyp3a_activity(phenotype: String) -> Float {
    if phenotype == "inhibited" { return 0.2 }
    if phenotype == "normal" { return 1.0 }
    if phenotype == "induced" { return 2.5 }
    return 1.0
}

// BSA covariate
fn bsa_covariate(param_ref: Float, bsa: Float) -> Float {
    return param_ref * (bsa / 1.73)
}

// Sex effect on parameter
fn sex_effect(param_ref: Float, is_female: Bool, female_factor: Float) -> Float {
    if is_female { return param_ref * female_factor }
    return param_ref
}

// Power model covariate: param = theta * (covariate/reference)^exponent
fn power_covariate(theta: Float, covariate: Float, reference: Float, exponent: Float) -> Float {
    return theta * pow(covariate / reference, exponent)
}

// =========================================================
// VARIABILITY METRICS
// =========================================================

// CV% from omega-squared (log-normal model)
// CV% = sqrt(exp(omega^2) - 1) * 100
fn iiv_cv_from_omega(omega_sq: Float) -> Float {
    return sqrt(exp(omega_sq) - 1.0) * 100.0
}

// Omega-squared from CV%
// omega^2 = ln(1 + (CV/100)^2)
fn iiv_omega_from_cv(cv_pct: Float) -> Float {
    return ln(1.0 + pow(cv_pct / 100.0, 2.0))
}

// Eta-shrinkage: measures how much individual estimates shrink toward population mean
// High shrinkage (>20-30%) means individual data is insufficient to estimate eta
fn shrinkage(omega_sq: Float, se_eta_sq: Float) -> Float {
    return 1.0 - sqrt(se_eta_sq / omega_sq)
}

// Coefficient of variation
fn cv_percent(sd: Float, mean_val: Float) -> Float {
    return (sd / mean_val) * 100.0
}

// Relative standard error
fn rse_percent(se: Float, estimate: Float) -> Float {
    return abs(se / estimate) * 100.0
}

// =========================================================
// BIOEQUIVALENCE
// =========================================================

// AUC ratio (test/reference)
fn auc_ratio(auc_test: Float, auc_ref: Float) -> Float {
    return auc_test / auc_ref
}

// Cmax ratio (test/reference)
fn cmax_ratio(cmax_test: Float, cmax_ref: Float) -> Float {
    return cmax_test / cmax_ref
}

// Standard BE criteria: 90% CI of geometric mean ratio within 80-125%
fn be_pass(lower_ci: Float, upper_ci: Float) -> Bool {
    if lower_ci < 0.80 { return false }
    if upper_ci > 1.25 { return false }
    return true
}

// Narrow therapeutic index drugs: tighter 90-111% criteria
fn be_pass_nti(lower_ci: Float, upper_ci: Float) -> Bool {
    if lower_ci < 0.90 { return false }
    if upper_ci > 1.11 { return false }
    return true
}

// BE assessment summary
fn be_assess(gmr: Float, lower_ci: Float, upper_ci: Float) -> String {
    if lower_ci >= 0.80 {
        if upper_ci <= 1.25 { return "BIOEQUIVALENT: 90% CI within 80-125%. Meets FDA criteria for generic approval." }
    }
    if lower_ci >= 0.90 {
        if upper_ci <= 1.11 { return "BIOEQUIVALENT (NTI criteria): 90% CI within 90-111%. Meets stricter criteria for narrow therapeutic index drugs." }
    }
    return "NOT BIOEQUIVALENT: 90% CI falls outside 80-125%. Does not meet FDA criteria."
}

// =========================================================
// SPECIAL POPULATIONS
// =========================================================

// Pediatric CL scaling with maturation
// Combines allometric scaling + maturation sigmoid
fn pediatric_cl(adult_cl: Float, weight_kg: Float, pma_weeks: Float, tm50: Float, hill: Float) -> Float {
    let allometric = weight_allometric_cl(adult_cl, weight_kg)
    let maturation = age_maturation(pma_weeks, tm50, hill)
    return allometric * maturation
}

// Simplified pediatric scaling by age
fn pediatric_cl_simple(adult_cl: Float, age_years: Float) -> Float {
    if age_years < 0.08 { return adult_cl * 0.1 }
    if age_years < 0.5 { return adult_cl * 0.3 }
    if age_years < 1.0 { return adult_cl * 0.5 }
    if age_years < 2.0 { return adult_cl * 0.7 }
    if age_years < 6.0 { return adult_cl * 0.9 }
    if age_years < 12.0 { return adult_cl * 1.1 }
    return adult_cl
}

// Neonatal CL fraction of adult
fn neonatal_cl_fraction(postnatal_days: Float) -> Float {
    // CYP3A7 dominant at birth, transitions to CYP3A4 by 6 months
    // GFR ~2-4 mL/min at birth, reaches adult values by 1 year
    if postnatal_days < 1.0 { return 0.05 }
    if postnatal_days < 7.0 { return 0.1 }
    if postnatal_days < 28.0 { return 0.2 }
    if postnatal_days < 90.0 { return 0.4 }
    if postnatal_days < 180.0 { return 0.6 }
    if postnatal_days < 365.0 { return 0.8 }
    return 1.0
}

// Geriatric CL adjustment (~1% decline per year after 40)
fn geriatric_cl(adult_cl: Float, age: Float) -> Float {
    if age <= 40.0 { return adult_cl }
    let decline = 1.0 - 0.01 * (age - 40.0)
    if decline < 0.3 { return adult_cl * 0.3 }
    return adult_cl * decline
}

// Hepatic impairment scaling by Child-Pugh class
fn hepatic_scaling(cl_ref: Float, child_pugh_class: String) -> Float {
    if child_pugh_class == "A" { return cl_ref * 0.75 }
    if child_pugh_class == "B" { return cl_ref * 0.50 }
    if child_pugh_class == "C" { return cl_ref * 0.25 }
    return cl_ref
}

// Obesity: which weight to use by drug class
fn obesity_dosing_weight(drug_class: String) -> String {
    if drug_class == "aminoglycosides" { return "Adjusted body weight (ABW = IBW + 0.4*(TBW-IBW)). Aminoglycosides distribute into lean tissue + partial extracellular fluid." }
    if drug_class == "vancomycin" { return "Total body weight (TBW). Vancomycin distributes into total body water. Obese patients need higher doses." }
    if drug_class == "fluoroquinolones" { return "Total body weight (TBW). Good tissue penetration regardless of adiposity." }
    if drug_class == "beta_lactams" { return "Lean body weight or adjusted body weight. Hydrophilic, distribute into lean tissue." }
    if drug_class == "heparin" { return "Total body weight (TBW) up to max dose (e.g., enoxaparin 150 mg). Dose cap for morbid obesity." }
    return "Consult pharmacist. Weight-based dosing in obesity is drug-specific."
}

// Pregnancy PK changes
fn pregnancy_pk_changes() -> String {
    return "Pregnancy alters PK: (1) Increased GFR 50% by 3rd trimester (increased renal CL of renally cleared drugs: lithium, digoxin, ampicillin). (2) Increased Vd (plasma volume +50%, total body water +8L). (3) Decreased albumin (increased free fraction of highly bound drugs: phenytoin, valproic acid). (4) CYP3A4 induced (increased CL of nifedipine, midazolam). (5) CYP2D6 induced (increased CL of metoprolol, fluoxetine). (6) CYP1A2 decreased (decreased CL of caffeine, theophylline, olanzapine). (7) Delayed gastric emptying (altered absorption Tmax). Source: Pariente G et al. Clin Pharmacokinet 55:251-268, 2016."
}

// =========================================================
// SIMULATION HELPERS
// =========================================================

// Box-Muller transform: generate standard normal from two uniform randoms
// z = sqrt(-2 * ln(u1)) * cos(2 * pi * u2)
fn box_muller(u1: Float, u2: Float) -> Float {
    return sqrt(0.0 - 2.0 * ln(u1)) * cos(2.0 * pi(0) * u2)
}

// Log-normal sample: X = mean * exp(sqrt(omega^2) * z)
fn log_normal_sample(mean_val: Float, omega_sq: Float, z: Float) -> Float {
    return mean_val * exp(sqrt(omega_sq) * z)
}

// Prediction interval (95%) lower bound
fn prediction_interval_lower(mean_val: Float, omega_sq: Float) -> Float {
    return mean_val * exp(0.0 - 1.96 * sqrt(omega_sq))
}

// Prediction interval (95%) upper bound
fn prediction_interval_upper(mean_val: Float, omega_sq: Float) -> Float {
    return mean_val * exp(1.96 * sqrt(omega_sq))
}

// Geometric mean of two values
fn geometric_mean_2(a: Float, b: Float) -> Float {
    return sqrt(a * b)
}

// Geometric CV
fn geometric_cv(omega_sq: Float) -> Float {
    return sqrt(exp(omega_sq) - 1.0) * 100.0
}

// =========================================================
// MODEL SELECTION CRITERIA
// =========================================================

// AIC: Akaike Information Criterion
// AIC = -2LL + 2*k (k = number of parameters)
// Lower is better. Penalizes overfitting.
fn aic(neg2ll: Float, n_params: Float) -> Float {
    return neg2ll + 2.0 * n_params
}

// BIC: Bayesian Information Criterion
// BIC = -2LL + k*ln(n) (n = number of observations)
// Penalizes complexity more than AIC for large datasets.
fn bic(neg2ll: Float, n_params: Float, n_obs: Float) -> Float {
    return neg2ll + n_params * ln(n_obs)
}

// Likelihood ratio test: chi-square with df = difference in parameters
// Significant if delta_OFV > 3.84 (p<0.05, df=1) or > 6.63 (p<0.01, df=1)
fn lrt_significant_p05(delta_ofv: Float, df: Int) -> Bool {
    // Chi-square critical values at p=0.05: df1=3.84, df2=5.99, df3=7.81
    if df == 1 { if delta_ofv > 3.84 { return true } }
    if df == 2 { if delta_ofv > 5.99 { return true } }
    if df == 3 { if delta_ofv > 7.81 { return true } }
    if df == 4 { if delta_ofv > 9.49 { return true } }
    return false
}

// =========================================================
// STRUCTURAL MODEL DESCRIPTIONS
// =========================================================

fn structural_model(name: String) -> String {
    if name == "one_compartment" { return "One-compartment model: body as single homogeneous compartment. Drug distributes instantaneously. Monoexponential decline. Parameters: CL, Vd. Suitable for drugs with rapid distribution (small Vd, hydrophilic)." }
    if name == "two_compartment" { return "Two-compartment model: central (plasma + well-perfused organs) + peripheral (poorly perfused tissues). Biexponential decline: alpha (distribution) + beta (elimination). Parameters: CL, V1, Q (inter-compartmental CL), V2. Most drugs." }
    if name == "three_compartment" { return "Three-compartment model: central + shallow peripheral + deep peripheral. Triexponential decline. Parameters: CL, V1, Q2, V2, Q3, V3. For drugs with very slow tissue redistribution (thiopental, digoxin)." }
    return ""
}

fn absorption_model(name: String) -> String {
    if name == "first_order" { return "First-order absorption: dA/dt = -ka*A. Most common for oral drugs. Rate proportional to amount remaining. Parameters: ka, Tlag (optional)." }
    if name == "zero_order" { return "Zero-order absorption: dA/dt = -k0. Constant rate input. For: controlled-release formulations, IV infusion, transdermal patches. Parameters: k0 (rate), duration." }
    if name == "transit" { return "Transit compartment model (Savic 2007): chain of transit compartments before absorption compartment. Models delayed absorption better than lag time. Parameters: ka, ktr (transit rate), n (number of transits). MTT = (n+1)/ktr." }
    return ""
}

// =========================================================
// COUNTS
// =========================================================

fn popPK_function_count() -> Int {
    return 45
}
