import sys, json, math
sys.path += ['/root/.local/lib/python3.13/site-packages','/usr/local/lib/python3.13/dist-packages','/usr/lib/python3/dist-packages']
import numpy as np
from scipy.spatial import cKDTree
SEC=json.load(open('/tmp/verify_daesan2_run/sections_final.json'))
FIT=json.load(open('/tmp/verify_daesan2_run/fits.json'))
C0=(14126.3,15305.0)
def pts(st, pred):
    out=[]
    for sid,v in SEC[st]['solids'].items():
        if pred(v['group'], v['id']):
            for poly in v['polys']: out.append(np.array(poly)[:,:2])
    return np.vstack(out)
def arc(cx,cy,r,a0,a1,step=1.0):
    if a1<a0: a1+=360
    n=max(2,int(r*math.radians(a1-a0)/step)); th=np.radians(np.linspace(a0,a1,n))
    return np.c_[cx+r*np.cos(th), cy+r*np.sin(th)]
def seg(x0,y0,x1,y1,step=1.0):
    n=max(2,int(math.hypot(x1-x0,y1-y0)/step)); t=np.linspace(0,1,n); return np.c_[x0+(x1-x0)*t,y0+(y1-y0)*t]
def circ(cx,cy,r,step=1.0):
    n=max(8,int(2*math.pi*r/step)); th=np.linspace(0,2*math.pi,n,endpoint=False); return np.c_[cx+r*np.cos(th),cy+r*np.sin(th)]
REF={
 'LIN_crown_outer': arc(*C0,7282,40,140), 'LIN_crown_inner': arc(*C0,6982,40,140),
 'LIN_haunchL_outer': arc(11959.9,17122.8,4454,140,204.67), 'LIN_haunchL_inner': arc(11621.3,17406.9,3712,140,195.56),
 'LIN_haunchR_outer': arc(16105.0,16965.3,4699,336.35,40), 'LIN_haunchR_inner': arc(16479.6,17279.7,3910,344.38,40),
 'LIN_footL': np.vstack([seg(8171.3,16410.9,8471.3,16410.9), seg(8471.3,16410.9,8621.3,16410.9), seg(8621.3,16410.9,8651.3,16000.9), seg(8651.3,16000.9,8771.3,15825.9), seg(8771.3,15825.9,8771.3,15750.9)]),
 'LIN_footR': np.vstack([seg(20181.3,16226.9,19781.3,16226.9), seg(19781.3,16226.9,19631.3,16226.9), seg(19631.3,16226.9,19601.3,15816.9), seg(19601.3,15816.9,19481.3,15641.9), seg(19481.3,15641.9,19481.3,15566.9)]),
 'EXC_crown': arc(*C0,7402,40,140), 'EXC_haunchL': arc(11959.9,17122.8,4574,140,203.98), 'EXC_haunchR': arc(16105.0,16965.3,4819,336.97,40),
 'PAV_top': seg(9526.3,15764.0,18726.3,15580.0), 'PAV_bot': seg(9526.3,15484.0,18726.3,15300.0), 'PAV_sideL': seg(9526.3,15764.0,9526.3,15484.0), 'PAV_sideR': seg(18726.3,15580.0,18726.3,15300.0),
 'BASE_top': seg(9526.3,15484.0,18726.3,15300.0), 'BASE_bot': seg(9526.3,15334.0,18726.3,15150.0),
 'DRN_wallL1': seg(9101.3,14864.0,9101.3,15695.5), 'DRN_wallL2': seg(9401.3,14864.0,9401.3,15701.5), 'DRN_botL': seg(9101.3,14864.0,9401.3,14864.0), 'DRN_pipeL': circ(9251.3,15064.0,150),
 'DRN_wallR1': seg(18851.3,14680.0,18851.3,15517.5), 'DRN_wallR2': seg(19151.3,14680.0,19151.3,15511.5), 'DRN_botR': seg(18851.3,14680.0,19151.3,14680.0), 'DRN_pipeR': circ(19001.3,14880.0,150),
}
def to_dxf(P, st, sig):
    uc,vc,_=FIT[st]['outer']
    return np.c_[C0[0]+sig*(P[:,0]-uc)*1000.0, C0[1]+(P[:,1]-vc)*1000.0]
