"""Independent standard-library check of cohort and diagnostic event counts.

Reads the selected data but uses no pandas or primary-audit functions. Outputs
aggregate counts only. No sample identifiers or clinical rows are saved.
"""
import collections,csv,gzip,json,math,pathlib,sys
ROOT=pathlib.Path(__file__).resolve().parents[1]
def number(v):
 try:return float(v)
 except (ValueError,TypeError):return math.nan
def table(p):
 op=gzip.open if str(p).endswith('.gz') else open
 with op(p,'rt') as f:
  rr=csv.reader(f,delimiter='\t');head=next(rr)[1:]
  rows={r[0]:[float(v) for v in r[1:]] for r in rr}
 return head,rows
with (ROOT/'data/metadata.csv').open() as f:meta=list(csv.DictReader(f))
dh,dr=table(sys.argv[1]);rh,rr=table(sys.argv[2]);th,tr=table(ROOT/'data/taxonomy.tsv.gz')
def mapping(headers,assay,suffix):
 byid={r['External ID']:r for r in meta if r['data_type']==assay}
 found=[(i,h[:-len(suffix)],byid[h[:-len(suffix)]]) for i,h in enumerate(headers) if h[:-len(suffix)] in byid]
 normal={r['site_sub_coll'] for _,h,r in found if not h.endswith('_P')}
 groups=collections.defaultdict(list)
 for i,h,r in found:
  if '_TR' not in h and not(h.endswith('_P') and r['site_sub_coll'] in normal):groups[r['site_sub_coll']].append((i,r))
 return {s:v[0] for s,v in groups.items() if s and len(v)==1}
dm=mapping(dh,'metagenomics','_level4ec');rm=mapping(rh,'metatranscriptomics','_level4ec');tm=mapping(th,'metagenomics','_profile')
species={k.rsplit('|',1)[-1][3:]:v for k,v in tr.items() if k.rsplit('|',1)[-1].startswith('s__')}
paired=[]
for s in sorted(set(dm)&set(rm)):
 di,d=dm[s];ri,r=rm[s]
 if d['diagnosis']!='CD' or not math.isfinite(number(d['week_num'])):continue
 assert (d['Participant ID'],d['week_num'],d['diagnosis'])==(r['Participant ID'],r['week_num'],r['diagnosis'])
 ti=tm[s][0] if s in tm else None
 valid=ti is not None and 99<=sum(v[ti] for v in species.values())<=101 and sum(v[ti]>0 for v in species.values())>1
 paired.append({'d':di,'r':ri,'t':ti,'m':d,'rm':r,'valid':valid})
filtered=[x for x in paired if number(x['m']['reads_filtered'])>=1e6 and number(x['rm']['reads_filtered'])>=1e6]
by=collections.defaultdict(list)
for x in filtered:
 if x['valid']:by[x['m']['Participant ID']].append(x)
edges=[]
for xs in by.values():
 xs.sort(key=lambda x:number(x['m']['week_num']))
 assert len({x['m']['week_num'] for x in xs})==len(xs)
 edges.extend((a,b) for a,b in zip(xs,xs[1:]) if 2<=number(b['m']['week_num'])-number(a['m']['week_num'])<=12)
res={'paired':len(paired),'participants':len({x['m']['Participant ID'] for x in filtered}),'filtered':len(filtered),'intervals':len(edges),'events':{}}
for ec in ('2.8.3.8','2.7.2.7','2.8.3.9','2.3.1.19'):
 rows={k:v for k,v in dr.items() if k.split(':',1)[0]==ec and '|g__' in k}
 count=0;fc=0;people=set()
 for a,b in edges:
  ds={k:v[a['d']] for k,v in rows.items() if v[a['d']]>0}
  if len(ds)<2:continue
  dom=max(ds,key=ds.get)
  if ds[dom]/sum(ds.values())<0.5:continue
  tx=species.get(dom.rsplit('s__',1)[-1])
  if tx is None or tx[a['t']]<0.1 or tx[b['t']]>0.5*tx[a['t']]:continue
  if rr.get(dom,[0]*len(rh))[a['r']]<=0:continue
  if not any(k!=dom and k in rr and rr[k][a['r']]>0 and rr[k][b['r']]>0 for k in ds):continue
  count+=1;people.add(a['m']['Participant ID'])
  if all(math.isfinite(number(x['m']['fecalcal'])) and number(x['m']['fecalcal'])>=0 and number(x['m']['fecalcal'])==number(x['rm']['fecalcal']) for x in (a,b)):fc+=1
 res['events'][ec]={'count':count,'participants':len(people),'fecalcal':fc}
primary=json.loads((ROOT/'results/audit.json').read_text());s=primary['reads_ge1million_each_assay']
checks=[res['paired']==primary['all_nonredundant']['paired_CD_specimens'],res['participants']==s['participants'],res['filtered']==s['paired_CD_specimens'],res['intervals']==s['consecutive_paired_specimen_intervals_2to12weeks']]
key='taxon_decline_ge50pct_baseline_ge0.1pct_alternative_RNA_both_times'
for ec,v in res['events'].items():
 e=s['enzymes'][ec]
 checks.extend([v['count']==e['events'].get(key,0),v['participants']==e['event_participants'].get(key,0),v['fecalcal']==e['events'].get(key+'_fecalcal_both_times',0)])
assert all(checks),'Independent eligibility calculation disagrees'
res['comparisons']=len(checks);res['all_agree']=all(checks)
(ROOT/'results/validation.json').write_text(json.dumps(res,indent=2)+'\n')
print(json.dumps({'independent_checks':len(checks),'all_agree':True}))