CLAUDE: 3-C rung 2 - corr2D_peak_eig (frozen 15-arg getMaxXYCmEig contract) +...

CLAUDE: 3-C rung 2 - corr2D_peak_eig (frozen 15-arg getMaxXYCmEig contract) + checkpoint-2 replay + synthetic numpy gate

Andrey's ruling 07/13/2026: the GPU peak contract = the 15-arg superset -
closed-form 2x2 eigen, recentering as a runtime param, envelope de-bias as
an optional array pointer (NULL = off).

- src/tp_peak.{h,cu}: corr2D_peak_eig, one block/tile - parallel load/NaN/
  argmax(first-occurrence)/8-connected grow, thread-0 serial DOUBLE centroid/
  covariance/eigen/recentering in the Java oracle's exact summation order.
- jna/tp_jna.cpp: kernel in the NVRTC module (24 kernels) +
  tp_proc_exec_corr2d_peak / tp_proc_get_peaks / tp_proc_set_peak_debias.
- test_pose_corr + test_pose_corr_jna: checkpoint-2 replay (in-place
  corr2D_normalize -> corr_pd tol tier, peaks vs the double Java oracle at
  --peak-tol) when the case carries it; pre-rung-2 cases skip with a note.
- make_peak_case.py + test_peak: synthetic gate vs an INDEPENDENT numpy
  float64 port (8 adversarial tiles + refine/debias case) - PASS at
  max|diff| 0 / 9.6e-17.
- Regressions: test_pose_corr_jna bit-exact @tol 0 on the existing v013 case,
  avg_td / avg_td_oob / convert_direct PASS.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

jna/test_pose_corr_jna.cu

@@ -8,6 +8,16 @@
 //   exec_corr2d_inter_td vs the resident center TD (set_clt ref cam 0) ->
 //   get_corr_indices/get_corr_td vs oracle, ORDER-INDEPENDENT by packed index.
 //
+// CHECKPOINT 2 (rung 2, 07/13/2026), when the case carries expected_pd_indices:
+//   exec_corr2d_normalize(combo=0) IN PLACE on the inter buffer (production
+//   path) -> get_corr2d vs expected_corr_pd AT TOL 0 (same NVRTC module ==
+//   the production capture) -> set_peak_debias (if the case has one) +
+//   exec_corr2d_peak (frozen 15-arg getMaxXYCmEig contract, manifest params)
+//   -> get_peaks vs the double-precision Java oracle at --peak-tol (default
+//   1e-6; double kernel math vs Java double - FMA/libm ULPs, not a defect
+//   tier). Once production runs the kernel and re-exports GPU-produced peaks,
+//   THAT compare becomes the tol-0 tier. Old cases skip with a note.
+//
 // Why this exists next to src/tests/test_pose_corr.cu (direct launches, -G,
 // cuda-gdb steppable): an offline nvcc build cannot be bit-exact against the
 // NVRTC-JIT production module (FMA/codegen divergence; measured: tasks match
@@ -36,6 +46,7 @@
 
 #include <cuda_runtime.h>
 #include "tests/tp_test_data.h"
+#include "tp_peak.h"        // PEAK_ROW_FLOATS (checkpoint 2) // By Claude on 07/13/2026
 
 extern "C" {
   void* tp_create_module(const char*, const char*);
@@ -59,6 +70,12 @@ extern "C" {
   int   tp_proc_get_corr_indices(void*, int*, int);
   int   tp_proc_get_corr_td(void*, float*);
   int   tp_proc_num_corr_tiles(void*);
+  // checkpoint 2 (rung 2) // By Claude on 07/13/2026
+  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_get_peaks(void*, float*);
   void  tp_proc_destroy(void*);
   void  tp_destroy_module(void*);
 }
@@ -91,14 +108,16 @@ int main(int argc, char** argv){
   std::string srcdir   = "/home/elphel/git/tile_processor_gpu/src";
   std::string devrt    = "/usr/local/cuda-12.8/lib64/libcudadevrt.a";
   float tol = 0.0f;
+  float peak_tol = 1e-6f; // Java-double-oracle tier for the peak kernel (header note) // By Claude on 07/13/2026
   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],"--peak-tol") && (i+1<argc)) peak_tol = (float)atof(argv[++i]);
     else if (!strcmp(argv[i],"--src")   && (i+1<argc)) srcdir = argv[++i];
     else if (!strcmp(argv[i],"--devrt") && (i+1<argc)) devrt = argv[++i];
     else if (!strcmp(argv[i],"--list"))                list_only = 1;
-    else { fprintf(stderr,"Usage: %s --data <dir> [--tol X] [--src <kernel src>] [--devrt <libcudadevrt.a>] [--list]\n", argv[0]); return EXIT_FAILURE; }
+    else { fprintf(stderr,"Usage: %s --data <dir> [--tol X] [--peak-tol X] [--src <kernel src>] [--devrt <libcudadevrt.a>] [--list]\n", argv[0]); return EXIT_FAILURE; }
   }
   printf("%s: data=%s tol=%g (NVRTC module from %s)\n", argv[0], data_dir.c_str(), tol, srcdir.c_str());
   TpTestData data(data_dir);
