Indexing validation

[baltzell]

I was looking for optimizations, in this case some smaller methods for potential inlining in #1258 in RungeKuttaDoca.computeCovMat, and was curious about the indexing here, particularly theC.get(2,0/1):

coatjava/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/RungeKuttaDoca.java

Lines 238 to 254 in 0e88d70

	cNext[3][0] = C.get(2,0)*sF.dtydtx0+C.get(3,0)*sF.dtydty0+C.get(4,0)*sF.dtydq0
	+ sF.dxdq0 *(C.get(2,4)*sF.dtydtx0+C.get(3,4)*sF.dtydty0+C.get(4,4)*sF.dtydq0)
	+ sF.dxdtx0*(C.get(2,2)*sF.dtydtx0+C.get(3,2)*sF.dtydty0+C.get(4,2)*sF.dtydq0)
	+ sF.dxdty0*(C.get(2,3)*sF.dtydtx0+C.get(3,3)*sF.dtydty0+C.get(4,3)*sF.dtydq0);
	cNext[3][1] = C.get(2,1)*sF.dtydtx0+C.get(3,1)*sF.dtydty0+C.get(4,1)*sF.dtydq0
	+ sF.dydq0 *(C.get(2,4)*sF.dtydtx0+C.get(3,4)*sF.dtydty0+C.get(4,4)*sF.dtydq0)
	+ sF.dydtx0*(C.get(2,2)*sF.dtydtx0+C.get(3,2)*sF.dtydty0+C.get(4,2)*sF.dtydq0)
	+ sF.dydty0*(C.get(2,3)*sF.dtydtx0+C.get(3,3)*sF.dtydty0+C.get(4,3)*sF.dtydq0);
	cNext[3][2] = sF.dtxdq0 *(C.get(2,4)*sF.dtydtx0+C.get(3,4)*sF.dtydty0+C.get(4,4)*sF.dtydq0)
	+ sF.dtxdtx0*(C.get(2,2)*sF.dtydtx0+C.get(3,2)*sF.dtydty0+C.get(4,2)*sF.dtydq0)
	+ sF.dtxdty0*(C.get(2,3)*sF.dtydtx0+C.get(3,3)*sF.dtydty0+C.get(4,3)*sF.dtydq0);
	cNext[3][3] = sF.dtydq0 *(C.get(2,4)*sF.dtydtx0+C.get(3,4)*sF.dtydty0+C.get(4,4)*sF.dtydq0)
	+ sF.dtydtx0*(C.get(2,2)*sF.dtydtx0+C.get(3,2)*sF.dtydty0+C.get(4,2)*sF.dtydq0)
	+ sF.dtydty0*(C.get(2,3)*sF.dtydtx0+C.get(3,3)*sF.dtydty0+C.get(4,3)*sF.dtydq0);

Is this snippet really correct relative to the previous 3 indices? Those do C.get(i,0/1).

[zieglerv]

The covariance matrix calculation is cNext = F * C * F^T and the state order is {x, y, tx, ty, q}.
The jacobian rows are:
F[0] = {1, 0, dxdtx0, dxdty0, dxdq0}
F[1] = {0, 1, dydtx0, dydty0, dydq0}
F[2] = {0, 0, dtxdtx0, dtxdty0, dtxdq0}
F[3] = {0, 0, dtydtx0, dtydty0, dtydq0}
F[4] = {0, 0, 0, 0, 1}

But in the proposed code change the doit method is used:
doit(0, cNext, C, sF);
doit(1, cNext, C, sF);
doit2(2, cNext, C, sF);
doit2(3, cNext, C, sF);

But doit(1, ...) still uses the x-row derivatives:
sF.dxdtx0
sF.dxdty0
sF.dxdq0

When it should use:
sF.dydtx0
sF.dydty0
sF.dydq0

Similarly, doit2(3, ...) still uses the tx-row derivatives:
sF.dtxdtx0
sF.dtxdty0
sF.dtxdq0

When it should use:
sF.dtydtx0
sF.dtydty0
sF.dtydq0

So these are incorrect:
cNext[1][0..4]
cNext[3][0..4]

Regarding the indexing question: The C.get(2,0/1) indexing is correct for the rows corresponding to the transported tx derivative row of the Jacobian.
For the transported tx row, the Jacobian row is F_tx = {0, 0, dtxdtx0, dtxdty0, dtxdq0}.
So when computing entries like Cnext[2][0] or Cnext[2][1], the left side of the multiplication starts from row/column index 2 of C, not from index 0 or 1, because tx has no direct derivative with respect to the initial x0 or y0.
That's why the expression begins with C.get(2,0) * sF.dtxdtx0 + C.get(3,0) * sF.dtxdty0 + C.get(4,0) * sF.dtxdq0 and similarly, for the first column with C.get(2,1) * sF.dtxdtx0 + C.get(3,1) * sF.dtxdty0 + C.get(4,1) * sF.dtxdq0.

