skip_on_cran()
skip_if_not_installed("lixoftConnectors")


initRsSimulx()
demo_path <- gsub("simulx", "monolix", .lixoftCall("getDemoPath"), fixed = TRUE)
skip_if(!dir.exists(demo_path), message = NULL)

project <- file.path(demo_path, "1.creating_and_using_models",
                     "1.1.libraries_of_models", "theophylline_project.mlxtran")
project <- get_project(project, runStdErrorsIfNeed = TRUE)

project_cov <- paste0(tools::file_path_sans_ext(project), "_cov.mlxtran")
.loadProject(project, software = "monolix")
.lixoftCall("setCovariateModel", list(list(V=c(SEX=T))))
.lixoftCall("saveProject", list(project_cov))
.lixoftCall("runScenario")

initRsSimulx(force = TRUE)

test_that("Simulx runs with parameters automatically generated by monolix", {
  res <- simulx(project = project,
                parameter = "mlx_Pop")
  param <- unlist(.lixoftCall("getPopulationElements")$mlx_Pop$data)
  expect_equal(res$population, param)

  res <- simulx(project = project,
                nrep=10,
                parameter = "mlx_PopUncertainLin")
  expect_equal(nrow(res$population), 10)
  
  res <- simulx(project = project,
                parameter = "mlx_PopIndiv")
  param <- unlist(.lixoftCall("getIndividualElements")$mlx_PopIndiv$data)
  param <- data.frame(
    id = as.factor(seq(1, 12)),
    original_id = seq(1, 12),
    ka = param[["ka"]],
    V = param[["V"]],
    Cl = param[["Cl"]]
  )
  expect_equal(subset(res$parameter, select = -c(WEIGHT, SEX)), param, tolerance = 1.e-5)

  # currently get warning (CONNECTORS-658)
  res <- simulx(project = project,
                parameter = "mlx_CondMean",
                settings = list(sharedIds = c("treatment", "individual", "covariate", "output")))
  param <- .renameColumns(.lixoftCall("getIndividualElements")$mlx_CondMean$data, "ID", "id")
  param$id <- factor(as.integer(param$id))
  expect_equal(subset(res$parameter, select = -c(original_id, WEIGHT, SEX)), param, tolerance = 1.e-5)
  
  # currently get warning (CONNECTORS-658)
  expect_message(
    res <- simulx(project = project,
                  parameter = "mean"),
    "'mode', 'mean', 'pop' are deprecated parameters."
  )
  param <- .renameColumns(.lixoftCall("getIndividualElements")$mlx_CondMean$data, "ID", "id")
  param$id <- factor(as.integer(param$id))
  expect_equal(subset(res$parameter, select = -c(original_id, WEIGHT, SEX)), param, tolerance = 1.e-5)

  # currently get warning (CONNECTORS-658)
  res <- simulx(project = project,
                parameter = "mlx_EBEs")
  param <- .renameColumns(.lixoftCall("getIndividualElements")$mlx_EBEs$data, "ID", "id")
  param$id <- factor(as.integer(param$id))
  expect_equal(subset(res$parameter, select = -c(original_id, WEIGHT, SEX)), param, tolerance = 1.e-5)

  # currently get warning (CONNECTORS-658)
  expect_message(
    res <- simulx(project = project,
                  parameter = "mode"),
    "'mode', 'mean', 'pop' are deprecated parameters."
  )
  param <- .renameColumns(.lixoftCall("getIndividualElements")$mlx_EBEs$data, "ID", "id")
  param$id <- factor(as.integer(param$id))
  expect_equal(subset(res$parameter, select = -c(original_id, WEIGHT, SEX)), param, tolerance = 1.e-5)

  # currently get warning (CONNECTORS-658)
  res <- simulx(project = project,
                parameter = "mlx_CondDistSample")
  param <- .renameColumns(.lixoftCall("getIndividualElements")$mlx_CondDistSample$data, "ID", "id")
  param$id <- factor(as.integer(param$id))
  expect_equal(subset(res$parameter, select = -c(original_id, WEIGHT, SEX)), param, tolerance = 1.e-5)
  
})

