Merge origin/main into parameter-support-v2 and add CLAUDE.md

Resolve merge conflicts in parameter.c and scalar_mult.c: keep
parameter-support-v2 features (parameter_expr, param_id, free_type_data)
while adopting _impl suffix naming from main. Fix old 2-arg
new_left_matmul calls in test_chain_rule_jacobian.h to use new 3-arg
signature with NULL param_node.

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

Author: Transurgeon

CLAUDE.md

@@ -0,0 +1,109 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Build & Test
+
+```bash
+# Build (from repo root)
+mkdir -p build && cd build && cmake .. -DCMAKE_BUILD_TYPE=Debug && make
+
+# Run all tests
+./build/all_tests
+
+# Run tests via CTest
+cd build && ctest --output-on-failure
+
+# Format check (requires clang-format)
+find src include -name '*.c' -o -name '*.h' | xargs clang-format --dry-run -Werror
+
+# Build with sanitizers (for debugging)
+cd build && cmake .. -DCMAKE_BUILD_TYPE=Debug -DCMAKE_C_FLAGS="-fsanitize=address,undefined" && make
+```
+
+There is no way to run a single test in isolation — `all_tests` runs the full suite. To test selectively, temporarily comment out `mu_run_test()` calls in `tests/all_tests.c`.
+
+## Architecture
+
+### Expr Node (Virtual Dispatch via Function Pointers)
+
+Every operation is an `expr` node (`include/expr.h`). Polymorphism is achieved through function pointers set at construction time via `init_expr()`:
+
+- `forward` — evaluate node value from input vector `u`
+- `jacobian_init` / `eval_jacobian` — allocate then fill sparse Jacobian (CSR)
+- `wsum_hess_init` / `eval_wsum_hess` — allocate then fill weighted-sum Hessian (CSR)
+- `is_affine` — used to skip zero Hessian computation
+- `local_jacobian` / `local_wsum_hess` — element-wise ops only: diagonal f'(g(x)) and f''(g(x))
+- `free_type_data` — cleanup for type-specific allocations (e.g., param_source pointers)
+
+Type-specific data uses C struct inheritance: cast the first field (`expr base`) to a subtype defined in `include/subexpr.h` (e.g., `parameter_expr`, `scalar_mult_expr`, `left_matmul_expr`).
+
+### Operation Categories
+
+| Directory | Description | Chain rule pattern |
+|-----------|-------------|-------------------|
+| `src/affine/` | Linear ops (variable, parameter, add, index, matmul, etc.) | J = structural combination of child Jacobians |
+| `src/elementwise_full_dom/` | Univariate ops (exp, sin, log, power, etc.) | J = diag(f'(g(x))) · J_child; common chain rule in `common.c` |
+| `src/elementwise_restricted_dom/` | Univariate ops with domain restrictions (log, entropy) | Same pattern as full_dom |
+| `src/bivariate_full_dom/` | Two-arg ops (elementwise multiply) | Product rule on two children |
+| `src/bivariate_restricted_dom/` | Two-arg restricted domain (quad_over_lin, rel_entr) | Custom chain rules |
+| `src/other/` | Special ops (prod, quad_form) | Custom implementations |
+
+### Init vs Eval Split
+
+Derivative computation is split into two phases:
+1. **Init** (`jacobian_init`, `wsum_hess_init`) — allocates CSR matrices and sets sparsity patterns. Called once per expression tree. Guarded by `jacobian_init()` / `wsum_hess_init()` wrappers in `expr.c` to handle DAGs where nodes are visited multiple times.
+2. **Eval** (`eval_jacobian`, `eval_wsum_hess`) — fills numerical values into pre-allocated structures. Called each time input changes.
+
+### Elementwise Common Chain Rule
+
+`src/elementwise_full_dom/common.c` and `src/elementwise_restricted_dom/common.c` implement shared chain rule logic for all univariate ops. Each op only provides `local_jacobian` (diagonal of f') and `local_wsum_hess` (diagonal of f''). The common code handles:
+- Variable child: Jacobian is diagonal
+- Composite child: J = diag(f') · J_child; H = J_child^T · diag(w·f'') · J_child + f' · H_child
+
+### Parameter System
+
+Parameters (`src/affine/parameter.c`) unify constants and updatable values:
+- `param_id == PARAM_FIXED` (-1): fixed constant, value set at construction
+- `param_id >= 0`: updatable parameter, indexed into a global theta array
+
+Problem-level API (`src/problem.c`):
+- `problem_register_params()` — register parameter nodes
+- `problem_update_params()` — update all parameter values from theta vector
+
+Operations that use parameters (scalar_mult, vector_mult, left/right_matmul) store a `param_source` pointer to the parameter node and read its value during forward/jacobian/hessian evaluation.
+
+### Workspace (`Expr_Work`)
+
+Each node has a `work` pointer for cached intermediate results: CSC conversion of Jacobian, local Jacobian diagonal, Hessian term1/term2 matrices. Allocated during init, reused across evals.
+
+## Adding a New Operation
+
+1. If needed, add a subtype struct in `include/subexpr.h`
+2. Declare the constructor in the appropriate header (`include/affine.h`, `include/elementwise_full_dom.h`, etc.)
+3. Implement in `src/<category>/<op_name>.c` with static functions for forward, jacobian_init_impl, eval_jacobian, wsum_hess_init_impl, eval_wsum_hess, is_affine
+4. For elementwise univariate ops: only implement `local_jacobian` and `local_wsum_hess`, then use `common_jacobian_init`/`common_eval_jacobian`/etc. from `common.c`
+5. Add tests in `tests/forward_pass/`, `tests/jacobian_tests/`, `tests/wsum_hess/` and register them in `tests/all_tests.c`
+6. Use `check_jacobian()` and `check_wsum_hess()` from `tests/numerical_diff.h` to validate against finite differences
+
+## Code Style
+
+C99. Enforced by `.clang-format`: Allman braces, 4-space indent, 85-column limit, right-aligned pointers. Run `clang-format -i <file>` before committing.
+
+## CI Workflows (`.github/workflows/`)
+
+- `cmake.yml` — build and test on Linux, macOS, Windows
+- `formatting.yml` — clang-format check
+- `sanitizer.yml` — ASan, MSan, UBSan
+- `valgrind.yml` — memory leak detection
+- `release.yml` — release packaging
+
+## Platform Notes
+
+- **macOS**: Links Accelerate for BLAS
+- **Linux**: Links OpenBLAS (`libopenblas-dev`)
+- **Windows**: Links OpenBLAS via vcpkg; no dense matmul in some configurations
+
+## Test Tolerances
+
+`ABS_TOL = 1e-6`, `REL_TOL = 1e-6` (defined in test headers). Numerical differentiation step size: `NUMERICAL_DIFF_DEFAULT_H = 1e-7`.

include/expr.h

@@ -73,8 +73,8 @@ typedef struct expr
     CSR_Matrix *jacobian;
     CSR_Matrix *wsum_hess;
     forward_fn forward;
-    jacobian_init_fn jacobian_init;
-    wsum_hess_init_fn wsum_hess_init;
+    jacobian_init_fn jacobian_init_impl;
+    wsum_hess_init_fn wsum_hess_init_impl;
     eval_jacobian_fn eval_jacobian;
     wsum_hess_fn eval_wsum_hess;
 
@@ -99,6 +99,11 @@ void init_expr(expr *node, int d1, int d2, int n_vars, forward_fn forward,
 
 void free_expr(expr *node);
 
+/* Guarded init: skips if already initialized (safe for DAGs
+ * where a node may be visited through multiple parents). */
+void jacobian_init(expr *node);
+void wsum_hess_init(expr *node);
+
 /* Initialize CSC form of the Jacobian from the CSR Jacobian.
  * Must be called after jacobian_init. */
 void jacobian_csc_init(expr *node);

src/affine/add.c

@@ -34,11 +34,11 @@ static void forward(expr *node, const double *u)
     }
 }
 
-static void jacobian_init(expr *node)
+static void jacobian_init_impl(expr *node)
 {
     /* initialize children's jacobians */
-    node->left->jacobian_init(node->left);
-    node->right->jacobian_init(node->right);
+    jacobian_init(node->left);
+    jacobian_init(node->right);
 
     /* we never have to store more than the sum of children's nnz */
     int nnz_max = node->left->jacobian->nnz + node->right->jacobian->nnz;
@@ -60,11 +60,11 @@ static void eval_jacobian(expr *node)
                                  node->jacobian);
 }
 
-static void wsum_hess_init(expr *node)
+static void wsum_hess_init_impl(expr *node)
 {
     /* initialize children's wsum_hess */
-    node->left->wsum_hess_init(node->left);
-    node->right->wsum_hess_init(node->right);
+    wsum_hess_init(node->left);
+    wsum_hess_init(node->right);
 
     /* we never have to store more than the sum of children's nnz */
     int nnz_max = node->left->wsum_hess->nnz + node->right->wsum_hess->nnz;
@@ -95,8 +95,8 @@ expr *new_add(expr *left, expr *right)
 {
     assert(left->d1 == right->d1 && left->d2 == right->d2);
     expr *node = (expr *) calloc(1, sizeof(expr));
-    init_expr(node, left->d1, left->d2, left->n_vars, forward, jacobian_init,
-              eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess, NULL);
+    init_expr(node, left->d1, left->d2, left->n_vars, forward, jacobian_init_impl,
+              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
     node->left = left;
     node->right = right;
     expr_retain(left);

src/affine/broadcast.c

@@ -66,10 +66,10 @@ static void forward(expr *node, const double *u)
     }
 }
 
-static void jacobian_init(expr *node)
+static void jacobian_init_impl(expr *node)
 {
     expr *x = node->left;
-    x->jacobian_init(x);
+    jacobian_init(x);
     broadcast_expr *bcast = (broadcast_expr *) node;
     int total_nnz;
 
@@ -185,10 +185,10 @@ static void eval_jacobian(expr *node)
     }
 }
 
-static void wsum_hess_init(expr *node)
+static void wsum_hess_init_impl(expr *node)
 {
     expr *x = node->left;
-    x->wsum_hess_init(x);
+    wsum_hess_init(x);
 
     /* Same sparsity as child - weights get summed */
     node->wsum_hess = new_csr_copy_sparsity(x->wsum_hess);
@@ -279,8 +279,8 @@ expr *new_broadcast(expr *child, int d1, int d2)
     // --------------------------------------------------------------------------
     //                  initialize the rest of the expression
     // --------------------------------------------------------------------------
-    init_expr(node, d1, d2, child->n_vars, forward, jacobian_init, eval_jacobian,
-              is_affine, wsum_hess_init, eval_wsum_hess, NULL);
+    init_expr(node, d1, d2, child->n_vars, forward, jacobian_init_impl,
+              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
     node->left = child;
     expr_retain(child);
     bcast->type = type;

src/affine/diag_vec.c

@@ -44,11 +44,11 @@ static void forward(expr *node, const double *u)
     }
 }
 
-static void jacobian_init(expr *node)
+static void jacobian_init_impl(expr *node)
 {
     expr *x = node->left;
     int n = x->size;
-    x->jacobian_init(x);
+    jacobian_init(x);
 
     CSR_Matrix *Jx = x->jacobian;
     CSR_Matrix *J = new_csr_matrix(node->size, node->n_vars, Jx->nnz);
@@ -92,12 +92,12 @@ static void eval_jacobian(expr *node)
     }
 }
 
