#' Compute Support Intervals
#'
#' A support interval contains only the values of the parameter that predict the observed data better
#' than average, by some degree \emph{k}; these are values of the parameter that are associated with an
#' updating factor greater or equal than \emph{k}. From the perspective of the Savage-Dickey Bayes factor, testing
#' against a point null hypothesis for any value within the support interval will yield a Bayes factor smaller
#' than \emph{1/k}.
#' \cr \cr
#' \strong{For more info, in particular on specifying correct priors for factors with more than 2 levels,
#' see \href{https://easystats.github.io/bayestestR/articles/bayes_factors.html}{the Bayes factors vignette}.}
#'
#' @param BF The amount of support required to be included in the support interval.
#' @inheritParams bayesfactor_parameters
#' @inheritParams hdi
#'
#' @details This method is used to compute support intervals based on prior and posterior distributions.
#' For the computation of support intervals, the model priors must be proper priors (at the very least
#' they should be \emph{not flat}, and it is preferable that they be \emph{informative} - note
#' that by default, \code{brms::brm()} uses flat priors for fixed-effects; see example below).
#'
#' \subsection{Choosing a value of \code{BF}}{
#' The choice of \code{BF} (the level of support) depends on what we want our interval to represent:
#' \itemize{
#'   \item A \code{BF} = 1 contains values whose credibility is not decreased by observing the data.
#'   \item A \code{BF} > 1 contains values who recived more impressive support from the data.
#'   \item A \code{BF} < 1 contains values whose credibility has \emph{not} been impressively decreased by observing the data.
#'   Testing against values outside this interval will produce a Bayes factor larger than 1/\code{BF} in support of
#'   the alternative. E.g., if an SI (BF = 1/3) excludes 0, the Bayes factor against the point-null will be larger than 3.
#' }
#' }
#'
#' @inheritSection bayesfactor_parameters Setting the correct \code{prior}
#'
#' @note There is also a \href{https://easystats.github.io/see/articles/bayestestR.html}{\code{plot()}-method} implemented in the \href{https://easystats.github.io/see/}{\pkg{see}-package}.
#'
#' @return
#' A data frame containing the lower and upper bounds of the SI.
#' \cr
#' Note that if the level of requested support is higher than observed in the data, the
#' interval will be \code{[NA,NA]}.
#'
#' @examples
#' library(bayestestR)
#'
#' prior <- distribution_normal(1000, mean = 0, sd = 1)
#' posterior <- distribution_normal(1000, mean = .5, sd = .3)
#'
#' si(posterior, prior)
#' \dontrun{
#' # rstanarm models
#' # ---------------
#' library(rstanarm)
#' contrasts(sleep$group) <- contr.bayes # see vingette
#' stan_model <- stan_lmer(extra ~ group + (1 | ID), data = sleep)
#' si(stan_model)
#' si(stan_model, BF = 3)
#'
#' # emmGrid objects
#' # ---------------
#' library(emmeans)
#' group_diff <- pairs(emmeans(stan_model, ~group))
#' si(group_diff, prior = stan_model)
#'
#' # brms models
#' # -----------
#' library(brms)
#' contrasts(sleep$group) <- contr.bayes # see vingette
#' my_custom_priors <-
#'   set_prior("student_t(3, 0, 1)", class = "b") +
#'   set_prior("student_t(3, 0, 1)", class = "sd", group = "ID")
#'
#' brms_model <- brm(extra ~ group + (1 | ID),
#'   data = sleep,
#'   prior = my_custom_priors
#' )
#' si(brms_model)
#' }
#' @references
#' Wagenmakers, E., Gronau, Q. F., Dablander, F., & Etz, A. (2018, November 22). The Support Interval. \doi{10.31234/osf.io/zwnxb}
#'
#' @export
si <- function(posterior, prior = NULL, BF = 1, verbose = TRUE, ...) {
  UseMethod("si")
}

#' @rdname si
#' @export
si.numeric <- function(posterior, prior = NULL, BF = 1, verbose = TRUE, ...) {
  if (is.null(prior)) {
    prior <- posterior
    if (verbose) {
      warning(
        "Prior not specified! ",
        "Please specify a prior (in the form 'prior = distribution_normal(1000, 0, 1)')",
        " to get meaningful results."
      )
    }
  }
  prior <- data.frame(X = prior)
  posterior <- data.frame(X = posterior)

  # Get SIs
  out <- si.data.frame(
    posterior = posterior, prior = prior,
    BF = BF, verbose = verbose, ...
  )
  out$Parameter <- NULL
  out
}

