#!/usr/bin/env python
"""Validate radio_class (size cut) vs CORNISH-S native SourceType; fix S fallback."""
import numpy as np
from astropy.table import Table
from astropy.coordinates import SkyCoord
import astropy.units as u
from collections import Counter
LOG=open('validate_s.log','w')
def p(*a):
    s=' '.join(str(x) for x in a);print(s);print(s,file=LOG);LOG.flush()

def sky(t,ra='RAJ2000',de='DEJ2000'):
    r=t[ra];d=t[de]
    if r.dtype.kind in ('S','U'):
        return SkyCoord([x.decode() if isinstance(x,bytes) else str(x) for x in r],
                        [x.decode() if isinstance(x,bytes) else str(x) for x in d],unit=(u.hourangle,u.deg))
    return SkyCoord(np.array(r,float)*u.deg,np.array(d,float)*u.deg)
PCAS=0.0048481
m=Table.read('master_phase1.fits')
mc=sky(m); dist=np.array(m['Dist'],float)
styp=np.array([str(x) for x in m['simple_type']])
cs=Table.read('cornish_south_catalog.fits')
ccs=sky(cs)
iS,dS,_=mc.match_to_catalog_sky(ccs); mS=dS<5*u.arcsec
Size=np.ma.filled(cs['Size'].astype(float),np.nan)[iS]
Stype=np.array([str(x) for x in cs['SourceType']])[iS]
# fixed S angular size: deconvolved Size if resolved, else beam 2.5" upper limit
BEAM=2.5
angS=np.where(np.isfinite(Size)&(Size>0),Size,BEAM)
Dpc=angS*dist*PCAS
cls=np.where(Dpc<=0.1,'UCHII',np.where(Dpc<=0.5,'compact','classical'))
# crosstab for matched S sources
sel=mS
p('Matched CORNISH-S:',int(sel.sum()))
p('\nSize-cut radio_class  vs  CORNISH-S SourceType (matched only):')
import itertools
types=['UCHII Region','HII Region','Diffuse HII Region','PN','MYSO','Infrared Quiet','Radio Star','--']
for rc in ['UCHII','compact','classical']:
    row=[]
    for st in types:
        row.append(int(((cls==rc)&(Stype==st)&sel).sum()))
    p(f'  {rc:10s}: '+' '.join(f'{t.split()[0][:6]}={n}' for t,n in zip(types,row)))
# how many of my UCHII (size) are CORNISH UCHII-labelled?
ucsel=sel&(cls=='UCHII')
p('\nOf my size-UCHII (S):',int(ucsel.sum()))
p('  CORNISH labels:',dict(Counter(Stype[ucsel]).most_common()))
# how many CORNISH "UCHII Region" matched RMS young, and their size class
cuc=sel&(Stype=='UCHII Region')
p('\nCORNISH-S UCHII Region matched to RMS young:',int(cuc.sum()))
p('  my size class:',dict(Counter(cls[cuc]).most_common()))
p('  RMS type:',dict(Counter(styp[cuc]).most_common()))
LOG.close()