"""Patient-report agreement development evaluator, never a clinical validity score.

References are repeated, deidentified reports; no optical causality is inferred
from agreement. Screen matching and recall are separate exploratory strata.
"""

from functools import lru_cache
import json
import math
from pathlib import Path

import numpy as np

VERSION="0.4.0-patient-draft"
KINDS={"circle","ellipse","arc","freehand"}


def finite(value,label,lo=-30.,hi=30.):
    if isinstance(value,bool) or not isinstance(value,(float,int)) or not math.isfinite(value) or not lo<=value<=hi:
        raise ValueError(f"{label}: finite number in [{lo},{hi}] required")
    return float(value)


def target_key(record):
    return "|".join(str(record[k]) for k in ("patient_id","visit_id","eye","condition_id"))


def validate_feature(feature,source_ids):
    if not isinstance(feature,dict):raise ValueError("feature must be an object")
    if feature.get("kind") not in KINDS:raise ValueError("unknown feature kind")
    if feature.get("source_id") not in source_ids:raise ValueError("feature source_id not in source list")
    finite(feature.get("width_deg"),"reported width",.005,10.)
    if feature["kind"]=="freehand":
        points=feature.get("points_deg",[])
        if not 2<=len(points)<=256:raise ValueError("freehand needs 2–256 points")
        for xy in points:
            if not isinstance(xy,list) or len(xy)!=2:raise ValueError("invalid contour point")
            for v in xy:finite(v,"contour coordinate")
    else:
        for key in ("cx_deg","cy_deg"):finite(feature.get(key),key)
        finite(feature.get("rx_deg"),"radius x",.005,30.)
        finite(feature.get("ry_deg",feature["rx_deg"]),"radius y",.005,30.)
        finite(feature.get("rotation_deg",0.),"rotation",-360.,360.)
        if feature["kind"]=="arc":
            finite(feature.get("start_deg",0.),"arc start",-360.,360.)
            finite(feature.get("sweep_deg"),"arc sweep",1.,360.)
    brightness=feature.get("brightness8")
    if brightness is not None:
        if not isinstance(brightness,list) or len(brightness)!=8:raise ValueError("brightness8 needs 8 entries or null")
        for v in brightness:
            if v is not None:finite(v,"ordinal brightness",0.,4.)


def validate_record(record):
    if record.get("schema_version")!=VERSION:raise ValueError("unsupported patient schema")
    for key in ("record_id","patient_id","visit_id","condition_id"):
        value=record.get(key)
        if not isinstance(value,str) or not 1<=len(value)<=80 or "|" in value:raise ValueError(f"invalid {key}")
    if type(record.get("repeat_id")) is not int or record["repeat_id"]<1:raise ValueError("positive integer repeat_id required")
    if record.get("eye") not in ("OD","OS","OU"):raise ValueError("OD/OS/OU required")
    if record.get("origin") not in ("synthetic","research"):raise ValueError("explicit origin required")
    if record.get("status") not in ("draft","complete"):raise ValueError("draft/complete status required")
    if record.get("method") not in ("external_standard_light","recall","screen_matching"):raise ValueError("explicit observation method required")
    frame=record.get("frame_deg")
    if not isinstance(frame,list) or len(frame)!=4:raise ValueError("four frame bounds required")
    for x in frame:finite(x,"frame coordinate")
    if not .5<=frame[1]-frame[0]<=30 or not .5<=frame[3]-frame[2]<=30:raise ValueError("unsupported frame extent")
    sources=record.get("sources",[])
    if not 1<=len(sources)<=4:raise ValueError("1–4 sources required")
    ids=[s.get("id") for s in sources]
    if any(not isinstance(v,str) or not v for v in ids) or len(set(ids))!=len(ids):raise ValueError("unique source ids required")
    for source in sources:
        finite(source.get("x_deg"),"source x");finite(source.get("y_deg"),"source y")
    q=record.get("questionnaire",{})
    if q.get("presence") not in ("present","absent","uncertain"):raise ValueError("explicit presence report required")
    for key in ("visibility","bother","impact7d"):
        if q.get(key) is not None:finite(q[key],key,0.,10.)
    if q.get("confidence") is not None:finite(q["confidence"],"confidence",0.,100.)
    if not isinstance(q.get("free_text",""),str) or len(q.get("free_text",""))>2000:raise ValueError("free report must be text up to2000 characters")
    features=record.get("features",[])
    if not isinstance(features,list) or len(features)>8:raise ValueError("at most 8 features")
    for feature in features:validate_feature(feature,ids)
    if q["presence"]=="absent" and features:raise ValueError("absence must have no annotated features")
    if record["status"]=="complete":
        if record.get("consent_documented") is not True:raise ValueError("completed record requires documented consent")
        if q["presence"]=="present" and not features:raise ValueError("present report needs geometry or remain a draft")
        if any(f.get("confirmed") is not True for f in features):raise ValueError("unconfirmed geometry cannot be complete")
    return record


