CLAUDE: test_pose_corr - lean pose checkpoint-1 replay (roadmap 3-C rung 1)

Standalone per-iteration replay of the CuasPoseRT.leanMeasure GPU chain
against the PoseCorrExport case (imagej-elphel curt.kernel_test=pose_corr):
calcReverseDistortionTable + calc_rot_deriv + calculate_tiles_offsets (NEW
under test; GPU-filled tasks compared bit-exact vs the oracle post-offsets
streams, header+xy fields - the disp_dist tail is not round-tripped by the
Java flatten), erase_clt_tiles(NaN) + convert_direct for both MB sets,
the v1 consolidation chain into the cam-0 slot (inputs = the oracle
post-offsets streams, the exact bytes production flattens), then
correlate2D_inter vs the resident center TD from gpu_clt_ref cam 0 (NEW
under test). Raw-TD compare is ORDER-INDEPENDENT (corr slot placement is
atomicAdd-nondeterministic): tiles matched by packed index, tol 0 goal;
results/corr_td_it<k> saved in oracle order for the Java round trip. All
four LPF/HPF __constant__ tables overwritten from the case via extern
__constant__ + cudaMemcpyToSymbol (correlate2D_inter reads LoG_corr and
lpf_rb_corr). Built by build_tests.sh (auto-discovered), cuda-gdb steppable.
Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

src/tests/test_pose_corr.cu

