#!/usr/bin/env python3
"""A09 reproducible NASA stimulus audit and prospective factorial-design tools.
No participant outcomes accompany this package. No network calls except fetch.
"""
from __future__ import annotations
import argparse, csv, hashlib, io, json, math, pathlib, sys, urllib.parse, urllib.request
from collections import Counter, defaultdict
from datetime import datetime, timezone
import numpy as np
from scipy import stats

QUERY = ('select pl_name,hostname,default_flag,soltype,discoverymethod,pl_controv_flag,sy_snum,cb_flag,'
         'pl_orbper,pl_orbpererr1,pl_orbpererr2,pl_orbperlim,pl_orbsmax,pl_orbsmaxerr1,pl_orbsmaxerr2,'
         'pl_orbsmaxlim,pl_orbeccen,pl_orbeccenerr1,pl_orbeccenerr2,pl_orbeccenlim,pl_orbincl,pl_refname,disc_year '
         'from ps where default_flag=1 order by pl_name')
URL = 'https://exoplanetarchive.ipac.caltech.edu/TAP/sync?' + urllib.parse.urlencode({'query':QUERY,'format':'csv'})
ALLOWED_METHODS = {'Transit','Radial Velocity','Transit Timing Variations'}
REQUIRED = ['pl_name','hostname','default_flag','soltype','discoverymethod','pl_controv_flag','sy_snum','cb_flag',
            'pl_orbper','pl_orbsmax','pl_orbeccen','pl_orbincl','pl_refname','disc_year'] + [x+y for x in ['pl_orbper','pl_orbsmax','pl_orbeccen'] for y in ['err1','err2','lim']]
CELLS = [(0,0),(0,1),(1,0),(1,1)]

def number(v):
    try:
        x=float(v)
        return x if math.isfinite(x) else None
    except (TypeError,ValueError): return None

def audit_csv(path):
    p=pathlib.Path(path); raw=p.read_bytes()
    reader=csv.DictReader(io.StringIO(raw.decode('utf-8-sig')))
    if not reader.fieldnames or not set(REQUIRED)<=set(reader.fieldnames): raise ValueError('Missing required NASA columns')
    rows=list(reader)
    names=[r['pl_name'] for r in rows]
    if any(not n for n in names) or len(set(names))!=len(names): raise ValueError('Default rows must have unique nonempty planet names')
    if any(number(r['default_flag'])!=1 for r in rows): raise ValueError('Query must contain default_flag=1 only')
    counts={k:sum(number(r[k]) is None for r in rows) for k in ['pl_orbper','pl_orbsmax','pl_orbeccen','pl_orbincl']}
    masks=[('publishedConfirmed',lambda r:r['soltype']=='Published Confirmed'),
           ('nonControversial',lambda r:number(r['pl_controv_flag'])==0),
           ('singleStarNotCircumbinary',lambda r:number(r['sy_snum'])==1 and number(r['cb_flag'])==0),
           ('nonProjectedDiscoveryMethods',lambda r:r['discoverymethod'] in ALLOWED_METHODS),
           ('finitePositivePeriodAndAxis',lambda r:all(number(r[k]) is not None and number(r[k])>0 for k in ['pl_orbper','pl_orbsmax'])),
           ('ellipticEccentricity',lambda r:number(r['pl_orbeccen']) is not None and 0<=number(r['pl_orbeccen'])<.95),
           ('allThreeLimitFlagsZero',lambda r:all(number(r[k+'lim'])==0 for k in ['pl_orbper','pl_orbsmax','pl_orbeccen'])),
           ('parameterReferencePresent',lambda r:bool(r['pl_refname'].strip()))]
    selected=rows; stages=[]
    for label,test in masks:
        before=len(selected);selected=[r for r in selected if test(r)]
        stages.append({'criterion':label,'remaining':len(selected),'excludedAtThisStage':before-len(selected)})
    def error_usable(r):
        for k in ['pl_orbper','pl_orbsmax','pl_orbeccen']:
            v,hi,lo=number(r[k]),number(r[k+'err1']),number(r[k+'err2'])
            if hi is None or lo is None or hi<=0 or lo>=0: return False
            if k!='pl_orbeccen' and v+lo<=0:return False
            if k=='pl_orbeccen' and (v+lo<0 or v+hi>=1):return False
        return True
    enhanced=[r for r in selected if error_usable(r)]
    qs=lambda xs:{'minimum':float(min(xs)),'median':float(np.median(xs)),'maximum':float(max(xs))} if xs else None
    return {'analysisKind':'Executed public astronomical stimulus-database feasibility audit; no educational participant data',
            'query':QUERY,'queryUrl':URL,'sha256':hashlib.sha256(raw).hexdigest(),'sourceBytes':len(raw),'rows':len(rows),
            'uniquePlanetNames':len(set(names)),'solutionTypes':dict(sorted(Counter(r['soltype'] for r in rows).items())),
            'missingOrNonfiniteByVariable':counts,'selectionStages':stages,'coreStimulusCandidates':len(selected),
            'uncertaintyChallengeCandidates':len(enhanced),
            'coreEccentricityZeroRows':sum(number(r['pl_orbeccen'])==0 for r in selected),
            'coreWithoutStrictTwoSidedUncertainty':len(selected)-len(enhanced),
            'enhancedParameterRanges':{k:qs([number(r[k]) for r in enhanced]) for k in ['pl_orbper','pl_orbsmax','pl_orbeccen']},
            'selectionRuleNote':'Sequential exclusions, not additive marginal missingness; imaging/microlensing projected separations excluded by method. Blank eccentricities remain missing. Exact-zero eccentricities are not certified measured circular orbits.',
            'rightsBoundary':'Public TAP access and archive citation instructions verified; no blanket CC license asserted. Raw source CSV, third-party text and images are not redistributed. Only aggregate audit counts, code and original instruments accompany this edition.',
            'educationalOutcomesPresent':False,'actualTrialParticipants':0}

