R/ui.R

Defines functions .nlmixrUpdate model ini .processLotri .getUif .getPipedData .clearPipedData .deparse1 .bodyDewrap .deparse nlmixrfindRhsLhs nlmixrfindLhs nlmixrfindRhs

Documented in ini model nlmixrfindLhs

nlmixrfindRhs <- function(x) {
  ## Modified from http://adv-r.had.co.nz/Expressions.html find_assign4
  if (is.atomic(x)) {
    character()
  } else if (is.name(x)) {
    return(as.character(x))
  } else if (is.call(x)) {
    if ((identical(x[[1]], quote(`<-`)) ||
      identical(x[[1]], quote(`=`))) &&
        is.name(x[[2]])) {
      unique(c(unlist(nlmixrfindRhs(x[[3]]))))
    } else {
      x1 <- x[-1]
      unique(unlist(lapply(x1, nlmixrfindRhs)))
    }
  } else if (is.pairlist(x)) {
    unique(unlist(lapply(x, nlmixrfindRhs)))
  } else {
    stop("Don't know how to handle type ", typeof(x),
         call. = FALSE
         )
  }
}

nlmixrfindLhs <- function(x) {
  ## Modified from http://adv-r.had.co.nz/Expressions.html find_assign4
  if (is.atomic(x) || is.name(x)) {
    character()
  } else if (is.call(x)) {
    if ((identical(x[[1]], quote(`<-`)) ||
      identical(x[[1]], quote(`=`))) &&
      is.name(x[[2]])) {
      lhs <- as.character(x[[2]])
    } else {
      lhs <- character()
    }
    unique(c(lhs, unlist(lapply(x, nlmixrfindLhs))))
  } else if (is.pairlist(x)) {
    unique(unlist(lapply(x, nlmixrfindLhs)))
  } else {
    stop("Don't know how to handle type ", typeof(x),
      call. = FALSE
    )
  }
}

nlmixrfindRhsLhs <- function(x) {
  .lhs <- nlmixrfindLhs(x)
  .rhs <- nlmixrfindRhs(x)
  .rhs <- setdiff(.rhs, .lhs)
  list(lhs=.lhs, rhs=.rhs)
}

.deparse <- function(expr) {
  deparse(expr, width.cutoff = 500, control = "useSource")
}

.bodyDewrap <- function(ret) {
  .ret <- ret
  if (length(.ret) > 1) {
    .ret[1] <- sub(rex::rex(
      start, or(group("function", any_spaces, "(", any_spaces, ")", any_spaces), ""),
      any_spaces, any_of("("), any_spaces, "{", any_spaces
    ), "", .ret[1], perl = TRUE)
    if (.ret[1] == "") .ret <- .ret[-1]
  }
  if (length(.ret) > 1) {
    .len <- length(.ret)
    .ret[.len] <- sub(rex::rex(any_spaces, "}", any_spaces, any_of(")"), any_spaces, end), "", .ret[.len])
    if (.ret[.len] == "") .ret <- .ret[-.len]
  }
  return(.ret)
}

.deparse1 <- function(expr) {
  return(.bodyDewrap(deparse(expr, width.cutoff = 500)))
}

.pipedData <- NULL
.clearPipedData <- function() {
  assignInMyNamespace(".pipedData", NULL)
}

.getPipedData <- function() {
  .pipedData
}

.getUif <- function(what) {
  if (inherits(what, "nlmixrFitCore")) {
    .uif <- what$uif
    assignInMyNamespace(".pipedData", getData(what))
  } else if (inherits(what, "nlmixrUI")) {
    .uif <- what
  } else if (inherits(what, "function")) {
    .uif <- nlmixrUI(what)
  } else {
    stop("Do not know how to handle object")
  }
  return(.uif)
}

.processLotri <- function(.code, .uif) {
  .mat <- try(eval(parse(text = sprintf("lotri::lotri(%s)", .code))), silent = TRUE)
  if (inherits(.mat, "matrix")) {
    .d <- dimnames(.mat)[[1]]
    .ini2 <- .as.data.frame(.uif$ini)
    .ini1 <- .ini2[is.na(.ini2$neta1), ]
    .ini2 <- .ini2[!is.na(.ini2$neta1), ]
    .ini4 <- .ini2[!is.na(.ini2$err), ]
    .ini2 <- .ini2[is.na(.ini2$err), ]
    .m2 <- as.vector(.mat)
    .l2 <- length(.d) * 2
    .df <- .data.frame(n1 = rep(.d, each = length(.d)), n2 = rep(.d, length(.d)), val = .m2)
    .dfI <- .as.data.frame(.uif$ini[, c("neta1", "neta2", "name")])
    .dfI <- .dfI[!is.na(.dfI$neta1), ]
    .dfI <- .dfI[which(.dfI$neta1 == .dfI$neta2), c("neta1", "name")]
    .d2 <- paste(.dfI$name)
    .diff <- setdiff(.d, .d2)
    if (length(.diff) > 0) {
      stop(sprintf(
        "trying to provide an estimate for non-existant eta: %s",
        paste(.diff, collapse = ", ")
      ))
    }
    names(.dfI)[2] <- "n1"
    .df$n1 <- paste(.df$n1)
    .dfI$n1 <- paste(.dfI$n1)
    .df <- merge(.df, .dfI, all.x = TRUE)
    names(.dfI) <- c("neta2", "n2")
    .df <- merge(.df, .dfI, all.x = TRUE)
    .df <- .df[order(.df$neta1, .df$neta2), ]
    ## Name becomes (eta.cl,eta.ka)
    .df$name <- ifelse(.df$neta1 == .df$neta2, .df$n1, paste0("(", .df$n1, ",", .df$n2, ")"))
    .df$name2 <- ifelse(.df$neta1 == .df$neta2, .df$n1, paste0("(", .df$n2, ",", .df$n1, ")"))
    .ini3 <- do.call("rbind", lapply(seq_along(.df$name), function(i) {
      .name1 <- paste(.df$name[i])
      .name2 <- paste(.df$name2[i])
      .w <- which(.name1 == .ini2$name)
      if (length(.w) == 1) {
        .ini2$est[.w] <<- .df$val[i]
        return(NULL)
      }
      .w <- which(.ini2$name == .name2)
      if (length(.w) == 1) {
        .ini2$est[.w] <<- .df$val[i]
        return(NULL)
      }
      if (.df$neta1[i] < .df$neta2[i]) {
        return(NULL)
      }
      .df2 <- nlmixrBoundsTemplate
      .df2$neta1 <- .df$neta1[i]
      .df2$neta2 <- .df$neta2[i]
      .df2$name <- .df$name[i]
      .df2$est <- .df$val[i]
      .df2$lower <- -Inf
      .df2$upper <- Inf
      .df2$fix <- FALSE
      .df2$condition <- "ID"
      return(.df2)
    }))
    .ini2 <- rbind(.ini2, .ini3)
    .ini2 <- .ini2[order(.ini2$neta1, .ini2$neta2), ]
    .ini2 <- rbind(.ini1, .ini2, .ini4)
    class(.ini2) <- c("nlmixrBounds", "data.frame")
    .uif$ini <- .ini2
    return(.uif)
  } else {
    stop(sprintf("invalid syntax: %s", .code))
  }
}

#' nlmixr ini block handling
#'
#' The ini block controls initial conditions for 'theta' (fixed effects),
#' 'omega' (random effects), and 'sigma' (residual error) elements of the model.
#'
#' 'theta' and 'sigma' can be set using either \code{<-} or \code{=} such as
#' \code{tvCL <- 1} or equivalently \code{tvCL = 1}.  'omega' can be set with a
#' \code{~}.
#'
#' Parameters can be named or unnamed (though named parameters are preferred).
#' A named parameter is set using the name on the left of the assignment while
#' unnamed parameters are set without an assignment operator.  \code{tvCL <- 1}
#' would set a named parameter of \code{tvCL} to \code{1}.  Unnamed parameters
#' are set using just the value, such as \code{1}.
#'
#' For some estimation methods, lower and upper bounds can be set for 'theta'
#' and 'sigma' values.  To set a lower and/or upper bound, use a vector of
#' values.  The vector is \code{c(lower, estimate, upper)}.  The vector may be
#' given with just the estimate (\code{c(estimate)}), the lower bound and
#' estimate (\code{c(lower, estimate)}), or all three (\code{c(lower, estimate,
#' upper)}).  To set an estimate and upper bound without a lower bound, set the
#' lower bound to \code{-Inf}, \code{c(-Inf, estimate, upper)}.  When an
#' estimation method does not support bounds, the bounds will be ignored with a
#' warning.
#'
#' 'omega' values can be set as a single value or as the values of a
#' lower-triangular matrix.  The values may be set as either a
#' variance-covariance matrix (the default) or as a correlation matrix for the
#' off-diagonals with the standard deviations on the diagonals.  Names may be
#' set on the left side of the \code{~}.  To set a variance-covariance matrix
#' with variance values of 2 and 3 and a covariance of -2.5 use \code{~c(2, 2.5,
#' 3)}.  To set the same matrix with names of \code{iivKa} and \code{iivCL}, use
#' \code{iivKa + iivCL~c(2, 2.5, 3)}.  To set a correlation matrix with standard
#' deviations on the diagonal, use \code{cor()} like \code{iivKa + iivCL~cor(2,
#' -0.5, 3)}.
#'
#' Values may be fixed (and therefore not estimated) using either the name
#' \code{fixed} at the end of the assignment or by calling \code{fixed()} as a
#' function for the value to fix.  For 'theta' and 'sigma', either the estimate
#' or the full definition (including lower and upper bounds) may be included in
#' the fixed setting.  For example, the following are all effectively equivalent
#' to set a 'theta' or 'sigma' to a fixed value (because the lower and upper
#' bounds are ignored for a fixed value): \code{tvCL <- fixed(1)}, \code{tvCL <-
#' fixed(0, 1)}, \code{tvCL <- fixed(0, 1, 2)}, \code{tvCL <- c(0, fixed(1),
#' 2)}, or \code{tvCL <- c(0, 1, fixed)}.  For 'omega' assignment, the full
#' block or none of the block must be set as \code{fixed}.  Examples of setting
#' an 'omega' value as fixed are: \code{iivKa~fixed(1)}, \code{iivKa +
#' iivCL~fixed(1, 2, 3)}, or \code{iivKa + iivCL~c(1, 2, 3, fixed)}.  Anywhere
#' that \code{fixed} is used, \code{FIX}, \code{FIXED}, or \code{fix} may be
#' used equivalently.
#'
#' For any value, standard mathematical operators or functions may be used to
#' define the value.  For example, \code{exp(2)} and \code{24*30} may be used to
#' define a value anywhere that a number can be used (e.g. lower bound,
#' estimate, upper bound, variance, etc.).
#'
#' Values may be labeled using the \code{label()} function after the assignment.
#' Labels are are used to make reporting easier by giving a human-readable
#' description of the parameter, but the labels do not have any effect on
#' estimation.  The typical way to set a label so that the parameter \code{tvCL}
#' has a label of "Typical Value of Clearance (L/hr)" is \code{tvCL <- 1;
#' label("Typical Value of Clearance (L/hr)")}.
#'
#' \code{nlmixr} will attempt to determine some back-transformations for the
#' user.  For example, \code{CL <- exp(tvCL)} will detect that \code{tvCL} must
#' be back-transformed by \code{exp()} for easier interpretation.  When you want
#' to control the back-transformation, you can specify the back-transformation
#' using \code{backTransform()} after the assignment.  For example, to set the
#' back-transformation to \code{exp()}, you can use \code{tvCL <- 1;
#' backTransform(exp())}.
#'
#' @param ini Ini block or nlmixr model object
#' @param ... Other arguments parsed by nlmixr
#' @return bounds expression or parsed ui object
#' @author Matthew L. Fidler
#' @export
ini <- function(ini, ...) {
  if (is(substitute(ini), "{")) {
    .f <- eval(parse(text = sprintf("function() %s", paste(.deparse(substitute(ini)), collapse = "\n"))))
    .fb <- nlmixrBounds(.f)
    if (!exists(".ini", parent.frame(1))) {
      assign(".ini", .f, parent.frame(1))
    } else {
      .ini <- get(".ini", parent.frame(1))
      if (is(.ini, "function")) {
        assign(".ini", .f, parent.frame(1))
      }
    }
    return(.fb)
  } else {
    .uif <- .getUif(ini)
    .ini <- .uif$ini
    .call <- match.call(expand.dots = TRUE)[-1]
    .call <- .call[-1]
    .ns <- names(.call)
    if (!is.null(.ns)) {
      .ini <- .uif$ini
      for (.n in .ns) {
        .w <- which(.n == .ini$name)
        if (length(.w) == 1) {
          if (any(.deparse(.call[[.n]]) == c("fix", "fixed", "FIX", "FIXED"))) {
            if (.uif$ini$fix[.w]) {
              warning("trying to fix '", n, "', but already fixed")
            } else {
              .uif$ini$fix[.w] <- TRUE
            }
          } else if (any(.deparse(.call[[.n]]) == c("unfix", "unfixed", "UNFIX", "UNFIXED"))) {
            if (.uif$ini$fix[.w]) {
              .uif$ini$fix[.w] <- FALSE
            } else {
              warning("trying to unfix '", .n, "', but not fixed")
            }
          } else if (regexpr(rex::rex(or(c("fix", "fixed", "FIX", "FIXED")), "(", anything, ")"), .deparse(.call[[.n]])) != -1) {
            .val <- eval(.call[[.n]][[2]])
            if (.uif$ini$fix[.w]) {
              warning("trying to fix '", .n, "', but already fixed, still assigned to '", .val, "'")
            } else {
              .uif$ini$fix[.w] <- TRUE
            }
            .uif$ini$est[.w] <- .val
          } else if (regexpr(rex::rex(or(c("unfix", "unfixed", "UNFIX", "UNFIXED")), "(", anything, ")"), .deparse(.call[[.n]])) != -1) {
            .val <- eval(.call[[.n]][[2]])
            if (.uif$ini$fix[.w]) {
              .uif$ini$fix[.w] <- FALSE
            } else {
              warning("trying to unfix '", .n, "', but not fixed, still assigned to '", .val, "'")
            }
            .uif$ini$est[.w] <- .val
          } else {
            .val <- try(eval(.call[[.n]]), silent = TRUE)
            .found <- TRUE
            if (inherits(.val, "try-error")) {
              .found <- FALSE
              .val <- as.character(.call[[.n]])
              .val2 <- .val
              .frames <- seq(1, sys.nframe())
              .frames <- .frames[.frames != 0]
              for (.f in .frames) {
                .env <- parent.frame(.f)
                if (exists(.val, envir = .env)) {
                  .val2 <- try(get(.val, envir = .env), silent = TRUE)
                  if (!inherits(.val2, "try-error")) {
                    .val <- .val2
                    .found <- TRUE
                    break
                  }
                }
              }
            }
            if (!.found) {
              stop(sprintf("object '%s' not found", .val))
            }
            if (length(.val) == 1) {
              .uif$ini$est[.w] <- .val
            } else if (length(.val) == 2) {
              .uif$ini$lower[.w] <- .val[1]
              .uif$ini$est[.w] <- .val[2]
              ## Warning here? The upper should be inf?
              .uif$ini$est[.w] <- Inf
              .uif$ini$fix[.w] <- FALSE
            } else if (length(.val) == 3) {
              .uif$ini$lower[.w] <- .val[1]
              .uif$ini$est[.w] <- .val[2]
              ## Warning here? The upper should be inf?
              .uif$ini$upper[.w] <- .val[3]
              .uif$ini$fix[.w] <- FALSE
            } else {
              stop(sprintf("Cannot figure out what you are trying to do to the '%s' estimate.", .n))
            }
          }
        } else {
          warning("the model does not have a parameter named '", .n, "', modification ignored")
        }
      }
    } else if (length(.call) == 1) {
      .lst <- eval(.call[[1]])
      .ns <- names(.lst)
      if (inherits(.lst, "list") && !is.null(.ns)) {
        for (.n in .ns) {
          .w <- which(.n == .ini$name)
          if (length(.w) == 1) {
            .val <- .lst[[.n]]
            if (length(.val) == 1) {
              .uif$ini$est[.w] <- .val
            } else if (length(.val) == 2) {
              .uif$ini$lower[.w] <- .val[1]
              .uif$ini$est[.w] <- .val[2]
              .uif$ini$upper[.w] <- Inf
              .uif$ini$fix[.w] <- FALSE
            } else if (length(.val) == 3) {
              .uif$ini$lower[.w] <- .val[1]
              .uif$ini$est[.w] <- .val[2]
              .uif$ini$upper[.w] <- .val[3]
              .uif$ini$fix[.w] <- FALSE
            } else {
              stop(sprintf("Cannot figure out what to do with '%s'", .n))
            }
          }
        }
      } else if (inherits(.lst, "numeric") && !is.null(.ns)) {
        for (.n in .ns) {
          .w <- which(.n == .ini$name)
          if (length(.w) == 1) {
            .uif$ini$est[.w] <- .lst[.n]
          }
        }
      } else {
        ## If this is a+b~c(...) parse using lotri
        ## FIXME handle conditional changes like a+b~c(...) | occ
        .code <- paste(deparse(.lst), collapse = " ")
        .uif <- .processLotri(.code, .uif)
      }
    } else {
      stop("Do not know what to do with the model's ini call.")
    }
    return(.uif)
  }
}


.thetamodelVars <- rex::rex(or("tv", "t", "pop", "POP", "Pop", "TV", "T", "cov", "err", "eff"))
.thetaModelReg <- rex::rex(or(
  group(start, .thetamodelVars),
  group(.thetamodelVars, end)
))
.etaParts <- c(
  "eta", "ETA", "Eta", "ppv", "PPV", "Ppv", "iiv", "Iiv", "bsv", "Bsv", "BSV",
  "bpv", "Bpv", "BPV", "psv", "PSV", "Psv"
)
.etaModelReg <- rex::rex(or(group(start, or(.etaParts)), group(or(.etaParts), end)))

##' nlmixr model block
##'
##' @param model Model specification
##' @param ... Other arguments to model object parsed by nlmixr
##' @param .lines This is an internal argument when code{model} is
##'     being called recursively and should not be used.
##' @return Parsed UI object
##' @author Matthew L. Fidler
##' @export
model <- function(model, ..., .lines = NULL) {
  if (is(substitute(model), "{")) {
    .f <- eval(parse(text = sprintf("function() %s", paste(.deparse(substitute(model)), collapse = "\n"))))
    if (!exists(".model", parent.frame(1))) {
      assign(".model", .f, parent.frame(1))
    } else {
      .model <- get(".model", parent.frame(1))
      if (is(.model, "function")) {
        assign(".model", .f, parent.frame(1))
      }
    }
    if (exists(".ini", parent.frame(1))) {
      .ini <- get(".model", parent.frame(1))
      .model <- get(".model", parent.frame(1))
      if (is(.ini, "function") & is(.model, "function")) {
        return(nlmixrUI(parent.frame(1)))
      }
    }
    return(.f)
  } else {
    .uif <- .getUif(model)
    .ini <- .as.data.frame(.uif$ini)
    .call <- match.call(expand.dots = TRUE)[-(1:2)]
    if (is.null(.lines)) {
      .lines <- .deparse(.call[[1]])
    }
    .f <- try(eval(parse(text = sprintf("function() %s", paste(.lines, collapse = "\n")))), silent = TRUE)
    .origLines <- strsplit(.uif$fun.txt, "\n")[[1]]
    .fNew <- c("function(){", .origLines, "}")
    .origLines <- .fNew
    if (!inherits(.f, "try-error")) {
      .fOrig <- eval(parse(text = paste(.fNew, collapse = "\n")))
      .lhsOrig <- nlmixrfindLhs(body(.fOrig))
      .lhsNew <- nlmixrfindLhs(body(.f))
      .shared <- intersect(.lhsNew, .lhsOrig)
      .distOrig <- .findDist(body(.fOrig))
      .distNew <- .findDist(body(.f))
      .sharedDist <- intersect(.distNew, .distOrig)
      if (length(.shared) == 0 && length(.sharedDist) == 0) {
        .code <- paste(.lines, collapse = " ")
        .processLotri(.code, .uif)
        .uif <- try(.processLotri(.code, .uif), silent = TRUE)
        if (inherits(.uif, "try-error")) {
          stop("Could not find a part of the model to modify.")
        }
        return(.uif)
      }
      if (length(.sharedDist) > 0) {
        for (.v in .sharedDist) {
          .regShared <- rex::rex(start, any_spaces, .v, any_spaces, "~")
          .newLine <- .lines[regexpr(.regShared, .lines) != -1]
          .oldLine <- .origLines[regexpr(.regShared, .origLines) != -1]
          .w <- which(regexpr(.regShared, .fNew) != -1)
          if (length(.w) > 1) stop("Cannot modify a multi-line definition")
          .fNew[.w] <- .newLine
          .oldVarsL <- allVars(body(eval(parse(text = sprintf("function(){%s}", .oldLine)))))
          .newVarsL <- allVars(body(eval(parse(text = sprintf("function(){%s}", .newLine)))))
          .rmVars <- setdiff(.oldVarsL, .newVarsL)
          .addVars <- setdiff(.newVarsL, .oldVarsL)
          for (.rm in .rmVars) {
            .wh <- .ini[.ini$name == .rm, ]
            if (length(.wh$name) == 1) {
              if (is.na(.wh$ntheta)) {
                ## removing eta
                .curEta <- .wh$neta1
                .maxEta <- max(.ini$neta1, na.rm = TRUE)
                .s <- seq(.curEta, .maxEta)
                .s <- .s[-length(.s)]
                .w <- unique(sort(c(
                  which(is.na(.ini$neta1)),
                  which(.ini$neta1 != .curEta),
                  which(.ini$neta2 != .curEta)
                )))
                .ini <- .ini[.w, ]
                for (.rmI in .s) {
                  .ini$neta1[.ini$neta1 == .rmI + 1] <- .rmI
                  .ini$neta2[.ini$neta2 == .rmI + 1] <- .rmI
                }
              } else {
                .curTheta <- .wh$ntheta
                .maxTheta <- max(.ini$ntheta, na.rm = TRUE)
                .s <- seq(.curTheta, .maxTheta)
                .s <- .s[-length(.s)]
                .w <- unique(sort(c(
                  which(is.na(.ini$ntheta)),
                  which(.ini$ntheta != .curTheta)
                )))
                .ini <- .ini[.w, ]
                for (.rmI in .s) {
                  .ini$ntheta[.ini$ntheta == .rmI + 1] <- .rmI
                }
              }
            }
          }
          ## Now add variables; Currently assume thetas; Perhaps something better?
          ## Something other than thetas are not currently supported by UI currently.
          for (.new in .addVars) {
            if (!any(.ini$name == .new)) {
              .maxTheta <- max(.ini$ntheta, na.rm = TRUE)
              .d2 <- nlmixrBoundsTemplate
              .d2$name <- .new
              .d2$ntheta <- .maxTheta + 1
              .d2$est <- 1
              .d2$lower <- -Inf
              .d2$upper <- Inf
              .d2$fix <- FALSE
              .ini <- rbind(.ini, .d2)
            }
          }
        }
      }
      if (length(.shared) > 0) {
        for (.v in .shared) {
          .regShared <- rex::rex(start, any_spaces, .v, any_spaces, or("=", "<-"))
          .newLine <- .lines[regexpr(.regShared, .lines) != -1]
          .oldLine <- .origLines[regexpr(.regShared, .origLines) != -1]
          .w <- which(regexpr(.regShared, .fNew) != -1)
          if (length(.w) > 1) stop("Cannot modify a multi-line definition")
          .fNew[.w] <- .newLine
          .oldVarsL <- allVars(body(eval(parse(text = sprintf("function(){%s}", .oldLine)))))
          .newVarsL <- allVars(body(eval(parse(text = sprintf("function(){%s}", .newLine)))))
          .rmVars <- setdiff(.oldVarsL, .newVarsL)
          .addVars <- setdiff(.newVarsL, .oldVarsL)
          for (.rm in .rmVars) {
            .wh <- .ini[.ini$name == .rm, ]
            if (length(.wh$name) == 1) {
              if (is.na(.wh$ntheta)) {
                ## removing eta
                .curEta <- .wh$neta1
                .maxEta <- max(.ini$neta1, na.rm = TRUE)
                .s <- seq(.curEta, .maxEta)
                .s <- .s[-length(.s)]
                .w <- unique(sort(c(
                  which(is.na(.ini$neta1)),
                  which(.ini$neta1 != .curEta),
                  which(.ini$neta2 != .curEta)
                )))
                .ini <- .ini[.w, ]
                for (.rmI in .s) {
                  .ini$neta1[.ini$neta1 == .rmI + 1] <- .rmI
                  .ini$neta2[.ini$neta2 == .rmI + 1] <- .rmI
                }
              } else {
                .curTheta <- .wh$ntheta
                .maxTheta <- max(.ini$ntheta, na.rm = TRUE)
                .s <- seq(.curTheta, .maxTheta)
                .s <- .s[-length(.s)]
                .w <- unique(sort(c(
                  which(is.na(.ini$ntheta)),
                  which(.ini$ntheta != .curTheta)
                )))
                .ini <- .ini[.w, ]
                for (.rmI in .s) {
                  .ini$ntheta[.ini$ntheta == .rmI + 1] <- .rmI
                }
              }
            }
          }
          ## Now add variables; Currently parsed based on variable name; Perhaps something better?
          for (.new in .addVars) {
            if (!any(.ini$name == .new)) {
              if (regexpr(.etaModelReg, .new) != -1) {
                .maxEta <- suppressWarnings(max(.ini$neta1, na.rm = TRUE))
                if (is.infinite(.maxEta)) .maxEta <- 0
                .maxTheta <- max(.ini$ntheta, na.rm = TRUE)
                .d2 <- nlmixrBoundsTemplate
                .d2$neta1 <- .maxEta + 1
                .d2$neta2 <- .maxEta + 1
                .d2$est <- 1
                .d2$lower <- -Inf
                .d2$upper <- Inf
                .d2$fix <- FALSE
                .d2$name <- .new
                .d2$condition <- "ID"
                .ini <- rbind(.ini, .d2)
              } else if (regexpr(.thetaModelReg, .new) != -1) {
                .maxTheta <- max(.ini$ntheta, na.rm = TRUE)
                .d2 <- nlmixrBoundsTemplate
                .d2$ntheta <- .maxTheta + 1
                .d2$est <- 1
                .d2$lower <- -Inf
                .d2$upper <- Inf
                .d2$fix <- FALSE
                .d2$name <- .new
                .ini <- rbind(.ini, .d2)
              }
            }
          }
        }
      }

      class(.ini) <- c("nlmixrBounds", "data.frame")
      .model <- eval(parse(text = paste(.fNew, collapse = "\n"), keep.source = TRUE))
      return(.finalizeUiModel(
        nlmixrUIModel(.model, .ini, NULL),
        as.list(.uif$meta)
      ))
    }
  }
}

.nlmixrUpdate <- function(object, ...) {
  .uif <- .getUif(object)
  .iniNames <- paste(.uif$ini$name)
  .call <- match.call(expand.dots = TRUE)[-c(1:2)]
  .callNames <- names(.call)
  if (is.null(.callNames)) {
    .callNames <- rep("", length(.call))
  }
  if (any(.callNames %in% c("data", "est", "control", "table", "save"))) {
    stop("Cannot update these arguments: 'data', 'est', 'control', 'table', or 'save'\nUse these options with nlmixr.")
  }
  ## Any time the parameters are named, use ini call
  .iniArgs <- which(.callNames %in% .iniNames)
  .iniNames2 <- .callNames[.iniArgs]
  if (length(.iniArgs) > 0) {
    .uif <- do.call(ini, c(list(.uif), setNames(as.list(.call[.iniArgs]), .iniNames2)))
    .iniI <- seq_along(.call)[-.iniArgs]
  } else {
    .iniI <- seq_along(.call)
  }
  for (.i in .iniI) {
    ## Call model() multiple times if needed.
    if (.callNames[.i] != "") {
      ## construct model(.uif,{name()<-expr()})
      .cn <- eval(parse(text = paste0("quote(", .callNames[.i], ")")))
      .uif <- eval(as.call(list(
        `model`, quote(.uif),
        as.call(list(
          quote(`{`),
          as.call(list(
            quote(`<-`), .cn,
            .call[[.i]]
          ))
        ))
      )))
    } else {
      .x <- .call[[.i]]
      .didCall <- FALSE
      if (is.call(.x) || is.pairlist(.x)) {
        if (identical(.x[[1]], quote(`c`)) ||
          identical(.x[[1]], quote(`list`))) {
          .uif <- eval(as.call(list(`ini`, quote(.uif), .x)))
          .didCall <- TRUE
        }
      } else if (is.name(.x)) {
        .cur <- 1
        .tmp <- eval(.x, parent.frame(1))
        if (is(.tmp, "list") || is(.tmp, "numeric") || is(.tmp, "integer")) {
          .uif <- eval(bquote(ini(.(.uif), .(.tmp))))
        } else {
          .uif <- model(.uif, .lines = .deparse(.tmp))
        }
        .didCall <- TRUE
      }
      if (!.didCall) {
        .uif <- model(.uif, .lines = .deparse(.x))
      }
    }
  }
  return(.uif)
}

##' @export
update.nlmixrUI <- .nlmixrUpdate

##' @export
update.nlmixrFitCore <- .nlmixrUpdate

##' @export
update.function <- .nlmixrUpdate

.finalizeUiModel <- function(fun2, meta) {
  class(fun2) <- "nlmixrUI"
  var.ini <- c(fun2$focei.names, fun2$eta.names)
  var.def <- fun2$all.vars
  diff <- setdiff(var.ini, var.def)
  if (length(diff) > 0) {
    stop(sprintf("Model error: initial estimates provided without variables being used: %s", paste(diff, collapse = ", ")))
  }
  ns <- fun2$name[which(!is.na(fun2$neta1) & !is.na(fun2$err))]
  if (length(ns) > 0) {
    stop("residual error component(s) need to be defined with assignment ('=' or '<-') in ini block (not '~'): ",
      paste(ns, collapse = ", "),
      call. = FALSE
    )
  }
  ns <- fun2$name[is.na(fun2$est)]
  if (length(ns) > 0) {
    ## Now handled by bounds (in theory), keeping here just in case.
    if (!all(is.na(ns))) {
      stop(sprintf("the following parameters initial estimates are NA: %s", paste(ns, collapse = ", ")),
        call. = FALSE
      ) # nocov
    }
  }
  fun2$meta <- list2env(meta, parent = emptyenv())
  .mv <- RxODE::rxModelVars(paste(c(.deparse1(body(fun2$focei.fun1)),
    .deparse1(body(fun2$pred)), fun2$extra
  ), collapse = "\n"))
  .tmp <- paste(fun2$nmodel$predDf$var)
  .testVars <- c(.mv$lhs, .mv$state)
  .tmp <- .tmp[!(.tmp %in% c(.testVars, "nlmixr_lincmt_pred"))]
  if (length(.tmp > 0)) {
    .predDf <- fun2$predDf
    if (length(.predDf$var) > 1) {
      .predDf <- .predDf[.predDf$var %in% .tmp, ]
      .predDf <- .predDf[.predDf$var == .predDf$cond, , drop = FALSE]
      if (length(.predDf$var) > 0) {
        stop("multiple compartment models with expressions need to be conditioned by `|`\n",
          "ie log(cp) ~ add(err) | cmt\n",
          "The following endpoints need to be corrected: ",
          paste(.predDf$var, collapse = ", "),
          call. = FALSE
        )
      }
    }
    .tmp <- lapply(.tmp, function(x) {
      .newPars <- RxODE::rxModelVars(paste0("nlmixr_tmp=", x))$params
      if (all(.newPars %in% .testVars)) {
        return(NA_character_)
      }
      return(.newPars)
    })
    .tmp <- stats::na.omit(unlist(.tmp))
    if (length(.tmp) > 0) {
      stop(sprintf(
        "Modeled responses need to be defined in the model; Add definition for: %s",
        paste(.tmp, collapse = ", ")
      ))
    }
  }
  return(fun2)
}

##' Prepares the UI function and returns a list.
##'
##' @param fun UI function
##' @return nlmixr UI function
##' @author Matthew L. Fidler
##' @keywords internal
##' @export
nlmixrUI <- function(fun) {
  if (is(fun, "function")) {
    lhs0 <- nlmixrfindLhs(body(fun))
    dum.fun <- function() {
      return(TRUE)
    }
    env.here <- environment(dum.fun)
    env <- new.env(parent = .GlobalEnv)
    assign("fun", fun, env)
    fun2 <- attr(fun, "srcref")
    if (is.null(fun2)) {
      fun2 <- deparse(fun)
      assign("fun", fun, env)
      cli::cli_alert_info("parameter labels from comments are typically ignored in non-interactive mode")
      cli::cli_alert_info("Need to run with the source intact to parse comments")
      ## message(sprintf("Cannot run this way in non-interactive mode.\nTry running:\n\nR -e 'source(\"script\", keep.source=TRUE)'\n\nFor Batch-mode type of running, you can use:\n\nR -e \"source('script.R', keep.source=TRUE, echo=TRUE)\" %s> script.Rout 2>&1", ifelse((.Platform$OS.type == "unix"), "", "1")))
      ## stop("option \"keep.source\" must be TRUE for nlmixr models.")
      ## return(eval(fun(), parent.frame(1)))
    } else {
      fun2 <- as.character(fun2, useSource = TRUE)
    }
    rg <- rex::rex("function", any_spaces, "(", anything, ")")
    w <- which(regexpr(rg, fun2) != -1)
    if (length(w) > 0) {
      w <- w[1]
      fun2[w] <- sub(rg, "", fun2[w])
    }
    fun2 <- gsub(rex::rex(boundary, "ini(", any_spaces, "{"), ".ini <- function()({", fun2)
    fun2 <- gsub(rex::rex(boundary, "model(", any_spaces, "{"), ".model <- function()({", fun2)
    if (fun2[length(fun2)] != "}") {
      fun2[length(fun2)] <- sub(rex::rex("}", end), "", fun2[length(fun2)])
      fun2[length(fun2) + 1] <- "}"
    }
    w <- which(regexpr(rex::rex(start, any_spaces, "#", anything), fun2) != -1)
    if (length(w) > 0 && all(lhs0 != "desc")) {
      w2 <- w[1]
      if (length(w) > 1) {
        for (i in 2:length(w)) {
          if (w[i] - 1 == w[i - 1]) {
            w2[i] <- w[i]
          } else {
            break
          }
        }
      }
      desc <- paste(gsub(
        rex::rex(any_spaces, end), "",
        gsub(rex::rex(start, any_spaces, any_of("#"), any_spaces), "", fun2[w2])
      ),
      collapse = " "
      )
      lhs0 <- c(lhs0, "desc")
    }
    .model <- .ini <- NULL ## Rcheck hack
    eval(parse(text = paste(fun2, collapse = "\n"), keep.source = TRUE))
  } else if (is(fun, "environment")) {
    env.here <- fun
    .model <- fun$.model
    .ini <- fun$.ini
    env <- new.env(parent = .GlobalEnv)
    lhs0 <- c("data", "desc", "ref", "imp", "est", "control", "table")
    warning("some information (like parameter labels) is lost by evaluating a nlmixr function")
    fun <- paste0(
      "function(){\n",
      paste(sapply(lhs0, function(var) {
        if (exists(var, envir = env.here)) {
          if (!inherits(get(var, envir = env.here), "function")) {
            return(sprintf("\n%s <- %s;", var, .deparse(get(var, envir = env.here))))
          } else {
            return("")
          }
        } else {
          return("")
        }
      }), collapse = ""),
      sprintf("\nini(%s)", paste(.deparse(body(.ini)), collapse = "\n")),
      sprintf("\nmodel(%s)", paste(.deparse(body(.model)), collapse = "\n")),
      "\n}"
    )
    fun <- eval(parse(text = fun))
    assign("fun", fun, env)
  }
  ini <- nlmixrBounds(.ini)
  meta <- list()
  for (var in lhs0) {
    if (!any(var == c(".ini", ".model")) &&
      exists(var, envir = env.here)) {
      meta[[var]] <- get(var, envir = env.here)
    }
  }
  if (inherits(.ini, "try-error")) {
    stop("Error parsing initial estimates.")
  }

  fun2 <- nlmixrUIModel(.model, ini, fun)
  ## if (inherits(fun2, "try-error")){
  ##     stop("Error parsing model.")
  ## }
  return(.finalizeUiModel(fun2, meta))
}
nlmixrUI.multipleEndpoint <- function(x) {
  if (length(x$nmodel$predDf$cond) > 1) {
    .info <- x$nmodel$predDf
    if (getOption("RxODE.combine.dvid", TRUE)) {
      .info <- .info[order(.info$dvid), ]
    }
    .info <- with(.info, .data.frame(
      variable = paste(var, "~", ifelse(use.utf(), "\u2026", "...")),
      cmt = paste0("cmt='", cond, "' or cmt=", cmt),
      "dvid*" = ifelse(is.na(dvid), "",
        paste0("dvid='", cond, "' or dvid=", dvid)
      ),
      check.names = FALSE
    ))
    if (!getOption("RxODE.combine.dvid", TRUE)) {
      .info <- .info[, names(.info) != "dvid*"]
    }
    if (requireNamespace("huxtable", quietly = TRUE)) {
      .hux <- huxtable::hux(.info) %>%
        huxtable::add_colnames() %>%
        huxtable::set_bold(row = 1, col = huxtable::everywhere, value = TRUE) %>%
        huxtable::set_position("center") %>%
        huxtable::set_all_borders(TRUE)
      if (getOption("RxODE.combine.dvid", TRUE)) {
        .hux <- .hux %>%
          huxtable::add_footnote("* If dvids are outside this range, all dvids are re-numered sequentially, ie 1,7, 10 becomes 1,2,3 etc")
      }
    } else {
      .hux <- .info
    }
    return(.hux)
  } else {
    return(NULL)
  }
}
##' Print UI function
##'
##' @param x  UI function
##' @param ... other arguments
##' @author Matthew L. Fidler
##' @return original object (invisibly)
##' @export
print.nlmixrUI <- function(x, ...) {
  cat(cli::rule(x$model.desc, line = "bar2"), "\n")
  cat(cli::rule(crayon::bold("Initialization:")), "\n")
  print(x$ini)
  if (length(x$all.covs) > 0) {
    cat("\n Covariates or Uninitialized Parameters ($all.covs)\n")
    print(x$all.covs)
  }
  if (length(x$predDf$cond) > 1) {
    cat(cli::rule(paste0(crayon::bold("Multiple Endpoint Model"), " (", crayon::bold$blue("$multipleEndpoint"), "):")), "\n")

    if (requireNamespace("huxtable", quietly = TRUE)) {
      x$multipleEndpoint %>%
        huxtable::print_screen(colnames = FALSE)
    } else {
      print(x$multipleEndpoint)
    }
    cat("\n")
  }
  .mu <- x$muRefTable
  if (length(.mu) > 0) {
    cat(cli::rule(paste0(
      crayon::bold(paste0(ifelse(use.utf(), "\u03bc", "mu"), "-referencing")),
      " (", crayon::bold$blue("$muRefTable"), "):"
    )), "\n")
    if (requireNamespace("huxtable", quietly = TRUE)) {
      .mu %>%
        huxtable::print_screen(colnames = FALSE)
    } else {
      print(.mu)
    }
    cat("\n")
  }
  cat(cli::rule(crayon::bold(sprintf("Model%s:", ifelse(inherits(x$rxode, "RxODE"), " (RxODE)", "")))), "\n")
  cat(x$fun.txt, "\n")
  cat(cli::rule(line = "bar2"), "\n")
  invisible(x)
}

## This is a list of supported distributions with the number of arguments they currently support.
dists <- list(
  "dpois" = 0,
  "dbinom" = 0:1,
  "dbern" = 0,
  "bern" = 0,
  "dbeta" = 2:3,
  ##
  ## "dnbinom"=2:3,  ## dnbinom is in R; FIXME: how does ot compare to dneg_binomial
  ## "dneg_binomial", ## not in base R (but in glnmm2)
  ##
  ## Available as external package http://ugrad.stat.ubc.ca/R/library/rmutil/html/BetaBinom.html
  ## "dbetabinomial", ## not in base R (but in glnmm2)
  "dt" = 1:2,
  "pois" = 0,
  "binom" = 0:1,
  "beta" = 2:3,
  "t" = 1:2,
  "add" = 1,
  "norm" = 1,
  "dnorm" = 1,
  "prop" = 1,
  "propT" = 1,
  "pow" = 2,
  "powT" = 2,
  "tbs" = 1,
  "boxCox" = 1,
  "tbsYj" = 1,
  "yeoJohnson" = 1,
  "logn" = 1,
  "dlogn" = 1,
  "logitNorm" = 1:3,
  "probitNorm" = 1:3,
  "lnorm" = 1,
  "dlnorm" = 1
)

distsPositive <- c("add", "norm", "dnorm", "prop", "propT", "pow", "powT", "logn", "dlogn", "lnorm", "dlnorm", "logitNorm", "probitNorm")

allVars <- function(x) {
  defined <- character()
  this.env <- environment()
  f <- function(x) {
    if (is.atomic(x)) {
      character()
    } else if (is.name(x)) {
      as.character(x)
    } else if (is.call(x) || is.pairlist(x)) {
      if (identical(x[[1]], quote(`~`)) ||
        identical(x[[1]], quote(`=`)) ||
        identical(x[[1]], quote(`<-`))) {
        if (is.call(x[[3]])) {
          ret <- unique(unlist(lapply(x[[3]][-1], f)))
        } else {
          ret <- unique(unlist(lapply(x[[3]], f)))
        }
        ret <- ret[!(ret %in% defined)]
        assign("defined", unique(c(defined, x[[2]])), this.env)
        return(ret)
      } else {
        children <- lapply(x[-1], f)
        unique(unlist(children))
      }
    } else {
      stop("Don't know how to handle type ", typeof(x),
        call. = FALSE
      )
    }
  }
  f(x)
}

allNames <- function(x) {
  if (is.atomic(x)) {
    character()
  } else if (is.name(x)) {
    as.character(x)
  } else if (is.call(x) || is.pairlist(x)) {
    children <- lapply(x[-1], allNames)
    unique(unlist(children))
  } else {
    stop("Don't know how to handle type ", typeof(x),
      call. = FALSE
    )
  }
}

allCalls <- function(x) {
  if (is.atomic(x) || is.name(x)) {
    character()
  } else if (is.call(x)) {
    fname <- as.character(x[[1]])
    children <- lapply(x[-1], allCalls)
    unique(c(fname, unlist(children)))
  } else if (is.pairlist(x)) {
    unique(unlist(lapply(x[-1], allCalls), use.names = FALSE))
  } else {
    stop("Don't know how to handle type ", typeof(x), call. = FALSE)
  }
}

##' Find the assignments in R expression
##'
##' @param x R expression
##' @return list of assigned parameters
##' @author Hadley Wickham and Matthew L. Fidler
##' @keywords internal
##' @export
nlmixrfindLhs <- function(x) {
  ## Modified from http://adv-r.had.co.nz/Expressions.html find_assign4
  if (is.atomic(x) || is.name(x)) {
    character()
  } else if (is.call(x)) {
    if ((identical(x[[1]], quote(`<-`)) ||
      identical(x[[1]], quote(`=`))) &&
      is.name(x[[2]])) {
      lhs <- as.character(x[[2]])
    } else {
      lhs <- character()
    }
    unique(c(lhs, unlist(lapply(x, nlmixrfindLhs))))
  } else if (is.pairlist(x)) {
    unique(unlist(lapply(x, nlmixrfindLhs)))
  } else {
    stop("Don't know how to handle type ", typeof(x),
      call. = FALSE
    )
  }
}

.findDist <- function(x) {
  ## Modified from http://adv-r.had.co.nz/Expressions.html find_assign4
  if (is.atomic(x) || is.name(x)) {
    character()
  } else if (is.call(x)) {
    if (identical(x[[1]], quote(`~`))) {
      if (is.name(x[[2]])) {
        lhs <- as.character(x[[2]])
      } else {
        if (is.call(x[[2]]) && identical(x[[2]][[1]], quote(`linCmt`))) {
          lhs <- "linCmt()"
        } else {
          lhs <- character()
        }
      }
    } else {
      lhs <- character()
    }
    unique(c(lhs, unlist(lapply(x, .findDist))))
  } else if (is.pairlist(x)) {
    unique(unlist(lapply(x, .findDist)))
  } else {
    stop("Don't know how to handle type ", typeof(x),
      call. = FALSE
    )
  }
}

