"""ST000284 serum metabolomics: prespecified, leakage-controlled reanalysis.
Offline default: python analysis.py --input inputs/ST000284_AN000452.json --output rerun
Full run uses 3x5 outer / 3 inner CV, 99 permutations, 1000 bootstraps.
No package installation or remote execution. See analysis_protocol.md.
"""
from pathlib import Path
import argparse, collections, datetime, hashlib, importlib.metadata, json, sys, time, warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.feature_selection import VarianceThreshold
from sklearn.linear_model import LogisticRegression
from sklearn.dummy import DummyClassifier
from sklearn.model_selection import StratifiedKFold, GridSearchCV
from sklearn.metrics import (roc_auc_score, average_precision_score, balanced_accuracy_score,
    recall_score, brier_score_loss, confusion_matrix, roc_curve)
from sklearn.exceptions import ConvergenceWarning
from threadpoolctl import threadpool_limits

SEEDS = [17, 29, 43]
C_GRID = [0.01, 0.1, 1.0, 10.0]
GROUPS = ["CRC", "Polyp", "Healthy"]
COLORS = {"CRC": "#B34A00", "Polyp": "#0072B2", "Healthy": "#007A59"}
INPUT_SHA256 = "cba77023eff08616c974c21d3e034468f601c31c1646bed12c6cc87b099ed7a9"

def write_json(path, obj):
    Path(path).write_text(json.dumps(obj, indent=2, default=lambda x: x.item() if isinstance(x, np.generic) else str(x)), encoding="utf-8")

def load_inputs(path):
    path = Path(path)
    checksum = hashlib.sha256(path.read_bytes()).hexdigest()
    if checksum != INPUT_SHA256:
        raise ValueError("Input checksum differs from the frozen deposit; review before analysis.")
    data = json.loads(path.read_text(encoding="utf-8"))
    metadata = pd.DataFrame([{"sample_id": r["Sample ID"], "subject_id": r["Subject ID"],
        "group": r["Factors"]["Patient group"], **r.get("Additional sample data", {})}
        for r in data["SUBJECT_SAMPLE_FACTORS"]]).set_index("sample_id")
    matrix = pd.DataFrame(data["MS_METABOLITE_DATA"]["Data"]).set_index("Metabolite").T
    assert matrix.index.is_unique and matrix.columns.is_unique and metadata.index.is_unique
    assert set(matrix.index) == set(metadata.index)
    matrix = matrix.loc[metadata.index].apply(pd.to_numeric, errors="coerce")
    assert matrix.shape == (224, 113) and metadata.group.value_counts().to_dict() == {"Healthy":84,"Polyp":76,"CRC":64}
    assert not matrix.duplicated().any(), "Duplicate profiles require a grouped protocol."
    assert metadata.subject_id.eq("-").all(), "Reassess grouping if identities change."
    return matrix, metadata, data

class PositiveFeatureFilter(BaseEstimator, TransformerMixin):
    def __init__(self, max_missing=0.20, floor_missing=False):
        self.max_missing = max_missing
        self.floor_missing = floor_missing
    def clean(self, X):
        arr = np.asarray(X, dtype=float).copy()
        bad = ~np.isfinite(arr) | (arr <= 0)
        if self.floor_missing:
            bad |= arr == 1
        arr[bad] = np.nan
        return arr
    def fit(self, X, y=None):
        arr = self.clean(X)
        self.keep_ = np.mean(np.isnan(arr), axis=0) <= self.max_missing
        if not self.keep_.any():
            raise ValueError("No features pass training-only missingness filter.")
        return self
    def transform(self, X):
        return self.clean(X)[:, self.keep_]

def make_pipeline(kind="metabolites"):
    if kind in ("metabolites", "floor_missing"):
        steps = [("valid", PositiveFeatureFilter(floor_missing=kind=="floor_missing")),
                 ("impute", SimpleImputer(strategy="median", keep_empty_features=True)),
                 ("log10", FunctionTransformer(np.log10))]
    else:
        steps = [("impute", SimpleImputer(strategy="median", keep_empty_features=True))]
    steps += [("variance", VarianceThreshold()), ("scale", StandardScaler()),
              ("model", LogisticRegression(C=1., solver="liblinear", max_iter=2000,
                                            tol=1e-6, random_state=20260929))]
    return Pipeline(steps)

