recovery_rates variable has been renamed as (disease) progression_rates.

[RJSheppard]

I've created two new functions to standardise drug and resistance calculations across treatment and mass treatments:

calculate_successful_treatments
and
update_mass_drug_admin

calculate_successful_treatments handles the durg efficacy and antimalarial resistance, generating key outputs for each intervention type and returning a bitset of those treated, a vector of the drugs used to treat those individuals and a vector of treated delays for those individual if resistance (SPC) is switched on. This function is called during calculate_treated, create_mda_listeners and create_pmc_process.

update_mass_drug_admin calls update_infection twice, assigning D and LM-detectable A infections to Tr, while non-LM detectable A infections, U and Tr infections to S. Those assigned to S have 0 infectivity and 0 dt, while those who are assigned to Tr are assigned a fraction of their infectivity in the current timestep and dt from the parameter set or calculated using antimalarial resistance.

[giovannic]

Great stuff. Thank you for this!

I've requested some minor documentation and performance improvements. But almost there for me.

[giovannic]

Suggested change

	#' @param target bitset for individuals who have been successfully treated
	#' @param target a list containing the successfully treated, the drug used and resistance parameters

[giovannic]

Suggested change

	dt_update_vector <- target$dt_spc_combined[target$successfully_treated$to_vector() %in% to_treated$to_vector()]
	dt_update_vector <- target$dt_spc_combined[
	target$successfully_treated$copy()$and(to_treated)$to_vector()
	]

[giovannic]

Suggested change

	#' @param renderer simulation renderer
	#' @noRd
	#' @param renderer simulation renderer
	#' @param int_name the intervention name to use for rendering, use "" for frontline treatment
	#' @noRd

[giovannic]

Suggested change

	treated <- bitset_at(individual::Bitset$new(parameters$human_population)$not(), target)
	treated <- individual::Bitset$new(parameters$human_population)$insert(target)

[giovannic]

Suggested change

	treated <- bitset_at(individual::Bitset$new(parameters$human_population)$not(), target)
	treated <- individual::Bitset$new(parameters$human_population)$insert(target)

[RJSheppard]

  # Set the size of the human population and an initial entomological inoculation rate (EIR)
  human_population <- 10000
  starting_EIR <- 10`
  
  simparams <- get_parameters(
    list(human_population = human_population)
  )
  
  # Use the set_equilibrium() function to update the individual-based model parameters to those 
  # required to match the specified initial EIR (starting_EIR) at equilibrium:
  simparams <- set_equilibrium(simparams, starting_EIR)
  
  # Update the parameter list with the default parameters for sulphadoxine-pyrimethamine
  # amodiaquine (SP-AQ)
  simparams <- set_drugs(simparams, list(AL_params))
  
  # Use set_mda() function to set the proportion of the population that the MDA reaches and the age 
  # ranges which can receive treatment.
  clinparams <- set_clinical_treatment(
    simparams,
    drug = 1,
    timesteps = 1:2 * 365,
    coverages = c(1,0)
  ) |> 
    set_equilibrium(, init_EIR = starting_EIR)
  
  clin_output <- run_simulation(365 * 3, clinparams)
  
  
  # Use set_mda() function to set the proportion of the population that the MDA reaches and the age 
  # ranges which can receive treatment.
  mdaparams <- set_mda(
    simparams,
    drug = 1,
    timesteps = 365,
    coverages = 0.99,
    min_ages = 0,
    max_ages = 200 * 365
  )
  
  mda_output <- run_simulation(365 * 3, mdaparams)
  
  # PMC 
  pmcparams <- set_drugs(simparams, list(SP_AQ_params))
  pmcparams <- set_pmc(
    pmcparams,
    drug = 1,
    timesteps = 1:2*365,
    coverages = c(0.99,0),
    ages = seq(30, 365*100, by = 30)
  )
  
  pmcparams$prevalence_rendering_min_ages <- 0
  pmcparams$prevalence_rendering_max_ages <- 365*100
  pmc_output <- run_simulation(365 * 3, pmcparams)
  
  plot(
    x = clin_output$timestep,
    y = clin_output$n_detect_lm_730_3650/clin_output$n_age_730_3650,
    xlab = "Time (days)",
    ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.8,
    ylim = c(0, 1),
    type = "l",
    main = paste0(name, ": clin")
  )
  lines(x = clin_output$timestep,
        y = clin_output$n_detect_pcr_730_3650/clin_output$n_age_730_3650, col = "blue")
  abline(v = 1:2*365, lty = 2)
  
  plot(
    x = mda_output$timestep,
    y = mda_output$n_detect_lm_730_3650/mda_output$n_age_730_3650,
    xlab = "Time (days)",
    ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.8,
    ylim = c(0, 1),
    type = "l",
    main = paste0(name, ": MDA")
  )
  lines(x = mda_output$timestep,
        y = mda_output$n_detect_pcr_730_3650/mda_output$n_age_730_3650, col = "blue")
  abline(v = 365, lty = 2)
  
  plot(
    x = pmc_output$timestep,
    y = pmc_output$n_detect_lm_0_36500/pmc_output$n_age_0_36500,
    xlab = "Time (days)",
    ylab = expression(paste(italic(Pf),"PR"[0-100])), cex = 0.8,
    ylim = c(0, 1),
    type = "l",
    main = paste0(name, ": PMC")
  )
  lines(x = pmc_output$timestep,
        y = pmc_output$n_detect_pcr_0_36500/pmc_output$n_age_0_36500, col = "blue")
  abline(v = 1:2*365, lty = 2)
  
}

par(mfrow = c(2,3))
# remove.packages("malariasimulation")
remotes::install_github("mrc-ide/malariasimulation@master")
treatment_testing_function("master")

remove.packages("malariasimulation")
remotes::install_github("mrc-ide/malariasimulation@progression_rates_mass_admin")
treatment_testing_function("progression_rates_mass_admin")

[RJSheppard]

[pwinskill]

Thanks @RJSheppard - looks good to me