@@ -153,6 +172,23 @@ int main(int argc, char** argv){
   const int num_pairs = num_cams*(num_cams-1)/2;
   tp_proc_setup_rbg_corr(p, num_pairs, 0,0,0,0, 1.f,1.f,1.f, 7);
 
+  // checkpoint 2 (rung 2) statics: FZ-norm in place + frozen peak contract // By Claude on 07/13/2026
+  const int  corr_rad = 7, pd_size = 15*15;
+  const bool have_ckpt2 = data.has("expected_pd_indices_it0");
+  const double fz_inter = data.param("fz_inter", 10000.0);
+  const float p_abs_min   = (float)data.param("eig_min_abs",   0.0);
+  const float p_rel_min   = (float)data.param("eig_min_rel",   0.0);
+  const float p_min_peak  = (float)data.param("min_str_sum",   0.0);
+  const float p_sub_frac  = (float)data.param("eig_sub_frac",  0.0);
+  const int   p_recenter  =        data.paramInt("pose_recenter", 0);
+  const float p_sub_frac1 = (float)data.param("eig_sub_frac1", 0.0);
+  const float p_scale_ax  = (float)data.param("eig_scale_axes",1.2);
+  const float p_inc_ax    = (float)data.param("eig_inc_axes",  1.0);
+  if (have_ckpt2){
+    if (data.has("debias")) tp_proc_set_peak_debias(p, data.hostFloat("debias"), (int)data.elements("debias"));
+    printf("checkpoint 2: fz=%g recenter=%d debias=%s\n", fz_inter, p_recenter, data.has("debias")?"ON":"off");
+  } else printf("NOTE: checkpoint-2 buffers absent (pre-rung-2 case vintage) - FZ-norm/peak replay skipped\n");
+
   int fail = 0;
   char nm[64];
   for (int k = 0; k < num_iters; k++){
@@ -228,6 +264,72 @@ int main(int argc, char** argv){
     snprintf(nm,sizeof(nm),"corr_td_it%d",k);
     data.saveResult(nm, scratch, {n_exp, chunk});
     cudaFree(scratch);
+
+    // ---- CHECKPOINT 2 (rung 2): FZ-norm IN PLACE + peak kernel, production entry points ----
+    // By Claude on 07/13/2026
+    if (have_ckpt2){
+      snprintf(nm,sizeof(nm),"expected_pd_indices_it%d",k);
+      const int* pd_idx = data.hostInt(nm);
+      const int  n_pd   = (int)data.elements(nm);
+      snprintf(nm,sizeof(nm),"expected_corr_pd_it%d",k);
+      const float* exp_pd = data.hostFloat(nm);
+      snprintf(nm,sizeof(nm),"expected_peaks_it%d",k);
+      const float* exp_pk = data.hostFloat(nm);
+      if(tp_proc_exec_corr2d_normalize(p, 0, fz_inter, corr_rad)){
+        printf("FAIL exec_corr2d_normalize: %s\n", tp_last_error()); return EXIT_FAILURE; }
+      std::vector<float> got_pd((size_t)num_corr*pd_size);
+      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);
+      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());
+      // 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; long pd_nan = 0, pk_nan = 0, pd_missing = 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);
+        if (it == by_index.end()){ pd_missing++; continue; }
+        const float* gpd = got_pd.data() + (size_t)it->second*pd_size;
+        const float* epd = exp_pd + (size_t)i*pd_size;
+        memcpy(res_pd.data() + (size_t)i*pd_size, gpd, pd_size*sizeof(float));
+        for (int j = 0; j < pd_size; j++){
+          const bool gn = std::isnan(gpd[j]), en = std::isnan(epd[j]);
+          if (gn || en){ if(gn!=en) pd_nan++; continue; }
+          const double d = std::fabs((double)gpd[j]-(double)epd[j]);
+          if (d > pd_maxd) pd_maxd = d;
+        }
+        const float* gpk = got_pk.data() + (size_t)it->second*PEAK_ROW_FLOATS;
+        const float* epk = exp_pk + (size_t)i*7;
+        memcpy(res_pk.data() + (size_t)i*7, gpk, 7*sizeof(float));
+        for (int j = 0; j < 7; j++){
+          const bool gn = std::isnan(gpk[j]), en = std::isnan(epk[j]);
+          if (gn || en){ if(gn!=en) pk_nan++; continue; }
+          const double d = std::fabs((double)gpk[j]-(double)epk[j]);
+          if (d > pk_maxd) pk_maxd = d;
+        }
+      }
+      const int pd_fail = (pd_maxd>tol)      || (pd_nan>0) || (pd_missing>0);
+      const int pk_fail = (pk_maxd>peak_tol) || (pk_nan>0) || (pd_missing>0);
+      printf("it%d: compare 'corr_pd' (%d tiles, in-place FZ-norm): max|diff|=%g NaN-mismatch=%ld missing=%ld -> %s\n",
+             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 NaN-mismatch=%ld -> %s (peak_tol %g)\n",
+             k, n_pd, pk_maxd, pk_nan, pk_fail?"FAIL":"PASS", peak_tol);
+      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);
+      snprintf(nm,sizeof(nm),"corr_pd_it%d",k);
+      data.saveResult(nm, scratch, {n_pd, pd_size});
+      cudaFree(scratch);
+      cudaMalloc((void**)&scratch, res_pk.size()*sizeof(float));
+      cudaMemcpy(scratch, res_pk.data(), res_pk.size()*sizeof(float), cudaMemcpyHostToDevice);
+      snprintf(nm,sizeof(nm),"peaks_it%d",k);
+      data.saveResult(nm, scratch, {n_pd, 7});
+      cudaFree(scratch);
+    }
   }
   printf("%s: %s\n", argv[0], fail ? "FAIL" : "PASS");
   tp_proc_destroy(p); tp_destroy_module(m);

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 23 funcs), // By Claude on 07/12/2026: +4 tp_consolidate kernels
+// extern "C": tp_create_module (NVRTC-compile kernels + cuLink(libcudadevrt) + load 24 funcs), // By Claude on 07/13/2026: +corr2D_peak_eig (rung 2)
 // tp_module_num_functions, tp_last_error, tp_destroy_module. By Claude on 2026-06-25.
 #include <cuda.h>
 #include <nvrtc.h>
