sparse matrix multiplication and associated fixes

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

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)

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

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

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, ...)"
 )

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

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))