test_that("multicollinearity vifcor", {
  skip_on_cran()
  set.seed(1)
  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- dplyr::select(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(div=bio1/bio12, prod=bio1*bio12, sub=bio12-bio1, soma=bio1+bio12)
  oc <- occurrences_sdm(occ, crs=6933)
  i <- input_sdm(oc, sa)
  th=0.8
  suppressWarnings(i <- multicollinearity_sdm(i, method = "vifcor", th=th))
  expect_true("variable_selection" %in% names(i$predictors))
  expect_true("vif" %in% names(i$predictors$variable_selection))
  expect_equal(i$predictors$variable_selection$vif$threshold, th)
  expect_true(any(get_predictor_names(i) %in% i$predictors$variable_selection$vif$selected_variables))
  expect_equal(selected_variables(i), c("div", "prod"))
  expect_snapshot(vif_summary(i))

  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- dplyr::select(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(div=bio1/bio12, prod=bio1*bio12, sub=bio12-bio1, soma=bio1+bio12)
  oc <- occurrences_sdm(occ, crs=6933)
  i <- input_sdm(oc, sa)
  expect_error(multicollinearity_sdm(i, method = "vifcor", th=NA))
  expect_error(multicollinearity_sdm(i, method = "vifcor", th="a"))
  expect_error(multicollinearity_sdm(i, method = "vifcor", variables_selected = 'bio99'))
  expect_error(multicollinearity_sdm(i$occurrences))
  expect_error(multicollinearity_sdm(i$predictors))
  suppressWarnings(i <- multicollinearity_sdm(i))
  expect_error(selected_variables(i))
  expect_error(vif_summary(i))

  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- dplyr::select(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(test=c(rep("a", 19), rep("b", 10)))
  oc <- occurrences_sdm(occ, crs=6933)
  oc <- join_area(oc, sa)
  i <- input_sdm(oc, sa)
  th=0.9
  suppressWarnings(i <- multicollinearity_sdm(i, method = "vifcor", th=th))
  expect_true("variable_selection" %in% names(i$predictors))
  expect_true("vif" %in% names(i$predictors$variable_selection))
  expect_equal(i$predictors$variable_selection$vif$threshold, th)
  expect_true(any(get_predictor_names(i) %in% i$predictors$variable_selection$vif$selected_variables))
  expect_equal(selected_variables(i), c("bio1", "bio12", "test"))
  expect_snapshot(vif_summary(i))
})

test_that("multicollinearity vifstep", {
  skip_on_cran()
  set.seed(1)
  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- dplyr::select(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(div=bio1/bio12, prod=bio1*bio12, sub=bio12-bio1, soma=bio1+bio12)
  oc <- occurrences_sdm(occ, crs=6933)
  i <- input_sdm(oc, sa)
  th=5
  suppressWarnings(i <- multicollinearity_sdm(i, method = "vifstep", th=th))
  expect_true("variable_selection" %in% names(i$predictors))
  expect_true("vif" %in% names(i$predictors$variable_selection))
  expect_equal(i$predictors$variable_selection$vif$threshold, th)
  expect_true(any(get_predictor_names(i) %in% i$predictors$variable_selection$vif$selected_variables))
  expect_equal(selected_variables(i), c("div", "prod"))
  expect_snapshot(vif_summary(i))

  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- dplyr::select(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(div=bio1/bio12, prod=bio1*bio12, sub=bio12-bio1, soma=bio1+bio12)
  oc <- occurrences_sdm(occ, crs=6933)
  i <- input_sdm(oc, sa)
  expect_error(multicollinearity_sdm(i, method = "vifstep", th=NA))
  expect_error(multicollinearity_sdm(i, method = "vifstep", th="a"))
  expect_error(multicollinearity_sdm(i, method = "vifstep", variables_selected = 'bio99'))
  expect_error(multicollinearity_sdm(i$occurrences))
  expect_error(multicollinearity_sdm(i$predictors))
  suppressWarnings(i <- multicollinearity_sdm(i))
  expect_error(selected_variables(i))
  expect_error(vif_summary(i))

  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- dplyr::select(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(test=c(rep("a", 19), rep("b", 10)))
  oc <- occurrences_sdm(occ, crs=6933)
  oc <- join_area(oc, sa)
  i <- input_sdm(oc, sa)
  th=5
  suppressWarnings(i <- multicollinearity_sdm(i, method = "vifstep", th=th))
  expect_true("variable_selection" %in% names(i$predictors))
  expect_true("vif" %in% names(i$predictors$variable_selection))
  expect_equal(i$predictors$variable_selection$vif$threshold, th)
  expect_true(any(get_predictor_names(i) %in% i$predictors$variable_selection$vif$selected_variables))
  expect_equal(selected_variables(i), c("bio1", "bio12", "test"))
  expect_snapshot(vif_summary(i))
})

test_that("multicollinearity pca", {
  skip_on_cran()
  set.seed(1)
  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- select_predictors(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(div=bio1/bio12, prod=bio1*bio12, sub=bio12-bio1, soma=bio1+bio12)
  sa <- add_scenarios(sa)
  oc <- occurrences_sdm(occ, crs=6933)
  i <- input_sdm(oc, sa)
  expect_error(get_pca_model(i))
  i <- multicollinearity_sdm(i, method = "pca", cumulative_proportion = 0.99)
  expect_true("variable_selection" %in% names(i$predictors))
  expect_true("pca" %in% names(i$predictors$variable_selection))
  expect_true(any(get_predictor_names(i)[-c(1:6)] %in% i$predictors$variable_selection$pca$selected_variables))
  expect_true(all(get_predictor_names(i)[-c(1:6)] %in% colnames(i$scenarios$data$current)))

  val <- pca_summary(i)
  expect_equal(val$importance[,1] |> as.numeric() |> round(4), c(573.37600, 0.99440, 0.99440))
  expect_equal(val$sdev |> round(3), c(573.376, 42.846, 0.259, 0.000, 0.000, 0.000))
  expect_equal(as.numeric(val$rotation[,1]) |> round(4), c(-0.00170, -0.02410, -0.00000, -0.99910, -0.02240, -0.02580))
  expect_equal(as.numeric(val$x[1,]) |> round(4), c(-19.95450, 73.23320, 0.45290, -0.00020, 0.00000, -0.00000))
  expect_equal(as.numeric(val$center) |> round(3), c(19.2340, 299.4870, 0.0650, 5734.4410, 280.2530, 318.7220))
  expect_equal(class(get_pca_model(i)), "prcomp")

  # No scenarios
  set.seed(1)
  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- select_predictors(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(div=bio1/bio12, prod=bio1*bio12, sub=bio12-bio1, soma=bio1+bio12)
  oc <- occurrences_sdm(occ, crs=6933)
  i <- input_sdm(oc, sa)
  expect_error(get_pca_model(i))
  i <- multicollinearity_sdm(i, method = "pca", cumulative_proportion = 0.99)
  expect_true("variable_selection" %in% names(i$predictors))
  expect_true("pca" %in% names(i$predictors$variable_selection))
  expect_true(any(get_predictor_names(i)[-c(1:6)] %in% i$predictors$variable_selection$pca$selected_variables))
  expect_true(all(get_predictor_names(i)[-c(1:6)] %in% colnames(i$scenarios$data$current)))

  val <- pca_summary(i)
  expect_equal(val$importance[,1] |> as.numeric() |> round(4), c(573.37600, 0.99440, 0.99440))
  expect_equal(val$sdev |> round(3), c(573.376, 42.846, 0.259, 0.000, 0.000, 0.000))
  expect_equal(as.numeric(val$rotation[,1]) |> round(4), c(-0.00170, -0.02410, -0.00000, -0.99910, -0.02240, -0.02580))
  expect_equal(as.numeric(val$x[1,]) |> round(4), c(-19.95450, 73.23320, 0.45290, -0.00020, 0.00000, -0.00000))
  expect_equal(as.numeric(val$center) |> round(3), c(19.2340, 299.4870, 0.0650, 5734.4410, 280.2530, 318.7220))
  expect_equal(class(get_pca_model(i)), "prcomp")
})

test_that("multicollinearity custom function", {
  skip_on_cran()
  set.seed(1)
  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- select_predictors(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(div=bio1/bio12, prod=bio1*bio12, sub=bio12-bio1, soma=bio1+bio12)
  sa <- add_scenarios(sa)
  oc <- occurrences_sdm(occ, crs=6933)
  i <- input_sdm(oc, sa)
  expect_error(get_pca_model(i))
  custom_function <- function(env_sf, occ_sf) {
    env_df <- dplyr::select(sf::st_drop_geometry(env_sf), -"cell_id")
    correlations <- cor(env_df)
    col <- caret::findCorrelation(correlations, cutoff = 0.7)
    selected <- colnames(correlations)[-col]
    return(selected)
  }
  i <- multicollinearity_sdm(i, method = custom_function)
  expect_true("variable_selection" %in% names(i$predictors))
  expect_true("custom_function" %in% names(i$predictors$variable_selection))
  expect_true(any(get_predictor_names(i) %in% i$predictors$variable_selection$custom_function$selected_variables))
  expect_true(all(get_predictor_names(i) %in% colnames(i$scenarios$data$current)))

  # No scenarios
  sa <- sdm_area(parana, 1)
  sa <- add_predictors(sa, bioc)
  sa <- select_predictors(sa, c("bio1", "bio12"))
  sa <- sa |> dplyr::mutate(div=bio1/bio12, prod=bio1*bio12, sub=bio12-bio1, soma=bio1+bio12)
  oc <- occurrences_sdm(occ, crs=6933)
  i <- input_sdm(oc, sa)
  expect_error(get_pca_model(i))
  i <- multicollinearity_sdm(i, method = custom_function)
  expect_true("variable_selection" %in% names(i$predictors))
  expect_true("custom_function" %in% names(i$predictors$variable_selection))
  expect_true(any(get_predictor_names(i) %in% i$predictors$variable_selection$custom_function$selected_variables))
})