New upstream version 0.5.4

# scikit-bio changelog

## Version 0.5.4 (2018-08-23)

### Features

* Added `FSVD`, an alternative fast heuristic method to perform Principal Coordinates Analysis, to `skbio.stats.ordination.pcoa`.

### Backward-incompatible changes [stable]

### Backward-incompatible changes [experimental]

### Performance enhancements

* Added optimized utility methods `f_matrix_inplace` and `e_matrix_inplace` which perform `f_matrix` and `e_matrix` computations in-place and are used by the new `center_distance_matrix` method in `skbio.stats.ordination`.

### Bug fixes

### Deprecated functionality [stable]

### Deprecated functionality [experimental]

### Miscellaneous

## Version 0.5.3 (2018-08-07)

### Features

           'TreeNode', 'nj', 'read', 'write', 'OrdinationResults']

__credits__ = "https://github.com/biocore/scikit-bio/graphs/contributors"

__version__ = "0.5.3"

__version__ = "0.5.4"

mottos = [

    # 03/15/2014

from ._canonical_correspondence_analysis import cca

from ._principal_coordinate_analysis import pcoa, pcoa_biplot

from ._ordination_results import OrdinationResults

from ._utils import (mean_and_std, scale, svd_rank, corr, e_matrix, f_matrix)

from ._utils import (mean_and_std, scale, svd_rank, corr, e_matrix, f_matrix,

                     center_distance_matrix)

__all__ = ['ca', 'rda', 'cca', 'pcoa', 'pcoa_biplot', 'OrdinationResults',

           'mean_and_std', 'scale', 'svd_rank', 'corr',

           'e_matrix', 'f_matrix']

           'e_matrix', 'f_matrix', 'center_distance_matrix']

test = TestRunner(__file__).test

# The full license is in the file COPYING.txt, distributed with this software.

# ----------------------------------------------------------------------------

from warnings import warn

import pandas as pd

import numpy as np

import pandas as pd

from numpy import dot, hstack

from numpy.linalg import qr, svd

from numpy.random import standard_normal

from scipy.linalg import eigh

from warnings import warn

from skbio.stats.distance import DistanceMatrix

from skbio.util._decorator import experimental

from ._ordination_results import OrdinationResults

from ._utils import e_matrix, f_matrix, scale

# - In cogent, after computing eigenvalues/vectors, the imaginary part

#   is dropped, if any. We know for a fact that the eigenvalues are

#   real, so that's not necessary, but eigenvectors can in principle

#   be complex (see for example

#   http://math.stackexchange.com/a/47807/109129 for details) and in

#   that case dropping the imaginary part means they'd no longer be

#   so, so I'm not doing that.

from ._utils import center_distance_matrix, scale

@experimental(as_of="0.4.0")

def pcoa(distance_matrix):

