#!/usr/bin/env python3
"""Conservation-by-construction physics substrate for the PD20 rebuild.

Every physical process moves mass ONLY through `transfer_fluid` (water + carried
solute) or `transfer_solute` (diffusive, solute only). Both are atomic paired
updates clamped to the source pool, so total water and total per-solute mmol are
invariant under any sequence of transfers — mass cannot be created or destroyed.
Each flux is derived from one constitutive law (citations below); none is tuned
to any expected outcome. See PHYSICS_REBUILD.md.

Run:  python3 physics_core.py   (executes the self-tests)
"""

import math                                # film-theory polarisation only (cp_beta)

SOLUTES = ['Sodium', 'Chloride', 'Glucose', 'Potassium', 'Urea',
           'Albumin', 'Globulin', 'Icodextrin', 'Bicarbonate', 'Lactate', 'Calcium', 'Magnesium', 'Phosphate', 'Sulfate',
           'Creatinine']
# GLOBULIN added 2026-08-09. Plasma protein is albumin PLUS globulin, and globulin is ~40% of the
# mass (normal 2.0-3.5 g/dL against albumin 3.5-5.0). The engine modelled albumin alone, which left
# the colloid osmotic pressure — the force that OPPOSES ultrafiltration — structurally too small.
# Globulins are also genuinely lost to the peritoneal effluent, so they cannot be represented as a
# fixed multiple of albumin: PD strips the smaller albumin preferentially, so the albumin FRACTION
# falls over time and a fixed ratio would freeze exactly the quantity that ought to move.
# Calcium and Magnesium added 2026-08-02. Every commercial PD fluid carries Ca 1.25-2.0 and
# Mg 0.25-0.75 mmol/L (Hashimoto & Kamijo, Life 2025, PMC11856993, Table 1); this engine had
# NEITHER, so a device circulating ~29-86 L/day of Ca-free Mg-free dialysate would have stripped
# both severely — a patient-safety gap invisible to the simulator because the species did not exist.

# Van't Hoff: π = R·T·C.  R=62.3637 L·mmHg/(mol·K), T=310.15 K, C in mmol/L→/1000.
R_T = 62.3637 * 310.15 / 1000.0  # SOURCE: textbook-standard (ideal gas constant, 62.3637 L·mmHg/(mol·K))          # = 19.34 mmHg per (mmol/L)

# Ionic VALENCE, needed for the Gibbs-Donnan equilibrium across the peritoneal membrane.
# Albumin is the dominant IMPERMEANT anion of plasma (sieving 0.01 here — it barely crosses),
# so it holds cations in and excludes diffusible anions: at equilibrium the serum/dialysate ratio
# is r for anions and 1/r for cations, not 1. Without this the model let chloride equilibrate
# 1:1 with the dialysate, which is what put serum Cl and HCO3 ~11 mEq/L above physiologic and
# collapsed the anion gap to ~2 (real uraemic 12-18). Added 2026-08-06.
# Albumin charge from Figge J, Mydosh T, Fencl V, J Lab Clin Med 1992;120:713 — the standard
# quantitative model: charge (mEq/L) = [alb g/L] x (0.123*pH - 0.631), i.e. 0.279 mEq per g/L at
# pH 7.4, which at MW 66.5 kDa is 18.6 negative charges per mmol.
VALENCE = {'Sodium': +1, 'Potassium': +1, 'Calcium': +2, 'Magnesium': +2,
           'Chloride': -1, 'Bicarbonate': -1, 'Lactate': -1,
           'Phosphate': -1.8,          # ~80% HPO4(2-) / 20% H2PO4(-) at pH 7.4
           'Sulfate': -2, 'Albumin': -18.6,
           # GLOBULIN carries no NET charge at pH 7.4 and is deliberately given valence 0. The
           # globulin fraction spans isoelectric points on both sides of physiologic pH — gamma
           # globulins (pI 6.5-8.5) are near-neutral to slightly cationic while alpha/beta globulins
           # are slightly anionic — so the fractions cancel. This is the standard treatment: the
           # Figge-Fencl quantitative acid-base model represents plasma protein charge by ALBUMIN
           # and PHOSPHATE only, and does not assign globulin a charge at all. Consequence in this
           # engine: globulin adds colloid osmotic pressure (it has mass) but contributes nothing to
           # the Gibbs-Donnan ratio, nothing to the anion gap, and nothing to conductivity.
            'Globulin': 0,
            'Glucose': 0, 'Urea': 0, 'Icodextrin': 0, 'Creatinine': 0}

# Peritoneal OSMOTIC reflection coefficients σ (three-pore; Rippe B, Venturoli D, Simonsen O, de Arteaga J, Perit Dial Int 2004;24(1):10-27, PMID 15104333).
# Small ions ~0.05, urea 0.02, glucose 0.03, albumin 0.84, icodextrin 0.5.
# Lactate is a small anion (MW 89) -> 0.05 (first-principles; the legacy code's
# 0.5 default for lactate was a bug — PHYSICS_REVIEW.md).
# GLUCOSE 0.03 -> 0.049 on 2026-08-09, from MEASUREMENT rather than from a quoted sigma. Oberg CM
# & Rippe B, Kidney Int Rep 2017;2(6):1128-38 (PMC5733752, Table 1), human PD patients of average
# transport type, report the OSMOTIC CONDUCTANCE TO GLUCOSE directly: LpS*sigma_g = 3.6 uL/min/mmHg,
# alongside LpS = 0.074 mL/min/mmHg. That fixes sigma_g = 3.6/74 = 0.049. The engine was producing
# 0.080*0.03 = 2.4 uL/min/mmHg, i.e. 33% BELOW the measured osmotic drive — and this is the single
# most ultrafiltration-relevant coefficient in the model, because it sets how hard dextrose pulls.
# GLOBULIN 0.90: reflected MORE than albumin because it is bigger. The globulin fraction averages
# ~150 kDa (IgG) against albumin's 66.5 kDa, and in the three-pore model sigma rises steeply with
# molecular radius through the large-pore term. 0.90 sits between albumin's measured 0.84 and the
# 1.0 of a solute the membrane excludes completely.
REFLECTION = {'Sodium': 0.05, 'Chloride': 0.05, 'Glucose': 0.049, 'Potassium': 0.05,
              'Urea': 0.02, 'Albumin': 0.84, 'Globulin': 0.90, 'Icodextrin': 0.5, 'Bicarbonate': 0.05,
              'Lactate': 0.05, 'Calcium': 0.05, 'Magnesium': 0.05, 'Phosphate': 0.05, 'Sulfate': 0.05,
              'Creatinine': 0.02}   # small neutral solute, like urea (113 Da)   # small divalent ions

