scipy/cluster/vq.py

""" K-means Clustering and Vector Quantization Module

    Provides routines for k-means clustering, generating code books
    from k-means models, and quantizing vectors by comparing them with
    centroids in a code book.

    The k-means algorithm takes as input the number of clusters to
    generate, k, and a set of observation vectors to cluster.  It
    returns a set of centroids, one for each of the k clusters.  An
    observation vector is classified with the cluster number or
    centroid index of the centroid closest to it.

    A vector v belongs to cluster i if it is closer to centroid i than
    any other centroids. If v belongs to i, we say centroid i is the
    dominating centroid of v. Common variants of k-means try to
    minimize distortion, which is defined as the sum of the distances
    between each observation vector and its dominating centroid.  Each
    step of the k-means algorithm refines the choices of centroids to
    reduce distortion. The change in distortion is often used as a
    stopping criterion: when the change is lower than a threshold, the
    k-means algorithm is not making sufficient progress and
    terminates.

    Since vector quantization is a natural application for k-means,
    information theory terminology is often used.  The centroid index
    or cluster index is also referred to as a "code" and the table
    mapping codes to centroids and vice versa is often referred as a
    "code book". The result of k-means, a set of centroids, can be
    used to quantize vectors. Quantization aims to find an encoding of
    vectors that reduces the expected distortion.

    For example, suppose we wish to compress a 24-bit color image
    (each pixel is represented by one byte for red, one for blue, and
    one for green) before sending it over the web.  By using a smaller
    8-bit encoding, we can reduce the amount of data by two
    thirds. Ideally, the colors for each of the 256 possible 8-bit
    encoding values should be chosen to minimize distortion of the
    color. Running k-means with k=256 generates a code book of 256
    codes, which fills up all possible 8-bit sequences.  Instead of
    sending a 3-byte value for each pixel, the 8-bit centroid index
    (or code word) of the dominating centroid is transmitted. The code
    book is also sent over the wire so each 8-bit code can be
    translated back to a 24-bit pixel value representation. If the
    image of interest was of an ocean, we would expect many 24-bit
    blues to be represented by 8-bit codes. If it was an image of a
    human face, more flesh tone colors would be represented in the
    code book.

    All routines expect obs to be a M by N array where the rows are
    the observation vectors. The codebook is a k by N array where the
    i'th row is the centroid of code word i. The observation vectors
    and centroids have the same feature dimension.

    whiten(obs) --
        Normalize a group of observations so each feature has unit
        variance.
    vq(obs,code_book) --
        Calculate code book membership of a set of observation
        vectors.
    kmeans(obs,k_or_guess,iter=20,thresh=1e-5) --
        Clusters a set of observation vectors. Learns centroids with
        the k-means algorithm, trying to minimize distortion.  A code
        book is generated that can be used to quantize vectors.
    kmeans2 --
        A different implementation of k-means with more methods for
        initializing centroids.  Uses maximum number of iterations as
        opposed to a distortion threshold as its stopping criterion.

"""
__docformat__ = 'restructuredtext'

__all__ = ['whiten', 'vq', 'kmeans', 'kmeans2']

# TODO:
#   - implements high level method for running several times k-means with
#   different initialialization
#   - warning: what happens if different number of clusters ? For now, emit a
#   warning, but it is not great, because I am not sure it really make sense to
#   succeed in this case (maybe an exception is better ?)
import warnings

from numpy.random import randint
from numpy import shape, zeros, sqrt, argmin, minimum, array, \
     newaxis, arange, compress, equal, common_type, single, double, take, \
     std, mean
import numpy as np

class ClusterError(Exception):
    pass

def whiten(obs):
    """ Normalize a group of observations on a per feature basis.

    Before running k-means, it is beneficial to rescale each feature
    dimension of the observation set with whitening. Each feature is
    divided by its standard deviation across all observations to give
    it unit variance.

    :Parameters:
        obs : ndarray
            Each row of the array is an observation.  The
            columns are the features seen during each observation.
            ::

                      #   f0    f1    f2
                obs = [[  1.,   1.,   1.],  #o0
                       [  2.,   2.,   2.],  #o1
                       [  3.,   3.,   3.],  #o2
                       [  4.,   4.,   4.]]) #o3

            XXX perhaps should have an axis variable here.

    :Returns:
        result : ndarray
            Contains the values in obs scaled by the standard devation
            of each column.

    Examples
    --------

    >>> from numpy import array
    >>> from scipy.cluster.vq import whiten
    >>> features  = array([[  1.9,2.3,1.7],
    ...                    [  1.5,2.5,2.2],
    ...                    [  0.8,0.6,1.7,]])
    >>> whiten(features)
    array([[ 3.41250074,  2.20300046,  5.88897275],
           [ 2.69407953,  2.39456571,  7.62102355],
           [ 1.43684242,  0.57469577,  5.88897275]])

    """
    std_dev = std(obs, axis=0)
    return obs / std_dev

