filling in an element of a matrix?

[bbolker]

I have an epidemiological model where I pre-compute the per capita flows from compartment i to j and store them in a matrix M: this way, the only term that needs to be updated over the course of the epidemic is the S -> E per capita flows: all the rest stay constant (eventually could relax this, e.g. due to time-varying contact-tracing effort which would move infected people to quarantine faster). This formulation also makes it easier to go back and forth between deterministic and stochastic formulations (e.g. the "flow matrix" M is set up nicely for pomp's reulermultinom function for stochastically moving discrete individuals on the basis of multiple competing stochastic processes).

I may be completely missing something, but I'm having trouble compiling the gradient function when I try to set a particular component of M. Is there a way to do this? (As a workaround, is there a way to generate the C code with the M[2,1] <- ... line commented out and then hack the code before compiling?)

Happy to submit this to some more appropriate forum if appropriate ...

remotes::install_github("bbolker/McMasterPandemic")
library(McMasterPandemic)
## setup, via McMasterPandemic
pars <- read_params("ICU1.csv")
state <- make_state(param=pars)
M <- make_ratemat(state, pars)
odin_gradfun <- odin::odin({
    s_init[] <- user()
    initial(s) <- s_init
    M[,] <- user()
    beta0 <- user()
    Ca <- user()
    Cp <- user()
    iso_m <- user()
    Cm <- user()
    iso_s <- user()
    Cs <- user()
    N <- user()
    M[2,1] <- beta0/N*(Ca*s[3]+Cp*s[4]+(1-iso_m)*Cm*s[5]+(1-iso_s)*Cs*s[6])
    flows <- M[i,j]*s[j]
    deriv(s) <- sum(M[i,])-sum(M[,j])
})

[richfitz]

Hi @bbolker - thanks for the message, and sorry for missing it (I declared github notification bankruptcy a long time ago but you can ping me with a notification or on the rOpenSci slack if you run into trouble).

I think you just need to disambiguate this a bit for odin. Here's a slightly smaller case, as the one above won't compile without the dim() calls. This fails with the same error as yours:

odin::odin({
  M[, ] <- user()
  M[1,2] <- 1
  dim(M) <- c(4, 4)
  initial(x) <- 0
  deriv(x) <- sum(M)
})
#> Error: user() may only be used on a single-line array
#>   M[, ] <- user() # (line 1)
#>   M[1, 2] <- 1 # (line 2)

however, this one will compile

gen <- odin::odin({
  m0 <- user()
  M[, ] <- m0
  M[1,2] <- 1
  dim(M) <- c(4, 4)
  initial(x) <- 0
  deriv(x) <- sum(M)
})

if your initial value for M is a matrix (sounds like it is) then that would be:

gen <- odin::odin({
  M0[, ] <- user()
  M[, ] <- M0[i, j]
  M[1,2] <- 1
  dim(M) <- c(4, 4)
  dim(M0) <- c(4, 4)
  initial(x) <- 0
  deriv(x) <- sum(M)
})

Creating the model and showing its internal variables:

mod <- gen(M0 = matrix(-seq_len(16), 4, 4))
mod$contents()
#> > mod <- gen(M0 = matrix(-seq_len(16), 4, 4))
#> > # This shows both M0 and M:
#> > mod$contents()
#> $dim_M
#> [1] 16
#> 
#> $dim_M_1
#> [1] 4
#> 
#> $dim_M_2
#> [1] 4
#> 
#> $dim_M0
#> [1] 16
#> 
#> $dim_M0_1
#> [1] 4
#> 
#> $dim_M0_2
#> [1] 4
#> 
#> $initial_x
#> [1] 0
#> 
#> $M
#>      [,1] [,2] [,3] [,4]
#> [1,]   -1    1   -9  -13
#> [2,]   -2   -6  -10  -14
#> [3,]   -3   -7  -11  -15
#> [4,]   -4   -8  -12  -16
#> 
#> $M0
#>      [,1] [,2] [,3] [,4]
#> [1,]   -1   -5   -9  -13
#> [2,]   -2   -6  -10  -14
#> [3,]   -3   -7  -11  -15
#> [4,]   -4   -8  -12  -16

here, because M can be computed without reference to anything time-varying M has been set before any integration happens. But if your value you insert into M does depend on time there will be a copy each time unfortunately:

gen <- odin::odin({
  M0[, ] <- user()
  M[, ] <- M0[i, j]
  M[1,2] <- 1 + (t * 0)
  dim(M) <- c(4, 4)
  dim(M0) <- c(4, 4)
  initial(x) <- 0
  deriv(x) <- sum(M)
})
mod <- gen(M0 = matrix(-seq_len(16), 4, 4))
mod$contents()$M
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    0    0    0
#> [2,]    0    0    0    0
#> [3,]    0    0    0    0
#> [4,]    0    0    0    0
mod$deriv(0, 0)
mod$contents()$M
#>      [,1] [,2] [,3] [,4]
#> [1,]   -1    1   -9  -13
#> [2,]   -2   -6  -10  -14
#> [3,]   -3   -7  -11  -15
#> [4,]   -4   -8  -12  -16

[richfitz]

I've opened a new ticket (#195) which would allow this to behave more efficiently as currently the generated code looks like:

  for (int i = 1; i <= internal->dim_M_1; ++i) {
    for (int j = 1; j <= internal->dim_M_2; ++j) {
      internal->M[i - 1 + internal->dim_M_1 * (j - 1)] = internal->M0[internal->dim_M0_1 * (j - 1) + i - 1];
    }
  }
  {
     int i = 1;
     int j = 2;
     internal->M[i - 1 + internal->dim_M_1 * (j - 1)] = 1 + (t * 0);
  }

with that whole top section being a waste of time. I doubt this will be a huge time sink though?