#' Optimize a function using a line search algorithm.
#'
#' \code{optim_LS} performs sequential grid search optimization of an arbitrary function with respect
#' to each of the parameters to be optimized over.
#' The function works for both discrete and continuous optimization parameters
#' and allows for box-constraints (by using the upper and lower function arguments) or sets of allowed values (by using the
#' allowed_values function argument) for all parameters, or on a parameter per parameter basis.
#'
#' @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 par A vector of initial values for the parameters to be optimized over.
#' @param fn A function to be minimized (or maximized),
#' with first argument the vector of parameters over which minimization is to take place.
#' It should return a scalar result.
#' @param upper Upper bounds on the parameters. A vector.
#' @param lower Lower bounds on the parameters. A vector.
#' @param line_length The number of different parameter values per parameter to evaluate. The values are
#' selected as an evenly spaced grid between the upper and lower bounds.
#' @param trace Should the algorithm output results intermittently.
#' @param maximize Should the function be maximized? Default is to minimize.
#' @param parallel Should we use parallel computations?
#' @param parallel_type Which type of parallelization should be used?
#' Can be "snow" or "multicore". "snow" works on Linux-like systems & Windows. "multicore" works only on
#' Linux-like systems. By default this is chosen for you depending on your operating system.
#' See \code{\link{start_parallel}}.
#' @param num_cores The number of cores to use in the parallelization. By default is set to the number
#' output from
#' \code{parallel::detectCores()}.
#' See \code{\link{start_parallel}}.
#' @param mrgsolve_model If the computations require a mrgsolve model and you
#' are using the "snow" method then you need to specify the name of the model
#' object created by \code{mread} or \code{mcode}.
#' @param seed The random seed to use in the algorithm,
#' @param allow_replicates Should the algorithm allow parameters to have the same value?
#' @param replicates_index A vector, the same length as the parameters.
#' If two values are the same in this vector then the parameters may not assume the same value in the optimization.
#' @param ofv_initial An initial objective function value (OFV). If not NULL then the initial design is not evaluated
#' and the OFV value is assumed to be this number.
#' @param closed_bounds Are the upper and lower limits open (boundaries not allowed) or closed (boundaries allowed) bounds?
#' @param ... Additional arguments passed to \code{fn} and \code{start_parallel}.
#'
#' @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
#'
#' @example tests/testthat/examples_fcn_doc/examples_optim_LS.R
#' @export
optim_LS <- function(par,
fn,
lower=NULL,
upper=NULL,
allowed_values=NULL,
#constraints=NULL,
#control=NULL,
#loc_fac=4,
#no_bounds_sd = par,
line_length=50,
trace = TRUE,
maximize=F, # default is to minimize
parallel=F, # T or F or a list of (type, n_cores)
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
ofv_initial=NULL,
closed_bounds=TRUE, # are upper and lower limits allowed values?
...){
# parallel with # T or F or a list of (type, n_cores)
# also change with optim_ARS and poped_optim and maybe change to optim_poped
#---------------- start trace
if((trace)){
tic(name=".ls_savedTime")
}
#--------------- checks
if((is.null(lower) || is.null(upper)) && is.null(allowed_values)){
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
itvector <- c()
dmfvector <- c()
counter = 0
best_changed = FALSE
# handle replicates
if(!allow_replicates){
replicates_index <- rep(1,length(par))
}
fn1 <- function(par) fn(par, ...)
par_opt <- par
if(is.null(ofv_initial)){
ofv_opt <- fn1(par_opt)
} else {
ofv_opt <- ofv_initial
}
if(!is.null(seed)) set.seed(seed)
# start parallel computing
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")))
}
#if(is.null(iter_chunk)) if(parallel) iter_chunk <- attr(parallel,"cores") else iter_chunk <- 1
# continuous and discrete parameters
par_type <- rep("cont",length(par))
if(!is.null(allowed_values)){
for(k in 1:length(par)){
if(!all(is.na(allowed_values[[k]])) && length(allowed_values[[k]]>0)){
par_type[k] <- "cat"
}
}
}
#resample <- function(x,size=1) x[sample.int(length(x),size=size)]
#if(!is.null(par)) iter = iter+1
par_tried <- list(par) # investigated parameter vectors
compare <-function(a,b) a < b
which_compare <- function(x) which.min(x)
if(maximize){
compare <-function(a,b) a > b
which_compare <- function(x) which.max(x)
}
# ------------ generate new parameter sets
gen_par_list <- function (par,i, ...) {
# generate new continuous parameters
if(par_type[i]=="cont"){
if(!is.finite(lower[i])) lower[i]=par[i]-1000*par[i]
if(!is.finite(upper[i])) upper[i]=par[i]+1000*par[i]
if(closed_bounds){
par_i_set <- seq(lower[i],upper[i],len=line_length)
} else {
par_i_set <- seq(lower[i],upper[i],len=(line_length+2))
par_i_set <- par_i_set[-c(1,length(par_i_set))]
}
} else {
par_i_set <- allowed_values[[i]]
}
# 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]])]
}
# create full parameter list
par_list <- lapply(par_i_set,function(x){ par[i] <- x; par})
# remove any design that has already been tested (works but is slow!)
#par_list <- par_list[unlist(lapply(par_list,function(x){!(list(x) %in% par_tried)}))]
return(par_list)
} # end function
# Randomize which parameter to investigate first
par_order <- sample.int(length(par))
#------------- evaluate new parameter sets
if(trace) cat("\n Initial parameters:",par_opt, "\n")
if(trace) cat(" Initial OFV:",ofv_opt, "\n\n")
for(i in par_order){
par_list <- gen_par_list(par_opt,i)
par_tried <- c(par_tried,par_list)
if(trace){
cat(" Searching parameter",i,"\n")
#pb <- txtProgressBar(min=0, max=n, initial=1, style=3)
#setTxtProgressBar(pb, i)
#close(pb)
}
if(length(par_list)==0){
if(trace) cat(paste(" No change","\n"))
next
}
if(parallel && (attr(parallel,"type")=="multicore")){
if(!is.null(seed)) set.seed(seed+i)
res <- parallel::mclapply(par_list,fn1,mc.cores=attr(parallel, "cores"))
} else if(parallel && (attr(parallel,"type")=="snow")){
res <- parallel::parLapply(attr(parallel, "cluster"),par_list,fn1)
} else {
res <- lapply(par_list,fn1)
}
res <- unlist(res)
best_index <- which_compare(res)
ofv <- res[best_index]
par <- par_list[[best_index]]
if((compare(ofv,ofv_opt) || is.null(ofv_opt))){
if(trace) cat(sprintf(" Changed from %g to %g",par_opt[i],par[i]))
if(trace) cat(sprintf(" ; OFV = %g \n",ofv))
par_opt <- par
ofv_opt <- ofv
} else {
if(trace) cat(paste(" No change","\n"))
}
}
#--------- Write results
if((trace)){
cat("\n ")
toc(name=".ls_savedTime")
cat("\n Final 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.