# Convective SIEVING S = fraction of solute swept through pores with UF water (solvent-drag
# sieving coefficient, three-pore). NOTE this is the CONVECTIVE coefficient and is INDEPENDENT
# of the OSMOTIC reflection σ above — they are two different three-pore quantities (osmotic
# water flow is aquaporin/small-pore dominated; solvent drag is small/large-pore dominated), so
# S ≠ 1−σ in general. Values (Nolph; Rippe) reproduce published Na removal & UF. Macromolecules ~0.
#
# SODIUM 0.72 -> 0.55 on 2026-08-09, ON THE INVENTOR'S INSTRUCTION to make the simulation as
# realistic and literature-supported as possible. This is the single most load-bearing constant in
# the engine and it had been sitting at the FAVOURABLE EDGE of its range.
#   The apparent sodium sieving coefficient is set by how much ultrafiltrate crosses by free water
#   (aquaporin, carrying NO sodium) versus through small pores (carrying sodium at nearly plasma
#   concentration), so S_Na is approximately 1 - fFWT. Two independent sources:
#     Bernardo 2012 measured free-water transport at 0.45 +/- 0.16 of ultrafiltration -> S_Na
#       0.39-0.71, centre 0.55.
#     Helman J, Wahlgren H, Andersson L, Morelle J, Oberg CM, 2024 (PMC10914194) state a fixed
#       40-50% of ultrafiltration is free water and adopt 50% -> S_Na 0.50-0.60.
#   The two agree on ~0.55. The old 0.72 implied only ~28% free water, BELOW every measured value,
#   and it flattered sodium removal while suppressing ultrafiltration.
#   MEASURED CONSEQUENCE (P4 anuric, 30 d, sieving_fwt.log): cumulative UF 40.3 L at 0.72 -> 57.0 L
#   at 0.55, with glucose absorbed 184 -> 205 g/day. Every published ultrafiltration figure moves.
# THE DISCOUNT IS APPLIED TO EVERY SMALL SOLUTE, NOT TO SODIUM ALONE. Free water crossing AQP1
# carries NO solute of any kind, so the (1 - fFWT) reduction is a property of the water pathway and
# therefore of the whole small-solute class. Rippe 2001 (PMID 11687943) states it exactly that way:
# "nearly 40% of the total osmotic water flow occurs through molecular water channels, termed
# aquaporin-1. This causes an inequality between 1 - sigma and the sieving coefficient FOR SMALL
# SOLUTES" — small solutes, not sodium.
# Discounting sodium alone was a real error, caught by adversarial review within the hour: it made
# the sieved ultrafiltrate carry a large net anion excess, stripping many more mEq of anion than of
# cation from the patient per litre.
# THE CHARGE NUMBERS BELOW WERE RE-DERIVED FROM THESE COEFFICIENTS ON 2026-08-09 and are no longer
# quoted from memory. A second adversarial pass challenged them, computed -22.52 mEq/L, and was
# WRONG — but it was right that they had never been re-checked after the coefficients moved, so the
# published figures were unverified. Evaluate sum(SIEVING[s]*VALENCE[s]*C[s]) against the engine's
# own t0 serum (physics_sim serum_c: Na 140, Cl 103, K 4, HCO3 24, lactate 1, Ca 1.20, Mg 0.85,
# phosphate 1.60, sulfate 2.1, albumin 0.70) and the answer is:
#     sodium discounted ALONE (the bug)      -19.2 mEq/L
#     uniform 0.764 discount (this dict)      +3.52 mEq/L
# Non-selective bulk removal at S=1 throughout carries 0.00 mEq/L — the starting serum is now
# ELECTRONEUTRAL (sulfate was 5.0 mmol/L = 10 mEq/L, double-counting the anion gap against albumin's
# full Figge charge; corrected to 2.1 mmol/L on 2026-08-22, S19/P5). The earlier -5.80 mEq/L figure
# was that pre-existing cation deficit, now gone.
# Measured against that reference, sieving distorts the removed fluid by +3.52 mEq/L relative to
# plasma. Both forms are given because the raw sum alone invites exactly the confusion above.
# Every value below is the previous one times 0.55/0.72 = 0.764, which preserves the measured
# RELATIVE restriction between species while setting the class reference to the free-water centre.
# Macromolecules are untouched: albumin and icodextrin barely cross at all, so the aquaporin
# argument does not apply to them.
SIEVING = {'Sodium': 0.550, 'Chloride': 0.588, 'Glucose': 0.382, 'Potassium': 0.535,
           'Urea': 0.688, 'Albumin': 0.01, 'Globulin': 0.005, 'Icodextrin': 0.02, 'Bicarbonate': 0.535,
           'Lactate': 0.611, 'Calcium': 0.458, 'Magnesium': 0.458, 'Phosphate': 0.497,
           'Sulfate': 0.458,
           'Creatinine': 0.65}   # 113 Da, between urea (0.688) and glucose (0.382)
# GLOBULIN 0.005, half of albumin's: convective transport through the large pores falls steeply with
# size, and clinically the effluent protein of a PD patient IS albumin-dominated (albumin is ~50-60%
# of a total protein loss of 5-15 g/day even though it is only ~60% of plasma protein by mass, i.e.
# it crosses preferentially). This is what makes the albumin FRACTION of plasma protein fall in a
# long-term PD patient, and it is why globulin is tracked separately rather than scaled off albumin.

# TEMPERATURE AND ACTIVITY, STATED PROPERLY (corrected 2026-08-22, literature audit finding #9).
# These are LIMITING (infinite-dilution) molar conductivities Λ°, tabulated at 25 °C (CRC). Two
# things were wrong with the old one-line label "~25 °C":
#   1. NO IONIC-STRENGTH CORRECTION. Real solutions are not at infinite dilution. By Kohlrausch,
#      at the physiologic I ≈ 0.15 M the actual Λm is ≈ 0.85 Λ°, so summing Λ° overstates a true
#      25 °C reading by ≈ 18 %.
#   2. Conductivity rises ≈ 2 %/°C, so a body-temperature (37 °C) fluid reads ≈ 24 % above its
#      25 °C value.
# Those two errors run in OPPOSITE directions and very nearly cancel: the number this function
# returns lands within ~5 % of an UNCOMPENSATED 37 °C probe reading. So the VALUE is defensible for
# a sensor sitting in body-temperature dialysate and reporting without temperature compensation —
# but the stated temperature was the wrong one, and that matters because THIS IS A LIVE CONTROLLER
# INPUT: the device's three conductivity sensors are its only chemical sensing. Any hardware
# specification, any calibration, and any published methods section must state which convention the
# probe uses. Computed with these values: dialysate ≈ 15.7 mS/cm, serum ≈ 17.0 mS/cm.
# CONVENTION OF THIS MODEL: an uncompensated probe in 37 °C fluid. Do not "fix" the 18 % activity
# gap without also removing the 24 % temperature gap, or the result will be wrong by both.
#
# Limiting molar ionic conductivity Λ° (S·cm²/mol, 25 °C table values; CRC). Only IONS carry current —
# glucose/urea/albumin/icodextrin are neutral and contribute nothing. Conductivity is the
# ONE fluid property the device measures on the dialysate/effluent (a bulk electrical
# reading, NOT any patient solute concentration): κ = Σ Λ°_i · c_i. It is dominated by
# Na+/Cl- (the major ions), so it tracks the sodium/osmotic state well and potassium only
# weakly — the device senses the ionic bath, never a specific serum concentration.
LAMBDA_ION = {'Sodium': 50.1, 'Potassium': 73.5, 'Chloride': 76.3,
              'Bicarbonate': 44.5, 'Lactate': 38.8,   # lactate was 40.0; CRC gives 38.8
              # Divalent cations: CRC limiting molar conductivity is quoted per EQUIVALENT
              # (1/2 Ca2+ = 59.5, 1/2 Mg2+ = 53.1), so PER MOLE they carry twice that.
              'Calcium': 119.0, 'Magnesium': 106.2,
              # Phosphate at pH 7.4 is ~80% HPO4(2-) / 20% H2PO4(-), mean valence ~ -1.8.
              # Lambda per mole scaled from the equivalent conductivity accordingly.
              'Phosphate': 98.4,  # recomputed from the code's OWN stated 80/20 mixture at
              #               pH 7.4: 0.8*114 (HPO4 2-) + 0.2*36 (H2PO4 -) = 98.4,
              #               not 103.0 (audit finding, +4.7% as written),
              # SULFATE stands for the whole retained uraemic anion pool — sulfate from sulfur
              # amino acids plus the organic acids (hippurate, indoxyl sulfate, p-cresyl sulfate)
              # that the failing kidney cannot excrete. It is what makes ESRD a HIGH-anion-gap
              # acidosis. Divalent: lambda per mole is 2x the per-equivalent 80.0.
              'Sulfate': 160.0}


