#!/usr/bin/env python3
"""TEST — plasma colloid osmotic pressure by Landis-Pappenheimer, with globulins.

WRITTEN BEFORE THE CHANGE, and it FAILS on the pre-change engine — that is the point.
The pre-change engine modelled plasma protein as albumin alone under van't Hoff, which is
LINEAR in concentration and therefore materially underestimates colloid osmotic pressure:
it produced a transperitoneal protein oncotic gradient of 11.2 mmHg against the 22 mmHg
measured in human PD patients (Öberg CM & Rippe B, Kidney Int Rep 2017;2(6):1128-38,
PMC5733752, Table 1: "Transperitoneal oncotic pressure gradient (Delta-pi-prot) 22 mm Hg").

TWO separate defects produced that gap and this file asserts both are gone:
  1. NO GLOBULINS. Globulins are ~40% of plasma protein mass and the model had none.
  2. VAN'T HOFF. Plasma protein is strongly non-ideal at physiologic concentration; the
     osmotic pressure rises faster than linearly. The standard empirical replacement is
     Landis EM & Pappenheimer JR, Handbook of Physiology, Sect 2 Vol II, 1963, p.975:
         COP (mmHg) = 2.1*C + 0.16*C^2 + 0.009*C^3,  C = TOTAL plasma protein in g/dL
     At a normal total protein of 7.0 g/dL that is 25.6 mmHg, against an ideal van't Hoff
     of ~15 — the non-ideal terms contribute nearly half the real pressure.

Note on what is NOT asserted here: the value 22 mmHg is Öberg's TRANSPERITONEAL GRADIENT,
quoted in that table as a raw colloid osmotic pressure and NOT pre-multiplied by a
reflection coefficient (their Starling term is RT*sum(phi_i*(C_P,i - C_D,i)), so phi is
applied separately). So the raw gradient is what is checked against 22 here, and the
engine then applies its own sigma on top, exactly as Rippe's equation does.

Run: OMP_NUM_THREADS=1 python3 test_oncotic_landis.py
"""
import os
os.environ.setdefault('OMP_NUM_THREADS', '1')
os.environ.setdefault('PARALLEL', '1')

import sys

import physics_core as PC
from physics_core import SOLUTES, REFLECTION, R_T, Pool
import physics_sim as PS

FAIL = []


def check(cond, name, detail=''):
    print('  %-4s %-52s %s' % ('OK' if cond else 'FAIL', name, detail))
    if not cond:
        FAIL.append(name)


print('=' * 100)
print('LANDIS-PAPPENHEIMER COLLOID OSMOTIC PRESSURE — with globulins')
print('=' * 100)

# ---------------------------------------------------------------- 1. the law itself
print('\n1. THE LAW  COP = 2.1C + 0.16C^2 + 0.009C^3   (C = g/dL total plasma protein)')
cop = PC.colloid_osmotic_pressure_g_dl
check(abs(cop(0.0)) < 1e-12, 'zero protein gives zero pressure', '%.3f mmHg' % cop(0.0))
c70 = cop(7.0)
check(25.0 < c70 < 26.5, 'normal total protein 7.0 g/dL -> ~25.6 mmHg',
      '%.2f mmHg (textbook plasma COP 25-28)' % c70)
# The published anchor, computed by hand so the test does not just re-run the implementation:
hand = 2.1 * 7.0 + 0.16 * 49.0 + 0.009 * 343.0
check(abs(c70 - hand) < 1e-9, 'matches the hand-evaluated cubic', '%.4f vs %.4f' % (c70, hand))
mono = all(cop(c) < cop(c + 0.1) for c in [x * 0.1 for x in range(0, 100)])
check(mono, 'monotone increasing over 0-10 g/dL')
# The cubic must NOT be extrapolated past its validated range — see physics_core. A dehydrating
# patient haemoconcentrates, and the adversary reached 11.97 g/dL in a dying 30-day arm.
PC.COP_EXTRAPOLATIONS[0] = 0
edge = PC.COP_VALID_MAX_G_DL
check(abs(cop(edge) - PC._cop_cubic(edge)) < 1e-12, 'the guard is exactly the cubic AT the boundary',
      '%.4f mmHg at %.1f g/dL' % (cop(edge), edge))
