#' Kernel density estimatior based on simplified vine copulas
#'
#' Implements the vine-copula based estimator of Nagler and Czado (2016). The
#' marginal densities are estimated by \code{\link{kde1d}}, the vine copula
#' density by \code{\link{kdevinecop}}. Discrete variables are convoluted with
#' the uniform distribution (see, Nagler, 2017). If a variable should be treated
#' as discrete, declare it as [ordered()]. Factors are expanded into binary
#' dummy codes.
#'
#' @param x (\eqn{n x d}) data matrix.
#' @param mult_1d numeric; all bandwidhts for marginal kernel density estimation
#'   are multiplied with \code{mult_1d}. Defaults to `log(1 + d)` where `d` is
#'   the number of variables after applying [cctools::expand_as_numeric()].
#' @param xmin numeric vector of length d; see \code{\link{kde1d}}.
#' @param xmax numeric vector of length d; see \code{\link{kde1d}}.
#' @param copula.type either \code{"kde"} (default) or \code{"parametric"} for
#'   kernel or parametric estimation of the vine copula.
#' @param ... further arguments passed to \code{\link{kde1d}} or
#'   \code{\link{kdevinecop}}.
#'
#' @return An object of class \code{kdevine}.
#'
#' @seealso \code{\link{dkdevine}} \code{\link{kde1d}} \code{\link{kdevinecop}}
#'
#' @references Nagler, T., Czado, C. (2016) *Evading the curse of
#'   dimensionality in nonparametric density estimation with simplified vine
#'   copulas.* Journal of Multivariate Analysis 151, 69-89
#'   (doi:10.1016/j.jmva.2016.07.003) \cr \cr
#'   Nagler, T. (2017). *A generic approach to nonparametric function
#'   estimation with mixed data.* [arXiv:1704.07457](https://arxiv.org/abs/1704.07457)

#'
#' @examples
#' # load data
#' data(wdbc, package = "kdecopula")
#' \dontshow{wdbc <- wdbc[1:30, ]}
#' # estimate density (use xmin to indicate positive support)
#' fit <- kdevine(wdbc[, 5:7], xmin = rep(0, 3))
#'
#' # evaluate density estimate
#' dkdevine(c(1000, 0.1, 0.1), fit)
#'
#' # plot simulated data
#' pairs(rkdevine(nrow(wdbc), fit))
#'
#' @importFrom VineCopula RVineStructureSelect RVineCopSelect
#' @importFrom cctools expand_vec
#' @export
kdevine <- function(x, mult_1d = NULL, xmin = NULL,
                    xmax = NULL, copula.type = "kde", ...) {
    if (missing(x) & !is.null(list(...)$data))  # for backwards compatibilitiy
        x <- list(...)$data
    if (!is.null(list(...)$mult.1d))  # for backwards compatibilitiy
        mult_1d <- list(...)$mult.1d
    x_cc <- cont_conv(x)
    if (NCOL(x_cc) == 1)
        stop("x must be multivariate or a factor.")
    d <- ncol(x_cc)

    ## sanity checks
    if (!is.null(xmin)) {
        xmin <- cctools::expand_vec(xmin, x)
    }
    if (!is.null(xmax)) {
        xmax <- cctools::expand_vec(xmax, x)
    }
    bw <- list(...)$bw
    if (!is.null(bw)) {
        bw <- cctools::expand_vec(bw, x)
    }

    ## estimation of the marginals
    i_disc <- attr(x_cc, "i_disc")
    if (is.null(mult_1d))
        mult_1d <- log(1 + ncol(x_cc))
    marg.dens <- as.list(numeric(d))
    for (k in 1:d) {
        marg.dens[[k]] <- kde1d(
            x_cc[, k],
            xmin = xmin[k],
            xmax = xmax[k],
            bw   = bw[k],
            mult = mult_1d,
            bw_min = ifelse(k %in% i_disc, 0.5 - attr(x_cc, "theta"), 0)
        )
    }
    res <- list(x_cc = x_cc, marg.dens = marg.dens)

    ## estimation of the R-vine copula (only if d > 1)
    if (d > 1) {
        # transform to copula data
        u <- sapply(1:d, function(k) pkde1d(x_cc[, k], marg.dens[[k]]))

        if (copula.type == "kde") {
            res$vine  <- suppressWarnings(
                kdevinecop(
                    u,
                    matrix      = list(...)$matrix,
                    method      = list(...)$method,
                    mult        = list(...)$mult,
                    info        = list(...)$info,
                    test.level  = list(...)$test.level,
                    trunc.level = list(...)$trunc.level,
                    treecrit    = list(...)$treecrit,
                    cores       = list(...)$cores
                )
            )
        } else if (copula.type == "parametric") {
            # get family and matrix if available
            fam <- list(...)$familyset
            if (is.null(fam))
                fam <- NA
            mat <- list(...)$Matrix

            # fit parametric vine
            res$vine <- if (is.null(mat) & d > 2) {
                # structure selection if no matrix is provided and d > 2
                RVineStructureSelect(u, familyset = fam)
            } else if (d == 2) {
                # select copula for default structure if d = 2
                RVineCopSelect(u,
                               familyset = fam,
                               Matrix = matrix(c(2, 1, 0, 1), 2, 2))
            } else {
                # or select copulas to provided structure
                RVineCopSelect(u,
                               familyset = fam,
                               Matrix = mat)
            }
        } else {
            stop("copula.type not implemented.")
        }
    }

    ## return results
    res$copula.type <- copula.type
    class(res) <- "kdevine"
    res
}

