"""Independent public-data feasibility audit; no RNA response or outcome tests.

Author R code was read, not executed. Implement its pilot/techrep exclusions
and strict site_sub_coll matching independently, rejecting ambiguities.
Inputs are supplied as CLI paths to avoid bundling individual clinical rows.
"""
import argparse, collections, hashlib, json, pathlib
import numpy as np
import pandas as pd

p = argparse.ArgumentParser()
p.add_argument('--metadata', required=True)
p.add_argument('--dna', required=True)
p.add_argument('--rna', required=True)
p.add_argument('--taxonomy', required=True)
p.add_argument('--output', required=True)
a = p.parse_args()
meta = pd.read_csv(a.metadata, dtype=str, keep_default_na=False)
frames = {
    'DNA': pd.read_csv(a.dna, sep='\t', index_col=0),
    'RNA': pd.read_csv(a.rna, sep='\t', index_col=0),
    'taxonomy': pd.read_csv(a.taxonomy, sep='\t', index_col=0),
}
audit = {}

def load(label, assay, ending):
    frame = frames[label]
    mm = meta[meta.data_type.eq(assay)].set_index('External ID')
    assert mm.index.is_unique
    records = []
    missing = 0
    for col in frame:
        assert col.endswith(ending)
        external = col[:-len(ending)]
        if external not in mm.index:
            missing += 1
            continue
        row = mm.loc[external]
        records.append(dict(external=external, column=col,
            specimen=row.site_sub_coll, participant=row['Participant ID'],
            week=pd.to_numeric(row.week_num, errors='coerce'),
            diagnosis=row.diagnosis,
            reads=pd.to_numeric(row.reads_filtered, errors='coerce'),
            fecalcal=pd.to_numeric(row.fecalcal, errors='coerce')))
    nonpilot = {r['specimen'] for r in records if not r['external'].endswith('_P')}
    tech = [r for r in records if '_TR' in r['external']]
    pilot = [r for r in records if r['external'].endswith('_P') and r['specimen'] in nonpilot]
    retained = [r for r in records if '_TR' not in r['external'] and not
        (r['external'].endswith('_P') and r['specimen'] in nonpilot)]
    grouped = collections.defaultdict(list)
    for r in retained:
        grouped[r['specimen']].append(r)
    ambiguous = {s for s, rr in grouped.items() if len(rr) != 1 or not s}
    clean = {s: rr[0] for s, rr in grouped.items() if s not in ambiguous}
    audit[label] = dict(headers=len(frame.columns), missing_metadata=missing,
        technical_replicates_removed=len(tech), duplicate_pilots_removed=len(pilot),
        ambiguous_specimens_excluded=len(ambiguous), retained_specimens=len(clean),
        retained_CD_specimens=sum(r['diagnosis']=='CD' for r in clean.values()),
        low_or_missing_read_count=sum(not np.isfinite(r['reads']) or r['reads']<1e6 for r in clean.values()))
    return clean

dna = load('DNA', 'metagenomics', '_level4ec')
rna = load('RNA', 'metatranscriptomics', '_level4ec')
tax = load('taxonomy', 'metagenomics', '_profile')
paired = []
inconsistent = 0
for s in sorted(dna.keys() & rna.keys()):
    d, r = dna[s], rna[s]
    if (d['participant'], d['week'], d['diagnosis']) != (r['participant'], r['week'], r['diagnosis']):
        inconsistent += 1
        continue
    if d['diagnosis'] != 'CD' or not np.isfinite(d['week']):
        continue
    paired.append(dict(specimen=s, participant=d['participant'], week=d['week'],
        dna=frames['DNA'][d['column']], rna=frames['RNA'][r['column']],
        qc=np.isfinite(d['reads']) and np.isfinite(r['reads']) and min(d['reads'],r['reads'])>=1e6,
        fecalcal=d['fecalcal'], fecalcal_agrees=(d['fecalcal']==r['fecalcal']),
        tax=frames['taxonomy'][tax[s]['column']] if s in tax else None))

# Species rows only; never sum nested genus/family/kingdom levels.
species_rows = {i.split('|')[-1][3:]: i for i in frames['taxonomy'].index
    if i.split('|')[-1].startswith('s__')}
assert len(species_rows) == sum(i.split('|')[-1].startswith('s__') for i in frames['taxonomy'].index)
tax_sums = frames['taxonomy'].loc[list(species_rows.values())].sum(axis=0)
for r in paired:
    values = r['tax'].reindex(list(species_rows.values())) if r['tax'] is not None else pd.Series(dtype=float)
    r['tax_valid'] = 99 <= values.sum() <= 101 and (values>0).sum()>1
targets = ['2.8.3.8','2.7.2.7','2.8.3.9','2.3.1.19']

