From e103ac9cc71757115baa6d403dfc7935717e7a90 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Sat, 12 Jul 2025 10:18:16 -0400 Subject: [PATCH 01/13] sparse matrix multiplication and associated fixes --- DESCRIPTION | 4 +- Makefile | 2 +- R/engine_functions.R | 32 +++- R/enum.R | 1 + R/parse_expr.R | 21 ++- man/engine_functions.Rd | 32 +++- misc/dev/dev.cpp | 220 ++++++++++++++-------------- src/macpan2.cpp | 220 ++++++++++++++-------------- tests/testthat/test-expr-list.R | 4 +- tests/testthat/test-expr-parser.R | 9 +- tests/testthat/test-sparse-matrix.R | 59 ++++++++ 11 files changed, 364 insertions(+), 240 deletions(-) create mode 100644 tests/testthat/test-sparse-matrix.R diff --git a/DESCRIPTION b/DESCRIPTION index b741282dd..e85fc8189 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: macpan2 Title: Fast and Flexible Compartmental Modelling -Version: 3.0.0 +Version: 3.1.0 Authors@R: c( person("Steve Walker", email="swalk@mcmaster.ca", role=c("cre", "aut")), person("Weiguang Guan", role="aut"), @@ -10,7 +10,7 @@ Authors@R: c( person("Irena Papst", role="ctb"), person("Michael Li", role="ctb"), person("Kevin Zhao", role="ctb")) -Description: Fast and flexible compartmental modelling with Template Model Builder. +Description: Tools for building and calibrating compartmental models of infectious disease. License: GPL-3 Depends: R (>= 4.1.0) diff --git a/Makefile b/Makefile index 39728d46b..723b409d0 100644 --- a/Makefile +++ b/Makefile @@ -34,7 +34,7 @@ full-install: # haven't updated macpan.cpp (but have perhaps modified dev.cpp) # and (3) do not require a roxygen update. quick-install: enum-update enum-meth-update - R CMD INSTALL --no-multiarch --install-tests . + R CMD INSTALL --no-multiarch --install-tests --configure-args="CFLAGS=-g CXXFLAGS=-g" . quick-doc-install: R/*.R misc/dev/dev.cpp diff --git a/R/engine_functions.R b/R/engine_functions.R index ac5363138..eb5e9e185 100644 --- a/R/engine_functions.R +++ b/R/engine_functions.R @@ -220,13 +220,22 @@ #' #' * `x %*% y` : Standard matrix multiplication. #' * `x %x% y` : Kronecker product +#' * `sparse_mat_mult(x, i, j, y, z)` : Matrix multiplication +#' when the left matrix is represented as a column vector, `x`, +#' of non-zero elements and integer vectors of row, `i`, and +#' column, `j`, indices. The right matrix and the resulting +#' matrix are not represented as sparse matrices. #' #' ### Arguments #' -#' * `x` : A matrix. For the standard product, `x` -#' must have as many columns as `y` has rows. -#' * `y` : A matrix. For standard product, `y` -#' must have as many rows as `x` has columns. +#' * `x` : A matrix. +#' * `y` : A matrix. +#' * `i` : Integer vector the same length as `x` giving +#' zero-based row indices for sparse matrix representation. +#' * `j` : Integer vector the same length as `x` giving +#' zero-based column indices for sparse matrix representation. +#' * `z` : A matrix with dimensions equal to the result of +#' the sparse matrix multiplication (see details). #' #' ### Return #' @@ -239,6 +248,20 @@ #' engine_eval(~ (1:10) %x% t(1:10)) #' ``` #' +#' ### Details +#' +#' For standard matrix multiplication, `x %*% y`, the number of +#' columns of `x` equals the number of rows of `y`. +#' +#' +#' Think about `sparse_mat_mult(x, i, j, y, z)` as similar to +#' `z ~ x %*% y`, where `x` is represented differently. In +#' particular, the argument `x` is a column vector containing the +#' non-zero elements of the left matrix, `i` contains the +#' zero-based row indices associated with each element in `x`, +#' and `j` contains the zero-based column indices associated with +#' each element in `x`. +#' #' ## Parenthesis #' #' The order of operations can be enforced in the usual @@ -1057,6 +1080,7 @@ #' @aliases invlogit #' @aliases logit #' @aliases cumsum +#' @aliases sparse_mat_mult #' @aliases assign #' @aliases unpack NULL diff --git a/R/enum.R b/R/enum.R index f655bb28c..878ce82ea 100644 --- a/R/enum.R +++ b/R/enum.R @@ -60,6 +60,7 @@ valid_func_sigs = c( , "fwrap: invlogit(x)" , "fwrap: logit(x)" , "fwrap: cumsum(x)" + , "fwrap: sparse_mat_mult(x, i, j, y, z)" , "fwrap: assign(x, i, j, v)" , "fwrap: unpack(x, ...)" ) diff --git a/R/parse_expr.R b/R/parse_expr.R index 3b835c8aa..3ea959115 100644 --- a/R/parse_expr.R +++ b/R/parse_expr.R @@ -8,6 +8,21 @@ is_na = function(x) { return(i) } +## TODO: Allow for no-op functions? +arg_placeholders = function(x) { + rep_0 = rep(0L, length(x) - 1L) + has_no_args = length(rep_0) == 0 + is_not_tilde = as.character(x[[1]]) != "~" + is_call = is.call(x) + if (has_no_args & is_call & is_not_tilde) { + stop( + "The following function was called without arguments:\n " + , as.character(x[[1L]]) + ) + } + return(rep_0) +} + #' Generate an Arithmetic Expression Parser #' #' @param parser_name Name of the parsing function as a character @@ -101,14 +116,16 @@ make_expr_parser = function( # the list of expressions x$x = append(x$x, as.list(x$x[[expr_id]])[-1L]) + # add placeholders for the i and n that are associated # with these new arguments -- these placeholders # will get updated at the next recursion if it is # discovered that they are also functions - rep_0 = rep(0L, length(x$x[[expr_id]]) - 1L) + rep_0 = arg_placeholders(x$x[[expr_id]]) x$n = append(x$n, rep_0) x$i = append(x$i, rep_0) - + + # update the current expression to reflect how it has # been reduced x$n[[expr_id]] = length(x$x[[expr_id]]) - 1L diff --git a/man/engine_functions.Rd b/man/engine_functions.Rd index c54633a84..8b9c3e912 100644 --- a/man/engine_functions.Rd +++ b/man/engine_functions.Rd @@ -62,6 +62,7 @@ \alias{invlogit} \alias{logit} \alias{cumsum} +\alias{sparse_mat_mult} \alias{assign} \alias{unpack} \title{Functions Available in the Simulation Engine} @@ -315,15 +316,24 @@ engine_eval(~ recycle(t(1:4), 3, 4)) \itemize{ \item \code{x \%*\% y} : Standard matrix multiplication. \item \code{x \%x\% y} : Kronecker product +\item \code{sparse_mat_mult(x, i, j, y, z)} : Matrix multiplication +when the left matrix is represented as a column vector, \code{x}, +of non-zero elements and integer vectors of row, \code{i}, and +column, \code{j}, indices. The right matrix and the resulting +matrix are not represented as sparse matrices. } } \subsection{Arguments}{ \itemize{ -\item \code{x} : A matrix. For the standard product, \code{x} -must have as many columns as \code{y} has rows. -\item \code{y} : A matrix. For standard product, \code{y} -must have as many rows as \code{x} has columns. +\item \code{x} : A matrix. +\item \code{y} : A matrix. +\item \code{i} : Integer vector the same length as \code{x} giving +zero-based row indices for sparse matrix representation. +\item \code{j} : Integer vector the same length as \code{x} giving +zero-based column indices for sparse matrix representation. +\item \code{z} : A matrix with dimensions equal to the result of +the sparse matrix multiplication (see details). } } @@ -340,6 +350,20 @@ engine_eval(~ (1:10) \%x\% t(1:10)) }\if{html}{\out{}} } +\subsection{Details}{ + +For standard matrix multiplication, \code{x \%*\% y}, the number of +columns of \code{x} equals the number of rows of \code{y}. + +Think about \code{sparse_mat_mult(x, i, j, y, z)} as similar to +\code{z ~ x \%*\% y}, where \code{x} is represented differently. In +particular, the argument \code{x} is a column vector containing the +non-zero elements of the left matrix, \code{i} contains the +zero-based row indices associated with each element in \code{x}, +and \code{j} contains the zero-based column indices associated with +each element in \code{x}. +} + } \subsection{Parenthesis}{ diff --git a/misc/dev/dev.cpp b/misc/dev/dev.cpp index 1f584d26c..92bf6aa93 100644 --- a/misc/dev/dev.cpp +++ b/misc/dev/dev.cpp @@ -136,8 +136,9 @@ enum macpan2_func , MP2_INVLOGIT = 58 // fwrap: invlogit(x) , MP2_LOGIT = 59 // fwrap: logit(x) , MP2_CUMSUM = 60 // fwrap: cumsum(x) - , MP2_ASSIGN = 61 // fwrap: assign(x, i, j, v) - , MP2_UNPACK = 62 // fwrap: unpack(x, ...) + , MP2_SPARSE_MAT_MULT = 61 // fwrap: sparse_mat_mult(x, i, j, y, z) + , MP2_ASSIGN = 62 // fwrap: assign(x, i, j, v) + , MP2_UNPACK = 63 // fwrap: unpack(x, ...) }; enum macpan2_meth @@ -216,7 +217,17 @@ Type FUN(const Type x, const int &index) { \ } \ -#define MP2_ERR(CODE, MSG, FUN) SetError(CODE, MSG, row, FUN, args.all_rows(), args.all_cols(), args.all_type_ints(), t); \ +#define MP2_ERR(CODE, MSG, FUN) \ +SetError(CODE \ + , MSG \ + , row \ + , FUN \ + , args.all_rows() \ + , args.all_cols() \ + , args.all_type_ints() \ + , t \ +); \ +return empty_matrix_for_errors; \ // ENGINE FUNCTIONS USED BY OTHER ENGINE FUNCTIONS @@ -530,23 +541,18 @@ int getNthMatIndex( } template -int CheckIndices(matrix x, const std::vector &row_indices, const std::vector &col_indices) -{ +int CheckIndices(matrix x, const std::vector &row_indices, const std::vector &col_indices) { int rows = x.rows(); int cols = x.cols(); - for (int row_index : row_indices) - { - if (row_index < 0 || row_index >= rows) - { + for (int row_index : row_indices) { + if (row_index < 0 || row_index >= rows) { return 1; // Indices are not valid } } - for (int col_index : col_indices) - { - if (col_index < 0 || col_index >= cols) - { + for (int col_index : col_indices) { + if (col_index < 0 || col_index >= cols) { return 1; // Indices are not valid } } @@ -876,6 +882,12 @@ class ArgList return return_val; } + + int check_n_args(int lower, int upper) const { + if (items_.size() < lower) return 1; + if (items_.size() > upper) return 2; + return 0; + } // Method to set an error code void set_error_code(int error) @@ -996,6 +1008,7 @@ class ExprEvaluator // in either the c++ or r sense. matrix m, m1, m2, m3, m4, m5; // return values std::vector v, v1, v2, v3, v4, v5; // integer vectors + matrix empty_matrix_for_errors; bool is_finite_mat; vector u; // FIXME: why not std::vector here?? matrix Y, X, A; @@ -1158,13 +1171,11 @@ class ExprEvaluator if (is_int_in(table_x[row] + 1, mp_math)) { if (args.all_matrices() == 0) { MP2_ERR(205, "All arguments to math functions must be numeric matrices, but at least one is an integer", table_x[row] + 1); - return m; } } if (is_int_in(table_x[row] + 1, mp_history)) { if (!mats_save_hist[index2mats[0]]) { MP2_ERR(205, "All arguments to functions that act on the simulation history must have a first argument that is a non-empty matrix with saved history", table_x[row] + 1); - return m; } } @@ -1178,19 +1189,15 @@ class ExprEvaluator switch (err_code) { case 201: MP2_ERR(201, "The two operands do not have the same number of columns", table_x[row] + 1); - return m; case 202: MP2_ERR(202, "The two operands do not have the same number of rows", table_x[row] + 1); - return m; case 203: MP2_ERR(203, "The two operands do not have the same number of columns or rows", table_x[row] + 1); - return m; } } else if (table_x[row] + 1 == MP2_MATRIX_MULTIPLY) { // %*% matrix multiplication if (args[0].cols() != args[1].rows()) { MP2_ERR(204, "The two operands are not compatible to do matrix multiplication", table_x[row] + 1); - return m; } } @@ -1457,7 +1464,6 @@ class ExprEvaluator to = args.get_as_int(1); if (from > to) { MP2_ERR(MP2_COLON, "Lower bound greater than upper bound in : operation", MP2_COLON); - return m; } m = matrix::Zero(to - from + 1, 1); for (int i = from; i <= to; i++) @@ -1483,7 +1489,6 @@ class ExprEvaluator by = args[2].coeff(0, 0); if (length <= 0) { MP2_ERR(MP2_SEQUENCE, "Sequence length is less than or equal to zero in seq operation", MP2_SEQUENCE); - return m; } m = matrix::Zero(length, 1); for (int i = 0; i < length; i++) @@ -1550,7 +1555,6 @@ class ExprEvaluator m = args[0]; if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_RECYCLE); - return m; } return m; @@ -1560,13 +1564,22 @@ class ExprEvaluator // #' // #' * `x %*% y` : Standard matrix multiplication. // #' * `x %x% y` : Kronecker product + // #' * `sparse_mat_mult(x, i, j, y, z)` : Matrix multiplication + // #' when the left matrix is represented as a column vector, `x`, + // #' of non-zero elements and integer vectors of row, `i`, and + // #' column, `j`, indices. The right matrix and the resulting + // #' matrix are not represented as sparse matrices. // #' // #' ### Arguments // #' - // #' * `x` : A matrix. For the standard product, `x` - // #' must have as many columns as `y` has rows. - // #' * `y` : A matrix. For standard product, `y` - // #' must have as many rows as `x` has columns. + // #' * `x` : A matrix. + // #' * `y` : A matrix. + // #' * `i` : Integer vector the same length as `x` giving + // #' zero-based row indices for sparse matrix representation. + // #' * `j` : Integer vector the same length as `x` giving + // #' zero-based column indices for sparse matrix representation. + // #' * `z` : A matrix with dimensions equal to the result of + // #' the sparse matrix multiplication (see details). // #' // #' ### Return // #' @@ -1579,6 +1592,11 @@ class ExprEvaluator // #' engine_eval(~ (1:10) %x% t(1:10)) // #' ``` // #' + // #' ### Details + // #' + // #' For standard matrix multiplication, `x %*% y`, the number of + // #' columns of `x` equals the number of rows of `y`. + // #' case MP2_MATRIX_MULTIPLY: // %*% return args[0] * args[1]; @@ -1593,6 +1611,52 @@ class ExprEvaluator } return m; + // #' + // #' Think about `sparse_mat_mult(x, i, j, y, z)` as similar to + // #' `z ~ x %*% y`, where `x` is represented differently. In + // #' particular, the argument `x` is a column vector containing the + // #' non-zero elements of the left matrix, `i` contains the + // #' zero-based row indices associated with each element in `x`, + // #' and `j` contains the zero-based column indices associated with + // #' each element in `x`. + // #' + case MP2_SPARSE_MAT_MULT: // sparse_mat_mult(x, i, j, y, z) + err_code = args.check_n_args(5, 5); + if (err_code == 1) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "Require exactly 5 arguments, but got fewer.", MP2_SPARSE_MAT_MULT); + } + if (err_code == 2) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "Require exactly 5 arguments, but got more.", MP2_SPARSE_MAT_MULT); + } + m = args[0]; // x + if (m.cols() != 1) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "First matrix must be represented as a vector with associated row and column indices.", MP2_SPARSE_MAT_MULT); + } + v1 = args.get_as_int_vec(1); // i + v2 = args.get_as_int_vec(2); // j + m1 = args[3]; // B + rows = args.rows(4); // nrow(Y) + cols = args.cols(4); // ncol(Y) + m2 = matrix::Zero(rows, cols); // Y + + v.push_back(m1.cols() - 1); + err_code = CheckIndices(m2, v1, v); + if (err_code) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "Number of rows in the output matrix must be greater than all of the row indices in the sparse matrix representation.", MP2_SPARSE_MAT_MULT); + } + err_code = CheckIndices(m1, v2, v); + if (err_code) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "Number of rows in the second (dense) matrix must be greater than all of the column indices in the sparse matrix representation.", MP2_SPARSE_MAT_MULT); + } + + for (int l = 0; l < m.rows(); l++) { + for (int j = 0; j < m1.cols(); j++) { + m2.coeffRef(v1[l], j) += m.coeff(l, 0) * m1.coeff(v2[l], j); + } + } + return m2; + + // #' ## Parenthesis // #' // #' The order of operations can be enforced in the usual @@ -1697,7 +1761,6 @@ class ExprEvaluator } else { MP2_ERR(MP2_CBIND, "Inconsistent size in cbind function", MP2_CBIND); - return m; } } return m; @@ -1711,26 +1774,20 @@ class ExprEvaluator int totrows, rowmarker; totrows = 0; rowmarker = 0; - for (int j = 0; j < n; j++) - { + for (int j = 0; j < n; j++) { totrows += args[j].rows(); } m = matrix::Zero(totrows, cols); - for (int i = 0; i < n; i++) - { - if (args[i].cols() == cols) - { + for (int i = 0; i < n; i++) { + if (args[i].cols() == cols) { rows_per_arg = args[i].rows(); - for (int k = 0; k < rows_per_arg; k++) - { + for (int k = 0; k < rows_per_arg; k++) { m.row(rowmarker + k) = args[i].row(k); } rowmarker += rows_per_arg; } - else - { + else { MP2_ERR(MP2_RBIND, "Inconsistent size in rbind function", MP2_RBIND); - return m; } } return m; @@ -1759,7 +1816,6 @@ class ExprEvaluator } else if (size_in < size_out) { if (size_out % size_in) { MP2_ERR(MP2_MATRIX, "The size of the input matrix is not compatible with the requested shape of the output matrix.", MP2_MATRIX); - return m1; } m = mp2_rep(m, size_out / size_in); m.resize(rows, cols); @@ -1903,12 +1959,11 @@ class ExprEvaluator m = args[0].colwise().sum().matrix(); return m; - case MP2_GROUPSUMS: // group_sums(x) + case MP2_GROUPSUMS: // group_sums(x, f, n) v1 = args.get_as_int_vec(1); err_code = args.check_indices(2, v1, {0}); if (err_code) { MP2_ERR(MP2_GROUPSUMS, "Group indexes are out of range.", MP2_GROUPSUMS); - return m; } m = args[0]; if (m.cols() != 1) { @@ -1918,7 +1973,6 @@ class ExprEvaluator m1 = matrix::Zero(rows, 1); if (v1.size() != m.rows()) { MP2_ERR(MP2_GROUPSUMS, "Number of rows in x must equal the number of indices in f in group_sums(x, f, n).", MP2_GROUPSUMS); - return m; } for (int i = 0; i < m.rows(); i++) { m1.coeffRef(v1[i], 0) += m.coeff(i, 0); @@ -1928,7 +1982,6 @@ class ExprEvaluator case MP2_MEAN: // mean(x) if (n != 1) { MP2_ERR(MP2_MEAN, "The mean function can only take a single matrix.", MP2_MEAN); - return m; } m = matrix::Zero(1, 1); sum = args.get_as_mat(0).mean(); @@ -1938,7 +1991,6 @@ class ExprEvaluator case MP2_SD: // sd(x) if (n != 1) { MP2_ERR(MP2_SD, "The sd function can only take a single matrix.", MP2_SD); - return m; } m = args.get_as_mat(0); m1 = matrix::Zero(1, 1); @@ -2058,7 +2110,6 @@ class ExprEvaluator err_code = args.check_indices(0, v1, v2); if (err_code) { MP2_ERR(MP2_SQUARE_BRACKET, "Illegal index to square bracket", MP2_SQUARE_BRACKET); - return m; } m = args[0]; m1 = matrix::Zero(nrow, ncol); @@ -2082,12 +2133,10 @@ class ExprEvaluator err_code = args.check_indices(0, v1, v2); if (err_code){ MP2_ERR(MP2_BLOCK, "Illegal starting index to block", MP2_BLOCK); - return m; } err_code = args.check_indices(0, v3, v4); if (err_code){ MP2_ERR(MP2_BLOCK, "Illegal index to block, requesting more elements than available in input", MP2_BLOCK); - return m; } return args[0].block(rowIndex, colIndex, rows, cols); @@ -2151,14 +2200,12 @@ class ExprEvaluator return m; // have not built up any previous iterations yet, so returning empty matrix if (t == 0) { MP2_ERR(154, "The simulation loop has not yet begun and so rbind_time (or rbind_lag) cannot be used", int_func); - return m; } matIndex = index2mats[0]; // m if ((matIndex < 0) | (index2what[0] != 0)) { MP2_ERR(MP2_RBIND_TIME, "Can only rbind_time (or rbind_lag) named matrices not expressions of matrices and not integer vectors", int_func); // Rcpp::Rcout << "return 2 " << std::endl; - return m; } if ((n == 1) & (!doing_lag)) { @@ -2184,7 +2231,6 @@ class ExprEvaluator // (usually a column vector) of initial values that take // us back into negative time steps. need to do the same // for convolution and anything else that looks backwards. - return m; } } // timeIndex = -timeIndex; @@ -2198,11 +2244,9 @@ class ExprEvaluator lowerTimeBound = args.get_as_int(2); if (lowerTimeBound < 0) { MP2_ERR(MP2_RBIND_TIME, "Lower time bound (third argument) is less than zero", int_func); - return m; } if (lowerTimeBound > t_max) { MP2_ERR(MP2_RBIND_TIME, "Lower time bound (third argument) is greater than the number of time steps", int_func); - return m; } } else if (doing_lag) { @@ -2249,7 +2293,6 @@ class ExprEvaluator { // Shall we allow inconsistent rows? MP2_ERR(MP2_RBIND_TIME, "Inconsistent rows or columns in rbind_time (or rbind_lag)", int_func); // Rcpp::Rcout << "return 5 " << std::endl; - return args[0]; } } @@ -2365,7 +2408,6 @@ class ExprEvaluator if (lag < 0) { MP2_ERR(MP2_TIME_STEP, "Time lag needs to be non-negative", MP2_TIME_STEP); - return m; } if (t > lag) { @@ -2376,7 +2418,6 @@ class ExprEvaluator case MP2_TIME_GROUP: // time_group(i, change_points) if (index2what[0] != 0) { MP2_ERR(MP2_TIME_GROUP, "First argument needs to be a matrix.", MP2_TIME_GROUP); - return m; } m = args[0]; off = args.get_as_int(0); @@ -2390,22 +2431,18 @@ class ExprEvaluator case MP2_TIME_VAR: // time_var(x, change_points) if (t == 0) { MP2_ERR(MP2_TIME_VAR, "Time variation is not allowed before the simulation loop begins.", MP2_TIME_VAR); - return m; } if (t > t_max) { MP2_ERR(MP2_TIME_VAR, "Time variation is not allowed after the simulation loop ends.", MP2_TIME_VAR); - return m; } v = args.get_as_int_vec(1); off = args.get_as_int(1); if (off < 0) { MP2_ERR(MP2_TIME_VAR, "The first element of the second argument must not be less than zero.", MP2_TIME_VAR); - return m; } if (off >= v.size()) { MP2_ERR(MP2_TIME_VAR, "The first element of the second argument must be less than the number of elements in the first.", MP2_TIME_VAR); - return m; } // first argument can have its rows indexed @@ -2427,12 +2464,10 @@ class ExprEvaluator } else { MP2_ERR(MP2_TIME_VAR, "Time variation pointers need to be length-1 integer vectors.", MP2_TIME_VAR); - return m; } } else if (cp < 0) { MP2_ERR(MP2_TIME_VAR, "Negative times are not allowed.", MP2_TIME_VAR); - return m; } } m = matrix::Zero(args[0].cols(), 1); @@ -2506,7 +2541,6 @@ class ExprEvaluator matIndex = index2mats[0]; // m if (matIndex == -1) { MP2_ERR(MP2_CONVOLUTION, "Can only convolve named matrices not expressions of matrices", MP2_CONVOLUTION); - return args[0]; } #ifdef MP_VERBOSE @@ -2549,10 +2583,8 @@ class ExprEvaluator return m; } - else - { + else { MP2_ERR(MP2_CONVOLUTION, "Either empty or non-column vector used as kernel in convolution", MP2_CONVOLUTION); - return m; } // #' ## Clamp @@ -2686,7 +2718,6 @@ class ExprEvaluator case MP2_POISSON_DENSITY: if (n < 2) { MP2_ERR(MP2_POISSON_DENSITY, "dpois needs two arguments: matrices with observed and expected values", MP2_POISSON_DENSITY); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2695,7 +2726,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_POISSON_DENSITY); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) { @@ -2708,7 +2738,6 @@ class ExprEvaluator case MP2_NEGBIN_DENSITY: if (n < 3) { MP2_ERR(MP2_NEGBIN_DENSITY, "dnbinom needs three arguments: matrices with observed values, expected values, and dispersion parameters", MP2_NEGBIN_DENSITY); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2718,7 +2747,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_NEGBIN_DENSITY); - return m; } // var ~ variance // mu ~ mean @@ -2739,7 +2767,6 @@ class ExprEvaluator case MP2_NORMAL_DENSITY: if (n < 3) { MP2_ERR(MP2_NORMAL_DENSITY, "dnorm needs three arguments: matrices with observed values, expected values, and standard deviation parameters", MP2_NORMAL_DENSITY); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2749,7 +2776,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_NORMAL_DENSITY); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) { @@ -2762,7 +2788,6 @@ class ExprEvaluator case MP2_BINOM_DENSITY: if (n < 3) { MP2_ERR(MP2_BINOM_DENSITY, "dbinom needs three arguments: matrices with observed values, numbers of trials, and probabilities", MP2_BINOM_DENSITY); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2772,7 +2797,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_BINOM_DENSITY); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) { @@ -2855,10 +2879,8 @@ class ExprEvaluator return m; case MP2_NEGBIN_SIM: - if (n < 2) - { + if (n < 2) { MP2_ERR(MP2_NEGBIN_SIM, "rnbinom needs two arguments: matrices with means and dispersion parameters", MP2_NEGBIN_SIM); - return m; } eps = 1e-8; rows = args[0].rows(); @@ -2866,10 +2888,8 @@ class ExprEvaluator v1.push_back(1); args = args.recycle_to_shape(v1, rows, cols); err_code = args.get_error_code(); - if (err_code != 0) - { + if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_NEGBIN_SIM); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) @@ -2888,10 +2908,8 @@ class ExprEvaluator return m; case MP2_NORMAL_SIM: - if (n < 2) - { + if (n < 2) { MP2_ERR(MP2_NORMAL_SIM, "rnorm needs two arguments: matrices with means and standard deviations", MP2_NORMAL_SIM); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2899,10 +2917,8 @@ class ExprEvaluator args = args.recycle_to_shape(v1, rows, cols); err_code = args.get_error_code(); // err_code = RecycleInPlace(args[1], rows, cols); - if (err_code != 0) - { + if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_NORMAL_SIM); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) @@ -2918,7 +2934,6 @@ class ExprEvaluator // rbinom(size, prob) if (n != 2) { MP2_ERR(MP2_BINOM_SIM, "rbinom needs two arguments: matrices with size and probability", MP2_BINOM_SIM); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2927,7 +2942,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(MP2_BINOM_SIM, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_BINOM_SIM); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) { @@ -2948,7 +2962,6 @@ class ExprEvaluator //Rcpp::Rcout << "++++++" << std::endl; //Rcpp::Rcout << args[0] << std::endl; MP2_ERR(MP2_EULER_MULTINOM_SIM, "The first 'size' argument must be scalar.", MP2_EULER_MULTINOM_SIM); - return m; } if (args[1].cols() != 1) { //Rcpp::Rcout << "------" << std::endl; @@ -3018,7 +3031,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_PGAMMA); - return m; } m1 = args.get_as_mat(0); // q m2 = args.get_as_mat(1); // shape @@ -3040,7 +3052,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_PGAMMA); - return m; } m1 = args.get_as_mat(0); // q m2 = args.get_as_mat(1); // mean @@ -3173,22 +3184,16 @@ class ExprEvaluator case MP2_ASSIGN: cols = args[1].cols(); - if (cols != 1) - { + if (cols != 1) { MP2_ERR(255, "Assignment index matrices must have a single column", MP2_ASSIGN); - return m; } cols = args[2].cols(); - if (cols != 1) - { + if (cols != 1) { MP2_ERR(255, "Assignment index matrices must have a single column", MP2_ASSIGN); - return m; } cols = args[3].cols(); - if (cols != 1) - { + if (cols != 1) { MP2_ERR(255, "Assignment value matrices must have a single column", MP2_ASSIGN); - return m; } // apparently we still need this check @@ -3198,10 +3203,8 @@ class ExprEvaluator // printIntVector(v2); err_code = args.check_indices(0, v1, v2); // CheckIndices(args[0], args[1], m1); // err_code = CheckIndices(args[0], args[1], args[2]); - if (err_code) - { + if (err_code) { MP2_ERR(MP2_ASSIGN, "Illegal index used in assign", MP2_ASSIGN); - return m; } rows = args[3].rows(); @@ -3211,7 +3214,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_ASSIGN); - return m; } for (int k = 0; k < rows; k++) { @@ -3220,7 +3222,6 @@ class ExprEvaluator matIndex = index2mats[0]; if (matIndex == -1) { MP2_ERR(MP2_ASSIGN, "Can only assign to named matrices not expressions of matrices", MP2_ASSIGN); - return args[0]; } valid_vars.m_matrices[matIndex].coeffRef(rowIndex, colIndex) = args[3].coeff(k, 0); } @@ -3279,19 +3280,15 @@ class ExprEvaluator m.resize(size, 1); start = 0; - for (int i = 1; i < n; i++) - { + for (int i = 1; i < n; i++) { sz = args[i].rows() * args[i].cols(); - if (size >= sz) - { + if (size >= sz) { m1 = m.block(start, 0, sz, 1); m1.resize(args[i].rows(), args[i].cols()); // Rcpp::Rcout << "MATRIX " << valid_vars.m_matrices[index2mats[i]] << std::endl << std::endl; matIndex = index2mats[i]; - if (matIndex == -1) - { + if (matIndex == -1) { MP2_ERR(MP2_ASSIGN, "Can only unpack into named matrices not expressions of matrices", MP2_UNPACK); - return args[0]; } valid_vars.m_matrices[matIndex] = m1; // args[i] = m1; @@ -3305,7 +3302,6 @@ class ExprEvaluator default: MP2_ERR(255, "invalid operator in arithmetic expression", -99); - return m; } } // switch (table_n[row]) }; diff --git a/src/macpan2.cpp b/src/macpan2.cpp index 2171a7fcd..fff47462c 100644 --- a/src/macpan2.cpp +++ b/src/macpan2.cpp @@ -137,8 +137,9 @@ enum macpan2_func , MP2_INVLOGIT = 58 // fwrap: invlogit(x) , MP2_LOGIT = 59 // fwrap: logit(x) , MP2_CUMSUM = 60 // fwrap: cumsum(x) - , MP2_ASSIGN = 61 // fwrap: assign(x, i, j, v) - , MP2_UNPACK = 62 // fwrap: unpack(x, ...) + , MP2_SPARSE_MAT_MULT = 61 // fwrap: sparse_mat_mult(x, i, j, y, z) + , MP2_ASSIGN = 62 // fwrap: assign(x, i, j, v) + , MP2_UNPACK = 63 // fwrap: unpack(x, ...) }; enum macpan2_meth @@ -217,7 +218,17 @@ Type FUN(const Type x, const int &index) { \ } \ -#define MP2_ERR(CODE, MSG, FUN) SetError(CODE, MSG, row, FUN, args.all_rows(), args.all_cols(), args.all_type_ints(), t); \ +#define MP2_ERR(CODE, MSG, FUN) \ +SetError(CODE \ + , MSG \ + , row \ + , FUN \ + , args.all_rows() \ + , args.all_cols() \ + , args.all_type_ints() \ + , t \ +); \ +return empty_matrix_for_errors; \ // ENGINE FUNCTIONS USED BY OTHER ENGINE FUNCTIONS @@ -531,23 +542,18 @@ int getNthMatIndex( } template -int CheckIndices(matrix x, const std::vector &row_indices, const std::vector &col_indices) -{ +int CheckIndices(matrix x, const std::vector &row_indices, const std::vector &col_indices) { int rows = x.rows(); int cols = x.cols(); - for (int row_index : row_indices) - { - if (row_index < 0 || row_index >= rows) - { + for (int row_index : row_indices) { + if (row_index < 0 || row_index >= rows) { return 1; // Indices are not valid } } - for (int col_index : col_indices) - { - if (col_index < 0 || col_index >= cols) - { + for (int col_index : col_indices) { + if (col_index < 0 || col_index >= cols) { return 1; // Indices are not valid } } @@ -877,6 +883,12 @@ class ArgList return return_val; } + + int check_n_args(int lower, int upper) const { + if (items_.size() < lower) return 1; + if (items_.size() > upper) return 2; + return 0; + } // Method to set an error code void set_error_code(int error) @@ -997,6 +1009,7 @@ class ExprEvaluator // in either the c++ or r sense. matrix m, m1, m2, m3, m4, m5; // return values std::vector v, v1, v2, v3, v4, v5; // integer vectors + matrix empty_matrix_for_errors; bool is_finite_mat; vector u; // FIXME: why not std::vector here?? matrix Y, X, A; @@ -1159,13 +1172,11 @@ class ExprEvaluator if (is_int_in(table_x[row] + 1, mp_math)) { if (args.all_matrices() == 0) { MP2_ERR(205, "All arguments to math functions must be numeric matrices, but at least one is an integer", table_x[row] + 1); - return m; } } if (is_int_in(table_x[row] + 1, mp_history)) { if (!mats_save_hist[index2mats[0]]) { MP2_ERR(205, "All arguments to functions that act on the simulation history must have a first argument that is a non-empty matrix with saved history", table_x[row] + 1); - return m; } } @@ -1179,19 +1190,15 @@ class ExprEvaluator switch (err_code) { case 201: MP2_ERR(201, "The two operands do not have the same number of columns", table_x[row] + 1); - return m; case 202: MP2_ERR(202, "The two operands do not have the same number of rows", table_x[row] + 1); - return m; case 203: MP2_ERR(203, "The two operands do not have the same number of columns or rows", table_x[row] + 1); - return m; } } else if (table_x[row] + 1 == MP2_MATRIX_MULTIPLY) { // %*% matrix multiplication if (args[0].cols() != args[1].rows()) { MP2_ERR(204, "The two operands are not compatible to do matrix multiplication", table_x[row] + 1); - return m; } } @@ -1458,7 +1465,6 @@ class ExprEvaluator to = args.get_as_int(1); if (from > to) { MP2_ERR(MP2_COLON, "Lower bound greater than upper bound in : operation", MP2_COLON); - return m; } m = matrix::Zero(to - from + 1, 1); for (int i = from; i <= to; i++) @@ -1484,7 +1490,6 @@ class ExprEvaluator by = args[2].coeff(0, 0); if (length <= 0) { MP2_ERR(MP2_SEQUENCE, "Sequence length is less than or equal to zero in seq operation", MP2_SEQUENCE); - return m; } m = matrix::Zero(length, 1); for (int i = 0; i < length; i++) @@ -1551,7 +1556,6 @@ class ExprEvaluator m = args[0]; if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_RECYCLE); - return m; } return m; @@ -1561,13 +1565,22 @@ class ExprEvaluator // #' // #' * `x %*% y` : Standard matrix multiplication. // #' * `x %x% y` : Kronecker product + // #' * `sparse_mat_mult(x, i, j, y, z)` : Matrix multiplication + // #' when the left matrix is represented as a column vector, `x`, + // #' of non-zero elements and integer vectors of row, `i`, and + // #' column, `j`, indices. The right matrix and the resulting + // #' matrix are not represented as sparse matrices. // #' // #' ### Arguments // #' - // #' * `x` : A matrix. For the standard product, `x` - // #' must have as many columns as `y` has rows. - // #' * `y` : A matrix. For standard product, `y` - // #' must have as many rows as `x` has columns. + // #' * `x` : A matrix. + // #' * `y` : A matrix. + // #' * `i` : Integer vector the same length as `x` giving + // #' zero-based row indices for sparse matrix representation. + // #' * `j` : Integer vector the same length as `x` giving + // #' zero-based column indices for sparse matrix representation. + // #' * `z` : A matrix with dimensions equal to the result of + // #' the sparse matrix multiplication (see details). // #' // #' ### Return // #' @@ -1580,6 +1593,11 @@ class ExprEvaluator // #' engine_eval(~ (1:10) %x% t(1:10)) // #' ``` // #' + // #' ### Details + // #' + // #' For standard matrix multiplication, `x %*% y`, the number of + // #' columns of `x` equals the number of rows of `y`. + // #' case MP2_MATRIX_MULTIPLY: // %*% return args[0] * args[1]; @@ -1594,6 +1612,52 @@ class ExprEvaluator } return m; + // #' + // #' Think about `sparse_mat_mult(x, i, j, y, z)` as similar to + // #' `z ~ x %*% y`, where `x` is represented differently. In + // #' particular, the argument `x` is a column vector containing the + // #' non-zero elements of the left matrix, `i` contains the + // #' zero-based row indices associated with each element in `x`, + // #' and `j` contains the zero-based column indices associated with + // #' each element in `x`. + // #' + case MP2_SPARSE_MAT_MULT: // sparse_mat_mult(x, i, j, y, z) + err_code = args.check_n_args(5, 5); + if (err_code == 1) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "Require exactly 5 arguments, but got fewer.", MP2_SPARSE_MAT_MULT); + } + if (err_code == 2) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "Require exactly 5 arguments, but got more.", MP2_SPARSE_MAT_MULT); + } + m = args[0]; // x + if (m.cols() != 1) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "First matrix must be represented as a vector with associated row and column indices.", MP2_SPARSE_MAT_MULT); + } + v1 = args.get_as_int_vec(1); // i + v2 = args.get_as_int_vec(2); // j + m1 = args[3]; // B + rows = args.rows(4); // nrow(Y) + cols = args.cols(4); // ncol(Y) + m2 = matrix::Zero(rows, cols); // Y + + v.push_back(m1.cols() - 1); + err_code = CheckIndices(m2, v1, v); + if (err_code) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "Number of rows in the output matrix must be greater than all of the row indices in the sparse matrix representation.", MP2_SPARSE_MAT_MULT); + } + err_code = CheckIndices(m1, v2, v); + if (err_code) { + MP2_ERR(MP2_SPARSE_MAT_MULT, "Number of rows in the second (dense) matrix must be greater than all of the column indices in the sparse matrix representation.", MP2_SPARSE_MAT_MULT); + } + + for (int l = 0; l < m.rows(); l++) { + for (int j = 0; j < m1.cols(); j++) { + m2.coeffRef(v1[l], j) += m.coeff(l, 0) * m1.coeff(v2[l], j); + } + } + return m2; + + // #' ## Parenthesis // #' // #' The order of operations can be enforced in the usual @@ -1698,7 +1762,6 @@ class ExprEvaluator } else { MP2_ERR(MP2_CBIND, "Inconsistent size in cbind function", MP2_CBIND); - return m; } } return m; @@ -1712,26 +1775,20 @@ class ExprEvaluator int totrows, rowmarker; totrows = 0; rowmarker = 0; - for (int j = 0; j < n; j++) - { + for (int j = 0; j < n; j++) { totrows += args[j].rows(); } m = matrix::Zero(totrows, cols); - for (int i = 0; i < n; i++) - { - if (args[i].cols() == cols) - { + for (int i = 0; i < n; i++) { + if (args[i].cols() == cols) { rows_per_arg = args[i].rows(); - for (int k = 0; k < rows_per_arg; k++) - { + for (int k = 0; k < rows_per_arg; k++) { m.row(rowmarker + k) = args[i].row(k); } rowmarker += rows_per_arg; } - else - { + else { MP2_ERR(MP2_RBIND, "Inconsistent size in rbind function", MP2_RBIND); - return m; } } return m; @@ -1760,7 +1817,6 @@ class ExprEvaluator } else if (size_in < size_out) { if (size_out % size_in) { MP2_ERR(MP2_MATRIX, "The size of the input matrix is not compatible with the requested shape of the output matrix.", MP2_MATRIX); - return m1; } m = mp2_rep(m, size_out / size_in); m.resize(rows, cols); @@ -1904,12 +1960,11 @@ class ExprEvaluator m = args[0].colwise().sum().matrix(); return m; - case MP2_GROUPSUMS: // group_sums(x) + case MP2_GROUPSUMS: // group_sums(x, f, n) v1 = args.get_as_int_vec(1); err_code = args.check_indices(2, v1, {0}); if (err_code) { MP2_ERR(MP2_GROUPSUMS, "Group indexes are out of range.", MP2_GROUPSUMS); - return m; } m = args[0]; if (m.cols() != 1) { @@ -1919,7 +1974,6 @@ class ExprEvaluator m1 = matrix::Zero(rows, 1); if (v1.size() != m.rows()) { MP2_ERR(MP2_GROUPSUMS, "Number of rows in x must equal the number of indices in f in group_sums(x, f, n).", MP2_GROUPSUMS); - return m; } for (int i = 0; i < m.rows(); i++) { m1.coeffRef(v1[i], 0) += m.coeff(i, 0); @@ -1929,7 +1983,6 @@ class ExprEvaluator case MP2_MEAN: // mean(x) if (n != 1) { MP2_ERR(MP2_MEAN, "The mean function can only take a single matrix.", MP2_MEAN); - return m; } m = matrix::Zero(1, 1); sum = args.get_as_mat(0).mean(); @@ -1939,7 +1992,6 @@ class ExprEvaluator case MP2_SD: // sd(x) if (n != 1) { MP2_ERR(MP2_SD, "The sd function can only take a single matrix.", MP2_SD); - return m; } m = args.get_as_mat(0); m1 = matrix::Zero(1, 1); @@ -2059,7 +2111,6 @@ class ExprEvaluator err_code = args.check_indices(0, v1, v2); if (err_code) { MP2_ERR(MP2_SQUARE_BRACKET, "Illegal index to square bracket", MP2_SQUARE_BRACKET); - return m; } m = args[0]; m1 = matrix::Zero(nrow, ncol); @@ -2083,12 +2134,10 @@ class ExprEvaluator err_code = args.check_indices(0, v1, v2); if (err_code){ MP2_ERR(MP2_BLOCK, "Illegal starting index to block", MP2_BLOCK); - return m; } err_code = args.check_indices(0, v3, v4); if (err_code){ MP2_ERR(MP2_BLOCK, "Illegal index to block, requesting more elements than available in input", MP2_BLOCK); - return m; } return args[0].block(rowIndex, colIndex, rows, cols); @@ -2152,14 +2201,12 @@ class ExprEvaluator return m; // have not built up any previous iterations yet, so returning empty matrix if (t == 0) { MP2_ERR(154, "The simulation loop has not yet begun and so rbind_time (or rbind_lag) cannot be used", int_func); - return m; } matIndex = index2mats[0]; // m if ((matIndex < 0) | (index2what[0] != 0)) { MP2_ERR(MP2_RBIND_TIME, "Can only rbind_time (or rbind_lag) named matrices not expressions of matrices and not integer vectors", int_func); // Rcpp::Rcout << "return 2 " << std::endl; - return m; } if ((n == 1) & (!doing_lag)) { @@ -2185,7 +2232,6 @@ class ExprEvaluator // (usually a column vector) of initial values that take // us back into negative time steps. need to do the same // for convolution and anything else that looks backwards. - return m; } } // timeIndex = -timeIndex; @@ -2199,11 +2245,9 @@ class ExprEvaluator lowerTimeBound = args.get_as_int(2); if (lowerTimeBound < 0) { MP2_ERR(MP2_RBIND_TIME, "Lower time bound (third argument) is less than zero", int_func); - return m; } if (lowerTimeBound > t_max) { MP2_ERR(MP2_RBIND_TIME, "Lower time bound (third argument) is greater than the number of time steps", int_func); - return m; } } else if (doing_lag) { @@ -2250,7 +2294,6 @@ class ExprEvaluator { // Shall we allow inconsistent rows? MP2_ERR(MP2_RBIND_TIME, "Inconsistent rows or columns in rbind_time (or rbind_lag)", int_func); // Rcpp::Rcout << "return 5 " << std::endl; - return args[0]; } } @@ -2366,7 +2409,6 @@ class ExprEvaluator if (lag < 0) { MP2_ERR(MP2_TIME_STEP, "Time lag needs to be non-negative", MP2_TIME_STEP); - return m; } if (t > lag) { @@ -2377,7 +2419,6 @@ class ExprEvaluator case MP2_TIME_GROUP: // time_group(i, change_points) if (index2what[0] != 0) { MP2_ERR(MP2_TIME_GROUP, "First argument needs to be a matrix.", MP2_TIME_GROUP); - return m; } m = args[0]; off = args.get_as_int(0); @@ -2391,22 +2432,18 @@ class ExprEvaluator case MP2_TIME_VAR: // time_var(x, change_points) if (t == 0) { MP2_ERR(MP2_TIME_VAR, "Time variation is not allowed before the simulation loop begins.", MP2_TIME_VAR); - return m; } if (t > t_max) { MP2_ERR(MP2_TIME_VAR, "Time variation is not allowed after the simulation loop ends.", MP2_TIME_VAR); - return m; } v = args.get_as_int_vec(1); off = args.get_as_int(1); if (off < 0) { MP2_ERR(MP2_TIME_VAR, "The first element of the second argument must not be less than zero.", MP2_TIME_VAR); - return m; } if (off >= v.size()) { MP2_ERR(MP2_TIME_VAR, "The first element of the second argument must be less than the number of elements in the first.", MP2_TIME_VAR); - return m; } // first argument can have its rows indexed @@ -2428,12 +2465,10 @@ class ExprEvaluator } else { MP2_ERR(MP2_TIME_VAR, "Time variation pointers need to be length-1 integer vectors.", MP2_TIME_VAR); - return m; } } else if (cp < 0) { MP2_ERR(MP2_TIME_VAR, "Negative times are not allowed.", MP2_TIME_VAR); - return m; } } m = matrix::Zero(args[0].cols(), 1); @@ -2507,7 +2542,6 @@ class ExprEvaluator matIndex = index2mats[0]; // m if (matIndex == -1) { MP2_ERR(MP2_CONVOLUTION, "Can only convolve named matrices not expressions of matrices", MP2_CONVOLUTION); - return args[0]; } #ifdef MP_VERBOSE @@ -2550,10 +2584,8 @@ class ExprEvaluator return m; } - else - { + else { MP2_ERR(MP2_CONVOLUTION, "Either empty or non-column vector used as kernel in convolution", MP2_CONVOLUTION); - return m; } // #' ## Clamp @@ -2687,7 +2719,6 @@ class ExprEvaluator case MP2_POISSON_DENSITY: if (n < 2) { MP2_ERR(MP2_POISSON_DENSITY, "dpois needs two arguments: matrices with observed and expected values", MP2_POISSON_DENSITY); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2696,7 +2727,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_POISSON_DENSITY); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) { @@ -2709,7 +2739,6 @@ class ExprEvaluator case MP2_NEGBIN_DENSITY: if (n < 3) { MP2_ERR(MP2_NEGBIN_DENSITY, "dnbinom needs three arguments: matrices with observed values, expected values, and dispersion parameters", MP2_NEGBIN_DENSITY); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2719,7 +2748,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_NEGBIN_DENSITY); - return m; } // var ~ variance // mu ~ mean @@ -2740,7 +2768,6 @@ class ExprEvaluator case MP2_NORMAL_DENSITY: if (n < 3) { MP2_ERR(MP2_NORMAL_DENSITY, "dnorm needs three arguments: matrices with observed values, expected values, and standard deviation parameters", MP2_NORMAL_DENSITY); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2750,7 +2777,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_NORMAL_DENSITY); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) { @@ -2763,7 +2789,6 @@ class ExprEvaluator case MP2_BINOM_DENSITY: if (n < 3) { MP2_ERR(MP2_BINOM_DENSITY, "dbinom needs three arguments: matrices with observed values, numbers of trials, and probabilities", MP2_BINOM_DENSITY); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2773,7 +2798,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_BINOM_DENSITY); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) { @@ -2856,10 +2880,8 @@ class ExprEvaluator return m; case MP2_NEGBIN_SIM: - if (n < 2) - { + if (n < 2) { MP2_ERR(MP2_NEGBIN_SIM, "rnbinom needs two arguments: matrices with means and dispersion parameters", MP2_NEGBIN_SIM); - return m; } eps = 1e-8; rows = args[0].rows(); @@ -2867,10 +2889,8 @@ class ExprEvaluator v1.push_back(1); args = args.recycle_to_shape(v1, rows, cols); err_code = args.get_error_code(); - if (err_code != 0) - { + if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_NEGBIN_SIM); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) @@ -2889,10 +2909,8 @@ class ExprEvaluator return m; case MP2_NORMAL_SIM: - if (n < 2) - { + if (n < 2) { MP2_ERR(MP2_NORMAL_SIM, "rnorm needs two arguments: matrices with means and standard deviations", MP2_NORMAL_SIM); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2900,10 +2918,8 @@ class ExprEvaluator args = args.recycle_to_shape(v1, rows, cols); err_code = args.get_error_code(); // err_code = RecycleInPlace(args[1], rows, cols); - if (err_code != 0) - { + if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_NORMAL_SIM); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) @@ -2919,7 +2935,6 @@ class ExprEvaluator // rbinom(size, prob) if (n != 2) { MP2_ERR(MP2_BINOM_SIM, "rbinom needs two arguments: matrices with size and probability", MP2_BINOM_SIM); - return m; } rows = args[0].rows(); cols = args[0].cols(); @@ -2928,7 +2943,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(MP2_BINOM_SIM, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_BINOM_SIM); - return m; } m = matrix::Zero(rows, cols); for (int i = 0; i < rows; i++) { @@ -2949,7 +2963,6 @@ class ExprEvaluator //Rcpp::Rcout << "++++++" << std::endl; //Rcpp::Rcout << args[0] << std::endl; MP2_ERR(MP2_EULER_MULTINOM_SIM, "The first 'size' argument must be scalar.", MP2_EULER_MULTINOM_SIM); - return m; } if (args[1].cols() != 1) { //Rcpp::Rcout << "------" << std::endl; @@ -3019,7 +3032,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_PGAMMA); - return m; } m1 = args.get_as_mat(0); // q m2 = args.get_as_mat(1); // shape @@ -3041,7 +3053,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_PGAMMA); - return m; } m1 = args.get_as_mat(0); // q m2 = args.get_as_mat(1); // mean @@ -3174,22 +3185,16 @@ class ExprEvaluator case MP2_ASSIGN: cols = args[1].cols(); - if (cols != 1) - { + if (cols != 1) { MP2_ERR(255, "Assignment index matrices must have a single column", MP2_ASSIGN); - return m; } cols = args[2].cols(); - if (cols != 1) - { + if (cols != 1) { MP2_ERR(255, "Assignment index matrices must have a single column", MP2_ASSIGN); - return m; } cols = args[3].cols(); - if (cols != 1) - { + if (cols != 1) { MP2_ERR(255, "Assignment value matrices must have a single column", MP2_ASSIGN); - return m; } // apparently we still need this check @@ -3199,10 +3204,8 @@ class ExprEvaluator // printIntVector(v2); err_code = args.check_indices(0, v1, v2); // CheckIndices(args[0], args[1], m1); // err_code = CheckIndices(args[0], args[1], args[2]); - if (err_code) - { + if (err_code) { MP2_ERR(MP2_ASSIGN, "Illegal index used in assign", MP2_ASSIGN); - return m; } rows = args[3].rows(); @@ -3212,7 +3215,6 @@ class ExprEvaluator err_code = args.get_error_code(); if (err_code != 0) { MP2_ERR(err_code, "cannot recycle rows and/or columns because the input is inconsistent with the recycling request", MP2_ASSIGN); - return m; } for (int k = 0; k < rows; k++) { @@ -3221,7 +3223,6 @@ class ExprEvaluator matIndex = index2mats[0]; if (matIndex == -1) { MP2_ERR(MP2_ASSIGN, "Can only assign to named matrices not expressions of matrices", MP2_ASSIGN); - return args[0]; } valid_vars.m_matrices[matIndex].coeffRef(rowIndex, colIndex) = args[3].coeff(k, 0); } @@ -3280,19 +3281,15 @@ class ExprEvaluator m.resize(size, 1); start = 0; - for (int i = 1; i < n; i++) - { + for (int i = 1; i < n; i++) { sz = args[i].rows() * args[i].cols(); - if (size >= sz) - { + if (size >= sz) { m1 = m.block(start, 0, sz, 1); m1.resize(args[i].rows(), args[i].cols()); // Rcpp::Rcout << "MATRIX " << valid_vars.m_matrices[index2mats[i]] << std::endl << std::endl; matIndex = index2mats[i]; - if (matIndex == -1) - { + if (matIndex == -1) { MP2_ERR(MP2_ASSIGN, "Can only unpack into named matrices not expressions of matrices", MP2_UNPACK); - return args[0]; } valid_vars.m_matrices[matIndex] = m1; // args[i] = m1; @@ -3306,7 +3303,6 @@ class ExprEvaluator default: MP2_ERR(255, "invalid operator in arithmetic expression", -99); - return m; } } // switch (table_n[row]) }; diff --git a/tests/testthat/test-expr-list.R b/tests/testthat/test-expr-list.R index 80c7846bf..497c6b392 100644 --- a/tests/testthat/test-expr-list.R +++ b/tests/testthat/test-expr-list.R @@ -20,10 +20,10 @@ test_that("formula validity is enforced", { ) }) -test_that("proper error message comes when output matrices are not initialized", { +test_that("proper error message when output matrices are not initialized", { m = macpan2:::TMBModel( init_mats = macpan2:::MatsList(a = 1), - expr_list =mp_tmb_expr_list(before = list(b ~ a)) + expr_list = mp_tmb_expr_list(before = list(b ~ a)) ) expect_error( m$data_arg(), diff --git a/tests/testthat/test-expr-parser.R b/tests/testthat/test-expr-parser.R index 9c4f43689..cedecfed7 100644 --- a/tests/testthat/test-expr-parser.R +++ b/tests/testthat/test-expr-parser.R @@ -48,7 +48,7 @@ test_that("parse_tables ...", { ) }) -test_that("argument names cannot be part of macpan2 expressions",{ +test_that("argument names cannot be part of macpan2 expressions", { expect_error(engine_eval(~ matrix(1:6,2,3,byrow = TRUE)) ,"Argument names cannot be used in engine expressions.") @@ -58,3 +58,10 @@ test_that("argument names cannot be part of macpan2 expressions",{ expect_equal(matrix(1:6,2,3,byrow = TRUE), engine_eval(~ x, x = matrix(1:6,2,3,byrow = TRUE))) }) + +test_that("no-op functions are not allowed", { + expect_error( + engine_eval(~ 1 + sin() + A, A = 1) + , regexp = "The following function was called without arguments" + ) +}) diff --git a/tests/testthat/test-sparse-matrix.R b/tests/testthat/test-sparse-matrix.R new file mode 100644 index 000000000..97ecb2028 --- /dev/null +++ b/tests/testthat/test-sparse-matrix.R @@ -0,0 +1,59 @@ +test_that("sparse_mat_mult works properly", { + make_args = function(dn, dm, tol) { + set.seed(1) + n = 10 + m = 5 + p = 3 + A = matrix(rnorm(n * (m + dm)) / 10, n, m + dm) + B = matrix(rnorm(m * p) / 10, m, p) + a = macpan2:::sparse_matrix_notation(A, tol = tol) + traj = simple_sims(list(Y ~ sparse_mat_mult(m, i, j, B, Y)) + , time_steps = 1 + , int_vecs = list( + i = a$row_index + , j = a$col_index + + ) + , mats = list( + m = a$values + , B = B + , Y = matrix(0, n + dn, p) + ) + ) + Y_mp = traj |> filter(matrix == "Y") |> with(matrix(value, max(row) + 1)) + Y_base_r = a$Msparse %*% B + nlist(Y_mp, Y_base_r) + } + + expect_error( + make_args(0, 1, 0.1) + , regexp = "Number of rows in the second" + ) + expect_error( + make_args(-2, 0, 0.1) + , regexp = "Number of rows in the output matrix" + ) + + args = make_args(-1, 0, 0.1) + expect_equal(args$Y_mp, args$Y_base_r[1:9, ]) + expect_equal(dim(args$Y_mp), c(9, 3)) + + args = make_args(0, 0, 0.1) + expect_equal(args$Y_mp, args$Y_base_r) + + args = make_args(0, 0, 0) + expect_equal(args$Y_mp, args$Y_base_r) + + expect_error( + engine_eval(~sparse_mat_mult(matrix(0, 2, 2), 0, 0, 0, 0)) + , "First matrix must be represented as a vector with associated row and column indices" + ) + expect_error( + engine_eval(~sparse_mat_mult(0, 0, 0)) + , "Require exactly 5 arguments, but got fewer" + ) + expect_error( + engine_eval(~sparse_mat_mult(0, 0, 0, 0, 0, 0)) + , "Require exactly 5 arguments, but got more" + ) +}) From 897a18887587415ff75c1f3ac37a2b30572a5b6b Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Mon, 14 Jul 2025 11:44:07 -0400 Subject: [PATCH 02/13] fix broken test --- tests/testthat/test-dot-layout.R | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-dot-layout.R b/tests/testthat/test-dot-layout.R index b4bda74f9..62a0b50d4 100644 --- a/tests/testthat/test-dot-layout.R +++ b/tests/testthat/test-dot-layout.R @@ -1,9 +1,11 @@ test_that("dot layouts produce appropriate errors and output classes", { if (require(Rgraphviz, quietly = TRUE)) { specs <- mp_tmb_entire_library() - no_mpflows <- c("lotka_volterra_competition", - "lotka_volterra_predator_prey", - "nfds") + n_flows = (specs + |> lapply(mp_flow_frame, topological_sort = FALSE) + |> vapply(nrow, integer(1L)) + ) + no_mpflows <- which(n_flows == 0L) |> names() for (s in no_mpflows) { expect_error(dot_layout(specs[[s]]), "was spec defined") } From 565e317277af0f4ce61b5186da71cc578d32c24f Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Mon, 14 Jul 2025 12:27:24 -0400 Subject: [PATCH 03/13] documenting sparse matrix multiplication --- NAMESPACE | 1 + R/index_matrices.R | 33 ++++- _pkgdown.yml | 2 + man/sparse_matrix_notation.Rd | 43 ++++++ vignettes/matrix_multiplication.Rmd | 206 ++++++++++++++++++++++++++++ 5 files changed, 284 insertions(+), 1 deletion(-) create mode 100644 man/sparse_matrix_notation.Rd create mode 100644 vignettes/matrix_multiplication.Rmd diff --git a/NAMESPACE b/NAMESPACE index 097559b05..3dd9d5a35 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -366,6 +366,7 @@ export(show_models) export(si_example_code) export(si_example_object) export(simple_sims) +export(sparse_matrix_notation) export(to_labels) export(to_name) export(to_name_pairs) diff --git a/R/index_matrices.R b/R/index_matrices.R index 7af715b7c..3f4a2ad5d 100644 --- a/R/index_matrices.R +++ b/R/index_matrices.R @@ -35,7 +35,38 @@ binary_matrix_notation <- function(M){ } - +#' Extract Sparse Matrix Notation from a Dense Matrix +#' +#' Converts a dense matrix to a sparse representation by extracting its non-zero +#' entries and their indices. Entries with absolute value less than or equal to +#' `tol` are treated as zeros. +#' +#' @param M A numeric matrix or object coercible to a matrix. +#' @param zero_based Logical; if `TRUE` (default), the returned row and column +#' indices are zero-based (starting at 0). If `FALSE`, indices are one-based +#' (as in standard R matrices). +#' @param tol Numeric tolerance (default `1e-4`). Entries with absolute value +#' less than or equal to `tol` are treated as zero. +#' +#' @return A named list with components: +#' \describe{ +#' \item{`row_index`}{Integer vector of row indices for non-zero entries +#' (adjusted by `zero_based`).} +#' \item{`col_index`}{Integer vector of column indices for non-zero entries +#' (adjusted by `zero_based`).} +#' \item{`values`}{Numeric vector of the non-zero entries.} +#' \item{`M`}{The original input matrix, coerced to a dense matrix.} +#' \item{`Msparse`}{A copy of `M` with near-zero entries (as determined by +#' `tol`) replaced by exact zeros.} +#' } +#' +#' @examples +#' M <- matrix(c(5, 0, 0, +#' 0, 0, 3, +#' 0, 2, 0), nrow = 3, byrow = TRUE) +#' sparse_matrix_notation(M) +#' +#' @export sparse_matrix_notation = function(M, zero_based = TRUE, tol = 1e-4) { M = as.matrix(M) non_zero_loc = abs(M) > tol diff --git a/_pkgdown.yml b/_pkgdown.yml index 12c5ce151..217b68be4 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -184,6 +184,7 @@ reference: - empty_trajectory - mp_zero_vector - rbf + - sparse_matrix_notation - mp_binary_operator - subtitle: Developer and Power-User Utilities @@ -252,6 +253,7 @@ articles: - calibration - real_data - options + - matrix_multiplication - FAQs - title: Specs desc: Specification Documents diff --git a/man/sparse_matrix_notation.Rd b/man/sparse_matrix_notation.Rd new file mode 100644 index 000000000..2ea2b060c --- /dev/null +++ b/man/sparse_matrix_notation.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/index_matrices.R +\name{sparse_matrix_notation} +\alias{sparse_matrix_notation} +\title{Extract Sparse Matrix Notation from a Dense Matrix} +\usage{ +sparse_matrix_notation(M, zero_based = TRUE, tol = 1e-04) +} +\arguments{ +\item{M}{A numeric matrix or object coercible to a matrix.} + +\item{zero_based}{Logical; if \code{TRUE} (default), the returned row and column +indices are zero-based (starting at 0). If \code{FALSE}, indices are one-based +(as in standard R matrices).} + +\item{tol}{Numeric tolerance (default \code{1e-4}). Entries with absolute value +less than or equal to \code{tol} are treated as zero.} +} +\value{ +A named list with components: +\describe{ +\item{\code{row_index}}{Integer vector of row indices for non-zero entries +(adjusted by \code{zero_based}).} +\item{\code{col_index}}{Integer vector of column indices for non-zero entries +(adjusted by \code{zero_based}).} +\item{\code{values}}{Numeric vector of the non-zero entries.} +\item{\code{M}}{The original input matrix, coerced to a dense matrix.} +\item{\code{Msparse}}{A copy of \code{M} with near-zero entries (as determined by +\code{tol}) replaced by exact zeros.} +} +} +\description{ +Converts a dense matrix to a sparse representation by extracting its non-zero +entries and their indices. Entries with absolute value less than or equal to +\code{tol} are treated as zeros. +} +\examples{ +M <- matrix(c(5, 0, 0, + 0, 0, 3, + 0, 2, 0), nrow = 3, byrow = TRUE) +sparse_matrix_notation(M) + +} diff --git a/vignettes/matrix_multiplication.Rmd b/vignettes/matrix_multiplication.Rmd new file mode 100644 index 000000000..c0e11ba8a --- /dev/null +++ b/vignettes/matrix_multiplication.Rmd @@ -0,0 +1,206 @@ +--- +title: "Multiplying a Sparse Matrix with a Dense Matrix" +output: + rmarkdown::html_vignette: + toc: true +vignette: > + %\VignetteIndexEntry{Multiplying a Sparse Matrix with a Dense Matrix} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} +editor_options: + chunk_output_type: console +--- + +[![status](https://img.shields.io/badge/status-mature%20draft-yellow)](https://canmod.github.io/macpan2/articles/vignette-status#mature-draft) + +```{r settings, echo = FALSE, eval = TRUE, message=FALSE, warning=FALSE, error=FALSE} +library(macpan2) +library(splines) +library(dplyr) +library(ggplot2) +theme_bw = function() ggplot2::theme_bw(base_size = 18) +``` + +## Overview + +When working with models in `macpan2`, it is often necessary to perform matrix multiplications where the left-hand matrix is sparse—meaning most of its entries are zero. Passing large sparse matrices as dense matrices with zeros into `macpan2` models can lead to unnecessary memory use and performance slowdowns. + +To address this, `macpan2` provides the engine function `sparse_mat_mult()`, which multiplies a sparse matrix (given in compressed form) by a dense matrix. Instead of passing the full sparse matrix, you pass: + +- `x` — the non-zero values of the sparse matrix, as a numeric vector +- `i` — the zero-based row indices of each value in `x` +- `j` — the zero-based column indices of each value in `x` +- `y` — the dense matrix on the right-hand side of the multiplication +- `z` — a matrix with the correct dimensions, which is passed both as an input and as the destination for the result + +This operation is equivalent to: +```r +z ~ LeftSparseMatrix %*% y +``` +but avoids ever constructing `LeftSparseMatrix` explicitly. You pass `z` as a pre-allocated matrix of the correct size, and `sparse_mat_mult()` fills it with the result of the multiplication. + +This method is particularly useful in `macpan2` for: + +- Reducing memory footprint when working with sparse structures +- Avoiding unnecessary computation with zero entries +- Embedding efficient matrix multiplications in a model specification + +We also plan to allow passing sparse matrices from the `Matrix` package directly into `macpan2` models, in the same way integer vectors can be passed today. The `Matrix` package works naturally with `TMB`, so this will allow us to leverage native sparse matrix support inside the `macpan2` engine. +Even with our current simple integer-vector-based approach, we have observed practically useful performance gains for large problems, especially when the left-hand matrix is highly sparse. + +## Extracting Sparse Representation + +To simplify the creation of sparse representations, `macpan2` provides the helper function `sparse_matrix_notation()`. +This function takes a dense matrix and returns: + +- `values` — the non-zero entries of the matrix +- `row_index` — the zero-based row indices of those entries +- `col_index` — the zero-based column indices of those entries +- `M` — the original matrix +- `Msparse` — a version of `M` with near-zero entries replaced by exact zeros + +By default, `sparse_matrix_notation()` applies a numerical tolerance to treat small values as zero and converts indices to zero-based, matching the requirements of `sparse_mat_mult()`. + +Example usage: +```{r} +library(macpan2) + +# Define matrix to convert +M <- matrix(c(5, 0, 0, + 0, 0, 3, + 0, 2, 0), nrow = 3, byrow = TRUE) + +# Extract sparse representation +sparse <- macpan2:::sparse_matrix_notation(M, zero_based = TRUE) +print(sparse) + +# Calculate sparsity percentage +total_entries <- length(M) +non_zero_entries <- length(sparse$values) +sparsity <- 100 * (1 - non_zero_entries / total_entries) +cat(sprintf("Matrix sparsity: %.1f%%\n", sparsity)) +``` + +With this simple example, the matrix is 66.7% sparse. +The performance benefit of using sparse matrix multiplication grows with both matrix size and sparsity percentage. + +## Example: Small Sparse Matrix Multiplication + +```{r} +library(macpan2) + +# Full matrix for validation +A <- matrix(c(5, 0, 0, + 0, 0, 3, + 0, 2, 0), nrow = 3, byrow = TRUE) + +# Extract sparse representation +sparse <- macpan2:::sparse_matrix_notation(A, zero_based = TRUE) + +# Dense matrix to multiply (3×2) +y <- matrix(c(1, 4, + 2, 5, + 3, 6), nrow = 3, byrow = TRUE) + +# Pre-allocate output matrix (3×2) +z <- matrix(0, nrow = 3, ncol = 2) + +# Model specification with sparse_mat_mult in before block +spec <- mp_tmb_model_spec( + before = z ~ sparse_mat_mult(x, i, j, y, z), + default = nlist(x = sparse$values, y, z), + integers = list(i = sparse$row_index, j = sparse$col_index) +) + +# Run simulation for zero time steps to execute 'before' block +result <- spec |> + mp_simulator(time_steps = 0, outputs = "z") |> + mp_final_list() + +# Result from macpan2 +result$z + +# Direct multiplication for validation +A %*% y +``` + +- `sparse_mat_mult()` multiplies a sparse matrix (given as values and indices) by a dense matrix without constructing the full sparse matrix. +- The helper `sparse_matrix_notation()` extracts the required sparse representation from any matrix, applying zero-based indexing and a numerical tolerance. +- The function writes the result into `z`, which must be allocated beforehand. +- You can use `mp_tmb_model_spec()` with `time_steps = 0` to evaluate this operation as part of a `macpan2` model workflow. + +## Example: Time-Varying Transmission in an SI Model + +We implement an SI model where the transmission rate $\beta(t)$ varies over time, modeled as a spline basis expansion: +$$ +\log \beta(t) = B(t) \cdot \theta +$$ +We use `sparse_matrix_notation()` to encode the spline basis sparsely and compute $\log \beta(t)$ using `sparse_mat_mult()` in the `before` block. In this case, the spline basis contains no exact zeros, but many near-zero values that contribute little to the result. The tolerance ensures these negligible entries are dropped, improving sparsity without meaningfully affecting the result. For basis matrices like these, choosing a reasonable tolerance (e.g., `1e-4`) helps strike a balance between computational efficiency and model fidelity. + +```{r, fig.height = 8, fig.width=4} +library(macpan2) +library(splines) +library(dplyr) +library(ggplot2) + +# Time points for simulation (100 steps) +n_steps <- 100 +time_points <- seq(0, 1, length.out = n_steps) + +# Create a cubic B-spline basis with 8 degrees of freedom +B_dense <- bs(time_points, df = 8, degree = 3, intercept = TRUE) + +# Convert basis to sparse form with small tolerance +B_sparse <- macpan2:::sparse_matrix_notation(B_dense, zero_based = TRUE, tol = 1e-4) + +# Print basis sparsity +total_entries <- length(B_dense) +non_zero_entries <- length(B_sparse$values) +sparsity <- 100 * (1 - non_zero_entries / total_entries) +cat(sprintf("Spline basis matrix sparsity (tol = 1e-4): %.1f%%\n", sparsity)) + +# Example spline coefficients for log(beta) +set.seed(12) +theta <- rnorm(8, -2, 1) + +# Pre-allocate beta_log[t] for 100 time steps +beta_log <- matrix(0, nrow = n_steps, ncol = 1) + +# Define SI model with time-varying beta +spec <- mp_tmb_model_spec( + before = beta_log ~ sparse_mat_mult(x, i, j, theta, beta_log), + during = list( + beta ~ exp(beta_log[time_step(1)]), + mp_per_capita_flow( + from = "S", + to = "I", + rate = "beta * I / N", + flow_name = "infection" + ) + ), + default = nlist(N = 100, S = 99, I = 1, x = B_sparse$values, theta, beta_log), + integers = nlist(i = B_sparse$row_index, j = B_sparse$col_index) +) + +# Simulate for 100 time steps +outputs = c("S", "I", "beta", "infection") +result <- spec |> + mp_simulator(time_steps = n_steps, outputs = outputs) |> + mp_trajectory() |> + mutate(variable = factor(matrix, levels = outputs)) + +# Plot the results +(result + |> ggplot() + + aes(time, value) + + geom_line() + + facet_wrap(~variable, scales = "free_y", ncol = 1) + + theme_bw() +) +``` + +- Transmission rate $\beta(t)$ is computed dynamically using `sparse_mat_mult()`. +- The spline basis is passed in sparse form, allowing efficient calculation of $\beta(t)$ at each time step. +- The SI model uses `mp_per_capita_flow()` for infection, with time-varying transmission. +- For this example, the spline basis (with tolerance $10^{-4}$) has approximately `r sprintf("%.1f%% sparsity", sparsity)`. + Larger models with higher sparsity will benefit even more from this approach. From 1ce47c6ea40c728e66619ce85bc0d1e138611868 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Mon, 14 Jul 2025 14:17:59 -0400 Subject: [PATCH 04/13] using sparse --- Makefile | 2 +- R/tmb_model_editors.R | 43 ++++++++++++++++++++++++----- _pkgdown.yml | 4 +-- man/mp_tmb_insert_glm_timevar.Rd | 6 ++-- tests/testthat/test-change-points.R | 2 ++ 5 files changed, 45 insertions(+), 12 deletions(-) diff --git a/Makefile b/Makefile index 723b409d0..aab6e5c98 100644 --- a/Makefile +++ b/Makefile @@ -34,7 +34,7 @@ full-install: # haven't updated macpan.cpp (but have perhaps modified dev.cpp) # and (3) do not require a roxygen update. quick-install: enum-update enum-meth-update - R CMD INSTALL --no-multiarch --install-tests --configure-args="CFLAGS=-g CXXFLAGS=-g" . + R CMD INSTALL --no-multiarch --install-tests --configure-args="CFLAGS=-g -O0 CXXFLAGS=-g -O0" . quick-doc-install: R/*.R misc/dev/dev.cpp diff --git a/R/tmb_model_editors.R b/R/tmb_model_editors.R index 41cf1fb3f..8313f525d 100644 --- a/R/tmb_model_editors.R +++ b/R/tmb_model_editors.R @@ -288,7 +288,8 @@ mp_tmb_insert_reports = function(model #' @param parameter_name Character string giving the name of the parameter #' to make time-varying. #' @param design_matrix Matrix of time variation. -#' @param timevar_coef Initial coefficient vector. +#' @param timevar_coef Initial coefficient matrix with the same number of rows +#' as there are columns of the `design_matrix`. #' @param link_function Link function given by functions like #' \code{\link{mp_log}}. #' @param matrix_coef_name Name of the vector containing values of the non-zero @@ -322,13 +323,30 @@ mp_tmb_insert_glm_timevar = function(model , timevar_coef_name = sprintf("time_var_%s", parameter_name) , time_index_name = sprintf("time_index_%s", parameter_name) , sparsity_tolerance = 0 + , engine_function = c("sparse_mat_mult", "group_sums") ) { + if (!parameter_name %in% names(model$default)) { + sprintf( + "Parameter '%s' not found in model$default. You must define it before adding time variation." + , parameter_name + ) |> stop() + } + + design_matrix = as.matrix(design_matrix) + if (nrow(design_matrix) < 1L) stop("The design matrix must have at least one row.") + if (ncol(design_matrix) != nrow(timevar_coef)) { + stop("The design_matrix and timevar_coef matrix are not compatible. The number of columns in the former must equal the number of rows in the latter.") + } + + sparse_matrix = sparse_matrix_notation(design_matrix, tol = sparsity_tolerance) + engine_function = match.arg(engine_function) + matrix_coefs = sparse_matrix$values matrix_row = sparse_matrix$row_index matrix_col = sparse_matrix$col_index - linear_pred = timeseries = numeric(nrow(design_matrix)) + linear_pred = timeseries = matrix(0, nrow(design_matrix), ncol(timevar_coef)) time_var = list(timevar_coef) |> setNames(timevar_coef_name) time_index = seq_along(linear_pred) @@ -342,11 +360,22 @@ mp_tmb_insert_glm_timevar = function(model , c(matrix_row_name, matrix_col_name, time_index_name) ) - matmult = sprintf("%s ~ %s + group_sums(%s * %s[%s], %s, %s)" - , linear_pred_name, link_function$ref(parameter_name) - , matrix_coef_name, timevar_coef_name - , matrix_col_name, matrix_row_name, linear_pred_name - ) |> as.formula() + if (engine_function == "group_sums") { + matmult = sprintf("%s ~ %s + group_sums(%s * %s[%s], %s, %s)" + , linear_pred_name, link_function$ref(parameter_name) + , matrix_coef_name, timevar_coef_name + , matrix_col_name, matrix_row_name, linear_pred_name + ) |> as.formula() + } else if(engine_function == "sparse_mat_mult") { + matmult = sprintf("%s ~ %s + sparse_mat_mult(%s, %s, %s, %s, %s)" + , linear_pred_name + , link_function$ref(parameter_name) + , matrix_coef_name, matrix_row_name, matrix_col_name + , timevar_coef_name + , linear_pred_name + ) |> as.formula() + } + linkfn = sprintf("%s ~ %s" , timeseries_name , link_function$ref_inv(linear_pred_name) diff --git a/_pkgdown.yml b/_pkgdown.yml index 217b68be4..af3f27a21 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -253,7 +253,6 @@ articles: - calibration - real_data - options - - matrix_multiplication - FAQs - title: Specs desc: Specification Documents @@ -263,8 +262,9 @@ articles: - title: Superuser Guides desc: Articles aimed at advanced users contents: - - calibration_advanced + - matrix_multiplication - time_varying_parameters_advanced + - calibration_advanced - engine_agnostic_grammar - title: Developer Guides desc: Articles aimed at package developers diff --git a/man/mp_tmb_insert_glm_timevar.Rd b/man/mp_tmb_insert_glm_timevar.Rd index 93e1954c2..ebdad2fbd 100644 --- a/man/mp_tmb_insert_glm_timevar.Rd +++ b/man/mp_tmb_insert_glm_timevar.Rd @@ -17,7 +17,8 @@ mp_tmb_insert_glm_timevar( timeseries_name = sprintf("timeseries_\%s", parameter_name), timevar_coef_name = sprintf("time_var_\%s", parameter_name), time_index_name = sprintf("time_index_\%s", parameter_name), - sparsity_tolerance = 0 + sparsity_tolerance = 0, + engine_function = c("sparse_mat_mult", "group_sums") ) } \arguments{ @@ -28,7 +29,8 @@ to make time-varying.} \item{design_matrix}{Matrix of time variation.} -\item{timevar_coef}{Initial coefficient vector.} +\item{timevar_coef}{Initial coefficient matrix with the same number of rows +as there are columns of the \code{design_matrix}.} \item{link_function}{Link function given by functions like \code{\link{mp_log}}.} diff --git a/tests/testthat/test-change-points.R b/tests/testthat/test-change-points.R index cb5f956d3..2edcc6ed6 100644 --- a/tests/testthat/test-change-points.R +++ b/tests/testthat/test-change-points.R @@ -26,3 +26,5 @@ test_that("time_var grabs entire matrix rows each time step", { answer = A[rep(1:3, 1:3), ] |> t() |> as.vector() expect_equal(result, answer) }) + + From 7f8685dcc118a453b1610dcc20da8f6689a77332 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Tue, 15 Jul 2025 16:31:19 -0400 Subject: [PATCH 05/13] roxygen and tests --- R/tmb_model_editors.R | 5 ++ man/mp_tmb_insert_glm_timevar.Rd | 6 ++ tests/testthat/test-glm-timevar.R | 142 ++++++++++++++++++++++++++++++ 3 files changed, 153 insertions(+) create mode 100644 tests/testthat/test-glm-timevar.R diff --git a/R/tmb_model_editors.R b/R/tmb_model_editors.R index 8313f525d..1d0d39284 100644 --- a/R/tmb_model_editors.R +++ b/R/tmb_model_editors.R @@ -308,6 +308,11 @@ mp_tmb_insert_reports = function(model #' parameter changes. #' @param sparsity_tolerance Make design matrix coefficients exactly zero #' when they are below this tolerance. +#' @param engine_function Which of two \code{\link{engine_functions}}, +#' \code{\link{sparse_mat_mult}} or \code{\link{group_sums}}, should be +#' used to compute the matrix multiplication. The only reason to use +#' anything other than the default is for back-compatibility with an +#' existing script. #' #' @export mp_tmb_insert_glm_timevar = function(model diff --git a/man/mp_tmb_insert_glm_timevar.Rd b/man/mp_tmb_insert_glm_timevar.Rd index ebdad2fbd..10a59d08a 100644 --- a/man/mp_tmb_insert_glm_timevar.Rd +++ b/man/mp_tmb_insert_glm_timevar.Rd @@ -58,6 +58,12 @@ parameter changes.} \item{sparsity_tolerance}{Make design matrix coefficients exactly zero when they are below this tolerance.} + +\item{engine_function}{Which of two \code{\link{engine_functions}}, +\code{\link{sparse_mat_mult}} or \code{\link{group_sums}}, should be +used to compute the matrix multiplication. The only reason to use +anything other than the default is for back-compatibility with an +existing script.} } \description{ Insert GLM Time Variation diff --git a/tests/testthat/test-glm-timevar.R b/tests/testthat/test-glm-timevar.R new file mode 100644 index 000000000..95ee5dd6b --- /dev/null +++ b/tests/testthat/test-glm-timevar.R @@ -0,0 +1,142 @@ +test_that("Basic insertion with identity link works", { + design_matrix <- diag(3) + timevar_coef <- matrix(0.1, nrow = 3, ncol = 1) + + expect_error( + mp_tmb_insert_glm_timevar( + mp_tmb_model_spec(), + parameter_name = "beta", + design_matrix = design_matrix, + timevar_coef = timevar_coef, + link_function = mp_identity + ) + , regexp = "not found in model" + ) + out = mp_tmb_insert_glm_timevar( + mp_tmb_model_spec(default = list(beta = 0.2)), + parameter_name = "beta", + design_matrix = design_matrix, + timevar_coef = timevar_coef, + link_function = mp_identity # Assuming identity link for simplicity + ) + + # Should have inserted default variables + expect_true("matrix_coef_beta" %in% names(out$default)) + expect_true("linear_pred_beta" %in% names(out$default)) + expect_true("time_var_beta" %in% names(out$default)) + expect_true("timeseries_beta" %in% names(out$default)) + + # Should have inserted integer variables + expect_true("matrix_row_beta" %in% names(out$integers)) + expect_true("matrix_col_beta" %in% names(out$integers)) + expect_true("time_index_beta" %in% names(out$integers)) + + # Should insert expressions into 'before' and 'during' phases + expect_true(any(grepl("linear_pred_beta", out$before))) + expect_true(any(grepl("timeseries_beta", out$before))) + expect_true(any(grepl("beta ~ time_var", out$during))) +}) + +test_that("Sparse matrix conversion zeros small elements", { + model <- mp_tmb_model_spec(default = list(gamma = 0.2)) + design_matrix <- matrix(c(0.001, 0.1, 0, 0), nrow = 2) + timevar_coef <- matrix(0.2, nrow = 2, ncol = 1) + + tol = 0.01 + out <- mp_tmb_insert_glm_timevar( + model, + parameter_name = "gamma", + design_matrix = design_matrix, + timevar_coef = timevar_coef, + sparsity_tolerance = tol + ) + + matrix_coefs <- out$default[["matrix_coef_gamma"]] + expect_false(any(abs(matrix_coefs) < tol)) +}) + +test_that("Handles empty design matrix gracefully", { + model <- mp_tmb_model_spec(default = list(empty = 1)) + + expect_error( + mp_tmb_insert_glm_timevar( + model, + "empty", + empty_matrix, + empty_matrix + ) + , regexp = "The design matrix must have at least one row" + ) + + expect_error( + mp_tmb_insert_glm_timevar( + model, + "empty", + matrix(0, 3, 2), + matrix(0, 3, 2) + ) + , regexp = "The design_matrix and timevar_coef matrix " + ) +}) + +test_that("Link function is correctly referenced", { + model <- mp_tmb_model_spec(default = list(phi = 1)) + design_matrix <- diag(1) + timevar_coef <- matrix(1, nrow = 1, ncol = 1) + + out <- mp_tmb_insert_glm_timevar( + model, + "phi", + design_matrix, + timevar_coef, + link_function = mp_log + ) + + expect_true(any(grepl("log\\(phi\\)", out$before))) +}) + +test_that("Time-varying beta affects SI model dynamics as expected", { + # Define time points and design matrix for linear trend + t_steps <- 1:50 + X <- cbind(1, t_steps) # Intercept + linear trend + + # Coefficients: intercept = log(0.05), slope = 0.02 + coef_mat <- matrix(c(log(0.1), 0.04), ncol = 1) + + # Build SI model spec with default beta + model <- mp_tmb_model_spec( + before = S ~ N - I + , during = mp_per_capita_flow("S", "I", "beta * I / N", "infection") + , default = list(S = 999, I = 1, N = 1000, beta = 1) + ) + + # Insert time-varying beta (log link) + model_tv <- mp_tmb_insert_glm_timevar(mp_hazard(model) + , parameter_name = "beta" + , design_matrix = X + , timevar_coef = coef_mat + , link_function = mp_log + ) + + # Run simulation + sims <- (model_tv + |> mp_simulator( + time_steps = max(t_steps) + , outputs = c("beta", "infection") + ) + |> mp_trajectory() + ) + + actual_beta = sims |> filter(matrix == "beta") |> pull(value) + + # Check that beta increases over time as per the design matrix + expected_beta <- exp(X %*% coef_mat) |> as.numeric() + expect_equal(actual_beta, expected_beta, tolerance = 1e-6) + + # (sims + # |> ggplot() + # + aes(time, value) + # + geom_line() + # + facet_wrap(~matrix, ncol = 1, scales = 'free') + # ) +}) From 9d0b36a89237794bfb9a0d77f996bc6a26ffa965 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Tue, 15 Jul 2025 16:34:47 -0400 Subject: [PATCH 06/13] [skip ci] bump version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index e85fc8189..82a1b0397 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: macpan2 Title: Fast and Flexible Compartmental Modelling -Version: 3.1.0 +Version: 3.2.0 Authors@R: c( person("Steve Walker", email="swalk@mcmaster.ca", role=c("cre", "aut")), person("Weiguang Guan", role="aut"), From add6b5a316506f9b9e854708037e23d45a4b6e46 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Wed, 16 Jul 2025 15:32:46 -0400 Subject: [PATCH 07/13] bug fixes and methods --- R/lists.R | 6 +++--- R/mp_tmb_model_spec.R | 6 ++++++ 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/R/lists.R b/R/lists.R index ce5de0bc7..93cccc428 100644 --- a/R/lists.R +++ b/R/lists.R @@ -119,8 +119,8 @@ melt_list_of_char_vecs = function(x) { clean_dimnames = function(dn) { if (!is.null(dn)) { - if (identical(as.character(dn[[2L]]), "0")) dn[[2L]] = NULL - if (identical(as.character(dn[[1L]]), "0")) dn[[1L]] = NULL + if (identical(as.character(dn[[2L]]), "0")) dn[2L] = list(NULL) + if (identical(as.character(dn[[1L]]), "0")) dn[1L] = list(NULL) } return(dn) } @@ -133,7 +133,7 @@ cast_default_matrix_list = function(x) { ncol = vapply(col_list, length, integer(1L)) mapply(matrix, val_list, nrow, ncol, dimnames = dimnames, SIMPLIFY = FALSE, USE.NAMES = TRUE) } -cast_default_matrix_list = memoise(cast_default_matrix_list) +#cast_default_matrix_list = memoise(cast_default_matrix_list) empty_named_list = function() list() |> setNames(character(0L)) diff --git a/R/mp_tmb_model_spec.R b/R/mp_tmb_model_spec.R index 280becdf1..63edf4318 100644 --- a/R/mp_tmb_model_spec.R +++ b/R/mp_tmb_model_spec.R @@ -54,6 +54,12 @@ TMBModelSpec = function( self$all_formula_vars = function() { self$expr_list()$all_formula_vars() } + self$all_lhs_vars = function() { + self$expr_list()$all_formula_vars("left") + } + self$all_rhs_vars = function() { + self$expr_list()$all_formula_vars("right") + } self$all_default_mats = function() { setdiff( self$all_default_vars() From f90b48790ab0f7abb476ed8417f8aa65ea1486a7 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Fri, 25 Jul 2025 10:46:31 -0400 Subject: [PATCH 08/13] kroncker work and associated fixes --- DESCRIPTION | 2 +- NAMESPACE | 2 + R/binary_operator.R | 59 +++ R/engine_functions.R | 102 +++-- R/enum.R | 3 +- R/mp_tmb_model_spec.R | 17 +- man/engine_functions.Rd | 112 +++-- misc/dev/dev.cpp | 204 +++++---- .../general_compartmental_model.Rmd | 2 +- src/macpan2.cpp | 204 +++++---- tests/testthat/test-clamp.R | 8 +- vignettes/kronecker.Rmd | 419 ++++++++++++++++++ vignettes/matrix_multiplication.Rmd | 328 +++++++++++--- vignettes/state_dependent_rates.Rmd | 191 +++++++- 14 files changed, 1316 insertions(+), 337 deletions(-) create mode 100644 vignettes/kronecker.Rmd diff --git a/DESCRIPTION b/DESCRIPTION index 82a1b0397..ed889704c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -65,5 +65,5 @@ Config/Needs/website: canmod/macpan2 Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.2.9000 +RoxygenNote: 7.3.2 URL: https://canmod.github.io/macpan2/, https://github.com/canmod/macpan2 diff --git a/NAMESPACE b/NAMESPACE index a1f2eb257..8cce770b4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -259,6 +259,7 @@ export(mp_inflow) export(mp_initial) export(mp_initial_list) export(mp_join) +export(mp_kronecker_operator) export(mp_labels) export(mp_layout_grid) export(mp_layout_paths) @@ -313,6 +314,7 @@ export(mp_simulator) export(mp_slices) export(mp_sqrt) export(mp_square) +export(mp_square_operator) export(mp_state_dependence_frame) export(mp_state_flow_vars) export(mp_state_vars) diff --git a/R/binary_operator.R b/R/binary_operator.R index 8b5451cf2..08471908d 100644 --- a/R/binary_operator.R +++ b/R/binary_operator.R @@ -59,6 +59,42 @@ BinaryOperator = function(operator) { } } +KroneckerOperator = function(operator) { + function(x, y) { + z = kronecker(as.matrix(x), as.matrix(y), FUN = operator, make.dimnames = TRUE) + + clean = function(dn) { + if (!is.null(dn)) { + if (identical(as.character(dn[[2L]]), ":")) dn[2L] = list(NULL) + if (identical(as.character(dn[[1L]]), ":")) dn[1L] = list(NULL) + } + return(dn) + } + ## handle difference of opinion about dimnames + + dimnames(z) = clean(dimnames(z)) + rownames(z) = gsub(":", ".", rownames(z)) + colnames(z) = gsub(":", ".", colnames(z)) + rownames(z) = sub("^\\.", "", rownames(z)) + rownames(z) = sub("\\.$", "", rownames(z)) + colnames(z) = sub("^\\.", "", colnames(z)) + colnames(z) = sub("\\.$", "", colnames(z)) + + if (ncol(z) == 1L) z = setNames(c(z), rownames(z)) + return(z) + } +} + +SquareOperator = function(operator) { + function(x) { + if (is.matrix(x)) stop("can only produce a square diagonal matrix for a vector") + if (length(x) == 1) return(x) + nms = names(x) + x = operator(x) + dimnames(x) = list(nms, nms) + return(x) + } +} #' Binary Operator #' @@ -85,4 +121,27 @@ BinaryOperator = function(operator) { #' @export mp_binary_operator = BinaryOperator +#' Kronecker Operator +#' +#' Convert a function that represents a scalar binary +#' operator into one that computes the \code{\link{kronecker}} +#' version, but with dimensions named in a way that is more +#' convenient for use with `macpan2`. +#' +#' @param operator A scalar binary operator. +#' @return A Kronecker operator convenient for use with `macpan2`. +#' +#' @export +mp_kronecker_operator = KroneckerOperator +#' Square Matrix Operator +#' +#' Convert a unary operator that takes a vector, into one with +#' dimensions named in a way that is more convenient for use with +#' `macpan2`. +#' +#' @param operator A unary operator of a vector. +#' @return An operator. +#' +#' @export +mp_square_operator = SquareOperator diff --git a/R/engine_functions.R b/R/engine_functions.R index eb5e9e185..9d4f48124 100644 --- a/R/engine_functions.R +++ b/R/engine_functions.R @@ -115,11 +115,15 @@ #' engine_eval(~ log(y), y = c(2, 0.5)) #' ``` #' -#' ## Proportions +#' ## Limiting Values #' +#' +#' #' ### Functions #' #' * `proportions(x, limit, eps)` +#' * `limit(x, limit, eps)` +#' * `clamp(x, eps, limit)` #' #' ### Arguments #' @@ -131,13 +135,64 @@ #' #' * matrix of `x / sum(x)` or `rep(limit, length(x))` if #' `sum(x) < eps`. -#' +#' #' ### Examples #' #' ``` #' engine_eval(~ proportions(y, 0.5, 1e-8), y = c(2, 0.5)) #' ``` #' +#' ### Details +#' +#' The return value depends on the function. +#' * `proportions` : matrix of `x / sum(x)` or `rep(limit, length(x))` if +#' `sum(x) < eps`. +#' * `limit` : +#' * `clamp` : +#' ### Details +#' +#' The divide_safe() function conducts elementwise division +#' of the first two arguments, and then allows two other +#' arguments: +#' * `limit` : Value to return if the `sqrt(denominator^2)` is +#' ### Details +#' +#' The `clamp` function smoothly clamps the elements of a +#' matrix so that they +#' do not get closer to 0 than a tolerance, `eps`, with +#' a default of 1e-12. This `clamp` function is the following +#' modification of the +#' [squareplus function](https://arxiv.org/abs/2112.11687). +#' +#' \deqn{f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}} +#' +#' Where the two parameters are defined as follows. +#' +#' \deqn{\epsilon_0 = f(0)} +#' +#' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)} +#' +#' This function is differentiable everywhere, monotonically +#' increasing, and \eqn{f(x) \approx x} if \eqn{x} is positive +#' and not too close to zero. By modifying the parameters, you +#' can control the distance between \eqn{f(x)} and the +#' horizontal axis at two 'places' -- \eqn{0} and \eqn{-\infty}. +#' [See issue #93](https://github.com/canmod/macpan2/issues/93). +#' for more information. +#' +#' For `clamp` the arguments specifically mean. +#' * `x` : A matrix with elements that should remain positive. +#' * `eps` : A small positive number, \eqn{\epsilon_0 = f(0)}, +#' giving the value of the function when the input is zero. +#' The default value is 1e-11 +#' * `limit` : A small positive number, +#' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)}, giving the +#' value of the function as the input goes to negative +#' infinity. The default is `limit = 1e-12`. This `limit` +#' should be chosen to be less than `eps` to ensure that +#' `clamp` is twice differentiable. +#' +#' #' ## Integer Sequences #' #' ### Functions @@ -673,48 +728,6 @@ #' ) #' ``` #' -#' ## Clamp -#' -#' Smoothly clamp the elements of a matrix so that they -#' do not get closer to 0 than a tolerance, `eps`, with -#' a default of 1e-12. This `clamp` function is the following -#' modification of the -#' [squareplus function](https://arxiv.org/abs/2112.11687). -#' -#' \deqn{f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}} -#' -#' Where the two parameters are defined as follows. -#' -#' \deqn{\epsilon_0 = f(0)} -#' -#' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)} -#' -#' This function is differentiable everywhere, monotonically -#' increasing, and \eqn{f(x) \approx x} if \eqn{x} is positive -#' and not too close to zero. By modifying the parameters, you -#' can control the distance between \eqn{f(x)} and the -#' horizontal axis at two 'places' -- \eqn{0} and \eqn{-\infty}. -#' [See issue #93](https://github.com/canmod/macpan2/issues/93). -#' for more information. -#' -#' ### Functions -#' -#' * `clamp(x, eps, limit)` -#' -#' ### Arguments -#' -#' * `x` : A matrix with elements that should remain positive. -#' * `eps` : A small positive number, \eqn{\epsilon_0 = f(0)}, -#' giving the value of the function when the input is zero. -#' The default value is 1e-11 -#' * `limit` : A small positive number, -#' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)}, giving the -#' value of the function as the input goes to negative -#' infinity. The default is `limit = 1e-12`. This `limit` -#' should be chosen to be less than `eps` to ensure that -#' `clamp` is twice differentiable. -#' -#' #' ## Probability Densities #' #' All probability densities have the same first two @@ -1081,6 +1094,7 @@ #' @aliases logit #' @aliases cumsum #' @aliases sparse_mat_mult +#' @aliases divide_safe #' @aliases assign #' @aliases unpack NULL diff --git a/R/enum.R b/R/enum.R index 878ce82ea..fc502f4c7 100644 --- a/R/enum.R +++ b/R/enum.R @@ -50,7 +50,7 @@ valid_func_sigs = c( , "fwrap: pgamma(q, shape, scale)" , "fwrap: mean(x)" , "fwrap: sd(x)" - , "fwrap: proportions(x)" + , "fwrap: proportions(x, limit, eps)" , "fwrap: last(x)" , "fwrap: check_finite(x)" , "fwrap: dbinom(observed, size, probability)" @@ -61,6 +61,7 @@ valid_func_sigs = c( , "fwrap: logit(x)" , "fwrap: cumsum(x)" , "fwrap: sparse_mat_mult(x, i, j, y, z)" + , "fwrap: divide_safe(x, eps, limit)" , "fwrap: assign(x, i, j, v)" , "fwrap: unpack(x, ...)" ) diff --git a/R/mp_tmb_model_spec.R b/R/mp_tmb_model_spec.R index 63edf4318..265608565 100644 --- a/R/mp_tmb_model_spec.R +++ b/R/mp_tmb_model_spec.R @@ -98,11 +98,22 @@ TMBModelSpec = function( ## an 'implied' integer vector for subsetting vectors in before, during, ## and after expressions by position name self$all_integers = function() { - ## TODO: make smarter so that only used integer vectors - ## are produced and maybe even check if an integer vector + ## TODO: maybe make smarter by checking if an integer vector ## is being used in the wrong numeric vector implied_integers = implied_position_vectors(self$default) - c(implied_integers, self$integers) + integers_we_need = intersect( + self$all_formula_vars() + , names(implied_integers) + ) + filtered_implied_integers = list() + nms = names(integers_we_need) + for (nm in unique(nms)) { + integers_nm = integers_we_need[nms == nm] + # if (sum(!duplicated(integers_nm)) != 1L) { + # } + filtered_implied_integers[[nm]] = integers_nm[[1L]] + } + c(filtered_implied_integers, self$integers) } self$empty_matrices = function() { diff --git a/man/engine_functions.Rd b/man/engine_functions.Rd index 8b9c3e912..76d0775ec 100644 --- a/man/engine_functions.Rd +++ b/man/engine_functions.Rd @@ -63,6 +63,7 @@ \alias{logit} \alias{cumsum} \alias{sparse_mat_mult} +\alias{divide_safe} \alias{assign} \alias{unpack} \title{Functions Available in the Simulation Engine} @@ -191,10 +192,12 @@ the results of applying the function to each element of \code{x}. } -\subsection{Proportions}{ +\subsection{Limiting Values}{ \subsection{Functions}{ \itemize{ \item \code{proportions(x, limit, eps)} +\item \code{limit(x, limit, eps)} +\item \code{clamp(x, eps, limit)} } } @@ -219,6 +222,67 @@ the results of applying the function to each element of \code{x}. }\if{html}{\out{}} } +\subsection{Details}{ + +The return value depends on the function. +\itemize{ +\item \code{proportions} : matrix of \code{x / sum(x)} or \code{rep(limit, length(x))} if +\code{sum(x) < eps}. +\item \code{limit} : +\item \code{clamp} : +} +} + +\subsection{Details}{ + +The divide_safe() function conducts elementwise division +of the first two arguments, and then allows two other +arguments: +\itemize{ +\item \code{limit} : Value to return if the \code{sqrt(denominator^2)} is +} +} + +\subsection{Details}{ + +The \code{clamp} function smoothly clamps the elements of a +matrix so that they +do not get closer to 0 than a tolerance, \code{eps}, with +a default of 1e-12. This \code{clamp} function is the following +modification of the +\href{https://arxiv.org/abs/2112.11687}{squareplus function}. + +\deqn{f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}} + +Where the two parameters are defined as follows. + +\deqn{\epsilon_0 = f(0)} + +\deqn{\epsilon_- = \lim_{x \to -\infty}f(x)} + +This function is differentiable everywhere, monotonically +increasing, and \eqn{f(x) \approx x} if \eqn{x} is positive +and not too close to zero. By modifying the parameters, you +can control the distance between \eqn{f(x)} and the +horizontal axis at two 'places' -- \eqn{0} and \eqn{-\infty}. +\href{https://github.com/canmod/macpan2/issues/93}{See issue #93}. +for more information. + +For \code{clamp} the arguments specifically mean. +\itemize{ +\item \code{x} : A matrix with elements that should remain positive. +\item \code{eps} : A small positive number, \eqn{\epsilon_0 = f(0)}, +giving the value of the function when the input is zero. +The default value is 1e-11 +\item \code{limit} : A small positive number, +\deqn{\epsilon_- = \lim_{x \to -\infty}f(x)}, giving the +value of the function as the input goes to negative +infinity. The default is \code{limit = 1e-12}. This \code{limit} +should be chosen to be less than \code{eps} to ensure that +\code{clamp} is twice differentiable. +} +} + } \subsection{Integer Sequences}{ @@ -838,52 +902,6 @@ calibrating) for time steps less than \eqn{m}. } -\subsection{Clamp}{ - -Smoothly clamp the elements of a matrix so that they -do not get closer to 0 than a tolerance, \code{eps}, with -a default of 1e-12. This \code{clamp} function is the following -modification of the -\href{https://arxiv.org/abs/2112.11687}{squareplus function}. - -\deqn{f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}} - -Where the two parameters are defined as follows. - -\deqn{\epsilon_0 = f(0)} - -\deqn{\epsilon_- = \lim_{x \to -\infty}f(x)} - -This function is differentiable everywhere, monotonically -increasing, and \eqn{f(x) \approx x} if \eqn{x} is positive -and not too close to zero. By modifying the parameters, you -can control the distance between \eqn{f(x)} and the -horizontal axis at two 'places' -- \eqn{0} and \eqn{-\infty}. -\href{https://github.com/canmod/macpan2/issues/93}{See issue #93}. -for more information. -\subsection{Functions}{ -\itemize{ -\item \code{clamp(x, eps, limit)} -} -} - -\subsection{Arguments}{ -\itemize{ -\item \code{x} : A matrix with elements that should remain positive. -\item \code{eps} : A small positive number, \eqn{\epsilon_0 = f(0)}, -giving the value of the function when the input is zero. -The default value is 1e-11 -\item \code{limit} : A small positive number, -\deqn{\epsilon_- = \lim_{x \to -\infty}f(x)}, giving the -value of the function as the input goes to negative -infinity. The default is \code{limit = 1e-12}. This \code{limit} -should be chosen to be less than \code{eps} to ensure that -\code{clamp} is twice differentiable. -} -} - -} - \subsection{Probability Densities}{ All probability densities have the same first two diff --git a/misc/dev/dev.cpp b/misc/dev/dev.cpp index a4068c73f..4bda09de8 100644 --- a/misc/dev/dev.cpp +++ b/misc/dev/dev.cpp @@ -126,7 +126,7 @@ enum macpan2_func , MP2_PGAMMA = 48 // fwrap: pgamma(q, shape, scale) , MP2_MEAN = 49 // fwrap: mean(x) , MP2_SD = 50 // fwrap: sd(x) - , MP2_PROPORTIONS = 51 // fwrap: proportions(x) + , MP2_PROPORTIONS = 51 // fwrap: proportions(x, limit, eps) , MP2_LAST = 52 // fwrap: last(x) , MP2_CHECK_FINITE = 53 // fwrap: check_finite(x) , MP2_BINOM_DENSITY = 54 // fwrap: dbinom(observed, size, probability) @@ -137,8 +137,9 @@ enum macpan2_func , MP2_LOGIT = 59 // fwrap: logit(x) , MP2_CUMSUM = 60 // fwrap: cumsum(x) , MP2_SPARSE_MAT_MULT = 61 // fwrap: sparse_mat_mult(x, i, j, y, z) - , MP2_ASSIGN = 62 // fwrap: assign(x, i, j, v) - , MP2_UNPACK = 63 // fwrap: unpack(x, ...) + , MP2_DIVIDE_SAFE = 62 // fwrap: divide_safe(x, eps, limit) + , MP2_ASSIGN = 63 // fwrap: assign(x, i, j, v) + , MP2_UNPACK = 64 // fwrap: unpack(x, ...) }; enum macpan2_meth @@ -1012,10 +1013,14 @@ class ExprEvaluator bool is_finite_mat; vector u; // FIXME: why not std::vector here?? matrix Y, X, A; - std::vector timeIndex; // for rbind_time and rbind_lag + std::vector timeIndex; // rbind_time and rbind_lag int doing_lag = 0; - Type sum, eps, limit, var, by, left_over, remaining_prop, p0; // intermediate scalars - Type delta_t; // for reulermultinom + Type zero = 0; + Type var; // variance + Type by; // seq + Type sum, eps, limit; // proportions, divide_safe, clamp + Type den, abs_den; // divide_safe + Type p0, remaining_prop, left_over, delta_t; // reulermultinom int rows, cols, lag, rowIndex, colIndex, matIndex, cp, off, size, times; unsigned int grpIndex; int size_in, size_out; @@ -1385,11 +1390,15 @@ class ExprEvaluator case MP2_LOGIT: return (-(1 / args[0].array() - 1).log()).matrix(); - // #' ## Proportions + // #' ## Limiting Values // #' + // #' + // #' // #' ### Functions // #' // #' * `proportions(x, limit, eps)` + // #' * `limit(x, limit, eps)` + // #' * `clamp(x, eps, limit)` // #' // #' ### Arguments // #' @@ -1401,13 +1410,20 @@ class ExprEvaluator // #' // #' * matrix of `x / sum(x)` or `rep(limit, length(x))` if // #' `sum(x) < eps`. - // #' + // #' // #' ### Examples // #' // #' ``` // #' engine_eval(~ proportions(y, 0.5, 1e-8), y = c(2, 0.5)) // #' ``` // #' + // #' ### Details + // #' + // #' The return value depends on the function. + // #' * `proportions` : matrix of `x / sum(x)` or `rep(limit, length(x))` if + // #' `sum(x) < eps`. + // #' * `limit` : + // #' * `clamp` : case MP2_PROPORTIONS: m = args.get_as_mat(0); m1 = matrix::Zero(1, 1); @@ -1423,7 +1439,104 @@ class ExprEvaluator } } return m2; - + + // #' ### Details + // #' + // #' The divide_safe() function conducts elementwise division + // #' of the first two arguments, and then allows two other + // #' arguments: + // #' * `limit` : Value to return if the `sqrt(denominator^2)` is + // + + case MP2_DIVIDE_SAFE: + args = args.recycle_for_bin_op(); + m = args.get_as_mat(0); + m1 = args.get_as_mat(1); + m2 = m.cwiseQuotient(m1); + limit = args.get_as_mat(2).coeff(0, 0); + eps = args.get_as_mat(3).coeff(0, 0); + for (int i = 0; i < m.rows(); i++) { + for (int j = 0; j < m.cols(); j++) { + den = m1.coeff(i, j); + abs_den = CppAD::CondExpGt(den, zero, den, -den); + m2.coeffRef(i, j) = CppAD::CondExpLt( + abs_den + , eps + , limit + , m2.coeffRef(i, j) + ); + } + } + return m2; + + // #' ### Details + // #' + // #' The `clamp` function smoothly clamps the elements of a + // #' matrix so that they + // #' do not get closer to 0 than a tolerance, `eps`, with + // #' a default of 1e-12. This `clamp` function is the following + // #' modification of the + // #' [squareplus function](https://arxiv.org/abs/2112.11687). + // #' + // #' \deqn{f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}} + // #' + // #' Where the two parameters are defined as follows. + // #' + // #' \deqn{\epsilon_0 = f(0)} + // #' + // #' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)} + // #' + // #' This function is differentiable everywhere, monotonically + // #' increasing, and \eqn{f(x) \approx x} if \eqn{x} is positive + // #' and not too close to zero. By modifying the parameters, you + // #' can control the distance between \eqn{f(x)} and the + // #' horizontal axis at two 'places' -- \eqn{0} and \eqn{-\infty}. + // #' [See issue #93](https://github.com/canmod/macpan2/issues/93). + // #' for more information. + // #' + // #' For `clamp` the arguments specifically mean. + // #' * `x` : A matrix with elements that should remain positive. + // #' * `eps` : A small positive number, \eqn{\epsilon_0 = f(0)}, + // #' giving the value of the function when the input is zero. + // #' The default value is 1e-11 + // #' * `limit` : A small positive number, + // #' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)}, giving the + // #' value of the function as the input goes to negative + // #' infinity. The default is `limit = 1e-12`. This `limit` + // #' should be chosen to be less than `eps` to ensure that + // #' `clamp` is twice differentiable. + // #' + // #' + case MP2_CLAMP: + eps = 1e-11; // default + if (n > 1) + eps = args[1].coeff(0, 0); + if (n == 3) { + limit = args[2].coeff(0, 0); + } else { + limit = 1e-12; // default + } + X = args[0]; + rows = X.rows(); + cols = X.cols(); + m = matrix::Zero(rows, cols); + + // https://github.com/canmod/macpan2/issues/93 + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + m.coeffRef(i, j) = limit + ( + ( + X.coeff(i, j) - limit + + sqrt( + pow(X.coeff(i, j) - limit, 2.0) + + pow(2.0 * eps - limit, 2.0) - + pow(limit, 2.0) + ) + ) / 2.0 + ); + } + } + return m; // #' ## Integer Sequences @@ -2612,79 +2725,6 @@ class ExprEvaluator MP2_ERR(MP2_CONVOLUTION, "Either empty or non-column vector used as kernel in convolution", MP2_CONVOLUTION); } - // #' ## Clamp - // #' - // #' Smoothly clamp the elements of a matrix so that they - // #' do not get closer to 0 than a tolerance, `eps`, with - // #' a default of 1e-12. This `clamp` function is the following - // #' modification of the - // #' [squareplus function](https://arxiv.org/abs/2112.11687). - // #' - // #' \deqn{f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}} - // #' - // #' Where the two parameters are defined as follows. - // #' - // #' \deqn{\epsilon_0 = f(0)} - // #' - // #' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)} - // #' - // #' This function is differentiable everywhere, monotonically - // #' increasing, and \eqn{f(x) \approx x} if \eqn{x} is positive - // #' and not too close to zero. By modifying the parameters, you - // #' can control the distance between \eqn{f(x)} and the - // #' horizontal axis at two 'places' -- \eqn{0} and \eqn{-\infty}. - // #' [See issue #93](https://github.com/canmod/macpan2/issues/93). - // #' for more information. - // #' - // #' ### Functions - // #' - // #' * `clamp(x, eps, limit)` - // #' - // #' ### Arguments - // #' - // #' * `x` : A matrix with elements that should remain positive. - // #' * `eps` : A small positive number, \eqn{\epsilon_0 = f(0)}, - // #' giving the value of the function when the input is zero. - // #' The default value is 1e-11 - // #' * `limit` : A small positive number, - // #' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)}, giving the - // #' value of the function as the input goes to negative - // #' infinity. The default is `limit = 1e-12`. This `limit` - // #' should be chosen to be less than `eps` to ensure that - // #' `clamp` is twice differentiable. - // #' - // #' - case MP2_CLAMP: - eps = 1e-11; // default - if (n > 1) - eps = args[1].coeff(0, 0); - if (n == 3) { - limit = args[2].coeff(0, 0); - } else { - limit = 1e-12; // default - } - X = args[0]; - rows = X.rows(); - cols = X.cols(); - m = matrix::Zero(rows, cols); - - // https://github.com/canmod/macpan2/issues/93 - for (int i = 0; i < rows; i++) { - for (int j = 0; j < cols; j++) { - m.coeffRef(i, j) = limit + ( - ( - X.coeff(i, j) - limit + - sqrt( - pow(X.coeff(i, j) - limit, 2.0) + - pow(2.0 * eps - limit, 2.0) - - pow(limit, 2.0) - ) - ) / 2.0 - ); - } - } - return m; - // #' ## Probability Densities // #' // #' All probability densities have the same first two diff --git a/misc/old-vignettes/general_compartmental_model.Rmd b/misc/old-vignettes/general_compartmental_model.Rmd index 23215bb73..142949248 100644 --- a/misc/old-vignettes/general_compartmental_model.Rmd +++ b/misc/old-vignettes/general_compartmental_model.Rmd @@ -304,7 +304,7 @@ magnitude of the flow. The user may choose from a number of flow types for each We present an example model that is just complex enough to illustrate the general model. -![](../misc/diagrams/model_matrices.png) +![](../diagrams/model_matrices.png) The state vector is given as follows. diff --git a/src/macpan2.cpp b/src/macpan2.cpp index e6654a9b3..0857e26e1 100644 --- a/src/macpan2.cpp +++ b/src/macpan2.cpp @@ -127,7 +127,7 @@ enum macpan2_func , MP2_PGAMMA = 48 // fwrap: pgamma(q, shape, scale) , MP2_MEAN = 49 // fwrap: mean(x) , MP2_SD = 50 // fwrap: sd(x) - , MP2_PROPORTIONS = 51 // fwrap: proportions(x) + , MP2_PROPORTIONS = 51 // fwrap: proportions(x, limit, eps) , MP2_LAST = 52 // fwrap: last(x) , MP2_CHECK_FINITE = 53 // fwrap: check_finite(x) , MP2_BINOM_DENSITY = 54 // fwrap: dbinom(observed, size, probability) @@ -138,8 +138,9 @@ enum macpan2_func , MP2_LOGIT = 59 // fwrap: logit(x) , MP2_CUMSUM = 60 // fwrap: cumsum(x) , MP2_SPARSE_MAT_MULT = 61 // fwrap: sparse_mat_mult(x, i, j, y, z) - , MP2_ASSIGN = 62 // fwrap: assign(x, i, j, v) - , MP2_UNPACK = 63 // fwrap: unpack(x, ...) + , MP2_DIVIDE_SAFE = 62 // fwrap: divide_safe(x, eps, limit) + , MP2_ASSIGN = 63 // fwrap: assign(x, i, j, v) + , MP2_UNPACK = 64 // fwrap: unpack(x, ...) }; enum macpan2_meth @@ -1013,10 +1014,14 @@ class ExprEvaluator bool is_finite_mat; vector u; // FIXME: why not std::vector here?? matrix Y, X, A; - std::vector timeIndex; // for rbind_time and rbind_lag + std::vector timeIndex; // rbind_time and rbind_lag int doing_lag = 0; - Type sum, eps, limit, var, by, left_over, remaining_prop, p0; // intermediate scalars - Type delta_t; // for reulermultinom + Type zero = 0; + Type var; // variance + Type by; // seq + Type sum, eps, limit; // proportions, divide_safe, clamp + Type den, abs_den; // divide_safe + Type p0, remaining_prop, left_over, delta_t; // reulermultinom int rows, cols, lag, rowIndex, colIndex, matIndex, cp, off, size, times; unsigned int grpIndex; int size_in, size_out; @@ -1386,11 +1391,15 @@ class ExprEvaluator case MP2_LOGIT: return (-(1 / args[0].array() - 1).log()).matrix(); - // #' ## Proportions + // #' ## Limiting Values // #' + // #' + // #' // #' ### Functions // #' // #' * `proportions(x, limit, eps)` + // #' * `limit(x, limit, eps)` + // #' * `clamp(x, eps, limit)` // #' // #' ### Arguments // #' @@ -1402,13 +1411,20 @@ class ExprEvaluator // #' // #' * matrix of `x / sum(x)` or `rep(limit, length(x))` if // #' `sum(x) < eps`. - // #' + // #' // #' ### Examples // #' // #' ``` // #' engine_eval(~ proportions(y, 0.5, 1e-8), y = c(2, 0.5)) // #' ``` // #' + // #' ### Details + // #' + // #' The return value depends on the function. + // #' * `proportions` : matrix of `x / sum(x)` or `rep(limit, length(x))` if + // #' `sum(x) < eps`. + // #' * `limit` : + // #' * `clamp` : case MP2_PROPORTIONS: m = args.get_as_mat(0); m1 = matrix::Zero(1, 1); @@ -1424,7 +1440,104 @@ class ExprEvaluator } } return m2; - + + // #' ### Details + // #' + // #' The divide_safe() function conducts elementwise division + // #' of the first two arguments, and then allows two other + // #' arguments: + // #' * `limit` : Value to return if the `sqrt(denominator^2)` is + // + + case MP2_DIVIDE_SAFE: + args = args.recycle_for_bin_op(); + m = args.get_as_mat(0); + m1 = args.get_as_mat(1); + m2 = m.cwiseQuotient(m1); + limit = args.get_as_mat(2).coeff(0, 0); + eps = args.get_as_mat(3).coeff(0, 0); + for (int i = 0; i < m.rows(); i++) { + for (int j = 0; j < m.cols(); j++) { + den = m1.coeff(i, j); + abs_den = CppAD::CondExpGt(den, zero, den, -den); + m2.coeffRef(i, j) = CppAD::CondExpLt( + abs_den + , eps + , limit + , m2.coeffRef(i, j) + ); + } + } + return m2; + + // #' ### Details + // #' + // #' The `clamp` function smoothly clamps the elements of a + // #' matrix so that they + // #' do not get closer to 0 than a tolerance, `eps`, with + // #' a default of 1e-12. This `clamp` function is the following + // #' modification of the + // #' [squareplus function](https://arxiv.org/abs/2112.11687). + // #' + // #' \deqn{f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}} + // #' + // #' Where the two parameters are defined as follows. + // #' + // #' \deqn{\epsilon_0 = f(0)} + // #' + // #' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)} + // #' + // #' This function is differentiable everywhere, monotonically + // #' increasing, and \eqn{f(x) \approx x} if \eqn{x} is positive + // #' and not too close to zero. By modifying the parameters, you + // #' can control the distance between \eqn{f(x)} and the + // #' horizontal axis at two 'places' -- \eqn{0} and \eqn{-\infty}. + // #' [See issue #93](https://github.com/canmod/macpan2/issues/93). + // #' for more information. + // #' + // #' For `clamp` the arguments specifically mean. + // #' * `x` : A matrix with elements that should remain positive. + // #' * `eps` : A small positive number, \eqn{\epsilon_0 = f(0)}, + // #' giving the value of the function when the input is zero. + // #' The default value is 1e-11 + // #' * `limit` : A small positive number, + // #' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)}, giving the + // #' value of the function as the input goes to negative + // #' infinity. The default is `limit = 1e-12`. This `limit` + // #' should be chosen to be less than `eps` to ensure that + // #' `clamp` is twice differentiable. + // #' + // #' + case MP2_CLAMP: + eps = 1e-11; // default + if (n > 1) + eps = args[1].coeff(0, 0); + if (n == 3) { + limit = args[2].coeff(0, 0); + } else { + limit = 1e-12; // default + } + X = args[0]; + rows = X.rows(); + cols = X.cols(); + m = matrix::Zero(rows, cols); + + // https://github.com/canmod/macpan2/issues/93 + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + m.coeffRef(i, j) = limit + ( + ( + X.coeff(i, j) - limit + + sqrt( + pow(X.coeff(i, j) - limit, 2.0) + + pow(2.0 * eps - limit, 2.0) - + pow(limit, 2.0) + ) + ) / 2.0 + ); + } + } + return m; // #' ## Integer Sequences @@ -2613,79 +2726,6 @@ class ExprEvaluator MP2_ERR(MP2_CONVOLUTION, "Either empty or non-column vector used as kernel in convolution", MP2_CONVOLUTION); } - // #' ## Clamp - // #' - // #' Smoothly clamp the elements of a matrix so that they - // #' do not get closer to 0 than a tolerance, `eps`, with - // #' a default of 1e-12. This `clamp` function is the following - // #' modification of the - // #' [squareplus function](https://arxiv.org/abs/2112.11687). - // #' - // #' \deqn{f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}} - // #' - // #' Where the two parameters are defined as follows. - // #' - // #' \deqn{\epsilon_0 = f(0)} - // #' - // #' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)} - // #' - // #' This function is differentiable everywhere, monotonically - // #' increasing, and \eqn{f(x) \approx x} if \eqn{x} is positive - // #' and not too close to zero. By modifying the parameters, you - // #' can control the distance between \eqn{f(x)} and the - // #' horizontal axis at two 'places' -- \eqn{0} and \eqn{-\infty}. - // #' [See issue #93](https://github.com/canmod/macpan2/issues/93). - // #' for more information. - // #' - // #' ### Functions - // #' - // #' * `clamp(x, eps, limit)` - // #' - // #' ### Arguments - // #' - // #' * `x` : A matrix with elements that should remain positive. - // #' * `eps` : A small positive number, \eqn{\epsilon_0 = f(0)}, - // #' giving the value of the function when the input is zero. - // #' The default value is 1e-11 - // #' * `limit` : A small positive number, - // #' \deqn{\epsilon_- = \lim_{x \to -\infty}f(x)}, giving the - // #' value of the function as the input goes to negative - // #' infinity. The default is `limit = 1e-12`. This `limit` - // #' should be chosen to be less than `eps` to ensure that - // #' `clamp` is twice differentiable. - // #' - // #' - case MP2_CLAMP: - eps = 1e-11; // default - if (n > 1) - eps = args[1].coeff(0, 0); - if (n == 3) { - limit = args[2].coeff(0, 0); - } else { - limit = 1e-12; // default - } - X = args[0]; - rows = X.rows(); - cols = X.cols(); - m = matrix::Zero(rows, cols); - - // https://github.com/canmod/macpan2/issues/93 - for (int i = 0; i < rows; i++) { - for (int j = 0; j < cols; j++) { - m.coeffRef(i, j) = limit + ( - ( - X.coeff(i, j) - limit + - sqrt( - pow(X.coeff(i, j) - limit, 2.0) + - pow(2.0 * eps - limit, 2.0) - - pow(limit, 2.0) - ) - ) / 2.0 - ); - } - } - return m; - // #' ## Probability Densities // #' // #' All probability densities have the same first two diff --git a/tests/testthat/test-clamp.R b/tests/testthat/test-clamp.R index 1b90b9f5e..81bab55ed 100644 --- a/tests/testthat/test-clamp.R +++ b/tests/testthat/test-clamp.R @@ -4,10 +4,10 @@ test_that("clamping function is modified squareplus", { clamp_base_r = function(x, eps = 1e-12, limit = eps) { limit + ( ( - x - limit + + x - limit + sqrt( - (x - limit)^2 + - (2 * eps - limit)^2 - + (x - limit)^2 + + (2 * eps - limit)^2 - limit^2 ) ) / 2 @@ -16,3 +16,5 @@ test_that("clamping function is modified squareplus", { clamp_macpan = function(x) engine_eval(~clamp(x), x = x) |> c() expect_equal(clamp_base_r(x), clamp_macpan(x)) }) + +engine_eval(~divide_safe(1e-9, 1e-9, 0, 1e-12)) diff --git a/vignettes/kronecker.Rmd b/vignettes/kronecker.Rmd new file mode 100644 index 000000000..f30765eb7 --- /dev/null +++ b/vignettes/kronecker.Rmd @@ -0,0 +1,419 @@ +--- +title: "Structured SI Model with Kronecker Products in macpan2" +output: + rmarkdown::html_vignette: + toc: true +vignette: > + %\VignetteIndexEntry{Structured SI Model with Kronecker Products in macpan2} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} +editor_options: + chunk_output_type: console +--- + +[![status](https://img.shields.io/badge/status-mature%20draft-yellow)](https://canmod.github.io/macpan2/articles/vignette-status#mature-draft) + +```{r settings, echo = FALSE, message=FALSE, warning=FALSE, error=FALSE} +knitr::opts_chunk$set(echo = TRUE) +library(macpan2); library(ggplot2); library(tidyr); library(dplyr); library(stringr) +theme_bw = function() ggplot2::theme_bw(base_size = 14) +``` + +## Introduction + +This document illustrates how to implement structured compartmental models in `macpan2` using Kronecker products to define state vectors, parameter matrices, and flows over stratified dimensions. As an example, we stratify an SI model along the following dimensions: age group, vaccination status, and symptom severity. + +```{r model-dimensions} +model_dimensions = c("age", "vax", "symp") +age = c("young", "old") +vax = c("unvax", "vax") +symp = c("mild", "severe") +``` + +## Utilities and Packages + +We load the necessary packages and define some utility functions that simplify the construction of named vectors and matrices over these stratification dimensions. + +```{r packages} +library(macpan2); library(ggplot2); library(tidyr); library(dplyr); library(stringr) +``` + +We overwrite the base-R Kronecker product operator, `%x%`, so that it preserves row and column names in a consistent manner. We will use this operator many times when defining model structure below. + +```{r operator-kron} +`%x%` <- mp_kronecker_operator(`*`) +``` + +The `ones_vec()` utility creates a named vector of ones over the cross-product of the input character vectors (strata). It uses Kronecker products to ensure the names follow the same ordering convention as Kronecker-structured vectors and matrices elsewhere in the model. + +```{r function-ones-vec} +ones_vec = function(...) { + nms = list(...) + vecs = lapply(nms, \(x) setNames(rep(1, length(x)), x)) + Reduce(`%x%`, vecs) +} +``` + +The `diagonal()` utility is much like the base-R `diag()` function, but takes its row and column names from the names of the input vector. +```{r diagonal-func} +diagonal = function(x) { + if (is.matrix(x)) + stop("can only produce a square diagonal matrix for a vector") + if (length(x) == 1) return(x) + nms = names(x) + x = diag(x) + dimnames(x) = list(nms, nms) + return(x) +} +``` + +The `identity()` function takes character vectors and returns an identity matrix over their Cartesian product. It uses `ones_vec()` and `diagonal()` to define the dimension names and place ones on the diagonal. + +```{r function-identity} +identity = function(...) ones_vec(...) |> diagonal() +``` + +The `ones_col()` (`ones_row()`) function takes character vectors and returns a column (row) vector of ones, with row (column) names based on the Cartesian product of the input strata. + +```{r function-ones-col} +ones_col = function(...) ones_vec(...) |> t() |> t() +ones_row = function(...) ones_col(...) |> t() +``` + + +Here are a few examples. + +```{r example-1-mats} +identity(age, vax) +ones_col(age, vax, symp) +ones_row(vax) +``` + + +The `ones_col()`, `ones_row()`, and `identity()` functions make it easier to build vectors and matrices that align with stratified model dimensions, enabling clear and consistent operations like summing over strata or selecting specific blocks. + + +## Initial State Vectors + +We now initialize the state vectors `S` (susceptible) and `I` (infectious) over the stratification dimensions. + +```{r initial-states} +S = c(young = 75, old = 25) %x% c(unvax = 1, vax = 0) +I = c(young = 1, old = 0) %x% c(unvax = 1, vax = 0) %x% c(mild = 1, severe = 0) +print(S) +print(I) +``` + +## Model Parameters + +We define default values for the following parameters: + +* Baseline transmission +* Vaccination rates +* Vaccine efficacy +* Contact matrix +* Symptom severity probabilities +* Multiplicative effects of stratification on susceptibility and infectivity. + +```{r parameters} +beta = 0.4 +contact_matrix = rbind( + young = c(young = 0.8, old = 0.2), + old = c(young = 0.2, old = 0.8) +) +symp_probs = rbind(mild = 0.6, severe = 0.4) +susceptibility_age = rbind(young = 0.9, old = 1) +infectivity_age = cbind(young = 1, old = 0.9) +infectivity_vax = cbind(unvax = 1, vax = 0.8) +infectivity_symp = cbind(mild = 0.95, severe = 1) +vax_rates = rbind(young = 0.1, old = 0.4) +VE = rbind(unvax = 0, vax = 0.8) +``` + +Note that instead of `c` for constructing vectors, we use `rbind` (for column vectors) and `cbind` (for row vectors). It is good practice to be explicit about the orientation of vectors when Kronecker products are involved. + +## Derived State Vectors + +To define flows and compute quantities like force of infection or vaccination rates, we often need to aggregate or restrict the original state vectors by strata. Derived state vectors let us isolate relevant subpopulations or repeat totals to match dimensional structure in downstream operations. + +```{r derived-states} +M_unvax = identity(age) %x% cbind(unvax = 1, vax = 0) +M_S = identity(age) %x% ones_row(vax) +M_I = identity(age) %x% ones_row(vax, symp) +M_N = ones_col(vax, symp) +S_unvax = M_unvax %*% S +N_age = (M_S %*% S) + (M_I %*% I) +N = N_age %x% M_N +print(S_unvax) +print(N_age) +print(N) +``` + +Note that `N` repeats the total population size within each age group to match the structure of the force of infection. `S_unvax` isolates the unvaccinated susceptible compartments, which are the only ones eligible for vaccination flows. + + +## Per-Capita Transmission Matrix + +We construct the full per-capita transmission matrix by combining component matrices using Kronecker products. + +```{r per-capita-transmission-matrix} +B_age = beta * contact_matrix * (susceptibility_age %x% infectivity_age) +B_vax = (1 - VE) %x% infectivity_vax +B_symp = infectivity_symp +B = B_age %x% B_vax %x% B_symp +print(B_age) +print(B_vax) +print(B_symp) +print(B) +``` + +The per-capita transmission matrix maps infectious compartments to susceptible compartments and determines the force of infection associated with each susceptible stratum. The row names give the susceptible strata, and the column names the infectious strata. + +## Per-Capita Force of Infection + +The force of infection is computed by multiplying the per-capita transmission matrix by the normalized infectious population. + +```{r force-of-infection} +foi = B %*% (I / N) +print(foi) +``` + +## Absolute Flows + +We compute total flows due to infection and vaccination by multiplying per-capita rates by the relevant state vectors. + +```{r absolute-flows} +infection = foi * S +vaccination = vax_rates * S_unvax +print(infection) +print(vaccination) +``` + +## Allocation Matrices + +In structured models, flows often act on families of compartments rather than individual compartments. This can obscure which specific sub-compartments are gaining or losing individuals. Allocation matrices resolve this by mapping flow vectors into change vectors that align with the structure of state vectors. + +They allow us to: + +* Distribute outflows from a compartment to multiple destination compartments (e.g., infections split by symptom severity). +* Apply signed changes within a state vector (e.g., moving individuals from unvaccinated to vaccinated compartments within the susceptible state vector). + +This ensures that each flow is properly aligned with the compartments it affects. + +```{r allocations} +A_infection = identity(age, vax) %x% symp_probs +A_vaccination = identity(age) %x% rbind(unvax = -1, vax = 1) +infection_S = infection +infection_I = A_infection %*% infection +vaccination_S = A_vaccination %*% vaccination +print(infection_S) +print(infection_I) +print(vaccination_S) +``` + +## Euler Step Update + +Using the change vectors produced using products of allocation matrices and state vectors, we update the original state vectors with a forward Euler step. + +```{r euler-step} +S_new = S - infection_S + vaccination_S +I_new = I + infection_I +print(S_new) +print(I_new) +``` + +## Model Specification + +We formalize the model specification with `mp_tmb_model_spec`, defining the derived quantities, flows, and update rules. + +```{r model-spec} +spec = mp_tmb_model_spec( + before = list( + B_age ~ beta * contact_matrix * (susceptibility_age %x% infectivity_age) + , B_vax ~ (1 - VE) %x% infectivity_vax + , B_symp ~ infectivity_symp + , B ~ B_age %x% B_vax %x% B_symp + + ), + during = list( + + # derived state vectors + N_age ~ (M_S %*% S) + (M_I %*% I) + , N ~ N_age %x% M_N + , S_unvax ~ M_unvax %*% S + + # state-dependent per-capita rates + , foi ~ B %*% (I / N) + + # absolute flow vectors + , infection ~ foi * S + , vaccination ~ vax_rates * S_unvax + + # absolute change vectors + , infection_S ~ -infection + , infection_I ~ A_infection %*% infection + , vaccination_S ~ A_vaccination %*% vaccination + + # state update (euler step) + , S ~ S + infection_S + vaccination_S + , I ~ I + infection_I + ), + + # initial state vectors + inits = nlist(S, I), + + default = nlist( + # parameter scalars + beta + + # parameter vectors + , vax_rates, VE + , susceptibility_age + , infectivity_age, infectivity_vax, infectivity_symp + + # parameter matrices + , contact_matrix + + # state vector transformation matrices + , M_S, M_I, M_N, M_unvax + + # allocation matrices + , A_infection, A_vaccination + + # absolute flow and change vectors + # - included here so that row names can be used + # in simulated trajectories + # - if removed from this list, the simulations + # will use generic integer IDs for each + # stratum + , infection, vaccination + , infection_S, infection_I, vaccination_S + ) +) +``` + +## Simulation and Plots + +We simulate the model over time and visualize the resulting trajectories of the state vectors and flow rates by stratification. The first step is to generate the simulations. + +```{r simulations} + +## state vectors +state = c("S", "I") + +## absolute flow rates +flow = c("infection", "vaccination") + +# absolute change vectors +change = c("infection_S", "infection_I", "vaccination_S") + +traj = (spec + |> mp_simulator(time_steps = 100, outputs = c(state, flow, change)) + |> mp_trajectory() + + ## separate row names into + ## names for each stratum + |> separate_wider_delim(row + , delim = "." + , names = model_dimensions + , too_few = "align_start" + ) + + ## make stratum names easier to read + |> mutate(vax = case_when( + vax == "unvax" ~ "Unvaccinated" + , vax == "vax" ~ "Vaccinated" + )) + |> mutate(age = str_to_title(age)) +) + +print(traj) +``` + +Here is a plot of the state vector trajectories. + +```{r simulation-plot, fig.height=7, fig.width=6} +(traj + |> filter(matrix %in% state) + + ## make names easier to read + |> mutate(state = case_when( + matrix == "S" ~ "Susceptible" + , matrix == "I" ~ "Infectious" + )) + |> mutate(state = if_else( + state == "Infectious" + , sprintf("%s\n(%s)", state, symp) + , state + )) + + ## plot + |> ggplot() + + aes(time, value) + + geom_line() + + facet_grid(state ~ age + vax, space = 'free', scales = 'free') + + scale_x_continuous(breaks = c(0, 40, 80)) + + theme_bw() +) +``` + +Structured SI models have the potential to become more interesting relative to the unstructured version. For a simple example in the above plot, two susceptible compartments display intermediate peaks that reflect a balance between vaccination and infection processes. In contrast, the susceptible compartment in the unstructured SI model can only decrease over time. You could produce this plot for different values of [the parameters](#model-parameters) (e.g., vaccine efficacy) to explore how these differences affect dynamics. + +We can make a similar plot of the absolute rates of flow (per time step) among compartments. These flow rates are stored in the `infection` and `vaccination` vectors. + +```{r simulation-plot-incidence, fig.height=7, fig.width=5} +(traj + |> filter(matrix %in% flow) + + |> mutate(matrix = str_to_title(matrix)) + |> mutate(flow = if_else( + matrix == "Infection" + , sprintf("%s\n(%s)", matrix, vax) + , matrix + )) + + |> ggplot() + + aes(time, value) + + geom_line() + + facet_grid(flow ~ age, scales = 'free') + + scale_x_continuous(breaks = c(0, 40, 80)) + + theme_bw() +) +``` + +In structured models, absolute flow rates translate into process-specific rates of change of state variables in a more complex way than in unstructured models. One example of this added complexity is that flows alone do not determine how individuals are distributed across symptom status compartments. To address this, allocation matrices -- `A_infection` and `A_vaccination` -- are used to map each process-specific flow onto the appropriate rate of change in the state variables. These matrices define how the flows are distributed across strata, such as symptom status. The plot below illustrates the resulting changes in the state vector. + +```{r simulation-plot-change, fig.height=7, fig.width=5} +(traj + |> filter(matrix %in% change) + + |> separate_wider_delim(matrix + , delim = "_" + , names = c("flow", "state") + ) + |> mutate(state = case_when( + state == "S" ~ "Susceptible" + , state == "I" ~ "Infectious" + )) + |> mutate(state = if_else( + state == "Infectious" + , sprintf("%s\n(%s)", state, symp) + , state + )) + |> mutate(flow = sprintf("Change due to\n%s", flow)) + + |> ggplot() + + aes(time, value, colour = vax) + + geom_line() + + facet_grid(flow + state ~ age, scales = 'free') + + scale_x_continuous(breaks = c(0, 40, 80)) + + theme_bw() + + guides(colour = guide_legend(position = "bottom", title = "")) +) +``` + +This plot reveals how susceptible individuals are allocated across symptom status compartments. By using color to represent vaccination status—instead of faceting into separate panels—it simultaneously highlights differences in rates of change across vaccination groups. + +## Limitations + +One cannot use [alternate state update methods](https://canmod.github.io/macpan2/reference/state_updates) with structured models, although this is [on the roadmap](https://github.com/canmod/macpan2/issues/288). diff --git a/vignettes/matrix_multiplication.Rmd b/vignettes/matrix_multiplication.Rmd index c0e11ba8a..1765534cc 100644 --- a/vignettes/matrix_multiplication.Rmd +++ b/vignettes/matrix_multiplication.Rmd @@ -1,10 +1,10 @@ --- -title: "Multiplying a Sparse Matrix with a Dense Matrix" +title: "Matrix Operations in macpan2 Models" output: rmarkdown::html_vignette: toc: true vignette: > - %\VignetteIndexEntry{Multiplying a Sparse Matrix with a Dense Matrix} + %\VignetteIndexEntry{Matrix Operations in macpan2 Models} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: @@ -13,82 +13,257 @@ editor_options: [![status](https://img.shields.io/badge/status-mature%20draft-yellow)](https://canmod.github.io/macpan2/articles/vignette-status#mature-draft) -```{r settings, echo = FALSE, eval = TRUE, message=FALSE, warning=FALSE, error=FALSE} +```{r settings, echo = FALSE, message=FALSE, warning=FALSE, error=FALSE} library(macpan2) -library(splines) library(dplyr) library(ggplot2) +library(splines) theme_bw = function() ggplot2::theme_bw(base_size = 18) ``` ## Overview -When working with models in `macpan2`, it is often necessary to perform matrix multiplications where the left-hand matrix is sparse—meaning most of its entries are zero. Passing large sparse matrices as dense matrices with zeros into `macpan2` models can lead to unnecessary memory use and performance slowdowns. +Matrix operations are central to `macpan2` models. The engine supports four ways to multiply two matrices together: -To address this, `macpan2` provides the engine function `sparse_mat_mult()`, which multiplies a sparse matrix (given in compressed form) by a dense matrix. Instead of passing the full sparse matrix, you pass: +- Standard matrix multiplication (`%*%`) +- Kronecker products (`%x%`) +- Elementwise multiplication with recycling (`*`) +- Left sparse matrix multiplication (`sparse_mat_mult()`) -- `x` — the non-zero values of the sparse matrix, as a numeric vector -- `i` — the zero-based row indices of each value in `x` -- `j` — the zero-based column indices of each value in `x` -- `y` — the dense matrix on the right-hand side of the multiplication -- `z` — a matrix with the correct dimensions, which is passed both as an input and as the destination for the result +This vignette introduces each operation and illustrates typical use cases. -This operation is equivalent to: -```r -z ~ LeftSparseMatrix %*% y +## Standard Matrix Multiplication (`%*%`) + +Dense matrix multiplication works in `macpan2` exactly as in base R. + +```{r} +A <- matrix(1:6, 3, 2) +B <- matrix(1:4, 2, 2) +A %*% B +macpan2::engine_eval(~A %*% B, A = A, B = B) +``` + +## Kronecker Product (`%x%`) + +The Kronecker product builds structured block matrices, exactly as in base R. + +```{r} +A <- matrix(c(1, 2), nrow = 1) +B <- matrix(c(0, 1, 1, 0), nrow = 2) +A %x% B +macpan2::engine_eval(~A %x% B, A = A, B = B) ``` -but avoids ever constructing `LeftSparseMatrix` explicitly. You pass `z` as a pre-allocated matrix of the correct size, and `sparse_mat_mult()` fills it with the result of the multiplication. -This method is particularly useful in `macpan2` for: +Kronecker products are useful when replicating structure across groups or compartments, as the follow section illustrates. + +### Example -- Structured SI Model with Kronecker products + +We describe an SI model structured by age group, vaccination status, and symptom status, entirely in terms of matrices and vectors so that Kronecker products can be used to consistently replicate structure across these groups. + +#### Step 1 — Stratifying State Variables + +We divide a population into: + +- \( n_{\text{age}} \) age groups (indexed by \( a = 1, \dots, n_{\text{age}} \)) +- \( n_{\text{vax}} \) vaccination categories (indexed by \( v = 1, \dots, n_{\text{vax}} \)) +- \( n_{\text{symp}} \) symptom categories (indexed by \( s = 1, \dots, n_{\text{symp}} \)) + +All population quantities are represented as column vectors, stacked in this order: + +- Age group (leading index) — varies slowest +- Vaccination status — varies next +- Symptom status — varies fastest (only for infectious individuals) + +This convention ensures that aggregation and mixing operations defined by Kronecker products naturally correspond to the structure of the population. + +Susceptible vector \( S \) with \( n_{\text{age}} = 2 \), \( n_{\text{vax}} = 2 \): +\[ +S = \begin{bmatrix} +S_{1,1} \\ +S_{1,2} \\ +S_{2,1} \\ +S_{2,2} +\end{bmatrix} +\] + +Infectious vector \( I \) with \( n_{\text{age}} = 2 \), \( n_{\text{vax}} = 2 \), \( n_{\text{symp}} = 3 \): +\[ +I = \begin{bmatrix} +I_{1,1,1} \\ +I_{1,1,2} \\ +I_{1,1,3} \\ +I_{1,2,1} \\ +I_{1,2,2} \\ +I_{1,2,3} \\ +I_{2,1,1} \\ +I_{2,1,2} \\ +I_{2,1,3} \\ +I_{2,2,1} \\ +I_{2,2,2} \\ +I_{2,2,3} +\end{bmatrix} +\] + +#### Step 2 — Total Population by Age Group + +We construct summation matrices to aggregate the population counts into totals by age group. Each state variable requires its own summation matrix because they are structured differently: + +- **\( S \)** is structured by age and vaccination status. +- **\( I \)** is structured by age, vaccination status, and symptom status. + +To define the summation matrices (and expansion matrices below), we use: + +- \( \mathbb{I}_n \) — the \( n \times n \) identity matrix (ones on the diagonal, zeros elsewhere). +- \( \mathbf{1}_{1 \times n} \) — a \( 1 \times n \) row vector of ones. +- \( \mathbf{1}_{n \times 1} \) — an \( n \times 1 \) column vector of ones. + + +These matrices are combined using Kronecker products, ensuring summation within the correct strata. + +- For **susceptibles** \( S \), we sum over vaccination status while preserving age: +\[ +M_S = \mathbb{I}_{n_{\text{age}}} \otimes \mathbf{1}_{1 \times n_{\text{vax}}} +\] + +- For **infectious individuals** \( I \), we sum over both vaccination status and symptom status while preserving age: +\[ +M_I = \mathbb{I}_{n_{\text{age}}} \otimes \mathbf{1}_{1 \times (n_{\text{vax}} \cdot n_{\text{symp}})} +\] + +We then compute the total population per age group by applying the summation matrices: +\[ +N_{\text{age}} = M_S \cdot S + M_I \cdot I +\] + +This yields a column vector of length \( n_{\text{age}} \), giving the total population in each age group. + +#### Step 3 — Expanded Total Population Vector + +We compute a population vector \( N \) with the same structure and length as the infectious vector \( I \), so that it can be used directly in the force of infection expression below. + +We obtain \( N \) by expanding the vector of total population by age group, \( N_{\text{age}} \), using the Kronecker product with a vector of ones over the vaccination and symptom strata: +\[ +N = N_{\text{age}} \otimes \mathbf{1}_{n_{\text{vax}} \cdot n_{\text{symp}} \times 1} +\] + +Each entry of \( N \) gives the total population of the corresponding age group, repeated across all vaccination and symptom categories. + +#### Step 4 — WAIFW Matrix Construction + +The above notation allows us to very nicely partition the effects of each stratum on the *Who Acquires Infection From Whom* (WAIFW) matrix for our structured SI model, using Kronecker products of component WAIFW matrices: + +\[ +W = C_{\text{age}} \otimes C_{\text{vax}} \otimes C_{\text{symp}} +\] + +This representation reflects the structure of our state variables: susceptible individuals are stratified by age and vaccination status, whereas infectious individuals are stratified by age, vaccination status, and symptom status. The dimensions and interpretations of each component matrix follow from this: + +- \( C_{\text{age}} \in \mathbb{R}^{n_{\text{age}} \times n_{\text{age}}} \): + A square matrix representing transmission potential between susceptible individuals in each age group (rows) and infectious individuals in each age group (columns). This is typically a contact matrix, where entry \((i, j)\) encodes the probability that an individual in age group \(j\) will contact an individual in age group \(i\). + +- \( C_{\text{vax}} \in \mathbb{R}^{n_{\text{vax}} \times n_{\text{vax}}} \): + This matrix modulates the impact of vaccination status on transmission. Each entry \((i, j)\) quantifies the transmission potential from an infectious individual with vaccination status \(j\) to a susceptible individual with status \(i\). Often, this captures relative reductions in susceptibility or infectiousness due to vaccination. -- Reducing memory footprint when working with sparse structures -- Avoiding unnecessary computation with zero entries -- Embedding efficient matrix multiplications in a model specification +- \( C_{\text{symp}} \in \mathbb{R}^{1 \times n_{\text{symp}}} \): + Unlike age and vaccination status, susceptibles do not have a symptom status—one might say their symptom status is implicitly "none". Consequently, \( C_{\text{symp}} \) is a 1-by-\(n_{\text{symp}}\) row vector, where each entry adjusts the infectiousness of individuals depending on their symptom status (e.g., asymptomatic, presymptomatic, mildly symptomatic). This can encode, for example, that asymptomatic individuals are half as infectious as symptomatic ones. -We also plan to allow passing sparse matrices from the `Matrix` package directly into `macpan2` models, in the same way integer vectors can be passed today. The `Matrix` package works naturally with `TMB`, so this will allow us to leverage native sparse matrix support inside the `macpan2` engine. -Even with our current simple integer-vector-based approach, we have observed practically useful performance gains for large problems, especially when the left-hand matrix is highly sparse. +The full WAIFW matrix \( W \) has dimension \((n_{\text{age}} n_{\text{vax}}) \times (n_{\text{age}} n_{\text{vax}} n_{\text{symp}})\) and defines how each susceptible subgroup (by age and vaccination) acquires infection from each infectious subgroup (by age, vaccination, and symptoms), assuming homogeneous mixing within each cell of this stratified structure. -## Extracting Sparse Representation +#### Step 5 — Force of Infection Vector -To simplify the creation of sparse representations, `macpan2` provides the helper function `sparse_matrix_notation()`. -This function takes a dense matrix and returns: +We can now write the force of infection vector, with elements corresponding to the elements of the $S$ vector. -- `values` — the non-zero entries of the matrix -- `row_index` — the zero-based row indices of those entries -- `col_index` — the zero-based column indices of those entries -- `M` — the original matrix -- `Msparse` — a version of `M` with near-zero entries replaced by exact zeros +\[ +\Lambda = \beta \cdot W \cdot I \cdot \text{diag}^{-1}(N) +\] -By default, `sparse_matrix_notation()` applies a numerical tolerance to treat small values as zero and converts indices to zero-based, matching the requirements of `sparse_mat_mult()`. +#### Step 6 — Inflows and Outflows with Allocation Matrices + +In the structured SI model, the force of infection determines the rate at which individuals leave the susceptible compartment. However, because the infectious compartment \( I \) includes an additional stratification (symptom status), the force of infection does not specify how new infections are distributed among the strata of \( I \). We introduce allocation matrices to define this distribution explicitly and ensure consistent tracking of flows. + +We define: + +- \( \text{Out}_{S} \) — the outflow vector from the susceptible compartment, with the same structure as \( S \). +- \( \text{In}_{I} \) — the inflow vector into the infectious compartment, with the same structure as \( I \). + +The outflow from \( S \) due to infection is computed elementwise: +\[ +\text{Out}_{S} = \Lambda \circ S +\] +where \( \circ \) denotes elementwise multiplication. + +To compute the inflow vector we introduce the allocation matrix \( A_{S \to I} \), a \( (n_{\text{age}} \cdot n_{\text{vax}} \cdot n_{\text{symp}}) \times (n_{\text{age}} \cdot n_{\text{vax}}) \) matrix such that: + +- All elements of $A_{S \to I}$ are $\geq 0$. +- All columns sum to one. + +The inflow to \( I \) is then given by: +\[ +\text{In}_{I} = A_{S \to I} \cdot \text{Out}_{S} +\] + +The required properties of \( A_{S \to I} \) ensure that the outflow from each susceptible stratum is completely allocated among the infectious strata. Rows may sum to any non-negative value, including zero, which means that some infectious strata may receive no inflow directly from \( S \). + +#### Structured SI Model Specification -Example usage: ```{r} -library(macpan2) -# Define matrix to convert +``` + +## Elementwise Multiplication with Recycling (`*`) + +Elementwise multiplication applies entrywise and supports matrix recycling. + +Given $z = x * y$, the entries of $z$ are defined by: + +$$ +z_{i,j} = \begin{cases} + x_{i,j} \times y_{i,j} & \text{if } \dim(x) = \dim(y) \\ + x_{i,j} \times y_{i,1} & \text{if } n_{\text{row}}(x) = n_{\text{row}}(y),\ n_{\text{col}}(y) = 1 \\ + x_{i,j} \times y_{1,j} & \text{if } n_{\text{col}}(x) = n_{\text{col}}(y),\ n_{\text{row}}(y) = 1 \\ + x_{i,1} \times y_{i,j} & \text{if } n_{\text{row}}(x) = n_{\text{row}}(y),\ n_{\text{col}}(x) = 1 \\ + x_{1,j} \times y_{i,j} & \text{if } n_{\text{col}}(x) = n_{\text{col}}(y),\ n_{\text{row}}(x) = 1 \\ + x_{1,1} \times y_{i,j} & \text{if } x \text{ is } 1 \times 1 \\ + x_{i,j} \times y_{1,1} & \text{if } y \text{ is } 1 \times 1 \\ +\end{cases} +$$ + +See [this article](https://canmod.github.io/macpan2/articles/elementwise_binary_operators) for a general treatment of elementwise binary operators with examples. + +## Multiplying a Sparse Matrix with a Dense Matrix + +It is often necessary to perform matrix multiplications where the left-hand matrix is sparse—meaning most of its entries are zero. Passing large sparse matrices as dense matrices with zeros into `macpan2` models can lead to unnecessary memory use and performance slowdowns. + +To address this, `macpan2` provides the engine function `sparse_mat_mult()`, which multiplies a sparse matrix (given in compressed form) by a dense matrix. Instead of passing the full sparse matrix, you pass: + +- `x` — non-zero values +- `i` — zero-based row indices +- `j` — zero-based column indices +- `y` — dense right-hand matrix +- `z` — pre-allocated output matrix + +### Extracting Sparse Representation + +To simplify the creation of sparse representations, `macpan2` provides the helper function `sparse_matrix_notation()`. This function takes a dense matrix and returns: + +```{r} M <- matrix(c(5, 0, 0, 0, 0, 3, 0, 2, 0), nrow = 3, byrow = TRUE) -# Extract sparse representation sparse <- macpan2:::sparse_matrix_notation(M, zero_based = TRUE) print(sparse) -# Calculate sparsity percentage total_entries <- length(M) non_zero_entries <- length(sparse$values) sparsity <- 100 * (1 - non_zero_entries / total_entries) cat(sprintf("Matrix sparsity: %.1f%%\n", sparsity)) ``` -With this simple example, the matrix is 66.7% sparse. -The performance benefit of using sparse matrix multiplication grows with both matrix size and sparsity percentage. +### Small Example of Sparse Matrix Multiplication -## Example: Small Sparse Matrix Multiplication ```{r} -library(macpan2) - # Full matrix for validation A <- matrix(c(5, 0, 0, 0, 0, 3, @@ -113,12 +288,13 @@ spec <- mp_tmb_model_spec( ) # Run simulation for zero time steps to execute 'before' block -result <- spec |> - mp_simulator(time_steps = 0, outputs = "z") |> - mp_final_list() +result <- (spec + |> mp_simulator(time_steps = 0, outputs = "z") + |> mp_final_list() +) # Result from macpan2 -result$z +print(result$z) # Direct multiplication for validation A %*% y @@ -129,20 +305,17 @@ A %*% y - The function writes the result into `z`, which must be allocated beforehand. - You can use `mp_tmb_model_spec()` with `time_steps = 0` to evaluate this operation as part of a `macpan2` model workflow. -## Example: Time-Varying Transmission in an SI Model - +### Example of Sparse Matrix Multiplication to Model Time-Varying Transmission Rate + We implement an SI model where the transmission rate $\beta(t)$ varies over time, modeled as a spline basis expansion: $$ \log \beta(t) = B(t) \cdot \theta $$ We use `sparse_matrix_notation()` to encode the spline basis sparsely and compute $\log \beta(t)$ using `sparse_mat_mult()` in the `before` block. In this case, the spline basis contains no exact zeros, but many near-zero values that contribute little to the result. The tolerance ensures these negligible entries are dropped, improving sparsity without meaningfully affecting the result. For basis matrices like these, choosing a reasonable tolerance (e.g., `1e-4`) helps strike a balance between computational efficiency and model fidelity. -```{r, fig.height = 8, fig.width=4} -library(macpan2) -library(splines) -library(dplyr) -library(ggplot2) +We compute $\log \beta(t)$ from a sparse spline basis. +```{r, fig.height = 8, fig.width = 4} # Time points for simulation (100 steps) n_steps <- 100 time_points <- seq(0, 1, length.out = n_steps) @@ -170,37 +343,50 @@ beta_log <- matrix(0, nrow = n_steps, ncol = 1) spec <- mp_tmb_model_spec( before = beta_log ~ sparse_mat_mult(x, i, j, theta, beta_log), during = list( - beta ~ exp(beta_log[time_step(1)]), - mp_per_capita_flow( - from = "S", - to = "I", - rate = "beta * I / N", - flow_name = "infection" - ) + beta ~ exp(beta_log[time_step(1)]), + mp_per_capita_flow( + from = "S", + to = "I", + rate = "beta * I / N", + flow_name = "infection" + ) ), default = nlist(N = 100, S = 99, I = 1, x = B_sparse$values, theta, beta_log), integers = nlist(i = B_sparse$row_index, j = B_sparse$col_index) ) # Simulate for 100 time steps -outputs = c("S", "I", "beta", "infection") -result <- spec |> - mp_simulator(time_steps = n_steps, outputs = outputs) |> - mp_trajectory() |> - mutate(variable = factor(matrix, levels = outputs)) +result <- (spec + |> mp_simulator(time_steps = n_steps, outputs = c("S", "I", "beta", "infection")) + |> mp_trajectory() + |> mutate(variable = factor(matrix, levels = c("S", "I", "beta", "infection"))) +) # Plot the results -(result - |> ggplot() - + aes(time, value) - + geom_line() - + facet_wrap(~variable, scales = "free_y", ncol = 1) - + theme_bw() -) +ggplot(result, aes(time, value)) + + geom_line() + + facet_wrap(~variable, scales = "free_y", ncol = 1) + + theme_bw() ``` - Transmission rate $\beta(t)$ is computed dynamically using `sparse_mat_mult()`. - The spline basis is passed in sparse form, allowing efficient calculation of $\beta(t)$ at each time step. - The SI model uses `mp_per_capita_flow()` for infection, with time-varying transmission. - For this example, the spline basis (with tolerance $10^{-4}$) has approximately `r sprintf("%.1f%% sparsity", sparsity)`. - Larger models with higher sparsity will benefit even more from this approach. +- Larger models with higher sparsity will benefit more from this approach. + +## Summary of Matrix Operations + +| Operation | Description | macpan2 Support | +|-----------|-------------|------------------| +| `%*%` | Dense matrix multiplication | Direct | +| `%x%` | Kronecker product | Direct | +| `*` | Elementwise multiplication with recycling | Direct | +| `sparse_mat_mult()` | Sparse × Dense multiplication | Direct | + +## When to Use Each + +- **Dense × Dense (`%*%`)**: Use for small/moderate dense matrices. +- **Kronecker (`%x%`)**: Use for structured expansions (e.g., age). +- **Elementwise (`*`)**: Use for scaling rates in structured models (e.g., susceptibility in force of infection). +- **Sparse × Dense (`sparse_mat_mult()`)**: Use for large, sparse left-hand matrices. diff --git a/vignettes/state_dependent_rates.Rmd b/vignettes/state_dependent_rates.Rmd index e74f46e41..8dfc66df3 100644 --- a/vignettes/state_dependent_rates.Rmd +++ b/vignettes/state_dependent_rates.Rmd @@ -20,11 +20,198 @@ knitr::opts_chunk$set( ) ``` + +Choose one susceptible and one infectious candidate from the population. + +* S-candidate is susceptible, $S_S$ +* S-candidate is in age group, $S_{A_i}$ +* S-candidate is in immunity status, $S_{V_j}$ +* I-candidate is infectious, $I_I$ +* I-candidate is in age group, $I_{A_i}$ +* I-candidate is in immunity status, $I_{V_j}$ + +We only consider $\text{pr}(S_S) = 1$. + +We let $q_{jk} = \text{Pr}(I_{V_j} | S_{V_k})$ + + + +In most compartmental models, flow rates depend not only on the compartments directly involved in the flow but also on other compartments. A typical example is infection: the rate at which individuals move from a susceptible compartment to an exposed compartment depends on the number of infectious individuals. We frame this setup generically, with \( S \) representing the source (“from”) compartments, \( E \) the destination (“to”) compartments, and \( I \) a set of other compartments whose states influence the flow. While infection is the motivating example, this structure also captures processes like vaccination, where uptake rates depend on vaccine supply, or treatment initiation, where rates depend on healthcare capacity. + +We define +$$ +S = (S_1, \ldots, S_n)^\top, \quad E = (E_1, \ldots, E_q)^\top, \quad I = (I_1, \ldots, I_m)^\top, +$$ +where each element represents the number of individuals in a specific stratum (e.g., by age, location, symptom status, or immunity level). The stratifications of \( S \), \( E \), and \( I \) may differ. For example, \( S \) and \( E \) might be stratified by age alone, while \( I \) is stratified jointly by age and symptom severity. As a result, \( n \), \( q \), and \( m \) are generally distinct. In special cases—such as omitting a latent period—one may simplify by setting \( E = I \). + +--- + +### Force of Infection and Flow Allocation + +The **force of infection** is given by +$$ +\Lambda = \beta I, +$$ +where \( \beta \) is an \( n \times m \) transmission matrix with entries \( \beta_{ij} \) representing the per-capita contribution of infectious individuals in stratum \( j \) to the infection risk of susceptibles in stratum \( i \). + +The vector of outflows from \( S \), denoted \( F^{\text{out}} \), is +$$ +F^{\text{out}} = \Lambda \circ S, +$$ +where \( \circ \) denotes elementwise multiplication. The vector \( F^{\text{out}} \) has length \( n \), matching \( S \). + +To account for different stratifications of \( S \) and \( E \), we define an **allocation matrix** \( M \) of size \( q \times n \), where \( M_{ki} \) is the proportion of individuals leaving \( S_i \) who enter \( E_k \). The inflow into \( E \), denoted \( F^{\text{in}} \), is given by +$$ +F^{\text{in}} = M F^{\text{out}}. +$$ +Here, \( F^{\text{in}} \) is a length-\( q \) vector matching \( E \). This is standard matrix multiplication: each entry of \( F^{\text{in}} \) is a weighted sum of the outflows from \( S \), with weights given by \( M \). By construction, the total flow is conserved: +$$ +\sum_{k=1}^q F^{\text{in}}_k = \sum_{i=1}^n F^{\text{out}}_i. +$$ + +--- + +### Stratified Transmission and Multiplicative Decomposition + +When compartments are stratified, heterogeneity in susceptibility, contact patterns, and infectivity can be captured by factorizing the transmission matrix \( \beta \). We write: +$$ +\beta = \operatorname{diag}(\sigma) \, P \, \operatorname{diag}(\tau), +$$ +where: +- \( \sigma \) is a length‑\( n \) vector of *susceptibilities* for each \( S \) stratum, +- \( P \) is an \( n \times m \) *contact matrix* describing contact rates between \( S \) and \( I \) strata, +- \( \tau \) is a length‑\( m \) vector of *infectivities* for each \( I \) stratum. + +The **force of infection** is then +$$ +\Lambda = \beta I, +$$ +and the absolute outflow from each \( S \) compartment is +$$ +F^{\text{out}} = \Lambda \circ S. +$$ + +This decomposition emphasizes how susceptibility, contact rates, and infectivity combine multiplicatively to determine transmission dynamics across strata. + +--- + +### Relation to the Standard SI Model + +In a standard unstratified SI model, there is a single susceptible compartment \( S \) and a single infectious compartment \( I \), with total population \( N = S + I \). The force of infection is typically written as +$$ +\Lambda = \beta \frac{I}{N}, +$$ +so the rate at which susceptibles become infected is \( \Lambda \cdot S = \beta \frac{I}{N} \cdot S \). + +In terms of our general framework: +- \( n = q = m = 1 \), +- \( \sigma = 1 \), +- \( P = 1 \), +- \( \tau = \beta / N \), + +so that +$$ +\beta = \sigma \cdot P \cdot \tau = 1 \cdot 1 \cdot \frac{\beta}{N} = \frac{\beta}{N}. +$$ + +The outflow from \( S \) is +$$ +F^{\text{out}} = \Lambda \cdot S = \frac{\beta}{N} \cdot I \cdot S, +$$ +matching the standard mass-action incidence. + +Thus, the classical SI model corresponds to the special case where there is no stratification, susceptibility is homogeneous, contact is uniform, and infectious individuals transmit at a per-contact rate scaled by \( 1/N \). + +### Product Models for Modular Stratification + +The framework above allows for arbitrary stratifications of the \( S \), \( E \), and \( I \) compartments. In practice, complex stratifications are often built by combining simpler ones — for example, stratifying by age, location, and symptom status simultaneously. + +We formalize this using **product models**, where the full stratification is represented as the Cartesian product of simpler factors (e.g., age groups \( \times \) locations \( \times \) symptom statuses). Each compartment corresponds to a unique combination of levels across these factors. + +For example, suppose we stratify by age (child, adult) and location (urban, rural). The susceptible vector \( S \) then has four entries: +$$ +S = \begin{pmatrix} S_{\text{child, urban}} \\ S_{\text{child, rural}} \\ S_{\text{adult, urban}} \\ S_{\text{adult, rural}} \end{pmatrix}. +$$ +Parameters like susceptibility or contact rates can vary by age alone, location alone, or both, depending on model design. + +Product models allow us to: +- Define parameters that vary along specific dimensions, +- Build large, structured transmission models from simpler components, +- Control complexity by selectively including or collapsing stratification dimensions. + +The force of infection and flow allocation machinery developed above applies directly to these multi-dimensional strata. + + +### Example: Building a Transmission Matrix with Kronecker Products + +When stratifications are defined by a product of factors, transmission matrices can often be constructed as Kronecker products of simpler matrices associated with each factor. + +For example, suppose: +- Contact patterns differ by **age** (with a \( 2 \times 2 \) contact matrix \( P_{\text{age}} \)), +- Individuals interact equally across **locations** (represented by a \( 2 \times 2 \) matrix of ones, \( 1_{2 \times 2} \)). + +The full contact matrix \( P \) for the combined age-location stratification (with 4 age-location groups) is given by the Kronecker product: +$$ +P = P_{\text{age}} \otimes 1_{2 \times 2}. +$$ +This expands \( P_{\text{age}} \) uniformly across locations, assuming location has no effect on contact rates. + +If location-specific effects were also important, we could replace \( 1_{2 \times 2} \) with a location contact matrix \( P_{\text{loc}} \), giving +$$ +P = P_{\text{age}} \otimes P_{\text{loc}}. +$$ + +This approach lets us modularly build structured transmission matrices by combining simple factors — mirroring the way stratifications are built in the product model. + +----- + +#### Product Expansion + +The full susceptibility vector \( \sigma \), contact matrix \( P \), and infectivity vector \( \tau \) are given by **Kronecker products** of their factor-specific components: +$$ +\sigma = \sigma^{(1)} \otimes \sigma^{(2)} \otimes \cdots \otimes \sigma^{(d)}, +$$ +$$ +P = P^{(1)} \otimes P^{(2)} \otimes \cdots \otimes P^{(d)}, +$$ +$$ +\tau = \tau^{(1)} \otimes \tau^{(2)} \otimes \cdots \otimes \tau^{(d)}. +$$ + +#### Transmission Matrix Construction + +The overall transmission matrix \( \beta \) is then constructed exactly as before: +$$ +\beta = \operatorname{diag}(\sigma) \, P \, \operatorname{diag}(\tau). +$$ + +This formulation means that: +- Susceptibility effects, +- Contact patterns, and +- Infectivity effects + +are each **structured by Kronecker products reflecting the model’s stratification factors** — and then combined multiplicatively through this decomposition. + +--- + +This triple-Kronecker structure allows the transmission process to respect the modular stratification of the product model while retaining the familiar separation of susceptibility, contact, and infectivity. + + + +------------ + ## Decomposition of State-Dependent Rates -The vast majority of compartmental models contain flows that depend on the states of compartments that are not directly involved in those flows. For example, the magnitude of the flow from compartment A to compartment B could depend on another compartment, C. The most important example of these dependencies is infection, and so we use the terminology from infection processes. However, our model will subsume a very large number of flows as special cases. +The vast majority of compartmental models contain flows that depend on the states of compartments that are not directly involved in those flows (i.e., the magnitude of the flow from compartment A to compartment B depends on another compartment, C). The most important example of these dependencies is infection, and so we use the terminology from infection processes. However, our model will subsume a very large number of flows as special cases. + +Let \( S = (S_1, \ldots, S_n)^\top \), \( E = (E_1, \ldots, E_q)^\top \), and \( I = (I_1, \ldots, I_m)^\top \) denote the vectors of numbers of susceptible, exposed, and infectious individuals, respectively. Each element of these vectors corresponds to a specific stratum—such as an age group, location, symptom status, or immunity level—and the stratifications may differ between \( S \), \( E \), and \( I \). Consequently, \( n \), \( q \), and \( m \) may be different. + +The force of infection vector, $\Lambda$, is $\beta I$, where beta is the transmission matrix (show the elements $\beta_ij$ out please). + +The absolute rate of flow out of each of the $S$ states is $\Lambda \circ S$ (ethat's elementwise mult). + +Because $n$ and $q$ are not necessarily the same, the absolute rate of flow out of -Let $s$ and $x$ be subsets of the state vector. The length-$n$ vector $s$ contains the from-compartments -- typically the compartments that are susceptible to infection -- for a set of flows. The length-$m$ vector $x$ contains compartments -- typically the infectious compartments -- that affect these flows. The components of $x$ do not need to affect all components of $y$ but they do need to affect some of them. Note that we do not keep track of the to-compartments in these flows because they do not matter (I think). We define the $n$-by-$m$ transmission matrix as the following decomposition. From c8f9598512a630ae066c72d44d837ae961bc9307 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Fri, 25 Jul 2025 23:53:21 -0400 Subject: [PATCH 09/13] readme --- README.Rmd | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/README.Rmd b/README.Rmd index 7f806047b..79d495a40 100644 --- a/README.Rmd +++ b/README.Rmd @@ -18,6 +18,7 @@ library(dplyr) theme_bw = function() ggplot2::theme_bw(base_size = 18) ``` + [![macpan2 status badge](https://canmod.r-universe.dev/badges/macpan2)](https://canmod.r-universe.dev/macpan2) [![R-CMD-check](https://github.com/canmod/macpan2/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/canmod/macpan2/actions/workflows/R-CMD-check.yaml) @@ -26,7 +27,8 @@ theme_bw = function() ggplot2::theme_bw(base_size = 18) [![commit activity](https://img.shields.io/github/commit-activity/m/canmod/macpan2)](https://github.com/canmod/macpan2/commits) [![contributors](https://img.shields.io/github/contributors/canmod/macpan2)](https://github.com/canmod/macpan2/graphs/contributors) -![logo](man/figures/logo.png) +![](man/figures/logo.png) + `macpan2` is a software platform for building and calibrating compartmental models of infectious disease. It supports flexible model specification and fast parameter calibration, making it easier for modellers to respond to emerging public health threats. Developed through collaboration with the [Public Health Agency of Canada (PHAC)](https://www.canada.ca/en/public-health.html), `macpan2` is being used to support responses to diseases such as mpox [@milwid2023mpox], measles [@phac2024measles, [interactive measles model](https://wzmli.shinyapps.io/two_pop_measles_shiny/)], and COVID-19 [@simmons2025cost;@miranda2024strategies], and internal work on [pandemic preparedness](https://phac-nml-phrsd.github.io/EPACmodel/). From afb25b2c04737de4425d83dfb221d1ea8475af3c Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Wed, 30 Jul 2025 15:39:42 -0400 Subject: [PATCH 10/13] fix tests --- R/mp_tmb_model_spec.R | 16 +- R/tmb_model.R | 15 - inst/starter_models/hiv/README.R | 4 +- inst/starter_models/macpan_base/README.R | 310 +++++++++++++++++- inst/starter_models/nfds/README.R | 2 +- inst/starter_models/shiver/README.R | 30 +- inst/starter_models/sir_mosquito/README.R | 2 + inst/starter_models/ww/README.R | 4 +- man/mp_kronecker_operator.Rd | 20 ++ man/mp_square_operator.Rd | 19 ++ tests/testthat.R | 7 - tests/testthat/setup.R | 2 + tests/testthat/test-c.R | 1 - tests/testthat/test-calibrate-bad-specs.R | 8 +- tests/testthat/test-calibrator-par.R | 1 - tests/testthat/test-calibrator-time.R | 1 - tests/testthat/test-calibrator-traj.R | 1 - tests/testthat/test-calibrator-tv.R | 1 - tests/testthat/test-change-model.R | 1 - tests/testthat/test-clamp.R | 1 - tests/testthat/test-convolution.R | 1 - tests/testthat/test-date-time-steps.R | 1 - .../testthat/test-dimnames-and-shape-change.R | 3 - tests/testthat/test-distributions.R | 1 - tests/testthat/test-dynamic-model.R | 3 - tests/testthat/test-index.R | 1 - tests/testthat/test-insert-state-vars.R | 12 +- tests/testthat/test-list-utils.R | 1 - tests/testthat/test-model-modifications.R | 1 - tests/testthat/test-model-structure.R | 2 - tests/testthat/test-mp-rk4.R | 1 - tests/testthat/test-partitions.R | 2 - tests/testthat/test-rbind-time-lag.R | 1 - tests/testthat/test-scalar-2-vector.R | 2 +- tests/testthat/test-state-flow-order.R | 3 - tests/testthat/test-tmb-model-spec.R | 1 - tests/testthat/test-tmb.R | 4 +- 37 files changed, 390 insertions(+), 96 deletions(-) create mode 100644 man/mp_kronecker_operator.Rd create mode 100644 man/mp_square_operator.Rd diff --git a/R/mp_tmb_model_spec.R b/R/mp_tmb_model_spec.R index 265608565..7ffff5911 100644 --- a/R/mp_tmb_model_spec.R +++ b/R/mp_tmb_model_spec.R @@ -100,20 +100,12 @@ TMBModelSpec = function( self$all_integers = function() { ## TODO: maybe make smarter by checking if an integer vector ## is being used in the wrong numeric vector - implied_integers = implied_position_vectors(self$default) - integers_we_need = intersect( + implied_integers = implied_position_vectors(self$default) ## named list of integers + integers_we_need = intersect( ## character vector of names of integers to add self$all_formula_vars() , names(implied_integers) - ) - filtered_implied_integers = list() - nms = names(integers_we_need) - for (nm in unique(nms)) { - integers_nm = integers_we_need[nms == nm] - # if (sum(!duplicated(integers_nm)) != 1L) { - # } - filtered_implied_integers[[nm]] = integers_nm[[1L]] - } - c(filtered_implied_integers, self$integers) + ) |> unique() + c(implied_integers[integers_we_need], self$integers) } self$empty_matrices = function() { diff --git a/R/tmb_model.R b/R/tmb_model.R index 8cb571cf9..2508dee43 100644 --- a/R/tmb_model.R +++ b/R/tmb_model.R @@ -52,21 +52,6 @@ mp_simulator.TMBModelSpec = function(model model$simulator_fresh(time_steps, outputs, default) } -# mp_simulator.TMBSimulator = function(model -# , time_steps -# , outputs -# , default = list() -# , inits = list() -# ) { -# stop("under construction") -# if (!missing(time_steps)) { -# model$simulator$replace$time_steps(time_steps) -# } -# ## TODO: -# ## set the params vector as the last best params vector -# ## update the outputs -# } - #' @export mp_simulator.TMBCalibrator = function(model , time_steps diff --git a/inst/starter_models/hiv/README.R b/inst/starter_models/hiv/README.R index 050fb3396..b4d74a84b 100644 --- a/inst/starter_models/hiv/README.R +++ b/inst/starter_models/hiv/README.R @@ -96,8 +96,8 @@ calibrator = (spec_for_cal |> mp_tmb_calibrator( data = simulated_data , traj = list( - treated = mp_poisson() - , untreated = mp_poisson() + treated = mp_pois() + , untreated = mp_pois() ) , par = c("log_lambda0", "logit_n") ) diff --git a/inst/starter_models/macpan_base/README.R b/inst/starter_models/macpan_base/README.R index b503f8275..6650ab4a2 100644 --- a/inst/starter_models/macpan_base/README.R +++ b/inst/starter_models/macpan_base/README.R @@ -2,7 +2,9 @@ knitr::opts_chunk$set( collapse = TRUE, comment = "#>", - fig.path = "./figures/" + fig.path = "./figures/", + message = FALSE, + warning = FALSE ) round_coef_tab = function(x, digits = 4) { id_cols = c("term", "mat", "row", "col", "type") @@ -24,6 +26,10 @@ library(ggplot2) library(dplyr) +## ----options------------------------------------------------------------------ +options(macpan2_verbose = FALSE) + + ## ----model_lib---------------------------------------------------------------- spec = mp_tmb_library( "starter_models" @@ -37,3 +43,305 @@ system.file("utils", "box-drawing.R", package = "macpan2") |> source() layout = mp_layout_paths(spec, x_gap = 0.1, y_gap = 0.1) plot_flow_diagram(layout) + +## ----get_data----------------------------------------------------------------- +ts_data = ("https://github.com/canmod/macpan2" + |> file.path("releases/download/macpan1.5_data/covid_on.RDS") + |> url() + |> readRDS() +) + + +## ----prep_ts------------------------------------------------------------------ +prepped_ts_data = (ts_data + |> rename(matrix = var) + # dates from base model calibration (Figure 4) + |> filter(date >= "2020-02-24" & date < "2020-08-31") + # create unique time identifier + |> arrange(date) + |> group_by(date) + |> mutate(time = cur_group_id()) + |> ungroup() + # one negative value for daily deaths (removing for now time==178) + # this explains negative values: + # https://github.com/ccodwg/Covid19Canada?tab=readme-ov-file#datasets + |> filter(matrix %in% c("death","report") & value >= 0) +) + + +## ----prep_mob----------------------------------------------------------------- +prepped_mobility_data = (ts_data + |> filter(var == "mobility_index") + # dates from base model calibration (Figure 4) + |> filter(date >= "2020-02-24" & date < "2020-08-31") + # create unique time identifier + |> arrange(date) + |> group_by(date) + |> mutate(time = cur_group_id()) + |> ungroup() + |> mutate(log_mobility_ind = log(value)) +) + + +## ----mob_breaks--------------------------------------------------------------- +mobility_breaks = (prepped_ts_data + |> filter(date %in% c("2020-04-01", "2020-08-07"), matrix == "report") + |> pull(time) +) + + +## ----time_steps--------------------------------------------------------------- +time_steps = nrow(prepped_ts_data %>% select(date) %>% unique()) +prepped_ts_data = select(prepped_ts_data, -date) + + +## ----mobility_matrix---------------------------------------------------------- +S_j = function(t, tau_j, s) 1/(1 + exp((t - tau_j) / s)) +X = cbind( + prepped_mobility_data$log_mobility_ind + , S_j(1:time_steps, mobility_breaks[1], 3) + , S_j(1:time_steps, mobility_breaks[1], 3) * prepped_mobility_data$log_mobility_ind + , S_j(1:time_steps, mobility_breaks[2], 3) + , S_j(1:time_steps, mobility_breaks[2], 3) * prepped_mobility_data$log_mobility_ind + ) %>% as.matrix() +X_sparse = macpan2:::sparse_matrix_notation(X, TRUE) +model_matrix_values = X_sparse$values +row_ind = X_sparse$row_index +col_ind = X_sparse$col_index + + +## ----focal_model-------------------------------------------------------------- +focal_model = (spec + # add variable transformations: + |> mp_tmb_insert(phase = "before" + , at = 1L + , expressions = list( + zeta ~ exp(log_zeta) + , beta0 ~ exp(log_beta0) + , nonhosp_mort ~ 1/(1+exp((-logit_nonhosp_mort))) + , E ~ exp(log_E) + ) + , default = list( + log_zeta = empty_matrix + , log_beta0 = empty_matrix + , logit_nonhosp_mort = empty_matrix + , log_E = empty_matrix + ) + + ) + + # add accumulator variables: + # death - new deaths each time step + |> mp_tmb_insert(phase = "during" + , at = Inf + , expressions = list(death ~ ICUd.D + Is.D) + ) + + # add phenomenological heterogeneity: + |> mp_tmb_update(phase = "during" + , at = 1L + , expressions = list( + mp_per_capita_flow( + "S", "E" + , "((S/N)^zeta) * (beta / N) * (Ia * Ca + Ip * Cp + Im * Cm * (1 - iso_m) + Is * Cs *(1 - iso_s))" + , "S.E" + )) + + ) + + # compute gamma-density delay kernel for convolution: + |> mp_tmb_insert_reports("S.E" + , report_prob = 0.1 + , mean_delay = 11 + , cv_delay = 0.25 + , reports_name = "report" + ) + + # add time-varying transmission with mobility data: + |> mp_tmb_insert(phase = "before" + , at = Inf + , expressions = list(relative_beta_values ~ group_sums(model_matrix_values * log_mobility_coefficients[col_ind], row_ind, model_matrix_values)) + , default = list(log_mobility_coefficients = empty_matrix, model_matrix_values = empty_matrix) + , integers = list(row_ind = row_ind, col_ind = col_ind) + ) + |> mp_tmb_insert(phase = "during" + , at = 1L + , expressions = list( + beta ~ exp(log_beta0 + relative_beta_values[time_step(1)]) + ) + ) + +) + + +## ----calibrator--------------------------------------------------------------- +focal_calib = mp_tmb_calibrator( + spec = focal_model |> mp_rk4() + , data = prepped_ts_data + , traj = c("report","death") + , par = c( + "log_zeta" + , "log_beta0" + , "logit_nonhosp_mort" + , "log_mobility_coefficients" + , "log_E" + # negative binomial dispersion parameters for reports and deaths get added + # automatically with options(macpan2_default_loss = "neg_bin") set above + ) + , outputs = c("death","report") + , default = list( + + # states + # Population of Ontario (2019) from: + # https://github.com/mac-theobio/macpan_base/blob/main/code/ontario_calibrate_comb.R + S = 14.57e6 - 5 + , log_E = log(5) + , Ia = 0 + , Ip = 0 + , Im = 0 + , Is = 0 + , R = 0 + , H = 0 + , ICUs = 0 + , ICUd = 0 + , H2 = 0 + , D = 0 + + , model_matrix_values = model_matrix_values + + # set initial parameter values for optimizer + + , log_beta0 = log(5) + , logit_nonhosp_mort = -0.5 + , log_mobility_coefficients = rep(0, 5) + , log_zeta = 1 + + ) +) +# converges +mp_optimize(focal_calib) + + +## ----get_trajectory----------------------------------------------------------- +fitted_data = mp_trajectory_sd(focal_calib, conf.int = TRUE) + + +## ----coefficients------------------------------------------------------------- +mp_tmb_coef(focal_calib, conf.int = TRUE) |> round_coef_tab() + + +## ----plot_fit----------------------------------------------------------------- +(ggplot(prepped_ts_data, aes(time,value)) + + geom_point() + + geom_line(aes(time, value) + , data = fitted_data |> filter(matrix %in% c("death","report")) + , colour = "red" + ) + + geom_ribbon(aes(time, ymin = conf.low, ymax = conf.high) + , data = fitted_data |> filter(matrix %in% c("death","report")) + , alpha = 0.2 + , colour = "red" + ) + + facet_wrap(vars(matrix),scales = 'free') + + theme_bw() +) + + +## ----cohort_model------------------------------------------------------------- +cohort_model_ph = ( + focal_model + %>% mp_tmb_update(phase="during" + , at = 2L + , expressions = list(S.E ~ (S^zeta) * (beta / (N^(zeta) * N_cohort)) * (Ia * Ca + Ip * Cp + Im * Cm * (1 - iso_m) + Is * Cs *(1 - iso_s))) + , default = list(S.E = 0)) +) + + +## ----simulate_cohort---------------------------------------------------------- +cohort_sim_ph = (mp_simulator(cohort_model_ph + , time_steps = 100L + , outputs = "S.E" + , default = list( + S = 14.57e6 - 1 + , log_E = log(1) + , Ia = 0 + , Ip = 0 + , Im = 0 + , Is = 0 + , R = 0 + , H = 0 + , ICUs = 0 + , ICUd = 0 + , H2 = 0 + , D = 0 + , N_cohort = 1 + , model_matrix_values = model_matrix_values + , qmax = 34 + , log_beta0=log(1) + , logit_nonhosp_mort = -0.5 + , log_mobility_coefficients = rep(0,5) + , log_zeta = 1 + ) + ) |> mp_trajectory() +) + + +## ----R0----------------------------------------------------------------------- +R0_ph = sum(cohort_sim_ph$value) +print(R0_ph) + + +## ----cohort_no_ph------------------------------------------------------------- +cohort_model = (focal_model + %>% mp_tmb_update(phase="during" + , at = 2L + , expressions = list(mp_per_capita_flow( + "S", "E" + , "(beta / N) * (Ia * Ca + Ip * Cp + Im * Cm * (1 - iso_m) + Is * Cs * (1 - iso_s))" + , "S.E" + ) + )) + %>% mp_tmb_insert(phase="during" + , at = Inf + , expressions = list(foi ~ (beta / N) * (Ia * Ca + Ip * Cp + Im * Cm * (1 - iso_m) + Is * Cs *(1 - iso_s))) + ) +) + +cohort_sim = (mp_simulator(cohort_model + , time_steps = 100L + , outputs = "foi" + , default = list( + S = 0 + , log_E = log(1) + , Ia = 0 + , Ip = 0 + , Im = 0 + , Is = 0 + , R = 0 + , H = 0 + , ICUs = 0 + , ICUd = 0 + , H2 = 0 + , D = 0 + , N = 1 + , model_matrix_values = model_matrix_values + , qmax = 34 + , log_beta0=log(1) + , logit_nonhosp_mort = -0.5 + , log_mobility_coefficients = rep(0,5) + , log_zeta = 1) # not actually used + ) |> mp_trajectory() +) + +R0 = sum(cohort_sim$value) +print(R0) + + +## ----check_R0----------------------------------------------------------------- +all.equal(R0_ph, R0) + + +## ----euler_lotka-------------------------------------------------------------- +euler_lotka = function(r) sum(cohort_sim$value * exp(-r * cohort_sim$time)) - 1 +uniroot(euler_lotka, c(0,10)) + diff --git a/inst/starter_models/nfds/README.R b/inst/starter_models/nfds/README.R index a779d256c..ad4d3f680 100644 --- a/inst/starter_models/nfds/README.R +++ b/inst/starter_models/nfds/README.R @@ -106,7 +106,7 @@ binary_matrix_notation <- function(M){ # get row indices row_index = as.integer(rep(1:nrow(M), times=rowSums(M))-1) - return(nlist(col_index,row_index)) + return(macpan2:::nlist(col_index,row_index)) } # get non-zero indices of G transpose (loci-genotype matrix) diff --git a/inst/starter_models/shiver/README.R b/inst/starter_models/shiver/README.R index 07c14232a..a1f477a03 100644 --- a/inst/starter_models/shiver/README.R +++ b/inst/starter_models/shiver/README.R @@ -318,11 +318,11 @@ print(reparameterized_spec) ## ----reparam_calib------------------------------------------------------------ prior_distributions = list( - log_beta = mp_uniform() - , log_E_I_ratio = mp_uniform() - , logit_p = mp_normal(qlogis(1/4), 8) - , sigma = mp_uniform() - , gamma_h = mp_uniform() + log_beta = mp_unif() + , log_E_I_ratio = mp_unif() + , logit_p = mp_norm(qlogis(1/4), 8) + , sigma = mp_unif() + , gamma_h = mp_unif() ) shiver_calibrator = mp_tmb_calibrator( spec = reparameterized_spec @@ -451,14 +451,14 @@ multi_traj_spec = (reparameterized_spec sd_par = 1 ## for convenience we give all parameters the same prior sd, for now sd_state = 4 ## extremely vague priors on state variables prior_distributions = list( - log_beta = mp_normal(log(0.2), sd_par) - , log_sigma = mp_normal(log(sigma), sd_par) - , log_gamma_h = mp_normal(log(0.07), sd_par) - , logit_report_prob = mp_normal(qlogis(0.1), sd_par) - , logit_p = mp_normal(qlogis(1/4), 4) - , log_E_I_ratio = mp_normal(0, sd_par) - , log_I = mp_normal(log(I0), sd_state) - , log_H = mp_normal(log(H0), sd_state) + log_beta = mp_norm(log(0.2), sd_par) + , log_sigma = mp_norm(log(sigma), sd_par) + , log_gamma_h = mp_norm(log(0.07), sd_par) + , logit_report_prob = mp_norm(qlogis(0.1), sd_par) + , logit_p = mp_norm(qlogis(1/4), 4) + , log_E_I_ratio = mp_norm(0, sd_par) + , log_I = mp_norm(log(I0), sd_state) + , log_H = mp_norm(log(H0), sd_state) ) ## put the data together @@ -473,8 +473,8 @@ shiver_calibrator = mp_tmb_calibrator( # (changed from negative binomial because apparently it is easier # to fit standard deviations than dispersion parameters) , traj = list( - H = mp_normal(sd = mp_fit(1)) - , reported_incidence = mp_normal(sd = mp_fit(1)) + H = mp_norm(sd = mp_fit(1)) + , reported_incidence = mp_norm(sd = mp_fit(1)) ) , par = prior_distributions # fit the transmission rate using five radial basis functions for diff --git a/inst/starter_models/sir_mosquito/README.R b/inst/starter_models/sir_mosquito/README.R index 921780da4..d2e2126c2 100644 --- a/inst/starter_models/sir_mosquito/README.R +++ b/inst/starter_models/sir_mosquito/README.R @@ -1,6 +1,8 @@ ## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, + warning = FALSE, + message = FALSE, comment = "#>", fig.path = "./figures/" ) diff --git a/inst/starter_models/ww/README.R b/inst/starter_models/ww/README.R index bcc008eb3..cdb6ababc 100644 --- a/inst/starter_models/ww/README.R +++ b/inst/starter_models/ww/README.R @@ -161,8 +161,8 @@ focal_calib = mp_tmb_calibrator( , data = obs_data , traj = list( # set likelihoods for trajectories we are fitting to - reported_incidence = mp_neg_bin(disp = 0.1) - , W = mp_normal(sd = 1) + reported_incidence = mp_nbinom(disp = 0.1) + , W = mp_norm(sd = 1) ) , par = c( # parameters to fit diff --git a/man/mp_kronecker_operator.Rd b/man/mp_kronecker_operator.Rd new file mode 100644 index 000000000..edc9e2d12 --- /dev/null +++ b/man/mp_kronecker_operator.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/binary_operator.R +\name{mp_kronecker_operator} +\alias{mp_kronecker_operator} +\title{Kronecker Operator} +\usage{ +mp_kronecker_operator(operator) +} +\arguments{ +\item{operator}{A scalar binary operator.} +} +\value{ +A Kronecker operator convenient for use with \code{macpan2}. +} +\description{ +Convert a function that represents a scalar binary +operator into one that computes the \code{\link{kronecker}} +version, but with dimensions named in a way that is more +convenient for use with \code{macpan2}. +} diff --git a/man/mp_square_operator.Rd b/man/mp_square_operator.Rd new file mode 100644 index 000000000..cf041ba64 --- /dev/null +++ b/man/mp_square_operator.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/binary_operator.R +\name{mp_square_operator} +\alias{mp_square_operator} +\title{Square Matrix Operator} +\usage{ +mp_square_operator(operator) +} +\arguments{ +\item{operator}{A unary operator of a vector.} +} +\value{ +An operator. +} +\description{ +Convert a unary operator that takes a vector, into one with +dimensions named in a way that is more convenient for use with +\code{macpan2}. +} diff --git a/tests/testthat.R b/tests/testthat.R index 2098c2921..22c989be3 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -6,11 +6,4 @@ # * https://r-pkgs.org/tests.html # * https://testthat.r-lib.org/reference/test_package.html#special-files -library(testthat) -library(macpan2) -library(dplyr) -library(tidyr) -library(ggplot2) -library(broom.mixed) - test_check("macpan2") diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index d9de5a949..ed95cd3d8 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -6,6 +6,8 @@ suppressPackageStartupMessages({ library(ggplot2) library(broom.mixed) library(deSolve) + library(lubridate) + library(oor) }) options(macpan2_verbose = FALSE) system.file("utils", "test-cache.R", package = "macpan2") |> source() diff --git a/tests/testthat/test-c.R b/tests/testthat/test-c.R index 1a831e0ec..2435f23ad 100644 --- a/tests/testthat/test-c.R +++ b/tests/testthat/test-c.R @@ -22,7 +22,6 @@ test_that("concatenation works with many different shapes of input", { }) test_that("String objects can be concatenated", { - library(testthat); library(macpan2) x = macpan2:::StringUndottedVector("S", "E", "I", "R") y = macpan2:::StringUndottedVector("D") z = macpan2:::StringUndottedVector("S", "E", "I", "R", "D") diff --git a/tests/testthat/test-calibrate-bad-specs.R b/tests/testthat/test-calibrate-bad-specs.R index 71a670f66..bcf5cad3b 100644 --- a/tests/testthat/test-calibrate-bad-specs.R +++ b/tests/testthat/test-calibrate-bad-specs.R @@ -1,6 +1,8 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) -test_that("vectors and vector elements cannot have the same name.", { - spec = mp_tmb_model_spec(default = list(a = 0, b = c(a = 0))) +test_that("vectors and vector elements cannot have the same name, if they are used", { + spec = mp_tmb_model_spec( + before = list(a ~ b[a]) + , default = list(a = 0, b = c(a = 0)) + ) expect_error( mp_tmb_calibrator(spec, empty_trajectory) , "The following names were used for one or more purposes" diff --git a/tests/testthat/test-calibrator-par.R b/tests/testthat/test-calibrator-par.R index 7e4753558..4104e4ade 100644 --- a/tests/testthat/test-calibrator-par.R +++ b/tests/testthat/test-calibrator-par.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("bad parameterizations give errors", { sir = mp_tmb_library("starter_models", "sir", package = "macpan2") sir_sim = test_cache_read("TRAJ-sir_5_infection.rds") diff --git a/tests/testthat/test-calibrator-time.R b/tests/testthat/test-calibrator-time.R index 294b12a71..f94c07127 100644 --- a/tests/testthat/test-calibrator-time.R +++ b/tests/testthat/test-calibrator-time.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("times can be supplied and assume daily time-step (with warning)", { sir_obs = ("starter_models/shiver/data/hospitalizations_ontario.csv" |> system.file(package = "macpan2") diff --git a/tests/testthat/test-calibrator-traj.R b/tests/testthat/test-calibrator-traj.R index 8a31733d7..cb57dd7f7 100644 --- a/tests/testthat/test-calibrator-traj.R +++ b/tests/testthat/test-calibrator-traj.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("bad outputs give warnings", { sir = mp_tmb_library("starter_models", "sir", package = "macpan2") expect_warning( diff --git a/tests/testthat/test-calibrator-tv.R b/tests/testthat/test-calibrator-tv.R index 21d6c6ec5..ef1d522df 100644 --- a/tests/testthat/test-calibrator-tv.R +++ b/tests/testthat/test-calibrator-tv.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("rbf basis dimension is the same size as initial weights", { spec = ( mp_tmb_library("starter_models", "sir", package = "macpan2") diff --git a/tests/testthat/test-change-model.R b/tests/testthat/test-change-model.R index ec11d26cc..87b73eb0a 100644 --- a/tests/testthat/test-change-model.R +++ b/tests/testthat/test-change-model.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr) spec = mp_tmb_model_spec( during = list( k1_S ~ 1 ## rk4 coefficient name to test for name conflict avoidance diff --git a/tests/testthat/test-clamp.R b/tests/testthat/test-clamp.R index 81bab55ed..c9210c16a 100644 --- a/tests/testthat/test-clamp.R +++ b/tests/testthat/test-clamp.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("clamping function is modified squareplus", { x = seq(-1e-11, 1e-11, by = 1e-14) clamp_base_r = function(x, eps = 1e-12, limit = eps) { diff --git a/tests/testthat/test-convolution.R b/tests/testthat/test-convolution.R index 9288a0240..f190141ec 100644 --- a/tests/testthat/test-convolution.R +++ b/tests/testthat/test-convolution.R @@ -1,5 +1,4 @@ ## not done by any stretch -library(macpan2) r = simple_sims( list( diff --git a/tests/testthat/test-date-time-steps.R b/tests/testthat/test-date-time-steps.R index 89d248ecb..2fdad4827 100644 --- a/tests/testthat/test-date-time-steps.R +++ b/tests/testthat/test-date-time-steps.R @@ -1,5 +1,4 @@ test_that("dates can be turned into time-steps", { - library(lubridate) base_date = as.Date("2024-01-25") expect_error( macpan2:::Weekly(start = base_date - days(30), end = base_date) diff --git a/tests/testthat/test-dimnames-and-shape-change.R b/tests/testthat/test-dimnames-and-shape-change.R index fc1119184..87a9cebb4 100644 --- a/tests/testthat/test-dimnames-and-shape-change.R +++ b/tests/testthat/test-dimnames-and-shape-change.R @@ -1,6 +1,3 @@ -library(macpan2) -library(testthat) - test_that("matrices with dimnames can change size", { m = macpan2:::TMBModel( init_mats = macpan2:::MatsList(x = c(a = 0), .mats_to_save = "x", .mats_to_return = "x"), diff --git a/tests/testthat/test-distributions.R b/tests/testthat/test-distributions.R index cc9ee20c1..8bf897f35 100644 --- a/tests/testthat/test-distributions.R +++ b/tests/testthat/test-distributions.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2); library(deSolve) test_that("distributions give appropriate variable assumption warnings", { # At this time the only distribution with variable assumptions is the diff --git a/tests/testthat/test-dynamic-model.R b/tests/testthat/test-dynamic-model.R index 579520bec..64b198c6d 100644 --- a/tests/testthat/test-dynamic-model.R +++ b/tests/testthat/test-dynamic-model.R @@ -4,9 +4,6 @@ ## from the engine_agnostic_grammar vignette. ## ----setup, echo = FALSE, message = FALSE, warning = FALSE-------------------- -library(macpan2) -library(ggplot2) -library(dplyr) ## ----SIR-starter, echo = FALSE------------------------------------------------ diff --git a/tests/testthat/test-index.R b/tests/testthat/test-index.R index f3b0af48d..64421eb61 100644 --- a/tests/testthat/test-index.R +++ b/tests/testthat/test-index.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("functions that act on index objects behave", { # SEIR epi compartments with vaccination status epi_states = c("S","E","I","R") diff --git a/tests/testthat/test-insert-state-vars.R b/tests/testthat/test-insert-state-vars.R index ee32a933d..817002e4a 100644 --- a/tests/testthat/test-insert-state-vars.R +++ b/tests/testthat/test-insert-state-vars.R @@ -18,12 +18,14 @@ test_that("inserted expressions are able to refer to single states/flows", { , to = c(1, 2) ) ) - s = mp_simulator(spec - , time_steps = 50 - , outputs = c("S", "I") + s = (spec + |> mp_tmb_insert(at = Inf + , expressions = list(test_ratio ~ state[I] / state[S]) + ) + |> mp_simulator(time_steps = 50 + , outputs = c("S", "I", "test_ratio") + ) ) - s$add$matrices(test_ratio = empty_matrix, .mats_to_save = "test_ratio", .mats_to_return = "test_ratio") - s$insert$expressions(test_ratio ~ state[I] / state[S], .at = Inf, .phase = "during") v = mp_trajectory(s) expect_equal( v[v$row == "I", "value"] / v[v$row == "S", "value"], diff --git a/tests/testthat/test-list-utils.R b/tests/testthat/test-list-utils.R index f1fb461c0..9186ab6f0 100644 --- a/tests/testthat/test-list-utils.R +++ b/tests/testthat/test-list-utils.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("nlist and list are equivalent when explicit names are used", { expect_identical(nlist(a = 1), list(a = 1)) }) diff --git a/tests/testthat/test-model-modifications.R b/tests/testthat/test-model-modifications.R index 693a50978..f13fdb6ee 100644 --- a/tests/testthat/test-model-modifications.R +++ b/tests/testthat/test-model-modifications.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("expressions can be deleted from model spec",{ sir = mp_tmb_library("starter_models","sir", package = "macpan2") # remove I to R flow diff --git a/tests/testthat/test-model-structure.R b/tests/testthat/test-model-structure.R index d102dcdef..382d8023d 100644 --- a/tests/testthat/test-model-structure.R +++ b/tests/testthat/test-model-structure.R @@ -1,5 +1,3 @@ -library(macpan2) - Epi = mp_index(Epi = c("S", "I", "R")) Age = mp_index(Age = c("young", "old")) diff --git a/tests/testthat/test-mp-rk4.R b/tests/testthat/test-mp-rk4.R index d5aa3f199..519981a8f 100644 --- a/tests/testthat/test-mp-rk4.R +++ b/tests/testthat/test-mp-rk4.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("mp_rk4() does not repeat preceding random computations to flows", { seir = mp_tmb_library("starter_models","seir",package = "macpan2") diff --git a/tests/testthat/test-partitions.R b/tests/testthat/test-partitions.R index 024ca1720..8ece225cb 100644 --- a/tests/testthat/test-partitions.R +++ b/tests/testthat/test-partitions.R @@ -1,5 +1,3 @@ -library(oor) -library(macpan2) r = CSVReader(system.file("model_library", "sir_vax", "variables.csv", package = "macpan2")) x = macpan2:::Partition(r$read()) x$name() diff --git a/tests/testthat/test-rbind-time-lag.R b/tests/testthat/test-rbind-time-lag.R index 1f990aab7..099e3579a 100644 --- a/tests/testthat/test-rbind-time-lag.R +++ b/tests/testthat/test-rbind-time-lag.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) test_that("a selection of the iterations in the simulation history of a matrix that doesn't change shape can be rbinded at the end", { steps = 10 x = matrix(1:12, 4, 3) diff --git a/tests/testthat/test-scalar-2-vector.R b/tests/testthat/test-scalar-2-vector.R index 03acf0f48..59dbe2799 100644 --- a/tests/testthat/test-scalar-2-vector.R +++ b/tests/testthat/test-scalar-2-vector.R @@ -1,3 +1,3 @@ -#library(macpan2) + #box_model = Model(ModelFiles("inst/starter_models/seir_symp_vax")) #scalar2vector = Scalar2Vector(box_model) diff --git a/tests/testthat/test-state-flow-order.R b/tests/testthat/test-state-flow-order.R index 35c7910ce..fe39c579b 100644 --- a/tests/testthat/test-state-flow-order.R +++ b/tests/testthat/test-state-flow-order.R @@ -1,6 +1,5 @@ test_that("the order of state variables and flow variables is defined by settings not code", { skip("test depends on old 'Compartmental' model approach") - library(macpan2) sir = Compartmental(system.file("starter_models", "sir", package = "macpan2")) N = 100 state = c(S = N - 1, I = 1, R = 0) @@ -25,8 +24,6 @@ test_that("the order of state variables and flow variables is defined by setting test_that("the order of state variables and flow variables in settings does not impact simulations", { skip("test depends on old 'Compartmental' model approach") - library(macpan2) - library(dplyr) sir = Compartmental(system.file("starter_models", "sir", package = "macpan2")) sir_mixup = Compartmental(system.file("testing_models", "sir", package = "macpan2")) N = 100 diff --git a/tests/testthat/test-tmb-model-spec.R b/tests/testthat/test-tmb-model-spec.R index e88b9f5f5..bbf8bde2b 100644 --- a/tests/testthat/test-tmb-model-spec.R +++ b/tests/testthat/test-tmb-model-spec.R @@ -1,4 +1,3 @@ -library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) si = mp_tmb_model_spec( before = list( I ~ 1 diff --git a/tests/testthat/test-tmb.R b/tests/testthat/test-tmb.R index a2024925c..2e27d6ff7 100644 --- a/tests/testthat/test-tmb.R +++ b/tests/testthat/test-tmb.R @@ -1,8 +1,6 @@ - -library(macpan2) m = macpan2:::TMBModel( macpan2:::MatsList(x = 0, y = 0, .mats_to_save = c("x", "y"), .mats_to_return = c("x", "y")), - mp_tmb_expr_list(during = list(x ~ x + 1)), + mp_tmb_expr_list(during = list(x ~ x + 1)), macpan2:::OptParamsList(0, par_id = 0, mat = "x", row_id = 0, col_id = 0), macpan2:::OptParamsList(), macpan2:::ObjectiveFunction(~x), From 7398137558ea31dd830e68589a547219a717f9b3 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Wed, 30 Jul 2025 16:08:41 -0400 Subject: [PATCH 11/13] cleaned up too much --- tests/testthat.R | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/testthat.R b/tests/testthat.R index 22c989be3..23c164c2e 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -6,4 +6,5 @@ # * https://r-pkgs.org/tests.html # * https://testthat.r-lib.org/reference/test_package.html#special-files +library(testthat) test_check("macpan2") From a1159b1bf31b7cb5da566e4e01082ecc61edc252 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Wed, 30 Jul 2025 16:32:55 -0400 Subject: [PATCH 12/13] tweaks --- vignettes/kronecker.Rmd | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/vignettes/kronecker.Rmd b/vignettes/kronecker.Rmd index f30765eb7..edcef22cf 100644 --- a/vignettes/kronecker.Rmd +++ b/vignettes/kronecker.Rmd @@ -95,11 +95,18 @@ The `ones_col()`, `ones_row()`, and `identity()` functions make it easier to bui ## Initial State Vectors -We now initialize the state vectors `S` (susceptible) and `I` (infectious) over the stratification dimensions. +We now use Kronecker products to initialize the state vectors `S` (susceptible) and `I` (infectious) over the stratification dimensions. ```{r initial-states} -S = c(young = 75, old = 25) %x% c(unvax = 1, vax = 0) -I = c(young = 1, old = 0) %x% c(unvax = 1, vax = 0) %x% c(mild = 1, severe = 0) +S = ( + c(young = 75, old = 25) + %x% c(unvax = 1, vax = 0) +) +I = ( + c(young = 1, old = 0) + %x% c(unvax = 1, vax = 0) + %x% c(mild = 1, severe = 0) +) print(S) print(I) ``` @@ -157,7 +164,11 @@ Note that `N` repeats the total population size within each age group to match t We construct the full per-capita transmission matrix by combining component matrices using Kronecker products. ```{r per-capita-transmission-matrix} -B_age = beta * contact_matrix * (susceptibility_age %x% infectivity_age) +B_age = ( + beta + * contact_matrix + * (susceptibility_age %x% infectivity_age) +) B_vax = (1 - VE) %x% infectivity_vax B_symp = infectivity_symp B = B_age %x% B_vax %x% B_symp @@ -229,7 +240,11 @@ We formalize the model specification with `mp_tmb_model_spec`, defining the deri ```{r model-spec} spec = mp_tmb_model_spec( before = list( - B_age ~ beta * contact_matrix * (susceptibility_age %x% infectivity_age) + B_age ~ ( + beta + * contact_matrix + * (susceptibility_age %x% infectivity_age) + ) , B_vax ~ (1 - VE) %x% infectivity_vax , B_symp ~ infectivity_symp , B ~ B_age %x% B_vax %x% B_symp From 94c0c173b9e0996b94be8e2040a0eded8dff7a0d Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Wed, 30 Jul 2025 16:34:26 -0400 Subject: [PATCH 13/13] [skip ci] tweaks --- vignettes/matrix_multiplication.Rmd | 155 +--------------------------- 1 file changed, 1 insertion(+), 154 deletions(-) diff --git a/vignettes/matrix_multiplication.Rmd b/vignettes/matrix_multiplication.Rmd index 1765534cc..c26a60059 100644 --- a/vignettes/matrix_multiplication.Rmd +++ b/vignettes/matrix_multiplication.Rmd @@ -54,161 +54,8 @@ A %x% B macpan2::engine_eval(~A %x% B, A = A, B = B) ``` -Kronecker products are useful when replicating structure across groups or compartments, as the follow section illustrates. +Kronecker products are useful when replicating structure across groups or compartments. -### Example -- Structured SI Model with Kronecker products - -We describe an SI model structured by age group, vaccination status, and symptom status, entirely in terms of matrices and vectors so that Kronecker products can be used to consistently replicate structure across these groups. - -#### Step 1 — Stratifying State Variables - -We divide a population into: - -- \( n_{\text{age}} \) age groups (indexed by \( a = 1, \dots, n_{\text{age}} \)) -- \( n_{\text{vax}} \) vaccination categories (indexed by \( v = 1, \dots, n_{\text{vax}} \)) -- \( n_{\text{symp}} \) symptom categories (indexed by \( s = 1, \dots, n_{\text{symp}} \)) - -All population quantities are represented as column vectors, stacked in this order: - -- Age group (leading index) — varies slowest -- Vaccination status — varies next -- Symptom status — varies fastest (only for infectious individuals) - -This convention ensures that aggregation and mixing operations defined by Kronecker products naturally correspond to the structure of the population. - -Susceptible vector \( S \) with \( n_{\text{age}} = 2 \), \( n_{\text{vax}} = 2 \): -\[ -S = \begin{bmatrix} -S_{1,1} \\ -S_{1,2} \\ -S_{2,1} \\ -S_{2,2} -\end{bmatrix} -\] - -Infectious vector \( I \) with \( n_{\text{age}} = 2 \), \( n_{\text{vax}} = 2 \), \( n_{\text{symp}} = 3 \): -\[ -I = \begin{bmatrix} -I_{1,1,1} \\ -I_{1,1,2} \\ -I_{1,1,3} \\ -I_{1,2,1} \\ -I_{1,2,2} \\ -I_{1,2,3} \\ -I_{2,1,1} \\ -I_{2,1,2} \\ -I_{2,1,3} \\ -I_{2,2,1} \\ -I_{2,2,2} \\ -I_{2,2,3} -\end{bmatrix} -\] - -#### Step 2 — Total Population by Age Group - -We construct summation matrices to aggregate the population counts into totals by age group. Each state variable requires its own summation matrix because they are structured differently: - -- **\( S \)** is structured by age and vaccination status. -- **\( I \)** is structured by age, vaccination status, and symptom status. - -To define the summation matrices (and expansion matrices below), we use: - -- \( \mathbb{I}_n \) — the \( n \times n \) identity matrix (ones on the diagonal, zeros elsewhere). -- \( \mathbf{1}_{1 \times n} \) — a \( 1 \times n \) row vector of ones. -- \( \mathbf{1}_{n \times 1} \) — an \( n \times 1 \) column vector of ones. - - -These matrices are combined using Kronecker products, ensuring summation within the correct strata. - -- For **susceptibles** \( S \), we sum over vaccination status while preserving age: -\[ -M_S = \mathbb{I}_{n_{\text{age}}} \otimes \mathbf{1}_{1 \times n_{\text{vax}}} -\] - -- For **infectious individuals** \( I \), we sum over both vaccination status and symptom status while preserving age: -\[ -M_I = \mathbb{I}_{n_{\text{age}}} \otimes \mathbf{1}_{1 \times (n_{\text{vax}} \cdot n_{\text{symp}})} -\] - -We then compute the total population per age group by applying the summation matrices: -\[ -N_{\text{age}} = M_S \cdot S + M_I \cdot I -\] - -This yields a column vector of length \( n_{\text{age}} \), giving the total population in each age group. - -#### Step 3 — Expanded Total Population Vector - -We compute a population vector \( N \) with the same structure and length as the infectious vector \( I \), so that it can be used directly in the force of infection expression below. - -We obtain \( N \) by expanding the vector of total population by age group, \( N_{\text{age}} \), using the Kronecker product with a vector of ones over the vaccination and symptom strata: -\[ -N = N_{\text{age}} \otimes \mathbf{1}_{n_{\text{vax}} \cdot n_{\text{symp}} \times 1} -\] - -Each entry of \( N \) gives the total population of the corresponding age group, repeated across all vaccination and symptom categories. - -#### Step 4 — WAIFW Matrix Construction - -The above notation allows us to very nicely partition the effects of each stratum on the *Who Acquires Infection From Whom* (WAIFW) matrix for our structured SI model, using Kronecker products of component WAIFW matrices: - -\[ -W = C_{\text{age}} \otimes C_{\text{vax}} \otimes C_{\text{symp}} -\] - -This representation reflects the structure of our state variables: susceptible individuals are stratified by age and vaccination status, whereas infectious individuals are stratified by age, vaccination status, and symptom status. The dimensions and interpretations of each component matrix follow from this: - -- \( C_{\text{age}} \in \mathbb{R}^{n_{\text{age}} \times n_{\text{age}}} \): - A square matrix representing transmission potential between susceptible individuals in each age group (rows) and infectious individuals in each age group (columns). This is typically a contact matrix, where entry \((i, j)\) encodes the probability that an individual in age group \(j\) will contact an individual in age group \(i\). - -- \( C_{\text{vax}} \in \mathbb{R}^{n_{\text{vax}} \times n_{\text{vax}}} \): - This matrix modulates the impact of vaccination status on transmission. Each entry \((i, j)\) quantifies the transmission potential from an infectious individual with vaccination status \(j\) to a susceptible individual with status \(i\). Often, this captures relative reductions in susceptibility or infectiousness due to vaccination. - -- \( C_{\text{symp}} \in \mathbb{R}^{1 \times n_{\text{symp}}} \): - Unlike age and vaccination status, susceptibles do not have a symptom status—one might say their symptom status is implicitly "none". Consequently, \( C_{\text{symp}} \) is a 1-by-\(n_{\text{symp}}\) row vector, where each entry adjusts the infectiousness of individuals depending on their symptom status (e.g., asymptomatic, presymptomatic, mildly symptomatic). This can encode, for example, that asymptomatic individuals are half as infectious as symptomatic ones. - -The full WAIFW matrix \( W \) has dimension \((n_{\text{age}} n_{\text{vax}}) \times (n_{\text{age}} n_{\text{vax}} n_{\text{symp}})\) and defines how each susceptible subgroup (by age and vaccination) acquires infection from each infectious subgroup (by age, vaccination, and symptoms), assuming homogeneous mixing within each cell of this stratified structure. - -#### Step 5 — Force of Infection Vector - -We can now write the force of infection vector, with elements corresponding to the elements of the $S$ vector. - -\[ -\Lambda = \beta \cdot W \cdot I \cdot \text{diag}^{-1}(N) -\] - -#### Step 6 — Inflows and Outflows with Allocation Matrices - -In the structured SI model, the force of infection determines the rate at which individuals leave the susceptible compartment. However, because the infectious compartment \( I \) includes an additional stratification (symptom status), the force of infection does not specify how new infections are distributed among the strata of \( I \). We introduce allocation matrices to define this distribution explicitly and ensure consistent tracking of flows. - -We define: - -- \( \text{Out}_{S} \) — the outflow vector from the susceptible compartment, with the same structure as \( S \). -- \( \text{In}_{I} \) — the inflow vector into the infectious compartment, with the same structure as \( I \). - -The outflow from \( S \) due to infection is computed elementwise: -\[ -\text{Out}_{S} = \Lambda \circ S -\] -where \( \circ \) denotes elementwise multiplication. - -To compute the inflow vector we introduce the allocation matrix \( A_{S \to I} \), a \( (n_{\text{age}} \cdot n_{\text{vax}} \cdot n_{\text{symp}}) \times (n_{\text{age}} \cdot n_{\text{vax}}) \) matrix such that: - -- All elements of $A_{S \to I}$ are $\geq 0$. -- All columns sum to one. - -The inflow to \( I \) is then given by: -\[ -\text{In}_{I} = A_{S \to I} \cdot \text{Out}_{S} -\] - -The required properties of \( A_{S \to I} \) ensure that the outflow from each susceptible stratum is completely allocated among the infectious strata. Rows may sum to any non-negative value, including zero, which means that some infectious strata may receive no inflow directly from \( S \). - -#### Structured SI Model Specification - -```{r} - -``` ## Elementwise Multiplication with Recycling (`*`)