|
| 1 | +/* |
| 2 | + * Copyright 2026 Daniel Cederberg and William Zhang |
| 3 | + * |
| 4 | + * This file is part of the SparseDiffEngine project. |
| 5 | + * |
| 6 | + * Licensed under the Apache License, Version 2.0 (the "License"); |
| 7 | + * you may not use this file except in compliance with the License. |
| 8 | + * You may obtain a copy of the License at |
| 9 | + * |
| 10 | + * http://www.apache.org/licenses/LICENSE-2.0 |
| 11 | + * |
| 12 | + * Unless required by applicable law or agreed to in writing, software |
| 13 | + * distributed under the License is distributed on an "AS IS" BASIS, |
| 14 | + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 15 | + * See the License for the specific language governing permissions and |
| 16 | + * limitations under the License. |
| 17 | + */ |
| 18 | +#include "atoms/affine.h" |
| 19 | +#include "subexpr.h" |
| 20 | +#include "utils/CSR_matrix.h" |
| 21 | +#include "utils/sparse_matrix.h" |
| 22 | +#include "utils/tracked_alloc.h" |
| 23 | +#include <assert.h> |
| 24 | +#include <string.h> |
| 25 | + |
| 26 | +/* Kronecker product Z = kron(A, B), where one operand is variable-free (held by |
| 27 | + * param_source) and the other (child = node->left) carries the variables. |
| 28 | + * |
| 29 | + * With column-major flattening, output entry OUT = (i*r + k) + (j*s + l)*(p*r) |
| 30 | + * equals A[i,j] * B[k,l], so each output entry depends on a single child entry: |
| 31 | + * |
| 32 | + * Z[OUT] = coeff[OUT] * vec(child)[child_row[OUT]] |
| 33 | + * J[OUT, :] = coeff[OUT] * J_child[child_row[OUT], :] |
| 34 | + * |
| 35 | + * where coeff[OUT] = param_source->value[coeff_idx[OUT]]. The constructors fill |
| 36 | + * child_row/coeff_idx only for the constant operand's active (nonzero) blocks; |
| 37 | + * the remaining rows keep child_row == -1 and are structurally zero. */ |
| 38 | + |
| 39 | +/* Pull current parameter values through any broadcast/promote wrappers. */ |
| 40 | +static void refresh_param_values(kron_expr *knode) |
| 41 | +{ |
| 42 | + if (!knode->base.needs_parameter_refresh) |
| 43 | + { |
| 44 | + return; |
| 45 | + } |
| 46 | + |
| 47 | + knode->param_source->forward(knode->param_source, NULL); |
| 48 | + knode->base.needs_parameter_refresh = false; |
| 49 | +} |
| 50 | + |
| 51 | +static void forward(expr *node, const double *u) |
| 52 | +{ |
| 53 | + expr *child = node->left; |
| 54 | + kron_expr *knode = (kron_expr *) node; |
| 55 | + |
| 56 | + refresh_param_values(knode); |
| 57 | + child->forward(child, u); |
| 58 | + |
| 59 | + const double *a = knode->param_source->value; |
| 60 | + const double *x = child->value; |
| 61 | + double *y = node->value; |
| 62 | + for (int out = 0; out < node->size; out++) |
| 63 | + { |
| 64 | + int cr = knode->child_row[out]; |
| 65 | + y[out] = (cr < 0) ? 0.0 : a[knode->coeff_idx[out]] * x[cr]; |
| 66 | + } |
| 67 | +} |
| 68 | + |
| 69 | +static void jacobian_init_impl(expr *node) |
| 70 | +{ |
| 71 | + expr *child = node->left; |
| 72 | + kron_expr *knode = (kron_expr *) node; |
| 73 | + |
| 74 | + jacobian_init(child); |
| 75 | + |
| 76 | + /* Row OUT of the result copies the sparsity of child row child_row[OUT] |
| 77 | + (with repetition); inactive rows are empty. */ |
| 78 | + CSR_matrix *Jc = child->jacobian->to_csr(child->jacobian); |
| 79 | + |
| 80 | + int total = 0; |
| 81 | + for (int out = 0; out < node->size; out++) |
| 82 | + { |
| 83 | + int cr = knode->child_row[out]; |
| 84 | + if (cr >= 0) |
| 85 | + { |
| 86 | + total += Jc->p[cr + 1] - Jc->p[cr]; |
| 87 | + } |
| 88 | + } |
| 89 | + |
| 90 | + CSR_matrix *Jk = new_CSR_matrix(node->size, node->n_vars, total); |
| 91 | + int idx = 0; |
| 92 | + Jk->p[0] = 0; |
| 93 | + for (int out = 0; out < node->size; out++) |
| 94 | + { |
| 95 | + int cr = knode->child_row[out]; |
| 96 | + if (cr >= 0) |
| 97 | + { |
| 98 | + int row_nnz = Jc->p[cr + 1] - Jc->p[cr]; |
| 99 | + memcpy(Jk->i + idx, Jc->i + Jc->p[cr], row_nnz * sizeof(int)); |
| 100 | + idx += row_nnz; |
| 101 | + } |
| 102 | + Jk->p[out + 1] = idx; |
| 103 | + } |
| 104 | + node->jacobian = new_sparse_matrix(Jk); |
| 105 | +} |
| 106 | + |
| 107 | +static void eval_jacobian(expr *node) |
| 108 | +{ |
| 109 | + expr *child = node->left; |
| 110 | + kron_expr *knode = (kron_expr *) node; |
| 111 | + |
| 112 | + child->eval_jacobian(child); |
| 113 | + |
| 114 | + /* Sparsity is fixed after jacobian_init, so the row offsets still align; |
| 115 | + refill active rows as scale * child-row-values. */ |
| 116 | + CSR_matrix *Jc = child->jacobian->to_csr(child->jacobian); |
| 117 | + CSR_matrix *Jk = node->jacobian->to_csr(node->jacobian); |
| 118 | + const double *a = knode->param_source->value; |
| 119 | + |
| 120 | + int idx = 0; |
| 121 | + for (int out = 0; out < node->size; out++) |
| 122 | + { |
| 123 | + int cr = knode->child_row[out]; |
| 124 | + if (cr < 0) |
| 125 | + { |
| 126 | + continue; |
| 127 | + } |
| 128 | + double scale = a[knode->coeff_idx[out]]; |
| 129 | + for (int t = Jc->p[cr]; t < Jc->p[cr + 1]; t++) |
| 130 | + { |
| 131 | + Jk->x[idx++] = scale * Jc->x[t]; |
| 132 | + } |
| 133 | + } |
| 134 | +} |
| 135 | + |
| 136 | +static void wsum_hess_init_impl(expr *node) |
| 137 | +{ |
| 138 | + expr *child = node->left; |
| 139 | + |
| 140 | + wsum_hess_init(child); |
| 141 | + node->wsum_hess = child->wsum_hess->copy_sparsity(child->wsum_hess); |
| 142 | + /* backprop workspace: one weight per child entry */ |
| 143 | + node->work->dwork = (double *) sp_malloc(child->size * sizeof(double)); |
| 144 | +} |
| 145 | + |
| 146 | +static void eval_wsum_hess(expr *node, const double *w) |
| 147 | +{ |
| 148 | + expr *child = node->left; |
| 149 | + kron_expr *knode = (kron_expr *) node; |
| 150 | + const double *a = knode->param_source->value; |
| 151 | + double *w_prime = node->work->dwork; |
| 152 | + |
| 153 | + /* kron is affine in the child, so we only push the weights back through the |
| 154 | + gather: w'[child_row[OUT]] += coeff[OUT] * w[OUT]. Many output rows map to |
| 155 | + one child entry, hence the accumulation. */ |
| 156 | + memset(w_prime, 0, child->size * sizeof(double)); |
| 157 | + for (int out = 0; out < node->size; out++) |
| 158 | + { |
| 159 | + int cr = knode->child_row[out]; |
| 160 | + if (cr >= 0) |
| 161 | + { |
| 162 | + w_prime[cr] += a[knode->coeff_idx[out]] * w[out]; |
| 163 | + } |
| 164 | + } |
| 165 | + |
| 166 | + child->eval_wsum_hess(child, w_prime); |
| 167 | + memcpy(node->wsum_hess->x, child->wsum_hess->x, |
| 168 | + node->wsum_hess->nnz * sizeof(double)); |
| 169 | +} |
| 170 | + |
| 171 | +static bool is_affine(const expr *node) |
| 172 | +{ |
| 173 | + return node->left->is_affine(node->left); |
| 174 | +} |
| 175 | + |
| 176 | +static void free_type_data(expr *node) |
| 177 | +{ |
| 178 | + kron_expr *knode = (kron_expr *) node; |
| 179 | + sp_free(knode->child_row); |
| 180 | + sp_free(knode->coeff_idx); |
| 181 | + free_expr(knode->param_source); |
| 182 | + |
| 183 | + knode->child_row = NULL; |
| 184 | + knode->coeff_idx = NULL; |
| 185 | + knode->param_source = NULL; |
| 186 | +} |
| 187 | + |
| 188 | +/* Allocate a kron node and its (all-inactive) index arrays. The left/right |
| 189 | + constructors then fill the active rows. */ |
| 190 | +static kron_expr *new_kron_common(expr *param_node, expr *child, int p, int q, int r, |
| 191 | + int s) |
| 192 | +{ |
| 193 | + int size_out = (p * r) * (q * s); |
| 194 | + |
| 195 | + kron_expr *knode = (kron_expr *) sp_calloc(1, sizeof(kron_expr)); |
| 196 | + expr *node = &knode->base; |
| 197 | + init_expr(node, p * r, q * s, child->n_vars, forward, jacobian_init_impl, |
| 198 | + eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, |
| 199 | + free_type_data); |
| 200 | + node->left = child; |
| 201 | + expr_retain(child); |
| 202 | + |
| 203 | + knode->param_source = param_node; |
| 204 | + expr_retain(param_node); |
| 205 | + |
| 206 | + knode->child_row = (int *) sp_malloc(size_out * sizeof(int)); |
| 207 | + knode->coeff_idx = (int *) sp_malloc(size_out * sizeof(int)); |
| 208 | + for (int out = 0; out < size_out; out++) |
| 209 | + { |
| 210 | + knode->child_row[out] = -1; /* inactive until an active block fills it */ |
| 211 | + } |
| 212 | + |
| 213 | + knode->base.needs_parameter_refresh = true; |
| 214 | + return knode; |
| 215 | +} |
| 216 | + |
| 217 | +/* Z = kron(A, B) with A = param_node (p x q) the constant, B = child (r x s) the |
| 218 | + variable. active_blocks holds column-major indices i + j*p of A's nonzeros. */ |
| 219 | +expr *new_left_kron(expr *param_node, expr *child, int p, int q, int r, int s, |
| 220 | + const int *active_blocks, int n_active) |
| 221 | +{ |
| 222 | + kron_expr *knode = new_kron_common(param_node, child, p, q, r, s); |
| 223 | + int n_rows = p * r; |
| 224 | + for (int b = 0; b < n_active; b++) |
| 225 | + { |
| 226 | + int bidx = active_blocks[b]; /* = i + j*p into A */ |
| 227 | + assert(0 <= bidx && bidx < p * q); |
| 228 | + int i = bidx % p; |
| 229 | + int j = bidx / p; |
| 230 | + for (int l = 0; l < s; l++) |
| 231 | + { |
| 232 | + for (int k = 0; k < r; k++) |
| 233 | + { |
| 234 | + int out = (i * r + k) + (j * s + l) * n_rows; |
| 235 | + knode->child_row[out] = k + l * r; /* col-major into B */ |
| 236 | + knode->coeff_idx[out] = bidx; /* col-major into A */ |
| 237 | + } |
| 238 | + } |
| 239 | + } |
| 240 | + return &knode->base; |
| 241 | +} |
| 242 | + |
| 243 | +/* Z = kron(A, B) with A = child (p x q) the variable, B = param_node (r x s) the |
| 244 | + constant. active_blocks holds column-major indices k + l*r of B's nonzeros. */ |
| 245 | +expr *new_right_kron(expr *param_node, expr *child, int p, int q, int r, int s, |
| 246 | + const int *active_blocks, int n_active) |
| 247 | +{ |
| 248 | + kron_expr *knode = new_kron_common(param_node, child, p, q, r, s); |
| 249 | + int n_rows = p * r; |
| 250 | + for (int b = 0; b < n_active; b++) |
| 251 | + { |
| 252 | + int bidx = active_blocks[b]; /* = k + l*r into B */ |
| 253 | + assert(0 <= bidx && bidx < r * s); |
| 254 | + int k = bidx % r; |
| 255 | + int l = bidx / r; |
| 256 | + for (int j = 0; j < q; j++) |
| 257 | + { |
| 258 | + for (int i = 0; i < p; i++) |
| 259 | + { |
| 260 | + int out = (i * r + k) + (j * s + l) * n_rows; |
| 261 | + knode->child_row[out] = i + j * p; /* col-major into A */ |
| 262 | + knode->coeff_idx[out] = bidx; /* col-major into B */ |
| 263 | + } |
| 264 | + } |
| 265 | + } |
| 266 | + return &knode->base; |
| 267 | +} |
0 commit comments