def vq(obs, code_book):
    """ Vector Quantization: assign codes from a code book to observations.

    Assigns a code from a code book to each observation. Each
    observation vector in the M by N obs array is compared with the
    centroids in the code book and assigned the code of the closest
    centroid.

    The features in obs should have unit variance, which can be
    acheived by passing them through the whiten function.  The code
    book can be created with the k-means algorithm or a different
    encoding algorithm.

    :Parameters:
        obs : ndarray
            Each row of the NxM array is an observation.  The columns are the
            "features" seen during each observation. The features must be
            whitened first using the whiten function or something equivalent.
        code_book : ndarray.
            The code book is usually generated using the k-means algorithm.
            Each row of the array holds a different code, and the columns are
            the features of the code.

            ::

                            #   f0    f1    f2   f3
                code_book = [[  1.,   2.,   3.,   4.],  #c0
                             [  1.,   2.,   3.,   4.],  #c1
                             [  1.,   2.,   3.,   4.]]) #c2

    :Returns:
        code : ndarray
            A length N array holding the code book index for each observation.
        dist : ndarray
            The distortion (distance) between the observation and its nearest
            code.

    Notes
    -----
    This currently forces 32-bit math precision for speed.  Anyone know
    of a situation where this undermines the accuracy of the algorithm?

    Examples
    --------
    >>> from numpy import array
    >>> from scipy.cluster.vq import vq
    >>> code_book = array([[1.,1.,1.],
    ...                    [2.,2.,2.]])
    >>> features  = array([[  1.9,2.3,1.7],
    ...                    [  1.5,2.5,2.2],
    ...                    [  0.8,0.6,1.7]])
    >>> vq(features,code_book)
    (array([1, 1, 0],'i'), array([ 0.43588989,  0.73484692,  0.83066239]))

    """
    try:
        import _vq
        ct = common_type(obs, code_book)
        c_obs = obs.astype(ct)
        c_code_book = code_book.astype(ct)
        if ct is single:
            results = _vq.vq(c_obs, c_code_book)
        elif ct is double:
            results = _vq.vq(c_obs, c_code_book)
        else:
            results = py_vq(obs, code_book)
    except ImportError:
        results = py_vq(obs, code_book)
    return results

def py_vq(obs, code_book):
    """ Python version of vq algorithm.

    The algorithm computes the euclidian distance between each
    observation and every frame in the code_book.

    :Parameters:
        obs : ndarray
            Expects a rank 2 array. Each row is one observation.
        code_book : ndarray
            Code book to use. Same format than obs. Should have same number of
            features (eg columns) than obs.

    :Note:
        This function is slower than the C version but works for
        all input types.  If the inputs have the wrong types for the
        C versions of the function, this one is called as a last resort.

        It is about 20 times slower than the C version.

    :Returns:
        code : ndarray
            code[i] gives the label of the ith obversation, that its code is
            code_book[code[i]].
        mind_dist : ndarray
            min_dist[i] gives the distance between the ith observation and its
            corresponding code.

    """
    # n = number of observations
    # d = number of features
    if np.ndim(obs) == 1:
        if not np.ndim(obs) == np.ndim(code_book):
            raise ValueError(
                    "Observation and code_book should have the same rank")
        else:
            return _py_vq_1d(obs, code_book)
    else:
        (n, d) = shape(obs)

    # code books and observations should have same number of features and same
    # shape
    if not np.ndim(obs) == np.ndim(code_book):
        raise ValueError("Observation and code_book should have the same rank")
    elif not d == code_book.shape[1]:
        raise ValueError("Code book(%d) and obs(%d) should have the same " \
                         "number of features (eg columns)""" %
                         (code_book.shape[1], d))

    code = zeros(n, dtype=int)
    min_dist = zeros(n)
    for i in range(n):
        dist = np.sum((obs[i] - code_book) ** 2, 1)
        code[i] = argmin(dist)
        min_dist[i] = dist[code[i]]

    return code, sqrt(min_dist)

