#' Inclusion Bayes Factors for testing predictors across Bayesian models
#'
#' The `bf_*` function is an alias of the main function.
#' \cr \cr
#' For more info, see [the Bayes factors vignette](https://easystats.github.io/bayestestR/articles/bayes_factors.html).
#'
#' @author Mattan S. Ben-Shachar
#' @param models An object of class [bayesfactor_models()] or `BFBayesFactor`.
#' @param match_models See details.
#' @param prior_odds Optional vector of prior odds for the models. See
#'   `BayesFactor::priorOdds<-`.
#' @param ... Arguments passed to or from other methods.
#'
#' @return a data frame containing the prior and posterior probabilities, and
#'   log(BF) for each effect (Use `as.numeric()` to extract the non-log Bayes
#'   factors; see examples).
#'
#' @details Inclusion Bayes factors answer the question: Are the observed data
#'   more probable under models with a particular effect, than they are under
#'   models without that particular effect? In other words, on average - are
#'   models with effect \eqn{X} more likely to have produced the observed data
#'   than models without effect \eqn{X}?
#'
#' \subsection{Match Models}{
#' If `match_models=FALSE` (default), Inclusion BFs are computed by comparing
#' all models with a term against all models without that term. If `TRUE`,
#' comparison is restricted to models that (1) do not include any interactions
#' with the term of interest; (2) for interaction terms, averaging is done only
#' across models that containe the main effect terms from which the interaction
#' term is comprised.
#' }
#'
#' @inheritSection bayesfactor_parameters Interpreting Bayes Factors
#'
#' @note Random effects in the `lmer` style are converted to interaction terms:
#' i.e., `(X|G)` will become the terms `1:G` and `X:G`.
#'
#' @seealso [weighted_posteriors()] for Bayesian parameter averaging.
#'
#' @examplesIf require("BayesFactor")
#' library(bayestestR)
#'
#' # Using bayesfactor_models:
#' # ------------------------------
#' mo0 <- lm(Sepal.Length ~ 1, data = iris)
#' mo1 <- lm(Sepal.Length ~ Species, data = iris)
#' mo2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)
#' mo3 <- lm(Sepal.Length ~ Species * Petal.Length, data = iris)
#'
#' BFmodels <- bayesfactor_models(mo1, mo2, mo3, denominator = mo0)
#' (bf_inc <- bayesfactor_inclusion(BFmodels))
#'
#' as.numeric(bf_inc)
#'
#' \donttest{
#' # BayesFactor
#' # -------------------------------
#' BF <- BayesFactor::generalTestBF(len ~ supp * dose, ToothGrowth, progress = FALSE)
#' bayesfactor_inclusion(BF)
#'
#' # compare only matched models:
#' bayesfactor_inclusion(BF, match_models = TRUE)
#' }
#'
#' @references
#' - Hinne, M., Gronau, Q. F., van den Bergh, D., and Wagenmakers, E. (2019, March 25).
#'   A conceptual introduction to Bayesian Model Averaging. \doi{10.31234/osf.io/wgb64}
#'
#' - Clyde, M. A., Ghosh, J., & Littman, M. L. (2011). Bayesian adaptive sampling
#'   for variable selection and model averaging. Journal of Computational and Graphical Statistics,
#'   20(1), 80-101.
#'
#' - Mathot, S. (2017). Bayes like a Baws: Interpreting Bayesian Repeated Measures in JASP.
#'   [Blog post](https://www.cogsci.nl/blog/interpreting-bayesian-repeated-measures-in-jasp).
#'
#' @export
bayesfactor_inclusion <- function(models,
                                  match_models = FALSE,
                                  prior_odds = NULL,
                                  ...) {
  UseMethod("bayesfactor_inclusion")
}

#' @rdname bayesfactor_inclusion
#' @export
bf_inclusion <- bayesfactor_inclusion

