/**
 **
 ** CuasTD.java - transform-domain utilities for the CUAS RT path:
 ** cross-sensor consolidation of per-sensor CLT data (CPU oracle for the
 ** future GPU kernel).
 **
 ** Copyright (C) 2026 Elphel, Inc.
 **
 ** -----------------------------------------------------------------------------**
 **
 **  CuasTD.java is free software: you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License as published by
 **  the Free Software Foundation, either version 3 of the License, or
 **  (at your option) any later version.
 **
 **  This program is distributed in the hope that it will be useful,
 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 **  GNU General Public License for more details.
 **
 **  You should have received a copy of the GNU General Public License
 **  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ** -----------------------------------------------------------------------------**
 **
 */
package com.elphel.imagej.cuas.rt;

import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.cuas.CuasMotion;
import com.elphel.imagej.gpu.GPUTileProcessor;
import com.elphel.imagej.gpu.GpuQuad;
import com.elphel.imagej.gpu.TpTask; // By Claude on 07/12/2026
import com.elphel.imagej.tileprocessor.QuadCLT;

import ij.ImagePlus;

/**
 * Transform-domain utilities for the CUAS RT pose/detection path.
 *
 * The central operation is the NaN-aware cross-sensor consolidation of the 16
 * per-sensor CLT channels into ONE averaged TD channel - the "average images
 * BEFORE correlation, multiply averages instead of averaging products" scheme
 * (better SNR: sensor-noise cross-terms drop before the nonlinear FZ-normalize,
 * and 1 correlation replaces 16). Not available on the GPU (combine_inter only
 * sums per-sensor correlation PRODUCTS, unweighted); this CPU implementation +
 * the GpuQuad get/setCltData D2H/H2D bridge is the phase-A2 path and the bit
 * oracle for the future clt_average_sensors CUDA kernel.
 *
 * TD granularity note: validity is per TILE per sensor (a sensor either has the
 * converted tile or it was erased to NaN / never written) - so the consolidation
 * unit is one [4][8][8] CLT tile (per color), and the weight is a per-tile
 * contributing-sensor COUNT, not per element.
 *
 * By Claude on 07/05/2026, from Andrey's design.
 */
public class CuasTD {
	public static final int TD_CHUNK = 4 * GPUTileProcessor.DTT_SIZE * GPUTileProcessor.DTT_SIZE; // one CLT tile (per color)
	// OOB source-coordinate priority margins (roadmap 2c, frozen spec Andrey 07/12/2026):
	// depth = distance of the per-sensor tile center to the nearest source-frame edge.
	// SOFT ~12 px = MCLT 24x24 support overhang (defined-but-synthetic replicated rows) ->
	// deprioritize: use only when no sensor has the tile deeper. HARD ~8 px = 16x16 core
	// overhang -> always skip. 0 = disabled (v0-compatible). By Claude on 07/12/2026.
	public static final double OOB_SOFT_DEFAULT = 12.0;
	public static final double OOB_HARD_DEFAULT =  8.0;