def _py_vq_1d(obs, code_book):
    """ Python version of vq algorithm for rank 1 only.

    :Parameters:
        obs : ndarray
            Expects a rank 1 array. Each item is one observation.
        code_book : ndarray
            Code book to use. Same format than obs. Should rank 1 too.

    :Returns:
        code : ndarray
            code[i] gives the label of the ith obversation, that its code is
            code_book[code[i]].
        mind_dist : ndarray
            min_dist[i] gives the distance between the ith observation and its
            corresponding code.

    """
    raise RuntimeError("_py_vq_1d buggy, do not use rank 1 arrays for now")
    n = obs.size
    nc = code_book.size
    dist = np.zeros((n, nc))
    for i in range(nc):
        dist[:, i] = np.sum(obs - code_book[i])
    print dist
    code = argmin(dist)
    min_dist = dist[code]

    return code, sqrt(min_dist)

def py_vq2(obs, code_book):
    """2nd Python version of vq algorithm.

    The algorithm simply computes the euclidian distance between each
    observation and every frame in the code_book/

    :Parameters:
        obs : ndarray
            Expect a rank 2 array. Each row is one observation.
        code_book : ndarray
            Code book to use. Same format than obs. Should have same number of
            features (eg columns) than obs.

    :Note:
        This could be faster when number of codebooks is small, but it
        becomes a real memory hog when codebook is large. It requires
        N by M by O storage where N=number of obs, M = number of
        features, and O = number of codes.

    :Returns:
        code : ndarray
            code[i] gives the label of the ith obversation, that its code is
            code_book[code[i]].
        mind_dist : ndarray
            min_dist[i] gives the distance between the ith observation and its
            corresponding code.

    """
    d = shape(obs)[1]

    # code books and observations should have same number of features
    if not d == code_book.shape[1]:
        raise ValueError("""
            code book(%d) and obs(%d) should have the same
            number of features (eg columns)""" % (code_book.shape[1], d))

    diff = obs[newaxis, :, :] - code_book[:,newaxis,:]
    dist = sqrt(np.sum(diff * diff, -1))
    code = argmin(dist, 0)
    min_dist = minimum.reduce(dist, 0) #the next line I think is equivalent
                                      #  - and should be faster
    #min_dist = choose(code,dist) # but in practice, didn't seem to make
                                  # much difference.
    return code, min_dist

def _kmeans(obs, guess, thresh=1e-5):
    """ "raw" version of k-means.

    :Returns:
        code_book :
            the lowest distortion codebook found.
        avg_dist :
            the average distance a observation is from a code in the book.
            Lower means the code_book matches the data better.

    :SeeAlso:
        - kmeans : wrapper around k-means

    XXX should have an axis variable here.

    Examples
    --------

    Note: not whitened in this example.

    >>> from numpy import array
    >>> from scipy.cluster.vq import _kmeans
    >>> features  = array([[ 1.9,2.3],
    ...                    [ 1.5,2.5],
    ...                    [ 0.8,0.6],
    ...                    [ 0.4,1.8],
    ...                    [ 1.0,1.0]])
    >>> book = array((features[0],features[2]))
    >>> _kmeans(features,book)
    (array([[ 1.7       ,  2.4       ],
           [ 0.73333333,  1.13333333]]), 0.40563916697728591)

    """

    code_book = array(guess, copy = True)
    avg_dist = []
    diff = thresh+1.
    while diff > thresh:
        nc = code_book.shape[0]
        #compute membership and distances between obs and code_book
        obs_code, distort = vq(obs, code_book)
        avg_dist.append(mean(distort, axis=-1))
        #recalc code_book as centroids of associated obs
        if(diff > thresh):
            has_members = []
            for i in arange(nc):
                cell_members = compress(equal(obs_code, i), obs, 0)
                if cell_members.shape[0] > 0:
                    code_book[i] = mean(cell_members, 0)
                    has_members.append(i)
            #remove code_books that didn't have any members
            code_book = take(code_book, has_members, 0)
        if len(avg_dist) > 1:
            diff = avg_dist[-2] - avg_dist[-1]
    #print avg_dist
    return code_book, avg_dist[-1]