def pcoa(distance_matrix, method="eigh", number_of_dimensions=0,

         inplace=False):

    r"""Perform Principal Coordinate Analysis.

    Principal Coordinate Analysis (PCoA) is a method similar to PCA

    that works from distance matrices, and so it can be used with

    ecologically meaningful distances like UniFrac for bacteria.

    Principal Coordinate Analysis (PCoA) is a method similar

    to Principal Components Analysis (PCA) with the difference that PCoA

    operates on distance matrices, typically with non-euclidian and thus

    ecologically meaningful distances like UniFrac in microbiome research.

    In ecology, the euclidean distance preserved by Principal

    Component Analysis (PCA) is often not a good choice because it

    species is present in two sites, that means that the sites are

    similar.).

    Note that the returned eigenvectors are not normalized to unit length.

    Parameters

    ----------

    distance_matrix : DistanceMatrix

        A distance matrix.

    method : str, optional

        Eigendecomposition method to use in performing PCoA.

        By default, uses SciPy's `eigh`, which computes exact

        eigenvectors and eigenvalues for all dimensions. The alternate

        method, `fsvd`, uses faster heuristic eigendecomposition but loses

        accuracy. The magnitude of accuracy lost is dependent on dataset.

    number_of_dimensions : int, optional

        Dimensions to reduce the distance matrix to. This number determines

        how many eigenvectors and eigenvalues will be returned.

        By default, equal to the number of dimensions of the distance matrix,

        as default eigendecomposition using SciPy's `eigh` method computes

        all eigenvectors and eigenvalues. If using fast heuristic

        eigendecomposition through `fsvd`, a desired number of dimensions

        should be specified. Note that the default eigendecomposition

        method `eigh` does not natively support a specifying number of

        dimensions to reduce a matrix to, so if this parameter is specified,

        all eigenvectors and eigenvalues will be simply be computed with no

        speed gain, and only the number specified by `number_of_dimensions`

        will be returned. Specifying a value of `0`, the default, will

        set `number_of_dimensions` equal to the number of dimensions of the

        specified `distance_matrix`.

    inplace : bool, optional

        If true, centers a distance matrix in-place in a manner that reduces

        memory consumption.

    Returns

    -------

    Notes

    -----

    It is sometimes known as metric multidimensional scaling or

    classical scaling.

    .. note::

       If the distance is not euclidean (for example if it is a

    .. note:: If the distance is not euclidean (for example if it is a

        semimetric and the triangle inequality doesn't hold),

        negative eigenvalues can appear. There are different ways

        to deal with that problem (see Legendre & Legendre 1998, \S

        9.2.3), but none are currently implemented here.

        However, a warning is raised whenever negative eigenvalues

        appear, allowing the user to decide if they can be safely

        ignored.

    """

    distance_matrix = DistanceMatrix(distance_matrix)

    E_matrix = e_matrix(distance_matrix.data)

    # If the used distance was euclidean, pairwise distances

    # needn't be computed from the data table Y because F_matrix =

    # Y.dot(Y.T) (if Y has been centred).

    F_matrix = f_matrix(E_matrix)

    # If the eigendecomposition ever became a bottleneck, it could

    # be replaced with an iterative version that computes the

    # largest k eigenvectors.

    eigvals, eigvecs = eigh(F_matrix)

    # eigvals might not be ordered, so we order them (at least one

    # is zero). cogent makes eigenvalues positive by taking the

    # Center distance matrix, a requirement for PCoA here

    matrix_data = center_distance_matrix(distance_matrix.data, inplace=inplace)

    # If no dimension specified, by default will compute all eigenvectors

    # and eigenvalues

    if number_of_dimensions == 0:

        if method == "fsvd" and matrix_data.shape[0] > 10:

            warn("FSVD: since no value for number_of_dimensions is specified, "

                 "PCoA for all dimensions will be computed, which may "

                 "result in long computation time if the original "

                 "distance matrix is large.", RuntimeWarning)

        # distance_matrix is guaranteed to be square

        number_of_dimensions = matrix_data.shape[0]

    elif number_of_dimensions < 0:

        raise ValueError('Invalid operation: cannot reduce distance matrix '

                         'to negative dimensions using PCoA. Did you intend '

                         'to specify the default value "0", which sets '

                         'the number_of_dimensions equal to the '

                         'dimensionality of the given distance matrix?')

    # Perform eigendecomposition

    if method == "eigh":

        # eigh does not natively support specifying number_of_dimensions, i.e.

        # there are no speed gains unlike in FSVD. Later, we slice off unwanted

        # dimensions to conform the result of eigh to the specified

        # number_of_dimensions.

        eigvals, eigvecs = eigh(matrix_data)

        long_method_name = "Principal Coordinate Analysis"

    elif method == "fsvd":

        eigvals, eigvecs = _fsvd(matrix_data, number_of_dimensions)

        long_method_name = "Approximate Principal Coordinate Analysis " \

                           "using FSVD"

    else:

        raise ValueError(

            "PCoA eigendecomposition method {} not supported.".format(method))

    # cogent makes eigenvalues positive by taking the

    # abs value, but that doesn't seem to be an approach accepted

    # by L&L to deal with negative eigenvalues. We raise a warning

    # in that case. First, we make values close to 0 equal to 0.

                                                  eigvals.max()),

            RuntimeWarning

        )

    # eigvals might not be ordered, so we first sort them, then analogously

    # sort the eigenvectors by the ordering of the eigenvalues too

    idxs_descending = eigvals.argsort()[::-1]

    eigvals = eigvals[idxs_descending]

    eigvecs = eigvecs[:, idxs_descending]

    # Scale eigenvalues to have lenght = sqrt(eigenvalue). This

    # works because np.linalg.eigh returns normalized

    # eigenvectors. Each row contains the coordinates of the

    # objects in the space of principal coordinates. Note that at

    # least one eigenvalue is zero because only n-1 axes are

    # needed to represent n points in an euclidean space.

    # If we return only the coordinates that make sense (i.e., that have a

    # corresponding positive eigenvalue), then Jackknifed Beta Diversity

    # won't work as it expects all the OrdinationResults to have the same

    eigvecs[:, num_positive:] = np.zeros(eigvecs[:, num_positive:].shape)

    eigvals[num_positive:] = np.zeros(eigvals[num_positive:].shape)

    if method == "fsvd":

        # Since the dimension parameter, hereafter referred to as 'd',

        # restricts the number of eigenvalues and eigenvectors that FSVD

        # computes, we need to use an alternative method to compute the sum

        # of all eigenvalues, used to compute the array of proportions

        # explained. Otherwise, the proportions calculated will only be

        # relative to d number of dimensions computed; whereas we want

        # it to be relative to the entire dimensionality of the

        # centered distance matrix.

        # An alternative method of calculating th sum of eigenvalues is by

        # computing the trace of the centered distance matrix.

        # See proof outlined here: https://goo.gl/VAYiXx

        sum_eigenvalues = np.trace(matrix_data)

    else:

        # Calculate proportions the usual way

        sum_eigenvalues = np.sum(eigvals)

    proportion_explained = eigvals / sum_eigenvalues

    # In case eigh is used, eigh computes all eigenvectors and -values.

    # So if number_of_dimensions was specified, we manually need to ensure

    # only the requested number of dimensions

    # (number of eigenvectors and eigenvalues, respectively) are returned.

    eigvecs = eigvecs[:, :number_of_dimensions]

    eigvals = eigvals[:number_of_dimensions]

    proportion_explained = proportion_explained[:number_of_dimensions]

    # Scale eigenvalues to have length = sqrt(eigenvalue). This

    # works because np.linalg.eigh returns normalized

    # eigenvectors. Each row contains the coordinates of the

    # objects in the space of principal coordinates. Note that at

    # least one eigenvalue is zero because only n-1 axes are

    # needed to represent n points in a euclidean space.

    coordinates = eigvecs * np.sqrt(eigvals)

    proportion_explained = eigvals / eigvals.sum()

    axis_labels = ['PC%d' % i for i in range(1, eigvals.size + 1)]

    axis_labels = ["PC%d" % i for i in range(1, number_of_dimensions + 1)]

    return OrdinationResults(

        short_method_name='PCoA',

        long_method_name='Principal Coordinate Analysis',

        short_method_name="PCoA",

        long_method_name=long_method_name,

        eigvals=pd.Series(eigvals, index=axis_labels),

        samples=pd.DataFrame(coordinates, index=distance_matrix.ids,

                             columns=axis_labels),

                                       index=axis_labels))

