search

epiverse-trace / finalsize

R package to calculate the final size of an SIR epidemic in populations with heterogeneity in social contacts and disease susceptibility
https://epiverse-trace.github.io/finalsize/
Other
11 stars 7 forks source link

Move finalsize specific element from solver function to main finalsize() function #68

Open Bisaloo opened 2 years ago

Bisaloo commented 2 years ago

It would be good to refactor the code to extract the parts specific to finalsize from the solver functions. In other words, the solver function should only contain the solver algorithm and nothing more.

https://github.com/epiverse-trace/finalsize/blob/185e99f1f14f1d885f78825540d6dd02e73736b4/R/iterative_solver.R#L30-L48

pratikunterwegs commented 2 years ago

That's a fair point - but I see for example that both solvers use the contact matrix in their algorithms (via helper functions , which is finalsize-specific, if I've understood what you mean. Will look into how this could be separated from the solver algorithms.

Bisaloo commented 2 years ago

Yes, generic formulation of these solvers usually take a function to optimize (with extra arguments passed as ...). I think all you will need to do is pass extra variables to these helper functions instead of relying on global variables.

pratikunterwegs commented 2 years ago

Hi @BlackEdder could you weigh in on this? That would help me get a sense of where this might fit in the current development pipeline, thanks!

pratikunterwegs commented 1 year ago

Could I get an opinion on using a solver from Boost instead of our own custom solvers? See https://www.boost.org/doc/libs/1_81_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/getting_started/short_example.html. This would require moving away from using Eigen::arrays to using std::vectors or Boost::arrays instead.

Using Boost would be straightforward via the {bh} package. The benefits include reducing the codebase we have to maintain, and instead using a widely used and high quality library (Boost). I could roll this into #146 if there is support for this idea. Thoughts @Bisaloo and @BlackEdder?

PS. Discussing this with Edwin in person in a short while.

Bisaloo commented 1 year ago

Would this make #150 obsolete? Do we expect other users or packages to use anything else than the solvers provided in finalsize?

pratikunterwegs commented 1 year ago

Just raised this point in the discussion with Edwin and Roz, and there may be some interest in using the epi preprocessing functions, or the functions that are passed to the solvers in other packages (although to be fair I don't know of such a case just yet - perhaps {epidemics}?).

Edwin points out that the solvers are optimised for edge cases of $R_0 \approx$ 1.0, where other solvers that he has tried (from stats::optim) did not work for these cases.

The course of action we agreed upon was to tackle #68 in #150, and later assess whether one of the Boost solvers would be a suitable replacement in a separate issue and potential PR.

Bisaloo commented 1 year ago

The course of action we agreed upon was to tackle https://github.com/epiverse-trace/finalsize/issues/68 in https://github.com/epiverse-trace/finalsize/pull/150, and later assess whether one of the Boost solvers would be a suitable replacement in a separate issue and potential PR.

Yes, this sounds like a solid strategy :+1:

pratikunterwegs commented 1 year ago

Hi @Bisaloo, I've moved the code section highlighted in the issue from the solver functions to the main R function, and this is part of #150. Could you say whether this is what you intended as the solution to the issue?

Bisaloo commented 1 year ago

No, to be able to use the solvers in other contexts, you would need to allow the users to pass your f (in the iterative solver) and your dx_f (in the newton solver) as arguments.

pratikunterwegs commented 1 year ago

Sure, so an implementation similar to deSolve::lsoda(). I think there's some debate as to whether passing functions that are going to be called in a loop as Rcpp::Function is a good idea, but can look into that.

pratikunterwegs commented 1 year ago

No, to be able to use the solvers in other contexts, you would need to allow the users to pass your f (in the iterative solver) and your dx_f (in the newton solver) as arguments.

Looking into this further, I'm not able to see a change to the solver functions that is compatible with this implementation, and also useful to potential users including ourselves in other packages such as {epidemics}.

I see a use case where a user is implementing an ODE model in Rcpp, and for each step of the model wants to calculate the expected final epidemic size from a contact matrix, demography vector, and matrix of susceptibility/p_susceptibility - which would already be in use for the ODE model. For this use case, exporting the solver functions in the finalsize C++ namepsace, with the appropriate epi processing of the contact matrix included within each function, is probably a useful functionality to have.

To implement changes to the solvers suggested here:

  1. f() from iterative_solver.h and dx_f() from newton_solver.h would have to be defined as functions that take the contact matrix, demography, susceptibility, and some initial conditions as arguments (do-able);
  2. Reformulate the solver functions to take a function as an argument (this is do-able with std::function).

Overall, this increases the complexity of the codebase.

And to use this in future Rcpp packages:

  1. Users would have to process the contact matrix correctly, since that step has now been moved to final_size.R following https://github.com/epiverse-trace/finalsize/pull/150/commits/11dea230b2f1afdaee11ba38cc52a134083f8cb6 and https://github.com/epiverse-trace/finalsize/pull/150/commits/daf50ffae254104806afaa313245460e2f96d434.
  2. Users would have to rely on the functions f() and dx_f() in the C++ namespace finalsize, and pass these to the solvers, also from the same namespace - I'm not sure that these functions can be used separately from the solvers they are associated with; or
  3. Users would have to write their own equivalents of f() or dx_f() to pass to one of the solvers.

If (2), which I see as the case for less advanced users, they would be better off using the solver and intermediate functions we provide, and if (3), for potentially more advanced users, they could be better off using an available solver from say, Boost.

Overall I would say that the goal of {finalsize} is not to replace or replicate {deSolve} or odeint. I would propose to drop this issue (and also point out that this was opened on the R-only version), and export the processing+solver functions as they were prior to this most recent discussion: this would require https://github.com/epiverse-trace/finalsize/pull/150/commits/11dea230b2f1afdaee11ba38cc52a134083f8cb6 and https://github.com/epiverse-trace/finalsize/pull/150/commits/daf50ffae254104806afaa313245460e2f96d434 to be reverted - thoughts welcome @BlackEdder @TimTaylor and @Bisaloo.

Bisaloo commented 1 year ago

I would argue that this actually decreases code complexity because it clearly separates conceptually distinct elements.

(3), for potentially more advanced users, they could be better off using an available solver from say, Boost.

It would indeed be interesting to determine how the custom solvers compare to solvers from Boost (which is the goal from ) but I do think this issue will have to be addressed in any case:

pratikunterwegs commented 1 year ago

Do you mean the separation of the initial matrix processing from further steps (including the intermediate functions), or both the processing, and the intermediate functions f() and dx_f(), and the further steps into 3 separate components?

pratikunterwegs commented 1 year ago

I've been looking into implementing the fix for this issue, and my evaluation is that implementing these changes will be technically challenging and time consuming, while delivering no additional benefit to users. One implementation I considered would even reduce performance as it requires accessing and copying list elements multiple times.

I think the changes in #150 that allow importing the C++ code into other packages are still useful, after reverting the final two commits. I'm appending the 'wontfix' and 'help wanted' tags here - if anybody with the required programming skill and time would like to take this up, I'm happy to review the PR.