From ff0110fc11ebf2a107e728202c32a8d88c706392 Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Wed, 8 Apr 2026 21:00:27 -0400 Subject: [PATCH 01/12] add initial attempt for fixing parameters and broadcasting --- src/atoms/affine/left_matmul.c | 8 ++++++++ src/atoms/affine/scalar_mult.c | 12 +++++++++++- src/atoms/affine/vector_mult.c | 12 +++++++++++- tests/all_tests.c | 6 ++++++ 4 files changed, 36 insertions(+), 2 deletions(-) diff --git a/src/atoms/affine/left_matmul.c b/src/atoms/affine/left_matmul.c index 9ff397e..8da9e3c 100644 --- a/src/atoms/affine/left_matmul.c +++ b/src/atoms/affine/left_matmul.c @@ -63,6 +63,14 @@ static void refresh_param_values(left_matmul_expr *lnode) static void forward(expr *node, const double *u) { left_matmul_expr *lnode = (left_matmul_expr *) node; + + /* Refresh param_source expression tree if parameters changed.*/ + if (lnode->param_source != NULL && lnode->base.needs_parameter_refresh) + { + /* pass NULL to forward: constant param_source never depends on u */ + lnode->param_source->forward(lnode->param_source, NULL); + } + refresh_param_values(lnode); expr *x = node->left; diff --git a/src/atoms/affine/scalar_mult.c b/src/atoms/affine/scalar_mult.c index baaf16e..57c3335 100644 --- a/src/atoms/affine/scalar_mult.c +++ b/src/atoms/affine/scalar_mult.c @@ -28,7 +28,17 @@ static void forward(expr *node, const double *u) { expr *child = node->left; - double a = ((scalar_mult_expr *) node)->param_source->value[0]; + scalar_mult_expr *snode = (scalar_mult_expr *) node; + + /* Refresh param_source expression tree if parameters changed.*/ + if (node->needs_parameter_refresh) + { + /* pass NULL to forward: constant param_source never depends on u */ + snode->param_source->forward(snode->param_source, NULL); + node->needs_parameter_refresh = false; + } + + double a = snode->param_source->value[0]; /* child's forward pass */ child->forward(child, u); diff --git a/src/atoms/affine/vector_mult.c b/src/atoms/affine/vector_mult.c index 116f012..d256106 100644 --- a/src/atoms/affine/vector_mult.c +++ b/src/atoms/affine/vector_mult.c @@ -28,7 +28,17 @@ static void forward(expr *node, const double *u) { expr *child = node->left; - const double *a = ((vector_mult_expr *) node)->param_source->value; + vector_mult_expr *vnode = (vector_mult_expr *) node; + + /* Refresh param_source expression tree if parameters changed.*/ + if (node->needs_parameter_refresh) + { + /* pass NULL to forward: constant param_source never depends on u */ + vnode->param_source->forward(vnode->param_source, NULL); + node->needs_parameter_refresh = false; + } + + const double *a = vnode->param_source->value; /* child's forward pass */ child->forward(child, u); diff --git a/tests/all_tests.c b/tests/all_tests.c index 3ad1b7c..ac69c16 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -49,6 +49,7 @@ #include "jacobian_tests/other/test_prod_axis_zero.h" #include "jacobian_tests/other/test_quad_form.h" #include "numerical_diff/test_numerical_diff.h" +#include "problem/test_param_broadcast.h" #include "problem/test_param_prob.h" #include "problem/test_problem.h" #include "utils/test_cblas.h" @@ -358,6 +359,11 @@ int main(void) mu_run_test(test_param_right_matmul_rectangular, tests_run); mu_run_test(test_param_shared_left_matmul_problem, tests_run); mu_run_test(test_param_fixed_skip_in_update, tests_run); + + printf("\n--- Parameter + Broadcast Tests ---\n"); + mu_run_test(test_param_broadcast_vector_mult, tests_run); + mu_run_test(test_param_sum_scalar_mult, tests_run); + mu_run_test(test_param_broadcast_left_matmul, tests_run); #endif /* PROFILE_ONLY */ #ifdef PROFILE_ONLY From b4fa176eac8d98115e2b0e8659cfad7fbeb2d0c9 Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Wed, 8 Apr 2026 21:05:10 -0400 Subject: [PATCH 02/12] add test for params with broadcast --- tests/problem/test_param_broadcast.h | 213 +++++++++++++++++++++++++++ 1 file changed, 213 insertions(+) create mode 100644 tests/problem/test_param_broadcast.h diff --git a/tests/problem/test_param_broadcast.h b/tests/problem/test_param_broadcast.h new file mode 100644 index 0000000..6d20ead --- /dev/null +++ b/tests/problem/test_param_broadcast.h @@ -0,0 +1,213 @@ +#ifndef TEST_PARAM_BROADCAST_H +#define TEST_PARAM_BROADCAST_H + +#include +#include +#include + +#include "atoms/affine.h" +#include "expr.h" +#include "minunit.h" +#include "problem.h" +#include "subexpr.h" +#include "test_helpers.h" + +/* + * Tests for param_source being a non-trivial expression tree + * (broadcast, sum, etc.) rather than a bare parameter node. + * + * The key invariant: after problem_update_params, consuming nodes must + * call param_source->forward() to propagate updated parameter values + * through any intermediate nodes before reading param_source->value. + */ + +/* ------------------------------------------------------------------ */ +/* vector_mult(broadcast(parameter), variable) */ +/* */ +/* p is (1,3), broadcast ROW to (2,3), x is (2,3) variable. */ +/* broadcast(p) = [p0,p0, p1,p1, p2,p2] in column-major. */ +/* constraint_k = broadcast(p)_k * x_k (elementwise). */ +/* Jacobian = diag(broadcast(p)). */ +/* ------------------------------------------------------------------ */ +const char *test_param_broadcast_vector_mult(void) +{ + int n = 6; + + expr *x = new_variable(2, 3, 0, n); + expr *objective = new_sum(x, -1); + + expr *p_param = new_parameter(1, 3, 0, n, NULL); + expr *p_bcast = new_broadcast(p_param, 2, 3); + expr *constraint = new_vector_mult(p_bcast, x); + + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, false); + + expr *param_nodes[1] = {p_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + double x_vals[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + + /* --- theta = [1, 2, 3] ---------------------------------------- */ + /* broadcast = [1,1, 2,2, 3,3] */ + /* constraint = [1*1, 1*2, 2*3, 2*4, 3*5, 3*6] */ + /* = [1, 2, 6, 8, 15, 18] */ + double theta1[3] = {1.0, 2.0, 3.0}; + problem_update_params(prob, theta1); + problem_constraint_forward(prob, x_vals); + + double fwd1[6] = {1.0, 2.0, 6.0, 8.0, 15.0, 18.0}; + mu_assert("bcast vmul fwd1", cmp_double_array(prob->constraint_values, fwd1, 6)); + + problem_jacobian(prob); + double jac1[6] = {1.0, 1.0, 2.0, 2.0, 3.0, 3.0}; + mu_assert("bcast vmul jac1", cmp_double_array(prob->jacobian->x, jac1, 6)); + + /* --- theta = [10, 20, 30] ------------------------------------- */ + /* broadcast = [10,10, 20,20, 30,30] */ + /* constraint = [10, 20, 60, 80, 150, 180] */ + double theta2[3] = {10.0, 20.0, 30.0}; + problem_update_params(prob, theta2); + problem_constraint_forward(prob, x_vals); + + double fwd2[6] = {10.0, 20.0, 60.0, 80.0, 150.0, 180.0}; + mu_assert("bcast vmul fwd2", cmp_double_array(prob->constraint_values, fwd2, 6)); + + problem_jacobian(prob); + double jac2[6] = {10.0, 10.0, 20.0, 20.0, 30.0, 30.0}; + mu_assert("bcast vmul jac2", cmp_double_array(prob->jacobian->x, jac2, 6)); + + free_problem(prob); + return 0; +} + +/* ------------------------------------------------------------------ */ +/* scalar_mult(sum(parameter), variable) */ +/* */ +/* p is (2,1), sum(p) is (1,1) scalar = p0+p1. */ +/* x is (3,1) variable. */ +/* constraint_i = (p0+p1) * x_i. */ +/* Jacobian = diag(p0+p1, p0+p1, p0+p1). */ +/* ------------------------------------------------------------------ */ +const char *test_param_sum_scalar_mult(void) +{ + int n = 3; + + expr *x = new_variable(3, 1, 0, n); + expr *objective = new_sum(x, -1); + + expr *p_param = new_parameter(2, 1, 0, n, NULL); + expr *p_sum = new_sum(p_param, -1); + expr *constraint = new_scalar_mult(p_sum, x); + + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, false); + + expr *param_nodes[1] = {p_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + double x_vals[3] = {1.0, 2.0, 3.0}; + + /* --- theta = [1, 2], sum(p) = 3 ------------------------------- */ + /* constraint = [3, 6, 9] */ + double theta1[2] = {1.0, 2.0}; + problem_update_params(prob, theta1); + problem_constraint_forward(prob, x_vals); + + double fwd1[3] = {3.0, 6.0, 9.0}; + mu_assert("sum smul fwd1", cmp_double_array(prob->constraint_values, fwd1, 3)); + + problem_jacobian(prob); + double jac1[3] = {3.0, 3.0, 3.0}; + mu_assert("sum smul jac1", cmp_double_array(prob->jacobian->x, jac1, 3)); + + /* --- theta = [5, 10], sum(p) = 15 ----------------------------- */ + /* constraint = [15, 30, 45] */ + double theta2[2] = {5.0, 10.0}; + problem_update_params(prob, theta2); + problem_constraint_forward(prob, x_vals); + + double fwd2[3] = {15.0, 30.0, 45.0}; + mu_assert("sum smul fwd2", cmp_double_array(prob->constraint_values, fwd2, 3)); + + problem_jacobian(prob); + double jac2[3] = {15.0, 15.0, 15.0}; + mu_assert("sum smul jac2", cmp_double_array(prob->jacobian->x, jac2, 3)); + + free_problem(prob); + return 0; +} + +/* ------------------------------------------------------------------ */ +/* left_matmul_dense(broadcast(parameter), variable) */ +/* */ +/* p is (1,2), broadcast ROW to (3,2). */ +/* broadcast stores column-major: [a,a,a, b,b,b]. */ +/* refresh_dense_left interprets this as A in column-major, giving */ +/* A = [[a,b],[a,b],[a,b]] (3x2). */ +/* x is (2,1) variable. */ +/* constraint = A @ x, a (3,1) vector. */ +/* Jacobian = A (each row = [a, b]). */ +/* ------------------------------------------------------------------ */ +const char *test_param_broadcast_left_matmul(void) +{ + int n = 2; + + expr *x = new_variable(2, 1, 0, n); + expr *objective = new_sum(x, -1); + + expr *p_param = new_parameter(1, 2, 0, n, NULL); + expr *p_bcast = new_broadcast(p_param, 3, 2); + + /* initial data (overwritten on first param refresh) */ + double A_init[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + expr *constraint = new_left_matmul_dense(p_bcast, x, 3, 2, A_init); + + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, false); + + expr *param_nodes[1] = {p_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + double x_vals[2] = {3.0, 4.0}; + int Ap[4] = {0, 2, 4, 6}; + int Ai[6] = {0, 1, 0, 1, 0, 1}; + + /* --- theta = [1, 2] ------------------------------------------- */ + /* A = [[1,2],[1,2],[1,2]], A@x = [11, 11, 11] */ + double theta1[2] = {1.0, 2.0}; + problem_update_params(prob, theta1); + problem_constraint_forward(prob, x_vals); + + double fwd1[3] = {11.0, 11.0, 11.0}; + mu_assert("bcast lmul fwd1", cmp_double_array(prob->constraint_values, fwd1, 3)); + + problem_jacobian(prob); + double jac1[6] = {1.0, 2.0, 1.0, 2.0, 1.0, 2.0}; + mu_assert("bcast lmul jac1", cmp_double_array(prob->jacobian->x, jac1, 6)); + mu_assert("bcast lmul rows1", cmp_int_array(prob->jacobian->p, Ap, 4)); + mu_assert("bcast lmul cols1", cmp_int_array(prob->jacobian->i, Ai, 6)); + + /* --- theta = [5, 10] ------------------------------------------ */ + /* A = [[5,10],[5,10],[5,10]], A@x = [55, 55, 55] */ + double theta2[2] = {5.0, 10.0}; + problem_update_params(prob, theta2); + problem_constraint_forward(prob, x_vals); + + double fwd2[3] = {55.0, 55.0, 55.0}; + mu_assert("bcast lmul fwd2", cmp_double_array(prob->constraint_values, fwd2, 3)); + + problem_jacobian(prob); + double jac2[6] = {5.0, 10.0, 5.0, 10.0, 5.0, 10.0}; + mu_assert("bcast lmul jac2", cmp_double_array(prob->jacobian->x, jac2, 6)); + mu_assert("bcast lmul rows2", cmp_int_array(prob->jacobian->p, Ap, 4)); + mu_assert("bcast lmul cols2", cmp_int_array(prob->jacobian->i, Ai, 6)); + + free_problem(prob); + return 0; +} + +#endif /* TEST_PARAM_BROADCAST_H */ From 3310a9195e766ec812a4a2f883627ccacb224c2b Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Wed, 8 Apr 2026 21:19:22 -0400 Subject: [PATCH 03/12] cleanup test to fit more with other tests style --- tests/problem/test_param_broadcast.h | 201 +++++++++++---------------- 1 file changed, 85 insertions(+), 116 deletions(-) diff --git a/tests/problem/test_param_broadcast.h b/tests/problem/test_param_broadcast.h index 6d20ead..be7c52a 100644 --- a/tests/problem/test_param_broadcast.h +++ b/tests/problem/test_param_broadcast.h @@ -12,199 +12,168 @@ #include "subexpr.h" #include "test_helpers.h" -/* - * Tests for param_source being a non-trivial expression tree - * (broadcast, sum, etc.) rather than a bare parameter node. - * - * The key invariant: after problem_update_params, consuming nodes must - * call param_source->forward() to propagate updated parameter values - * through any intermediate nodes before reading param_source->value. - */ - -/* ------------------------------------------------------------------ */ -/* vector_mult(broadcast(parameter), variable) */ -/* */ -/* p is (1,3), broadcast ROW to (2,3), x is (2,3) variable. */ -/* broadcast(p) = [p0,p0, p1,p1, p2,p2] in column-major. */ -/* constraint_k = broadcast(p)_k * x_k (elementwise). */ -/* Jacobian = diag(broadcast(p)). */ -/* ------------------------------------------------------------------ */ const char *test_param_broadcast_vector_mult(void) { int n = 6; + /* minimize sum(x) subject to broadcast(p) ∘ x, with p parameter */ expr *x = new_variable(2, 3, 0, n); expr *objective = new_sum(x, -1); - expr *p_param = new_parameter(1, 3, 0, n, NULL); expr *p_bcast = new_broadcast(p_param, 2, 3); expr *constraint = new_vector_mult(p_bcast, x); - expr *constraints[1] = {constraint}; problem *prob = new_problem(objective, constraints, 1, false); + /* register parameters and fill sparsity patterns */ expr *param_nodes[1] = {p_param}; problem_register_params(prob, param_nodes, 1); problem_init_derivatives(prob); + /* point for evaluating */ double x_vals[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; - /* --- theta = [1, 2, 3] ---------------------------------------- */ - /* broadcast = [1,1, 2,2, 3,3] */ - /* constraint = [1*1, 1*2, 2*3, 2*4, 3*5, 3*6] */ - /* = [1, 2, 6, 8, 15, 18] */ - double theta1[3] = {1.0, 2.0, 3.0}; - problem_update_params(prob, theta1); + /* test 1: p=[1,2,3] */ + double theta[3] = {1.0, 2.0, 3.0}; + problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); - - double fwd1[6] = {1.0, 2.0, 6.0, 8.0, 15.0, 18.0}; - mu_assert("bcast vmul fwd1", cmp_double_array(prob->constraint_values, fwd1, 6)); - + double constrs[6] = {1.0, 2.0, 6.0, 8.0, 15.0, 18.0}; problem_jacobian(prob); - double jac1[6] = {1.0, 1.0, 2.0, 2.0, 3.0, 3.0}; - mu_assert("bcast vmul jac1", cmp_double_array(prob->jacobian->x, jac1, 6)); - - /* --- theta = [10, 20, 30] ------------------------------------- */ - /* broadcast = [10,10, 20,20, 30,30] */ - /* constraint = [10, 20, 60, 80, 150, 180] */ - double theta2[3] = {10.0, 20.0, 30.0}; - problem_update_params(prob, theta2); + double jac_x[6] = {1.0, 1.0, 2.0, 2.0, 3.0, 3.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + + /* test 2: p=[10,20,30] */ + theta[0] = 10.0; + theta[1] = 20.0; + theta[2] = 30.0; + problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); - - double fwd2[6] = {10.0, 20.0, 60.0, 80.0, 150.0, 180.0}; - mu_assert("bcast vmul fwd2", cmp_double_array(prob->constraint_values, fwd2, 6)); - problem_jacobian(prob); - double jac2[6] = {10.0, 10.0, 20.0, 20.0, 30.0, 30.0}; - mu_assert("bcast vmul jac2", cmp_double_array(prob->jacobian->x, jac2, 6)); + constrs[0] = 10.0; + constrs[1] = 20.0; + constrs[2] = 60.0; + constrs[3] = 80.0; + constrs[4] = 150.0; + constrs[5] = 180.0; + jac_x[0] = 10.0; + jac_x[1] = 10.0; + jac_x[2] = 20.0; + jac_x[3] = 20.0; + jac_x[4] = 30.0; + jac_x[5] = 30.0; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); free_problem(prob); return 0; } -/* ------------------------------------------------------------------ */ -/* scalar_mult(sum(parameter), variable) */ -/* */ -/* p is (2,1), sum(p) is (1,1) scalar = p0+p1. */ -/* x is (3,1) variable. */ -/* constraint_i = (p0+p1) * x_i. */ -/* Jacobian = diag(p0+p1, p0+p1, p0+p1). */ -/* ------------------------------------------------------------------ */ const char *test_param_sum_scalar_mult(void) { int n = 3; + /* minimize sum(x) subject to sum(p) * x, with p parameter */ expr *x = new_variable(3, 1, 0, n); expr *objective = new_sum(x, -1); - expr *p_param = new_parameter(2, 1, 0, n, NULL); expr *p_sum = new_sum(p_param, -1); expr *constraint = new_scalar_mult(p_sum, x); - expr *constraints[1] = {constraint}; problem *prob = new_problem(objective, constraints, 1, false); + /* register parameters and fill sparsity patterns */ expr *param_nodes[1] = {p_param}; problem_register_params(prob, param_nodes, 1); problem_init_derivatives(prob); + /* point for evaluating */ double x_vals[3] = {1.0, 2.0, 3.0}; - /* --- theta = [1, 2], sum(p) = 3 ------------------------------- */ - /* constraint = [3, 6, 9] */ - double theta1[2] = {1.0, 2.0}; - problem_update_params(prob, theta1); + /* test 1: p=[1,2], sum(p)=3 */ + double theta[2] = {1.0, 2.0}; + problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); - - double fwd1[3] = {3.0, 6.0, 9.0}; - mu_assert("sum smul fwd1", cmp_double_array(prob->constraint_values, fwd1, 3)); - + double constrs[3] = {3.0, 6.0, 9.0}; problem_jacobian(prob); - double jac1[3] = {3.0, 3.0, 3.0}; - mu_assert("sum smul jac1", cmp_double_array(prob->jacobian->x, jac1, 3)); - - /* --- theta = [5, 10], sum(p) = 15 ----------------------------- */ - /* constraint = [15, 30, 45] */ - double theta2[2] = {5.0, 10.0}; - problem_update_params(prob, theta2); + double jac_x[3] = {3.0, 3.0, 3.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 3)); + + /* test 2: p=[5,10], sum(p)=15 */ + theta[0] = 5.0; + theta[1] = 10.0; + problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); - - double fwd2[3] = {15.0, 30.0, 45.0}; - mu_assert("sum smul fwd2", cmp_double_array(prob->constraint_values, fwd2, 3)); - problem_jacobian(prob); - double jac2[3] = {15.0, 15.0, 15.0}; - mu_assert("sum smul jac2", cmp_double_array(prob->jacobian->x, jac2, 3)); + constrs[0] = 15.0; + constrs[1] = 30.0; + constrs[2] = 45.0; + jac_x[0] = 15.0; + jac_x[1] = 15.0; + jac_x[2] = 15.0; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 3)); free_problem(prob); return 0; } -/* ------------------------------------------------------------------ */ -/* left_matmul_dense(broadcast(parameter), variable) */ -/* */ -/* p is (1,2), broadcast ROW to (3,2). */ -/* broadcast stores column-major: [a,a,a, b,b,b]. */ -/* refresh_dense_left interprets this as A in column-major, giving */ -/* A = [[a,b],[a,b],[a,b]] (3x2). */ -/* x is (2,1) variable. */ -/* constraint = A @ x, a (3,1) vector. */ -/* Jacobian = A (each row = [a, b]). */ -/* ------------------------------------------------------------------ */ const char *test_param_broadcast_left_matmul(void) { int n = 2; + /* minimize sum(x) subject to broadcast(p)@x, with p parameter */ expr *x = new_variable(2, 1, 0, n); expr *objective = new_sum(x, -1); - expr *p_param = new_parameter(1, 2, 0, n, NULL); expr *p_bcast = new_broadcast(p_param, 3, 2); - - /* initial data (overwritten on first param refresh) */ - double A_init[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; - expr *constraint = new_left_matmul_dense(p_bcast, x, 3, 2, A_init); - + double Ax[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + expr *constraint = new_left_matmul_dense(p_bcast, x, 3, 2, Ax); expr *constraints[1] = {constraint}; problem *prob = new_problem(objective, constraints, 1, false); + /* register parameters and fill sparsity patterns */ expr *param_nodes[1] = {p_param}; problem_register_params(prob, param_nodes, 1); problem_init_derivatives(prob); + /* point for evaluating and utilities for test */ double x_vals[2] = {3.0, 4.0}; int Ap[4] = {0, 2, 4, 6}; int Ai[6] = {0, 1, 0, 1, 0, 1}; - /* --- theta = [1, 2] ------------------------------------------- */ - /* A = [[1,2],[1,2],[1,2]], A@x = [11, 11, 11] */ - double theta1[2] = {1.0, 2.0}; - problem_update_params(prob, theta1); + /* test 1: p=[1,2] */ + double theta[2] = {1.0, 2.0}; + problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); - - double fwd1[3] = {11.0, 11.0, 11.0}; - mu_assert("bcast lmul fwd1", cmp_double_array(prob->constraint_values, fwd1, 3)); - + double constrs[3] = {11.0, 11.0, 11.0}; problem_jacobian(prob); - double jac1[6] = {1.0, 2.0, 1.0, 2.0, 1.0, 2.0}; - mu_assert("bcast lmul jac1", cmp_double_array(prob->jacobian->x, jac1, 6)); - mu_assert("bcast lmul rows1", cmp_int_array(prob->jacobian->p, Ap, 4)); - mu_assert("bcast lmul cols1", cmp_int_array(prob->jacobian->i, Ai, 6)); - - /* --- theta = [5, 10] ------------------------------------------ */ - /* A = [[5,10],[5,10],[5,10]], A@x = [55, 55, 55] */ - double theta2[2] = {5.0, 10.0}; - problem_update_params(prob, theta2); + double jac_x[6] = {1.0, 2.0, 1.0, 2.0, 1.0, 2.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 4)); + mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 6)); + + /* test 2: p=[5,10] */ + theta[0] = 5.0; + theta[1] = 10.0; + problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); - - double fwd2[3] = {55.0, 55.0, 55.0}; - mu_assert("bcast lmul fwd2", cmp_double_array(prob->constraint_values, fwd2, 3)); - problem_jacobian(prob); - double jac2[6] = {5.0, 10.0, 5.0, 10.0, 5.0, 10.0}; - mu_assert("bcast lmul jac2", cmp_double_array(prob->jacobian->x, jac2, 6)); - mu_assert("bcast lmul rows2", cmp_int_array(prob->jacobian->p, Ap, 4)); - mu_assert("bcast lmul cols2", cmp_int_array(prob->jacobian->i, Ai, 6)); + constrs[0] = 55.0; + constrs[1] = 55.0; + constrs[2] = 55.0; + jac_x[0] = 5.0; + jac_x[1] = 10.0; + jac_x[2] = 5.0; + jac_x[3] = 10.0; + jac_x[4] = 5.0; + jac_x[5] = 10.0; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 4)); + mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 6)); free_problem(prob); return 0; From 6176be6863c3897f3266426138e174ecdadadf77 Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Thu, 9 Apr 2026 13:45:53 -0400 Subject: [PATCH 04/12] some progress on supporting backwards compatible constants --- src/atoms/affine/parameter.c | 11 +- src/atoms/affine/scalar_mult.c | 11 +- src/atoms/affine/vector_mult.c | 11 +- src/problem.c | 3 + tests/all_tests.c | 4 +- tests/problem/test_param_broadcast.h | 145 ++++----------------------- tests/problem/test_param_prob.h | 47 +++++---- 7 files changed, 65 insertions(+), 167 deletions(-) diff --git a/src/atoms/affine/parameter.c b/src/atoms/affine/parameter.c index fa96377..a951aae 100644 --- a/src/atoms/affine/parameter.c +++ b/src/atoms/affine/parameter.c @@ -18,6 +18,7 @@ #include "atoms/affine.h" #include "subexpr.h" #include "utils/tracked_alloc.h" +#include #include #include @@ -65,11 +66,13 @@ expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *valu is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL); pnode->param_id = param_id; - - if (values != NULL) + + if (values == NULL) { - memcpy(node->value, values, node->size * sizeof(double)); + 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; } diff --git a/src/atoms/affine/scalar_mult.c b/src/atoms/affine/scalar_mult.c index 57c3335..ef3f6bc 100644 --- a/src/atoms/affine/scalar_mult.c +++ b/src/atoms/affine/scalar_mult.c @@ -30,13 +30,9 @@ static void forward(expr *node, const double *u) expr *child = node->left; scalar_mult_expr *snode = (scalar_mult_expr *) node; - /* Refresh param_source expression tree if parameters changed.*/ - if (node->needs_parameter_refresh) - { - /* pass NULL to forward: constant param_source never depends on u */ - snode->param_source->forward(snode->param_source, NULL); - node->needs_parameter_refresh = false; - } + /* call forward for param_source expr tree + ex: broadcast(param) or promote(const)*/ + snode->param_source->forward(snode->param_source, NULL); double a = snode->param_source->value[0]; @@ -128,6 +124,7 @@ expr *new_scalar_mult(expr *param_node, expr *child) mult_node->param_source = param_node; expr_retain(param_node); + //node->needs_parameter_refresh = true; return node; } diff --git a/src/atoms/affine/vector_mult.c b/src/atoms/affine/vector_mult.c index d256106..ed40795 100644 --- a/src/atoms/affine/vector_mult.c +++ b/src/atoms/affine/vector_mult.c @@ -30,13 +30,9 @@ static void forward(expr *node, const double *u) expr *child = node->left; vector_mult_expr *vnode = (vector_mult_expr *) node; - /* Refresh param_source expression tree if parameters changed.*/ - if (node->needs_parameter_refresh) - { - /* pass NULL to forward: constant param_source never depends on u */ - vnode->param_source->forward(vnode->param_source, NULL); - node->needs_parameter_refresh = false; - } + /* call forward for param_source expr tree + ex: broadcast(param) or promote(const)*/ + vnode->param_source->forward(vnode->param_source, NULL); const double *a = vnode->param_source->value; @@ -138,6 +134,7 @@ expr *new_vector_mult(expr *param_node, expr *child) vnode->param_source = param_node; expr_retain(param_node); + //node->needs_parameter_refresh = true; return node; } diff --git a/src/problem.c b/src/problem.c index 0419063..28448fe 100644 --- a/src/problem.c +++ b/src/problem.c @@ -369,6 +369,9 @@ void problem_register_params(problem *prob, expr **param_nodes, int n_param_node prob->total_parameter_size = 0; for (int i = 0; i < n_param_nodes; i++) { + // TODO do we need to skip fixed params? maybe we adopt the convention + // that we don't ever register fixed params? + if (((parameter_expr *) param_nodes[i])->param_id == PARAM_FIXED) continue; prob->total_parameter_size += param_nodes[i]->size; } } diff --git a/tests/all_tests.c b/tests/all_tests.c index ac69c16..fb1f9a6 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -361,9 +361,7 @@ int main(void) mu_run_test(test_param_fixed_skip_in_update, tests_run); printf("\n--- Parameter + Broadcast Tests ---\n"); - mu_run_test(test_param_broadcast_vector_mult, tests_run); - mu_run_test(test_param_sum_scalar_mult, tests_run); - mu_run_test(test_param_broadcast_left_matmul, tests_run); + mu_run_test(test_constant_broadcast_vector_mult, tests_run); #endif /* PROFILE_ONLY */ #ifdef PROFILE_ONLY diff --git a/tests/problem/test_param_broadcast.h b/tests/problem/test_param_broadcast.h index be7c52a..3532c88 100644 --- a/tests/problem/test_param_broadcast.h +++ b/tests/problem/test_param_broadcast.h @@ -12,30 +12,24 @@ #include "subexpr.h" #include "test_helpers.h" -const char *test_param_broadcast_vector_mult(void) +const char *test_constant_broadcast_vector_mult(void) { int n = 6; - /* minimize sum(x) subject to broadcast(p) ∘ x, with p parameter */ + /* minimize sum(x) subject to broadcast(c) ∘ x, with c constant */ expr *x = new_variable(2, 3, 0, n); expr *objective = new_sum(x, -1); - expr *p_param = new_parameter(1, 3, 0, n, NULL); - expr *p_bcast = new_broadcast(p_param, 2, 3); - expr *constraint = new_vector_mult(p_bcast, x); + double c_vals[3] = {1.0, 2.0, 3.0}; + expr *c = new_parameter(1, 3, PARAM_FIXED, n, c_vals); + expr *c_bcast = new_broadcast(c, 2, 3); + expr *constraint = new_vector_mult(c_bcast, x); expr *constraints[1] = {constraint}; problem *prob = new_problem(objective, constraints, 1, false); - - /* register parameters and fill sparsity patterns */ - expr *param_nodes[1] = {p_param}; - problem_register_params(prob, param_nodes, 1); problem_init_derivatives(prob); /* point for evaluating */ double x_vals[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; - /* test 1: p=[1,2,3] */ - double theta[3] = {1.0, 2.0, 3.0}; - problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); double constrs[6] = {1.0, 2.0, 6.0, 8.0, 15.0, 18.0}; problem_jacobian(prob); @@ -43,140 +37,37 @@ const char *test_param_broadcast_vector_mult(void) mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); - /* test 2: p=[10,20,30] */ - theta[0] = 10.0; - theta[1] = 20.0; - theta[2] = 30.0; - problem_update_params(prob, theta); - problem_constraint_forward(prob, x_vals); - problem_jacobian(prob); - constrs[0] = 10.0; - constrs[1] = 20.0; - constrs[2] = 60.0; - constrs[3] = 80.0; - constrs[4] = 150.0; - constrs[5] = 180.0; - jac_x[0] = 10.0; - jac_x[1] = 10.0; - jac_x[2] = 20.0; - jac_x[3] = 20.0; - jac_x[4] = 30.0; - jac_x[5] = 30.0; - mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); - mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); - free_problem(prob); return 0; } -const char *test_param_sum_scalar_mult(void) +const char *test_constant_promote_vector_mult(void) { - int n = 3; + int n = 6; - /* minimize sum(x) subject to sum(p) * x, with p parameter */ - expr *x = new_variable(3, 1, 0, n); + /* minimize sum(x) subject to promote(c) ∘ x, with c constant */ + expr *x = new_variable(2, 3, 0, n); expr *objective = new_sum(x, -1); - expr *p_param = new_parameter(2, 1, 0, n, NULL); - expr *p_sum = new_sum(p_param, -1); - expr *constraint = new_scalar_mult(p_sum, x); + double c_vals = 3.0; + expr *c = new_parameter(1, 1, PARAM_FIXED, n, &c_vals); + expr *c_bcast = new_promote(c, 2, 3); + expr *constraint = new_vector_mult(c_bcast, x); expr *constraints[1] = {constraint}; problem *prob = new_problem(objective, constraints, 1, false); - /* register parameters and fill sparsity patterns */ - expr *param_nodes[1] = {p_param}; - problem_register_params(prob, param_nodes, 1); problem_init_derivatives(prob); /* point for evaluating */ - double x_vals[3] = {1.0, 2.0, 3.0}; - - /* test 1: p=[1,2], sum(p)=3 */ - double theta[2] = {1.0, 2.0}; - problem_update_params(prob, theta); - problem_constraint_forward(prob, x_vals); - double constrs[3] = {3.0, 6.0, 9.0}; - problem_jacobian(prob); - double jac_x[3] = {3.0, 3.0, 3.0}; - mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); - mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 3)); - - /* test 2: p=[5,10], sum(p)=15 */ - theta[0] = 5.0; - theta[1] = 10.0; - problem_update_params(prob, theta); - problem_constraint_forward(prob, x_vals); - problem_jacobian(prob); - constrs[0] = 15.0; - constrs[1] = 30.0; - constrs[2] = 45.0; - jac_x[0] = 15.0; - jac_x[1] = 15.0; - jac_x[2] = 15.0; - mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); - mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 3)); - - free_problem(prob); - return 0; -} - -const char *test_param_broadcast_left_matmul(void) -{ - int n = 2; - - /* minimize sum(x) subject to broadcast(p)@x, with p parameter */ - expr *x = new_variable(2, 1, 0, n); - expr *objective = new_sum(x, -1); - expr *p_param = new_parameter(1, 2, 0, n, NULL); - expr *p_bcast = new_broadcast(p_param, 3, 2); - double Ax[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; - expr *constraint = new_left_matmul_dense(p_bcast, x, 3, 2, Ax); - expr *constraints[1] = {constraint}; - problem *prob = new_problem(objective, constraints, 1, false); - - /* register parameters and fill sparsity patterns */ - expr *param_nodes[1] = {p_param}; - problem_register_params(prob, param_nodes, 1); - problem_init_derivatives(prob); - - /* point for evaluating and utilities for test */ - double x_vals[2] = {3.0, 4.0}; - int Ap[4] = {0, 2, 4, 6}; - int Ai[6] = {0, 1, 0, 1, 0, 1}; - - /* test 1: p=[1,2] */ - double theta[2] = {1.0, 2.0}; - problem_update_params(prob, theta); - problem_constraint_forward(prob, x_vals); - double constrs[3] = {11.0, 11.0, 11.0}; - problem_jacobian(prob); - double jac_x[6] = {1.0, 2.0, 1.0, 2.0, 1.0, 2.0}; - mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); - mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); - mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 4)); - mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 6)); + double x_vals[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; - /* test 2: p=[5,10] */ - theta[0] = 5.0; - theta[1] = 10.0; - problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); + double constrs[6] = {3.0, 6.0, 9.0, 12.0, 15.0, 18.0}; problem_jacobian(prob); - constrs[0] = 55.0; - constrs[1] = 55.0; - constrs[2] = 55.0; - jac_x[0] = 5.0; - jac_x[1] = 10.0; - jac_x[2] = 5.0; - jac_x[3] = 10.0; - jac_x[4] = 5.0; - jac_x[5] = 10.0; - mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + double jac_x[6] = {3.0, 3.0, 3.0, 3.0, 3.0, 3.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); - mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 4)); - mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 6)); free_problem(prob); return 0; } - #endif /* TEST_PARAM_BROADCAST_H */ diff --git a/tests/problem/test_param_prob.h b/tests/problem/test_param_prob.h index f5141e5..59c1c7a 100644 --- a/tests/problem/test_param_prob.h +++ b/tests/problem/test_param_prob.h @@ -20,7 +20,8 @@ const char *test_param_scalar_mult_problem(void) /* minimize a * sum(log(x)), with a parameter */ expr *x = new_variable(2, 1, 0, n); expr *log_x = new_log(x); - expr *a_param = new_parameter(1, 1, 0, n, NULL); + double theta[1] = {3.0}; + expr *a_param = new_parameter(1, 1, 0, n, theta); expr *scaled = new_scalar_mult(a_param, log_x); expr *objective = new_sum(scaled, -1); problem *prob = new_problem(objective, NULL, 0, false); @@ -34,8 +35,6 @@ const char *test_param_scalar_mult_problem(void) double x_vals[2] = {1.0, 2.0}; /* test 1: a=3 */ - double theta[1] = {3.0}; - problem_update_params(prob, theta); double obj_val = problem_objective_forward(prob, x_vals); problem_gradient(prob); double expected_obj = 3.0 * log(2.0); @@ -66,7 +65,8 @@ const char *test_param_vector_mult_problem(void) /* minimize sum(x) subject to p ∘ x, with p parameter */ expr *x = new_variable(2, 1, 0, n); expr *objective = new_sum(x, -1); - expr *p_param = new_parameter(2, 1, 0, n, NULL); + double theta[2] = {3.0, 4.0}; + expr *p_param = new_parameter(2, 1, 0, n, theta); expr *constraint = new_vector_mult(p_param, x); expr *constraints[1] = {constraint}; problem *prob = new_problem(objective, constraints, 1, false); @@ -83,8 +83,6 @@ const char *test_param_vector_mult_problem(void) double Ax[2] = {3.0, 4.0}; /* test 1: p=[3,4] */ - double theta[2] = {3.0, 4.0}; - problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); double constrs[2] = {3.0, 8.0}; problem_jacobian(prob); @@ -121,7 +119,8 @@ const char *test_param_left_matmul_problem(void) /* minimize sum(x) subject to Ax = ?, with A parameter */ expr *x = new_variable(2, 1, 0, n); expr *objective = new_sum(x, -1); - expr *A_param = new_parameter(2, 2, 0, n, NULL); + double theta[4] = {1.0, 0.0, 0.0, 1.0}; + expr *A_param = new_parameter(2, 2, 0, n, theta); /* dense 2x2 matrix */ double Ax[4] = {2.0, 0.0, 0.0, 1.0}; @@ -149,7 +148,7 @@ const char *test_param_left_matmul_problem(void) mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 4)); /* test 2: A = [[1,2],[3,4]] (column-major [1,3,2,4]) */ - double theta[4] = {1.0, 3.0, 2.0, 4.0}; + theta[0] = 1.0; theta[1] = 3.0; theta[2] = 2.0; theta[3] = 4.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); @@ -177,7 +176,8 @@ const char *test_param_right_matmul_problem(void) /* minimize sum(x) subject to xA = ?, with A parameter */ expr *x = new_variable(1, 2, 0, n); expr *objective = new_sum(x, -1); - expr *A_param = new_parameter(2, 2, 0, n, NULL); + double theta[4] = {1.0, 0.0, 0.0, 1.0}; + expr *A_param = new_parameter(2, 2, 0, n, theta); /* dense 2x2 matrix */ double Ax[4] = {2.0, 0.0, 0.0, 1.0}; @@ -205,7 +205,10 @@ const char *test_param_right_matmul_problem(void) mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 4)); /* test 2: A = [[1,2],[3,4]] (column-major [1,3,2,4]) */ - double theta[4] = {1.0, 3.0, 2.0, 4.0}; + theta[0] = 1.0; + theta[1] = 3.0; + theta[2] = 2.0; + theta[3] = 4.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); @@ -238,7 +241,8 @@ const char *test_param_fixed_skip_in_update(void) expr *a_log = new_scalar_mult(a_param, log_x); expr *sum_a_log = new_sum(a_log, -1); - expr *b_param = new_parameter(1, 1, 0, n, NULL); + double b_val = 0.0; + expr *b_param = new_parameter(1, 1, 0, n, &b_val); expr *b_x = new_scalar_mult(b_param, x); expr *sum_b_x = new_sum(b_x, -1); @@ -246,8 +250,8 @@ const char *test_param_fixed_skip_in_update(void) problem *prob = new_problem(objective, NULL, 0, false); /* register parameters and fill sparsity patterns */ - expr *param_nodes[2] = {a_param, b_param}; - problem_register_params(prob, param_nodes, 2); + expr *param_nodes[1] = {b_param}; + problem_register_params(prob, param_nodes, 1); problem_init_derivatives(prob); /* point for evaluating */ @@ -288,7 +292,8 @@ const char *test_param_left_matmul_rectangular(void) /* minimize sum(x) subject to Ax = ?, with A parameter (3x2) */ expr *x = new_variable(2, 1, 0, n); expr *objective = new_sum(x, -1); - expr *A_param = new_parameter(3, 2, 0, n, NULL); + double theta[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + expr *A_param = new_parameter(3, 2, 0, n, theta); /* dense 3x2 matrix */ double Ax[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; @@ -316,7 +321,8 @@ const char *test_param_left_matmul_rectangular(void) mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 6)); /* test 2: A = [[7,8],[9,10],[11,12]] (column-major [7,9,11,8,10,12]) */ - double theta[6] = {7.0, 9.0, 11.0, 8.0, 10.0, 12.0}; + theta[0] = 7.0; theta[1] = 9.0; theta[2] = 11.0; + theta[3] = 8.0; theta[4] = 10.0; theta[5] = 12.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); @@ -347,7 +353,8 @@ const char *test_param_right_matmul_rectangular(void) /* minimize sum(x) subject to xA = ?, with A parameter (2x3) */ expr *x = new_variable(1, 2, 0, n); expr *objective = new_sum(x, -1); - expr *A_param = new_parameter(2, 3, 0, n, NULL); + double theta[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + expr *A_param = new_parameter(2, 3, 0, n, theta); /* dense 2x3 matrix */ double Ax[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; @@ -376,7 +383,8 @@ const char *test_param_right_matmul_rectangular(void) mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); /* test 2: A = [[7,8,9],[10,11,12]] (column-major [7,10,8,11,9,12]) */ - double theta[6] = {7.0, 10.0, 8.0, 11.0, 9.0, 12.0}; + theta[0] = 7.0; theta[1] = 10.0; theta[2] = 8.0; + theta[3] = 11.0; theta[4] = 9.0; theta[5] = 12.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); @@ -408,7 +416,8 @@ const char *test_param_shared_left_matmul_problem(void) expr *x = new_variable(2, 1, 0, n); expr *y = new_variable(2, 1, 2, n); expr *objective = new_sum(x, -1); - expr *A_param = new_parameter(2, 2, 0, n, NULL); + double theta[4] = {1.0, 0.0, 0.0, 1.0}; + expr *A_param = new_parameter(2, 2, 0, n, theta); /* dense 2x2 identity */ double Ax[4] = {1.0, 0.0, 0.0, 1.0}; @@ -438,7 +447,7 @@ const char *test_param_shared_left_matmul_problem(void) mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 8)); /* test 2: A = [[1,2],[3,4]] (column-major [1,3,2,4]) */ - double theta[4] = {1.0, 3.0, 2.0, 4.0}; + theta[0] = 1.0; theta[1] = 3.0; theta[2] = 2.0; theta[3] = 4.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); From 9153a88c7d074a8496ed7e43d65a84faee54c81c Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Thu, 9 Apr 2026 14:02:00 -0400 Subject: [PATCH 05/12] add some parameter broadcast tests as well --- tests/all_tests.c | 3 + tests/problem/test_param_broadcast.h | 84 ++++++++++++++++++++++++++++ 2 files changed, 87 insertions(+) diff --git a/tests/all_tests.c b/tests/all_tests.c index fb1f9a6..96a6ff7 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -362,6 +362,9 @@ int main(void) printf("\n--- Parameter + Broadcast Tests ---\n"); mu_run_test(test_constant_broadcast_vector_mult, tests_run); + mu_run_test(test_constant_promote_vector_mult, tests_run); + mu_run_test(test_param_broadcast_vector_mult, tests_run); + mu_run_test(test_param_promote_vector_mult, tests_run); #endif /* PROFILE_ONLY */ #ifdef PROFILE_ONLY diff --git a/tests/problem/test_param_broadcast.h b/tests/problem/test_param_broadcast.h index 3532c88..4457160 100644 --- a/tests/problem/test_param_broadcast.h +++ b/tests/problem/test_param_broadcast.h @@ -70,4 +70,88 @@ const char *test_constant_promote_vector_mult(void) free_problem(prob); return 0; } + +const char *test_param_broadcast_vector_mult(void) +{ + int n = 6; + + /* minimize sum(x) subject to broadcast(p) ∘ x, with p parameter */ + expr *x = new_variable(2, 3, 0, n); + expr *objective = new_sum(x, -1); + double c_vals[3] = {1.0, 2.0, 3.0}; + expr *c = new_parameter(1, 3, 0, n, c_vals); + expr *c_bcast = new_broadcast(c, 2, 3); + expr *constraint = new_vector_mult(c_bcast, x); + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, false); + + expr *param_nodes[1] = {c}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* point for evaluating */ + double x_vals[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + + problem_constraint_forward(prob, x_vals); + double constrs[6] = {1.0, 2.0, 6.0, 8.0, 15.0, 18.0}; + problem_jacobian(prob); + double jac_x[6] = {1.0, 1.0, 2.0, 2.0, 3.0, 3.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + + /* second iteration after updating parameter */ + double theta[3] = {5.0, 4.0, 3.0}; + problem_update_params(prob, theta); + problem_constraint_forward(prob, x_vals); + problem_jacobian(prob); + double updated_constrs[6] = {5.0, 10.0, 12.0, 16.0, 15.0, 18.0}; + double updated_jac_x[6] = {5.0, 5.0, 4.0, 4.0, 3.0, 3.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, updated_constrs, 6)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, updated_jac_x, 6)); + + free_problem(prob); + return 0; +} + +const char *test_param_promote_vector_mult(void) +{ + int n = 6; + + /* minimize sum(x) subject to promote(p) ∘ x, with p parameter */ + expr *x = new_variable(2, 3, 0, n); + expr *objective = new_sum(x, -1); + double c_vals = 3.0; + expr *c = new_parameter(1, 1, 0, n, &c_vals); + expr *c_bcast = new_promote(c, 2, 3); + expr *constraint = new_vector_mult(c_bcast, x); + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, false); + + expr *param_nodes[1] = {c}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* point for evaluating */ + double x_vals[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + + problem_constraint_forward(prob, x_vals); + double constrs[6] = {3.0, 6.0, 9.0, 12.0, 15.0, 18.0}; + problem_jacobian(prob); + double jac_x[6] = {3.0, 3.0, 3.0, 3.0, 3.0, 3.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + + /* second iteration after updating parameter */ + double theta = 5.0; + problem_update_params(prob, &theta); + problem_constraint_forward(prob, x_vals); + problem_jacobian(prob); + double updated_constrs[6] = {5.0, 10.0, 15.0, 20.0, 25.0, 30.0}; + double updated_jac_x[6] = {5.0, 5.0, 5.0, 5.0, 5.0, 5.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, updated_constrs, 6)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, updated_jac_x, 6)); + + free_problem(prob); + return 0; +} #endif /* TEST_PARAM_BROADCAST_H */ From d0ed832207030bfe96ef77e3c555d95bbfb6f37d Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Thu, 9 Apr 2026 14:10:57 -0400 Subject: [PATCH 06/12] cleanup left matmul as well --- src/atoms/affine/left_matmul.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/atoms/affine/left_matmul.c b/src/atoms/affine/left_matmul.c index 8da9e3c..b6f1788 100644 --- a/src/atoms/affine/left_matmul.c +++ b/src/atoms/affine/left_matmul.c @@ -64,10 +64,11 @@ static void forward(expr *node, const double *u) { left_matmul_expr *lnode = (left_matmul_expr *) node; - /* Refresh param_source expression tree if parameters changed.*/ - if (lnode->param_source != NULL && lnode->base.needs_parameter_refresh) + /* Always call forward on param_source if it exists */ + /* Should we also adopt a convention that left_matmul always + points to a param_source, even if its constant? */ + if (lnode->param_source != NULL) { - /* pass NULL to forward: constant param_source never depends on u */ lnode->param_source->forward(lnode->param_source, NULL); } From 860dd5a3113a8c2bf1ad2249dc0ae649a4375023 Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Thu, 9 Apr 2026 14:14:36 -0400 Subject: [PATCH 07/12] some very minor cleanups --- src/atoms/affine/parameter.c | 2 +- src/atoms/affine/scalar_mult.c | 1 - src/atoms/affine/vector_mult.c | 1 - 3 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/atoms/affine/parameter.c b/src/atoms/affine/parameter.c index a951aae..e83f87b 100644 --- a/src/atoms/affine/parameter.c +++ b/src/atoms/affine/parameter.c @@ -66,7 +66,7 @@ expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *valu is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL); pnode->param_id = param_id; - + if (values == NULL) { fprintf(stderr, "Parameter values should always be set, this is a bug and" diff --git a/src/atoms/affine/scalar_mult.c b/src/atoms/affine/scalar_mult.c index ef3f6bc..bb1c80e 100644 --- a/src/atoms/affine/scalar_mult.c +++ b/src/atoms/affine/scalar_mult.c @@ -124,7 +124,6 @@ expr *new_scalar_mult(expr *param_node, expr *child) mult_node->param_source = param_node; expr_retain(param_node); - //node->needs_parameter_refresh = true; return node; } diff --git a/src/atoms/affine/vector_mult.c b/src/atoms/affine/vector_mult.c index ed40795..5ed4242 100644 --- a/src/atoms/affine/vector_mult.c +++ b/src/atoms/affine/vector_mult.c @@ -134,7 +134,6 @@ expr *new_vector_mult(expr *param_node, expr *child) vnode->param_source = param_node; expr_retain(param_node); - //node->needs_parameter_refresh = true; return node; } From 36b8ed9c0aae607980b5be92180a7b1219c60d2a Mon Sep 17 00:00:00 2001 From: dance858 Date: Thu, 9 Apr 2026 16:42:05 -0700 Subject: [PATCH 08/12] some error checks and numerical diff to tests --- src/problem.c | 14 ++++ tests/all_tests.c | 1 + .../jacobian_tests/affine/test_left_matmul.h | 2 +- .../composite/test_chain_rule_jacobian.h | 26 +++---- .../composite/test_composite_exp.h | 2 +- tests/numerical_diff.c | 2 +- tests/numerical_diff.h | 2 +- tests/numerical_diff/test_numerical_diff.h | 2 +- tests/problem/test_param_broadcast.h | 25 ++++++- tests/problem/test_param_prob.h | 73 +++++++++++++++++-- 10 files changed, 123 insertions(+), 26 deletions(-) diff --git a/src/problem.c b/src/problem.c index 28448fe..cec796d 100644 --- a/src/problem.c +++ b/src/problem.c @@ -369,6 +369,13 @@ void problem_register_params(problem *prob, expr **param_nodes, int n_param_node prob->total_parameter_size = 0; for (int i = 0; i < n_param_nodes; i++) { + + if (((parameter_expr *) param_nodes[i])->param_id == PARAM_FIXED) + { + fprintf(stderr, "can this ever happen? \n"); + exit(1); + } + // TODO do we need to skip fixed params? maybe we adopt the convention // that we don't ever register fixed params? if (((parameter_expr *) param_nodes[i])->param_id == PARAM_FIXED) continue; @@ -390,6 +397,13 @@ void problem_update_params(problem *prob, const double *theta) { expr *pnode = prob->param_nodes[i]; parameter_expr *param = (parameter_expr *) pnode; + + if (param->param_id == PARAM_FIXED) + { + fprintf(stderr, "can this ever happen? \n"); + exit(1); + } + if (param->param_id == PARAM_FIXED) continue; int offset = param->param_id; memcpy(pnode->value, theta + offset, pnode->size * sizeof(double)); diff --git a/tests/all_tests.c b/tests/all_tests.c index 96a6ff7..09a991c 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -359,6 +359,7 @@ int main(void) mu_run_test(test_param_right_matmul_rectangular, tests_run); mu_run_test(test_param_shared_left_matmul_problem, tests_run); mu_run_test(test_param_fixed_skip_in_update, tests_run); + mu_run_test(test_param_scalar_mult_problem_with_constant, tests_run); printf("\n--- Parameter + Broadcast Tests ---\n"); mu_run_test(test_constant_broadcast_vector_mult, tests_run); diff --git a/tests/jacobian_tests/affine/test_left_matmul.h b/tests/jacobian_tests/affine/test_left_matmul.h index 41cccc4..a0c1385 100644 --- a/tests/jacobian_tests/affine/test_left_matmul.h +++ b/tests/jacobian_tests/affine/test_left_matmul.h @@ -138,7 +138,7 @@ const char *test_jacobian_left_matmul_exp_composite(void) expr *A_exp_Bx = new_left_matmul(NULL, exp_Bx, A); mu_assert("check_jacobian failed", - check_jacobian(A_exp_Bx, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(A_exp_Bx, x_vals, NUMERICAL_DIFF_DEFAULT_H)); free_csr_matrix(A); free_csr_matrix(B); diff --git a/tests/jacobian_tests/composite/test_chain_rule_jacobian.h b/tests/jacobian_tests/composite/test_chain_rule_jacobian.h index 829ac4d..bb7a23c 100644 --- a/tests/jacobian_tests/composite/test_chain_rule_jacobian.h +++ b/tests/jacobian_tests/composite/test_chain_rule_jacobian.h @@ -17,7 +17,7 @@ const char *test_jacobian_exp_sum(void) expr *exp_sum_x = new_exp(sum_x); mu_assert("check_jacobian failed", - check_jacobian(exp_sum_x, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(exp_sum_x, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(exp_sum_x); return 0; @@ -34,7 +34,7 @@ const char *test_jacobian_exp_sum_mult(void) expr *exp_sum_xy = new_exp(sum_xy); mu_assert("check_jacobian failed", - check_jacobian(exp_sum_xy, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(exp_sum_xy, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(exp_sum_xy); return 0; @@ -49,7 +49,7 @@ const char *test_jacobian_sin_cos(void) expr *sin_cos_x = new_sin(cos_x); mu_assert("check_jacobian failed", - check_jacobian(sin_cos_x, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(sin_cos_x, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(sin_cos_x); return 0; @@ -68,7 +68,7 @@ const char *test_jacobian_cos_sin_multiply(void) expr *multiply = new_elementwise_mult(sum, sin_y); mu_assert("check_jacobian failed", - check_jacobian(multiply, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(multiply, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(multiply); return 0; @@ -87,7 +87,7 @@ const char *test_jacobian_Ax_Bx_multiply(void) expr *multiply = new_elementwise_mult(Ax, Bx); mu_assert("check_jacobian failed", - check_jacobian(multiply, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(multiply, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(multiply); free_csr_matrix(A); @@ -107,7 +107,7 @@ const char *test_jacobian_AX_BX_multiply(void) expr *multiply = new_elementwise_mult(new_sin(AX), new_cos(BX)); mu_assert("check_jacobian failed", - check_jacobian(multiply, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(multiply, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(multiply); free_csr_matrix(A); @@ -137,7 +137,7 @@ const char *test_jacobian_quad_form_Ax(void) expr *node = new_quad_form(sin_Ax, Q); mu_assert("check_jacobian failed", - check_jacobian(node, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(node, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(node); free_csr_matrix(A); @@ -164,7 +164,7 @@ const char *test_jacobian_quad_form_exp(void) expr *node = new_quad_form(exp_x, Q); mu_assert("check_jacobian failed", - check_jacobian(node, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(node, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(node); free_csr_matrix(Q); @@ -183,7 +183,7 @@ const char *test_jacobian_matmul_exp_exp(void) expr *Z = new_matmul(exp_X, exp_Y); mu_assert("check_jacobian failed", - check_jacobian(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(Z); return 0; @@ -201,7 +201,7 @@ const char *test_jacobian_matmul_sin_cos(void) expr *Z = new_matmul(sin_X, cos_Y); mu_assert("check_jacobian failed", - check_jacobian(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(Z); return 0; @@ -222,7 +222,7 @@ const char *test_jacobian_matmul_Ax_By(void) expr *Z = new_matmul(AX, BY); /* 3x2 */ mu_assert("check_jacobian failed", - check_jacobian(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(Z); free_csr_matrix(A); @@ -246,7 +246,7 @@ const char *test_jacobian_matmul_sin_Ax_cos_Bx(void) expr *Z = new_matmul(sin_AX, cos_BX); /* 2x2 */ mu_assert("check_jacobian failed", - check_jacobian(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(Z); free_csr_matrix(A); @@ -263,7 +263,7 @@ const char *test_jacobian_matmul_X_X(void) expr *Z = new_matmul(X, X); /* 2x2 */ mu_assert("check_jacobian failed", - check_jacobian(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(Z); return 0; diff --git a/tests/jacobian_tests/composite/test_composite_exp.h b/tests/jacobian_tests/composite/test_composite_exp.h index 725e153..32671d1 100644 --- a/tests/jacobian_tests/composite/test_composite_exp.h +++ b/tests/jacobian_tests/composite/test_composite_exp.h @@ -71,7 +71,7 @@ const char *test_jacobian_composite_exp_add(void) expr *sum = new_add(exp_Ax, exp_By); mu_assert("check_jacobian failed", - check_jacobian(sum, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(sum, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(sum); free_csr_matrix(A); diff --git a/tests/numerical_diff.c b/tests/numerical_diff.c index e018998..9c9812a 100644 --- a/tests/numerical_diff.c +++ b/tests/numerical_diff.c @@ -61,7 +61,7 @@ double *numerical_jacobian(expr *node, const double *u, double h) return J; } -int check_jacobian(expr *node, const double *u, double h) +int check_jacobian_num(expr *node, const double *u, double h) { int m = node->size; int n = node->n_vars; diff --git a/tests/numerical_diff.h b/tests/numerical_diff.h index e303e51..60008fa 100644 --- a/tests/numerical_diff.h +++ b/tests/numerical_diff.h @@ -13,7 +13,7 @@ double *numerical_jacobian(expr *node, const double *u, double h); /* Evaluate analytical Jacobian, compute numerical Jacobian, * and compare. Returns 1 on match, 0 on mismatch. * Prints diagnostic on first failing entry. */ -int check_jacobian(expr *node, const double *u, double h); +int check_jacobian_num(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. diff --git a/tests/numerical_diff/test_numerical_diff.h b/tests/numerical_diff/test_numerical_diff.h index f77e3e9..8af633d 100644 --- a/tests/numerical_diff/test_numerical_diff.h +++ b/tests/numerical_diff/test_numerical_diff.h @@ -22,7 +22,7 @@ const char *test_check_jacobian_composite_exp(void) expr *exp_node = new_exp(Au); mu_assert("check_jacobian failed", - check_jacobian(exp_node, u_vals, NUMERICAL_DIFF_DEFAULT_H)); + check_jacobian_num(exp_node, u_vals, NUMERICAL_DIFF_DEFAULT_H)); free_expr(exp_node); free_csr_matrix(A); diff --git a/tests/problem/test_param_broadcast.h b/tests/problem/test_param_broadcast.h index 4457160..55c18eb 100644 --- a/tests/problem/test_param_broadcast.h +++ b/tests/problem/test_param_broadcast.h @@ -8,6 +8,7 @@ #include "atoms/affine.h" #include "expr.h" #include "minunit.h" +#include "numerical_diff.h" #include "problem.h" #include "subexpr.h" #include "test_helpers.h" @@ -37,6 +38,9 @@ const char *test_constant_broadcast_vector_mult(void) mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + free_problem(prob); return 0; } @@ -67,6 +71,9 @@ const char *test_constant_promote_vector_mult(void) mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + free_problem(prob); return 0; } @@ -99,6 +106,9 @@ const char *test_param_broadcast_vector_mult(void) mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + /* second iteration after updating parameter */ double theta[3] = {5.0, 4.0, 3.0}; problem_update_params(prob, theta); @@ -106,9 +116,13 @@ const char *test_param_broadcast_vector_mult(void) problem_jacobian(prob); double updated_constrs[6] = {5.0, 10.0, 12.0, 16.0, 15.0, 18.0}; double updated_jac_x[6] = {5.0, 5.0, 4.0, 4.0, 3.0, 3.0}; - mu_assert("vals fail", cmp_double_array(prob->constraint_values, updated_constrs, 6)); + mu_assert("vals fail", + cmp_double_array(prob->constraint_values, updated_constrs, 6)); mu_assert("vals fail", cmp_double_array(prob->jacobian->x, updated_jac_x, 6)); + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + free_problem(prob); return 0; } @@ -141,6 +155,9 @@ const char *test_param_promote_vector_mult(void) mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 6)); mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + /* second iteration after updating parameter */ double theta = 5.0; problem_update_params(prob, &theta); @@ -148,9 +165,13 @@ const char *test_param_promote_vector_mult(void) problem_jacobian(prob); double updated_constrs[6] = {5.0, 10.0, 15.0, 20.0, 25.0, 30.0}; double updated_jac_x[6] = {5.0, 5.0, 5.0, 5.0, 5.0, 5.0}; - mu_assert("vals fail", cmp_double_array(prob->constraint_values, updated_constrs, 6)); + mu_assert("vals fail", + cmp_double_array(prob->constraint_values, updated_constrs, 6)); mu_assert("vals fail", cmp_double_array(prob->jacobian->x, updated_jac_x, 6)); + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + free_problem(prob); return 0; } diff --git a/tests/problem/test_param_prob.h b/tests/problem/test_param_prob.h index 59c1c7a..40350d6 100644 --- a/tests/problem/test_param_prob.h +++ b/tests/problem/test_param_prob.h @@ -58,6 +58,53 @@ const char *test_param_scalar_mult_problem(void) return 0; } +const char *test_param_scalar_mult_problem_with_constant(void) +{ + int n = 2; + + /* minimize a * (sum(log(x)) + c), with a parameter, c = 4 constant */ + expr *x = new_variable(2, 1, 0, n); + expr *log_x = new_log(x); + double theta[1] = {3.0}; + double c_val = 4.0; + expr *a_param = new_parameter(1, 1, 0, n, theta); + expr *c = new_parameter(1, 1, PARAM_FIXED, n, &c_val); + expr *scaled = new_scalar_mult(a_param, new_add(new_sum(log_x, -1), c)); + expr *objective = new_sum(scaled, -1); + problem *prob = new_problem(objective, NULL, 0, false); + + /* register parameters and fill sparsity patterns */ + expr *param_nodes[1] = {a_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* point for evaluating */ + double x_vals[2] = {1.0, 2.0}; + + /* test 1: a=3 */ + double obj_val = problem_objective_forward(prob, x_vals); + problem_gradient(prob); + double expected_obj = 3.0 * (log(2.0) + 4.0); + mu_assert("vals fail", fabs(obj_val - expected_obj) < 1e-10); + double grad[2] = {3.0, 1.5}; + mu_assert("vals fail", cmp_double_array(prob->gradient_values, grad, 2)); + + /* test 2: a=5 */ + theta[0] = 5.0; + problem_update_params(prob, theta); + obj_val = problem_objective_forward(prob, x_vals); + problem_gradient(prob); + expected_obj = 5.0 * (log(2.0) + 4.0); + mu_assert("vals fail", fabs(obj_val - expected_obj) < 1e-10); + grad[0] = 5.0; + grad[1] = 2.5; + mu_assert("vals fail", cmp_double_array(prob->gradient_values, grad, 2)); + + free_problem(prob); + + return 0; +} + const char *test_param_vector_mult_problem(void) { int n = 2; @@ -148,7 +195,10 @@ const char *test_param_left_matmul_problem(void) mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 4)); /* test 2: A = [[1,2],[3,4]] (column-major [1,3,2,4]) */ - theta[0] = 1.0; theta[1] = 3.0; theta[2] = 2.0; theta[3] = 4.0; + theta[0] = 1.0; + theta[1] = 3.0; + theta[2] = 2.0; + theta[3] = 4.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); @@ -321,8 +371,12 @@ const char *test_param_left_matmul_rectangular(void) mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 6)); /* test 2: A = [[7,8],[9,10],[11,12]] (column-major [7,9,11,8,10,12]) */ - theta[0] = 7.0; theta[1] = 9.0; theta[2] = 11.0; - theta[3] = 8.0; theta[4] = 10.0; theta[5] = 12.0; + theta[0] = 7.0; + theta[1] = 9.0; + theta[2] = 11.0; + theta[3] = 8.0; + theta[4] = 10.0; + theta[5] = 12.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); @@ -383,8 +437,12 @@ const char *test_param_right_matmul_rectangular(void) mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 6)); /* test 2: A = [[7,8,9],[10,11,12]] (column-major [7,10,8,11,9,12]) */ - theta[0] = 7.0; theta[1] = 10.0; theta[2] = 8.0; - theta[3] = 11.0; theta[4] = 9.0; theta[5] = 12.0; + theta[0] = 7.0; + theta[1] = 10.0; + theta[2] = 8.0; + theta[3] = 11.0; + theta[4] = 9.0; + theta[5] = 12.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); @@ -447,7 +505,10 @@ const char *test_param_shared_left_matmul_problem(void) mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 8)); /* test 2: A = [[1,2],[3,4]] (column-major [1,3,2,4]) */ - theta[0] = 1.0; theta[1] = 3.0; theta[2] = 2.0; theta[3] = 4.0; + theta[0] = 1.0; + theta[1] = 3.0; + theta[2] = 2.0; + theta[3] = 4.0; problem_update_params(prob, theta); problem_constraint_forward(prob, x_vals); problem_jacobian(prob); From 5fc185f10bd22b57c3a6c45b8658edbd1b25bc51 Mon Sep 17 00:00:00 2001 From: dance858 Date: Thu, 9 Apr 2026 18:57:56 -0700 Subject: [PATCH 09/12] we don't always have to run forward of parameter in left matmul --- src/atoms/affine/left_matmul.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atoms/affine/left_matmul.c b/src/atoms/affine/left_matmul.c index b6f1788..e56561d 100644 --- a/src/atoms/affine/left_matmul.c +++ b/src/atoms/affine/left_matmul.c @@ -67,7 +67,7 @@ static void forward(expr *node, const double *u) /* Always call forward on param_source if it exists */ /* Should we also adopt a convention that left_matmul always points to a param_source, even if its constant? */ - if (lnode->param_source != NULL) + if (lnode->param_source != NULL && lnode->base.needs_parameter_refresh) { lnode->param_source->forward(lnode->param_source, NULL); } From 73679141551dfedbe5ab33aff652faab86c7aa73 Mon Sep 17 00:00:00 2001 From: dance858 Date: Thu, 9 Apr 2026 19:10:18 -0700 Subject: [PATCH 10/12] we don't always have to call forward for parameter node in vector mult --- src/atoms/affine/left_matmul.c | 4 +--- src/atoms/affine/vector_mult.c | 14 +++++++++++--- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/atoms/affine/left_matmul.c b/src/atoms/affine/left_matmul.c index e56561d..fdda2ae 100644 --- a/src/atoms/affine/left_matmul.c +++ b/src/atoms/affine/left_matmul.c @@ -64,9 +64,7 @@ static void forward(expr *node, const double *u) { left_matmul_expr *lnode = (left_matmul_expr *) node; - /* Always call forward on param_source if it exists */ - /* Should we also adopt a convention that left_matmul always - points to a param_source, even if its constant? */ + /* 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); diff --git a/src/atoms/affine/vector_mult.c b/src/atoms/affine/vector_mult.c index 5ed4242..b1663e6 100644 --- a/src/atoms/affine/vector_mult.c +++ b/src/atoms/affine/vector_mult.c @@ -30,9 +30,14 @@ static void forward(expr *node, const double *u) expr *child = node->left; vector_mult_expr *vnode = (vector_mult_expr *) node; - /* call forward for param_source expr tree - ex: broadcast(param) or promote(const)*/ - vnode->param_source->forward(vnode->param_source, NULL); + /* call forward for param_source expr tree (this extra logic is needed + in case the parameter is a broadcast or promote node which needs to refresh + its values) */ + if (vnode->base.needs_parameter_refresh) + { + vnode->param_source->forward(vnode->param_source, NULL); + vnode->base.needs_parameter_refresh = false; + } const double *a = vnode->param_source->value; @@ -135,5 +140,8 @@ expr *new_vector_mult(expr *param_node, expr *child) vnode->param_source = param_node; expr_retain(param_node); + /* special case for handling broadcasting of constants correctly */ + vnode->base.needs_parameter_refresh = true; + return node; } From 0e7a58afc906ef11d3d33eae7522ee28efe63358 Mon Sep 17 00:00:00 2001 From: dance858 Date: Thu, 9 Apr 2026 19:11:36 -0700 Subject: [PATCH 11/12] comment out forward parameter pass in scalar mult because it is not needed, I think --- src/atoms/affine/scalar_mult.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atoms/affine/scalar_mult.c b/src/atoms/affine/scalar_mult.c index bb1c80e..5b4cd67 100644 --- a/src/atoms/affine/scalar_mult.c +++ b/src/atoms/affine/scalar_mult.c @@ -32,7 +32,7 @@ static void forward(expr *node, const double *u) /* call forward for param_source expr tree ex: broadcast(param) or promote(const)*/ - snode->param_source->forward(snode->param_source, NULL); + // snode->param_source->forward(snode->param_source, NULL); double a = snode->param_source->value[0]; From 45776043dc67340d819c8f30f3a017db1a828b77 Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Fri, 10 Apr 2026 15:28:01 -0400 Subject: [PATCH 12/12] add test for scalar case to be consistent with vector mult --- src/atoms/affine/parameter.c | 3 +- src/atoms/affine/scalar_mult.c | 14 ++++- tests/all_tests.c | 2 + tests/problem/test_param_broadcast.h | 82 ++++++++++++++++++++++++++++ 4 files changed, 97 insertions(+), 4 deletions(-) diff --git a/src/atoms/affine/parameter.c b/src/atoms/affine/parameter.c index e83f87b..49bade7 100644 --- a/src/atoms/affine/parameter.c +++ b/src/atoms/affine/parameter.c @@ -64,7 +64,8 @@ expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *valu 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); - + + // TODO we should assert that the values array has the correct size. pnode->param_id = param_id; if (values == NULL) diff --git a/src/atoms/affine/scalar_mult.c b/src/atoms/affine/scalar_mult.c index 5b4cd67..b10d17f 100644 --- a/src/atoms/affine/scalar_mult.c +++ b/src/atoms/affine/scalar_mult.c @@ -30,9 +30,14 @@ static void forward(expr *node, const double *u) expr *child = node->left; scalar_mult_expr *snode = (scalar_mult_expr *) node; - /* call forward for param_source expr tree - ex: broadcast(param) or promote(const)*/ - // snode->param_source->forward(snode->param_source, NULL); + /* call forward for param_source expr tree (this extra logic is needed + in case the parameter is a broadcast or promote node which needs to refresh + its values) */ + if (snode->base.needs_parameter_refresh) + { + snode->param_source->forward(snode->param_source, NULL); + snode->base.needs_parameter_refresh = false; + } double a = snode->param_source->value[0]; @@ -125,5 +130,8 @@ expr *new_scalar_mult(expr *param_node, expr *child) mult_node->param_source = param_node; expr_retain(param_node); + /* special case for handling broadcasting of constants correctly */ + mult_node->base.needs_parameter_refresh = true; + return node; } diff --git a/tests/all_tests.c b/tests/all_tests.c index 09a991c..6eeab3c 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -366,6 +366,8 @@ int main(void) mu_run_test(test_constant_promote_vector_mult, tests_run); mu_run_test(test_param_broadcast_vector_mult, tests_run); mu_run_test(test_param_promote_vector_mult, tests_run); + mu_run_test(test_const_sum_scalar_mult, tests_run); + mu_run_test(test_param_sum_scalar_mult, tests_run); #endif /* PROFILE_ONLY */ #ifdef PROFILE_ONLY diff --git a/tests/problem/test_param_broadcast.h b/tests/problem/test_param_broadcast.h index 55c18eb..a985a36 100644 --- a/tests/problem/test_param_broadcast.h +++ b/tests/problem/test_param_broadcast.h @@ -175,4 +175,86 @@ const char *test_param_promote_vector_mult(void) free_problem(prob); return 0; } + +const char *test_const_sum_scalar_mult(void) +{ + int n = 6; + + /* minimize sum(x) subject to sum(c) * x, with c constant */ + expr *x = new_variable(1, 1, 0, n); + expr *objective = new_sum(x, -1); + double c_vals[3] = {1.0, 2.0, 3.0}; + expr *c = new_parameter(1, 3, PARAM_FIXED, n, c_vals); + expr *c_sum = new_sum(c, -1); + expr *constraint = new_scalar_mult(c_sum, x); + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, false); + + problem_init_derivatives(prob); + + /* point for evaluating */ + double x_vals[1] = {4.0}; + + problem_constraint_forward(prob, x_vals); + double constrs[1] = {6.0 * 4.0}; + problem_jacobian(prob); + double jac_x[1] = {6.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 1)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 1)); + + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + free_problem(prob); + return 0; +} + +const char *test_param_sum_scalar_mult(void) +{ + int n = 6; + + /* minimize sum(x) subject to sum(p) * x, with p parameter */ + expr *x = new_variable(1, 1, 0, n); + expr *objective = new_sum(x, -1); + double c_vals[3] = {1.0, 2.0, 3.0}; + expr *c = new_parameter(1, 3, 0, n, c_vals); + expr *c_sum = new_sum(c, -1); + expr *constraint = new_scalar_mult(c_sum, x); + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, false); + + expr *param_nodes[1] = {c}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* point for evaluating */ + double x_vals[1] = {4.0}; + + problem_constraint_forward(prob, x_vals); + double constrs[1] = {6.0 * 4.0}; + problem_jacobian(prob); + double jac_x[1] = {6.0}; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 1)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, jac_x, 1)); + + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + + /* second iteration after updating parameter */ + double theta[3] = {5.0, 4.0, 3.0}; + problem_update_params(prob, theta); + problem_constraint_forward(prob, x_vals); + problem_jacobian(prob); + double updated_constrs[1] = {12.0 * 4.0}; + double updated_jac_x[1] = {12.0}; + mu_assert("vals fail", + cmp_double_array(prob->constraint_values, updated_constrs, 1)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, updated_jac_x, 1)); + + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + + free_problem(prob); + return 0; +} + #endif /* TEST_PARAM_BROADCAST_H */