-static void wsum_hess_init(expr *node)
+static void wsum_hess_init_impl(expr *node)
 {
     expr *x = node->left;
 
     /* initialize child's wsum_hess */
-    x->wsum_hess_init(x);
+    wsum_hess_init(x);
 
     /* workspace for extracting diagonal weights */
     node->work->dwork = (double *) calloc(x->size, sizeof(double));
@@ -137,8 +137,8 @@ expr *new_diag_vec(expr *child)
     /* n is the number of elements (works for both row and column vectors) */
     int n = child->size;
     expr *node = (expr *) calloc(1, sizeof(expr));
-    init_expr(node, n, n, child->n_vars, forward, jacobian_init, eval_jacobian,
-              is_affine, wsum_hess_init, eval_wsum_hess, NULL);
+    init_expr(node, n, n, child->n_vars, forward, jacobian_init_impl, eval_jacobian,
+              is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
     node->left = child;
     expr_retain(child);
 

src/affine/hstack.c

@@ -42,7 +42,7 @@ static void forward(expr *node, const double *u)
     }
 }
 
-static void jacobian_init(expr *node)
+static void jacobian_init_impl(expr *node)
 {
     hstack_expr *hnode = (hstack_expr *) node;
 
@@ -51,7 +51,7 @@ static void jacobian_init(expr *node)
     for (int i = 0; i < hnode->n_args; i++)
     {
         assert(hnode->args[i] != NULL);
-        hnode->args[i]->jacobian_init(hnode->args[i]);
+        jacobian_init(hnode->args[i]);
         nnz += hnode->args[i]->jacobian->nnz;
     }
 
@@ -100,14 +100,14 @@ static void eval_jacobian(expr *node)
     }
 }
 