	/**
	 * CPU oracle of the consolidation TILE/SENSOR selection - the contract of the
	 * future index_consolidate (pass 1) + OOB priority (pass 2) GPU kernels
	 * (roadmap 2c). Inputs mirror the kernel: TWO task sets (motion-blur pair,
	 * tp_tasks[1] may be null = no-MB mode - same convention as
	 * GpuQuad.setInterTasksMotionBlur()), sets index-matched by tile txy (looked
	 * up, not by array position). Per tile:
	 *   pass 1: the tile must be present (task != 0) in EVERY provided set,
	 *           else mask = 0 (tile dropped);
	 *   pass 2: per sensor, depth = min over provided sets of the distance of
	 *           the sensor's tile center (TpTask.xy) to the nearest frame edge;
	 *           depth < hard_margin -> sensor always skipped; if any sensor has
	 *           depth >= soft_margin only those are kept, else the hard
	 *           survivors serve as the last resort. NaN/undefined xy -> sensor
	 *           excluded (when filtering is enabled). Margin 0 = that tier off;
	 *           both 0 -> all-sensor mask (v0-compatible).
	 * By Claude on 07/12/2026, from Andrey's design.
	 *
	 * @param tp_tasks    [2][]: task set 0 (required) + MB pair set or null
	 * @param num_sens    number of sensors (mask width, LSB = sensor 0)
	 * @param tilesX      tile grid width  (tasks may be full-grid or compacted)
	 * @param tilesY      tile grid height
	 * @param width       source sensor frame width, px
	 * @param height      source sensor frame height, px
	 * @param soft_margin ~12.0 (support overhang), 0 = disabled
	 * @param hard_margin ~8.0 (core overhang), 0 = disabled
	 * @return int[tilesY*tilesX]: bit s set = average sensor s; 0 = tile poisoned
	 */
	public static int [] oobSensorMasks(
			final TpTask [][] tp_tasks,
			final int         num_sens,
			final int         tilesX,
			final int         tilesY,
			final int         width,
			final int         height,
			final double      soft_margin,
			final double      hard_margin) {
		final int num_tiles = tilesY * tilesX;
		final TpTask [][] by_nt = new TpTask [2][];
		for (int nset = 0; nset < 2; nset++) {
			if ((tp_tasks == null) || (nset >= tp_tasks.length) || (tp_tasks[nset] == null)) continue;
			by_nt[nset] = new TpTask [num_tiles];
			for (TpTask t : tp_tasks[nset]) if (t != null) {
				by_nt[nset][t.getTileY() * tilesX + t.getTileX()] = t;
			}
		}
		final boolean filter = (soft_margin > 0) || (hard_margin > 0);
		final int all_mask = (num_sens >= 32) ? -1 : ((1 << num_sens) - 1);
		final int [] masks = new int [num_tiles];
		int no_xy_tiles = 0; // tiles that could not be geometry-filtered (missing xy)
		for (int nt = 0; nt < num_tiles; nt++) {
			boolean present = (by_nt[0] != null) || (by_nt[1] != null);
			for (int nset = 0; nset < 2; nset++) if (by_nt[nset] != null) {
				final TpTask t = by_nt[nset][nt];
				present &= (t != null) && (t.getTask() != 0);
			}
			if (!present) continue; // mask stays 0 - tile dropped (pass-1 semantics)
			if (!filter) {
				masks[nt] = all_mask;
				continue;
			}
			boolean have_xy = true;
			for (int nset = 0; nset < 2; nset++) if ((by_nt[nset] != null) && (by_nt[nset][nt].getXY() == null)) {
				have_xy = false;
			}
			if (!have_xy) { // no geometry - cannot filter, pass all (warn once below)
				no_xy_tiles++;
				masks[nt] = all_mask;
				continue;
			}
			int soft_mask = 0, hard_mask = 0;
			for (int s = 0; s < num_sens; s++) {
				double depth = Double.POSITIVE_INFINITY; // min over the provided sets: OOB in EITHER MB pass = OOB
				for (int nset = 0; nset < 2; nset++) if (by_nt[nset] != null) {
					final float [] xy = by_nt[nset][nt].getXY()[s];
					final double d = Math.min(Math.min(xy[0], xy[1]),
							Math.min((width - 1) - xy[0], (height - 1) - xy[1]));
					depth = Math.min(depth, Double.isNaN(d) ? Double.NEGATIVE_INFINITY : d);
				}
				if ((hard_margin <= 0) || (depth >= hard_margin)) {
					hard_mask |= 1 << s;
					if ((soft_margin > 0) && (depth >= soft_margin)) {
						soft_mask |= 1 << s;
					}
				}
			}
			// priority, not rejection: deep sensors when any exist, else the soft-zone
			// survivors as the last resort; nobody above hard -> 0 = poisoned tile
			masks[nt] = ((soft_margin > 0) && (soft_mask != 0)) ? soft_mask : hard_mask;
		}
		if (no_xy_tiles > 0) {
			System.out.println("CuasTD.oobSensorMasks(): WARNING - "+no_xy_tiles+
					" tiles have no per-sensor xy (geometry filter skipped for them)");
		}
		return masks;
	}

	/**
	 * Consolidate per-sensor TD (CLT) data into a single averaged channel:
	 * per tile - add the tiles of the sensors that HAVE it (skip NaN tiles),
	 * count them, divide by the count. A tile is considered absent for a sensor
	 * if its first element is NaN (erase_clt_tiles fills whole tiles); any stray
	 * NaN inside an otherwise-valid tile poisons the sum and the result tile is
	 * set to all-NaN (with the count kept) - fail-visible, never fail-silent.
	 *
	 * @param fclt   per-sensor CLT as returned by GpuQuad.getCltData():
	 *               [num_sensors][num_tiles*num_colors*TD_CHUNK]
	 * @param counts null or int[num_tiles*num_colors] - filled with the number of
	 *               contributing sensors per tile (the weight plane)
	 * @return averaged CLT channel, same length as fclt[0]; all-NaN tiles where
	 *               no sensor contributed
	 */
	public static float [] consolidateSensorsTD(
			final float [][] fclt,
			final int []     counts) {
		return consolidateSensorsTD(fclt, counts, null, 1); // v1 with no mask == bit-exact v0 // By Claude on 07/12/2026
	}

