diff --git a/src/atoms/affine/left_matmul.c b/src/atoms/affine/left_matmul.c index 9ff397e..fdda2ae 100644 --- a/src/atoms/affine/left_matmul.c +++ b/src/atoms/affine/left_matmul.c @@ -63,6 +63,13 @@ 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; + + /* 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; diff --git a/src/atoms/affine/parameter.c b/src/atoms/affine/parameter.c index fa96377..49bade7 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 @@ -63,13 +64,16 @@ 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) + 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 baaf16e..b10d17f 100644 --- a/src/atoms/affine/scalar_mult.c +++ b/src/atoms/affine/scalar_mult.c @@ -28,7 +28,18 @@ 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; + + /* 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]; /* child's forward pass */ child->forward(child, u); @@ -119,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/src/atoms/affine/vector_mult.c b/src/atoms/affine/vector_mult.c index 116f012..b1663e6 100644 --- a/src/atoms/affine/vector_mult.c +++ b/src/atoms/affine/vector_mult.c @@ -28,7 +28,18 @@ 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; + + /* 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; /* child's forward pass */ child->forward(child, u); @@ -129,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; } diff --git a/src/problem.c b/src/problem.c index 0419063..cec796d 100644 --- a/src/problem.c +++ b/src/problem.c @@ -369,6 +369,16 @@ 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; prob->total_parameter_size += param_nodes[i]->size; } } @@ -387,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 3ad1b7c..6eeab3c 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,15 @@ 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); + 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/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 new file mode 100644 index 0000000..a985a36 --- /dev/null +++ b/tests/problem/test_param_broadcast.h @@ -0,0 +1,260 @@ +#ifndef TEST_PARAM_BROADCAST_H +#define TEST_PARAM_BROADCAST_H + +#include +#include +#include + +#include "atoms/affine.h" +#include "expr.h" +#include "minunit.h" +#include "numerical_diff.h" +#include "problem.h" +#include "subexpr.h" +#include "test_helpers.h" + +const char *test_constant_broadcast_vector_mult(void) +{ + int n = 6; + + /* 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); + 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); + 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)); + + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + + free_problem(prob); + return 0; +} + +const char *test_constant_promote_vector_mult(void) +{ + int n = 6; + + /* 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); + 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); + + 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)); + + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + + 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)); + + 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[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)); + + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + + 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)); + + 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); + 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)); + + mu_assert("check_jacobian failed", + check_jacobian_num(constraint, x_vals, NUMERICAL_DIFF_DEFAULT_H)); + + 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 */ diff --git a/tests/problem/test_param_prob.h b/tests/problem/test_param_prob.h index f5141e5..40350d6 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); @@ -59,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; @@ -66,7 +112,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 +130,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 +166,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 +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]) */ - 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 +226,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 +255,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 +291,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 +300,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 +342,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 +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]) */ - 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 +407,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 +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]) */ - 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 +474,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 +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]) */ - 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);