As written doit() / doit2() are problematic because the output row index and the Jacobian row are not the same thing. An alternate version that makes each transported row explicit and preserves the correct indexing is:

/** Compute the final covariance matrix. covMat = F C F^T. */
private double[][] computeCovMat(Matrix C, RK4Vec sF) {
double[][] cNext = new double[5][5];

final double c00 = C.get(0,0);
final double c01 = C.get(0,1);
final double c02 = C.get(0,2);
final double c03 = C.get(0,3);
final double c04 = C.get(0,4);

final double c10 = C.get(1,0);
final double c11 = C.get(1,1);
final double c12 = C.get(1,2);
final double c13 = C.get(1,3);
final double c14 = C.get(1,4);

final double c20 = C.get(2,0);
final double c21 = C.get(2,1);
final double c22 = C.get(2,2);
final double c23 = C.get(2,3);
final double c24 = C.get(2,4);

final double c30 = C.get(3,0);
final double c31 = C.get(3,1);
final double c32 = C.get(3,2);
final double c33 = C.get(3,3);
final double c34 = C.get(3,4);

final double c40 = C.get(4,0);
final double c41 = C.get(4,1);
final double c42 = C.get(4,2);
final double c43 = C.get(4,3);
final double c44 = C.get(4,4);

/*
 * State order:
 *
 *   0 = x
 *   1 = y
 *   2 = tx
 *   3 = ty
 *   4 = q
 *
 * Transport Jacobian rows:
 *
 *   Fx  = {1, 0, dxdtx0,  dxdty0,  dxdq0}
 *   Fy  = {0, 1, dydtx0,  dydty0,  dydq0}
 *   Ftx = {0, 0, dtxdtx0, dtxdty0, dtxdq0}
 *   Fty = {0, 0, dtydtx0, dtydty0, dtydq0}
 *   Fq  = {0, 0, 0,       0,       1}
 */

final double x_tx  = sF.dxdtx0;
final double x_ty  = sF.dxdty0;
final double x_q   = sF.dxdq0;

final double y_tx  = sF.dydtx0;
final double y_ty  = sF.dydty0;
final double y_q   = sF.dydq0;

final double tx_tx = sF.dtxdtx0;
final double tx_ty = sF.dtxdty0;
final double tx_q  = sF.dtxdq0;

final double ty_tx = sF.dtydtx0;
final double ty_ty = sF.dtydty0;
final double ty_q  = sF.dtydq0;

/*
 * Compute C * F_j^T for each transported row.
 *
 * The tx and ty rows start at covariance index 2 and 3 because
 * those Jacobian rows have no x0 or y0 component.
 */

// C * Fx^T
final double cx0 = c00 + c02*x_tx + c03*x_ty + c04*x_q;
final double cx1 = c10 + c12*x_tx + c13*x_ty + c14*x_q;
final double cx2 = c20 + c22*x_tx + c23*x_ty + c24*x_q;
final double cx3 = c30 + c32*x_tx + c33*x_ty + c34*x_q;
final double cx4 = c40 + c42*x_tx + c43*x_ty + c44*x_q;

// C * Fy^T
final double cy0 = c01 + c02*y_tx + c03*y_ty + c04*y_q;
final double cy1 = c11 + c12*y_tx + c13*y_ty + c14*y_q;
final double cy2 = c21 + c22*y_tx + c23*y_ty + c24*y_q;
final double cy3 = c31 + c32*y_tx + c33*y_ty + c34*y_q;
final double cy4 = c41 + c42*y_tx + c43*y_ty + c44*y_q;

// C * Ftx^T
final double ctx2 = c22*tx_tx + c23*tx_ty + c24*tx_q;
final double ctx3 = c32*tx_tx + c33*tx_ty + c34*tx_q;
final double ctx4 = c42*tx_tx + c43*tx_ty + c44*tx_q;

// C * Fty^T
final double cty2 = c22*ty_tx + c23*ty_ty + c24*ty_q;
final double cty3 = c32*ty_tx + c33*ty_ty + c34*ty_q;
final double cty4 = c42*ty_tx + c43*ty_ty + c44*ty_q;

/*
 * Upper triangle of F C F^T.
 *
 * The lower triangle is filled by symmetry.
 */

final double o00 = cx0 + x_tx*cx2 + x_ty*cx3 + x_q*cx4;
final double o01 = cy0 + x_tx*cy2 + x_ty*cy3 + x_q*cy4;
final double o02 =       x_tx*ctx2 + x_ty*ctx3 + x_q*ctx4;
final double o03 =       x_tx*cty2 + x_ty*cty3 + x_q*cty4;
final double o04 = c04 + x_tx*c24 + x_ty*c34 + x_q*c44;

final double o11 = cy1 + y_tx*cy2 + y_ty*cy3 + y_q*cy4;
final double o12 =       y_tx*ctx2 + y_ty*ctx3 + y_q*ctx4;
final double o13 =       y_tx*cty2 + y_ty*cty3 + y_q*cty4;
final double o14 = c14 + y_tx*c24 + y_ty*c34 + y_q*c44;

final double o22 = tx_tx*ctx2 + tx_ty*ctx3 + tx_q*ctx4;
final double o23 = tx_tx*cty2 + tx_ty*cty3 + tx_q*cty4;
final double o24 = tx_tx*c24 + tx_ty*c34 + tx_q*c44;

final double o33 = ty_tx*cty2 + ty_ty*cty3 + ty_q*cty4;
final double o34 = ty_tx*c24 + ty_ty*c34 + ty_q*c44;

final double o44 = c44;

cNext[0][0] = o00;
cNext[0][1] = o01;
cNext[0][2] = o02;
cNext[0][3] = o03;
cNext[0][4] = o04;

cNext[1][0] = o01;
cNext[1][1] = o11;
cNext[1][2] = o12;
cNext[1][3] = o13;
cNext[1][4] = o14;

cNext[2][0] = o02;
cNext[2][1] = o12;
cNext[2][2] = o22;
cNext[2][3] = o23;
cNext[2][4] = o24;

cNext[3][0] = o03;
cNext[3][1] = o13;
cNext[3][2] = o23;
cNext[3][3] = o33;
cNext[3][4] = o34;

cNext[4][0] = o04;
cNext[4][1] = o14;
cNext[4][2] = o24;
cNext[4][3] = o34;
cNext[4][4] = o44;

return cNext;

}

