r"""Composition Statistics (:mod:`skbio.stats.composition`)
=======================================================

.. currentmodule:: skbio.stats.composition

This module provides functions for compositional data analysis.

Many omics datasets are inherently compositional -- meaning that they are best
interpreted as proportions or percentages rather than absolute counts.

Formally, sample :math:`x` is a composition if :math:`\sum_{i=0}^D x_{i} = c`
and :math:`x_{i} > 0`, :math:`1 \leq i \leq D` and :math:`c` is a real-valued
constant and there are :math:`D` components (features) for this composition.
In this module :math:`c=1`. Compositional data can be analyzed using
**Aitchison geometry** [1]_.

However, in this framework, standard real Euclidean operations such as addition
and multiplication no longer apply. Only operations such as perturbation and
power can be used to manipulate this data.

This module allows two styles of manipulation of compositional data.
Compositional data can be analyzed using perturbation and power operations,
which can be useful for simulation studies. The alternative strategy is to
transform compositional data into the real space. Right now, the centre log
ratio transform (clr) and the isometric log ratio transform (ilr) [2]_ can be
used to accomplish this. This transform can be useful for performing standard
statistical methods such as parametric hypothesis testing, regression and more.

The major caveat of using this framework is dealing with zeros. In Aitchison
geometry, only compositions with non-zero components can be considered.
The multiplicative replacement technique [3]_ can be used to substitute these
zeros with small pseudocounts without introducing major distortions to the
data.

Functions
---------

.. autosummary::
   :toctree:

   closure
   multi_replace
   multiplicative_replacement
   perturb
   perturb_inv
   power
   inner
   clr
   clr_inv
   ilr
   ilr_inv
   alr
   alr_inv
   centralize
   vlr
   pairwise_vlr
   tree_basis
   ancom
   sbp_basis
   dirmult_ttest

References
----------
.. [1] V. Pawlowsky-Glahn, J. J. Egozcue, R. Tolosana-Delgado (2015),
   Modeling and Analysis of Compositional Data, Wiley, Chichester, UK

.. [2] J. J. Egozcue.,  "Isometric Logratio Transformations for
   Compositional Data Analysis" Mathematical Geology, 35.3 (2003)

.. [3] J. A. Martin-Fernandez,  "Dealing With Zeros and Missing Values in
   Compositional Data Sets Using Nonparametric Imputation",
   Mathematical Geology, 35.3 (2003)

"""  # noqa: D205, D415

# ----------------------------------------------------------------------------
# Copyright (c) 2013--, scikit-bio development team.
#
# Distributed under the terms of the Modified BSD License.
#
# The full license is in the file LICENSE.txt, distributed with this software.
# ----------------------------------------------------------------------------

from warnings import warn, simplefilter

import numpy as np
import pandas as pd
import scipy.stats
from scipy.sparse import coo_matrix
from scipy.stats import t, gmean
from statsmodels.stats.weightstats import CompareMeans

from skbio.stats.distance import DistanceMatrix
from skbio.util import find_duplicates
from skbio.util._misc import get_rng
from skbio.util._warning import _warn_deprecated
from statsmodels.stats.multitest import multipletests as sm_multipletests


def closure(mat):
    """Perform closure to ensure that all elements add up to 1.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        The matrix where all of the values are non-zero and each composition
        (row) adds up to 1.

    Raises
    ------
    ValueError
        If any values are negative.
    ValueError
        If the matrix has more than two dimensions.
    ValueError
        If there is a row that has all zeros.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import closure
    >>> X = np.array([[2, 2, 6], [4, 4, 2]])
    >>> closure(X)
    array([[ 0.2,  0.2,  0.6],
           [ 0.4,  0.4,  0.2]])

    """
    mat = np.atleast_2d(mat)
    if np.any(mat < 0):
        raise ValueError("Cannot have negative proportions")
    if mat.ndim > 2:
        raise ValueError("Input matrix can only have two dimensions or less")
    if np.all(mat == 0, axis=1).sum() > 0:
        raise ValueError("Input matrix cannot have rows with all zeros")
    mat = mat / mat.sum(axis=1, keepdims=True)
    return mat.squeeze()


def multi_replace(mat, delta=None):
    r"""Replace all zeros with small non-zero values.

    It uses the multiplicative replacement strategy [1]_, replacing zeros with
    a small positive :math:`\delta` and ensuring that the compositions still
    add up to 1.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    delta : float, optional
        A small number to be used to replace zeros. If not specified, the
        default value is :math:`\delta = \frac{1}{N^2}` where :math:`N` is the
        number of components.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        The matrix where all of the values are non-zero and each composition
        (row) adds up to 1.

    Raises
    ------
    ValueError
        If negative proportions are created due to a large ``delta``.

    Notes
    -----
    This method will result in negative proportions if a large delta is chosen.

    References
    ----------
    .. [1] J. A. Martin-Fernandez. "Dealing With Zeros and Missing Values in
           Compositional Data Sets Using Nonparametric Imputation"

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import multi_replace
    >>> X = np.array([[.2, .4, .4, 0],[0, .5, .5, 0]])
    >>> multi_replace(X)
    array([[ 0.1875,  0.375 ,  0.375 ,  0.0625],
           [ 0.0625,  0.4375,  0.4375,  0.0625]])

    """
    mat = closure(mat)
    z_mat = mat == 0

    num_feats = mat.shape[-1]
    tot = z_mat.sum(axis=-1, keepdims=True)

    if delta is None:
        delta = (1.0 / num_feats) ** 2

    zcnts = 1 - tot * delta
    if np.any(zcnts) < 0:
        raise ValueError(
            "The multiplicative replacement created negative "
            "proportions. Consider using a smaller `delta`."
        )
    mat = np.where(z_mat, delta, zcnts * mat)
    return mat.squeeze()


def multiplicative_replacement(mat, delta=None):
    r"""Replace all zeros with small non-zero values.

    This function is an alias for ``multi_replace``.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    delta : float, optional
        A small number to be used to replace zeros. If not specified, the
        default value is :math:`\delta = \frac{1}{N^2}` where :math:`N` is the
        number of components.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        The matrix where all of the values are non-zero and each composition
        (row) adds up to 1.

    Raises
    ------
    ValueError
        If negative proportions are created due to a large ``delta``.

    Warnings
    --------
    ``multiplicative_replacement`` is deprecated as of ``0.6.0`` in favor of
    ``multi_replace``.

    See Also
    --------
    multi_replace

    """
    _warn_deprecated(multiplicative_replacement, "0.6.0")
    return multi_replace(mat, delta)


def perturb(x, y):
    r"""Perform the perturbation operation.

    This operation is defined as:

    .. math::
        x \oplus y = C[x_1 y_1, \ldots, x_D y_D]

    :math:`C[x]` is the closure operation defined as:

    .. math::
        C[x] = \left[\frac{x_1}{\sum_{i=1}^{D} x_i},\ldots,
                     \frac{x_D}{\sum_{i=1}^{D} x_i} \right]

    for some :math:`D` dimensional real vector :math:`x` and
    :math:`D` is the number of components for every composition.

    Parameters
    ----------
    x : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    y : array_like of shape (n_compositions, n_components)
        A matrix of proportions.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
       A matrix of proportions where all of the values are non-zero and each
       composition (row) adds up to 1.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import perturb

    Consider a very simple environment with only three species. The species in
    the environment are evenly distributed and their proportions are equal:

    >>> before = np.array([1/3, 1/3, 1/3])

    Suppose that an antibiotic kills off half of the population for the first
    two species, but doesn't harm the third species. Then the perturbation
    vector would be as follows:

    >>> after = np.array([1/2, 1/2, 1])

    And the resulting perturbation would be:

    >>> perturb(before, after)
    array([ 0.25,  0.25,  0.5 ])

    """
    x, y = closure(x), closure(y)
    return closure(x * y)