def conductivity(pool, eff_volume=None):
    """Bulk conductivity of a pool's fluid (mS/cm-scale proxy) = Σ Λ°_i·[ion]_i. Accepts a
    Pool (uses its concentrations) or a plain {solute: conc} dict. Neutral solutes ignored."""
    if hasattr(pool, 'conc'):
        if isinstance(eff_volume, dict):
            c = {s: pool.conc(s, eff_volume.get(s)) for s in LAMBDA_ION}
        else:
            c = {s: pool.conc(s, eff_volume) for s in LAMBDA_ION}
    else:
        c = pool
    return sum(LAMBDA_ION[s] * c.get(s, 0.0) for s in LAMBDA_ION) / 1000.0


class Pool:
    """A fluid compartment: water volume (L) + solute amounts (mmol)."""
    def __init__(self, volume=0.0, solutes=None):
        self.volume = float(volume)
        self.solutes = {s: 0.0 for s in SOLUTES}
        if solutes:
            for s, v in solutes.items():
                self.solutes[s] = float(v)

    def conc(self, s, eff_volume=None):
        v = self.volume if eff_volume is None else eff_volume
        return self.solutes[s] / v if v > 1e-12 else 0.0


def transfer_fluid(src, dst, volume, carried_conc=None):
    """Move `volume` L water src->dst. Carried solute per L = carried_conc[s] if
    given, else src's own concentration (mmol / src.volume). Clamped to src.volume
    and to each solute amount. Atomic paired update => conserves water+solute.
    Returns volume actually moved. (carried_conc lets UF carry SIEVED solute and
    serum use per-solute Vd concentrations; default = bulk carry, used by lymph.)"""
    if volume <= 0.0:
        return 0.0
    volume = min(volume, src.volume)
    for s in SOLUTES:
        if carried_conc is None:
            c = src.solutes[s] / src.volume if src.volume > 1e-12 else 0.0
        else:
            c = carried_conc.get(s, 0.0)
        m = min(c * volume, src.solutes[s])      # never move more solute than exists
        src.solutes[s] -= m
        dst.solutes[s] += m
    src.volume -= volume
    dst.volume += volume
    return volume


def transfer_solute(src, dst, s, mmol):
    """Diffusive solute transfer (no water), signed: + = src->dst. Clamped so
    neither pool goes negative. Atomic paired update => conserves solute s."""
    if mmol >= 0.0:
        mmol = min(mmol, src.solutes[s])
    else:
        mmol = -min(-mmol, dst.solutes[s])
    src.solutes[s] -= mmol
    dst.solutes[s] += mmol
    return mmol


def total_water(pools):
    return sum(p.volume for p in pools)


def total_solute(pools, s):
    return sum(p.solutes[s] for p in pools)


# ---------------------------------------------------------------------------
# COLLOID osmotic pressure of plasma protein — LANDIS-PAPPENHEIMER, not van't Hoff.
#
# Van't Hoff (π = R·T·C) is the ideal-solution law and it is simply WRONG for plasma protein at
# physiologic concentration: protein solutions are strongly non-ideal (excluded volume, protein-
# protein interaction, and the Donnan contribution of the protein's own charge), so the real
# osmotic pressure rises much faster than linearly with concentration. The standard empirical
# replacement, in every physiology text since:
#
#     COP (mmHg) = 2.1*C + 0.16*C^2 + 0.009*C^3          C = TOTAL plasma protein, g/dL
#
#     Landis EM & Pappenheimer JR, "Exchange of substances through the capillary walls",
#     Handbook of Physiology, Section 2: Circulation, Vol II, American Physiological Society,
#     1963, p.975.
#
# At a normal total protein of 7.0 g/dL: 14.7 + 7.84 + 3.09 = 25.6 mmHg, which is the textbook
# plasma colloid osmotic pressure. The ideal van't Hoff value on the same protein mass is ~15,
# so the NON-IDEAL terms supply about 40% of the real pressure. That missing 40%, together with
# the absent globulin fraction, is the entire reason this engine produced a transperitoneal
# protein oncotic gradient of 11.2 mmHg where human PD patients measure 22
# (Öberg CM & Rippe B, Kidney Int Rep 2017;2(6):1128-38, PMC5733752, Table 1, "Transperitoneal
# oncotic pressure gradient (Δπ_prot) 22 mm Hg").
#
# WHY THIS MATTERS PHYSICALLY: the protein oncotic gradient is the force holding water INSIDE the
# circulation, i.e. it OPPOSES ultrafiltration. At half its true strength the engine ultrafiltrated
# too easily, and that error had been partly hidden by two others pulling the other way (a glucose
# reflection coefficient 33% too low and an intraperitoneal pressure law 40% too low). Those two
# were corrected earlier the same day; this is the third of the three, and they had to be fixed
# together or not at all, because the model only looked plausible through their cancellation.
#
# Molecular weights are used ONLY to convert the engine's mmol/L pools into the g/dL the law wants.
PROTEIN_MW_G_PER_MMOL = {'Albumin': 66.5,      # human serum albumin, 66.5 kDa
                         'Globulin': 150.0}    # mass-weighted mean of the globulin fractions,
#                                                dominated by IgG (150 kDa); alpha/beta/gamma span
#                                                ~50 kDa to ~720 kDa (alpha-2-macroglobulin).
PROTEINS = ('Albumin', 'Globulin')


