metrumresearchgroup / bbr

R interface for model and project management
https://metrumresearchgroup.github.io/bbr/
Other
22 stars 2 forks source link

Function to run NMTRAN on a model object #650

Open barrettk opened 5 months ago

barrettk commented 5 months ago

This idea has been brought up numerous times, but the core idea is to provide a method for validating a bbr model object via NMTRAN before submitting the model. There are other applications for this as well, such as validating a recently tweaked model, or just ensuring NONMEM syntax is correct in general.

Limitations

Below is a prototype scaffold of how we can determine the NMTRAN executable path:

Prototype functions ```r locate_nmtran <- function(.mod, .config_path = NULL, nmtran_exe = NULL){ if(is.null(nmtran_exe)){ model_dir <- get_model_working_directory(.mod) config_path <- .config_path %||% file.path(model_dir, "bbi.yaml") if(!file_exists(config_path)){ stop(paste("No bbi configuration was found in the execution directory.", "Please run `bbi_init()` with the appropriate directory to continue.")) } if (!is.null(.config_path)) { config_path <- normalizePath(.config_path) } bbi_config <- yaml::read_yaml(config_path) nm_config <- bbi_config$nonmem # look for default nonmem installation default_nm <- purrr::keep(nm_config, function(nm_ver){ !is.null(nm_ver$default) }) # Set nonmem path if(length(default_nm) > 0){ default_nm <- default_nm[[1]] }else{ # If no default, use the last one (likely higher version) default_nm <- nm_config[[length(nm_config)]] } # Set NMTRAN executable path nm_path <- default_nm$home nmtran_exe <- file.path(nm_path, "tr", "NMTRAN.exe") } if(!file_exists(nmtran_exe)){ stop(glue("Could not find an NMTRAN executable at `{nmtran_exe}`")) } return(nmtran_exe) } run_nmtran <- function(.mod, .config_path = NULL, nmtran_exe = NULL){ nmtran_exe <- locate_nmtran(.mod, .config_path, nmtran_exe) } ```
barrettk commented 5 months ago

FYI I finalized a working prototype. Can open up a scratch branch/draft PR, though I imagine there might be some discussion around A) supporting windows, and B) whether or not we need to copy the model to a temporary location like I have. It's worth noting I took quite a bit from NMProject (here & here). The main differences are A) how we determine the NMTRAN executable path (using bbi.yaml) and B) using bbr and nmrec for getting model paths and modifying/fetching data paths.