-static void wsum_hess_init(expr *node)
+static void wsum_hess_init_impl(expr *node)
 {
     /* initialize children's hessians */
     hstack_expr *hnode = (hstack_expr *) node;
     int nnz = 0;
     for (int i = 0; i < hnode->n_args; i++)
     {
-        hnode->args[i]->wsum_hess_init(hnode->args[i]);
+        wsum_hess_init(hnode->args[i]);
         nnz += hnode->args[i]->wsum_hess->nnz;
     }
 
@@ -187,8 +187,9 @@ expr *new_hstack(expr **args, int n_args, int n_vars)
     /* Allocate the type-specific struct */
     hstack_expr *hnode = (hstack_expr *) calloc(1, sizeof(hstack_expr));
     expr *node = &hnode->base;
-    init_expr(node, args[0]->d1, d2, n_vars, forward, jacobian_init, eval_jacobian,
-              is_affine, wsum_hess_init, wsum_hess_eval, free_type_data);
+    init_expr(node, args[0]->d1, d2, n_vars, forward, jacobian_init_impl,
+              eval_jacobian, is_affine, wsum_hess_init_impl, wsum_hess_eval,
+              free_type_data);
 
     /* Set type-specific fields (deep copy args array) */
     hnode->args = (expr **) calloc(n_args, sizeof(expr *));

src/affine/index.c

@@ -57,11 +57,11 @@ static void forward(expr *node, const double *u)
     }
 }
 