# The Landis-Pappenheimer fit is an EMPIRICAL curve over the range plasma actually occupies. Above
# ~10 g/dL total protein it is extrapolation, and because the dominant term out there is CUBIC, the
# extrapolation diverges fast and confidently: at 20 g/dL the cubic asserts 178 mmHg on no evidence
# whatever. That state is reachable in this engine — not in a normal patient, but by HAEMO-
# CONCENTRATION in a dying one. Adversarial review measured total protein reaching 11.97 g/dL in a
# 30-day arm whose patient was dehydrating to the fatal -7 L trip: serum volume had collapsed from
# 17.1 L to 8.2 L, so the protein AMOUNT was normal and only the volume was gone.
# Beyond the validated ceiling the curve therefore continues along its own TANGENT at that point:
# C1-continuous (no step, no kink in the force), monotone, and honest about being a straight-line
# extension rather than a measured law. The two agree to within 3% just past the boundary (62.0 vs
# 63.8 mmHg at 12 g/dL), so this changes nothing in any physiologic state — it only stops an
# unvalidated cubic from driving the Starling balance in an extreme one.
# Occurrences are COUNTED, never silently swallowed: if this fires during a run that is a fact about
# the run, and test_oncotic_landis asserts the counter stays at zero in normal operation.
COP_VALID_MAX_G_DL = 10.0  # SOURCE: INFERRED/UNVERIFIED (extrapolation guard)
COP_EXTRAPOLATIONS = [0]          # list so it is mutable from module scope; index 0 = count


def _cop_cubic(c):
    return 2.1 * c + 0.16 * c * c + 0.009 * c * c * c


def _cop_slope(c):
    return 2.1 + 0.32 * c + 0.027 * c * c


def colloid_osmotic_pressure_g_dl(c_total_g_dl):
    """Landis-Pappenheimer colloid osmotic pressure (mmHg) of a plasma-protein solution at
    total protein concentration `c_total_g_dl` (g/dL). The published cubic over its validated
    0-10 g/dL range; its own tangent beyond that (see above)."""
    c = max(0.0, c_total_g_dl)
    if c <= COP_VALID_MAX_G_DL:
        return _cop_cubic(c)
    COP_EXTRAPOLATIONS[0] += 1
    return (_cop_cubic(COP_VALID_MAX_G_DL)
            + _cop_slope(COP_VALID_MAX_G_DL) * (c - COP_VALID_MAX_G_DL))


def protein_g_dl_conc(conc):
    """Total protein (g/dL) from a {solute: mmol/L} mapping."""
    return sum(conc.get(p, 0.0) * PROTEIN_MW_G_PER_MMOL[p] for p in PROTEINS) / 10.0


def protein_cop_conc(conc):
    """Colloid osmotic pressure (mmHg) of a {solute: mmol/L} mapping."""
    return colloid_osmotic_pressure_g_dl(protein_g_dl_conc(conc))


def protein_sigma_conc(conc):
    """Effective osmotic reflection coefficient of the protein MIXTURE, mass-weighted between
    albumin and globulin. Not a constant: as peritoneal dialysis strips the smaller albumin
    preferentially, the surviving protein is globulin-richer and therefore better reflected, so
    sigma drifts up over a long treatment. Falls back to albumin's value on an empty solution."""
    num = den = 0.0
    for p in PROTEINS:
        g = conc.get(p, 0.0) * PROTEIN_MW_G_PER_MMOL[p]
        num += REFLECTION[p] * g
        den += g
    return num / den if den > 1e-12 else REFLECTION['Albumin']


def _pool_conc(pool, vd):
    return {s: (pool.conc(s) if vd is None else pool.conc(s, vd[s])) for s in SOLUTES}


def protein_cop(pool, vd=None):
    return protein_cop_conc(_pool_conc(pool, vd))


def protein_sigma(pool, vd=None):
    return protein_sigma_conc(_pool_conc(pool, vd))


def protein_oncotic_diff(layer, serum, serum_vd):
    """σ_prot·(COP_layer − COP_serum), mmHg, signed like the rest of Δπ (layer minus serum).
    Negative in every realistic case, because the dialysate carries essentially no protein."""
    lc = _pool_conc(layer, None)
    sc = _pool_conc(serum, serum_vd)
    return _protein_term(lc, sc)


def _protein_term(lc, sc):
    """The protein contribution to Δπ from two concentration mappings. σ is evaluated on the
    SERUM side, which is the side that has protein — evaluating it on the near-protein-free
    dialysate would make the coefficient meaningless (0/0) and its fallback would then decide
    the answer."""
    return protein_sigma_conc(sc) * (protein_cop_conc(lc) - protein_cop_conc(sc))


def osmotic_pressure_diff_conc(lc, sc):
    """Δπ (mmHg) across the peritoneal membrane, layer minus serum, from two {solute: mmol/L}
    mappings. Positive => layer hypertonic => UF water moves serum->layer.

    CRYSTALLOIDS follow van't Hoff, Σ σ_s·R·T·(C_layer − C_serum), which is correct for small
    solutes at these concentrations. PROTEIN does not: it is evaluated by Landis-Pappenheimer
    on each side and differenced (see above). This is the single shared implementation — the
    per-layer loop in physics_sim.step calls it rather than keeping its own copy, so the two
    cannot drift apart, which they previously could."""
    dpi = 0.0
    for s in SOLUTES:
        if s in PROTEINS:
            continue
        dpi += REFLECTION[s] * R_T * (lc[s] - sc[s])
    return dpi + _protein_term(lc, sc)


def osmotic_pressure_diff(layer, serum, serum_vd):
    """Δπ (mmHg) across the peritoneal membrane, layer minus serum. Pool-based wrapper around
    osmotic_pressure_diff_conc; the layer uses its actual volume, serum its per-solute Vd."""
    return osmotic_pressure_diff_conc(_pool_conc(layer, None), _pool_conc(serum, serum_vd))


# ---------------------------------------------------------------------------
# Flux 4 — lymphatic absorption  (layer -> serum), carries solute @ layer conc.
#   lymph_rate·(V_i/V_perit)·dt   [Flessner, Am J Physiol 1991;261:H1040]
# ---------------------------------------------------------------------------
def lymph_flux(layer, serum, lymph_rate, vol_frac, dt):
    return transfer_fluid(layer, serum, lymph_rate * vol_frac * dt)


# ---------------------------------------------------------------------------
# Flux 3+6 — peritoneal ultrafiltration + convective solute drag (serum<->layer)
#   UF_i = Kf_i·(ΔP_hydro_i + Δπ_i)·dt   [Rippe 3-pore; Flessner]
#   convective solute carried at SIEVING[s]·C_serum (forward) / layer conc (back).
# ---------------------------------------------------------------------------
def uf_flux(serum, layer, Kf_layer, dP_hydro, dt, serum_vd):
    dpi = osmotic_pressure_diff(layer, serum, serum_vd)
    uf = Kf_layer * (dP_hydro + dpi) * dt
    if uf >= 0.0:   # serum -> layer, sieved
        carried = {s: SIEVING[s] * serum.conc(s, serum_vd[s]) for s in SOLUTES}
        return transfer_fluid(serum, layer, uf, carried_conc=carried)
    else:           # back-filtration layer -> serum, carries at layer conc
        return -transfer_fluid(layer, serum, -uf)