def dist_stats(M, R, label):
    """M model pts (dxf frame), R reference pts. Returns d(ref->model) max, d(model->ref) max within band."""
    if len(M)==0 or len(R)==0: return None
    tM=cKDTree(M); tR=cKDTree(R)
    dR,_=tM.query(R); dM,_=tR.query(M)
    return dict(ref_to_model_max=float(dR.max()), ref_to_model_mean=float(dR.mean()), model_to_ref_max=float(dM.max()), model_to_ref_mean=float(dM.mean()), n_ref=len(R), n_model=len(M))
def run(sig, st, feats):
    res={}
    L=to_dxf(pts(st,lambda g,i:g=='라이닝'),st,sig); E=to_dxf(pts(st,lambda g,i:g=='굴착면'),st,sig)
    P=to_dxf(pts(st,lambda g,i:g=='노면(계획고)'),st,sig); B=to_dxf(pts(st,lambda g,i:g=='기층·인버트(가정)'),st,sig)
    D=to_dxf(pts(st,lambda g,i:g=='배수(측구)'),st,sig)
    def side_clip(R, M, side):
        # keep reference points on given side (L: x<C0x, R: x>C0x) and above model's lowest y on that side minus 50mm
        sm = M[(M[:,0]<C0[0])] if side=='L' else M[(M[:,0]>C0[0])]
        sr = R[(R[:,0]<C0[0])] if side=='L' else R[(R[:,0]>C0[0])]
        if len(sm)==0: return sr, 0, None
        ymin=sm[:,1].min()-0.05*1000
        keep=sr[sr[:,1]>=ymin]
        dropped=len(sr)-len(keep)
        rng=None if dropped==0 else [round(float(sr[sr[:,1]<ymin,1].min())),round(float(sr[sr[:,1]<ymin,1].max()))]
        return keep, dropped, rng
    out={}
    # Lining: arcs and foot lines; clip below model extent per side
    lin_ref=np.vstack([REF[k] for k in ['LIN_crown_outer','LIN_crown_inner','LIN_haunchL_outer','LIN_haunchL_inner','LIN_haunchR_outer','LIN_haunchR_inner']])
    keepL,dL,rL=side_clip(lin_ref[lin_ref[:,0]<C0[0]],L,'L'); keepR,dR_,rR=side_clip(lin_ref[lin_ref[:,0]>C0[0]],L,'R')
    crown=REF['LIN_crown_outer'];
    out['LIN_arcs']=dist_stats(L, np.vstack([crown, REF['LIN_crown_inner'], keepL, keepR]), 'LIN')
    out['LIN_arcs']['clipped_below_model_extent']={'L':dL,'R':dR_,'yrange_L':rL,'yrange_R':rR}
    out['LIN_footlines']=dist_stats(L, np.vstack([REF['LIN_footL'],REF['LIN_footR']]), 'LIN_foot')
    # per-arc max deviations (ref->model), by arc name
    perarc={}
    for k in ['LIN_crown_outer','LIN_crown_inner','LIN_haunchL_outer','LIN_haunchL_inner','LIN_haunchR_outer','LIN_haunchR_inner']:
        R=REF[k]; tM=cKDTree(L); d,_=tM.query(R); perarc[k]=round(float(d.max())*1,1)
    out['LIN_perarc_ref2model_max_mm']=perarc
    # Excavation
    exc_ref=np.vstack([REF['EXC_crown'],REF['EXC_haunchL'],REF['EXC_haunchR']])
    keepE_L,eL,reL=side_clip(exc_ref,E,'L'); keepE_R,eR,reR=side_clip(exc_ref,E,'R')
    out['EXC_arcs']=dist_stats(E, np.vstack([REF['EXC_crown'],keepE_L,keepE_R]),'EXC')
    out['EXC_arcs']['clipped_below_model_extent']={'L':eL,'R':eR,'yrange_L':reL,'yrange_R':reR}
    tE=cKDTree(E); d,_=tE.query(REF['EXC_crown']); out['EXC_crown_ref2model_max_mm']=round(float(d.max()),1)
    # Pavement: top/bottom/sides
    pav_ref=np.vstack([REF['PAV_top'],REF['PAV_bot'],REF['PAV_sideL'],REF['PAV_sideR']])
    out['PAV']=dist_stats(P[(P[:,0]>9526.3-150)&(P[:,0]<18726.3+150)], pav_ref,'PAV')
    # Base: model points restricted to road width (drawn span)
    Bw=B[(B[:,0]>=9526.3-1)&(B[:,0]<=18726.3+1)]
    out['BASE_drawn_span']=dist_stats(Bw, np.vstack([REF['BASE_top'],REF['BASE_bot']]),'BASE')
    out['BASE_model_extent_x']=[round(float(B[:,0].min())),round(float(B[:,0].max()))]
    out['BASE_model_ymin_ymax']=[round(float(B[:,1].min())),round(float(B[:,1].max()))]
    # Drains: both ditches vs both references
    dref=np.vstack([REF[k] for k in ['DRN_wallL1','DRN_wallL2','DRN_botL','DRN_pipeL','DRN_wallR1','DRN_wallR2','DRN_botR','DRN_pipeR']])
    dm=D[(D[:,0]>8900)&(D[:,0]<19300)]
    out['DRN']=dist_stats(dm, dref,'DRN')
    # per-ditch: ref walls/bottom/pipe to model
    tD=cKDTree(dm)
    out['DRN_perref_ref2model_max_mm']={k: round(float(tD.query(REF[k])[0].max()),1) for k in ['DRN_wallL1','DRN_wallL2','DRN_botL','DRN_pipeL','DRN_wallR1','DRN_wallR2','DRN_botR','DRN_pipeR']}
    out['DRN_model_xrange']=[round(float(dm[:,0].min())),round(float(dm[:,0].max()))]
    out['DRN_model_yrange']=[round(float(dm[:,1].min())),round(float(dm[:,1].max()))]
    return out