-static void jacobian_init(expr *node)
+static void jacobian_init_impl(expr *node)
 {
     expr *x = node->left;
     index_expr *idx = (index_expr *) node;
-    x->jacobian_init(x);
+    jacobian_init(x);
 
     CSR_Matrix *Jx = x->jacobian;
     CSR_Matrix *J = new_csr_matrix(node->size, node->n_vars, Jx->nnz);
@@ -96,12 +96,12 @@ static void eval_jacobian(expr *node)
     }
 }
 
-static void wsum_hess_init(expr *node)
+static void wsum_hess_init_impl(expr *node)
 {
     expr *x = node->left;
 
     /* initialize child's wsum_hess */
-    x->wsum_hess_init(x);
+    wsum_hess_init(x);
 
     /* for setting weight vector to evaluate hessian of child */
     node->work->dwork = (double *) calloc(x->size, sizeof(double));
@@ -166,8 +166,9 @@ expr *new_index(expr *child, int d1, int d2, const int *indices, int n_idxs)
     index_expr *idx = (index_expr *) calloc(1, sizeof(index_expr));
     expr *node = &idx->base;
 
-    init_expr(node, d1, d2, child->n_vars, forward, jacobian_init, eval_jacobian,
-              is_affine, wsum_hess_init, eval_wsum_hess, free_type_data);
+    init_expr(node, d1, d2, child->n_vars, forward, jacobian_init_impl,
+              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
+              free_type_data);
 
     node->left = child;
     expr_retain(child);

src/affine/left_matmul.c

@@ -104,13 +104,13 @@ static void free_type_data(expr *node)
     lnode->param_source = NULL;
 }
 
-static void jacobian_init(expr *node)
+static void jacobian_init_impl(expr *node)
 {
     expr *x = node->left;
     left_matmul_expr *lnode = (left_matmul_expr *) node;
 
     /* initialize child's jacobian and precompute sparsity of its CSC */
-    x->jacobian_init(x);
+    jacobian_init(x);
     lnode->Jchild_CSC = csr_to_csc_fill_sparsity(x->jacobian, node->work->iwork);
 
     /* precompute sparsity of this node's jacobian in CSC and CSR */
@@ -136,11 +136,11 @@ static void eval_jacobian(expr *node)
     csc_to_csr_fill_values(J_CSC, node->jacobian, lnode->csc_to_csr_work);
 }
 
-static void wsum_hess_init(expr *node)
+static void wsum_hess_init_impl(expr *node)
 {
     /* initialize child's hessian */
     expr *x = node->left;
-    x->wsum_hess_init(x);
+    wsum_hess_init(x);
 
     /* allocate this node's hessian with the same sparsity as child's */
     node->wsum_hess = new_csr_copy_sparsity(x->wsum_hess);
@@ -220,8 +220,8 @@ expr *new_left_matmul(expr *param_node, expr *u, const CSR_Matrix *A)
     left_matmul_expr *lnode =
         (left_matmul_expr *) calloc(1, sizeof(left_matmul_expr));
     expr *node = &lnode->base;
-    init_expr(node, d1, d2, u->n_vars, forward, jacobian_init, eval_jacobian,
-              is_affine, wsum_hess_init, eval_wsum_hess, free_type_data);
+    init_expr(node, d1, d2, u->n_vars, forward, jacobian_init_impl, eval_jacobian,
+              is_affine, wsum_hess_init_impl, eval_wsum_hess, free_type_data);
     node->left = u;
     expr_retain(u);
 
@@ -274,8 +274,8 @@ expr *new_left_matmul_dense(expr *param_node, expr *u, int m, int n,
     left_matmul_expr *lnode =
         (left_matmul_expr *) calloc(1, sizeof(left_matmul_expr));
     expr *node = &lnode->base;
-    init_expr(node, d1, d2, u->n_vars, forward, jacobian_init, eval_jacobian,
-              is_affine, wsum_hess_init, eval_wsum_hess, free_type_data);
+    init_expr(node, d1, d2, u->n_vars, forward, jacobian_init_impl, eval_jacobian,
+              is_affine, wsum_hess_init_impl, eval_wsum_hess, free_type_data);
     node->left = u;
     expr_retain(u);