diff --git a/include/atoms/affine.h b/include/atoms/affine.h index c51bc64..3fa076d 100644 --- a/include/atoms/affine.h +++ b/include/atoms/affine.h @@ -31,7 +31,7 @@ expr *new_vstack(expr **args, int n_args, int n_vars); expr *new_promote(expr *child, int d1, int d2); expr *new_trace(expr *child); -expr *new_constant(int d1, int d2, int n_vars, const double *values); +expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *values); expr *new_variable(int d1, int d2, int var_id, int n_vars); expr *new_index(expr *child, int d1, int d2, const int *indices, int n_idxs); @@ -40,26 +40,27 @@ expr *new_broadcast(expr *child, int target_d1, int target_d2); expr *new_diag_vec(expr *child); expr *new_transpose(expr *child); -/* Left matrix multiplication: A @ f(x) where A is a constant sparse - * matrix */ -expr *new_left_matmul(expr *u, const CSR_Matrix *A); +/* Left matrix multiplication: A @ f(x) where A is a constant or parameter + * sparse matrix. param_node is NULL for fixed constants. */ +expr *new_left_matmul(expr *param_node, expr *u, const CSR_Matrix *A); -/* Left matrix multiplication: A @ f(x) where A is a constant dense - * matrix (row-major, m x n). Uses CBLAS for efficient computation. */ -expr *new_left_matmul_dense(expr *u, int m, int n, const double *data); +/* Left matrix multiplication: A @ f(x) where A is a constant or parameter + * dense matrix (row-major, m x n). Uses CBLAS for efficient computation. */ +expr *new_left_matmul_dense(expr *param_node, expr *u, int m, int n, + const double *data); -/* Right matrix multiplication: f(x) @ A where A is a constant - * matrix */ -expr *new_right_matmul(expr *u, const CSR_Matrix *A); +/* Right matrix multiplication: f(x) @ A where A is a constant or parameter + * matrix. */ +expr *new_right_matmul(expr *param_node, expr *u, const CSR_Matrix *A); -expr *new_right_matmul_dense(expr *u, int m, int n, const double *data); +expr *new_right_matmul_dense(expr *param_node, expr *u, int m, int n, + const double *data); -/* Constant scalar multiplication: a * f(x) where a is a constant - * double */ -expr *new_const_scalar_mult(double a, expr *child); +/* Scalar multiplication: a * f(x) where a comes from param_node */ +expr *new_scalar_mult(expr *param_node, expr *child); -/* Constant vector elementwise multiplication: a . f(x) where a is - * constant */ -expr *new_const_vector_mult(const double *a, expr *child); +/* Vector elementwise multiplication: a . f(x) where a comes from + * param_node */ +expr *new_vector_mult(expr *param_node, expr *child); #endif /* AFFINE_H */ diff --git a/include/problem.h b/include/problem.h index 459a6a4..f5775ad 100644 --- a/include/problem.h +++ b/include/problem.h @@ -49,6 +49,11 @@ typedef struct problem int n_vars; int total_constraint_size; + /* parameter support */ + expr **param_nodes; + int n_param_nodes; + int total_parameter_size; + /* allocated by new_problem */ double *constraint_values; double *gradient_values; @@ -79,6 +84,9 @@ void problem_init_jacobian_coo(problem *prob); void problem_init_hessian_coo_lower_triangular(problem *prob); void free_problem(problem *prob); +void problem_register_params(problem *prob, expr **param_nodes, int n_param_nodes); +void problem_update_params(problem *prob, const double *theta); + double problem_objective_forward(problem *prob, const double *u); void problem_constraint_forward(problem *prob, const double *u); void problem_gradient(problem *prob); diff --git a/include/subexpr.h b/include/subexpr.h index 44a964b..2b638da 100644 --- a/include/subexpr.h +++ b/include/subexpr.h @@ -26,8 +26,22 @@ /* Forward declaration */ struct int_double_pair; +/* Parameter ID for fixed constants (not updatable) */ +#define PARAM_FIXED -1 + /* Type-specific expression structures that "inherit" from expr */ +/* Unified constant/parameter node. Constants use param_id == PARAM_FIXED. + * Updatable parameters use param_id >= 0 (offset into global theta). */ +typedef struct parameter_expr +{ + expr base; + int param_id; + /* Set to true by problem_update_params(), cleared by + refresh_param_values() after propagating new values. */ + bool needs_refresh; +} parameter_expr; + /* Linear operator: y = A * x + b * The matrix A is stored as node->jacobian (CSR). */ typedef struct linear_op_expr @@ -118,18 +132,24 @@ typedef struct left_matmul_expr CSC_Matrix *Jchild_CSC; CSC_Matrix *J_CSC; int *csc_to_csr_work; + expr *param_source; + void (*refresh_param_values)(struct left_matmul_expr *); } left_matmul_expr; -/* Right matrix multiplication: y = f(x) * A where f(x) is an expression. - * f(x) has shape p x n, A has shape n x q, output y has shape p x q. - * Uses vec(y) = B * vec(f(x)) where B = A^T kron I_p. */ -typedef struct right_matmul_expr +/* Scalar multiplication: y = a * child where a comes from param_source */ +typedef struct scalar_mult_expr +{ + expr base; + expr *param_source; +} scalar_mult_expr; + +/* Vector elementwise multiplication: y = a \circ child where a comes from + * param_source */ +typedef struct vector_mult_expr { expr base; - CSR_Matrix *B; /* B = A^T kron I_p */ - CSR_Matrix *BT; /* B^T for backpropagating Hessian weights */ - CSC_Matrix *CSC_work; -} right_matmul_expr; + expr *param_source; +} vector_mult_expr; /* Bivariate matrix multiplication: Z = f(u) @ g(u) where both children * may be composite expressions. */ @@ -153,20 +173,6 @@ typedef struct matmul_expr int *idx_map_Hg; } matmul_expr; -/* Constant scalar multiplication: y = a * child where a is a constant double */ -typedef struct const_scalar_mult_expr -{ - expr base; - double a; -} const_scalar_mult_expr; - -/* Constant vector elementwise multiplication: y = a \circ child for constant a */ -typedef struct const_vector_mult_expr -{ - expr base; - double *a; /* length equals node->size */ -} const_vector_mult_expr; - /* Index/slicing: y = child[indices] where indices is a list of flat positions */ typedef struct index_expr { diff --git a/include/utils/dense_matrix.h b/include/utils/dense_matrix.h index 8121020..edcea66 100644 --- a/include/utils/dense_matrix.h +++ b/include/utils/dense_matrix.h @@ -17,4 +17,6 @@ Matrix *new_dense_matrix(int m, int n, const double *data); /* Transpose helper */ Matrix *dense_matrix_trans(const Dense_Matrix *self); +void A_transpose(double *AT, const double *A, int m, int n); + #endif /* DENSE_MATRIX_H */ diff --git a/include/utils/matrix.h b/include/utils/matrix.h index 478fabd..96538aa 100644 --- a/include/utils/matrix.h +++ b/include/utils/matrix.h @@ -31,6 +31,7 @@ typedef struct Matrix const CSC_Matrix *J, int p); void (*block_left_mult_values)(const struct Matrix *self, const CSC_Matrix *J, CSC_Matrix *C); + void (*update_values)(struct Matrix *self, const double *new_values); void (*free_fn)(struct Matrix *self); } Matrix; diff --git a/src/atoms/affine/left_matmul.c b/src/atoms/affine/left_matmul.c index c8f8e19..bce4201 100644 --- a/src/atoms/affine/left_matmul.c +++ b/src/atoms/affine/left_matmul.c @@ -49,16 +49,34 @@ #include "utils/tracked_alloc.h" #include "utils/utils.h" +static void refresh_param_values(left_matmul_expr *lnode) +{ + if (lnode->param_source == NULL) + { + return; + } + parameter_expr *param = (parameter_expr *) lnode->param_source; + if (!param->needs_refresh) + { + return; + } + param->needs_refresh = false; + lnode->refresh_param_values(lnode); +} + static void forward(expr *node, const double *u) { + left_matmul_expr *lnode = (left_matmul_expr *) node; + refresh_param_values(lnode); + expr *x = node->left; /* child's forward pass */ node->left->forward(node->left, u); /* y = A_kron @ vec(f(x)) */ - Matrix *A = ((left_matmul_expr *) node)->A; - int n_blocks = ((left_matmul_expr *) node)->n_blocks; + Matrix *A = lnode->A; + int n_blocks = lnode->n_blocks; A->block_left_mult_vec(A, x->value, node->value, n_blocks); } @@ -75,11 +93,16 @@ static void free_type_data(expr *node) free_csc_matrix(lnode->Jchild_CSC); free_csc_matrix(lnode->J_CSC); free(lnode->csc_to_csr_work); + if (lnode->param_source != NULL) + { + free_expr(lnode->param_source); + } lnode->A = NULL; lnode->AT = NULL; lnode->Jchild_CSC = NULL; lnode->J_CSC = NULL; lnode->csc_to_csr_work = NULL; + lnode->param_source = NULL; } static void jacobian_init_impl(expr *node) @@ -99,8 +122,8 @@ static void jacobian_init_impl(expr *node) static void eval_jacobian(expr *node) { - expr *x = node->left; left_matmul_expr *lnode = (left_matmul_expr *) node; + expr *x = node->left; CSC_Matrix *Jchild_CSC = lnode->Jchild_CSC; CSC_Matrix *J_CSC = lnode->J_CSC; @@ -131,9 +154,11 @@ static void wsum_hess_init_impl(expr *node) static void eval_wsum_hess(expr *node, const double *w) { + left_matmul_expr *lnode = (left_matmul_expr *) node; + /* compute A^T w*/ - Matrix *AT = ((left_matmul_expr *) node)->AT; - int n_blocks = ((left_matmul_expr *) node)->n_blocks; + Matrix *AT = lnode->AT; + int n_blocks = lnode->n_blocks; AT->block_left_mult_vec(AT, w, node->work->dwork, n_blocks); node->left->eval_wsum_hess(node->left, node->work->dwork); @@ -141,7 +166,30 @@ static void eval_wsum_hess(expr *node, const double *w) node->wsum_hess->nnz * sizeof(double)); } -expr *new_left_matmul(expr *u, const CSR_Matrix *A) +static void refresh_sparse_left(left_matmul_expr *lnode) +{ + (void) lnode; + fprintf(stderr, + "Error in refresh_sparse_left: parameter for a sparse matrix not " + "supported \n"); + exit(1); +} + +static void refresh_dense_left(left_matmul_expr *lnode) +{ + Dense_Matrix *dm_A = (Dense_Matrix *) lnode->A; + Dense_Matrix *dm_AT = (Dense_Matrix *) lnode->AT; + int m = dm_A->base.m; + int n = dm_A->base.n; + + /* The parameter represents the A in left_matmul_dense(A, x) in column-major. + In this diffengine, we store A in row-major order. Hence, param->vals + actually corresponds to the transpose of A, and we transpose AT to get A. */ + memcpy(dm_AT->x, lnode->param_source->value, m * n * sizeof(double)); + A_transpose(dm_A->x, dm_AT->x, n, m); +} + +expr *new_left_matmul(expr *param_node, expr *u, const CSR_Matrix *A) { /* We expect u->d1 == A->n. However, numpy's broadcasting rules allow users to do A @ u where u is (n, ) which in C is actually (1, n). In that case @@ -188,10 +236,24 @@ expr *new_left_matmul(expr *u, const CSR_Matrix *A) lnode->AT = sparse_matrix_trans((const Sparse_Matrix *) lnode->A, node->work->iwork); + /* parameter support */ + lnode->param_source = param_node; + if (param_node != NULL) + { + + fprintf(stderr, "Error in new_left_matmul: parameter for a sparse matrix " + "not supported \n"); + exit(1); + + expr_retain(param_node); + lnode->refresh_param_values = refresh_sparse_left; + } + return node; } -expr *new_left_matmul_dense(expr *u, int m, int n, const double *data) +expr *new_left_matmul_dense(expr *param_node, expr *u, int m, int n, + const double *data) { int d1, d2, n_blocks; if (u->d1 == n) @@ -227,5 +289,13 @@ expr *new_left_matmul_dense(expr *u, int m, int n, const double *data) lnode->A = new_dense_matrix(m, n, data); lnode->AT = dense_matrix_trans((const Dense_Matrix *) lnode->A); + /* parameter support */ + lnode->param_source = param_node; + if (param_node != NULL) + { + expr_retain(param_node); + lnode->refresh_param_values = refresh_dense_left; + } + return node; } diff --git a/src/atoms/affine/constant.c b/src/atoms/affine/parameter.c similarity index 71% rename from src/atoms/affine/constant.c rename to src/atoms/affine/parameter.c index 607b019..7cd2909 100644 --- a/src/atoms/affine/constant.c +++ b/src/atoms/affine/parameter.c @@ -16,39 +16,37 @@ * limitations under the License. */ #include "atoms/affine.h" +#include "subexpr.h" #include "utils/tracked_alloc.h" #include #include static void forward(expr *node, const double *u) { - /* Constants don't depend on u; values are already set */ + /* Parameters/constants don't depend on u; values are already set */ (void) node; (void) u; } static void jacobian_init_impl(expr *node) { - /* Constant jacobian is all zeros: size x n_vars with 0 nonzeros. - * new_csr_matrix uses calloc for row pointers, so they're already 0. */ + /* Zero jacobian: size x n_vars with 0 nonzeros. */ node->jacobian = new_csr_matrix(node->size, node->n_vars, 0); } static void eval_jacobian(expr *node) { - /* Constant jacobian never changes - nothing to evaluate */ (void) node; } static void wsum_hess_init_impl(expr *node) { - /* Constant Hessian is all zeros: n_vars x n_vars with 0 nonzeros. */ + /* Zero Hessian: n_vars x n_vars with 0 nonzeros. */ node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, 0); } static void eval_wsum_hess(expr *node, const double *w) { - /* Constant Hessian is always zero - nothing to compute */ (void) node; (void) w; } @@ -59,12 +57,20 @@ static bool is_affine(const expr *node) return true; } -expr *new_constant(int d1, int d2, int n_vars, const double *values) +expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *values) { - expr *node = (expr *) SP_CALLOC(1, sizeof(expr)); + parameter_expr *pnode = (parameter_expr *) SP_CALLOC(1, sizeof(parameter_expr)); + expr *node = &pnode->base; init_expr(node, d1, d2, n_vars, forward, jacobian_init_impl, eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL); - memcpy(node->value, values, node->size * sizeof(double)); + + pnode->param_id = param_id; + pnode->needs_refresh = false; + + if (values != NULL) + { + memcpy(node->value, values, node->size * sizeof(double)); + } return node; } diff --git a/src/atoms/affine/right_matmul.c b/src/atoms/affine/right_matmul.c index b49a1a9..d0f9082 100644 --- a/src/atoms/affine/right_matmul.c +++ b/src/atoms/affine/right_matmul.c @@ -16,16 +16,53 @@ * limitations under the License. */ #include "atoms/affine.h" - +#include "subexpr.h" #include "utils/CSR_Matrix.h" +#include "utils/dense_matrix.h" #include "utils/tracked_alloc.h" +#include #include /* This file implements the atom 'right_matmul' corresponding to the operation y = f(x) @ A, where A is a given matrix and f(x) is an arbitrary expression. We implement this by expressing right matmul in terms of left matmul and - transpose: f(x) @ A = (A^T @ f(x)^T)^T. */ -expr *new_right_matmul(expr *u, const CSR_Matrix *A) + transpose: f(x) @ A = (A^T @ f(x)^T)^T. + + For the parameter case: + - param_source stores A values in CSR data order + - inner left_matmul stores AT as its A-matrix and A as its AT-matrix + - on refresh: update AT (inner's AT, the original A) from param_source, + then recompute A^T (inner's A) from the updated A. */ + +/* Refresh for sparse right_matmul: param stores A in CSR data order. + Inner left_matmul: lnode->A = AT (transposed), lnode->AT = A (original). + So: update lnode->AT from param values, then recompute lnode->A. */ +static void refresh_sparse_right(left_matmul_expr *lnode) +{ + (void) lnode; + fprintf(stderr, + "Error in refresh_sparse_right: parameter for a sparse matrix not " + "supported \n"); + exit(1); +} + +static void refresh_dense_right(left_matmul_expr *lnode) +{ + /* This left_matmul_expr node corresponds to left multiplication with B = AT, + where A is the original (m x n) matrix given to the right_matmul function. + Furthermore, lnode->param_source->value corresponds to the column-major + version of A, which is BT (an m x n matrix) */ + + Dense_Matrix *B = (Dense_Matrix *) lnode->AT; + Dense_Matrix *BT = (Dense_Matrix *) lnode->A; + int m = B->base.n; + int n = B->base.m; + + memcpy(BT->x, lnode->param_source->value, m * n * sizeof(double)); + A_transpose(B->x, BT->x, m, n); +} + +expr *new_right_matmul(expr *param_node, expr *u, const CSR_Matrix *A) { /* We can express right matmul using left matmul and transpose: u @ A = (A^T @ u^T)^T. */ @@ -33,7 +70,23 @@ expr *new_right_matmul(expr *u, const CSR_Matrix *A) CSR_Matrix *AT = transpose(A, work_transpose); expr *u_transpose = new_transpose(u); - expr *left_matmul = new_left_matmul(u_transpose, AT); + expr *left_matmul = new_left_matmul(NULL, u_transpose, AT); + + /* If parameterized, attach param_source and custom refresh to inner + left_matmul */ + if (param_node != NULL) + { + + fprintf(stderr, "Error in new_right_matmul: parameter for a sparse matrix " + "not supported \n"); + exit(1); + + left_matmul_expr *lnode = (left_matmul_expr *) left_matmul; + lnode->param_source = param_node; + expr_retain(param_node); + lnode->refresh_param_values = refresh_sparse_right; + } + expr *node = new_transpose(left_matmul); free_csr_matrix(AT); @@ -41,21 +94,25 @@ expr *new_right_matmul(expr *u, const CSR_Matrix *A) return node; } -expr *new_right_matmul_dense(expr *u, int m, int n, const double *data) +expr *new_right_matmul_dense(expr *param_node, expr *u, int m, int n, + const double *data) { - /* We express: u @ A = (A^T @ u^T)^T - A is m x n, so A^T is n x m. */ + /* We express: u @ A = (A^T @ u^T)^T. A is m x n, so A^T is n x m. */ double *AT = (double *) SP_MALLOC(n * m * sizeof(double)); - for (int i = 0; i < m; i++) + A_transpose(AT, data, m, n); + + expr *u_transpose = new_transpose(u); + expr *left_matmul_node = new_left_matmul_dense(NULL, u_transpose, n, m, AT); + + /* If parameterized, attach param_source and custom refresh */ + if (param_node != NULL) { - for (int j = 0; j < n; j++) - { - AT[j * m + i] = data[i * n + j]; - } + left_matmul_expr *lnode = (left_matmul_expr *) left_matmul_node; + lnode->param_source = param_node; + expr_retain(param_node); + lnode->refresh_param_values = refresh_dense_right; } - expr *u_transpose = new_transpose(u); - expr *left_matmul_node = new_left_matmul_dense(u_transpose, n, m, AT); expr *node = new_transpose(left_matmul_node); free(AT); diff --git a/src/atoms/affine/const_scalar_mult.c b/src/atoms/affine/scalar_mult.c similarity index 76% rename from src/atoms/affine/const_scalar_mult.c rename to src/atoms/affine/scalar_mult.c index f410d4a..baaf16e 100644 --- a/src/atoms/affine/const_scalar_mult.c +++ b/src/atoms/affine/scalar_mult.c @@ -23,17 +23,17 @@ #include #include -/* Constant scalar multiplication: y = a * child where a is a constant double */ +/* Scalar multiplication: y = a * child where a comes from param_source */ static void forward(expr *node, const double *u) { expr *child = node->left; + double a = ((scalar_mult_expr *) node)->param_source->value[0]; /* child's forward pass */ child->forward(child, u); /* local forward pass: multiply each element by scalar a */ - double a = ((const_scalar_mult_expr *) node)->a; for (int i = 0; i < node->size; i++) { node->value[i] = a * child->value[i]; @@ -54,7 +54,7 @@ static void jacobian_init_impl(expr *node) static void eval_jacobian(expr *node) { expr *child = node->left; - double a = ((const_scalar_mult_expr *) node)->a; + double a = ((scalar_mult_expr *) node)->param_source->value[0]; /* evaluate child */ child->eval_jacobian(child); @@ -82,7 +82,7 @@ static void eval_wsum_hess(expr *node, const double *w) expr *x = node->left; x->eval_wsum_hess(x, w); - double a = ((const_scalar_mult_expr *) node)->a; + double a = ((scalar_mult_expr *) node)->param_source->value[0]; for (int j = 0; j < x->wsum_hess->nnz; j++) { node->wsum_hess->x[j] = a * x->wsum_hess->x[j]; @@ -91,21 +91,33 @@ static void eval_wsum_hess(expr *node, const double *w) static bool is_affine(const expr *node) { - /* Affine iff the child is affine */ return node->left->is_affine(node->left); } -expr *new_const_scalar_mult(double a, expr *child) +static void free_type_data(expr *node) { - const_scalar_mult_expr *mult_node = - (const_scalar_mult_expr *) SP_CALLOC(1, sizeof(const_scalar_mult_expr)); + scalar_mult_expr *snode = (scalar_mult_expr *) node; + if (snode->param_source != NULL) + { + free_expr(snode->param_source); + snode->param_source = NULL; + } +} + +expr *new_scalar_mult(expr *param_node, expr *child) +{ + scalar_mult_expr *mult_node = + (scalar_mult_expr *) SP_CALLOC(1, sizeof(scalar_mult_expr)); expr *node = &mult_node->base; init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init_impl, - eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL); + eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, + free_type_data); node->left = child; - mult_node->a = a; expr_retain(child); + mult_node->param_source = param_node; + expr_retain(param_node); + return node; } diff --git a/src/atoms/affine/const_vector_mult.c b/src/atoms/affine/vector_mult.c similarity index 80% rename from src/atoms/affine/const_vector_mult.c rename to src/atoms/affine/vector_mult.c index 8d965a4..116f012 100644 --- a/src/atoms/affine/const_vector_mult.c +++ b/src/atoms/affine/vector_mult.c @@ -22,12 +22,13 @@ #include #include -/* Constant vector elementwise multiplication: y = a \circ child */ +/* Vector elementwise multiplication: y = a \circ child + * where a comes from param_source */ static void forward(expr *node, const double *u) { expr *child = node->left; - const double *a = ((const_vector_mult_expr *) node)->a; + const double *a = ((vector_mult_expr *) node)->param_source->value; /* child's forward pass */ child->forward(child, u); @@ -53,7 +54,7 @@ static void jacobian_init_impl(expr *node) static void eval_jacobian(expr *node) { expr *x = node->left; - const double *a = ((const_vector_mult_expr *) node)->a; + const double *a = ((vector_mult_expr *) node)->param_source->value; /* evaluate x */ x->eval_jacobian(x); @@ -84,7 +85,7 @@ static void wsum_hess_init_impl(expr *node) static void eval_wsum_hess(expr *node, const double *w) { expr *x = node->left; - const double *a = ((const_vector_mult_expr *) node)->a; + const double *a = ((vector_mult_expr *) node)->param_source->value; /* scale weights w by a */ for (int i = 0; i < node->size; i++) @@ -100,20 +101,23 @@ static void eval_wsum_hess(expr *node, const double *w) static void free_type_data(expr *node) { - const_vector_mult_expr *vnode = (const_vector_mult_expr *) node; - free(vnode->a); + vector_mult_expr *vnode = (vector_mult_expr *) node; + if (vnode->param_source != NULL) + { + free_expr(vnode->param_source); + vnode->param_source = NULL; + } } static bool is_affine(const expr *node) { - /* Affine iff the child is affine */ return node->left->is_affine(node->left); } -expr *new_const_vector_mult(const double *a, expr *child) +expr *new_vector_mult(expr *param_node, expr *child) { - const_vector_mult_expr *vnode = - (const_vector_mult_expr *) SP_CALLOC(1, sizeof(const_vector_mult_expr)); + vector_mult_expr *vnode = + (vector_mult_expr *) SP_CALLOC(1, sizeof(vector_mult_expr)); expr *node = &vnode->base; init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init_impl, @@ -122,9 +126,8 @@ expr *new_const_vector_mult(const double *a, expr *child) node->left = child; expr_retain(child); - /* copy a vector */ - vnode->a = (double *) SP_MALLOC(child->size * sizeof(double)); - memcpy(vnode->a, a, child->size * sizeof(double)); + vnode->param_source = param_node; + expr_retain(param_node); return node; } diff --git a/src/problem.c b/src/problem.c index 399079b..e436e82 100644 --- a/src/problem.c +++ b/src/problem.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "problem.h" +#include "subexpr.h" #include "utils/CSR_sum.h" #include "utils/tracked_alloc.h" #include "utils/utils.h" @@ -335,6 +336,9 @@ void free_problem(problem *prob) print_end_message(&prob->stats); } + /* Free param_nodes array (weak refs, don't free the nodes) */ + free(prob->param_nodes); + /* Free allocated arrays */ free(prob->constraint_values); free(prob->gradient_values); @@ -356,6 +360,35 @@ void free_problem(problem *prob) free(prob); } +void problem_register_params(problem *prob, expr **param_nodes, int n_param_nodes) +{ + prob->n_param_nodes = n_param_nodes; + prob->param_nodes = (expr **) SP_MALLOC(n_param_nodes * sizeof(expr *)); + memcpy(prob->param_nodes, param_nodes, n_param_nodes * sizeof(expr *)); + + prob->total_parameter_size = 0; + for (int i = 0; i < n_param_nodes; i++) + { + prob->total_parameter_size += param_nodes[i]->size; + } +} + +void problem_update_params(problem *prob, const double *theta) +{ + for (int i = 0; i < prob->n_param_nodes; i++) + { + expr *pnode = prob->param_nodes[i]; + parameter_expr *param = (parameter_expr *) pnode; + if (param->param_id == PARAM_FIXED) continue; + int offset = param->param_id; + memcpy(pnode->value, theta + offset, pnode->size * sizeof(double)); + param->needs_refresh = true; + } + + /* Force re-evaluation of affine Jacobians on next call */ + prob->jacobian_called = false; +} + double problem_objective_forward(problem *prob, const double *u) { Timer timer; diff --git a/src/utils/dense_matrix.c b/src/utils/dense_matrix.c index a973abe..445d36e 100644 --- a/src/utils/dense_matrix.c +++ b/src/utils/dense_matrix.c @@ -42,6 +42,12 @@ static void dense_block_left_mult_vec(const Matrix *A, const double *x, double * n, 0.0, y, m); } +static void dense_update_values(Matrix *self, const double *new_values) +{ + Dense_Matrix *dm = (Dense_Matrix *) self; + memcpy(dm->x, new_values, dm->base.m * dm->base.n * sizeof(double)); +} + static void dense_free(Matrix *A) { Dense_Matrix *dm = (Dense_Matrix *) A; @@ -58,6 +64,7 @@ Matrix *new_dense_matrix(int m, int n, const double *data) dm->base.block_left_mult_vec = dense_block_left_mult_vec; dm->base.block_left_mult_sparsity = I_kron_A_alloc; dm->base.block_left_mult_values = I_kron_A_fill_values; + dm->base.update_values = dense_update_values; dm->base.free_fn = dense_free; dm->x = (double *) SP_MALLOC(m * n * sizeof(double)); memcpy(dm->x, data, m * n * sizeof(double)); @@ -71,15 +78,20 @@ Matrix *dense_matrix_trans(const Dense_Matrix *A) int n = A->base.n; double *AT_x = (double *) SP_MALLOC(m * n * sizeof(double)); + A_transpose(AT_x, A->x, m, n); + + Matrix *result = new_dense_matrix(n, m, AT_x); + free(AT_x); + return result; +} + +void A_transpose(double *AT, const double *A, int m, int n) +{ for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { - AT_x[j * m + i] = A->x[i * n + j]; + AT[j * m + i] = A[i * n + j]; } } - - Matrix *result = new_dense_matrix(n, m, AT_x); - free(AT_x); - return result; } diff --git a/src/utils/sparse_matrix.c b/src/utils/sparse_matrix.c index e327886..e413f85 100644 --- a/src/utils/sparse_matrix.c +++ b/src/utils/sparse_matrix.c @@ -19,6 +19,7 @@ #include "utils/matrix.h" #include "utils/tracked_alloc.h" #include +#include static void sparse_block_left_mult_vec(const Matrix *self, const double *x, double *y, int p) @@ -41,6 +42,12 @@ static void sparse_block_left_mult_values(const Matrix *self, const CSC_Matrix * block_left_multiply_fill_values(sm->csr, J, C); } +static void sparse_update_values(Matrix *self, const double *new_values) +{ + Sparse_Matrix *sm = (Sparse_Matrix *) self; + memcpy(sm->csr->x, new_values, sm->csr->nnz * sizeof(double)); +} + static void sparse_free(Matrix *self) { Sparse_Matrix *sm = (Sparse_Matrix *) self; @@ -56,6 +63,7 @@ Matrix *new_sparse_matrix(const CSR_Matrix *A) sm->base.block_left_mult_vec = sparse_block_left_mult_vec; sm->base.block_left_mult_sparsity = sparse_block_left_mult_sparsity; sm->base.block_left_mult_values = sparse_block_left_mult_values; + sm->base.update_values = sparse_update_values; sm->base.free_fn = sparse_free; sm->csr = new_csr(A); return &sm->base; @@ -70,6 +78,7 @@ Matrix *sparse_matrix_trans(const Sparse_Matrix *self, int *iwork) sm->base.block_left_mult_vec = sparse_block_left_mult_vec; sm->base.block_left_mult_sparsity = sparse_block_left_mult_sparsity; sm->base.block_left_mult_values = sparse_block_left_mult_values; + sm->base.update_values = sparse_update_values; sm->base.free_fn = sparse_free; sm->csr = AT; return &sm->base; diff --git a/tests/all_tests.c b/tests/all_tests.c index 2bce8f3..c0bb1b9 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -13,7 +13,7 @@ #include "forward_pass/affine/test_neg.h" #include "forward_pass/affine/test_promote.h" #include "forward_pass/affine/test_sum.h" -#include "forward_pass/affine/test_variable_constant.h" +#include "forward_pass/affine/test_variable_parameter.h" #include "forward_pass/affine/test_vstack.h" #include "forward_pass/bivariate_full_dom/test_matmul.h" #include "forward_pass/composite/test_composite.h" @@ -23,17 +23,17 @@ #include "forward_pass/other/test_prod_axis_one.h" #include "forward_pass/other/test_prod_axis_zero.h" #include "jacobian_tests/affine/test_broadcast.h" -#include "jacobian_tests/affine/test_const_scalar_mult.h" -#include "jacobian_tests/affine/test_const_vector_mult.h" #include "jacobian_tests/affine/test_hstack.h" #include "jacobian_tests/affine/test_index.h" #include "jacobian_tests/affine/test_left_matmul.h" #include "jacobian_tests/affine/test_neg.h" #include "jacobian_tests/affine/test_promote.h" #include "jacobian_tests/affine/test_right_matmul.h" +#include "jacobian_tests/affine/test_scalar_mult.h" #include "jacobian_tests/affine/test_sum.h" #include "jacobian_tests/affine/test_trace.h" #include "jacobian_tests/affine/test_transpose.h" +#include "jacobian_tests/affine/test_vector_mult.h" #include "jacobian_tests/affine/test_vstack.h" #include "jacobian_tests/bivariate_full_dom/test_elementwise_mult.h" #include "jacobian_tests/bivariate_full_dom/test_matmul.h" @@ -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_prob.h" #include "problem/test_problem.h" #include "utils/test_cblas.h" #include "utils/test_coo_matrix.h" @@ -59,15 +60,15 @@ #include "utils/test_linalg_utils_matmul_chain_rule.h" #include "utils/test_matrix.h" #include "wsum_hess/affine/test_broadcast.h" -#include "wsum_hess/affine/test_const_scalar_mult.h" -#include "wsum_hess/affine/test_const_vector_mult.h" #include "wsum_hess/affine/test_hstack.h" #include "wsum_hess/affine/test_index.h" #include "wsum_hess/affine/test_left_matmul.h" #include "wsum_hess/affine/test_right_matmul.h" +#include "wsum_hess/affine/test_scalar_mult.h" #include "wsum_hess/affine/test_sum.h" #include "wsum_hess/affine/test_trace.h" #include "wsum_hess/affine/test_transpose.h" +#include "wsum_hess/affine/test_vector_mult.h" #include "wsum_hess/affine/test_vstack.h" #include "wsum_hess/bivariate_full_dom/test_matmul.h" #include "wsum_hess/bivariate_full_dom/test_multiply.h" @@ -147,10 +148,10 @@ int main(void) mu_run_test(test_jacobian_matmul_sin_Ax_cos_Bx, tests_run); mu_run_test(test_jacobian_matmul_X_X, tests_run); mu_run_test(test_jacobian_composite_exp_add, tests_run); - mu_run_test(test_jacobian_const_scalar_mult_log_vector, tests_run); - mu_run_test(test_jacobian_const_scalar_mult_log_matrix, tests_run); - mu_run_test(test_jacobian_const_vector_mult_log_vector, tests_run); - mu_run_test(test_jacobian_const_vector_mult_log_matrix, tests_run); + mu_run_test(test_jacobian_scalar_mult_log_vector, tests_run); + mu_run_test(test_jacobian_scalar_mult_log_matrix, tests_run); + mu_run_test(test_jacobian_vector_mult_log_vector, tests_run); + mu_run_test(test_jacobian_vector_mult_log_matrix, tests_run); mu_run_test(test_jacobian_rel_entr_vector_args_1, tests_run); mu_run_test(test_jacobian_rel_entr_vector_args_2, tests_run); mu_run_test(test_jacobian_rel_entr_matrix_args, tests_run); @@ -251,10 +252,10 @@ int main(void) mu_run_test(test_wsum_hess_quad_over_lin_xy, tests_run); mu_run_test(test_wsum_hess_quad_over_lin_yx, tests_run); mu_run_test(test_wsum_hess_quad_form, tests_run); - mu_run_test(test_wsum_hess_const_scalar_mult_log_vector, tests_run); - mu_run_test(test_wsum_hess_const_scalar_mult_log_matrix, tests_run); - mu_run_test(test_wsum_hess_const_vector_mult_log_vector, tests_run); - mu_run_test(test_wsum_hess_const_vector_mult_log_matrix, tests_run); + mu_run_test(test_wsum_hess_scalar_mult_log_vector, tests_run); + mu_run_test(test_wsum_hess_scalar_mult_log_matrix, tests_run); + mu_run_test(test_wsum_hess_vector_mult_log_vector, tests_run); + mu_run_test(test_wsum_hess_vector_mult_log_matrix, tests_run); mu_run_test(test_wsum_hess_multiply_linear_ops, tests_run); mu_run_test(test_wsum_hess_multiply_sparse_random, tests_run); mu_run_test(test_wsum_hess_multiply_1, tests_run); @@ -347,6 +348,15 @@ int main(void) mu_run_test(test_problem_jacobian_multi, tests_run); mu_run_test(test_problem_constraint_forward, tests_run); mu_run_test(test_problem_hessian, tests_run); + + printf("\n--- Parameter Tests ---\n"); + mu_run_test(test_param_scalar_mult_problem, tests_run); + mu_run_test(test_param_vector_mult_problem, tests_run); + mu_run_test(test_param_left_matmul_problem, tests_run); + mu_run_test(test_param_right_matmul_problem, tests_run); + mu_run_test(test_param_left_matmul_rectangular, tests_run); + mu_run_test(test_param_right_matmul_rectangular, tests_run); + mu_run_test(test_param_fixed_skip_in_update, tests_run); #endif /* PROFILE_ONLY */ #ifdef PROFILE_ONLY diff --git a/tests/forward_pass/affine/test_add.h b/tests/forward_pass/affine/test_add.h index 0c4e41f..1f8a61a 100644 --- a/tests/forward_pass/affine/test_add.h +++ b/tests/forward_pass/affine/test_add.h @@ -12,7 +12,7 @@ const char *test_addition(void) double u[2] = {3.0, 4.0}; double c[2] = {1.0, 2.0}; expr *var = new_variable(2, 1, 0, 2); - expr *const_node = new_constant(2, 1, 0, c); + expr *const_node = new_parameter(2, 1, PARAM_FIXED, 0, c); expr *sum = new_add(var, const_node); sum->forward(sum, u); double expected[2] = {4.0, 6.0}; diff --git a/tests/forward_pass/affine/test_left_matmul_dense.h b/tests/forward_pass/affine/test_left_matmul_dense.h index e911ba5..e8c8eb1 100644 --- a/tests/forward_pass/affine/test_left_matmul_dense.h +++ b/tests/forward_pass/affine/test_left_matmul_dense.h @@ -24,7 +24,7 @@ const char *test_left_matmul_dense(void) double A_data[9] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}; /* Build expression Z = A @ X */ - expr *Z = new_left_matmul_dense(X, 3, 3, A_data); + expr *Z = new_left_matmul_dense(NULL, X, 3, 3, A_data); /* Variable values in column-major order */ double u[9] = {1.0, 2.0, 3.0, /* first column */ diff --git a/tests/forward_pass/affine/test_sum.h b/tests/forward_pass/affine/test_sum.h index 74614a9..5c91bf2 100644 --- a/tests/forward_pass/affine/test_sum.h +++ b/tests/forward_pass/affine/test_sum.h @@ -18,7 +18,7 @@ const char *test_sum_axis_neg1(void) Stored as: [1, 2, 3, 4, 5, 6] */ double values[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; - expr *const_node = new_constant(3, 2, 0, values); + expr *const_node = new_parameter(3, 2, PARAM_FIXED, 0, values); expr *log_node = new_log(const_node); expr *sum_node = new_sum(log_node, -1); sum_node->forward(sum_node, NULL); @@ -43,7 +43,7 @@ const char *test_sum_axis_0(void) Stored as: [1, 2, 3, 4, 5, 6] */ double values[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; - expr *const_node = new_constant(3, 2, 0, values); + expr *const_node = new_parameter(3, 2, PARAM_FIXED, 0, values); expr *log_node = new_log(const_node); expr *sum_node = new_sum(log_node, 0); sum_node->forward(sum_node, NULL); @@ -70,7 +70,7 @@ const char *test_sum_axis_1(void) Stored as: [1, 2, 3, 4, 5, 6] */ double values[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; - expr *const_node = new_constant(3, 2, 0, values); + expr *const_node = new_parameter(3, 2, PARAM_FIXED, 0, values); expr *log_node = new_log(const_node); expr *sum_node = new_sum(log_node, 1); sum_node->forward(sum_node, NULL); diff --git a/tests/forward_pass/affine/test_variable_constant.h b/tests/forward_pass/affine/test_variable_parameter.h similarity index 91% rename from tests/forward_pass/affine/test_variable_constant.h rename to tests/forward_pass/affine/test_variable_parameter.h index 2878019..7ef80c5 100644 --- a/tests/forward_pass/affine/test_variable_constant.h +++ b/tests/forward_pass/affine/test_variable_parameter.h @@ -21,7 +21,7 @@ const char *test_constant(void) { double c[2] = {5.0, 10.0}; double u[2] = {0.0, 0.0}; - expr *const_node = new_constant(2, 1, 0, c); + expr *const_node = new_parameter(2, 1, PARAM_FIXED, 0, c); const_node->forward(const_node, u); mu_assert("Constant test failed", cmp_double_array(const_node->value, c, 2)); free_expr(const_node); diff --git a/tests/forward_pass/composite/test_composite.h b/tests/forward_pass/composite/test_composite.h index b710aac..a465fd4 100644 --- a/tests/forward_pass/composite/test_composite.h +++ b/tests/forward_pass/composite/test_composite.h @@ -17,7 +17,7 @@ const char *test_composite(void) /* Build tree: log(exp(x) + c) */ expr *var = new_variable(2, 1, 0, 2); expr *exp_node = new_exp(var); - expr *const_node = new_constant(2, 1, 0, c); + expr *const_node = new_parameter(2, 1, PARAM_FIXED, 0, c); expr *sum = new_add(exp_node, const_node); expr *log_node = new_log(sum); diff --git a/tests/forward_pass/other/test_prod_axis_one.h b/tests/forward_pass/other/test_prod_axis_one.h index 663f1f8..b22b102 100644 --- a/tests/forward_pass/other/test_prod_axis_one.h +++ b/tests/forward_pass/other/test_prod_axis_one.h @@ -16,7 +16,7 @@ const char *test_forward_prod_axis_one(void) Stored as: [1, 2, 3, 4, 5, 6] */ double values[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; - expr *const_node = new_constant(2, 3, 0, values); + expr *const_node = new_parameter(2, 3, PARAM_FIXED, 0, values); expr *prod_node = new_prod_axis_one(const_node); prod_node->forward(prod_node, NULL); diff --git a/tests/forward_pass/other/test_prod_axis_zero.h b/tests/forward_pass/other/test_prod_axis_zero.h index fe6b67a..5fffce2 100644 --- a/tests/forward_pass/other/test_prod_axis_zero.h +++ b/tests/forward_pass/other/test_prod_axis_zero.h @@ -16,7 +16,7 @@ const char *test_forward_prod_axis_zero(void) Stored as: [1, 2, 3, 4, 5, 6] */ double values[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; - expr *const_node = new_constant(2, 3, 0, values); + expr *const_node = new_parameter(2, 3, PARAM_FIXED, 0, values); expr *prod_node = new_prod_axis_zero(const_node); prod_node->forward(prod_node, NULL); diff --git a/tests/jacobian_tests/affine/test_broadcast.h b/tests/jacobian_tests/affine/test_broadcast.h index e3d79d2..f41082d 100644 --- a/tests/jacobian_tests/affine/test_broadcast.h +++ b/tests/jacobian_tests/affine/test_broadcast.h @@ -141,7 +141,7 @@ const char *test_double_broadcast(void) /* form the expression x + b */ expr *x = new_variable(5, 1, 0, 5); - expr *b = new_constant(1, 5, 5, b_vals); + expr *b = new_parameter(1, 5, PARAM_FIXED, 5, b_vals); expr *bcast_x = new_broadcast(x, 5, 5); expr *bcast_b = new_broadcast(b, 5, 5); expr *sum = new_add(bcast_x, bcast_b); diff --git a/tests/jacobian_tests/affine/test_left_matmul.h b/tests/jacobian_tests/affine/test_left_matmul.h index c4968d6..41cccc4 100644 --- a/tests/jacobian_tests/affine/test_left_matmul.h +++ b/tests/jacobian_tests/affine/test_left_matmul.h @@ -42,7 +42,7 @@ const char *test_jacobian_left_matmul_log(void) memcpy(A->x, A_x, 7 * sizeof(double)); expr *log_x = new_log(x); - expr *A_log_x = new_left_matmul(log_x, A); + expr *A_log_x = new_left_matmul(NULL, log_x, A); A_log_x->forward(A_log_x, x_vals); jacobian_init(A_log_x); @@ -86,7 +86,7 @@ const char *test_jacobian_left_matmul_log_matrix(void) memcpy(A->x, A_x, 7 * sizeof(double)); expr *log_x = new_log(x); - expr *A_log_x = new_left_matmul(log_x, A); + expr *A_log_x = new_left_matmul(NULL, log_x, A); A_log_x->forward(A_log_x, x_vals); jacobian_init(A_log_x); @@ -135,7 +135,7 @@ const char *test_jacobian_left_matmul_exp_composite(void) expr *Bx = new_linear(x, B, NULL); expr *exp_Bx = new_exp(Bx); - expr *A_exp_Bx = new_left_matmul(exp_Bx, A); + 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)); diff --git a/tests/jacobian_tests/affine/test_right_matmul.h b/tests/jacobian_tests/affine/test_right_matmul.h index e9f40da..8c2fa04 100644 --- a/tests/jacobian_tests/affine/test_right_matmul.h +++ b/tests/jacobian_tests/affine/test_right_matmul.h @@ -27,7 +27,7 @@ const char *test_jacobian_right_matmul_log(void) memcpy(A->x, A_x, 4 * sizeof(double)); expr *log_x = new_log(x); - expr *log_x_A = new_right_matmul(log_x, A); + expr *log_x_A = new_right_matmul(NULL, log_x, A); log_x_A->forward(log_x_A, x_vals); jacobian_init(log_x_A); @@ -76,7 +76,7 @@ const char *test_jacobian_right_matmul_log_vector(void) memcpy(A->x, A_x, 4 * sizeof(double)); expr *log_x = new_log(x); - expr *log_x_A = new_right_matmul(log_x, A); + expr *log_x_A = new_right_matmul(NULL, log_x, A); log_x_A->forward(log_x_A, x_vals); jacobian_init(log_x_A); diff --git a/tests/jacobian_tests/affine/test_const_scalar_mult.h b/tests/jacobian_tests/affine/test_scalar_mult.h similarity index 82% rename from tests/jacobian_tests/affine/test_const_scalar_mult.h rename to tests/jacobian_tests/affine/test_scalar_mult.h index ca27f94..a16a6fa 100644 --- a/tests/jacobian_tests/affine/test_const_scalar_mult.h +++ b/tests/jacobian_tests/affine/test_scalar_mult.h @@ -5,11 +5,12 @@ #include "atoms/elementwise_restricted_dom.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" -/* Test: y = a * log(x) where a is a scalar constant */ +/* Test: y = a * log(x) where a is a scalar parameter */ -const char *test_jacobian_const_scalar_mult_log_vector(void) +const char *test_jacobian_scalar_mult_log_vector(void) { /* Create variable x: [1.0, 2.0, 4.0] with 3 elements */ double u_vals[3] = {1.0, 2.0, 4.0}; @@ -19,8 +20,9 @@ const char *test_jacobian_const_scalar_mult_log_vector(void) expr *log_node = new_log(x); /* Create scalar mult node: y = 2.5 * log(x) */ - double a = 2.5; - expr *y = new_const_scalar_mult(a, log_node); + double a_val = 2.5; + expr *a_param = new_parameter(1, 1, PARAM_FIXED, 3, &a_val); + expr *y = new_scalar_mult(a_param, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -45,7 +47,7 @@ const char *test_jacobian_const_scalar_mult_log_vector(void) return 0; } -const char *test_jacobian_const_scalar_mult_log_matrix(void) +const char *test_jacobian_scalar_mult_log_matrix(void) { /* Create variable x as 2x2 matrix: [[1.0, 2.0], [4.0, 8.0]] */ double u_vals[4] = {1.0, 2.0, 4.0, 8.0}; @@ -55,8 +57,9 @@ const char *test_jacobian_const_scalar_mult_log_matrix(void) expr *log_node = new_log(x); /* Create scalar mult node: y = 3.0 * log(x) */ - double a = 3.0; - expr *y = new_const_scalar_mult(a, log_node); + double a_val = 3.0; + expr *a_param = new_parameter(1, 1, PARAM_FIXED, 4, &a_val); + expr *y = new_scalar_mult(a_param, log_node); /* Forward pass */ y->forward(y, u_vals); diff --git a/tests/jacobian_tests/affine/test_transpose.h b/tests/jacobian_tests/affine/test_transpose.h index 4419315..a02e22d 100644 --- a/tests/jacobian_tests/affine/test_transpose.h +++ b/tests/jacobian_tests/affine/test_transpose.h @@ -21,7 +21,7 @@ const char *test_jacobian_transpose(void) // X = [1 2; 3 4] (columnwise: x = [1 3 2 4]) expr *X = new_variable(2, 2, 0, 4); - expr *AX = new_left_matmul(X, A); + expr *AX = new_left_matmul(NULL, X, A); expr *transpose_AX = new_transpose(AX); double u[4] = {1, 3, 2, 4}; transpose_AX->forward(transpose_AX, u); diff --git a/tests/jacobian_tests/affine/test_const_vector_mult.h b/tests/jacobian_tests/affine/test_vector_mult.h similarity index 77% rename from tests/jacobian_tests/affine/test_const_vector_mult.h rename to tests/jacobian_tests/affine/test_vector_mult.h index fbee7df..cd6fa1a 100644 --- a/tests/jacobian_tests/affine/test_const_vector_mult.h +++ b/tests/jacobian_tests/affine/test_vector_mult.h @@ -5,11 +5,12 @@ #include "atoms/elementwise_restricted_dom.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" -/* Test: y = a ∘ log(x) where a is a constant vector */ +/* Test: y = a ∘ log(x) where a is a parameter vector */ -const char *test_jacobian_const_vector_mult_log_vector(void) +const char *test_jacobian_vector_mult_log_vector(void) { /* Create variable x: [1.0, 2.0, 4.0] */ double u_vals[3] = {1.0, 2.0, 4.0}; @@ -20,7 +21,8 @@ const char *test_jacobian_const_vector_mult_log_vector(void) /* Create vector mult node: y = [2.0, 3.0, 4.0] ∘ log(x) */ double a[3] = {2.0, 3.0, 4.0}; - expr *y = new_const_vector_mult(a, log_node); + expr *a_param = new_parameter(3, 1, PARAM_FIXED, 3, a); + expr *y = new_vector_mult(a_param, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -29,11 +31,6 @@ const char *test_jacobian_const_vector_mult_log_vector(void) jacobian_init(y); y->eval_jacobian(y); - /* Expected jacobian (row-wise scaling): - * Row 0: 2.0 * [1/1] = [2.0] - * Row 1: 3.0 * [1/2] = [1.5] - * Row 2: 4.0 * [1/4] = [1.0] - */ double expected_x[3] = {2.0, 1.5, 1.0}; int expected_p[4] = {0, 1, 2, 3}; int expected_i[3] = {0, 1, 2}; @@ -49,7 +46,7 @@ const char *test_jacobian_const_vector_mult_log_vector(void) return 0; } -const char *test_jacobian_const_vector_mult_log_matrix(void) +const char *test_jacobian_vector_mult_log_matrix(void) { /* Create variable x as 2x2 matrix: [[1.0, 2.0], [4.0, 8.0]] */ double u_vals[4] = {1.0, 2.0, 4.0, 8.0}; @@ -60,7 +57,8 @@ const char *test_jacobian_const_vector_mult_log_matrix(void) /* Create vector mult node: y = [1.5, 2.5, 3.5, 4.5] ∘ log(x) */ double a[4] = {1.5, 2.5, 3.5, 4.5}; - expr *y = new_const_vector_mult(a, log_node); + expr *a_param = new_parameter(4, 1, PARAM_FIXED, 4, a); + expr *y = new_vector_mult(a_param, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -69,12 +67,6 @@ const char *test_jacobian_const_vector_mult_log_matrix(void) jacobian_init(y); y->eval_jacobian(y); - /* Expected jacobian (row-wise scaling): - * Row 0: 1.5 * [1/1] = [1.5] - * Row 1: 2.5 * [1/2] = [1.25] - * Row 2: 3.5 * [1/4] = [0.875] - * Row 3: 4.5 * [1/8] = [0.5625] - */ double expected_x[4] = {1.5, 1.25, 0.875, 0.5625}; int expected_p[5] = {0, 1, 2, 3, 4}; int expected_i[4] = {0, 1, 2, 3}; diff --git a/tests/jacobian_tests/composite/test_chain_rule_jacobian.h b/tests/jacobian_tests/composite/test_chain_rule_jacobian.h index 8139a33..829ac4d 100644 --- a/tests/jacobian_tests/composite/test_chain_rule_jacobian.h +++ b/tests/jacobian_tests/composite/test_chain_rule_jacobian.h @@ -82,8 +82,8 @@ const char *test_jacobian_Ax_Bx_multiply(void) CSR_Matrix *A = new_csr_random(2, 2, 1.0); CSR_Matrix *B = new_csr_random(2, 2, 1.0); expr *x = new_variable(2, 1, 1, 4); - expr *Ax = new_left_matmul(x, A); - expr *Bx = new_left_matmul(x, B); + expr *Ax = new_left_matmul(NULL, x, A); + expr *Bx = new_left_matmul(NULL, x, B); expr *multiply = new_elementwise_mult(Ax, Bx); mu_assert("check_jacobian failed", @@ -102,8 +102,8 @@ const char *test_jacobian_AX_BX_multiply(void) CSR_Matrix *A = new_csr_random(2, 2, 1.0); CSR_Matrix *B = new_csr_random(2, 2, 1.0); expr *X = new_variable(2, 2, 0, 4); - expr *AX = new_left_matmul(X, A); - expr *BX = new_left_matmul(X, B); + expr *AX = new_left_matmul(NULL, X, A); + expr *BX = new_left_matmul(NULL, X, B); expr *multiply = new_elementwise_mult(new_sin(AX), new_cos(BX)); mu_assert("check_jacobian failed", @@ -132,7 +132,7 @@ const char *test_jacobian_quad_form_Ax(void) memcpy(Q->p, Qp, 4 * sizeof(int)); expr *x = new_variable(4, 1, 1, 6); - expr *Ax = new_left_matmul(x, A); + expr *Ax = new_left_matmul(NULL, x, A); expr *sin_Ax = new_sin(Ax); expr *node = new_quad_form(sin_Ax, Q); @@ -215,11 +215,11 @@ const char *test_jacobian_matmul_Ax_By(void) CSR_Matrix *A = new_csr_random(3, 2, 1.0); CSR_Matrix *B = new_csr_random(2, 3, 1.0); - expr *X = new_variable(2, 2, 0, 10); /* 2x2, vars 0-3 */ - expr *Y = new_variable(3, 2, 4, 10); /* 3x2, vars 4-9 */ - expr *AX = new_left_matmul(X, A); /* 3x2 */ - expr *BY = new_left_matmul(Y, B); /* 2x2 */ - expr *Z = new_matmul(AX, BY); /* 3x2 */ + expr *X = new_variable(2, 2, 0, 10); /* 2x2, vars 0-3 */ + expr *Y = new_variable(3, 2, 4, 10); /* 3x2, vars 4-9 */ + expr *AX = new_left_matmul(NULL, X, A); /* 3x2 */ + expr *BY = new_left_matmul(NULL, Y, B); /* 2x2 */ + expr *Z = new_matmul(AX, BY); /* 3x2 */ mu_assert("check_jacobian failed", check_jacobian(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); @@ -238,12 +238,12 @@ const char *test_jacobian_matmul_sin_Ax_cos_Bx(void) CSR_Matrix *A = new_csr_random(2, 3, 1.0); CSR_Matrix *B = new_csr_random(2, 3, 1.0); - expr *X = new_variable(3, 2, 0, 6); /* 3x2, vars 0-5 */ - expr *AX = new_left_matmul(X, A); /* 2x2 */ - expr *BX = new_left_matmul(X, B); /* 2x2 */ - expr *sin_AX = new_sin(AX); /* 2x2 */ - expr *cos_BX = new_cos(BX); /* 2x2 */ - expr *Z = new_matmul(sin_AX, cos_BX); /* 2x2 */ + expr *X = new_variable(3, 2, 0, 6); /* 3x2, vars 0-5 */ + expr *AX = new_left_matmul(NULL, X, A); /* 2x2 */ + expr *BX = new_left_matmul(NULL, X, B); /* 2x2 */ + expr *sin_AX = new_sin(AX); /* 2x2 */ + expr *cos_BX = new_cos(BX); /* 2x2 */ + expr *Z = new_matmul(sin_AX, cos_BX); /* 2x2 */ mu_assert("check_jacobian failed", check_jacobian(Z, u_vals, NUMERICAL_DIFF_DEFAULT_H)); diff --git a/tests/problem/test_param_prob.h b/tests/problem/test_param_prob.h new file mode 100644 index 0000000..76b24e3 --- /dev/null +++ b/tests/problem/test_param_prob.h @@ -0,0 +1,403 @@ +#ifndef TEST_PARAM_PROB_H +#define TEST_PARAM_PROB_H + +#include +#include +#include + +#include "atoms/affine.h" +#include "atoms/elementwise_restricted_dom.h" +#include "expr.h" +#include "minunit.h" +#include "problem.h" +#include "subexpr.h" +#include "test_helpers.h" + +const char *test_param_scalar_mult_problem(void) +{ + int n = 2; + + /* 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); + expr *scaled = new_scalar_mult(a_param, log_x); + 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 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); + 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); + 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; + + /* 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); + expr *constraint = new_vector_mult(p_param, 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 and utilities for test */ + double x_vals[2] = {1.0, 2.0}; + int Ap[3] = {0, 1, 2}; + int Ai[2] = {0, 1}; + 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); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, Ax, 2)); + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 3)); + mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 2)); + + /* test 2: p=[5,6] */ + theta[0] = 5.0; + theta[1] = 6.0; + problem_update_params(prob, theta); + problem_constraint_forward(prob, x_vals); + problem_jacobian(prob); + constrs[0] = 5.0; + constrs[1] = 12.0; + Ax[0] = 5.0; + Ax[1] = 6.0; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, Ax, 2)); + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 3)); + mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 2)); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + + free_problem(prob); + + return 0; +} + +const char *test_param_left_matmul_problem(void) +{ + int n = 2; + + /* 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); + + /* dense 2x2 matrix */ + double Ax[4] = {2.0, 0.0, 0.0, 1.0}; + expr *constraint = new_left_matmul_dense(A_param, x, 2, 2, Ax); + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, 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 and utilities for test */ + double x_vals[2] = {1.0, 2.0}; + int Ap[3] = {0, 2, 4}; + int Ai[4] = {0, 1, 0, 1}; + + /* test 1: initial jacobian */ + problem_constraint_forward(prob, x_vals); + double constrs[2] = {2.0, 2.0}; + problem_jacobian(prob); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, Ax, 4)); + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 3)); + 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}; + problem_update_params(prob, theta); + problem_constraint_forward(prob, x_vals); + problem_jacobian(prob); + constrs[0] = 5; + constrs[1] = 11; + Ax[0] = 1.0; + Ax[1] = 2.0; + Ax[2] = 3.0; + Ax[3] = 4.0; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, Ax, 4)); + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 3)); + mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 4)); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + + free_problem(prob); + + return 0; +} + +const char *test_param_right_matmul_problem(void) +{ + int n = 2; + + /* 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); + + /* dense 2x2 matrix */ + double Ax[4] = {2.0, 0.0, 0.0, 1.0}; + expr *constraint = new_right_matmul_dense(A_param, x, 2, 2, Ax); + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, 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 and utilities for test */ + double x_vals[2] = {1.0, 2.0}; + int Ap[3] = {0, 2, 4}; + int Ai[4] = {0, 1, 0, 1}; + + /* test 1: initial jacobian */ + problem_constraint_forward(prob, x_vals); + double constrs[2] = {2.0, 2.0}; + problem_jacobian(prob); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, Ax, 4)); + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 3)); + 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}; + problem_update_params(prob, theta); + problem_constraint_forward(prob, x_vals); + problem_jacobian(prob); + constrs[0] = 7; + constrs[1] = 10; + Ax[0] = 1.0; + Ax[1] = 3.0; + Ax[2] = 2.0; + Ax[3] = 4.0; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, Ax, 4)); + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 3)); + mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 4)); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 2)); + + free_problem(prob); + + return 0; +} + +const char *test_param_fixed_skip_in_update(void) +{ + int n = 2; + + /* minimize a * sum(log(x)) + b * sum(x), a fixed, b updatable */ + expr *x = new_variable(2, 1, 0, n); + expr *log_x = new_log(x); + double a_val = 2.0; + expr *a_param = new_parameter(1, 1, PARAM_FIXED, n, &a_val); + 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); + expr *b_x = new_scalar_mult(b_param, x); + expr *sum_b_x = new_sum(b_x, -1); + + expr *objective = new_add(sum_a_log, sum_b_x); + 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); + problem_init_derivatives(prob); + + /* point for evaluating */ + double x_vals[2] = {1.0, 2.0}; + + /* test 1: b=3, a stays 2 */ + double theta[1] = {3.0}; + problem_update_params(prob, theta); + mu_assert("vals fail", fabs(a_param->value[0] - 2.0) < 1e-10); + double obj_val = problem_objective_forward(prob, x_vals); + problem_gradient(prob); + double expected_obj = 2.0 * log(2.0) + 9.0; + mu_assert("vals fail", fabs(obj_val - expected_obj) < 1e-10); + double grad[2] = {5.0, 4.0}; + mu_assert("vals fail", cmp_double_array(prob->gradient_values, grad, 2)); + + /* test 2: b=5, a stays 2 */ + theta[0] = 5.0; + problem_update_params(prob, theta); + mu_assert("vals fail", fabs(a_param->value[0] - 2.0) < 1e-10); + obj_val = problem_objective_forward(prob, x_vals); + problem_gradient(prob); + expected_obj = 2.0 * log(2.0) + 15.0; + mu_assert("vals fail", fabs(obj_val - expected_obj) < 1e-10); + grad[0] = 7.0; + grad[1] = 6.0; + mu_assert("vals fail", cmp_double_array(prob->gradient_values, grad, 2)); + + free_problem(prob); + + return 0; +} + +const char *test_param_left_matmul_rectangular(void) +{ + int n = 2; + + /* 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); + + /* dense 3x2 matrix */ + double Ax[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + expr *constraint = new_left_matmul_dense(A_param, 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] = {A_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* point for evaluating and utilities for test */ + double x_vals[2] = {1.0, 2.0}; + int Ap[4] = {0, 2, 4, 6}; + int Ai[6] = {0, 1, 0, 1, 0, 1}; + + /* test 1: initial jacobian */ + problem_constraint_forward(prob, x_vals); + double constrs[3] = {5.0, 11.0, 17.0}; + problem_jacobian(prob); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, Ax, 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: 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}; + problem_update_params(prob, theta); + problem_constraint_forward(prob, x_vals); + problem_jacobian(prob); + constrs[0] = 23.0; + constrs[1] = 29.0; + constrs[2] = 35.0; + Ax[0] = 7.0; + Ax[1] = 8.0; + Ax[2] = 9.0; + Ax[3] = 10.0; + Ax[4] = 11.0; + Ax[5] = 12.0; + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + mu_assert("vals fail", cmp_double_array(prob->jacobian->x, Ax, 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)); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + + free_problem(prob); + + return 0; +} + +const char *test_param_right_matmul_rectangular(void) +{ + int n = 2; + + /* 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); + + /* dense 2x3 matrix */ + double Ax[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + expr *constraint = new_right_matmul_dense(A_param, x, 2, 3, Ax); + expr *constraints[1] = {constraint}; + problem *prob = new_problem(objective, constraints, 1, 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 and utilities for test */ + double x_vals[2] = {1.0, 2.0}; + int Ap[4] = {0, 2, 4, 6}; + int Ai[6] = {0, 1, 0, 1, 0, 1}; + + /* test 1: initial jacobian */ + problem_constraint_forward(prob, x_vals); + double constrs[3] = {9.0, 12.0, 15.0}; + problem_jacobian(prob); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + 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 jac_x[6] = {1.0, 4.0, 2.0, 5.0, 3.0, 6.0}; + 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}; + problem_update_params(prob, theta); + problem_constraint_forward(prob, x_vals); + problem_jacobian(prob); + constrs[0] = 27.0; + constrs[1] = 30.0; + constrs[2] = 33.0; + jac_x[0] = 7.0; + jac_x[1] = 10.0; + jac_x[2] = 8.0; + jac_x[3] = 11.0; + jac_x[4] = 9.0; + jac_x[5] = 12.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)); + mu_assert("vals fail", cmp_double_array(prob->constraint_values, constrs, 3)); + + free_problem(prob); + + return 0; +} + +#endif /* TEST_PARAM_PROB_H */ diff --git a/tests/profiling/profile_left_matmul.h b/tests/profiling/profile_left_matmul.h index df80a90..4644a41 100644 --- a/tests/profiling/profile_left_matmul.h +++ b/tests/profiling/profile_left_matmul.h @@ -31,7 +31,7 @@ const char *profile_left_matmul(void) } A->p[n] = n * n; - expr *AX = new_left_matmul(X, A); + expr *AX = new_left_matmul(NULL, X, A); double *x_vals = (double *) malloc(n * n * sizeof(double)); for (int i = 0; i < n * n; i++) diff --git a/tests/wsum_hess/affine/test_left_matmul.h b/tests/wsum_hess/affine/test_left_matmul.h index f5441b9..313624a 100644 --- a/tests/wsum_hess/affine/test_left_matmul.h +++ b/tests/wsum_hess/affine/test_left_matmul.h @@ -63,7 +63,7 @@ const char *test_wsum_hess_left_matmul(void) memcpy(A->x, A_x, 7 * sizeof(double)); expr *log_x = new_log(x); - expr *A_log_x = new_left_matmul(log_x, A); + expr *A_log_x = new_left_matmul(NULL, log_x, A); A_log_x->forward(A_log_x, x_vals); jacobian_init(A_log_x); @@ -118,7 +118,7 @@ const char *test_wsum_hess_left_matmul_exp_composite(void) expr *Bx = new_linear(x, B, NULL); expr *exp_Bx = new_exp(Bx); - expr *A_exp_Bx = new_left_matmul(exp_Bx, A); + expr *A_exp_Bx = new_left_matmul(NULL, exp_Bx, A); mu_assert("check_wsum_hess failed", check_wsum_hess(A_exp_Bx, x_vals, w, NUMERICAL_DIFF_DEFAULT_H)); @@ -170,7 +170,7 @@ const char *test_wsum_hess_left_matmul_matrix(void) memcpy(A->x, A_x, 7 * sizeof(double)); expr *log_x = new_log(x); - expr *A_log_x = new_left_matmul(log_x, A); + expr *A_log_x = new_left_matmul(NULL, log_x, A); A_log_x->forward(A_log_x, x_vals); jacobian_init(A_log_x); diff --git a/tests/wsum_hess/affine/test_right_matmul.h b/tests/wsum_hess/affine/test_right_matmul.h index a74df5c..2f51b9b 100644 --- a/tests/wsum_hess/affine/test_right_matmul.h +++ b/tests/wsum_hess/affine/test_right_matmul.h @@ -33,7 +33,7 @@ const char *test_wsum_hess_right_matmul(void) memcpy(A->x, A_x, 4 * sizeof(double)); expr *log_x = new_log(x); - expr *log_x_A = new_right_matmul(log_x, A); + expr *log_x_A = new_right_matmul(NULL, log_x, A); log_x_A->forward(log_x_A, x_vals); jacobian_init(log_x_A); @@ -83,7 +83,7 @@ const char *test_wsum_hess_right_matmul_vector(void) memcpy(A->x, A_x, 4 * sizeof(double)); expr *log_x = new_log(x); - expr *log_x_A = new_right_matmul(log_x, A); + expr *log_x_A = new_right_matmul(NULL, log_x, A); log_x_A->forward(log_x_A, x_vals); jacobian_init(log_x_A); diff --git a/tests/wsum_hess/affine/test_const_scalar_mult.h b/tests/wsum_hess/affine/test_scalar_mult.h similarity index 72% rename from tests/wsum_hess/affine/test_const_scalar_mult.h rename to tests/wsum_hess/affine/test_scalar_mult.h index 08e6d86..ba6bae9 100644 --- a/tests/wsum_hess/affine/test_const_scalar_mult.h +++ b/tests/wsum_hess/affine/test_scalar_mult.h @@ -6,11 +6,12 @@ #include "atoms/elementwise_restricted_dom.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" -/* Test: y = a * log(x) where a is a scalar constant */ +/* Test: y = a * log(x) where a is a scalar parameter */ -const char *test_wsum_hess_const_scalar_mult_log_vector(void) +const char *test_wsum_hess_scalar_mult_log_vector(void) { /* Create variable x: [1.0, 2.0, 4.0] */ double u_vals[3] = {1.0, 2.0, 4.0}; @@ -20,8 +21,9 @@ const char *test_wsum_hess_const_scalar_mult_log_vector(void) expr *log_node = new_log(x); /* Create scalar mult node: y = 2.5 * log(x) */ - double a = 2.5; - expr *y = new_const_scalar_mult(a, log_node); + double a_val = 2.5; + expr *a_param = new_parameter(1, 1, PARAM_FIXED, 3, &a_val); + expr *y = new_scalar_mult(a_param, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -32,15 +34,6 @@ const char *test_wsum_hess_const_scalar_mult_log_vector(void) double w[3] = {1.0, 0.5, 0.25}; y->eval_wsum_hess(y, w); - /* For y = a * log(x), the Hessian is: - * H = a * H_log = a * diag([-1/x_i^2]) - * With weights w and scalar a: - * H_weighted = a * diag([-w_i / x_i^2]) - * - * Expected diagonal: 2.5 * [-1/1^2, -0.5/2^2, -0.25/4^2] - * = 2.5 * [-1, -0.125, -0.015625] - * = [-2.5, -0.3125, -0.0390625] - */ double expected_x[3] = {-2.5, -0.3125, -0.0390625}; int expected_p[4] = {0, 1, 2, 3}; int expected_i[3] = {0, 1, 2}; @@ -56,7 +49,7 @@ const char *test_wsum_hess_const_scalar_mult_log_vector(void) return 0; } -const char *test_wsum_hess_const_scalar_mult_log_matrix(void) +const char *test_wsum_hess_scalar_mult_log_matrix(void) { /* Create variable x as 2x2 matrix: [[1.0, 2.0], [4.0, 8.0]] */ double u_vals[4] = {1.0, 2.0, 4.0, 8.0}; @@ -66,8 +59,9 @@ const char *test_wsum_hess_const_scalar_mult_log_matrix(void) expr *log_node = new_log(x); /* Create scalar mult node: y = 3.0 * log(x) */ - double a = 3.0; - expr *y = new_const_scalar_mult(a, log_node); + double a_val = 3.0; + expr *a_param = new_parameter(1, 1, PARAM_FIXED, 4, &a_val); + expr *y = new_scalar_mult(a_param, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -78,10 +72,6 @@ const char *test_wsum_hess_const_scalar_mult_log_matrix(void) double w[4] = {1.0, 1.0, 1.0, 1.0}; y->eval_wsum_hess(y, w); - /* Expected diagonal: 3.0 * [-1/1^2, -1/2^2, -1/4^2, -1/8^2] - * = 3.0 * [-1, -0.25, -0.0625, -0.015625] - * = [-3.0, -0.75, -0.1875, -0.046875] - */ double expected_x[4] = {-3.0, -0.75, -0.1875, -0.046875}; int expected_p[5] = {0, 1, 2, 3, 4}; int expected_i[4] = {0, 1, 2, 3}; diff --git a/tests/wsum_hess/affine/test_const_vector_mult.h b/tests/wsum_hess/affine/test_vector_mult.h similarity index 73% rename from tests/wsum_hess/affine/test_const_vector_mult.h rename to tests/wsum_hess/affine/test_vector_mult.h index d265c71..1d75ad3 100644 --- a/tests/wsum_hess/affine/test_const_vector_mult.h +++ b/tests/wsum_hess/affine/test_vector_mult.h @@ -6,11 +6,12 @@ #include "atoms/elementwise_restricted_dom.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" -/* Test: y = a ∘ log(x) where a is a constant vector */ +/* Test: y = a ∘ log(x) where a is a parameter vector */ -const char *test_wsum_hess_const_vector_mult_log_vector(void) +const char *test_wsum_hess_vector_mult_log_vector(void) { /* Create variable x: [1.0, 2.0, 4.0] */ double u_vals[3] = {1.0, 2.0, 4.0}; @@ -21,7 +22,8 @@ const char *test_wsum_hess_const_vector_mult_log_vector(void) /* Create vector mult node: y = [2.0, 3.0, 4.0] ∘ log(x) */ double a[3] = {2.0, 3.0, 4.0}; - expr *y = new_const_vector_mult(a, log_node); + expr *a_param = new_parameter(3, 1, PARAM_FIXED, 3, a); + expr *y = new_vector_mult(a_param, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -32,13 +34,6 @@ const char *test_wsum_hess_const_vector_mult_log_vector(void) double w[3] = {1.0, 0.5, 0.25}; y->eval_wsum_hess(y, w); - /* For y = a ∘ log(x), the weighted Hessian is: - * H = diag([a_i * (-w_i / x_i^2)]) - * - * Expected diagonal: [2.0 * (-1 / 1^2), 3.0 * (-0.5 / 2^2), 4.0 * (-0.25 / 4^2)] - * = [2.0 * (-1), 3.0 * (-0.125), 4.0 * (-0.015625)] - * = [-2.0, -0.375, -0.0625] - */ double expected_x[3] = {-2.0, -0.375, -0.0625}; int expected_p[4] = {0, 1, 2, 3}; int expected_i[3] = {0, 1, 2}; @@ -54,7 +49,7 @@ const char *test_wsum_hess_const_vector_mult_log_vector(void) return 0; } -const char *test_wsum_hess_const_vector_mult_log_matrix(void) +const char *test_wsum_hess_vector_mult_log_matrix(void) { /* Create variable x as 2x2 matrix: [[1.0, 2.0], [4.0, 8.0]] */ double u_vals[4] = {1.0, 2.0, 4.0, 8.0}; @@ -65,7 +60,8 @@ const char *test_wsum_hess_const_vector_mult_log_matrix(void) /* Create vector mult node: y = [1.5, 2.5, 3.5, 4.5] ∘ log(x) */ double a[4] = {1.5, 2.5, 3.5, 4.5}; - expr *y = new_const_vector_mult(a, log_node); + expr *a_param = new_parameter(4, 1, PARAM_FIXED, 4, a); + expr *y = new_vector_mult(a_param, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -76,11 +72,6 @@ const char *test_wsum_hess_const_vector_mult_log_matrix(void) double w[4] = {1.0, 1.0, 1.0, 1.0}; y->eval_wsum_hess(y, w); - /* Expected diagonal: [1.5 * (-1 / 1^2), 2.5 * (-1 / 2^2), - * 3.5 * (-1 / 4^2), 4.5 * (-1 / 8^2)] - * = [1.5 * (-1), 2.5 * (-0.25), 3.5 * (-0.0625), 4.5 * - * (-0.015625)] = [-1.5, -0.625, -0.21875, -0.0703125] - */ double expected_x[4] = {-1.5, -0.625, -0.21875, -0.0703125}; int expected_p[5] = {0, 1, 2, 3, 4}; int expected_i[4] = {0, 1, 2, 3}; diff --git a/tests/wsum_hess/composite/test_chain_rule_wsum_hess.h b/tests/wsum_hess/composite/test_chain_rule_wsum_hess.h index ca414f7..c5c1a60 100644 --- a/tests/wsum_hess/composite/test_chain_rule_wsum_hess.h +++ b/tests/wsum_hess/composite/test_chain_rule_wsum_hess.h @@ -126,8 +126,8 @@ const char *test_wsum_hess_Ax_Bx_multiply(void) CSR_Matrix *A = new_csr_random(2, 2, 1.0); CSR_Matrix *B = new_csr_random(2, 2, 1.0); expr *x = new_variable(2, 1, 1, 4); - expr *Ax = new_left_matmul(x, A); - expr *Bx = new_left_matmul(x, B); + expr *Ax = new_left_matmul(NULL, x, A); + expr *Bx = new_left_matmul(NULL, x, B); expr *multiply = new_elementwise_mult(Ax, Bx); mu_assert("check_wsum_hess failed", @@ -162,8 +162,8 @@ const char *test_wsum_hess_AX_BX_multiply(void) CSR_Matrix *A = new_csr_random(2, 2, 1.0); CSR_Matrix *B = new_csr_random(2, 2, 1.0); expr *X = new_variable(2, 2, 0, 4); - expr *AX = new_left_matmul(X, A); - expr *BX = new_left_matmul(X, B); + expr *AX = new_left_matmul(NULL, X, A); + expr *BX = new_left_matmul(NULL, X, B); expr *multiply = new_elementwise_mult(new_sin(AX), new_cos(BX)); mu_assert("check_wsum_hess failed", @@ -184,8 +184,8 @@ const char *test_wsum_hess_multiply_deep_composite(void) CSR_Matrix *B = new_csr_random(2, 2, 1.0); expr *X = new_variable(2, 2, 0, 8); expr *Y = new_variable(2, 2, 0, 8); - expr *AX = new_left_matmul(X, A); - expr *BY = new_left_matmul(Y, B); + expr *AX = new_left_matmul(NULL, X, A); + expr *BY = new_left_matmul(NULL, Y, B); expr *sin_AX = new_sin(AX); expr *cos_BY = new_cos(BY); expr *sin_AX_mult_sin_AX = new_elementwise_mult(sin_AX, sin_AX); @@ -217,7 +217,7 @@ const char *test_wsum_hess_quad_form_Ax(void) memcpy(Q->p, Qp, 4 * sizeof(int)); expr *x = new_variable(4, 1, 1, 6); - expr *Ax = new_left_matmul(x, A); + expr *Ax = new_left_matmul(NULL, x, A); expr *node = new_quad_form(Ax, Q); mu_assert("check_wsum_hess failed", @@ -246,7 +246,7 @@ const char *test_wsum_hess_quad_form_sin_Ax(void) memcpy(Q->p, Qp, 4 * sizeof(int)); expr *x = new_variable(4, 1, 1, 6); - expr *Ax = new_left_matmul(x, A); + expr *Ax = new_left_matmul(NULL, x, A); expr *sin_Ax = new_sin(Ax); expr *node = new_quad_form(sin_Ax, Q); @@ -308,9 +308,9 @@ const char *test_wsum_hess_matmul_Ax_By(void) expr *X = new_variable(2, 2, 0, 10); expr *Y = new_variable(3, 2, 4, 10); - expr *AX = new_left_matmul(X, A); /* 3x2 */ - expr *BY = new_left_matmul(Y, B); /* 2x2 */ - expr *Z = new_matmul(AX, BY); /* 3x2 */ + expr *AX = new_left_matmul(NULL, X, A); /* 3x2 */ + expr *BY = new_left_matmul(NULL, Y, B); /* 2x2 */ + expr *Z = new_matmul(AX, BY); /* 3x2 */ mu_assert("check_wsum_hess failed", check_wsum_hess(Z, u_vals, w, NUMERICAL_DIFF_DEFAULT_H)); @@ -331,8 +331,8 @@ const char *test_wsum_hess_matmul_sin_Ax_cos_Bx(void) CSR_Matrix *B = new_csr_random(2, 3, 1.0); expr *X = new_variable(3, 2, 0, 6); - expr *AX = new_left_matmul(X, A); /* 2x2 */ - expr *BX = new_left_matmul(X, B); /* 2x2 */ + expr *AX = new_left_matmul(NULL, X, A); /* 2x2 */ + expr *BX = new_left_matmul(NULL, X, B); /* 2x2 */ expr *sin_AX = new_sin(AX); expr *cos_BX = new_cos(BX); expr *Z = new_matmul(sin_AX, cos_BX); /* 2x2 */