Merge pull request #28 from bmreilly/main

added example for debugging NMupdateSizes in devel/NMsim_NWPRI_sizes_testing.R

R/NMsim_NWPRI.R

@@ -101,7 +101,7 @@ NMsim_NWPRI <- function(file.sim,file.mod,data.sim,PLEV=0.999){
     cov.l <- addParType(cov.l,suffix="i")
     cov.l <- addParType(cov.l,suffix="j")
 
-    # Identify fixed thetas and set their diagonal in the cov matrix to be 1.0 for compatibility with nm7.60
+    # Identify fixed thetas and set their diagonal in the $THETAPV cov matrix to be 1.0 for compatibility with nm7.60
     cov.l = merge.data.table(x = cov.l, y = pars[par.type=="THETA",.(i, parameter.i=parameter,FIX)], by = c("i", "parameter.i"))
     cov.l[i==j&FIX==1, value := 1.0]    
 
@@ -140,6 +140,15 @@ NMsim_NWPRI <- function(file.sim,file.mod,data.sim,PLEV=0.999){
 ### add TRUE=PRIOR to $SIMULATION
     lines.sim <- NMdata:::NMwriteSectionOne(lines=lines.sim, section="SIMULATION", location="after", newlines="TRUE=PRIOR", backup=FALSE, quiet=TRUE)
 
+### update $SIZES LTH and LVR to reflect the parameters in NWPRI (not resized automatically like other subroutines)
+    # add 10 to both numbers per Bob Bauer (doesn't hurt to have slightly more memory/size)
+    # LTH = number of thetas = $THETA + $THETAP + $OMEGAPD + $SIGMAPD
+    # LVR = number of diagonal omegas + sigmas = $OMEGA + $OMEGAP + $SIGMA + $SIGMAP + $THETAPV
+    lth = 2*nrow(pars[par.type=="THETA"]) + nrow(nwpri_df) + 10 
+    lvr = 2*nrow(pars[(par.type=="OMEGA"|par.type=="SIGMA")&i==j]) + nrow(pars[par.type=="THETA"]) + 10
+
+    lines.sim = NMsim::NMupdateSizes(file.mod=NULL, newfile=NULL,lines = lines.sim, LTH = lth, LVR = lvr)
+
 ### update the simulation control stream
     ## if(return.text){
     ##     return(lines.sim)

R/NMupdateSizes.R

@@ -73,12 +73,30 @@ NMupdateSizes <- function(file.mod=NULL,newfile=file.mod,lines=NULL,wipe=FALSE,w
     ),collapse=" "
     )
 
-    if(!is.null(sizes.old)){
-        lines <- NMwriteSection(files=file.mod,newfile=newfile,section="SIZES",newlines="",location="replace",write=FALSE)
-    } else if(!is.null(file.mod)) {
+    # create lines from file.mod if not passed in as 'lines'
+    if(!is.null(file.mod) & is.null(lines)) {
         lines <- readLines(file.mod,warn=FALSE)
+    } 
+
+    # if model already has $SIZES and we want to overwrite or append to it.
+    if(!is.null(sizes.old)) {
+        # if we want to overwrite it or append to it, sizes.new and lines.new will already be modified for that, so replace
+        textlines = NMdata:::NMwriteSectionOne(lines = lines,section="SIZES",newlines=lines.new,location="replace")
+
+    # else if model does not have $SIZES, create it
+    } else if (is.null(NMreadSizes(lines=lines))) {
+        textlines = NMdata:::NMwriteSectionOne(lines=lines,section="SIZES",newlines=lines.new,location="first")
+
+    # else if it had $SIZES but we wiped it, and we want to replace/append it:
+    } else if (!is.null(NMreadSizes(lines=lines)) & wipe) {
+        textlines = NMdata:::NMwriteSectionOne(lines=lines,section="SIZES",newlines=lines.new,location="replace")
+    }
+
+    if(write & !is.null(newfile)) {
+        writeTextFile(textlines, newfile)
+    } else {
+        return(textlines)
     }
-    NMdata:::NMwriteSectionOne(lines=lines,newfile=newfile,section="SIZES",newlines=lines.new,location="first",write=write)
 
 }
 

devel/NMsim_NWPRI_sizes_testing.R

@@ -0,0 +1,141 @@
+#### testing NMsim::NMsim_NWPRI after accounting for $SIZES LVR and LTH
+
+unloadNamespace("NMsim")
+unloadNamespace("NMdata")
+library(devtools)
+library(here)
+library(scales)
+library(tidyverse)
+library(stringr)
+
+
+### lets try to make the directories in the script only depend on wdirs. And maybe one more for setwd if needed
+wdirs <- "/data/sandbox/trunk/analysis/NMsim/wdirs"
+setwd("/data/sandbox/trunk/analysis/NMsim/wdirs")
+# wdirs <- "~/wdirs"
+load_all(file.path(wdirs,"NMdata"))
+load_all(file.path(wdirs,"NMsim"))
+
+
+NMdataConf(path.nonmem="/data/sandbox/trunk/analysis/NMsim/nonmem_pre_releases/nm760/run/nmfe76",
+           as.fun="data.table")
+# NMdataConf(path.nonmem="/opt/NONMEM/nm75/run/nmfe75",
+#            as.fun="data.table")
+
+
+printSection <- function(section) {
+  if(!is.list(section)) section <- list(section)
+  res <- lapply(section,function(x){
+    x <- x[!grepl("^ *$",section)]
+    cat(paste(paste(x,collapse="\n"),"\n\n"))
+
+  })
+}
+
+##### Choose a model
+
+### The model shown in the NMsim-ParamUncertain vignette
+file.project <- function(...)file.path(system.file("examples",package="NMsim"),...)
+file.mod <- file.project("nonmem/xgxr032.mod")
+## dat.sim <- read_fst(here::here("wdirs/NMsim/vignettes/simulate-results/dat_sim.fst"))
+data.sim <- read_fst(file.path(wdirs,"NMsim/vignettes/simulate-results/dat_sim.fst"))
+
+
+pars <- NMreadExt(file.mod,return="pars",as.fun="data.table")
+pars[value!=0]
+pars[, parlab := stringr::str_remove_all(string = parameter, pattern = "\\(|\\)") %>% stringr::str_replace(",", "\\_")]
+
+NMsim(
+  file.mod=file.mod,              ## Path to estimation input control stream
+  data=data.sim                    ## simulation input data
+  ,dir.sims="NMsim/devel/simtmp_NWPRI_nm76"
+  # ,modify.model = list(
+  #   ## name the THETAS/OMEGAS/SIGMAS in $ERROR so we can compare
+  #   ## the distributions to other methods.
+  #   ERROR = add( paste0(pars$parlab, " = ", pars$par.name) )
+  # )  
+  # ,table.vars = paste0( "PRED IPRED Y ", paste0(pars$parlab, collapse = " ") )   
+  ,method.sim=NMsim_default
+  , execute = FALSE
+  ## Var-Cov parameter sampling
+  ,name.sim="NWPRI_sizes"               ## a recognizable directory name
+  ,subproblems=250                 ## sampling multiple models
+  ,sge=FALSE                      ## run simulations in parallel please
+  ,nmquiet=T
+  ,reuse.results=F
+)
+
+file.sim = "/data/sandbox/trunk/analysis/NMsim/wdirs/NMsim/devel/simtmp_NWPRI_nm76/xgxr032_NWPRI_sizes/xgxr032_NWPRI_sizes.mod"
+
+# a file with no SIZES line
+lines = readLines(file.sim)
+NMsim:::NMreadSizes(file.mod = file.sim) # works
+NMsim::NMupdateSizes(lines =  lines, LTV = 50) # works
+NMsim::NMupdateSizes(file.mod=file.sim, LTV = 50, write=F) # works
+
+NMsim::NMupdateSizes(file.mod = file.sim, LTV=50) #works
+NMsim::NMupdateSizes(file.mod = file.sim, newfile = fnAppend(file.sim, "updateSize"), LTV=5000) # works
+
+# a file with SIZES
+lines = c("$SIZES LIM1=10000 LTH=40 LVR=20", lines)
+NMsim::NMupdateSizes(lines = lines, LTH=50) # works
+
+
+# 2822
+file.mod = "/data/prod_vx548_phase3_analysis/trunk/analysis/NDA/models/PK/2822/2822.mod"
+data.sim = NMdata::NMscanInput(file.mod,translate = T,recover.cols = FALSE, apply.filters = TRUE) %>% 
+  # dplyr::filter(USUBJIDN==10001) %>% 
+  # dplyr::select(USUBJIDN:WTBLI) %>% 
+  rename(ID=USUBJIDN, TIME=AFRLT) %>% 
+  filter(ID %in% sample(ID, size = 10))
+
+pars <- NMreadExt(file.mod,return="pars",as.fun="data.table")
+pars[, parlab := stringr::str_remove_all(string = parameter, pattern = "\\(|\\)") %>% stringr::str_replace(",", "\\_")]
+
+# # first generate the empty sim control stream:
+NMsim::NMsim(
+    file.mod=file.mod,              ## Path to estimation input control stream
+    data=data.sim                    ## simulation input data
+    ,dir.sims="NMsim/devel/simtmp_NWPRI_nm76" ## where to store temporary simulation files
+    ,modify.model = list(
+        ## name the THETAS/OMEGAS/SIGMAS in $ERROR so we can compare
+        ## the distributions to other methods.
+        ERROR = add( paste0(pars$parlab, " = ", pars$par.name) )
+    )
+    ,table.vars = paste0( "PRED IPRED Y ", paste0(pars$parlab, collapse = " ") )
+    ,method.sim=NMsim_NWPRI         ## Var-Cov parameter sampling
+    , execute = FALSE
+    ,name.sim="NWPRI_sizes"               ## a recognizable directory name
+    ,subproblems=250                 ## sampling multiple models
+    ,sge=FALSE                      ## run simulations in parallel please
+    ,nmquiet=F
+    ,reuse.results=F
+)
+file.sim="NMsim/devel/simtmp_NWPRI_nm76/2822_NWPRI_sizes/2822_NWPRI_sizes.mod"
+
+# check if model looks right.
+readLines(file.sim)
+
+
+NMsim::NMsim(
+  file.mod=file.mod,              ## Path to estimation input control stream
+  data=data.sim                    ## simulation input data
+  ,dir.sims="NMsim/devel/simtmp_NWPRI_nm76" ## where to store temporary simulation files
+  ##   ,dir.res="simulate-results"      ## where to store simulation results files
+  ,modify.model = list(
+    ## name the THETAS/OMEGAS/SIGMAS in $ERROR so we can compare
+    ## the distributions to other methods.
+    ERROR = add( paste0(pars$parlab, " = ", pars$par.name) )
+  )
+  ,table.vars = paste0( "PRED IPRED Y ", paste0(pars$parlab, collapse = " ") )
+  ,method.sim=NMsim_NWPRI         ## Var-Cov parameter sampling
+  ,name.sim="NWPRI_sizes"               ## a recognizable directory name
+  ,subproblems=50                 ## sampling multiple models
+  ,sge=FALSE                      ## run simulations in parallel please
+  ,nmquiet=F
+  ,reuse.results=F
+  ,wait = T 
+)
+
+
+

tests/testthat/test_NMsim_NWPRI.R

@@ -2,7 +2,7 @@
 # library(testthat)
 # withr::with_libpaths(new = file.path("../../../../renv/local/"), devtools::install_github("philipdelff/NMdata@v0.1.8" ))
 # library("NMdata",lib.loc = "../../../../renv/local")
-
+# setwd("NMsim/tests/testthat")
 # devtools::load_all("../../../NMdata")
 # devtools::load_all("../../../NMsim")
 context("NMsim_NWPRI.R")
@@ -68,6 +68,8 @@ test_that("NMsim_NWPRI",{
         mod$SIMULATION
 
 compareCols(ref,mod)
+        ref$SIZES
+        mod$SIZES
 
         ref$OMEGAP
         mod$OMEGAP