Merge remote-tracking branch 'upstream/master'

Merge in more PopED upstream bug fixes/improvements
(e.g. diff equation solver defaults)

NAMESPACE

@@ -27,6 +27,7 @@ export(diag_matlab)
 export(downsizing_general_design)
 export(ed_laplace_ofv)
 export(ed_mftot)
+export(efficiency)
 export(evaluate.e.ofv.fim)
 export(evaluate.fim)
 export(feps.add)

R/blockfinal.R

@@ -126,9 +126,11 @@ blockfinal <- function(fn,fmf,dmf,groupsize,ni,xt,x,a,model_switch,bpop,d,docc,s
     #           ofv_criterion(dmf,npar,poped.db)/ofv_criterion(dmf_init,npar,poped.db))
     # }
 
-    fprintf(fn,'\nEfficiency (Final/Initial): %0.5g\n',
-            ofv_criterion(dmf,npar,poped.db)/ofv_criterion(dmf_init,npar,poped.db),both=TRUE)
-
+    eff <- efficiency(dmf_init, dmf, poped.db)
+    fprintf(fn,"\nEfficiency: \n  (%s) = %.5g\n",attr(eff,"description"),eff,both=TRUE)
+    # fprintf(fn,'\nEfficiency (Final/Initial): %0.5g\n',
+    #         ofv_criterion(dmf,npar,poped.db)/ofv_criterion(dmf_init,npar,poped.db),both=TRUE)
+    # 
     #fprintf(fn,'\nEfficiency criterion: det(FIM)^(1/npar) = %g\n',dmf^(1/length(params)))
     #fprintf(fn,'\nEfficiency (final_design/initial_design): %g\n',(dmf^(1/length(params)))/(dmf_init^(1/length(params))))
     #if(fn!="") fprintf('\nEfficiency (final_design/initial_design): %g\n',(dmf^(1/length(params)))/(dmf_init^(1/length(params))))

R/create.poped.database.R

@@ -220,7 +220,7 @@
 #' @param hle Step length of derivative of model w.r.t. sigma
 #' @param AbsTol The absolute tolerance for the diff equation solver
 #' @param RelTol The relative tolerance for the diff equation solver
-#' @param iDiffSolverMethod The diff equation solver method, 0, no other option 
+#' @param iDiffSolverMethod The diff equation solver method, NULL as default.
 #' @param bUseMemorySolver If the differential equation results should be stored in memory (1) or not (0)
 #' @param rsit Number of Random search iterations 
 #' @param sgit Number of stochastic gradient iterations
@@ -576,11 +576,11 @@ create.poped.database <-
            ## -- Step length of derivative of model w.r.t. sigma --
            hle=poped.choose(popedInput$hle,0.00001),
            ## -- The absolute tolerance for the diff equation solver --
-           AbsTol=poped.choose(popedInput$AbsTol,0.00001),
+           AbsTol=poped.choose(popedInput$AbsTol,0.000001),
            ## -- The relative tolerance for the diff equation solver --
-           RelTol=poped.choose(popedInput$RelTol,0.00001),
+           RelTol=poped.choose(popedInput$RelTol,0.000001),
            ## -- The diff equation solver method, 0, no other option --
-           iDiffSolverMethod=poped.choose(popedInput$iDiffSolverMethod,0),
+           iDiffSolverMethod=poped.choose(popedInput$iDiffSolverMethod,NULL),
            ## -- If the differential equation results should be stored in memory (1) or not (0) --
            bUseMemorySolver=poped.choose(popedInput$bUseMemorySolver,FALSE),
            ## -- Number of Random search iterations --

R/create_design.R

@@ -79,7 +79,11 @@ create_design <- function(
       #all_cov_names <- unique(unlist(sapply(a,names)))
 
       #a <- as.matrix(plyr::rbind.fill(lapply(a,function(x){data.frame(rbind(unlist(x)))})))
-      a <- as.matrix(dplyr::rbind_all(lapply(a,function(x){data.frame(rbind(unlist(x)))})))
+      if(packageVersion("dplyr") >= "0.5.0"){
+        a <- as.matrix(dplyr::bind_rows(lapply(a,function(x){data.frame(rbind(unlist(x)))})))
+      } else {
+        a <- as.matrix(dplyr::rbind_all(lapply(a,function(x){data.frame(rbind(unlist(x)))})))
+      }
     }
     colnam <- names(a)
     if(is.null(colnam)) colnam <- colnames(a)
@@ -93,7 +97,13 @@ create_design <- function(
 
   ## for x ----------
   if(!is.null(x)){
-    if(is.list(x)) x <- as.matrix(dplyr::rbind_all(lapply(x,function(x){data.frame(rbind(unlist(x)))})))
+    if(is.list(x)){
+      if(packageVersion("dplyr") >= "0.5.0"){
+        x <- as.matrix(dplyr::bind_rows(lapply(x,function(x){data.frame(rbind(unlist(x)))})))
+      } else {
+        x <- as.matrix(dplyr::rbind_all(lapply(x,function(x){data.frame(rbind(unlist(x)))})))
+      }
+    } 
     colnam <- names(x)
     if(is.null(colnam)) colnam <- colnames(x)
     if(size(x,1)==1 && m!=1) x <- matrix(rep(x,m),ncol=length(x),nrow=m,byrow=T,dimnames=list(paste("grp_",1:m,sep=""),colnam))

R/create_design_space.R

@@ -348,7 +348,11 @@ create_design_space <- function(
     ## for a ---------
     if(!is.null(maxa)){
       if(is.list(maxa)){
-        maxa <- as.matrix(dplyr::rbind_all(lapply(maxa,function(x){data.frame(rbind(unlist(x)))})))
+        if(packageVersion("dplyr") >= "0.5.0"){
+          maxa <- as.matrix(dplyr::bind_rows(lapply(maxa,function(x){data.frame(rbind(unlist(x)))})))
+        } else {
+          maxa <- as.matrix(dplyr::rbind_all(lapply(maxa,function(x){data.frame(rbind(unlist(x)))})))
+        }
       }
       if(size(maxa,1)==1 && m!=1) maxa <- matrix(rep(maxa,m),ncol=length(maxa),nrow=m,byrow=T)
       if(!is.matrix(maxa)) maxa  <- rbind(maxa)
@@ -360,7 +364,11 @@ create_design_space <- function(
 
     if(!is.null(mina)){
       if(is.list(mina)){
-        mina <- as.matrix(dplyr::rbind_all(lapply(mina,function(x){data.frame(rbind(unlist(x)))})))
+        if(packageVersion("dplyr") >= "0.5.0"){
+          mina <- as.matrix(dplyr::bind_rows(lapply(mina,function(x){data.frame(rbind(unlist(x)))})))
+        } else {
+          mina <- as.matrix(dplyr::rbind_all(lapply(mina,function(x){data.frame(rbind(unlist(x)))})))
+        }
       }
       if(size(mina,1)==1 && m!=1) mina <- matrix(rep(mina,m),ncol=length(mina),nrow=m,byrow=T)
       if(!is.matrix(mina)) mina  <- rbind(mina)

R/model_prediction.R

@@ -292,7 +292,11 @@ model_prediction <- function(poped.db=NULL,
               if(length(unique(tmp.df[nam]))==1) dose.df[nam] <- tmp.df[1,nam]  
             }
             dose.df$dose_record_tmp <- 1
-            tmp.df <- dplyr::rbind_list(dose.df,tmp.df)
+            if(packageVersion("dplyr") >= "0.5.0"){
+              tmp.df <- dplyr::bind_rows(dose.df,tmp.df)
+            } else {
+              tmp.df <- dplyr::rbind_list(dose.df,tmp.df)
+            }
             tmp.df <- tmp.df[order(tmp.df$Time,tmp.df$dose_record_tmp),]
             tmp.df$dose_record_tmp <- NULL
           }

R/plot_efficiency_of_windows.R

@@ -79,9 +79,9 @@ plot_efficiency_of_windows <- function(poped.db,
 
   if(y_eff){
     eff = zeros(1,iNumSimulations)
-    if(ofv_calc_type==4) {
-      d_eff = zeros(1,iNumSimulations)
-    }
+    # if(ofv_calc_type==4) {
+    #   d_eff = zeros(1,iNumSimulations)
+    # }
   }
   if(y_rse) rse <- zeros(iNumSimulations,p)
 
@@ -206,13 +206,13 @@ plot_efficiency_of_windows <- function(poped.db,
       tmp2 <- ofv_criterion(ofv_fim(ref_fmf,poped.db,ofv_calc_type=ofv_calc_type),
                             p,poped.db,ofv_calc_type=ofv_calc_type) 
       eff[1,i] = tmp1/tmp2  
-      if(ofv_calc_type==4) {
-        tmp1 <- ofv_criterion(ofv_fim(fmf,poped.db,ofv_calc_type=1),
-                              p,poped.db,ofv_calc_type=1)
-        tmp2 <- ofv_criterion(ofv_fim(ref_fmf,poped.db,ofv_calc_type=1),
-                              p,poped.db,ofv_calc_type=1) 
-        d_eff[1,i] = tmp1/tmp2  
-      }
+      # if(ofv_calc_type==4) {
+      #   tmp1 <- ofv_criterion(ofv_fim(fmf,poped.db,ofv_calc_type=1),
+      #                         p,poped.db,ofv_calc_type=1)
+      #   tmp2 <- ofv_criterion(ofv_fim(ref_fmf,poped.db,ofv_calc_type=1),
+      #                         p,poped.db,ofv_calc_type=1) 
+      #   d_eff[1,i] = tmp1/tmp2  
+      # }
     }
     if(y_rse){
       rse_tmp <- get_rse(fmf,poped.db)
@@ -324,16 +324,16 @@ plot_efficiency_of_windows <- function(poped.db,
   ind <- NULL
   if(y_eff){
     efficiency <- eff[1,]*100
-    if(ofv_calc_type==4) {
-      d_efficiency <-  d_eff[1,]*100
-    }
+    # if(ofv_calc_type==4) {
+    #   d_efficiency <-  d_eff[1,]*100
+    # }
   }
   df <- data.frame(sample=c(1:iNumSimulations))
   if(y_eff){
     df$Efficiency <- efficiency
-    if(ofv_calc_type==4) {
-      df["D-Efficiency"] <- d_efficiency
-    }
+    # if(ofv_calc_type==4) {
+    #   df["D-Efficiency"] <- d_efficiency
+    # }
   }
   if(y_rse){
     rse_df <- data.frame(rse)
@@ -345,11 +345,11 @@ plot_efficiency_of_windows <- function(poped.db,
   names(df_stack) <- c("sample","values","ind")
   if(y_eff){
     levs <- levels(df_stack$ind)
-    if(ofv_calc_type==4) {
-      df_stack$ind <- factor(df_stack$ind,levels=c("D-Efficiency","Efficiency",levs[-c(grep("Efficiency",levs))]))        
-    } else {
-      df_stack$ind <- factor(df_stack$ind,levels=c("Efficiency",levs[-c(grep("Efficiency",levs))]))      
-    }
+    # if(ofv_calc_type==4) {
+    #   df_stack$ind <- factor(df_stack$ind,levels=c("D-Efficiency","Efficiency",levs[-c(grep("Efficiency",levs))]))        
+    # } else {
+    df_stack$ind <- factor(df_stack$ind,levels=c("Efficiency",levs[-c(grep("Efficiency",levs))]))      
+    # }
     #levels(df_stack$ind) <- c("Efficiency",levs[-c(grep("Efficiency",levs))])
   }
 

R/poped_optim.R

@@ -548,9 +548,11 @@ poped_optim <- function(poped.db,
       fprintf("Relative difference in OFV:  %.3g%%\n",rel_diff*100)
 
       # efficiency
-      p = get_fim_size(poped.db)
-      eff = ofv_criterion(output$ofv,p,poped.db)/ofv_criterion(ofv_init,p,poped.db)
-      cat("Efficiency: ",sprintf("%.5g",eff), "\n")
+
+
+      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")
@@ -708,3 +710,63 @@ poped_optim <- function(poped.db,
   return(invisible(list( ofv= output$ofv, FIM=FIM, poped.db = poped.db ))) 
 }
 
+#' Compute efficiency.
+#' 
+#' Efficiency calculation between two designs.
+#' 
+#' 
+#' @param ofv_init An initial objective function
+#' @param ofv_final A final objective function.
+#' @param npar The number of parameters to use for normalization.
+#' @param poped_db a poped database
+#' @inheritParams ofv_fim
+#' @inheritParams poped_optim
+#' @inheritParams create.poped.database
+#' 
+#' @return The specified efficiency value depending on the ofv_calc_type.  
+#' The attribute "description" tells you how the calculation was made 
+#' \code{attr(return_vale,"description")}
+#' 
+#' @family FIM
+#' 
+#' 
+## @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
+## @example tests/testthat/examples_fcn_doc/examples_ofv_criterion.R
+#' 
+#' @export
+
+efficiency <- function(ofv_init, ofv_final, poped_db,
+                       npar = get_fim_size(poped_db),
+                       ofv_calc_type=poped_db$settings$ofv_calc_type,
+                       ds_index = poped_db$parameters$ds_index) {
+
+
+  eff = ofv_final/ofv_init
+  attr(eff,"description") <- "ofv_final / ofv_init"
+
+  if((ofv_calc_type==1) ){#D-Optimal Design
+    eff = eff^(1/npar)
+    attr(eff,"description") <- "(ofv_final / ofv_init)^(1/n_parameters)"
+  }
+  if((ofv_calc_type==4) ){#lnD-Optimal Design
+    eff = (exp(ofv_final)/exp(ofv_init))^(1/npar)
+    attr(eff,"description") <- "(exp(ofv_final) / exp(ofv_init))^(1/n_parameters)"
+
+  }
+  # if((ofv_calc_type==2) ){#A-Optimal Design
+  #   eff=ofv_f/npar
+  # }
+  # 
+  # if((ofv_calc_type==3) ){#S-Optimal Design
+  #   stop(sprintf('Criterion for S-optimal design not implemented yet'))
+  # }
+  # 
+  if((ofv_calc_type==6) ){#Ds-Optimal design
+    eff = eff^(1/sum(ds_index))
+    attr(eff,"description") <- "(ofv_final / ofv_init)^(1/sum(interesting_parameters))"
+
+  }   
+
+  return( eff ) 
+}
+