class DistanceMatrix(Matrix):
    def __init__(self):
        self.df = pd.DataFrame()

    def __eq__(self, other):
        if (np.abs(self.sort().values - other.sort().values) < 1e-10).all():
            return True

    @classmethod
    def from_csv(cls, filename, **kwargs):
        with fileIO.freader(filename) as handle:
            new_instance = cls()
            if not 'index_col' in kwargs:
                kwargs['index_col'] = 0  # new behaviour in pandas
            new_instance.df = pd.read_csv(handle, **kwargs)
        return new_instance

    def to_csv(self, filename, **kwargs):
        with fileIO.fwriter(filename) as handle:
            self.df.to_csv(handle)

    @classmethod
    def from_array(cls, array, names=None):
        nrow, ncol = array.shape
        if nrow != ncol:
            errmsg = 'Not a square array'
            raise ValueError(errmsg)
        new_instance = cls()
        new_instance.df = pd.DataFrame(array)
        try:
            new_instance.set_names(names)
            return new_instance
        except ValueError as err:
            sys.stderr.write(str(err))
            return new_instance

    def get_names(self):
        return [str(x) for x in self.df.index]

    def set_names(self, names):
        if names is None:
            return
        if len(names) != self.df.index.size:
            errmsg = 'Expected {} names, got {}'.format(self.df.index.size, len(names))
            raise ValueError(errmsg)
        self.df.index = self.df.columns = names

    def add_noise(self, std=0.0001):
        ix = np.triu_indices(len(self.df), 1)
        rev_ix = ix[::-1]
        noise = np.random.normal(0, std, len(ix[0]))
        noisy = self.to_array().copy()
        noisy[ix] += noise
        noisy[rev_ix] += noise
        noisy[noisy < 0] = np.abs(noisy[noisy < 0])
        return self.__class__(noisy, self.df.index)

    def affinity(self, mask=None, scale=None):
        return affinity(self.to_array(), mask, scale)

    def double_centre(self, square_input=True):
        return double_centre(self.to_array(), square_input)

    def embedding(self, dimensions, method, **kwargs):
        """
        Embeds the distance matrix in a coordinate space. Implemented methods are:
            cmds: Classical MultiDimensional Scaling
            kpca: Kernel Principal Components Analysis
            mmds: Metric MultiDimensional Scaling
            nmmds: Non-Metric MultiDimensional Scaling
            spectral: Spectral decomposition of Laplacian matrix
            tsne: t-distributed Stochastic Neighbour Embedding

        Valid kwargs:
            kpca: affinity_matrix - a precomputed array of affinities
                  sigma - the value of sigma to use when computing the affinity matrix via
                          the Radial Basis Function
            nmmds: initial_coords - a set of coordinates to refine. NMMDS works very badly
                                    without this
            spectral: affinity_matrix, sigma
                      unit_length - scale the coordinates to unit length, so points sit
                                    on the surface of the unit sphere
        :param dimensions: (int) number of coordinate axes to use
        :param method: (string) one of cmds, kpca, mmds, nmmds, spectral
        :param kwargs: unit_length (bool), affinity_matrix (np.array), sigma (float), initial_coords (np.array)
        :return: coordinate matrix (np.array)
        """
        errors.optioncheck(method, ['cmds', 'kpca', 'mmds', 'nmmds', 'spectral', 'tsne'])
        if method == 'cmds':
            array =  _embedding_classical_mds(self.to_array(), dimensions, **kwargs)
        elif method == 'kpca':
            array = _embedding_kernel_pca(self.to_array(), dimensions, **kwargs)
        elif method == 'mmds':
            array = _embedding_metric_mds(self.to_array(), dimensions)
        elif method == 'nmmds':
            array = _embedding_nonmetric_mds(self.to_array(), dimensions, **kwargs)
        elif method == 'spectral':
            array = _embedding_spectral(self.to_array(), dimensions, **kwargs)
        elif method == 'tsne':
            array = _embedding_tsne(self.to_array(), dimensions, **kwargs)

        return CoordinateMatrix(array, names=self.df.index)

    def reorder(self, new_order):
        reordered_df = self.df.reindex(columns=new_order, index=new_order)
        reordered_names = reordered_df.columns
        newobj = self.__class__()
        newobj.df = reordered_df
        return newobj

    def sort(self):
        order = self.df.index.argsort()
        return self.reorder(order) def rbf(matrix, sigma=1):
    """ Returns the Radial Basis Function kernel from the distance matrix
    """
    return np.exp(-matrix ** 2 / (2 * sigma ** 2)) def binsearch_mask(matrix, logic='or'):
    mink = 1
    maxk = matrix.shape[0]
    guessk = int(np.log(maxk).round())
    result = [guessk, None]
    while maxk - mink != 1:
        test_dist = kdists(matrix, k=guessk)
        test_mask = kmask(matrix, dists=test_dist, logic=logic)
        if isconnected(test_mask) and (test_dist > 1e-6).all():
            maxk = guessk  # either correct or too high
            result = (guessk, test_mask)
            guessk = mink + old_div((guessk - mink), 2)  # try a lower number
        else:
            mink = guessk  # too low
            guessk += old_div((maxk - guessk), 2)

    if result[0] == guessk + 1:
        k, mask = result
        scale = kscale(matrix, k)
        return k, mask, scale
    else:
        k = guessk + 1
        mask = kmask(matrix, k=guessk + 1, logic=logic)
        scale = kscale(matrix, k)
    return k, mask, scale def kmask(matrix, k=7, dists=None, logic='or'):
    """ Creates a boolean mask to include points within k nearest
    neighbours, and exclude the rest.
    Logic can be OR or AND. OR gives the k-nearest-neighbour mask,
    AND gives the mutual k-nearest-neighbour mask."""

    dists = (kdists(matrix, k=k) if dists is None else dists)
    mask = (matrix <= dists)
    if logic == 'or' or logic == '|':
        return mask | mask.T
    elif logic == 'and' or logic == '&':
        return mask & mask.T
    return mask def kscale(matrix, k=7, dists=None):
    """ Returns the local scale based on the k-th nearest neighbour """
    dists = (kdists(matrix, k=k) if dists is None else dists)
    scale = dists.dot(dists.T)
    return scale def affinity(
        matrix,
        mask=None,
        scale=None,
):
    """
    Mask is a 2d boolean matrix. Scale is a 2d local scale matrix,
    as output by kscale(). It's the outer product of the kdists column
    vector produced by kdists.
    """

    mask = (mask if mask is not None else np.ones(matrix.shape, dtype=bool))
    assert isconnected(mask)
    scale = (scale if scale is not None else np.ones(matrix.shape))

    ix = np.where(np.logical_not(mask))
    scaled_matrix = -matrix ** 2 / scale
    # inputs where distance = 0 and scale = 0 result in NaN:
    # the next line replaces NaNs with -1.0
    scaled_matrix[np.where(np.isnan(scaled_matrix))] = -1.0
    affinity_matrix = np.exp(scaled_matrix)
    affinity_matrix[ix] = 0.  # mask
    affinity_matrix.flat[::len(affinity_matrix) + 1] = 0.  # diagonal
    return affinity_matrix def laplace(affinity_matrix, shi_malik_type=False):
    """ Converts affinity matrix into normalised graph Laplacian,
    for spectral clustering.
    (At least) two forms exist:

    L = (D^-0.5).A.(D^-0.5) - default

    L = (D^-1).A - `Shi-Malik` type, from Shi Malik paper"""

    diagonal = affinity_matrix.sum(axis=1) - affinity_matrix.diagonal()
    zeros = diagonal <= 1e-10
    diagonal[zeros] = 1
    if (diagonal <= 1e-10).any():  # arbitrarily small value
        raise ZeroDivisionError
    if shi_malik_type:
        inv_d = np.diag(1 / diagonal)
        return inv_d.dot(affinity_matrix)
    diagonal = np.sqrt(diagonal)
    return affinity_matrix / diagonal / diagonal[:, np.newaxis] def eigen(matrix):
    """ Calculates the eigenvalues and eigenvectors of the input matrix.
    Returns a tuple of (eigenvalues, eigenvectors, cumulative percentage of
    variance explained). Eigenvalues and eigenvectors are sorted in order of
    eigenvalue magnitude, high to low """

    (vals, vecs) = np.linalg.eigh(matrix)
    ind = vals.argsort()[::-1]
    vals = vals[ind]
    vecs = vecs[:, ind]
    vals_ = vals.copy()
    vals_[vals_ < 0] = 0.
    cum_var_exp = np.cumsum(vals_ / vals_.sum())
    return Decomp(matrix.copy(), vals, vecs, cum_var_exp) def double_centre(matrix, square_input=True):
    """ Double-centres the input matrix: From each element: Subtract the row
    mean Subtract the column mean Add the grand mean Divide by -2 Method
    from: Torgerson, W S (1952). Multidimensional scaling: I. Theory and
    method. Alternatively M = -0.5 * (I - 1/n)D[^2](I - 1/n) """

    m = matrix.copy()
    if square_input:
        m **= 2

    (rows, cols) = m.shape
    cm = np.mean(m, axis=0)  # column means
    rm = np.mean(m, axis=1).reshape((rows, 1))  # row means
    gm = np.mean(cm)  # grand mean
    m -= rm + cm - gm
    m /= -2
    return m def normalise_rows(matrix):
    """ Scales all rows to length 1. Fails when row is 0-length, so it
    leaves these unchanged """

    lengths = np.apply_along_axis(np.linalg.norm, 1, matrix)
    if not (lengths > 0).all():
        # raise ValueError('Cannot normalise 0 length vector to length 1')
        # print(matrix)
        lengths[lengths == 0] = 1
    return matrix / lengths[:, np.newaxis] class CoordinateMatrix(Matrix):
    def __init__(self, array, names=None):
        nrow = array.shape[0]
        if names is not None and len(names) != nrow:
            msg = 'Length of names array doesn\'t match the number of rows in the coordinate matrix ({}, {})'.format(
                len(names), nrow)
            raise ValueError(msg)

        self.df = pd.DataFrame(array, index=names) class Partition(object):
    """ Class to store clustering information """

    def __init__(self, partition_vector):
        self._partition_vector = None
        self.partition_vector = partition_vector
        self.membership = self.get_membership()

    def __str__(self):
        return str(self.partition_vector)

    def __repr__(self):
        return self.__class__.__name__ + '({0})'.format(str(self))

    def __len__(self):
        """
        This gives the number of groups in the partition, rather than
        the number of elements in the data
        """
        return max(self.partition_vector) + 1

    def __getitem__(self, index):
        return self.membership[index]

    def __eq__(self, other):
        return self.partition_vector == other.partition_vector

    @property
    def partition_vector(self):
        return self._partition_vector

    @partition_vector.setter
    def partition_vector(self, vec):
        self._partition_vector = self._restricted_growth_notation(vec)

    def is_minimal(self):
        """ The partition describes all members being in the same cluster
        :return: boolean
        """
        return self.num_groups() == 1

    def is_maximal(self):
        """ The partition describes every member being in its own group
        :return: boolean
        """
        return self.num_groups() == self.num_elements()

    @staticmethod
    def _restricted_growth_notation(l):
        """ The clustering returned by the hcluster module gives group
        membership without regard for numerical order This function preserves
        the group membership, but sorts the labelling into numerical order """

        list_length = len(l)

        d = defaultdict(list)
        for (i, element) in enumerate(l):
            d[element].append(i)

        l2 = [None] * list_length

        for (name, index_list) in enumerate(sorted(d.values(), key=min)):
            for index in index_list:
                l2[index] = name

        return tuple(l2)

    @classmethod
    def random(cls, alpha, size):
        """
        Generate a random start using expected proportions, alpha.
        These are used to parameterise a random draw from a Dirichlet
        distribution.
        An example, to split a dataset of 20 items into 3 groups of [10,
        6, 4] items:
         - alpha = [10, 6, 4],
         - alpha = [100, 60, 40],
         - alpha = [5, 3, 2],
        would all work. Variance is inversely related to sum(alpha)
        """
        props = np.concatenate([[0], (scipy.stats.dirichlet.rvs(alpha) * size).cumsum().round().astype(int)])
        indices = np.array(list(range(size)))
        random.shuffle(indices)
        x = []
        for i in range(len(props)-1):
            ix = indices[props[i]:props[i+1]]
            x.append(ix)
        return cls.from_membership(x)

    def num_elements(self):
        return len(self.partition_vector)

    def num_groups(self):
        return max(self.partition_vector) + 1

    @classmethod
    def from_membership(cls, membership):
        tmp = {}
        for group, members in enumerate(membership):
            for member in members:
                tmp[member] = group
        return cls([tmp[k] for k in sorted(tmp)])

    def get_membership(self):
        """
        Alternative representation of group membership -
        creates a list with one tuple per group; each tuple contains
        the indices of its members

        Example:
        partition  = (0,0,0,1,0,1,2,2)
        membership = [(0,1,2,4), (3,5), (6,7)]

        :return: list of tuples giving group memberships by index
        """
        result = defaultdict(list)
        for (position, value) in enumerate(self.partition_vector):
            result[value].append(position)
        return sorted([tuple(x) for x in result.values()])

    @classmethod
    def read(cls, filename):
        with open(filename) as reader:
            s = reader.read()
            t = tuple(s.rstrip().split(','))
        return cls(t)

    def normalised_mutual_information(self, other):
        (entropy_1, entropy_2, mut_inf) = entropies(self, other)

        return 2 * mut_inf / (entropy_1 + entropy_2)

    def variation_of_information(self, other):
        """ calculates Variation of Information Metric between two clusterings
        of the same data - SEE Meila, M. (2007). Comparing clusterings: an
        information based distance. Journal of Multivariate Analysis, 98(5),
        873-895. doi:10.1016/j.jmva.2006.11.013"""
        (entropy_1, entropy_2, mut_inf) = entropies(self, other)

        return entropy_1 + entropy_2 - 2 * mut_inf

    def write(self, filename):
        with open(filename, 'w') as writer:
            writer.write(','.join(str(x) for x in self.partition_vector) + '\n') def enum(*sequential, **named):
    """creates an Enum type with given values"""
    enums = dict(list(zip(sequential, list(range(len(sequential))))), **named)
    enums['reverse'] = dict((value, key) for key, value in list(enums.items()))
    return type('Enum', (object, ), enums) class OptionError(Exception):
    """ Reports on disallowed options passed to functions """

    def __init__(self, option, choices):
        self.value = option
        self.choices = choices

    def __str__(self):
        return '\'{0}\' is not a valid option. Valid options are {1}'.format(self.value,
                                                                             self.choices) def isnumbercheck(n):
    if not isinstance(n, numbers.Number):
        raise ValueError('{0} is not a number.')
    return n def rangecheck(n, lower, upper):
    if n < lower or n > upper:
        raise RangeError(n, lower, upper)
    return n
