# Product Models

This is a working document for defining product models, which is a framework that we are considering for macpan.

1. Worden & Porco, 2017
* Multy-layer compartmental models can be constructed by taking a modified Cartesian product of single-layer models
* There are 4 different products discussed, they are the Cartesian product, the modified Cartesian product, the strong product, and the modified strong product
2. Friston Et Al, 2020
* Transition rates between compartments are specified via conditional probabilities
* Discribes a single compartmental model as opposed to a structure for creating them.

## Product Definitions
Suppose there are 2 single-layer models A and B. Let A have $m$ compartments $a_1, a_2, ..., a_m$ and let B have $n$ compartments $b1, b2, ..., b_n$. Thyen, in all cases, the compartments of the product model will be the pairs $(a_i, b_j)$ where $i \in {1, ..., m}$ and $j\in{1, ..., n}$. The differences between the products is a matter of how the flows are defined.

#### Two single-layer compartmental models A(sir) and B(uvp)

 

### Unmodified products
#### Cartesian Product
When taking the Cartesian product of A and B (which we denote here by $A*B$)the flows in the product model are defined acording to the following two rules:
* If A has a flow from $a_k$ to $a_l$ then for every $j\in{1, ..., n}$ $A*B$ has a flow from $(a_k, b_j)$ to $(a_l, b_j)$. The symetric rule applies for the flows of B.
* if A has a flow from $a_k$ to $a_l$ and the flow rate is given by an expression $f(a_{x1}, a_{x2}, ..., a_{xs})$ where the $a_{xi}$ denote any compartment of A that is NOT $a_k$. Then for every $j\in{1, ..., n}$ the flow in $A*B$ from $(a_k, b_j)$ to $(a_l, b_j)$ is given by $f((a_{x1}, b_j), ..., (a_{xs}, b_j))$.

#### The unmodified Cartesian product of A and B

#### Strong Product
The Strong Product between A and B (which we denote by $A**B$) is constructed in the same way as the Cartesian product but with the additional two rules:
* If A has a flow from $a_k$ to $a_l$ and B has a flow from $b_u$ to $b_v$ then $A**B$ has a flow from $(a_k, b_u)$ to $(a_l, b_v)$. Note that this is in addition to the flows inherited from the Cartesian product.
* Question: How are the 'diagonal' flow rates defined?
* The Worden paper doesn't appear to address this directly. In their explanatory example the use a symmetry where the flow from $a_k$ to $a_l$ is given by $f(a_k, a_l)$ and the flow from $b_u$ to $b_v$ is $f(b_u, b_v)$. They then show the flow from $(a_k, b_u)$ to $(a_l, b_v)$ as $f((a_k, b_u), (a_l, b_v))$.
* Idea: If the flow from $a_k$ to $a_l$ is given by a function $f(a_k, a_l)$ and the flow from $b_u$ to $b_v$ is given by a function $g(b_u, b_v)$ can the flow from $(a_k, b_u)$ to $(a_l, b_v)$ be given by $fg((a_k, b_u), (a_l, b_v))$? At least as a default/preliminary definition.

#### The unmodified strong product of A and B

### Modified Products
The modified versions of the Cartesian and Strong products will produce the same compartments and flows as their unmodified analogs. The only difference will be in the flow rates. In particular if the unmodified product model has a flow from $(a_k, b_i)$ to $(a_l, b_j)$ with the flow rate defined as $f((a_{x1}, b_{y1}), ..., (a_{xs}, b_{yt}))$ then the modified flow rate will be given by $\sum\limits_{x1}...\sum\limits_{xs}f((a_{x1}, b_{y1}), ..., (a_{xs}, b_{yt}))$. That is to say if the unmodified flow rate depends on the population of a single compartment $(a_1, b_y)$ then the modified flow rate will depend on the population of all compartments $(a_x, b_y)$ where $x\in{1, ..., m}$.

The above definition is an attempt to allow for the flow rate to depend on the population of an arbitrary number of compartments. For the purpose of explanation its propably easier to consider the most likely case which is that the flow rate depends only on the population of the destination compartment. In that case if the unmodified product model has a flow rate $f((a_l, b_j))$ then the modified flow rate will be $\sum\limits_zf((a_z, b_j))$ where $a_z$ ranges over all compartments in model A.

