Unify Constant and Parameter into single parameter type

Merge the separate Constant and Parameter leaf nodes into a unified
parameter_expr with PARAM_FIXED sentinel (-1) for constants. This
eliminates duplicate code paths and consolidates 7 bivariate constructors
into 3 unified ones:

- new_const_scalar_mult / new_param_scalar_mult -> new_scalar_mult
- new_const_vector_mult / new_param_vector_mult -> new_vector_mult
- new_left_matmul (CSR) / new_left_param_matmul -> new_left_matmul (param node)

Key changes:
- Add PARAM_FIXED define and extend new_parameter() to accept initial values
- Delete constant.c (absorbed by parameter.c)
- Remove direct value storage (double a, double *a) from scalar/vector mult
  structs; always read from param_source
- left_matmul builds sparse CSR for fixed params (preserving sparsity) and
  dense CSR for updatable params
- right_matmul internally creates a fixed parameter node from transposed A
- problem_update_params skips PARAM_FIXED nodes
- Update all test callers to use new_parameter with PARAM_FIXED

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: Transurgeon

include/affine.h

@@ -32,9 +32,8 @@ expr *new_hstack(expr **args, int n_args, int n_vars);
 expr *new_promote(expr *child, int d1, int d2);
 expr *new_trace(expr *child);
 
-expr *new_constant(int d1, int d2, int n_vars, const double *values);
 expr *new_variable(int d1, int d2, int var_id, int n_vars);
-expr *new_parameter(int d1, int d2, int param_id, int n_vars);
+expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *values);
 
 expr *new_index(expr *child, int d1, int d2, const int *indices, int n_idxs);
 expr *new_reshape(expr *child, int d1, int d2);

include/bivariate.h

@@ -30,25 +30,16 @@ expr *new_rel_entr_second_arg_scalar(expr *left, expr *right);
 /* Matrix multiplication: Z = X @ Y */
 expr *new_matmul(expr *x, expr *y);
 
-/* Left matrix multiplication: A @ f(x) where A is a constant matrix */
-expr *new_left_matmul(expr *u, const CSR_Matrix *A);
+/* Left matrix multiplication: A @ f(x) where A comes from a parameter node */
+expr *new_left_matmul(expr *param_node, expr *child);
 
 /* Right matrix multiplication: f(x) @ A where A is a constant matrix */
 expr *new_right_matmul(expr *u, const CSR_Matrix *A);
 
-/* Constant scalar multiplication: a * f(x) where a is a constant double */
-expr *new_const_scalar_mult(double a, expr *child);
+/* Scalar multiplication: a * f(x) where a comes from a parameter node */
+expr *new_scalar_mult(expr *param_node, expr *child);
 
-/* Constant vector elementwise multiplication: a ∘ f(x) where a is constant */
-expr *new_const_vector_mult(const double *a, expr *child);
-
-/* Left matrix multiplication: P @ f(x) where P is a parameter */
-expr *new_left_param_matmul(expr *param_node, expr *child);
-
-/* Parameter scalar multiplication: p * f(x) where p is a parameter */
-expr *new_param_scalar_mult(expr *param_node, expr *child);
-
-/* Parameter vector elementwise multiplication: p ∘ f(x) where p is a parameter */
-expr *new_param_vector_mult(expr *param_node, expr *child);
+/* Vector elementwise multiplication: a ∘ f(x) where a comes from a parameter node */
+expr *new_vector_mult(expr *param_node, expr *child);
 
 #endif /* BIVARIATE_H */

include/subexpr.h

@@ -25,12 +25,17 @@
 /* Forward declaration */
 struct int_double_pair;
 
-/* Parameter node: like constant but with updatable values via problem_update_params
+/* param_id value for fixed (constant) parameters */
+#define PARAM_FIXED -1
+
+/* Parameter node: unified leaf for constants and updatable parameters.
+ * Constants use param_id == PARAM_FIXED and have values set at creation.
+ * Updatable parameters have param_id >= 0 and are updated via problem_update_params.
  */
 typedef struct parameter_expr
 {
     expr base;
-    int param_id; /* offset into global theta vector */
+    int param_id; /* offset into global theta vector, or PARAM_FIXED */
 } parameter_expr;
 
 /* Type-specific expression structures that "inherit" from expr */
@@ -136,20 +141,18 @@ typedef struct right_matmul_expr
     CSC_Matrix *CSC_work;
 } right_matmul_expr;
 
-/* Constant scalar multiplication: y = a * child where a is a constant double */
+/* Scalar multiplication: y = a * child where a comes from a parameter node */
 typedef struct const_scalar_mult_expr
 {
     expr base;
-    double a;
-    expr *param_source; /* if non-NULL, read a from param_source->value[0] */
+    expr *param_source; /* always set; read a from param_source->value[0] */
 } const_scalar_mult_expr;
 
-/* Constant vector elementwise multiplication: y = a \circ child for constant a */
+/* Vector elementwise multiplication: y = a \circ child where a comes from a parameter node */
 typedef struct const_vector_mult_expr
 {
     expr base;
-    double *a;          /* length equals node->size */
-    expr *param_source; /* if non-NULL, use param_source->value instead of a */
+    expr *param_source; /* always set; read a from param_source->value */
 } const_vector_mult_expr;
 
 /* Index/slicing: y = child[indices] where indices is a list of flat positions */

src/affine/constant.c

@@ -16,10 +16,14 @@
  * limitations under the License.
  */
 