def factorial_estimates(block_cells):
    """Each row must be ordered 00,01,10,11 (immersion,revision).
    Teacher-block mean contrasts, across-block conservative variance; approximate t CI.
    """
    x=np.asarray(block_cells,dtype=float)
    if x.ndim!=2 or x.shape[1]!=4 or len(x)<3 or not np.isfinite(x).all(): raise ValueError('Need at least three complete four-cell instructor blocks')
    contrasts={'interaction':np.array([1,-1,-1,1]),'revisionMain':np.array([-.5,.5,-.5,.5]),'immersionMain':np.array([-.5,-.5,.5,.5])}
    out={}
    for label,c in contrasts.items():
        d=x@c; estimate=float(d.mean());se=float(d.std(ddof=1)/math.sqrt(len(d)))
        if se<=0:raise ValueError('Zero-variance contrast cannot support the declared t inference')
        critical=float(stats.t.ppf(.975,len(d)-1))
        out[label]={'estimate':estimate,'standardError':se,'degreesOfFreedom':len(d)-1,
                    'ci95':[estimate-critical*se,estimate+critical*se],
                    'twoSidedP':float(2*stats.t.sf(abs(estimate/se),len(d)-1))}
    return out

def sensitivity(seed=20261016,repetitions=10000):
    if repetitions<100:raise ValueError('At least 100 simulation replicates required')
    rng=np.random.default_rng(seed);out=[]
    for blocks in [16,32,48]:
        for rho in [.05,.10,.15]:
            v=rho+(1-rho)/24
            noise=rng.normal(0,math.sqrt(v),(repetitions,blocks,4))
            # Standardized class means; common instructor intercept cancels in every contrast.
            d=noise@np.array([1,-1,-1,1]);se=d.std(axis=1,ddof=1)/math.sqrt(blocks)
            crit=stats.t.ppf(.975,blocks-1)
            null=np.mean(abs(d.mean(axis=1)/se)>crit)
            power=np.mean(abs((d.mean(axis=1)+.4)/se)>crit)
            out.append({'instructorBlocks':blocks,'schools':blocks,'courseSections':4*blocks,'studentsPerSection':24,'plannedStudents':96*blocks,
                        'classICC':rho,'assumedInteractionSD':.4,'power':float(power),'powerMonteCarloSE':math.sqrt(power*(1-power)/repetitions),
                        'nullRejection':float(null),'nullMonteCarloSE':math.sqrt(null*(1-null)/repetitions),
                        'normalApproxInteractionMDE80':float((stats.norm.ppf(.975)+stats.norm.ppf(.8))*math.sqrt(4*v/blocks))})
    return {'analysisKind':'Executed artificial-data design sensitivity; not observed student outcomes','seed':seed,'replicatesPerScenario':repetitions,
            'numpyVersion':np.__version__,'assumptions':{'individualMarginalVariance':1,'equalClassSize':24,'completeOutcomes':True,'independentClassErrorsWithinInstructor':True,
            'commonInstructorInterceptCancels':True,'noInterference':True,'GaussianWorkingModel':True,'boundedRubricFloorCeilingNotModeled':True},'scenarios':out}