def eligibility(record):
    reasons=[]
    if record.get("status")!="complete":reasons.append("draft")
    if record.get("consent_documented") is not True:reasons.append("consent_not_documented")
    if record.get("method")!="external_standard_light":reasons.append("exploratory_observation_method")
    if record.get("questionnaire",{}).get("presence")=="uncertain":reasons.append("uncertain_presence")
    stimulus=record.get("stimulus",{})
    if not stimulus.get("protocol_id"):reasons.append("missing_stimulus_protocol")
    if stimulus.get("geometry_verified") is not True:reasons.append("source_geometry_unverified")
    if stimulus.get("photometry_verified") is not True:reasons.append("source_photometry_unverified")
    try:
        for key,lo,hi in (("source_luminance_cd_m2",.000001,1e9),("ambient_lux",0,1e6),("adaptation_minutes",0,180),
                          ("source_distance_mm",1,1e9),("source_diameter_mm",.000001,1e6)):
            finite(stimulus.get(key),key,lo,hi)
    except ValueError:reasons.append("stimulus_measurements_missing")
    calibration=record.get("calibration",{})
    if calibration.get("status")!="verified":reasons.append("angular_geometry_unverified")
    else:
        try:
            finite(calibration.get("ruler_css_px"),"ruler pixels",20,2000)
            finite(calibration.get("ruler_mm"),"ruler millimetres",5,1000)
            finite(calibration.get("view_distance_mm"),"viewing distance",100,3000)
        except ValueError:reasons.append("calibration_values_missing")
    for key in ("frame_clipped","source_localization_uncertain","calibration_changed"):
        if record.get("flags",{}).get(key) is True:reasons.append(key)
    xmin,xmax,ymin,ymax=record["frame_deg"]
    for feature in record["features"]:
        points=contour(feature,1024);half=feature["width_deg"]/2
        if np.any(points[:,0]-half<xmin) or np.any(points[:,0]+half>xmax) or np.any(points[:,1]-half<ymin) or np.any(points[:,1]+half>ymax):
            reasons.append("annotation_outside_frame");break
    if any(not (xmin<=s["x_deg"]<=xmax and ymin<=s["y_deg"]<=ymax) for s in record["sources"]):reasons.append("source_outside_frame")
    return reasons


def contour(feature,count=128):
    if feature["kind"]=="freehand":
        points=np.asarray(feature["points_deg"],dtype=float)
        length=np.r_[0.,np.cumsum(np.linalg.norm(np.diff(points,axis=0),axis=1))]
        keep=np.r_[True,np.diff(length)>1e-12];points=points[keep];length=length[keep]
        if length[-1]<=1e-12:return np.repeat(points[:1],count,axis=0)
        samples=np.linspace(0,length[-1],count)
        return np.column_stack([np.interp(samples,length,points[:,axis]) for axis in (0,1)])
    kind=feature["kind"]
    start=feature.get("start_deg",0.) if kind=="arc" else 0.
    sweep=feature.get("sweep_deg",360.) if kind=="arc" else 360.
    theta=np.radians(np.linspace(start,start+sweep,count,endpoint=kind=="arc"))
    rx=feature["rx_deg"];ry=feature.get("ry_deg",rx) if kind=="ellipse" else rx
    x=rx*np.cos(theta);y=ry*np.sin(theta)
    rot=math.radians(feature.get("rotation_deg",0.))
    return np.column_stack((feature["cx_deg"]+x*np.cos(rot)-y*np.sin(rot),feature["cy_deg"]+x*np.sin(rot)+y*np.cos(rot)))