# ---------------------------------------------------------------------------
# Flux 5 — diffusive membrane transport (serum<->layer), solute only.
#   J_s = MTAC_s·(C_serum − C_layer)·dt    [Smit 2003; Rippe]
# ---------------------------------------------------------------------------
def diffusion_flux(serum, layer, mtac, dt, serum_vd):
    for s in SOLUTES:
        j = mtac[s] * (serum.conc(s, serum_vd[s]) - layer.conc(s)) * dt
        transfer_solute(serum, layer, s, j)


# ---------------------------------------------------------------------------
# Flux 7 — serum<->tissue water (Starling), crystalloid osmotic drive.
#   Jv = Kf_t·Δπ_crys·dt   [Starling 1896; Levick & Michel 2010]
# ---------------------------------------------------------------------------
def starling_flux(serum, tissue, Kf_t, dt, serum_vd, tissue_vd):
    """SUPERSEDED — NOT LIVE PHYSICS. Kept only as a reference implementation, exactly like
    fickian_flux above, and pinned as DIVERGENT by test_engine_equivalence2.check_starling_divergent
    so this dead reference is never mistaken for the model.

    THE BARRIER IS THE CELL MEMBRANE, NOT A CAPILLARY. `tissue` is the ICF compartment
    (physics_sim: tissue0 = 2/3 of total body water, seeded with 140 mM potassium), so serum<->tissue
    is ECF<->ICF. A flat sigma = 0.05 for every solute is wrong physiology there: sodium and glucose
    are EFFECTIVE osmoles across a cell membrane (sigma ~ 1, they drive tonicity) and urea is
    INEFFECTIVE (sigma ~ 0, it equilibrates and shifts no water). The 0.05 is the PERITONEAL
    small-ion value from REFLECTION, borrowed for a membrane it does not describe; physics_sim
    line 191 records that borrowing as the "H3 fix" and it was corrected there, not here.

    WHAT THE LIVE ENGINE ACTUALLY DOES (physics_sim step, "Transcellular water PARTITION"): it
    places serum and ICF at ALGEBRAIC iso-tonic equilibrium every step on EFFECTIVE osmoles only —
    A = 2*Na + glucose against the tracked ICF osmoles — moving pure water that carries no solute
    and preserving serum+tissue volume. Urea is absent from A by design, which is the correct
    treatment and is why dialysis does not shift cell water on urea alone. Being algebraic rather
    than rate-based, it cannot oscillate or drain the ICF.

    HISTORY OF THIS DOCSTRING, since it is a lesson about reviewing: on 2026-08-16 I annotated this
    function as crossing "the systemic capillary wall" and flagged its 0.05 as a live uncited
    coefficient needing an inventor decision. BOTH WERE WRONG. It is the cell membrane, the function
    is dead, and the engine had already found and fixed the exact error I thought I had found. I had
    read physics_core in isolation without tracing whether anything calls it — the same
    incomplete-trace mistake that produced the earlier RO-pressure error, in a different costume.
    """
    dpi = 0.0
    for s in SOLUTES:
        if s in ('Albumin', 'Globulin', 'Icodextrin'):
            continue
        dpi += 0.05 * R_T * (tissue.conc(s, tissue_vd[s]) - serum.conc(s, serum_vd[s]))
    jv = Kf_t * dpi * dt
    if jv >= 0.0:   # serum -> tissue
        carried = {s: serum.conc(s, serum_vd[s]) for s in SOLUTES}
        return transfer_fluid(serum, tissue, jv, carried_conc=carried)
    else:
        carried = {s: tissue.conc(s, tissue_vd[s]) for s in SOLUTES}
        return -transfer_fluid(tissue, serum, -jv, carried_conc=carried)


# ---------------------------------------------------------------------------
# Flux 8 — inter-layer Fickian mixing (layer_i <-> layer_j), solute only.
#   J = D_eff·A·ΔC/Δz·dt  [Fick]. NOT capped at a fraction of the pool — the only bound is
#   transfer_solute's clamp to the source amount, which prevents a negative pool but does not
#   bound the step for stability. (The comment here used to claim a 0.5·pool cap that the code
#   has never implemented.) This law is SUPERSEDED: the live engine models the cavity as
#   well-mixed and does not call it — test_engine_equivalence2.check_fickian_superseded asserts
#   exactly that, so it must not be described as live physics anywhere.
# ---------------------------------------------------------------------------
def fickian_flux(layer_lo, layer_hi, D_eff_cm2s, area_cm2, dz_cm, dt_min):
    dt_s = dt_min * 60.0
    for s in SOLUTES:
        # ΔC in mmol/L = mmol/(1000 cm^3); flux mmol = D[cm2/s]*A[cm2]*ΔC/dz[cm]*t[s]
        c_lo = layer_lo.conc(s); c_hi = layer_hi.conc(s)
        j = D_eff_cm2s * area_cm2 * ((c_hi - c_lo) / 1000.0) / dz_cm * dt_s  # hi->lo if positive
        transfer_solute(layer_hi, layer_lo, s, j)


# ---------------------------------------------------------------------------
# Flux 9+10 — RO dialysate regeneration. The drained spent dialysate ("feed")
# is split into PERMEATE (recovered clean water, re-infused) + REJECT (waste,
# leaves the system). Solution-diffusion: R_s = 1−B_s/(B_s+Jw), Jw = A·(ΔP−Δπ);
# flux-feasible cap zeroes permeate when Δπ ≥ ΔP (v62s). Mass-conserving by
# construction: feed is fully partitioned into permeate + reject (feed ends empty).
# [Wijmans & Baker, J Membr Sci 1995;107:1]
# ---------------------------------------------------------------------------
# Osmotic weighting used when estimating the RO feed's osmotic pressure. Small solutes count 1
# osmole per mole; the MACROMOLECULES are discounted by their molecular weight in kDa. Values for
# icodextrin and albumin are UNCHANGED from the original inline expression this constant replaced,
# so no shipped number moves; globulin follows the same convention at 150 kDa. FLAGGED, not fixed:
# per van't Hoff a mole of macromolecule contributes one osmole regardless of its mass, so this
# discount is dimensionally odd. It is numerically inert either way — the RO feed is spent
# DIALYSATE, whose protein content is ~2e-3 mmol/L, i.e. four orders below the salts — so changing
# it would move nothing and would only break bit-comparability with every published run.
# ICODEXTRIN USES THE NUMBER-AVERAGE MW (fixed 2026-08-22, literature audit finding #8).
# Was 1/16.2, from icodextrin's WEIGHT-average MW of 16.2 kDa. Colligative properties — osmotic
# pressure among them — depend on the NUMBER of particles, so the correct basis is the
# NUMBER-average Mn, which for clinical icodextrin is ~5-6.5 kDa. Using Mw undercounted icodextrin
# osmoles by roughly 2.5-3x. Icodextrin is inert in the shipped configuration (it starts at 0
# everywhere) so nothing published moves, but it is the SECOND osmotic agent in this model and any
# icodextrin arm would have been wrong by that factor — which is exactly the kind of latent error
# that surfaces only once someone runs the experiment. 5.75 kDa = midpoint of the published range.
# The albumin/globulin entries remain dimensionally odd (van't Hoff counts moles, not kilodaltons)
# and are numerically inert here: the RO feed is spent DIALYSATE, protein ~2e-3 mmol/L, four orders
# below the salts. Left as-is deliberately rather than silently changing a term nothing exercises.
MACRO_OSM_FACTOR = {'Icodextrin': 1 / 5.75, 'Albumin': 1 / 66.5, 'Globulin': 1 / 150.0}


