Parallelising rTPC

[padpadpadpad]

I have used the new purrr to parallelise fitting of TPCs, but have also changed fundamentally how multiple models are fit. Instead of fitting each model using a separate map call, we write a custom function of all the models we want to fit that we then pass to purrr. This makes it much easier to play ball with how purrr wants to parallelise. And I think should be the default way we try and fit going forwards. Maybe having an argument for models to fit and then returning a dataframe with them all. It may of course allow people to think less, but this seems to be the best way going forwards and avoids code repeats. Probably quite similar to #77.

Proof of principle code is here. It goes from 48s to 14s, a pretty massive speed up when you have LOADS of curves. This approach could replace all of the vignettes we have, at least the model fitting ones.

Open to thoughts @fwimp

# if librarian is not installed, install it
if (!requireNamespace("librarian", quietly = TRUE)) {
  install.packages("librarian")
}
# load packages
librarian::shelf(
  tidyverse,
  rTPC,
  nls.multstart,
  progressr,
  microbenchmark,
  mirai
)

## ---------------------------

# load in data
data("chlorella_tpc")

d <- chlorella_tpc

# just some of the models, edit as wanted
d2 <- filter(d, curve_id <= 60)

# ever more complex uses of parallelisation of purrr and usign rtpc

#----------------------------------------------#
# 1. just fit a single model to the dataset ####
#----------------------------------------------#

# write a function to fit a single model

# fit single rTPC
fit_gaussian <- function(d, ...) {
  # get start values and fit model
  start_vals <- rTPC::get_start_vals(
    d$temp,
    d$rate,
    model_name = 'gaussian_1987'
  )
  lower <- rTPC::get_lower_lims(d$temp, d$rate, model_name = 'gaussian_1987')
  upper <- rTPC::get_upper_lims(d$temp, d$rate, model_name = 'gaussian_1987')

  # fit model
  mod <- nls.multstart::nls_multstart(
    rate ~ rTPC::gaussian_1987(temp = temp, rmax, topt, a),
    data = d,
    iter = c(4, 4, 4),
    start_lower = start_vals - 10,
    start_upper = start_vals + 10,
    lower = rTPC::get_lower_lims(d$temp, d$rate, model_name = 'gaussian_1987'),
    upper = rTPC::get_upper_lims(d$temp, d$rate, model_name = 'gaussian_1987'),
    supp_errors = 'Y',
    convergence_count = FALSE
  )

  return(mod)
}

# fit gaussian without parallelisation
start_normal <- Sys.time()

d_mods <- d2 %>%
  nest(data = c(temp, rate)) %>%
  mutate(
    mods = map(
      data,
      ~ fit_gaussian(d = .x)
    ),
    .progress = TRUE
  )

end_normal <- Sys.time()

print(end_normal - start_normal)

# create a parallel plan
start_time_parallel <- Sys.time()

# try parallelisation
daemons(3)

d_mods1 <- d2 %>%
  nest(data = c(temp, rate)) %>%
  mutate(
    mods = map(
      data,
      in_parallel(
        \(x) fit_gaussian(d = x),
        fit_gaussian = fit_gaussian
      ),
      .progress = TRUE
    )
  )


daemons(0)

end_time_parallel <- Sys.time()

# print time taken
print(end_time_parallel - start_time_parallel)

#------------------------------------------#
# 2. fit multiple models to the dataset ####
#------------------------------------------#

# instead of running multiple map calls within a single mutate call, we can include a single model fitting function to fit the models we want

