Matrix-Based Models¶
Intro¶
We want to define and integrate a set of \(N_c\) ordinary differential equations (ODEs) as
To this end, epipack.numeric_matrix_epi_models.MatrixEpiModel
contains three
parameter-carrying linear algebra objects, the numpy array birth_rates
that is equal to the vector \(\gamma\),
the scipy sparse matrix linear_rates
that is equal to the matrix \(\beta\) and a list
of scipy sparse matrices quadratic_rates
that contains a scipy sparse matrix in each of
its entries i, which makes it correspond to the tensor \(\alpha\).
Choice of Data Structures¶
We chose scipy sparse matrices because their API is simple and tailored for
fast execution of dot products. The ODEs are computed using a current-state vector
\(Y\) and dot products between the operators linear_rates
and quadratic_rates
.
One might wonder about the specific choice of sparse matrices instead of arrays. The reasoning
here is that these models can be potentially used to set up reaction-diffusion systems where a spatial
dimension is introduced using a (graph) Laplacian. In order to efficiently expand the system
to evolve on every node on a network (discretized point in space), the operators birth_rates
,
linear_rates
, and quadratic_rates
can be set system-wide using the Kronecker product
(see scipy.sparse.kron).
Setting Rates¶
There are two ways to set rates
Directly. To this end, use the methods
model.set_linear_rates
andmodel.set_quadratic_rates
. Rates are encoded as a list of tuples that contain the coupling compartments, the affected compartment, and the rate value. Methods whose name begins withset
will override previously set rates.Through proccesses. To this end, use the methods
model.set_processes
ormodel.add_*_processes
. Processes are encoded as a list of tuples that contain the coupling compartments, a rate value, and the affected compartment. Methods whose name begins withset
will override previously set rates, while methods whose name begins withadd
will not override previously set processes.
When rates for quadratic processes are set, the components that are affected by quadratic processes are
saved in model.affected_by_quadratic_process
as integers. This saves time in situations where there are
many compartments but only few transmission processes.
Integrating ODEs¶
The momenta \(dY_i/dt\) are automatically computed by means of dot products. They are integrated using
a Runge-Kutta 4(5) method with adaptive step size (Dormand-Prince) (see epipack.integrators.integrate_dopri5()
).