def _fsvd(centered_distance_matrix, number_of_dimensions=10):

    """

    Performs singular value decomposition, or more specifically in

    this case eigendecomposition, using fast heuristic algorithm

    nicknamed "FSVD" (FastSVD), adapted and optimized from the algorithm

    described by Halko et al (2011).

    Parameters

    ----------

    centered_distance_matrix : np.array

       Numpy matrix representing the distance matrix for which the

       eigenvectors and eigenvalues shall be computed

    number_of_dimensions : int

       Number of dimensions to keep. Must be lower than or equal to the

       rank of the given distance_matrix.

    Returns

    -------

    np.array

       Array of eigenvectors, each with number_of_dimensions length.

    np.array

       Array of eigenvalues, a total number of number_of_dimensions.

    Notes

    -----

    The algorithm is based on 'An Algorithm for the Principal

    Component analysis of Large Data Sets'

    by N. Halko, P.G. Martinsson, Y. Shkolnisky, and M. Tygert.

    Original Paper: https://arxiv.org/abs/1007.5510

    Ported from MATLAB implementation described here:

    https://stats.stackexchange.com/a/11934/211065

    """

    m, n = centered_distance_matrix.shape

    # Number of levels of the Krylov method to use.

    # For most applications, num_levels=1 or num_levels=2 is sufficient.

    num_levels = 1

    # Changes the power of the spectral norm, thus minimizing the error).

    use_power_method = False

    # Note: a (conjugate) transpose is removed for performance, since we

    # only expect square matrices.

    if m != n:

        raise ValueError('FSVD expects square distance matrix')

    if number_of_dimensions > m or number_of_dimensions > n:

        raise ValueError('FSVD: number_of_dimensions cannot be larger than'

                         ' the dimensionality of the given distance matrix.')

    if number_of_dimensions < 0:

        raise ValueError('Invalid operation: cannot reduce distance matrix '

                         'to negative dimensions using PCoA. Did you intend '

                         'to specify the default value "0", which sets '

                         'the number_of_dimensions equal to the '

                         'dimensionality of the given distance matrix?')

    k = number_of_dimensions + 2

    # Form a real nxl matrix G whose entries are independent, identically

    # distributed Gaussian random variables of zero mean and unit variance

    G = standard_normal(size=(n, k))

    if use_power_method:

        # use only the given exponent

        H = dot(centered_distance_matrix, G)

        for x in range(2, num_levels + 2):

            # enhance decay of singular values

            # note: distance_matrix is no longer transposed, saves work

            # since we're expecting symmetric, square matrices anyway

            # (Daniel McDonald's changes)

            H = dot(centered_distance_matrix, dot(centered_distance_matrix, H))

    else:

        # compute the m x l matrices H^{(0)}, ..., H^{(i)}

        # Note that this is done implicitly in each iteration below.

        H = dot(centered_distance_matrix, G)

        # to enhance performance

        H = hstack(

            (H,

             dot(centered_distance_matrix, dot(centered_distance_matrix, H))))

        for x in range(3, num_levels + 2):

            tmp = dot(centered_distance_matrix,

                      dot(centered_distance_matrix, H))

            H = hstack(

                (H, dot(centered_distance_matrix,

                        dot(centered_distance_matrix, tmp))))

    # Using the pivoted QR-decomposition, form a real m * ((i+1)l) matrix Q

    # whose columns are orthonormal, s.t. there exists a real

    # ((i+1)l) * ((i+1)l) matrix R for which H = QR

    Q, R = qr(H)

    # Compute the n * ((i+1)l) product matrix T = A^T Q

    T = dot(centered_distance_matrix, Q)  # step 3

    # Form an SVD of T

    Vt, St, W = svd(T, full_matrices=False)

    W = W.transpose()

    # Compute the m * ((i+1)l) product matrix

    Ut = dot(Q, W)

    U_fsvd = Ut[:, :number_of_dimensions]

    S = St[:number_of_dimensions]

    # drop imaginary component, if we got one

    # Note:

    #   In cogent, after computing eigenvalues/vectors, the imaginary part

    #   is dropped, if any. We know for a fact that the eigenvalues are

    #   real, so that's not necessary, but eigenvectors can in principle

    #   be complex (see for example

    #   http://math.stackexchange.com/a/47807/109129 for details)

    eigenvalues = S.real

    eigenvectors = U_fsvd.real

    return eigenvalues, eigenvectors

