"""Session 29: fixed-inventory adjacency baseline and section heterogeneity.
Python + NumPy; run beside lp_core.py. No network or plaintext ranking.
The conditional allocation model is explicitly a Poisson approximation.
"""
import hashlib, itertools, json, math
from lp_core import ROOT, LET, lp, vals, np

SEED=3301029
PERMUTATIONS=9999
ALLOCATIONS=99999

def repeats(c): return int(np.count_nonzero(c[:-1]==c[1:]))
def moments(c):
    n=len(c);f=np.bincount(c,minlength=29).astype(float)
    assert n>=4
    a=f*(f-1);p2=float(a.sum()/(n*(n-1)))
    p3=float((a*(f-2)).sum()/(n*(n-1)*(n-2)))
    p22=float((a.sum()**2-(a*a).sum()+(a*(f-2)*(f-3)).sum())/(n*(n-1)*(n-2)*(n-3)))
    m=n-1;overlap=n-2;disjoint=m*(m-1)//2-overlap
    variance=m*p2*(1-p2)+2*overlap*(p3-p2*p2)+2*disjoint*(p22-p2*p2)
    return float(a.sum()/n),float(variance)

def controls():
    rows=[]
    for c in ([0,0,1,1],[0,0,0,1,2],[0,0,1,1,2],[0,1,2,3],[0,0,0,0]):
        unique=set(itertools.permutations(c))
        d=np.array([repeats(np.array(x)) for x in unique])
        mu,var=moments(np.array(c))
        assert np.isclose(mu,d.mean()) and np.isclose(var,d.var()),(c,mu,var)
        rows.append({'inventory':c,'unique_orders':len(unique),'mean':mu,'variance':var,'passed':True})
    return rows

def pearson(x,expected):return ((x-expected)**2/expected).sum(axis=-1)
def allocation_comparison(obs,weights,rng):
    n=int(obs.sum());prob=weights/weights.sum();expected=n*prob
    simulated=rng.multinomial(n,prob,size=ALLOCATIONS)
    observed=float(pearson(obs,expected));statistics=pearson(simulated,expected)
    ge=int(np.count_nonzero(statistics>=observed))
    # Same first-four / last-four contrast used in Session 25, including uncertainty.
    k=int(obs[:4].sum());p=float(prob[:4].sum())
    probabilities=np.array([math.comb(n,j)*p**j*(1-p)**(n-j) for j in range(n+1)])
    two_sided=float(probabilities[probabilities<=probabilities[k]+1e-15].sum())
    return dict(expected=expected.tolist(),pearson_statistic=observed,
        simulated_at_or_above=ge,monte_carlo_upper_tail=(1+ge)/(ALLOCATIONS+1),
        first_four_contrast=dict(observed=k,expected=n*p,conditional_binomial_two_sided=two_sided),
        model='Independent Poisson repeat counts with section means rho * weight; conditioning on total repeats removes fitted rho and yields a multinomial. Approximation for overlapping adjacency events, not an exact fixed-inventory cipher model.')

def main():
    checks=controls()
    pieces=[(i,vals(s)) for i,s in enumerate('%'.join(lp[:56]).split('$')) if len(vals(s))>100]
    assert [i for i,c in pieces]==[0,1,2,4,5,6,7,8]
    rows=[];sims=[];rng=np.random.default_rng(SEED)
    for i,c in pieces:
        mu,var=moments(c);o=repeats(c)
        a=np.array([repeats(rng.permutation(c)) for _ in range(PERMUTATIONS)])
        le=int(np.count_nonzero(a<=o));sims.append(a)
        rows.append(dict(section=i,runes=len(c),pairs=len(c)-1,observed=o,
            counts=np.bincount(c,minlength=29).tolist(),analytic_expected=mu,analytic_variance=var,
            observed_over_expected=o/mu,observed_rate=o/(len(c)-1),
            permutation=dict(mean=float(a.mean()),variance=float(a.var()),
                central_95_percent=np.quantile(a,[.025,.975]).tolist(),minimum=int(a.min()),
                at_or_below=le,lower_tail=(1+le)/(PERMUTATIONS+1),
                eight_section_bonferroni=min(1,8*(1+le)/(PERMUTATIONS+1)))))
    obs=np.array([r['observed'] for r in rows]);expected=np.array([r['analytic_expected'] for r in rows])
    totals=sum(sims);o=int(obs.sum());mu=float(expected.sum());var=sum(r['analytic_variance'] for r in rows)
    assert o==86 and sum(r['pairs'] for r in rows)==12939
    rate=o/mu;heterogeneity=allocation_comparison(obs,expected,np.random.default_rng(SEED+1))
    # Sensitivity only: replace frequency-adjusted exposure by the number of pairs.
    unadjusted=allocation_comparison(obs,np.array([r['pairs'] for r in rows]),np.random.default_rng(SEED+2))
    out=dict(seed=SEED,permutations=PERMUTATIONS,allocations=ALLOCATIONS,labels=LET,
        scope='Eight substantial LP2 sections across pages 0..55. Excludes nine-rune segment and cross-section edges. All non-rune material consumes no positions; adjacent means rune-stream adjacency, not necessarily neighboring printed glyphs.',
        permutation_null='Uniformly random ordering of each section with its exact rune inventory fixed; independent sections. Does not preserve per-page inventories, word content or ordered-pair counts.',
        expectation_formula='E[D_s] = sum_a n_sa (n_sa - 1) / N_s',
        sections=rows,total=dict(runes=sum(r['runes'] for r in rows),pairs=sum(r['pairs'] for r in rows),observed=o,
            analytic_expected=mu,analytic_standard_deviation=math.sqrt(var),
            observed_over_expected=rate,fewer_than_expected_fraction=1-rate,
            permutation_mean=float(totals.mean()),central_95_percent=np.quantile(totals,[.025,.975]).tolist(),
            minimum=int(totals.min()),at_or_below=int(np.count_nonzero(totals<=o)),
            lower_tail=(1+int(np.count_nonzero(totals<=o)))/(PERMUTATIONS+1)),
        frequency_adjusted_heterogeneity=heterogeneity,pair_exposure_sensitivity=unadjusted,
        exhaustive_moment_controls=checks,
        transcription_sha256=hashlib.sha256((ROOT/'liber-primus__transcription--master.txt').read_bytes()).hexdigest(),
        limits=['The suppression is already documented in Sessions 4, 5 and 25; this is not a new discovery or independent validation dataset.',
            'Monte Carlo resolution is 1/(B+1); no claim that a zero exceedance count means probability zero.',
            'Global fixed-inventory randomization is an explicit comparison model, not a description of all ciphers.',
            'The fitted 86/E ratio measures this ciphertext relative to its shuffle baseline, not removed plaintext letters or a key parameter.',
            'Section heterogeneity uses a Poisson approximation and has limited information with 86 events; nonsignificance does not establish identical algorithms.',
            'No fresh full-page image transcription, decryption, source attribution or key search.'])
    (ROOT/'repeat_frequency_results.json').write_text(json.dumps(out,indent=2)+'\n')
    print(json.dumps({k:out[k] for k in ['total','frequency_adjusted_heterogeneity','pair_exposure_sensitivity']},indent=2))
    for r in rows:print(r['section'],r['runes'],r['observed'],round(r['analytic_expected'],3),r['permutation'])
if __name__=='__main__':main()