	/**
	 * v1 (roadmap 2c): same NaN-aware consolidation gated by per-TILE sensor masks
	 * (see {@link #oobSensorMasks}): only sensors with their mask bit set are
	 * accumulated; mask 0 = whole tile poisoned (all-NaN, count 0 - the v0
	 * whole-tile-poison semantics, never a synthetic average). sens_masks == null
	 * reproduces v0 bit-exactly (identical op order). By Claude on 07/12/2026.
	 *
	 * @param sens_masks null, or int[tilesY*tilesX] per-tile sensor mask
	 * @param colors     chunk-to-tile mapping (chunk = tile*colors + color); 1 for LWIR
	 */
	public static float [] consolidateSensorsTD(
			final float [][] fclt,
			final int []     counts,
			final int []     sens_masks,
			final int        colors) {
		final int len = fclt[0].length;
		final int num_chunks = len / TD_CHUNK;
		final float [] avg = new float [len];
		for (int chunk = 0; chunk < num_chunks; chunk++) {
			final int offs = chunk * TD_CHUNK;
			final int mask = (sens_masks != null) ? sens_masks[chunk / colors] : -1; // By Claude on 07/12/2026
			int n = 0;
			boolean first = true;
			for (int nsens = 0; nsens < fclt.length; nsens++) {
				if (((mask >> nsens) & 1) == 0) continue; // sensor OOB-excluded / tile dropped // By Claude on 07/12/2026
				final float [] f = fclt[nsens];
				if (Float.isNaN(f[offs])) continue; // sensor does not have this tile
				n++;
				if (first) {
					System.arraycopy(f, offs, avg, offs, TD_CHUNK);
					first = false;
				} else {
					for (int i = 0; i < TD_CHUNK; i++) avg[offs + i] += f[offs + i];
				}
			}
			if (n == 0) {
				java.util.Arrays.fill(avg, offs, offs + TD_CHUNK, Float.NaN);
			} else if (n > 1) {
				final float scale = 1.0f / n;
				boolean poisoned = false;
				for (int i = 0; i < TD_CHUNK; i++) {
					avg[offs + i] *= scale;
					poisoned |= Float.isNaN(avg[offs + i]);
				}
				if (poisoned) { // stray NaN inside a "valid" tile - make the whole tile NaN
					java.util.Arrays.fill(avg, offs, offs + TD_CHUNK, Float.NaN);
				}
			}
			if (counts != null) counts[chunk] = n;
		}
		return avg;
	}

	private static boolean gpu_bridge_reported = false; // one-line console note on the first GPU-chain run // By Claude on 07/12/2026