# RO PRESSURE ACTUATOR ENVELOPE (bar). Lives HERE, in the shared base, because from 2026-08-22 the
# controller may COMMAND pressure (inventor: "controller also needs to control the pressure") and
# both the plant (physics_sim, which enforces the envelope) and the policy (physics_controller,
# which must not ask for what cannot be built) need the same two numbers. Defining it in either one
# would force the other to import it and create a cycle, or — worse — a second copy that drifts.
# CEILING: the highest REAL element class, Dow XUS180808 ultra-high-pressure at 120 bar; the
# seawater SW30 class tops out at 82.7. The network's TRAINING span deliberately runs to 200 bar
# (_RNG_SPAN['pressure']) so no shipped point sits at the edge of what it has seen — but a training
# span is not a hardware rating, and the PLANT is never commanded past what can be built.
# FLOOR: well under the 120 bar shipped point, so the loop has somewhere to go when water is slack.
RO_PRESSURE_MIN_BAR, RO_PRESSURE_MAX_BAR = 40.0, 120.0


def ro_nominal_perm(ro_wp=5.0, nominal_rejection=None, ro_pressure=15.0):  # SOURCE: design choice (reference operating point)
    nr = nominal_rejection or {'Sodium': 0.995, 'Chloride': 0.995, 'Glucose': 0.99,
        # Monovalent salt rejection raised to the COMMERCIALLY-AVAILABLE BEST (2026-08-24, inventor:
        # "don't assume cheap membrane, do best there is commercially available"). Brackish-water
        # polyamide TFC elements reject NaCl at 99.0-99.8% at standard conditions (Dow FILMTEC BW30
        # nominal 99.5%; Hydranautics CPA 99.7%). The old Na 0.93 / Cl 0.94 was a deliberate
        # cheap-membrane assumption that was NOT conservative: it recycled ~7% of the feed salt back
        # to the patient. At ~99.5% the permeate is essentially pure water, so the blind mixer's
        # "added = concentration" holds for sodium exactly as it already did for glucose.
        # Glucose 0.20 -> 0.99 (2026-07-17). A real RO membrane rejects sugars ~96-99% (measured
        # polyamide TFC: 96.5% glucose @10 bar; RO is used industrially to CONCENTRATE sugars), and a
        # neutral 180 Da solute is rejected MORE than Na+, not 4x less — the old 0.20 was inverted and
        # unphysical. It also forced the mixer to read permeate.conc() (an oracle the device cannot
        # have): blind dosing at 0.20 lands +47% over target tonicity. At 0.99 permeate ~ pure water,
        # so "added = concentration" holds and the mixer needs no knowledge of current dialysate.
        # Urea stays 0.50: RO genuinely rejects urea poorly (~40-70%) — that is real physics, not a bug.
        'Potassium': 0.99, 'Urea': 0.50, 'Albumin': 0.99, 'Globulin': 0.995, 'Icodextrin': 0.99,
        'Bicarbonate': 0.95, 'Lactate': 0.92,
        # Divalent cations are rejected MORE than monovalent by a polyamide RO membrane (higher
        # charge -> stronger Donnan exclusion; manufacturer data consistently shows Ca/Mg rejection
        # above Na/K, ~99.5%+ commercially). Raised with the monovalents so the ordering is preserved.
        'Calcium': 0.996, 'Magnesium': 0.996, 'Phosphate': 0.99, 'Sulfate': 0.998,
        'Creatinine': 0.50}   # small neutral polar, like urea (rejected poorly)
    # Reference water flux = ro_wp·(applied pressure − reference feed osmotic pressure ~8 bar).
    # Tracks ro_pressure so the flux-feasible cap (ff = jw_raw/ref_jw) stays valid when the applied
    # pressure is retuned (e.g. a power/energy study). Default 15 bar reproduces the prior value 35.
    ref_jw = max(ro_wp * (ro_pressure - 8.0), 0.1)
    return {s: ref_jw * (1.0 - min(0.999, max(0.001, nr.get(s, 0.5))))
            / min(0.999, max(0.001, nr.get(s, 0.5)))
            for s in SOLUTES}, ref_jw


def cp_beta(osm_bulk, ro_wp, ro_pressure, ref_jw, cp_beta_ref):
    """CONCENTRATION POLARISATION modulus β = c_wall/c_bulk, film theory: β = exp(Jw/k)
    [Wijmans & Baker; Baker, Membrane Technology and Applications].

    Solute rejected at the membrane accumulates in the boundary layer, so the wall sees a HIGHER
    concentration than the bulk. Two adverse consequences, both previously absent from this model:
      1. salt passage rises — permeate concentration is set by the WALL, c_perm = β·c_bulk·(1−R);
      2. flux falls — the opposing osmotic pressure is the WALL value, so the driving force is
         (ΔP − β·π_bulk), which in turn lowers Jw and hence β. The two are solved together.

    k is back-solved from `cp_beta_ref`, the modulus at the REFERENCE flux: k = ref_jw/ln(β_ref).
    cp_beta_ref = 1.0 (the DEFAULT) gives k → ∞, β ≡ 1, and reproduces the pre-polarisation engine
    BYTE-IDENTICALLY — so enabling this is opt-in and no published number moves until someone sets it.

    β_ref is NOT VERIFIED against a source: it could not be obtained in the 2026-07-19 session
    (PubMed returned a CFD paper that explicitly avoids film theory; MDPI 403'd). Typical RO moduli
    are quoted ~1.1-1.4. exp_polarization.py sweeps that whole range and the conclusion it draws is
    invariant across it, so nothing currently depends on pinning the value."""
    if cp_beta_ref <= 1.0 + 1e-12:
        return 1.0
    k = ref_jw / math.log(cp_beta_ref)
    jw = ro_wp * (ro_pressure - osm_bulk)          # start from the unpolarised flux
    for _ in range(12):                            # fixed-point: β raises π_wall, which lowers Jw
        b = math.exp(max(0.0, jw) / k)
        jw = ro_wp * (ro_pressure - b * osm_bulk)
    return math.exp(max(0.0, jw) / k)


