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

[Bisaloo]

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.

• it will avoid redundant code across solvers:

finalsize/R/iterative_solver.R

Lines 30 to 48 in 185e99f

	# count demography groups
	n_dim <- length(demography_vector)
	# make vector of zeros
	zeros <- rep(0.0, n_dim)
	# make vector of initial final size guesses = 0.5
	epi_final_size <- rep(0.5, n_dim)
	# set contact_matrix values to zero if there are no contacts among
	# demography groups, or if demography groups are empty
	i_here <- demography_vector == 0 | susceptibility == 0 |
	rowSums(contact_matrix) == 0
	zeros[i_here] <- 1.0
	epi_final_size[i_here] <- 0.0
	# matrix filled by columns
	contact_matrix_ <- contact_matrix * demography_vector %o% susceptibility
	contact_matrix_[i_here, i_here] <- 0.0

• it will make it easier to plug external solvers, such as solvers provided by base R (e.goptim()) or other packages.

[pratikunterwegs]

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]

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]

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]

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]

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

[pratikunterwegs]

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]

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.

Yes, this sounds like a solid strategy 👍

[pratikunterwegs]

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]

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]

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]

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 11dea23 and daf50ff.
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 11dea23 and daf50ff to be reverted - thoughts welcome @BlackEdder @TimTaylor and @Bisaloo.

[Bisaloo]

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:

• either solvers from Boost are better and they should be used in finalsize, which as far as I can tell, will require solving the present issue first
• or solvers from finalsize are better (even if only in some contexts) and then there is value for external users to use them, which requires solving the present issue.

[pratikunterwegs]

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]

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.