"""
PHOENIX research — 12-channel candidate precision, CONFOUND-KILLED (v3).

v2 showed (powered, intensity-matched, placebo 0.50): full spectrum 0.79 > derived 0.69 AUC.
Open confound: real fires occur at specific times/illumination, so visible/NIR bands could encode
TIME-OF-DAY, not combustion. v3 kills it two ways:
  (1) SAME-FRAME matched negatives: each negative is drawn from the SAME frame as its positive,
      intensity-matched (MIR +/-2K, MIR-TIR +/-2K). Identical time, sun angle, atmosphere -> any
      spectral separation is the fire pixel vs a non-fire pixel under identical conditions.
  (2) Ablation: SPECTRAL_all vs SPECTRAL_THERMAL_ONLY (drop vis006/vis008/nir016 + their diffs).
      If the win survives thermal-only, it is combustion spectral structure, not solar illumination.
Honest spatial GroupKFold CV + placebo.
"""
import numpy as np, glob, os, re, sqlite3, time, json
from datetime import datetime, timezone
from collections import defaultdict
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_v3_results.json"
def log(m): print("[v3 %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']
VIS={'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'}
R=re.compile(r'-NA-(\d{14})')
def fts(n):
    m=R.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
    for f in ('%Y-%m-%dT%H:%M:%S','%a, %d %b %Y %H:%M:%S'):
        try: return datetime.strptime(ts.strip().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"%T)
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 positives=%d"%int(Y.sum()))
MIRf=F[:,mir];TIRf=F[:,tir];doz=MIRf-TIRf
def szf(x):
    out=np.zeros_like(x)
    for i in range(x.shape[0]):
        v=x[i];m=v[v>0]
        if m.size: out[i]=np.nan_to_num((v-m.mean())/(m.std()+1e-6),nan=0.0)
    return out
doz_z=szf(doz);mir_z=szf(MIRf)
def feats(ti,hi,wi,thermal_only=False):
    idxs=[k for k in range(C) if not(thermal_only and CH[k] in VIS)]
    raw=np.stack([F[ti,k,hi,wi] for k in idxs],1)
    pairs=[(a,b) for ai,a in enumerate(idxs) for b in idxs[ai+1:]]
    diffs=np.stack([F[ti,a,hi,wi]-F[ti,b,hi,wi] for (a,b) in pairs],1)
    return np.concatenate([raw,diffs],1)
def derived(ti,hi,wi):
    return np.stack([MIRf[ti,hi,wi],TIRf[ti,hi,wi],doz[ti,hi,wi],mir_z[ti,hi,wi],doz_z[ti,hi,wi]],1)
pti,phi,pwi=np.where(Y==1);nP=len(pti)
# same-frame intensity-matched negatives
rng=np.random.RandomState(0);negs=[];pos_keep=[]
for i in range(nP):
    fi=pti[i];pm=MIRf[fi,phi[i],pwi[i]];pd=doz[fi,phi[i],pwi[i]]
    fr=F[fi];m_ok=(Y[fi]==0)&(np.abs(MIRf[fi]-pm)<=2.0)&(np.abs(doz[fi]-pd)<=2.0)
    cand=np.argwhere(m_ok)
    if len(cand)==0: continue
    take=cand[rng.choice(len(cand),size=min(5,len(cand)),replace=False)]
    for (hh,ww) in take: negs.append((fi,hh,ww))
    pos_keep.append(i)
pk=np.array(pos_keep);negs=np.array(negs)
log("same-frame design: pos=%d neg=%d"%(len(pk),len(negs)))
pti2,phi2,pwi2=pti[pk],phi[pk],pwi[pk];nti,nhi,nwi=negs[:,0],negs[:,1],negs[:,2]
y=np.concatenate([np.ones(len(pk),int),np.zeros(len(negs),int)])
lat=np.concatenate([LAT[pti2,phi2,pwi2],LAT[nti,nhi,nwi]]);lon=np.concatenate([LON[pti2,phi2,pwi2],LON[nti,nhi,nwi]])
groups=(np.floor(lat/0.5).astype(int)*1000+np.floor(lon/0.5).astype(int))
Xder=np.concatenate([derived(pti2,phi2,pwi2),derived(nti,nhi,nwi)])
Xall=np.concatenate([feats(pti2,phi2,pwi2,False),feats(nti,nhi,nwi,False)])
Xtho=np.concatenate([feats(pti2,phi2,pwi2,True),feats(nti,nhi,nwi,True)])
Xint=Xder[:,[0,2]]
def cv(X,y,g,sh=False,sd=0):
    yy=np.random.RandomState(sd).permutation(y) if sh else y.copy()
    k=min(5,len(np.unique(g)))
    if k<2: return None,None,k
    A=[]
    for tr,te in GroupKFold(n_splits=k).split(X,yy,g):
        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=sd)
        c.fit(X[tr],yy[tr])
        try:A.append(roc_auc_score(yy[te],c.predict_proba(X[te])[:,1]))
        except:pass
    return (float(np.mean(A)),float(np.std(A)),k) if A else (None,None,k)
res={"design":"same-frame intensity-matched (kills time/illumination confound)","n_pos":int(len(pk)),"n_neg":int(len(negs)),"groups":int(len(np.unique(groups)))}
for nm,X in [("intensity_only",Xint),("DERIVED_0037",Xder),("SPECTRAL_all",Xall),("SPECTRAL_thermal_only",Xtho)]:
    m,s,k=cv(X,y,groups);res[nm]={"cv_auc":m,"cv_std":s,"folds":k,"n_feat":X.shape[1]};log("%s AUC=%s +/-%s feat=%d"%(nm,m,s,X.shape[1]))
pm,ps,_=cv(Xall,y,groups,sh=True,sd=7);res["PLACEBO"]={"cv_auc":pm,"cv_std":ps};log("PLACEBO=%s"%pm)
json.dump(res,open(OUT,'w'),indent=1);print("RESULT_JSON_BEGIN");print(json.dumps(res,indent=1));print("RESULT_JSON_END")