def summarize_metrics(y, probability):
    predicted = np.asarray(probability) >= 0.5
    tn, fp, fn, tp = confusion_matrix(y, predicted, labels=[0,1]).ravel()
    return {"roc_auc":roc_auc_score(y, probability), "average_precision":average_precision_score(y, probability),
            "balanced_accuracy":balanced_accuracy_score(y, predicted),
            "sensitivity":tp/(tp+fn), "specificity":tn/(tn+fp),
            "brier":brier_score_loss(y, probability), "tn":int(tn), "fp":int(fp), "fn":int(fn), "tp":int(tp)}

def run_nested(X, y, name, kind="metabolites", seeds=SEEDS, record=True):
    X = X.copy()
    y = np.asarray(y, dtype=int)
    assert X.index.is_unique and len(X)==len(y)
    predictions, metrics, settings, inner_records = [], [], [], []
    with threadpool_limits(limits=1):
        for rep, seed in enumerate(seeds):
            outer = StratifiedKFold(n_splits=5, shuffle=True, random_state=seed)
            coverage = np.zeros(len(y), dtype=int)
            repeat_p = np.full(len(y), np.nan)
            dummy_p = np.full(len(y), np.nan)
            for fold, (train, test) in enumerate(outer.split(X, y)):
                assert len(set(train) & set(test)) == 0
                assert set(X.index[train]).isdisjoint(X.index[test])
                inner_seed = 1000 + 10*rep + fold
                inner = list(StratifiedKFold(n_splits=3, shuffle=True, random_state=inner_seed).split(X.iloc[train],y[train]))
                if record:
                    for inner_fold, (itr, ite) in enumerate(inner):
                        assert set(itr).isdisjoint(ite)
                        for idx in train[ite]:
                            inner_records.append({"analysis":name,"repeat":rep+1,"outer_fold":fold+1,
                                                  "inner_fold":inner_fold+1,"sample_id":X.index[idx]})
                search = GridSearchCV(make_pipeline(kind), {"model__C":C_GRID}, scoring="roc_auc",
                                      cv=inner, n_jobs=1, refit=True, error_score="raise")
                with warnings.catch_warnings(record=True) as caught:
                    warnings.simplefilter("always", ConvergenceWarning)
                    search.fit(X.iloc[train], y[train])
                convergence = [str(w.message) for w in caught if issubclass(w.category, ConvergenceWarning)]
                p = search.predict_proba(X.iloc[test])[:,1]
                dummy = DummyClassifier(strategy="prior").fit(X.iloc[train],y[train])
                dp = dummy.predict_proba(X.iloc[test])[:,1]
                assert np.isfinite(p).all() and ((p>=0)&(p<=1)).all()
                repeat_p[test], dummy_p[test] = p, dp
                coverage[test] += 1
                if record:
                    settings.append({"analysis":name,"repeat":rep+1,"outer_fold":fold+1,"seed":seed,
                         "inner_seed":inner_seed,"n_train":len(train),"n_test":len(test),
                         "C":search.best_params_["model__C"],"inner_auc":search.best_score_,
                         "features_retained":int(search.best_estimator_.named_steps["variance"].get_support().sum()),
                         "convergence_warning_count":len(convergence), "warnings":" | ".join(convergence)})
                    for idx, pred, base in zip(test,p,dp):
                        predictions.append({"analysis":name,"repeat":rep+1,"outer_fold":fold+1,
                            "sample_id":X.index[idx],"y":int(y[idx]),"probability":pred,"dummy_probability":base})
            assert (coverage==1).all() and np.isfinite(repeat_p).all()
            metrics.append({"analysis":name,"repeat":rep+1,"model":kind,"n":len(y),
                            "n_crc":int(y.sum()),"prevalence":y.mean(),**summarize_metrics(y,repeat_p)})
            metrics.append({"analysis":name,"repeat":rep+1,"model":"dummy","n":len(y),
                            "n_crc":int(y.sum()),"prevalence":y.mean(),**summarize_metrics(y,dummy_p)})
    return {"predictions":pd.DataFrame(predictions),"metrics":pd.DataFrame(metrics),
            "settings":pd.DataFrame(settings),"inner_splits":pd.DataFrame(inner_records)}

