CLAUDE: peak kernel precision two-step - templated body, float instantiation...

CLAUDE: peak kernel precision two-step - templated body, float instantiation corr2D_peak_eig_f (Andrey's methodology)

Step 1 (double vs double Java oracle) was the validated state; the body is
now a template on the accumulator type with two extern-C instantiations:
corr2D_peak_eig (double, the oracle-matched reference) and corr2D_peak_eig_f
(float, ~64x cheaper FLOPs + half the shared footprint - kept adequate for
full-frame, Andrey's ruling). tp_proc_exec_corr2d_peak grows a use_float arg;
NVRTC module = 25 kernels. All three tests run BOTH and report float-vs-double
grouped stats (dx/dy px, strength rel, eigvec abs, lambda rel).

Step 2 measured:
- synthetic: dx/dy <= 4.8e-7 px, lambda <= 9.4e-7 rel; eigvec 0.031 ONLY on
  the (near-)isotropic tiles where the direction is ill-defined (float can't
  fire the oracle's e=1e-12 branch) AND unused (setEigenTransform degenerates
  to isotropic weighting when lambda0~lambda1) - documented in tp_peak.cu.
- REAL case (4 iterations x 150 tiles): dx/dy <= 7.2e-7 px, strength
  IDENTICAL, eigvec <= 8.9e-5, lambda <= 1.2e-6 rel - 4+ orders below the
  ~0.01 px measurement noise. Float is adequate for production; double stays
  as the validation/oracle tier. Double regression unchanged (synthetic 0 /
  9.6e-17; real case ALL PASS).
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

jna/test_pose_corr_jna.cu

@@ -74,7 +74,7 @@ extern "C" {
   int   tp_proc_exec_corr2d_normalize(void*, int, double, int);
   int   tp_proc_get_corr2d(void*, float*, int);
   int   tp_proc_set_peak_debias(void*, const float*, int);
-  int   tp_proc_exec_corr2d_peak(void*, int, float,float,float,float, int, float,float,float);
+  int   tp_proc_exec_corr2d_peak(void*, int, float,float,float,float, int, float,float,float, int);
   int   tp_proc_get_peaks(void*, float*);
   void  tp_proc_destroy(void*);
   void  tp_destroy_module(void*);
@@ -282,14 +282,22 @@ int main(int argc, char** argv){
       tp_proc_get_corr2d(p, got_pd.data(), corr_rad);
       const int n_peaks = tp_proc_exec_corr2d_peak(p, 0xff,
           p_abs_min, p_rel_min, p_min_peak, p_sub_frac,
-          p_recenter, p_sub_frac1, p_scale_ax, p_inc_ax);
+          p_recenter, p_sub_frac1, p_scale_ax, p_inc_ax, 0); // double = the PASS-gated step-1 kernel
       if (n_peaks < 0){ printf("FAIL exec_corr2d_peak: %s\n", tp_last_error()); return EXIT_FAILURE; }
       std::vector<float> got_pk((size_t)n_peaks*PEAK_ROW_FLOATS);
       tp_proc_get_peaks(p, got_pk.data());
+      // precision step 2 (informational): the FLOAT instantiation on the same data
+      std::vector<float> got_pkf((size_t)n_peaks*PEAK_ROW_FLOATS);
+      if (tp_proc_exec_corr2d_peak(p, 0xff, p_abs_min, p_rel_min, p_min_peak, p_sub_frac,
+          p_recenter, p_sub_frac1, p_scale_ax, p_inc_ax, 1) == n_peaks){
+        tp_proc_get_peaks(p, got_pkf.data());
+      } else { printf("WARN: float peak kernel failed: %s\n", tp_last_error()); got_pkf = got_pk; }
       // match by packed index (dense tile -> (tile<<8)|0xff sum slot), oracle row order out
       std::vector<float> res_pd((size_t)n_pd*pd_size, NAN);
       std::vector<float> res_pk((size_t)n_pd*7, NAN);
       double pd_maxd = 0, pk_maxd = 0, pk_maxulp = 0; long pd_nan = 0, pk_nan = 0, pd_missing = 0, pk_viol = 0;
+      // float-vs-double per-component-group maxima (precision step 2, informational)
+      double f_dxy = 0, f_str = 0, f_eigv = 0, f_lam = 0; long f_nan = 0;
       for (int i = 0; i < n_pd; i++){
         const int packed = (pd_idx[i] << 8) | 0xff; // CORR_NTILE_SHIFT
         std::map<int,int>::const_iterator it = by_index.find(packed);
@@ -322,6 +330,18 @@ int main(int argc, char** argv){
           if (du > pk_maxulp) pk_maxulp = du;
           if ((d > peak_tol) && (du > PEAK_TOL_ULP)) pk_viol++;
         }
+        // float-vs-double, grouped: dx/dy absolute (px), strength relative,
+        // eigenvector components absolute (unit vector), lambdas relative
+        const float* gpf = got_pkf.data() + (size_t)it->second*PEAK_ROW_FLOATS;
+        for (int j = 0; j < 7; j++){
+          const bool fn = std::isnan(gpf[j]), dn = std::isnan(gpk[j]);
+          if (fn || dn){ if(fn!=dn) f_nan++; continue; }
+          const double d = std::fabs((double)gpf[j]-(double)gpk[j]);
+          if      (j < 2)  { if (d > f_dxy) f_dxy = d; }
+          else if (j == 2) { const double r = d/std::max(1e-30,std::fabs((double)gpk[j])); if (r > f_str) f_str = r; }
+          else if (j < 5)  { if (d > f_eigv) f_eigv = d; }
+          else             { const double r = d/std::max(1e-30,std::fabs((double)gpk[j])); if (r > f_lam) f_lam = r; }
+        }
       }
       const int pd_fail = (pd_maxd>tol) || (pd_nan>0) || (pd_missing>0);
       const int pk_fail = (pk_viol>0)   || (pk_nan>0) || (pd_missing>0);
@@ -329,6 +349,8 @@ int main(int argc, char** argv){
              k, n_pd, pd_maxd, pd_nan, pd_missing, pd_fail?"FAIL":"PASS");
       printf("it%d: compare 'peaks' (%d tiles vs the double Java oracle): max|diff|=%g max=%.2f ULP viol=%ld NaN-mismatch=%ld -> %s (abs %g | %d ULP)\n",
              k, n_pd, pk_maxd, pk_maxulp, pk_viol, pk_nan, pk_fail?"FAIL":"PASS", peak_tol, PEAK_TOL_ULP);
+      printf("it%d: float-vs-double peaks (informational, step 2): max |d dx,dy|=%.3g px, str rel=%.3g, |d eigvec|=%.3g, lambda rel=%.3g, NaN-mismatch=%ld\n",
+             k, f_dxy, f_str, f_eigv, f_lam, f_nan);
       if (pd_fail || pk_fail) fail = 1;
       cudaMalloc((void**)&scratch, res_pd.size()*sizeof(float));
       cudaMemcpy(scratch, res_pd.data(), res_pd.size()*sizeof(float), cudaMemcpyHostToDevice);

jna/tp_jna.cpp

 // Stage-0b: the probe refactored into a JNA-callable shared lib (libtileproc.so).
-// extern "C": tp_create_module (NVRTC-compile kernels + cuLink(libcudadevrt) + load 24 funcs), // By Claude on 07/13/2026: +corr2D_peak_eig (rung 2)
+// extern "C": tp_create_module (NVRTC-compile kernels + cuLink(libcudadevrt) + load 25 funcs), // By Claude on 07/14/2026: +corr2D_peak_eig{,_f} (rung 2)
 // tp_module_num_functions, tp_last_error, tp_destroy_module. By Claude on 2026-06-25.
 #include <cuda.h>
 #include <nvrtc.h>
@@ -47,7 +47,7 @@ static const char* KERNELS[] = {
   "gen_texture_list","clear_texture_rbga","textures_accumulate","create_nonoverlap_list","erase_clt_tiles",
   "calculate_tiles_offsets","calc_rot_deriv","calcReverseDistortionTable",
   "clt_average_sensors","index_consolidate","consolidate_oob","clt_average_sensors_list", // By Claude on 07/12/2026
-  "corr2D_peak_eig"}; // By Claude on 07/13/2026
+  "corr2D_peak_eig","corr2D_peak_eig_f"}; // double + float instantiations (precision two-step) // By Claude on 07/13-14/2026
 static const int N_KERNELS = sizeof(KERNELS)/sizeof(KERNELS[0]);
 
 struct TpModule { CUcontext ctx; CUmodule mod; int nfun; };
@@ -985,12 +985,15 @@ int tp_proc_set_peak_debias(TpProc* p, const float* data, int len){
 // (gpu_corrs, i.e. AFTER tp_proc_exec_corr2d_normalize combo=0) - one row of
 // PEAK_ROW_FLOATS {dx,dy,str,e0x,e0y,l0,l1,valid} per resident corr tile.
 // pair_select: process only slots with this pair code (lean sum = 0xff; -1 = all).
-// Returns the number of corr tiles (rows) or <0 on error.
+// use_float: 0 = double accumulators (the oracle-matched step-1 kernel),
+// !=0 = float instantiation (precision step 2). Returns the row count or <0.
 int tp_proc_exec_corr2d_peak(TpProc* p, int pair_select,
                              float abs_min, float rel_min, float min_peak, float eig_sub_frac,
-                             int refine, float eig_sub_frac1, float scale_axes, float inc_axes){
+                             int refine, float eig_sub_frac1, float scale_axes, float inc_axes,
+                             int use_float){ // By Claude on 07/14/2026 (precision two-step)
   if(!p||!p->have_corr){ seterr("exec_corr2d_peak: no corr buffers"); return -1; } cuCtxSetCurrent(p->mod->ctx);
-  CUfunction f=getfun(p->mod,"corr2D_peak_eig"); if(!f){ seterr("corr2D_peak_eig missing"); return -2; }
+  CUfunction f=getfun(p->mod, use_float ? "corr2D_peak_eig_f" : "corr2D_peak_eig");
+  if(!f){ seterr("corr2D_peak_eig%s missing", use_float?"_f":""); return -2; }
   int nct=p->last_num_corr_tiles;
   if(nct<=0){ seterr("exec_corr2d_peak: no corr tiles (%d)",nct); return -3; }
   if(nct>p->peak_capacity){

src/tests/test_peak.cu

@@ -90,6 +90,16 @@ int main(int argc, char **argv)
 	checkCudaErrors(cudaGetLastError());
 	checkCudaErrors(cudaDeviceSynchronize());
 	std::vector<float> got = gpu_peaks.download();
+	// precision step 2 (informational): the FLOAT instantiation on the same tiles
+	corr2D_peak_eig_f<<<num_corr, 256>>>(
+			num_corr, gpu_idx.dev(), 0xff,
+			gpu_pd.pitchFloats(), gpu_pd.dev(), corr_rad,
+			p_abs_min, p_rel_min, p_min_peak, p_sub_frac,
+			p_recenter, p_sub_frac1, p_scale_ax, p_inc_ax,
+			d_debias, gpu_peaks.dev());
+	checkCudaErrors(cudaGetLastError());
+	checkCudaErrors(cudaDeviceSynchronize());
+	std::vector<float> gotf = gpu_peaks.download();
 
 	const int * pd_idx = data.hostInt("expected_pd_indices");
 	const int   n_exp  = (int) data.elements("expected_pd_indices");
@@ -131,6 +141,26 @@ int main(int argc, char **argv)
 	if ((maxd > tol) || (nan_mm > 0)) fail = 1;
 	printf("compare 'peaks' (%d tiles): max|diff|=%g NaN-mismatch=%ld tol=%g -> %s\n",
 			n_exp, maxd, nan_mm, tol, ((maxd > tol) || nan_mm) ? "FAIL" : "PASS");
+	// float-vs-double summary (informational, precision step 2)
+	{
+		double f_dxy = 0, f_str = 0, f_eigv = 0, f_lam = 0;
+		long f_nan = 0;
+		for (int i = 0; i < num_corr; i++) {
+			const float * gd = got.data()  + (size_t) i * PEAK_ROW_FLOATS;
+			const float * gf = gotf.data() + (size_t) i * PEAK_ROW_FLOATS;
+			for (int j = 0; j < 7; j++) {
+				const bool fn = std::isnan(gf[j]), dn = std::isnan(gd[j]);
+				if (fn || dn) { if (fn != dn) f_nan++; continue; }
+				const double d = std::fabs((double) gf[j] - (double) gd[j]);
+				if      (j < 2)  { if (d > f_dxy) f_dxy = d; }
+				else if (j == 2) { const double r = d / std::max(1e-30, std::fabs((double) gd[j])); if (r > f_str) f_str = r; }
+				else if (j < 5)  { if (d > f_eigv) f_eigv = d; }
+				else             { const double r = d / std::max(1e-30, std::fabs((double) gd[j])); if (r > f_lam) f_lam = r; }
+			}
+		}
+		printf("float-vs-double peaks (informational): max |d dx,dy|=%.3g px, str rel=%.3g, |d eigvec|=%.3g, lambda rel=%.3g, NaN-mismatch=%ld\n",
+				f_dxy, f_str, f_eigv, f_lam, f_nan);
+	}
 	printf("%s: %s\n", argv[0], fail ? "FAIL" : "PASS");
 	return fail ? EXIT_FAILURE : EXIT_SUCCESS;
 }

src/tests/test_pose_corr.cu

@@ -503,12 +503,27 @@ int main(int argc, char **argv)
 			std::vector<float> got_pk((size_t) num_corr * PEAK_ROW_FLOATS);
 			checkCudaErrors(cudaMemcpy(got_pk.data(), gpu_peaks.dev(),
 					got_pk.size() * sizeof(float), cudaMemcpyDeviceToHost));
+			// precision step 2 (informational): the FLOAT instantiation on the same PD
+			corr2D_peak_eig_f<<<num_corr, 256>>>(
+					num_corr, gpu_corr_indices.dev(), 0xff,
+					gpu_corrs_pd.pitchFloats(), gpu_corrs_pd.dev(), corr_rad,
+					p_abs_min, p_rel_min, p_min_peak, p_sub_frac,
+					p_recenter, p_sub_frac1, p_scale_ax, p_inc_ax,
+					d_debias, gpu_peaks.dev());
+			checkCudaErrors(cudaGetLastError());
+			checkCudaErrors(cudaDeviceSynchronize());
+			std::vector<float> got_pkf((size_t) num_corr * PEAK_ROW_FLOATS);
+			checkCudaErrors(cudaMemcpy(got_pkf.data(), gpu_peaks.dev(),
+					got_pkf.size() * sizeof(float), cudaMemcpyDeviceToHost));
 			// match by packed index (dense tile -> (tile<<shift)|0xff sum slot), reusing
 			// the checkpoint-1 map; save results in ORACLE row order for the Java round trip
 			std::vector<float> res_pd((size_t) n_pd * pd_size, NAN);
 			std::vector<float> res_pk((size_t) n_pd * 7, NAN);
 			double pd_maxd = 0, pk_maxd = 0, pk_maxulp = 0;
 			long pd_nan = 0, pk_nan = 0, pd_missing = 0, pk_viol = 0;
+			// float-vs-double per-component-group maxima (precision step 2, informational)
+			double f_dxy = 0, f_str = 0, f_eigv = 0, f_lam = 0;
+			long f_nan = 0;
 			for (int i = 0; i < n_pd; i++) {
 				const int packed = (pd_idx[i] << CORR_NTILE_SHIFT) | 0xff;
 				const std::map<int, int>::const_iterator it = got_by_index.find(packed);
@@ -538,6 +553,18 @@ int main(int argc, char **argv)
 					if (du > pk_maxulp) pk_maxulp = du;
 					if ((d > peak_tol) && (du > PEAK_TOL_ULP)) pk_viol++;
 				}
+				// float-vs-double, grouped: dx/dy absolute (px), strength relative,
+				// eigenvector components absolute (unit vector), lambdas relative
+				const float * gpf = got_pkf.data() + (size_t) it->second * PEAK_ROW_FLOATS;
+				for (int j = 0; j < 7; j++) {
+					const bool fn = std::isnan(gpf[j]), dn = std::isnan(gpk[j]);
+					if (fn || dn) { if (fn != dn) f_nan++; continue; }
+					const double d = std::fabs((double) gpf[j] - (double) gpk[j]);
+					if      (j < 2)  { if (d > f_dxy) f_dxy = d; }
+					else if (j == 2) { const double r = d / std::max(1e-30, std::fabs((double) gpk[j])); if (r > f_str) f_str = r; }
+					else if (j < 5)  { if (d > f_eigv) f_eigv = d; }
+					else             { const double r = d / std::max(1e-30, std::fabs((double) gpk[j])); if (r > f_lam) f_lam = r; }
+				}
 			}
 			const int pd_fail = (pd_maxd > tol) || (pd_nan > 0) || (pd_missing > 0);
 			const int pk_fail = (pk_viol > 0)   || (pk_nan > 0) || (pd_missing > 0);
@@ -546,6 +573,9 @@ int main(int argc, char **argv)
 			printf("it%d: compare 'peaks' (%d tiles vs the double Java oracle): max|diff|=%g max=%.2f ULP "
 					"viol=%ld NaN-mismatch=%ld -> %s (abs %g | %d ULP)\n",
 					k, n_pd, pk_maxd, pk_maxulp, pk_viol, pk_nan, pk_fail ? "FAIL" : "PASS", peak_tol, PEAK_TOL_ULP);
+			printf("it%d: float-vs-double peaks (informational, step 2): max |d dx,dy|=%.3g px, str rel=%.3g, "
+					"|d eigvec|=%.3g, lambda rel=%.3g, NaN-mismatch=%ld\n",
+					k, f_dxy, f_str, f_eigv, f_lam, f_nan);
 			if (pd_fail || pk_fail) fail = 1;
 			checkCudaErrors(cudaMemcpy(gpu_result_scratch.dev(), res_pd.data(),
 					res_pd.size() * sizeof(float), cudaMemcpyHostToDevice));

src/tp_peak.h

@@ -20,11 +20,15 @@
  * NOT in the contract: the oracle's fpn_mask/ignore_border branch — the lean
  * chain always passes null (input is FPN-subtracted); host-side policy.
  *
- * All accumulation and the eigen/refine math run in DOUBLE with the SAME
- * summation order as the Java oracle (serial linescan), so the float-input /
- * double-math results track the double oracle to rounding (offline-nvcc vs
- * NVRTC FMA/codegen caveats per the rung-1 lesson: the tol-0 verdict is
- * NVRTC-replay-vs-production-capture, Java-oracle compare is an epsilon tier).
+ * PRECISION — TWO VARIANTS, same contract (Andrey's two-step methodology
+ * 07/14/2026): corr2D_peak_eig accumulates in DOUBLE with the SAME summation
+ * order as the Java oracle (serial linescan) — step 1, matched bit-tight
+ * against the oracle (1-2 float-ULP cast flips). corr2D_peak_eig_f is the
+ * FLOAT instantiation of the same templated body — step 2, measured against
+ * both the oracle and the double variant on real data to pick the production
+ * precision (as convert_direct/correlations: float GPU vs double Java).
+ * Offline-nvcc vs NVRTC FMA/codegen caveats per the rung-1 lesson: the tol-0
+ * verdict is NVRTC-replay-vs-production-capture.
  *
  * Java oracle: Correlation2d.getMaxXYCmEig (15-arg, eig_fast2x2=true);
  * case export: PoseCorrExport (curt.kernel_test = pose_corr, checkpoint 2).
@@ -65,4 +69,23 @@ extern "C" __global__ void corr2D_peak_eig(
 		float       * gpu_peaks);       // [num_corr_tiles][PEAK_ROW_FLOATS]; rejected/skipped ->
 		                                // 7x NaN + valid 0.0f; measured -> values + valid 1.0f
 
+// FLOAT-accumulator instantiation of the same body/contract (precision step 2)
+extern "C" __global__ void corr2D_peak_eig_f(
+		int           num_corr_tiles,
+		const int   * gpu_corr_indices,
+		int           pair_select,
+		const size_t  corr_stride,
+		const float * gpu_corrs,
+		int           corr_radius,
+		float         abs_min,
+		float         rel_min,
+		float         min_peak,
+		float         eig_sub_frac,
+		int           refine,
+		float         eig_sub_frac1,
+		float         scale_axes,
+		float         inc_axes,
+		const float * gpu_debias,
+		float       * gpu_peaks);
+
 #endif /* SRC_TP_PEAK_H_ */