"""
PHOENIX research — 12-channel candidate precision, POWERED + CONFOUND-CONTROLLED (v2).

v1 finding (kept, it is the headline): a contextual MIR-TIR front-end emits 48,534 candidates
over the 2,591-frame corpus and only 10 are real FIRMS fires (precision 0.0002) — the SEVIRI floor
means real fires (767 FIRMS-matched pixels) almost never produce a detectable thermal anomaly
(10/767 = 1.3% enter the candidate set). With n_real=10 the apparent spectral AUC (0.95) is a
small-sample mirage; the placebo scored 0.65 (not 0.5) = leakage. So v1 is not a usable answer.

v2 asks the question fairly and with power: among pixels at the SAME thermal-anomaly level, does the
FULL raw 11-band spectral vector separate real fire from non-fire BEYOND the derived intensity
features (entry 0037)? Positives = all 767 FIRMS-matched fire pixels. Negatives = non-fire pixels
sampled to MATCH each positive's (MIR, MIR-TIR) so intensity cannot be a shortcut. Honest spatial
GroupKFold CV + placebo (shuffled labels must be ~0.5).
"""
import numpy as np, glob, os, re, sqlite3, time, json
from datetime import datetime, timezone
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import GroupKFold
from sklearn.metrics import roc_auc_score

FRAMEDIR="/workspace/seviri_frames"; DB="/workspace/data/ground_truth.sqlite"
OUT="/workspace/spec12_v2_results.json"; FEAT_OUT="/workspace/spec12_v2_features.npz"
def log(m): print("[v2 %s] %s"%(time.strftime("%H:%M:%S"),m),flush=True)
CH=['bt_ir039','bt_wv062','bt_wv073','bt_ir087','bt_ir097','bt_ir108','bt_ir120','bt_ir134',
    'bt_vis006','bt_vis008','bt_nir016']
MIR='bt_ir039'; TIR='bt_ir108'
FIRMS={'firms_viirs_noaa20','firms_viirs_noaa21','firms_viirs_snpp','firms_modis_nrt',
       'slstr_frp_s3a','slstr_frp_s3b','mtg_af_l2','viirs_dnb_fire',
       'firms_viirs_noaa20_nrt_archive','firms_viirs_noaa21_nrt_archive','firms_viirs_snpp_nrt_archive',
       'firms_modis_nrt_archive','firms_viirs_noaa20_global'}
FN_TS=re.compile(r'-NA-(\d{14})')
def fts(n):
    m=FN_TS.search(n); return datetime.strptime(m.group(1),'%Y%m%d%H%M%S').replace(tzinfo=timezone.utc) if m else None
def pts(ts):
    if not ts: return None
    s=ts.strip()
    for f in ('%Y-%m-%dT%H:%M:%S','%a, %d %b %Y %H:%M:%S'):
        try: return datetime.strptime(s.replace('Z','').split('.')[0].split('+')[0].strip(),f).replace(tzinfo=timezone.utc)
        except: pass
    return None
files=sorted(glob.glob(os.path.join(FRAMEDIR,'*.npz')))
fr_,ft_,fla_,flo_=[],[],[],[]
for p in files:
    t=fts(os.path.basename(p));
    if t is None: continue
    z=np.load(p,allow_pickle=True)
    if MIR not in z or TIR not in z: continue
    try: cube=np.stack([np.nan_to_num(z[c].astype(np.float32),nan=0.0) for c in CH],0)
    except: continue
    if cube.shape[1:]!=(57,100): continue
    fr_.append(cube);ft_.append(t);fla_.append(z['lat'].astype(np.float32));flo_.append(z['lon'].astype(np.float32))
F=np.stack(fr_,0);LAT=np.stack(fla_,0);LON=np.stack(flo_,0);T,C,H,W=F.shape
o=np.argsort([t.timestamp() for t in ft_]);F=F[o];LAT=LAT[o];LON=LON[o];ft_=[ft_[i] for i in o]
mir=CH.index(MIR);tir=CH.index(TIR)
log("frames=%d cube=%s"%(T,(C,H,W)))
con=sqlite3.connect(DB);cur=con.cursor()
fires=[]
for (s,la,lo,ts,frp) in cur.execute("SELECT source,lat,lng,ts,frp_w FROM external_fires"):
    if s not in FIRMS: continue
    t=pts(ts)
    if t and la is not None and lo is not None: fires.append((la,lo,t.timestamp()))