def audit_data(X, meta, out):
    out=Path(out); out.mkdir(parents=True,exist_ok=True)
    finite=np.isfinite(X); invalid=~finite | X.le(0); floor=X.eq(1)
    valid=X.mask(invalid)
    log=np.log10(valid)
    medlog=log.median(axis=1)
    mad=float(np.median(np.abs(medlog-np.median(medlog))))
    robust_z=(medlog-np.median(medlog))/(1.4826*mad) if mad>0 else medlog*0
    sample=pd.DataFrame({"group":meta.group,"invalid_n":invalid.sum(axis=1),
        "exact_one_n":floor.sum(axis=1),"median_log10_cps":medlog,"total_cps":valid.sum(axis=1),
        "median_log_robust_z":robust_z,"outlier_flag_abs_z_gt_3_5":robust_z.abs()>3.5})
    sample.index.name="sample_id"; sample.to_csv(out/"sample_quality.csv")
    feature=pd.DataFrame({"invalid_n":invalid.sum(),"exact_one_n":floor.sum(),"unique_values":X.nunique(),
                         "minimum":X.min(),"median":X.median(),"maximum":X.max(),"log10_sd":log.std()})
    feature.index.name="metabolite"; feature.to_csv(out/"feature_quality.csv")
    meta.to_csv(out/"sample_metadata.csv")
    X.to_csv(out/"metabolite_matrix.csv",index_label="sample_id")
    ledger=pd.DataFrame({"group":meta.group,"included_primary":meta.group.isin(["CRC","Polyp"]),
            "included_healthy_comparison":meta.group.isin(["CRC","Healthy"]),
            "included_pooled_comparison":True,"technical_exclusion":False,
            "primary_reason":np.where(meta.group.eq("Healthy"),"Healthy excluded from primary disease-control comparison","Included")})
    ledger.to_csv(out/"cohort_ledger.csv",index_label="sample_id")
    clean=meta.replace({"-":np.nan,"":np.nan})
    clean.groupby("group").count().to_csv(out/"metadata_nonmissing_counts.csv")
    numeric=[]
    for column in ["Age at  Consent","BMI [kg/m²]","Height [cm]","Weight [kg]"]:
        tab=pd.to_numeric(meta[column],errors="coerce").groupby(meta.group).agg(["count","mean","std","median","min","max"]).reset_index()
        tab.insert(0,"variable",column);numeric.append(tab)
    pd.concat(numeric).to_csv(out/"clinical_numeric_summary.csv",index=False)
    for column,short in [("Gender","sex"),("Stage","stage"),("Smoking condition","smoking"),("Alcohol consumption","alcohol")]:
        pd.crosstab(meta.group,meta[column]).to_csv(out/(short+"_counts.csv"))
    summary={"n_samples":len(X),"n_features":X.shape[1],"groups":meta.group.value_counts().to_dict(),
        "invalid_cells":int(invalid.sum().sum()),"duplicate_samples":int(X.duplicated(keep=False).sum()),
        "duplicate_features":int(X.T.duplicated(keep=False).sum()),"constant_features":int((X.nunique()<=1).sum()),
        "exact_one_cells":int(floor.sum().sum()),"exact_one_samples":int(floor.any(axis=1).sum()),
        "exact_one_features":int(floor.any(axis=0).sum()),"outlier_sample_ids":list(X.index[robust_z.abs()>3.5]),
        "subject_ids_available":False,"batch_site_dates_available":False,"excluded_samples":0,
        "upstream_imputation_normalization_feature_selection":"Not fully verified; exact 1 values may encode a floor.",
        "paper_vs_deposit":{"paper_samples":234,"deposited_samples":224,"paper_monitored_metabolites":158,"deposited_metabolites":113}}
    write_json(out/"quality_summary.json",summary)
    with plt.rc_context({"font.size":10,"axes.spines.top":False,"axes.spines.right":False}):
        fig,axes=plt.subplots(2,2,figsize=(10,7),layout="constrained")
        rng=np.random.default_rng(721)
        fields=[medlog, floor.sum(axis=1),pd.to_numeric(meta["Age at  Consent"],errors="coerce"),pd.to_numeric(meta["BMI [kg/m²]"],errors="coerce")]
        labels=["Median log10 metabolite intensity (cps)","Count of exact-1 entries","Age (years)","BMI (kg/m²)"]
        for ax,field,label,letter in zip(axes.flat,fields,labels,"ABCD"):
            for j,g in enumerate(GROUPS):
                vals=field[meta.group.eq(g)].dropna()
                ax.scatter(j+rng.uniform(-.17,.17,len(vals)),vals,s=17,color=COLORS[g],alpha=.65,marker=["o","s","^"][j],linewidths=0)
                ax.plot([j-.22,j+.22],[vals.median()]*2,color="black",lw=2)
            ticks=[f"{g}\nn={int(field[meta.group.eq(g)].notna().sum())}" for g in GROUPS]
            ax.set(xticks=range(3),xticklabels=ticks,ylabel=label,title=f"{letter}  {label}")
            if letter=="B":ax.set_ylim(bottom=-.15)
        fig.suptitle("ST000284 quality and cohort audit — all 224 deposited samples",fontsize=13)
        fig.savefig(out/"data_quality.png",dpi=180,facecolor="white")
        plt.close(fig)
    return summary

