[WIP] Unify constants and parameters into single parameter_expr type (#53)

* Unify constants and parameters into single parameter_expr type

Constants and updatable parameters are now represented by a single
parameter_expr node. Constants use param_id == PARAM_FIXED (-1),
while updatable parameters use param_id >= 0. Bivariate ops
(left_matmul, right_matmul, scalar_mult, vector_mult) accept an
optional param_source node and refresh their internal data on
forward/jacobian/hessian evaluation. Adds problem_register_params
and problem_update_params for problem-level parameter management.
Also adds update_values to the Matrix interface for parameter refresh.

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

* Address review comments on parameter-support-v2

- Remove redundant refresh_param_values calls from eval_jacobian and
  eval_wsum_hess in left_matmul (forward always runs first)
- Use memcpy in problem_register_params for pointer array copy
- Add PARAM_FIXED guard in problem_update_params to skip fixed constants
- Remove unused right_matmul_expr struct from subexpr.h
- Add test_param_fixed_skip_in_update covering mixed fixed/updatable params
- Add CLAUDE.md for Claude Code guidance

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

* Run clang-format and remove CLAUDE.md

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

* Set has_been_refreshed to false in parameter constructor

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

* Remove CLAUDE.md from PR

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

* remove unused constructors

* run formatter and push it up

* Update includes in test_param_prob.h to use atoms/ prefix

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

* Fix column-major parameter ordering in parameterized matmul (#72)

* Fix column-major parameter ordering in parameterized matmul

CVXPY sends parameter values in Fortran (column-major) order, but the
matmul refresh functions assumed row-major/CSR order via raw memcpy.
This produced incorrect matrix values for non-symmetric matrices.

For sparse matrices, iterate the CSR pattern and index into the
column-major source array. For dense matrices, exploit the fact that
column-major A is row-major A^T to memcpy directly into AT, then
transpose to get A.

Also fixes a latent bug where sparse update_values would blindly copy
the first nnz values from the full d1*d2 parameter array, which is
wrong for matrices with structural zeros.

Adds tests for rectangular (3x2) and sparse (3x3 with zeros) cases.

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

* Run clang-format on changed files

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

* introduce explicit transpose function for dense matrix

* clean up refresh dense right

* clean up tests...

* one more test

---------

Co-authored-by: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Co-authored-by: dance858 <danielcederberg1@gmail.com>

* add test that fails with current appraoch

* new parameter abstraction that fixess the fialing test

* remove dead code

* reset csc

* add raise error to catch bugs

* add initial attempt for fixing parameters and broadcasting (#73)

* add initial attempt for fixing parameters and broadcasting

* add test for params with broadcast

* cleanup test to fit more with other tests style

* some progress on supporting backwards compatible constants

* add some parameter broadcast tests as well

* cleanup left matmul as well

* some very minor cleanups

* some error checks and numerical diff to tests

* we don't always have to run forward of parameter in left matmul

* we don't always have to call forward for parameter node in vector mult

* comment out forward parameter pass in scalar mult because it is not needed, I think

* add test for scalar case to be consistent with vector mult

---------

Co-authored-by: dance858 <danielcederberg1@gmail.com>

* error message and format

---------

Co-authored-by: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Co-authored-by: dance858 <danielcederberg1@gmail.com>

Author: 3 people

include/atoms/affine.h

@@ -31,7 +31,7 @@ expr *new_vstack(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_parameter(int d1, int d2, int param_id, int n_vars, const double *values);
 expr *new_variable(int d1, int d2, int var_id, int n_vars);
 
 expr *new_index(expr *child, int d1, int d2, const int *indices, int n_idxs);
@@ -40,26 +40,27 @@ expr *new_broadcast(expr *child, int target_d1, int target_d2);
 expr *new_diag_vec(expr *child);
 expr *new_transpose(expr *child);
 
-/* Left matrix multiplication: A @ f(x) where A is a constant sparse
- * matrix */
-expr *new_left_matmul(expr *u, const CSR_Matrix *A);
+/* Left matrix multiplication: A @ f(x) where A is a constant or parameter
+ * sparse matrix. param_node is NULL for fixed constants. */
+expr *new_left_matmul(expr *param_node, expr *u, const CSR_Matrix *A);
 
-/* Left matrix multiplication: A @ f(x) where A is a constant dense
- * matrix (row-major, m x n). Uses CBLAS for efficient computation. */
-expr *new_left_matmul_dense(expr *u, int m, int n, const double *data);
+/* Left matrix multiplication: A @ f(x) where A is a constant or parameter
+ * dense matrix (row-major, m x n). Uses CBLAS for efficient computation. */
+expr *new_left_matmul_dense(expr *param_node, expr *u, int m, int n,
+                            const double *data);
 
-/* Right matrix multiplication: f(x) @ A where A is a constant
- * matrix */
-expr *new_right_matmul(expr *u, const CSR_Matrix *A);
+/* Right matrix multiplication: f(x) @ A where A is a constant or parameter
+ * matrix. */
+expr *new_right_matmul(expr *param_node, expr *u, const CSR_Matrix *A);
 
-expr *new_right_matmul_dense(expr *u, int m, int n, const double *data);
+expr *new_right_matmul_dense(expr *param_node, expr *u, int m, int n,
+                             const double *data);
 
-/* 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 param_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);
+/* Vector elementwise multiplication: a . f(x) where a comes from
+ * param_node */
+expr *new_vector_mult(expr *param_node, expr *child);
 
 #endif /* AFFINE_H */

include/expr.h

@@ -86,6 +86,12 @@ typedef struct expr
     local_wsum_hess_fn local_wsum_hess; /* used by elementwise univariate atoms*/
     free_type_data_fn free_type_data;   /* Cleanup for type-specific fields */
     Expr_Work *work;                    /* derivative workspace */
+    /* Set to true on all nodes by problem_update_params() via
+       expr_set_needs_refresh(). Atoms that cache parameter data
+       (e.g. left_matmul_dense) check this flag before their forward
+       pass: if true, they refresh their cached matrices from
+       param_source->value and clear the flag to false. */
+    bool needs_parameter_refresh;
 
     // name of node just for debugging - should be removed later
     char name[32];
@@ -108,6 +114,9 @@ void wsum_hess_init(expr *node);
  * Must be called after jacobian_init. */
 void jacobian_csc_init(expr *node);
 
+/* Recursively set needs_parameter_refresh on node and all children */
+void expr_set_needs_refresh(expr *node);
+
 /* Reference counting helpers */
 void expr_retain(expr *node);
 

include/problem.h

@@ -49,6 +49,11 @@ typedef struct problem
     int n_vars;
     int total_constraint_size;
 
+    /* parameter support */
+    expr **param_nodes;
+    int n_param_nodes;
+    int total_parameter_size;
+
     /* allocated by new_problem */
     double *constraint_values;
     double *gradient_values;
@@ -79,6 +84,9 @@ void problem_init_jacobian_coo(problem *prob);
 void problem_init_hessian_coo_lower_triangular(problem *prob);
 void free_problem(problem *prob);
 
+void problem_register_params(problem *prob, expr **param_nodes, int n_param_nodes);
+void problem_update_params(problem *prob, const double *theta);
+
 double problem_objective_forward(problem *prob, const double *u);
 void problem_constraint_forward(problem *prob, const double *u);
 void problem_gradient(problem *prob);

include/subexpr.h

@@ -26,8 +26,19 @@
 /* Forward declaration */
 struct int_double_pair;
 
+/* Parameter ID for fixed constants (not updatable) */
+#define PARAM_FIXED -1
+
 /* Type-specific expression structures that "inherit" from expr */
 
+/* Unified constant/parameter node. Constants use param_id == PARAM_FIXED.
+ * Updatable parameters use param_id >= 0 (offset into global theta). */
+typedef struct parameter_expr
+{
+    expr base;
+    int param_id;
+} parameter_expr;
+
 /* Linear operator: y = A * x + b
  * The matrix A is stored as node->jacobian (CSR). */
 typedef struct linear_op_expr
@@ -118,18 +129,24 @@ typedef struct left_matmul_expr
     CSC_Matrix *Jchild_CSC;
     CSC_Matrix *J_CSC;
     int *csc_to_csr_work;
+    expr *param_source;
+    void (*refresh_param_values)(struct left_matmul_expr *);
 } left_matmul_expr;
 
-/* Right matrix multiplication: y = f(x) * A where f(x) is an expression.
- * f(x) has shape p x n, A has shape n x q, output y has shape p x q.
- * Uses vec(y) = B * vec(f(x)) where B = A^T kron I_p. */
-typedef struct right_matmul_expr
+/* Scalar multiplication: y = a * child where a comes from param_source */
+typedef struct scalar_mult_expr
+{
+    expr base;
+    expr *param_source;
+} scalar_mult_expr;
+
+/* Vector elementwise multiplication: y = a \circ child where a comes from
+ * param_source */
+typedef struct vector_mult_expr
 {
     expr base;
-    CSR_Matrix *B;  /* B = A^T kron I_p */
-    CSR_Matrix *BT; /* B^T for backpropagating Hessian weights */
-    CSC_Matrix *CSC_work;
-} right_matmul_expr;
+    expr *param_source;
+} vector_mult_expr;
 
 /* Bivariate matrix multiplication: Z = f(u) @ g(u) where both children
  * may be composite expressions. */
@@ -153,20 +170,6 @@ typedef struct matmul_expr
     int *idx_map_Hg;
 } matmul_expr;
 
-/* Constant scalar multiplication: y = a * child where a is a constant double */
-typedef struct const_scalar_mult_expr
-{
-    expr base;
-    double a;
-} const_scalar_mult_expr;
-
-/* Constant vector elementwise multiplication: y = a \circ child for constant a */
-typedef struct const_vector_mult_expr
-{
-    expr base;
-    double *a; /* length equals node->size */
-} const_vector_mult_expr;
-
 /* Index/slicing: y = child[indices] where indices is a list of flat positions */
 typedef struct index_expr
 {

include/utils/dense_matrix.h

@@ -17,4 +17,6 @@ Matrix *new_dense_matrix(int m, int n, const double *data);
 /* Transpose helper */
 Matrix *dense_matrix_trans(const Dense_Matrix *self);
 
+void A_transpose(double *AT, const double *A, int m, int n);
+
 #endif /* DENSE_MATRIX_H */

include/utils/matrix.h

@@ -31,6 +31,7 @@ typedef struct Matrix
                                             const CSC_Matrix *J, int p);
     void (*block_left_mult_values)(const struct Matrix *self, const CSC_Matrix *J,
                                    CSC_Matrix *C);
+    void (*update_values)(struct Matrix *self, const double *new_values);
     void (*free_fn)(struct Matrix *self);
 } Matrix;
 