@@ -21,6 +21,7 @@
 #include "tp_utils.h"   // copyalloc_kernel_gpu / copyalloc_pointers_gpu / copyalloc_image_gpu / alloc_kernel_gpu
 #include "tp_files.h"   // readFloatsFromFile / readAllFloatsFromFile
 #include "tp_consolidate.h" // TD_TILE_FLOATS (host use in tp_proc_exec_consolidate) // By Claude on 07/12/2026
+#include "tp_peak.h"        // PEAK_ROW_FLOATS/PEAK_MAX_PIX (host use in tp_proc_exec_corr2d_peak) // By Claude on 07/13/2026
 
 static char g_err[8192] = "";
 static void seterr(const char* fmt, ...){ va_list ap; va_start(ap,fmt); vsnprintf(g_err,sizeof(g_err),fmt,ap); va_end(ap); }
@@ -38,13 +39,15 @@ static const char* TP_DEFINES =
   "#define LIST_TEXTURE_BIT 8\n#define TEXT_NTILE_SHIFT 9\n#define FAT_ZERO_WEIGHT 0.0001\n"
   "#define THREADS_DYNAMIC_BITS 5\n#define RBYRDIST_LEN 5001\n#define RBYRDIST_STEP 0.0004\n#define TILES_PER_BLOCK_GEOM 2\n";
 static const char* SRC_FILES[] = {"dtt8x8.h","dtt8x8.cu","geometry_correction.h","geometry_correction.cu","TileProcessor.h","TileProcessor.cu",
-  "tp_consolidate.h","tp_consolidate.cu"}; // TD consolidation chain (roadmap 2c) // By Claude on 07/12/2026
+  "tp_consolidate.h","tp_consolidate.cu", // TD consolidation chain (roadmap 2c) // By Claude on 07/12/2026
+  "tp_peak.h","tp_peak.cu"}; // peak measurement kernel (roadmap 3-C rung 2) // By Claude on 07/13/2026
 static const char* KERNELS[] = {
   "convert_direct","imclt_rbg_all","correlate2D","correlate2D_inter","corr2D_combine","corr2D_normalize",
   "textures_nonoverlap","generate_RBGA","clear_texture_list","mark_texture_tiles","mark_texture_neighbor_tiles",
   "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
+  "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
 static const int N_KERNELS = sizeof(KERNELS)/sizeof(KERNELS[0]);
 
 struct TpModule { CUcontext ctx; CUmodule mod; int nfun; };
@@ -611,6 +614,12 @@ struct TpProc {
   int cons_capacity;
   float *cons_ftasks0, *cons_ftasks1, *cons_pairs, *cons_tasks_out, *cons_td_sens, *cons_td_avg;
   int   *cons_counters; // [0] pairs, [1] surviving tiles, [2] misaligned
+  // peak measurement (tp_peak.cu, roadmap 3-C rung 2) — lazily allocated in
+  // tp_proc_exec_corr2d_peak; gpu_peak_debias uploaded once via
+  // tp_proc_set_peak_debias (len 0 = off -> kernel gets NULL). By Claude on 07/13/2026.
+  int    peak_capacity;   // corr tiles gpu_peaks is sized for (PEAK_ROW_FLOATS each)
+  float *gpu_peaks;
+  float *gpu_peak_debias; int peak_debias_len;
 };
 
 extern "C" {
@@ -630,6 +639,7 @@ TpProc* tp_proc_create(TpModule* m){
   p->gpu_tex_indices=p->gpu_tex_pnum=nullptr; p->last_num_tex_tiles=0;
   p->cons_capacity=0; p->cons_ftasks0=p->cons_ftasks1=p->cons_pairs=p->cons_tasks_out=p->cons_td_sens=p->cons_td_avg=nullptr;
   p->cons_counters=nullptr; // By Claude on 07/12/2026
+  p->peak_capacity=0; p->gpu_peaks=nullptr; p->gpu_peak_debias=nullptr; p->peak_debias_len=0; // By Claude on 07/13/2026
   return p;
 }
 
@@ -955,6 +965,56 @@ int tp_proc_exec_corr2d_normalize(TpProc* p, int combo, double fat_zero, int cor
     return launch1(f,1,1,1,1,1,1,a,"corr2D_normalize_pair")?-3:0;
   }
 }
+
+// --- peak measurement (tp_peak.cu corr2D_peak_eig, roadmap 3-C rung 2) --- // By Claude on 07/13/2026
+// Optional envelope de-bias array (CuasPoseRT.windowAutocorrDebias, (2*corr_rad+1)^2
+// floats): upload once per config; len 0 (or null data) clears it -> kernel gets NULL.
+int tp_proc_set_peak_debias(TpProc* p, const float* data, int len){
+  if(!p){ seterr("set_peak_debias: null"); return -1; } cuCtxSetCurrent(p->mod->ctx);
+  if(!data || (len<=0)){ cudaFree(p->gpu_peak_debias); p->gpu_peak_debias=nullptr; p->peak_debias_len=0; return 0; }
+  if(len>PEAK_MAX_PIX){ seterr("set_peak_debias: len %d > %d",len,PEAK_MAX_PIX); return -2; }
+  if(len!=p->peak_debias_len){ cudaFree(p->gpu_peak_debias); p->gpu_peak_debias=nullptr;
+    if(cudaMalloc((void**)&p->gpu_peak_debias,(size_t)len*sizeof(float))!=cudaSuccess){
+      p->peak_debias_len=0; seterr("set_peak_debias: alloc failed"); return -3; } }
+  if(cudaMemcpy(p->gpu_peak_debias,data,(size_t)len*sizeof(float),cudaMemcpyHostToDevice)!=cudaSuccess){
+    seterr("set_peak_debias: HtoD failed"); return -3; }
+  p->peak_debias_len=len;
+  return 0;
+}
+// Run the frozen 15-arg getMaxXYCmEig contract on the pixel-domain correlation buffer
+// (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.
+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){
+  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; }
+  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){
+    cudaFree(p->gpu_peaks); p->gpu_peaks=nullptr;
+    if(cudaMalloc((void**)&p->gpu_peaks,(size_t)nct*PEAK_ROW_FLOATS*sizeof(float))!=cudaSuccess){
+      p->peak_capacity=0; seterr("exec_corr2d_peak: peaks alloc failed (%d tiles)",nct); return -4; }
+    p->peak_capacity=nct;
+  }
+  size_t sp=p->dstride_corr/sizeof(float); int crad=p->corr_out_rad;
+  void* a[]={ &nct,&p->gpu_corr_indices,&pair_select,&sp,&p->gpu_corrs,&crad,
+              &abs_min,&rel_min,&min_peak,&eig_sub_frac,
+              &refine,&eig_sub_frac1,&scale_axes,&inc_axes,
+              &p->gpu_peak_debias,&p->gpu_peaks }; // gpu_peak_debias nullptr = de-bias off
+  if(launch1(f,nct,1,1,256,1,1,a,"corr2D_peak_eig")) return -5; // 256 = PEAK_THREADS (tp_peak.cu)
+  return nct;
+}
+// D2H the peak rows (last_num_corr_tiles x PEAK_ROW_FLOATS floats). Returns the count.
+int tp_proc_get_peaks(TpProc* p, float* out){
+  if(!p||!p->gpu_peaks){ seterr("get_peaks: no peaks buffer"); return -1; } cuCtxSetCurrent(p->mod->ctx);
+  int nct=p->last_num_corr_tiles;
+  if(cudaMemcpy(out,p->gpu_peaks,(size_t)nct*PEAK_ROW_FLOATS*sizeof(float),cudaMemcpyDeviceToHost)!=cudaSuccess){
+    seterr("get_peaks: DtoH failed"); return -2; }
+  return nct;
+}
 // de-pitch per-pair gpu_corrs (last_num_corr_tiles x (2*corr_rad+1)^2). returns count.
 int tp_proc_get_corr2d(TpProc* p, float* out, int corr_rad){ if(!p||!p->have_corr)return -1; cuCtxSetCurrent(p->mod->ctx);
   int per=(2*corr_rad+1)*(2*corr_rad+1);
@@ -1161,6 +1221,7 @@ void tp_proc_destroy(TpProc* p){
     cudaFree(p->gpu_tex_pnum); cudaFree(p->gpu_tex_color_weights); cudaFree(p->gpu_tex_rbga_params); }
   cudaFree(p->cons_ftasks0); cudaFree(p->cons_ftasks1); cudaFree(p->cons_pairs); cudaFree(p->cons_tasks_out); // By Claude on 07/12/2026
   cudaFree(p->cons_td_sens); cudaFree(p->cons_td_avg); cudaFree(p->cons_counters); // cudaFree(nullptr) is a no-op // By Claude on 07/12/2026
+  cudaFree(p->gpu_peaks); cudaFree(p->gpu_peak_debias); // By Claude on 07/13/2026
   delete p;
 }
 

