R/Doptim.R

In PopED: Population (and Individual) Optimal Experimental Design

#' D-family optimization function
#' 
#' Optimize the objective function. There are 4 different optimization
#' algorithms used in this function \enumerate{ \item Adaptive random search.
#' See \code{\link{RS_opt}}. \item Stochastic gradient. \item A Broyden Fletcher
#' Goldfarb Shanno (BFGS) method for nonlinear minimization with box
#' constraints. \item A line search. See \code{\link{a_line_search}}. } The
#' optimization algorithms run in series, taking as input the output from the
#' previous method. The stopping rule used is to test if the line search
#' algorithm fids a better optimum then its initial value. If so, then the chain
#' of algorithms is run again.  If line search is not used then the argument
#' \code{iter_tot} defines the number of times the chain of algorithms is run. 
#' 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 evaluate.fim
#' @inheritParams create.poped.database
#' @param bpopdescr Matrix defining the fixed effects, per row (row number =
#'   parameter_number) we should have: \itemize{ \item column 1 the type of the
#'   distribution for E-family designs (0 = Fixed, 1 = Normal, 2 = Uniform, 3 =
#'   User Defined Distribution, 4 = lognormal and 5 = truncated normal) \item
#'   column 2  defines the mean. \item column 3 defines the variance of the
#'   distribution (or length of uniform distribution). }
#' @param ddescr Matrix defining the diagonals of the IIV (same logic as for
#'   the \code{bpopdescr}).
#' @param fmf The initial value of the FIM. If set to zero then it is computed.
#' @param dmf The initial OFV. If set to zero then it is computed.
#' @param trflag Should the optimization be output to the screen and to a file?
#' @param ls_step_size Number of grid points in the line search.
#' @param iter_tot Number of iterations to use if line search is not used. Must
#'   be less than \code{iter_max} to be used.
#' @param iter_max If line search is used then the algorithm tests if line
#'   search (always run at the end of the optimization iteration) changes the 
#'   design in any way.  If not, the algorithm stops.  If yes, then a new
#'   iteration is run unless \code{iter_max} iterations have already been run.
#'   
#' @references \enumerate{ \item M. Foracchia, A.C. Hooker, P. Vicini and A.
#'   Ruggeri, "PopED, a software for 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/warfarin_optimize.R
#' @example tests/testthat/examples_fcn_doc/examples_Doptim.R
#' @export
#' @keywords internal
#' 
## Function translated using 'matlab.to.r()'
## Then manually adjusted to make work
## Author: Andrew Hooker