# fit single rTPC
fit_TPCs <- function(d, ...) {
  # get start values and fit model
  start_vals <- rTPC::get_start_vals(
    d$temp,
    d$rate,
    model_name = 'gaussian_1987'
  )

  # fit model
  mod <- nls.multstart::nls_multstart(
    rate ~ rTPC::gaussian_1987(temp = temp, rmax, topt, a),
    data = d,
    iter = c(4, 4, 4),
    start_lower = start_vals - 10,
    start_upper = start_vals + 10,
    lower = rTPC::get_lower_lims(d$temp, d$rate, model_name = 'gaussian_1987'),
    upper = rTPC::get_upper_lims(d$temp, d$rate, model_name = 'gaussian_1987'),
    supp_errors = 'Y',
    convergence_count = FALSE,
    ...
  )

  # get start values and fit model
  start_vals <- rTPC::get_start_vals(
    d$temp,
    d$rate,
    model_name = 'boatman_2017'
  )

  # fit model
  mod1 <- nls.multstart::nls_multstart(
    rate ~ rTPC::boatman_2017(temp = temp, rmax, tmin, tmax, a, b),
    data = d,
    iter = c(4, 4, 4, 4, 4),
    start_lower = start_vals - 10,
    start_upper = start_vals + 10,
    lower = rTPC::get_lower_lims(d$temp, d$rate, model_name = 'boatman_2017'),
    upper = rTPC::get_upper_lims(d$temp, d$rate, model_name = 'boatman_2017'),
    supp_errors = 'Y',
    convergence_count = FALSE,
    ...
  )

  # get start values and fit model
  start_vals <- rTPC::get_start_vals(
    d$temp,
    d$rate,
    model_name = 'sharpeschoolhigh_1981'
  )

  # fit model
  mod2 <- nls.multstart::nls_multstart(
    rate ~
      rTPC::sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 20),
    data = d,
    iter = c(3, 3, 3, 3),
    start_lower = start_vals - 10,
    start_upper = start_vals + 10,
    lower = rTPC::get_lower_lims(
      d$temp,
      d$rate,
      model_name = 'sharpeschoolhigh_1981'
    ),
    upper = rTPC::get_upper_lims(
      d$temp,
      d$rate,
      model_name = 'sharpeschoolhigh_1981'
    ),
    supp_errors = 'Y',
    convergence_count = FALSE,
    ...
  )

  return(tibble::tibble(
    gaussian = list(mod),
    boatman = list(mod1),
    sharpeschoolhigh = list(mod2)
  ))
}

# fit all models without parallelisation
start_normal <- Sys.time()

d_mods <- d2 %>%
  nest(data = c(temp, rate)) %>%
  mutate(
    mods = map(
      data,
      ~ fit_TPCs(d = .x),
      .progress = TRUE
    )
  ) %>%
  unnest(mods)

end_normal <- Sys.time()

print(end_normal - start_normal)

# create a parallel plan
start_time_parallel <- Sys.time()

# try parallelisation
daemons(6)

d_mods1 <- d2 %>%
  nest(data = c(temp, rate)) %>%
  mutate(
    mods = map(
      data,
      in_parallel(
        \(x) fit_TPCs(d = x),
        fit_TPCs = fit_TPCs
      ),
      .progress = TRUE
    )
  ) %>%
  unnest(mods)

daemons(0)

end_time_parallel <- Sys.time()

print(end_time_parallel - start_time_parallel)

[fwimp]

Oh this is awesome! I think with a bit of clever design, a lot of this can leverage the new funcs from #77 combined with some of the fun code I originally wrote in #59 (comment) .

Shall I try implementing this on the branch where I'm working on #77?

In terms of the "allowing people not to think" I believe that, as in #77, providing some basic ability to get things done in a sensible way is a good thing, however if people want more flexibility then they will definitely want to be writing their own custom fitting funcs anyway. We should provide not only an easy general solution, but break down how it's put together in the vignette so people can easily build their own in a similar style.

Also that speedup is WILD!

[padpadpadpad]

Yes this did remind me of some of the code you wrote in that comment!

Implementing in that branch sounds great.

Having a single map call in this way also allows the native progress bar that exists in purrr now to work out the box which is so much nicer

Also obviously happy to help where I can @fwimp. Let me know if there are things I can do. Other than work on a paper to actually get recognition for all these changes...

[fwimp]

Having a go at it now.

Interesting question just popped up for me: how do you feel about being able to specify bare model functions as well as string representations? There are basically 3 levels of this with varying levels of complexity under the hood to implement:

1. Just take models as strings (e.g. c("ashrafi1_2018", "briere1_1999", "briere2_1999") - very easy to do, requires no changes.
2. Handle bare functions (e.g. c(ashrafi1_2018, briere1_1999, briere2_1999)) - a bit trickier, but much more friendly for users that are using autocomplete.
3. Handle called functions (e.g. c(ashrafi1_2018(), briere1_1999(), briere2_1999()) - requires a whole pass over all major functions and some non-standard function behaviour (returning itself called with no args). Allows for really easy autocomplete stuff.

Personally I'm happy with either 1 or 2, but I think 3 would cause more issues than it'd solve day-to-day

Edit: I've actually worked out a nice function to do this, so implementing 2 is really pretty easy:

Edit2: ignore that, it's harder than I wanted. Switching back to 1 for now.

[fwimp]

Right, looks like that internal fit_TPCs() function can be reworked to use the quickfit_tpc() function like so:

.quickfit_multi <- function(data, model_names, temp, trait, start_adjusts = 0, iter = 150, lhstype = "none", gridstart=FALSE, force=FALSE) {
  # Coerce args (basically make df of params)
  df <- tibble::tibble(model=model_names, temp, trait, start_adjusts, iter, lhstype , gridstart, force)
  
  # Map tpcs using quickfit
  out <- pmap(df, \(...) list(quickfit_tpc(data, ...)))
  
  # Get some lovely names
  names(out) <- model_names
  
  # Cast to tibble and return
  return(tibble::as.tibble(out))
}

That's the first step, then I'll just need to adapt this code to use .quickfit_multi() and I think we're off to the races!

d_mods1 <- d2 %>%
  nest(data = c(temp, rate)) %>%
  mutate(
    mods = map(
      data,
      in_parallel(
        \(x) fit_TPCs(d = x),
        fit_TPCs = fit_TPCs
      ),
      .progress = TRUE
    )
  ) %>%
  unnest(mods)

Of course, that's always easier said than done. :)

[fwimp]

Final update of the night: I now have a somewhat working version of this.

Currently it's a little hacky (in order to use functions that are not currently part of the installable package [i.e. quickfit_tpc() and get_tpc_as_formula()] within a parallel map, you need to redefine those functions within the fitting function [i.e. .quickfit_multi()]).

You also need to do some explicit scope and environment management to make the stuff play nicely, but that's all okay.

I'll look at this tomorrow to see if I can tidy it up a bit. But writing the final code will likely mean I need to test by constantly actually installing the package...

(Also uncovered a bug in purrr, which I need to raise with them)

[padpadpadpad]

Coding at midnight is not advised, but this looks great.

I am not sure there is anyway around defining your functions within quickfit_multi() if you are using any custom function.

And I agree 1. seems fine. 2 would be ncie given how much people use the tidyverse and take bare unquoted functions for granted, and 3 seems not as useful.

[fwimp]

I think that if the functions are in the installed version of the package then they work if you reference them with their full :: name, however this does not include ad-hoc loading of the dev package with devtools. A hacky fix for development might be to force the workers to load the package using devtools. However I don't love that. I'll play with properly installing the dev version and see how that goes.

And agreed, I'll put 2 on my todo list as I have it nearly sorted, but not quite there yet.

[padpadpadpad]

I dont understand what you mean by devtools. If you have installed the development package you should still be able to :: functions and it all work inside the fit_TPC() or whatever we call it. It is if you have custom functions not in the package that it will start having issues, is my (admittedly rudimentary) understanding of the new purrr::map()

[fwimp]

So the way I work with packages when developing is constantly reloading with devtools::load_all(".") which loads the package temporarily. These functions are not propagated to the mirai demons, which instead need to (by default) run with the installed packages (rather than any new ones you've added since you last installed the package).

I just found that mirai have actually thought of this and so included the mirai::everywhere() command, which allows you to make every daemon run some code. So now just working on getting everything pulled together.

So for testing purposes you'll want something like

if (!mirai::daemons_set()) {
    mirai::daemons(6)
}
# TESTING ONLY
mirai::everywhere(devtools::load_all("."))

# Other stuff...

[fwimp]

One thing that's quite an important decision: mirai recommends that the daemons calls are all handled by the user. Should we follow that paradigm and just assume(/check) that mirai daemons are running?

[fwimp]

@padpadpadpad I've got all of these parts now implemented and (as far as I can tell) working.

I can put in a PR should you wish, or if you'd like to test it out a bit yourself then feel free to do that!

Just as a note as well: I went through the suggests field and removed any packages that were specified on their own AND as a part of tidyverse (as that's a suggestion) and also bumped the working version up to 1.1.0 (I think we've done enough since 1.0.4 that 1.1.0 really makes sense here).

[padpadpadpad]

Awesome @fwimp. Do you want to submit a PR and we can then workout what vignettes might need changing - think most need a big old rewrite if we're not having multiple map() inside a mutate() anymore. I can make a start on those.