import sys
exec(open('/tmp/verify_daesan2_run/scripts/hd2.py').read().split('allres={}')[0])
from scipy.spatial import cKDTree
import numpy as np
extra={
 'PAV_sideL': seg(9526.3,15764.0,9526.3,15484.0), 'PAV_sideR': seg(18726.3,15580.0,18726.3,15300.0),
 'DRN_topL': seg(9026.3,15694.0,9476.3,15703.0), 'DRN_topR': seg(18776.3,15519.0,19226.3,15510.0),
}
REF.update(extra)
def dd(A,B):
    d1,_=cKDTree(A).query(B); d2,_=cKDTree(B).query(A); return d1,d2
for st in ['4745.0005','5100.0','5469.9995']:
    print("=== station", st, "sigma=+1 slope(dxf)", round(FIT[st]['slope_top'],5))
    L=to_dxf(pts(st,lambda g,i:g=='라이닝'),st,1); E=to_dxf(pts(st,lambda g,i:g=='굴착면'),st,1)
    P=to_dxf(pts(st,lambda g,i:g=='노면(계획고)'),st,1); B=to_dxf(pts(st,lambda g,i:g=='기층·인버트(가정)'),st,1)
    D=to_dxf(pts(st,lambda g,i:g=='배수(측구)'),st,1)
    # LIN full arcs (no clipping) and clipped to model extent
    arcs=['LIN_crown_outer','LIN_crown_inner','LIN_haunchL_outer','LIN_haunchL_inner','LIN_haunchR_outer','LIN_haunchR_inner']
    for k in arcs:
        d,_=cKDTree(L).query(REF[k]); print(f"  LIN {k:18s} max={d.max():7.1f} mm  (at y_min_ref={REF[k][:,1].min():.0f})")
    ymin_L=L[L[:,0]<C0[0]][:,1].min(); ymin_R=L[L[:,0]>C0[0]][:,1].min()
    lr=np.vstack([REF[k] for k in arcs]); d,_=cKDTree(L).query(lr)
    inside=((lr[:,0]<C0[0])&(lr[:,1]>=ymin_L))|((lr[:,0]>=C0[0])&(lr[:,1]>=ymin_R))
    print(f"  LIN within model extent (y>={ymin_L:.0f}/{ymin_R:.0f}) max={d[inside].max():.1f} mm; drawn part below model extent (not modelled) n={int((~inside).sum())} pts, y range {lr[~inside,1].min():.0f}-{lr[~inside,1].max():.0f}")
    fl=np.vstack([REF['LIN_footL'],REF['LIN_footR']]); a,b=dd(L,fl); print(f"  LIN foot lines (CS-CONC-MAJR wall/foot) ref->model max={a.max():.1f} model->ref(near band) ")
    # EXC
    ea=['EXC_crown','EXC_haunchL','EXC_haunchR']
    eymL=E[E[:,0]<C0[0]][:,1].min(); eymR=E[E[:,0]>C0[0]][:,1].min()
    er=np.vstack([REF[k] for k in ea]); d,_=cKDTree(E).query(er)
    ins=((er[:,0]<C0[0])&(er[:,1]>=eymL))|((er[:,0]>=C0[0])&(er[:,1]>=eymR))
    print(f"  EXC arcs within model extent (y>={eymL:.0f}/{eymR:.0f}) max={d[ins].max():.1f} mm; crown max={cKDTree(E).query(REF['EXC_crown'])[0].max():.1f}; below model extent n={int((~ins).sum())} y {er[~ins,1].min():.0f}-{er[~ins,1].max():.0f}")
    # PAV model vs drawn top/bottom/sides
    pref=np.vstack([REF['PAV_top'],REF['PAV_bot'],REF['PAV_sideL'],REF['PAV_sideR']])
    Pm=P[(P[:,0]>9526.3-150)&(P[:,0]<18726.3+150)]
    a,b=dd(Pm,pref); print(f"  PAV ref->model max={a.max():.1f} mm ; model->ref max={b.max():.1f} mm (top/bottom/sides)")
    # BASE: top envelope vs base top line; bottom envelope vs base bottom line, within road span
    bw=B[(B[:,0]>=9526.3)&(B[:,0]<=18726.3)]
    xs=np.linspace(9526.3,18726.3,921)
    top_line=15484.0+(xs-9526.3)*(15300-15484)/9200.0; bot_line=15334.0+(xs-9526.3)*(15150-15334)/9200.0
    # model envelope per x-bin (1 m bins for robust): max y (top) and min y (bottom)
    tops=[];bots=[]
    for x0 in np.arange(9526.3,18726.3,200):
        m=(bw[:,0]>=x0)&(bw[:,0]<x0+200)
        if m.sum()==0: continue
        xm=bw[m,0].mean(); tops.append((xm,bw[m,1].max())); bots.append((xm,bw[m,1].min()))
    tops=np.array(tops); bots=np.array(bots)
    et=np.abs(tops[:,1]-np.interp(tops[:,0],xs,top_line)); eb=np.abs(bots[:,1]-np.interp(bots[:,0],xs,bot_line))
    print(f"  BASE top envelope vs drawn top (15484->15300): max={et.max():.1f} mm, mean={et.mean():.1f}; bottom envelope vs drawn bottom (15334->15150): max={eb.max():.1f} mm, mean={eb.mean():.1f}")
    print(f"  BASE model bottom envelope at span ends: x={bots[0,0]:.0f} y={bots[0,1]:.1f} / x={bots[-1,0]:.0f} y={bots[-1,1]:.1f}; top envelope ends y={tops[0,1]:.1f}/{tops[-1,1]:.1f}")
    print(f"  BASE model extent outside road span: x {B[:,0].min():.0f}..{B[:,0].max():.0f}, y {B[:,1].min():.0f}..{B[:,1].max():.0f}")
    # DRN
    dref=np.vstack([REF[k] for k in ['DRN_wallL1','DRN_wallL2','DRN_botL','DRN_pipeL','DRN_wallR1','DRN_wallR2','DRN_botR','DRN_pipeR','DRN_topL','DRN_topR']])
    dm=D[(D[:,0]>8900)&(D[:,0]<19300)]
    a,b=dd(dm,dref); print(f"  DRN ref->model max={a.max():.1f} mm ; model->ref max={b.max():.1f} mm (walls, bottoms, pipes, tops)")
    # vertical datum check at tunnel centre and road centre
    fit=FIT[st]; uc,vc,_=fit['outer']
    top_c=fit['slope_top']*(0-uc)+0 ; 
    # top envelope at u for DXF tunnel centre (u=uc) and road centre (x=9776.3 -> u=uc-4.35)
    uu=P[:,0]; 
    def top_env(xdxf):
        um=uc+(xdxf-C0[0])/1000.0
        m=np.abs(pts(st,lambda g,i:g=='노면(계획고)')[:,0]-um)<0.05
        return pts(st,lambda g,i:g=='노면(계획고)')[m,1].max()
    Pm_raw=pts(st,lambda g,i:g=='노면(계획고)')
    for lab,xdxf,ydxf in [('tunnel centre x=14126.3',14126.3,15672.0),('road centre x=9776.3',9776.3,15759.0)]:
        um=uc+(xdxf-C0[0])/1000.0
        m=np.abs(Pm_raw[:,0]-um)<0.05
        yt=(Pm_raw[m,1].max()-vc)*1000; yb=(Pm_raw[m,1].min()-vc)*1000
        print(f"  datum {lab}: model top={yt:.1f} mm, bottom={yb:.1f} mm above lining centre | DXF top at this x = {ydxf-15305:.1f} mm above lining centre")