def fit_apparent(X,y):
    with threadpool_limits(limits=1):
        cv=StratifiedKFold(3,shuffle=True,random_state=1000)
        search=GridSearchCV(make_pipeline(),{"model__C":C_GRID},scoring="roc_auc",cv=cv,n_jobs=1,error_score="raise")
        with warnings.catch_warnings(record=True) as caught:
            warnings.simplefilter("always",ConvergenceWarning); search.fit(X,y)
        p=search.predict_proba(X)[:,1]
    return {"chosen_C":search.best_params_["model__C"],
            "convergence_warning_count":sum(issubclass(w.category,ConvergenceWarning) for w in caught),
            **summarize_metrics(y,p)}

def evaluate_all(X,meta,out,primary_only=False):
    out=Path(out);out.mkdir(parents=True,exist_ok=True)
    specs=[("CRC_vs_Polyp",["CRC","Polyp"],"metabolites")]
    if not primary_only:
        specs += [("CRC_vs_Healthy",["CRC","Healthy"],"metabolites"),
                  ("CRC_vs_Pooled",GROUPS,"metabolites"),
                  ("clinical_only",["CRC","Polyp"],"clinical"),
                  ("floor_as_missing",["CRC","Polyp"],"floor_missing"),
                  ("floor_indicator_only",["CRC","Polyp"],"indicators")]
    results=[]
    for name,groups,kind in specs:
        keep=meta.group.isin(groups); xm=X.loc[keep]; mm=meta.loc[keep]; y=mm.group.eq("CRC").astype(int).to_numpy()
        if kind=="clinical":
            xm=pd.DataFrame({"age":pd.to_numeric(mm["Age at  Consent"],errors="coerce"),
                             "male":mm.Gender.map({"M":1.,"F":0.}),
                             "bmi":pd.to_numeric(mm["BMI [kg/m²]"],errors="coerce")},index=mm.index)
        elif kind=="indicators":
            xm=xm.eq(1).astype(float)
        started=time.monotonic()
        res=run_nested(xm,y,name,kind)
        results.append(res)
        print(name,res["metrics"].loc[res["metrics"].model.ne("dummy"),["repeat","roc_auc","balanced_accuracy"]].round(4).to_dict("records"),f"seconds={time.monotonic()-started:.1f}",flush=True)
    combined={key:pd.concat([r[key] for r in results],ignore_index=True) for key in results[0]}
    for key,value in combined.items():value.to_csv(out/(key+".csv"),index=False)
    mask=meta.group.isin(["CRC","Polyp"]);y=meta.loc[mask].group.eq("CRC").astype(int).to_numpy()
    apparent=fit_apparent(X.loc[mask],y)
    write_json(out/"apparent_training_metrics.json",apparent)
    # Save stage/control-specific behavior for first repeat without retraining.
    first=combined["predictions"].query("repeat == 1").copy()
    first["group"]=first.sample_id.map(meta.group)
    first["stage"]=first.sample_id.map(meta.Stage)
    subgroup=[]
    for (name,g), frame in first.groupby(["analysis","group"]):
        subgroup.append({"analysis":name,"group":g,"n":len(frame),"fraction_predicted_CRC":float((frame.probability>=.5).mean())})
    for stage,frame in first.query("analysis == 'CRC_vs_Polyp' and y == 1").groupby("stage"):
        subgroup.append({"analysis":"CRC_vs_Polyp","group":stage,"n":len(frame),"fraction_predicted_CRC":float((frame.probability>=.5).mean())})
    pd.DataFrame(subgroup).to_csv(out/"subgroup_threshold_behavior.csv",index=False)
    return combined,apparent