src/tests/build_tests.sh

@@ -25,7 +25,7 @@ fi
 
 # Kernel + host sources shared by every test ("load all kernels"):
 COMMON="$SRC/TileProcessor.cu $SRC/dtt8x8.cu $SRC/geometry_correction.cu \
-        $SRC/tp_consolidate.cu $SRC/tp_utils.cu $SRC/tp_files.cu tp_test_data.cu"
+        $SRC/tp_consolidate.cu $SRC/tp_peak.cu $SRC/tp_utils.cu $SRC/tp_files.cu tp_test_data.cu"
 
 TESTS="${@:-$(ls test_*.cu | sed 's/\.cu$//')}"
 for t in $TESTS; do

src/tests/test_peak.cu

+/*
+ * test_peak.cu
+ *
+ * Standalone test of corr2D_peak_eig (tp_peak.cu, roadmap 3-C rung 2) against
+ * the SYNTHETIC numpy reference case (src/tests/make_peak_case.py): float64
+ * port of the frozen 15-arg getMaxXYCmEig contract, same summation order, so
+ * the double-math kernel must agree to rounding. NaN rows (rejects) compare
+ * NaN==NaN; interleaved non-sum slots verify the pair_select skip.
+ *
+ * This is the KERNEL-LOGIC gate that runs with no Java/GPU-export dependency
+ * (flood-fill connectivity, argmax tie-break, eigen, recentering, de-bias);
+ * the REAL-DATA gates are test_pose_corr / test_pose_corr_jna checkpoint 2
+ * on a rung-2 pose_corr case export.
+ *
+ * Usage: test_peak --data testdata/synth4/peak [--tol 1e-9] [--list]
+ *
+ *  Created on: Jul 13, 2026
+ *      Author: Claude (Anthropic), for Elphel
+ */
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+#include <string>
+#include <vector>
+
+#include <cuda_runtime.h>
+#include <helper_cuda.h>
+
+#include "tp_defines.h"            // CORR_NTILE_SHIFT
+#include "tp_peak.h"
+#include "tests/tp_test_data.h"
+#include "tests/tp_gpu_buf.h"
+
+int main(int argc, char **argv)
+{
+	std::string data_dir = "testdata/synth4/peak";
+	float tol = 1e-9f;             // double-math kernel vs the float64 numpy reference
+	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_corr  = data.paramInt("num_corr_tiles", 0);
+	const int corr_rad  = data.paramInt("corr_radius", 7);
+	const int pd_size   = (2 * corr_rad + 1) * (2 * corr_rad + 1);
+	const float p_abs_min   = (float) data.param("eig_min_abs",   0.0);
+	const float p_rel_min   = (float) data.param("eig_min_rel",   0.0);
+	const float p_min_peak  = (float) data.param("min_str_sum",   0.0);
+	const float p_sub_frac  = (float) data.param("eig_sub_frac",  0.0);
+	const int   p_recenter  =         data.paramInt("pose_recenter", 0);
+	const float p_sub_frac1 = (float) data.param("eig_sub_frac1", 0.0);
+	const float p_scale_ax  = (float) data.param("eig_scale_axes",1.2);
+	const float p_inc_ax    = (float) data.param("eig_inc_axes",  1.0);
+	printf("tiles=%d rad=%d abs/rel/min=%g/%g/%g sub=%g/%g recenter=%d axes=%gx+%g debias=%s\n",
+			num_corr, corr_rad, p_abs_min, p_rel_min, p_min_peak, p_sub_frac, p_sub_frac1,
+			p_recenter, p_scale_ax, p_inc_ax, data.has("debias") ? "ON" : "off");
+
+	// PD tiles into a pitched buffer (the kernel takes a stride, as production)
+	GpuPitchedBuf<float> gpu_pd(pd_size, num_corr);
+	gpu_pd.upload(data.hostFloat("corr_pd"));
+	GpuBuf<int> gpu_idx(num_corr);
+	gpu_idx.upload(data.hostInt("corr_indices"));
+	GpuBuf<float> gpu_peaks((size_t) num_corr * PEAK_ROW_FLOATS);
+	GpuBuf<float> gpu_debias(pd_size);
+	float * d_debias = nullptr;
+	if (data.has("debias")) {
+		gpu_debias.upload(data.hostFloat("debias"));
+		d_debias = gpu_debias.dev();
+	}
+
+	corr2D_peak_eig<<<num_corr, 256>>>(
+			num_corr, gpu_idx.dev(),
+			0xff,                                   // sum slots only (interleaved 0x00 slots skipped)
+			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> got = gpu_peaks.download();
+
+	const int * pd_idx = data.hostInt("expected_pd_indices");
+	const int   n_exp  = (int) data.elements("expected_pd_indices");
+	const float * exp_pk = data.hostFloat("expected_peaks");
+	const int * idx = data.hostInt("corr_indices");
+	int fail = 0;
+	// non-sum slots must be skipped (NaN + valid 0)
+	for (int i = 0; i < num_corr; i++) {
+		if ((idx[i] & ((1 << CORR_NTILE_SHIFT) - 1)) != 0xff) {
+			const float * g = got.data() + (size_t) i * PEAK_ROW_FLOATS;
+			if (!std::isnan(g[0]) || (g[PEAK_ROW_FLOATS - 1] != 0.0f)) {
+				printf("slot %d (pair 0x%02x): NOT skipped -> FAIL\n", i, idx[i] & 0xff);
+				fail = 1;
+			}
+		}
+	}
+	double maxd = 0;
+	long nan_mm = 0;
+	for (int i = 0; i < n_exp; i++) {
+		// find the sum slot of this dense tile
+		const int packed = (pd_idx[i] << CORR_NTILE_SHIFT) | 0xff;
+		int row = -1;
+		for (int j = 0; j < num_corr; j++) if (idx[j] == packed) { row = j; break; }
+		if (row < 0) { printf("tile %d: sum slot missing -> FAIL\n", pd_idx[i]); fail = 1; continue; }
+		const float * g = got.data() + (size_t) row * PEAK_ROW_FLOATS;
+		const float * e = exp_pk + (size_t) i * 7;
+		double tile_maxd = 0;
+		for (int j = 0; j < 7; j++) {
+			const bool gn = std::isnan(g[j]), en = std::isnan(e[j]);
+			if (gn || en) { if (gn != en) nan_mm++; continue; }
+			const double d = std::fabs((double) g[j] - (double) e[j]);
+			if (d > tile_maxd) tile_maxd = d;
+		}
+		if (tile_maxd > maxd) maxd = tile_maxd;
+		printf("tile %2d: got {%9.6f,%9.6f,%8.5f | %8.5f,%8.5f,%9.4f,%9.4f} v=%g max|diff|=%g%s\n",
+				pd_idx[i], g[0], g[1], g[2], g[3], g[4], g[5], g[6], g[7], tile_maxd,
+				(tile_maxd > tol) ? " <- OVER TOL" : "");
+	}
+	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");
+	printf("%s: %s\n", argv[0], fail ? "FAIL" : "PASS");
+	return fail ? EXIT_FAILURE : EXIT_SUCCESS;
+}

src/tests/test_pose_corr.cu

@@ -45,7 +45,19 @@
  * --tol (~1 ULP at scale) or against goldens blessed from THIS binary;
  * use test_pose_corr_jna for the production-equivalence check.
  *
- * Usage: test_pose_corr --data <dir> [--tol X] [--list]
+ * CHECKPOINT 2 (rung 2, added 07/13/2026), when the case carries the
+ * expected_pd_indices_it<k> buffers (re-export with the rung-2 Java):
+ *   corr2D_normalize IN PLACE on the inter TD buffer (weights NULL == the
+ *   production uniform 1.0) -> compare pixel-domain tiles vs expected_corr_pd
+ *   (dense-tile-indexed, matched via packed (tile<<8)|0xff, --tol tier) ->
+ *   corr2D_peak_eig (the frozen 15-arg getMaxXYCmEig contract, params from
+ *   the manifest incl. optional debias array) -> compare {dx,dy,str,eigen}
+ *   vs the double-precision Java oracle rows at --peak-tol (default 1e-6:
+ *   the kernel computes in double but offline-nvcc FMA/codegen and the float
+ *   PD input make bit-parity with Java unattainable - the tol-0 verdict is
+ *   the NVRTC replay, test_pose_corr_jna). Old cases skip with a note.
+ *
+ * Usage: test_pose_corr --data <dir> [--tol X] [--peak-tol X] [--list]
  *
  *  Created on: Jul 13, 2026
  *      Author: Claude (Anthropic), for Elphel
@@ -67,6 +79,7 @@
 #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 "tp_peak.h"               // corr2D_peak_eig (checkpoint 2, rung 2) // By Claude on 07/13/2026
 #include "tests/tp_test_data.h"
 #include "tests/tp_gpu_buf.h"
 
@@ -121,17 +134,19 @@ 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
+	float peak_tol  = 1e-6f;  // Java-double-oracle tier for the peak kernel (header note)
 	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], "--peak-tol") && (i + 1 < argc)) peak_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]);