def ro_split(feed, permeate, reject_waste, recovery, ro_pressure, ro_wp, ro_sp, ref_jw,  # SOURCE: design choice (solution-diffusion split)
             cp_beta_ref=1.0):
    osm = 0.0258 * sum(feed.conc(s) * MACRO_OSM_FACTOR.get(s, 1.0) for s in SOLUTES)
    beta = cp_beta(osm, ro_wp, ro_pressure, ref_jw, cp_beta_ref)   # 1.0 by default => no-op
    jw_raw = ro_wp * (ro_pressure - beta * osm)             # wall osmotic pressure opposes ΔP
    jw = max(jw_raw, 0.01)
    flux_feasible = max(0.0, min(1.0, jw_raw / ref_jw))     # v62s cap
    permeate_vol = recovery * flux_feasible * feed.volume
    carried = {s: beta * feed.conc(s) * (1.0 - max(0.0, min(1.0, 1.0 - ro_sp[s] / (ro_sp[s] + jw))))
               for s in SOLUTES}                            # permeate conc = β·feed·(1−rejection)
    transfer_fluid(feed, permeate, permeate_vol, carried_conc=carried)
    transfer_fluid(feed, reject_waste, feed.volume)         # all remaining feed -> waste
    return permeate_vol


# ---------------------------------------------------------------------------
# Fluxes 11-15 — external exchange. Each is serum<->an external pool (so the
# external pools count in the conservation total). Urine/insensible/respiratory/
# stool leave; intake/GI enter. All paired transfers => conserved.
# ---------------------------------------------------------------------------
def renal_clearance(serum, urine, vol_rate, gfr_frac, dt, serum_vd):
    """Water vol_rate·dt to urine; solute filtered at gfr_frac·serum conc."""
    carried = {s: gfr_frac * serum.conc(s, serum_vd[s]) for s in SOLUTES}
    return transfer_fluid(serum, urine, vol_rate * dt, carried_conc=carried)

def insensible_loss(serum, lost, vol_rate, dt):
    """Pure water loss (skin/respiratory) — carries no solute."""
    return transfer_fluid(serum, lost, vol_rate * dt, carried_conc={s: 0.0 for s in SOLUTES})

def oral_intake(intake_src, serum, vol, conc):
    """Drinking/food water + solute into serum."""
    return transfer_fluid(intake_src, serum, vol, carried_conc=conc)


# ---------------------------------------------------------------------------
# Flux 16 — Cori reaction (in serum): 1 Lactate -> 1 Bicarbonate. NOT a transfer
# (species conversion). Conserves total osmoles + charge (both −1 anions). Tracked
# so the per-species conservation invariant accounts for it. [Cori cycle]
# ---------------------------------------------------------------------------
def cori_react(serum, amount_mmol):
    amount_mmol = max(0.0, min(amount_mmol, serum.solutes['Lactate']))
    serum.solutes['Lactate'] -= amount_mmol
    serum.solutes['Bicarbonate'] += amount_mmol
    return amount_mmol


# ===========================================================================
# SELF-TESTS  (each asserts conservation EVERY step + flux matches its own law)
# ===========================================================================
def _approx(a, b, tol=1e-9):  # SOURCE: design choice (numerical tolerance)
    return abs(a - b) <= tol * max(1.0, abs(a), abs(b))


def _assert_conserved(pools, W0, tot0, label):
    assert _approx(total_water(pools), W0), f"{label}: water not conserved"
    for s in SOLUTES:
        assert _approx(total_solute(pools, s), tot0[s]), f"{label}: {s} not conserved"


def test_transfer_fluid_conserves():
    src = Pool(2.0, {'Sodium': 280.0, 'Urea': 40.0}); dst = Pool(1.0, {'Sodium': 140.0})
    W0 = total_water([src, dst]); tot0 = {s: total_solute([src, dst], s) for s in SOLUTES}
    transfer_fluid(src, dst, 0.5)
    _assert_conserved([src, dst], W0, tot0, "transfer_fluid")
    assert _approx(dst.solutes['Sodium'], 140.0 + 70.0)   # 0.5 L * 140 mM
    return "conserves water+solute; carries at src conc"


def test_transfer_clamps_to_source():
    src = Pool(0.3, {'Sodium': 42.0}); dst = Pool(1.0)
    moved = transfer_fluid(src, dst, 5.0)
    assert _approx(moved, 0.3) and _approx(src.volume, 0.0) and src.solutes['Sodium'] >= -1e-12
    return "clamps to source, no negative pools"


def test_transfer_solute_bidirectional():
    a = Pool(1.0, {'Urea': 10.0}); b = Pool(1.0, {'Urea': 2.0}); U0 = 12.0
    transfer_solute(a, b, 'Urea', 100.0); assert _approx(a.solutes['Urea'], 0.0) and _approx(b.solutes['Urea'], U0)
    transfer_solute(a, b, 'Urea', -100.0); assert _approx(b.solutes['Urea'], 0.0) and _approx(a.solutes['Urea'], U0)
    return "clamped both directions, conserves"


def test_lymph_isolation():
    layer = Pool(0.5, {'Sodium': 70.0, 'Urea': 10.0}); serum = Pool(14.0, {'Sodium': 1960.0, 'Urea': 280.0})
    pools = [layer, serum]; W0 = total_water(pools); tot0 = {s: total_solute(pools, s) for s in SOLUTES}
    moved = 0.0
    for _ in range(100):
        moved += lymph_flux(layer, serum, 0.001, 1.0, 1.0)
        _assert_conserved(pools, W0, tot0, "lymph")
    assert _approx(moved, min(0.001 * 1.0 * 100, 0.5))
    return f"conserved every step; moved {moved:.4f} L = formula"


def test_uf_isolation():
    """Hypertonic layer -> UF pulls water serum->layer; verify direction, sieving,
    conservation, and that UF reverses (back-filtration) when layer is hypotonic."""
    serum_vd = {s: 14.0 for s in SOLUTES}                 # flat Vd for the test
    serum = Pool(14.0, {'Sodium': 14.0 * 140})            # 140 mM
    layer = Pool(1.0, {'Sodium': 1.0 * 200})              # 200 mM -> hypertonic
    pools = [serum, layer]; W0 = total_water(pools); Na0 = total_solute(pools, 'Sodium')
    v_layer0 = layer.volume
    for _ in range(50):
        uf_flux(serum, layer, Kf_layer=8e-5, dP_hydro=0.0, dt=1.0, serum_vd=serum_vd)
        assert _approx(total_water(pools), W0), "UF broke water conservation"
        assert _approx(total_solute(pools, 'Sodium'), Na0), "UF broke Na conservation"
    assert layer.volume > v_layer0, "UF should move water INTO hypertonic layer"
    # hypotonic layer => back-filtration (water leaves layer)
    layer2 = Pool(1.0, {'Sodium': 1.0 * 100}); serum2 = Pool(14.0, {'Sodium': 14.0 * 140})
    uf_flux(serum2, layer2, 8e-5, 0.0, 1.0, serum_vd)
    assert layer2.volume < 1.0, "hypotonic layer should back-filtrate (water out)"
    return "UF moves water toward hypertonic side, sieved, conserved; reverses correctly"