-/* Parameter leaf node: behaviorally identical to constant (zero derivatives
-   w.r.t. variables), but its values are updatable via problem_update_params.
-   This allows re-solving with different parameter values without rebuilding
-   the expression tree. */
+/* Unified parameter/constant leaf node.
+ *
+ * When param_id == PARAM_FIXED, this is a constant whose values are set at
+ * creation and never change. When param_id >= 0, values are updated via
+ * problem_update_params.
+ *
+ * In both cases the derivative behavior is identical: zero Jacobian and
+ * Hessian with respect to variables (always affine). */
 
 #include "affine.h"
 #include "subexpr.h"
@@ -28,26 +32,26 @@
 
 static void forward(expr *node, const double *u)
 {
-    /* Values are set by problem_update_params, not by forward pass */
+    /* Values are set at creation (constants) or by problem_update_params */
     (void) node;
     (void) u;
 }
 
 static void jacobian_init(expr *node)
 {
-    /* Parameter jacobian is all zeros: size x n_vars with 0 nonzeros */
+    /* Parameter/constant jacobian is all zeros: size x n_vars with 0 nonzeros */
     node->jacobian = new_csr_matrix(node->size, node->n_vars, 0);
 }
 
 static void eval_jacobian(expr *node)
 {
-    /* Parameter jacobian never changes */
+    /* Jacobian never changes */
     (void) node;
 }
 
 static void wsum_hess_init(expr *node)
 {
-    /* Parameter Hessian is all zeros */
+    /* Hessian is all zeros */
     node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, 0);
 }
 
@@ -63,12 +67,19 @@ static bool is_affine(const expr *node)
     return true;
 }
 
-expr *new_parameter(int d1, int d2, int param_id, int n_vars)
+expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *values)
 {
     parameter_expr *pnode = (parameter_expr *) calloc(1, sizeof(parameter_expr));
     init_expr(&pnode->base, d1, d2, n_vars, forward, jacobian_init, eval_jacobian,
               is_affine, wsum_hess_init, eval_wsum_hess, NULL);
     pnode->param_id = param_id;
-    /* values will be populated by problem_update_params */
+
+    /* If values provided (fixed constant), copy them now */
+    if (values != NULL)
+    {
+        memcpy(pnode->base.value, values, pnode->base.size * sizeof(double));
+    }
+    /* Otherwise values will be populated by problem_update_params */
+
     return &pnode->base;
 }

src/affine/parameter.c

@@ -22,7 +22,7 @@
 #include <stdlib.h>
 #include <string.h>
 
-/* Constant scalar multiplication: y = a * child where a is a constant double */
+/* Scalar multiplication: y = a * child where a comes from a parameter node */
 
 static void forward(expr *node, const double *u)
 {
@@ -33,7 +33,7 @@ static void forward(expr *node, const double *u)
 
     /* local forward pass: multiply each element by scalar a */
     const_scalar_mult_expr *sn = (const_scalar_mult_expr *) node;
-    double a = sn->param_source ? sn->param_source->value[0] : sn->a;
+    double a = sn->param_source->value[0];
     for (int i = 0; i < node->size; i++)
     {
         node->value[i] = a * child->value[i];
@@ -57,7 +57,7 @@ static void eval_jacobian(expr *node)
 {
     expr *child = node->left;
     const_scalar_mult_expr *sn = (const_scalar_mult_expr *) node;
-    double a = sn->param_source ? sn->param_source->value[0] : sn->a;
+    double a = sn->param_source->value[0];
 
     /* evaluate child */
     child->eval_jacobian(child);
@@ -88,7 +88,7 @@ static void eval_wsum_hess(expr *node, const double *w)
     x->eval_wsum_hess(x, w);
 
     const_scalar_mult_expr *sn = (const_scalar_mult_expr *) node;
-    double a = sn->param_source ? sn->param_source->value[0] : sn->a;
+    double a = sn->param_source->value[0];
     for (int j = 0; j < x->wsum_hess->nnz; j++)
     {
         node->wsum_hess->x[j] = a * x->wsum_hess->x[j];
@@ -101,23 +101,7 @@ static bool is_affine(const expr *node)
     return node->left->is_affine(node->left);
 }
 
-expr *new_const_scalar_mult(double a, expr *child)
-{
-    const_scalar_mult_expr *mult_node =
-        (const_scalar_mult_expr *) calloc(1, sizeof(const_scalar_mult_expr));
-    expr *node = &mult_node->base;
-
-    init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init,
-              eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess, NULL);
-    node->left = child;
-    mult_node->a = a;
-    mult_node->param_source = NULL;
-    expr_retain(child);
-
-    return node;
-}
-
-static void free_param_type_data(expr *node)
+static void free_type_data(expr *node)
 {
     const_scalar_mult_expr *sn = (const_scalar_mult_expr *) node;
     if (sn->param_source)
@@ -126,17 +110,16 @@ static void free_param_type_data(expr *node)
     }
 }
 
-expr *new_param_scalar_mult(expr *param_node, expr *child)
+expr *new_scalar_mult(expr *param_node, expr *child)
 {
     const_scalar_mult_expr *mult_node =
         (const_scalar_mult_expr *) calloc(1, sizeof(const_scalar_mult_expr));
     expr *node = &mult_node->base;
 
     init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init,
               eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess,
-              free_param_type_data);
+              free_type_data);
     node->left = child;
-    mult_node->a = param_node->value[0]; /* initial value */
     mult_node->param_source = param_node;
     expr_retain(child);
     expr_retain(param_node);