	/**
	 * The consolidation bridge shared by the RT call sites (lean pose measurement,
	 * virtual render "merged"): consolidate the 16-sensor TD into the cam-0 CLT slot.
	 * GPU-first: {@link GpuQuad#execConsolidateSensorsTD} runs the on-device v1 chain
	 * (index_consolidate -> consolidate_oob -> clt_average_sensors_list, JNA backend
	 * only - the per-sensor TD never leaves the device). CPU fallback (the pre-07/12
	 * bridge, now OOB-mask-capable): getCltData D2H -> oobSensorMasks + masked
	 * consolidateSensorsTD -> setCltData(0) H2D - used on the JCuda backend (Andrey's
	 * 07/12/2026 ruling: no JCuda launch code for new kernels) or if the native call
	 * fails. Margins 0/0 = v0 unfiltered (also when tp_tasks == null).
	 * By Claude on 07/12/2026.
	 *
	 * @param gpuQuad     GPU with a current per-sensor TD (after convert_direct)
	 * @param tp_tasks    [2][]: task set 0 + MB pair set ([1] may be null = no-MB),
	 *                    the sets used for the convert; null = no geometry filtering
	 * @param soft_margin OOB soft margin, px (curt.oob_soft; 0 = tier disabled)
	 * @param hard_margin OOB hard margin, px (curt.oob_hard; 0 = tier disabled)
	 */
	public static void consolidateToSlot0(
			final GpuQuad     gpuQuad,
			final TpTask [][] tp_tasks,
			final double      soft_margin,
			final double      hard_margin) {
		if (gpuQuad.execConsolidateSensorsTD(tp_tasks, soft_margin, hard_margin)) {
			if (!gpu_bridge_reported) {
				gpu_bridge_reported = true;
				System.out.println("CuasTD.consolidateToSlot0(): GPU consolidation chain active (soft/hard = "+
						soft_margin+"/"+hard_margin+" px, no D2H/H2D)");
			}
			return;
		}
		final float [][] fclt16 = gpuQuad.getCltData(false);
		int [] masks = null;
		if ((tp_tasks != null) && ((soft_margin > 0) || (hard_margin > 0))) {
			final int [] wh = gpuQuad.getWH(false);
			final int tilesX = wh[0] / GPUTileProcessor.DTT_SIZE;
			final int tilesY = wh[1] / GPUTileProcessor.DTT_SIZE;
			masks = oobSensorMasks(tp_tasks, gpuQuad.num_cams, tilesX, tilesY,
					wh[0], wh[1], soft_margin, hard_margin);
		}
		final float [] avg_td = consolidateSensorsTD(fclt16, null, masks, gpuQuad.num_colors);
		gpuQuad.setCltData(0, avg_td, false);
	}

	/**
	 * Validate the consolidation against the linearity oracle: imclt is linear, so
	 * imclt(TD-average of sensors) must equal the pixel-domain average of the
	 * per-sensor imclt renders (up to float summation order). Renders both through
	 * the SAME GPU imclt, compares, prints max|diff| / RMS, optionally saves a
	 * 3-slice debug stack, and restores the original per-sensor TD.
	 *
	 * @param gpuQuad   GPU with a current 16-sensor TD (e.g. after convert_direct)
	 * @param save_qclt null or QuadCLT whose model directory receives the debug
	 *                  stack (-CUAS-TDAVG-CHECK: [pixel-avg, imclt(TD-avg), diff])
	 * @param debugLevel debug level
	 * @return max abs difference (finite pixels); expect float-epsilon scale
	 */
	public static double validateConsolidation(
			final GpuQuad gpuQuad,
			final QuadCLT save_qclt,
			final int     debugLevel) {
		final int num_sens = gpuQuad.getNumSensors();
		final float [][] orig = gpuQuad.getCltData(false);
		final float [][] imgs = CuasMotion.perSensorImagesFromTD(gpuQuad, false); // renders BEFORE overwrite
		final int width  = gpuQuad.getImageWidth();
		final int height = gpuQuad.getImageHeight();
		final int [] counts = new int [orig[0].length / TD_CHUNK];
		final float [] avg_td = consolidateSensorsTD(orig, counts);
		int cmin = Integer.MAX_VALUE, cmax = 0; long csum = 0; int cnz = 0;
		for (int c : counts) if (c > 0) { cmin = Math.min(cmin, c); cmax = Math.max(cmax, c); csum += c; cnz++; }
		System.out.println("CuasTD.validateConsolidation(): contributing sensors per tile: min="+
				(cnz > 0 ? cmin : 0)+", max="+cmax+", mean="+String.format("%.2f", (cnz > 0)? ((double) csum)/cnz : 0.0)+
				", tiles with data="+cnz+" of "+counts.length);
		// upload the averaged channel to ALL sensor slots, render, take sensor 0
		for (int ncam = 0; ncam < num_sens; ncam++) gpuQuad.setCltData(ncam, avg_td, false);
		final float [][] imgs_avg = CuasMotion.perSensorImagesFromTD(gpuQuad, false);
		final float [] render_of_avg = imgs_avg[0];
		// pixel-domain NaN-aware average of the per-sensor renders
		final float [] pix_avg = new float [width * height];
		final float [] diff =    new float [width * height];
		double max_diff = 0, sum2 = 0; int n_cmp = 0;
		for (int i = 0; i < pix_avg.length; i++) {
			double s = 0; int n = 0;
			for (int ncam = 0; ncam < num_sens; ncam++) {
				final float v = imgs[ncam][i];
				if (!Float.isNaN(v)) { s += v; n++; }
			}
			pix_avg[i] = (n > 0) ? (float) (s / n) : Float.NaN;
			diff[i] = pix_avg[i] - render_of_avg[i];
			if (!Float.isNaN(diff[i])) {
				final double ad = Math.abs(diff[i]);
				if (ad > max_diff) max_diff = ad;
				sum2 += diff[i] * diff[i];
				n_cmp++;
			}
		}
		System.out.println(String.format(
				"CuasTD.validateConsolidation(): imclt(TD-avg) vs pixel-avg(per-sensor imclt): max|diff|=%.6g, RMS=%.6g over %d pixels",
				max_diff, (n_cmp > 0) ? Math.sqrt(sum2 / n_cmp) : Double.NaN, n_cmp));
		// restore the original per-sensor TD
		for (int ncam = 0; ncam < num_sens; ncam++) gpuQuad.setCltData(ncam, orig[ncam], false);
		if (save_qclt != null) {
			final ImagePlus imp = ShowDoubleFloatArrays.makeArrays(
					new float [][] {pix_avg, render_of_avg, diff}, width, height,
					save_qclt.getImageName()+"-CUAS-TDAVG-CHECK",
					new String [] {"pixel_avg", "imclt_td_avg", "diff"});
			save_qclt.saveImagePlusInModelDirectory(null, imp);
		}
		return max_diff;
	}

