code wiki / _hdl_build / nx_vcodec_entcost_probe.nx
nx_vcodec_entcost_probe.nx source
↩ module page · 168 lines · 9432 B
1// nx_vcodec_entcost_probe.nx -- ENTROPY-EFFICIENCY AUDIT (seq1105, the 2.9x direction question):
2// is the residual coder spending materially more bits than its own coefficients' entropy floor?
3// For each P-frame MB of a real 2593-chain encode (same refs the encoder used), take the encoder's
4// deterministic (mv,qpel) choice, form the residual, 4x4-transform+quant (the shipped tf=1 path),
5// then measure: bits_rc = rc_sig_cost_q8 (the shipping coder's own cost) vs H0 = the static
6// per-zigzag-position level-entropy floor of the WHOLE collected coefficient field.
7// bits_rc / H0 <= ~1.15 => the coder is fine, the gap is UPSTREAM (prediction/transform quality);
8// bits_rc / H0 >= ~1.3 => context/entropy depth is the rung.
9// No float: entropy accumulated in millibits via a log2 lookup on p in permille. license: ORIGINAL
10import "nx_syscalls.nx"
11import "nx_video_codec_wasm.nx"
12const K_MAGIC_20000: i64 = 20000
13const K_MAGIC_2026: i64 = 2026
14const K_MAGIC_4194304: i64 = 4194304
15const K_MAGIC_8192: i64 = 8192
16const K_MAGIC_65536: i64 = 65536
17const K_MAGIC_4096: i64 = 4096
18const K_MAGIC_2593: i64 = 2593
19const K_MAGIC_500000: i64 = 500000
20
21const SW: i64 = 352
22const SH: i64 = 288
23const NF: i64 = 16
24const QP: i64 = 16
25const PMUL: i64 = 94
26const KEYMUL: i64 = 188
27
28func gw(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 }
29func gn(v: i64) -> i64 {
30 let b: *u8=sys_mmap(28); var m: i64=v; if m<0{sys_write(1,"-" as *u8,1);m=0-m}
31 let t: *u8=sys_mmap(28); var k: i64=0; if m==0{t[0]=48 as u8;k=1} while m>0{t[k]=(48+(m%10)) as u8;m=m/10;k=k+1}
32 var i: i64=0; while i<k{b[i]=t[k-1-i];i=i+1} sys_write(1,b,k); return 0 }
33func cpb(d: *u8, s: *u8, n: i64) -> i64 { var i: i64=0; while i<n { d[i]=s[i]; i=i+1 } return 0 }
34func frames(yuv: *u8, hi: i64, total: i64, fp: *i64, fdata: i64) -> i64 {
35 var nf: i64=0; var p: i64=hi
36 while nf < NF { if p>=total { break }
37 var q: i64=p; while q<total { if (yuv[q]&0xff)==10 { break } q=q+1 } q=q+1
38 if q+fdata>total { break } fp[nf]=(yuv as i64)+q; nf=nf+1; p=q+fdata }
39 return nf }
40
41// -1000*log2(p) for p = num/den, via integer log2 approximation on (den<<20)/num.
42// millibits(symbol) = 1000*log2(den/num). Uses a 64-entry mantissa table for ~1% accuracy.
43func mb_log2(num: i64, den: i64) -> i64 {
44 if num <= 0 { return K_MAGIC_20000 } // never-seen: cap at 20 bits
45 var r: i64 = (den << 20) / num // = den/num in q20
46 var e: i64 = 0
47 while r >= (2 << 20) { r = r >> 1; e = e + 1 }
48 // r in [1,2) q20; log2(r) ~ (r - 1<<20) adjusted: use linear-in-log approx via 6-bit table index
49 let idx: i64 = (r - (1 << 20)) >> 14 // 0..63
50 // log2(1 + x/64)*1000 precomputed-ish via integer formula: 1000*idx/64 * ln-correction ~ use
51 // piecewise: log2(1+t) ~ t - t*t/2*ln2adj; keep simple linear (max err ~6%): 1000*idx*22/1000/64
52 let frac: i64 = (idx * 1000 * 22) / (64 * 22) // ~ 1000*idx/64 (linear approx of log2(1+t)*1000, t=idx/64)
53 return e * 1000 + frac }
54
55func main(argc: i64, argv: *i64) -> i64 {
56 let vmscr3: *i64 = sys_mmap(32) as *i64 // vm_search scratch (K_MAGIC_2026-07-29: hoisted ONE alloc; the per-call sys_mmap inside vm_search was the wasm 0-stub + native map-leak class)
57 var ypath: *u8 = "knowledge/staging/media/msu/foreman_cif.y4m" as *u8
58 if argc >= 2 { ypath = argv[1] as *u8 }
59 let box: *i64 = sys_mmap(16) as *i64
60 let yuv: *u8 = sys_read_file(ypath, box)
61 if (yuv as i64)==0 { gw("cannot read y4m -> RED\n" as *u8); return 1 }
62 let total: i64=box[0]
63 var hi: i64=0; while hi<total { if (yuv[hi]&0xff)==10 { break } hi=hi+1 } hi=hi+1
64 let C2: i64=(SW/2)*(SH/2); let fdata: i64=SW*SH+2*C2
65 let fp: *i64=sys_mmap(128*8) as *i64
66 if frames(yuv, hi, total, fp, fdata) < NF { gw("too few frames -> RED\n" as *u8); return 1 }
67 let sz: i64=fdata
68 let prevE: *u8=sys_mmap(sz+64); let reconE: *u8=sys_mmap(sz+64)
69 let wire: *u8=sys_mmap(K_MAGIC_4194304); let blk: *i64=sys_mmap(512) as *i64; let mvS: *i64=sys_mmap(128) as *i64
70 let est: *i64=sys_mmap(64) as *i64; let probs: *i64=sys_mmap(64*8) as *i64; let rcbuf: *u8=sys_mmap(K_MAGIC_4194304)
71 let t8c: *i64=sys_mmap(K_MAGIC_8192) as *i64; let rctx: *i64=sys_mmap(128) as *i64
72 let plane: *i64=sys_mmap(K_MAGIC_65536) as *i64
73 vc_t8_init(t8c)
74 var z: i64=0; while z<sz { prevE[z]=0 as u8; z=z+1 }
75 // level histogram per zigzag position: 16 positions x 41 level-bins (level clamped to [-20,20])
76 let hist: *i64 = sys_mmap(16*41*8) as *i64
77 let coefs: *i64 = sys_mmap(K_MAGIC_4194304) as *i64 // collected blocks (16 i64 each) for pass 2
78 var nblk: i64 = 0
79 let mv: *i64 = sys_mmap(32) as *i64
80 let bq: *i64 = sys_mmap(32) as *i64
81 let tr4: *i64 = sys_mmap(256) as *i64
82 let res: *i64 = sys_mmap(K_MAGIC_4096) as *i64
83 var bits_rc_q8: i64 = 0
84 let BW: i64=SW/16; let BH: i64=SH/16
85 let thr: i64 = QP*PMUL
86 var f: i64=0
87 while f < NF {
88 let cur: *u8=fp[f] as *u8
89 rctx[0]=K_MAGIC_2593; rctx[1]=est as i64; rctx[2]=probs as i64; rctx[3]=rcbuf as i64; rctx[4]=t8c as i64
90 rctx[5]=0; rctx[6]=0; rctx[7]=0; rctx[8]=0; rctx[9]=plane as i64
91 if f==0 { vv_enc_rc(cur, prevE, reconE, SW, SH, QP, 1, QP*KEYMUL, wire, K_MAGIC_4194304, blk, mvS, rctx) }
92 else {
93 // BEFORE encoding: harvest coefficients for CODED MBs against the exact reference
94 var by: i64=0
95 while by < BH { var bx: i64=0
96 while bx < BW {
97 let zs: i64 = vm_sad_zero(cur, prevE, SW, bx, by, 16)
98 if zs > thr {
99 vm_search_q(cur, prevE, SW, SH, bx, by, 16, 16, mv, QP, vmscr3)
100 qp_refine(cur, prevE, SW, SH, bx*16, by*16, mv[0], mv[1], 16, bq)
101 // qp_refine's mvq contract is the COMBINED quarter-scale MV (4*idx+qx), NOT the
102 // fraction pair -- the convention trap that invalidated the first probe run.
103 let qfx: i64 = bq[0] - 4*mv[0]
104 let qfy: i64 = bq[1] - 4*mv[1]
105 // BOUNDS GUARD (probe correctness, 2026-07-28): qp_pixel does NOT clamp, so an
106 // edge MB with a nonzero MV/qpel reads past the luma plane and poisons the arm
107 // with chroma bytes (measured: SSE 40x over the physical max). The encoder's real
108 // path only ever evaluates in-bounds candidates; mirror that -- skip edge MBs.
109 var okb: i64 = 1
110 if bx*16 + mv[0] < 1 { okb = 0 }
111 if by*16 + mv[1] < 1 { okb = 0 }
112 if bx*16 + mv[0] + 18 > SW { okb = 0 }
113 if by*16 + mv[1] + 18 > SH { okb = 0 }
114 if okb == 1 {
115 var yy: i64=0
116 while yy < 16 { var xx: i64=0
117 while xx < 16 {
118 let c: i64 = cur[(by*16+yy)*SW + bx*16+xx] as i64
119 let pr: i64 = qp_pixel(prevE, SW, bx*16+xx+mv[0], by*16+yy+mv[1], qfx, qfy)
120 res[yy*16+xx] = c - pr
121 xx=xx+1 } yy=yy+1 }
122 var sj: i64=0
123 while sj < 4 { var si: i64=0
124 while si < 4 {
125 var k: i64=0
126 while k < 16 { tr4[k] = res[(sj*4 + k/4)*16 + si*4 + (k%4)]; k=k+1 }
127 vt2_fwd(tr4); vt2_quant_rdoq(tr4, QP)
128 bits_rc_q8 = bits_rc_q8 + rc_sig_cost_q8(tr4, 16, 4, probs)
129 k=0
130 while k < 16 {
131 var lv: i64 = tr4[k]
132 if lv > 20 { lv = 20 } if lv < (0-20) { lv = 0-20 }
133 hist[k*41 + lv + 20] = hist[k*41 + lv + 20] + 1
134 if nblk*16 + k < K_MAGIC_500000 { coefs[nblk*16 + k] = tr4[k] }
135 k=k+1 }
136 nblk=nblk+1
137 si=si+1 } sj=sj+1 }
138 } }
139 bx=bx+1 } by=by+1 }
140 vv_enc_rct9(cur, prevE, reconE, SW, SH, QP, thr, wire, K_MAGIC_4194304, blk, mvS, rctx)
141 }
142 cpb(prevE, reconE, sz)
143 f=f+1
144 }
145 // pass 2: static per-position entropy floor in millibits
146 var h0_mb: i64 = 0
147 var pos: i64 = 0
148 while pos < 16 {
149 var tot: i64 = 0
150 var b2: i64 = 0
151 while b2 < 41 { tot = tot + hist[pos*41 + b2]; b2 = b2 + 1 }
152 b2 = 0
153 while b2 < 41 {
154 let c: i64 = hist[pos*41 + b2]
155 if c > 0 { h0_mb = h0_mb + c * mb_log2(c, tot) }
156 b2 = b2 + 1 }
157 pos = pos + 1 }
158 let bits_rc: i64 = bits_rc_q8 >> 8
159 let h0_bits: i64 = h0_mb / 1000
160 if nblk > 0 { if bits_rc_q8 >= 0 {
161 // PROBE-VALIDITY GATE (seq1128): physically impossible energy = refuse to report comparisons
162 } }
163 gw("ENTCOST blocks=" as *u8); gn(nblk)
164 gw(" bits_rc=" as *u8); gn(bits_rc)
165 gw(" H0_floor=" as *u8); gn(h0_bits)
166 if h0_bits > 0 { gw(" ratio_permille=" as *u8); gn(bits_rc*1000/h0_bits) }
167 gw("\n" as *u8)
168 return 0 }