src/atoms/affine/left_matmul.c

@@ -49,16 +49,37 @@
 #include "utils/tracked_alloc.h"
 #include "utils/utils.h"
 
+static void refresh_param_values(left_matmul_expr *lnode)
+{
+    if (lnode->param_source == NULL || !lnode->base.needs_parameter_refresh)
+    {
+        return;
+    }
+
+    lnode->base.needs_parameter_refresh = false;
+    lnode->refresh_param_values(lnode);
+}
+
 static void forward(expr *node, const double *u)
 {
+    left_matmul_expr *lnode = (left_matmul_expr *) node;
+
+    /* call forward on param_source if it exists and needs refresh */
+    if (lnode->param_source != NULL && lnode->base.needs_parameter_refresh)
+    {
+        lnode->param_source->forward(lnode->param_source, NULL);
+    }
+
+    refresh_param_values(lnode);
+
     expr *x = node->left;
 
     /* child's forward pass */
     node->left->forward(node->left, u);
 
     /* y = A_kron @ vec(f(x)) */
-    Matrix *A = ((left_matmul_expr *) node)->A;
-    int n_blocks = ((left_matmul_expr *) node)->n_blocks;
+    Matrix *A = lnode->A;
+    int n_blocks = lnode->n_blocks;
     A->block_left_mult_vec(A, x->value, node->value, n_blocks);
 }
 
