Upgrade upper_tri jacobian and hessian tests to 4x4

The 3x3 tests produced identical indices for both row-major and
column-major orderings, so they didn't actually verify the fix.
4x4 indices [4,8,12,9,13,14] differ from column-major [4,8,9,12,13,14].

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: Transurgeon

tests/jacobian_tests/affine/test_upper_tri.h

@@ -8,55 +8,64 @@
 
 const char *test_upper_tri_jacobian_variable(void)
 {
-    /* upper_tri of a 3x3 variable (9 vars total)
-     * Upper tri flat indices: [3, 6, 7]
-     * Jacobian is 3x9 CSR: row k has col indices[k] */
-    double u[9] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
-    expr *var = new_variable(3, 3, 0, 9);
+    /* upper_tri of a 4x4 variable (16 vars total)
+     * Row-major upper tri indices: [4, 8, 12, 9, 13, 14]
+     * Jacobian is 6x16 CSR: row k has a single 1.0 at col indices[k] */
+    double u[16];
+    for (int k = 0; k < 16; k++)
+    {
+        u[k] = (double) (k + 1);
+    }
+    expr *var = new_variable(4, 4, 0, 16);
     expr *ut = new_upper_tri(var);
 
     ut->forward(ut, u);
     jacobian_init(ut);
     ut->eval_jacobian(ut);
 
-    double expected_x[3] = {1.0, 1.0, 1.0};
-    int expected_p[4] = {0, 1, 2, 3};
-    int expected_i[3] = {3, 6, 7};
+    double expected_x[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
+    int expected_p[7] = {0, 1, 2, 3, 4, 5, 6};
+    int expected_i[6] = {4, 8, 12, 9, 13, 14};
 
     mu_assert("upper_tri jac vals",
-              cmp_double_array(ut->jacobian->x, expected_x, 3));
-    mu_assert("upper_tri jac p", cmp_int_array(ut->jacobian->p, expected_p, 4));
-    mu_assert("upper_tri jac i", cmp_int_array(ut->jacobian->i, expected_i, 3));
+              cmp_double_array(ut->jacobian->x, expected_x, 6));
+    mu_assert("upper_tri jac p",
+              cmp_int_array(ut->jacobian->p, expected_p, 7));
+    mu_assert("upper_tri jac i",
+              cmp_int_array(ut->jacobian->i, expected_i, 6));
 
     free_expr(ut);
     return 0;
 }
 
 const char *test_upper_tri_jacobian_of_log(void)
 {
-    /* upper_tri(log(X)) where X is 3x3 variable
-     * Upper tri flat indices: [3, 6, 7]
-     * X values at those positions: x[3]=4, x[6]=7, x[7]=8
-     * d/dx log at those positions:
-     * Row 0: 1/4 = 0.25 at col 3
-     * Row 1: 1/7 at col 6
-     * Row 2: 1/8 = 0.125 at col 7 */
-    double u[9] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
-    expr *var = new_variable(3, 3, 0, 9);
+    /* upper_tri(log(X)) where X is 4x4 variable
+     * Row-major upper tri indices: [4, 8, 12, 9, 13, 14]
+     * Values at those positions: u[4]=5, u[8]=9, u[12]=13,
+     *                            u[9]=10, u[13]=14, u[14]=15
+     * d/dx log at those positions: 1/5, 1/9, 1/13, 1/10, 1/14, 1/15 */
+    double u[16];
+    for (int k = 0; k < 16; k++)
+    {
+        u[k] = (double) (k + 1);
+    }
+    expr *var = new_variable(4, 4, 0, 16);
     expr *log_node = new_log(var);
     expr *ut = new_upper_tri(log_node);
 
     ut->forward(ut, u);
     jacobian_init(ut);
     ut->eval_jacobian(ut);
 
-    double expected_x[3] = {0.25, 1.0 / 7.0, 0.125};
-    int expected_i[3] = {3, 6, 7};
+    double expected_x[6] = {0.2, 1.0 / 9.0, 1.0 / 13.0,
+                            0.1, 1.0 / 14.0, 1.0 / 15.0};
+    int expected_i[6] = {4, 8, 12, 9, 13, 14};
 
     mu_assert("upper_tri log jac vals",
-              cmp_double_array(ut->jacobian->x, expected_x, 3));
+              cmp_double_array(ut->jacobian->x, expected_x, 6));
     mu_assert("upper_tri log jac cols",
-              cmp_int_array(ut->jacobian->i, expected_i, 3));
+              cmp_int_array(ut->jacobian->i, expected_i, 6));
 
     free_expr(ut);
     return 0;