	// ===================== OOB oracle self-test (headless, no GPU) =====================
	// java -cp target/classes:<deps> com.elphel.imagej.cuas.rt.CuasTD
	// Synthetic validation of oobSensorMasks() + masked consolidateSensorsTD() against
	// the frozen 2c spec BEFORE any GPU run (the CuasLogFold pattern). Exit 0 = PASS.
	// By Claude on 07/12/2026.
	private static int test_fails = 0;
	private static void check(boolean ok, String what) {
		System.out.println((ok ? "PASS: " : "FAIL: ") + what);
		if (!ok) test_fails++;
	}
	/** task with all sensors' tile centers at the given frame-edge depths */
	private static TpTask testTask(int num_sens, int tx, int ty, int width, int height, double [] depths) {
		final TpTask t = new TpTask(num_sens, tx, ty);
		t.task = 1; // any nonzero = active (the practical GPU contract)
		final float [][] xy = new float [num_sens][2];
		for (int s = 0; s < num_sens; s++) { // depth from the LEFT edge, y centered (depth-neutral)
			xy[s][0] = (float) depths[s];
			xy[s][1] = height / 2.0f;
		}
		t.xy = xy;
		return t;
	}
	public static void main(String [] args) {
		final int NS = 4, TX = 4, TY = 2, W = 640, H = 512;
		final double SOFT = OOB_SOFT_DEFAULT, HARD = OOB_HARD_DEFAULT;
		final int ALL = (1 << NS) - 1;
		// grid tiles: 0 all-deep | 1 mixed | 2 all-soft-zone | 3 all-below-hard |
		//             4 missing in MB set | 5 task==0 in set 0 | 6 NaN xy sensor | 7 absent everywhere
		final TpTask [] set0 = new TpTask [7];
		final TpTask [] set1 = new TpTask [6]; // no entries for tiles 4 (dropped) and 7
		set0[0] = testTask(NS, 0, 0, W, H, new double [] {100, 200, 300, 400});
		set0[1] = testTask(NS, 1, 0, W, H, new double [] {100,  50,  10,   5}); // deep,deep,soft,below-hard
		set0[2] = testTask(NS, 2, 0, W, H, new double [] {  9,  10,  11,  11.9}); // all soft-zone
		set0[3] = testTask(NS, 3, 0, W, H, new double [] {  3,   5,   0,   7.9}); // all below hard
		set0[4] = testTask(NS, 0, 1, W, H, new double [] {100, 100, 100, 100});
		set0[5] = testTask(NS, 1, 1, W, H, new double [] {100, 100, 100, 100});
		set0[5].task = 0; // inactive in set 0
		set0[6] = testTask(NS, 2, 1, W, H, new double [] {100, 100, 100, 100});
		set0[6].getXY()[1][0] = Float.NaN; // sensor 1 projection undefined
		for (int i = 0; i < 4; i++) { // MB pair for tiles 0-3: same geometry shifted 1 px deeper
			final TpTask s = set0[i];
			final double [] d = new double [NS];
			for (int ss = 0; ss < NS; ss++) d[ss] = s.getXY()[ss][0] + 1;
			set1[i] = testTask(NS, s.getTileX(), s.getTileY(), W, H, d);
		}
		// NO set-1 entry for tile 4 (0,1) - the MB-drop case
		set1[4] = testTask(NS, 1, 1, W, H, new double [] {100, 100, 100, 100}); // tile 5: active here, task==0 in set 0
		set1[5] = testTask(NS, 2, 1, W, H, new double [] {100, 100, 100, 100}); // tile 6: MB entry with clean xy (NaN only in set 0)
		// --- masks, MB mode ---
		int [] m = oobSensorMasks(new TpTask [][] {set0, set1}, NS, TX, TY, W, H, SOFT, HARD);
		check(m[0] == ALL,            "tile 0: all sensors deep -> full mask (got "+m[0]+")");
		check(m[1] == 0b0011,         "tile 1: soft tier wins, soft-zone + below-hard dropped (got "+m[1]+")");
		check(m[2] == ALL,            "tile 2: nobody deep -> soft-zone last resort keeps all (got "+m[2]+")");
		check(m[3] == 0,              "tile 3: all below hard -> poisoned (got "+m[3]+")");
		check(m[4] == 0,              "tile 4: missing in MB set -> dropped (got "+m[4]+")");
		check(m[5] == 0,              "tile 5: task==0 in set 0 -> dropped (got "+m[5]+")");
		check(m[6] == (ALL & ~0b10),  "tile 6: NaN xy sensor excluded (got "+m[6]+")");
		check(m[7] == 0,              "tile 7: absent everywhere -> 0 (got "+m[7]+")");
		// --- no-MB mode (set 1 null): tile 4 comes back, tile 5 still inactive ---
		m = oobSensorMasks(new TpTask [][] {set0, null}, NS, TX, TY, W, H, SOFT, HARD);
		check(m[4] == ALL,            "no-MB: tile 4 present again (got "+m[4]+")");
		check(m[5] == 0,              "no-MB: tile 5 still task==0 (got "+m[5]+")");
		// --- margins disabled = v0-compatible: every present tile gets the full mask ---
		m = oobSensorMasks(new TpTask [][] {set0, set1}, NS, TX, TY, W, H, 0, 0);
		check((m[1] == ALL) && (m[2] == ALL) && (m[3] == ALL), "margins 0: no geometry filtering");
		check(m[4] == 0,              "margins 0: pass-1 MB drop still applies (got "+m[4]+")");
		// --- masked consolidation: sensor s constant s+1 per present tile ---
		final int num_tiles = TX * TY;
		final float [][] fclt = new float [NS][num_tiles * TD_CHUNK];
		for (int s = 0; s < NS; s++) java.util.Arrays.fill(fclt[s], s + 1.0f);
		java.util.Arrays.fill(fclt[2], 5 * TD_CHUNK, 6 * TD_CHUNK, Float.NaN); // sensor 2 absent on tile 5
		final int [] masks = new int [] {ALL, 0b0011, ALL, 0, ALL, 0b0100, ALL, 0};
		final int [] counts = new int [num_tiles];
		final float [] avg = consolidateSensorsTD(fclt, counts, masks, 1);
		check(Math.abs(avg[0]          - 2.5f) < 1e-6, "avg tile 0 = mean(1..4) = 2.5 (got "+avg[0]+")");
		check(Math.abs(avg[TD_CHUNK]   - 1.5f) < 1e-6, "avg tile 1 = mean(1,2) = 1.5 (got "+avg[TD_CHUNK]+")");
		check(Float.isNaN(avg[3 * TD_CHUNK]) && (counts[3] == 0), "tile 3 poisoned NaN, count 0");
		check(Float.isNaN(avg[5 * TD_CHUNK]) && (counts[5] == 0), "tile 5: only allowed sensor absent -> NaN, count 0");
		check(counts[0] == 4 && counts[1] == 2, "counts follow allowed+present sensors");
		// --- null masks == v0 bit-exact ---
		final float [] v0 = consolidateSensorsTD(fclt, null);
		final float [] v1 = consolidateSensorsTD(fclt, null, null, 1);
		boolean same = true;
		for (int i = 0; i < v0.length; i++) same &= (v0[i] == v1[i]) || (Float.isNaN(v0[i]) && Float.isNaN(v1[i]));
		check(same, "null masks reproduce v0 bit-exactly");
		System.out.println(test_fails == 0 ? "CuasTD OOB oracle self-test: ALL PASSED"
				: "CuasTD OOB oracle self-test: " + test_fails + " FAILED");
		System.exit((test_fails == 0) ? 0 : 1);
	}
}