allres={}
for sig in [1,-1]:
    for st in ['4745.0005','5100.0','5469.9995']:
        r=run(sig,st,None); allres[f"{sig}|{st}"]=r
        print(f"=== sigma={sig} st={st}  slope={FIT[st]['slope_top']*sig:.5f}")
        for k in ['LIN_arcs','EXC_arcs','PAV','BASE_drawn_span','DRN']:
            v=r[k]; 
            if v is None: print('  ',k,'none'); continue
            print(f"  {k:16s} ref->model max={v['ref_to_model_max']*1000:7.1f}mm mean={v['ref_to_model_mean']*1000:6.1f} | model->ref max={v['model_to_ref_max']*1000:7.1f}mm mean={v['model_to_ref_mean']*1000:6.1f}")
            if 'clipped_below_model_extent' in v: print('     clipped', v['clipped_below_model_extent'])
        print('   LIN foot lines model->ref / ref->model:', None if r['LIN_footlines'] is None else (round(r['LIN_footlines']['model_to_ref_max']*1000,1), round(r['LIN_footlines']['ref_to_model_max']*1000,1)))
        print('   LIN per-arc ref->model max mm:', r['LIN_perarc_ref2model_max_mm'])
        print('   EXC crown ref->model max mm:', r['EXC_crown_ref2model_max_mm'])
        print('   BASE model x-range', r['BASE_model_extent_x'], 'y-range', r['BASE_model_ymin_ymax'])
        print('   DRN per-ref ref->model max mm:', r['DRN_perref_ref2model_max_mm'], 'model x', r['DRN_model_xrange'], 'y', r['DRN_model_yrange'])
json.dump(allres, open('/tmp/verify_daesan2_run/hd2_results.json','w'), indent=1, default=str)