import logging
import numpy as np
import fastcluster
import sklearn.manifold
import pkg_resources
        import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import fcluster, dendrogram
from scipy.spatial.distance import squareform
    from Bio.Cluster import kmedoids
from sklearn.cluster import AffinityPropagation, DBSCAN, KMeans
    from sklearn.mixture import GaussianMixture
    from sklearn.mixture import GMM
from .distance_matrix import DistanceMatrix, rbf, binsearch_mask, kmask, kscale, affinity, laplace, eigen, \
    double_centre, normalise_rows, CoordinateMatrix
from .partition import Partition
from .utils import enum
from .errors import OptionError, isnumbercheck, rangecheck
        br_bottom = 0
    else:
        br_bottom = linkmat[linkmat_size - nclusters][2]
    threshold = 0.5 * (br_top + br_bottom)
    return threshold

def _hclust(linkmat, nclusters):
    threshold = _get_threshold(linkmat, nclusters)
    t = fcluster(linkmat, threshold, criterion='distance')
    return Partition(t)


class ClusteringManager(object):
    """
    Clustering manager base class
    """
    def __init__(self, dm):
        if isinstance(dm, np.ndarray):
            dm = DistanceMatrix.from_array(dm)

        if not isinstance(dm, DistanceMatrix):
            raise ValueError('Distance matrix should be a numpy array or treeCl.DistanceMatrix')
        self.dm = dm

    def __str__(self):
        return str(self.dm)

    def get_dm(self, noise):
        return self.dm.add_noise().values if noise else self.dm.values