run_nmtran & helper functions ```r run_nmtran <- function( .mod, .config_path = NULL, nmtran_exe = NULL, delete_on_exit = TRUE, ... ){ nmtran_exe <- locate_nmtran(.mod, .config_path, nmtran_exe) nm_ver <- attr(nmtran_exe, "nonmem_version") mod_path <- get_model_path(.mod) data_path <- get_data_path_from_ctl(.mod) # make temporary directory in current directory mod_name <- fs::path_ext_remove(basename(mod_path)) tempdir0 <- paste0("nmtran_", mod_name, "_", basename(tempdir())) dir.create(tempdir0) if(isTRUE(delete_on_exit)){ on.exit(unlink(tempdir0, recursive = TRUE, force = TRUE)) } # copy model and dataset file.copy(mod_path, tempdir0) file.copy(data_path, tempdir0) # overwrite $DATA of new model modify_data_path_ctl( mod_path = file.path(tempdir0, basename(mod_path)), data_path = basename(data_path) ) # Get command & append the control file name cmd <- stringr::str_glue(nmtran_exe, .envir = list(ctl_name = basename(mod_path)), .na = NULL) cmd <- paste(cmd, "<", basename(mod_path)) # Run NMTRAN if(!is.null(nm_ver)){ message(glue("Running NMTRAN with NONMEM version `{nm_ver}`")) } system_nm(cmd, dir = tempdir0, wait = TRUE, ...) if(isFALSE(delete_on_exit)){ return(tempdir0) } } locate_nmtran <- function(.mod, .config_path = NULL, nmtran_exe = NULL){ if(is.null(nmtran_exe)){ model_dir <- get_model_working_directory(.mod) config_path <- .config_path %||% file.path(model_dir, "bbi.yaml") if(!file_exists(config_path)){ stop(paste("No bbi configuration was found in the execution directory.", "Please run `bbi_init()` with the appropriate directory to continue.")) } if (!is.null(.config_path)) { config_path <- normalizePath(.config_path) } bbi_config <- yaml::read_yaml(config_path) nm_config <- bbi_config$nonmem # look for default nonmem installation default_nm <- purrr::keep(nm_config, function(nm_ver){ !is.null(nm_ver$default) }) # Set nonmem path if(length(default_nm) > 0){ default_nm <- default_nm[[1]] }else{ # If no default, use the last one (likely higher version) default_nm <- nm_config[[length(nm_config)]] } # Set NMTRAN executable path nm_path <- default_nm$home nmtran_exe <- file.path(nm_path, "tr", "NMTRAN.exe") # If executable found via bbi.yaml, append NONMEM version as attribute attr(nmtran_exe, "nonmem_version") <- basename(default_nm$home) } if(!file_exists(nmtran_exe)){ stop(glue("Could not find an NMTRAN executable at `{nmtran_exe}`")) } return(nmtran_exe) } system_nm_default <- function(cmd, ...) { if (.Platform$OS.type == "windows") { local_env_vars <- Sys.getenv() stdout_unit_vars <- local_env_vars[grepl("STDOUT_UNIT|STDERR_UNIT", names(local_env_vars))] for (i in seq_along(stdout_unit_vars)) { Sys.unsetenv(names(stdout_unit_vars)[i]) } on.exit({ if (length(stdout_unit_vars) > 0) { do.call(Sys.setenv, as.list(stdout_unit_vars)) } }) args <- list(...) if (!"wait" %in% names(args)) wait <- FALSE else wait <- args$wait if (wait == FALSE) { shell(paste("START CMD /C", cmd), ...) } else { shell(cmd, ...) } } else { system(cmd, ...) } } system_nm <- function(cmd, dir = NULL, ...) { if (is.null(dir) || !file.exists(dir)) dir <- "." if (file.exists(dir)) { currentwd <- getwd() setwd(dir) on.exit(setwd(currentwd)) } else { stop(paste0("Directory \"", dir, "\" doesn't exist.")) } system_nm_default(cmd, ...) } get_data_path_from_ctl <- function(.mod){ mod_path <- get_model_path(.mod) ctl <- nmrec::read_ctl(mod_path) data_rec <- nmrec::select_records(ctl, "data")[[1]] data_path <- nmrec::get_record_option(data_rec, "filename")$value data_path_norm <- fs::path_norm(file.path(mod_path, data_path)) if(!fs::file_exists(data_path_norm)){ stop(glue("Could not find data at {data_path_norm}")) } return(data_path_norm) } modify_data_path_ctl <- function(mod_path, data_path){ # Get data record ctl <- nmrec::read_ctl(mod_path) data_rec <- nmrec::select_records(ctl, "data")[[1]] data_rec$parse() # Overwrite 'filename' option # TODO: confirm this works with .mod extensions data_rec$values <- purrr::map(data_rec$values, function(data_opt){ if(inherits(data_opt, "nmrec_option_pos") && data_opt$name == "filename"){ data_opt$value <- data_path } data_opt }) # Write out modified ctl nmrec::write_ctl(ctl, mod_path) } ```

Example output

> run_nmtran(MOD1)
Running NMTRAN with NONMEM version `nm75`

 WARNINGS AND ERRORS (IF ANY) FOR PROBLEM    1

 (WARNING  2) NM-TRAN INFERS THAT THE DATA ARE POPULATION.

Note: Analytical 2nd Derivatives are constructed in FSUBS but are never used.
      You may insert $ABBR DERIV2=NO after the first $PROB to save FSUBS construction and compilation time

One unfortunate thing is that NMTRAN does require a valid data path to properly execute. While this makes sense, and is there for good reason, this does prevent us from testing this out on a few of our 'complex' models (within bbr), as we do not have the necessary dataset as part of this package.

Something to think about/look into:

barrettk commented 5 months ago

