Commit 814f38e0 authored by Andrey Filippov's avatar Andrey Filippov

CLAUDE: rung C1 - resident prepareLMA wired: measured inputs never leave the GPU

Design 3-A4i rung C1 (native tile_processor_gpu 7db7cf6):
- IntersceneLma: PoseLmaPrepareProvider + residentPrepare() (mirrors the
  getFxDerivs head: ERS pokes + setupERS + camera/centers capture; per-cycle
  H2D = camera state + pose vector + centers + 6 policy floats). Production
  lean prepare skips setSamplesWeights / both fx passes / WJtJ+reg+normalize /
  y-build entirely; weights/y stay GPU-resident; initial RMS seeds from the
  first resident step's packed[19,20] (kills the runLma first-step
  re-linearization = C0's 'LMA CPU remainder'). Invalid/unavailable falls
  back loudly to the legacy Java prepare.
- One-shot oracle (pose_lma_debug>=1, first prepare): the capture cycle runs
  FULLY legacy, then the resident prepare captures its buffers and the new
  IntersceneLmaFloat.prepareResidentOracle (serial float clone) must match
  BIT-EXACTLY ('resident CUDA vs Java-float PREPARE: ... mismatches=0').
- lmaStep: prepared-resident steps pass NULL weights/y/eigen (no per-step
  H2D); first-block linearization skipped when prepared.
- CuasPoseRT: provider wiring (slot 0xff + CORR_NTILE_SHIFT = host policy),
  new marker 'resident CUDA prepareLMA active'; cycle_rms_meas moved after
  runLma reading getInitialRms() (same value on both paths).