check(PC.COP_EXTRAPOLATIONS[0] == 0, 'and evaluating AT the boundary is not an extrapolation')
_h = 1e-6
_sl_in = (cop(edge) - cop(edge - _h)) / _h
_sl_out = (cop(edge + _h) - cop(edge)) / _h
check(abs(_sl_in - _sl_out) < 1e-3, 'C1-continuous across the boundary — no kink in the force',
      'slope %.6f in, %.6f out (mmHg per g/dL)' % (_sl_in, _sl_out))
check(cop(20.0) < PC._cop_cubic(20.0), 'the tangent extension is BELOW the runaway cubic',
      '%.1f vs the cubic\'s %.1f mmHg at 20 g/dL — the cubic asserts that on no evidence'
      % (cop(20.0), PC._cop_cubic(20.0)))
check(PC.COP_EXTRAPOLATIONS[0] > 0, 'and extrapolations are COUNTED, never silently swallowed',
      '%d recorded' % PC.COP_EXTRAPOLATIONS[0])
check(all(cop(c) < cop(c + 0.1) for c in [x * 0.5 for x in range(0, 60)]),
      'still monotone across the boundary and out to 30 g/dL')
PC.COP_EXTRAPOLATIONS[0] = 0
# Non-ideality must be POSITIVE and large: van't Hoff on the same protein MASS underestimates.
# The reference has to be the real plasma protein MIXTURE, not albumin alone. Van't Hoff counts
# MOLES, so pretending all 7 g/dL were the smallest protein maximises the mole count and therefore
# flatters the ideal law — my first version of this check did exactly that (66.5 kDa throughout,
# giving 20.4 mmHg) and failed a comparison that was wrong in the test, not in the engine. The
# physiologic 60/40 albumin/globulin split has a mass-weighted mean MW of ~85 kDa.
_f_alb = 0.60
_mw_mix = 1.0 / (_f_alb / PC.PROTEIN_MW_G_PER_MMOL['Albumin']
                 + (1.0 - _f_alb) / PC.PROTEIN_MW_G_PER_MMOL['Globulin'])
vh = R_T * (7.0 * 10.0 / _mw_mix)
check(c70 > vh * 1.5, 'non-ideal COP far exceeds ideal van\'t Hoff',
      'LP %.2f vs van\'t Hoff %.2f mmHg on the same mass (mean MW %.1f kDa)' % (c70, vh, _mw_mix))
# Öberg's 22 mmHg must correspond to a PLAUSIBLE PD-patient total protein, not a fantasy one.
lo, hi = 0.0, 12.0
for _ in range(80):
    mid = 0.5 * (lo + hi)
    if cop(mid) < 22.0:
        lo = mid
    else:
        hi = mid
c22 = 0.5 * (lo + hi)
check(5.8 < c22 < 7.0, 'Oberg 22 mmHg <-> a real PD total protein',
      '%.2f g/dL (PD patients run 6.0-6.8, mildly hypoproteinaemic)' % c22)

# ---------------------------------------------------------------- 2. globulins exist
print('\n2. GLOBULINS ARE MODELLED')
check('Globulin' in SOLUTES, 'Globulin is a tracked solute species')
check(PC.VALENCE.get('Globulin', None) == 0, 'globulin carries no net charge at pH 7.4',
      'Figge-Fencl models albumin and phosphate only; globulin pI spread straddles 7.4')
check('Globulin' not in PC.LAMBDA_ION, 'globulin contributes nothing to conductivity',
      'it is uncharged, so the device cannot see it')
check(REFLECTION['Globulin'] > REFLECTION['Albumin'], 'globulin is reflected MORE than albumin',
      'sigma %.2f vs %.2f — 150 kDa vs 66.5 kDa' % (REFLECTION['Globulin'], REFLECTION['Albumin']))
check(PC.SIEVING['Globulin'] < PC.SIEVING['Albumin'], 'globulin sieves LESS than albumin',
      'S %.3f vs %.3f' % (PC.SIEVING['Globulin'], PC.SIEVING['Albumin']))