class EMMixin(object):
    """
    Provide methods to do kmeans and GMM estimation
    """
    @staticmethod
    def kmeans(nclusters, coords):
        est = KMeans(n_clusters=nclusters, n_init=50, max_iter=500)
        est.fit(coords)
        return Partition(est.labels_)

    @staticmethod
    def gmm(nclusters, coords, n_init=50, n_iter=500):
        if USE_GAUSSIAN_MIXTURE:
            est = GaussianMixture(n_components=nclusters, n_init=n_init, max_iter=n_iter)
        else:
            est = GMM(n_components=nclusters, n_init=n_init, n_iter=n_iter)
        est.fit(coords)
        return Partition(est.predict(coords))

def _check_val(opt, min_, max_):
    isnumbercheck(opt)
    rangecheck(opt, min_, max_)


class Spectral(ClusteringManager, EMMixin):
    """
    Manager for spectral clustering and Kernel PCA clustering
    """
    def __init__(self, dm, 
                 pruning_option=options.PRUNING_NONE, 
                 scale_option=options.LOCAL_SCALE_MEDIAN, 
                 manual_pruning=None, 
                 manual_scale=None,
                 verbosity=0):
        super(Spectral, self).__init__(dm)
        try:
            options.reverse[pruning_option]
        except KeyError:
            raise OptionError(pruning_option, list(options.reverse.values()))
        try:
            options.reverse[scale_option]
        except KeyError:
            raise OptionError(scale_option, list(options.reverse.values()))

        if pruning_option == options.PRUNING_MANUAL:
            _check_val(manual_pruning, 2, self.dm.df.shape[0])

        if scale_option == options.LOCAL_SCALE_MANUAL:
            _check_val(manual_scale, 2, self.dm.df.shape[0])

        self._pruning_option = pruning_option
        self._scale_option = scale_option
        self._manual_pruning = manual_pruning
        self._manual_scale = manual_scale
        self._verbosity = verbosity
        self._affinity = self.decompose()

    def __str__(self):
        return ('Spectral Clustering with local scaling:\n'
                'Pruning option: {}\n'
                'Scaling option: {}'
                .format(options.reverse[self._pruning_option], options.reverse[self._scale_option]))

    def decompose(self,
            noise=False,
            verbosity=0,
            logic='or',
            **kwargs):
        """ Use prune to remove links between distant points:

        prune is None: no pruning
        prune={int > 0}: prunes links beyond `prune` nearest neighbours
        prune='estimate': searches for the smallest value that retains a fully
        connected graph

        """
        matrix = self.get_dm(noise)

        # get local scale estimate
        est_scale = None

        # ADJUST MASK
        if self._pruning_option == options.PRUNING_NONE:  
            # Set kp to max value
            kp = len(matrix) - 1
            mask = np.ones(matrix.shape, dtype=bool)
        elif self._pruning_option == options.PRUNING_MANUAL:
            # Manually set value of kp
            kp = self._manual_pruning

            mask = kmask(matrix, self._manual_pruning, logic=logic)