def kmeans(obs, k_or_guess, iter=20, thresh=1e-5):
    """Performs k-means on a set of observation vectors forming k
       clusters. This yields a code book mapping centroids to codes
       and vice versa. The k-means algorithm adjusts the centroids
       until sufficient progress cannot be made, i.e. the change in
       distortion since the last iteration is less than some
       threshold.

    :Parameters:
        obs : ndarray
            Each row of the M by N array is an observation vector. The
            columns are the features seen during each observation.
            The features must be whitened first with the whiten
            function.

        k_or_guess : int or ndarray
            The number of centroids to generate. A code is assigned to
            each centroid, which is also the row index of the centroid
            in the code_book matrix generated.

            The initial k centroids are chosen by randomly selecting
            observations from the observation matrix. Alternatively,
            passing a k by N array specifies the initial k centroids.

        iter : int
            The number of times to run k-means, returning the codebook
            with the lowest distortion. This argument is ignored if
            initial centroids are specified with an array for the
            k_or_guess paramter. This parameter does not represent the
            number of iterations of the k-means algorithm.

        thresh : float
            Terminates the k-means algorithm if the change in
            distortion since the last k-means iteration is less than
            thresh.

    :Returns:
        codebook : ndarray
            A k by N array of k centroids. The i'th centroid
            codebook[i] is represented with the code i. The centroids
            and codes generated represent the lowest distortion seen,
            not necessarily the globally minimal distortion.

        distortion : float
           The distortion between the observations passed and the
           centroids generated.

    :SeeAlso:
        - kmeans2: a different implementation of k-means clustering
          with more methods for generating initial centroids but without
          using a distortion change threshold as a stopping criterion.
        - whiten: must be called prior to passing an observation matrix
          to kmeans.

    Examples
    --------

    >>> from numpy import array
    >>> from scipy.cluster.vq import vq, kmeans, whiten
    >>> features  = array([[ 1.9,2.3],
    ...                    [ 1.5,2.5],
    ...                    [ 0.8,0.6],
    ...                    [ 0.4,1.8],
    ...                    [ 0.1,0.1],
    ...                    [ 0.2,1.8],
    ...                    [ 2.0,0.5],
    ...                    [ 0.3,1.5],
    ...                    [ 1.0,1.0]])
    >>> whitened = whiten(features)
    >>> book = array((whitened[0],whitened[2]))
    >>> kmeans(whitened,book)
    (array([[ 2.3110306 ,  2.86287398],
           [ 0.93218041,  1.24398691]]), 0.85684700941625547)

    >>> from numpy import random
    >>> random.seed((1000,2000))
    >>> codes = 3
    >>> kmeans(whitened,codes)
    (array([[ 2.3110306 ,  2.86287398],
           [ 1.32544402,  0.65607529],
           [ 0.40782893,  2.02786907]]), 0.5196582527686241)

    """
    if int(iter) < 1:
        raise ValueError, 'iter must be >= to 1.'
    if type(k_or_guess) == type(array([])):
        guess = k_or_guess
        if guess.size < 1:
            raise ValueError("Asked for 0 cluster ? initial book was %s" % \
                             guess)
        result = _kmeans(obs, guess, thresh = thresh)
    else:
        #initialize best distance value to a large value
        best_dist = 100000
        No = obs.shape[0]
        k = k_or_guess
        if k < 1:
            raise ValueError("Asked for 0 cluster ? ")
        for i in range(iter):
            #the intial code book is randomly selected from observations
            guess = take(obs, randint(0, No, k), 0)
            book, dist = _kmeans(obs, guess, thresh = thresh)
            if dist < best_dist:
                best_book = book
                best_dist = dist
        result = best_book, best_dist
    return result

def _kpoints(data, k):
    """Pick k points at random in data (one row = one observation).

    This is done by taking the k first values of a random permutation of 1..N
    where N is the number of observation.

    :Parameters:
        data : ndarray
            Expect a rank 1 or 2 array. Rank 1 are assumed to describe one
            dimensional data, rank 2 multidimensional data, in which case one
            row is one observation.
        k : int
            Number of samples to generate.

    """
    if data.ndim > 1:
        n = data.shape[0]
    else:
        n = data.size

    p = np.random.permutation(n)
    x = data[p[:k], :].copy()

    return x

def _krandinit(data, k):
    """Returns k samples of a random variable which parameters depend on data.

    More precisely, it returns k observations sampled from a Gaussian random
    variable which mean and covariances are the one estimated from data.

    :Parameters:
        data : ndarray
            Expect a rank 1 or 2 array. Rank 1 are assumed to describe one
            dimensional data, rank 2 multidimensional data, in which case one
            row is one observation.
        k : int
            Number of samples to generate.

    """
    def init_rank1(data):
        mu  = np.mean(data)
        cov = np.cov(data)
        x = np.random.randn(k)
        x *= np.sqrt(cov)
        x += mu
        return x
    def init_rankn(data):
        mu  = np.mean(data, 0)
        cov = np.atleast_2d(np.cov(data, rowvar = 0))

        # k rows, d cols (one row = one obs)
        # Generate k sample of a random variable ~ Gaussian(mu, cov)
        x = np.random.randn(k, mu.size)
        x = np.dot(x, np.linalg.cholesky(cov).T) + mu
        return x

    nd = np.ndim(data)
    if nd == 1:
        return init_rank1(data)
    else:
        return init_rankn(data)