#' Evaluate the density of a kdevine object
#'
#' @param x (\eqn{m x d}) matrix of evaluation points (or vector of length \eqn{d}).
#' @param obj a \code{kdevine} object.
#'
#' @return The density estimate evaluated at \code{x}.
#'
#' @seealso
#' \code{\link{kdevine}}
#'
#' @examples
#' # load data
#' data(wdbc)
#' \dontshow{wdbc <- wdbc[1:30, ]}
#' # estimate density (use xmin to indicate positive support)
#' fit <- kdevine(wdbc[, 5:7], xmin = rep(0, 3))
#'
#' # evaluate density estimate
#' dkdevine(c(1000, 0.1, 0.1), fit)
#'
#' @importFrom VineCopula RVinePDF
#'
#' @export
dkdevine <- function(x, obj) {
    stopifnot(class(obj) == "kdevine")
    if (NCOL(x) == 1)
        x <- t(x)
    nms <- colnames(x)
    # must be numeric, factors are expanded
    x <- expand_as_numeric(x)
    # variables must be in same order
    if (!is.null(nms))
        x <- x[, colnames(obj$x_cc), drop = FALSE]

    ## evaluate marginals
    d <- ncol(x)
    margvals <- u <- x
    for (k in 1:d) {
        x_k <- x[, k]
        if (k %in% attr(obj$x_cc, "i_disc")) {
            # use normalization if discrete
            attr(x_k, "i_disc") <- 1
            obj$marg.dens[[k]]$levels <- attr(obj$x_cc, "levels")[[k]]
        }
        margvals[, k] <- dkde1d(x_k, obj$marg.dens[[k]])
    }

    if (!is.null(obj$vine)) {
        # PIT to copula level
        for (k in 1:d) {
            if (k %in% attr(obj$x_cc, "i_disc")) {
                # use continuous variant for PIT
                attr(x_k, "i_disc") <- integer(0)
                obj$marg.dens[[k]]$levels <- NULL
            }
            u[, k] <- pkde1d(x[, k], obj$marg.dens[[k]])
        }
        if (inherits(obj$vine, "kdevinecop")) {
            vinevals <- dkdevinecop(u, obj = obj$vine, stable = TRUE)
        } else if (inherits(obj$vine, "RVineMatrix")) {
            vinevals <- RVinePDF(u, obj$vine)
        } else {
            stop("vine has incompatible type")
        }
    } else {
        vinevals <- rep(1, nrow(x))
    }

    ## final density estimate is product of marginals and copula density
    apply(cbind(margvals, vinevals), 1, prod)
}


#' Simulate from a kdevine object
#'
#' @param n number of observations.
#' @param obj a \code{kdevine} object.
#' @param quasi logical; the default (\code{FALSE}) returns pseudo-random
#' numbers, use \code{TRUE} for quasi-random numbers (generalized Halton, only
#' works for fully nonparametric fits).
#'
#' @return An \eqn{n x d} matrix of simulated data from the \code{kdevine}
#' object.
#'
#' @seealso
#' \code{\link{kdevine}},
#' \code{\link{rkdevinecop}},
#' \code{\link{rkde1d}}
#'
#' @examples
#' # load and plot data
#' data(wdbc)
#' \dontshow{wdbc <- wdbc[1:30, ]}
#' # estimate density
#' fit <- kdevine(wdbc[, 5:7], xmin = rep(0, 3))
#'
#' # plot simulated data
#' pairs(rkdevine(nrow(wdbc), fit))
#'
#' @importFrom VineCopula pobs RVineSim
#' @export
rkdevine <- function(n, obj, quasi = FALSE) {
    # simulate from copula
    usim <- switch(obj$copula.type,
                   "kde" = rkdevinecop(n, obj$vine, quasi = quasi),
                   "parametric" = RVineSim(n, obj$vine))
    # use quantile transformation for marginals
    sapply(seq_len(ncol(usim)), function(i)
        qkde1d(usim[, i], obj$marg.dens[[i]])
    )
}