+			fprintf(stderr, "Usage: %s --data <dir> [--tol X] [--peak-tol X] [--list]\n", argv[0]);
 			return EXIT_FAILURE;
 		}
 	}
-	printf("%s: data=%s tol=%g\n", argv[0], data_dir.c_str(), tol);
+	printf("%s: data=%s tol=%g peak_tol=%g\n", argv[0], data_dir.c_str(), tol, peak_tol);
 
 	findCudaDevice(argc, (const char **) argv);
 	TpTestData data(data_dir);
@@ -252,6 +267,37 @@ int main(int argc, char **argv)
 	GpuBuf<float> gpu_result_scratch(max_corr * chunk);      // reordered result upload for saveResult
 	const size_t corr_stride = gpu_corrs_td.pitchFloats();
 
+	// ---- checkpoint 2 (rung 2) statics: FZ-norm in place + frozen peak contract ----
+	// By Claude on 07/13/2026
+	const int  corr_rad  = DTT_SIZE - 1;                                 // 7 -> 15x15
+	const int  pd_size   = (2 * corr_rad + 1) * (2 * corr_rad + 1);      // 225
+	const bool have_ckpt2 = data.has("expected_pd_indices_it0");
+	const float fz_inter = (float) data.param("fz_inter", 10000.0);
+	const float fz2 = (float) ((double) fz_inter * (double) fz_inter);   // GpuQuad: (float)(fat_zero*fat_zero)
+	const float p_abs_min   = (float) data.param("eig_min_abs",   0.0);
+	const float p_rel_min   = (float) data.param("eig_min_rel",   0.0);
+	const float p_min_peak  = (float) data.param("min_str_sum",   0.0);
+	const float p_sub_frac  = (float) data.param("eig_sub_frac",  0.0);
+	const int   p_recenter  =         data.paramInt("pose_recenter", 0);
+	const float p_sub_frac1 = (float) data.param("eig_sub_frac1", 0.0);
+	const float p_scale_ax  = (float) data.param("eig_scale_axes",1.2);
+	const float p_inc_ax    = (float) data.param("eig_inc_axes",  1.0);
+	GpuPitchedBuf<float> gpu_corrs_pd(pd_size, (int) max_corr);
+	GpuBuf<float> gpu_peaks(max_corr * PEAK_ROW_FLOATS);
+	GpuBuf<float> gpu_debias(pd_size);
+	float * d_debias = nullptr;                                          // NULL = de-bias off (the contract)
+	if (have_ckpt2 && data.has("debias")) {
+		gpu_debias.upload(data.hostFloat("debias"));
+		d_debias = gpu_debias.dev();
+	}
+	if (have_ckpt2) {
+		printf("checkpoint 2: fz=%g recenter=%d debias=%s peak params abs/rel/min=%g/%g/%g sub=%g/%g axes=%gx+%g\n",
+				fz_inter, p_recenter, d_debias ? "ON" : "off",
+				p_abs_min, p_rel_min, p_min_peak, p_sub_frac, p_sub_frac1, p_scale_ax, p_inc_ax);
+	} else {
+		printf("NOTE: checkpoint-2 buffers absent (pre-rung-2 case vintage) - FZ-norm/peak replay skipped\n");
+	}
+
 	int fail = 0;
 	const int IDX_THREADS = 256;
 	for (int k = 0; k < num_iters; k++) {
@@ -420,6 +466,88 @@ int main(int argc, char **argv)
 				reordered.size() * sizeof(float), cudaMemcpyHostToDevice));
 		snprintf(nm, sizeof(nm), "corr_td_it%d", k);
 		data.saveResult(nm, gpu_result_scratch.dev(), {n_exp, chunk});