#### The modified Cartesian product of A and B

#### The modified strong product of A and B

## Tentative framework for macpan

* Define $n$ single layer compartmental models
* Lable them $m_1$ through $m_n$
* Each should be able to run on its own
* Select from our models those that are "epi like"
* What does it mean for a model to be "epi-like"? Obviously anything that models a contagion but is there anything else?
* label them $m_{i1}$ through $m_{ik}$
* label the non-epi-like models $m_{j1}$ through $m_{jl}$
* Note $k+l=n$
* Create a new, multi-layer model $M_{epi} = m_{i1}**m_{i2}**...**m_{ik}$
* Using the modified strong product
* Call this model the Epi model
* Does order of operations matter here?
* Create our final model $M = M_{epi}*m_{j1}*m_{j2}*...*m_{jl}$
* Using the modified Cartesian product
* Does order of operations matter here?
* Provide "sculpting tools" so user can introduce irregularities

## Types of Rates in Models with Structure

flows in an n-d sub-model have rates that get replicated across m_1 other dimensions and rates that vary across m_2 other dimensions, such that the total dimension of the product model is N = n + m_1 + m_2

the replication case is easy in that if you specify one rate for every flow in the n-d sub-model, then you are done specifying for the full model

1. flows in a n-d sub-model with rates that get replicated across every other dimension
2. flows in a nb-d sub-model with rates that vary across the other dimensions ...
- where the variation needs to be exhaustively described
- according to some specific pattern
-

### Proposal

1. **States of Models are Products of Sub-Models --** Every model has `n` sub-models that divide the population into groups (e.g. SEIR, unvax/vax/protected, young/old), and the state of each individual is described by their state in each of these sub-models (e.g. S_unvax_young)
* Question: How do we order the naming convention? e.g. why S_unvax_young instead of unvax_S_young? Do we allow any order? user specified order?
3. **Specify Transitions within a Focal Sub-Model --** At least when initially setting up the model, every from-to pair should specify a single transition within a single focal sub-model (e.g. `S->E`, `unvax->vax`), with the default assumption that this transition occurs at every level of the other non-focal sub-models (e.g. `S->E` implies that `S_unvax_young->E_unvax_young`, `S_vax_young->E_vax_young`, ...)
* Note that we should consider allowing the user to select only certain states in the non-focal sub-models, which will be helpful for the case of vaccination that only applies to certain epidemiological categories
* On the other hand, we could handle this in #7 below with editing
4. **Sub-Models have Parameter Arrays --** Every sub-model has its own list of `n-1`-dimensional parameter arrays, where these dimensions allow for (but do not require) variation across the states of the `n-1` non-focal sub-models
5. **Sub-Models can have Structure Arrays --** Sub-models can optionally have a set of `n`-dimensional arrays that group parameters and/or state variables together, such that the focal dimension of these arrays identifies particular aspects of the focal model (e.g. one might define an `I` array describing a vector of different types of infected individuals for each combination of the non-focal model `I_unvax_young = [Ia_unvax_young, Ip_unvax_young, Im_unvax_young, Is_unvax_young]`)
6. **Specify Rates as Scalars --** Every transition-rate expression should be written to evaluate to a scalar, with the assumption that these scalar-valued rates will vary across the states of the non-focal sub-models in response to variation in parameters and state variables described in #3 above (e.g. the S->E rate is different for different vaccination states because vaccine efficacy varies among these states)
7. **Elements of Structure Vectors are (Implicitly or Explicitly?) Summed --** To ensure rate expressions evaluate to scalars, the elements of structure vectors (or combinations of structure vectors) are summed when they occur in transition-rate expressions (e.g. continuing the example from #4 above, `S->E: (1-VE) * beta0 * (1/N) * I` implies that `I` in this expression corresponds to the sum of the elements in the `I` vector defined in #4)
* Note that if we decide to go with this implicit version it seems (at least to me) to be connected with the Einstein summation convention for physical tensors
8. **Edit Rates --** If you really need to, you can add, delete, or overwrite individual transition-rates within the full model

### Nice things that we may not be able to have all at the same time

1. Every rate expression evaluates to a scalar
2. Freely allowed to make whatever matrix/vector you want (i.e. the struc stuff)
3. Ability to mix 1 and 2 -- sounds nice because it blends simplicity of 1 with the flexibility of 2, but is it actually any better?
4. Grouping values that work for both states and parameters (e.g. is it `I` that gets grouped by `a`, `p`, `m`, and `s` or is it `Ia`, `Ip`, `Im`, `Is`?)
5.

### Vaccination Example

#### Repeated rates

Here are some examples where the same rate gets repeated for all states of the non-focal model.

Ia->R ~ gamma_a
Is->H ~ (1 - nonhosp_mort) * gamma_s * phi1

In both of these cases the non-focal model is the vaccination model and the focal model is the SEIR+ model.

add_rate(
    from = "Is", 
    to = "H", 
    rate = ~ (1 - nonhosp_mort) * gamma_s * phi1
)

Here is another case where the vaccination model is focal.

unvax->vaxdose1 
    ~ vax_prop_first_dose 
    * vax_doses_per_day 
    * (1 / asymp_unvax_N))