def uncertainty(X,meta,predictions,out,n_perm=99,n_boot=1000,cutoff=None):
    out=Path(out)
    primary=predictions.query("analysis == 'CRC_vs_Polyp' and repeat == 1").set_index("sample_id")
    mask=meta.group.isin(["CRC","Polyp"]); xp=X.loc[mask]; yp=meta.loc[mask].group.eq("CRC").astype(int).to_numpy()
    pp=primary.loc[xp.index,"probability"].to_numpy()
    assert np.array_equal(primary.loc[xp.index,"y"],yp)
    observed=summarize_metrics(yp,pp)
    rng=np.random.default_rng(20260929);idxs=[np.flatnonzero(yp==k) for k in [0,1]]
    boot=[]
    for b in range(n_boot):
        ids=np.concatenate([rng.choice(i,size=len(i),replace=True) for i in idxs])
        boot.append({"bootstrap":b+1,**summarize_metrics(yp[ids],pp[ids])})
    boot=pd.DataFrame(boot);boot.to_csv(out/"bootstrap_metrics.csv",index=False)
    intervals=pd.DataFrame([{"metric":m,"estimate":observed[m],"lower_95":boot[m].quantile(.025),"upper_95":boot[m].quantile(.975),
         "method":"1000 stratified sample bootstraps of fixed first-repeat OOF predictions; conditional approximate"}
         for m in ["roc_auc","average_precision","balanced_accuracy","sensitivity","specificity","brier"]])
    intervals.to_csv(out/"conditional_intervals.csv",index=False)
    rng=np.random.default_rng(20260930);perms=[]
    for pidx in range(n_perm):
        if cutoff and datetime.datetime.now(datetime.timezone.utc)>=datetime.datetime.fromisoformat(cutoff):
            break
        py=rng.permutation(yp)
        result=run_nested(xp,py,"permutation","metabolites",seeds=[17],record=False)
        score=float(result["metrics"].query("model == 'metabolites'").roc_auc.iloc[0])
        perms.append({"permutation":pidx+1,"roc_auc":score})
        pd.DataFrame(perms).to_csv(out/"permutation_auc.csv",index=False)
        if (pidx+1)%10==0: print(f"Permutations completed: {pidx+1}/{n_perm}",flush=True)
    perm_df=pd.DataFrame(perms,columns=["permutation","roc_auc"])
    perm_df.to_csv(out/"permutation_auc.csv",index=False)
    complete=len(perms)==n_perm
    null_ge=int((perm_df.roc_auc>=observed["roc_auc"]).sum()) if len(perms) else 0
    summary={"observed_first_repeat_auc":observed["roc_auc"],"planned_permutations":n_perm,"completed_permutations":len(perms),
             "permutation_complete":complete,"null_at_least_observed":null_ge,
             "permutation_p_plus_one":(1+null_ge)/(1+len(perms)) if complete and len(perms)>0 else None,
             "planned_bootstraps":n_boot,"completed_bootstraps":len(boot),
             "independence":"Sample independence and exchangeability assumed; participant/group identities absent.",
             "interval_limit":"Conditional OOF bootstrap does not include complete training-process uncertainty or external transport uncertainty."}
    write_json(out/"uncertainty_summary.json",summary)
    bins=pd.qcut(pp,q=5,duplicates="drop")
    calibration=pd.DataFrame({"probability":pp,"y":yp,"bin":bins}).groupby("bin",observed=True).agg(n=("y","size"),mean_prediction=("probability","mean"),observed_fraction=("y","mean"))
    calibration.to_csv(out/"descriptive_calibration.csv")
    return intervals,perm_df,summary

