diff --git a/include/affine.h b/include/affine.h index ada58f1..2210dd2 100644 --- a/include/affine.h +++ b/include/affine.h @@ -32,8 +32,8 @@ expr *new_hstack(expr **args, int n_args, int n_vars); expr *new_promote(expr *child, int d1, int d2); expr *new_trace(expr *child); -expr *new_constant(int d1, int d2, int n_vars, const double *values); expr *new_variable(int d1, int d2, int var_id, int n_vars); +expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *values); expr *new_index(expr *child, int d1, int d2, const int *indices, int n_idxs); expr *new_reshape(expr *child, int d1, int d2); diff --git a/include/bivariate.h b/include/bivariate.h index aa005ed..a832d14 100644 --- a/include/bivariate.h +++ b/include/bivariate.h @@ -30,16 +30,17 @@ expr *new_rel_entr_second_arg_scalar(expr *left, expr *right); /* Matrix multiplication: Z = X @ Y */ expr *new_matmul(expr *x, expr *y); -/* Left matrix multiplication: A @ f(x) where A is a constant matrix */ -expr *new_left_matmul(expr *u, const CSR_Matrix *A); +/* Left matrix multiplication: A @ f(x) where A comes from a parameter node. + Only the forward pass possibly updates the parameter. */ +expr *new_left_matmul(expr *param_node, expr *child, const CSR_Matrix *A); -/* 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 comes from a parameter node. */ +expr *new_right_matmul(expr *param_node, expr *u, const CSR_Matrix *A); -/* Constant scalar multiplication: a * f(x) where a is a constant double */ -expr *new_const_scalar_mult(double a, expr *child); +/* Scalar multiplication: a * f(x) where a comes from a parameter node */ +expr *new_scalar_mult(expr *param_node, expr *child); -/* Constant vector elementwise multiplication: a ∘ f(x) where a is constant */ -expr *new_const_vector_mult(const double *a, expr *child); +/* Vector elementwise multiplication: a ∘ f(x) where a comes from a parameter node */ +expr *new_vector_mult(expr *param_node, expr *child); #endif /* BIVARIATE_H */ diff --git a/include/memory_wrappers.h b/include/memory_wrappers.h new file mode 100644 index 0000000..00ed6bd --- /dev/null +++ b/include/memory_wrappers.h @@ -0,0 +1,13 @@ +#ifndef MEMORY_WRAPPERS_H +#define MEMORY_WRAPPERS_H + +#include + +#define FREE_AND_NULL(p) \ + do \ + { \ + free(p); \ + (p) = NULL; \ + } while (0) + +#endif /* MEMORY_WRAPPERS_H */ diff --git a/include/problem.h b/include/problem.h index 2462ffd..c84be7d 100644 --- a/include/problem.h +++ b/include/problem.h @@ -59,6 +59,11 @@ typedef struct problem * hessian are called */ bool jacobian_called; + /* Parameter tracking for fast parameter updates. */ + expr **param_nodes; /* weak references to parameter nodes in tree */ + int n_param_nodes; + int total_parameter_size; + /* Statistics for performance measurement */ Diff_engine_stats stats; bool verbose; @@ -78,4 +83,8 @@ void problem_gradient(problem *prob); void problem_jacobian(problem *prob); void problem_hessian(problem *prob, double obj_w, const double *w); +/* Parameter support */ +void problem_register_params(problem *prob, expr **param_nodes, int n_param_nodes); +void problem_update_params(problem *prob, const double *theta); + #endif diff --git a/include/subexpr.h b/include/subexpr.h index b3f4170..fdb57e4 100644 --- a/include/subexpr.h +++ b/include/subexpr.h @@ -25,6 +25,20 @@ /* Forward declaration */ struct int_double_pair; +/* param_id value for fixed (constant) parameters */ +#define PARAM_FIXED -1 + +/* Parameter node: unified leaf for constants and updatable parameters. + * Constants use param_id == PARAM_FIXED and have values set at creation. + * Updatable parameters have param_id >= 0 and are updated via problem_update_params. + */ +typedef struct parameter_expr +{ + expr base; + int param_id; /* offset into global theta vector, or PARAM_FIXED */ + bool has_been_refreshed; /* tracks whether parameter has been refreshed */ +} parameter_expr; + /* Type-specific expression structures that "inherit" from expr */ /* Linear operator: y = A * x + b */ @@ -113,6 +127,9 @@ typedef struct left_matmul_expr CSC_Matrix *Jchild_CSC; CSC_Matrix *J_CSC; int *csc_to_csr_workspace; + int *AT_iwork; /* work for computing AT values from A */ + expr *param_source; /* parameter node; A/AT values are refreshed from this */ + void (*refresh_param_values)(struct left_matmul_expr *lin_node); } left_matmul_expr; /* Right matrix multiplication: y = f(x) * A where f(x) is an expression. @@ -126,19 +143,20 @@ typedef struct right_matmul_expr CSC_Matrix *CSC_work; } right_matmul_expr; -/* Constant scalar multiplication: y = a * child where a is a constant double */ -typedef struct const_scalar_mult_expr +/* Scalar multiplication: y = a * child where a comes from a parameter node */ +typedef struct scalar_mult_expr { expr base; - double a; -} const_scalar_mult_expr; + expr *param_source; /* always set; read a from param_source->value[0] */ +} scalar_mult_expr; -/* Constant vector elementwise multiplication: y = a \circ child for constant a */ -typedef struct const_vector_mult_expr +/* Vector elementwise multiplication: y = a \circ child where a comes from a + * parameter node */ +typedef struct vector_mult_expr { expr base; - double *a; /* length equals node->size */ -} const_vector_mult_expr; + expr *param_source; /* always set; read a from param_source->value */ +} 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/Vec_macros.h b/include/utils/Vec_macros.h index def45e4..4f0c8ba 100644 --- a/include/utils/Vec_macros.h +++ b/include/utils/Vec_macros.h @@ -19,6 +19,7 @@ #ifndef VEC_MACROS_H #define VEC_MACROS_H +#include "memory_wrappers.h" #include #include #include @@ -48,7 +49,7 @@ vec->data = (TYPE *) malloc(capacity * sizeof(TYPE)); \ if (vec->data == NULL) \ { \ - free(vec); \ + FREE_AND_NULL(vec); \ return NULL; \ } \ \ @@ -59,8 +60,9 @@ \ static inline void TYPE_NAME##Vec_free(TYPE_NAME##Vec *vec) \ { \ - free(vec->data); \ - free(vec); \ + if (!vec) return; \ + FREE_AND_NULL(vec->data); \ + FREE_AND_NULL(vec); \ } \ \ static inline void TYPE_NAME##Vec_clear_no_resize(TYPE_NAME##Vec *vec) \ diff --git a/src/affine/hstack.c b/src/affine/hstack.c index e2235d6..39b6057 100644 --- a/src/affine/hstack.c +++ b/src/affine/hstack.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "affine.h" +#include "memory_wrappers.h" #include "utils/CSR_sum.h" #include #include @@ -165,13 +166,10 @@ static void free_type_data(expr *node) for (int i = 0; i < hnode->n_args; i++) { free_expr(hnode->args[i]); - hnode->args[i] = NULL; } free_csr_matrix(hnode->CSR_work); - hnode->CSR_work = NULL; - free(hnode->args); - hnode->args = NULL; + FREE_AND_NULL(hnode->args); } expr *new_hstack(expr **args, int n_args, int n_vars) diff --git a/src/affine/index.c b/src/affine/index.c index 9577a05..c646d9c 100644 --- a/src/affine/index.c +++ b/src/affine/index.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "affine.h" +#include "memory_wrappers.h" #include "subexpr.h" #include #include @@ -38,7 +39,7 @@ static bool check_for_duplicates(const int *indices, int n_idxs, int max_idx) } seen[indices[i]] = true; } - free(seen); + FREE_AND_NULL(seen); return has_dup; } @@ -154,11 +155,7 @@ static bool is_affine(const expr *node) static void free_type_data(expr *node) { index_expr *idx = (index_expr *) node; - if (idx->indices) - { - free(idx->indices); - idx->indices = NULL; - } + FREE_AND_NULL(idx->indices); } expr *new_index(expr *child, int d1, int d2, const int *indices, int n_idxs) diff --git a/src/affine/linear_op.c b/src/affine/linear_op.c index a8d2863..1eb252e 100644 --- a/src/affine/linear_op.c +++ b/src/affine/linear_op.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "affine.h" +#include "memory_wrappers.h" #include #include #include @@ -55,18 +56,12 @@ static void free_type_data(expr *node) if (!node->jacobian) { free_csr_matrix(lin_node->A_csr); + lin_node->A_csr = NULL; } free_csc_matrix(lin_node->A_csc); - - if (lin_node->b != NULL) - { - free(lin_node->b); - lin_node->b = NULL; - } - - lin_node->A_csr = NULL; lin_node->A_csc = NULL; + FREE_AND_NULL(lin_node->b); } static void jacobian_init(expr *node) diff --git a/src/affine/constant.c b/src/affine/parameter.c similarity index 56% rename from src/affine/constant.c rename to src/affine/parameter.c index 59da3e2..7564827 100644 --- a/src/affine/constant.c +++ b/src/affine/parameter.c @@ -15,39 +15,48 @@ * See the License for the specific language governing permissions and * limitations under the License. */ + +/* Unified parameter/constant leaf node. + * + * When param_id == PARAM_FIXED, this is a constant whose values are set at + * creation and never change. When param_id >= 0, values are updated via + * problem_update_params. + * + * In both cases the derivative behavior is identical: zero Jacobian and + * Hessian with respect to variables (always affine). */ + #include "affine.h" +#include "subexpr.h" #include #include static void forward(expr *node, const double *u) { - /* Constants don't depend on u; values are already set */ + /* Values are set at creation (constants) or by problem_update_params */ (void) node; (void) u; } static void jacobian_init(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. */ + /* Parameter/constant jacobian is all zeros: size x n_vars with 0 nonzeros */ node->jacobian = new_csr_matrix(node->size, node->n_vars, 0); } static void eval_jacobian(expr *node) { - /* Constant jacobian never changes - nothing to evaluate */ + /* Jacobian never changes */ (void) node; } static void wsum_hess_init(expr *node) { - /* Constant Hessian is all zeros: n_vars x n_vars with 0 nonzeros. */ + /* Hessian is all zeros */ 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; } @@ -58,12 +67,21 @@ 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 *) calloc(1, sizeof(expr)); + parameter_expr *pnode = (parameter_expr *) calloc(1, sizeof(parameter_expr)); + expr *node = &pnode->base; init_expr(node, d1, d2, n_vars, forward, jacobian_init, eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess, NULL); - memcpy(node->value, values, node->size * sizeof(double)); + pnode->param_id = param_id; + pnode->has_been_refreshed = false; + + /* If values provided (fixed constant), copy them now. + Otherwise values will be populated by problem_update_params. */ + if (values != NULL) + { + memcpy(node->value, values, node->size * sizeof(double)); + } return node; } diff --git a/src/affine/sum.c b/src/affine/sum.c index d1fec14..ce30d70 100644 --- a/src/affine/sum.c +++ b/src/affine/sum.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "affine.h" +#include "memory_wrappers.h" #include "utils/CSR_sum.h" #include "utils/int_double_pair.h" #include "utils/mini_numpy.h" @@ -175,7 +176,7 @@ static bool is_affine(const expr *node) static void free_type_data(expr *node) { sum_expr *snode = (sum_expr *) node; - free(snode->idx_map); + FREE_AND_NULL(snode->idx_map); } expr *new_sum(expr *child, int axis) diff --git a/src/affine/trace.c b/src/affine/trace.c index 1732c40..5a25bfb 100644 --- a/src/affine/trace.c +++ b/src/affine/trace.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "affine.h" +#include "memory_wrappers.h" #include "utils/CSR_sum.h" #include "utils/int_double_pair.h" #include "utils/utils.h" @@ -139,7 +140,7 @@ static void free_type_data(expr *node) if (node) { trace_expr *tnode = (trace_expr *) node; - free(tnode->idx_map); + FREE_AND_NULL(tnode->idx_map); } } diff --git a/src/bivariate/left_matmul.c b/src/bivariate/left_matmul.c index 31d8251..1cecfc1 100644 --- a/src/bivariate/left_matmul.c +++ b/src/bivariate/left_matmul.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "bivariate.h" +#include "memory_wrappers.h" #include "subexpr.h" #include "utils/Timer.h" #include "utils/linalg_sparse_matmuls.h" @@ -23,10 +24,10 @@ #include #include -/* This file implement the atom 'left_matmul' corresponding to the operation y = - A @ f(x), where A is a given matrix and f(x) is an arbitrary expression. - Here, f(x) can be a vector-valued expression and a matrix-valued - expression. The dimensions are A - m x n, f(x) - n x p, y - m x p. +/* This file implements the atom 'left_matmul' corresponding to the operation y = + A @ f(x), where A is a given matrix (from a parameter node) and f(x) is an + arbitrary expression. Here, f(x) can be a vector-valued expression and a + matrix-valued expression. The dimensions are A - m x n, f(x) - n x p, y - m x p. Note that here A does not have global column indices but it is a local matrix. This is an important distinction compared to linear_op_expr. @@ -45,24 +46,45 @@ Working in terms of A_kron unifies the implementation of f(x) being vector-valued or matrix-valued. - - I (dance858) think we can get additional big speedups when A is dense by - introducing a dense matrix class. */ #include "utils/utils.h" +#include +#include + +/* Refresh A and AT values from param_source. + A is the small m x n matrix (NOT block-diagonal). + No-op when param_source is NULL (fixed constant — values already in A). */ +static void refresh_param_values(left_matmul_expr *lin_node) +{ + parameter_expr *param = (parameter_expr *) lin_node->param_source; + + if (!param || param->has_been_refreshed) return; + param->has_been_refreshed = true; + + assert(param->param_id != PARAM_FIXED); + + /* update values of A */ + memcpy(lin_node->A->x, lin_node->param_source->value, + lin_node->A->nnz * sizeof(double)); + + /* update values of AT */ + AT_fill_values(lin_node->A, lin_node->AT, lin_node->AT_iwork); +} static void forward(expr *node, const double *u) { expr *x = node->left; + left_matmul_expr *lin_node = (left_matmul_expr *) node; + + /* possibly refresh A and AT */ + if (lin_node->refresh_param_values) lin_node->refresh_param_values(lin_node); /* child's forward pass */ node->left->forward(node->left, u); /* y = A_kron @ vec(f(x)) */ - CSR_Matrix *A = ((left_matmul_expr *) node)->A; - int n_blocks = ((left_matmul_expr *) node)->n_blocks; - block_left_multiply_vec(A, x->value, node->value, n_blocks); + block_left_multiply_vec(lin_node->A, x->value, node->value, lin_node->n_blocks); } static bool is_affine(const expr *node) @@ -77,12 +99,9 @@ static void free_type_data(expr *node) free_csr_matrix(lin_node->AT); free_csc_matrix(lin_node->Jchild_CSC); free_csc_matrix(lin_node->J_CSC); - free(lin_node->csc_to_csr_workspace); - lin_node->A = NULL; - lin_node->AT = NULL; - lin_node->Jchild_CSC = NULL; - lin_node->J_CSC = NULL; - lin_node->csc_to_csr_workspace = NULL; + FREE_AND_NULL(lin_node->csc_to_csr_workspace); + FREE_AND_NULL(lin_node->AT_iwork); + free_expr(lin_node->param_source); } static void jacobian_init(expr *node) @@ -137,7 +156,7 @@ static void wsum_hess_init(expr *node) static void eval_wsum_hess(expr *node, const double *w) { - /* compute A^T w*/ + /* compute A^T w */ CSR_Matrix *AT = ((left_matmul_expr *) node)->AT; int n_blocks = ((left_matmul_expr *) node)->n_blocks; block_left_multiply_vec(AT, w, node->dwork, n_blocks); @@ -147,20 +166,21 @@ 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) +expr *new_left_matmul(expr *param_node, expr *child, const CSR_Matrix *A) { - /* We expect u->d1 == A->n. However, numpy's broadcasting rules allow users + /* Dimension logic: handle numpy broadcasting (1, n) as (n, )/ + 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 the result of A @ u is (m, ), which is (1, m) according to broadcasting rules. We therefore check if this is the case. */ int d1, d2, n_blocks; - if (u->d1 == A->n) + if (child->d1 == A->n) { d1 = A->m; - d2 = u->d2; - n_blocks = u->d2; + d2 = child->d2; + n_blocks = child->d2; } - else if (u->d2 == A->n && u->d1 == 1) + else if (child->d2 == A->n && child->d1 == 1) { d1 = 1; d2 = A->m; @@ -168,7 +188,7 @@ expr *new_left_matmul(expr *u, const CSR_Matrix *A) } else { - fprintf(stderr, "Error in new_left_matmul: dimension mismatch \n"); + fprintf(stderr, "Error in new_left_matmul: dimension mismatch\n"); exit(1); } @@ -176,22 +196,25 @@ expr *new_left_matmul(expr *u, const CSR_Matrix *A) left_matmul_expr *lin_node = (left_matmul_expr *) calloc(1, sizeof(left_matmul_expr)); expr *node = &lin_node->base; - init_expr(node, d1, d2, u->n_vars, forward, jacobian_init, eval_jacobian, + init_expr(node, d1, d2, child->n_vars, forward, jacobian_init, eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess, free_type_data); - node->left = u; - expr_retain(u); - - /* allocate workspace. iwork is used for transposing A (requiring size A->n) - and for converting J_child csr to csc (requring size node->n_vars). - csc_to_csr_workspace is used for converting J_CSC to CSR (requring node->size) - */ - node->iwork = (int *) malloc(MAX(A->n, node->n_vars) * sizeof(int)); + node->left = child; + expr_retain(child); + + /* Store small A (NOT block-diagonal) — block functions handle the rest + Allocate workspace. iwork is used for converting J_child csr to csc + (requring size node->n_vars). csc_to_csr_workspace is used for + converting J_CSC to CSR (requring node->size) */ + node->iwork = (int *) malloc(node->n_vars * sizeof(int)); + lin_node->AT_iwork = (int *) malloc(A->n * sizeof(int)); lin_node->csc_to_csr_workspace = (int *) malloc(node->size * sizeof(int)); lin_node->n_blocks = n_blocks; - - /* store A and AT */ lin_node->A = new_csr(A); - lin_node->AT = transpose(lin_node->A, node->iwork); + lin_node->AT = transpose(lin_node->A, lin_node->AT_iwork); + lin_node->param_source = param_node; + lin_node->refresh_param_values = refresh_param_values; + + if (param_node) expr_retain(param_node); return node; } diff --git a/src/bivariate/multiply.c b/src/bivariate/multiply.c index d8836ae..dceea8d 100644 --- a/src/bivariate/multiply.c +++ b/src/bivariate/multiply.c @@ -27,7 +27,7 @@ // ------------------------------------------------------------------------------ // Implementation of elementwise multiplication when both arguments are vectors. // If one argument is a scalar variable, the broadcasting should be represented -// as a linear operator child node? How to treat if one variable is a constant? +// as a linear operator child node. // ------------------------------------------------------------------------------ static void forward(expr *node, const double *u) { diff --git a/src/bivariate/quad_over_lin.c b/src/bivariate/quad_over_lin.c index 6b781d9..4da4135 100644 --- a/src/bivariate/quad_over_lin.c +++ b/src/bivariate/quad_over_lin.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "bivariate.h" +#include "memory_wrappers.h" #include "subexpr.h" #include "utils/CSC_Matrix.h" #include @@ -102,7 +103,7 @@ static void jacobian_init(expr *node) } assert(nonzero_cols == node->jacobian->nnz); - free(col_nz); + FREE_AND_NULL(col_nz); /* insert y variable index at correct position */ insert_idx(y->var_id, node->jacobian->i, node->jacobian->nnz); diff --git a/src/bivariate/right_matmul.c b/src/bivariate/right_matmul.c index 6ec593b..166ed72 100644 --- a/src/bivariate/right_matmul.c +++ b/src/bivariate/right_matmul.c @@ -17,27 +17,53 @@ */ #include "affine.h" #include "bivariate.h" +#include "memory_wrappers.h" #include "subexpr.h" #include "utils/CSR_Matrix.h" #include "utils/linalg_sparse_matmuls.h" +#include #include +#include + +/* Refresh AT and A values from param_source for right matmul. + param_source stores values in CSR order for the original A. */ +static void refresh_param_values(left_matmul_expr *lin_node) +{ + parameter_expr *param = (parameter_expr *) lin_node->param_source; + + if (!param || param->has_been_refreshed) return; + param->has_been_refreshed = true; + + assert(param->param_id != PARAM_FIXED); + + /* update values of original A (stored in lin_node->AT) */ + memcpy(lin_node->AT->x, lin_node->param_source->value, + lin_node->AT->nnz * sizeof(double)); + + /* update values of A^T (stored in lin_node->A) */ + AT_fill_values(lin_node->AT, lin_node->A, lin_node->AT_iwork); +} /* 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. + 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) +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. */ int *work_transpose = (int *) malloc(A->n * sizeof(int)); CSR_Matrix *AT = transpose(A, work_transpose); - expr *u_transpose = new_transpose(u); - expr *left_matmul = new_left_matmul(u_transpose, AT); - expr *node = new_transpose(left_matmul); + expr *left_matmul_node = new_left_matmul(param_node, u_transpose, AT); + expr *node = new_transpose(left_matmul_node); + + /* functionality for parameter */ + left_matmul_expr *left_matmul_data = (left_matmul_expr *) left_matmul_node; + FREE_AND_NULL(left_matmul_data->AT_iwork); + left_matmul_data->AT_iwork = work_transpose; + left_matmul_data->refresh_param_values = refresh_param_values; free_csr_matrix(AT); - free(work_transpose); return node; } diff --git a/src/bivariate/const_scalar_mult.c b/src/bivariate/scalar_mult.c similarity index 80% rename from src/bivariate/const_scalar_mult.c rename to src/bivariate/scalar_mult.c index 4898389..415687e 100644 --- a/src/bivariate/const_scalar_mult.c +++ b/src/bivariate/scalar_mult.c @@ -22,7 +22,7 @@ #include #include -/* Constant scalar multiplication: y = a * child where a is a constant double */ +/* Scalar multiplication: y = a * child where a comes from a parameter node */ static void forward(expr *node, const double *u) { @@ -32,7 +32,7 @@ static void forward(expr *node, const double *u) child->forward(child, u); /* local forward pass: multiply each element by scalar a */ - double a = ((const_scalar_mult_expr *) node)->a; + double a = ((scalar_mult_expr *) node)->param_source->value[0]; for (int i = 0; i < node->size; i++) { node->value[i] = a * child->value[i]; @@ -55,7 +55,7 @@ static void jacobian_init(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); @@ -85,7 +85,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]; @@ -98,17 +98,25 @@ static bool is_affine(const expr *node) 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 *) calloc(1, sizeof(const_scalar_mult_expr)); + scalar_mult_expr *sn = (scalar_mult_expr *) node; + free_expr(sn->param_source); +} + +expr *new_scalar_mult(expr *param_node, expr *child) +{ + scalar_mult_expr *mult_node = + (scalar_mult_expr *) calloc(1, sizeof(scalar_mult_expr)); expr *node = &mult_node->base; init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init, - eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess, NULL); + eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess, + free_type_data); node->left = child; - mult_node->a = a; + mult_node->param_source = param_node; expr_retain(child); + expr_retain(param_node); return node; } diff --git a/src/bivariate/const_vector_mult.c b/src/bivariate/vector_mult.c similarity index 83% rename from src/bivariate/const_vector_mult.c rename to src/bivariate/vector_mult.c index 65823a7..6956aec 100644 --- a/src/bivariate/const_vector_mult.c +++ b/src/bivariate/vector_mult.c @@ -21,12 +21,13 @@ #include #include -/* Constant vector elementwise multiplication: y = a \circ child */ +/* Vector elementwise multiplication: y = a \circ child + * where a comes from a parameter node */ 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); @@ -54,7 +55,7 @@ static void jacobian_init(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); @@ -87,7 +88,7 @@ static void wsum_hess_init(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++) @@ -101,22 +102,22 @@ static void eval_wsum_hess(expr *node, const double *w) memcpy(node->wsum_hess->x, x->wsum_hess->x, x->wsum_hess->nnz * sizeof(double)); } -static void free_type_data(expr *node) -{ - const_vector_mult_expr *vnode = (const_vector_mult_expr *) node; - free(vnode->a); -} - 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) +static void free_type_data(expr *node) +{ + vector_mult_expr *vnode = (vector_mult_expr *) node; + free_expr(vnode->param_source); +} + +expr *new_vector_mult(expr *param_node, expr *child) { - const_vector_mult_expr *vnode = - (const_vector_mult_expr *) calloc(1, sizeof(const_vector_mult_expr)); + vector_mult_expr *vnode = + (vector_mult_expr *) calloc(1, sizeof(vector_mult_expr)); expr *node = &vnode->base; init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init, @@ -125,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 *) 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/expr.c b/src/expr.c index 99a7175..d093064 100644 --- a/src/expr.c +++ b/src/expr.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "expr.h" +#include "memory_wrappers.h" #include "utils/int_double_pair.h" #include #include @@ -61,19 +62,16 @@ void free_expr(expr *node) } /* free value array and jacobian */ - free(node->value); + FREE_AND_NULL(node->value); free_csr_matrix(node->jacobian); free_csr_matrix(node->wsum_hess); - free(node->dwork); - free(node->iwork); - node->value = NULL; + FREE_AND_NULL(node->dwork); + FREE_AND_NULL(node->iwork); node->jacobian = NULL; node->wsum_hess = NULL; - node->dwork = NULL; - node->iwork = NULL; /* free the node itself */ - free(node); + FREE_AND_NULL(node); } void expr_retain(expr *node) diff --git a/src/other/prod_axis_one.c b/src/other/prod_axis_one.c index cd43229..8a3e02b 100644 --- a/src/other/prod_axis_one.c +++ b/src/other/prod_axis_one.c @@ -1,3 +1,4 @@ +#include "memory_wrappers.h" #include "other.h" #include #include @@ -371,9 +372,9 @@ static bool is_affine(const expr *node) static void free_type_data(expr *node) { prod_axis *pnode = (prod_axis *) node; - free(pnode->num_of_zeros); - free(pnode->zero_index); - free(pnode->prod_nonzero); + FREE_AND_NULL(pnode->num_of_zeros); + FREE_AND_NULL(pnode->zero_index); + FREE_AND_NULL(pnode->prod_nonzero); } expr *new_prod_axis_one(expr *child) diff --git a/src/other/prod_axis_zero.c b/src/other/prod_axis_zero.c index ad85b6e..128e1a0 100644 --- a/src/other/prod_axis_zero.c +++ b/src/other/prod_axis_zero.c @@ -1,3 +1,4 @@ +#include "memory_wrappers.h" #include "other.h" #include #include @@ -330,9 +331,9 @@ static bool is_affine(const expr *node) static void free_type_data(expr *node) { prod_axis *pnode = (prod_axis *) node; - free(pnode->num_of_zeros); - free(pnode->zero_index); - free(pnode->prod_nonzero); + FREE_AND_NULL(pnode->num_of_zeros); + FREE_AND_NULL(pnode->zero_index); + FREE_AND_NULL(pnode->prod_nonzero); } /* TODO: refactor to remove diagonal entry as nonzero since it's always zero */ diff --git a/src/other/quad_form.c b/src/other/quad_form.c index d03c70e..fb7a97f 100644 --- a/src/other/quad_form.c +++ b/src/other/quad_form.c @@ -1,3 +1,4 @@ +#include "memory_wrappers.h" #include "other.h" #include "subexpr.h" #include "utils/CSC_Matrix.h" @@ -129,7 +130,7 @@ static void jacobian_init(expr *node) } } assert(nonzero_cols == node->jacobian->nnz); - free(col_nz); + FREE_AND_NULL(col_nz); node->jacobian->p[0] = 0; node->jacobian->p[1] = node->jacobian->nnz; @@ -174,7 +175,6 @@ static void free_type_data(expr *node) { quad_form_expr *qnode = (quad_form_expr *) node; free_csr_matrix(qnode->Q); - qnode->Q = NULL; } static bool is_affine(const expr *node) diff --git a/src/problem.c b/src/problem.c index 3413858..294db92 100644 --- a/src/problem.c +++ b/src/problem.c @@ -16,6 +16,8 @@ * limitations under the License. */ #include "problem.h" +#include "memory_wrappers.h" +#include "subexpr.h" #include "utils/CSR_sum.h" #include "utils/utils.h" #include @@ -66,6 +68,9 @@ problem *new_problem(expr *objective, expr **constraints, int n_constraints, prob->stats.nnz_affine = 0; prob->stats.nnz_nonlinear = 0; + prob->param_nodes = NULL; + prob->n_param_nodes = 0; + prob->verbose = verbose; return prob; @@ -231,7 +236,7 @@ void problem_init_hessian(problem *prob) prob->hess_idx_map = (int *) malloc(nnz * sizeof(int)); int *iwork = (int *) malloc(MAX(nnz, prob->n_vars) * sizeof(int)); problem_lagrange_hess_fill_sparsity(prob, iwork); - free(iwork); + FREE_AND_NULL(iwork); clock_gettime(CLOCK_MONOTONIC, &timer.end); prob->stats.time_init_derivatives += GET_ELAPSED_SECONDS(timer); @@ -282,11 +287,14 @@ void free_problem(problem *prob) if (prob == NULL) return; /* Free allocated arrays */ - free(prob->constraint_values); - free(prob->gradient_values); + FREE_AND_NULL(prob->constraint_values); + FREE_AND_NULL(prob->gradient_values); free_csr_matrix(prob->jacobian); free_csr_matrix(prob->lagrange_hessian); - free(prob->hess_idx_map); + FREE_AND_NULL(prob->hess_idx_map); + + /* Free parameter node array (weak references, not owned) */ + FREE_AND_NULL(prob->param_nodes); /* Release expression references (decrements refcount) */ free_expr(prob->objective); @@ -294,7 +302,7 @@ void free_problem(problem *prob) { free_expr(prob->constraints[i]); } - free(prob->constraints); + FREE_AND_NULL(prob->constraints); if (prob->verbose) { @@ -302,7 +310,7 @@ void free_problem(problem *prob) } /* Free problem struct */ - free(prob); + FREE_AND_NULL(prob); } double problem_objective_forward(problem *prob, const double *u) @@ -440,3 +448,30 @@ void problem_hessian(problem *prob, double obj_w, const double *w) clock_gettime(CLOCK_MONOTONIC, &timer.end); prob->stats.time_eval_hessian += GET_ELAPSED_SECONDS(timer); } + +void problem_register_params(problem *prob, expr **param_nodes, int n_param_nodes) +{ + prob->n_param_nodes = n_param_nodes; + prob->param_nodes = (expr **) 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++) + { + parameter_expr *p = (parameter_expr *) prob->param_nodes[i]; + if (p->param_id == PARAM_FIXED) continue; + memcpy(p->base.value, theta + p->param_id, p->base.size * sizeof(double)); + p->has_been_refreshed = false; + } + + /* force re-evaluation of affine Jacobians on next call */ + prob->jacobian_called = false; +} diff --git a/src/utils/CSC_Matrix.c b/src/utils/CSC_Matrix.c index d017a62..c445917 100644 --- a/src/utils/CSC_Matrix.c +++ b/src/utils/CSC_Matrix.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "utils/CSC_Matrix.h" +#include "memory_wrappers.h" #include "utils/iVec.h" #include #include @@ -32,10 +33,10 @@ CSC_Matrix *new_csc_matrix(int m, int n, int nnz) if (!matrix->p || !matrix->i || !matrix->x) { - free(matrix->p); - free(matrix->i); - free(matrix->x); - free(matrix); + FREE_AND_NULL(matrix->p); + FREE_AND_NULL(matrix->i); + FREE_AND_NULL(matrix->x); + FREE_AND_NULL(matrix); return NULL; } @@ -50,10 +51,10 @@ void free_csc_matrix(CSC_Matrix *matrix) { if (matrix) { - free(matrix->p); - free(matrix->i); - free(matrix->x); - free(matrix); + FREE_AND_NULL(matrix->p); + FREE_AND_NULL(matrix->i); + FREE_AND_NULL(matrix->x); + FREE_AND_NULL(matrix); } } @@ -105,7 +106,7 @@ CSR_Matrix *ATA_alloc(const CSC_Matrix *A) symmetrize_csr(Cp, Ci->data, n, C); /* free workspace */ - free(Cp); + FREE_AND_NULL(Cp); iVec_free(Ci); return C; @@ -211,7 +212,7 @@ CSC_Matrix *csr_to_csc(const CSR_Matrix *A) } } - free(count); + FREE_AND_NULL(count); return C; } @@ -395,7 +396,7 @@ CSR_Matrix *BTA_alloc(const CSC_Matrix *A, const CSC_Matrix *B) memcpy(C->i, Ci->data, nnz * sizeof(int)); /* free workspace */ - free(Cp); + FREE_AND_NULL(Cp); iVec_free(Ci); return C; diff --git a/src/utils/CSR_Matrix.c b/src/utils/CSR_Matrix.c index 513a457..0778804 100644 --- a/src/utils/CSR_Matrix.c +++ b/src/utils/CSR_Matrix.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "utils/CSR_Matrix.h" +#include "memory_wrappers.h" #include "utils/int_double_pair.h" #include "utils/utils.h" #include @@ -49,10 +50,10 @@ void free_csr_matrix(CSR_Matrix *matrix) { if (matrix) { - free(matrix->p); - free(matrix->i); - free(matrix->x); - free(matrix); + FREE_AND_NULL(matrix->p); + FREE_AND_NULL(matrix->i); + FREE_AND_NULL(matrix->x); + FREE_AND_NULL(matrix); } } @@ -448,5 +449,5 @@ void symmetrize_csr(const int *Ap, const int *Ai, int m, CSR_Matrix *C) } } - free(counts); + FREE_AND_NULL(counts); } diff --git a/src/utils/int_double_pair.c b/src/utils/int_double_pair.c index 6b49021..1c0d01b 100644 --- a/src/utils/int_double_pair.c +++ b/src/utils/int_double_pair.c @@ -16,6 +16,7 @@ * limitations under the License. */ #include "utils/int_double_pair.h" +#include "memory_wrappers.h" #include static int compare_int_double_pair(const void *a, const void *b) @@ -45,7 +46,7 @@ void set_int_double_pair_array(int_double_pair *pair, int *ints, double *doubles void free_int_double_pair_array(int_double_pair *array) { - free(array); + FREE_AND_NULL(array); } void sort_int_double_pair_array(int_double_pair *array, int size) diff --git a/src/utils/linalg_sparse_matmuls.c b/src/utils/linalg_sparse_matmuls.c index 82c03e0..0bbae7e 100644 --- a/src/utils/linalg_sparse_matmuls.c +++ b/src/utils/linalg_sparse_matmuls.c @@ -15,6 +15,7 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include "memory_wrappers.h" #include "utils/CSC_Matrix.h" #include "utils/CSR_Matrix.h" #include "utils/iVec.h" @@ -183,7 +184,7 @@ CSC_Matrix *block_left_multiply_fill_sparsity(const CSR_Matrix *A, memcpy(C->i, Ci->data, Ci->len * sizeof(int)); /* Clean up workspace */ - free(Cp); + FREE_AND_NULL(Cp); iVec_free(Ci); return C; @@ -312,7 +313,7 @@ CSR_Matrix *csr_csc_matmul_alloc(const CSR_Matrix *A, const CSC_Matrix *B) CSR_Matrix *C = new_csr_matrix(m, p, nnz); memcpy(C->p, Cp, (m + 1) * sizeof(int)); memcpy(C->i, Ci->data, nnz * sizeof(int)); - free(Cp); + FREE_AND_NULL(Cp); iVec_free(Ci); return C; diff --git a/tests/all_tests.c b/tests/all_tests.c index 6aa699a..e652965 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -11,7 +11,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/composite/test_composite.h" #include "forward_pass/elementwise/test_exp.h" #include "forward_pass/elementwise/test_log.h" @@ -20,8 +20,6 @@ #include "forward_pass/test_prod_axis_zero.h" #include "jacobian_tests/test_broadcast.h" #include "jacobian_tests/test_composite.h" -#include "jacobian_tests/test_const_scalar_mult.h" -#include "jacobian_tests/test_const_vector_mult.h" #include "jacobian_tests/test_elementwise_mult.h" #include "jacobian_tests/test_hstack.h" #include "jacobian_tests/test_index.h" @@ -39,9 +37,12 @@ #include "jacobian_tests/test_rel_entr_scalar_vector.h" #include "jacobian_tests/test_rel_entr_vector_scalar.h" #include "jacobian_tests/test_right_matmul.h" +#include "jacobian_tests/test_scalar_mult.h" #include "jacobian_tests/test_sum.h" #include "jacobian_tests/test_trace.h" #include "jacobian_tests/test_transpose.h" +#include "jacobian_tests/test_vector_mult.h" +#include "problem/test_param_prob.h" #include "problem/test_problem.h" #include "utils/test_csc_matrix.h" #include "utils/test_csr_csc_conversion.h" @@ -56,8 +57,6 @@ #include "wsum_hess/elementwise/test_trig.h" #include "wsum_hess/elementwise/test_xexp.h" #include "wsum_hess/test_broadcast.h" -#include "wsum_hess/test_const_scalar_mult.h" -#include "wsum_hess/test_const_vector_mult.h" #include "wsum_hess/test_hstack.h" #include "wsum_hess/test_index.h" #include "wsum_hess/test_left_matmul.h" @@ -72,9 +71,11 @@ #include "wsum_hess/test_rel_entr_scalar_vector.h" #include "wsum_hess/test_rel_entr_vector_scalar.h" #include "wsum_hess/test_right_matmul.h" +#include "wsum_hess/test_scalar_mult.h" #include "wsum_hess/test_sum.h" #include "wsum_hess/test_trace.h" #include "wsum_hess/test_transpose.h" +#include "wsum_hess/test_vector_mult.h" #endif /* PROFILE_ONLY */ #ifdef PROFILE_ONLY @@ -90,7 +91,7 @@ int main(void) #ifndef PROFILE_ONLY printf("--- Forward Pass Tests ---\n"); mu_run_test(test_variable, tests_run); - mu_run_test(test_constant, tests_run); + mu_run_test(test_fixed_parameter, tests_run); mu_run_test(test_addition, tests_run); mu_run_test(test_linear_op, tests_run); mu_run_test(test_neg_forward, tests_run); @@ -280,6 +281,10 @@ 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); + 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); #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 7ce859c..a087c9a 100644 --- a/tests/forward_pass/affine/test_add.h +++ b/tests/forward_pass/affine/test_add.h @@ -5,6 +5,7 @@ #include "affine.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_addition() @@ -12,7 +13,7 @@ const char *test_addition() 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_sum.h b/tests/forward_pass/affine/test_sum.h index 2e3d6ab..e700ebd 100644 --- a/tests/forward_pass/affine/test_sum.h +++ b/tests/forward_pass/affine/test_sum.h @@ -6,6 +6,7 @@ #include "elementwise_univariate.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_sum_axis_neg1() @@ -17,7 +18,7 @@ const char *test_sum_axis_neg1() 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); @@ -42,7 +43,7 @@ const char *test_sum_axis_0() 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); @@ -69,7 +70,7 @@ const char *test_sum_axis_1() 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 70% rename from tests/forward_pass/affine/test_variable_constant.h rename to tests/forward_pass/affine/test_variable_parameter.h index c964654..2cfea7b 100644 --- a/tests/forward_pass/affine/test_variable_constant.h +++ b/tests/forward_pass/affine/test_variable_parameter.h @@ -5,6 +5,7 @@ #include "affine.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_variable() @@ -17,13 +18,14 @@ const char *test_variable() return 0; } -const char *test_constant() +const char *test_fixed_parameter() { 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)); + mu_assert("Fixed parameter test failed", + cmp_double_array(const_node->value, c, 2)); free_expr(const_node); return 0; } diff --git a/tests/forward_pass/composite/test_composite.h b/tests/forward_pass/composite/test_composite.h index 92074a0..253aa6a 100644 --- a/tests/forward_pass/composite/test_composite.h +++ b/tests/forward_pass/composite/test_composite.h @@ -6,6 +6,7 @@ #include "elementwise_univariate.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_composite() @@ -16,7 +17,7 @@ const char *test_composite() /* 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/test_prod_axis_one.h b/tests/forward_pass/test_prod_axis_one.h index 6f4b3bb..499214e 100644 --- a/tests/forward_pass/test_prod_axis_one.h +++ b/tests/forward_pass/test_prod_axis_one.h @@ -6,6 +6,7 @@ #include "expr.h" #include "minunit.h" #include "other.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_forward_prod_axis_one() @@ -16,7 +17,7 @@ const char *test_forward_prod_axis_one() 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/test_prod_axis_zero.h b/tests/forward_pass/test_prod_axis_zero.h index 6504502..30b8cdc 100644 --- a/tests/forward_pass/test_prod_axis_zero.h +++ b/tests/forward_pass/test_prod_axis_zero.h @@ -6,6 +6,7 @@ #include "expr.h" #include "minunit.h" #include "other.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_forward_prod_axis_zero() @@ -16,7 +17,7 @@ const char *test_forward_prod_axis_zero() 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/test_broadcast.h b/tests/jacobian_tests/test_broadcast.h index 612e5cf..da99597 100644 --- a/tests/jacobian_tests/test_broadcast.h +++ b/tests/jacobian_tests/test_broadcast.h @@ -5,6 +5,7 @@ #include "affine.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_broadcast_row_jacobian() @@ -141,7 +142,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/test_left_matmul.h b/tests/jacobian_tests/test_left_matmul.h index f79eee7..feda189 100644 --- a/tests/jacobian_tests/test_left_matmul.h +++ b/tests/jacobian_tests/test_left_matmul.h @@ -1,10 +1,12 @@ #include #include +#include #include "bivariate.h" #include "elementwise_univariate.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_jacobian_left_matmul_log() @@ -30,7 +32,7 @@ const char *test_jacobian_left_matmul_log() double x_vals[3] = {1.0, 2.0, 3.0}; expr *x = new_variable(3, 1, 0, 3); - /* Create sparse matrix A in CSR format */ + /* A is 4x3: [1, 0, 2; 3, 0, 4; 5, 0, 6; 7, 0, 0] */ CSR_Matrix *A = new_csr_matrix(4, 3, 7); int A_p[5] = {0, 2, 4, 6, 7}; int A_i[7] = {0, 2, 0, 2, 0, 2, 0}; @@ -40,7 +42,8 @@ const char *test_jacobian_left_matmul_log() 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); + free_csr_matrix(A); A_log_x->forward(A_log_x, x_vals); A_log_x->jacobian_init(A_log_x); @@ -63,7 +66,6 @@ const char *test_jacobian_left_matmul_log() mu_assert("cols fail", cmp_int_array(A_log_x->jacobian->i, expected_Ai, 7)); mu_assert("rows fail", cmp_int_array(A_log_x->jacobian->p, expected_Ap, 5)); - free_csr_matrix(A); free_expr(A_log_x); return 0; } @@ -74,7 +76,7 @@ const char *test_jacobian_left_matmul_log_matrix() double x_vals[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; expr *x = new_variable(3, 2, 0, 6); - /* Create sparse matrix A in CSR format (4x3) */ + /* A is 4x3: [1, 0, 2; 3, 0, 4; 5, 0, 6; 7, 0, 0] */ CSR_Matrix *A = new_csr_matrix(4, 3, 7); int A_p[5] = {0, 2, 4, 6, 7}; int A_i[7] = {0, 2, 0, 2, 0, 2, 0}; @@ -84,7 +86,8 @@ const char *test_jacobian_left_matmul_log_matrix() 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); + free_csr_matrix(A); A_log_x->forward(A_log_x, x_vals); A_log_x->jacobian_init(A_log_x); @@ -102,7 +105,6 @@ const char *test_jacobian_left_matmul_log_matrix() mu_assert("cols fail", cmp_int_array(A_log_x->jacobian->i, expected_Ai, 14)); mu_assert("rows fail", cmp_int_array(A_log_x->jacobian->p, expected_Ap, 9)); - free_csr_matrix(A); free_expr(A_log_x); return 0; } @@ -145,7 +147,7 @@ const char *test_jacobian_left_matmul_log_composite() memcpy(B->i, B_i, 9 * sizeof(int)); memcpy(B->x, B_x, 9 * sizeof(double)); - /* Create A matrix */ + /* A is 4x3: [1, 0, 2; 3, 0, 4; 5, 0, 6; 7, 0, 0] */ CSR_Matrix *A = new_csr_matrix(4, 3, 7); int A_p[5] = {0, 2, 4, 6, 7}; int A_i[7] = {0, 2, 0, 2, 0, 2, 0}; @@ -156,7 +158,8 @@ const char *test_jacobian_left_matmul_log_composite() expr *Bx = new_linear(x, B, NULL); expr *log_Bx = new_log(Bx); - expr *A_log_Bx = new_left_matmul(log_Bx, A); + expr *A_log_Bx = new_left_matmul(NULL, log_Bx, A); + free_csr_matrix(A); A_log_Bx->forward(A_log_Bx, x_vals); A_log_Bx->jacobian_init(A_log_Bx); @@ -176,7 +179,6 @@ const char *test_jacobian_left_matmul_log_composite() mu_assert("cols fail", cmp_int_array(A_log_Bx->jacobian->i, expected_Ai, 12)); mu_assert("rows fail", cmp_int_array(A_log_Bx->jacobian->p, expected_Ap, 5)); - free_csr_matrix(A); free_csr_matrix(B); free_expr(A_log_Bx); return 0; diff --git a/tests/jacobian_tests/test_right_matmul.h b/tests/jacobian_tests/test_right_matmul.h index 931df44..a181d60 100644 --- a/tests/jacobian_tests/test_right_matmul.h +++ b/tests/jacobian_tests/test_right_matmul.h @@ -27,7 +27,7 @@ const char *test_jacobian_right_matmul_log() 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); log_x_A->jacobian_init(log_x_A); @@ -76,7 +76,7 @@ const char *test_jacobian_right_matmul_log_vector() 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); log_x_A->jacobian_init(log_x_A); diff --git a/tests/jacobian_tests/test_const_scalar_mult.h b/tests/jacobian_tests/test_scalar_mult.h similarity index 90% rename from tests/jacobian_tests/test_const_scalar_mult.h rename to tests/jacobian_tests/test_scalar_mult.h index 3143929..6e1cfd0 100644 --- a/tests/jacobian_tests/test_const_scalar_mult.h +++ b/tests/jacobian_tests/test_scalar_mult.h @@ -1,9 +1,11 @@ #include +#include "affine.h" #include "bivariate.h" #include "elementwise_univariate.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 */ @@ -19,7 +21,8 @@ const char *test_jacobian_const_scalar_mult_log_vector() /* Create scalar mult node: y = 2.5 * log(x) */ double a = 2.5; - expr *y = new_const_scalar_mult(a, log_node); + expr *a_node = new_parameter(1, 1, PARAM_FIXED, 3, &a); + expr *y = new_scalar_mult(a_node, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -55,7 +58,8 @@ const char *test_jacobian_const_scalar_mult_log_matrix() /* Create scalar mult node: y = 3.0 * log(x) */ double a = 3.0; - expr *y = new_const_scalar_mult(a, log_node); + expr *a_node = new_parameter(1, 1, PARAM_FIXED, 4, &a); + expr *y = new_scalar_mult(a_node, log_node); /* Forward pass */ y->forward(y, u_vals); diff --git a/tests/jacobian_tests/test_transpose.h b/tests/jacobian_tests/test_transpose.h index 4fd22d5..0870d20 100644 --- a/tests/jacobian_tests/test_transpose.h +++ b/tests/jacobian_tests/test_transpose.h @@ -3,25 +3,29 @@ #define TEST_TRANSPOSE_H #include "affine.h" +#include "bivariate.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" #include #include +#include const char *test_jacobian_transpose() { - // A = [1 2; 3 4] + /* A = [1 2; 3 4] as dense 2x2 CSR */ CSR_Matrix *A = new_csr_matrix(2, 2, 4); - int A_p[3] = {0, 2, 4}; - int A_i[4] = {0, 1, 0, 1}; - double A_x[4] = {1, 2, 3, 4}; - memcpy(A->p, A_p, 3 * sizeof(int)); - memcpy(A->i, A_i, 4 * sizeof(int)); - memcpy(A->x, A_x, 4 * sizeof(double)); + int Ap[3] = {0, 2, 4}; + int Ai[4] = {0, 1, 0, 1}; + double Ax[4] = {1.0, 2.0, 3.0, 4.0}; + memcpy(A->p, Ap, 3 * sizeof(int)); + memcpy(A->i, Ai, 4 * sizeof(int)); + memcpy(A->x, Ax, 4 * sizeof(double)); // 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); + free_csr_matrix(A); expr *transpose_AX = new_transpose(AX); double u[4] = {1, 3, 2, 4}; transpose_AX->forward(transpose_AX, u); @@ -40,7 +44,6 @@ const char *test_jacobian_transpose() mu_assert("jacobian col idx fail", cmp_int_array(transpose_AX->jacobian->i, expected_i, 8)); free_expr(transpose_AX); - free_csr_matrix(A); return 0; } diff --git a/tests/jacobian_tests/test_const_vector_mult.h b/tests/jacobian_tests/test_vector_mult.h similarity index 91% rename from tests/jacobian_tests/test_const_vector_mult.h rename to tests/jacobian_tests/test_vector_mult.h index 4658dd4..3122e2e 100644 --- a/tests/jacobian_tests/test_const_vector_mult.h +++ b/tests/jacobian_tests/test_vector_mult.h @@ -1,9 +1,11 @@ #include +#include "affine.h" #include "bivariate.h" #include "elementwise_univariate.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 */ @@ -19,7 +21,8 @@ const char *test_jacobian_const_vector_mult_log_vector() /* 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_node = new_parameter(3, 1, PARAM_FIXED, 3, a); + expr *y = new_vector_mult(a_node, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -59,7 +62,8 @@ const char *test_jacobian_const_vector_mult_log_matrix() /* 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_node = new_parameter(2, 2, PARAM_FIXED, 4, a); + expr *y = new_vector_mult(a_node, log_node); /* Forward pass */ y->forward(y, u_vals); diff --git a/tests/problem/test_param_prob.h b/tests/problem/test_param_prob.h new file mode 100644 index 0000000..29b9c1b --- /dev/null +++ b/tests/problem/test_param_prob.h @@ -0,0 +1,351 @@ +#ifndef TEST_PARAM_PROB_H +#define TEST_PARAM_PROB_H + +#include +#include +#include + +#include "affine.h" +#include "bivariate.h" +#include "elementwise_univariate.h" +#include "expr.h" +#include "minunit.h" +#include "problem.h" +#include "test_helpers.h" + +/* + * Test 1: param_scalar_mult in objective + * + * Problem: minimize a * sum(log(x)), no constraints, x size 2 + * a is a scalar parameter (param_id=0) + * + * At x=[1,2], a=3: + * obj = 3*(log(1)+log(2)) = 3*log(2) + * gradient = [3/1, 3/2] = [3.0, 1.5] + * + * After update a=5: + * obj = 5*log(2) + * gradient = [5.0, 2.5] + */ +const char *test_param_scalar_mult_problem(void) +{ + int n_vars = 2; + + /* Build tree: sum(a * log(x)) */ + expr *x = new_variable(2, 1, 0, n_vars); + expr *log_x = new_log(x); + expr *a_param = new_parameter(1, 1, 0, n_vars, NULL); + expr *scaled = new_scalar_mult(a_param, log_x); + expr *objective = new_sum(scaled, -1); + + /* Create problem (no constraints) */ + problem *prob = new_problem(objective, NULL, 0, true); + + /* Register parameter */ + expr *param_nodes[1] = {a_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* Set a=3 and evaluate at x=[1,2] */ + double theta[1] = {3.0}; + problem_update_params(prob, theta); + + double u[2] = {1.0, 2.0}; + double obj_val = problem_objective_forward(prob, u); + problem_gradient(prob); + + double expected_obj = 3.0 * log(2.0); + mu_assert("obj wrong (a=3)", fabs(obj_val - expected_obj) < 1e-10); + + double expected_grad[2] = {3.0, 1.5}; + mu_assert("gradient wrong (a=3)", + cmp_double_array(prob->gradient_values, expected_grad, 2)); + + /* Update a=5 and re-evaluate */ + theta[0] = 5.0; + problem_update_params(prob, theta); + + obj_val = problem_objective_forward(prob, u); + problem_gradient(prob); + + expected_obj = 5.0 * log(2.0); + mu_assert("obj wrong (a=5)", fabs(obj_val - expected_obj) < 1e-10); + + double expected_grad2[2] = {5.0, 2.5}; + mu_assert("gradient wrong (a=5)", + cmp_double_array(prob->gradient_values, expected_grad2, 2)); + + free_problem(prob); + + return 0; +} + +/* + * Test 2: param_vector_mult in constraint + * + * Problem: minimize sum(x), subject to p ∘ x, x size 2 + * p is a vector parameter of size 2 (param_id=0) + * + * At x=[1,2], p=[3,4]: + * constraint_values = [3, 8] + * jacobian = diag([3, 4]) + * + * After update p=[5,6]: + * constraint_values = [5, 12] + * jacobian = diag([5, 6]) + */ +const char *test_param_vector_mult_problem(void) +{ + int n_vars = 2; + + /* Objective: sum(x) */ + expr *x_obj = new_variable(2, 1, 0, n_vars); + expr *objective = new_sum(x_obj, -1); + + /* Constraint: p ∘ x */ + expr *x_con = new_variable(2, 1, 0, n_vars); + expr *p_param = new_parameter(2, 1, 0, n_vars, NULL); + expr *constraint = new_vector_mult(p_param, x_con); + + expr *constraints[1] = {constraint}; + + /* Create problem */ + problem *prob = new_problem(objective, constraints, 1, true); + + expr *param_nodes[1] = {p_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* Set p=[3,4] and evaluate at x=[1,2] */ + double theta[2] = {3.0, 4.0}; + problem_update_params(prob, theta); + + double u[2] = {1.0, 2.0}; + problem_constraint_forward(prob, u); + problem_jacobian(prob); + + double expected_cv[2] = {3.0, 8.0}; + mu_assert("constraint values wrong (p=[3,4])", + cmp_double_array(prob->constraint_values, expected_cv, 2)); + + CSR_Matrix *jac = prob->jacobian; + mu_assert("jac rows wrong", jac->m == 2); + mu_assert("jac cols wrong", jac->n == 2); + + int expected_p[3] = {0, 1, 2}; + mu_assert("jac->p wrong (p=[3,4])", cmp_int_array(jac->p, expected_p, 3)); + + int expected_i[2] = {0, 1}; + mu_assert("jac->i wrong (p=[3,4])", cmp_int_array(jac->i, expected_i, 2)); + + double expected_x[2] = {3.0, 4.0}; + mu_assert("jac->x wrong (p=[3,4])", cmp_double_array(jac->x, expected_x, 2)); + + /* Update p=[5,6] and re-evaluate */ + double theta2[2] = {5.0, 6.0}; + problem_update_params(prob, theta2); + + problem_constraint_forward(prob, u); + problem_jacobian(prob); + + double expected_cv2[2] = {5.0, 12.0}; + mu_assert("constraint values wrong (p=[5,6])", + cmp_double_array(prob->constraint_values, expected_cv2, 2)); + + double expected_x2[2] = {5.0, 6.0}; + mu_assert("jac->x wrong (p=[5,6])", cmp_double_array(jac->x, expected_x2, 2)); + + free_problem(prob); + + return 0; +} + +/* + * Test 3: left_param_matmul in constraint + * + * Problem: minimize sum(x), subject to A @ x, x size 2, A is 2x2 + * A is a 2x2 matrix parameter (param_id=0, size=4, CSR data order) + * A = [[1,2],[3,4]] → CSR data order theta = [1,2,3,4] + * + * At x=[1,2]: + * constraint_values = [1*1+2*2, 3*1+4*2] = [5, 11] + * jacobian = [[1,2],[3,4]] + * + * After update A = [[5,6],[7,8]] → theta = [5,6,7,8]: + * constraint_values = [5*1+6*2, 7*1+8*2] = [17, 23] + * jacobian = [[5,6],[7,8]] + */ +const char *test_param_left_matmul_problem(void) +{ + int n_vars = 2; + + /* Objective: sum(x) */ + expr *x_obj = new_variable(2, 1, 0, n_vars); + expr *objective = new_sum(x_obj, -1); + + /* Constraint: A @ x */ + expr *x_con = new_variable(2, 1, 0, n_vars); + expr *A_param = new_parameter(2, 2, 0, n_vars, NULL); + + /* Dense 2x2 CSR with placeholder zeros (values refreshed from A_param) */ + CSR_Matrix *A = new_csr_matrix(2, 2, 4); + int Ap[3] = {0, 2, 4}; + int Ai[4] = {0, 1, 0, 1}; + double Ax[4] = {0.0, 0.0, 0.0, 0.0}; + memcpy(A->p, Ap, 3 * sizeof(int)); + memcpy(A->i, Ai, 4 * sizeof(int)); + memcpy(A->x, Ax, 4 * sizeof(double)); + + expr *constraint = new_left_matmul(A_param, x_con, A); + free_csr_matrix(A); + + expr *constraints[1] = {constraint}; + + /* Create problem */ + problem *prob = new_problem(objective, constraints, 1, true); + + expr *param_nodes[1] = {A_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* Set A = [[1,2],[3,4]], CSR data order: [1,2,3,4] */ + double theta[4] = {1.0, 2.0, 3.0, 4.0}; + problem_update_params(prob, theta); + + double u[2] = {1.0, 2.0}; + problem_constraint_forward(prob, u); + problem_jacobian(prob); + + double expected_cv[2] = {5.0, 11.0}; + mu_assert("constraint values wrong (A1)", + cmp_double_array(prob->constraint_values, expected_cv, 2)); + + CSR_Matrix *jac = prob->jacobian; + mu_assert("jac rows wrong", jac->m == 2); + mu_assert("jac cols wrong", jac->n == 2); + + /* Dense jacobian = [[1,2],[3,4]], CSR: row 0 → cols 0,1 vals 1,2; + * row 1 → cols 0,1 vals 3,4 */ + int expected_p[3] = {0, 2, 4}; + mu_assert("jac->p wrong (A1)", cmp_int_array(jac->p, expected_p, 3)); + + int expected_i[4] = {0, 1, 0, 1}; + mu_assert("jac->i wrong (A1)", cmp_int_array(jac->i, expected_i, 4)); + + double expected_x[4] = {1.0, 2.0, 3.0, 4.0}; + mu_assert("jac->x wrong (A1)", cmp_double_array(jac->x, expected_x, 4)); + + /* Update A = [[5,6],[7,8]], CSR data order: [5,6,7,8] */ + double theta2[4] = {5.0, 6.0, 7.0, 8.0}; + problem_update_params(prob, theta2); + + problem_constraint_forward(prob, u); + problem_jacobian(prob); + + double expected_cv2[2] = {17.0, 23.0}; + mu_assert("constraint values wrong (A2)", + cmp_double_array(prob->constraint_values, expected_cv2, 2)); + + double expected_x2[4] = {5.0, 6.0, 7.0, 8.0}; + mu_assert("jac->x wrong (A2)", cmp_double_array(jac->x, expected_x2, 4)); + + free_problem(prob); + + return 0; +} + +/* + * Test 4: right_param_matmul in constraint + * + * Problem: minimize sum(x), subject to x @ A, x size 1x2, A is 2x2 + * A is a 2x2 matrix parameter (param_id=0, size=4, CSR data order) + * A = [[1,2],[3,4]] → CSR data order theta = [1,2,3,4] + * + * At x=[1,2]: + * constraint_values = [1*1+2*3, 1*2+2*4] = [7, 10] + * jacobian = [[1,3],[2,4]] = A^T + * + * After update A = [[5,6],[7,8]] → theta = [5,6,7,8]: + * constraint_values = [1*5+2*7, 1*6+2*8] = [19, 22] + * jacobian = [[5,7],[6,8]] = A^T + */ +const char *test_param_right_matmul_problem(void) +{ + int n_vars = 2; + + /* Objective: sum(x) */ + expr *x_obj = new_variable(1, 2, 0, n_vars); + expr *objective = new_sum(x_obj, -1); + + /* Constraint: x @ A */ + expr *x_con = new_variable(1, 2, 0, n_vars); + expr *A_param = new_parameter(2, 2, 0, n_vars, NULL); + + /* Dense 2x2 CSR with placeholder zeros (values refreshed from A_param) */ + CSR_Matrix *A = new_csr_matrix(2, 2, 4); + int Ap[3] = {0, 2, 4}; + int Ai[4] = {0, 1, 0, 1}; + double Ax[4] = {0.0, 0.0, 0.0, 0.0}; + memcpy(A->p, Ap, 3 * sizeof(int)); + memcpy(A->i, Ai, 4 * sizeof(int)); + memcpy(A->x, Ax, 4 * sizeof(double)); + + expr *constraint = new_right_matmul(A_param, x_con, A); + free_csr_matrix(A); + + expr *constraints[1] = {constraint}; + + /* Create problem */ + problem *prob = new_problem(objective, constraints, 1, true); + + expr *param_nodes[1] = {A_param}; + problem_register_params(prob, param_nodes, 1); + problem_init_derivatives(prob); + + /* Set A = [[1,2],[3,4]], CSR data order: [1,2,3,4] */ + double theta[4] = {1.0, 2.0, 3.0, 4.0}; + problem_update_params(prob, theta); + + double u[2] = {1.0, 2.0}; + problem_constraint_forward(prob, u); + problem_jacobian(prob); + + double expected_cv[2] = {7.0, 10.0}; + mu_assert("constraint values wrong (A1)", + cmp_double_array(prob->constraint_values, expected_cv, 2)); + + CSR_Matrix *jac = prob->jacobian; + mu_assert("jac rows wrong", jac->m == 2); + mu_assert("jac cols wrong", jac->n == 2); + + /* Dense jacobian = [[1,3],[2,4]] = A^T, CSR: row 0 → cols 0,1 vals 1,3; + * row 1 → cols 0,1 vals 2,4 */ + int expected_p[3] = {0, 2, 4}; + mu_assert("jac->p wrong (A1)", cmp_int_array(jac->p, expected_p, 3)); + + int expected_i[4] = {0, 1, 0, 1}; + mu_assert("jac->i wrong (A1)", cmp_int_array(jac->i, expected_i, 4)); + + double expected_x[4] = {1.0, 3.0, 2.0, 4.0}; + mu_assert("jac->x wrong (A1)", cmp_double_array(jac->x, expected_x, 4)); + + /* Update A = [[5,6],[7,8]], CSR data order: [5,6,7,8] */ + double theta2[4] = {5.0, 6.0, 7.0, 8.0}; + problem_update_params(prob, theta2); + + problem_constraint_forward(prob, u); + problem_jacobian(prob); + + double expected_cv2[2] = {19.0, 22.0}; + mu_assert("constraint values wrong (A2)", + cmp_double_array(prob->constraint_values, expected_cv2, 2)); + + double expected_x2[4] = {5.0, 7.0, 6.0, 8.0}; + mu_assert("jac->x wrong (A2)", cmp_double_array(jac->x, expected_x2, 4)); + + 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 a8e819b..38a4bb3 100644 --- a/tests/profiling/profile_left_matmul.h +++ b/tests/profiling/profile_left_matmul.h @@ -7,31 +7,38 @@ #include "bivariate.h" #include "elementwise_univariate.h" #include "expr.h" +#include "memory_wrappers.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" #include "utils/Timer.h" const char *profile_left_matmul() { - /* A @ X where A is 50 x 50 dense stored in CSR and X is 50 x 50 variable */ + /* A @ X where A is 100 x 100 dense (all ones) and X is 100 x 100 variable */ int n = 100; expr *X = new_variable(n, n, 0, n * n); - CSR_Matrix *A = new_csr_matrix(n, n, n * n); - for (int i = 0; i < n * n; i++) - { - A->x[i] = 1.0; /* dense matrix of all ones */ - } - for (int row = 0; row < n; row++) + + /* Build dense n x n CSR (all ones) */ + int nnz = n * n; + CSR_Matrix *A = new_csr_matrix(n, n, nnz); { - A->p[row] = row * n; - for (int col = 0; col < n; col++) + int idx = 0; + for (int row = 0; row < n; row++) { - A->i[row * n + col] = col; + A->p[row] = idx; + for (int col = 0; col < n; col++) + { + A->i[idx] = col; + A->x[idx] = 1.0; + idx++; + } } + A->p[n] = idx; } - A->p[n] = n * n; - expr *AX = new_left_matmul(X, A); + expr *AX = new_left_matmul(NULL, X, A); + free_csr_matrix(A); double *x_vals = (double *) malloc(n * n * sizeof(double)); for (int i = 0; i < n * n; i++) @@ -55,8 +62,7 @@ const char *profile_left_matmul() printf("left_matmul jacobian eval time: %8.3f seconds\n", GET_ELAPSED_SECONDS(timer)); - free(x_vals); - free_csr_matrix(A); + FREE_AND_NULL(x_vals); free_expr(AX); return 0; } diff --git a/tests/utils/test_csr_csc_conversion.h b/tests/utils/test_csr_csc_conversion.h index c7daeb6..ba6f64a 100644 --- a/tests/utils/test_csr_csc_conversion.h +++ b/tests/utils/test_csr_csc_conversion.h @@ -3,6 +3,7 @@ #include #include +#include "memory_wrappers.h" #include "minunit.h" #include "test_helpers.h" #include "utils/CSC_Matrix.h" @@ -46,7 +47,7 @@ const char *test_csr_to_csc_split() mu_assert("C vals incorrect", cmp_double_array(C->x, Cx_correct, 5)); - free(iwork); + FREE_AND_NULL(iwork); free_csr_matrix(A); free_csc_matrix(C); @@ -90,7 +91,7 @@ const char *test_csc_to_csr_sparsity() mu_assert("C dimensions incorrect", C->m == 4 && C->n == 5); mu_assert("C nnz incorrect", C->nnz == 5); - free(iwork); + FREE_AND_NULL(iwork); free_csc_matrix(A); free_csr_matrix(C); @@ -123,7 +124,7 @@ const char *test_csc_to_csr_values() mu_assert("C vals incorrect", cmp_double_array(C->x, Cx_correct, 5)); - free(iwork); + FREE_AND_NULL(iwork); free_csc_matrix(A); free_csr_matrix(C); @@ -161,8 +162,8 @@ const char *test_csr_csc_csr_roundtrip() mu_assert("Round-trip: col indices incorrect", cmp_int_array(C->i, Ai, 8)); mu_assert("Round-trip: row pointers incorrect", cmp_int_array(C->p, Ap, 4)); - free(iwork_csc); - free(iwork_csr); + FREE_AND_NULL(iwork_csc); + FREE_AND_NULL(iwork_csr); free_csr_matrix(A); free_csc_matrix(B); free_csr_matrix(C); diff --git a/tests/utils/test_csr_matrix.h b/tests/utils/test_csr_matrix.h index 09b2e7d..c2d4f9c 100644 --- a/tests/utils/test_csr_matrix.h +++ b/tests/utils/test_csr_matrix.h @@ -2,6 +2,7 @@ #include #include +#include "memory_wrappers.h" #include "minunit.h" #include "test_helpers.h" #include "utils/CSR_Matrix.h" @@ -432,7 +433,7 @@ const char *test_AT_alloc_and_fill() free_csr_matrix(A); free_csr_matrix(AT); - free(iwork); + FREE_AND_NULL(iwork); return 0; } diff --git a/tests/wsum_hess/test_left_matmul.h b/tests/wsum_hess/test_left_matmul.h index 28d1cec..5722f25 100644 --- a/tests/wsum_hess/test_left_matmul.h +++ b/tests/wsum_hess/test_left_matmul.h @@ -8,6 +8,7 @@ #include "elementwise_univariate.h" #include "expr.h" #include "minunit.h" +#include "subexpr.h" #include "test_helpers.h" const char *test_wsum_hess_left_matmul() @@ -52,7 +53,7 @@ const char *test_wsum_hess_left_matmul() expr *x = new_variable(3, 1, 0, 3); - /* Create sparse matrix A in CSR format */ + /* A is 4x3: [1, 0, 2; 3, 0, 4; 5, 0, 6; 7, 0, 0] */ CSR_Matrix *A = new_csr_matrix(4, 3, 7); int A_p[5] = {0, 2, 4, 6, 7}; int A_i[7] = {0, 2, 0, 2, 0, 2, 0}; @@ -62,7 +63,8 @@ const char *test_wsum_hess_left_matmul() 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); + free_csr_matrix(A); A_log_x->forward(A_log_x, x_vals); A_log_x->jacobian_init(A_log_x); @@ -84,7 +86,6 @@ const char *test_wsum_hess_left_matmul() mu_assert("cols incorrect", cmp_int_array(A_log_x->wsum_hess->i, expected_i, 3)); mu_assert("rows incorrect", cmp_int_array(A_log_x->wsum_hess->p, expected_p, 4)); - free_csr_matrix(A); free_expr(A_log_x); return 0; } @@ -150,7 +151,7 @@ const char *test_wsum_hess_left_matmul_composite() memcpy(B->i, B_i, 9 * sizeof(int)); memcpy(B->x, B_x, 9 * sizeof(double)); - /* Create A matrix */ + /* A is 4x3: [1, 0, 2; 3, 0, 4; 5, 0, 6; 7, 0, 0] */ CSR_Matrix *A = new_csr_matrix(4, 3, 7); int A_p[5] = {0, 2, 4, 6, 7}; int A_i[7] = {0, 2, 0, 2, 0, 2, 0}; @@ -161,7 +162,8 @@ const char *test_wsum_hess_left_matmul_composite() expr *Bx = new_linear(x, B, NULL); expr *log_Bx = new_log(Bx); - expr *A_log_Bx = new_left_matmul(log_Bx, A); + expr *A_log_Bx = new_left_matmul(NULL, log_Bx, A); + free_csr_matrix(A); A_log_Bx->forward(A_log_Bx, x_vals); A_log_Bx->jacobian_init(A_log_Bx); @@ -186,7 +188,6 @@ const char *test_wsum_hess_left_matmul_composite() mu_assert("rows incorrect", cmp_int_array(A_log_Bx->wsum_hess->p, expected_p, 4)); - free_csr_matrix(A); free_csr_matrix(B); free_expr(A_log_Bx); return 0; @@ -223,7 +224,7 @@ const char *test_wsum_hess_left_matmul_matrix() expr *x = new_variable(3, 2, 0, 6); - /* Create sparse matrix A in CSR format */ + /* A is 4x3: [1, 0, 2; 3, 0, 4; 5, 0, 6; 7, 0, 0] */ CSR_Matrix *A = new_csr_matrix(4, 3, 7); int A_p[5] = {0, 2, 4, 6, 7}; int A_i[7] = {0, 2, 0, 2, 0, 2, 0}; @@ -233,7 +234,8 @@ const char *test_wsum_hess_left_matmul_matrix() 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); + free_csr_matrix(A); A_log_x->forward(A_log_x, x_vals); A_log_x->jacobian_init(A_log_x); @@ -257,7 +259,6 @@ const char *test_wsum_hess_left_matmul_matrix() mu_assert("cols incorrect", cmp_int_array(A_log_x->wsum_hess->i, expected_i, 6)); mu_assert("rows incorrect", cmp_int_array(A_log_x->wsum_hess->p, expected_p, 7)); - free_csr_matrix(A); free_expr(A_log_x); return 0; } diff --git a/tests/wsum_hess/test_right_matmul.h b/tests/wsum_hess/test_right_matmul.h index ca109f8..575f38c 100644 --- a/tests/wsum_hess/test_right_matmul.h +++ b/tests/wsum_hess/test_right_matmul.h @@ -33,7 +33,7 @@ const char *test_wsum_hess_right_matmul() 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); log_x_A->jacobian_init(log_x_A); @@ -83,7 +83,7 @@ const char *test_wsum_hess_right_matmul_vector() 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); log_x_A->jacobian_init(log_x_A); diff --git a/tests/wsum_hess/test_const_scalar_mult.h b/tests/wsum_hess/test_scalar_mult.h similarity index 92% rename from tests/wsum_hess/test_const_scalar_mult.h rename to tests/wsum_hess/test_scalar_mult.h index 7172be4..b4654d5 100644 --- a/tests/wsum_hess/test_const_scalar_mult.h +++ b/tests/wsum_hess/test_scalar_mult.h @@ -1,10 +1,12 @@ #include #include +#include "affine.h" #include "bivariate.h" #include "elementwise_univariate.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 */ @@ -20,7 +22,8 @@ const char *test_wsum_hess_const_scalar_mult_log_vector() /* Create scalar mult node: y = 2.5 * log(x) */ double a = 2.5; - expr *y = new_const_scalar_mult(a, log_node); + expr *a_node = new_parameter(1, 1, PARAM_FIXED, 3, &a); + expr *y = new_scalar_mult(a_node, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -65,7 +68,8 @@ const char *test_wsum_hess_const_scalar_mult_log_matrix() /* Create scalar mult node: y = 3.0 * log(x) */ double a = 3.0; - expr *y = new_const_scalar_mult(a, log_node); + expr *a_node = new_parameter(1, 1, PARAM_FIXED, 4, &a); + expr *y = new_scalar_mult(a_node, log_node); /* Forward pass */ y->forward(y, u_vals); diff --git a/tests/wsum_hess/test_const_vector_mult.h b/tests/wsum_hess/test_vector_mult.h similarity index 92% rename from tests/wsum_hess/test_const_vector_mult.h rename to tests/wsum_hess/test_vector_mult.h index 88bc127..6bee2b8 100644 --- a/tests/wsum_hess/test_const_vector_mult.h +++ b/tests/wsum_hess/test_vector_mult.h @@ -1,10 +1,12 @@ #include #include +#include "affine.h" #include "bivariate.h" #include "elementwise_univariate.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 */ @@ -20,7 +22,8 @@ const char *test_wsum_hess_const_vector_mult_log_vector() /* 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_node = new_parameter(3, 1, PARAM_FIXED, 3, a); + expr *y = new_vector_mult(a_node, log_node); /* Forward pass */ y->forward(y, u_vals); @@ -63,7 +66,8 @@ const char *test_wsum_hess_const_vector_mult_log_matrix() /* 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_node = new_parameter(2, 2, PARAM_FIXED, 4, a); + expr *y = new_vector_mult(a_node, log_node); /* Forward pass */ y->forward(y, u_vals);