CLAUDE: CuasLogFold - LoG-into-kernel fold + synthetic impulse/grating...

CLAUDE: CuasLogFold - LoG-into-kernel fold + synthetic impulse/grating self-test (roadmap RT-seed step 1)

Findings (all validated to stated precision by the self-test, exit 0):
- Per-bin TD quad-multiply == pixel convolution with cos(pi*dx/16)*cos(pi*dy/16)
  ATTENUATED taps (matches norm_sym_weights in calculateCLTKernel; shift-tap
  error exactly 1-cos(pi/16) before compensation). getLogTd() therefore
  pre-compensates LoG taps by 1/cosAtten so the EFFECTIVE filter is the true
  pixel LoG (machine-exact at non-bin gratings).
- Residual floors: symmetric-extension wrap ~0.5-2.5% (inherent to the lapped
  structure); folded-vs-pixel-LoG A/B cross-term ~7% on a pessimistic spread
  kernel (test 5 - predicts the full-run A/B floor).
- clt_symmetrize_kernel quadrant-3 (AA) sign pattern looks like a latent
  copy-paste of the AS expression: pure-AA kernel test - math sign 2.6e-3,
  code sign 1.9e-1. Production kernel files likely carry the wrong AA sign.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

src/main/java/com/elphel/imagej/cuas/rt/CuasLogFold.java