Doptim <- function(poped.db,ni, xt, model_switch, x, a, bpopdescr, 
                   ddescr, maxxt, minxt,maxa,mina,fmf=0,dmf=0,
                   trflag=TRUE,
                   bUseRandomSearch=poped.db$settings$bUseRandomSearch,
                   bUseStochasticGradient=poped.db$settings$bUseStochasticGradient,
                   bUseBFGSMinimizer=poped.db$settings$bUseBFGSMinimizer,
                   bUseLineSearch=poped.db$settings$bUseLineSearch,
                   sgit=poped.db$settings$sgit,
                   ls_step_size=poped.db$settings$ls_step_size,
                   BFGSConvergenceCriteriaMinStep=poped.db$settings$BFGSConvergenceCriteriaMinStep,
                   BFGSProjectedGradientTol=poped.db$settings$BFGSProjectedGradientTol,
                   BFGSTolerancef=poped.db$settings$BFGSTolerancef,
                   BFGSToleranceg=poped.db$settings$BFGSToleranceg,
                   BFGSTolerancex=poped.db$settings$BFGSTolerancex,
                   iter_tot=poped.db$settings$iNumSearchIterationsIfNotLineSearch,
                   iter_max=10,
                   ...){
  
  
  ## 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)))    }
  }
  iter=0 #Search iterations
  test_change=TRUE
  
  trflag = trflag
  
  # ----------------- type of optimization determination
  axt=poped.db$settings$optsw[2]*poped.db$settings$cfaxt*matrix(1,poped.db$design$m,max(poped.db$design_space$maxni))
  aa=poped.db$settings$optsw[4]*poped.db$settings$cfaa*matrix(1,poped.db$design$m,size(poped.db$design$a,2))
  optxt=poped.db$settings$optsw[2]
  optx=poped.db$settings$optsw[3]
  opta=poped.db$settings$optsw[4]
  bfgs_init=matrix(0,0,0)
  # ----------------- initialization of size variables
  m=size(ni,1)
  maxni=size(xt,2)
  
  iMaxSearchIterations = iter_tot
  
  # ----------------- initialization of model parameters
  bpop=bpopdescr[,2,drop=F]
  d=getfulld(ddescr[,2,drop=F],poped.db$parameters$covd)
  docc_full = getfulld(poped.db$parameters$docc[,2,drop=F],poped.db$parameters$covdocc)
  
  Engine = list(Type=1,Version=version$version.string)
  
  if((dmf==0) ){#Only first time
    returnArgs <-  mftot(model_switch,poped.db$design$groupsize,ni,xt,x,a,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
    fmf <- returnArgs[[1]]
    poped.db <- returnArgs[[2]]
    dmf=ofv_fim(fmf,poped.db)
    fmf_init <- fmf
    dmf_init <- dmf
  }
  
  if((optxt || optx || opta)){
    # ------------------ Write output file header
    #fn=blockheader(FALSE,iter,poped.db)
    #if((trflag)){
    fn=blockheader(poped.db,name='D_cont_opt',iter=1,
                     opt_xt=optxt,opt_a=opta,opt_x=optx,
                     opt_inds=F,opt_samps=F,
                     fmf=fmf_init,dmf=dmf_init,
                     bpop=bpopdescr,d=ddescr,docc=poped.db$parameters$docc,sigma=poped.db$parameters$sigma,
                     trflag=trflag,
                     ...)
    #       param_vars_init=diag_matlab(inv(fmf))
    #       returnArgs <-  get_cv(param_vars_init,bpop=bpopdescr,d=ddescr,docc=poped.db$parameters$docc,sigma=poped.db$parameters$sigma,poped.db) 
    #       params_init <- returnArgs[[1]]
    #       param_cvs_init <- returnArgs[[2]]
    #}
    while(test_change==TRUE && iMaxSearchIterations>0 && iter < iter_max){
      iter=iter+1
      
      
      
      
      if(((bUseRandomSearch || bUseStochasticGradient) && poped.db$settings$bShowGraphs)){
        #figure(1)
        #if((poped.db$settings$Engine$Type==1)){
        #    set(1,'Name','Random Search and Stochastic Gradient')
        #}
      }
      
      itvector <- c()
      dmfvector <- c()
      
      # ----------------- initialization of optimum variables (and old)
      xtopt=xt
      xopt=x
      aopt=a
      xtopto=xt
      xopto=x
      aopto=a
      lgxto=zeros(m*maxni)
      lgao=zeros(m*size(poped.db$design$a,2))
      
      # ----------------- RS variables initialization
      dxt=(maxxt-minxt)/poped.db$settings$rslxt
      da=(maxa-mina)/poped.db$settings$rsla
      xtoptn=xtopt
      xoptn=xopt
      aoptn=aopt
      nullit=1
      ff=1
      
      # ----------------- SG variables initialization
      kitxt=1
      kita=1
      mnormxt=zeros(m,maxni)
      mnorma=zeros(m,size(poped.db$design$a,2))
      normgxt=zeros(m,maxni)
      normga=zeros(m,size(poped.db$design$a,2))
      
      # ----------------- initialization of trace support variables
      odmf=0
      inversionxt=FALSE
      inversiona=FALSE
      
      if((bUseRandomSearch) ){#If we want to perform random search
        # ----------------- RANDOM SEARCH BEGINS HERE
        #save for graphical output
        if((poped.db$settings$parallel$bParallelLS == 0)){
          tic()              
        }
        itvector[1] = 0
        dmfvector[1] = dmf
        
        if((trflag)){
          Dtrace(fn,0,ni,xtopt,xopt,aopt,matrix(0,0,0),matrix(0,0,0),dmf,matrix(0,0,0),matrix(0,0,0),matrix(0,0,0),itvector,dmfvector,poped.db)
        }
        
        if((poped.db$settings$parallel$bParallelRS)){
          
          # Generate the input designs
          
          designsin = cell(1,0)
          for(it in 1:poped.db$settings$rsit){
            if((optxt==TRUE)){
              if((poped.db$design_space$bUseGrouped_xt)){
                xtoptn=grouped_rand(poped.db$design_space$G_xt,xtopt,dxt,ff,axt)
              } else {
                xtoptn=xtopt+dxt/ff*randn(m,maxni)*(axt>0)
              }
              xtoptn=xtoptn-((xtoptn>maxxt)*(xtoptn-maxxt))
              xtoptn=xtoptn+((xtoptn<minxt)*(minxt-xtoptn))
            }
            if((optx==TRUE)){
              xoptn=get_discrete_x(poped.db$design_space$G_x,poped.db$design_space$discrete_x,poped.db$design_space$bUseGrouped_x)
            }
            if((opta==TRUE)){
              if((poped.db$design_space$bUseGrouped_a)){
                aoptn=grouped_rand_a(poped.db$design_space$G_a,aopt,da,ff,aa)
              } else {
                aoptn=aopt+da/ff*randn(m,size(poped.db$design$a,2))*(aa>0)
              }
              aoptn=aoptn-((aoptn>maxa)*(aoptn-maxa))
              aoptn=aoptn+((aoptn<mina)*(mina-aoptn))
            }
            designsin = update_designinlist(designsin,poped.db$design$groupsize,ni,xtoptn,xoptn,aoptn,-1,0)
          }
          
          stop("Parallel execution not yet implemented in R")
          designsout = designsin
          #designsout = execute_parallel(designsin,poped.db)
          #Store the optimal design
          for(it in 1:poped.db$settings$rsit){
            if((designsout[[it]]$ofv>dmf) ){##ok<USENS>
              if((optxt==TRUE)){
                xtopt=designsin[[it]]$xt
              }
              if((optx==TRUE)){
                xopt=designsin[[it]]$x
              }
              if((opta==TRUE)){
                aopt=designsin[[it]]$a
              }
              dmf=designsout[[it]]$ofv
              fmf=designsout[[it]]$FIM
            }
            
            if((trflag && ((it %% poped.db$settings$rsit_output)==0 || it==poped.db$settings$rsit))){
              itvector[ceiling(it/poped.db$settings$rsit_output)+1]=it
              dmfvector[ceiling(it/poped.db$settings$rsit_output)+1]=dmf
              Dtrace(fn,it,ni,xtopt,xopt,aopt,matrix(0,0,0),matrix(0,0,0),dmf,matrix(0,0,0),matrix(0,0,0),matrix(0,0,0),itvector,dmfvector,poped.db)
            }
          }
        } else {
          for(it in 1:poped.db$settings$rsit){
            
            if((optxt==TRUE)){
              if((poped.db$design_space$bUseGrouped_xt)){
                xtoptn=grouped_rand(poped.db$design_space$G_xt,xtopt,dxt,ff,axt)
              } else {
                xtoptn=xtopt+dxt/ff*randn(m,maxni)*(axt>0)
              }
              xtoptn=xtoptn-((xtoptn>maxxt)*(xtoptn-maxxt))
              xtoptn=xtoptn+((xtoptn<minxt)*(minxt-xtoptn))
            }
            if((optx==TRUE)){
              xoptn=get_discrete_x(poped.db$design_space$G_x,poped.db$design_space$discrete_x,poped.db$design_space$bUseGrouped_x)
            }
            if((opta==TRUE)){
              if((poped.db$design_space$bUseGrouped_a)){
                aoptn=grouped_rand_a(poped.db$design_space$G_a,aopt,da,ff,aa)
              } else {
                aoptn=aopt+da/ff*randn(m,size(poped.db$design$a,2))*(aa>0)
              }
              aoptn=aoptn-((aoptn>maxa)*(aoptn-maxa))
              aoptn=aoptn+((aoptn<mina)*(mina-aoptn))
            }
            returnArgs <- mftot(model_switch,poped.db$design$groupsize,ni,xtoptn,xoptn,aoptn,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
            nfmf <- returnArgs[[1]]
            poped.db <- returnArgs[[2]]
            ndmf=ofv_fim(nfmf,poped.db)
            if((ndmf>dmf)){
              if((optxt==TRUE)){
                xtopt=xtoptn
              }
              if((optx==TRUE)){
                xopt=xoptn
              }
              if((opta==TRUE)){
                aopt=aoptn
              }
              dmf=ndmf
              fmf=nfmf
              nullit=1
              ff=1
            } else {
              nullit=nullit+1
            }
            if((nullit==poped.db$settings$maxrsnullit)){
              ff=ff+1
              nullit=1
            }
            
            if((!isempty(poped.db$settings$strIterationFileName))){
              write_iterationfile('Random Search',it,xtopt,aopt,xopt,ni,poped.db$design$groupsize,fmf,dmf,poped.db)
            }
            
            if((trflag && ((it %% poped.db$settings$rsit_output)==0 || it==poped.db$settings$rsit))){
              itvector[ceiling(it/poped.db$settings$rsit_output)+1]=it
              dmfvector[ceiling(it/poped.db$settings$rsit_output)+1]=dmf
              Dtrace(fn,it,ni,xtopt,xopt,aopt,matrix(0,0,0),matrix(0,0,0),dmf,matrix(0,0,0),matrix(0,0,0),matrix(0,0,0),itvector,dmfvector,poped.db)
            }
          }  
        }
        
        if((poped.db$settings$parallel$bParallelLS == 0)){
          timeLS = toc(echo=FALSE)
          fprintf('Run time for random search: %g seconds\n\n',timeLS)
          #if(trflag) fprintf(fn,'Elapsed time for Serial Random search with %g iterations: %g seconds\n',it,timeLS)
        }
        # ----------------- RANDOM SEARCH ENDS HERE 
      }
      # ----------------- initialization of best optimum variables
      
      bestxt=xtopt
      bestx=xopt
      besta=aopt
      
      diff=poped.db$settings$convergence_eps+1
      
      if((bUseBFGSMinimizer )){
        if((Engine$Type==2)){
          stop(sprintf('BFGS optimization can not be used with Freemat in this version!'))
        }
        
        f_name <- 'calc_ofv_and_grad' 
        f_options <- list(x,optxt, opta, model_switch,aa,axt,poped.db$design$groupsize,ni,xtopt,xopt,aopt,bpop,d,poped.db$parameters$sigma,docc_full,poped.db)
        x_k=matrix(0,0,0)
        lb=matrix(0,0,0)
        ub=matrix(0,0,0)
        options=list('factr'=BFGSConvergenceCriteriaMinStep,'pgtol'=BFGSProjectedGradientTol,'ftol'=BFGSTolerancef,'gtol'=BFGSToleranceg,'xtol'=BFGSTolerancex)
        if((optxt==TRUE)){
          index=t(1:numel(xtopt))
          if(poped.db$design_space$bUseGrouped_xt){
            returnArgs <- unique(poped.db$design_space$G_xt) 
            temp <- returnArgs[[1]]
            index <- returnArgs[[2]]
            temp2 <- returnArgs[[3]]
          }
          index=index[minxt!=maxxt]
          x_k=t(t(xtopt[index]))
          lb=t(t(minxt[index]))
          ub=t(t(maxxt[index]))
        }
        if((opta==TRUE)){
          index=t(1:numel(aopt))
          if(poped.db$design_space$bUseGrouped_a){
            returnArgs <- unique(poped.db$design_space$G_a) 
            temp1 <- returnArgs[[1]]
            index <- returnArgs[[2]]
            temp2 <- returnArgs[[3]]
          }
          index=index[mina!=maxa]
          x_k=t(t(c(x_k,aopt[index])))
          lb=t(t(c(lb,mina[index])))
          ub=t(t(c(ub,maxa[index])))
          
          #x_k(end+index)=aopt[index]
          #lb(end+index)=mina[index]
          #ub(end+index)=maxa[index]
        }
        if((any(x_k<lb))){
          x_k[x_k<lb]=lb[x_k<lb]
        }
        if((isempty(bfgs_init) || any(x_k!=bfgs_init))){
          bfgs_init=x_k
          fprintf('Starting BGFS minimization with OFV of %g \n', dmf)
          returnArgs <- bfgsb_min(f_name,f_options, x_k,lb,ub,options) 
          x_opt  <- returnArgs[[1]]
          f_k <- returnArgs[[2]]
          B_k <- returnArgs[[3]]
          
          if(optxt){
            notfixed=minxt!=maxxt
            if(poped.db$design_space$bUseGrouped_xt){
              xtopt[notfixed]=x_opt[poped.db$design_space$G_xt[notfixed]]
              x_opt <- x_opt[-(1:numel(unique(poped.db$design_space$G_xt[notfixed])))]
            } else {
              xtopt[notfixed]=x_opt[1:numel(xtopt[notfixed])]
              x_opt <- x_opt[-(1:numel(xtopt[notfixed]))]
            }
          }
          
          if(opta){
            notfixed=mina!=maxa
            
            if(poped.db$design_space$bUseGrouped_a){
              aopt[notfixed]=x_opt(poped.db$design_space$G_a[notfixed])
            } else {
              aopt[notfixed]=x_opt
            }
          }
          returnArgs <-  mftot(model_switch,poped.db$design$groupsize,ni,xtopt,xopt,aopt,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
          nfmf <- returnArgs[[1]]
          poped.db <- returnArgs[[2]]
          ndmf=ofv_fim(nfmf,poped.db)
          
          fprintf('BFGS minimization finished. New OFV: %g \n', ndmf)
          if((ndmf>dmf)){
            dmf=ndmf
            fmf=nfmf
            bestxt=xtopt
            bestx=xopt
            besta=aopt
          }
        }
      }
      
      if((bUseStochasticGradient)){
        tic()
        # ----------------- SG AUTO-FOCUS BEGINS HERE
        if((optxt==TRUE)){
          returnArgs <-  gradofv_xt(model_switch,axt,poped.db$design$groupsize,ni,xtopt,xopt,aopt,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
          gradxt <- returnArgs[[1]]
          poped.db <- returnArgs[[2]]
          normgxt=sign(gradxt)*(maxxt-minxt)
          xtopto=xtopt
          returnArgs <-  calc_autofocus(m,ni,dmf,xtopt,xtopto,maxxt,minxt,gradxt,normgxt,axt,model_switch,poped.db$design$groupsize,xtopt,xopt,aopt,ni,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
          axt <- returnArgs[[1]]
          poped.db <- returnArgs[[2]]
        }
        
        if((optx==TRUE)){
          # No gradient optimization for discrete values
          xopto=xopt
        }
        
        if((opta==TRUE)){
          returnArgs <- gradofv_a(model_switch,aa,poped.db$design$groupsize,ni,xtopt,xopt,aopt,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
          grada <- returnArgs[[1]]
          poped.db <- returnArgs[[2]]
          normga=grada/abs(grada)*(maxa-mina)
          aopto=aopt
          na_vector = matrix(1,1,m)*size(poped.db$design$a,2)
          returnArgs <- calc_autofocus(m,na_vector,dmf,aopt,aopto,maxa,mina,grada,normga,aa,model_switch,poped.db$design$groupsize,xtopt,xopt,aopt,ni,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
          aa <- returnArgs[[1]]
          poped.db <- returnArgs[[2]]
        }
        
        xtopt=xtopto
        xopt=xopto
        aopt=aopto
        
        if((all(axt==0) && all(aa==0))){
          diff=0
        }
        
        # ----------------- SG AUTO-FOCUS ENDS HERE
        #
        # ----------------- STOCHASTIC GRADIENT BEGINS HERE
        it=1
        dmfvector <- c()
        itvector <- c()
        while((it<=sgit) && (abs(diff)>poped.db$settings$convergence_eps)){
          if((optxt==TRUE)){
            returnArgs <- gradofv_xt(model_switch,axt,poped.db$design$groupsize,ni,xtopto,xopto,aopto,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
            gradxt <- returnArgs[[1]]
            poped.db <- returnArgs[[2]]
            returnArgs <- sg_search(gradxt,mnormxt,axt,maxxt,minxt,xtopto,lgxto,kitxt,it,m,maxni) 
            lgxto <- returnArgs[[1]]
            kitxt <- returnArgs[[2]]
            inversionxt <- returnArgs[[3]]
            xtopt <- returnArgs[[4]]
          }
          if((opta==TRUE)){
            returnArgs <- gradofv_a(model_switch,aa,poped.db$design$groupsize,ni,xtopto,xopto,aopto,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
            grada <- returnArgs[[1]]
            poped.db <- returnArgs[[2]]
            returnArgs <- sg_search(grada,mnorma,aa,maxa,mina,aopto,lgao,kita,it,m,size(poped.db$design$a,2)) 
            lgao <- returnArgs[[1]]
            kita <- returnArgs[[2]]
            inversiona <- returnArgs[[3]]
            aopt <- returnArgs[[4]]
          }
          xtopto=xtopt
          xopto=xopt
          aopto=aopt
          
          
          if((sum(sum(is.nan(xtopt)))>0 || sum(sum(is.nan(xopt)))>0 || sum(sum(is.nan(aopt)))>0)){
            break
          }
          
          
          returnArgs <-  mftot(model_switch,poped.db$design$groupsize,ni,xtopto,xopto,aopto,bpop,d,poped.db$parameters$sigma,docc_full,poped.db) 
          nfmf <- returnArgs[[1]]
          poped.db <- returnArgs[[2]]
          ndmf=ofv_fim(nfmf,poped.db)
          
          if((ndmf>dmf)){
            dmf=ndmf
            fmf=nfmf
            bestxt=xtopt
            bestx=xopt
            besta=aopt
          }
          
          if((!isempty(poped.db$settings$strIterationFileName))){
            write_iterationfile('Stochastic Gradient',it,bestxt,besta,bestx,ni,poped.db$design$groupsize,fmf,dmf,poped.db)
          }
          
          if((poped.db$settings$convergence_eps!=0 || trflag==TRUE)){
            if((it>1)){
              if((odmf!=0)){
                diff=abs((ndmf-odmf))/abs(odmf)
              }
            }
            odmf=ndmf
          }
          if((trflag==TRUE && ((it %% poped.db$settings$sgit_output)==0 || abs(diff)<poped.db$settings$convergence_eps || it==sgit))){
            itvector[ceiling(it/poped.db$settings$sgit_output)]=it
            dmfvector[ceiling(it/poped.db$settings$sgit_output)]=dmf
            ga_tmp=0
            if(opta) ga_tmp = grada/dmf
            gxt_tmp=0
            if(optxt) gxt_tmp = gradxt/dmf
            Dtrace(fn,poped.db$settings$rsit+it,ni,bestxt,bestx,besta,gxt_tmp,ga_tmp,dmf,diff,inversionxt,inversiona,itvector,dmfvector,poped.db)
            #Dtrace(fn,poped.db$settings$rsit+it,ni,bestxt,bestx,besta,normgxt,normga,dmf,diff,inversionxt,inversiona,itvector,dmfvector,poped.db)
            inversionxt=FALSE
            inversiona=FALSE
          }
          it=it+1
        }
        # ----------------- STOCHASTIC GRADIENT ENDS HERE
        sg_time=toc(echo=FALSE)
        fprintf('Stochastic gradient run time: %g seconds\n\n',sg_time)
        #if(trflag) fprintf(fn,'Elapsed time for stochastic gradient search with %g iterations: %g seconds\n',sgit,sg_time)
        
        
      }
      xtopt=bestxt
      xopt=bestx
      aopt=besta
      xt=xtopt
      x=xopt
      a=aopt
      
      #poped.db$gxt = xtopt
      #poped.db$gx = xopt
      #poped.db$design$a = aopt
      
      poped.db$design$xt <- xtopt
      poped.db$design$x <-xopt
      poped.db$design$a <-aopt
      
      if((bUseLineSearch)){
        #------------------------------- LINE SEARCH optimization START HERE
        strLineSearchFile=sprintf('%s_LS_%g%s',poped.db$settings$strOutputFileName,iter,poped.db$settings$strOutputFileExtension)
        if (!is.character(poped.db$settings$strOutputFilePath)) poped.db$settings$strOutputFilePath = '.'
        strLineSearchFile = file.path(poped.db$settings$strOutputFilePath,strLineSearchFile)
        #returnArgs <- a_line_search(strLineSearchFile,FALSE,0,fmf,dmf,poped.db) 
        returnArgs <- a_line_search(poped.db,fn,FALSE,0,fmf,dmf) 
        fmf <- returnArgs[[1]]
        dmf <- returnArgs[[2]]
        test_change <- returnArgs[[3]]
        xt <- returnArgs[[4]]
        x <- returnArgs[[5]]
        a <- returnArgs[[6]]
        poped.db <- returnArgs[[7]]
        #------------------------------- LINE SEARCH optimization ENDS HERE
      } else {
        iMaxSearchIterations = iMaxSearchIterations-1
      }
      
      
    }
    #--------- Write results
    #     if((trflag)){
    blockfinal(fn,fmf,dmf,poped.db$design$groupsize,ni,xt,x,a,model_switch,
                 bpopdescr[,2],ddescr,poped.db$parameters$docc,poped.db$parameters$sigma,poped.db,
                 opt_xt=optxt,opt_a=opta,opt_x=optx,fmf_init=fmf_init,
                 dmf_init=dmf_init,trflag=trflag,...)
    #}
    #blockfinal(fn,fmf,dmf,poped.db$design$groupsize,ni,xt,x,a,bpopdescr,ddescr,poped.db$parameters$docc,poped.db$parameters$sigma,m,poped.db)
    
    #close(fn)
    
  } else {
    fprintf('No Continuous Optimization Performed \n')
  } # matches -- if (optxt|optx|opta)
  return(list( xt= xt,x=x,a=a,fmf=fmf,dmf=dmf,poped.db =poped.db )) 
}

#' Compute the autofocus portion of the stochastic gradient routine
#' 
#' Compute the autofocus portion of the stochastic gradient routine
#'    
#' @return A list containing:
#' \item{navar}{The autofocus parameter.}
#' \item{poped.db}{PopED database.}
#' 
#' 
#' @inheritParams RS_opt
#' @inheritParams evaluate.fim
#' @inheritParams Doptim
#' @inheritParams create.poped.database
#' @param ni_var The ni_var.
#' @param varopt The varopt.
#' @param varopto The varopto.
#' @param maxvar The maxvar.
#' @param minvar The minvar.
#' @param gradvar The gradvar.
#' @param normgvar The normgvar.
#' @param avar The avar.
#' @param xtopt The optimal sampling times matrix.
#' @param xopt The optimal discrete design variables matrix.
#' @param aopt The optimal continuous design variables matrix.
#' @family Optimize
#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
#' @example tests/testthat/examples_fcn_doc/examples_calc_autofocus.R
#' @export
#' @keywords internal
calc_autofocus <- function(m,ni_var,dmf,varopt,varopto,maxvar,minvar,gradvar,normgvar,
                           avar,model_switch,groupsize,xtopt,xopt,aopt,ni,bpop,d,sigma,docc,poped.db){
  navar = avar
  for(i in 1:m){
    for(ct1 in 1:ni_var[i]){
      if((varopt[i,ct1]==maxvar[i,ct1] && gradvar[i,ct1]>0)){
        avar[i,ct1]=0
      }
      if((varopt[i,ct1]==minvar[i,ct1] && gradvar[i,ct1]<0)){
        avar[i,ct1]=0
      }
      if((avar[i,ct1]!=0)){
        varopt=varopto
        tavar=avar[i,ct1]
        varopt[i,ct1]=varopto[i,ct1]+tavar*normgvar[i,ct1]
        varopt[i,ct1]=min(maxvar[i,ct1], max(varopt[i,ct1], minvar[i,ct1]))
        returnArgs <-  mftot(model_switch,groupsize,ni,xtopt,xopt,aopt,bpop,d,sigma,docc,poped.db) 
        mf_tmp <- returnArgs[[1]]
        poped.db <- returnArgs[[2]]
        ndmf=ofv_fim(mf_tmp,poped.db)
        ct2=1
        while(ct2<=10 && ndmf<dmf){
          tavar=tavar/2
          varopt[i,ct1]=varopto[i,ct1]+tavar*normgvar[i,ct1]
          varopt[i,ct1]=min(maxvar[i,ct1], max(varopt[i,ct1], minvar[i,ct1]))
          returnArgs <- mftot(model_switch,groupsize,ni,xtopt,xopt,aopt,bpop,d,sigma,docc,poped.db) 
          mf_tmp <- returnArgs[[1]]
          poped.db <- returnArgs[[2]]
          ndmf=ofv_fim(mf_tmp,poped.db)
          ct2=ct2+1
        }
        if((ndmf<dmf)){
          navar[i,ct1]=0
        } else {
          navar[i,ct1]=tavar
        }
      }
    }
  }
  return(list( navar= navar,poped.db=poped.db)) 
}

sg_search <- function(graddetvar,mnormvar,avar,maxvar,minvar,varopto,lgvaro,oldkitvar,it,m,numvar){
  mnormvar=((mnormvar*(it-1))+(graddetvar^2))/it
  normgvar=(avar/oldkitvar)*graddetvar/sqrt(mnormvar)*(maxvar-minvar)
  varopt=varopto+normgvar
  varopt=varopt-((varopt>maxvar)*(varopt-maxvar))
  varopt=varopt+((varopt<minvar)*(minvar-varopt))
  lgvar=varopt-varopto
  dim(lgvar) = c(m*numvar,1)
  if(any(t(lgvar)%*%lgvaro<0)){
    kitvar=oldkitvar+1
    inversionvar=TRUE
  } else {
    inversionvar=FALSE
    kitvar=oldkitvar
  }
  lgvaro=lgvar
  return(list( lgvaro= lgvaro, kitvar= kitvar, inversionvar= inversionvar, varopt= varopt)) 
}

# sign <- function(x){
# # s = zeros(size(x))
# s = x/abs(x)*(x!=0)
# return( s) 
# }


#' Compute an objective function and gradient 
#' 
#' Compute an objective function and gradient with respect to the optimization parameters.
#' This function can be passed to the Broyden Fletcher Goldfarb Shanno (BFGS) 
#' method for nonlinear minimization with box constraints implemented in \code{\link{bfgsb_min}}.
#'   
#' @return A list containing:
#' \item{f}{The objective function.}
#' \item{g}{The gradient.}
#' 
#' 
#' @inheritParams RS_opt
#' @inheritParams evaluate.fim
#' @inheritParams Doptim
#' @inheritParams create.poped.database
#' @inheritParams calc_autofocus
#' @param aa The aa value
#' @param axt the axt value
#' @param xtopto the xtopto value
#' @param xopto the xopto value
#' @param optxt If sampling times are optimized
#' @param opta If continuous design variables are optimized
#' @param aopto the aopto value
#' @param only_fim Should the gradient be calculated?
#' @family Optimize
#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
#' @example tests/testthat/examples_fcn_doc/examples_calc_ofv_and_grad.R
#' @export
#' @keywords internal
calc_ofv_and_grad <- function(x,optxt,opta, model_switch,aa,axt,groupsize,ni,xtopto,
                              xopto,aopto,bpop,d,sigma,docc_full,poped.db,only_fim=FALSE){
  if(optxt){
    notfixed <- poped.db$design_space$minxt!=poped.db$design_space$maxxt
    if(poped.db$design_space$bUseGrouped_xt){
      xtopto[notfixed]=x[poped.db$design_space$G_xt[notfixed]]
      ##x[1:numel(unique(poped.db$design_space$G_xt[notfixed]))]=matrix(0,0,0)
      x=x[-c(1:numel(unique(poped.db$design_space$G_xt[notfixed])))]
    } else {
      xtopto[notfixed]=x[1:numel(xtopto[notfixed])]
      x=x[-c(1:numel(xtopto[notfixed]))]
      ##x[1:numel(xtopto[notfixed])]=matrix(0,0,0)
    }
  }
  if(opta){
    notfixed <- poped.db$design_space$mina!=poped.db$design_space$maxa
    if(poped.db$design_space$bUseGrouped_a){
      aopto[notfixed]=x[poped.db$design_space$G_a[notfixed]]
    } else {
      aopto[notfixed]=x
    }
  }
  returnArgs <-  mftot(model_switch,groupsize,ni,xtopto,xopto,aopto,bpop,d,sigma,docc_full,poped.db) 
  nfmf <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  f=-ofv_fim(nfmf,poped.db)
  gradxt=matrix(0,0,0)
  grada=matrix(0,0,0)
  if((only_fim)){
    g=0
    return(list( f= f,g=g))
  }
  if((optxt==TRUE)){
    notfixed=poped.db$design_space$minxt!=poped.db$design_space$maxxt
    returnArgs <- gradofv_xt(model_switch,axt,groupsize,ni,xtopto,xopto,aopto,bpop,d,sigma,docc_full,poped.db) 
    gradxt <- returnArgs[[1]]
    poped.db <- returnArgs[[2]]
    gradxt=gradxt[notfixed]
    if(poped.db$design_space$bUseGrouped_xt){
      returnArgs <- unique(poped.db$design_space$G_xt) 
      temp1 <- returnArgs[[1]]
      index <- returnArgs[[2]]
      temp2 <- returnArgs[[3]]
      gradxt=gradxt[index]
    }
  }
  if((opta==TRUE)){
    notfixed=poped.db$design_space$mina!=poped.db$design_space$maxa
    returnArgs <- gradofv_a(model_switch,aa,groupsize,ni,xtopto,xopto,aopto,bpop,d,sigma,docc_full,poped.db) 
    grada <- returnArgs[[1]]
    poped.db <- returnArgs[[2]]
    grada=grada[notfixed]
    if(poped.db$design_space$bUseGrouped_a){
      returnArgs <- unique(poped.db$design_space$G_a) 
      temp1 <- returnArgs[[1]]
      index <- returnArgs[[2]]
      temp2 <- returnArgs[[3]]
      grada=grada[index]
    }
  }
  g=-matrix(c(gradxt,grada),ncol=1,byrow=T)
  return(list( f= f,g=g)) 
}