def perturb_inv(x, y):
    r"""Perform the inverse perturbation operation.

    This operation is defined as:

    .. math::
        x \ominus y = C[x_1 y_1^{-1}, \ldots, x_D y_D^{-1}]

    :math:`C[x]` is the closure operation defined as:

    .. math::
        C[x] = \left[\frac{x_1}{\sum_{i=1}^{D} x_i},\ldots,
                     \frac{x_D}{\sum_{i=1}^{D} x_i} \right]

    for some :math:`D` dimensional real vector :math:`x` and :math:`D` is the
    number of components for every composition.

    Parameters
    ----------
    x : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    y : array_like of shape (n_compositions, n_components)
        A matrix of proportions.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        A matrix of proportions where all of the values are non-zero and each
        composition (row) adds up to 1.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import perturb_inv
    >>> x = np.array([.1, .3, .4, .2])
    >>> y = np.array([1/6, 1/6, 1/3, 1/3])
    >>> perturb_inv(x, y)
    array([ 0.14285714,  0.42857143,  0.28571429,  0.14285714])

    """
    x, y = closure(x), closure(y)
    return closure(x / y)


def power(x, a):
    r"""Perform the power operation.

    This operation is defined as follows:

    .. math::
        `x \odot a = C[x_1^a, \ldots, x_D^a]

    :math:`C[x]` is the closure operation defined as:

    .. math::
        C[x] = \left[\frac{x_1}{\sum_{i=1}^{D} x_i},\ldots,
                     \frac{x_D}{\sum_{i=1}^{D} x_i} \right]

    for some :math:`D` dimensional real vector :math:`x` and
    :math:`D` is the number of components for every composition.

    Parameters
    ----------
    x : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    a : float
        A scalar exponent.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
       The matrix where all of the values are non-zero and each composition
       (row) adds up to 1.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import power
    >>> x = np.array([.1, .3, .4, .2])
    >>> power(x, .1)
    array([ 0.23059566,  0.25737316,  0.26488486,  0.24714631])

    """
    x = closure(x)
    return closure(x**a).squeeze()


def inner(x, y):
    r"""Calculate the Aitchson inner product.

    This inner product is defined as follows:

    .. math::
        \langle x, y \rangle_a =
        \frac{1}{2D} \sum\limits_{i=1}^{D} \sum\limits_{j=1}^{D}
        \ln\left(\frac{x_i}{x_j}\right) \ln\left(\frac{y_i}{y_j}\right)

    Parameters
    ----------
    x : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    y : array_like of shape (n_compositions, n_components)
        A matrix of proportions.

    Returns
    -------
    ndarray or scalar of shape (n_compositions, n_compositions)
        Inner product result.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import inner
    >>> x = np.array([.1, .3, .4, .2])
    >>> y = np.array([.2, .4, .2, .2])
    >>> inner(x, y)  # doctest: +ELLIPSIS
    0.2107852473...

    """
    x = closure(x)
    y = closure(y)
    a, b = clr(x), clr(y)
    return a.dot(b.T)


def clr(mat):
    r"""Perform centre log ratio transformation.

    This function transforms compositions from Aitchison geometry to the real
    space. The :math:`clr` transform is both an isometry and an isomorphism
    defined on the following spaces:

    .. math::
        clr: S^D \rightarrow U

    where :math:`U=
    \{x :\sum\limits_{i=1}^D x = 0 \; \forall x \in \mathbb{R}^D\}`

    It is defined for a composition :math:`x` as follows:

    .. math::
        clr(x) = \ln\left[\frac{x_1}{g_m(x)}, \ldots, \frac{x_D}{g_m(x)}\right]

    where :math:`g_m(x) = (\prod\limits_{i=1}^{D} x_i)^{1/D}` is the geometric
    mean of :math:`x`.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        Clr-transformed matrix.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import clr
    >>> x = np.array([.1, .3, .4, .2])
    >>> clr(x)
    array([-0.79451346,  0.30409883,  0.5917809 , -0.10136628])

    """
    mat = closure(mat)
    lmat = np.log(mat)
    gm = lmat.mean(axis=-1, keepdims=True)
    return (lmat - gm).squeeze()


def clr_inv(mat):
    r"""Perform inverse centre log ratio transformation.

    This function transforms compositions from the real space to Aitchison
    geometry. The :math:`clr^{-1}` transform is both an isometry, and an
    isomorphism defined on the following spaces:

    .. math::
        clr^{-1}: U \rightarrow S^D

    where :math:`U=
    \{x :\sum\limits_{i=1}^D x = 0 \; \forall x \in \mathbb{R}^D\}`

    This transformation is defined as follows:

    .. math::
        clr^{-1}(x) = C[\exp( x_1, \ldots, x_D)]

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of clr-transformed data.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        Inverse clr-transformed matrix.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import clr_inv
    >>> x = np.array([.1, .3, .4, .2])
    >>> clr_inv(x)
    array([ 0.21383822,  0.26118259,  0.28865141,  0.23632778])

    """
    # for numerical stability (aka softmax trick)
    mat = np.atleast_2d(mat)
    emat = np.exp(mat - mat.max(axis=-1, keepdims=True))
    return closure(emat)


def ilr(mat, basis=None, check=True):
    r"""Perform isometric log ratio transformation.

    This function transforms compositions from Aitchison simplex to the real
    space. The :math:`ilr` transform is both an isometry, and an isomorphism
    defined on the following spaces:

    .. math::
        ilr: S^D \rightarrow \mathbb{R}^{D-1}

    The ilr transformation is defined as follows:

    .. math::
        ilr(x) =
        [\langle x, e_1 \rangle_a, \ldots, \langle x, e_{D-1} \rangle_a]

    where :math:`[e_1,\ldots,e_{D-1}]` is an orthonormal basis in the simplex.

    If an orthornormal basis isn't specified, the J. J. Egozcue orthonormal
    basis derived from Gram-Schmidt orthogonalization will be used by default.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    basis : ndarray or sparse matrix, optional
        Orthonormal basis for Aitchison simplex. Defaults to J. J. Egozcue
        orthonormal basis.
    check : bool
        Check to see if basis is orthonormal.

    Returns
    -------
    ndarray of shape (n_compositions, n_components - 1)
        Ilr-transformed matrix.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import ilr
    >>> x = np.array([.1, .3, .4, .2])
    >>> ilr(x)
    array([-0.7768362 , -0.68339802,  0.11704769])

    Notes
    -----
    If the ``basis`` parameter is specified, it is expected to be a basis in
    the Aitchison simplex. If there are :math:`D - 1` elements specified in
    ``mat``, then the dimensions of the basis needs be :math:`(D-1) \times D`,
    where rows represent basis vectors, and the columns represent proportions.

    """
    mat = closure(mat)
    if basis is None:
        d = mat.shape[-1]
        basis = _gram_schmidt_basis(d)  # dimension (d-1) x d
    else:
        if len(basis.shape) != 2:
            raise ValueError(
                "Basis needs to be a 2D matrix, "
                "not a %dD matrix." % (len(basis.shape))
            )
        if check:
            _check_orthogonality(basis)

    return clr(mat) @ basis.T


def ilr_inv(mat, basis=None, check=True):
    r"""Perform inverse isometric log ratio transform.

    This function transforms compositions from the real space to Aitchison
    geometry. The :math:`ilr^{-1}` transform is both an isometry, and an
    isomorphism defined on the following spaces:

    .. math::
        ilr^{-1}: \mathbb{R}^{D-1} \rightarrow S^D

    The inverse ilr transformation is defined as follows:

    .. math::
        ilr^{-1}(x) = \bigoplus\limits_{i=1}^{D-1} x \odot e_i

    where :math:`[e_1,\ldots, e_{D-1}]` is an orthonormal basis in the simplex.

    If an orthonormal basis isn't specified, the J. J. Egozcue orthonormal
    basis derived from Gram-Schmidt orthogonalization will be used by
    default.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components - 1)
        A matrix of ilr-transformed data.
    basis : ndarray or sparse matrix, optional
        Orthonormal basis for Aitchison simplex. Defaults to J. J. Egozcue
        orthonormal basis.
    check : bool
        Check to see if basis is orthonormal.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        Inverse ilr-transformed matrix.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import ilr
    >>> x = np.array([.1, .3, .6,])
    >>> ilr_inv(x)
    array([ 0.34180297,  0.29672718,  0.22054469,  0.14092516])

    Notes
    -----
    If the ``basis`` parameter is specified, it is expected to be a basis in
    the Aitchison simplex. If there are :math:`D - 1` elements specified in
    ``mat``, then the dimensions of the basis needs be :math:`(D-1) \times D`,
    where rows represent basis vectors, and the columns represent proportions.

    """
    mat = np.atleast_2d(mat)
    if basis is None:
        # dimension d-1 x d basis
        basis = _gram_schmidt_basis(mat.shape[-1] + 1)
    else:
        if len(basis.shape) != 2:
            raise ValueError(
                "Basis needs to be a 2D matrix, "
                "not a %dD matrix." % (len(basis.shape))
            )
        if check:
            _check_orthogonality(basis)
        # this is necessary, since the clr function
        # performs np.squeeze()
        basis = np.atleast_2d(basis)

    return clr_inv(mat @ basis)


def alr(mat, denominator_idx=0):
    r"""Perform additive log ratio transformation.

    This function transforms compositions from a D-part Aitchison simplex to
    a non-isometric real space of D-1 dimensions. The argument
    ``denominator_col`` defines the index of the column used as the common
    denominator. The :math:`alr` transformed data are amenable to multivariate
    analysis as long as statistics don't involve distances.

    .. math::
        alr: S^D \rightarrow \mathbb{R}^{D-1}

    The alr transformation is defined as follows

    .. math::
        alr(x) = \left[ \ln \frac{x_1}{x_D}, \ldots,
        \ln \frac{x_{D-1}}{x_D} \right]

    where :math:`D` is the index of the part used as common denominator.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    denominator_idx : int
        The index of the column (2-D matrix) or position (vector) of ``mat``
        which should be used as the reference composition. Default is 0 which
        specifies the first column or position.

    Returns
    -------
    ndarray of shape (n_compositions, n_components - 1)
        Alr-transformed data projected in a non-isometric real space of
        :math:`D - 1` dimensions for a *D*-parts composition.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import alr
    >>> x = np.array([.1, .3, .4, .2])
    >>> alr(x)
    array([ 1.09861229,  1.38629436,  0.69314718])

    """
    mat = closure(mat)
    if mat.ndim == 2:
        mat_t = mat.T
        numerator_idx = list(range(0, mat_t.shape[0]))
        del numerator_idx[denominator_idx]
        lr = np.log(mat_t[numerator_idx, :] / mat_t[denominator_idx, :]).T
    elif mat.ndim == 1:
        numerator_idx = list(range(0, mat.shape[0]))
        del numerator_idx[denominator_idx]
        lr = np.log(mat[numerator_idx] / mat[denominator_idx])
    else:
        raise ValueError("mat must be either 1D or 2D")
    return lr


def alr_inv(mat, denominator_idx=0):
    r"""Perform inverse additive log ratio transform.

    This function transforms compositions from the non-isometric real space of
    alrs to Aitchison geometry.

    .. math::
        alr^{-1}: \mathbb{R}^{D-1} \rightarrow S^D

    The inverse alr transformation is defined as follows:

    .. math::
         alr^{-1}(x) = C[exp([y_1, y_2, ..., y_{D-1}, 0])]

    where :math:`C[x]` is the closure operation defined as

    .. math::
        C[x] = \left[\frac{x_1}{\sum_{i=1}^{D} x_i},\ldots,
                     \frac{x_D}{\sum_{i=1}^{D} x_i} \right]

    for some :math:`D` dimensional real vector :math:`x` and
    :math:`D` is the number of components for every composition.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components - 1)
        A matrix of alr-transformed data.
    denominator_idx : int
        The index of the column (2-D matrix) or position (vector) of ``mat``
        which should be used as the reference composition. Default is 0 which
        specifies the first column or position.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        Inverse alr-transformed matrix or vector where rows sum to 1.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import alr, alr_inv
    >>> x = np.array([.1, .3, .4, .2])
    >>> alr_inv(alr(x))
    array([ 0.1,  0.3,  0.4,  0.2])

    """
    mat = np.array(mat)
    if mat.ndim == 2:
        mat_idx = np.insert(mat, denominator_idx, np.repeat(0, mat.shape[0]), axis=1)
        comp = np.zeros(mat_idx.shape)
        comp[:, denominator_idx] = 1 / (np.exp(mat).sum(axis=1) + 1)
        numerator_idx = list(range(0, comp.shape[1]))
        del numerator_idx[denominator_idx]
        for i in numerator_idx:
            comp[:, i] = comp[:, denominator_idx] * np.exp(mat_idx[:, i])
    elif mat.ndim == 1:
        mat_idx = np.insert(mat, denominator_idx, 0, axis=0)
        comp = np.zeros(mat_idx.shape)
        comp[denominator_idx] = 1 / (np.exp(mat).sum(axis=0) + 1)
        numerator_idx = list(range(0, comp.shape[0]))
        del numerator_idx[denominator_idx]
        for i in numerator_idx:
            comp[i] = comp[denominator_idx] * np.exp(mat_idx[i])
    else:
        raise ValueError("mat must be either 1D or 2D")
    return comp


def centralize(mat):
    r"""Center data around its geometric average.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.

    Returns
    -------
    ndarray of shape (n_compositions, n_components)
        Centered composition matrix.

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import centralize
    >>> X = np.array([[.1, .3, .4, .2], [.2, .2, .2, .4]])
    >>> centralize(X)
    array([[ 0.17445763,  0.30216948,  0.34891526,  0.17445763],
           [ 0.32495488,  0.18761279,  0.16247744,  0.32495488]])

    """
    mat = closure(mat)
    cen = gmean(mat, axis=0)
    return perturb_inv(mat, cen)


def _vlr(x, y, ddof):
    r"""Calculate variance log ratio.

    Parameters
    ----------
    x : array_like of shape (n_components,)
        A vector of proportions.
    y : array_like of shape (n_components,)
        A vector of proportions.
    ddof : int
        Degrees of freedom.

    Returns
    -------
    float
        Variance log ratio value.

    """
    # Log transformation
    x = np.log(x)
    y = np.log(y)

    # Variance log ratio
    return np.var(x - y, ddof=ddof)


def _robust_vlr(x, y, ddof):
    r"""Calculate variance log ratio while masking zeros.

    Parameters
    ----------
    x : array_like of shape (n_components,)
        A vector of proportions.
    y : array_like of shape (n_components,)
        A vector of proportions.
    ddof : int
        Degrees of freedom.

    Returns
    -------
    float
        Variance log ratio value.

    """
    # Mask zeros
    x = np.ma.masked_array(x, mask=x == 0)
    y = np.ma.masked_array(y, mask=y == 0)

    # Log transformation
    x = np.ma.log(x)
    y = np.ma.log(y)

    # Variance log ratio
    return np.ma.var(x - y, ddof=ddof)


def vlr(x, y, ddof=1, robust=False):
    r"""Calculate variance log ratio.

    Parameters
    ----------
    x : array_like of shape (n_components,)
        A vector of proportions.
    y : array_like of shape (n_components,)
        A vector of proportions.
    ddof : int
        Degrees of freedom.
    robust : bool
        Whether to mask zeros at the cost of performance.

    Returns
    -------
    float
        Variance log ratio value.

    Notes
    -----
    Variance log ratio was described in [1]_ and [2]_.

    References
    ----------
    .. [1] V. Lovell D, Pawlowsky-Glahn V, Egozcue JJ, Marguerat S,
           Bähler J (2015) Proportionality: A Valid Alternative to
           Correlation for Relative Data. PLoS Comput Biol 11(3): e1004075.
           https://doi.org/10.1371/journal.pcbi.1004075
    .. [2] Erb, I., Notredame, C.
           How should we measure proportionality on relative gene
           expression data?. Theory Biosci. 135, 21-36 (2016).
           https://doi.org/10.1007/s12064-015-0220-8

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import vlr
    >>> x = np.exp([1, 2, 3])
    >>> y = np.exp([2, 3, 4])
    >>> vlr(x, y)  # no zeros
    0.0

    """
    # Convert array_like to numpy array
    x = closure(x)
    y = closure(y)

    # Set up input and parameters
    kwargs = {
        "x": x,
        "y": y,
        "ddof": ddof,
    }

    # Run backend function
    if robust:
        return _robust_vlr(**kwargs)
    else:
        return _vlr(**kwargs)


def _pairwise_vlr(mat, ddof):
    r"""Perform pairwise variance log ratio transformation.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    ddof : int
        Degrees of freedom.

    Returns
    -------
    ndarray of shape (n_compositions, n_compositions)
        Distance matrix of variance log ratio values.

    """
    # Log Transform
    X_log = np.log(mat)

    # Variance Log Ratio
    covariance = np.cov(X_log.T, ddof=ddof)
    diagonal = np.diagonal(covariance)
    vlr_data = -2 * covariance + diagonal[:, np.newaxis] + diagonal
    return vlr_data


def pairwise_vlr(mat, ids=None, ddof=1, robust=False, validate=True):
    r"""Perform pairwise variance log ratio transformation.

    Parameters
    ----------
    mat : array_like of shape (n_compositions, n_components)
        A matrix of proportions.
    ids : array_like of str of shape (n_components,)
        Component names.
    ddof : int
        Degrees of freedom.
    robust : bool
        Whether to mask zeros at the cost of performance.
    validate : bool
        Whether to validate the distance matrix after construction.

    Returns
    -------
    skbio.DistanceMatrix if validate=True
        Distance matrix of variance log ratio values.
    skbio.DissimilarityMatrix if validate=False
        Dissimilarity matrix of variance log ratio values.

    Notes
    -----
    Pairwise variance log ratio transformation was described in [1]_ and [2]_.

    References
    ----------
    .. [1] V. Lovell D, Pawlowsky-Glahn V, Egozcue JJ, Marguerat S,
           Bähler J (2015) Proportionality: A Valid Alternative to
           Correlation for Relative Data. PLoS Comput Biol 11(3): e1004075.
           https://doi.org/10.1371/journal.pcbi.1004075
    .. [2] Erb, I., Notredame, C.
           How should we measure proportionality on relative gene
           expression data?. Theory Biosci. 135, 21-36 (2016).
           https://doi.org/10.1007/s12064-015-0220-8

    Examples
    --------
    >>> import numpy as np
    >>> from skbio.stats.composition import pairwise_vlr
    >>> mat = np.array([np.exp([1, 2, 2]),
    ...                 np.exp([2, 3, 6]),
    ...                 np.exp([2, 3, 12])]).T
    >>> dism = pairwise_vlr(mat)
    >>> dism.redundant_form()
    array([[  0.,   3.,  27.],
           [  3.,   0.,  12.],
           [ 27.,  12.,   0.]])

    """
    # Mask zeros
    mat = closure(mat.astype(np.float64))

    # Set up input and parameters
    kwargs = {
        "mat": mat,
        "ddof": ddof,
    }

    # Variance log ratio
    if robust:
        raise NotImplementedError("Pairwise version of robust VLR not implemented.")
    else:
        vlr_data = _pairwise_vlr(**kwargs)

    # Return distance matrix
    if validate:
        vlr_data = 0.5 * (vlr_data + vlr_data.T)
        return DistanceMatrix(vlr_data, ids=ids)

    # Return dissimilarity matrix
    else:
        return DistanceMatrix(vlr_data, ids=ids, validate=False)


def tree_basis(tree):
    r"""Calculate the sparse representation of an ilr basis from a tree.

    This computes an orthonormal basis specified from a bifurcating tree.

    Parameters
    ----------
    tree : skbio.TreeNode
        Input bifurcating tree. Must be strictly bifurcating (i.e. every
        internal node needs to have exactly two children). This is used to
        specify the ilr basis.

    Returns
    -------
    scipy.sparse.coo_matrix
        The ilr basis required to perform the ilr_inv transform. This is also
        known as the sequential binary partition. Note that this matrix is
        represented in clr coordinates.
    list of str
        List of tree node names indicating the ordering in the basis.

    Raises
    ------
    ValueError
        If the tree doesn't contain two branches.

    Examples
    --------
    >>> from skbio import TreeNode
    >>> tree = u"((b,c)a, d)root;"
    >>> t = TreeNode.read([tree])
    >>> basis, nodes = tree_basis(t)
    >>> basis.toarray()
    array([[-0.40824829, -0.40824829,  0.81649658],
           [-0.70710678,  0.70710678,  0.        ]])

    """
    # Specifies which child is numerator and denominator
    # within any given node in a tree.
    NUMERATOR = 1
    DENOMINATOR = 0

    # this is inspired by @wasade in
    # https://github.com/biocore/gneiss/pull/8
    t = tree.copy()
    D = len(list(tree.tips()))

    # calculate number of tips under each node
    for n in t.postorder(include_self=True):
        if n.is_tip():
            n._tip_count = 1
        else:
            if len(n.children) == 2:
                left, right = (
                    n.children[NUMERATOR],
                    n.children[DENOMINATOR],
                )
            else:
                raise ValueError("Not a strictly bifurcating tree.")
            n._tip_count = left._tip_count + right._tip_count

    # calculate k, r, s, t coordinate for each node
    left, right = (
        t.children[NUMERATOR],
        t.children[DENOMINATOR],
    )
    t._k, t._r, t._s, t._t = 0, left._tip_count, right._tip_count, 0
    for n in t.preorder(include_self=False):
        if n.is_tip():
            n._k, n._r, n._s, n._t = 0, 0, 0, 0

        elif n == n.parent.children[NUMERATOR]:
            n._k = n.parent._k
            n._r = n.children[NUMERATOR]._tip_count
            n._s = n.children[DENOMINATOR]._tip_count
            n._t = n.parent._s + n.parent._t
        elif n == n.parent.children[DENOMINATOR]:
            n._k = n.parent._r + n.parent._k
            n._r = n.children[NUMERATOR]._tip_count
            n._s = n.children[DENOMINATOR]._tip_count
            n._t = n.parent._t
        else:
            raise ValueError("Tree topology is not correct.")

    # navigate through tree to build the basis in a sparse matrix form
    value = []
    row, col = [], []
    nodes = []
    i = 0

    for n in t.levelorder(include_self=True):
        if n.is_tip():
            continue

        for j in range(n._k, n._k + n._r):
            row.append(i)
            # consider tips in reverse order. May want to rethink
            # this orientation in the future.
            col.append(D - 1 - j)
            A = np.sqrt(n._s / (n._r * (n._s + n._r)))

            value.append(A)

        for j in range(n._k + n._r, n._k + n._r + n._s):
            row.append(i)
            col.append(D - 1 - j)
            B = -np.sqrt(n._r / (n._s * (n._s + n._r)))

            value.append(B)
        i += 1
        nodes.append(n.name)

    basis = coo_matrix((value, (row, col)), shape=(D - 1, D))

    return basis, nodes


def _calc_p_adjust(name, p):
    """
    Calculate the p-value adjustment for a given method.

    Parameters
    -------
        name : str
            The name of the *p*-value correction function.
            This should match one of the method names available
            in `statsmodels.stats.multitest.multipletests`.
        p : ndarray of shape (n_tests,)
            Original *p*-values.

    Returns
    -------
        p : ndarray of shape (n_tests,)
            Corrected *p*-values.

    Raises
    -------
        ValueError: If the given method name is not available.

    See Also
    --------
    statsmodels.stats.multitest.multipletests

    References
    ----------
    .. [1] https://www.statsmodels.org/dev/generated/statsmodels.stats.multitest.multipletests.html

    """
    name_ = name.lower()

    # Original options are kept for backwards compatibility
    if name_ in ("holm", "holm-bonferroni"):
        name_ = "holm"
    if name_ in ("bh", "fdr_bh", "benjamini-hochberg"):
        name_ = "fdr_bh"

    try:
        res = sm_multipletests(pvals=p, alpha=0.05, method=name_)
    except ValueError as e:
        if "method not recognized" in str(e):
            raise ValueError(f"{name} is not an available FDR correction method.")
        else:
            raise ValueError(f"Cannot perform FDR correction using the {name} method.")
    else:
        return res[1]


def ancom(
    table,
    grouping,
    alpha=0.05,
    tau=0.02,
    theta=0.1,
    p_adjust="holm",
    significance_test="f_oneway",
    percentiles=(0.0, 25.0, 50.0, 75.0, 100.0),
    multiple_comparisons_correction="holm-bonferroni",
):
    r"""Perform a differential abundance test using ANCOM.

    Analysis of composition of microbiomes (ANCOM) is done by calculating
    pairwise log ratios between all features and performing a significance
    test to determine if there is a significant difference in feature ratios
    with respect to the variable of interest.

    In an experiment with only two treatments, this tests the following
    hypothesis for feature :math:`i`:

    .. math::

        H_{0i}: \mathbb{E}[\ln(u_i^{(1)})] = \mathbb{E}[\ln(u_i^{(2)})]

    where :math:`u_i^{(1)}` is the mean abundance for feature :math:`i` in the
    first group and :math:`u_i^{(2)}` is the mean abundance for feature
    :math:`i` in the second group.

    Parameters
    ----------
    table : pd.DataFrame
        A 2-D matrix of strictly positive values (i.e. counts or proportions)
        where the rows correspond to samples and the columns correspond to
        features.
    grouping : pd.Series
        Vector indicating the assignment of samples to groups. For example,
        these could be strings or integers denoting which group a sample
        belongs to. It must be the same length as the samples in `table`.
        The index must be the same on `table` and `grouping` but need not be
        in the same order.
    alpha : float, optional
        Significance level for each of the statistical tests. This can can be
        anywhere between 0 and 1 exclusive.
    tau : float, optional
        A constant used to determine an appropriate cutoff. A value close to
        zero indicates a conservative cutoff. This can can be anywhere between
        0 and 1 exclusive.
    theta : float, optional
        Lower bound for the proportion for the *W*-statistic. If all *W*-
        statistics are lower than theta, then no features will be detected to
        be significantly different. This can can be anywhere between 0 and 1
        exclusive.
    p_adjust : str or None, optional
        Method to correct *p*-values for multiple comparisons. Options are Holm-
        Boniferroni ("holm" or "holm-bonferroni") (default), Benjamini-
        Hochberg ("bh", "fdr_bh" or "benjamini-hochberg"), or any method supported
        by statsmodels' ``multipletests`` function. Case-insensitive. If None, no
        correction will be performed.

        .. versionchanged:: 0.6.0

            Replaces ``multiple_comparisons_correction`` for conciseness.

    significance_test : str or callable, optional
        A function to test for significance between classes. It must be able to
        accept at least two vectors of floats and returns a test statistic and
        a *p*-value. Functions under ``scipy.stats`` can be directly specified
        by name. The default is one-way ANOVA ("f_oneway").

        .. versionchanged:: 0.6.0

            Accepts test names in addition to functions.

    percentiles : iterable of floats, optional
        Percentile abundances to return for each feature in each group. By
        default, will return the minimum, 25th percentile, median, 75th
        percentile, and maximum abundances for each feature in each group.
    multiple_comparisons_correction : str or None, optional
        Alias for ``p_adjust``. For backward compatibility. Deprecated.

    Returns
    -------
    pd.DataFrame
        A table of features, their *W*-statistics and whether the null
        hypothesis is rejected.

        - ``W``: *W*-statistic, or the number of features that the current
          feature is tested to be significantly different against.

        - ``Reject null hypothesis``: Whether the feature is differentially
          abundant across groups (``True``) or not (``False``).

    pd.DataFrame
        A table of features and their percentile abundances in each group. If
        ``percentiles`` is empty, this will be an empty ``pd.DataFrame``. The
        rows in this object will be features, and the columns will be a
        multi-index where the first index is the percentile, and the second
        index is the group.

    See Also
    --------
    multi_replace
    scipy.stats.ttest_ind
    scipy.stats.f_oneway
    scipy.stats.wilcoxon
    scipy.stats.kruskal

    Warnings
    --------
    ``multiple_comparisons_correction`` is deprecated as of ``0.6.0``. It has
    been renamed to ``p_adjust``.

    ``significance_test=None`` is deprecated as of ``0.6.0``. The default value
    is now "f_oneway".

    Notes
    -----
    The developers of ANCOM recommend the following significance tests ([1]_,
    Supplementary File 1, top of page 11):

    - If there are two groups, use the standard parametric *t*-test
      (``ttest_ind``) or the non-parametric Mann-Whitney rank test
      (``mannwhitneyu``).

    - For paired samples, use the parametric paired *t*-test (``ttest_rel``) or
      the non-parametric Wilcoxon signed-rank test (``wilcoxon``).

    - If there are more than two groups, use the parametric one-way ANOVA
      (``f_oneway``) or the non-parametric Kruskal-Wallis test (``kruskal``).

    - If there are multiple measurements obtained from the individuals, use a
      Friedman test (``friedmanchisquare``).

    Because one-way ANOVA is equivalent to the standard *t*-test when the
    number of groups is two, we default to ``f_oneway`` here, which can be used
    when there are two or more groups.

    Users should refer to the documentation of these tests in SciPy to
    understand the assumptions made by each test.

    This method cannot handle any zero counts as input, since the logarithm
    of zero cannot be computed.  While this is an unsolved problem, many
    studies, including [1]_, have shown promising results by adding
    pseudocounts to all values in the matrix. In [1]_, a pseudocount of 0.001
    was used, though the authors note that a pseudocount of 1.0 may also be
    useful. Zero counts can also be addressed using the ``multi_replace`` method.

    References
    ----------
    .. [1] Mandal et al. "Analysis of composition of microbiomes: a novel
       method for studying microbial composition", Microbial Ecology in Health
       & Disease, (2015), 26.

    Examples
    --------
    >>> from skbio.stats.composition import ancom
    >>> import pandas as pd

    Let's load in a DataFrame with six samples and seven features (e.g., these
    may be bacterial taxa):

    >>> table = pd.DataFrame([[12, 11, 10, 10, 10, 10, 10],
    ...                       [9,  11, 12, 10, 10, 10, 10],
    ...                       [1,  11, 10, 11, 10, 5,  9],
    ...                       [22, 21, 9,  10, 10, 10, 10],
    ...                       [20, 22, 10, 10, 13, 10, 10],
    ...                       [23, 21, 14, 10, 10, 10, 10]],
    ...                      index=['s1', 's2', 's3', 's4', 's5', 's6'],
    ...                      columns=['b1', 'b2', 'b3', 'b4', 'b5', 'b6',
    ...                               'b7'])

    Then create a grouping vector. In this example, there is a treatment group
    and a placebo group.

    >>> grouping = pd.Series(['treatment', 'treatment', 'treatment',
    ...                       'placebo', 'placebo', 'placebo'],
    ...                      index=['s1', 's2', 's3', 's4', 's5', 's6'])

    Now run ``ancom`` to determine if there are any features that are
    significantly different in abundance between the treatment and the placebo
    groups. The first DataFrame that is returned contains the ANCOM test
    results, and the second contains the percentile abundance data for each
    feature in each group.

    >>> ancom_df, percentile_df = ancom(table, grouping)
    >>> ancom_df['W'] # doctest: +ELLIPSIS
    b1    0
    b2    4
    b3    0
    b4    1
    b5    1
    b6    0
    b7    1
    Name: W, dtype: ...

    The *W*-statistic is the number of features that a single feature is tested
    to be significantly different against. In this scenario, ``b2`` was
    detected to have significantly different abundances compared to four of the
    other features. To summarize the results from the *W*-statistic, let's take
    a look at the results from the hypothesis test. The ``Reject null
    hypothesis`` column in the table indicates whether the null hypothesis was
    rejected, and that a feature was therefore observed to be differentially
    abundant across the groups.

    >>> ancom_df['Reject null hypothesis']
    b1    False
    b2     True
    b3    False
    b4    False
    b5    False
    b6    False
    b7    False
    Name: Reject null hypothesis, dtype: bool

    From this we can conclude that only ``b2`` was significantly different in
    abundance between the treatment and the placebo. We still don't know, for
    example, in which group ``b2`` was more abundant. We therefore may next be
    interested in comparing the abundance of ``b2`` across the two groups.
    We can do that using the second DataFrame that was returned. Here we
    compare the median (50th percentile) abundance of ``b2`` in the treatment
    and placebo groups:

    >>> percentile_df[50.0].loc['b2']
    Group
    placebo      21.0
    treatment    11.0
    Name: b2, dtype: float64

    We can also look at a full five-number summary for ``b2`` in the treatment
    and placebo groups:

    >>> percentile_df.loc['b2'] # doctest: +NORMALIZE_WHITESPACE
    Percentile  Group
    0.0         placebo      21.0
    25.0        placebo      21.0
    50.0        placebo      21.0
    75.0        placebo      21.5
    100.0       placebo      22.0
    0.0         treatment    11.0
    25.0        treatment    11.0
    50.0        treatment    11.0
    75.0        treatment    11.0
    100.0       treatment    11.0
    Name: b2, dtype: float64

    Taken together, these data tell us that ``b2`` is present in significantly
    higher abundance in the placebo group samples than in the treatment group
    samples.

    """
    if not isinstance(table, pd.DataFrame):
        raise TypeError(
            "`table` must be a `pd.DataFrame`, " "not %r." % type(table).__name__
        )
    if not isinstance(grouping, pd.Series):
        raise TypeError(
            "`grouping` must be a `pd.Series`," " not %r." % type(grouping).__name__
        )

    if np.any(table <= 0):
        raise ValueError(
            "Cannot handle zeros or negative values in `table`. "
            "Use pseudocounts or ``multi_replace``."
        )

    if not 0 < alpha < 1:
        raise ValueError("`alpha`=%f is not within 0 and 1." % alpha)

    if not 0 < tau < 1:
        raise ValueError("`tau`=%f is not within 0 and 1." % tau)

    if not 0 < theta < 1:
        raise ValueError("`theta`=%f is not within 0 and 1." % theta)

    # @deprecated
    if multiple_comparisons_correction != "holm-bonferroni":
        _warn_deprecated(ancom, "0.6.0")
        p_adjust = multiple_comparisons_correction

    if (grouping.isnull()).any():
        raise ValueError("Cannot handle missing values in `grouping`.")

    if (table.isnull()).any().any():
        raise ValueError("Cannot handle missing values in `table`.")

    percentiles = list(percentiles)
    for percentile in percentiles:
        if not 0.0 <= percentile <= 100.0:
            raise ValueError(
                "Percentiles must be in the range [0, 100], %r "
                "was provided." % percentile
            )

    duplicates = find_duplicates(percentiles)
    if duplicates:
        formatted_duplicates = ", ".join(repr(e) for e in duplicates)
        raise ValueError(
            "Percentile values must be unique. The following"
            " value(s) were duplicated: %s." % formatted_duplicates
        )

    groups = np.unique(grouping)
    num_groups = len(groups)

    if num_groups == len(grouping):
        raise ValueError(
            "All values in `grouping` are unique. This method cannot "
            "operate on a grouping vector with only unique values (e.g., "
            "there are no 'within' variance because each group of samples "
            "contains only a single sample)."
        )

    if num_groups == 1:
        raise ValueError(
            "All values the `grouping` are the same. This method cannot "
            "operate on a grouping vector with only a single group of samples"
            "(e.g., there are no 'between' variance because there is only a "
            "single group)."
        )

    # @deprecated
    if significance_test is None:
        significance_test = "f_oneway"

    table_index_len = len(table.index)
    grouping_index_len = len(grouping.index)
    mat, cats = table.align(grouping, axis=0, join="inner")
    if len(mat) != table_index_len or len(cats) != grouping_index_len:
        raise ValueError("`table` index and `grouping` " "index must be consistent.")

    n_feat = mat.shape[1]

    _logratio_mat = _log_compare(mat.values, cats.values, significance_test)
    logratio_mat = _logratio_mat + _logratio_mat.T

    # Multiple comparisons
    if p_adjust is not None:
        logratio_mat = np.apply_along_axis(
            lambda arr: _calc_p_adjust(p_adjust, arr), 1, logratio_mat
        )

    np.fill_diagonal(logratio_mat, 1)
    W = (logratio_mat < alpha).sum(axis=1)
    c_start = W.max() / n_feat
    if c_start < theta:
        reject = np.zeros_like(W, dtype=bool)
    else:
        # Select appropriate cutoff
        cutoff = c_start - np.linspace(0.05, 0.25, 5)
        prop_cut = np.array([(W > n_feat * cut).mean() for cut in cutoff])
        dels = np.abs(prop_cut - np.roll(prop_cut, -1))
        dels[-1] = 0

        if (dels[0] < tau) and (dels[1] < tau) and (dels[2] < tau):
            nu = cutoff[1]
        elif (dels[0] >= tau) and (dels[1] < tau) and (dels[2] < tau):
            nu = cutoff[2]
        elif (dels[1] >= tau) and (dels[2] < tau) and (dels[3] < tau):
            nu = cutoff[3]
        else:
            nu = cutoff[4]
        reject = W >= nu * n_feat

    feat_ids = mat.columns
    ancom_df = pd.DataFrame(
        {
            "W": pd.Series(W, index=feat_ids),
            "Reject null hypothesis": pd.Series(reject, index=feat_ids),
        }
    )

    if len(percentiles) == 0:
        return ancom_df, pd.DataFrame()
    else:
        data = []
        columns = []
        for group in groups:
            feat_dists = mat[cats == group]
            for percentile in percentiles:
                columns.append((percentile, group))
                data.append(np.percentile(feat_dists, percentile, axis=0))
        columns = pd.MultiIndex.from_tuples(columns, names=["Percentile", "Group"])
        percentile_df = pd.DataFrame(
            np.asarray(data).T, columns=columns, index=feat_ids
        )
        return ancom_df, percentile_df


def _log_compare(mat, cats, test="ttest_ind"):
    """Calculate pairwise log ratios and perform a significance test.

    Calculate pairwise log ratios between all features and perform a
    significance test (i.e. *t*-test) to determine if there is a significant
    difference in feature ratios with respect to the variable of interest.

    Parameters
    ----------
    mat : array_like of shape (n_samples, n_features)
        A matrix of proportions.
    cats : array_like of shape (n_samples,)
        A vector of categories.
    test : str or callable
        Statistical test to run.

    Returns
    -------
    log_ratio : ndarray
        Log ratio *p*-value matrix.

    Raises
    ------
    ValueError
        If specified test name is not a function under ``scipy.stats``.

    """
    c = mat.shape[1]
    log_ratio = np.zeros((c, c))
    log_mat = np.log(mat)
    cs = np.unique(cats)

    if isinstance(test, str):
        try:
            test = getattr(scipy.stats, test)
        except AttributeError:
            raise ValueError(f'Function "{test}" does not exist under scipy.stats.')

    def func(x):
        return test(*[x[cats == k] for k in cs])

    for i in range(c - 1):
        ratio = (log_mat[:, i].T - log_mat[:, i + 1 :].T).T
        _, p = np.apply_along_axis(func, axis=0, arr=ratio)
        log_ratio[i, i + 1 :] = np.squeeze(np.array(p.T))
    return log_ratio


def _gram_schmidt_basis(n):
    """Build clr-transformed basis derived from Gram-Schmidt orthogonalization.

    Parameters
    ----------
    n : int
        Dimension of the Aitchison simplex.

    Returns
    -------
    basis : array_like of shape (n - 1, n)
        Basis matrix.

    """
    basis = np.zeros((n, n - 1))
    for j in range(n - 1):
        i = j + 1
        e = np.array([(1 / i)] * i + [-1] + [0] * (n - i - 1)) * np.sqrt(i / (i + 1))
        basis[:, j] = e
    return basis.T


def sbp_basis(sbp):
    r"""Build an orthogonal basis from a sequential binary partition (SBP).

    A SBP is a hierarchical collection of binary divisions of compositional
    parts ([1]_). The child groups are divided again until all groups contain a
    single part. The SBP can be encoded in a :math:`(D - 1) \times D` matrix
    where, for each row, parts can be grouped by -1 and +1 tags, and 0 for
    excluded parts. The *i*-th balance is computed as follows:

    .. math::
        b_i = \sqrt{ \frac{r_i s_i}{r_i+s_i} }
        \ln \left( \frac{g(x_{r_i})}{g(x_{s_i})} \right)

    where :math:`b_i` is the *i*-th balance corresponding to the *i*-th row in
    the SBP, :math:`r_i` and :math:`s_i` and the number of respectively ``+1``
    and ``-1`` labels in the *i*-th row of the SBP and where :math:`g(x) =
    (\prod\limits_{i=1}^{D} x_i)^{1/D}` is the geometric mean of :math:`x`.

    Parameters
    ----------
    sbp : array_like of shape (n_partitions, n_features)
        A contrast matrix, also known as a sequential binary partition, where
        every row represents a partition between two groups of features. A part
        labelled ``+1`` would correspond to that feature being in the numerator
        of the given row partition, a part labelled ``-1`` would correspond to
        features being in the denominator of that given row partition, and
        ``0`` would correspond to features excluded in the row partition.

    Returns
    -------
    ndarray of shape (n_partitions, n_features)
        An orthonormal basis in the Aitchison simplex.

    Notes
    -----
    The ``sbp_basis`` method was derived from the ``gsi.buildilrBase()``
    function implemented in the R package "compositions" [2]_.

    Examples
    --------
    >>> import numpy as np
    >>> sbp = np.array([[1, 1,-1,-1,-1],
    ...                 [1,-1, 0, 0, 0],
    ...                 [0, 0, 1,-1,-1],
    ...                 [0, 0, 0, 1,-1]])
    ...
    >>> sbp_basis(sbp)
    array([[ 0.54772256,  0.54772256, -0.36514837, -0.36514837, -0.36514837],
           [ 0.70710678, -0.70710678,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.81649658, -0.40824829, -0.40824829],
           [ 0.        ,  0.        ,  0.        ,  0.70710678, -0.70710678]])

    References
    ----------
    .. [1] Parent, S.É., Parent, L.E., Egozcue, J.J., Rozane, D.E.,
       Hernandes, A., Lapointe, L., Hébert-Gentile, V., Naess, K.,
       Marchand, S., Lafond, J., Mattos, D., Barlow, P., Natale, W., 2013.
       The plant ionome revisited by the nutrient balance concept.
       Front. Plant Sci. 4, 39.
    .. [2] van den Boogaart, K. Gerald, Tolosana-Delgado, Raimon and Bren,
       Matevz, 2014. `compositions`: Compositional Data Analysis. R package
       version 1.40-1. https://CRAN.R-project.org/package=compositions.

    """
    n_pos = (sbp == 1).sum(axis=1)
    n_neg = (sbp == -1).sum(axis=1)
    psi = np.zeros(sbp.shape)
    for i in range(0, sbp.shape[0]):
        psi[i, :] = sbp[i, :] * np.sqrt(
            (n_neg[i] / n_pos[i]) ** sbp[i, :] / np.sum(np.abs(sbp[i, :]))
        )
    return psi


def _check_orthogonality(basis):
    r"""Check to see if basis is truly orthonormal in the Aitchison simplex.

    Parameters
    ----------
    basis : ndarray
        Basis in the Aitchison simplex of dimension :math:`(D - 1) \times D`.

    """
    basis = np.atleast_2d(basis)
    if not np.allclose(basis @ basis.T, np.identity(len(basis)), rtol=1e-4, atol=1e-6):
        raise ValueError("Basis is not orthonormal.")


def _welch_ttest(a, b):
    r"""Perform Welch's *t*-test on two samples of unequal variances.

    Parameters
    ----------
    a, b : 1-D array_like
        Samples to test.

    Returns
    -------
    pd.DataFrame
        Test result. Columns are: T statistic, df, pvalue, Difference, CI(2.5),
        CI(97.5).

    See Also
    --------
    scipy.stats.ttest_ind
    statsmodels.stats.weightstats.CompareMeans

    Notes
    -----
    Compared with ``scipy.stats.ttest_ind`` with ``equal_var=False``, this
    function additionally returns confidence intervals. This implementation
    uses the ``CompareMeans`` class from ``statsmodels.stats.weightstats``.

    """
    # See https://stats.stackexchange.com/a/475345
    # See https://www.statsmodels.org/dev/generated/statsmodels.stats.weightstats.CompareMeans.html

    # Creating a CompareMeans object to perform Welch's t-test
    statsmodel_cm_object = CompareMeans.from_data(
        data1=a, data2=b, weights1=None, weights2=None
    )

    # Performing Welch's t-test using the object to obtain tstat, pvalue, and df
    ttest_cm_result = statsmodel_cm_object.ttest_ind(
        alternative="two-sided", usevar="unequal", value=0
    )

    tstat = ttest_cm_result[0]
    p = ttest_cm_result[1]
    df = ttest_cm_result[2]

    # Calculating difference between the two means
    m1 = np.mean(a)
    m2 = np.mean(b)

    delta = m1 - m2

    # Calculating confidence intervals using the aformentioned CompareMeans object
    conf_int = statsmodel_cm_object.tconfint_diff(
        alpha=0.05, alternative="two-sided", usevar="unequal"
    )

    lb = conf_int[0]
    ub = conf_int[1]

    return pd.DataFrame(
        np.array([tstat, df, p, delta, lb, ub]).reshape(1, -1),
        columns=["T statistic", "df", "pvalue", "Difference", "CI(2.5)", "CI(97.5)"],
    )


def dirmult_ttest(
    table,
    grouping,
    treatment,
    reference,
    pseudocount=0.5,
    draws=128,
    p_adjust="holm",
    seed=None,
):
    r"""*T*-test using Dirichlet-multinomial distribution.

    The Dirichlet-multinomial distribution is a compound distribution that
    combines a Dirichlet distribution over the probabilities of a multinomial
    distribution. This distribution is used to model the distribution of
    species abundances in a community.

    To perform the *t*-test, we first fit a Dirichlet-multinomial distribution
    for each sample, and then we compute the fold change and *p*-value for each
    feature. The fold change is computed as the difference between the
    samples of the two groups. *t*-tests are then performed on the posterior
    samples, drawn from each Dirichlet-multinomial distribution. The
    log-fold changes as well as their credible intervals, the *p*-values and
    the multiple comparison corrected *p*-values are reported.

    This process mirrors the approach performed by the R package "ALDEx2" [1]_.

    Parameters
    ----------
    table : pd.DataFrame
        Contingency table of counts where rows are features and columns are samples.
    grouping : pd.Series
        Vector indicating the assignment of samples to groups. For example,
        these could be strings or integers denoting which group a sample
        belongs to. It must be the same length as the samples in ``table``.
        The index must be the same on ``table`` and ``grouping`` but need not be
        in the same order. The *t*-test is computed between the ``treatment``
        group and the ``reference`` group specified in the ``grouping`` vector.
    treatment : str
        Name of the treatment group.
    reference : str
        Name of the reference group.
    pseudocount : float, optional
        A non-zero value added to the input counts to ensure that all of the
        estimated abundances are strictly greater than zero.
    draws : int, optional
        The number of draws from the Dirichilet-multinomial posterior distribution
        More draws provide higher uncertainty surrounding the estimated
        log-fold changes and *p*-values.
    p_adjust : str or None, optional
        Method to correct *p*-values for multiple comparisons. Options are Holm-
        Boniferroni ("holm" or "holm-bonferroni") (default), Benjamini-
        Hochberg ("bh", "fdr_bh" or "benjamini-hochberg"), or any method supported
        by statsmodels' ``multipletests`` function. Case-insensitive. If None, no
        correction will be performed.
    seed : int or np.random.Generator, optional
        A user-provided random seed or random generator instance.

    Returns
    -------
    pd.DataFrame
        A table of features, their log-fold changes and other relevant statistics.

        ``T statistic`` is the *t*-statistic outputted from the *t*-test. *t*-statistics
        are generated from each posterior draw.  The reported ``T statistic`` is the
        average across all of the posterior draws.

        ``df`` is the degrees of freedom from the *t*-test.

        ``Log2(FC)`` is the expected log2-fold change. Within each posterior draw
        the log2 fold-change is computed as the difference between the mean
        log-abundance the ``treatment`` group and the ``reference`` group. All log2
        fold changes are expressed in clr coordinates. The reported ``Log2(FC)``
        is the average of all of the log2-fold changes computed from each of the
        posterior draws.

        ``CI(2.5)`` is the 2.5% quantile of the log2-fold change. The reported
        ``CI(2.5)`` is the 2.5% quantile of all of the log2-fold changes computed
        from each of the posterior draws.

        ``CI(97.5)`` is the 97.5% quantile of the log2-fold change. The
        reported ``CI(97.5)`` is the 97.5% quantile of all of the log2-fold
        changes computed from each of the posterior draws.

        ``pvalue`` is the *p*-value of the *t*-test. The reported values are the
        average of all of the *p*-values computed from the *t*-tests calculated
        across all of the posterior draws.

        ``qvalue`` is the *p*-value of the *t*-test after performing multiple
        comparison correction.

        ``Reject null hypothesis`` indicates if feature is differentially
        abundant across groups (``True``) or not (``False``). In order for a
        feature to be differentially abundant, the qvalue needs to be significant
        (i.e. <0.05) and the confidence intervals reported by ``CI(2.5)`` and
        ``CI(97.5)`` must not overlap with zero.

    See Also
    --------
    scipy.stats.ttest_ind

    Notes
    -----
    The confidence intervals are computed using the mininum 2.5% and maximum
    97.5% bounds computed across all of the posterior draws.

    The reference frame here is the geometric mean. Extracting absolute log
    fold changes from this test assumes that the average feature abundance
    between the ``treatment`` and the ``reference`` groups are the same. If this
    assumption is violated, then the log-fold changes will be biased, and the
    *p*-values will not be reliable. However, the bias is the same across each
    feature, as a result the ordering of the log-fold changes can still be useful.

    One benefit of using the Dirichlet-multinomial distribution is that the
    statistical power increases with regards to the abundance magnitude. More counts
    per sample will shrink the size of the confidence intervals, and can result in
    lower *p*-values.

    References
    ----------
    .. [1] Fernandes et al. "Unifying the analysis of
       high-throughput sequencing datasets: characterizing RNA-seq,
       16S rRNA gene sequencing and selective growth experiments by
       compositional data analysis." Microbiome (2014).

    Examples
    --------
    >>> import pandas as pd
    >>> from skbio.stats.composition import dirmult_ttest
    >>> table = pd.DataFrame([[20,  110, 100, 101, 100, 103, 104],
    ...                       [33,  110, 120, 100, 101, 100, 102],
    ...                       [12,  110, 100, 110, 100, 50,  90],
    ...                       [202, 201, 9,  10, 10, 11, 11],
    ...                       [200, 202, 10, 10, 13, 10, 10],
    ...                       [203, 201, 14, 10, 10, 13, 12]],
    ...                      index=['s1', 's2', 's3', 's4', 's5', 's6'],
    ...                      columns=['b1', 'b2', 'b3', 'b4', 'b5', 'b6',
    ...                               'b7'])
    >>> grouping = pd.Series(['treatment', 'treatment', 'treatment',
    ...                       'placebo', 'placebo', 'placebo'],
    ...                      index=['s1', 's2', 's3', 's4', 's5', 's6'])
    >>> lfc_result = dirmult_ttest(table, grouping, 'treatment', 'placebo',
    ...                            seed=0)
    >>> lfc_result[["Log2(FC)", "CI(2.5)", "CI(97.5)", "qvalue"]]
        Log2(FC)   CI(2.5)  CI(97.5)    qvalue
    b1 -4.991987 -7.884498 -2.293463  0.020131
    b2 -2.533729 -3.594590 -1.462339  0.007446
    b3  1.627677 -1.048219  4.750792  0.068310
    b4  1.707221 -0.467481  4.164998  0.065613
    b5  1.528243 -1.036910  3.978387  0.068310
    b6  1.182343 -0.702656  3.556061  0.068310
    b7  1.480232 -0.601277  4.043888  0.068310

    """
    rng = get_rng(seed)
    if not isinstance(table, pd.DataFrame):
        raise TypeError(
            "`table` must be a `pd.DataFrame`, " "not %r." % type(table).__name__
        )
    if not isinstance(grouping, pd.Series):
        raise TypeError(
            "`grouping` must be a `pd.Series`," " not %r." % type(grouping).__name__
        )

    if np.any(table < 0):
        raise ValueError("Cannot handle negative values in `table`. ")

    if (grouping.isnull()).any():
        raise ValueError("Cannot handle missing values in `grouping`.")

    if (table.isnull()).any().any():
        raise ValueError("Cannot handle missing values in `table`.")

    table_index_len = len(table.index)
    grouping_index_len = len(grouping.index)
    mat, cats = table.align(grouping, axis=0, join="inner")
    if len(mat) != table_index_len or len(cats) != grouping_index_len:
        raise ValueError("`table` index and `grouping` " "index must be consistent.")

    trt_group = grouping.loc[grouping == treatment]
    ref_group = grouping.loc[grouping == reference]
    posterior = [
        rng.dirichlet(table.values[i] + pseudocount) for i in range(table.shape[0])
    ]
    dir_table = pd.DataFrame(clr(posterior), index=table.index, columns=table.columns)
    res = [
        _welch_ttest(
            np.array(dir_table.loc[trt_group.index, x].values),
            np.array(dir_table.loc[ref_group.index, x].values),
        )
        for x in table.columns
    ]
    res = pd.concat(res)
    for i in range(1, draws):
        posterior = [
            rng.dirichlet(table.values[i] + pseudocount) for i in range(table.shape[0])
        ]
        dir_table = pd.DataFrame(
            clr(posterior), index=table.index, columns=table.columns
        )

        ires = [
            _welch_ttest(
                np.array(dir_table.loc[trt_group.index, x].values),
                np.array(dir_table.loc[ref_group.index, x].values),
            )
            for x in table.columns
        ]
        ires = pd.concat(ires)
        # online average to avoid holding all of the results in memory
        res["Difference"] = (i * res["Difference"] + ires["Difference"]) / (i + 1)
        res["pvalue"] = (i * res["pvalue"] + ires["pvalue"]) / (i + 1)
        res["CI(2.5)"] = np.minimum(res["CI(2.5)"], ires["CI(2.5)"])
        res["CI(97.5)"] = np.maximum(res["CI(97.5)"], ires["CI(97.5)"])
        res["T statistic"] = (i * res["T statistic"] + ires["T statistic"]) / (i + 1)

    res.index = table.columns
    # convert all log fold changes to base 2
    res["Difference"] = res["Difference"] / np.log(2)
    res["CI(2.5)"] = res["CI(2.5)"] / np.log(2)
    res["CI(97.5)"] = res["CI(97.5)"] / np.log(2)

    # multiple comparison
    if p_adjust is not None:
        qval = _calc_p_adjust(p_adjust, res["pvalue"])
    else:
        qval = res["pvalue"].values

    # test to see if confidence interval includes 0.
    sig = np.logical_or(
        np.logical_and(res["CI(2.5)"] > 0, res["CI(97.5)"] > 0),
        np.logical_and(res["CI(2.5)"] < 0, res["CI(97.5)"] < 0),
    )

    reject = np.logical_and(qval[0], sig)

    res = res.rename(columns={"Difference": "Log2(FC)"})
    res["qvalue"] = qval
    res["Reject null hypothesis"] = reject

    col_order = [
        "T statistic",
        "df",
        "Log2(FC)",
        "CI(2.5)",
        "CI(97.5)",
        "pvalue",
        "qvalue",
        "Reject null hypothesis",
    ]
    return res[col_order]