def test_diffusion_isolation():
    """Urea diffuses serum->layer down its gradient and approaches equilibrium;
    conserved every step."""
    serum_vd = {s: 14.0 for s in SOLUTES}
    serum = Pool(14.0, {'Urea': 14.0 * 20})               # 20 mM
    layer = Pool(2.0, {'Urea': 0.0})                      # 0 mM
    pools = [serum, layer]; W0 = total_water(pools); U0 = total_solute(pools, 'Urea')
    mtac = {s: 0.0 for s in SOLUTES}; mtac['Urea'] = 0.02   # L/min
    for _ in range(2000):
        diffusion_flux(serum, layer, mtac, 1.0, serum_vd)
        assert _approx(total_water(pools), W0) and _approx(total_solute(pools, 'Urea'), U0), "diffusion not conserved"
    # near equilibrium: serum conc (Vd) ≈ layer conc
    assert abs(serum.conc('Urea', 14.0) - layer.conc('Urea')) < 1.0, "should approach equilibrium"
    assert layer.solutes['Urea'] > 0.0, "urea should have diffused into layer"
    return f"urea diffused down gradient to equilibrium; conserved"


def test_starling_isolation():
    serum_vd = {s: 14.0 for s in SOLUTES}; tissue_vd = {s: 28.0 for s in SOLUTES}
    serum = Pool(14.0, {'Sodium': 14.0 * 140}); tissue = Pool(28.0, {'Sodium': 28.0 * 150})  # tissue hypertonic
    pools = [serum, tissue]; W0 = total_water(pools); Na0 = total_solute(pools, 'Sodium')
    v_t0 = tissue.volume
    for _ in range(50):
        starling_flux(serum, tissue, 0.001, 1.0, serum_vd, tissue_vd)
        assert _approx(total_water(pools), W0) and _approx(total_solute(pools, 'Sodium'), Na0), "starling not conserved"
    assert tissue.volume > v_t0, "water should move serum->tissue (tissue hypertonic)"
    return "Starling water flux toward hypertonic tissue, conserved"


def test_fickian_isolation():
    lo = Pool(1.0, {'Glucose': 0.0}); hi = Pool(1.0, {'Glucose': 100.0})   # gradient
    pools = [lo, hi]; W0 = total_water(pools); G0 = total_solute(pools, 'Glucose')
    dC0 = abs(hi.conc('Glucose') - lo.conc('Glucose'))
    prev = dC0
    for _ in range(2000):
        fickian_flux(lo, hi, 7.3e-4, 100.0, 3.0, 1.0)
        assert _approx(total_water(pools), W0) and _approx(total_solute(pools, 'Glucose'), G0), "fickian not conserved"
        dC = abs(hi.conc('Glucose') - lo.conc('Glucose'))
        assert dC <= prev + 1e-12, "gradient must shrink monotonically (down-gradient)"
        prev = dC
    assert dC < 0.05 * dC0, "should approach equilibrium given enough time"
    return "inter-layer Fickian mixing: monotone down-gradient to equilibrium, conserved"


def test_ro_isolation():
    """Feed split into permeate+reject; conserved; permeate cleaner & reject
    concentrated; flux-feasible cap zeroes permeate for hypertonic feed."""
    ro_sp, ref_jw = ro_nominal_perm()
    feed = Pool(1.0, {s: v for s, v in {'Sodium': 140.0, 'Chloride': 110.0, 'Urea': 20.0,
                'Glucose': 5.0, 'Bicarbonate': 30.0}.items()})
    perm = Pool(0.0); rej = Pool(0.0)
    W0 = total_water([feed, perm, rej]); tot0 = {s: total_solute([feed, perm, rej], s) for s in SOLUTES}
    cNa_feed = feed.conc('Sodium')
    ro_split(feed, perm, rej, recovery=0.85, ro_pressure=15.0, ro_wp=5.0, ro_sp=ro_sp, ref_jw=ref_jw)
    _assert_conserved([feed, perm, rej], W0, tot0, "RO")
    assert _approx(feed.volume, 0.0), "feed must be fully partitioned"
    assert perm.volume > 0 and perm.conc('Sodium') < cNa_feed, "permeate should be cleaner than feed"
    assert rej.conc('Sodium') > cNa_feed, "reject should be concentrated"
    # flux-feasible: very hypertonic feed (osmotic > applied 15 bar) => no permeate
    feed2 = Pool(1.0, {'Sodium': 400.0, 'Chloride': 400.0, 'Urea': 80.0, 'Glucose': 60.0})
    p2 = Pool(0.0); r2 = Pool(0.0)
    pv = ro_split(feed2, p2, r2, 0.85, 15.0, 5.0, ro_sp, ref_jw)
    assert _approx(pv, 0.0) and _approx(p2.volume, 0.0), "hypertonic feed must yield ~0 permeate (flux-feasible)"
    return "RO partitions feed->permeate+reject, conserved; permeate clean, reject concentrated; cap works"


def test_external_isolation():
    serum_vd = {s: 14.0 for s in SOLUTES}
    serum = Pool(14.0, {'Sodium': 14.0 * 140, 'Urea': 14.0 * 20})
    urine = Pool(0.0); lost = Pool(0.0); intake_src = Pool(1.0, {'Sodium': 100.0})
    pools = [serum, urine, lost, intake_src]
    W0 = total_water(pools); tot0 = {s: total_solute(pools, s) for s in SOLUTES}
    for _ in range(60):
        renal_clearance(serum, urine, 0.0003, 0.6, 1.0, serum_vd)   # ~0.43 L/day
        insensible_loss(serum, lost, 0.0005, 1.0)
        oral_intake(intake_src, serum, 0.005, {'Sodium': 100.0})
        _assert_conserved(pools, W0, tot0, "external")
    assert urine.volume > 0 and lost.volume > 0 and serum.solutes['Sodium'] != tot0['Sodium']
    assert _approx(lost.solutes['Sodium'], 0.0), "insensible loss carries no solute"
    return "renal/insensible/intake all conserved across body+external pools"


def test_cori_reaction():
    serum = Pool(14.0, {'Lactate': 30.0, 'Bicarbonate': 20.0})
    pair0 = serum.solutes['Lactate'] + serum.solutes['Bicarbonate']
    converted = cori_react(serum, 12.0)
    assert _approx(converted, 12.0) and _approx(serum.solutes['Lactate'], 18.0) and _approx(serum.solutes['Bicarbonate'], 32.0)
    assert _approx(serum.solutes['Lactate'] + serum.solutes['Bicarbonate'], pair0), "Cori conserves total anion"
    over = cori_react(serum, 999.0); assert _approx(over, 18.0) and _approx(serum.solutes['Lactate'], 0.0)
    return "Cori converts 1 Lactate->1 HCO3, clamped, conserves total osmoles+charge"


if __name__ == '__main__':
    tests = [test_transfer_fluid_conserves, test_transfer_clamps_to_source,
             test_transfer_solute_bidirectional, test_lymph_isolation,
             test_uf_isolation, test_diffusion_isolation, test_starling_isolation,
             test_fickian_isolation, test_ro_isolation, test_external_isolation,
             test_cori_reaction]
    ok = True
    for t in tests:
        try:
            print(f"  [PASS] {t.__name__}: {t()}")
        except AssertionError as e:
            ok = False; print(f"  [FAIL] {t.__name__}: {e}")
    print("RESULT:", "ALL PASS" if ok else "FAILURES")
    import sys
    sys.exit(0 if ok else 1)