@lru_cache(maxsize=256)
def _mask_cached(feature_json,frame,n):
    feature=json.loads(feature_json)
    xmin,xmax,ymin,ymax=frame
    x,y=np.meshgrid(np.linspace(xmin,xmax,n,endpoint=False)+(xmax-xmin)/n/2,np.linspace(ymin,ymax,n,endpoint=False)+(ymax-ymin)/n/2)
    if feature["kind"] in ("circle","arc"):
        xx=x-feature["cx_deg"];yy=y-feature["cy_deg"]
        distance=np.abs(np.hypot(xx,yy)-feature["rx_deg"])
        if feature["kind"]=="arc" and feature["sweep_deg"]<360:
            theta=(np.degrees(np.arctan2(yy,xx))-feature.get("start_deg",0.)-feature.get("rotation_deg",0.))%360
            valid=theta<=feature["sweep_deg"]
            endpoints=contour(feature,2)
            caps=np.minimum(np.hypot(x-endpoints[0,0],y-endpoints[0,1]),np.hypot(x-endpoints[1,0],y-endpoints[1,1]))
            distance=np.where(valid,distance,caps)
    else:
        points=np.asarray(feature["points_deg"],dtype=float) if feature["kind"]=="freehand" else contour(feature)
        if feature["kind"]=="ellipse":points=np.vstack((points,points[0]))
        distance=np.full_like(x,np.inf)
        for a,b in zip(points[:-1],points[1:]):
            delta=b-a;length=float(delta@delta)
            t=np.clip(((x-a[0])*delta[0]+(y-a[1])*delta[1])/max(length,1e-20),0,1)
            distance=np.minimum(distance,np.hypot(x-a[0]-t*delta[0],y-a[1]-t*delta[1]))
    return distance<=feature["width_deg"]/2


def mask(feature,frame,n):
    geometry={key:value for key,value in feature.items() if key not in ("id","confirmed","brightness8","confidence")}
    return _mask_cached(json.dumps(geometry,sort_keys=True),tuple(frame),n)


def iou(a,b):
    denominator=np.count_nonzero(a|b)
    return np.count_nonzero(a&b)/denominator if denominator else 1.


def optimal_pairs(candidates,references,cmasks,rmasks):
    # <=8 features: exact assignment, unmatched candidate/reference each cost1.
    @lru_cache(None)
    def solve(index,used):
        if index==len(candidates):return float(len(references)-used.bit_count()),()
        tail,pairs=solve(index+1,used);best=(1.+tail,pairs)
        for j,reference in enumerate(references):
            if used&(1<<j) or candidates[index]["source_id"]!=reference["source_id"]:continue
            tail,pairs=solve(index+1,used|(1<<j))
            proposal=(1.-iou(cmasks[index],rmasks[j])+tail,((index,j),)+pairs)
            if proposal[0]<best[0]:best=proposal
        return best
    return list(solve(0,0)[1])


