Commit 27d944d7 authored by Andrey Filippov's avatar Andrey Filippov

CODEX: Gate float CUDA pose LMA candidate

Co-authored-by: 's avatarCodex <codex@elphel.com>
parent e824a17c
......@@ -836,6 +836,8 @@ public class CuasPoseRT {
private static boolean lma_products_path_reported = false;
// once-per-program console note for production getFxDerivs() on CUDA
private static boolean pose_fx_path_reported = false;
// once-per-program console note for the fixed float runLma candidate path
private static boolean pose_lma_step_path_reported = false;
// one raw-pixel/J comparison of the float oracle against the current double geometry
private static boolean float_jacobian_oracle_reported = false;
......@@ -1032,6 +1034,16 @@ public class CuasPoseRT {
}
return result;
}); // first production gate: replace only getFxDerivs per-tile work; callers remain unchanged // By Codex on 07/15/2026
intersceneLma.setPoseLmaStepProvider((lambda, weights, jt, ymfxWeighted, vector) -> {
final IntersceneLmaFloat.LmaStepResult result = lmaGpu.execPoseLmaStep(
lambda, weights, jt, ymfxWeighted, vector);
if ((result != null) && !pose_lma_step_path_reported) {
pose_lma_step_path_reported = true;
System.out.println("CuasPoseRT: CUDA float runLma candidate active "+
"(one thread, per-call H2D/D2H; Java RMS/acceptance retained)");
}
return result;
}); // validation rung: only fixed 3x3 candidate calculation moves to CUDA // By Codex on 07/15/2026
intersceneLma.setNormalEquationProvider((weights, jt, ymfxWeighted) -> {
final long profileStart = (rtProfile != null) ? rtProfile.start() : 0L;
final double [] products = lmaGpu.execLmaNormalProducts(weights, jt, ymfxWeighted);
......
......@@ -2923,6 +2923,20 @@ public class GpuQuad{ // quad camera description
return null;
}
/**
* Fixed three-angle float LMA candidate gate. The JNA backend computes raw
* H/b and the explicit damped 3x3 solve in one CUDA thread; base/JCuda returns
* null so the Java-double solver remains the fallback.
*/
public IntersceneLmaFloat.LmaStepResult execPoseLmaStep(
float lambda,
float [] weights,
float [] jt,
float [] ymfx_weighted,
float [] vector) {
return null;
}
public void execRBGA(
double [] color_weights,
boolean is_lwir,
......
package com.elphel.imagej.gpu.jna;
import java.util.Arrays;
import com.elphel.imagej.gpu.GpuQuad;
import com.elphel.imagej.gpu.GPUTileProcessor;
import com.elphel.imagej.gpu.TpTask;
......@@ -532,6 +534,34 @@ public class GpuQuadJna extends GpuQuad {
numTiles, fx, calculate_jacobian ? nativeJt : null, valid, numValid);
}
// ---- fixed float pose LMA candidate (tp_lma.cu, roadmap rung 3-A4d) ---- // By Codex on 07/15/2026
@Override public IntersceneLmaFloat.LmaStepResult execPoseLmaStep(
final float lambda,
final float[] weights,
final float[] jt,
final float[] ymfxWeighted,
final float[] vector) {
if ((weights == null) || (weights.length == 0) ||
(jt == null) || (jt.length != IntersceneLmaFloat.NUM_PARAMS * weights.length) ||
(ymfxWeighted == null) || (ymfxWeighted.length != weights.length) ||
(vector == null) || (vector.length != IntersceneLmaFloat.NUM_PARAMS)) {
throw new IllegalArgumentException("GpuQuadJna.execPoseLmaStep: inconsistent inputs");
}
final float[] packed = new float[IntersceneLmaFloat.LMA_RESULT_FLOATS];
final int rc = lib.tp_proc_exec_pose_lma_step(
proc, lambda, weights, jt, ymfxWeighted, vector, weights.length, packed);
if (rc != 0) {
throw new IllegalStateException("GpuQuadJna.execPoseLmaStep rc="+rc+": "+lib.tp_last_error());
}
final int deltaOffset = 1 + IntersceneLmaFloat.LMA_PRODUCTS;
final int candidateOffset = deltaOffset + IntersceneLmaFloat.NUM_PARAMS;
return new IntersceneLmaFloat.LmaStepResult(
packed[0] != 0.0f,
Arrays.copyOfRange(packed, 1, deltaOffset),
Arrays.copyOfRange(packed, deltaOffset, candidateOffset),
Arrays.copyOfRange(packed, candidateOffset, IntersceneLmaFloat.LMA_RESULT_FLOATS));
}
// ---- LMA normal-equation products (tp_lma.cu, roadmap rung 3) ---- // By Codex on 07/14/2026
@Override public double[] execLmaNormalProducts(
double[] weights, double[][] jt, double[] ymfxWeighted) {
......
......@@ -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 = 29; // 19 base + 2 conditioning + 4 consolidate + 2 peak + pose fx/J + products
private static final int EXPECTED_KERNELS = 30; // 19 base + 2 conditioning + 4 consolidate + 2 peak + pose fx/J + pose step + products
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";
......
......@@ -12,8 +12,8 @@ public interface TpJna extends Library {
/** NVRTC-compile the kernels in srcdir (+ getTpDefines), cuLink(libcudadevrt), load the module.
* Returns an opaque module handle, or null on failure (see tp_last_error()). */
Pointer tp_create_module(String srcdir, String libcudadevrt);
/** Number of kernel functions resolved in the module (29 expected:
* 19 base + 2 conditioning + 4 consolidation + 2 peak + pose-fx/J + LMA-products),
/** Number of kernel functions resolved in the module (30 expected:
* 19 base + 2 conditioning + 4 consolidation + 2 peak + pose-fx/J + pose-step + LMA-products),
* or -1 if handle is null. */
int tp_module_num_functions(Pointer module);
/** Last error message (NVRTC/CUDA log), empty if none. */
......@@ -148,6 +148,12 @@ public interface TpJna extends Library {
float[] poseVectors, float[] centers, byte[] selection,
int numTiles, int calculateJacobian,
float[] fx, float[] jt, byte[] valid);
/** One-thread fixed three-angle float candidate calculation. Result layout is
* {valid, raw H[9], raw b[3], delta[3], candidate[3]}. */
int tp_proc_exec_pose_lma_step(Pointer proc,
float lambda, float[] weights, float[] jt,
float[] ymfxWeighted, float[] vector,
int numValues, float[] result);
/** 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. */
......
......@@ -100,6 +100,18 @@ public class IntersceneLma {
boolean calculateJacobian);
}
private PoseFxProvider poseFxProvider = null;
/** Optional one-thread float backend for the frozen three-angle LMA candidate. */
@FunctionalInterface
public interface PoseLmaStepProvider {
IntersceneLmaFloat.LmaStepResult calculate(
float lambda,
float [] weights,
float [] jt,
float [] ymfxWeighted,
float [] vector);
}
private PoseLmaStepProvider poseLmaStepProvider = null;
private static boolean poseLmaStepOracleReported = false;
private IntersceneLmaFloat.Camera poseFloatReference = null;
private IntersceneLmaFloat.Camera poseFloatScene = null;
private float [] poseFloatCenters = null;
......@@ -125,6 +137,10 @@ public class IntersceneLma {
public void setPoseFxProvider(final PoseFxProvider provider) {
this.poseFxProvider = provider;
}
/** Install the fixed three-angle float candidate backend. Null keeps the Java-double solve. */
public void setPoseLmaStepProvider(final PoseLmaStepProvider provider) {
this.poseLmaStepProvider = provider;
}
private void resetPoseFxInputs() {
poseFloatReference = null;
poseFloatScene = null;
......@@ -895,22 +911,23 @@ public class IntersceneLma {
}
Matrix wjtjlambda;
Matrix jty;
double [] normalProducts = null;
if (normalEquationProvider != null) {
final int numPars = this.last_jt.length;
final double [] products = normalEquationProvider.products(
normalProducts = normalEquationProvider.products(
this.weights, this.last_jt, this.last_ymfx);
if ((products == null) || (products.length != numPars * numPars + numPars)) {
if ((normalProducts == null) || (normalProducts.length != numPars * numPars + numPars)) {
throw new IllegalStateException("IntersceneLma: normal-equation provider returned "+
((products == null) ? "null" : products.length)+" values, expected "+
((normalProducts == null) ? "null" : normalProducts.length)+" values, expected "+
(numPars * numPars + numPars));
}
final double [][] h = new double [numPars][numPars];
for (int i = 0; i < numPars; i++) {
System.arraycopy(products, i * numPars, h[i], 0, numPars);
System.arraycopy(normalProducts, i * numPars, h[i], 0, numPars);
h[i][i] += h[i][i] * lambda; // same LMA convention as getWJtJlambda()
}
final double [] b = new double [numPars];
System.arraycopy(products, numPars * numPars, b, 0, numPars);
System.arraycopy(normalProducts, numPars * numPars, b, 0, numPars);
wjtjlambda = new Matrix(h);
jty = new Matrix(b, b.length);
} else {
......@@ -959,6 +976,44 @@ public class IntersceneLma {
for (int i = 0; i < parameters_vector.length; i++) {
new_vector[i] += scale * delta[i];
}
if ((poseLmaStepProvider != null) && isPoseLmaStepShape()) {
final float [] floatWeights = toFloatArray(weights);
final float [] floatJt = flattenJtFloat(last_jt);
final float [] floatYmfx = toFloatArray(last_ymfx);
final float [] floatVector = toFloatArray(parameters_vector);
final IntersceneLmaFloat.LmaStepResult poseStep = poseLmaStepProvider.calculate(
(float) lambda, floatWeights, floatJt, floatYmfx, floatVector);
if ((debug_level >= 1) && !poseLmaStepOracleReported && (poseStep != null)) {
poseLmaStepOracleReported = true;
final IntersceneLmaFloat.LmaStepResult javaFloatStep = IntersceneLmaFloat.lmaStep(
(float) lambda, floatWeights, floatJt, floatYmfx, floatVector);
final IntersceneLmaFloat.LmaStepComparison comparison =
IntersceneLmaFloat.compareLmaSteps(poseStep, javaFloatStep);
System.out.println("IntersceneLma CUDA vs Java-float runLma candidate: "+
comparison.format());
float maxProducts = 0.0f;
if ((normalProducts != null) && (normalProducts.length == poseStep.products.length)) {
for (int i = 0; i < normalProducts.length; i++) {
maxProducts = Math.max(maxProducts,
Math.abs(poseStep.products[i] - (float) normalProducts[i]));
}
}
float maxDelta = 0.0f;
float maxCandidate = 0.0f;
for (int i = 0; i < IntersceneLmaFloat.NUM_PARAMS; i++) {
maxDelta = Math.max(maxDelta, Math.abs(poseStep.delta[i] - (float) delta[i]));
maxCandidate = Math.max(maxCandidate,
Math.abs(poseStep.candidate[i] - (float) new_vector[i]));
}
System.out.println(String.format(
"IntersceneLma CUDA-float vs Java-double runLma candidate: valid=%s, max|H,b|=%g, max|delta|=%g, max|candidate|=%g",
poseStep.valid, maxProducts, maxDelta, maxCandidate));
}
if ((poseStep != null) && poseStep.valid &&
(poseStep.candidate.length == new_vector.length)) {
for (int i = 0; i < new_vector.length; i++) new_vector[i] = poseStep.candidate[i];
}
}
double [] fx = getFxDerivs(
new_vector, // double [] vector,
......@@ -1778,6 +1833,36 @@ public class IntersceneLma {
return true;
}
private boolean isPoseLmaStepShape() {
if ((num_components != 2) || (par_indices == null) ||
(par_indices.length != IntersceneLmaFloat.NUM_PARAMS) ||
(last_jt == null) || (last_jt.length != IntersceneLmaFloat.NUM_PARAMS)) return false;
for (int par = 0; par < IntersceneLmaFloat.NUM_PARAMS; par++) {
if (par_indices[par] != ErsCorrection.DP_DSAZ + par) return false;
}
return true;
}
private static float [] toFloatArray(final double [] source) {
final float [] result = new float [source.length];
for (int i = 0; i < source.length; i++) result[i] = (float) source[i];
return result;
}
private static float [] flattenJtFloat(final double [][] source) {
final int length = source[0].length;
final float [] result = new float [source.length * length];
for (int par = 0; par < source.length; par++) {
if ((source[par] == null) || (source[par].length != length)) {
throw new IllegalArgumentException("IntersceneLma: inconsistent Jacobian rows");
}
for (int value = 0; value < length; value++) {
result[par * length + value] = (float) source[par][value];
}
}
return result;
}
private double [][] getWJtJlambda( // USED in lwir
final double lambda,
final double [][] jt)
......
......@@ -40,6 +40,8 @@ import Jama.Matrix;
public final class IntersceneLmaFloat {
public static final int NUM_PARAMS = 3;
public static final int NUM_COMPONENTS = 2;
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 float INFINITY_DISPARITY = 0.01f;
static final int ERS_XYZ = 0;
......@@ -156,6 +158,136 @@ public final class IntersceneLmaFloat {
}
}
/** One fixed three-angle LMA candidate: raw H/b, delta and updated parameter vector. */
public static final class LmaStepResult {
public final boolean valid;
public final float [] products; // row-major H[3][3], then b[3], before damping
public final float [] delta;
public final float [] candidate;
public LmaStepResult(
final boolean valid,
final float [] products,
final float [] delta,
final float [] candidate) {
this.valid = valid;
this.products = products;
this.delta = delta;
this.candidate = candidate;
}
}
/** Statistics for checking the one-thread CUDA step against the Java float clone. */
public static final class LmaStepComparison {
public boolean validityMatch;
public boolean bitExact;
public float maxAbsProducts;
public float maxAbsDelta;
public float maxAbsCandidate;
public String format() {
return String.format(
"valid-match=%s, bit-exact=%s, max|H,b|=%g, max|delta|=%g, max|candidate|=%g",
validityMatch, bitExact, maxAbsProducts, maxAbsDelta, maxAbsCandidate);
}
}
/**
* Primitive-float clone of the fixed three-angle candidate calculation in
* {@code IntersceneLma.lmaStep()}. It deliberately keeps all reductions and
* the explicit 3x3 solve serial so the matching CUDA kernel can use one
* thread; parallel tile preparation/reduction is a later residency layer.
*/
public static LmaStepResult lmaStep(
final float lambda,
final float [] weights,
final float [] jt,
final float [] ymfxWeighted,
final float [] vector) {
if ((weights == null) || (weights.length == 0) ||
(jt == null) || (jt.length != NUM_PARAMS * weights.length) ||
(ymfxWeighted == null) || (ymfxWeighted.length != weights.length) ||
(vector == null) || (vector.length != NUM_PARAMS)) {
throw new IllegalArgumentException("inconsistent fixed pose LMA arrays");
}
final int numValues = weights.length;
final float [] products = new float [LMA_PRODUCTS];
for (int row = 0; row < NUM_PARAMS; row++) {
final int ri = row * numValues;
for (int column = row; column < NUM_PARAMS; column++) {
final int ci = column * numValues;
float sum = 0.0f;
for (int value = 0; value < numValues; value++) {
sum += (weights[value] * jt[ri + value]) * jt[ci + value];
}
products[row * NUM_PARAMS + column] = sum;
products[column * NUM_PARAMS + row] = sum;
}
float sum = 0.0f;
for (int value = 0; value < numValues; value++) {
sum += jt[ri + value] * ymfxWeighted[value];
}
products[NUM_PARAMS * NUM_PARAMS + row] = sum;
}
final float [] damped = Arrays.copyOf(products, NUM_PARAMS * NUM_PARAMS);
for (int i = 0; i < NUM_PARAMS; i++) {
final int diagonal = i * NUM_PARAMS + i;
damped[diagonal] += damped[diagonal] * lambda;
}
final float [] inverse = inverse3x3(damped);
final float [] delta = new float [] {Float.NaN, Float.NaN, Float.NaN};
final float [] candidate = vector.clone();
if (inverse == null) return new LmaStepResult(false, products, delta, candidate);
final int bi = NUM_PARAMS * NUM_PARAMS;
for (int row = 0; row < NUM_PARAMS; row++) {
delta[row] = inverse[row * NUM_PARAMS] * products[bi] +
inverse[row * NUM_PARAMS + 1] * products[bi + 1] +
inverse[row * NUM_PARAMS + 2] * products[bi + 2];
candidate[row] += delta[row];
}
return new LmaStepResult(allFinite(delta) && allFinite(candidate), products, delta, candidate);
}
public static LmaStepComparison compareLmaSteps(
final LmaStepResult actual,
final LmaStepResult expected) {
if ((actual == null) || (expected == null) ||
(actual.products == null) || (expected.products == null) ||
(actual.products.length != expected.products.length) ||
(actual.delta == null) || (expected.delta == null) ||
(actual.delta.length != expected.delta.length) ||
(actual.candidate == null) || (expected.candidate == null) ||
(actual.candidate.length != expected.candidate.length)) {
throw new IllegalArgumentException("incompatible float LMA step results");
}
final LmaStepComparison comparison = new LmaStepComparison();
comparison.validityMatch = actual.valid == expected.valid;
comparison.bitExact = comparison.validityMatch;
comparison.maxAbsProducts = compareFloatArrays(actual.products, expected.products, comparison);
comparison.maxAbsDelta = compareFloatArrays(actual.delta, expected.delta, comparison);
comparison.maxAbsCandidate = compareFloatArrays(actual.candidate, expected.candidate, comparison);
return comparison;
}
private static float compareFloatArrays(
final float [] actual,
final float [] expected,
final LmaStepComparison comparison) {
float maximum = 0.0f;
for (int i = 0; i < actual.length; i++) {
final boolean sameBits = Float.floatToRawIntBits(actual[i]) ==
Float.floatToRawIntBits(expected[i]);
if (!sameBits) comparison.bitExact = false;
if (finite(actual[i]) && finite(expected[i])) {
maximum = Math.max(maximum, Math.abs(actual[i] - expected[i]));
} else if (!sameBits) {
maximum = Float.POSITIVE_INFINITY;
}
}
return maximum;
}
/** Convert the existing sparse center representation to CUDA-style xyz triples. */
public static float [] flattenCenters(final double [][] centers) {
final float [] flat = new float [3 * centers.length];
......
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