#' Proportionality Distance

#' 

#' \code{proportionality} is a wrapper that compute proportionality distance for 

#' a clustering result (\code{pca}, \code{spca}, \code{pls}, \code{spls}, \code{block.pls}, \code{block.spls}).

#' and it performs a u-test to compare the median within a cluster to the median of the entire background set.

#' 

#' @param X an object of the class: \code{pca}, \code{spca}, \code{pls}, \code{spls}, \code{block.pls} or \code{block.spls}

#' 

#' @return 

#' Return a list containing the following components:

#'   \item{propr.distance}{Square matrix with proportionality distance between pairs of features}

#'   \item{propr.distance.w.cluster}{distance between pairs with cluster label}

#'   \item{pvalue}{Wilcoxon U-test p-value comparing the medians within clusters and with the entire background set}

#' 

#' 

#' @references 

#' Lovell, D., Pawlowsky-Glahn, V., Egozcue, J. J., Marguerat, S., Bähler, J. (2015). Proportionality: a valid alternative to correlation for relative data. PLoS Comput. Biol. 11, e1004075. doi: 10.1371/journal.pcbi.1004075

#' 

#' Quinn, T. P., Richardson, M. F., Lovell, D., Crowley, T. M. (2017). propr: an r-package for identifying proportionally abundant features using compositional data analysis. Sci. Rep. 7, 16252. doi: 10.1038/s41598-017-16520-0

#' 

#' @examples 

#' demo <- suppressWarnings(get_demo_cluster())

#' 

#' # pca

#' X <- demo$pca

#' propr.res <- proportionality(X)

#' plot(propr.res)

#' 

#' # pls

#' X <- demo$spls

#' propr.res <- proportionality(X)

#' plot(propr.res)

#' 

#' # block.pls

#' X <- demo$block.spls

#' propr.res <- proportionality(X)

#' plot(propr.res)

#' 

#' @importFrom dplyr mutate filter rename left_join

#' @importFrom magrittr %>%

#' @importFrom tibble rownames_to_column

#' @importFrom tidyr pivot_longer

#' 

#' @export

proportionality <- function(X){

    #stopifnot(is(X, c("pca", "spca", "mixo_pls", "mixo_spls", "block.pls", "block.spls")))

    # modif bioc 3.16

    stopifnot(any(class(X) %in% c("pca", "spca", "mixo_pls", "mixo_spls", "block.pls", "block.spls")))

    # 1. get cluster

    cluster.info <- getCluster(X) %>%

        dplyr::select(molecule, cluster)

    # 2. extract data and add cluster

    # pca / spca

    if(any(class(X) %in% c("pca", "spca"))){

        # unscaling + positive value

        data <- unscale(X$X) %>% `+`(abs(min(.))) 

    } else if(any(class(X) %in% c("mixo_pls", "mixo_spls"))){

        # unscale X, unscale Y, cat

        data.X <- unscale(X$X) %>%

            `+`(abs(min(.)))

        data.Y <- unscale(X$Y) %>%

            `+`(abs(min(.)))

        data <- cbind(data.X, data.Y)

    } else if(any(class(X) %in% c("block.pls", "block.spls"))){

        # if(is.null(X$Y)){  ## no need: Y is passed to X

        data <- lapply(X$X, function(x) x %>%

                           unscale %>%

                           `+`(abs(min(.)))) %>%

            do.call(what="cbind")

        # } else {

        #     data.X <- lapply(X$X, function(x) x %>% unscale %>% `+`(abs(min(.)))) %>%

        #         do.call(what="cbind")

        #     data.Y <- unscale(X$Y) %>% `+`(abs(min(.)))

        #     data <- cbind(data.X, data.Y)

        # }

    }

    # 3. compute phi_s

    # update 16 march 2023: propr deprecated;

    # implementation of phi_s only based on Quinn paper

    # data.propr <- suppressMessages(as.data.frame(propr::propr(data, metric = "phs")@matrix))

    # gather, add cluster info, compute basic stats

    data.propr <- get_phi_s(data)

    data.propr.gather <- data.propr %>%

        tibble::rownames_to_column("feature1") %>%

        tidyr::pivot_longer(-feature1, names_to = "feature2", values_to = 'value') %>%

        dplyr::filter(feature1 %in% cluster.info$molecule) %>%

        dplyr::left_join(cluster.info, by = c("feature1"="molecule")) %>%

        dplyr::rename("cluster1"="cluster") %>%

        dplyr::filter(feature2 %in% cluster.info$molecule) %>%

        dplyr::left_join(cluster.info, by = c("feature2"="molecule")) %>%

        dplyr::rename("cluster2"="cluster") %>%

        dplyr::mutate(insideout = ifelse(cluster1 == cluster2, "inside", "outside"))

    # compute stat (median, u-test pval, adj.pval)

    res.stat <- stat_median(data.propr.gather) %>%

        na.omit()

    res <- list()

    res[["propr.distance"]] <- data.propr

    res[["propr.distance.w.cluster"]] <- data.propr.gather

    res[["pvalue"]] <- res.stat

    class(res) <- "proportionality"

    return(res)

}