def compare_report(prediction,record,n=256):
    if not isinstance(prediction,dict):raise ValueError("prediction must be an object")
    if not isinstance(prediction.get("symptoms",{}),dict):raise ValueError("symptoms must be an object when supplied")
    sources={s["id"] for s in record["sources"]}
    p=finite(prediction.get("presence_probability"),"presence probability",0.,1.)
    proposed=prediction.get("features",[])
    if not isinstance(proposed,list) or len(proposed)>8:raise ValueError("at most8 predicted features")
    for f in proposed:validate_feature(f,sources)
    candidates=proposed if p>=.5 else []
    references=record["features"]
    frame=record["frame_deg"]
    cm=[mask(f,frame,n) for f in candidates];rm=[mask(f,frame,n) for f in references]
    if any(not m.any() for m in cm+rm):raise ValueError("feature has no pixels in evaluation frame; review clipping/resolution")
    c_union=np.logical_or.reduce(cm) if cm else np.zeros((n,n),bool)
    r_union=np.logical_or.reduce(rm) if rm else np.zeros((n,n),bool)
    spatial=iou(c_union,r_union)
    pairs=optimal_pairs(candidates,references,cm,rm)
    count_f1=2*len(pairs)/(len(candidates)+len(references)) if candidates or references else 1.
    distances=[];profile=[];profile_possible=0;intervals=[];parameters=[]
    for ci,ri in pairs:
        c,r=candidates[ci],references[ri]
        a,b=contour(c),contour(r)
        ds=np.linalg.norm(a[:,None,:]-b[None,:,:],axis=2)
        distances.append(.5*(float(ds.min(axis=0).mean())+float(ds.min(axis=1).mean())))
        if c["kind"]!="freehand" and r["kind"]!="freehand":
            parameters.append({"center_error_deg":math.hypot(c["cx_deg"]-r["cx_deg"],c["cy_deg"]-r["cy_deg"]),
                               "relative_source_offset_error_deg":math.hypot(c["cx_deg"]-r["cx_deg"],c["cy_deg"]-r["cy_deg"]),
                               "radius_x_error_deg":abs(c["rx_deg"]-r["rx_deg"]),
                               "radius_y_error_deg":abs(c.get("ry_deg",c["rx_deg"])-r.get("ry_deg",r["rx_deg"])),
                               "reported_width_error_deg":abs(c["width_deg"]-r["width_deg"])})
        rb=r.get("brightness8");cb=c.get("brightness8")
        if rb is not None:
            profile_possible+=sum(v is not None for v in rb)
            if cb is not None:profile.extend(abs(x-y)/4 for x,y in zip(cb,rb) if x is not None and y is not None)
        interval=c.get("radius_interval_deg")
        if interval is not None and c["kind"]==r["kind"]=="circle":
            if not isinstance(interval,list) or len(interval)!=2:raise ValueError("radius interval needs two bounds")
            lo=finite(interval[0],"interval lower",0,30);hi=finite(interval[1],"interval upper",0,30)
            if not lo<=c["rx_deg"]<=hi:raise ValueError("radius point estimate must lie inside interval")
            value=r["rx_deg"];alpha=.2
            interval_score=hi-lo+2/alpha*max(lo-value,0)+2/alpha*max(value-hi,0)
            intervals.append({"covered":lo<=value<=hi,"width_deg":hi-lo,"interval_score_deg":interval_score,"nominal_coverage":.8})
    contour_scores=[math.exp(-distance/.25) for distance in distances]
    denominator=max(len(candidates),len(references))
    contour_score=sum(contour_scores)/denominator if denominator else 1.
    y=1. if record["questionnaire"]["presence"]=="present" else 0.
    brier=(p-y)**2
    score=100*(.2*(1-brier)+.4*spatial+.2*contour_score+.2*count_f1)
    symptoms={}
    for key in ("visibility","bother","impact7d"):
        predicted=prediction.get("symptoms",{}).get(key);observed=record["questionnaire"].get(key)
        if predicted is not None:
            finite(predicted,key,0,10)
            if observed is not None:symptoms[key+"_absolute_error"]=abs(predicted-observed)
    return {"score":score,"presence_brier":brier,"region_iou":spatial,"contour_score":contour_score,
            "count_source_f1":count_f1,"matched_count":len(pairs),"candidate_count":len(candidates),"reference_count":len(references),
            "contour_errors_deg":distances,"structured_parameter_errors":parameters,
            "ordinal_brightness_mae":float(np.mean(profile)) if profile else None,
            "ordinal_profile_coverage":len(profile)/profile_possible if profile_possible else None,
            "radius_intervals":intervals,"symptom_errors":symptoms}