@experimental(as_of="0.5.3")

def pcoa_biplot(ordination, y):

    """Compute the projection of descriptors into a PCoA matrix

    col_means = E_matrix.mean(axis=0, keepdims=True)

    matrix_mean = E_matrix.mean()

    return E_matrix - row_means - col_means + matrix_mean

def center_distance_matrix(distance_matrix, inplace=False):

    """

    Centers a distance matrix.

    Note: If the used distance was euclidean, pairwise distances

    needn't be computed from the data table Y because F_matrix =

    Y.dot(Y.T) (if Y has been centered).

    But since we're expecting distance_matrix to be non-euclidian,

    we do the following computation as per

    Numerical Ecology (Legendre & Legendre 1998).

    Parameters

    ----------

    distance_matrix : 2D array_like

        Distance matrix.

    inplace : bool, optional

        Whether or not to center the given distance matrix in-place, which

        is more efficient in terms of memory and computation.

    """

    if inplace:

        return _f_matrix_inplace(_e_matrix_inplace(distance_matrix))

    else:

        return f_matrix(e_matrix(distance_matrix))

def _e_matrix_inplace(distance_matrix):

    """

    Compute E matrix from a distance matrix inplace.

    Squares and divides by -2 the input element-wise. Eq. 9.20 in

    Legendre & Legendre 1998.

    Modified from :func:`skbio.stats.ordination.e_matrix` function,

    performing row-wise operations to avoid excessive memory allocations.

    Parameters

    ----------

    distance_matrix : 2D array_like

        Distance matrix.

    """

    distance_matrix = distance_matrix.astype(np.float)

    for i in np.arange(len(distance_matrix)):

        distance_matrix[i] = (distance_matrix[i] * distance_matrix[i]) / -2

    return distance_matrix

def _f_matrix_inplace(e_matrix):

    """

    Compute F matrix from E matrix inplace.

    Centering step: for each element, the mean of the corresponding

    row and column are subtracted, and the mean of the whole

    matrix is added. Eq. 9.21 in Legendre & Legendre 1998.

    Modified from :func:`skbio.stats.ordination.f_matrix` function,

    performing row-wise operations to avoid excessive memory allocations.

    Parameters

    ----------

    e_matrix : 2D array_like

        A matrix representing the "E matrix" as described above.

    """

    e_matrix = e_matrix.astype(np.float)

    row_means = np.zeros(len(e_matrix), dtype=float)

    col_means = np.zeros(len(e_matrix), dtype=float)

    matrix_mean = 0.0

    for i in np.arange(len(e_matrix)):

        row_means[i] = e_matrix[i].mean()

        matrix_mean += e_matrix[i].sum()

        col_means += e_matrix[i]

    matrix_mean /= len(e_matrix) ** 2

    col_means /= len(e_matrix)

    for i in np.arange(len(e_matrix)):

        v = e_matrix[i]

        v -= row_means[i]

        v -= col_means

        v += matrix_mean

        e_matrix[i] = v

    return e_matrix