#' Optimize a function using adaptive random search.
#'
#' Optimize an objective function using an adaptive random search algorithm.
#' The function works for both discrete and continuous optimization parameters
#' and allows for box-constraints and sets of allowed values.
#'
#' @param loc_fac Locality factor for determining the standard deviation of the sampling distribution around the current
#' position of the parameters. The initial standard deviation is normally calculated as \code{(upper - lower)/loc_fac}
#' except in cases when there are no upper or lower limits (e.g. when \code{upper=Inf} or \code{lower=-Inf}).
#' @param no_bounds_sd The standard deviation of the sampling distribution around the current
#' position of the parameters when there are no upper or lower limits (e.g. when \code{upper=Inf} or \code{lower=-Inf}).
#' @param adapt_scale The scale for adapting the size of the sampling distribution. The adaptation of the
#' standard deviation of the sampling distribution around the current
#' position of the parameters is done after \code{iter_adapt} iteration with no change in the best objective function.
#' When adapting, the standard deviation of the sampling distribution is calculated as
#' \code{(upper - lower)/(loc_fac*ff*adapt_scale)} where ff starts at 1 and increases by 1 for each adaptation.
#' @param allowed_values A list containing allowed values for each parameter \code{list(par1=c(2,3,4,5,6),par2=c(5,6,7,8))}.
#' A vector containing allowed values for all parameters is also allowed \code{c(2,3,4,5,6)}.
#' @param iter The number of iterations for the algorithm to perform (this is a maximum number, it could be less).
#' @param iter_adapt The number of iterations before adapting (shrinking) the parameter search space.
#' @param max_run The maximum number of iterations to run without a change in the best parameter estimates.
#' @param trace_iter How many iterations between each update to the screen about the result of the search.
#' @param new_par_max_it The algorithm randomly chooses samples based on the current best set of parameters. If when drawing
#' these samples the new parameter set has already been tested then a new draw is performed. After \code{new_par_max_it} draws, with
#' no new parameter sets, then the algorithm stops.
#' @param allow_replicates Should the algorithm allow parameters to have the same value?
#' @param generator A user-defined function that generates new parameter sets to try in the algorithm. See examples below.
#'
#'
#' @references \enumerate{
#' \item M. Foracchia, A.C. Hooker, P. Vicini and A. Ruggeri, "PopED, a software fir optimal
#' experimental design in population kinetics", Computer Methods and Programs in Biomedicine, 74, 2004.
#' \item J. Nyberg, S. Ueckert, E.A. Stroemberg, S. Hennig, M.O. Karlsson and A.C. Hooker, "PopED: An extended,
#' parallelized, nonlinear mixed effects models optimal design tool",
#' Computer Methods and Programs in Biomedicine, 108, 2012.
#' }
#'
#' @family Optimize
#' @inheritParams optim_LS
#'
#' @example tests/testthat/examples_fcn_doc/examples_optim_ARS.R
#' @export
optim_ARS <- function(par,
fn,
lower=NULL,
upper=NULL,
allowed_values=NULL,
#constraints=NULL,
#control=NULL,
loc_fac=4,
no_bounds_sd = par,
iter=400,
iter_adapt = 50,# number of iterations with no change before adapting search space
adapt_scale = 1,
max_run=200,
trace = TRUE,
trace_iter=5,
#par_grouping=NULL,
new_par_max_it = 200,
maximize=F, # default is to minimize
parallel=F,
parallel_type=NULL,
num_cores = NULL,
mrgsolve_model=NULL,
seed=round(runif(1,0,10000000)),
allow_replicates=TRUE,
replicates_index=seq(1,length(par)), # same value, parameters can not be the same value
generator=NULL,
...){
# constratints
# uniform instead of normal sampling
# start trace ----------------
if((trace)){
tic(name=".ars_savedTime")
wd_iter <- nchar(iter)
}
# checks---------------
if((is.null(lower) || is.null(upper)) && is.null(allowed_values) && is.null(generator)){
stop("At least 'lower' and 'upper' or 'allowed_values' or 'generator' should be supplied.")
}
if(length(lower)==1 && length(par)>1) lower <- rep(lower,length(par))
if(length(upper)==1 && length(par)>1) upper <- rep(upper,length(par))
if(!is.null(allowed_values)){
if(!is.list(allowed_values)) allowed_values <- list(allowed_values)
if(length(allowed_values) == 1 && length(par)>1) allowed_values <- rep(allowed_values,length(par))
}
if(length(replicates_index)!=length(par)){
msg <- stringr::str_glue(
'The number of parameters ({length(par)}) ',
'is not the same as the the length of replicates_index ({length(replicates_index)})!'
)
stop(msg)
}
# initialization -----------------
nullit=1
runs <- 1
ff=1
ff_scale <- 1
par0=par
fn1 <- function(par) fn(par, ...)
#npar <- max(c(length(lower),length(upper),length(allowed_values)))
npar <- length(par)
ofv_opt <- NULL
par_opt <- par
dpar=(upper-lower)/loc_fac
dpar[is.infinite(dpar)] <- no_bounds_sd
if(!is.null(seed)) set.seed(seed)
par_vec <- cell(iter,1) # investigated parameter vectors
# continuous and discrete parameters
par_type <- rep("cont",npar)
if(!is.null(allowed_values)){
for(k in 1:npar){
if(!all(is.na(allowed_values[[k]])) && length(allowed_values[[k]]>0)){
par_type[k] <- "cat"
}
}
}
# handle replicates
if(!allow_replicates){
replicates_index <- rep(1,length(par))
}
handle_replicates <- FALSE
if(length(unique(replicates_index))!=length(par)) handle_replicates <- TRUE
# parallelization
if(parallel){
parallel <- start_parallel(parallel,seed=seed,
parallel_type=parallel_type,
num_cores=num_cores,
mrgsolve_model=mrgsolve_model,
...)
on.exit(if(parallel && (attr(parallel,"type")=="snow"))
parallel::stopCluster(attr(parallel,"cluster")))
}
iter_chunk = NULL
if(is.null(iter_chunk)) if(parallel) iter_chunk <- attr(parallel,"cores") else iter_chunk <- 1
# Functions ---------------------------------------------------------------
resample <- function(x,size=1) x[sample.int(length(x),size=size)]
compare <- function(a,b) a < b
if(maximize) compare <- function(a,b) a > b
# generate parameter sets
generate_par <- function(par,par_type,allowed_values,...){
par[par_type=="cont"]=par[par_type=="cont"] +
dpar[par_type=="cont"]/ff_scale*rnorm(length(par[par_type=="cont"]))
#par[par_type=="cont" && is.infinite(par)] <- par[par_type=="cont"]+1e10/ff*rnorm(length(par[par_type=="cont"]))
# generate new discrete parameters
par[par_type=="cat"] <- vapply(allowed_values[par_type=="cat"],resample,c(0))
return(par)
}
# ------------ generate and evaluate new parameter sets
gen_par_ofv <- function (it, par, ...) {
need_new_par <- TRUE
new_par_it <- 0
#---------------- generate new parameter set
while(need_new_par && (new_par_it < new_par_max_it)){
# compute OFV for initial par if supplied
if(it==1 && !is.null(par)){
need_new_par <- FALSE
break
}
# generate new continuous parameters if not supplied
if(is.null(par[par_type=="cont"])){
par[par_type=="cont"] <- lower[par_type=="cont"]+(upper[par_type=="cont"] - lower[par_type=="cont"])/2
par[is.nan(par) && par_type=="cont"] <- 0
}
if(!is.null(generator)){
par <- do.call(generator,list(par,...))
} else {
par <- generate_par(par,par_type,allowed_values)
}
need_new_par <- FALSE
# browser()
# handle replicates
# if(!allow_replicates){
# while(any(duplicated(par))) {
# cur_pars <- par[!duplicated(par)]
# tmp_allowed <- lapply(allowed_values[duplicated(par)],setdiff,cur_pars)
# #if(any(sapply(tmp_allowed,length)==0)) stop("Not enough potential values for all parameters to have different values.")
# par[duplicated(par)] <- generate_par(par[duplicated(par)],par_type[duplicated(par)],tmp_allowed)
# }
# }
# (par <- c(rep(2,5),rep(3,4)))
# (replicates_index <-c(rep(1,5),rep(2,4)))
# handle_replicates <- T
if(handle_replicates){
for(ii in unique(replicates_index)){
#ii <- unique(replicates_index)[[1]] # for testing
par_tmp <- par[replicates_index==ii]
if(any(duplicated(par_tmp))){
while(any(duplicated(par_tmp))) {
duped <- duplicated(par_tmp)
allowed_values_tmp <- allowed_values[replicates_index==ii]
par_type_tmp <- par_type[replicates_index==ii]
cur_pars <- par_tmp[!duped]
tmp_allowed <- lapply(allowed_values_tmp[duped],setdiff,cur_pars)
if(any(sapply(tmp_allowed,length)==0)){
need_new_par <- TRUE
break
} else {
par_tmp[duped] <- generate_par(par_tmp[duped],par_type_tmp[duped],tmp_allowed)
}
}
par[replicates_index==ii] <- par_tmp
}
}
}
# # handle replicates
# if(length(unique(replicates_index))!=length(par)){
# par_i_set <- par_i_set[!(par_i_set %in% par[replicates_index==replicates_index[i]])]
# }
# Group samples that should be grouped
# if(!is.null(par_grouping)){
# for(k in unique(par_grouping)){
# par[par_grouping==k] <- resample(par[par_grouping==k])
# }
# }
# set to boundary if beyond the boundary
if(!all(is.null(upper))) par[par>upper] <- upper[par>upper]
#par=par-((par>upper)*(par-upper))
if(!all(is.null(lower))) par[par<lower] <- lower[par<lower]
#par=par+((par<lower)*(lower-par))
# check if design has already been tested
#need_new_par <- any(sapply(sapply(par_vec, length)!=0,function(x) all(x == par)))
#par_vec[sapply(par_vec, length)!=0,]
if(!need_new_par) need_new_par <- any(sapply(par_vec[sapply(par_vec, length)!=0,],function(x) all(x==par)))
#need_new_par <- any(sapply(par_vec[sapply(par_vec, length)!=0,],function(x) all(x==par)))
#browser()
new_par_it <- new_par_it + 1
} # end while
#------------- evaluate new parameter sets
if(need_new_par){
ofv <- NULL
par <- NULL
} else {
ofv <- fn1(par)
}
return(list(ofv=ofv,par=par,need_new_par=need_new_par))
} # end function
# body --------------------------------------------------------------------
for(it in 1:ceiling(iter/iter_chunk)){
start_it <- (it-1)*iter_chunk+1
stop_it <- min(it*iter_chunk,iter)
it_seq <- start_it:stop_it
# generate new parameters and OFVs
if(parallel && (attr(parallel,"type")=="multicore")){
if(!is.null(seed)) set.seed(seed+it)
res <- mapply(c,parallel::mclapply(it_seq,gen_par_ofv,par_opt,mc.cores=attr(parallel, "cores")))
} else if(parallel && (attr(parallel,"type")=="snow")){
res <- mapply(c,parallel::parLapply(attr(parallel, "cluster"),it_seq,gen_par_ofv,par_opt))
# library(doParallel)
# foreach::foreach(i=it_seq) %dopar% gen_par_ofv(i,par_opt)
} else {
res <- mapply(c,lapply(it_seq,gen_par_ofv,par_opt))
}
# handle error messages and NULL values in OFV calculation
#print(res)
res2 <- res
if(is.null(dim(res2))) res2 <- mapply(c,res2[sapply(res2,is.list)])
res2 <- res2[,!sapply(res2["ofv",],is.null),drop=F]
if(maximize){
out <- res2[,which.max(res2["ofv",])]
} else {
out <- res2[,which.min(res2["ofv",])]
# x <- res2["ofv",]$ofv
# which(x == min(x, na.rm = F))
}
# check if a unique new parameter vector was generated
if(!is.null(out$need_new_par)){
if(out$need_new_par){
cat(paste0("Maximum number of duplicate parameter samples reached\n(new_par_max_it=",
new_par_max_it,"), optimization stopped.\n"))
break
}
}
ofv <- out$ofv
par <- out$par
par_vec[it_seq] <- res["par",] # save new parameter vectors in a list
if(
(
(compare(ofv,ofv_opt) && !is.null(ofv))
||
(is.null(ofv_opt) && !is.null(ofv))
)
){
par_opt <- par
ofv_opt <- ofv
nullit=1
runs <- 1
#ff=1
} else {
nullit=nullit+(length(it_seq))
runs <- runs+(length(it_seq))
}
if((nullit>=iter_adapt) ){# when to make the search area smaller
ff=ff+1
ff_scale <- ff*adapt_scale
nullit=1
}
if(trace && start_it==1){
cat(sprintf("Initial OFV = %g",res2[["ofv",1]]))
if(trace==2 | trace==3) cat(" | par = ",res2[["par",1]])
cat("\n")
}
if((trace && any((it_seq %% trace_iter)==0))){
if(length(it_seq)==1){
cat(sprintf(paste0("It. %",wd_iter,"i"),start_it))
} else {
cat(sprintf(paste0("It. %",wd_iter,"i to %",wd_iter,"i"),start_it,stop_it))
}
cat(sprintf(" | OFV = %g",ofv_opt))
if(trace==2) cat(" | opt. par. = ",par_opt)
#cat(" | runs = ",runs)
#cat(" | nullit = ",nullit)
if(trace==3) cat(" | best par tried = ",par)
cat("\n")
}
if(runs>=max_run){
cat(paste0("Maximum number of identical optimal values reached (max_run=",max_run,"), optimization stopped.\n"))
break
}
}
# Write results / end trace ---------
if((trace)){
cat("\nTotal iterations:",stop_it,"\n")
toc(name=".ars_savedTime")
cat("\nFinal OFV = ", ofv_opt, "\n")
cat("Parameters:",par_opt, "\n\n")
}
return(list( par=par_opt,ofv=ofv_opt))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.