import struct, os

d = r"C:\Project\MAG160C\analysis\pairs_recheck_20260811"
npix = 19200
nsegs = 3

def load16(name, n):
    b = open(os.path.join(d, name), "rb").read()
    return struct.unpack_from("<%dH" % n, b, 0)

def loadi16(name, n):
    b = open(os.path.join(d, name), "rb").read()
    return struct.unpack_from("<%dh" % n, b, 0)

f20 = load16("pair_000.f20", npix)
ref = load16("pair_000.ref", npix)
raw = load16("pair_000.raw", npix)

# working tables (from build-artifacts, verified = interp of EP2/EP3)
gain = load16(r"C:\Project\MAG160C\build-artifacts\mag160c_official_nuc_gain.bin", npix*nsegs*2)
thr = loadi16(r"C:\Project\MAG160C\build-artifacts\mag160c_official_nuc_thr.bin", npix*(nsegs-1))

# blind records from DDT (all endpoints identical; use EP2)
ddt = open(r"C:\Project\MAG160C\analysis\ida\Core160043865.ddt", "rb").read()
thr_ep0 = 0x80
block = (nsegs-1)*npix*2 + 1024 + nsegs*npix*4
blind_off = thr_ep0 + 6*block
recs = []
for i in range(32):
    r = blind_off + i*40
    target, typ = struct.unpack_from("<II", ddt, r)
    neigh = struct.unpack_from("<8I", ddt, r+8)
    recs.append((target, typ, neigh))
print("blind recs:", len(recs), "first:", recs[0])

def blind_apply(buf, recs):
    out = list(buf)
    for (t, typ, ng) in recs:
        if typ == 3:
            v = out[ng[0]] + out[ng[1]] + out[ng[2]]
            out[t] = v // 3
        elif typ == 4:
            v = out[ng[0]] + out[ng[1]] + out[ng[2]] + out[ng[3]]
            out[t] = v >> 2
        elif typ == 5:
            v = out[ng[0]] + out[ng[1]] + out[ng[2]] + out[ng[3]] + out[ng[4]]
            out[t] = v // 5
        elif typ == 6:
            v = sum(out[ng[j]] for j in range(6))
            out[t] = v // 6
        elif typ == 7:
            v = sum(out[ng[j]] for j in range(7))
            out[t] = v // 7
        elif typ == 8:
            v = sum(out[ng[j]] for j in range(8))
            out[t] = v >> 3
    return out

# NUC with exact official math
nuc = [0]*npix
for i in range(npix):
    d2 = (int(f20[i]) - int(ref[i])) >> 1
    seg = 0
    if d2 > thr[i*2]:
        seg = 1
        if d2 > thr[i*2+1]:
            seg = 2
    p = (seg*npix + i)*2
    v = (gain[p] * d2) >> 12 + gain[p+1]  # CAREFUL: precedence! official: off + ((gain*d2)>>12)
    v = gain[p+1] + ((gain[p] * d2) >> 12)
    if v < 0: v = 0
    if v > 65535: v = 65535
    nuc[i] = v

# stats before blind
diff = sum(1 for i in range(npix) if nuc[i] != raw[i])
sad = sum(abs(nuc[i]-raw[i]) for i in range(npix))
print(f"NUC only: diff={diff}/{npix} MAE={sad/npix:.4f}")

# apply blind comp
out = blind_apply(nuc, recs)
diff2 = sum(1 for i in range(npix) if out[i] != raw[i])
sad2 = sum(abs(out[i]-raw[i]) for i in range(npix))
print(f"NUC+blind: diff={diff2}/{npix} MAE={sad2/npix:.4f}")

# list mismatch pixels
mism = [(i, nuc[i], out[i], raw[i]) for i in range(npix) if out[i] != raw[i]]
print("mismatches:", len(mism))
for m in mism[:15]:
    print("  pix", m)

# blind target pixels check
for (t, typ, ng) in recs:
    print(f"  rec target={t} type={typ} neigh={list(ng[:typ-2])}")