sim = PS.PhysicsSim(weight=72.0, N=10)
alb_g = sim.serum.conc('Albumin', sim.serum_vd['Albumin']) * PC.PROTEIN_MW_G_PER_MMOL['Albumin'] / 10.0
glb_g = sim.serum.conc('Globulin', sim.serum_vd['Globulin']) * PC.PROTEIN_MW_G_PER_MMOL['Globulin'] / 10.0
check(2.0 < glb_g < 3.6, 'serum globulin is physiologic', '%.2f g/dL (reference 2.0-3.5)' % glb_g)
check(5.8 < alb_g + glb_g < 8.2, 'serum TOTAL protein is physiologic',
      '%.2f g/dL = albumin %.2f + globulin %.2f (reference 6.0-8.3)' % (alb_g + glb_g, alb_g, glb_g))
check(0.50 < alb_g / (alb_g + glb_g) < 0.68, 'albumin fraction of total protein is physiologic',
      '%.1f%% (reference 55-65%%)' % (100.0 * alb_g / (alb_g + glb_g)))

# ---------------------------------------------------------------- 3. the gradient
print('\n3. THE TRANSPERITONEAL PROTEIN ONCOTIC GRADIENT')
for _ in range(60):
    sim.step(1.0)
lay = sim.layers[0]
raw = PC.protein_cop(sim.serum, sim.serum_vd) - PC.protein_cop(lay, None)
check(17.0 <= raw <= 28.0, 'raw gradient brackets Oberg 22 mmHg',
      '%.2f mmHg (engine protein %.2f g/dL; was 13.4 pre-change, albumin-only van\'t Hoff)'
      % (raw, alb_g + glb_g))
sig = PC.protein_sigma(sim.serum, sim.serum_vd)
check(0.84 <= sig <= 0.90, 'effective protein sigma is the mass-weighted mix',
      '%.4f between albumin %.2f and globulin %.2f' % (sig, REFLECTION['Albumin'], REFLECTION['Globulin']))
# and it must OPPOSE ultrafiltration (negative contribution to dpi = layer - serum)
dpi_prot = PC.protein_oncotic_diff(lay, sim.serum, sim.serum_vd)
check(dpi_prot < 0.0, 'protein oncotic term OPPOSES ultrafiltration', '%+.2f mmHg' % dpi_prot)
check(abs(dpi_prot + sig * raw) < 1e-9, 'the term is exactly sigma x the raw gradient',
      '%.6f vs %.6f' % (-dpi_prot, sig * raw))
# NOT a regression on the old number: it must be materially BIGGER than the 11.2 it replaces.
check(abs(dpi_prot) > 15.0, 'strictly larger than the albumin-only van\'t Hoff term it replaces',
      '%.2f mmHg vs the pre-change 11.21' % abs(dpi_prot))

# ---------------------------------------------------------------------- 3b. THE REAL TEST
# The check above is a one-hour snapshot of a patient who still has a HEALTHY person's protein
# (7.46 g/dL), and Öberg's 22 mmHg is measured in PD PATIENTS, who are mildly hypoproteinaemic.
# So the claim worth asserting is the STEADY STATE one: put the patient on the treatment, let
# peritoneal protein loss do what it does, and the gradient must land on the measured value.
# Nothing here is fitted to 22 — globulin was set to its normal reference concentration and the
# transport coefficients come from the size argument, not from this number.
print('\n3b. STEADY STATE — the gradient a PD PATIENT actually runs at (Oberg measured 22 mmHg)')
for _ in range(4 * 1440):
    sim.step(1.0)
alb_ss = sim.serum.conc('Albumin', sim.serum_vd['Albumin']) * PC.PROTEIN_MW_G_PER_MMOL['Albumin'] / 10.0
glb_ss = sim.serum.conc('Globulin', sim.serum_vd['Globulin']) * PC.PROTEIN_MW_G_PER_MMOL['Globulin'] / 10.0
raw_ss = PC.protein_cop(sim.serum, sim.serum_vd) - PC.protein_cop(sim.layers[0], None)
check(5.9 <= alb_ss + glb_ss <= 6.9, 'total protein settles in the PD-patient range',
      '%.2f g/dL = albumin %.2f + globulin %.2f (started 7.46; PD reference 6.0-6.8)'
      % (alb_ss + glb_ss, alb_ss, glb_ss))
