R/UpdateMDLfrom_dmptxt.R

In Certara.NLME8: Utilities for Certara's Nonlinear Mixed-Effects Modeling Engine

#' Update Model text file from NLME output File
#'
#' This function updates a model file with parameter estimates obtained from a
#' dmp file (R structure format of output generated by NLME) text file. The
#' updated model file includes the estimated fixed effects, error terms and
#' random effects values.
#'
#' @param dmpfile The path to the DMP text file.
#' @param SharedWorkingDir The working directory. Used if `dmpfile`,
#'   `model_file`, `output_file` are given without path.
#' @param model_file The name of the model file to be updated (with optional
#'   full path).
#' @param compile A logical value indicating whether to compile the updated
#'   model file into NLME executable. Default is `TRUE`, it also overwrites
#'   `model_file` with updated estimates (i.e. making the same as
#'   `output_file`.)
#' @param output_file The name of the new model file with updated estimates.
#'
#' @details `TDL5` executable from NLME Engine is used. NLME engine
#' location is identified by `INSTALLDIR` environment variable. The current
#' function will give an error if `TDL5` cannot be executed.
#'
#' @return The path to the updated model file.
#'
#' @md
#' @export
UpdateMDLfrom_dmptxt <-
  function(dmpfile = "dmp.txt",
           SharedWorkingDir = getwd(),
           model_file = "test.mdl",
           compile = TRUE,
           output_file = "test.mdx") {
    if (dirname(model_file) == ".") {
      model_file <- file.path(SharedWorkingDir, model_file)
    }

    model_file <- normalizePath(model_file)

    # dmp has thetas omegas and sigma
    if (dirname(dmpfile) == ".") {
      dmpfile <- file.path(SharedWorkingDir, dmpfile)
    }

    dmpfile <- normalizePath(dmpfile)

    if (dirname(output_file) == ".") {
      output_file <- file.path(SharedWorkingDir, output_file)
    }

    output_file <- normalizePath(output_file, mustWork = FALSE)

    # trying to get the header
    mdlfile <- paste(readLines(model_file,
                               n = -1,
                               warn = FALSE),
                     collapse = "\n")
    # remove comments from header
    OneLine2Slash <- "(?:\\/\\/(?:\\\\\\n|[^\\n])*(?=$|\\n))"
    OneLineSharp <- "(?:#(?:\\\\\\n|[^\\n])*(?=$|\\n))"
    Asterisk <- "(?:\\/\\*[\\s\\S]*?\\*\\/)"
    pattern <- paste(OneLine2Slash,
                     OneLineSharp,
                     Asterisk,
                     sep = "|",
                     collapse = "|")

    fileWOComments <- gsub(pattern, "\n", mdlfile, perl = TRUE)
    mdlheader <-
      substr(fileWOComments, 1, regexpr("{", fileWOComments, fixed = TRUE)[1])


    # find forzen fixefs and enables
    # if hash is given use it
    NLME_HASH <- as.integer(Sys.getenv("NLME_HASH"))
    if (!is.na(NLME_HASH) && NLME_HASH > 0) {
      NLME_HASH <- paste("-hash", NLME_HASH)
    } else {
      NLME_HASH <- character(0)
    }

    dirnameModel <- dirname(model_file)

    TDL5_run <-
      system2(
        path.expand(file.path(
          Sys.getenv("INSTALLDIR"),
          Sys.getenv("PML_BIN_DIR") ,
          "TDL5"
        )),
        args = c(
          NLME_HASH,
          "-i",
          shQuote(model_file, type = "cmd"),
          shQuote(dirnameModel, type = "cmd")
        ),
        stdout = TRUE
      )

    ModelInfoPath <- file.path(dirnameModel, "ModelInfo.txt")

    if (!file.exists(ModelInfoPath)) {
      stop(
        paste0(
          "Cannot create file for mapping:\n",
          ModelInfoPath,
          "\nTDL5 output:\n",
          paste0(TDL5_run, collapse = "\n")
        )
      )
    }

    ModelInfo <- readLines(ModelInfoPath, warn = FALSE)
    fixefLine <-
      ModelInfo[grepl("^\\(fixedeffects ", ModelInfo)]
    fixefStrings <-
      unlist(strsplit(fixefLine, split = "(\\(fixedeffects\\W*\\()|(\\)\\W*\\()|(\\)\\))"))
    fixefStrings <- fixefStrings[fixefStrings != ""]
    splittedFixefsStrings <-
      strsplit(fixefStrings, split = " ")
    fixefInfoDF <- data.frame(splittedFixefsStrings)
    colnames(fixefInfoDF) <- fixefInfoDF[1, ]
    fixefInfoDF <- fixefInfoDF[-c(1), ]
    rownames(fixefInfoDF) <- c("enable", "freeze", "value")

    # get rid of residuals and read results
    dmpLines <- shrinkDmpDotTxt(dmpfile)
    source(textConnection(dmpLines), local = TRUE)

    # write model with override
    model_file_override <- paste0(tools::file_path_sans_ext(model_file),
                                "Override.",
                                tools::file_ext(model_file))
    file.copy(from = model_file,
              to = model_file_override,
              overwrite = TRUE)

    cat(
      "override",
      mdlheader,
      "\n",
      file = model_file_override,
      sep = "\n",
      append = TRUE
    )

    epsilon <- dmp.txt$sigma
    fixedEffects <- dmp.txt$coefficients$fixed

    if (length(epsilon) > 0) {
      for (epsilonName in colnames(epsilon)) {
        cat(
          "   \nerror(",
          epsilonName,
          " = ",
          sprintf("%.15g", fixedEffects[epsilonName]),
          ")\n",
          file = model_file_override,
          append = TRUE
        )
      }
    }

    # if there is a sigma, get rid of it since it is written as epsilon
    if (length(epsilon) > 0) {
      head(fixedEffects, -1)
    }

    fixedEffectNames <- names(fixedEffects)
    if (length(fixedEffects) > 0) {
      fixefString <- ""
      for (ifixef in seq_along(fixedEffects)) {
        EnableFreezeFlag <- character(0)
        FreezeFlag <-
          as.numeric(fixefInfoDF["freeze", fixedEffectNames[ifixef]])
        EnableFlag <-
          as.numeric(fixefInfoDF["enable", fixedEffectNames[ifixef]])
        if (fixedEffectNames[ifixef] %in% colnames(fixefInfoDF)) {
          if (FreezeFlag > 0) {
            # frozen
            EnableFreezeFlag <- "(freeze)"
          } else if (EnableFlag >= 0) {
            EnableFreezeFlag <- paste0("(enable = c(", EnableFlag, "))")
          }
        }

        fixefString <- paste0(
          fixefString,
          "   \nfixef(",
          fixedEffectNames[ifixef],
          EnableFreezeFlag,
          " = c(,",
          sprintf("%.15g", fixedEffects[ifixef]),
          ",))\n"
        )
      }

      cat(fixefString, file = model_file_override, append = TRUE)
    }

    randomEffects <- data.frame(dmp.txt$omega)
    randomEffectsNames <- colnames(randomEffects)
    if (length(randomEffects) > 0) {
      ranefString <-
        paste0("   \nranef(block(",
               paste(randomEffectsNames, collapse = ", "),
               ") = c(")

      for (iRanefRow in 1:length(dimnames(dmp.txt$omega)[[1]])) {
        for (iRanefCol in 1:iRanefRow) {
          if (iRanefCol + iRanefRow > 2) {
            ranefString <- paste0(ranefString, ", ")
          }

          ranefString <- paste0(ranefString,
                                sprintf("%.8g", dmp.txt$omega[iRanefRow, iRanefCol]))
        }
      }

      ranefString <- paste0(ranefString, "))")
      cat(ranefString, file = model_file_override, append = TRUE)
    }

    cat("}", file = model_file_override, append = TRUE)

    # mdx is stored initially in the same dir as original model
    TDL5_run <-
      system2(
        path.expand(file.path(
          Sys.getenv("INSTALLDIR"),
          Sys.getenv("PML_BIN_DIR"),
          "TDL5"
        )),
        args = c(
          NLME_HASH,
          "-r",
          shQuote(model_file_override, type = "cmd"),
          shQuote(dirnameModel, type = "cmd")
        ),
        stdout = TRUE
      )

    newModelFile <- gsub("\\w$", "x", model_file_override)
    if (!file.exists(newModelFile)) {
      stop(paste0(
        "Cannot update estimates in model file.",
        "\nTDL5 output:\n",
        paste0(TDL5_run, collapse = "\n")
      ))
    }

    if (newModelFile != output_file) {
      file.copy(from = newModelFile,
                to = output_file,
                overwrite = TRUE)
    }

    # compiling the resulted model if requested
    if (compile) {
      file.rename(model_file, paste0(model_file, "_old"))
      file.rename(newModelFile, model_file)
      newModelFile <- model_file
      generateNLMEScriptAndRun("COMPILE")
    }

    output_file
  }