#' @importFrom dplyr mutate filter pull

#' @importFrom purrr set_names

#' @importFrom magrittr %>%

stat_median <- function(res.phs.X){

    i = 1

    res.pval <- matrix(ncol = 4, nrow = 4) %>%

        as.data.frame() %>%

        purrr::set_names("cluster", "median_inside", "median_outside", "Pvalue")

    for(clu in unique(res.phs.X$cluster1)){

        inside <- res.phs.X %>%

            filter(cluster1 == clu) %>% 

            dplyr::filter(cluster2==clu) %>% 

            dplyr::pull(value)

        outside <- res.phs.X %>% 

            filter(cluster1 == clu) %>% 

            dplyr::filter(cluster2!=clu) %>% 

            dplyr::pull(value)

        #ttest.pval <- t.test(inside, outside)$p.value

        utest.pval <- stats::wilcox.test(inside, outside)$p.value

        res.pval[i,] <- c(clu, round(median(inside), digits = 2),

                          round(median(outside), digits = 2), utest.pval)

        i = i+1

    }

    as.data.frame(res.pval) %>% 

        dplyr::mutate("Adj.Pvalue" = stats::p.adjust(Pvalue, method = "fdr")) %>%

        na.omit

    return(res.pval)

}

#' @export

#' @import ggplot2

#' @importFrom mixOmics color.mixo

plot.proportionality <- function(x, ...){

    ggplot2::ggplot(data = x$propr.distance.w.cluster, 

                    aes(x=as.factor(cluster1), y=value, col=insideout)) + 

        geom_boxplot() + theme_bw() + 

        xlab("Cluster ID") + 

        ylab("Proportionality distance") + 

        labs(color = "Proportionality distance") +

        scale_color_manual(values = mixOmics::color.mixo(1:2))

}

# based on the works of Thom Quinn (https://github.com/tpq/propr)

get_phi_s <- function(counts){

    ct <- counts

    if (any(as.matrix(counts) == 0)) {

        #message("Alert: Replacing 0s with next smallest value.")

        zeros <- ct == 0

        ct[zeros] <- min(ct[!zeros])

    }

    #ivar = "clr"

    #use <- propr::ivar2index(ct, ivar)

    use <- 1:ncol(counts)

    logX <- log(ct)

    logSet <- logX[, use, drop = FALSE]

    ref <- rowMeans(logSet)

    lr <- sweep(logX, 1, ref, "-")

    mat <- as.data.frame(to_lr2phs(lr))

    colnames(mat) <- colnames(lr)

    rownames(mat) <- colnames(lr)

    return(mat)

}

to_lr2phs <- function(lr){

    # Calculate phs = var(a-b)/var(a+b)

    nfeats <- ncol(lr);

    mat <- matrix(nrow = nfeats, ncol = nfeats)

    for(i in 1:nfeats){

        for(j in 1:nfeats){

            if(i == j){

                mat[i,j] <- 0

            } else {

                a <- lr[,i]; b <- lr[,j]

                mat[i,j] <- var(a-b)/var(a+b)

            }

        }

    }

    # return the same matrix as propr

    return(mat)

}