mrc-ide / odin

ᚩ A DSL for describing and solving differential equations in R
https://mrc-ide.github.io/odin
Other
103 stars 13 forks source link

Handing products for multiplicative interventions #316

Open adamkucharski opened 1 month ago

adamkucharski commented 1 month ago

Thanks for all your great work on this package. We’ve been looking at combining the {epidemics} package interface (which makes contact matrices, interventions etc. easy for user to define, and can be reused in other packages) with {odin} model implementionss (which are easier to edit and can be passed to {mcstate} etc). One thing that came up is how to implement multiple interventions that multiplicatively reduce the contact rates (e.g. school closures on top of workplace closures).

It's my understanding that prod isn't a valid operation in an odin model, so we've been using log transforms as a work around (full draft vignette here):

  # Contact matrix with time-dependent interventions
  contact_reduction[,] <- if (t > intervention_start[i] && t < intervention_end[i]) (1.0 - intervention_effect[i,j]) else 1
  contact_reduction_log[,] <- log(contact_reduction[i,j] + 1e-10) # Avoid log(0)
  contact_reduction_total_log[] <- sum(contact_reduction_log[,i]) 
  contact_reduction_total[] <- exp(contact_reduction_total_log[i]) 

  # Specify how transmission varies over time
  lambda_prod[, ] <- C[i, j] * I[j] * beta * contact_reduction_total[j] * contact_reduction_total[i]
  lambda[] <- sum(lambda_prod[i, ])

Wanted to check that we weren't missing something more efficient, and to highlight a use case for value of having a prod operation if this is potentially in the pipeline.

richfitz commented 1 month ago

I have no objections to supporting prod() by extension of sum(), and you are right that this is not easily supported at present. When I saw the title I thought you were looking for matrix multiplication (%*%) which is also something we'd like to support but has proven hard to reason about (#134, #291, and other issues).

Progress is likely to be slow at the moment as we're taking the opportunity to dig the stack out of some pandemic-induced technical debt, but the good news is I'll get this implemented once going through and doing sum() in the new version.