tests/wsum_hess/affine/test_upper_tri.h

@@ -11,24 +11,29 @@
 
 const char *test_wsum_hess_upper_tri_log(void)
 {
-    /* upper_tri(log(X)) where X is 3x3, w = [1, 1, 1]
-     * X (column-major): [1, 2, 3, 4, 5, 6, 7, 8, 9]
-     * Upper tri flat indices (row-major): [3, 6, 7]
-     * Hessian of log is diag(-1/x^2)
-     * Weights scatter: parent_w[3]=1, parent_w[6]=1, parent_w[7]=1
-     * All other parent_w entries = 0
+    /* upper_tri(log(X)) where X is 4x4, w = [1, 1, 1, 1, 1, 1]
+     * X (column-major): [1..16]
+     * Row-major upper tri indices: [4, 8, 12, 9, 13, 14]
+     * Weights scatter to parent_w: 1.0 at positions
+     * {4, 8, 9, 12, 13, 14}, zero elsewhere.
      *
+     * Hessian of log is diag(-1/x^2).
      * H[k,k] = parent_w[k] * (-1/x[k]^2):
-     * H[0,0] = 0, H[1,1] = 0, H[2,2] = 0
-     * H[3,3] = -1/16 = -0.0625
-     * H[4,4] = 0, H[5,5] = 0
-     * H[6,6] = -1/49
-     * H[7,7] = -1/64 = -0.015625
-     * H[8,8] = 0 */
-    double u[9] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
-    double w[3] = {1.0, 1.0, 1.0};
+     *   H[4,4]   = -1/25  = -0.04
+     *   H[8,8]   = -1/81
+     *   H[9,9]   = -1/100 = -0.01
+     *   H[12,12] = -1/169
+     *   H[13,13] = -1/196
+     *   H[14,14] = -1/225
+     *   All others = 0 */
+    double u[16];
+    for (int k = 0; k < 16; k++)
+    {
+        u[k] = (double) (k + 1);
+    }
+    double w[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
 
-    expr *var = new_variable(3, 3, 0, 9);
+    expr *var = new_variable(4, 4, 0, 16);
     expr *log_node = new_log(var);
     expr *ut = new_upper_tri(log_node);
 
@@ -37,17 +42,22 @@ const char *test_wsum_hess_upper_tri_log(void)
     wsum_hess_init(ut);
     ut->eval_wsum_hess(ut, w);
 
-    double expected_x[9] = {0.0, 0.0,         0.0,       -0.0625, 0.0,
-                            0.0, -1.0 / 49.0, -0.015625, 0.0};
-    int expected_p[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
-    int expected_i[9] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+    double expected_x[16] = {
+        0.0,          0.0,          0.0,          0.0,
+        -1.0 / 25.0,  0.0,          0.0,          0.0,
+        -1.0 / 81.0,  -1.0 / 100.0, 0.0,          0.0,
+        -1.0 / 169.0, -1.0 / 196.0, -1.0 / 225.0, 0.0};
+    int expected_p[17] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
+                          11, 12, 13, 14, 15, 16};
+    int expected_i[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
+                          11, 12, 13, 14, 15};
 
     mu_assert("upper_tri log hess vals",
-              cmp_double_array(ut->wsum_hess->x, expected_x, 9));
+              cmp_double_array(ut->wsum_hess->x, expected_x, 16));
     mu_assert("upper_tri log hess p",
-              cmp_int_array(ut->wsum_hess->p, expected_p, 10));
+              cmp_int_array(ut->wsum_hess->p, expected_p, 17));
     mu_assert("upper_tri log hess i",
-              cmp_int_array(ut->wsum_hess->i, expected_i, 9));
+              cmp_int_array(ut->wsum_hess->i, expected_i, 16));
 
     free_expr(ut);
     return 0;