"""Independent numerical reference. Requires numpy, pandas and scipy.

Usage: python reproduce.py /path/to/Puromycin.csv
Input values must match the documented source checksum. No Gradiance imports.
"""
import argparse
import hashlib
import json
from pathlib import Path

import numpy as np
import pandas as pd
import scipy
from scipy.optimize import curve_fit

parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("source", type=Path, help="Original Puromycin.csv from the article")
args = parser.parse_args()
source = args.source
if not source.is_file():
    parser.error(f"Input file not found: {source}")
checksum = hashlib.sha256(source.read_bytes()).hexdigest()
if checksum != "389be46927d2217b13d48b5b9ec56374007a87ea607c8e6f004c17683c427c70":
    parser.error("Input SHA-256 does not match the published Puromycin source. Download the original CSV linked in the article.")
data = pd.read_csv(source)
results = {}


def michaelis_menten(concentration, vmax, km):
    return vmax * concentration / (km + concentration)


for state, observations in data.groupby("state"):
    x = observations.conc.to_numpy()
    y = observations.rate.to_numpy()
    fitted, _ = curve_fit(
        michaelis_menten, x, y, p0=[200, 0.05], bounds=(0, np.inf), maxfev=10000
    )
    residual = y - michaelis_menten(x, *fitted)
    results[state] = {
        "n": len(observations),
        "vmax": float(fitted[0]),
        "km": float(fitted[1]),
        "r_squared": float(1 - (residual**2).sum() / ((y - y.mean())**2).sum()),
        "rmse": float(np.sqrt(np.mean(residual**2))),
    }
print(json.dumps({"sha256": checksum, "scipy": scipy.__version__, "numpy": np.__version__, "results": results}, indent=2))