def orbit_xy(a,e,mean_anomaly):
    """Pedagogical two-body ellipse, AU; phase is stipulated, not observational ephemeris."""
    if not all(math.isfinite(float(v)) for v in [a,e,mean_anomaly]) or a<=0 or not 0<=e<.95:raise ValueError('Require a>0 and 0<=eccentricity<0.95, finite inputs')
    m=float(mean_anomaly)%(2*math.pi);lo=0.;hi=2*math.pi
    for _ in range(80):
        mid=(lo+hi)/2
        if mid-e*math.sin(mid)<m:lo=mid
        else:hi=mid
    E=(lo+hi)/2
    return {'xAU':a*(math.cos(E)-e),'yAU':a*math.sqrt(1-e*e)*math.sin(E),'eccentricAnomalyRadians':E,
            'assumption':'Two-body orbit; phase and viewing orientation are stipulated task variables, not actual sky positions.'}

def randomize(rows,seed,data_kind):
    if data_kind not in ['synthetic','study']:raise ValueError('Declare synthetic or study data')
    by=defaultdict(list);ids=[]; school_to_instructor={}; instructor_to_school={}
    for r in rows:
        if not r.get('school_id') or not r.get('instructor_id') or not r.get('section_id'):raise ValueError('Missing roster identifiers')
        school,ins=r['school_id'],r['instructor_id']
        if (school in school_to_instructor and school_to_instructor[school]!=ins) or (ins in instructor_to_school and instructor_to_school[ins]!=school):raise ValueError('Require one instructor block per school')
        school_to_instructor[school]=ins;instructor_to_school[ins]=school
        by[r['instructor_id']].append(r);ids.append(r['section_id'])
    if len(ids)!=len(set(ids)) or any(len(v)!=4 for v in by.values()):raise ValueError('Require four unique sections per instructor')
    if data_kind=='study' and len(by)!=32:raise ValueError('Declared study requires 32 instructor blocks; amend prospectively for another design')
    rng=np.random.default_rng(seed);out=[]
    for instructor in sorted(by):
        sections=sorted(by[instructor],key=lambda x:x['section_id'])
        order=rng.permutation(4)
        for row,index in zip(sections,order):out.append({**row,'immersion':CELLS[index][0],'revision':CELLS[index][1]})
    return {'seed':seed,'dataKind':data_kind,'assignmentUnit':'course section','blocks':len(by),'assignments':out,
            'notice':'A local randomization utility, not evidence of ethics approval, consent, recruitment or actual assignment.'}

