PopED: Population (and Individual) Optimal Experimental Design

Documented in poped_optim_2

#' 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)) 
}