diff --git a/DESCRIPTION b/DESCRIPTION index befe80e..3f6b32a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: cali Title: Good vibes and model calibration -Version: 0.2.1 +Version: 0.2.3 Authors@R: c( person("Pete", "Winskill", email = "p.winskill@imperial.ac.uk", role = c("aut", "cre")) ) diff --git a/R/calibrate.R b/R/calibrate.R index d6731af..b352574 100644 --- a/R/calibrate.R +++ b/R/calibrate.R @@ -6,10 +6,10 @@ #' @param target Values of target variable to calibrate to. #' @param summary_function A function that take the raw model output as an argument and #' produces a vector of the target variable. -#' @param tolerance The routine will complete when the absolute sum of the weighted difference between the target variable +#' @param tolerance The routine will complete when the average absolute weighted difference between the target variable #' and the target values falls below this value. Tolerance is specified in #' the units of the target variable (e.g. if my target variable is prevalence, and my tolerance is 0.05, -#' then the absolute sum of the weighted difference between the target variable and the target values must be <5% for the routine to succeed). +#' then the average absolute sum of the weighted difference between the target variable and the target values must be <5% for the routine to succeed). #' @param weights A numerical vector of weights the same length as target giving the weights to use for elements of target. #' @param elimination_check Check transmission is maintained for all target points with ongoing transmission before exiting early. #' @param low Lower boundof EIRs @@ -26,48 +26,18 @@ calibrate <- function(target, elimination_check = TRUE, maxiter = 20, low = 0.001, high = 2000){ - - x <- proposal(low, high) - - for(i in 1:maxiter){ - if(low == high) break - - y <- objective(x = x, - parameters = parameters, - summary_function = summary_function) - - difference <- y - target - weighted_difference <- difference * weights - print(signif(rbind(y, target, difference, weighted_difference)), 3) - diff <- sum(weighted_difference) - abs_diff <- sum(abs(weighted_difference)) - - # Can stop early if close enough and transmission maintained - within_tol <- abs_diff < tolerance - transmisison <- TRUE - if(elimination_check){ - transmission <- all(y[target > 0] > 0) - } - if(within_tol & transmission) break - - if(diff < 0){ - low <- x - x <- proposal(x, high) - } - if(diff > 0){ - high <- x - x <- proposal(low, x) - } - } - return(x) -} -#' Propose new EIR, moving on log scale -#' -#' @param a lower EIR -#' @param b upper EIR -#' -#' @return EIR -proposal <- function(a, b){ - exp(log(a) + (log(b) - log(a)) / 2) -} \ No newline at end of file + x <- stats::uniroot( + f = objective, + lower = low, + upper = high, + maxiter = maxiter, + parameters = parameters, + summary_function = summary_function, + target = target, + weights = weights, + tolerance = tolerance, + elimination_check = elimination_check + ) + return(x$root) +} diff --git a/R/objective.R b/R/objective.R index 2e06988..3916fdc 100644 --- a/R/objective.R +++ b/R/objective.R @@ -6,12 +6,33 @@ #' @inheritParams calibrate #' #' @return Difference between output and target. -objective <- function(x, parameters, summary_function){ +objective <- function(x, parameters, summary_function, target, weights, tolerance, elimination_check){ message("\nTrying EIR: ", signif(x, 5)) p <- malariasimulation::set_equilibrium(parameters, init_EIR = x) raw_output <- malariasimulation::run_simulation(timesteps = p$timesteps, parameters = p) model_output <- summary_function(raw_output) - - return(model_output) + difference <- model_output - target + weighted_difference <- difference * weights + + print(signif(rbind(model_output, target, difference, weighted_difference)), 3) + + if(elimination_check){ + bad_eliminated <- model_output == 0 & target != 0 + if(sum(bad_eliminated) > 0){ + message("Unwanted elimination") + weighted_difference[bad_eliminated] <- -1e6 + } + } + + mean_weighted_difference <- mean(weighted_difference) + + message("\nSum squared weighted difference: ", signif(mean_weighted_difference, 5)) + + if(mean(abs(weighted_difference)) < tolerance){ + message("Mean absolute difference < tolerance") + mean_weighted_difference <- 0 + } + + return(mean_weighted_difference) } \ No newline at end of file diff --git a/man/calibrate.Rd b/man/calibrate.Rd index 20c0325..be811f8 100644 --- a/man/calibrate.Rd +++ b/man/calibrate.Rd @@ -24,10 +24,10 @@ produces a vector of the target variable.} \item{parameters}{Other malariasimulation model parameters} -\item{tolerance}{The routine will complete when the absolute sum of the weighted difference between the target variable +\item{tolerance}{The routine will complete when the average absolute weighted difference between the target variable and the target values falls below this value. Tolerance is specified in the units of the target variable (e.g. if my target variable is prevalence, and my tolerance is 0.05, -then the absolute sum of the weighted difference between the target variable and the target values must be <5\% for the routine to succeed).} +then the average absolute sum of the weighted difference between the target variable and the target values must be <5\% for the routine to succeed).} \item{weights}{A numerical vector of weights the same length as target giving the weights to use for elements of target.} diff --git a/man/objective.Rd b/man/objective.Rd index 5ceef1f..e108481 100644 --- a/man/objective.Rd +++ b/man/objective.Rd @@ -4,7 +4,15 @@ \alias{objective} \title{Objective function} \usage{ -objective(x, parameters, summary_function) +objective( + x, + parameters, + summary_function, + target, + weights, + tolerance, + elimination_check +) } \arguments{ \item{x}{Input parameter to be varied. For malariasimulation this is EIR.} @@ -13,6 +21,17 @@ objective(x, parameters, summary_function) \item{summary_function}{A function that take the raw model output as an argument and produces a vector of the target variable.} + +\item{target}{Values of target variable to calibrate to.} + +\item{weights}{A numerical vector of weights the same length as target giving the weights to use for elements of target.} + +\item{tolerance}{The routine will complete when the average absolute weighted difference between the target variable +and the target values falls below this value. Tolerance is specified in +the units of the target variable (e.g. if my target variable is prevalence, and my tolerance is 0.05, +then the average absolute sum of the weighted difference between the target variable and the target values must be <5\% for the routine to succeed).} + +\item{elimination_check}{Check transmission is maintained for all target points with ongoing transmission before exiting early.} } \value{ Difference between output and target. diff --git a/man/proposal.Rd b/man/proposal.Rd deleted file mode 100644 index 2e8aeb8..0000000 --- a/man/proposal.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/calibrate.R -\name{proposal} -\alias{proposal} -\title{Propose new EIR, moving on log scale} -\usage{ -proposal(a, b) -} -\arguments{ -\item{a}{lower EIR} - -\item{b}{upper EIR} -} -\value{ -EIR -} -\description{ -Propose new EIR, moving on log scale -} diff --git a/tests/testthat/test-objective.R b/tests/testthat/test-objective.R index 2e94a11..c7a9ee0 100644 --- a/tests/testthat/test-objective.R +++ b/tests/testthat/test-objective.R @@ -4,6 +4,38 @@ test_that("objective works", { p <- malariasimulation::get_parameters(list(human_population = 1000)) mockery::stub(objective, "malariasimulation::run_simulation", x) expect_equal(objective(x = 1, parameters = p, - summary_function = summary_mean_pfpr_2_10), 0.5) + summary_function = summary_mean_pfpr_2_10, + target = 0.5, + weights = 1, + tolerance = 0, + elimination_check = FALSE), 0) mockery::stub(objective, "malariasimulation::run_simulation", x) + expect_equal(objective(x = 1, parameters = p, + summary_function = summary_mean_pfpr_2_10, + target = 0, + weights = 1, + tolerance = 0, + elimination_check = FALSE), 0.5) + + x <- data.frame(n_detect_730_3650 = c(0, 0, 0), + n_730_3650 = 100) + mockery::stub(objective, "malariasimulation::run_simulation", x) + expect_equal(objective(x = 1, + parameters = p, + summary_function = summary_mean_pfpr_2_10, + target = 0.5, + weights = 1, + tolerance = 0.02, + elimination_check = TRUE), -1e6) + + x <- data.frame(n_detect_730_3650 = c(0, 0, 0), + n_730_3650 = 100) + mockery::stub(objective, "malariasimulation::run_simulation", x) + expect_equal(objective(x = 1, + parameters = p, + summary_function = summary_mean_pfpr_2_10, + target = 0.001, + weights = 1, + tolerance = 0.02, + elimination_check = TRUE), -1e6) })