## cauchy = t w/df=1; Support?
## dgeom is a special negative binomial; Support?
## fixme
unsupported.dists <- c(
  "dchisq", "chisq", "dexp", "df", "f", "dgeom", "geom",
  "dhyper", "hyper", "dunif", "unif",
  "dweibull", "weibull",
  ## for testing...
  "nlmixrDist"
)
add.dists <- c("add", "prop", "propT", "norm", "pow", "powT", "dnorm", "logn", "lnorm", "dlnorm", "tbs", "tbsYj", "boxCox", "yeoJohnson", "logitNorm", "probitNorm")
nlmixrUIModel <- function(fun, ini = NULL, bigmodel = NULL) {
  .md5 <- digest::digest(list(
    .deparse1(body(fun)), .deparse(ini), .deparse(bigmodel),
    ## Should give different models for different nlmixr versions
    sessionInfo()$otherPkgs$nlmixr$Version
  ))
  .wd <- RxODE::rxTempDir()
  if (.wd == "") {
    warning("rxTempDir did not work")
    .wd <- tempfile()
    dir.create(.wd, recursive = TRUE)
  } else {
    ## This makes all the parsing files, cpp and so in their own
    ## directory.  No collisions.
    suppressWarnings({
      dir.create(.wd, recursive = TRUE)
    })
  }
  .uiFile <- file.path(.wd, paste0("ui-", .md5))
  .uiLock <- paste0(.uiFile, ".lock")
  .uiBad <- paste0(.uiFile, ".bad")
  .uiFile <- paste0(.uiFile, ".uid")
  if (file.exists(.uiLock)) {
    message(sprintf("Waiting for UI to parse on another thread (%s)", .uiLock), appendLF = FALSE)
    while (file.exists(.uiLock)) {
      Sys.sleep(0.5)
      message(".", appendLF = FALSE)
    }
    message("")
    if (file.exists(.uiBad)) {
      load(.uiBad)
      for (.w in .lst$ws) {
        warning(.w)
      }
      stop(paste(.lst$em, "(another thread)"))
    }
    load(file = .uiFile)
    return(ret)
  } else if (file.exists(.uiBad)) {
    load(.uiBad)
    for (.w in .lst$ws) {
      warning(.w)
    }
    stop(sprintf("%s\nbad parsed model (cached %s)", .lst$em, .uiBad))
  } else if (file.exists(.uiFile)) {
    load(file = .uiFile)
    return(ret)
  } else {
    sink(.uiLock)
    cat("")
    sink()
    on.exit(
      {
        unlink(.uiLock)
      },
      add = TRUE
    )
    sink(.uiBad)
    cat("")
    sink()
    .thisEnv <- environment()
    .thisEnv$ws <- character(0)
    ret <- tryCatch(suppressWarnings(withCallingHandlers(.nlmixrUIModel(fun, ini, bigmodel),
      warning = function(w) {
        assign("ws", unique(c(w$message, .thisEnv$ws)), .thisEnv)
      }
    )), error = function(e) {
      assign("em", e$message, .thisEnv)
    })
    .lst <- list()
    if (exists("ws", envir = .thisEnv)) {
      .lst$ws <- .thisEnv$ws
    }
    if (exists("em", envir = .thisEnv)) {
      .lst$em <- .thisEnv$em
      save(.lst, file = .uiBad)
      ## stop(.lst$em);
      ## Let it error out on its own for the backtrace...
      ret <- .nlmixrUIModel(fun, ini, bigmodel)
    }
    save(ret, .lst, file = .uiFile)
    unlink(.uiBad)
    return(ret)
  }
}
.nlmixrUIModel <- function(fun, ini = NULL, bigmodel = NULL) {
  ## Parses the UI function to extract predictions and errors, and the other model specification.
  .doAssign <- TRUE
  .assign <- function(...) {
    if (.doAssign) assign(...)
  }
  .fun000 <- fun
  rxode <- FALSE
  all.names <- allNames(body(fun))
  .diff <- setdiff(paste(ini$name), all.names)
  if (length(.diff) > 0) {
    .diff <- .diff[regexpr("[(]", .diff) == -1]
    if (length(.diff) > 0) {
      stop(sprintf("The following parameter(s) were in the ini block but not in the model block: %s", paste(.diff, collapse = ", ")))
    }
  }
  if (any(regexpr(rex::rex(start, or("rx_", "nlmixr_")), paste(all.names)) != -1)) {
    stop("Parameters/States/Variables cannot start with `rx_` or `nlmixr_`")
  }
  all.vars <- allVars(body(fun))
  all.funs <- allCalls(body(fun))
  all.lhs <- nlmixrfindLhs(body(fun))
  errs.specified <- c()
  add.prop.errs <- .data.frame(y = character(), add = logical(), prop = logical())
  bounds <- ini
  theta.names <- c()
  theta.ord <- c()
  eta.names <- c()
  .mu.ref <- list()
  .oneTheta <- c()
  cov.ref <- list()
  cov.theta <- c()
  log.theta <- c()
  logit.theta <- c()
  logit.theta.low <- c()
  logit.theta.hi <- c()
  probit.theta <- c()
  probit.theta.low <- c()
  probit.theta.hi <- c()
  log.eta <- c()
  logit.eta <- c()
  logit.eta.low <- c()
  logit.eta.hi <- c()
  probit.eta <- c()
  probit.eta.low <- c()
  probit.eta.hi <- c()
  this.env <- environment()
  if (!is.null(ini)) {
    unnamed.thetas <- ini$ntheta[(!is.na(ini$ntheta) & is.na(ini$name))]
    if (length(unnamed.thetas) > 0) {
      stop(sprintf("The following THETAs are unnamed: %s", paste(sprintf("THETA[%d]", unnamed.thetas), collapse = ", ")))
    }
    unnamed.etas <- ini$neta1[!is.na(ini$neta1) & (ini$neta1 == ini$neta2) & is.na(ini$name)]
    if (length(unnamed.etas) > 0) {
      stop(sprintf("The following ETAs are unnamed: %s", paste(sprintf("ETA[%d]", unnamed.etas), collapse = ", ")))
    }
    theta.names <- ini$theta.names
    eta.names <- ini$eta.names
  }
  errn <- 0

  any.theta.names <- function(what, th.names) {
    return(any(what[1] == th.names))
  }

  find.theta <- function(x) {
    if (length(x) == 0) {
    } else if (is.name(x) && length(x) == 1 && any.theta.names(as.character(x), theta.names)) {
      return(as.character(x))
    } else if (is.call(x) || is.pairlist(x)) {
      if (length(x) == 0) {
      } else if (length(x) == 1 && any.theta.names(as.character(x), theta.names)) {
        return(as.character(x))
      } else if (identical(x[[1]], quote(`+`))) {
        th <- c()
        if (length(x) >= 3) {
          if (any.theta.names(as.character(x[[3]]), theta.names)) {
            th <- as.character(x[[3]])
          }
        }
        if (length(x) >= 2) {
          if (length(x[[2]]) > 1) {
            return(c(th, find.theta(x[[2]])))
          } else {
            if (any.theta.names(as.character(x[[2]]), theta.names)) {
              th <- c(th, as.character(x[[2]]))
            }
            return(th)
          }
        }
      }
    }
  }
  .regPar <- rex::rex("nlmixr_", capture(anything), "_par")
  .doDist <- function(distName, distArgs, curCond = NULL) {
    if (any(regexpr(.regPar, distArgs) != -1)) {
      .tmp <- distArgs[regexpr(.regPar, distArgs) != -1]
      distArgs <- gsub(.regPar, "\\1", distArgs)
    }
    if (!any(names(dists) == distName)) {
      stop(sprintf("the %s distribution is currently unsupported", distName))
    }
    .nargs <- dists[[distName]]
    if (length(.nargs) == 1) {
      if (length(distArgs) != .nargs) {
        stop(sprintf("the %s distribution requires %s argument%s", distName, .nargs, ifelse(.nargs == 1, "", "s")))
      }
    } else {
      .minNargs <- min(.nargs)
      .maxNargs <- max(.nargs)
      if (length(distArgs) < .minNargs | length(distArgs) > .maxNargs) {
        stop(sprintf("the %s distribution requires %s-%s arguments", distName, min(.nargs), max(.nargs)))
      }
    }
    .trLow <- -Inf
    .trHi <- Inf
    if (distName == "logitNorm") {
      if (length(distArgs) >= 2) {
        .trLow <- suppressWarnings(as.numeric(distArgs[2]))
        if (is.na(.trLow)) {
          stop("logitNorm lower bound must be numeric")
        }
        .trHi <- suppressWarnings(as.numeric(distArgs[3]))
        if (length(distArgs) == 3) {
          if (is.na(suppressWarnings(.trHi))) {
            stop("logitNorm upper bound must be numeric")
          }
        }
      }
      distArgs <- distArgs[1]
    }
    if (distName == "probitNorm") {
      if (length(distArgs) >= 2) {
        .trLow <- suppressWarnings(as.numeric(distArgs[2]))
        if (is.na(.trLow)) {
          stop("probitNorm lower bound must be numeric")
        }
        .trHi <- suppressWarnings(as.numeric(distArgs[3]))
        if (length(distArgs) == 3) {
          if (is.na(.trHi)) {
            stop("probitNorm upper bound must be numeric")
          }
        }
      }
      distArgs <- distArgs[1]
    }
    .est0 <- -Inf
    if (length(distArgs) >= 1) {
      if (distArgs[1] == "NA") {
        if (distName %in% c("dlnorm", "lnorm", "logitNorm", "probitNorm", "dlogn")) {
          distArgs <- c()
          .est0 <- NA_real_
        }
      }
    }
    if (length(distArgs) == 0L) {
      .tmp <- .as.data.frame(bounds)
      .tmp1 <- .tmp[1, ]
      .tmp1[1, ] <- NA
      .tmp1[, "err"] <- distName
      .tmp1[, "est"] <- NA
      if (!is.null(curCond)) {
        .tmp1[, "condition"] <- .deparse(curCond)
      } else {
        .tmp1[, "condition"] <- ""
      }
      .tmp1[, "trHi"] <- .trHi
      .tmp1[, "trLow"] <- .trLow
      .tmp <- rbind(.tmp, .tmp1)
      class(.tmp) <- c("nlmixrBounds", "data.frame")
      .assign("bounds", .tmp, this.env)
      return(errn)
    }
    for (.i in seq_along(distArgs)) {
      .tmp <- suppressWarnings(as.numeric(distArgs[.i]))
      errn <- errn + 1
      if (!is.na(.tmp)) {
        ## FIXME: allow numeric estimates...?
        stop("distribution parameters cannot be numeric, but need to be estimated")
      }
      .w <- which(bounds$name == distArgs[.i])
      if (length(.w) == 0) {
        stop("residual distribution parameter(s) estimates were not found in ini block")
      }
      .tmp <- .as.data.frame(bounds)
      .tmp$err[.w] <- ifelse(.i == 1, distName, paste0(distName, .i))
      if (any(distName == distsPositive)) {
        if (any(is.na(.tmp$lower[.w])) || any(is.infinite(.tmp$lower[.w]))) {
          .tmp$lower[.w] <- 0
        }
        if ((.tmp$lower[.w] < 0) || .tmp$est[.w] < 0 || .tmp$upper[.w] < 0) {
          stop(sprintf("The distribution '%s' must have positive parameter estimates", distName))
        }
      }
      if (!is.null(curCond)) {
        .tmp$condition[.w] <- sub(rex::rex(or("cmt", "CMT"), any_spaces, "==", any_spaces), "", .deparse(curCond))
      } else {
        .tmp$condition[.w] <- ""
      }
      .tmp$trHi[.w] <- .trHi
      .tmp$trLow[.w] <- .trLow
      class(.tmp) <- c("nlmixrBounds", "data.frame")
      .assign("bounds", .tmp, this.env)
    }
    return(errn)
  }
  .predDf <- NULL
  .doDist1 <- function(err1, err1.v, err1.args, x2, x3, curCond = NULL) {
    .bCond <- is.null(curCond)
    if (is.null(curCond)) curCond <- x2
    if (!is.na(suppressWarnings(as.numeric(err1.v)))) {
      stop("Distribution parameters cannot be numeric, but need to be estimated.")
    }
    .assign("errs.specified", unique(errs.specified, err1), this.env)
    if (any(do.pred == c(2, 4, 5))) {
      return(quote(nlmixrIgnore()))
    }
    else if (any(do.pred == c(1, 4, 5))) {
      .assign(".predDf", rbind(
        .predDf,
        .data.frame(
          cond = ifelse(.bCond, "", sub(rex::rex(or("cmt", "CMT"), any_spaces, "==", any_spaces), "", .deparse(curCond))),
          var = .deparse(x2)
        )
      ), this.env)
      return(bquote(if (CMT == .(curCond)) {
        nlmixr_pred <- .(x2)
      }))
    } else if (do.pred == 3) {
      .doDist(err1, err1.args, curCond)
      tmp <- bounds
      if ((any(paste(tmp$err) == "add") || any(paste(tmp$err) == "norm") || any(paste(tmp$err) == "dnorm") ||
        any(paste(tmp$err) == "lnorm") || any(paste(tmp$err) == "dlnorm") || any(paste(tmp$err) == "logn") ||
        any(paste(tmp$err) == "dlogn")) &&
        (any(paste(tmp$err) == "prop") || any(paste(tmp$err) == "propT"))) {
        .assign("errn", errn + 1, this.env)
        .assign("add.prop.errs", rbind(
          add.prop.errs,
          .data.frame(y = sprintf("Y%02d", errn), add = TRUE, prop = TRUE)
        ), this.env)
      } else if ((any(paste(tmp$err) == "prop") || any(paste(tmp$err) == "propT"))) {
        .assign("errn", errn + 1, this.env)
        .assign("add.prop.errs", rbind(
          add.prop.errs,
          .data.frame(y = sprintf("Y%02d", errn), add = FALSE, prop = TRUE)
        ), this.env)
      }
      return(bquote(return(.(sprintf("Y%02d", errn)))))
    } else {
      return(bquote(if (CMT == .(curCond)) {
        return(.(x3))
      }))
    }
  }
  .doDist2 <- function(err1, err1.v, err1.args, err2, err2.v, err2.args, x2, x3, curCond = NULL) {
    .bCond <- is.null(curCond)
    if (is.null(curCond)) curCond <- x2
    if (!is.na(suppressWarnings(as.numeric(err1.v)))) {
      stop("Distribution parameters cannot be numeric, but need to be estimated.")
    }
    if (!is.na(suppressWarnings(as.numeric(err2.v)))) {
      stop("Distribution parameters cannot be numeric, but need to be estimated.")
    }
    if (any(err1 == add.dists) &&
      any(err2 == add.dists)) {
      tmp <- paste(sort(c(err1, err2)), collapse = "+")
      .assign("errs.specified", unique(errs.specified, tmp), this.env)
      if (any(do.pred == c(2, 4, 5))) {
        return(quote(nlmixrIgnore()))
      }
      else if (any(do.pred == c(1, 4, 5))) {
        .assign(".predDf", rbind(
          .predDf,
          .data.frame(
            cond = ifelse(.bCond, "", sub(rex::rex(or("cmt", "CMT"), any_spaces, "==", any_spaces), "", .deparse(curCond))),
            var = .deparse(x2)
          )
        ), this.env)
        return(bquote(if (CMT == .(curCond)) {
          nlmixr_pred <- .(x2)
        }))
      } else if (do.pred == 3) {
        .doDist(err1, err1.args, curCond)
        .doDist(err2, err2.args, curCond)
        tmp <- bounds
        if ((any(paste(tmp$err) == "add") ||
          any(paste(tmp$err) == "norm") ||
          any(paste(tmp$err) == "dnorm") ||
          any(paste(tmp$err) == "lnorm") ||
          any(paste(tmp$err) == "dlnorm") ||
          any(paste(tmp$err) == "logn") ||
          any(paste(tmp$err) == "dlogn")
        ) && (any(paste(tmp$err) == "prop") || any(paste(tmp$err) == "propT"))) {
          .assign("errn", errn + 1, this.env)
          .assign("add.prop.errs", rbind(
            add.prop.errs,
            .data.frame(y = sprintf("Y%02d", errn), add = TRUE, prop = TRUE)
          ), this.env)
        }
        return(bquote(return(.(sprintf("Y%02d", errn)))))
      } else {
        return(bquote(if (CMT == .(curCond)) {
          return(.(x3))
        }))
      }
    } else {
      stop(sprintf(
        "the %s and %s distributions cannot be combined\ncurrently can combine: %s",
        as.character(x3[[2]][[1]]), as.character(x3[[3]][[1]]),
        paste(add.dists, collapse = ", ")
      ))
    }
  }
  .doDist3 <- function(err1, err1.v, err1.args, err2, err2.v, err2.args, err3, err3.v, err3.args, x2, x3, curCond = NULL) {
    .bCond <- is.null(curCond)
    if (is.null(curCond)) curCond <- x2
    if (!is.na(suppressWarnings(as.numeric(err1.v)))) {
      stop("Distribution parameters cannot be numeric, but need to be estimated.")
    }
    if (!is.na(suppressWarnings(as.numeric(err2.v)))) {
      stop("Distribution parameters cannot be numeric, but need to be estimated.")
    }
    if (!is.na(suppressWarnings(as.numeric(err3.v)))) {
      stop("Distribution parameters cannot be numeric, but need to be estimated.")
    }
    if (any(err1 == add.dists) &&
      any(err2 == add.dists) &&
      any(err3 == add.dists)) {
      tmp <- paste(sort(c(err1, err2, err3)), collapse = "+")
      .assign("errs.specified", unique(errs.specified, tmp), this.env)
      if (any(do.pred == c(2, 4, 5))) {
        return(quote(nlmixrIgnore()))
      }
      else if (any(do.pred == c(1, 4, 5))) {
        .assign(".predDf", rbind(
          .predDf,
          .data.frame(
            cond = ifelse(.bCond, "", sub(rex::rex(or("cmt", "CMT"), any_spaces, "==", any_spaces), "", .deparse(curCond))),
            var = .deparse(x2)
          )
        ), this.env)
        return(bquote(if (CMT == .(curCond)) {
          nlmixr_pred <- .(x2)
        }))
      } else if (do.pred == 3) {
        .doDist(err1, err1.args, curCond)
        .doDist(err2, err2.args, curCond)
        .doDist(err3, err3.args, curCond)
        tmp <- bounds
        if ((any(paste(tmp$err) == "add") ||
          any(paste(tmp$err) == "dnorm") ||
          any(paste(tmp$err) == "norm") ||
          any(paste(tmp$err) == "lnorm") ||
          any(paste(tmp$err) == "dlnorm") ||
          any(paste(tmp$err) == "logn") ||
          any(paste(tmp$err) == "dlogn")
        ) && (any(paste(tmp$err) == "prop") || any(paste(tmp$err) == "propT"))) {
          .assign("errn", errn + 1, this.env)
          .assign("add.prop.errs", rbind(
            add.prop.errs,
            .data.frame(y = sprintf("Y%02d", errn), add = TRUE, prop = TRUE)
          ), this.env)
        }
        return(bquote(return(.(sprintf("Y%02d", errn)))))
      } else {
        return(bquote(if (CMT == .(curCond)) {
          return(.(x3))
        }))
      }
    } else {
      stop(sprintf(
        "the %s, %s and %s distributions cannot be combined\ncurrently can combine: %s",
        err1, err2, err3, paste(add.dists, collapse = ", ")
      ))
    }
  }
  .ff <- function(x) {
    if (is.name(x)) {
      .x <- as.character(x)
      .x <- sub("^nlmixr_(.*)_par$", "\\1", .x)
      if (any.theta.names(.x, theta.names)) {
        ## print(as.character(x))
        .assign("theta.ord", unique(c(theta.ord, .x)), this.env)
      }
    } else if (is.call(x)) {
      if ((identical(x[[1]], quote(`=`)) || identical(x[[1]], quote(`<-`)))) {
        .x2 <- deparse1(x[[2]])
        .x3 <- deparse1(x[[3]])
        if (paste0("nlmixr_", .x3, "_par") == .x2) {
          return(NULL)
        }
      }
      lapply(x, .ff)
    }
  }
  f <- function(x) {
    if (is.name(x)) {
      return(x)
    } else if (is.call(x)) {
      .isTilde <- identical(x[[1]], quote(`~`))
      if (.isTilde &&
        as.character(x[[3]][[1]]) == "|" &&
        any(as.character(x[[3]][[2]][[1]]) == c(names(dists), unsupported.dists))) {
        ch.dist <- as.character(x[[3]][[2]])
        if (length(x[[3]][[3]]) == 1) {
          curCond <- sprintf("%s", as.character(x[[3]][[3]]))
        } else {
          curCond <- .deparse(x[[3]][[3]])
        }
        curCond <- eval(parse(text = sprintf("quote(%s)", curCond)))
        if (length(ch.dist) == 1) {
          err1 <- as.character(x[[3]][[2]][[1]])
          err1.v <- err1
          err1.args <- character(0)
          .doDist1(err1, err1.v, err1.args, x[[2]], x[[3]][[2]], curCond = curCond)
        } else {
          err1 <- as.character(x[[3]][[2]][[1]])
          err1.v <- as.character(x[[3]][[2]][[2]])
          err1.args <- as.character(x[[3]][[2]][-1])
          .doDist1(err1, err1.v, err1.args, x[[2]], x[[3]][[2]], curCond = curCond)
        }
      } else if (.isTilde &&
        as.character(x[[3]][[1]]) == "|" &&
        identical(x[[3]][[2]][[1]], quote(`+`)) &&
        length(as.character(x[[3]][[2]])) == 3 &&
        any(as.character(x[[3]][[2]][[2]][[1]]) == c(names(dists), unsupported.dists)) &&
        any(as.character(x[[3]][[2]][[3]][[1]]) == c(names(dists), unsupported.dists))) {
        err1 <- as.character(x[[3]][[2]][[2]][[1]])
        err1.v <- as.character(x[[3]][[2]][[2]][[2]])
        err1.args <- as.character(x[[3]][[2]][[2]][-1])
        err2 <- as.character(x[[3]][[2]][[3]][[1]])
        err2.v <- as.character(x[[3]][[2]][[3]][[2]])
        err2.args <- as.character(x[[3]][[2]][[3]][-1])
        if (length(x[[3]][[3]]) == 1) {
          curCond <- sprintf("%s", as.character(x[[3]][[3]]))
        } else {
          curCond <- .deparse(x[[3]][[3]])
        }
        curCond <- eval(parse(text = sprintf("quote(%s)", curCond)))
        .doDist2(err1, err1.v, err1.args, err2, err2.v, err2.args, x[[2]], x[[3]][[2]], curCond = curCond)
      } else if (.isTilde &&
        as.character(x[[3]][[1]]) == "|" &&
        identical(x[[3]][[2]][[1]], quote(`+`)) &&
        length(as.character(x[[3]][[2]])) == 3 &&
        any(as.character(x[[3]][[2]][[2]][[2]][[1]]) == c(names(dists), unsupported.dists))) {
        err1 <- as.character(x[[3]][[2]][[2]][[2]][[1]])
        err1.v <- as.character(x[[3]][[2]][[2]][[2]][[2]])
        err1.args <- as.character(x[[3]][[2]][[2]][[2]][-1])
        err2 <- as.character(x[[3]][[2]][[3]][[1]])
        err2.v <- as.character(x[[3]][[2]][[3]][[2]])
        err2.args <- as.character(x[[3]][[2]][[3]][-1])
        err3 <- as.character(x[[3]][[2]][[2]][[3]][[1]])
        err3.v <- as.character(x[[3]][[2]][[2]][[3]][[2]])
        err3.args <- as.character(x[[3]][[2]][[2]][[3]][-1])
        if (length(x[[3]][[3]]) == 1) {
          curCond <- sprintf("%s", as.character(x[[3]][[3]]))
        } else {
          curCond <- .deparse(x[[3]][[3]])
        }
        curCond <- eval(parse(text = sprintf("quote(%s)", curCond)))
        .doDist3(err1, err1.v, err1.args, err2, err2.v, err2.args, err3, err3.v, err3.args, x[[2]], x[[3]][[2]], curCond = curCond)
      } else if (.isTilde &&
        any(as.character(x[[3]][[1]]) == c(names(dists), unsupported.dists))) {
        ch.dist <- as.character(x[[3]])
        if (length(ch.dist) == 1) {
          err1 <- as.character(x[[3]][[1]])
          err1.v <- err1
          err1.args <- character(0)
          .doDist1(err1, err1.v, err1.args, x[[2]], x[[3]], curCond = NULL)
        } else {
          err1 <- as.character(x[[3]][[1]])
          err1.v <- as.character(x[[3]][[2]])
          err1.args <- as.character(x[[3]][-1])
          .doDist1(err1, err1.v, err1.args, x[[2]], x[[3]], curCond = NULL)
        }
      } else if (.isTilde && ## Arg parsing 4 should be the last....
        identical(x[[3]][[1]], quote(`+`)) &&
        identical(x[[3]][[2]][[1]], quote(`+`)) &&
        identical(x[[3]][[2]][[2]][[1]], quote(`+`)) &&
        any(as.character(x[[3]][[2]][[2]][[2]][[1]]) == c(names(dists), unsupported.dists))) {
        stop(sprintf(
          "only 3 distributions can be combined.\ncurrently can combine: %s",
          paste(add.dists, collapse = ", ")
        ))
      } else if (.isTilde &&
        identical(x[[3]][[1]], quote(`+`)) &&
        identical(x[[3]][[2]][[1]], quote(`+`)) &&
        any(as.character(x[[3]][[2]][[2]][[1]]) == c(names(dists), unsupported.dists))) {
        err1 <- as.character(x[[3]][[2]][[2]][[1]])
        err1.v <- as.character(x[[3]][[2]][[2]][[2]])
        err1.args <- as.character(x[[3]][[2]][[2]][-1])
        err2 <- as.character(x[[3]][[3]][[1]])
        err2.v <- as.character(x[[3]][[3]][[2]])
        err2.args <- as.character(x[[3]][[3]][-1])
        err3 <- as.character(x[[3]][[2]][[3]][[1]])
        err3.v <- as.character(x[[3]][[2]][[3]][[2]])
        err3.args <- as.character(x[[3]][[2]][[3]][-1])
        .doDist3(err1, err1.v, err1.args, err2, err2.v, err2.args, err3, err3.v, err3.args, x[[2]], x[[3]], curCond = NULL)
      } else if (.isTilde &&
        identical(x[[3]][[1]], quote(`+`)) &&
        length(as.character(x[[3]])) == 3 &&
        any(as.character(x[[3]][[2]][[1]]) == c(names(dists), unsupported.dists)) &&
        any(as.character(x[[3]][[3]][[1]]) == c(names(dists), unsupported.dists))) {
        err1 <- as.character(x[[3]][[2]][[1]])
        err1.v <- as.character(x[[3]][[2]][[2]])
        err1.args <- as.character(x[[3]][[2]][-1])
        err2 <- as.character(x[[3]][[3]][[1]])
        err2.v <- as.character(x[[3]][[3]][[2]])
        err2.args <- as.character(x[[3]][[3]][-1])
        .doDist2(err1, err1.v, err1.args, err2, err2.v, err2.args, x[[2]], x[[3]])
      } else if (.isTilde && (do.pred != 2)) {
        return(quote(nlmixrIgnore()))
      } else if (identical(x[[1]], quote(`<-`)) && !any(do.pred == c(2, 4, 5))) {
        .ff(x)
        return(quote(nlmixrIgnore()))
      } else if (identical(x[[1]], quote(`=`)) && !any(do.pred == c(2, 4, 5))) {
        .ff(x)
        return(quote(nlmixrIgnore()))
      } else if (identical(x[[1]], quote(`<-`)) && do.pred == 4) {
        .ff(x)
        ## SAEM requires = instead of <-
        x[[1]] <- quote(`=`)
        return(as.call(lapply(x, f)))
      } else if (identical(x[[1]], quote(`probitInv`)) && any(do.pred == c(4, 5))) {
        .low <- 0
        .hi <- 1
        if (length(x) >= 3) {
          .tmp <- try(eval(x[[3]]), silent = TRUE)
          if (inherits(.tmp, "numeric")) {
            .low <- .tmp
          } else {
            .low <- NA_real_
          }
          if (length(x) >= 4) {
            .tmp <- try(eval(x[[4]]), silent = TRUE)
            if (inherits(.tmp, "numeric")) {
              .hi <- .tmp
            } else {
              .hi <- NA_real_
            }
          }
        }
        find.probit <- function(x) {
          if (is.atomic(x) || is.name(x)) {
            if (any.theta.names(as.character(x), theta.names)) {
              .assign("probit.theta", unique(c(probit.theta, as.character(x))), this.env)
              if (!(as.character(x) %in% names(probit.theta.hi))) {
                .assign("probit.theta.hi", c(probit.theta.hi, setNames(.hi, as.character(x))), this.env)
              }
              if (!(as.character(x) %in% names(probit.theta.low))) {
                .assign("probit.theta.low", c(probit.theta.low, setNames(.low, as.character(x))), this.env)
              }
            } else if (any.theta.names(as.character(x), eta.names)) {
              .assign("probit.eta", unique(c(probit.eta, as.character(x))), this.env)
              if (!(as.character(x) %in% names(probit.eta.hi))) {
                .assign("probit.eta.hi", c(probit.eta.hi, setNames(.hi, as.character(x))), this.env)
              }
              if (!(as.character(x) %in% names(probit.eta.low))) {
                .assign("probit.eta.low", c(probit.eta.low, setNames(.low, as.character(x))), this.env)
              }
            }
            return(x)
          } else if (is.pairlist(x)) {
            return(lapply(x, find.probit))
          } else if (is.call(x)) {
            return(lapply(x, find.probit))
          } else {
            stop("Don't know how to handle type ", typeof(x),
              call. = FALSE
            )
          }
        }
        find.probit(x[[2]])
        if (length(x[[2]]) == 1 && any.theta.names(as.character(x[[2]]), theta.names)) {
          tmp <- as.character(x[[2]])
          .assign(".oneTheta", unique(c(.oneTheta, tmp)), this.env)
        }
        return(as.call(lapply(x, f)))
      } else if (identical(x[[1]], quote(`expit`)) && any(do.pred == c(4, 5))) {
        .low <- 0
        .hi <- 1
        if (length(x) >= 3) {
          .tmp <- try(eval(x[[3]]), silent = TRUE)
          if (inherits(.tmp, "numeric")) {
            .low <- .tmp
          } else {
            .low <- NA_real_
          }
          if (length(x) >= 4) {
            .tmp <- try(eval(x[[4]]), silent = TRUE)
            if (inherits(.tmp, "numeric")) {
              .hi <- .tmp
            } else {
              .hi <- NA_real_
            }
          }
        }
        find.logit <- function(x) {
          if (is.atomic(x) || is.name(x)) {
            if (any.theta.names(as.character(x), theta.names)) {
              .assign("logit.theta", unique(c(logit.theta, as.character(x))), this.env)
              if (!(as.character(x) %in% names(logit.theta.hi))) {
                .assign("logit.theta.hi", c(logit.theta.hi, setNames(.hi, as.character(x))), this.env)
              }
              if (!(as.character(x) %in% names(logit.theta.low))) {
                .assign("logit.theta.low", c(logit.theta.low, setNames(.low, as.character(x))), this.env)
              }
            } else if (any.theta.names(as.character(x), eta.names)) {
              .assign("logit.eta", unique(c(logit.eta, as.character(x))), this.env)
              if (!(as.character(x) %in% names(logit.eta.hi))) {
                .assign("logit.eta.hi", c(logit.eta.hi, setNames(.hi, as.character(x))), this.env)
              }
              if (!(as.character(x) %in% names(logit.eta.low))) {
                .assign("logit.eta.low", c(logit.eta.low, setNames(.low, as.character(x))), this.env)
              }
            }
            return(x)
          } else if (is.pairlist(x)) {
            return(lapply(x, find.logit))
          } else if (is.call(x)) {
            return(lapply(x, find.logit))
          } else {
            stop("Don't know how to handle type ", typeof(x),
              call. = FALSE
            )
          }
        }
        find.logit(x[[2]])
        if (length(x[[2]]) == 1 && any.theta.names(as.character(x[[2]]), theta.names)) {
          tmp <- as.character(x[[2]])
          .assign(".oneTheta", unique(c(.oneTheta, tmp)), this.env)
        }
        return(as.call(lapply(x, f)))
      } else if (identical(x[[1]], quote(`exp`)) && any(do.pred == c(4, 5))) {
        ## Need traverse the parsing tree to get log theta/eta
        ## parameters.
        find.log <- function(x) {
          if (is.atomic(x) || is.name(x)) {
            if (any.theta.names(as.character(x), theta.names)) {
              .assign("log.theta", unique(c(log.theta, as.character(x))), this.env)
            } else if (any.theta.names(as.character(x), eta.names)) {
              .assign("log.eta", unique(c(log.eta, as.character(x))), this.env)
            }
            return(x)
          } else if (is.pairlist(x)) {
            return(lapply(x, find.log))
          } else if (is.call(x)) {
            return(lapply(x, find.log))
          } else {
            stop("Don't know how to handle type ", typeof(x),
              call. = FALSE
            )
          }
        }
        find.log(x[[2]])
        if (length(x[[2]]) == 1 && any.theta.names(as.character(x[[2]]), theta.names)) {
          tmp <- as.character(x[[2]])
          .assign(".oneTheta", unique(c(.oneTheta, tmp)), this.env)
        }
        return(as.call(lapply(x, f)))
      } else if (identical(x[[1]], quote(`+`)) ||
        identical(x[[1]], quote(`-`))) {
        isMinus <- identical(x[[1]], quote(`-`))
        ## print(as.character(x))
        if (any(do.pred == c(4, 5))) {
          if (length(x) >= 3) {
            ## message("---")
            ## print(x[[1]])
            ## print(x[[2]]);
            ## print(x[[3]]);
            wm <- NULL
            muRef <- FALSE
            if (length(x[[3]]) == 3) {
              if (identical(x[[3]][[1]], quote(`*`))) {
                wm <- 3
                wm0 <- 2
                muRef <- TRUE
              } else if (identical(x[[3]][[1]], quote(`/`)) ||
                identical(x[[3]][[1]], quote(`^`)) ||
                identical(x[[3]][[1]], quote(`**`))) {
                wm <- 3
                wm0 <- 2
              }
            }
            if (length(x[[2]]) == 3) {
              if (identical(x[[2]][[1]], quote(`*`))) {
                wm <- 2
                wm0 <- 3
                muRef <- TRUE
              } else if (identical(x[[2]][[1]], quote(`/`)) ||
                identical(x[[2]][[1]], quote(`^`)) ||
                identical(x[[2]][[1]], quote(`**`))) {
                wm <- 3
                wm0 <- 2
              }
            }
            if (!is.null(wm)) {
              cur <- 2
              th <- 3
              w <- try(which(all.covs == as.character(x[[wm]][[2]])[1]), silent = TRUE)
              if (inherits(w, "try-error")) w <- integer(0)
              if (length(w) == 0) {
                cur <- 3
                th <- 2
                w <- try(which(all.covs == as.character(x[[wm]][[3]])[1]), silent = TRUE)
                if (inherits(w, "try-error")) w <- integer(0)
              }
              if (length(w) == 1) {
                cov <- all.covs[w]
                th <- as.character(x[[wm]][[th]])
                th0 <- find.theta(x[[wm0]])
                if (muRef && !isMinus && length(th0) == 1 && do.pred == 4) {
                  tmp <- get("cov.ref", this.env)
                  tmp[[cov]] <- c(tmp[[cov]], structure(th0, .Names = th))
                  .assign("cov.ref", tmp, this.env)
                  return(f(x[[wm0]]))
                }
                tmp <- get("cov.theta", this.env)
                .assign("cov.theta", unique(c(tmp, th)), this.env)
              }
            }
            if (any.theta.names(as.character(x[[2]]), eta.names) &&
              any.theta.names(as.character(x[[3]]), theta.names)) {
              ## Found ETA+THETA
              tmp <- .mu.ref
              tmp[[as.character(x[[2]])]] <- as.character(x[[3]])
              ## .assign("mu.ref", tmp, this.env);
              .assign(".mu.ref", tmp, this.env)
              tmp <- as.character(x[[3]])
              .assign(".oneTheta", unique(c(.oneTheta, tmp)), this.env)
              ## Collapse to THETA
              return(x[[3]])
            } else if (any.theta.names(as.character(x[[3]]), eta.names) &&
              any.theta.names(as.character(x[[2]]), theta.names)) {
              ## Found THETA+ETA
              tmp <- .mu.ref
              tmp[[as.character(x[[3]])]] <- as.character(x[[2]])
              ## .assign(".mu.ref", tmp, this.env)
              .assign(".mu.ref", tmp, this.env)
              tmp <- as.character(x[[2]])
              .assign(".oneTheta", unique(c(.oneTheta, tmp)), this.env)
              ## Collapse to THETA
              ## model$omega=diag(c(1,1,0))
              ## 0 is not estimated.
              ## inits$omega has the initial estimate
              ## a+b*f
              ## mod$ares = initial estimate of res
              ## mod$bres = initial estimate of power/prop
              ## mod$cres = initial estimate of exponent
              ## mod$lres = initial estimate of lambda transformation (either box-cox or yeo-johnson)
              return(x[[2]])
            } else if (any.theta.names(as.character(x[[3]]), eta.names) &&
              length(x[[2]]) > 1) {
              ## This allows 123 + Cl + 123 + eta.Cl + 123
              ## And collapses to 123 + Cl + 123 + 123
              ## Useful for covariates...
              eta <- as.character(x[[3]])
              th <- find.theta(x[[2]])
              if (length(th) == 1) {
                tmp <- .mu.ref
                tmp[[eta]] <- th
                ## .assign("tmp", .mu.ref, this.env)
                .assign(".mu.ref", tmp, this.env)
                tmp <- as.character(th)
                .assign(".oneTheta", unique(c(.oneTheta, tmp)), this.env)
                return(f(as.call(x[[2]])))
              }
            } else if (length(x) < 3) {
            } else if (any.theta.names(as.character(x[[3]]), theta.names) &&
              length(x[[2]]) > 1) {
              ## This allows 123 + eta.Cl + 123 + Cl + 123
              ## And collapses to 123  + 123 + Cl + 123
              ## Useful for covariates...
              theta <- as.character(x[[3]])
              .etas <- c()
              find.etas <- function(x) {
                if (is.atomic(x) || is.name(x)) {
                  return(x)
                } else if (is.pairlist(x)) {
                  return(lapply(x, find.etas))
                } else if (is.call(x)) {
                  if (identical(x[[1]], quote(`+`)) &&
                    any.theta.names(as.character(x[[3]]), eta.names)) {
                    .assign(".etas", c(.etas, as.character(x[[3]])), this.env)
                    return(x[[2]])
                  }
                  return(as.call(lapply(x, find.etas)))
                } else {
                  stop("Don't know how to handle type ", typeof(x),
                    call. = FALSE
                  )
                }
              }
              new <- find.etas(x[[2]])
              if (length(.etas) == 1) {
                tmp <- .mu.ref
                tmp[[.etas]] <- theta
                ## .assign("tmp", .mu.ref, this.env);
                .assign(".mu.ref", tmp, this.env)
                tmp <- as.character(theta)
                assign(".oneTheta", unique(c(.oneTheta, tmp)), this.env)
                x[[2]] <- new
              }
            }
          }
        }
        return(as.call(lapply(x, f)))
      } else {
        return(as.call(lapply(x, f)))
      }
    } else if (is.pairlist(x)) {
      as.pairlist(lapply(x, f))
    } else if (is.atomic(x)) {
      return(x)
    } else {
      stop("Don't know how to handle type ", typeof(x),
        call. = FALSE
      )
    }
  }
  rm.empty <- function(x) {
    ## empty if/else
    ## First remove if () followed by nlmixrIgnore()
    ignoreit <- rex::rex(any_spaces, "nlmixrIgnore()", any_spaces)
    w1 <- which(regexpr(rex::rex(start, any_spaces, or(group("if", any_spaces, "(", anything, ")"), "else"), any_spaces, end), x) != -1)
    if (length(w1) > 0) {
      w2 <- which(regexpr(ignoreit, x[w1 + 1]) != -1)
      if (length(w2) > 0) {
        x <- x[-w1[w2]]
      }
    }
    x <- x[regexpr(ignoreit, x, perl = TRUE) == -1]
    w1 <- which(regexpr(rex::rex(start, any_spaces, or("if", "else"), anything, "{", end), x) != -1)
    if (length(w1) > 0) {
      w2 <- w1 + 1
      w3 <- which(regexpr(rex::rex(start, any_spaces, "}", end), x[w2]) != -1)
      if (length(w3) > 0) {
        w1 <- w1[w3]
        w2 <- w2[w3]
        return(x[-c(w1, w2)])
      } else {
        return(x)
      }
    } else {
      return(x)
    }
  }
  new.fn <- function(x) {
    x <- rm.empty(x)
    if (do.pred == 2) {
      .assign("rxode", any(regexpr(rex::rex(start, any_spaces, "d/dt(", anything, ")", any_spaces, or("=", "<-")), x) != -1), this.env)
    }
    x <- c("function(){", .bodyDewrap(x), "}")
    x <- eval(parse(text = paste(x, collapse = "\n")))
    return(x)
  }
  .fun00 <- function(funTxt) {
    .fun0 <- function(reg0 = rex::rex(
                        capture(any_spaces),
                        capture(except_any_of("()\n; ")),
                        capture(any_spaces)
                      ),
                      reg00 = "", reg01 = "", repE = "expr") {
      if (any(regexpr(.regRx, funTxt[1:w]) != -1)) {
        .rxBegin <- funTxt[1:w]
        .re <- rex::rex(
          start, capture(any_spaces), reg0,
          capture(any_spaces, or("=", "~", "<-")),
          capture(anything)
        )
        .w <- which(regexpr(.re, .rxBegin) != -1)
        .lines <- funTxt[.w]
        ## Extract any estimated parameter only expressions and promote them.
        .w2 <- which(regexpr(rex::rex(anything, or("=", "~", "<-"), .lhsReg), .lines) != -1)
        .doIt <- TRUE
        if (length(.w2) > 0) {
          ## Remove any parameters that depend on prior values and not covariates/expressions.
          .w <- .w[-.w2]
          if (length(.w) == 0) {
            .doIt <- FALSE
          } else {
            .lines <- funTxt[.w]
          }
        }
        if (.doIt) {
          ## Now remove any variables that are duplicated.  This often happens when
          ## if (x){y=a;}else{y=b;} blocks.
          .rxBegin <- funTxt[.w]
          .dups <- funTxt[regexpr(rex::rex(or("=", "~", "<-")), funTxt) != -1]
          .dups <- gsub(rex::rex(
            start, any_spaces, reg0, any_spaces,
            or("=", "~", "<-"), anything
          ), "\\2", .dups)
          .dups <- unique(.dups[duplicated(.dups)])
          .w2 <- which(regexpr(rex::rex(
            start, any_spaces,
            reg00, or(.dups), reg01, any_spaces,
            or("=", "~", "<-"), anything
          ), .lines) != -1)
          if (length(.w2) > 0) {
            .w <- .w[-.w2]
            if (length(.w) == 0) {
              .doIt <- FALSE
            } else {
              .lines <- funTxt[.w]
            }
          }
          if (.doIt) {
            .rxBegin <- funTxt
            .rxBegin[.w] <- gsub(.re, paste0("\\1\\2\\3\\4\\5nlmixr_\\3_", repE),
              .rxBegin[.w],
              perl = TRUE
            )
            .lines <- gsub(
              rex::rex(.re, anything),
              paste0("\\1nlmixr_\\3_", repE, "\\5\\6"), .lines
            )
            ## .w <- min(which(regexpr(.regRx, .rxBegin, perl=TRUE) != -1))-1;
            ## return(c(.rxBegin[1:.w],.lines,.rxBegin[-(1:.w)]));
            return(c(.lines, .rxBegin))
          }
        }
      }
      return(funTxt)
    }
    w <- which(regexpr(reg, funTxt, perl = TRUE) != -1)
    w <- max(w)
    funTxt <- .fun0(
      reg0 = rex::rex(capture(or("f(", "F(")), capture(except_any_of("()\n; ")), capture(")")),
      reg00 = rex::rex(or("f(", "F(")), reg01 = ")", repE = "F"
    )
    w <- which(regexpr(reg, funTxt, perl = TRUE) != -1)
    w <- max(w)
    funTxt <- .fun0(
      reg0 = rex::rex(capture(or("dur(", "d(")), capture(except_any_of("()\n; ")), capture(")")),
      reg00 = rex::rex(or("dur(", "d(")), reg01 = ")", repE = "dur"
    )
    w <- which(regexpr(reg, funTxt, perl = TRUE) != -1)
    w <- max(w)
    funTxt <- .fun0(
      reg0 = rex::rex(capture(or("lag(", "alag(")), capture(except_any_of("()\n; ")), capture(")")),
      reg00 = rex::rex(or("lag(", "alag(")), reg01 = ")", repE = "lag"
    )
    w <- which(regexpr(reg, funTxt, perl = TRUE) != -1)
    w <- max(w)
    funTxt <- .fun0(
      reg0 = rex::rex(capture(or("r(", "rate(")), capture(except_any_of("()\n; ")), capture(")")),
      reg00 = rex::rex(or("r(", "rate(")), reg01 = ")", repE = "rate"
    )
    w <- which(regexpr(reg, funTxt, perl = TRUE) != -1)
    w <- max(w)

    funTxt <- .fun0(
      reg0 = rex::rex(capture(any_spaces), capture(except_any_of("()\n; ")), capture("(0)")),
      reg00 = "", reg01 = "(0)", repE = "ini"
    )
    w <- which(regexpr(reg, funTxt, perl = TRUE) != -1)
    w <- max(w)

    funTxt <- .fun0()
    w <- which(regexpr(reg, funTxt, perl = TRUE) != -1)
    w <- max(w)

    .finalFix <- function() {
      if (any(regexpr(.regRx, funTxt[1:w],
        perl = TRUE
      ) != -1)) {
        ## There are still mixed PK parameters and ODEs
        w <- which(regexpr(.regRx, funTxt, perl = TRUE) != -1)
        w <- min(w) - 1
        if (w > 0) {
          .env <- new.env(parent = emptyenv())
          .assign("extra", NULL, .env)
          .subs <- function(x) {
            if (is.atomic(x)) {
              x
            } else if (is.name(x)) {
              if (any(as.character(x) == ini$name)) {
                .assign("extra", unique(c(.env$extra, as.character(x))), .env)
                return(eval(parse(text = sprintf("quote(nlmixr_%s_par)", as.character(x)))))
              } else {
                return(x)
              }
            } else if (is.call(x)) {
              as.call(lapply(x, .subs))
            } else if (is.pairlist(x)) {
              as.pairlist(lapply(x, .subs))
            } else {
              return(x)
            }
          }
          .x <- .deparse(.subs(body(eval(parse(text = paste("function(){\n", paste(funTxt[-(1:w)], collapse = "\n"), "\n}"))))))
          .x <- .x[-1]
          .x <- .x[-length(.x)]
          return(c(funTxt[1:w], paste0("nlmixr_", .env$extra, "_par <- ", .env$extra), .x))
        }
        return(funTxt)
      }
      return(funTxt)
    }
    funTxt <- .finalFix()
    return(funTxt)
  }
  .tmp <- .deparse1(body(fun))
  ## Assume ~ is boundaries
  reg <- rex::rex(or("=", "<-"), anything, boundary, or(ini$name), boundary)
  .regRx <- rex::rex(start, any_spaces, or(
    "d/dt(", "f(", "F(", "dur(", "d(",
    "lag(", "alag(", "r(", "rate(",
    group(anything, any_spaces, "(0)", any_spaces, or("=", "<-"))
  ))
  .lhs <- nlmixrfindLhs(body(
    eval(parse(text = paste(
      "function(){",
      paste(.tmp, collapse = "\n"),
      "}"
    )))
  ))
  .lhsReg <- rex::rex(boundary, or(.lhs), boundary)
  fun <- eval(parse(text = paste(c("function()({", .fun00(.tmp), "})"), collapse = "\n")))
  fun.all <- eval(parse(text = paste(c("function()({", .tmp, "})"), collapse = "\n")))
  all.covs <- character()
  do.pred <- 1
  pred.txt <- .deparse(f(body(fun)))
  .doAssign <- FALSE
  pred.txt.all <- .deparse(f(body(fun)))
  .doAssign <- TRUE
  .reg <- rex::rex(or(
    group(any_spaces, "(", any_spaces, "{", any_spaces),
    group(any_spaces, "}", any_spaces, ")", any_spaces),
    group(any_spaces, "nlmixrIgnore()", any_spaces),
    group(any_spaces, "(", any_spaces, "{", any_spaces, "nlmixrIgnore()", any_spaces),
    group("nlmixrIgnore()", any_spaces, "}", any_spaces, ")")
  ))
  .pred <- sapply(pred.txt, function(x) {
    regexpr(.reg, x) != -1
  })
  if (all(.pred)) {
    stop("There must be at least one prediction in the model({}) block.  Use `~` for predictions")
  }
  pred <- new.fn(pred.txt)
  do.pred <- 0
  err <- new.fn(.deparse(f(body(fun))))
  .doAssign <- FALSE
  pred.all <- new.fn(pred.txt.all)
  err.all <- new.fn(.deparse(f(body(fun.all))))
  .doAssign <- TRUE
  do.pred <- 2
  rest.txt <- .deparse(f(body(fun)))
  .doAssign <- FALSE
  rest.txt.all <- .deparse(f(body(fun.all)))
  rest0 <- new.fn(rest.txt.all)
  .doAssign <- TRUE
  rest <- new.fn(rest.txt)
  rest.funs <- allCalls(body(rest))
  rest.vars <- allVars(body(rest))
  all.covs <- setdiff(rest.vars, paste0(bounds$name))
  do.pred <- 3
  grp.fn <- new.fn(.deparse(f(body(fun))))
  do.pred <- 4
  saem.pars <- try(.deparse(f(body(fun))), silent = TRUE)
  .doAssign <- FALSE
  saem.pars.all <- try(.deparse(f(body(fun.all))), silent = TRUE)
  .doAssign <- TRUE
  nlme.mu.fun2 <- NULL
  saem.fun0 <- NULL
  if (inherits(saem.pars, "try-error")) {
    saem.pars <- NULL
  } else {
    .doAssign <- FALSE
    saem.fun0 <- new.fn(saem.pars.all)
    .doAssign <- TRUE
  }
  do.pred <- 5
  nlme.mu.fun <- try(.deparse(f(body(fun))), silent = TRUE)
  if (inherits(nlme.mu.fun, "try-error")) {
    nlme.mu.fun <- NULL
  }
  .pred <- FALSE
  if (!rxode) {
    rxode <- TRUE
    .pred <- TRUE
  }
  .linCmt <- FALSE
  if (any(regexpr(rex::rex("limCmt("), .deparse(body(fun))) != -1)) {
    stop("You used `limCmt`, did you mean `linCmt`?")
  }
  if (any(regexpr(rex::rex("linCmt("), .deparse(body(fun))) != -1)) {
    .linCmt <- TRUE
    .hasLinCmt <- any(regexpr(rex::rex("linCmt("), .deparse(body(rest))) == -1)
    rx.txt <- .deparse1(body(rest))
    .regLin <- rex::rex(
      start, any_spaces,
      capture(
        or(
          group(one_of("Kk"), some_of("AaEe0123456789")),
          group(
            one_of("V", "v"),
            any_of("c", "C", "P", "p", "T", "t", "S", "s", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9")
          ),
          group(one_of("Qq"), any_of("0":"9", "p", "c")),
          group(one_of("Cc"), one_of("Ll"), any_of("Dd2")),
          "aob", "AOB", "alpha", "ALPHA", "Alpha",
          "beta", "BETA", "Beta"
        )
      ), any_spaces, or("=", "<-"),
      capture(anything)
    )
    .pred <- gsub(
      .regLin,
      "nlmixr_lincmt_\\1 <- \\2", rx.txt
    )
    .pred <- .pred[regexpr(rex::rex(
      start, any_spaces, or("dur", "d", "rate", "r", "lag", "alag", "f", "F"),
      any_spaces, "(",
      any_spaces, or("depot", "central"), any_spaces,
      ")", any_spaces, or("=", "<-"), any_spaces, capture(anything)
    ), .pred) == -1]
    .vars <- gsub(.regLin, "\\1", rx.txt[regexpr(.regLin, rx.txt) != 0])
    .rx <- gsub(
      rex::rex(
        start, any_spaces, capture(or(.vars)), any_spaces, or("=", "<-"),
        capture(anything)
      ),
      "\\1 <- nlmixr_lincmt_\\1", rx.txt
    )
    .pred <- .pred[regexpr(rex::rex("linCmt("), .pred) == -1]
    if (any(regexpr(rex::rex("linCmt("), .rx) != -1)) .hasLinCmt <- FALSE
    rest <- eval(parse(text = paste(c(
      "function(){",
      .pred, .rx,
      ifelse(.hasLinCmt, "nlmixr_lincmt_pred <- linCmt()\n}", "}")
    ),
    collapse = "\n"
    )))
    if (!is.null(saem.pars)) {
      saem.pars <- gsub(
        .regLin,
        "nlmixr_lincmt_\\1 = \\2", saem.pars
      )
    }
    if (!is.null(nlme.mu.fun)) {
      nlme.mu.fun <- gsub(
        .regLin,
        "nlmixr_lincmt_\\1 = \\2", nlme.mu.fun
      )
    }
    pred.txt <- gsub(
      rex::rex("linCmt(", any_spaces, ")"),
      "nlmixr_lincmt_pred", pred.txt
    )
    rest.txt <- gsub(.regLin, "nlmixr_lincmt_\\1 = \\2", rest.txt)
    .tmp <- .bodyDewrap(as.character(attr(fun, "srcref"), useSource = TRUE))
    .tmp <- .tmp[regexpr("~", .tmp) != -1]
    .tmp <- gsub(rex::rex("linCmt(", any_spaces, ")"), "nlmixr_lincmt_pred", .tmp)
    .tmp <- c("function()({", .pred, .rx, ifelse(.hasLinCmt, "nlmixr_lincmt_pred <- linCmt()", ""), "})")
    fun <- eval(parse(text = paste(.tmp, collapse = "\n"), keep.source = TRUE))
    rxode <- TRUE
    .pred <- FALSE
  }
  if (rxode) {
    rx.txt <- .deparse1(body(rest))
    w <- which(regexpr(reg, rx.txt, perl = TRUE) != -1)
    if (length(w) == 0) {
      stop("Error parsing model -- no parameters found.")
    }
    w <- which(regexpr(reg, rx.txt, perl = TRUE) != -1)
    w <- max(w)
    if (any(regexpr(.regRx, rx.txt[1:w], perl = TRUE) != -1)) {
      if (length(rx.txt) == w) {
        tmp <- RxODE::rxGetModel(paste(rx.txt, collapse = "\n"))
        tmp <- paste(tmp$params, "=", tmp$params)
        w <- length(tmp)
        rx.txt <- c(tmp, rx.txt)
      } else {
        stop("Mixed estimation types and ODEs.")
      }
    }
    rx.ode <- rx.txt[-(1:w)]
    rx.pred <- try(eval(parse(text = paste(c("function() {", rx.txt[1:w], "}"), collapse = "\n"))), silent = TRUE)
    if (inherits(rx.pred, "try-error")) {
      .w <- which(regexpr(rex::rex("}"), rx.ode) != -1)
      if (length(.w) > 0) {
        .i <- 1
        .w0 <- w + .w[.i]
        rx.pred <- rx.txt[1:.w0]
        rx.ode <- rx.txt[-(1:.w0)]
        if (regexpr(rex::rex("else"), rx.ode[1]) != -1) {
          .i <- .i + 1
          .w0 <- w + .w[.i]
          rx.pred <- rx.txt[1:.w0]
          rx.ode <- rx.txt[-(1:.w0)]
          w <- .w0
          if (regexpr(rex::rex("else"), rx.ode[1]) != -1) {
            stop("else if are not supported in nlmixr")
          }
        }
        if (any(regexpr(.regRx, rx.pred, perl = TRUE) != -1)) {
          stop("Mixed estimation types and ODEs #2.")
        }
      } else {
        stop("Mixed estimation types and ODEs #3.")
      }
    }
    ## Now separate out parameters for SAEM.
    .tmp <- saem.pars
    nlme.mu.fun2 <- saem.pars
    ## Now separate out for nlme.mu.fun
    rx.ode <- rx.txt[-(1:w)]
    rx.pred <- eval(parse(text = paste(c("function() {", rx.txt[1:w], "}"), collapse = "\n")))
    ## Now separate out parameters for SAEM.
    .no.mu.etas <- c()
    if (!is.null(saem.pars)) {
      w <- max(which(regexpr(reg, saem.pars, perl = TRUE) != -1))
      .tmp <- try(parse(text = paste(c(saem.pars[1:w], "})"), collapse = "\n")), silent = TRUE)
      if (inherits(.tmp, "try-error")) {
        .saemPars2 <- saem.pars[-(1:w)]
        .w <- which(regexpr(rex::rex("}"), .saemPars2) != -1)
        if (length(.w) > 0) {
          .i <- 1
          .w0 <- w + .w[.i]
          .saemPars2 <- saem.pars[-(1:.w0)]
          if (regexpr(rex::rex("else"), .saemPars2[1]) != -1) {
            .i <- .i + 1
            .w0 <- w + .w[.i]
            w <- .w0
          }
        }
      }
      saem.pars <- c(saem.pars[1:w], "")
      .no.mu.etas <- intersect(
        paste(ini$name[!is.na(ini$neta1)]),
        allVars(eval(parse(text = paste0("quote({", paste(saem.pars[-1], collapse = "\n"), "})"))))
      )
      if (length(.no.mu.etas) > 0) {
        saem.pars <- NULL
        nlme.mu.fun2 <- NULL
        nlme.mu.fun <- NULL
      } else {
        ## sapply(saem.pars, message)
        nlme.mu.fun2 <- saem.pars
        w <- max(which(regexpr(reg, nlme.mu.fun, perl = TRUE) != -1))
        .tmp <- try(parse(text = paste(c(nlme.mu.fun[1:w], "})"), collapse = "\n")), silent = TRUE)
        if (inherits(.tmp, "try-error")) {
          .saemPars2 <- nlme.mu.fun[-(1:w)]
          .w <- which(regexpr(rex::rex("}"), .saemPars2) != -1)
          if (length(.w) > 0) {
            .i <- 1
            .w0 <- w + .w[.i]
            .saemPars2 <- nlme.mu.fun[-(1:.w0)]
            if (regexpr(rex::rex("else"), .saemPars2[1]) != -1) {
              .i <- .i + 1
              .w0 <- w + .w[.i]
              w <- .w0
            }
          }
        }
        nlme.mu.fun <- c(nlme.mu.fun[1:w], "")
      }
    }

    rxode <- paste(rx.ode, collapse = "\n")
    if (.linCmt) {
      rxode <- RxODE::rxNorm(RxODE::rxGetLin(rxode))
    }
    rest <- rx.pred
    all.vars <- all.vars[!(all.vars %in% RxODE::rxState(rxode))]
    rest.vars <- rest.vars[!(rest.vars %in% RxODE::rxState(rxode))]
    all.covs <- setdiff(rest.vars, paste0(bounds$name))
    all.covs <- all.covs[!(all.covs %in% RxODE::rxLhs(rxode))]
    all.covs <- all.covs[!(all.covs %in% RxODE::rxState(rxode))]
    all.covs <- setdiff(all.covs, c(
      "t", "time", "podo", "M_E", "M_LOG2E", "M_LOG10E", "M_LN2", "M_LN10", "M_PI", "M_PI_2", "M_PI_4", "M_1_PI",
      "M_2_PI", "M_2_SQRTPI", "M_SQRT2", "M_SQRT1_2", "M_SQRT_3", "M_SQRT_32", "M_LOG10_2", "M_2PI", "M_SQRT_PI",
      "M_1_SQRT_2PI", "M_SQRT_2dPI", "M_LN_SQRT_PI", "M_LN_SQRT_2PI", "M_LN_SQRT_PId2", "pi",
      nlmixrfindLhs(body(rest))
    ))
    .mv <- RxODE::rxModelVars(rxode)
    .state <- c(.mv$state, .mv$stateExtra)
    .tmp <- .data.frame(cmt = seq_along(.state), cond = .state)
    .predDf$dvid <- seq_along(.predDf$var)
    .predDf <- merge(.predDf, .tmp, all.x = TRUE, by = "cond")
    if (any(is.na(.predDf$cmt))) {
      .predDfNa <- .predDf[is.na(.predDf$cmt), names(.predDf) != "cmt"]
      .predDf <- .predDf[!is.na(.predDf$cmt), ]
      .tmp <- .data.frame(cmt = seq_along(.state), var = .state)
      .predDf <- rbind(
        .predDf,
        merge(.predDfNa, .tmp, all.x = TRUE, by = "var")[names(.predDf)]
      )
    }
    .w <- which(is.na(.predDf$cmt))
    .predDf$cond <- paste(.predDf$cond)
    .cmtEndpoints <- integer(0)
    if (length(.w) > 0) {
      .cmtEndpoints <- paste(.predDf$cmt[-.w])
    }
    .predDf$cmt[.w] <- length(.state) + seq_along(.w)
    .w <- .predDf$cond == "" | .predDf$cond != .predDf$var
    .predDf$cond <- paste(.predDf$cond)
    .predDf$cond[.predDf$cond == ""] <- paste(.predDf$var[.predDf$cond == ""])
    .predDf$cmt[.predDf$cond != .predDf$var] <- NA
    .extra <- paste(.predDf$cond[.w])
    .extra <- .extra[regexpr("[()+-]", .extra) == -1]
    if (length(.extra) > 0) {
      .extra <- paste(paste0("cmt(", .extra, ");\n"), collapse = "\n")
      rxode <- paste0(rxode, ";\n", .extra)
    } else {
      .extra <- ""
    }
    if (any(is.na(.predDf$cmt))) {
      .predDfNa <- .predDf[is.na(.predDf$cmt), names(.predDf) != "cmt"]
      .predDf <- .predDf[!is.na(.predDf$cmt), ]
      .mv <- RxODE::rxModelVars(rxode)
      .state <- c(.mv$state, .mv$stateExtra)
      .tmp <- .data.frame(cmt = seq_along(.state), cond = .state)
      .predDf <- rbind(
        .predDf,
        merge(.predDfNa, .tmp, all.x = TRUE, by = "cond")[names(.predDf)]
      )
    }
    .dvid <- ""
    if (length(.predDf$cond) > 1) {
      .dvid <- paste0("dvid(", paste(.predDf[order(.predDf$dvid), "cmt"], collapse = ","), ")")
      rxode <- paste0(rxode, ";\n", .dvid, ";\n")
    }
    if (length(.predDf$cmt) > 1L) {
      if (length(.w) > 0L) {
        .nums <- suppressWarnings(as.numeric(paste(.predDf[.w, "cond"])))
        .w2 <- which(!is.na(.nums))
        if (length(.w2) > 0L) {
          .predDf[.w[.w2], "cmt"] <- .nums[.w2]
          .w <- which(is.na(.predDf$cmt))
        }
      }
    }
  } else {
    .predDf$cmt <- -1
    all.covs <- setdiff(rest.vars, paste0(bounds$name))
    nlme.mu.fun2 <- saem.pars
    rxode <- NULL
  }
  fun2 <- c("function(){", .bodyDewrap(as.character(attr(fun, "srcref"), useSource = TRUE)), "}")
  fun2 <- eval(parse(text = paste0(fun2, collapse = "\n"), keep.source = TRUE))
  fun2 <- as.character(attr(fun2, "srcref"), useSource = TRUE)
  fun3 <- c(.bodyDewrap(as.character(attr(.fun000, "srcref"), useSource = TRUE)))
  fun3 <- paste0(fun3, collapse = "\n")

  misplaced.dists <- intersect(rest.funs, c(names(dists), unsupported.dists))
  if (length(misplaced.dists) == 1) {
    if (misplaced.dists == "dt") {
      if (!any(regexpr("[^/]\\bdt[(]", .deparse(rest), perl = TRUE) != -1)) {
        misplaced.dists <- character()
      }
    }
    if (length(misplaced.dists) == 1) {
      if (misplaced.dists == "f") {
        if (!any(regexpr(rex::rex(one_of("Ff"), any_spaces, "(", except_some_of(")\n"), ")", any_spaces, or("<-", "=")), .deparse(rest), perl = TRUE) != -1)) {
          misplaced.dists <- character()
        }
      }
    }
  } else if (length(misplaced.dists) == 2) {
    .tmp <- order(misplaced.dists)
    .tmp <- misplaced.dists[.tmp]
    if (all(misplaced.dists == c("dt", "f"))) {
      if (!any(regexpr("[^/]\\bdt[(]", .deparse(rest), perl = TRUE) != -1)) {
        misplaced.dists <- character()
      }
      if (!any(regexpr(rex::rex(one_of("Ff"), any_spaces, "(", except_some_of(")\n"), ")", any_spaces, or("<-", "=")), .deparse(rest), perl = TRUE) != -1)) {
        misplaced.dists <- character()
      }
    }
  }
  if (length(misplaced.dists) > 0) {
    stop(sprintf("distributions need to be on residual model lines (like f ~ add(add.err)).\nmisplaced distribution(s): %s", paste(misplaced.dists, collapse = ", ")))
  }
  tmp <- gsub(rex::rex("linCmt(", any_spaces, ")"), "nlmixr_lincmt_pred", .deparse(pred))
  pred <- eval(parse(text = paste(tmp, collapse = "\n")))
  tmp <- tmp[regexpr(rex::rex("nlmixr_pred <- "), tmp) != -1]
  if (!is.null(saem.pars)) {
    saem.pars <- new.fn(saem.pars)
    saem.theta.trans <- rep(NA, length(theta.names))
  } else {
    saem.pars <- NULL
    saem.theta.trans <- NULL
  }
  if (!is.null(nlme.mu.fun)) {
    nlme.mu.fun <- new.fn(nlme.mu.fun)
  }
  if (!is.null(nlme.mu.fun2)) {
    nlme.mu.fun2 <- new.fn(nlme.mu.fun2)
  }

  cov.theta.pars <- gsub(rex::rex(or(all.covs), "."), "", names(unlist(cov.ref)))
  for (i in seq_along(theta.names)) {
    if (!any(theta.names[i] == cov.theta.pars)) {
      w <- which(theta.names[i] == theta.ord)
      if (length(w) == 1) {
        if (!any(theta.names[i] == cov.theta.pars)) {
          saem.theta.trans[i] <- w
        }
      }
    }
  }
  cur <- 1
  if (!all(is.na(saem.theta.trans))) {
    while (cur <= max(saem.theta.trans, na.rm = TRUE)) {
      while (!any(saem.theta.trans[!is.na(abs(saem.theta.trans))] == cur)) {
        w <- which(saem.theta.trans > cur)
        saem.theta.trans[w] <- saem.theta.trans[w] - 1
      }
      cur <- cur + 1
    }
  }
  env <- new.env(parent = emptyenv())
  env$sum.prod <- FALSE
  ## Split out inPars
  saem.all.covs <- all.covs[all.covs %in% names(cov.ref)]
  saem.inPars <- all.covs[!(all.covs %in% names(cov.ref))]
  .subs <- function(x) {
    if (is.call(x)) {
      if (identical(x[[1]], quote(`==`)) &&
        all(as.character(x[[3]]) == .what)) {
        if (x[[2]] != "CMT") {
          stop("Multiple endpoints can only be defined in terms of CMT")
        }
        x[[3]] <- eval(parse(text = sprintf("quote(%s)", .with)))
      }
      as.call(lapply(x, .subs))
    } else if (is.pairlist(x)) {
      as.pairlist(lapply(x, .subs))
    } else {
      return(x)
    }
  }
  if (length(.predDf$cond) > 1) {
    for (i in seq_along(.predDf$cond)) {
      .what <- (.predDf$cond[i])
      .with <- paste(.predDf$cmt[i])
      body(pred) <- .subs(body(pred))
      body(err) <- .subs(body(err))
    }
  } else {
    ## Strip if clauses out.
    .strip <- function(x) {
      if (is.call(x)) {
        if (identical(x[[1]], quote(`if`))) {
          return(x[[3]][[2]])
        }
        as.call(lapply(x, .strip))
      } else if (is.pairlist(x)) {
        as.pairlist(lapply(x, .strip))
      } else {
        return(x)
      }
    }
    body(pred) <- .strip(body(pred))
    body(err) <- .strip(body(err))
  }
  .predDf <- .predDf[order(.predDf$cmt), ]
  .w <- which(is.na(.predDf$cond))
  if (length(.w) > 0) .predDf$cond[.w] <- ""
  .predDf$var <- gsub(rex::rex("linCmt(", any_spaces, ")"), "nlmixr_lincmt_pred", paste(.predDf$var))
  .predSaem <- eval(parse(text = sprintf("function(){\n%s;\n}", paste(paste(.predDf$var), collapse = ";\n"))))
  .w <- which(!is.na(bounds$err))
  .tmpErr <- paste(bounds$name[.w])
  .errReg <- rex::rex("nlmixr_", or(.tmpErr), "_par", any_spaces, or("<-", "="), any_spaces, or(.tmpErr))
  .subs <- function(x) {
    if (is.atomic(x)) {
      x
    } else if (is.name(x)) {
      .w <- which(as.character(x) == paste0("nlmixr_", .tmpErr, "_par"))
      if (length(.w) == 1) {
        return(eval(parse(text = sprintf("quote(%s)", .tmpErr[.w]))))
      } else {
        return(x)
      }
    } else if (is.call(x)) {
      as.call(lapply(x, .subs))
    } else if (is.pairlist(x)) {
      as.pairlist(lapply(x, .subs))
    } else {
      return(x)
    }
  }
  ## Now fix the errors again...
  fun2 <- fun2[regexpr(.errReg, fun2) == -1]
  fun2 <- .deparse(.subs(body(eval(parse(text = paste(fun2, collapse = "\n"))))))
  fun2[1] <- paste("function(){")
  body(err) <- .subs(body(err))
  .tmp <- .deparse(body(rest))
  .tmp <- .tmp[regexpr(.errReg, .tmp) == -1]
  .tmp[1] <- "function(){"
  rest <- eval(parse(text = paste(.tmp, collapse = "\n")))
  if (!is.null(saem.pars)) {
    .tmp <- .deparse(body(saem.pars))
    .tmp <- .tmp[regexpr(.errReg, .tmp) == -1]
    .tmp[1] <- "function(){"
    saem.pars <- eval(parse(text = paste(.tmp, collapse = "\n")))
  }
  if (!is.null(nlme.mu.fun)) {
    .tmp <- .deparse(body(nlme.mu.fun))
    .tmp <- .tmp[regexpr(.errReg, .tmp) == -1]
    .tmp[1] <- "function(){"
    nlme.mu.fun <- eval(parse(text = paste(.tmp, collapse = "\n")))
  }
  if (!is.null(nlme.mu.fun2)) {
    .tmp <- .deparse(body(nlme.mu.fun2))
    .tmp <- .tmp[regexpr(.errReg, .tmp) == -1]
    .tmp[1] <- "function(){"
    nlme.mu.fun2 <- eval(parse(text = paste(.tmp, collapse = "\n")))
  }
  .saemErr <- ""
  if (length(intersect(c("t", "time"), allVars(body(rest)))) != 0) {
    .saemErr <- "Initial parameters defined based on ini({}) block variables cannot be defined in terms of time;\nTo fix:\n1. Define ini({}) variables/relationships first.\n2. Use the new variables in an time-based if/else clause later in the `model({})`"
  }
  if (length(.predDf$cond) == 1) {
    .w <- which(bounds$condition == "")
    bounds$condition[.w] <- paste(.predDf$cond)
    saem.theta.trans <- setdiff(saem.theta.trans, bounds$ntheta[which(!is.na(bounds$err) & !is.na(bounds$ntheta))])
  }
  ## Take out error terms (if present)
  ret <- list(
    ini = bounds, model = bigmodel,
    nmodel = list(
      fun = fun2, fun.txt = fun3, pred = pred, error = err, rest = rest, rxode = rxode,
      all.vars = all.vars, rest.vars = rest.vars, all.names = all.names,
      all.funs = all.funs, all.lhs = all.lhs,
      all.covs = all.covs, saem.all.covs = saem.all.covs,
      saem.inPars = saem.inPars,
      errs.specified = errs.specified, add.prop.errs = add.prop.errs,
      grp.fn = grp.fn, mu.ref = .mu.ref, cov.ref = cov.ref, cov.theta = cov.theta,
      saem.pars = saem.pars, nlme.mu.fun = nlme.mu.fun, nlme.mu.fun2 = nlme.mu.fun2,
      log.theta = log.theta,
      log.eta = log.eta, theta.ord = theta.ord, saem.theta.trans = saem.theta.trans,
      predDf = .predDf, predSaem = .predSaem, env = env, predSys = .pred,
      noMuEtas = .no.mu.etas,
      saemErr = .saemErr, cmtEndpoints = .cmtEndpoints,
      oneTheta = .oneTheta, extra = .extra,
      logit.theta = logit.theta, logit.theta.hi = logit.theta.hi, logit.theta.low = logit.theta.low,
      logit.eta = logit.eta, logit.eta.hi = logit.eta.hi, logit.eta.low = logit.eta.low,
      probit.theta = probit.theta, probit.theta.hi = probit.theta.hi, probit.theta.low = probit.theta.low,
      probit.eta = probit.eta, probit.eta.hi = probit.eta.hi, probit.eta.low = probit.eta.low,
      saem.fun1 = saem.fun0, focei.fun1 = rest0, extra = paste0(.extra, ifelse(.extra == "", "", ";\n"), .dvid, ifelse(.dvid == "", "", ";\n"))
    )
  )
  if (.linCmt) {
    ret$nmodel$lin.solved <- TRUE
  } else {
    ret$nmodel$lin.solved <- NULL
  }
  return(ret)
}

##' Create the nlme specs list for nlmixr nlme solving
##' @inheritParams nlmixrUI.nlmefun
##' @param mu.type is the mu-referencing type of model hat nlme will be using.
##' @return specs list for nlme
##' @author Matthew L. Fidler
nlmixrUI.nlme.specs <- function(object, mu.type = c("thetas", "covariates", "none")) {
  mu.type <- match.arg(mu.type)
  if (mu.type == "thetas") {
    return(list(
      fixed = object$fixed.form,
      random = object$random.mu,
      start = object$theta
    ))
  } else if (mu.type == "covariates") {
    theta <- names(object$theta)
    cov.ref <- object$cov.ref
    cov.theta <- unique(as.vector(unlist(cov.ref)))
    cov.base <- theta[!(theta %in% cov.theta)]
    cov.base <- cov.base[!(cov.base %in% unlist(lapply(names(cov.ref), function(x) {
      names(cov.ref[[x]])
    })))]
    cov.lst <- list()
    new.theta <- cov.base
    for (n in names(cov.ref)) {
      cov.base <- cov.base[!(cov.base %in% (names(cov.ref[[n]])))]
      cur <- cov.ref[[n]]
      for (i in seq_along(cur)) {
        m <- cur[i]
        cov.lst[[m]] <- c(cov.lst[[m]], n)
        new.theta <- c(new.theta, as.vector(m), names(m))
      }
    }
    e1 <- paste(paste(cov.base, collapse = "+"), "~ 1")
    fixed.form <- paste(c(e1, sapply(names(cov.lst), function(x) {
      paste(x, "~", paste(cov.lst[[x]], collapse = "+"))
    })), collapse = ", ")
    fixed.form <- eval(parse(text = sprintf("list(%s)", fixed.form)))
    if (length(cov.base) == 0) {
      fixed.form <- fixed.form[-1]
    }
    theta <- theta[new.theta]
    return(list(
      fixed = fixed.form,
      random = object$random.mu,
      start = object$theta
    ))
  } else {
    return(list(
      fixed = object$fixed.form,
      random = object$random,
      start = object$theta
    ))
  }
}
##' Create the nlme parameter transform function from the UI object.
##'
##' @param object UI object
##' @param mu Is the model mu referenced?
##' \itemize{
##'
##' \item With the "thetas" only the population parameters are
##' mu-referenced; All covariates are included in the model parameter
##' function.  The between subject variability pieces are specified in
##' the \code{random} specs parameter.
##'
##' \item With the "covariates" option, the population parameters are
##' mu referenced and covariates are removed from the model function.
##' The covariates will be specified used in the fixed effects
##' parameterization of nlme, like \code{list(lKA+lCL~1, lV~WT)}
##'
##' \item With the "none" option, the model function is given to nlme
##' without any modification.
##'
##' }
##' @return Parameter function for nlme
##' @author Matthew L. Fidler
##' @keywords internal
nlmixrUI.nlmefun <- function(object, mu.type = c("thetas", "covariates", "none")) {
  ## create nlme function
  mu.type <- match.arg(mu.type)
  if (mu.type == "thetas") {
    if (length(object$mu.ref) == 0L) {
      return(NULL)
    }
    .all <- object$ini$name[which(object$ini$neta1 == object$ini$neta2)] %in% names(object$mu.ref)
    if (!all(.all)) {
      return(NULL)
    }
    fn <- eval(parse(text = sprintf("function(%s) NULL", paste(unique(c(names(object$ini$theta), object$all.covs)), collapse = ", "))))
    body(fn) <- body(object$nlme.mu.fun)
  } else if (mu.type == "covariates") {
    if (length(object$mu.ref) == 0L) {
      return(NULL)
    }
    vars <- unique(c(unlist(object$mu.ref), unlist(object$cov.ref)))
    fn <- eval(parse(text = sprintf("function(%s) NULL", paste(vars, collapse = ", "))))
    body(fn) <- body(object$nlme.mu.fun2)
    vars2 <- allVars(body(fn))
    if (length(vars) != length(vars2)) {
      return(NULL)
    }
  }