Merge branch 'master' of https://github.com/andrewhooker/PopED

rikardn · Jun 19, 2016 · 5a46b87 · 5a46b87
2 parents 14c45dd + 2ef2fdc
commit 5a46b87
Show file tree

Hide file tree

Showing 16 changed files with 387 additions and 56 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: PopED
 Type: Package
 Title: Population (and Individual) Optimal Experimental Design
-Version: 0.3.0.9005
+Version: 0.3.0.9006
 Depends:
     R (>= 2.14)
 Imports:

diff --git a/NAMESPACE b/NAMESPACE
@@ -50,6 +50,7 @@ export(inv)
 export(isempty)
 export(isfield)
 export(log_prior_pdf)
+export(mc_mean)
 export(median_hilow_poped)
 export(mf)
 export(mf3)

diff --git a/NEWS b/NEWS
@@ -9,6 +9,10 @@ PopED 0.3.0.9000
 
 * Update to error handling for Laplace approximation to ED objective function. 
 
+* Update to more easily allow discrete optimization of xt and a variables.
+
+* Added more intuitive cost function input.  See examples in `?poped_optim`
+
 
 PopED 0.3.0
 ===========

diff --git a/R/blockfinal.R b/R/blockfinal.R
@@ -140,10 +140,14 @@ blockfinal <- function(fn,fmf,dmf,groupsize,ni,xt,x,a,model_switch,bpop,d,docc,s
 
 
     if(is.null(param_cvs_init) && !is.null(fmf_init) && is.matrix(fmf_init) && compute_inv){
+      if(is.finite(dmf_init)) {
+
+
       param_vars_init=diag_matlab(inv(fmf_init))
       returnArgs <-  get_cv(param_vars_init,bpop,d,docc,sigma,poped.db) 
       params_init <- returnArgs[[1]]
       param_cvs_init <- returnArgs[[2]]
+      }
     }
 
     if(compute_inv){

diff --git a/R/blockheader.R b/R/blockheader.R
@@ -124,7 +124,7 @@ blockheader <- function(poped.db,name="Default",iter=NULL,
     }
   }
 
-  if(is.matrix(fmf) && compute_inv){
+  if(is.matrix(fmf) && compute_inv && is.finite(dmf)){
     param_vars=diag_matlab(inv(fmf))
     returnArgs <-  get_cv(param_vars,bpop,d,docc,sigma,poped.db) 
     params <- returnArgs[[1]]

diff --git a/R/calc_ofv_and_fim.R b/R/calc_ofv_and_fim.R
@@ -14,6 +14,7 @@
 #' @param laplace.fim Should an E(FIM) be calculated when computing the Laplace approximated E(OFV).  Typically
 #' the FIM does not need to be computed and, if desired,  this calculation
 #' is done usng the standard MC integration technique, so can be slow. 
+#' @param evaluate_fim Should the FIM be calculated?
 #' 
 #' @return A list containing the FIM and OFV(FIM) or the E(FIM) and E(OFV(FIM)) according  to the function arguments.
 #' 
@@ -43,32 +44,49 @@ calc_ofv_and_fim <- function (poped.db,
                               fim.calc.type=poped.db$settings$iFIMCalculationType, 
                               use_laplace=poped.db$settings$iEDCalculationType, 
                               laplace.fim=FALSE, 
+                              ofv_fun = poped.db$settings$ofv_fun,
+                              evaluate_fim = TRUE,
                               ...) {
 
   ## compute the OFV
   if((ofv==0)){
     if(d_switch){ 
-      if(!is.matrix(fim)){ 
-        fmf <- evaluate.fim(poped.db,
-                            bpop.val=bpop,
-                            d_full=d,
-                            docc_full=docc_full,
-                            sigma_full=poped.db$parameters$sigma,
-                            model_switch=model_switch,
-                            ni=ni,
-                            xt=xt,
-                            x=x,
-                            a=a,
-                            groupsize=poped.db$design$groupsize,
-                            fim.calc.type=fim.calc.type,
-                            ...)
-
-        #     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]]
+      if(is.null(ofv_fun)){
+        if(!is.matrix(fim)){ 
+          fmf <- evaluate.fim(poped.db,
+                              bpop.val=bpop,
+                              d_full=d,
+                              docc_full=docc_full,
+                              sigma_full=poped.db$parameters$sigma,
+                              model_switch=model_switch,
+                              ni=ni,
+                              xt=xt,
+                              x=x,
+                              a=a,
+                              groupsize=poped.db$design$groupsize,
+                              fim.calc.type=fim.calc.type,
+                              ...)
+
+          #     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,...)
+      } else {
+        ## 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]])))
+            if(!is.null(eval(parse(text=paste(i))))) eval(parse(text=paste(capture.output(default_args[[i]]),"<-",i)))
+          }
+        }
+        dmf <- do.call(ofv_fun,list(poped.db,...))
+        fmf <- NULL
       }
-      dmf=ofv_fim(fmf,poped.db,...)
-    } else {
+    } else { # e-family
+      if(is.null(ofv_fun)){
       output <- evaluate.e.ofv.fim(poped.db,
                                    fim.calc.type=fim.calc.type,
                                    bpop=bpopdescr,
@@ -87,7 +105,22 @@ calc_ofv_and_fim <- function (poped.db,
                                    ...)
       dmf <- output$E_ofv
       fmf <- output$E_fim 
+      } else {
+        ## 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]])))
+            if(!is.null(eval(parse(text=paste(i))))) eval(parse(text=paste(capture.output(default_args[[i]]),"<-",i)))
+          }
+        }
+
+        dmf <- mc_mean(ofv_fun,poped.db,...)
+        fmf <- NULL
+      }
     }
+
     ofv <- dmf
     if(!is.matrix(fim)) fim <- fmf
   }
@@ -100,7 +133,7 @@ calc_ofv_and_fim <- function (poped.db,
     calc_fim <- TRUE
   }    
   if(use_laplace && !laplace.fim) calc_fim <- FALSE
-
+  if(!evaluate_fim) calc_fim <- FALSE
 
   if(calc_fim){
     if(d_switch){ 

diff --git a/R/create.poped.database.R b/R/create.poped.database.R
@@ -126,6 +126,9 @@
 #' \code{ofv_calc_type=6}.
 #' @param strEDPenaltyFile Penalty function name or path and filename, empty string means no penalty.  
 #' User defined criterion can be defined this way.
+#' @param ofv_fun User defined function used to compute the objective function. The function must have a poped database object as its first
+#' argument and have "..." in its argument list.  Can be referenced as a function or as a file name where the funciton defined in the file has the same name as the file.
+#' e.g. "cost.txt" has a function named "cost" in it.
 
 #' @param iEDCalculationType  \itemize{
 #' \item \bold{******START OF E-FAMILY CRITERION SPECIFICATION OPTIONS**********}}
@@ -415,6 +418,7 @@ create.poped.database <-
            ds_index=popedInput$CriterionOptions$ds_index,   
            ## -- Penalty function, empty string means no penalty.  User defined criterion --
            strEDPenaltyFile=poped.choose(popedInput$strEDPenaltyFile,''),
+           ofv_fun = NULL,
 
 
 
@@ -949,6 +953,42 @@ create.poped.database <-
       }
     } 
 
+    # if(is.null(ofv_fun) || is.function(ofv_fun)){
+    #   poped.db$settings$ofv_fun = ofv_fun
+    # } else {
+    #   stop("ofv_fun must be a function or NULL")
+    # }
+
+    if(is.null(ofv_fun) || is.function(ofv_fun)){
+      poped.db$settings$ofv_fun <- 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))
+        poped.db$settings$ofv_fun <- 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.function(ofv_fun)){
+    #   poped.db$settings$ofv_fun = ofv_fun 
+    # } else if(exists(ofv_fun)){
+    #   poped.db$settings$ofv_fun = ofv_fun 
+    # } else {
+    #   source(ofv_fun)
+    #   returnArgs <-  fileparts(ofv_fun) 
+    #   strffModelFilePath <- returnArgs[[1]]
+    #   strffModelFilename  <- returnArgs[[2]]
+    #   ## if (~strcmp(strffModelFilePath,''))
+    #   ##    cd(strffModelFilePath);
+    #   ## end
+    #   poped.db$settings$ofv_fun = strffModelFilename
+    # }
+
+
     poped.db$model$auto_pointer=zeros(1,0)
     if((!as.character(strAutoCorrelationFile)=='')){
       if(exists(strAutoCorrelationFile)){

diff --git a/R/mc_mean.R b/R/mc_mean.R
@@ -0,0 +1,70 @@
+#' Compute the monte-carlo mean of a function
+#'
+#'Function computes the monte-carlo mean of a function by varying the parameter inputs to the function
+#'
+#' @param ofv_fcn A function with poped.db as the first input
+#' @inheritParams evaluate.fim
+#' @inheritParams Doptim
+#' @inheritParams create.poped.database
+#' @param doccdescr Matrix defining the IOV.
+#' 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 of the variance.
+#' \item column 3 defines the variance of the distribution (or length of uniform distribution).
+#' }
+#' @param user_distribution_pointer Function name for user defined distributions for E-family designs 
+#' @return The mean of the function evaluated at different parameter values.
+#' @export
+#'
+# @examples
+
+mc_mean <- function(
+  ofv_fcn,
+  poped.db,
+  bpopdescr=poped.db$parameters$bpop, 
+  ddescr=poped.db$parameters$d,
+  doccdescr=poped.db$parameters$d,
+  user_distribution_pointer=poped.db$model$user_distribution_pointer,
+  ED_samp_size=poped.db$settings$ED_samp_size,
+  bLHS=poped.db$settings$bLHS,
+  ...) 
+{
+  ofv_sum=0
+
+
+  bpop_gen  <-  pargen(bpopdescr,
+                       user_distribution_pointer,
+                       ED_samp_size,
+                       bLHS,
+                       NULL,#zeros(1,0),
+                       poped.db)
+
+  d_gen <- pargen(ddescr,
+                  user_distribution_pointer,
+                  ED_samp_size,
+                  bLHS,
+                  NULL,
+                  poped.db)
+
+  docc_gen <- pargen(doccdescr,
+                     user_distribution_pointer,
+                     ED_samp_size,
+                     bLHS,
+                     NULL,
+                     poped.db)
+
+  poped.db_tmp <- poped.db
+  for(ct in 1:ED_samp_size){
+    poped.db_tmp$parameters$bpop[,2] <- bpop_gen[ct,]
+    poped.db_tmp$parameters$d[,2] <- d_gen[ct,]
+    poped.db_tmp$parameters$docc[,2] <- docc_gen[ct,]
+
+    dmf_tmp <- do.call(ofv_fcn,list(poped.db_tmp,...))
+
+    ofv_sum=ofv_sum+dmf_tmp
+  }
+  ofv_mean=ofv_sum/poped.db$settings$ED_samp_size
+  return(ofv_mean)
+}
diff --git a/R/ofv_criterion.R b/R/ofv_criterion.R
@@ -27,6 +27,10 @@ ofv_criterion <- function(ofv_f,
 
   criterion_value = 0
 
+  if(ofv_calc_type==0){
+    criterion_value=ofv_f
+  }
+
   if((ofv_calc_type==1 || ofv_calc_type==4) ){#D-Optimal Design
     criterion_value = ofv_f^(1/num_parameters)
   }