Usage of integrate method with discrete solver (using stochastic model) confusing

[stijnvanhoey]

When checking the current models in the models.py, two out of the three models (COVID19_SEIRD_sto and COVID19_SEIRD_sto_spatial) will only work with the discrete=True option and do have an integrate method whereas it does not do any integration, but an update (cfr. as if the solver is explicit Euler).

It is well-explained in the example notebook (nice!):

Initialize the model, the only difference with the stochastic model is the discrete flag

However, I do think this is very confusing to other users/developers as one is not only required to use the discrete=True, but the setup of the integrate/update function changes completely.

Although solver-choice should in fact not matter, I do think that the current implementation of the stochastis models are probably only possible with the discrete step? As such, maybe it makes more sense to have a BaseModel which contains code that is independent from the solver-choice and have two subclasses:

1. DeterministicModel: solve_ivp (scipy ode solvers) oriented, which requires a subclass with an integrate method
2. StochasticModel: solve_discrete oriented, as used by the stochastic/discrete models with an update method

COVID19_SEIRD will subclass from DeterministicModel, whereas COVID19_SEIRD_sto and COVID19_SEIRD_sto_spatial from StochasticModel?

@twallema, @jorisvandenbossche and @JennaVergeynst any thoughts?

Maybe important to consider: which models are used for the research? @twallema and @mrollier, does it make sense to invest more time into the integration with the scipy solver or are you only using the discrete versions of the model?

[JennaVergeynst]

Two separate subclasses seems a good idea to me.

Considering model use: for predictions we are currently using the continuous version (COVID19_SEIRD), so this one is still important. To facilitate calibration, I put all calibration steps in a function (full_calibration_wave1 in covid19model/optimization/run_optimization.py). The idea is to have such a function for each wave (assuming that we have one ideal calibration sequence per wave that we are all using). @mrollier will probably need an adapted version of this function for calibration of the spatial model (say e.g. spatial_calibration_wave1). Does this seem a good approach to you?