def evaluate_patient_reports(reference_file,prediction_file,grid_n=256):
    if type(grid_n) is not int or not 64<=grid_n<=512:raise ValueError("evaluation grid64–512 required")
    records=json.loads(Path(reference_file).read_text())["records"]
    predictions=json.loads(Path(prediction_file).read_text())["predictions"]
    if len({r.get("origin") for r in records})>1:raise ValueError("score synthetic and research records in separate runs")
    lookup={}
    for prediction in predictions:
        if not isinstance(prediction,dict):raise ValueError("prediction must be an object")
        key=prediction.get("target_id")
        if not isinstance(key,str) or key in lookup:raise ValueError("unique prediction target_ids required")
        lookup[key]=prediction
    groups={};excluded=[];seen=set()
    for record in records:
        validate_record(record)
        identity=(target_key(record),record["repeat_id"])
        if identity in seen:raise ValueError("duplicate repeat within target")
        seen.add(identity)
        reasons=eligibility(record)
        if not reasons and any(not mask(f,record["frame_deg"],grid_n).any() for f in record["features"]):
            reasons.append("reference_below_grid_resolution")
        if reasons:excluded.append({"record_id":record["record_id"],"reasons":reasons});continue
        groups.setdefault(target_key(record),[]).append(record)
    reports=[]
    for key,observations in groups.items():
        if len(observations)<2:
            excluded.extend({"record_id":r["record_id"],"reasons":["fewer_than_two_eligible_repeats"]} for r in observations);continue
        first=observations[0]
        if any(r["frame_deg"]!=first["frame_deg"] or r["sources"]!=first["sources"] or r["origin"]!=first["origin"] or r.get("stimulus")!=first.get("stimulus") for r in observations):
            excluded.extend({"record_id":r["record_id"],"reasons":["inconsistent_condition_metadata"]} for r in observations);continue
        prediction=lookup.get(key)
        rows=[];invalid=None
        if prediction is None:invalid="missing prediction"
        else:
            try:rows=[compare_report(prediction,r,grid_n) for r in observations]
            except (ValueError,KeyError,TypeError) as exc:invalid=str(exc)
        # Observed repeat agreement is contextual, not a score ceiling or weight.
        repeat_agreement=[]
        for i,a in enumerate(observations):
            for b in observations[i+1:]:
                q={"presence_probability":1. if a["questionnaire"]["presence"]=="present" else 0.,"features":a["features"]}
                repeat_agreement.append(compare_report(q,b,grid_n)["score"])
        reports.append({"target_id":key,"patient_id":first["patient_id"],"origin":first["origin"],"repeat_count":len(observations),
                        "score":float(np.mean([r["score"] for r in rows])) if not invalid else 0.,"invalid":invalid,
                        "repeat_agreement_context":float(np.mean(repeat_agreement)),"observations":rows})
    patients={}
    for report in reports:patients.setdefault(report["patient_id"],[]).append(report["score"])
    macro={key:float(np.mean(values)) for key,values in patients.items()}
    interval=None
    if len(macro)>=2:
        rng=np.random.default_rng(20261027);scores=np.asarray(list(macro.values()))
        draws=rng.choice(scores,size=(2000,len(scores)),replace=True).mean(axis=1)
        interval=np.quantile(draws,[.025,.975]).tolist()
    raw_intervals=[item for report in reports for row in report["observations"] for item in row["radius_intervals"]]
    uncertainty={"evaluated_count":len(raw_intervals),"nominal_coverage":.8,
                 "observed_repeat_coverage":float(np.mean([r["covered"] for r in raw_intervals])) if raw_intervals else None,
                 "mean_width_deg":float(np.mean([r["width_deg"] for r in raw_intervals])) if raw_intervals else None,
                 "mean_interval_score_deg":float(np.mean([r["interval_score_deg"] for r in raw_intervals])) if raw_intervals else None}
    all_ids={target_key(r) for r in records}
    return {"schema_version":VERSION,"origin":records[0]["origin"] if records else None,"patient_agreement_development_score":float(np.mean(list(macro.values()))) if macro else None,
            "clinical_score":None,"mechanistic_understanding_score":None,"validated_instrument":False,
            "patient_count":len(macro),"target_count":len(reports),"eligible_record_count":sum(r["repeat_count"] for r in reports),
            "excluded_records":excluded,"invalid_target_count":sum(r["invalid"] is not None for r in reports),
            "unknown_prediction_targets":sorted(set(lookup)-all_ids),"patient_macro_scores":macro,
            "patient_cluster_bootstrap_95_interval":interval,"interval_prediction_diagnostics":uncertainty,
            "grid_n":grid_n,"targets":reports,
            "notes":["Draft development weights/scales, not clinical acceptance thresholds.",
                     "Repeated reports are not independent patients. Agreement is not proof of optical causality.",
                     "Recalls, unverified geometry/photometry, screen matching, uncertain/censored reports are excluded from primary scores, never relabeled as absence.",
                     "Research status and calibration confirmations are investigator attestations; the browser/scorer cannot verify apparatus or consent."]}