+package com.elphel.imagej.cuas.rt;
+
+import com.elphel.imagej.tileprocessor.DttRad2;
+import com.elphel.imagej.tileprocessor.ImageDttCPU;
+import com.elphel.imagej.tileprocessor.ImageDttParameters;
+
+/**
+ * Fold the CUAS LoG pre-filter into the aberration-correction CLT kernels
+ * ("LoG for free" in convert_direct) + a headless synthetic self-test.
+ *
+ * Math reference: attic/imagej-elphel-internal/handoffs/
+ * 2026-07-05_mclt_dtt_convolution_reference.md (esp. section 4). The LoG is
+ * even in both axes, so of its 4 DTT-III symmetry classes only the SS quadrant
+ * survives; folding is a per-bin REAL scalar multiply of all 4 quadrants of
+ * every stored kernel (rq = cq * L0), which commutes with offsetKernelSensor
+ * and the convert-time fractional shifts. CltExtra offsets are not touched,
+ * and folded kernels are NEVER re-normalized (LoG sum ~ 0).
+ *
+ * Self-test (synthetic impulse/grating validation, roadmap "RT seed" step 1,
+ * run BEFORE any full-system A/B):
+ *   java -cp target/classes:<deps> com.elphel.imagej.cuas.rt.CuasLogFold
+ * Uses only public CPU-oracle methods: ImageDttCPU.clt_2d (MCLT fold+DTT-IV),
+ * convolve_tile (Java twin of the GPU convolveTiles), iclt_2d (inverse +
+ * unfold + overlap-add). Exit code 0 = all tests passed.
+ *
+ * By Claude on 07/06/2026.
+ */
+public class CuasLogFold {
+	public static final int    TRANSFORM_SIZE = 8;
+	public static final int    KERNEL_WIDTH   = 2 * TRANSFORM_SIZE - 1;              // 15
+	public static final double TD_SCALE       = 2 * TRANSFORM_SIZE;                  // 16, see floatGetCltLoGFd / clt_normalize_kernel
+	public static final double LOG_Q_LEAK_MAX = 1e-12; // max allowed |quadrant 1..3| of the (even-even) LoG spectrum
+
+	/**
+	 * Build the per-bin transform-domain response L0[64] of the CUAS LoG
+	 * (CuasRTUtils.getLoGKenel, sigma = psf_radius/sqrt(2), flux-normalized),
+	 * scaled so an identity (delta) filter maps to all-ones (verified by the
+	 * self-test). This is the multiplier for foldKernels().
+	 * @param psf_radius curt.psf_radius (optical PSF radius, pix)
+	 * @param n_sigma    curt.n_sigma (LoG support cutoff, sigmas)
+	 * @return L0[TRANSFORM_SIZE*TRANSFORM_SIZE] per-bin real response
+	 */
+	public static double [] getLogTd(
+			double psf_radius,
+			double n_sigma) {
+		int radius = CuasRTUtils.getRadiusLoGKernel(psf_radius, n_sigma);
+		if (radius > TRANSFORM_SIZE - 1) {
+			throw new IllegalArgumentException("CuasLogFold.getLogTd(): LoG radius "+radius+
+					" exceeds kernel support "+(TRANSFORM_SIZE-1)+" (psf_radius="+psf_radius+", n_sigma="+n_sigma+")");
+		}
+		double [] log_pix = CuasRTUtils.getLoGKenel(psf_radius, radius);
+		// Pre-compensate taps by 1/(cos(pi*dx/16)*cos(pi*dy/16)): a per-bin TD
+		// quad-multiply is equivalent to pixel convolution with cos-ATTENUATED
+		// taps (verified to machine precision by the self-test; same weights as
+		// norm_sym_weights in QuadCLTCPU.calculateCLTKernel:8339). Dividing here
+		// makes the EFFECTIVE applied filter equal the true pixel LoG.
+		int lw = 2 * radius + 1;
+		double [] log_comp = new double [log_pix.length];
+		for (int dy = -radius; dy <= radius; dy++) {
+			for (int dx = -radius; dx <= radius; dx++) {
+				log_comp[(dy + radius) * lw + (dx + radius)] =
+						log_pix[(dy + radius) * lw + (dx + radius)] / cosAtten(dx, dy);
+			}
+		}
+		double [] k15 = embedCentered(log_comp, lw);
+		double [][] quad = symmetrizeDtt3(k15, true);
+		// even-even filter: quadrants 1..3 must vanish (guards the embedding/centering)
+		for (int q = 1; q < 4; q++) {
+			for (int i = 0; i < quad[q].length; i++) {
+				if (Math.abs(quad[q][i]) > LOG_Q_LEAK_MAX) {
+					throw new IllegalStateException("CuasLogFold.getLogTd(): LoG quadrant "+q+
+							" leak "+quad[q][i]+" at bin "+i+" - kernel not centered/even?");
+				}
+			}
+		}
+		return quad[0];
+	}
+
+	/**
+	 * Fold a per-bin real response into CLT convolution kernels: multiply all 4
+	 * quadrants of every tile/sensor kernel by log_td, copy everything else
+	 * (CltExtra offsets/derivatives) verbatim. Input is NOT modified - the
+	 * original kernels stay available for A/B re-upload mid-run.
+	 * @param clt_kernels [ncam][color][tilesY][tilesX][4+extra][64 or extra len]
+	 * @param log_td      per-bin response from getLogTd()
+	 * @return new kernel array of the same shape
+	 */
+	public static double [][][][][][] foldKernels(
+			double [][][][][][] clt_kernels,
+			double []           log_td) {
+		double [][][][][][] folded = new double [clt_kernels.length][][][][][];
+		for (int ncam = 0; ncam < clt_kernels.length; ncam++) {
+			if (clt_kernels[ncam] == null) continue;
+			folded[ncam] = new double [clt_kernels[ncam].length][][][][];
+			for (int col = 0; col < clt_kernels[ncam].length; col++) {
+				if (clt_kernels[ncam][col] == null) continue;
+				folded[ncam][col] = new double [clt_kernels[ncam][col].length][][][];
+				for (int ty = 0; ty < clt_kernels[ncam][col].length; ty++) {
+					if (clt_kernels[ncam][col][ty] == null) continue;
+					folded[ncam][col][ty] = new double [clt_kernels[ncam][col][ty].length][][];
+					for (int tx = 0; tx < clt_kernels[ncam][col][ty].length; tx++) {
+						double [][] ktile = clt_kernels[ncam][col][ty][tx];
+						if (ktile == null) continue;
+						double [][] ftile = new double [ktile.length][];
+						for (int p = 0; p < ktile.length; p++) {
+							ftile[p] = ktile[p].clone();
+							if (p < 4) {
+								for (int i = 0; i < ftile[p].length; i++) {
+									ftile[p][i] *= log_td[i];
+								}
+							}
+						}
+						folded[ncam][col][ty][tx] = ftile;
+					}
+				}
+			}
+		}
+		return folded;
+	}
+
+	/** Inherent per-tap attenuation of a TD per-bin multiply: pixel-domain
+	 *  equivalent kernel taps are scaled by cos(pi*dx/(2*TS))*cos(pi*dy/(2*TS))
+	 *  (the lapped-window/convolution non-commutation term; same weights as
+	 *  norm_sym_weights in QuadCLTCPU.calculateCLTKernel). */
+	public static double cosAtten(int dx, int dy) {
+		return Math.cos(Math.PI * dx / (2 * TRANSFORM_SIZE)) *
+		       Math.cos(Math.PI * dy / (2 * TRANSFORM_SIZE));
+	}
+
+	/** Embed a centered odd-sized pixel kernel into the centered KERNEL_WIDTH^2 support. */
+	public static double [] embedCentered(
+			double [] kernel,
+			int       kwidth) {
+		if (kwidth > KERNEL_WIDTH) {
+			throw new IllegalArgumentException("embedCentered(): kernel "+kwidth+" > "+KERNEL_WIDTH);
+		}
+		double [] k15 = new double [KERNEL_WIDTH * KERNEL_WIDTH];
+		int off = (KERNEL_WIDTH - kwidth) / 2;
+		for (int i = 0; i < kwidth; i++) {
+			System.arraycopy(kernel, i * kwidth, k15, (i + off) * KERNEL_WIDTH + off, kwidth);
+		}
+		return k15;
+	}
+
+	/**
+	 * Centered pixel kernel (KERNEL_WIDTH^2) -> 4 symmetry classes -> DTT-III
+	 * spectra * TD_SCALE. Mirrors ImageDttCPU.clt_symmetrize_kernel +
+	 * clt_dtt3_kernel (package-private there).
+	 * @param aa_math_sign true = mathematically-derived AA projection
+	 *        (+k00 -k01 -k10 +k11); false = bit-copy of the existing
+	 *        clt_symmetrize_kernel quadrant 3 (which textually repeats the AS
+	 *        expression - see the self-test that discriminates the two).
+	 */
+	public static double [][] symmetrizeDtt3(
+			double [] k15,
+			boolean   aa_math_sign) {
+		int ts = TRANSFORM_SIZE;
+		int in_size = KERNEL_WIDTH;
+		int m1 = ts - 1;
+		int center = m1 * in_size + m1;
+		double [][] sym = new double [4][ts * ts];
+		for (int i = 0; i < ts; i++) {
+			for (int j = 0; j < ts; j++) {
+				double kmm = k15[center - i * in_size - j]; // k(-i,-j)
+				double kmp = k15[center - i * in_size + j]; // k(-i,+j)
+				double kpm = k15[center + i * in_size - j]; // k(+i,-j)
+				double kpp = k15[center + i * in_size + j]; // k(+i,+j)
+				sym[0][i * ts + j] = 0.25 * ( kmm + kmp + kpm + kpp);
+				if (j > 0)            sym[1][ i    * ts + j - 1] = 0.25 * (-kmm + kmp - kpm + kpp);
+				if (i > 0)            sym[2][(i-1) * ts + j    ] = 0.25 * (-kmm - kmp + kpm + kpp);
+				if ((i > 0) && (j > 0)) {
+					sym[3][(i-1) * ts + (j-1)] = aa_math_sign?
+							(0.25 * ( kmm - kmp - kpm + kpp)) : // mathematical AA projection
+							(0.25 * (-kmm + kmp - kpm + kpp));  // as in clt_symmetrize_kernel:10531
+				}
+			}
+		}
+		DttRad2 dtt = new DttRad2(ts);
+		for (int q = 0; q < 4; q++) {
+			sym[q] = dtt.dttt_iiie(sym[q], q, ts);
+			for (int i = 0; i < sym[q].length; i++) {
+				sym[q][i] *= TD_SCALE;
+			}
+		}
+		return sym;
+	}
+
+	// ------------------------------------------------------------------
+	// Synthetic self-test below (headless, no GUI classes touched)
+	// ------------------------------------------------------------------
+
+	/** Plain zero-padded pixel-domain convolution with a centered odd-width kernel.
+	 *  convention: out(x,y) = sum k(dy,dx) * in(y - dy, x - dx) (true convolution) */
+	static double [] pixConvolve(
+			double [] img,
+			int       width,
+			double [] kernel,
+			int       kwidth,
+			boolean   mirror) { // mirror=true: correlation orientation in(y+dy, x+dx)
+		int height = img.length / width;
+		int kr = kwidth / 2;
+		double [] out = new double [img.length];
+		for (int y = 0; y < height; y++) {
+			for (int x = 0; x < width; x++) {
+				double s = 0.0;
+				for (int dy = -kr; dy <= kr; dy++) {
+					int yy = mirror? (y + dy) : (y - dy);
+					if ((yy < 0) || (yy >= height)) continue;
+					for (int dx = -kr; dx <= kr; dx++) {
+						int xx = mirror? (x + dx) : (x - dx);
+						if ((xx < 0) || (xx >= width)) continue;
+						s += kernel[(dy + kr) * kwidth + (dx + kr)] * img[yy * width + xx];
+					}
+				}
+				out[y * width + x] = s;
+			}
+		}
+		return out;
+	}
+
+	/** MCLT -> per-tile convolve_tile(ktile) -> IMCLT. Output is (width-8)x(height-8),
+	 *  out(x,y) corresponds to in(x+4, y+4) (iclt_2d half-tile offset). */
+	static double [] cltConvolve(
+			ImageDttCPU imageDtt,
+			double []   img,
+			int         width,
+			double [][] ktile) {
+		final int window_type = 1; // sine (Princen-Bradley / TDAC)
+		double [][][][] clt = imageDtt.clt_2d(
+				img, width, window_type,
+				0, 0,     // shiftX, shiftY
+				-1, -1,   // debug_tileX, debug_tileY (none)
+				0,        // debug_mode
+				8,        // threadsMax
+				0);       // globalDebugLevel
+		for (int ty = 0; ty < clt.length; ty++) {
+			for (int tx = 0; tx < clt[ty].length; tx++) {
+				imageDtt.convolve_tile(clt[ty][tx], ktile, false);
+			}
+		}
+		return imageDtt.iclt_2d(clt, window_type, 15, 0, 8, 0); // all 4 quadrants
+	}
+
+	/** max |a-b| over the interior / max |b| over the interior; a is the (width-8) clt image,
+	 *  b the full-size reference (b sampled at +4,+4). */
+	static double relErrInterior(
+			double [] a, // (width-8) x (height-8)
+			double [] b, // width x height reference
+			int       width,
+			int       margin) {
+		int aw = width - TRANSFORM_SIZE;
+		int ah = a.length / aw;
+		int half = TRANSFORM_SIZE / 2;
+		double max_diff = 0.0, max_ref = 0.0;
+		for (int y = margin; y < ah - margin; y++) {
+			for (int x = margin; x < aw - margin; x++) {
+				double ref = b[(y + half) * width + (x + half)];
+				double d = Math.abs(a[y * aw + x] - ref);
+				if (d > max_diff) max_diff = d;
+				double r = Math.abs(ref);
+				if (r > max_ref) max_ref = r;
+			}
+		}
+		return (max_ref > 0)? (max_diff / max_ref) : max_diff;
+	}
+
+	static double [] makeImpulses(int width, int height) {
+		double [] img = new double [width * height];
+		// smooth background so relative error is meaningful everywhere
+		for (int y = 0; y < height; y++) {
+			for (int x = 0; x < width; x++) {
+				img[y * width + x] = 2.0 + Math.cos(2 * Math.PI * (0.7 * x + 0.4 * y) / 64.0);
+			}
+		}
+		int [][] pts = {{53, 61}, {64, 64}, {58, 71}, {37, 44}, {90, 83}}; // y,x - off tile centers
+		double [] amps = {10.0, -7.0, 5.0, 8.0, -6.0};
+		for (int p = 0; p < pts.length; p++) {
+			img[pts[p][0] * width + pts[p][1]] += amps[p];
+		}
+		return img;
+	}
+
+	static double [] makeGrating(int width, int height, double fx, double fy, double phase) {
+		double [] img = new double [width * height];
+		for (int y = 0; y < height; y++) {
+			for (int x = 0; x < width; x++) {
+				img[y * width + x] = Math.cos(2 * Math.PI * (fx * x + fy * y) + phase);
+			}
+		}
+		return img;
+	}
+
+	static boolean report(String name, double rel_err, double tol) {
+		boolean pass = rel_err <= tol;
+		System.out.println(String.format("%-58s rel_err=%9.3e  tol=%7.1e  %s",
+				name, rel_err, tol, pass? "PASS" : "FAIL"));
+		return pass;
+	}
+
+	public static void main(String[] args) {
+		final int width = 128, height = 128;
+		final int margin = 24;
+		final double tol = 1e-9;      // machine-precision expectations
+		// symmetric-extension wrap residual of per-tile TD convolution
+		// (inherent to the lapped structure, grows with tap offset; measured
+		// 4.5e-3 for a 1-px shift tap, ~1e-2 for spread kernels)
+		final double WRAP_TOL = 1e-2;
+		boolean all_pass = true;
+
+		ImageDttCPU imageDtt = new ImageDttCPU(
+				16, TRANSFORM_SIZE, new ImageDttParameters(),
+				true,  // aux
+				true,  // mono
+				true,  // lwir
+				1.0);  // scale_strengths
+
+		// --- test 0: delta filter -> all-ones TD response (pins the TD_SCALE chain)
+		double [] delta = new double [KERNEL_WIDTH * KERNEL_WIDTH];
+		delta[(KERNEL_WIDTH * KERNEL_WIDTH) / 2] = 1.0;
+		double [][] id_ktile = symmetrizeDtt3(delta, true);
+		double max_dev = 0.0;
+		for (int i = 0; i < id_ktile[0].length; i++) {
+			max_dev = Math.max(max_dev, Math.abs(id_ktile[0][i] - 1.0));
+		}
+		for (int q = 1; q < 4; q++) {
+			for (int i = 0; i < id_ktile[q].length; i++) {
+				max_dev = Math.max(max_dev, Math.abs(id_ktile[q][i]));
+			}
+		}
+		all_pass &= report("0: delta kernel -> identity TD response", max_dev, 1e-12);
+
+		// --- test 1: identity kernel round trip (pins MCLT scaffold + iclt offset)
+		double [] impulses = makeImpulses(width, height);
+		double [] rt = cltConvolve(imageDtt, impulses, width, id_ktile);
+		all_pass &= report("1: identity kernel: MCLT->conv->IMCLT == input",
+				relErrInterior(rt, impulses, width, margin), tol);
+
+		// --- test 1b: single off-center tap. The TD path must equal pixel
+		//     convolution with the cos-ATTENUATED tap (the model this class is
+		//     built on): shift-by-1 scaled by cos(pi/16). Raw (unattenuated)
+		//     comparison is printed as reference - its error IS 1-cos(pi/16).
+		double [] shift1 = new double [KERNEL_WIDTH * KERNEL_WIDTH];
+		shift1[(KERNEL_WIDTH * KERNEL_WIDTH) / 2 + 1] = 1.0; // k(0,+1)
+		{
+			double [][] sh_ktile = symmetrizeDtt3(shift1, true);
+			double [] sh = cltConvolve(imageDtt, impulses, width, sh_ktile);
+			double e_atten = relErrInterior(sh,
+					pixConvolve(impulses, width, attenuate(shift1, KERNEL_WIDTH), KERNEL_WIDTH, false), width, margin);
+			double e_raw = relErrInterior(sh,
+					pixConvolve(impulses, width, shift1, KERNEL_WIDTH, false), width, margin);
+			System.out.println(String.format(
+					"   1b raw-tap err=%9.3e (should be 1-cos(pi/16)=%9.3e)", e_raw, 1.0 - Math.cos(Math.PI/16)));
+			all_pass &= report("1b: single-tap shift == cos-attenuated pixel conv", e_atten, WRAP_TOL);
+		}
+
+		// --- the LoG (curt defaults: psf_radius=1.0, n_sigma=4)
+		final double psf_radius = 1.0, n_sigma = 4.0;
+		double [] log_td = getLogTd(psf_radius, n_sigma);
+		int log_rad = CuasRTUtils.getRadiusLoGKernel(psf_radius, n_sigma);
+		int log_w = 2 * log_rad + 1;
+		double [] log_pix = CuasRTUtils.getLoGKenel(psf_radius, log_rad);
+		System.out.println("LoG: radius="+log_rad+" ("+log_w+"x"+log_w+
+				"), TD response min/max = "+min(log_td)+" / "+max(log_td));
+
+		// --- test 2: folded identity == pixel-domain LoG (impulses)
+		double [][] log_ktile = mulQuad(id_ktile, log_td);
+		double [] folded_img = cltConvolve(imageDtt, impulses, width, log_ktile);
+		double [] oracle_img = pixConvolve(impulses, width, log_pix, log_w, false);
+		all_pass &= report("2: folded identity == pixel LoG (impulses)",
+				relErrInterior(folded_img, oracle_img, width, margin), WRAP_TOL);
+
+		// --- test 3: grating sweep (pins per-bin scale across frequency).
+		//     Gated on error RELATIVE TO THE UNIT GRATING AMPLITUDE: near the
+		//     LoG null the filtered reference is tiny and ref-relative error
+		//     inflates without meaning. Non-bin frequencies are machine-exact
+		//     (the OLA-exact regime); bin frequencies carry the wrap residual.
+		double [][] freqs = {{1.0/32, 0}, {1.0/16, 0}, {0, 1.0/16}, {1.0/16, 1.0/16},
+				{3.0/32, 1.0/32}, {1.0/8, 1.0/8}, {3.0/16, 0}, {7.0/32, 5.0/32}};
+		double worst_grating = 0.0;
+		for (double [] f : freqs) {
+			double [] g = makeGrating(width, height, f[0], f[1], 0.6);
+			double [] gf = cltConvolve(imageDtt, g, width, log_ktile);
+			double [] go = pixConvolve(g, width, log_pix, log_w, false);
+			double e_rel = relErrInterior(gf, go, width, margin);
+			double e_amp = absErrInterior(gf, go, width, margin); // grating amplitude = 1
+			System.out.println(String.format(
+					"   grating fx=%7.4f fy=%7.4f  err/amp=%9.3e  (ref-relative %9.3e)",
+					f[0], f[1], e_amp, e_rel));
+			worst_grating = Math.max(worst_grating, e_amp);
+		}
+		all_pass &= report("3: grating sweep folded == pixel LoG (err/amplitude)", worst_grating, 3e-2);
+
+		// --- test 4: asymmetric kernel (all 4 quadrants nonzero) vs the
+		//     cos-attenuated pixel oracle - machine-precision MODEL test that
+		//     also discriminates the two AA sign variants of clt_symmetrize_kernel
+		double [] asym = new double [KERNEL_WIDTH * KERNEL_WIDTH];
+		int c = (KERNEL_WIDTH * KERNEL_WIDTH) / 2;
+		asym[c + 1 * KERNEL_WIDTH + 2] = 0.6;  // k(+1,+2)
+		asym[c - 2 * KERNEL_WIDTH + 1] = 0.3;  // k(-2,+1)
+		asym[c]                        = 0.4;  // k(0,0)
+		asym[c + 2 * KERNEL_WIDTH - 3] = -0.2; // k(+2,-3)
+		double [] asym_atten = attenuate(asym, KERNEL_WIDTH);
+		double [] best = new double [2]; // best rel_err per AA variant (0=math, 1=code)
+		for (int variant = 0; variant < 2; variant++) {
+			boolean aa_math = (variant == 0);
+			double [][] asym_ktile = symmetrizeDtt3(asym, aa_math);
+			double [] path_a = cltConvolve(imageDtt, impulses, width, asym_ktile);
+			double e_conv = relErrInterior(path_a,
+					pixConvolve(impulses, width, asym_atten, KERNEL_WIDTH, false), width, margin);
+			double e_corr = relErrInterior(path_a,
+					pixConvolve(impulses, width, asym_atten, KERNEL_WIDTH, true), width, margin);
+			best[variant] = Math.min(e_conv, e_corr);
+			System.out.println(String.format(
+					"   4%s AA %s sign: conv=%9.3e  corr=%9.3e  (best orientation: %s)",
+					aa_math? "a":"b", aa_math? "math":"code", e_conv, e_corr,
+					(e_conv <= e_corr)? "convolution" : "correlation"));
+		}
+		// exactly one AA sign variant should reproduce (attenuated) pixel-domain
+		// convolution; which one documents the file/GPU convention (and whether
+		// clt_symmetrize_kernel quadrant 3 is a latent copy-paste bug)
+		all_pass &= report(String.format(
+				"4: asym kernel model test (winner: AA %s sign)",
+				(best[0] <= best[1])? "math" : "code"),
+				Math.min(best[0], best[1]), 2e-2);
+		// --- test 4c: PURE-AA kernel - the sharpest discriminator. The "code"
+		//     variant projects AA content onto the (already-populated) AS class,
+		//     i.e. loses it entirely for a pure-AA kernel: expect O(1) error
+		//     there and wrap-floor error for the math variant.
+		double [] pure_aa = new double [KERNEL_WIDTH * KERNEL_WIDTH];
+		pure_aa[c + 1 * KERNEL_WIDTH + 1] =  0.25;
+		pure_aa[c + 1 * KERNEL_WIDTH - 1] = -0.25;
+		pure_aa[c - 1 * KERNEL_WIDTH + 1] = -0.25;
+		pure_aa[c - 1 * KERNEL_WIDTH - 1] =  0.25;
+		pure_aa[c] = 1.0; // carrier so output is not ~0
+		double [] pure_aa_atten = attenuate(pure_aa, KERNEL_WIDTH);
+		double [] aa_err = new double [2];
+		for (int variant = 0; variant < 2; variant++) {
+			double [][] aa_ktile = symmetrizeDtt3(pure_aa, variant == 0);
+			double [] out_aa = cltConvolve(imageDtt, impulses, width, aa_ktile);
+			aa_err[variant] = Math.min(
+					relErrInterior(out_aa, pixConvolve(impulses, width, pure_aa_atten, KERNEL_WIDTH, false), width, margin),
+					relErrInterior(out_aa, pixConvolve(impulses, width, pure_aa_atten, KERNEL_WIDTH, true),  width, margin));
+			System.out.println(String.format("   4c AA %s sign (pure-AA kernel): rel_err=%9.3e",
+					(variant == 0)? "math":"code", aa_err[variant]));
+		}
+		all_pass &= report("4c: pure-AA kernel, math sign", aa_err[0], 2e-2);
+		if (aa_err[1] < 10 * aa_err[0]) {
+			System.out.println("   WARNING: code AA sign not clearly refuted (expected O(1) error)");
+		}
+
+		// --- test 5: full-run A/B PREDICTION. Production-like aberration kernel
+		//     K (its taps are TD-effective by calibration), folded K*L0 vs the
+		//     oracle pixelLoG(TD-render with K). Per-bin folding composes as
+		//     (k (*) log/c)*c while the oracle is (k*c) (*) log - the difference
+		//     is a second-order cross term this test quantifies (REPORTED, loose
+		//     tolerance; if too big, the exact route is a per-tile pixel-domain
+		//     composite: q = ((k*c) (*) log)/c, host-side at upload).
+		double [] kab = new double [KERNEL_WIDTH * KERNEL_WIDTH];
+		double s_eff = 0.0;
+		for (int dy = -KERNEL_WIDTH/2; dy <= KERNEL_WIDTH/2; dy++) {
+			for (int dx = -KERNEL_WIDTH/2; dx <= KERNEL_WIDTH/2; dx++) {
+				// asymmetric decaying blur: Gaussian sigma 1.2 centered (+0.4,-0.3)
+				double g = Math.exp(-((dx-0.4)*(dx-0.4) + (dy+0.3)*(dy+0.3)) / (2*1.2*1.2));
+				kab[(dy + KERNEL_WIDTH/2) * KERNEL_WIDTH + (dx + KERNEL_WIDTH/2)] = g;
+				s_eff += g * cosAtten(dx, dy); // production normalization convention
+			}
+		}
+		for (int i = 0; i < kab.length; i++) kab[i] /= s_eff;
+		double [][] kab_ktile = symmetrizeDtt3(kab, true);
+		double [] render_k   = cltConvolve(imageDtt, impulses, width, kab_ktile);          // current pipeline
+		double [] render_kl  = cltConvolve(imageDtt, impulses, width, mulQuad(kab_ktile, log_td)); // folded
+		// oracle: pixel LoG applied to the current-pipeline render (interior only;
+		// render_k is already the cropped clt-size image)
+		int aw = width - TRANSFORM_SIZE;
+		double [] oracle_ab = pixConvolve(render_k, aw, log_pix, log_w, false);
+		double max_d = 0.0, max_r = 0.0;
+		for (int y = margin; y < (render_kl.length/aw) - margin; y++) {
+			for (int x = margin; x < aw - margin; x++) {
+				max_d = Math.max(max_d, Math.abs(render_kl[y*aw+x] - oracle_ab[y*aw+x]));
+				max_r = Math.max(max_r, Math.abs(oracle_ab[y*aw+x]));
+			}
+		}
+		double e5 = max_d / max_r;
+		System.out.println(String.format(
+				"   5: fold-vs-pixel-LoG cross-term residual = %9.3e (predicts full-run A/B floor)", e5));
+		all_pass &= report("5: A/B prediction: folded K*L0 vs pixelLoG(render(K))", e5, 1e-1);
+
+		System.out.println(all_pass? "ALL PASSED" : "FAILURES PRESENT");
+		System.exit(all_pass? 0 : 1);
+	}
+
+	/** max |a-b| over the interior (no normalization); a is the clt-cropped image. */
+	static double absErrInterior(
+			double [] a,
+			double [] b,
+			int       width,
+			int       margin) {
+		int aw = width - TRANSFORM_SIZE;
+		int ah = a.length / aw;
+		int half = TRANSFORM_SIZE / 2;
+		double max_diff = 0.0;
+		for (int y = margin; y < ah - margin; y++) {
+			for (int x = margin; x < aw - margin; x++) {
+				max_diff = Math.max(max_diff, Math.abs(a[y * aw + x] - b[(y + half) * width + (x + half)]));
+			}
+		}
+		return max_diff;
+	}
+
+	/** Copy of a centered odd-width tap kernel with each tap scaled by cosAtten. */
+	static double [] attenuate(double [] k, int kwidth) {
+		int kr = kwidth / 2;
+		double [] out = new double [k.length];
+		for (int dy = -kr; dy <= kr; dy++) {
+			for (int dx = -kr; dx <= kr; dx++) {
+				out[(dy + kr) * kwidth + (dx + kr)] = k[(dy + kr) * kwidth + (dx + kr)] * cosAtten(dx, dy);
+			}
+		}
+		return out;
+	}
+
+	static double [][] mulQuad(double [][] ktile, double [] l0) {
+		double [][] out = new double [ktile.length][];
+		for (int q = 0; q < ktile.length; q++) {
+			out[q] = ktile[q].clone();
+			if (q < 4) {
+				for (int i = 0; i < out[q].length; i++) {
+					out[q][i] *= l0[i];
+				}
+			}
+		}
+		return out;
+	}
+
+	static double min(double [] a) { double m = a[0]; for (double v : a) m = Math.min(m, v); return m; }
+	static double max(double [] a) { double m = a[0]; for (double v : a) m = Math.max(m, v); return m; }
+}