fa=np.array([t.timestamp() for t in ft_])
Y=np.zeros((T,H,W),np.uint8)
for (la,lo,tt) in fires:
    dt=np.abs(fa-tt); fi=int(dt.argmin())
    if dt[fi]>20*60: continue
    d2=(LAT[fi]-la)**2+(LON[fi]-lo)**2; idx=int(d2.argmin()); h2,w2=np.unravel_index(idx,(H,W))
    if d2[h2,w2]<=(0.05**2): Y[fi,h2,w2]=1
log("FIRMS positive pixels=%d"%int(Y.sum()))

def spatial_z(x):
    out=np.zeros_like(x)
    for i in range(x.shape[0]):
        v=x[i];m=v[v>0];mu=m.mean() if m.size else 0;sd=m.std()+1e-6;out[i]=(v-mu)/sd
    return out
MIRf=F[:,mir];TIRf=F[:,tir];doz=MIRf-TIRf
doz_z=spatial_z(doz);mir_z=spatial_z(MIRf)
dprev=np.zeros_like(MIRf);dprev[1:]=MIRf[1:]-MIRf[:-1];tdiff_z=spatial_z(dprev)
def local_bg(x):
    s=np.zeros_like(x);sq=np.zeros_like(x);c=0
    for dh in(-1,0,1):
        for dw in(-1,0,1):
            if dh==0 and dw==0: continue
            sh=np.roll(np.roll(x,dh,axis=1),dw,axis=2);s+=sh;sq+=sh*sh;c+=1
    mu=s/c;var=np.clip(sq/c-mu*mu,0,None);return mu,np.sqrt(var)+1e-6
bg_mu,bg_sd=local_bg(MIRf)

def feats(ti,hi,wi):
    raw=np.stack([F[ti,k,hi,wi] for k in range(C)],1)
    pairs=[(a,b) for a in range(C) for b in range(a+1,C)]
    diffs=np.stack([raw[:,a]-raw[:,b] for (a,b) in pairs],1)
    spec=np.concatenate([raw,diffs],1)
    der=np.stack([MIRf[ti,hi,wi],TIRf[ti,hi,wi],doz[ti,hi,wi],mir_z[ti,hi,wi],doz_z[ti,hi,wi],
                  tdiff_z[ti,hi,wi],bg_mu[ti,hi,wi],bg_sd[ti,hi,wi],
                  (MIRf[ti,hi,wi]-bg_mu[ti,hi,wi])/bg_sd[ti,hi,wi]],1)
    return der,spec

# positives
pti,phi,pwi=np.where(Y==1)
nP=len(pti); log("positives=%d"%nP)
# the floor story: where do real fires sit in anomaly space?
log("real-fire doz_z pct: p50=%.2f p90=%.2f p99=%.2f ; MIR p50=%.1fK"%(
    np.percentile(doz_z[pti,phi,pwi],50),np.percentile(doz_z[pti,phi,pwi],90),
    np.percentile(doz_z[pti,phi,pwi],99),np.percentile(MIRf[pti,phi,pwi],50)))
# negative pool: non-fire, warm-ish
negmask=(Y==0)&(MIRf>275.0)
nti,nhi,nwi=np.where(negmask)
log("negative pool=%d"%len(nti))
# intensity-matched sampling: bin by (round MIR to 1K, round doz to 1K)
pmir=np.round(MIRf[pti,phi,pwi]).astype(int); pdoz=np.round(doz[pti,phi,pwi]).astype(int)
nmir=np.round(MIRf[nti,nhi,nwi]).astype(int); ndoz=np.round(doz[nti,nhi,nwi]).astype(int)
from collections import defaultdict
binmap=defaultdict(list)
for j in range(len(nti)): binmap[(nmir[j],ndoz[j])].append(j)
rng=np.random.RandomState(0); R=5
seln=[]; unmatched=0
for i in range(nP):
    key=(pmir[i],pdoz[i]); pool=binmap.get(key,[])
    if len(pool)==0:
        # widen by +/-1K
        pool=[]
        for dk in (-1,0,1):
            for dd in (-1,0,1):
                pool+=binmap.get((pmir[i]+dk,pdoz[i]+dd),[])
    if len(pool)==0: unmatched+=1; continue
    take=rng.choice(pool,size=min(R,len(pool)),replace=False); seln+=list(take)