check(19.0 <= raw_ss <= 25.0, 'STEADY-STATE gradient reproduces Oberg\'s measured 22 mmHg',
      '%.2f mmHg — this is the headline claim' % raw_ss)
check(alb_ss < alb_g, 'peritoneal dialysis lowered serum albumin, as it does clinically',
      '%.2f -> %.2f g/dL' % (alb_g, alb_ss))
check(glb_ss < glb_g, 'and lowered globulin too — it is genuinely lost to the effluent',
      '%.2f -> %.2f g/dL' % (glb_g, glb_ss))
check(PC.COP_EXTRAPOLATIONS[0] == 0,
      'no out-of-range extrapolation anywhere in 4 days of normal operation',
      '%d evaluations above %.0f g/dL' % (PC.COP_EXTRAPOLATIONS[0], PC.COP_VALID_MAX_G_DL))

# ---------------------------------------------------------------- 4. the two code paths agree
print('\n4. THE TWO CODE PATHS AGREE (physics_core vs the inline loop in physics_sim)')
sc = {s: sim.serum.conc(s, sim.serum_vd[s]) for s in SOLUTES}
lc = {s: lay.conc(s) for s in SOLUTES}
core = PC.osmotic_pressure_diff(lay, sim.serum, sim.serum_vd)
inline = PC.osmotic_pressure_diff_conc(lc, sc)
check(abs(core - inline) < 1e-12, 'pool-based and concentration-based forms are identical',
      '%.12f vs %.12f mmHg' % (core, inline))
# and the concentration form is what physics_sim actually calls
import inspect
src = inspect.getsource(PS.PhysicsSim.step)
check('osmotic_pressure_diff_conc' in src, 'physics_sim.step calls the shared function',
      'no second hand-inlined copy of the osmotic law can drift out of sync')
check('REFLECTION[s] * R_T * (lc[i][s] - sc[s])' not in src,
      'the old hand-inlined van\'t Hoff sum is gone from physics_sim')

# ---------------------------------------------------------------- 5. conservation
print('\n5. GLOBULIN IS CONSERVED LIKE EVERY OTHER SOLUTE')
sim2 = PS.PhysicsSim(weight=70.0, enable={'external': False, 'apparatus': False, 'cori': False})
pools = [sim2.serum, sim2.tissue] + list(sim2.layers)
g0 = sum(p.solutes['Globulin'] for p in pools)
for _ in range(240):
    sim2.step(1.0)
g1 = sum(p.solutes['Globulin'] for p in pools)
check(abs(g1 - g0) < 1e-9 * max(1.0, abs(g0)), 'closed system conserves globulin to FP precision',
      'start %.9f -> end %.9f mmol (drift %.3e)' % (g0, g1, g1 - g0))
check(g0 > 0.0, 'and the pool is non-empty, so that is not a vacuous check', '%.4f mmol' % g0)

# ---------------------------------------------------------------- 6. blindness preserved
print('\n6. THE DEVICE STILL CANNOT SEE PROTEIN')
kappa_with = PC.conductivity(sim.serum, sim.serum_vd['Sodium'])
save = sim.serum.solutes['Globulin']
sim.serum.solutes['Globulin'] = save * 3.0
kappa_without = PC.conductivity(sim.serum, sim.serum_vd['Sodium'])
sim.serum.solutes['Globulin'] = save
check(abs(kappa_with - kappa_without) < 1e-12, 'tripling globulin does not move conductivity',
      'the one fluid property the device reads is untouched by protein')

print('\n' + '=' * 100)
if FAIL:
    print('FAILED %d: %s' % (len(FAIL), '; '.join(FAIL)))
    sys.exit(1)
print('ALL PASS — plasma colloid osmotic pressure is Landis-Pappenheimer on albumin + globulin.')