One challenge here is that this particular case does not require a transition for every state in the non-focal sub-model.

add_rate(
    from = "unvax",
    to = "vaxdose1",
    rate = ~ vax_prop_first_dose 
           * vax_doses_per_day 
           * (1 / asymp_unvax_N),
    filter = list(
        epi = c("S", "E", "Ia", "Ip", "R")
    )
)

#### Varying rates across one dimension

E->Ia ~ alpha * sigma
E->Ip ~ (1 - alpha) * sigma

#### Varying rates across one dimension with rates depending on inner products

What Worden would do (we think):

S->E 
    ~ beta0 
    * (1/N) 
    * (Ia * Ca + Ip * Cp + Im * Cm + Is * Cs)

Here the assumption is that `beta0` varies with vaccination status, and in this way encorporates vaccine efficacy.

What we would do without model structure arrays (see above):

S->E 
    ~ beta0 
    * (1 - VE)
    * (1/N) 
    * (Ia * Ca + Ip * Cp + Im * Cm + Is * Cs)

The difference here is that we can handle vaccine efficacy.

What we would do with model structure arrays:

S->E ~ (1 - VE) * beta0 * (1/N) * I * C

The `I * C` piece implies summation over the products of the elements of `I` and `C` -- i.e. this is an inner product.

What the macpan vaccination model would do:

kronecker(vax_trans_red, t(baseline_trans_rates)) %*% Istate

Where `vax_trans_red` is a matrix of expressions involving parameters, `baseline_trans_rates` is a vector of parameters and `Istate` is a vector of infected states. For more detail see [here](https://github.com/mac-theobio/McMasterPandemic/blob/a1c3a4d527c134824643a777e128b3316a3ed8f1/R/models.R#L126).

# Scratch

Base Parameter Object

| parameter | model | default |
| -------- | ----- | -------- |
| alpha | epi | 0.5 |
| sigma | epi | 0.9 |

Base State Table

| structured state object | epi dim | vax dim|
| -------- | -----| --- |
| I | Ia | NA
| I | Ip |NA
| I | Is |NA
| I | Im |NA

Grouping Map

|group key| group value |epi dim | vax dim |
| --- | ---- | --- | --- |
| prot | unprot | NA | unvax |
| prot | unprot | NA | vax |
| prot | prot | NA | prot |
| sympt | asymp | Ia | NA
| sympt | pre | Ip | NA
| sympt | mild | Im | NA
| sympt | severe | Is | NA

Parameter Variation

| base parameter | prot status | sympt_status |
| -------- | -------- | -------- |
| alpha | unprot | NA
| alpha | unprot | NA
| alpha | prot | NA
| sigma | NA | NA
| gamma_a | NA | NA
| gamma_s | NA. | NA | NA
| nonhosp_mort | NA. | NA | NA
| ve | unprot | NA
| ve | prot | NA
| C | NA | asympt |
| C | NA | pre |
| C | NA | mild |
| C | NA | severe |

Product Model Parameter Object

| parameter | defaut |
| -------- | -------- |
| alpha_unprot | 0.5 |
| alpha_prot | 0.9 |