|
| 1 | +#ifndef TEST_PARAM_SOURCE_REFRESH_H |
| 2 | +#define TEST_PARAM_SOURCE_REFRESH_H |
| 3 | + |
| 4 | +/* Composite param_source refresh across problem_update_params. |
| 5 | + * |
| 6 | + * A param_source subtree is evaluated by the owning atom's gated recursion, |
| 7 | + * not the main forward walk, so gated nodes *inside* a composite source |
| 8 | + * (promote, nested mults, cached matmuls) are only re-evaluated if the owner |
| 9 | + * marks the subtree before forwarding it. These tests pin the shapes that |
| 10 | + * served stale values before that mark existed: |
| 11 | + * |
| 12 | + * 1. left_matmul whose coefficient is p (.) A ("(p*A) @ x") |
| 13 | + * 2. quad_form whose matrix is g (.) Sigma ("quad_form(x, g*Sig)") |
| 14 | + * 3. two gated levels: coefficient (c*p) (.) A ("((2p)*A) @ x") |
| 15 | + * 4. kron whose left operand is p (.) A ("kron(p*A, X)") |
| 16 | + * |
| 17 | + * Bare-parameter sources are covered by test_param_prob.h; these are the |
| 18 | + * composite counterparts. |
| 19 | + */ |
| 20 | + |
| 21 | +#include <math.h> |
| 22 | +#include <stdio.h> |
| 23 | + |
| 24 | +#include "atoms/affine.h" |
| 25 | +#include "atoms/non_elementwise_full_dom.h" |
| 26 | +#include "expr.h" |
| 27 | +#include "minunit.h" |
| 28 | +#include "problem.h" |
| 29 | +#include "subexpr.h" |
| 30 | +#include "test_helpers.h" |
| 31 | + |
| 32 | +/* (p * A) @ x: the elementwise product of a promoted scalar parameter with a |
| 33 | + constant matrix feeds left_matmul as its param_source. */ |
| 34 | +const char *test_composite_source_left_matmul(void) |
| 35 | +{ |
| 36 | + int n = 2; |
| 37 | + |
| 38 | + expr *x = new_variable(2, 1, 0, n); |
| 39 | + expr *objective = new_sum(x, -1); |
| 40 | + |
| 41 | + double theta[1] = {2.0}; |
| 42 | + expr *p = new_parameter(1, 1, 0, n, theta); |
| 43 | + /* column-major A = [[1,2],[3,4]] */ |
| 44 | + double A_vals[4] = {1.0, 3.0, 2.0, 4.0}; |
| 45 | + expr *A_const = new_parameter(2, 2, PARAM_FIXED, n, A_vals); |
| 46 | + expr *coeff = new_vector_mult(new_promote(p, 2, 2), A_const); |
| 47 | + |
| 48 | + expr *constraint = new_left_matmul_dense(coeff, x, 2, 2, NULL); |
| 49 | + expr *constraints[1] = {constraint}; |
| 50 | + problem *prob = new_problem(objective, constraints, 1, false); |
| 51 | + |
| 52 | + expr *param_nodes[1] = {p}; |
| 53 | + problem_register_params(prob, param_nodes, 1); |
| 54 | + problem_init_derivatives(prob); |
| 55 | + |
| 56 | + double x_vals[2] = {1.0, 2.0}; |
| 57 | + int Ap[3] = {0, 2, 4}; |
| 58 | + int Ai[4] = {0, 1, 0, 1}; |
| 59 | + |
| 60 | + /* p = 2: effective A is 2*A, so A@x = 2*[5,11] */ |
| 61 | + problem_constraint_forward(prob, x_vals); |
| 62 | + problem_jacobian(prob); |
| 63 | + double constrs[2] = {10.0, 22.0}; |
| 64 | + double Ax[4] = {2.0, 4.0, 6.0, 8.0}; |
| 65 | + mu_assert("initial vals fail", |
| 66 | + cmp_double_array(prob->constraint_values, constrs, 2)); |
| 67 | + mu_assert("initial jac fail", cmp_double_array(prob->jacobian->x, Ax, 4)); |
| 68 | + mu_assert("rows fail", cmp_int_array(prob->jacobian->p, Ap, 3)); |
| 69 | + mu_assert("cols fail", cmp_int_array(prob->jacobian->i, Ai, 4)); |
| 70 | + |
| 71 | + /* p = 10: the promote inside the source must re-evaluate */ |
| 72 | + theta[0] = 10.0; |
| 73 | + problem_update_params(prob, theta); |
| 74 | + problem_constraint_forward(prob, x_vals); |
| 75 | + problem_jacobian(prob); |
| 76 | + constrs[0] = 50.0; |
| 77 | + constrs[1] = 110.0; |
| 78 | + Ax[0] = 10.0; |
| 79 | + Ax[1] = 20.0; |
| 80 | + Ax[2] = 30.0; |
| 81 | + Ax[3] = 40.0; |
| 82 | + mu_assert("stale constraint values after update", |
| 83 | + cmp_double_array(prob->constraint_values, constrs, 2)); |
| 84 | + mu_assert("stale jacobian values after update", |
| 85 | + cmp_double_array(prob->jacobian->x, Ax, 4)); |
| 86 | + |
| 87 | + free_problem(prob); |
| 88 | + |
| 89 | + return 0; |
| 90 | +} |
| 91 | + |
| 92 | +/* quad_form(x, g * Sigma): the scaled constant matrix feeds quad_form's |
| 93 | + cached Q through its param_source. */ |
| 94 | +const char *test_composite_source_quad_form(void) |
| 95 | +{ |
| 96 | + int n = 2; |
| 97 | + |
| 98 | + expr *x = new_variable(2, 1, 0, n); |
| 99 | + |
| 100 | + double theta[1] = {1.0}; |
| 101 | + expr *g = new_parameter(1, 1, 0, n, theta); |
| 102 | + /* Sigma = diag(2, 3); symmetric, so major order is irrelevant */ |
| 103 | + double S_vals[4] = {2.0, 0.0, 0.0, 3.0}; |
| 104 | + expr *S_const = new_parameter(2, 2, PARAM_FIXED, n, S_vals); |
| 105 | + expr *Q_src = new_vector_mult(new_promote(g, 2, 2), S_const); |
| 106 | + |
| 107 | + expr *objective = new_quad_form_dense(x, 2, NULL, Q_src); |
| 108 | + problem *prob = new_problem(objective, NULL, 0, false); |
| 109 | + |
| 110 | + expr *param_nodes[1] = {g}; |
| 111 | + problem_register_params(prob, param_nodes, 1); |
| 112 | + problem_init_derivatives(prob); |
| 113 | + |
| 114 | + double x_vals[2] = {1.0, 2.0}; |
| 115 | + |
| 116 | + /* g = 1: x'(g Sigma)x = 2*1 + 3*4 = 14; gradient 2(g Sigma)x = [4, 12] */ |
| 117 | + double obj_val = problem_objective_forward(prob, x_vals); |
| 118 | + problem_gradient(prob); |
| 119 | + double grad[2] = {4.0, 12.0}; |
| 120 | + mu_assert("initial obj fail", fabs(obj_val - 14.0) < 1e-10); |
| 121 | + mu_assert("initial grad fail", cmp_double_array(prob->gradient_values, grad, 2)); |
| 122 | + |
| 123 | + /* g = 5: Q's cached values must be re-copied from the re-evaluated source */ |
| 124 | + theta[0] = 5.0; |
| 125 | + problem_update_params(prob, theta); |
| 126 | + obj_val = problem_objective_forward(prob, x_vals); |
| 127 | + problem_gradient(prob); |
| 128 | + grad[0] = 20.0; |
| 129 | + grad[1] = 60.0; |
| 130 | + mu_assert("stale objective after update", fabs(obj_val - 70.0) < 1e-10); |
| 131 | + mu_assert("stale gradient after update", |
| 132 | + cmp_double_array(prob->gradient_values, grad, 2)); |
| 133 | + |
| 134 | + free_problem(prob); |
| 135 | + |
| 136 | + return 0; |
| 137 | +} |
| 138 | + |
| 139 | +/* ((c * p) * A) @ x: two gated levels inside the source (scalar_mult under |
| 140 | + promote under vector_mult) — the mark must recurse through all of them. */ |
| 141 | +const char *test_composite_source_nested_gates(void) |
| 142 | +{ |
| 143 | + int n = 2; |
| 144 | + |
| 145 | + expr *x = new_variable(2, 1, 0, n); |
| 146 | + expr *objective = new_sum(x, -1); |
| 147 | + |
| 148 | + double theta[1] = {1.0}; |
| 149 | + expr *p = new_parameter(1, 1, 0, n, theta); |
| 150 | + double c_val = 2.0; |
| 151 | + expr *c_const = new_parameter(1, 1, PARAM_FIXED, n, &c_val); |
| 152 | + expr *scaled = new_scalar_mult(p, c_const); /* = c * p, a gated node */ |
| 153 | + /* column-major A = [[1,2],[3,4]] */ |
| 154 | + double A_vals[4] = {1.0, 3.0, 2.0, 4.0}; |
| 155 | + expr *A_const = new_parameter(2, 2, PARAM_FIXED, n, A_vals); |
| 156 | + expr *coeff = new_vector_mult(new_promote(scaled, 2, 2), A_const); |
| 157 | + |
| 158 | + expr *constraint = new_left_matmul_dense(coeff, x, 2, 2, NULL); |
| 159 | + expr *constraints[1] = {constraint}; |
| 160 | + problem *prob = new_problem(objective, constraints, 1, false); |
| 161 | + |
| 162 | + expr *param_nodes[1] = {p}; |
| 163 | + problem_register_params(prob, param_nodes, 1); |
| 164 | + problem_init_derivatives(prob); |
| 165 | + |
| 166 | + double x_vals[2] = {1.0, 2.0}; |
| 167 | + |
| 168 | + /* p = 1: effective A is 2*A */ |
| 169 | + problem_constraint_forward(prob, x_vals); |
| 170 | + problem_jacobian(prob); |
| 171 | + double constrs[2] = {10.0, 22.0}; |
| 172 | + double Ax[4] = {2.0, 4.0, 6.0, 8.0}; |
| 173 | + mu_assert("initial vals fail", |
| 174 | + cmp_double_array(prob->constraint_values, constrs, 2)); |
| 175 | + mu_assert("initial jac fail", cmp_double_array(prob->jacobian->x, Ax, 4)); |
| 176 | + |
| 177 | + /* p = 5: effective A is 10*A; both gated levels must re-evaluate */ |
| 178 | + theta[0] = 5.0; |
| 179 | + problem_update_params(prob, theta); |
| 180 | + problem_constraint_forward(prob, x_vals); |
| 181 | + problem_jacobian(prob); |
| 182 | + constrs[0] = 50.0; |
| 183 | + constrs[1] = 110.0; |
| 184 | + Ax[0] = 10.0; |
| 185 | + Ax[1] = 20.0; |
| 186 | + Ax[2] = 30.0; |
| 187 | + Ax[3] = 40.0; |
| 188 | + mu_assert("stale constraint values after update", |
| 189 | + cmp_double_array(prob->constraint_values, constrs, 2)); |
| 190 | + mu_assert("stale jacobian values after update", |
| 191 | + cmp_double_array(prob->jacobian->x, Ax, 4)); |
| 192 | + |
| 193 | + free_problem(prob); |
| 194 | + |
| 195 | + return 0; |
| 196 | +} |
| 197 | + |
| 198 | +/* kron(p * A, X): the scaled constant operand feeds the kron atom's |
| 199 | + coefficient values through its param_source. */ |
| 200 | +const char *test_composite_source_kron(void) |
| 201 | +{ |
| 202 | + int n = 4; |
| 203 | + |
| 204 | + expr *X = new_variable(2, 2, 0, n); |
| 205 | + |
| 206 | + double theta[1] = {2.0}; |
| 207 | + expr *p = new_parameter(1, 1, 0, n, theta); |
| 208 | + /* column-major A = [[1,2],[3,4]]; all four blocks active */ |
| 209 | + double A_vals[4] = {1.0, 3.0, 2.0, 4.0}; |
| 210 | + expr *A_const = new_parameter(2, 2, PARAM_FIXED, n, A_vals); |
| 211 | + expr *coeff = new_vector_mult(new_promote(p, 2, 2), A_const); |
| 212 | + int active[4] = {0, 1, 2, 3}; |
| 213 | + |
| 214 | + expr *Z = new_left_kron(coeff, X, 2, 2, 2, 2, active, 4); |
| 215 | + expr *objective = new_sum(Z, -1); |
| 216 | + problem *prob = new_problem(objective, NULL, 0, false); |
| 217 | + |
| 218 | + expr *param_nodes[1] = {p}; |
| 219 | + problem_register_params(prob, param_nodes, 1); |
| 220 | + problem_init_derivatives(prob); |
| 221 | + |
| 222 | + double x_vals[4] = {1.0, 1.0, 1.0, 1.0}; |
| 223 | + |
| 224 | + /* p = 2: sum(kron(2A, ones)) = 2 * sum(A) * 4 = 80; |
| 225 | + gradient: each X entry sees sum(2A) = 20 */ |
| 226 | + double obj_val = problem_objective_forward(prob, x_vals); |
| 227 | + problem_gradient(prob); |
| 228 | + double grad[4] = {20.0, 20.0, 20.0, 20.0}; |
| 229 | + mu_assert("initial obj fail", fabs(obj_val - 80.0) < 1e-10); |
| 230 | + mu_assert("initial grad fail", cmp_double_array(prob->gradient_values, grad, 4)); |
| 231 | + |
| 232 | + /* p = 10: the promote and mult inside the operand must re-evaluate */ |
| 233 | + theta[0] = 10.0; |
| 234 | + problem_update_params(prob, theta); |
| 235 | + obj_val = problem_objective_forward(prob, x_vals); |
| 236 | + problem_gradient(prob); |
| 237 | + grad[0] = grad[1] = grad[2] = grad[3] = 100.0; |
| 238 | + mu_assert("stale objective after update", fabs(obj_val - 400.0) < 1e-10); |
| 239 | + mu_assert("stale gradient after update", |
| 240 | + cmp_double_array(prob->gradient_values, grad, 4)); |
| 241 | + |
| 242 | + free_problem(prob); |
| 243 | + |
| 244 | + return 0; |
| 245 | +} |
| 246 | + |
| 247 | +#endif /* TEST_PARAM_SOURCE_REFRESH_H */ |
0 commit comments