+
+		// ---- CHECKPOINT 2 (rung 2): FZ-norm IN PLACE + peak kernel ----
+		// By Claude on 07/13/2026
+		if (have_ckpt2) {
+			snprintf(nm, sizeof(nm), "expected_pd_indices_it%d", k);
+			const int * pd_idx = data.hostInt(nm);
+			const int   n_pd   = (int) data.elements(nm);
+			snprintf(nm, sizeof(nm), "expected_corr_pd_it%d", k);
+			const float * exp_pd = data.hostFloat(nm);
+			snprintf(nm, sizeof(nm), "expected_peaks_it%d", k);
+			const float * exp_pk = data.hostFloat(nm);
+			// normalize ALL resident slots where they sit (per-tile independent math;
+			// NULL weights == the production uniform 1.0 - see TDCorrTile.convertTDtoPDInPlace)
+			corr2D_normalize<<<1, 1>>>(
+					num_corr, corr_stride, gpu_corrs_td.dev(),
+					(float *) NULL,
+					gpu_corrs_pd.pitchFloats(), gpu_corrs_pd.dev(),
+					fz2, corr_rad);
+			checkCudaErrors(cudaGetLastError());
+			checkCudaErrors(cudaDeviceSynchronize());
+			corr2D_peak_eig<<<num_corr, 256>>>(
+					num_corr, gpu_corr_indices.dev(),
+					0xff,                                   // lean processes the sum slot only
+					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_pd((size_t) num_corr * pd_size);
+			checkCudaErrors(cudaMemcpy2D(got_pd.data(), (size_t) pd_size * sizeof(float),
+					gpu_corrs_pd.dev(), gpu_corrs_pd.pitchBytes(),
+					(size_t) pd_size * sizeof(float), num_corr, cudaMemcpyDeviceToHost));
+			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));
+			// 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;
+			long pd_nan = 0, pk_nan = 0, pd_missing = 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);
+				if (it == got_by_index.end()) { pd_missing++; continue; }
+				const float * gpd = got_pd.data() + (size_t) it->second * pd_size;
+				const float * epd = exp_pd + (size_t) i * pd_size;
+				memcpy(res_pd.data() + (size_t) i * pd_size, gpd, pd_size * sizeof(float));
+				for (int j = 0; j < pd_size; j++) {
+					const bool gn = std::isnan(gpd[j]), en = std::isnan(epd[j]);
+					if (gn || en) { if (gn != en) pd_nan++; continue; }
+					const double d = std::fabs((double) gpd[j] - (double) epd[j]);
+					if (d > pd_maxd) pd_maxd = d;
+				}
+				const float * gpk = got_pk.data() + (size_t) it->second * PEAK_ROW_FLOATS;
+				const float * epk = exp_pk + (size_t) i * 7;
+				memcpy(res_pk.data() + (size_t) i * 7, gpk, 7 * sizeof(float));
+				for (int j = 0; j < 7; j++) {
+					const bool gn = std::isnan(gpk[j]), en = std::isnan(epk[j]);
+					if (gn || en) { if (gn != en) pk_nan++; continue; }
+					const double d = std::fabs((double) gpk[j] - (double) epk[j]);
+					if (d > pk_maxd) pk_maxd = d;
+				}
+			}
+			const int pd_fail = (pd_maxd > tol)      || (pd_nan > 0) || (pd_missing > 0);
+			const int pk_fail = (pk_maxd > peak_tol) || (pk_nan > 0) || (pd_missing > 0);
+			printf("it%d: compare 'corr_pd' (%d tiles, in-place FZ-norm): max|diff|=%g NaN-mismatch=%ld "
+					"missing=%ld -> %s\n", 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 "
+					"NaN-mismatch=%ld -> %s (peak_tol %g)\n",
+					k, n_pd, pk_maxd, pk_nan, pk_fail ? "FAIL" : "PASS", peak_tol);
+			if (pd_fail || pk_fail) fail = 1;
+			checkCudaErrors(cudaMemcpy(gpu_result_scratch.dev(), res_pd.data(),
+					res_pd.size() * sizeof(float), cudaMemcpyHostToDevice));
+			snprintf(nm, sizeof(nm), "corr_pd_it%d", k);
+			data.saveResult(nm, gpu_result_scratch.dev(), {n_pd, pd_size});
+			checkCudaErrors(cudaMemcpy(gpu_result_scratch.dev(), res_pk.data(),
+					res_pk.size() * sizeof(float), cudaMemcpyHostToDevice));
+			snprintf(nm, sizeof(nm), "peaks_it%d", k);
+			data.saveResult(nm, gpu_result_scratch.dev(), {n_pd, 7});
+		}
 	}
 	printf("%s: %s\n", argv[0], fail ? "FAIL" : "PASS");
 	return fail ? EXIT_FAILURE : EXIT_SUCCESS;