seln=np.array(sorted(set(seln)))
log("matched negatives=%d (unmatched positives=%d)"%(len(seln),unmatched))
# build matrices
Pder,Pspec=feats(pti,phi,pwi)
Nder,Nspec=feats(nti[seln],nhi[seln],nwi[seln])
Xder=np.concatenate([Pder,Nder]); Xspec=np.concatenate([Pspec,Nspec]); Xcomb=np.concatenate([Xder,Xspec],1)
y=np.concatenate([np.ones(nP,int),np.zeros(len(seln),int)])
plat=np.concatenate([LAT[pti,phi,pwi],LAT[nti[seln],nhi[seln],nwi[seln]]])
plon=np.concatenate([LON[pti,phi,pwi],LON[nti[seln],nhi[seln],nwi[seln]]])
groups=(np.floor(plat/0.5).astype(int)*1000+np.floor(plon/0.5).astype(int))
log("design: pos=%d neg=%d (matched on MIR+MIR-TIR); groups=%d"%(int(y.sum()),int((1-y).sum()),len(np.unique(groups))))
# sanity: can intensity alone separate? (should be ~0.5 after matching)
def cvauc(X,y,groups,shuffle=False,seed=0):
    yy=y.copy()
    if shuffle: yy=np.random.RandomState(seed).permutation(yy)
    k=min(5,len(np.unique(groups)));
    if k<2: return None,None,k
    gkf=GroupKFold(n_splits=k);A=[]
    for tr,te in gkf.split(X,yy,groups):
        if yy[tr].sum()<5 or yy[te].sum()<2 or (1-yy[te]).sum()<2: continue
        c=HistGradientBoostingClassifier(max_iter=300,learning_rate=0.06,max_depth=4,l2_regularization=1.0,random_state=seed)
        c.fit(X[tr],yy[tr]);p=c.predict_proba(X[te])[:,1]
        try:A.append(roc_auc_score(yy[te],p))
        except:pass
    return (float(np.mean(A)),float(np.std(A)),k) if A else (None,None,k)
res={"n_frames":int(T),"n_positives":int(nP),"n_matched_neg":int(len(seln)),
     "real_fire_doz_z_p90":float(np.percentile(doz_z[pti,phi,pwi],90)),
     "v1_operating_point":{"candidates":48534,"real":10,"precision":10/48534,
        "note":"contextual front-end recovers only 10/767 FIRMS fires = SEVIRI floor"}}
for name,X in [("intensity_only_MIR_doz",Xder[:,[0,2]]),("DERIVED_0037",Xder),
               ("SPECTRAL_full_raw",Xspec),("COMBINED",Xcomb)]:
    m,s,k=cvauc(X,y,groups);res[name]={"cv_auc":m,"cv_std":s,"folds":k,"n_feat":X.shape[1]}
    log("%s  AUC=%s +/- %s (k=%s,feat=%d)"%(name,m,s,k,X.shape[1]))
pm,ps,pk=cvauc(Xspec,y,groups,shuffle=True,seed=7);res["PLACEBO_shuffled"]={"cv_auc":pm,"cv_std":ps}
log("PLACEBO=%s (must be ~0.5)"%pm)
json.dump(res,open(OUT,'w'),indent=1)
np.savez_compressed(FEAT_OUT,Xder=Xder,Xspec=Xspec,y=y,lat=plat,lon=plon,groups=groups)
print("RESULT_JSON_BEGIN");print(json.dumps(res,indent=1));print("RESULT_JSON_END")