private void covStep(RK4Vec next, double qv) {
double Csq = Csq(next.x, next.y);
double C = C(Csq);
double Ax = Ax(C, next.x, next.y);
double Ay = Ay(C, next.x, next.y);

next.tx = qv * Ax;
next.ty = qv * Ay;

double dAx_dtx = dAx_dtx(C, Csq, Ax, next.x, next.y);
double dAx_dty = dAx_dty(C, Csq, Ax, next.x, next.y);
double dAy_dtx = dAy_dtx(C, Csq, Ay, next.x, next.y);
double dAy_dty = dAy_dty(C, Csq, Ay, next.x, next.y);

next.dtxdtx0 = dtx_dtx0(qv, dAx_dtx, dAx_dty,
                        next.dxdtx0, next.dydtx0);

next.dtxdty0 = dtx_dty0(qv, dAx_dtx, dAx_dty,
                        next.dxdty0, next.dydty0);

next.dtxdq0  = dtx_dq0(qv, Ax, dAx_dtx, dAx_dty,
                       next.dxdq0, next.dydq0);

next.dtydtx0 = dty_dtx0(qv, dAy_dtx, dAy_dty,
                        next.dxdtx0, next.dydtx0);

next.dtydty0 = dty_dty0(qv, dAy_dtx, dAy_dty,
                        next.dxdty0, next.dydty0);

next.dtydq0  = dty_dq0(qv, Ay, dAy_dtx, dAy_dty,
                       next.dxdq0, next.dydq0);

}

Also: private static double dtx_dty0(double qv,
double dAx_dtx, double dAx_dty,
double dxdty0, double dydty0) {
return qv * (dAx_dtydtx_dty0 + dAx_dtydty_dty0);
}
Should be:
private static double dtx_dty0(double qv,
double dAx_dtx, double dAx_dty,
double dxdty0, double dydty0) {
return qv * (dAx_dtxdxdty0 + dAx_dtydydty0);
}

Same bug fix:
private static double dty_dty0(double qv,
double dAy_dtx, double dAy_dty,
double dxdty0, double dydty0) {
return qv * (dAy_dtxdxdty0 + dAy_dtydydty0);
}

[baltzell]

What proposed code? I'm referring to the main branch, see link in above.