_valid_init_meth = {'random': _krandinit, 'points': _kpoints}

def _missing_warn():
    """Print a warning when called."""
    warnings.warn("One of the clusters is empty. "
                 "Re-run kmean with a different initialization.")

def _missing_raise():
    """raise a ClusterError when called."""
    raise ClusterError, "One of the clusters is empty. "\
                        "Re-run kmean with a different initialization."

_valid_miss_meth = {'warn': _missing_warn, 'raise': _missing_raise}

def kmeans2(data, k, iter = 10, thresh = 1e-5, minit = 'random',
        missing = 'warn'):
    """Classify a set of observations into k clusters using the k-means
       algorithm.

    The algorithm attempts to minimize the Euclidian distance between
    observations and centroids. Several initialization methods are
    included.

    :Parameters:
        data : ndarray
            A M by N array of M observations in N dimensions or a length
            M array of M one-dimensional observations.
        k : int or ndarray
            The number of clusters to form as well as the number of
            centroids to generate. If minit initialization string is
            'matrix', or if a ndarray is given instead, it is
            interpreted as initial cluster to use instead.
        iter : int
            Number of iterations of the k-means algrithm to run. Note
            that this differs in meaning from the iters parameter to
            the kmeans function.
        thresh : float
            (not used yet).
        minit : string
            Method for initialization. Available methods are 'random',
            'points', 'uniform', and 'matrix':

            'random': generate k centroids from a Gaussian with mean and
            variance estimated from the data.

            'points': choose k observations (rows) at random from data for
            the initial centroids.

            'uniform': generate k observations from the data from a uniform
            distribution defined by the data set (unsupported).

            'matrix': interpret the k parameter as a k by M (or length k
            array for one-dimensional data) array of initial centroids.

    :Returns:
        centroid : ndarray
            A k by N array of centroids found at the last iteration of
            k-means.
        label : ndarray
            label[i] is the code or index of the centroid the
            i'th observation is closest to.
    """
    if missing not in _valid_miss_meth.keys():
        raise ValueError("Unkown missing method: %s" % str(missing))
    # If data is rank 1, then we have 1 dimension problem.
    nd  = np.ndim(data)
    if nd == 1:
        d = 1
        #raise ValueError("Input of rank 1 not supported yet")
    elif nd == 2:
        d = data.shape[1]
    else:
        raise ValueError("Input of rank > 2 not supported")

    if np.size(data) < 1:
        raise ValueError("Input has 0 items.")

    # If k is not a single value, then it should be compatible with data's
    # shape
    if np.size(k) > 1 or minit == 'matrix':
        if not nd == np.ndim(k):
            raise ValueError("k is not an int and has not same rank than data")
        if d == 1:
            nc = len(k)
        else:
            (nc, dc) = k.shape
            if not dc == d:
                raise ValueError("k is not an int and has not same rank than\
                        data")
        clusters = k.copy()
    else:
        try:
            nc = int(k)
        except TypeError:
            raise ValueError("k (%s) could not be converted to an integer " % str(k))

        if nc < 1:
            raise ValueError("kmeans2 for 0 clusters ? (k was %s)" % str(k))

        if not nc == k:
            warnings.warn("k was not an integer, was converted.")
        try:
            init = _valid_init_meth[minit]
        except KeyError:
            raise ValueError("unknown init method %s" % str(minit))
        clusters = init(data, k)

    assert not iter == 0
    return _kmeans2(data, clusters, iter, nc, _valid_miss_meth[missing])

def _kmeans2(data, code, niter, nc, missing):
    """ "raw" version of kmeans2. Do not use directly.

    Run k-means with a given initial codebook.  """
    for i in range(niter):
        # Compute the nearest neighbour for each obs
        # using the current code book
        label = vq(data, code)[0]
        # Update the code by computing centroids using the new code book
        for j in range(nc):
            mbs = np.where(label==j)
            if mbs[0].size > 0:
                code[j] = np.mean(data[mbs], axis=0)
            else:
                missing()

    return code, label

if __name__  == '__main__':
    pass
    #import _vq
    #a = np.random.randn(4, 2)
    #b = np.random.randn(2, 2)

    #print _vq.vq(a, b)
    #print _vq.vq(np.array([[1], [2], [3], [4], [5], [6.]]),
    #             np.array([[2.], [5.]]))
    #print _vq.vq(np.array([1, 2, 3, 4, 5, 6.]), np.array([2., 5.]))
    #_vq.vq(a.astype(np.float32), b.astype(np.float32))
    #_vq.vq(a, b.astype(np.float32))
    #_vq.vq([0], b)