src/tp_peak.h

+/*
+ * tp_peak.h
+ *
+ * Per-tile correlation peak measurement kernel — the GPU side of the lean
+ * pose-chain CHECKPOINT 2 second half (roadmap 3-C rung 2).
+ *
+ * FROZEN CONTRACT (Andrey's ruling 07/13/2026): the 15-arg
+ * Correlation2d.getMaxXYCmEig superset —
+ *   argmax (first occurrence in linescan order) -> weak-peak reject
+ *   -> 8-connected region above min(abs_min, rel_min*max) grown from the
+ *      argmax seed -> pedestal-subtracted centroid + covariance
+ *   -> closed-form 2x2 eigen (getEigen2x2 port, eig_fast2x2 == the contract)
+ *   -> OPTIONAL recentering passes (refine > 0): cosine-mask window over the
+ *      eigen-ellipse (scale_axes*sqrt(lambda)+inc_axes half-axes),
+ *      re-centroid each pass
+ *   -> OPTIONAL envelope de-bias: gpu_debias (NULL = off) is multiplied into
+ *      the tile before ANY peak math (the window-autocorr envelope
+ *      under-reports offsets by ~0.75-0.9x; array from
+ *      CuasPoseRT.windowAutocorrDebias, single Java source).
+ * 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).
+ *
+ * Java oracle: Correlation2d.getMaxXYCmEig (15-arg, eig_fast2x2=true);
+ * case export: PoseCorrExport (curt.kernel_test = pose_corr, checkpoint 2).
+ *
+ *  Created on: Jul 13, 2026
+ *      Author: Claude (Anthropic), for Elphel
+ */
+
+#ifndef SRC_TP_PEAK_H_
+#define SRC_TP_PEAK_H_
+
+#ifndef JCUDA // NVRTC concatenation (tp_jna.cpp TP_DEFINES) has no include paths
+#include <cuda_runtime.h> // __global__ in the prototype below, also for non-nvcc parsers
+#endif
+
+// correlation tile side bound: corr_radius <= DTT_SIZE-1 = 7 -> 15x15
+#define PEAK_MAX_WIDTH  15
+#define PEAK_MAX_PIX    (PEAK_MAX_WIDTH * PEAK_MAX_WIDTH)
+// output row stride (floats): {dx, dy, strength, e0x, e0y, lambda0, lambda1, valid}
+#define PEAK_ROW_FLOATS 8
+
+extern "C" __global__ void corr2D_peak_eig(
+		int           num_corr_tiles,   // resident correlation tiles (one block each)
+		const int   * gpu_corr_indices, // packed (tile << CORR_NTILE_SHIFT) | pair (0xff = sum)
+		int           pair_select,      // process only slots with this pair code (lean: 0xff); -1 = all
+		const size_t  corr_stride,      // in floats: pixel-domain correlation buffer stride
+		const float * gpu_corrs,        // FZ-normalized pixel-domain correlation tiles
+		int           corr_radius,      // 7 -> 15x15
+		float         abs_min,          // imp.eig_min_abs
+		float         rel_min,          // imp.eig_min_rel
+		float         min_peak,         // imp.min_str_sum
+		float         eig_sub_frac,     // pedestal subtract fraction (centroid/covariance)
+		int           refine,           // recentering passes (curt.pose_recenter; 0 = single-pass)
+		float         eig_sub_frac1,    // pedestal subtract during recentering
+		float         scale_axes,       // eigen-ellipse half-axis scale (imp.eig_scale_axes)
+		float         inc_axes,         // eigen-ellipse half-axis increment, pix (imp.eig_inc_axes)
+		const float * gpu_debias,       // NULL or [width*width] multiplicative envelope de-bias
+		float       * gpu_peaks);       // [num_corr_tiles][PEAK_ROW_FLOATS]; rejected/skipped ->
+		                                // 7x NaN + valid 0.0f; measured -> values + valid 1.0f
+
+#endif /* SRC_TP_PEAK_H_ */