#' @importFrom insight clean_parameters
#' @rdname si
#' @export
si.stanreg <- function(posterior, prior = NULL,
                       BF = 1, verbose = TRUE,
                       effects = c("fixed", "random", "all"),
                       component = c("conditional", "location", "zi", "zero_inflated", "all", "smooth_terms", "sigma", "distributional", "auxiliary"),
                       parameters = NULL,
                       ...) {
  cleaned_parameters <- insight::clean_parameters(posterior)
  effects <- match.arg(effects)
  component <- match.arg(component)

  samps <- .clean_priors_and_posteriors(posterior, prior,
    verbose = verbose,
    effects = effects, component = component,
    parameters = parameters
  )

  # Get SIs
  temp <- si.data.frame(
    posterior = samps$posterior, prior = samps$prior,
    BF = BF, verbose = verbose, ...
  )

  out <- .prepare_output(temp, cleaned_parameters, inherits(posterior, "stanmvreg"))

  attr(out, "object_name") <- .safe_deparse(substitute(posterior))
  class(out) <- class(temp)
  attr(out, "plot_data") <- attr(temp, "plot_data")

  out
}


#' @rdname si
#' @export
si.brmsfit <- si.stanreg


#' @rdname si
#' @export
si.emmGrid <- function(posterior, prior = NULL,
                       BF = 1, verbose = TRUE, ...) {
  samps <- .clean_priors_and_posteriors(posterior, prior,
    verbose = verbose
  )

  # Get SIs
  out <- si.data.frame(
    posterior = samps$posterior, prior = samps$prior,
    BF = BF, verbose = verbose, ...
  )

  attr(out, "object_name") <- .safe_deparse(substitute(posterior))
  out
}

#' @export
si.emm_list <- si.emmGrid


#' @rdname si
#' @export
si.data.frame <- function(posterior, prior = NULL, BF = 1, verbose = TRUE, ...) {
  if (length(BF) > 1) {
    SIs <- lapply(BF, function(i) {
      si(posterior, prior = prior, BF = i, verbose = verbose, ...)
    })
    out <- do.call(rbind, SIs)

    attr(out, "plot_data") <- attr(SIs[[1]], "plot_data")
    class(out) <- unique(c("bayestestR_si", "bayestestR_ci", class(out)))
    return(out)
  }

  if (is.null(prior)) {
    prior <- posterior
    warning(
      "Prior not specified! ",
      "Please specify priors (with column order matching 'posterior')",
      " to get meaningful results."
    )
  }

  sis <- matrix(NA, nrow = ncol(posterior), ncol = 2)
  for (par in seq_along(posterior)) {
    sis[par, ] <- .si(posterior[[par]],
      prior[[par]],
      BF = BF, ...
    )
  }

  out <- data.frame(
    Parameter = colnames(posterior),
    CI = BF,
    CI_low = sis[, 1],
    CI_high = sis[, 2],
    stringsAsFactors = FALSE
  )
  class(out) <- unique(c("bayestestR_si", "see_si", "bayestestR_ci", "see_ci", class(out)))
  attr(out, "plot_data") <- .make_BF_plot_data(posterior, prior, 0, 0, ...)$plot_data

  out
}


#' @export
si.stanfit <- function(posterior, prior = NULL, BF = 1, verbose = TRUE, effects = c("fixed", "random", "all"), ...) {
  si(insight::get_parameters(posterior, effects = effects))
}


#' @importFrom stats median mad na.omit
#' @keywords internal
.si <- function(posterior, prior, BF = 1, extend_scale = 0.05, precision = 2^8, ...) {
  if (!requireNamespace("logspline")) {
    stop("Package \"logspline\" needed for this function to work. Please install it.")
  }

  if (isTRUE(all.equal(prior, posterior))) {
    return(c(NA, NA))
  }

  x <- c(prior, posterior)
  x_range <- range(x)
  x_rangex <- stats::median(x) + 7 * stats::mad(x) * c(-1, 1)
  x_range <- c(
    max(c(x_range[1], x_rangex[1])),
    min(c(x_range[2], x_rangex[2]))
  )

  extension_scale <- diff(x_range) * extend_scale
  x_range <- x_range + c(-1, 1) * extension_scale

  x_axis <- seq(x_range[1], x_range[2], length.out = precision)

  f_prior <- .logspline(prior, ...)
  f_posterior <- .logspline(posterior, ...)
  d_prior <- logspline::dlogspline(x_axis, f_prior)
  d_posterior <- logspline::dlogspline(x_axis, f_posterior)

  relative_d <- d_posterior / d_prior

  crit <- relative_d >= BF

  cp <- rle(c(stats::na.omit(crit)))
  if (length(cp$lengths) > 3) {
    warning("More than 1 SI detected. Plot the result to investigate.", call. = FALSE)
  }

  x_supported <- stats::na.omit(x_axis[crit])
  if (length(x_supported) < 2) {
    return(c(NA, NA))
  } else {
    range(x_supported)
  }
}