def summarize(records):
    by_week = collections.defaultdict(list)
    for r in records: by_week[(r['participant'],r['week'])].append(r)
    ties = {k for k,v in by_week.items() if len(v)>1}
    longitudinal = [r for r in records if (r['participant'],r['week']) not in ties and r['tax_valid']]
    by_person = collections.defaultdict(list)
    for r in longitudinal: by_person[r['participant']].append(r)
    edges = []
    for person, rs in by_person.items():
        rs.sort(key=lambda r:r['week'])
        for first, second in zip(rs,rs[1:]):
            if 2 <= second['week']-first['week'] <= 12: edges.append((first,second))
    enzymes = {}
    for ec in targets:
        base = next(i for i in frames['DNA'].index if i.split(':')[0]==ec and '|' not in i)
        dr = [i for i in frames['DNA'].index if i.startswith(base+'|g__')]
        rr = [i for i in frames['RNA'].index if i.startswith(base+'|g__')]
        missing_names = {i.split('s__')[-1] for i in set(dr+rr) if i.split('s__')[-1] not in species_rows}
        counts = collections.Counter()
        persons = collections.defaultdict(set)
        coverage = collections.Counter()
        for r in records:
            dd = r['dna'].reindex(dr,fill_value=0)
            rd = r['rna'].reindex(rr,fill_value=0)
            if (dd>0).sum()>=2: coverage['two_DNA_contributors']+=1
            if (rd>0).sum()>=2: coverage['two_RNA_contributors']+=1
        for first, second in edges:
            person = first['participant']
            def add(k): counts[k]+=1; persons[k].add(person)
            dd = first['dna'].reindex(dr,fill_value=0)
            if (dd>0).sum()<2 or dd.sum()<=0: continue
            dominant = dd.idxmax()
            if dd[dominant]/dd.sum()<0.5: continue
            add('baseline_two_DNA_contributors_and_dominant_share_ge50pct')
            species = dominant.split('s__')[-1]
            if species not in species_rows:
                add('dominant_species_missing_taxonomy_annotation')
                continue
            tr = species_rows[species]
            t0,t1 = first['tax'][tr],second['tax'][tr]
            for threshold in [0.1,0.5,1.0]:
                if t0 < threshold or t1 > 0.5*t0: continue
                key = f'taxon_decline_ge50pct_baseline_ge{threshold:g}pct'
                add(key)
                # Existing alternative genes at baseline; detectable transcripts
                # required at both ends to support a later paired-response test.
                alternatives = [i for i in dd[dd>0].index if i != dominant]
                if first['rna'].get(dominant,0)<=0: continue
                add(key+'_initial_dominant_RNA')
                eligible = [i for i in alternatives if first['rna'].get(i,0)>0 and second['rna'].get(i,0)>0]
                if eligible:
                    add(key+'_alternative_RNA_both_times')
                    if all(np.isfinite(z['fecalcal']) and z['fecalcal']>=0 and z['fecalcal_agrees'] for z in [first,second]):
                        add(key+'_alternative_RNA_both_times_fecalcal_both_times')
        enzymes[ec] = dict(coverage=dict(coverage),
            functional_species_names_absent_from_taxonomy=len(missing_names),
            events=dict(counts), event_participants={k:len(v) for k,v in persons.items()})
    fc_observed = [r for r in records if np.isfinite(r['fecalcal']) and r['fecalcal']>=0 and r['fecalcal_agrees']]
    fc_by_person = collections.defaultdict(list)
    for r in fc_observed: fc_by_person[r['participant']].append(r)
    fc_only_intervals = 0
    for rs in fc_by_person.values():
        rs.sort(key=lambda r:r['week'])
        fc_only_intervals += sum(2 <= second['week']-first['week'] <= 12 for first,second in zip(rs,rs[1:]))
    return dict(paired_CD_specimens=len(records),participants=len({r['participant'] for r in records}),
        paired_taxonomy_specimens=sum(r['tax'] is not None for r in records),
        invalid_or_single_species_taxonomy_profiles_excluded=sum(not r['tax_valid'] for r in records),
        tied_participant_week_groups_excluded=len(ties),
        consecutive_paired_specimen_intervals_2to12weeks=len(edges),
        consecutive_intervals_with_numeric_fecalcal_both_times=sum(all(np.isfinite(z['fecalcal']) and z['fecalcal']>=0 and z['fecalcal_agrees'] for z in e) for e in edges),
        fecalcal_numeric_same_value_in_both_assay_metadata=len(fc_observed),
        fecalcal_participants=len({r['participant'] for r in fc_observed}),
        intervals_consecutive_only_among_fecalcal_observations_2to12weeks=fc_only_intervals,
        enzymes=enzymes)

out = dict(scope='Independent feasibility calculations only. No compensation, RNA response, inflammation association, metabolite production or causal effect tested. Species declines are compositional relative-abundance declines, not absolute population losses. Positive entries are detections, not a validated expression limit. EC annotations remain broader than validated butyrate-specific genes.',
    matching='Known header endings removed; exact assay-specific External ID metadata lookup. Exclude _TR and duplicate _P pilot profiles using author code conventions; strict specimen site_sub_coll matching; reject ambiguous specimens and participant/week/diagnosis conflicts. Do not average pilot/techrep abundances. Do not impute weeks. Exclude tied weeks for event ordering.',
    input_sha256={name:hashlib.sha256(pathlib.Path(path).read_bytes()).hexdigest() for name,path in [('metadata',a.metadata),('DNA_selected_rows',a.dna),('RNA_selected_rows',a.rna),('taxonomy',a.taxonomy)]},
    input_audit=audit, conflicting_pairs=inconsistent,
    species_relative_abundance_sum_min=float(tax_sums.min()),species_relative_abundance_sum_max=float(tax_sums.max()),
    all_nonredundant=summarize(paired),
    reads_ge1million_each_assay=summarize([r for r in paired if r['qc']]))
serialized = json.dumps(out,indent=2,default=lambda v:v.item())
pathlib.Path(a.output).write_text(serialized+'\n')
print(serialized)