- Stage0 kernel count 33 -> 36 (3 new prepare kernels).
Gates: mvn package + test PASS; Stage0 36/36; full native suite in
tile_processor_gpu 7db7cf6 (direct prepare BIT-EXACT, sanitizer 0 errors).
Co-Authored-By: 's avatarClaude Fable 5 <noreply@anthropic.com>
parent b6a074a1
......@@ -870,6 +870,7 @@ public class CuasPoseRT {
private static boolean pose_fx_path_reported = false;
// once-per-program console note for resident float preparation/reduction/candidate
private static boolean pose_lma_step_path_reported = false;
private static boolean pose_prepare_path_reported = false; // rung C1 // By Claude on 07/17/2026
// one raw-pixel/J comparison of the float oracle against the current double geometry
private static boolean float_jacobian_oracle_reported = false;
......@@ -1081,6 +1082,27 @@ public class CuasPoseRT {
}
return result;
}); // production candidate: resident raw fx/J -> parallel prepare/reduce -> fixed 3x3 tail // By Codex on 07/15/2026
// Rung C1 (design 2026-07-16): resident prepareLMA - assemble/normalize on the
// GPU from the resident peak buffers; per-cycle H2D = camera state + pose
// vector + centers + policy scalars; D2H = the 8-float compact result. The
// packed-index slot/shift are CUAS host policy. By Claude on 07/17/2026.
intersceneLma.setPoseLmaPrepareProvider((refCam, sceneCam, refXyz, refAtr,
sceneXyz, sceneAtr, centers, useEigen, eigMinSqrt, eigMaxSqrt,
minConfidence, sameWeights, pullTargets, paramRegweights, capture) -> {
final IntersceneLmaFloat.PrepareResidentResult result =
lmaGpu.execPoseLmaPrepareResident(
refCam, sceneCam, refXyz, refAtr, sceneXyz, sceneAtr, centers,
0xff, GPUTileProcessor.CORR_NTILE_SHIFT,
useEigen, eigMinSqrt, eigMaxSqrt, minConfidence, sameWeights,
pullTargets, paramRegweights, capture);
if ((result != null) && !pose_prepare_path_reported) {
pose_prepare_path_reported = true;
System.out.println("CuasPoseRT: resident CUDA prepareLMA active "+
"(assemble+normalize on-device from resident peaks; steps run with "+
"no per-step H2D; initial RMS from the first step)");
}
return result;
});
// Rung C0 (design 2026-07-16): route IntersceneLma sub-stage samples into the
// RT profiler - attributes prepareLMA's budget and the cross-cutting fx
// boundaries (JNA roundtrip, setupERS) before C1 goes resident. Null profiler
......@@ -1216,7 +1238,11 @@ public class CuasPoseRT {
cuda_comparison.format("CUDA", "float"));
}
}
cycle_rms_meas[nlma] = intersceneLma.getLastRms()[1]; // re-measured pure RMS at this cycle's pose // By Claude on 07/12/2026
// cycle_rms_meas moved after runLma (rung C1): the resident prepare defers
// the initial RMS to the first step's packed[19,20] - getInitialRms() holds
// the same "re-measured at this cycle's pose" value on BOTH paths (the
// legacy prepare seeds it, the resident first step seeds it).
// By Claude on 07/17/2026 (was getLastRms()[1] here, before runLma).
// pose_lma_debug >= 2: per-inner-step solver lines; header attributes them to
// their outer cycle (the solver lines carry no scene/cycle context).
// By Claude on 07/12/2026, from Andrey's request.
......@@ -1252,6 +1278,7 @@ public class CuasPoseRT {
last_diff_atr = diffs_atr[0]; // By Claude on 07/12/2026
cycle_datr[nlma] = diffs_atr[0]; // By Claude on 07/12/2026
cycle_inner[nlma] = lmaResult; // inner runLma iterations this cycle // By Claude on 07/12/2026
cycle_rms_meas[nlma] = intersceneLma.getInitialRms()[1]; // re-measured pure RMS at this cycle's pose (moved after runLma, rung C1) // By Claude on 07/17/2026
cycle_rms_solved[nlma] = intersceneLma.getLastRms()[1]; // model-predicted pure RMS after the solve // By Claude on 07/12/2026
ncyc = nlma + 1; // By Claude on 07/12/2026
if (rtProfile != null) rtProfile.addElapsed(RtPoseProfile.POSE_UPDATE, profileStart);
......
......@@ -2955,6 +2955,33 @@ public class GpuQuad{ // quad camera description
return null;
}
/**
* Rung C1 (design 2026-07-16): resident prepareLMA - assemble y/weights/eigen
* from the resident peak buffers, project fx/J at the current pose (no
* readback), normalize weights on-device. Base/JCuda returns null so the Java
* prepare remains the fallback. By Claude on 07/17/2026.
*/
public IntersceneLmaFloat.PrepareResidentResult execPoseLmaPrepareResident(
IntersceneLmaFloat.Camera reference,
IntersceneLmaFloat.Camera scene,
float [] reference_xyz,
float [] reference_atr,
float [] scene_xyz,
float [] scene_atr,
float [] centers,
int corr_slot,
int ntile_shift,
boolean use_eigen,
float eig_min_sqrt,
float eig_max_sqrt,
float min_confidence,
boolean same_weights,
float [] pull_targets,
float [] param_regweights,
boolean capture) {
return null;
}
public void execRBGA(
double [] color_weights,
boolean is_lwir,
......
......@@ -556,6 +556,69 @@ public class GpuQuadJna extends GpuQuad {
return unpackPoseLmaStep(packed);
}
// ---- rung C1: resident prepareLMA (design 2026-07-16) ---- // By Claude on 07/17/2026
@Override public IntersceneLmaFloat.PrepareResidentResult execPoseLmaPrepareResident(
final IntersceneLmaFloat.Camera reference,
final IntersceneLmaFloat.Camera scene,
final float[] reference_xyz,
final float[] reference_atr,
final float[] scene_xyz,
final float[] scene_atr,
final float[] centers,
final int corrSlot,
final int ntileShift,
final boolean useEigen,
final float eigMinSqrt,
final float eigMaxSqrt,
final float minConfidence,
final boolean sameWeights,
final float[] pullTargets,
final float[] paramRegweights,
final boolean capture) {
if ((reference == null) || (scene == null) || (centers == null) ||
((centers.length % 3) != 0) ||
(pullTargets == null) || (pullTargets.length != IntersceneLmaFloat.NUM_PARAMS) ||
(paramRegweights == null) || (paramRegweights.length != IntersceneLmaFloat.NUM_PARAMS)) {
throw new IllegalArgumentException("GpuQuadJna.execPoseLmaPrepareResident: inconsistent inputs");
}
final int numTiles = centers.length / 3;
final int numValues = IntersceneLmaFloat.NUM_COMPONENTS * numTiles +
IntersceneLmaFloat.NUM_PARAMS;
final float[] poseVectors = new float[12];
copyPoseVector3(reference_xyz, poseVectors, 0);
copyPoseVector3(reference_atr, poseVectors, 3);
copyPoseVector3(scene_xyz, poseVectors, 6);
copyPoseVector3(scene_atr, poseVectors, 9);
final float[] prepResult = new float[IntersceneLmaFloat.PREP_RESULT_FLOATS];
final float[] outY = capture ? new float[numValues] : null;
final float[] outWeights = capture ? new float[numValues] : null;
final float[] outEigen = capture ? new float[4 * numTiles] : null;
final int rc = lib.tp_proc_exec_pose_lma_prepare_resident(
proc,
poseCameraMeta(reference), reference.radial,
reference.rByRDist, reference.rByRDist.length, reference.ers, reference.ers.length,
poseCameraMeta(scene), scene.radial,
scene.rByRDist, scene.rByRDist.length, scene.ers, scene.ers.length,
poseVectors, centers,
numTiles, IntersceneLmaFloat.NUM_COMPONENTS,
corrSlot, ntileShift,
useEigen ? 1 : 0, eigMinSqrt, eigMaxSqrt,
minConfidence, sameWeights ? 1 : 0,
pullTargets, paramRegweights,
prepResult, capture ? 1 : 0, outY, outWeights, outEigen);
if (rc != 0) {
throw new IllegalStateException("GpuQuadJna.execPoseLmaPrepareResident rc="+rc+": "+
lib.tp_last_error());
}
return new IntersceneLmaFloat.PrepareResidentResult(
numTiles,
prepResult[0] != 0.0f,
(int) prepResult[1],
prepResult[2], prepResult[3], prepResult[4],
new float[] {prepResult[5], prepResult[6], prepResult[7]},
outY, outWeights, outEigen);
}
// ---- resident float pose preparation/reduction continuation (roadmap rung 3-A4e) ----
@Override public IntersceneLmaFloat.ResidentLmaStepResult execPoseLmaResidentStep(
final float lambda,
......@@ -569,9 +632,13 @@ public class GpuQuadJna extends GpuQuad {
final boolean capturePrepared) {
final int numValues = IntersceneLmaFloat.NUM_COMPONENTS * numTiles +
IntersceneLmaFloat.NUM_PARAMS;
if ((numTiles <= 0) || (weights == null) || (weights.length != numValues) ||
// Rung C1: all-null weights/y/eigen = use the resident buffers assembled by
// execPoseLmaPrepareResident (no per-step H2D). By Claude on 07/17/2026.
final boolean usePrepared = (weights == null) && (y == null) && (eigen == null);
if ((numTiles <= 0) ||
(!usePrepared && ((weights == null) || (weights.length != numValues) ||
(y == null) || (y.length != numValues) || (eigen == null) ||
(eigen.length != 4 * numTiles) || (vector == null) ||
(eigen.length != 4 * numTiles))) || (vector == null) ||
(vector.length != IntersceneLmaFloat.NUM_PARAMS) ||
!Float.isFinite(pureWeight) || !(pureWeight > 0.0f) ||
!Float.isFinite(rmsDiff)) {
......
......@@ -10,7 +10,7 @@ import com.sun.jna.Pointer;
* com.elphel.imagej.gpu.jna.Stage0 [kernel_src_dir] [libcudadevrt.a]
*/
public class Stage0 {
private static final int EXPECTED_KERNELS = 33; // 19 base + 2 conditioning + 4 consolidate + 2 peak + pose fx/J + 4 pose LMA + products
private static final int EXPECTED_KERNELS = 36; // 19 base + 2 conditioning + 4 consolidate + 2 peak + pose fx/J + 4 pose LMA + products + 3 resident prepare (rung C1) // By Claude on 07/17/2026
public static void main(String[] args) {
String srcdir = (args.length > 0) ? args[0] : "/home/elphel/git/tile_processor_gpu/src";
String devrt = (args.length > 1) ? args[1] : "/usr/local/cuda/targets/x86_64-linux/lib/libcudadevrt.a";
......
......@@ -162,6 +162,25 @@ public interface TpJna extends Library {
float[] eigen, float[] vector, float pureWeight,
float rmsDiff, int numTiles, int capturePrepared, float[] result,
float[] preparedJt, float[] preparedYmfx);
/** Rung C1 (design 2026-07-16): resident prepareLMA. Assembles y/weights/eigen
* from the RESIDENT peak/corr-index buffers, projects fx/J at the current pose
* (no readback) and normalizes weights on-device. prepResult layout:
* {valid, numMeas, sumWRaw, pureWeight, fullWeight, regW0, regW1, regW2}.
* outY/outWeights/outEigen nullable unless capture != 0 (one-shot oracle).
* By Claude on 07/17/2026. */
int tp_proc_exec_pose_lma_prepare_resident(Pointer proc,
float[] refMeta, float[] refRadial,
float[] refRbr, int refRbrLen, float[] refErs, int refErsLen,
float[] sceneMeta, float[] sceneRadial,
float[] sceneRbr, int sceneRbrLen, float[] sceneErs, int sceneErsLen,
float[] poseVectors, float[] centers,
int numTiles, int numComponents,
int corrSlot, int ntileShift,
int useEigen, float eigMinSqrt, float eigMaxSqrt,
float minConfidence, int sameWeights,
float[] pullTargets, float[] paramRegweights,
float[] prepResult, int capture,
float[] outY, float[] outWeights, float[] outEigen);
/** Deterministic double normal-equation products for the lean pose LMA.
* jt is parameter-major [numParams][numValues]; out is row-major H followed by b.
* Damping and the small solve remain in Java. 0 on success. */
......
......@@ -43,6 +43,7 @@ public final class IntersceneLmaFloat {
public static final int LMA_PRODUCTS = NUM_PARAMS * NUM_PARAMS + NUM_PARAMS;
public static final int LMA_RESULT_FLOATS = 1 + LMA_PRODUCTS + 2 * NUM_PARAMS;
public static final int LMA_RESIDENT_RESULT_FLOATS = LMA_RESULT_FLOATS + 6;
public static final int PREP_RESULT_FLOATS = 8; // rung C1 resident prepare compact result // By Claude on 07/17/2026
public static final float INFINITY_DISPARITY = 0.01f;
static final int ERS_XYZ = 0;
......@@ -213,6 +214,207 @@ public final class IntersceneLmaFloat {
}
/** Compact resident step plus optional prepared arrays for a one-shot oracle gate. */
/** Rung C1 (design 2026-07-16): compact result of the resident prepareLMA.
* capturedY/Weights/Eigen are non-null only on the one-shot oracle call.
* By Claude on 07/17/2026. */
public static final class PrepareResidentResult {
public final int numTiles;
public final boolean valid;
public final int numMeas;
public final float sumWRaw;
public final float pureWeight;
public final float fullWeight;
public final float [] regWeights;
public final float [] capturedY;
public final float [] capturedWeights;
public final float [] capturedEigen;
public PrepareResidentResult(
final int numTiles,
final boolean valid,
final int numMeas,
final float sumWRaw,
final float pureWeight,
final float fullWeight,
final float [] regWeights,
final float [] capturedY,
final float [] capturedWeights,
final float [] capturedEigen) {
this.numTiles = numTiles;
this.valid = valid;
this.numMeas = numMeas;
this.sumWRaw = sumWRaw;
this.pureWeight = pureWeight;
this.fullWeight = fullWeight;
this.regWeights = regWeights;
this.capturedY = capturedY;
this.capturedWeights = capturedWeights;
this.capturedEigen = capturedEigen;
}
}
/** Serial float clone of the GPU resident prepare - the bit-exact oracle for
* the captured buffers (assemble scatter is order-free; both tails are
* serial in ascending order on the GPU by design). By Claude on 07/17/2026. */
public static final class PrepareResidentOracle {
public final float [] y;
public final float [] weights;
public final float [] eigen;
public final int numMeas;
public final float sumWRaw;
public final float pureWeight;
public final float fullWeight;
public final float [] regWeights;
PrepareResidentOracle(final float [] y, final float [] weights, final float [] eigen,
final int numMeas, final float sumWRaw, final float pureWeight,
final float fullWeight, final float [] regWeights) {
this.y = y;
this.weights = weights;
this.eigen = eigen;
this.numMeas = numMeas;
this.sumWRaw = sumWRaw;
this.pureWeight = pureWeight;
this.fullWeight = fullWeight;
this.regWeights = regWeights;
}
}
public static PrepareResidentOracle prepareResidentOracle(
final double [][] vectorXYSDS, // exact float widenings of the peak dx/dy/str
final double [][] eigenSrc, // {e0x,e0y,l0,l1} per tile, or null = identity mode
final Result raw, // raw float fx/J/valid at the current pose (readback)
final float eigMinSqrt,
final float eigMaxSqrt,
final float minConfidence,
final boolean sameWeights,
final float [] pullTargets,
final float [] paramRegweights) {
final int numTiles = raw.numTiles;
final int ns = NUM_COMPONENTS * numTiles;
final int numValues = ns + NUM_PARAMS;
final float [] y = new float [numValues];
final float [] weights = new float [numValues];
final float [] eigen = new float [4 * numTiles];
// assemble (per-tile, order-free scatter on the GPU)
final boolean useEigen = (eigenSrc != null);
for (int tile = 0; tile < numTiles; tile++) {
if (vectorXYSDS[tile] == null) continue; // unmeasured / kernel-rejected
final float str = (float) vectorXYSDS[tile][2];
if (!(str >= minConfidence)) continue; // NaN fails here too
float e0x = 0f, e0y = 0f, l0 = 0f, l1 = 0f;
if (useEigen) {
if (eigenSrc[tile] == null) continue; // null transform skips the tile
e0x = (float) eigenSrc[tile][0];
e0y = (float) eigenSrc[tile][1];
l0 = (float) eigenSrc[tile][2];
l1 = (float) eigenSrc[tile][3];
if (!Float.isFinite(e0x) || !Float.isFinite(e0y) ||
!Float.isFinite(l0) || !Float.isFinite(l1)) continue;
}
if (!(str > 0.0f)) continue;
final float w = sameWeights ? 1.0f : str;
final int vi = NUM_COMPONENTS * tile;
weights[vi] = w;
y[vi] = (float) vectorXYSDS[tile][0];
y[vi + 1] = (float) vectorXYSDS[tile][1];
final int ei = 4 * tile;
if (useEigen) {
final float k0inv = 1.0f / eigMaxSqrt;
final float k0 = Math.max(0.0f,
1.0f / Math.max(eigMinSqrt, (float) Math.sqrt(l0)) - k0inv);
final float k1 = Math.max(0.0f,
1.0f / Math.max(eigMinSqrt, (float) Math.sqrt(l1)) - k0inv);
eigen[ei] = k0 * e0x;
eigen[ei + 1] = k0 * e0y;
eigen[ei + 2] = k1 * e0y;
eigen[ei + 3] = k1 * -e0x;
} else {
eigen[ei] = 1.0f;
eigen[ei + 3] = 1.0f;
}
}
// tail 1: sample-weight normalization (serial ascending, GPU order)
float sum = 0.0f;
int numMeas = 0;
for (int tile = 0; tile < numTiles; tile++) {
final float w = weights[NUM_COMPONENTS * tile];
sum += w + w;
if (w > 0.0f) numMeas++;
}
final float sNorm = 1.0f / sum;
for (int tile = 0; tile < numTiles; tile++) {
final int vi = NUM_COMPONENTS * tile;
final float w = weights[vi] * sNorm;
weights[vi] = w;
weights[vi + 1] = w;
}
for (int par = 0; par < NUM_PARAMS; par++) {
y[ns + par] = pullTargets[par];
}
// tail 2: y += fx (measured+valid), H diagonal, reg weights, full normalize
for (int tile = 0; tile < numTiles; tile++) {
final int vi = NUM_COMPONENTS * tile;
if ((weights[vi] > 0.0f) && (raw.valid[tile] != 0)) {
y[vi] = y[vi] + raw.fx[vi];
y[vi + 1] = y[vi + 1] + raw.fx[vi + 1];
}
}
final float [] h = new float [NUM_PARAMS];
for (int tile = 0; tile < numTiles; tile++) {
final int vi = NUM_COMPONENTS * tile;
final float w = weights[vi];
if (!(w > 0.0f)) continue;
final boolean haveRaw = raw.valid[tile] != 0;
final int ei = 4 * tile;
for (int par = 0; par < NUM_PARAMS; par++) {
final int rawI = par * ns + vi;
final float jx = haveRaw ? raw.jt[rawI] : 0.0f;
final float jy = haveRaw ? raw.jt[rawI + 1] : 0.0f;
final float tj0 = eigen[ei] * jx + eigen[ei + 1] * jy;
final float tj1 = eigen[ei + 2] * jx + eigen[ei + 3] * jy;
h[par] = h[par] + w * (tj0 * tj0);
h[par] = h[par] + w * (tj1 * tj1);
}
}
final float [] regWeights = new float [NUM_PARAMS];
for (int par = 0; par < NUM_PARAMS; par++) {
regWeights[par] = (h[par] > 0.0f) ?
(paramRegweights[par] / (float) Math.sqrt(h[par])) : 0.0f;
weights[ns + par] = regWeights[par];
}
float pure = 0.0f;
for (int i = 0; i < ns; i++) pure += weights[i];
float full = pure;
for (int par = 0; par < NUM_PARAMS; par++) full += weights[ns + par];
final float scale = 1.0f / full;
for (int i = 0; i < numValues; i++) weights[i] *= scale;
return new PrepareResidentOracle(y, weights, eigen,
numMeas, sum, pure / full, full, regWeights);
}
/** One-shot comparison of the captured GPU prepare buffers vs the serial
* float clone: exact float mismatch counts (0 expected). By Claude on 07/17/2026. */
public static String comparePrepareResident(
final PrepareResidentResult gpu,
final PrepareResidentOracle oracle) {
int badY = 0, badW = 0, badE = 0;
for (int i = 0; i < oracle.y.length; i++) {
if (Float.floatToIntBits(gpu.capturedY[i]) != Float.floatToIntBits(oracle.y[i])) badY++;
}
for (int i = 0; i < oracle.weights.length; i++) {
if (Float.floatToIntBits(gpu.capturedWeights[i]) != Float.floatToIntBits(oracle.weights[i])) badW++;
}
for (int i = 0; i < oracle.eigen.length; i++) {
if (Float.floatToIntBits(gpu.capturedEigen[i]) != Float.floatToIntBits(oracle.eigen[i])) badE++;
}
return String.format(
"y mismatches=%d, weight mismatches=%d, eigen mismatches=%d; "+
"num_meas %d/%d, pure_weight %g/%g, full_weight %g/%g",
badY, badW, badE, gpu.numMeas, oracle.numMeas,
gpu.pureWeight, oracle.pureWeight, gpu.fullWeight, oracle.fullWeight);
}
public static final class ResidentLmaStepResult {
public final LmaStepResult step;
public final PreparedLma prepared;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment