utils::globalVariables(c("predicted", "std"))

#' Fit nonlinear mixed models with SAEM
#'
#' This function uses [saemix::saemix()] as a backend for fitting nonlinear mixed
#' effects models created from [mmkin] row objects using the Stochastic Approximation
#' Expectation Maximisation algorithm (SAEM).
#'
#' An mmkin row object is essentially a list of mkinfit objects that have been
#' obtained by fitting the same model to a list of datasets using [mkinfit].
#'
#' Starting values for the fixed effects (population mean parameters, argument
#' psi0 of [saemix::saemixModel()] are the mean values of the parameters found
#' using [mmkin].
#'
#' @importFrom utils packageVersion
#' @importFrom saemix saemix
#' @param object An [mmkin] row object containing several fits of the same
#' [mkinmod] model to different datasets
#' @param verbose Should we print information about created objects of
#' type [saemix::SaemixModel] and [saemix::SaemixData]?
#' @param transformations Per default, all parameter transformations are done
#' in mkin. If this argument is set to 'saemix', parameter transformations
#' are done in 'saemix' for the supported cases, i.e. (as of version 1.1.2)
#' SFO, FOMC, DFOP and HS without fixing `parent_0`, and SFO or DFOP with
#' one SFO metabolite.
#' @param error_model Possibility to override the error model used in the mmkin object
#' @param degparms_start Parameter values given as a named numeric vector will
#' be used to override the starting values obtained from the 'mmkin' object.
#' @param test_log_parms If TRUE, an attempt is made to use more robust starting
#' values for population parameters fitted as log parameters in mkin (like
#' rate constants) by only considering rate constants that pass the t-test
#' when calculating mean degradation parameters using [mean_degparms].
#' @param conf.level Possibility to adjust the required confidence level
#' for parameter that are tested if requested by 'test_log_parms'.
#' @param solution_type Possibility to specify the solution type in case the
#' automatic choice is not desired
#' @param no_random_effect Character vector of degradation parameters for
#' which there should be no variability over the groups. Only used
#' if the covariance model is not explicitly specified.
#' @param covariance.model Will be passed to [saemix::saemixModel()]. Per
#' default, uncorrelated random effects are specified for all degradation
#' parameters.
#' @param omega.init Will be passed to [saemix::saemixModel()]. If using
#' mkin transformations and the default covariance model with optionally
#' excluded random effects, the variances of the degradation parameters
#' are estimated using [mean_degparms], with testing of untransformed
#' log parameters for significant difference from zero. If not using
#' mkin transformations or a custom covariance model, the default
#' initialisation of [saemix::saemixModel] is used for omega.init.
#' @param covariates A data frame with covariate data for use in
#' 'covariate_models', with dataset names as row names.
#' @param covariate_models A list containing linear model formulas with one explanatory
#' variable, i.e. of the type 'parameter ~ covariate'. Covariates must be available
#' in the 'covariates' data frame.
#' @param error.init Will be passed to [saemix::saemixModel()].
#' @param quiet Should we suppress the messages saemix prints at the beginning
#' and the end of the optimisation process?
#' @param nbiter.saemix Convenience option to increase the number of
#' iterations
#' @param control Passed to [saemix::saemix].
#' @param \dots Further parameters passed to [saemix::saemixModel].
#' @return An S3 object of class 'saem.mmkin', containing the fitted
#' [saemix::SaemixObject] as a list component named 'so'. The
#' object also inherits from 'mixed.mmkin'.
#' @seealso [summary.saem.mmkin] [plot.mixed.mmkin]
#' @examples
#' \dontrun{
#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
#'  function(x) subset(x$data[c("name", "time", "value")]))
#' names(ds) <- paste("Dataset", 6:10)
#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds,
#'   state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
#' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
#'
#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
#' f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
#' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
#' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
#' anova(f_saem_sfo, f_saem_fomc, f_saem_dfop)
#' anova(f_saem_sfo, f_saem_dfop, test = TRUE)
#' illparms(f_saem_dfop)
#' f_saem_dfop_red <- update(f_saem_dfop, no_random_effect = "g_qlogis")
#' anova(f_saem_dfop, f_saem_dfop_red, test = TRUE)
#'
#' anova(f_saem_sfo, f_saem_fomc, f_saem_dfop)
#' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
#' # functions from saemix
#' library(saemix)
#' compare.saemix(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)
#' plot(f_saem_fomc$so, plot.type = "convergence")
#' plot(f_saem_fomc$so, plot.type = "individual.fit")
#' plot(f_saem_fomc$so, plot.type = "npde")
#' plot(f_saem_fomc$so, plot.type = "vpc")
#'
#' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
#' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
#' anova(f_saem_fomc, f_saem_fomc_tc, test = TRUE)
#'
#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
#'   A1 = mkinsub("SFO"))
#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
#'   A1 = mkinsub("SFO"))
#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
#'   A1 = mkinsub("SFO"))
#' # The following fit uses analytical solutions for SFO-SFO and DFOP-SFO,
#' # and compiled ODEs for FOMC that are much slower
#' f_mmkin <- mmkin(list(
#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
#'   ds, quiet = TRUE)
#' # saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds
#' # each on this system, as we use analytical solutions written for saemix.
#' # When using the analytical solutions written for mkin this took around
#' # four minutes
#' f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
#' f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
#' # We can use print, plot and summary methods to check the results
#' print(f_saem_dfop_sfo)
#' plot(f_saem_dfop_sfo)
#' summary(f_saem_dfop_sfo, data = TRUE)
#'
#' # The following takes about 6 minutes
#' f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
#'   nbiter.saemix = c(200, 80))
#'
#' #anova(
#' #  f_saem_dfop_sfo,
#' #  f_saem_dfop_sfo_deSolve))
#'
#' # If the model supports it, we can also use eigenvalue based solutions, which
#' # take a similar amount of time
#' #f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen",
#' #  control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
#' }
#' @export
saem <- function(object, ...) UseMethod("saem")

#' @rdname saem
#' @export
saem.mmkin <- function(object,
  transformations = c("mkin", "saemix"),
  error_model = "auto",
  degparms_start = numeric(),
  test_log_parms = TRUE,
  conf.level = 0.6,
  solution_type = "auto",
  covariance.model = "auto",
  omega.init = "auto",
  covariates = NULL,
  covariate_models = NULL,
  no_random_effect = NULL,
  error.init = c(1, 1),
  nbiter.saemix = c(300, 100),
  control = list(displayProgress = FALSE, print = FALSE,
    nbiter.saemix = nbiter.saemix,
    save = FALSE, save.graphs = FALSE),
  verbose = FALSE, quiet = FALSE, ...)
{
  call <- match.call()
  transformations <- match.arg(transformations)
  m_saemix <- saemix_model(object, verbose = verbose,
    error_model = error_model,
    degparms_start = degparms_start,
    test_log_parms = test_log_parms, conf.level = conf.level,
    solution_type = solution_type,
    transformations = transformations,
    covariance.model = covariance.model,
    omega.init = omega.init,
    covariates = covariates,
    covariate_models = covariate_models,
    error.init = error.init,
    no_random_effect = no_random_effect,
    ...)
  d_saemix <- saemix_data(object, covariates = covariates, verbose = verbose)

  fit_failed <- FALSE
  FIM_failed <- NULL
  fit_time <- system.time({
    utils::capture.output(f_saemix <- try(saemix(m_saemix, d_saemix, control)), split = !quiet)
    if (inherits(f_saemix, "try-error")) fit_failed <- TRUE
  })

  return_data <- nlme_data(object)

  if (!fit_failed) {
    if (any(is.na(f_saemix@results@se.fixed))) FIM_failed <- c(FIM_failed, "fixed effects")
    if (any(is.na(c(f_saemix@results@se.omega, f_saemix@results@se.respar)))) {
      FIM_failed <- c(FIM_failed, "random effects and error model parameters")
    }

    transparms_optim <- f_saemix@results@fixed.effects
    names(transparms_optim) <- f_saemix@results@name.fixed

    if (transformations == "mkin") {
      bparms_optim <- backtransform_odeparms(transparms_optim,
        object[[1]]$mkinmod,
        object[[1]]$transform_rates,
        object[[1]]$transform_fractions)
    } else {
      bparms_optim <- transparms_optim
    }

    saemix_data_ds <- f_saemix@data@data$ds
    mkin_ds_order <- as.character(unique(return_data$ds))
    saemix_ds_order <- unique(saemix_data_ds)

    psi <- saemix::psi(f_saemix)
    rownames(psi) <- saemix_ds_order
    return_data$predicted <- f_saemix@model@model(
      psi = psi[mkin_ds_order, ],
      id = as.numeric(return_data$ds),
      xidep = return_data[c("time", "name")])

    return_data <- transform(return_data,
      residual = value - predicted,
      std = sigma_twocomp(predicted,
        f_saemix@results@respar[1], f_saemix@results@respar[2]))
    return_data <- transform(return_data,
      standardized = residual / std)
  }

  result <- list(
    mkinmod = object[[1]]$mkinmod,
    mmkin = object,
    solution_type = object[[1]]$solution_type,
    transformations = transformations,
    transform_rates = object[[1]]$transform_rates,
    transform_fractions = object[[1]]$transform_fractions,
    covariates = covariates,
    covariate_models = covariate_models,
    sm = m_saemix,
    so = f_saemix,
    call = call,
    time = fit_time,
    FIM_failed = FIM_failed,
    mean_dp_start = attr(m_saemix, "mean_dp_start"),
    bparms.fixed = object[[1]]$bparms.fixed,
    data = return_data,
    err_mod = object[[1]]$err_mod,
    date.fit = date(),
    saemixversion = as.character(utils::packageVersion("saemix")),
    mkinversion = as.character(utils::packageVersion("mkin")),
    Rversion = paste(R.version$major, R.version$minor, sep=".")
  )

  if (!fit_failed) {
    result$mkin_ds_order <- mkin_ds_order
    result$saemix_ds_order <- saemix_ds_order
    result$bparms.optim <- bparms_optim
  }

  class(result) <- c("saem.mmkin", "mixed.mmkin")
  return(result)
}

#' @export
#' @rdname saem
#' @param x An saem.mmkin object to print
#' @param digits Number of digits to use for printing
print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
  cat( "Kinetic nonlinear mixed-effects model fit by SAEM" )
  cat("\nStructural model:\n")
  diffs <- x$mmkin[[1]]$mkinmod$diffs
  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs)
  writeLines(strwrap(nice_diffs, exdent = 11))
  cat("\nData:\n")
  cat(nrow(x$data), "observations of",
    length(unique(x$data$name)), "variable(s) grouped in",
    length(unique(x$data$ds)), "datasets\n")

  if (inherits(x$so, "try-error")) {
    cat("\nFit did not terminate successfully\n")
  } else {
    cat("\nLikelihood computed by importance sampling\n")
    ll <- try(logLik(x$so, type = "is"), silent = TRUE)
    if (inherits(ll, "try-error")) {
      cat("Not available\n")
    } else {
      print(data.frame(
          AIC = AIC(x$so, type = "is"),
          BIC = BIC(x$so, type = "is"),
          logLik = logLik(x$so, type = "is"),
          row.names = " "), digits = digits)
    }

    cat("\nFitted parameters:\n")
    conf.int <- parms(x, ci = TRUE)
    print(conf.int, digits = digits)
  }

  invisible(x)
}

#' @rdname saem
#' @return An [saemix::SaemixModel] object.
#' @export
saemix_model <- function(object, solution_type = "auto",
  transformations = c("mkin", "saemix"), error_model = "auto",
  degparms_start = numeric(),
  covariance.model = "auto", no_random_effect = NULL,
  omega.init = "auto",
  covariates = NULL, covariate_models = NULL,
  error.init = numeric(),
  test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...)
{
  if (nrow(object) > 1) stop("Only row objects allowed")

  mkin_model <- object[[1]]$mkinmod

  if (length(mkin_model$spec) > 1 & solution_type[1] == "analytical") {
    stop("mkin analytical solutions not supported for more thane one observed variable")
  }

  degparms_optim <-  mean_degparms(object, test_log_parms = test_log_parms)
  na_degparms <- names(which(is.na(degparms_optim)))
  if (length(na_degparms) > 0) {
    message("Did not find valid starting values for ", paste(na_degparms, collapse = ", "), "\n",
      "Now trying with test_log_parms = FALSE")
    degparms_optim <-  mean_degparms(object, test_log_parms = FALSE)
  }
  if (transformations == "saemix") {
    degparms_optim <- backtransform_odeparms(degparms_optim,
      object[[1]]$mkinmod,
      object[[1]]$transform_rates,
      object[[1]]$transform_fractions)
  }
  degparms_fixed <- object[[1]]$bparms.fixed

  # Transformations are done in the degradation function by default
  # (transformations = "mkin")
  transform.par = rep(0, length(degparms_optim))

  odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
  odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE)

  odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names)
  odeparms_fixed <- degparms_fixed[odeparms_fixed_names]

  odeini_fixed <- degparms_fixed[odeini_fixed_parm_names]
  names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names)

  model_function <- FALSE

  # Model functions with analytical solutions
  # Fixed parameters, use_of_ff = "min" and turning off sinks currently not supported here
  # In general, we need to consider exactly how the parameters in mkinfit were specified,
  # as the parameters are currently mapped by position in these solutions
  sinks <- sapply(mkin_model$spec, function(x) x$sink)
  if (length(odeparms_fixed) == 0 & mkin_model$use_of_ff == "max" & all(sinks)) {
    # Parent only
    if (length(mkin_model$spec) == 1) {
      parent_type <- mkin_model$spec[[1]]$type
      if (length(odeini_fixed) == 1 && !grepl("_bound$", names(odeini_fixed))) {
        if (transformations == "saemix") {
          stop("saemix transformations are not supported for parent fits with fixed initial parent value")
        }
        if (parent_type == "SFO") {
          stop("saemix needs at least two parameters to work on.")
        }
        if (parent_type == "FOMC") {
          model_function <- function(psi, id, xidep) {
            odeini_fixed / (xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 1])
          }
        }
        if (parent_type == "DFOP") {
          model_function <- function(psi, id, xidep) {
            g <- plogis(psi[id, 3])
            t <- xidep[, "time"]
            odeini_fixed * (g * exp(- exp(psi[id, 1]) * t) +
              (1 - g) * exp(- exp(psi[id, 2]) * t))
          }
        }
        if (parent_type == "HS") {
          model_function <- function(psi, id, xidep) {
            tb <- exp(psi[id, 3])
            t <- xidep[, "time"]
            k1 = exp(psi[id, 1])
            odeini_fixed * ifelse(t <= tb,
              exp(- k1 * t),
              exp(- k1 * tb) * exp(- exp(psi[id, 2]) * (t - tb)))
          }
        }
      } else {
        if (length(odeini_fixed) == 2) {
          stop("SFORB with fixed initial parent value is not supported")
        }
        if (parent_type == "SFO") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"])
            }
          } else {
            model_function <- function(psi, id, xidep) {
              psi[id, 1] * exp( - psi[id, 2] * xidep[, "time"])
            }
            transform.par = c(0, 1)
          }
        }
        if (parent_type == "FOMC") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2])
            }
          } else {
            model_function <- function(psi, id, xidep) {
              psi[id, 1] / (xidep[, "time"]/psi[id, 3] + 1)^psi[id, 2]
            }
            transform.par = c(0, 1, 1)
          }
        }
        if (parent_type == "DFOP") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              g <- plogis(psi[id, 4])
              t <- xidep[, "time"]
              psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) +
                (1 - g) * exp(- exp(psi[id, 3]) * t))
            }
          } else {
            model_function <- function(psi, id, xidep) {
              g <- psi[id, 4]
              t <- xidep[, "time"]
              psi[id, 1] * (g * exp(- psi[id, 2] * t) +
                (1 - g) * exp(- psi[id, 3] * t))
            }
            transform.par = c(0, 1, 1, 3)
          }
        }
        if (parent_type == "SFORB") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              k_12 <- exp(psi[id, 3])
              k_21 <- exp(psi[id, 4])
              k_1output <- exp(psi[id, 2])
              t <- xidep[, "time"]

              sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 - k_1output * k_21)
              b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
              b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp

              psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
                ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
            }
          } else {
            model_function <- function(psi, id, xidep) {
              k_12 <- psi[id, 3]
              k_21 <- psi[id, 4]
              k_1output <- psi[id, 2]
              t <- xidep[, "time"]

              sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 - k_1output * k_21)
              b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
              b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp

              psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
                ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
            }
            transform.par = c(0, 1, 1, 1)
          }
        }
        if (parent_type == "HS") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              tb <- exp(psi[id, 4])
              t <- xidep[, "time"]
              k1 <- exp(psi[id, 2])
              psi[id, 1] * ifelse(t <= tb,
                exp(- k1 * t),
                exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb)))
            }
          } else {
            model_function <- function(psi, id, xidep) {
              tb <- psi[id, 4]
              t <- xidep[, "time"]
              psi[id, 1] * ifelse(t <= tb,
                exp(- psi[id, 2] * t),
                exp(- psi[id, 2] * tb) * exp(- psi[id, 3] * (t - tb)))
            }
            transform.par = c(0, 1, 1, 1)
          }
        }
      }
    }

    # Parent with one metabolite
    # Parameter names used in the model functions are as in
    # https://nbviewer.jupyter.org/urls/jrwb.de/nb/Symbolic%20ODE%20solutions%20for%20mkin.ipynb
    types <- unname(sapply(mkin_model$spec, function(x) x$type))
    if (length(mkin_model$spec) == 2 &! "SFORB" %in% types ) {
      # Initial value for the metabolite (n20) must be fixed
      if (names(odeini_fixed) == names(mkin_model$spec)[2]) {
        n20 <- odeini_fixed
        parent_name <- names(mkin_model$spec)[1]
        if (identical(types, c("SFO", "SFO"))) {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              t <- xidep[, "time"]
              n10 <- psi[id, 1]
              k1 <- exp(psi[id, 2])
              k2 <- exp(psi[id, 3])
              f12 <- plogis(psi[id, 4])
              ifelse(xidep[, "name"] == parent_name,
                n10 * exp(- k1 * t),
                (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
                  (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
              )
            }
          } else {
            model_function <- function(psi, id, xidep) {
              t <- xidep[, "time"]
              n10 <- psi[id, 1]
              k1 <- psi[id, 2]
              k2 <- psi[id, 3]
              f12 <- psi[id, 4]
              ifelse(xidep[, "name"] == parent_name,
                n10 * exp(- k1 * t),
                (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
                  (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
              )
            }
            transform.par = c(0, 1, 1, 3)
          }
        }
        if (identical(types, c("DFOP", "SFO"))) {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              t <- xidep[, "time"]
              n10 <- psi[id, 1]
              k2 <- exp(psi[id, 2])
              f12 <- plogis(psi[id, 3])
              l1 <- exp(psi[id, 4])
              l2 <- exp(psi[id, 5])
              g <- plogis(psi[id, 6])
              ifelse(xidep[, "name"] == parent_name,
                n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
                ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
                  (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
                  ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
                      ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
                        f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
                  ((l1 - k2) * l2 - k2 * l1 + k2^2)
              )
            }
          } else {
            model_function <- function(psi, id, xidep) {
              t <- xidep[, "time"]
              n10 <- psi[id, 1]
              k2 <- psi[id, 2]
              f12 <- psi[id, 3]
              l1 <- psi[id, 4]
              l2 <- psi[id, 5]
              g <- psi[id, 6]
              ifelse(xidep[, "name"] == parent_name,
                n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
                ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
                  (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
                  ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
                      ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
                        f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
                  ((l1 - k2) * l2 - k2 * l1 + k2^2)
              )
            }
            transform.par = c(0, 1, 3, 1, 1, 3)
          }
        }
      }
    }
  }

  if (is.function(model_function) & solution_type == "auto") {
    solution_type = "analytical saemix"
  } else {

    if (transformations == "saemix") {
      stop("Using saemix transformations is only supported if an analytical solution is implemented for saemix")
    }

    if (solution_type == "auto")
      solution_type <- object[[1]]$solution_type

    # Define some variables to avoid function calls in model function
    transparms_optim_names <- names(degparms_optim)
    odeini_optim_names <- gsub('_0$', '', odeini_optim_parm_names)
    diff_names <- names(mkin_model$diffs)
    ode_transparms_optim_names <- setdiff(transparms_optim_names, odeini_optim_parm_names)
    transform_rates <- object[[1]]$transform_rates
    transform_fractions <- object[[1]]$transform_fractions

    # Get native symbol info for speed
    use_symbols = FALSE
    if (solution_type == "deSolve" & !is.null(mkin_model$cf)) {
      mkin_model$symbols <- try(deSolve::checkDLL(
        dllname = mkin_model$dll_info[["name"]],
        func = "diffs", initfunc = "initpar",
        jacfunc = NULL, nout = 0, outnames = NULL))
      if (!inherits(mkin_model$symbols, "try-error")) {
        use_symbols = TRUE
      }
    }

    # Define the model function
    model_function <- function(psi, id, xidep) {

      uid <- unique(id)

      res_list <- lapply(uid, function(i) {

        transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict
        names(transparms_optim) <- transparms_optim_names

        odeini_optim <- transparms_optim[odeini_optim_parm_names]
        names(odeini_optim) <- odeini_optim_names

        odeini <- c(odeini_optim, odeini_fixed)[diff_names]

        odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model,
          transform_rates = transform_rates,
          transform_fractions = transform_fractions)
        odeparms <- c(odeparms_optim, odeparms_fixed)

        xidep_i <- xidep[which(id == i), ]

        if (solution_type[1] == "analytical") {
          out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms)
        } else {

          i_time <- xidep_i$time
          i_name <- xidep_i$name

          out_wide <- mkinpredict(mkin_model,
            odeparms = odeparms, odeini = odeini,
            solution_type = solution_type,
            outtimes = sort(unique(i_time)),
            na_stop = FALSE
          )

          out_index <- cbind(as.character(i_time), as.character(i_name))
          out_values <- out_wide[out_index]
        }
        return(out_values)
      })
      res <- unlist(res_list)
      return(res)
    }
  }

  if (identical(error_model, "auto")) {
    error_model = object[[1]]$err_mod
  }
  error.model <- switch(error_model,
    const = "constant",
    tc = "combined",
    obs = "constant")

  if (error_model == "obs") {
    warning("The error model 'obs' (variance by variable) can currently not be transferred to an saemix model")
  }

  degparms_psi0 <- degparms_optim
  degparms_psi0[names(degparms_start)] <- degparms_start
  psi0_matrix <- matrix(degparms_psi0, nrow = 1,
    dimnames = list("(Intercept)", names(degparms_psi0)))

  if (covariance.model[1] == "auto") {
    covariance_diagonal <- rep(1, length(degparms_optim))
    if (!is.null(no_random_effect)) {
      degparms_no_random <- which(names(degparms_psi0) %in% no_random_effect)
      covariance_diagonal[degparms_no_random] <- 0
    }
    covariance.model = diag(covariance_diagonal)
  }

  if (omega.init[1] == "auto") {
    if (transformations == "mkin") {
      degparms_eta_ini <- as.numeric( # remove names
        mean_degparms(object,
          random = TRUE, test_log_parms = TRUE)$eta)

      omega.init <- 2 * diag(degparms_eta_ini^2)
    } else {
      omega.init <- matrix(nrow = 0, ncol = 0)
    }
  }

  if (is.null(covariate_models)) {
    covariate.model <- matrix(nrow = 0, ncol = 0) # default in saemixModel()
  } else {
    degparms_dependent <- sapply(covariate_models, function(x) as.character(x[[2]]))
    covariates_in_models = unique(unlist(lapply(
      covariate_models, function(x)
        colnames(attr(terms(x), "factors"))
      )))
    covariates_not_available <- setdiff(covariates_in_models, names(covariates))
    if (length(covariates_not_available) > 0) {
      stop("Covariate(s) ", paste(covariates_not_available, collapse = ", "),
        " used in the covariate models not available in the covariate data")
    }
    psi0_matrix <- rbind(psi0_matrix,
      matrix(0, nrow = length(covariates), ncol = ncol(psi0_matrix),
        dimnames = list(names(covariates), colnames(psi0_matrix))))
    covariate.model <- matrix(0, nrow = length(covariates),
      ncol = ncol(psi0_matrix),
      dimnames = list(
        covariates = names(covariates),
        degparms = colnames(psi0_matrix)))
    if (transformations == "saemix") {
      stop("Covariate models with saemix transformations currently not supported")
    }
    parms_trans <- as.data.frame(t(sapply(object, parms, transformed = TRUE)))
    for (covariate_model in covariate_models) {
      covariate_name <- as.character(covariate_model[[2]])
      model_data <- cbind(parms_trans, covariates)
      ini_model <- lm(covariate_model, data = model_data)
      ini_coef <- coef(ini_model)
      psi0_matrix[names(ini_coef), covariate_name] <- ini_coef
      covariate.model[names(ini_coef)[-1], covariate_name] <- as.numeric(as.logical(ini_coef[-1]))
    }
  }

  res <- saemix::saemixModel(model_function,
    psi0 = psi0_matrix,
    "Mixed model generated from mmkin object",
    transform.par = transform.par,
    error.model = error.model,
    verbose = verbose,
    covariance.model = covariance.model,
    covariate.model = covariate.model,
    omega.init = omega.init,
    error.init = error.init,
    ...
  )
  attr(res, "mean_dp_start") <- degparms_optim
  return(res)
}

#' @rdname saem
#' @importFrom rlang !!!
#' @return An [saemix::SaemixData] object.
#' @export
saemix_data <- function(object, covariates = NULL, verbose = FALSE, ...) {
  if (nrow(object) > 1) stop("Only row objects allowed")
  ds_names <- colnames(object)

  ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")])
  names(ds_list) <- ds_names
  ds_saemix_all <- vctrs::vec_rbind(!!!ds_list, .names_to = "ds")
  ds_saemix <- data.frame(ds = ds_saemix_all$ds,
    name = as.character(ds_saemix_all$variable),
    time = ds_saemix_all$time,
    value = ds_saemix_all$observed,
    stringsAsFactors = FALSE)
  if (!is.null(covariates)) {
    name.covariates <- names(covariates)
    covariates$ds <- rownames(covariates)
    ds_saemix <- merge(ds_saemix, covariates, sort = FALSE)
  } else {
    name.covariates <- character(0)
  }

  res <- saemix::saemixData(ds_saemix,
    name.group = "ds",
    name.predictors = c("time", "name"),
    name.response = "value",
    name.covariates = name.covariates,
    verbose = verbose,
    ...)
  return(res)
}

#' logLik method for saem.mmkin objects
#'
#' @param object The fitted [saem.mmkin] object
#' @param \dots Passed to [saemix::logLik.SaemixObject]
#' @param method Passed to [saemix::logLik.SaemixObject]
#' @export
logLik.saem.mmkin <- function(object, ..., method = c("is", "lin", "gq")) {
  method <- match.arg(method)
  return(logLik(object$so, method = method))
}

