#' Optimization main module for PopED
#'
#' Optimize the objective function. The function works for both discrete and
#' continuous optimization variables. If more than one optimization method is
#' specified then the methods are run in series. If \code{loop_methods=TRUE}
#' then the series of optimization methods will be run for \code{iter_max}
#' iterations, or until the efficiency of the design after the current series
#' (compared to the start of the series) is less than \code{stop_crit_eff}.
#'
#' This function takes information from the PopED database supplied as an
#' argument. The PopED database supplies information about the the model,
#' parameters, design and methods to use. Some of the arguments coming from the
#' PopED database can be overwritten; if they are supplied then they are used
#' instead of the arguments from the PopED database.
#'
#' @inheritParams RS_opt
#' @inheritParams Doptim
#' @inheritParams create.poped.database
#' @inheritParams Dtrace
#' @inheritParams calc_ofv_and_fim
#' @inheritParams optim_LS
#' @param ... arguments passed to other functions.
#' @param control Contains control arguments for each method specified.
#' @param method A vector of optimization methods to use in a sequential
#' fashion. Options are \code{c("ARS","BFGS","LS","GA")}. \code{c("ARS")} is
#' for Adaptive Random Search \code{\link{optim_ARS}}. \code{c("LS")} is for
#' Line Search \code{\link{optim_LS}}. \code{c("BFGS")} is for Method
#' "L-BFGS-B" from \code{\link[stats]{optim}}. \code{c("GA")} is for the
#' genetic algorithm from \code{\link[GA]{ga}}.
#' @param out_file Save output from the optimization to a file.
#' @param loop_methods Should the optimization methods be looped for
#' \code{iter_max} iterations, or until the efficiency of the design after the
#' current series (compared to the start of the series) is less than, or equal to,
#' \code{stop_crit_eff}?
#' @param stop_crit_eff If \code{loop_methods==TRUE}, the looping will stop if the
#' efficiency of the design after the current series (compared to the start of
#' the series) is less than, or equal to, \code{stop_crit_eff} (if \code{maximize==FALSE} then 1/stop_crit_eff is the cut
#' off and the efficiency must be greater than or equal to this value to stop the looping).
#' @param stop_crit_diff If \code{loop_methods==TRUE}, the looping will stop if the
#' difference in criterion value of the design after the current series (compared to the start of
#' the series) is less than, or equal to, \code{stop_crit_diff} (if \code{maximize==FALSE} then -stop_crit_diff is the cut
#' off and the difference in criterion value must be greater than or equal to this value to stop the looping).
#' @param stop_crit_rel If \code{loop_methods==TRUE}, the looping will stop if the
#' relative difference in criterion value of the design after the current series (compared to the start of
#' the series) is less than, or equal to, \code{stop_crit_rel} (if \code{maximize==FALSE} then -stop_crit_rel is the cut
#' off and the relative difference in criterion value must be greater than or equal to this value to stop the looping).
#' @param maximize Should the objective function be maximized or minimized?
#' @param transform_parameters Should we transform the parameters before optimization?
#'
#'
#'
#' @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
#'
#' @keywords internal
# @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
# @example tests/testthat/examples_fcn_doc/examples_poped_optim.R
# @export
# uses create_ofv
poped_optim_2 <- function(poped.db,
opt_xt=poped.db$settings$optsw[2],
opt_a=poped.db$settings$optsw[4],
opt_x=poped.db$settings$optsw[3],
opt_samps=poped.db$settings$optsw[1],
opt_inds=poped.db$settings$optsw[5],
method=c("ARS","BFGS","LS"),
control=list(),
trace = TRUE,
fim.calc.type=poped.db$settings$iFIMCalculationType,
ofv_calc_type=poped.db$settings$ofv_calc_type,
approx_type=poped.db$settings$iApproximationMethod,
d_switch=poped.db$settings$d_switch,
ED_samp_size = poped.db$settings$ED_samp_size,
bLHS=poped.db$settings$bLHS,
use_laplace=poped.db$settings$iEDCalculationType,
out_file="",
parallel=F,
parallel_type=NULL,
num_cores = NULL,
loop_methods=ifelse(length(method)>1,TRUE,FALSE),
iter_max = 10,
stop_crit_eff = 1.001,
stop_crit_diff = NULL,
stop_crit_rel = NULL,
ofv_fun = poped.db$settings$ofv_fun,
maximize=T,
transform_parameters = F,
...){
#------------ update poped.db with options supplied in function
called_args <- match.call()
default_args <- formals()
for(i in names(called_args)[-1]){
if(length(grep("^poped\\.db\\$",capture.output(default_args[[i]])))==1) {
#eval(parse(text=paste(capture.output(default_args[[i]]),"<-",called_args[[i]])))
eval(parse(text=paste(capture.output(default_args[[i]]),"<-",i)))
}
}
#------------- initialization
fmf = 0 #The best FIM so far
dmf = 0 #The best ofv of FIM so far
#output <-calc_ofv_and_fim(poped.db,...)
#--------------------- Header information
## only for header information
if(is.null(ofv_fun) || is.function(ofv_fun)){
ofv_fun_user <- ofv_fun
} else {
# source explicit file
# here I assume that function in file has same name as filename minus .txt and pathnames
if(file.exists(as.character(ofv_fun))){
source(as.character(ofv_fun))
ofv_fun_user <- eval(parse(text=fileparts(ofv_fun)[["filename"]]))
} else {
stop("ofv_fun is not a function or NULL, and no file with that name was found")
}
}
if(!is.null(ofv_fun)){
poped.db$settings$ofv_calc_type = 0
}
output <-calc_ofv_and_fim(poped.db,d_switch=d_switch,
ED_samp_size=ED_samp_size,
bLHS=bLHS,
use_laplace=use_laplace,
ofv_calc_type=ofv_calc_type,
fim.calc.type=fim.calc.type,
ofv_fun = ofv_fun_user,
...)
fmf <- output$fim
dmf <- output$ofv
fmf_init <- fmf
dmf_init <- dmf
poped.db_init <- poped.db
if(is.nan(dmf_init)) stop("Objective function of initial design is NaN")
fn=blockheader(poped.db,name="optim",e_flag=!d_switch,
fmf=fmf_init,dmf=dmf_init,
out_file=out_file,
trflag=trace,
...)
#---------------------- ofv for optimization
my_ofv <- create_ofv(poped.db=poped.db,
opt_xt=opt_xt,
opt_a=opt_a,
opt_x=opt_x,
opt_samps=opt_samps,
opt_inds=opt_inds,
fim.calc.type=fim.calc.type,
ofv_calc_type=ofv_calc_type,
approx_type=approx_type,
d_switch=d_switch,
ED_samp_size = ED_samp_size,
bLHS=bLHS,
use_laplace=use_laplace,
ofv_fun = ofv_fun,
transform_parameters = transform_parameters,
...)
par <- my_ofv$par
ofv_fun <- my_ofv$fun
back_transform_par_fun <- my_ofv$back_transform_par
lower=my_ofv$space[["lower"]]
upper=my_ofv$space[["upper"]]
allowed_values=my_ofv$space[["allowed_values"]]
par_cat_cont=my_ofv$space[["par_cat_cont"]]
par_fixed_index=my_ofv$space[["par_fixed_index"]]
par_not_fixed_index=my_ofv$space[["par_not_fixed_index"]]
par_df_unique=my_ofv$space[["par_df_unique"]]
par_df=my_ofv$space[["par_df"]]
par_dim=my_ofv$space[["par_dim"]]
par_type=my_ofv$space[["par_type"]]
#------------ optimize
if(!(fn=="")) sink(fn, append=TRUE, split=TRUE)
iter <- 0
stop_crit <- FALSE
output$ofv <- dmf_init
while(stop_crit==FALSE && iter < iter_max){
ofv_init <- output$ofv
iter=iter+1
method_loop <- method
if(loop_methods){
cat("************* Iteration",iter," for all optimization methods***********************\n\n")
}
while(length(method_loop)>0){
cur_meth <- method_loop[1]
method_loop <- method_loop[-1]
if(cur_meth=="ARS"){
cat("*******************************************\n")
cat("Running Adaptive Random Search Optimization\n")
cat("*******************************************\n")
# handle control arguments
con <- list(trace = trace,
parallel=parallel,
parallel_type=parallel_type,
num_cores = num_cores)
nmsC <- names(con)
con[(namc <- names(control$ARS))] <- control$ARS
#if (length(noNms <- namc[!namc %in% nmsC])) warning("unknown names in control: ", paste(noNms, collapse = ", "))
output <- do.call(optim_ARS,c(list(par=par,
fn=ofv_fun,
lower=lower,
upper=upper,
allowed_values = allowed_values,
maximize=maximize
#par_df_full=par_df
),
#par_grouping=par_grouping),
con,
...))
}
if(cur_meth=="LS"){
cat("*******************************************\n")
cat("Running Line Search Optimization\n")
cat("*******************************************\n")
# handle control arguments
con <- list(trace = trace,
parallel=parallel,
parallel_type=parallel_type,
num_cores = num_cores)
nmsC <- names(con)
con[(namc <- names(control$LS))] <- control$LS
#if (length(noNms <- namc[!namc %in% nmsC])) warning("unknown names in control: ", paste(noNms, collapse = ", "))
output <- do.call(optim_LS,c(list(par=par,
fn=ofv_fun,
lower=lower,
upper=upper,
allowed_values = allowed_values,
maximize=maximize
#par_df_full=par_df
),
#par_grouping=par_grouping),
con,
...))
}
if(cur_meth=="BFGS"){
cat("*******************************************\n")
cat("Running L-BFGS-B Optimization\n")
cat("*******************************************\n")
if(all(par_cat_cont=="cat")){
cat("\nNo continuous variables to optimize, L-BFGS-B Optimization skipped\n\n")
next
}
if(trace) trace_optim=3
if(is.numeric(trace)) trace_optim = trace
#if(trace==2) trace_optim = 4
#if(trace==3) trace_optim = 5
#if(trace==4) trace_optim = 6
# handle control arguments
con <- list(trace=trace_optim)
nmsC <- names(con)
con[(namc <- names(control$BFGS))] <- control$BFGS
fnscale=-1
if(!maximize) fnscale=1
if(is.null(con[["fnscale"]])) con$fnscale <- fnscale
#if (length(noNms <- namc[!namc %in% nmsC])) warning("unknown names in control: ", paste(noNms, collapse = ", "))
par_full <- par
output <- optim(par=par[par_cat_cont=="cont"],
fn=ofv_fun,
gr=NULL,
#par_full=par_full,
only_cont=T,
lower=lower[par_cat_cont=="cont"],
upper=upper[par_cat_cont=="cont"],
method = "L-BFGS-B",
control=con)
output$ofv <- output$value
par_tmp <- output$par
output$par <- par_full
output$par[par_cat_cont=="cont"] <- par_tmp
fprintf('\n')
if(fn!="") fprintf(fn,'\n')
}
if(cur_meth=="GA"){
cat("*******************************************\n")
cat("Running Genetic Algorithm (GA) Optimization\n")
cat("*******************************************\n")
if (!requireNamespace("GA", quietly = TRUE)) {
stop("GA package needed for this function to work. Please install it.",
call. = FALSE)
}
if(all(par_cat_cont=="cat")){
cat("\nNo continuous variables to optimize, GA Optimization skipped\n\n")
next
}
# handle control arguments
parallel_ga <- parallel
if(!is.null(num_cores)) parallel_ga <- num_cores
if(!is.null(parallel_type)) parallel_ga <- parallel_type
con <- list(parallel=parallel_ga)
dot_vals <- dots(...)
#if(is.null(dot_vals[["monitor"]]) && packageVersion("GA")>="3.0.2") con$monitor <- GA::gaMonitor2
nmsC <- names(con)
con[(namc <- names(control$GA))] <- control$GA
#if (length(noNms <- namc[!namc %in% nmsC])) warning("unknown names in control: ", paste(noNms, collapse = ", "))
par_full <- par
ofv_fun_2 <- ofv_fun
if(!maximize) {
ofv_fun_2 <- function(par,only_cont=F,...){
-ofv_fun(par,only_cont=F,...)
}
}
output_ga <- do.call(GA::ga,c(list(type = "real-valued",
fitness = ofv_fun_2,
#par_full=par_full,
only_cont=T,
min=lower[par_cat_cont=="cont"],
max=upper[par_cat_cont=="cont"],
suggestions=par[par_cat_cont=="cont"]),
#allowed_values = allowed_values),
con,
...))
output$ofv <- output_ga@fitnessValue
if(!maximize) output$ofv <- -output$ofv
output$par <- output_ga@solution
par_tmp <- output$par
output$par <- par_full
output$par[par_cat_cont=="cont"] <- par_tmp
fprintf('\n')
if(fn!="") fprintf(fn,'\n')
}
par <- output$par
}
if(!loop_methods){
stop_crit <- TRUE
} else {
cat("*******************************************\n")
cat("Stopping criteria testing\n")
cat("(Compare between start of iteration and end of iteration)\n")
cat("*******************************************\n")
# relative difference
rel_diff <- (output$ofv - ofv_init)/abs(ofv_init)
abs_diff <- (output$ofv - ofv_init)
fprintf("Difference in OFV: %.3g\n",abs_diff)
fprintf("Relative difference in OFV: %.3g%%\n",rel_diff*100)
# efficiency
eff <- efficiency(ofv_init, output$ofv, poped.db)
fprintf("Efficiency: \n (%s) = %.5g\n",attr(eff,"description"),eff)
#cat("Efficiency: \n ", attr(eff,"description"), sprintf("%.5g",eff), "\n")
#if(eff<=stop_crit_eff) stop_crit <- TRUE
#cat("eff: ",sprintf("%.3g",(output$ofv - ofv_init)/p), "\n")
#cat("eff: ",sprintf("%.3g",(exp(output$ofv)/exp(ofv_init))^(1/p)), "\n")
compare <-function(crit,crit_stop,maximize,inv=FALSE,neg=FALSE,text=""){
if(is.null(crit_stop)){
#cat(" Stopping criteria not defined\n")
return(FALSE)
}
cat("\n",text,"\n")
if(is.nan(crit)){
fprintf(" Stopping criteria using 'NaN' as a comparitor cannot be used\n")
return(FALSE)
}
comparitor <- "<="
if(!maximize) comparitor <- ">="
if(inv) crit_stop <- 1/crit_stop
if(neg) crit_stop <- -crit_stop
fprintf(" Is (%0.5g %s %0.5g)? ",crit,comparitor,crit_stop)
res <- do.call(comparitor,list(crit,crit_stop))
if(res) cat(" Yes.\n Stopping criteria achieved.\n")
if(!res) cat(" No.\n Stopping criteria NOT achieved.\n")
return(res)
#if(maximize) cat("Efficiency stopping criteria (lower limit) = ",crit_stop, "\n")
#if(!maximize) cat("Efficiency stopping criteria (upper limit) = ",1/crit_stop, "\n")
#if(maximize) return(crit <= crit_stop)
#if(!maximize) return(crit >= 1/crit_stop)
}
if(all(is.null(c(stop_crit_eff,stop_crit_rel,stop_crit_diff)))){
cat("No stopping criteria defined")
} else {
stop_eff <- compare(eff,stop_crit_eff,maximize,inv=!maximize,
text="Efficiency stopping criteria:")
stop_abs <- compare(abs_diff,stop_crit_diff,maximize,neg=!maximize,
text="OFV difference stopping criteria:")
stop_rel <- compare(rel_diff,stop_crit_rel,maximize,neg=!maximize,
text="Relative OFV difference stopping criteria:")
if(stop_eff || stop_rel || stop_abs) stop_crit <- TRUE
if(stop_crit){
cat("\nStopping criteria achieved.\n")
} else {
cat("\nStopping criteria NOT achieved.\n")
}
cat("\n")
}
}
} # end of total loop
if(!(fn=="")) sink()
# back transform parameters
par <- back_transform_par_fun(par)
# add the results into a poped database
# expand results to full size
if(length(par_fixed_index)!=0){
par_df_2[par_not_fixed_index,"par"] <- par
par <- par_df_2$par
}
if(!is.null(par_df_unique)){
par_df_unique$par <- par
for(j in par_df_unique$par_grouping){
par_df[par_df$par_grouping==j,"par"] <- par_df_unique[par_df_unique$par_grouping==j,"par"]
}
} else {
par_df$par <- par
}
#poped.db$design$ni <- ni
#if(opt_xt) poped.db$design$xt[,]=matrix(par_df[par_type=="xt","par"],par_dim$xt)
if(opt_xt){
xt <- zeros(par_dim$xt)
par_xt <- par_df[par_type=="xt","par"]
for(i in 1:poped.db$design$m){
if((poped.db$design$ni[i]!=0 && poped.db$design$groupsize[i]!=0)){
xt[i,1:poped.db$design$ni[i]] <- par_xt[1:poped.db$design$ni[i]]
par_xt <- par_xt[-c(1:poped.db$design$ni[i])]
}
}
poped.db$design$xt[,]=xt[,]
}
if(opt_a) poped.db$design$a[,]=matrix(par_df[par_type=="a","par"],par_dim$a)
# if((!isempty(x))){
# poped.db$design$x[1:size(x,1),1:size(x,2)]=x
# #poped.db$design$x <- x
# }
# if((!isempty(a))){
# poped.db$design$a[1:size(a,1),1:size(a,2)]=a
# #poped.db$design$a <- a
# }
#--------- Write results
#if((trflag)){
# if(footer_flag){
#FIM <- evaluate.fim(poped.db,...)
# if(d_switch){
# FIM <- evaluate.fim(poped.db,...)
# } else{
# out <-calc_ofv_and_fim(poped.db,d_switch=d_switch,
# ED_samp_size=ED_samp_size,
# bLHS=bLHS,
# use_laplace=use_laplace,
# ofv_calc_type=ofv_calc_type,
# fim.calc.type=fim.calc.type,
# ...)
#
# FIM <- out$fim
# }
FIM <-calc_ofv_and_fim(poped.db,
ofv=output$ofv,
fim=0,
d_switch=d_switch,
ED_samp_size=ED_samp_size,
bLHS=bLHS,
use_laplace=use_laplace,
ofv_calc_type=ofv_calc_type,
fim.calc.type=fim.calc.type,
ofv_fun = ofv_fun_user,
...)[["fim"]]
time_value <-
blockfinal(fn=fn,fmf=FIM,
dmf=output$ofv,
groupsize=poped.db$design$groupsize,
ni=poped.db$design$ni,
xt=poped.db$design$xt,
x=poped.db$design$x,
a=poped.db$design$a,
model_switch=poped.db$design$model_switch,
poped.db$parameters$param.pt.val$bpop,
poped.db$parameters$param.pt.val$d,
poped.db$parameters$docc,
poped.db$parameters$param.pt.val$sigma,
poped.db,
fmf_init=fmf_init,
dmf_init=dmf_init,
...)
# }
#}
results <- list( ofv= output$ofv, FIM=FIM, initial=list(ofv=dmf_init, FIM=fmf_init, poped.db=poped.db_init),
run_time=time_value, poped.db = poped.db )
class(results) <- "poped_optim"
return(invisible(results))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.