Numerical derivative checker (#56)

* wsum checker

* formatter

Author: dance858

tests/all_tests.c

@@ -300,6 +300,7 @@ int main(void)
 
     printf("\n--- Numerical Diff Tests ---\n");
     mu_run_test(test_check_jacobian_composite_log, tests_run);
+    mu_run_test(test_check_wsum_hess_log_composite, tests_run);
 
     printf("\n--- Problem Struct Tests ---\n");
     mu_run_test(test_problem_new_free, tests_run);

tests/numerical_diff.c

@@ -97,3 +97,98 @@ int check_jacobian(expr *node, const double *u, double h)
     free(J_num);
     return result;
 }
+
+/* Compute g = J^T w where J is CSR (m x n) and w has m entries.
+ * Result written into g (size n), which must be zero-initialized. */
+static void csr_transpose_mult_vec(const CSR_Matrix *J, const double *w, double *g)
+{
+    for (int row = 0; row < J->m; row++)
+    {
+        for (int idx = J->p[row]; idx < J->p[row + 1]; idx++)
+        {
+            g[J->i[idx]] += J->x[idx] * w[row];
+        }
+    }
+}
+
+double *numerical_wsum_hess(expr *node, const double *u, const double *w, double h)
+{
+    int n = node->n_vars;
+    double inv_2h = 1.0 / (2.0 * h);
+
+    /* Initialize jacobian sparsity once, then forward */
+    node->jacobian_init(node);
+    node->forward(node, u);
+
+    double *H = calloc((size_t) n * n, sizeof(double));
+    double *u_work = malloc(n * sizeof(double));
+    double *g_plus = malloc(n * sizeof(double));
+    double *g_minus = malloc(n * sizeof(double));
+
+    memcpy(u_work, u, n * sizeof(double));
+
+    for (int j = 0; j < n; j++)
+    {
+        /* g(u + h*e_j) */
+        u_work[j] = u[j] + h;
+        node->forward(node, u_work);
+        node->eval_jacobian(node);
+        memset(g_plus, 0, n * sizeof(double));
+        csr_transpose_mult_vec(node->jacobian, w, g_plus);
+
+        /* g(u - h*e_j) */
+        u_work[j] = u[j] - h;
+        node->forward(node, u_work);
+        node->eval_jacobian(node);
+        memset(g_minus, 0, n * sizeof(double));
+        csr_transpose_mult_vec(node->jacobian, w, g_minus);
+
+        u_work[j] = u[j];
+
+        for (int i = 0; i < n; i++)
+        {
+            H[i * n + j] = (g_plus[i] - g_minus[i]) * inv_2h;
+        }
+    }
+
+    free(g_minus);
+    free(g_plus);
+    free(u_work);
+    return H;
+}
+
+int check_wsum_hess(expr *node, const double *u, const double *w, double h)
+{
+    int n = node->n_vars;
+
+    /* Compute numerical first (does its own jacobian_init) */
+    double *H_num = numerical_wsum_hess(node, u, w, h);
+
+    /* Now compute analytical (reuses jacobian from numerical) */
+    node->wsum_hess_init(node);
+    node->forward(node, u);
+    node->eval_jacobian(node);
+    node->eval_wsum_hess(node, w);
+
+    double *H_ana = calloc((size_t) n * n, sizeof(double));
+    csr_to_dense(node->wsum_hess, H_ana);
+
+    int result = 1;
+    for (int i = 0; i < n * n; i++)
+    {
+        if (!is_close(H_ana[i], H_num[i]))
+        {
+            int row = i / n;
+            int col = i % n;
+            printf("  check_wsum_hess FAILED at (%d, %d):"
+                   " analytical=%e, numerical=%e\n",
+                   row, col, H_ana[i], H_num[i]);
+            result = 0;
+            break;
+        }
+    }
+
+    free(H_ana);
+    free(H_num);
+    return result;
+}

tests/numerical_diff.h

@@ -15,4 +15,15 @@ double *numerical_jacobian(expr *node, const double *u, double h);
  * Prints diagnostic on first failing entry. */
 int check_jacobian(expr *node, const double *u, double h);
 
+/* Compute dense numerical weighted-sum Hessian via central
+ * differences on the gradient g(u) = J(u)^T w.
+ * Returns malloc'd row-major array (n_vars x n_vars).
+ * Caller must free(). */
+double *numerical_wsum_hess(expr *node, const double *u, const double *w, double h);
+
+/* Evaluate analytical wsum_hess, compute numerical wsum_hess,
+ * and compare. Returns 1 on match, 0 on mismatch.
+ * Prints diagnostic on first failing entry. */
+int check_wsum_hess(expr *node, const double *u, const double *w, double h);
+
 #endif /* NUMERICAL_DIFF_H */

tests/numerical_diff/test_numerical_diff.h

@@ -28,3 +28,27 @@ const char *test_check_jacobian_composite_log(void)
     free_csr_matrix(A);
     return 0;
 }
+
+const char *test_check_wsum_hess_log_composite(void)
+{
+    double u_vals[5] = {1, 2, 3, 4, 5};
+    double w[3] = {-1, -2, -3};
+    double Ax[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
+    int Ai[] = {0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4};
+    int Ap[] = {0, 5, 10, 15};
+    CSR_Matrix *A_csr = new_csr_matrix(3, 5, 15);
+    memcpy(A_csr->x, Ax, 15 * sizeof(double));
+    memcpy(A_csr->i, Ai, 15 * sizeof(int));
+    memcpy(A_csr->p, Ap, 4 * sizeof(int));
+
+    expr *x = new_variable(5, 1, 0, 5);
+    expr *Ax_node = new_linear(x, A_csr, NULL);
+    expr *log_node = new_log(Ax_node);
+
+    mu_assert("check_wsum_hess failed",
+              check_wsum_hess(log_node, u_vals, w, NUMERICAL_DIFF_DEFAULT_H));
+
+    free_expr(log_node);
+    free_csr_matrix(A_csr);
+    return 0;
+}