Though I dont imagine we would include this function in bbr, I wrote this helper function for comparing two models (similar to @kyleam's bash script):

compare_nmtran ```r compare_nmtran <- function( .mod, .mod_compare, .config_path = NULL, nmtran_exe = NULL ){ nmtran_exe <- locate_nmtran(.mod, .config_path, nmtran_exe) nmtran_exe2 <- locate_nmtran(.mod_compare, .config_path, nmtran_exe) if(nmtran_exe != nmtran_exe2){ rlang::warn( c( "!" = "Found two separate NMTRAN executables:", " " = paste("-", nmtran_exe), " " = paste("-", nmtran_exe2), "i" = "Defaulting to the first one" ) ) } # Run NMTRAN on each model (only message once) nmtran_mod <- run_nmtran( empty_prob_statement(.mod), delete_on_exit = FALSE, intern = TRUE ) nmtran_compare <- run_nmtran( empty_prob_statement(.mod_compare), delete_on_exit = FALSE, intern = TRUE ) %>% suppressMessages() # Force delete folders at the end on.exit(unlink(nmtran_mod, recursive = TRUE, force = TRUE)) on.exit(unlink(nmtran_compare, recursive = TRUE, force = TRUE), add = TRUE) # Compare FCON files nmtran_mod_fcon <- file.path(nmtran_mod, "FCON") nmtran_compare_fcon <- file.path(nmtran_compare, "FCON") cmd <- paste("cmp", nmtran_mod_fcon, nmtran_compare_fcon) # Warnings occur when files are different output <- system_nm(cmd, intern = TRUE, input = tempdir()) %>% suppressWarnings() if(length(output) == 0){ message("\nNo differences found") }else{ message("\nModels are not equivalent") } return(output) } # This function is used to remove problem statement differences introduced via `copy_model_from` empty_prob_statement <- function(.mod){ mod_path <- get_model_path(.mod) ctl <- nmrec::read_ctl(mod_path) prob_rec <- nmrec::select_records(ctl, "prob")[[1]] prob_rec$parse() # Overwrite 'text' option prob_rec$values <- purrr::map(prob_rec$values, function(prob_rec){ if(inherits(prob_rec, "nmrec_option_pos") && prob_rec$name == "text"){ prob_rec$value <- "" } prob_rec }) # Write out modified ctl nmrec::write_ctl(ctl, mod_path) return(.mod) } ```

Example

# Starting model - set a reference
open_model_file(MOD1)

# Make new model
MOD_COMPARE <- copy_model_from(MOD1)
# Make a change
open_model_file(MOD_COMPARE)

# compare NMTRAN evaluation
compare_nmtran(MOD1, MOD_COMPARE)

# delete new model at the end
delete_models(MOD_COMPARE, .tags = NULL, .force = TRUE)

Example output

> compare_nmtran(MOD1, MOD_COMPARE)
Running NMTRAN with NONMEM version `nm75`

Models are not equivalent
[1] "nmtran_1_RtmpdmX0QB/FCON nmtran_2_RtmpdmX0QB/FCON differ: byte 3595, line 46"
attr(,"status")
[1] 1

Real scenario example

I wanted to test whether diagonal matrices could specify standard deviation for one value, and variance for another

The reference model had this block:

$OMEGA
0.05 STANDARD   ; iiv CL
0.2     ; iiv V2

The new model had this block:

$OMEGA
0.05 STANDARD   ; iiv CL
0.2 VAR    ; iiv V2

I found no differences here, meaning that adding VAR to the second ETA value had no impact. This was an important observation, as it means I can assume one flag each for diagonal and off-diagonal values, per diagonal matrix-type record.

> compare_nmtran(MOD1, MOD_COMPARE)
Running NMTRAN with NONMEM version `nm75`

No differences found
character(0)
kyleam commented 5 months ago

[ just a drive-by comment ... ]

I know we initially discussed using processx, but system seems to work pretty well, and allowed me to port over more code from NMProject.

Please use processx. system() and friends invoke the command through the shell, which means you need to worry about things like path quoting and shell injection. Even if the latter isn't a realistic security concern for the given context, it's better to just avoid it entirely.

kyleam commented 5 months ago

[ another quick comment ]

Based on how nmfe* invokes NMTRAN.exe...

$ grep tr/NMTRAN /opt/NONMEM/nm75/run/nmfe75
# $dir/tr/NMTRAN.exe $prdefault $tprdefault $maxlim < $1 >& FMSG
  $dir/util/nmtran_presort < tempzzzz1 | $dir/tr/NMTRAN.exe $prdefault $tprdefault $maxlim $do2test >& FMSG
  $dir/util/nmtran_presort < tempzzzz1 | $dir/tr/NMTRAN.exe $prdefault $tprdefault $maxlim >& FMSG

... I think it's worth considering these questions:

The answer may be that none of that matters in the context of this check, but that'd still be good to document here.

barrettk commented 5 months ago

Please use processx. system() and friends invoke the command through the shell, which means you need to worry about things like path quoting and shell injection. Even if the latter isn't a realistic security concern for the given context, it's better to just avoid it entirely

FWIW I had a lot of trouble trying to use processx. An example command could be: /opt/NONMEM/nm75/tr/NMTRAN.exe < 1.ctl

When experimenting with various processx::run() configurations, I couldn't get it to not quote the arguments. Here are a couple examples I would return:

/opt/NONMEM/nm75/tr/NMTRAN.exe '< 1.ctl' # when treating it as a single arg (args = '<  1.ctl')
/opt/NONMEM/nm75/tr/NMTRAN.exe '<' 1.ctl # when making '<' and '1.ctl' two separate args (args = c('<', '1.ctl'))
'/opt/NONMEM/nm75/tr/NMTRAN.exe <' 1.ctl # when making '<' part of the command
/opt/NONMEM/nm75/tr/NMTRAN.exe < '1.ctl' # I think I produced this one as well, but am blanking on the configuration

In this case, < isn't really an argument, so I think that's part of the problem (though I couldnt make it part of the command either).

I agree with the motivation to use processx, but do you have any suggestions for how I could either A) reproduce the command /opt/NONMEM/nm75/tr/NMTRAN.exe < 1.ctl (with no quotes), or B) find an alternative to the < command?

kyleam commented 5 months ago

prog <foo is shell syntax that sends the contents of foo as standard input to prog. processx doesn't go through a shell. See processx's stdin arg.

barrettk commented 5 months ago

See processx's stdin arg.

I had looked into that, but will give it another look. Thanks