+/*
+ * test_pose_corr.cu
+ *
+ * Standalone test of the lean pose-chain CHECKPOINT 1 (roadmap 3-C rung 1):
+ * replay the per-LMA-iteration GPU chain of CuasPoseRT.leanMeasure for one
+ * real scene against the Java-oracle case exported by PoseCorrExport
+ * (curt.kernel_test = pose_corr):
+ *
+ *   per iteration (motion-blur pair):
+ *     set 0: upload tasks (pre-offsets) -> calcReverseDistortionTable +
+ *            calc_rot_deriv + calculate_tiles_offsets  [NEWLY EXERCISED]
+ *            -> compare GPU-filled tasks vs oracle -> erase_clt_tiles(NaN)
+ *            -> convert_direct (SET, positive scale)
+ *     set 1: same, no erase (SUBTRACT, negative scale in the task data)
+ *     consolidation: index_consolidate -> consolidate_oob ->
+ *            clt_average_sensors_list -> average into the cam-0 CLT slot
+ *            (inputs = the oracle post-offsets streams, the exact bytes the
+ *            production bridge flattens - the offsets check above already
+ *            validated that this test's own streams match them)
+ *     correlate2D_inter (sel_sensors=1) vs the resident center TD
+ *            (gpu_clt_ref cam 0, uploaded from the case) [NEWLY EXERCISED]
+ *     -> compare raw TD correlation tiles vs the oracle, ORDER-INDEPENDENT:
+ *        the corr slot placement comes from an atomicAdd (nondeterministic
+ *        run-to-run), so tiles are matched by their packed index
+ *        (tile << CORR_NTILE_SHIFT | sensor, 0xff = sum) and the data
+ *        compared at tol (0 = bit-exact goal). results/corr_td_it<k> are
+ *        saved re-ordered into the ORACLE index order, so the Java round
+ *        trip (PoseCorrExport.compareResults) compares element-wise.
+ *
+ * Everything is data-driven from the manifest (testdata_format.md); all four
+ * LPF/HPF __constant__ arrays are overwritten from the case (correlate2D_inter
+ * reads LoG_corr and lpf_rb_corr). Deterministic single pass - set breakpoints
+ * in calculate_tiles_offsets / correlate2D_inter_inner and step in
+ * NSight/cuda-gdb.
+ *
+ * Usage: test_pose_corr --data <dir> [--tol X] [--list]
+ *
+ *  Created on: Jul 13, 2026
+ *      Author: Claude (Anthropic), for Elphel
+ */
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+#include <map>
+#include <string>
+#include <vector>
+
+#include <cuda_runtime.h>
+#include <helper_cuda.h>           // checkCudaErrors, findCudaDevice
+
+#include "tp_defines.h"            // TASK_* bits, RBYRDIST_LEN, CORR_NTILE_SHIFT
+#include "dtt8x8.h"
+#include "geometry_correction.h"   // struct gc/corr_vector/trot_deriv, geometry kernels, TP_TASK_*_OFFSET
+#include "TileProcessor.h"         // convert_direct, correlate2D_inter, erase_clt_tiles
+#include "tp_consolidate.h"        // index_consolidate, consolidate_oob, clt_average_sensors_list
+#include "tests/tp_test_data.h"
+#include "tests/tp_gpu_buf.h"
+
+// LPF/HPF filter tables live as __constant__ definitions in TileProcessor.cu
+// (external linkage under -rdc=true); the case carries the exact payloads the
+// Java chain uploads, so the baked-in defaults never leak into the replay.
+extern __constant__ float lpf_data[4][64];
+extern __constant__ float lpf_corr[64];
+extern __constant__ float lpf_rb_corr[64];
+extern __constant__ float LoG_corr[64];
+
+// bit-exact (NaN==NaN) compare of the task fields the GPU/oracle fill:
+// 6 header floats + 2*num_cams per-sensor xy; the disp_dist tail is never
+// round-tripped by the Java flatten (asFloatArray) and is excluded.
+static int compareTasks(
+		const float * got,
+		const float * expected,
+		int           n_entries,
+		int           task_size,
+		int           num_cams,
+		const char  * label)
+{
+	const int cmp_floats = 6 + 2 * num_cams;
+	int bad = 0;
+	long first_bad = -1;
+	double maxd = 0;
+	for (int i = 0; i < n_entries; i++) {
+		for (int j = 0; j < cmp_floats; j++) {
+			const size_t idx = (size_t) i * task_size + j;
+			const float g = got[idx], e = expected[idx];
+			const bool gn = std::isnan(g), en = std::isnan(e);
+			if (gn && en) continue;
+			if (gn != en) { bad++; if (first_bad < 0) first_bad = (long) idx; continue; }
+			const double d = std::fabs((double) g - (double) e);
+			if (d > 0) {
+				bad++;
+				if (d > maxd) maxd = d;
+				if (first_bad < 0) first_bad = (long) idx;
+			}
+		}
+	}
+	if (bad) {
+		printf("compare '%s': %d of %d entries x %d floats mismatch, max|diff|=%g, first @%ld -> FAIL\n",
+				label, bad, n_entries, cmp_floats, maxd, first_bad);
+	} else {
+		printf("compare '%s': %d entries, bit-exact -> PASS\n", label, n_entries);
+	}
+	return bad;
+}
+
+int main(int argc, char **argv)
+{
+	std::string data_dir = "testdata/pose_corr";
+	float tol       = 0.0f;   // the whole chain is expected bit-exact vs the oracle
+	int   list_only = 0;
+	for (int i = 1; i < argc; i++) {
+		if      (!strcmp(argv[i], "--data") && (i + 1 < argc)) data_dir = argv[++i];
+		else if (!strcmp(argv[i], "--tol")  && (i + 1 < argc)) tol = (float) atof(argv[++i]);
+		else if (!strcmp(argv[i], "--list"))                   list_only = 1;
+		else {
+			fprintf(stderr, "Usage: %s --data <dir> [--tol X] [--list]\n", argv[0]);
+			return EXIT_FAILURE;
+		}
+	}
+	printf("%s: data=%s tol=%g\n", argv[0], data_dir.c_str(), tol);
+
+	findCudaDevice(argc, (const char **) argv);
+	TpTestData data(data_dir);
+	if (list_only) {
+		data.list();
+		return EXIT_SUCCESS;
+	}
+
+	const int num_cams     = data.paramInt("num_sensors",  16);
+	const int num_colors   = data.paramInt("colors",        1);
+	const int tilesx       = data.paramInt("tilesx",       80);
+	const int tilesy       = data.paramInt("tilesy",       64);
+	const int img_width    = data.paramInt("img_width",   640);
+	const int img_height   = data.paramInt("img_height",  512);
+	const int kernels_hor  = data.paramInt("kernels_hor",  82);
+	const int kernels_vert = data.paramInt("kernels_vert", 66);
+	const int task_size    = data.paramInt("task_size", 6 + 6 * num_cams);
+	const int num_iters    = data.paramInt("num_iters",     1);
+	const int has_mb       = data.paramInt("has_mb",        1);
+	const int sel_sensors  = data.paramInt("sel_sensors",   1);
+	const float soft_margin = (float) data.param("soft_margin", 12.0);
+	const float hard_margin = (float) data.param("hard_margin",  8.0);
+	const float scale0 = (float) data.param("scale0", 1.0);
+	const float scale1 = (float) data.param("scale1", 0.0);
+	const float scale2 = (float) data.param("scale2", 0.0);
+	const int chunk = 4 * DTT_SIZE * DTT_SIZE;                        // 256 floats per TD tile
+	const size_t slice = (size_t) tilesy * tilesx * num_colors * chunk; // one cam's CLT
+	const int num_sel = __builtin_popcount(sel_sensors);
+	printf("cams=%d colors=%d frame %dx%d tiles %dx%d kernels %dx%d task_size=%d iters=%d mb=%d "
+			"sel=0x%x soft/hard=%g/%g scales=(%g,%g,%g)\n",
+			num_cams, num_colors, img_width, img_height, tilesx, tilesy, kernels_hor, kernels_vert,
+			task_size, num_iters, has_mb, sel_sensors, soft_margin, hard_margin, scale0, scale1, scale2);
+
+	// ---- static inputs ----
+	// __constant__ filter tables (correlate2D_inter reads LoG_corr + lpf_rb_corr;
+	// lpf_data/lpf_corr set too so the module state matches the Java run exactly)
+	checkCudaErrors(cudaMemcpyToSymbol(lpf_data,    data.hostFloat("lpf_data"),    4 * 64 * sizeof(float)));
+	checkCudaErrors(cudaMemcpyToSymbol(lpf_corr,    data.hostFloat("lpf_corr"),    64 * sizeof(float)));
+	checkCudaErrors(cudaMemcpyToSymbol(lpf_rb_corr, data.hostFloat("lpf_rb_corr"), 64 * sizeof(float)));
+	checkCudaErrors(cudaMemcpyToSymbol(LoG_corr,    data.hostFloat("log_corr"),    64 * sizeof(float)));
+
+	// geometry: raw float images of the Java-side structs (toFloatArray layouts);
+	// buffer sizes mirror tp_proc_setup (gc/cv 512 floats - generous fixed slots)
+	GpuBuf<float> gpu_gc (512);
+	GpuBuf<float> gpu_cv (512);
+	GpuBuf<float> gpu_rbr(RBYRDIST_LEN);
+	GpuBuf<float> gpu_rot(5 * NUM_CAMS * 3 * 3);
+	checkCudaErrors(cudaMemset(gpu_gc.dev(), 0, 512 * sizeof(float)));
+	checkCudaErrors(cudaMemset(gpu_cv.dev(), 0, 512 * sizeof(float)));
+	checkCudaErrors(cudaMemcpy(gpu_gc.dev(), data.hostFloat("geometry_correction"),
+			data.elements("geometry_correction") * sizeof(float), cudaMemcpyHostToDevice));
+	checkCudaErrors(cudaMemcpy(gpu_cv.dev(), data.hostFloat("correction_vector"),
+			data.elements("correction_vector") * sizeof(float), cudaMemcpyHostToDevice));
+
+	// per-camera images (pitched) / kernels / offsets + device pointer arrays
+	const float * h_images  = data.hostFloat("images");
+	const float * h_kernels = data.hostFloat("kernels");
+	const float * h_offsets = data.hostFloat("kernel_offsets");
+	const int kern_elems = (int) (data.elements("kernels") / num_cams);
+	const int offs_elems = (int) (data.elements("kernel_offsets") / num_cams);
+	std::vector<GpuPitchedBuf<float>> images;
+	std::vector<GpuBuf<float>> kernels, offsets;
+	std::vector<float *> images_p, kernels_p, offsets_p, clt_p, clt_ref_p;
+	for (int ncam = 0; ncam < num_cams; ncam++) {
+		images.emplace_back(img_width, img_height);
+		images.back().upload(h_images + (size_t) ncam * img_width * img_height);
+		kernels.emplace_back(kern_elems);
+		kernels.back().upload(h_kernels + (size_t) ncam * kern_elems);
+		offsets.emplace_back(offs_elems);
+		offsets.back().upload(h_offsets + (size_t) ncam * offs_elems);
+		images_p.push_back (images.back().dev());
+		kernels_p.push_back(kernels.back().dev());
+		offsets_p.push_back(offsets.back().dev());
+	}
+	// CLT buffers CONTIGUOUS across cams (the consolidation kernels index one
+	// [num_cams][slice] block - production DtoD-gathers, here direct); the
+	// per-cam pointer arrays for convert/erase/correlate point into the block
+	GpuBuf<float> clt_all    ((size_t) num_cams * slice);
+	GpuBuf<float> clt_ref_all((size_t) num_cams * slice);
+	for (int ncam = 0; ncam < num_cams; ncam++) {
+		clt_p.push_back    (clt_all.dev()     + (size_t) ncam * slice);
+		clt_ref_p.push_back(clt_ref_all.dev() + (size_t) ncam * slice);
+	}
+	GpuBuf<float *> gpu_images (num_cams); gpu_images.upload (images_p);
+	GpuBuf<float *> gpu_kernels(num_cams); gpu_kernels.upload(kernels_p);
+	GpuBuf<float *> gpu_offsets(num_cams); gpu_offsets.upload(offsets_p);
+	GpuBuf<float *> gpu_clt    (num_cams); gpu_clt.upload    (clt_p);
+	GpuBuf<float *> gpu_clt_ref(num_cams); gpu_clt_ref.upload(clt_ref_p);
+	const size_t dstride_floats = images[0].pitchFloats();
+
+	// resident virtual-center TD: cam-0 slice of gpu_clt_ref (the only slice
+	// sensor mask 1 reads); the other ref cams stay zero
+	checkCudaErrors(cudaMemset(clt_ref_all.dev(), 0, (size_t) num_cams * slice * sizeof(float)));
+	checkCudaErrors(cudaMemcpy(clt_ref_all.dev(), data.hostFloat("clt_ref_cam0"),
+			slice * sizeof(float), cudaMemcpyHostToDevice));
+
+	// sizing scratch by the largest iteration
+	size_t max_entries = 0;
+	for (int k = 0; k < num_iters; k++) {
+		char name[64];
+		snprintf(name, sizeof(name), "ftasks0_it%d", k);
+		size_t n = data.elements(name) / task_size;
+		if (n > max_entries) max_entries = n;
+	}
+	GpuBuf<float> gpu_ftasks   (max_entries * task_size);   // the live task stream (offsets/convert/correlate)
+	GpuBuf<float> gpu_cons0    (max_entries * task_size);   // consolidation inputs (oracle post-offsets streams)
+	GpuBuf<float> gpu_cons1    (max_entries * task_size);
+	GpuBuf<float> gpu_pairs    (max_entries * 2 * task_size);
+	GpuBuf<float> gpu_tasks_out(max_entries * task_size);
+	GpuBuf<int>   gpu_counters (3);                          // [0] pairs, [1] surviving, [2] misaligned
+	GpuBuf<float> gpu_td_avg   (slice);
+	GpuBuf<int>   gpu_active_tiles(tilesx * tilesy);
+	GpuBuf<int>   gpu_num_active(1);
+	const size_t max_corr = max_entries * (num_sel + 1);
+	GpuPitchedBuf<float> gpu_corrs_td(chunk, (int) max_corr);
+	GpuBuf<int> gpu_corr_indices((size_t) max_corr);
+	GpuBuf<int> gpu_num_corr(1);
+	GpuBuf<float> gpu_result_scratch(max_corr * chunk);      // reordered result upload for saveResult
+	const size_t corr_stride = gpu_corrs_td.pitchFloats();
+
+	int fail = 0;
+	const int IDX_THREADS = 256;
+	for (int k = 0; k < num_iters; k++) {
+		char nm[64];
+		// ---- the two motion-blur task sets: offsets (NEW) + convert ----
+		int n_set[2] = {0, 0};
+		const int num_sets = has_mb ? 2 : 1;
+		for (int nset = 0; nset < num_sets; nset++) {
+			snprintf(nm, sizeof(nm), "ftasks%d_it%d", nset, k);
+			const int n = (int) (data.elements(nm) / task_size);
+			n_set[nset] = n;
+			checkCudaErrors(cudaMemcpy(gpu_ftasks.dev(), data.hostFloat(nm),
+					(size_t) n * task_size * sizeof(float), cudaMemcpyHostToDevice));
+			// geometry chain, exactly as tp_proc_exec_geometry (launch dims included):
+			// reverse distortion table + rotation matrices/derivatives + per-sensor XY
+			calcReverseDistortionTable<<<16, dim3(3, 3, 3)>>>(
+					(struct gc *) gpu_gc.dev(), gpu_rbr.dev());
+			checkCudaErrors(cudaGetLastError());
+			checkCudaErrors(cudaDeviceSynchronize());
+			calc_rot_deriv<<<num_cams, dim3(3, 3, 3)>>>(
+					num_cams, (struct corr_vector *) gpu_cv.dev(), (trot_deriv *) gpu_rot.dev());
+			checkCudaErrors(cudaGetLastError());
+			checkCudaErrors(cudaDeviceSynchronize());
+			calculate_tiles_offsets<<<1, 1>>>(
+					0,                                    // uniform_grid: 0 = use provided centers
+					num_cams,
+					gpu_ftasks.dev(),
+					n,
+					(struct gc *) gpu_gc.dev(),
+					(struct corr_vector *) gpu_cv.dev(),
+					gpu_rbr.dev(),
+					(trot_deriv *) gpu_rot.dev());
+			checkCudaErrors(cudaGetLastError());
+			checkCudaErrors(cudaDeviceSynchronize());
+			// checkpoint: GPU-filled tasks vs the oracle post-offsets stream
+			{
+				std::vector<float> got((size_t) n * task_size);
+				checkCudaErrors(cudaMemcpy(got.data(), gpu_ftasks.dev(),
+						got.size() * sizeof(float), cudaMemcpyDeviceToHost));
+				char lbl[64];
+				snprintf(lbl, sizeof(lbl), "tasks%d_it%d", nset, k);
+				snprintf(nm, sizeof(nm), "expected_tasks%d_it%d", nset, k);
+				if (compareTasks(got.data(), data.hostFloat(nm), n, task_size, num_cams, lbl)) fail = 1;
+			}
+			if (nset == 0) { // erase before the FIRST (SET) convert only, as interCorrTDMotionBlur
+				erase_clt_tiles<<<1, 1>>>(num_cams, num_colors, tilesx, tilesy, gpu_clt.dev(), NAN);
+				checkCudaErrors(cudaGetLastError());
+				checkCudaErrors(cudaDeviceSynchronize());
+			}
+			convert_direct<<<1, 1>>>(
+					num_cams, num_colors,
+					gpu_offsets.dev(), gpu_kernels.dev(), gpu_images.dev(),
+					gpu_ftasks.dev(), gpu_clt.dev(),
+					dstride_floats, n,
+					0,                                    // lpf_mask - always 0
+					img_width, img_height,
+					kernels_hor, kernels_vert,            // deconvolution applied (production)
+					gpu_active_tiles.dev(), gpu_num_active.dev(), tilesx);
+			checkCudaErrors(cudaGetLastError());
+			checkCudaErrors(cudaDeviceSynchronize());
+		}
+		// ---- consolidation (the proven v1 chain) into the cam-0 slot ----
+		snprintf(nm, sizeof(nm), "expected_tasks0_it%d", k);
+		checkCudaErrors(cudaMemcpy(gpu_cons0.dev(), data.hostFloat(nm),
+				(size_t) n_set[0] * task_size * sizeof(float), cudaMemcpyHostToDevice));
+		float * d_cons1 = nullptr;
+		if (has_mb) {
+			snprintf(nm, sizeof(nm), "expected_tasks1_it%d", k);
+			checkCudaErrors(cudaMemcpy(gpu_cons1.dev(), data.hostFloat(nm),
+					(size_t) n_set[1] * task_size * sizeof(float), cudaMemcpyHostToDevice));
+			d_cons1 = gpu_cons1.dev();
+		}
+		checkCudaErrors(cudaMemset(gpu_counters.dev(), 0, 3 * sizeof(int)));
+		checkCudaErrors(cudaMemset(gpu_td_avg.dev(), 0xFF, slice * sizeof(float))); // NaN poison prefill
+		index_consolidate<<<(n_set[0] + IDX_THREADS - 1) / IDX_THREADS, IDX_THREADS>>>(
+				task_size, gpu_cons0.dev(), d_cons1, n_set[0],
+				gpu_pairs.dev(), gpu_counters.dev() + 0, gpu_counters.dev() + 2);
+		checkCudaErrors(cudaGetLastError());
+		int num_pairs = 0, num_out = 0, num_mis = 0;
+		checkCudaErrors(cudaMemcpy(&num_pairs, gpu_counters.dev(), sizeof(int), cudaMemcpyDeviceToHost));
+		if (num_pairs > 0) {
+			consolidate_oob<<<(num_pairs + IDX_THREADS - 1) / IDX_THREADS, IDX_THREADS>>>(
+					task_size, num_cams, gpu_pairs.dev(), num_pairs,
+					(float) img_width, (float) img_height, soft_margin, hard_margin,
+					gpu_tasks_out.dev(), gpu_counters.dev() + 1);
+			checkCudaErrors(cudaGetLastError());
+		}
+		checkCudaErrors(cudaMemcpy(&num_out, gpu_counters.dev() + 1, sizeof(int), cudaMemcpyDeviceToHost));
+		if (num_out > 0) {
+			float * counts = nullptr; // production bridge passes no counts plane
+			clt_average_sensors_list<<<dim3(num_out, num_colors), TD_TILE_FLOATS>>>(
+					num_cams, clt_all.dev(), tilesx, tilesy, num_colors,
+					gpu_tasks_out.dev(), task_size, num_out,
+					gpu_td_avg.dev(), counts);
+			checkCudaErrors(cudaGetLastError());
+		}
+		checkCudaErrors(cudaDeviceSynchronize());
+		checkCudaErrors(cudaMemcpy(&num_mis, gpu_counters.dev() + 2, sizeof(int), cudaMemcpyDeviceToHost));
+		// the average -> the cam-0 slot (what the single conj-multiply correlates)
+		checkCudaErrors(cudaMemcpy(clt_all.dev(), gpu_td_avg.dev(), slice * sizeof(float),
+				cudaMemcpyDeviceToDevice));
+		printf("it%d: sets %d+%d tasks -> %d pairs -> %d consolidated tiles (%d misaligned)\n",
+				k, n_set[0], n_set[1], num_pairs, num_out, num_mis);
+		if (num_mis) fail = 1;
+		// ---- the single conj-multiply vs the resident center TD (NEW) ----
+		// device ftasks still holds the LAST uploaded set (set 1 in MB mode) -
+		// exactly the state execCorr2D_inter_TD sees in production
+		const int n_corr_tasks = has_mb ? n_set[1] : n_set[0];
+		correlate2D_inter<<<1, 1>>>(
+				num_cams, sel_sensors,
+				gpu_clt.dev(), gpu_clt_ref.dev(),
+				num_colors, scale0, scale1, scale2,
+				gpu_ftasks.dev(), n_corr_tasks, tilesx,
+				gpu_corr_indices.dev(), gpu_num_corr.dev(),
+				corr_stride, gpu_corrs_td.dev());
+		checkCudaErrors(cudaGetLastError());
+		checkCudaErrors(cudaDeviceSynchronize());
+		int num_corr = 0;
+		checkCudaErrors(cudaMemcpy(&num_corr, gpu_num_corr.dev(), sizeof(int), cudaMemcpyDeviceToHost));
+		// ---- checkpoint 1 compare: raw TD correlation, order-independent ----
+		snprintf(nm, sizeof(nm), "expected_corr_indices_it%d", k);
+		const int * exp_idx = data.hostInt(nm);
+		const int n_exp = (int) data.elements(nm);
+		snprintf(nm, sizeof(nm), "expected_corr_td_it%d", k);
+		const float * exp_td = data.hostFloat(nm);
+		std::vector<int> got_idx(num_corr);
+		checkCudaErrors(cudaMemcpy(got_idx.data(), gpu_corr_indices.dev(),
+				(size_t) num_corr * sizeof(int), cudaMemcpyDeviceToHost));
+		std::vector<float> got_td((size_t) num_corr * chunk);
+		checkCudaErrors(cudaMemcpy2D(got_td.data(), (size_t) chunk * sizeof(float),
+				gpu_corrs_td.dev(), gpu_corrs_td.pitchBytes(),
+				(size_t) chunk * sizeof(float), num_corr, cudaMemcpyDeviceToHost));
+		if (num_corr != n_exp) {
+			printf("it%d: num_corr_tiles %d != expected %d -> FAIL\n", k, num_corr, n_exp);
+			fail = 1;
+		}
+		std::map<int, int> got_by_index; // packed (tile<<CORR_NTILE_SHIFT)|sensor(0xff=sum) -> row
+		for (int i = 0; i < num_corr; i++) got_by_index[got_idx[i]] = i;
+		std::vector<float> reordered((size_t) n_exp * chunk, NAN); // oracle order, for results/ + Java round trip
+		double maxd = 0;
+		long nan_mm = 0, missing = 0;
+		long first_bad = -1;
+		for (int i = 0; i < n_exp; i++) {
+			const std::map<int, int>::const_iterator it = got_by_index.find(exp_idx[i]);
+			if (it == got_by_index.end()) { missing++; continue; }
+			const float * g = got_td.data() + (size_t) it->second * chunk;
+			const float * e = exp_td + (size_t) i * chunk;
+			memcpy(reordered.data() + (size_t) i * chunk, g, chunk * sizeof(float));
+			for (int j = 0; j < chunk; j++) {
+				const bool gn = std::isnan(g[j]), en = std::isnan(e[j]);
+				if (gn || en) {
+					if (gn != en) { nan_mm++; if (first_bad < 0) first_bad = (size_t) i * chunk + j; }
+					continue;
+				}
+				const double d = std::fabs((double) g[j] - (double) e[j]);
+				if (d > maxd) { maxd = d; if ((d > tol) && (first_bad < 0)) first_bad = (size_t) i * chunk + j; }
+			}
+		}
+		const int it_fail = (maxd > tol) || (nan_mm > 0) || (missing > 0) || (num_corr != n_exp);
+		printf("it%d: compare 'corr_td' (%d tiles, by packed index): max|diff|=%g NaN-mismatch=%ld "
+				"missing=%ld first@%ld -> %s\n",
+				k, n_exp, maxd, nan_mm, missing, first_bad, it_fail ? "FAIL" : "PASS");
+		if (it_fail) fail = 1;
+		// save the reordered result for the Java round trip
+		checkCudaErrors(cudaMemcpy(gpu_result_scratch.dev(), reordered.data(),
+				reordered.size() * sizeof(float), cudaMemcpyHostToDevice));
+		snprintf(nm, sizeof(nm), "corr_td_it%d", k);
+		data.saveResult(nm, gpu_result_scratch.dev(), {n_exp, chunk});
+	}
+	printf("%s: %s\n", argv[0], fail ? "FAIL" : "PASS");
+	return fail ? EXIT_FAILURE : EXIT_SUCCESS;
+}