def validation_figure(results,intervals,permutations,apparent,out):
    out=Path(out);pred=results["predictions"];metrics=results["metrics"]
    primary=pred.query("analysis == 'CRC_vs_Polyp' and repeat == 1")
    with plt.rc_context({"font.size":10,"axes.spines.top":False,"axes.spines.right":False}):
        fig,axes=plt.subplots(2,2,figsize=(11,8),layout="constrained")
        ax=axes[0,0]
        for name,color,style in [("CRC_vs_Polyp","#B34A00","-"),("clinical_only","#0072B2","--"),("floor_indicator_only","#007A59",":")]:
            f=pred.query("analysis == @name and repeat == 1")
            if len(f):
                fpr,tpr,_=roc_curve(f.y,f.probability)
                ax.plot(fpr,tpr,color=color,ls=style,label=f"{name.replace('_',' ')}: {roc_auc_score(f.y,f.probability):.3f}")
        ax.plot([0,1],[0,1],color="gray",lw=1,ls="--")
        ax.set(xlim=(0,1),ylim=(0,1.02),xlabel="False positive rate",ylabel="True positive rate",title="A  First-repeat held-out ROC (CRC vs polyp)")
        ax.legend(fontsize=8,loc="lower right")
        ax=axes[0,1]
        names=["CRC_vs_Polyp","CRC_vs_Healthy","CRC_vs_Pooled","clinical_only","floor_as_missing","floor_indicator_only"]
        present=[n for n in names if n in metrics.analysis.values]
        for j,name in enumerate(present):
            vals=metrics.loc[metrics.analysis.eq(name)&metrics.model.ne("dummy"),"roc_auc"].to_numpy()
            ax.scatter(vals,j+np.linspace(-.12,.12,len(vals)),marker="o",s=24,color="#0072B2")
        ci=intervals.query("metric == 'roc_auc'").iloc[0]
        ax.errorbar(ci.estimate,0,xerr=[[ci.estimate-ci.lower_95],[ci.upper_95-ci.estimate]],fmt="s",color="#B34A00",capsize=3,label="Primary conditional 95% interval")
        ax.axvline(.5,color="gray",ls="--",lw=1)
        ax.set(xlim=(0,1),yticks=range(len(present)),yticklabels=[n.replace("_"," ") for n in present],xlabel="ROC AUC",title="B  Three repeat estimates (not independent)")
        ax.invert_yaxis();ax.legend(fontsize=7,loc="upper left")
        ax=axes[1,0]
        if len(permutations):ax.hist(permutations.roc_auc,bins=np.linspace(0,1,21),color="#8BA9BA",edgecolor="white",label=f"{len(permutations)} permutations")
        ax.axvline(ci.estimate,color="#B34A00",lw=2,label=f"Observed {ci.estimate:.3f}")
        ax.set(xlim=(0,1),xlabel="First-repeat nested-CV ROC AUC",ylabel="Count",title="C  Full-procedure label-permutation null")
        ax.legend(fontsize=8)
        ax=axes[1,1];rng=np.random.default_rng(93)
        for j,(label,group) in enumerate([(0,"Polyp"),(1,"CRC")]):
            vals=primary.loc[primary.y.eq(label),"probability"]
            ax.scatter(j+rng.uniform(-.16,.16,len(vals)),vals,s=20,alpha=.65,color=COLORS[group],marker=["s","o"][j])
        ax.axhline(.5,color="gray",ls="--",lw=1)
        ax.set(xticks=[0,1],xticklabels=["Polyp (n=76)","CRC (n=64)"],ylim=(-.02,1.02),ylabel="Out-of-fold predicted CRC probability",title="D  Predictions at the fixed 0.5 threshold")
        fig.suptitle(f"ST000284: internal validation only | apparent training AUC {apparent['roc_auc']:.3f}",fontsize=13)
        fig.savefig(out/"validation_results.png",dpi=180,facecolor="white");plt.close(fig)

def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--input",default="inputs/ST000284_AN000452.json")
    parser.add_argument("--output",default="rerun")
    parser.add_argument("--permutations",type=int,default=99)
    parser.add_argument("--cutoff",default=None,help="Optional UTC ISO timestamp; incomplete permutation tests are flagged.")
    args=parser.parse_args()
    out=Path(args.output)
    out.mkdir(parents=True,exist_ok=False)
    X,meta,data=load_inputs(args.input)
    audit_data(X,meta,out)
    result,apparent=evaluate_all(X,meta,out)
    intervals,permutations,summary=uncertainty(X,meta,result["predictions"],out,n_perm=args.permutations,cutoff=args.cutoff)
    validation_figure(result,intervals,permutations,apparent,out)
    write_json(out/"runtime.json",{"python":sys.version,"packages":{p:importlib.metadata.version(p) for p in ["numpy","pandas","scipy","scikit-learn","matplotlib","threadpoolctl"]},"input_sha256":INPUT_SHA256})
    print(json.dumps(summary,indent=2))

if __name__ == "__main__":
    main()
