diff --git a/DESCRIPTION b/DESCRIPTION index c81decd0c..ed889704c 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"), @@ -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) @@ -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 cd0bcc9cf..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) @@ -363,6 +365,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/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 ac5363138..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 @@ -220,13 +275,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 +303,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 @@ -650,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 @@ -1057,6 +1093,8 @@ #' @aliases invlogit #' @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 f655bb28c..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)" @@ -60,6 +60,8 @@ valid_func_sigs = c( , "fwrap: invlogit(x)" , "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/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/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..7ffff5911 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() @@ -92,11 +98,14 @@ 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) + 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) + ) |> unique() + c(implied_integers[integers_we_need], self$integers) } self$empty_matrices = function() { 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/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/R/tmb_model_editors.R b/R/tmb_model_editors.R index 41cf1fb3f..1d0d39284 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 @@ -307,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 @@ -322,13 +328,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 +365,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/README.Rmd b/README.Rmd index d3258c7d7..80b359fb9 100644 --- a/README.Rmd +++ b/README.Rmd @@ -17,6 +17,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) diff --git a/_pkgdown.yml b/_pkgdown.yml index ec85a2454..38478c6ee 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 @@ -261,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/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/engine_functions.Rd b/man/engine_functions.Rd index c54633a84..76d0775ec 100644 --- a/man/engine_functions.Rd +++ b/man/engine_functions.Rd @@ -62,6 +62,8 @@ \alias{invlogit} \alias{logit} \alias{cumsum} +\alias{sparse_mat_mult} +\alias{divide_safe} \alias{assign} \alias{unpack} \title{Functions Available in the Simulation Engine} @@ -190,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)} } } @@ -218,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}{ @@ -315,15 +380,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 +414,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}{ @@ -814,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/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/man/mp_tmb_insert_glm_timevar.Rd b/man/mp_tmb_insert_glm_timevar.Rd index 93e1954c2..10a59d08a 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}}.} @@ -56,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/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/misc/dev/dev.cpp b/misc/dev/dev.cpp index 179d1d4ae..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) @@ -136,8 +136,10 @@ 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_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 @@ -216,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 @@ -530,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 } } @@ -876,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) @@ -996,13 +1009,18 @@ 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; - 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; @@ -1158,13 +1176,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 +1194,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; } } @@ -1378,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 // #' @@ -1394,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); @@ -1416,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 @@ -1457,7 +1577,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 +1602,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 +1668,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 +1677,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 +1705,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 +1724,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 +1874,6 @@ class ExprEvaluator } else { MP2_ERR(MP2_CBIND, "Inconsistent size in cbind function", MP2_CBIND); - return m; } } return m; @@ -1711,26 +1887,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 +1929,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 +2072,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 +2086,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 +2095,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 +2104,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 +2223,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 +2246,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 +2313,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 +2344,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 +2357,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 +2406,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 +2521,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 +2531,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 +2544,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; } // off is the 'initial offset', which lets the user @@ -2452,12 +2602,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); @@ -2531,7 +2679,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 @@ -2574,85 +2721,10 @@ 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 - // #' - // #' 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 @@ -2711,7 +2783,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(); @@ -2720,7 +2791,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++) { @@ -2733,7 +2803,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(); @@ -2743,7 +2812,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 @@ -2764,7 +2832,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(); @@ -2774,7 +2841,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++) { @@ -2787,7 +2853,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(); @@ -2797,7 +2862,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++) { @@ -2880,10 +2944,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(); @@ -2891,10 +2953,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++) @@ -2913,10 +2973,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(); @@ -2924,10 +2982,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++) @@ -2943,7 +2999,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(); @@ -2952,7 +3007,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++) { @@ -2973,7 +3027,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; @@ -3043,7 +3096,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 @@ -3065,7 +3117,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 @@ -3198,22 +3249,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 @@ -3223,10 +3268,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(); @@ -3236,7 +3279,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++) { @@ -3245,7 +3287,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); } @@ -3304,19 +3345,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; @@ -3330,7 +3367,6 @@ class ExprEvaluator default: MP2_ERR(255, "invalid operator in arithmetic expression", -99); - return m; } } // switch (table_n[row]) }; 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 54ae64d59..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) @@ -137,8 +137,10 @@ 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_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 @@ -217,7 +219,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 +543,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 +884,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,13 +1010,18 @@ 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; - 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; @@ -1159,13 +1177,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 +1195,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; } } @@ -1379,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 // #' @@ -1395,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); @@ -1417,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 @@ -1458,7 +1578,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 +1603,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 +1669,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 +1678,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 +1706,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 +1725,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 +1875,6 @@ class ExprEvaluator } else { MP2_ERR(MP2_CBIND, "Inconsistent size in cbind function", MP2_CBIND); - return m; } } return m; @@ -1712,26 +1888,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 +1930,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 +2073,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 +2087,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 +2096,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 +2105,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 +2224,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 +2247,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 +2314,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 +2345,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 +2358,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 +2407,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 +2522,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 +2532,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 +2545,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; } // off is the 'initial offset', which lets the user @@ -2453,12 +2603,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); @@ -2532,7 +2680,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 @@ -2575,85 +2722,10 @@ 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 - // #' - // #' 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 @@ -2712,7 +2784,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(); @@ -2721,7 +2792,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++) { @@ -2734,7 +2804,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(); @@ -2744,7 +2813,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 @@ -2765,7 +2833,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(); @@ -2775,7 +2842,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++) { @@ -2788,7 +2854,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(); @@ -2798,7 +2863,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++) { @@ -2881,10 +2945,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(); @@ -2892,10 +2954,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++) @@ -2914,10 +2974,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(); @@ -2925,10 +2983,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++) @@ -2944,7 +3000,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(); @@ -2953,7 +3008,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++) { @@ -2974,7 +3028,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; @@ -3044,7 +3097,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 @@ -3066,7 +3118,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 @@ -3199,22 +3250,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 @@ -3224,10 +3269,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(); @@ -3237,7 +3280,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++) { @@ -3246,7 +3288,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); } @@ -3305,19 +3346,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; @@ -3331,7 +3368,6 @@ class ExprEvaluator default: MP2_ERR(255, "invalid operator in arithmetic expression", -99); - return m; } } // switch (table_n[row]) }; diff --git a/tests/testthat.R b/tests/testthat.R index 2098c2921..23c164c2e 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -7,10 +7,4 @@ # * 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-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) }) + + diff --git a/tests/testthat/test-clamp.R b/tests/testthat/test-clamp.R index 1b90b9f5e..c9210c16a 100644 --- a/tests/testthat/test-clamp.R +++ b/tests/testthat/test-clamp.R @@ -1,13 +1,12 @@ -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) { 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 +15,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/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-expr-list.R b/tests/testthat/test-expr-list.R index e4a1c4c22..63d19b3fb 100644 --- a/tests/testthat/test-expr-list.R +++ b/tests/testthat/test-expr-list.R @@ -16,10 +16,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-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') + # ) +}) 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-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" + ) +}) 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), diff --git a/vignettes/kronecker.Rmd b/vignettes/kronecker.Rmd new file mode 100644 index 000000000..edcef22cf --- /dev/null +++ b/vignettes/kronecker.Rmd @@ -0,0 +1,434 @@ +--- +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 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) +) +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 new file mode 100644 index 000000000..c26a60059 --- /dev/null +++ b/vignettes/matrix_multiplication.Rmd @@ -0,0 +1,239 @@ +--- +title: "Matrix Operations in macpan2 Models" +output: + rmarkdown::html_vignette: + toc: true +vignette: > + %\VignetteIndexEntry{Matrix Operations in macpan2 Models} + %\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} +library(macpan2) +library(dplyr) +library(ggplot2) +library(splines) +theme_bw = function() ggplot2::theme_bw(base_size = 18) +``` + +## Overview + +Matrix operations are central to `macpan2` models. The engine supports four ways to multiply two matrices together: + +- Standard matrix multiplication (`%*%`) +- Kronecker products (`%x%`) +- Elementwise multiplication with recycling (`*`) +- Left sparse matrix multiplication (`sparse_mat_mult()`) + +This vignette introduces each operation and illustrates typical use cases. + +## 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) +``` + +Kronecker products are useful when replicating structure across groups or compartments. + + +## 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) + +sparse <- macpan2:::sparse_matrix_notation(M, zero_based = TRUE) +print(sparse) + +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)) +``` + +### Small Example of Sparse Matrix Multiplication + + +```{r} +# 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 +print(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 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. + +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) + +# 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 +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 +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 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.