#' @export
bayesfactor_inclusion.bayesfactor_models <- function(models,
                                                     match_models = FALSE,
                                                     prior_odds = NULL,
                                                     ...) {
  if (isTRUE(attr(models, "unsupported_models"))) {
    insight::format_error(
      "Can not compute inclusion Bayes factors - passed models are not (yet) supported."
    )
  }

  # Build Models Table #
  df.model <- .get_model_table(models, priorOdds = prior_odds)
  effnames <- colnames(df.model)[-(1:3)]

  # Build Interaction Matrix #
  if (isTRUE(match_models)) {
    effects.matrix <- as.matrix(df.model[, -(1:3)])

    df.interaction <- data.frame(effnames, stringsAsFactors = FALSE)

    for (eff in effnames) {
      df.interaction[, eff] <- sapply(effnames, .includes_interaction, effnames = eff)
    }
    rownames(df.interaction) <- effnames
    df.interaction <- as.matrix(df.interaction[, -1])
  }

  # Build Effect Table #
  df.effect <- data.frame(
    effnames,
    Pinc = rep(NA, length(effnames)),
    PincD = rep(NA, length(effnames)),
    log_BF = rep(NA, length(effnames)),
    stringsAsFactors = FALSE
  )

  for (eff in effnames) {
    if (isTRUE(match_models)) {
      idx1 <- df.interaction[eff, ]
      idx2 <- df.interaction[, eff]

      has_not_high_order_interactions <- !apply(effects.matrix[, idx1, drop = FALSE], 1, any)

      ind_include <- has_not_high_order_interactions & effects.matrix[, eff]

      ind_exclude <- apply(effects.matrix[, idx2, drop = FALSE], 1, all) &
        has_not_high_order_interactions &
        !effects.matrix[, eff]

      df.model_temp <- df.model[ind_include | ind_exclude, , drop = FALSE]
    } else {
      df.model_temp <- df.model
    }

    # models with effect
    mwith <- which(df.model_temp[[eff]])
    mwithprior <- sum(df.model_temp[mwith, "priorProbs"])
    mwithpost <- sum(df.model_temp[mwith, "postProbs"])

    # models without effect
    mwithoutprior <- sum(df.model_temp[-mwith, "priorProbs"])
    mwithoutpost <- sum(df.model_temp[-mwith, "postProbs"])

    # Save results
    df.effect$Pinc[effnames == eff] <- mwithprior
    df.effect$PincD[effnames == eff] <- mwithpost
    df.effect$log_BF[effnames == eff] <- (log(mwithpost) - log(mwithoutpost)) - (log(mwithprior) - log(mwithoutprior))
  }

  df.effect <- df.effect[, -1, drop = FALSE]
  colnames(df.effect) <- c("p_prior", "p_posterior", "log_BF")
  rownames(df.effect) <- effnames

  class(df.effect) <- c("bayesfactor_inclusion", class(df.effect))
  attr(df.effect, "matched") <- match_models
  attr(df.effect, "priorOdds") <- prior_odds

  return(df.effect)
}


#' @export
bayesfactor_inclusion.BFBayesFactor <- function(models,
                                                match_models = FALSE,
                                                prior_odds = NULL,
                                                ...) {
  models <- bayesfactor_models.BFBayesFactor(models)
  bayesfactor_inclusion.bayesfactor_models(models,
    match_models = match_models,
    prior_odds = prior_odds
  )
}


#' @keywords internal
.includes_interaction <- function(eff, effnames) {
  eff_b <- strsplit(eff, ":", fixed = TRUE)
  effnames_b <- strsplit(effnames, ":", fixed = TRUE)

  is_int <- vapply(effnames_b, function(x) length(x) > 1, TRUE)

  temp <- logical(length(effnames))

  for (rr in seq_along(effnames)) {
    if (is_int[rr]) {
      temp[rr] <- all(eff_b[[1]] %in% effnames_b[[rr]]) &
        !all(effnames_b[[rr]] %in% eff_b[[1]])
    }
  }

  temp
}