@@ -75,11 +96,16 @@ static void free_type_data(expr *node)
     free_csc_matrix(lnode->Jchild_CSC);
     free_csc_matrix(lnode->J_CSC);
     free(lnode->csc_to_csr_work);
+    if (lnode->param_source != NULL)
+    {
+        free_expr(lnode->param_source);
+    }
     lnode->A = NULL;
     lnode->AT = NULL;
     lnode->Jchild_CSC = NULL;
     lnode->J_CSC = NULL;
     lnode->csc_to_csr_work = NULL;
+    lnode->param_source = NULL;
 }
 
 static void jacobian_init_impl(expr *node)
@@ -99,8 +125,8 @@ static void jacobian_init_impl(expr *node)
 
 static void eval_jacobian(expr *node)
 {
-    expr *x = node->left;
     left_matmul_expr *lnode = (left_matmul_expr *) node;
+    expr *x = node->left;
 
     CSC_Matrix *Jchild_CSC = lnode->Jchild_CSC;
     CSC_Matrix *J_CSC = lnode->J_CSC;
@@ -131,17 +157,33 @@ static void wsum_hess_init_impl(expr *node)
 
 static void eval_wsum_hess(expr *node, const double *w)
 {
+    left_matmul_expr *lnode = (left_matmul_expr *) node;
+
     /* compute A^T w*/
-    Matrix *AT = ((left_matmul_expr *) node)->AT;
-    int n_blocks = ((left_matmul_expr *) node)->n_blocks;
+    Matrix *AT = lnode->AT;
+    int n_blocks = lnode->n_blocks;
     AT->block_left_mult_vec(AT, w, node->work->dwork, n_blocks);
 
     node->left->eval_wsum_hess(node->left, node->work->dwork);
     memcpy(node->wsum_hess->x, node->left->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
 }
 
-expr *new_left_matmul(expr *u, const CSR_Matrix *A)
+static void refresh_dense_left(left_matmul_expr *lnode)
+{
+    Dense_Matrix *dm_A = (Dense_Matrix *) lnode->A;
+    Dense_Matrix *dm_AT = (Dense_Matrix *) lnode->AT;
+    int m = dm_A->base.m;
+    int n = dm_A->base.n;
+
+    /* The parameter represents the A in left_matmul_dense(A, x) in column-major.
+       In this diffengine, we store A in row-major order. Hence, param->vals
+       actually corresponds to the transpose of A, and we transpose AT to get A. */
+    memcpy(dm_AT->x, lnode->param_source->value, m * n * sizeof(double));
+    A_transpose(dm_A->x, dm_AT->x, n, m);
+}
+
+expr *new_left_matmul(expr *param_node, expr *u, const CSR_Matrix *A)
 {
     /* We expect u->d1 == A->n. However, numpy's broadcasting rules allow users
        to do A @ u where u is (n, ) which in C is actually (1, n). In that case
@@ -188,10 +230,20 @@ expr *new_left_matmul(expr *u, const CSR_Matrix *A)
     lnode->AT =
         sparse_matrix_trans((const Sparse_Matrix *) lnode->A, node->work->iwork);
 
+    /* parameter support */
+    lnode->param_source = param_node;
+    if (param_node != NULL)
+    {
+        fprintf(stderr, "Error in new_left_matmul: parameter for a sparse matrix "
+                        "not supported \n");
+        exit(1);
+    }
+
     return node;
 }
 
-expr *new_left_matmul_dense(expr *u, int m, int n, const double *data)
+expr *new_left_matmul_dense(expr *param_node, expr *u, int m, int n,
+                            const double *data)
 {
     int d1, d2, n_blocks;
     if (u->d1 == n)
@@ -227,5 +279,13 @@ expr *new_left_matmul_dense(expr *u, int m, int n, const double *data)
     lnode->A = new_dense_matrix(m, n, data);
     lnode->AT = dense_matrix_trans((const Dense_Matrix *) lnode->A);
 
+    /* parameter support */
+    lnode->param_source = param_node;
+    if (param_node != NULL)
+    {
+        expr_retain(param_node);
+        lnode->refresh_param_values = refresh_dense_left;
+    }
+
     return node;
 }

src/atoms/affine/constant.c src/atoms/affine/parameter.csrc/atoms/affine/constant.c renamed to src/atoms/affine/parameter.c

@@ -16,39 +16,38 @@
  * limitations under the License.
  */
 #include "atoms/affine.h"
+#include "subexpr.h"
 #include "utils/tracked_alloc.h"
+#include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 
 static void forward(expr *node, const double *u)
 {
-    /* Constants don't depend on u; values are already set */
+    /* Parameters/constants don't depend on u; values are already set */
     (void) node;
     (void) u;
 }
 
 static void jacobian_init_impl(expr *node)
 {
-    /* Constant jacobian is all zeros: size x n_vars with 0 nonzeros.
-     * new_csr_matrix uses calloc for row pointers, so they're already 0. */
+    /* Zero jacobian: size x n_vars with 0 nonzeros. */
     node->jacobian = new_csr_matrix(node->size, node->n_vars, 0);
 }
 
 static void eval_jacobian(expr *node)
 {
-    /* Constant jacobian never changes - nothing to evaluate */
     (void) node;
 }
 
 static void wsum_hess_init_impl(expr *node)
 {
-    /* Constant Hessian is all zeros: n_vars x n_vars with 0 nonzeros. */
+    /* Zero Hessian: n_vars x n_vars with 0 nonzeros. */
     node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, 0);
 }
 
 static void eval_wsum_hess(expr *node, const double *w)
 {
-    /* Constant Hessian is always zero - nothing to compute */
     (void) node;
     (void) w;
 }
@@ -59,12 +58,22 @@ static bool is_affine(const expr *node)
     return true;
 }
 
-expr *new_constant(int d1, int d2, int n_vars, const double *values)
+expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *values)
 {
-    expr *node = (expr *) SP_CALLOC(1, sizeof(expr));
+    parameter_expr *pnode = (parameter_expr *) SP_CALLOC(1, sizeof(parameter_expr));
+    expr *node = &pnode->base;
     init_expr(node, d1, d2, n_vars, forward, jacobian_init_impl, eval_jacobian,
               is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
-    memcpy(node->value, values, node->size * sizeof(double));
 
+    // TODO we should assert that the values array has the correct size.
+    pnode->param_id = param_id;
+
+    if (values == NULL)
+    {
+        fprintf(stderr, "Parameter values should always be set, this is a bug and"
+                        " should be reported\n");
+        exit(1);
+    }
+    memcpy(node->value, values, node->size * sizeof(double));
     return node;
 }