#' @export
update.saem.mmkin <- function(object, ..., evaluate = TRUE) {
  call <- object$call
  # For some reason we get saem.mmkin in the call when using mhmkin
  # so we need to fix this so we do not have to export saem.mmkin in
  # addition to the S3 method
  call[[1]] <- saem

  # We also need to provide the mmkin object in the call, so it
  # will also be found when called by testthat or pkgdown
  call[[2]] <- object$mmkin

  update_arguments <- match.call(expand.dots = FALSE)$...

  if (length(update_arguments) > 0) {
    update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
  }

  for (a in names(update_arguments)[update_arguments_in_call]) {
    call[[a]] <- update_arguments[[a]]
  }

  update_arguments_not_in_call <- !update_arguments_in_call
  if(any(update_arguments_not_in_call)) {
    call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
    call <- as.call(call)
  }
  if(evaluate) eval(call, parent.frame())
  else call
}

#' @export
#' @rdname parms
#' @param ci Should a matrix with estimates and confidence interval boundaries
#' be returned? If FALSE (default), a vector of estimates is returned if no
#' covariates are given, otherwise a matrix of estimates is returned, with
#' each column corresponding to a row of the data frame holding the covariates
#' @param covariates A data frame holding covariate values for which to
#' return parameter values. Only has an effect if 'ci' is FALSE.
parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, ...) {
  cov.mod <- object$sm@covariance.model
  n_cov_mod_parms <- sum(cov.mod[upper.tri(cov.mod, diag = TRUE)])
  n_parms <- length(object$sm@name.modpar) +
    n_cov_mod_parms +
    length(object$sm@name.sigma)

  if (inherits(object$so, "try-error")) {
    conf.int <- matrix(rep(NA, 3 * n_parms), ncol = 3)
    colnames(conf.int) <- c("estimate", "lower", "upper")
  } else {
    conf.int <- object$so@results@conf.int[c("estimate", "lower", "upper")]
    rownames(conf.int) <- object$so@results@conf.int[["name"]]
    conf.int.var <- grepl("^Var\\.", rownames(conf.int))
    conf.int <- conf.int[!conf.int.var, ]
    conf.int.cov <- grepl("^Cov\\.", rownames(conf.int))
    conf.int <- conf.int[!conf.int.cov, ]
  }
  estimate <- conf.int[, "estimate"]

  names(estimate) <- rownames(conf.int)

  if (ci) {
    return(conf.int)
  } else {
    if (is.null(covariates)) {
      return(estimate)
    } else {
      est_for_cov <- matrix(NA,
        nrow = length(object$sm@name.modpar), ncol = nrow(covariates),
        dimnames = (list(object$sm@name.modpar, rownames(covariates))))
      covmods <- object$covariate_models
      names(covmods) <- sapply(covmods, function(x) as.character(x[[2]]))
      for (deg_parm_name in rownames(est_for_cov)) {
        if (deg_parm_name %in% names(covmods)) {
          covariate <- covmods[[deg_parm_name]][[3]]
          beta_degparm_name <- paste0("beta_", covariate,
            "(", deg_parm_name, ")")
          est_for_cov[deg_parm_name, ] <- estimate[deg_parm_name] +
            covariates[[covariate]] * estimate[beta_degparm_name]
        } else {
          est_for_cov[deg_parm_name, ] <- estimate[deg_parm_name]
        }
      }
      return(est_for_cov)
    }
  }
}