def analyze_complete(rows,data_kind):
    grouped=defaultdict(list);assign={};seen=set();school_to_instructor={};instructor_to_school={}
    for r in rows:
        sid=r.get('student_id')
        if not r.get('section_id') or not sid or sid in seen:raise ValueError('Missing section/student or globally repeated student ID')
        seen.add(sid)
        ins=r.get('instructor_id');cell=(number(r.get('immersion')),number(r.get('revision')))
        school=r.get('school_id')
        if not school or not ins or cell not in CELLS:raise ValueError('Invalid school, instructor or binary assignment')
        if (school in school_to_instructor and school_to_instructor[school]!=ins) or (ins in instructor_to_school and instructor_to_school[ins]!=school):raise ValueError('Require one instructor block per school')
        school_to_instructor[school]=ins;instructor_to_school[ins]=school
        score=number(r.get('delayed_score'))
        if score is None or not 0<=score<=20:raise ValueError('Complete-cohort tool rejects missing/out-of-range scores; use prespecified missing-data analysis separately')
        sec=r['section_id']
        if sec in assign and assign[sec]!=(ins,cell):raise ValueError('Inconsistent section assignment')
        assign[sec]=(ins,cell);grouped[sec].append(score)
    if any(len(v)!=24 for v in grouped.values()):raise ValueError('Requires all 24 baseline-cohort delayed outcomes per section')
    blocks=defaultdict(dict)
    for sec,scores in grouped.items():
        ins,cell=assign[sec]
        if cell in blocks[ins]:raise ValueError('Duplicate treatment cell in instructor block')
        blocks[ins][cell]=np.mean(scores)
    if any(set(v)!=set(CELLS) for v in blocks.values()):raise ValueError('Incomplete instructor factorial block')
    if data_kind=='study' and len(blocks)!=32:raise ValueError('Declared study requires 32 instructor blocks')
    if data_kind not in ['synthetic','study']:raise ValueError('Declare data kind')
    result=factorial_estimates([[blocks[k][c] for c in CELLS] for k in sorted(blocks)])
    return {'analysisKind':'Complete-cohort blocked factorial estimator','dataKind':data_kind,'blocks':len(blocks),
            'students':len(rows),'contrasts':result,'primaryContrast':'interaction','outcome':'delayed 0–20 transfer score',
            'scope':'Approximate across-instructor t inference; complete cohorts only. No imputation, mediation, equivalence test or participant certification implemented.'}

def main():
    p=argparse.ArgumentParser(description=__doc__);sp=p.add_subparsers(dest='command',required=True)
    for cmd in ['audit','design','fetch','randomize','analyze-complete']:
        s=sp.add_parser(cmd);s.add_argument('--output',required=True)
        if cmd in ['audit','randomize','analyze-complete']:s.add_argument('--input',required=True)
        if cmd in ['randomize','analyze-complete']:s.add_argument('--data-kind',required=True,choices=['synthetic','study'])
        if cmd in ['randomize','design']:s.add_argument('--seed',type=int,default=20261016)
        if cmd=='design':s.add_argument('--replicates',type=int,default=10000)
    a=p.parse_args()
    if a.command=='fetch':
        req=urllib.request.Request(URL,headers={'User-Agent':'JournalOS-A09-research-reproduction/1.0'})
        raw=urllib.request.urlopen(req,timeout=90).read();pathlib.Path(a.output).write_bytes(raw)
        print(json.dumps({'url':URL,'sha256':hashlib.sha256(raw).hexdigest(),'bytes':len(raw),'retrievedAt':datetime.now(timezone.utc).isoformat()}));return
    if a.command=='audit':out=audit_csv(a.input)
    elif a.command=='design':out=sensitivity(a.seed,a.replicates)
    else:
        rows=list(csv.DictReader(open(a.input,newline='')))
        out=randomize(rows,a.seed,a.data_kind) if a.command=='randomize' else analyze_complete(rows,a.data_kind)
    pathlib.Path(a.output).write_text(json.dumps(out,indent=2,ensure_ascii=False)+'\n')
if __name__=='__main__':main()