test_that("Groups with parameters automatically generated by monolix", {
  res <- simulx(project = project,
                nrep=2,
                group = list(list(parameter = "mlx_PopUncertainLin"),
                             list(parameter = "mlx_Pop")))
  mlxpop <- .lixoftCall("getPopulationElements")$mlx_Pop$data
  pop_g1 <- res$population[res$population$group == 1,]
  rownames(pop_g1) <- NULL
  pop_g2 <- res$population[res$population$group == 2,]
  rownames(pop_g2) <- NULL
  expect_equal(nrow(res$population), 4)
  expect_equal(nrow(unique(pop_g1[names(pop_g1) != "rep"])), 2)
  expect_equal(nrow(unique(pop_g2[names(pop_g2) != "rep"])), 1)
  expect_equal(unique(pop_g2[! names(pop_g2) %in% c("rep", "group")]), mlxpop)

  # currently get warning (CONNECTORS-658)
  res <- simulx(project = project,
                group = list(list(parameter = "mlx_PopIndiv"),
                             list(parameter = "mlx_CondMean"),
                             list(parameter = "mlx_EBEs"),
                             list(parameter = "mlx_CondDistSample")),
                settings = list(sharedIds = c("treatment", "individual", "covariate", "output")))
  popindiv_param <- .lixoftCall("getIndividualElements")$mlx_PopIndiv$data
  condmean_param <- .renameColumns(.lixoftCall("getIndividualElements")$mlx_CondMean$data, "ID", "id")
  condmean_param$original_id <- factor(as.integer(condmean_param$id))
  ebes_param <- .renameColumns(.lixoftCall("getIndividualElements")$mlx_EBEs$data, "ID", "id")
  ebes_param$original_id <- factor(as.integer(ebes_param$id))
  conddist_param <- .renameColumns(.lixoftCall("getIndividualElements")$mlx_CondDistSample$data, "ID", "id")
  conddist_param$original_id <- factor(as.integer(conddist_param$id))
  
  p1 <- subset(res$parameter[res$parameter$group == 1,], select=-c(group))
  p2 <- subset(res$parameter[res$parameter$group == 2,], select=-c(group))
  p3 <- subset(res$parameter[res$parameter$group == 3,], select=-c(group))
  p4 <- subset(res$parameter[res$parameter$group == 4,], select=-c(group))
  rownames(p1) <- rownames(p2) <- rownames(p3) <- rownames(p4) <- NULL

  expect_equal(nrow(unique(subset(p1, select=-c(id, original_id, WEIGHT, SEX)))), 1)
  expect_equal(unique(subset(p1, select=-c(id, original_id, WEIGHT, SEX))), popindiv_param, tolerance = 1.e-4)
  expect_equal(subset(p2, select=-c(id, original_id, WEIGHT, SEX)), subset(condmean_param, select=-c(id, original_id)), tolerance = 1.e-4)
  expect_equal(subset(p3, select=-c(id, original_id, WEIGHT, SEX)), subset(ebes_param, select=-c(id, original_id)), tolerance = 1.e-4)
  expect_equal(subset(p4, select=-c(id, original_id, WEIGHT, SEX)), subset(conddist_param, select=-c(id, original_id)), tolerance = 1.e-4)

})

test_that("Simulx runs with covariates automatically generated by monolix", {
  res <- simulx(project = project_cov,
                covariate = "mlx_Cov")
  cov <- .renameColumns(.lixoftCall("getCovariateElements")$mlx_Cov$data, "ID", "id")
  cov$id <- factor(as.numeric(cov$id))
  expect_equal(res$parameter[c("id", "WEIGHT", "SEX")], cov)
  
  res <- simulx(project = project_cov,
                group = list(size=5000),
                covariate = "mlx_CovDist")
  cov <- .renameColumns(.lixoftCall("getCovariateElements")$mlx_CovDist$data, "ID", "id")
  prob <- table(res$parameter$SEX) / nrow(res$parameter)
  expect_equal(prob[["F"]], cov$probabilities[which(cov$modalities == "F")], tolerance = 5e-2)
  
})

test_that("Groups with covariates automatically generated by monolix", {
  res <- simulx(project = project_cov,
                group = list(list(covariate="mlx_Cov"), list(size=5000, covariate="mlx_CovDist")))
  cov <- .renameColumns(.lixoftCall("getCovariateElements")$mlx_Cov$data, "ID", "id")
  covdist <- .lixoftCall("getCovariateElements")$mlx_CovDist$data
  
  c1 <- res$parameter[res$parameter$group == 1, c("id", "SEX")]
  c2 <- res$parameter[res$parameter$group == 2, c("id", "SEX")]
  row.names(c1) <- row.names(c2) <- NULL
  prob2 <- table(c2$SEX) / nrow(c2)
  
  expect_equal(subset(c1, select=SEX), subset(cov, select=SEX))
  expect_equal(prob2[["F"]], covdist$probabilities[which(covdist$modalities == "F")], tolerance = 5e-2)
  
  res <- simulx(project = project_cov,
                group = list(list(size=2000, covariate="mlx_CovDist"), list(size=5000, covariate="mlx_CovDist")))
  covdist <- .lixoftCall("getCovariateElements")$mlx_CovDist$data
  c1 <- res$parameter[res$parameter$group == 1, c("id", "SEX")]
  c2 <- res$parameter[res$parameter$group == 2, c("id", "SEX")]
  row.names(c1) <- row.names(c2) <- NULL
  prob1 <- table(c1$SEX) / nrow(c1)
  prob2 <- table(c2$SEX) / nrow(c2)

  expect_equal(prob1[["F"]], covdist$probabilities[which(covdist$modalities == "F")], tolerance = 5e-2)
  expect_equal(prob2[["F"]], covdist$probabilities[which(covdist$modalities == "F")], tolerance = 5e-2)

})

test_that("Simulx runs with treatments automatically generated by monolix", {
  res <- simulx(project = project,
                treatment = "mlx_Adm1")
  trt <- .renameColumns(.lixoftCall("getTreatmentElements")$mlx_Adm1$data, "ID", "id")
  trt$id <- factor(as.numeric(trt$id))
  trt$washout <- as.numeric(trt$washout)
  trt$admtype <- 1
  trt <- trt[c("id", "time", "amount", "admtype", "washout")]
  res$treatment$original_id <- NULL
  expect_equal(res$treatment, trt)

})

test_that("Define Groups with different Treatments", {
  initRsSimulx(force = TRUE)

  # defining two different dosing
  adm1 <- list(amount = 25*120, time=0, rate = 25)
  adm2 <- list(amount = 50*120, time=0, rate = 50)
  
  # defining the number of individuals for each group
  N <- 100
  
  # defining the groups
  g1 <- list(size=N, treatment=adm1)
  g2 <- list(size=N, treatment=adm2)
  
  res  <- simulx(project = project,
                 group = list(g1,g2))
  expect_true("group" %in% names(res$treatment))
  expect_equal(levels(res$treatment$group), c("1", "2"))
  expect_equal(nrow(res$treatment[res$treatment$group == 1,]), N)
  expect_equal(nrow(res$treatment[res$treatment$group == 2,]), N)
})