# Weighted 4PL standard curve and interpolation of samples. Python 3.13, NumPy, SciPy, pandas.
import numpy as np, pandas as pd
from scipy import optimize

d = pd.read_csv("plate.csv")                         # columns: type, id, conc_pg_ml, od; one row per well
st = d[d.type == "standard"]; x, y = st.conc_pg_ml.to_numpy(), st.od.to_numpy()
def f4(x, bottom, top, ec50, slope): return top + (bottom - top) / (1 + (x / ec50) ** slope)
def inv(yv, bottom, top, ec50, slope): return ec50 * ((bottom - top) / (yv - top) - 1) ** (1 / slope)
p, _ = optimize.curve_fit(f4, x, y, p0=[0.05, 3, 100, 1], sigma=y, maxfev=20000, xtol=1e-12, ftol=1e-12)   # sigma = y: weights 1/OD^2
print("Weighted 4PL: bottom {:.4f}, top {:.4f}, EC50 {:.2f} pg/mL, slope {:.4f}".format(*p))

lev = st.groupby("conc_pg_ml").od.mean().sort_index(ascending=False)
lloq, uloq = lev.index.min(), lev.index.max()
cv = st.groupby("conc_pg_ml").od.agg(lambda v: 100 * v.std(ddof=1) / v.mean())
print("Back-calculated standards (mean of duplicate wells), % of nominal; CV of the duplicate ODs:")
for c, od in lev.items():
    acc = 100 * inv(od, *p) / c; lim = 25 if c in (lloq, uloq) else 20
    print(f"  {c:8.2f} pg/mL: {acc:6.1f}%  {'ok' if abs(acc - 100) <= lim else 'OUTSIDE'} (limit +/-{lim}%); CV {cv[c]:.1f}%")

b1, b0 = np.polyfit(x, y, 1)                         # for comparison: a straight line through all standards
print("Straight line through all standards, back-calculated % of nominal:")
print("  " + ", ".join(f"{c:g}: {100 * ((od - b0) / b1) / c:.0f}%" for c, od in lev.items()))

lo_od, hi_od = f4(lloq, *p), f4(uloq, *p)
flat = p[0] + 0.75 * (p[1] - p[0])                  # example rule set in advance: upper quarter of the OD range = plateau
print(f"Samples (mean OD of duplicate wells), reportable range {lloq:g} to {uloq:g} pg/mL:")
for k, od in d[d.type == "sample"].groupby("id").od.mean().items():
    if od < lo_od: print(f"  {k}: OD {od:.4f} -> below the LLOQ (report as < {lloq:g} pg/mL)")
    elif od > hi_od: print(f"  {k}: OD {od:.4f} -> above the ULOQ (dilute and re-assay)")
    elif od > flat: print(f"  {k}: OD {od:.4f} -> {inv(od, *p):.1f} pg/mL, but on the upper plateau (dilute and re-assay)")
    else: print(f"  {k}: OD {od:.4f} -> {inv(od, *p):.1f} pg/mL")