try: 
    import psyco
    psyco.full()
except ImportError:
    print "Impossible to import psyco library"
    pass

try:
    import numpy
    import numpy.linalg
    import scipy
except ImportError:
    raise ImportError, "Impossible to import numpy/scipy module"

import math
import random
import copy
import itertools
import operator
import heapq

try:
    import statfuncs
except ImportError:
    raise ImportError, "Impossible to import statfuncs module"




__author__ = "Andrea Maffezzoli"
__version__ = "$Revision: 4.98 - 28 Oct 2010$"
__name__ = "kmeans"
__all__ = ["minsmaxs", "Euclide_distance", \
           "scatter_old", "_cov_tot", "_cov_between", \
           "_cov_within", "_cov_within_no_samplesnr", \
           "scattertrace", "scatterbetween", "scatterwilks", \
           "scatterwithin_no_samplesnr", "scatterwithin", \
           "kmeans", "fuzzy_kmeans", "kmeans_replica"]



def minsmaxs(L):
    """This function returns the min and max of each variable for samples in L"""

    if not(hasattr(L[0], '__iter__')): return [(min(L), max(L))]
    return [(i[0], i[-1]) for i in map(sorted, map(None, *L))]


def Euclide_distance(m, p):
    """This function calculates the Euclidean distance between the samples m and p
        For other techniques to compute the distance see:
        http://docs.scipy.org/doc/scipy/reference/spatial.distance.html
    """
    
    _pow = math.pow
    #################################################################################
    # When working with only one sample of one variable or the I Principal Component
    #################################################################################
    if isinstance(p, int) | isinstance(p, float): return math.fabs(m - p)
    elif hasattr(p, '__iter__'):
        #####################################################################################
        # When p comes from only the I Principal Component, whereas m is a list of variables 
        #####################################################################################
        if len(p) == 1:
            if hasattr(m, '__iter__'): return math.fabs(m[0] - p[0])
            return math.fabs(m - p[0]) 
        ###########################################
        # When working with more than one variable
        ###########################################
        return math.sqrt(sum((i-j)**2 for i, j in zip(m, p)))


def scatter_old(M, clusters_indexes, distance = Euclide_distance):

    k = len(set(clusters_indexes))
    if k <= 1: raise Exception, "The number of clusters has to be at least 2"
    clusters_data = [[i for i in M if clusters_indexes[M.index(i)] == z] \
                      for z in set(clusters_indexes)]
    mnval, centr = statfuncs.meanval, statfuncs.centroid
    sct = scatter
    scatter_within = mnval([sct(c) for c in clusters_data])
    scatter_between = statfuncs.meanval([statfuncs.meanval([distance(centr(clusters_data[z]), \
                            centr(clusters_data[i])) for i in range(z + 1, k)]) for z in range(k - 1)]) 
    return scatter_within / scatter_between # / (scatter_within + scatter_between)


def _cov_between(W):

    tuples_centroids = [(len(S), numpy.mean(S, axis = 0)) for S in W]
    global_centroid = numpy.mean(numpy.concatenate(W), axis = 0)
    ##########################################################################
    # Don't use numpy.dot to calculate _cov_between_temp, but use numpy.outer
    ##########################################################################
    _cov_between_temp =  sum([c[0] * numpy.outer(c[1] - global_centroid, \
                                                 numpy.transpose(c[1] - global_centroid)) \
                              for c in tuples_centroids])
    if len(W[0][0]) == 1: return _cov_between_temp[0][0]
    else: return _cov_between_temp


#####################################################################################
# For the following functions '_cov_within', '_cov_within_no_samplesnr', '_cov_tot' 
# use option 'bias=1' is necessary to give numpy.cov for clusters of only one sample
#####################################################################################
def _cov_within(W):

    return sum([len(S) * numpy.cov(S, rowvar = 0, bias = 1) for S in W])


def _cov_within_no_samplesnr(W):
    """See also equation (2) in reference:
       Chris Ding and Xiaofeng He. "K-means Clustering via Principal Component Analysis".
       Proc. of Int'l Conf. Machine Learning (ICML 2004), pp 225-232. July 2004.
    """

    _covariance = lambda L: numpy.cov(L, rowvar = 0, bias = 1)
    return sum(map(_covariance, W))


def _cov_tot(W):

    W_total = numpy.concatenate(W)
    return len(W_total) * numpy.cov(W_total, rowvar = 0, bias = 1)


def scattertrace(W):

    return numpy.trace(_cov_within(W))


def scatterdet(W):

    return numpy.linalg.det(_cov_within(W))


def scatterbetween(W):

    #########################
    # With only one variable
    #########################
    if len(W[0][0]) == 1: return numpy.cov(map(numpy.concatenate, W), rowvar = 0, bias = 1)
    else: return numpy.linalg.det(_cov_between(W))


def scatterwithin_no_samplesnr(W):

    #########################
    # With only one variable
    #########################
    if len(W[0][0]) == 1: return sum(map(numpy.var, map(numpy.concatenate, W)))
    #    return sum(numpy.cov(m, bias=1) for m in map(numpy.concatenate, W))
    #    return sum(numpy.cov(m) if len(m) > 1 else 0. for m in map(numpy.concatenate, W))
    _cw = _cov_within_no_samplesnr(W)
    if not(getattr(_cw, 'shape', False)): _cw = numpy.array([_cw])
    if _cw.shape in ((1,1), (1,)): return 0.
    return numpy.linalg.det(_cw)
    
def scatterwithin(W):

    #########################
    # With only one variable
    #########################
    if len(W[0][0]) == 1: return numpy.cov(map(numpy.concatenate, W), bias = 1)
    else:
        __cov_within = _cov_within(W)
        try:
            return numpy.linalg.det(__cov_within)
        #####################################################################
        # When W is only one sample and _cov_within(W) is only a float value
        #####################################################################
        except numpy.linalg.linalg.LinAlgError:
            print 'Sum of Within Covariance matrices is not an array: it is equal to', __cov_within
            return __cov_within


def scatterwilks(W):

    ############################################################
    # When using only one variable or the I Principal Component
    ############################################################
    if len(W[0][0]) == 1:
        return scatterwithin(map(numpy.concatenate, W)) / _cov_tot(W)
    else:
        try:
            return scatterwithin(W) / numpy.linalg.det(_cov_tot(W))
        except ZeroDivisionError:
            print 'Covariance matrix\n', _cov_tot(W), '\nof input data\n', W, '\nis singular!!!'
            print 'NB: Covariance matrix is singular when "samples nr.<=variables nr"!!!\n'
            raise



######################################################################################
##########                                                                  ##########
##########                              K-MEANS                             ##########
##########                                                                  ##########
######################################################################################
class kmeans(object):
    """This class performs the k-means algorithm.
       The reference is the book:
       - R.O. Duda, P.E. Hart, and D.G. Stork. Pattern Classification. 
         Wiley-Interscience, October 2000 (II edition).
       As authors say, it can be considered an iterative optimization clustering procedure. 
       It is useful to simplify the performing of clustering computations, e.g respect to
       Normal mixtures method, specially in the case when multivariate input data are not
       pairwisely correlated and with approx. the same variance on each variable.
       Here the means iteratively move in order to minimize a squared-error criterion: 
       it can be viewed as an approximate way to get maximum-likelihood estimated for the
       means.
       The matrix in input has as many lists as the number of samples, each one of such 
       lists with a length equal to the nr of variables/dimensions of the problem.
       Input data matrix D has to be a matrix of elements, such algorithm doesn't work 
       with lists.
       In the following functions input data W or D (or D.values() if D is a dictionary) 
       is shaped as in the following:
        D=[
            [
                [var1, var2, ...var-mth] for sample 1, 
                [var1, var2, ...var-mth] for sample 2, 
                [...],
                [var1, var2, ...var-mth] for sample n, 
            ] for cluster 1,
            [
                [var1, var2, ...var-mth] for sample 1, 
                [...]
            ] for cluster 2,
            [   
                [var1, var2, ...var-mth] for sample 1, 
                [...]
            ] for cluster m-th,
        ]
    """
    
    def __init__(self, D, clusters_nr,
                 start_centroids = {}, 
                 m_init = 'uniform',
                 distance = Euclide_distance,
                 scaling_down = True,
                 scatter_calc = True,
                 scatter_method = scatterwithin_no_samplesnr,
                 max_nr_iter = 10**3):
        """m_init is the method for clusters centroids initialization:
            - uniform if using a uniform pdf
            - gauss if using a Gaussian centered on the all samples matrix data centroid
              and with standard deviation equal to 2
            - samples if using a number of data samples picked at random
        """

        M = dict(enumerate(D)) if not(isinstance(D, dict)) else D
        if not(M.itervalues()): raise Exception, "Data matrix can't be empty!"
        centr = statfuncs.centroid
        if clusters_nr > len(M) and not(scaling_down): 
            print """Please give attention:
the nr of clusters > nr of samples and also parameter 'scaling_down' is False.
So now switching automatically the value of parameter 'scaling_down' to True. 
"""
            scaling_down = True

        #########################
        # Initializing centroids
        #########################
        if not(start_centroids):
            if m_init is 'uniform':
                mins_maxs = minsmaxs(M.values())
                centroids = dict((p, [random.uniform(m[0], m[1]) for m in mins_maxs]) for p in xrange(clusters_nr))
            elif m_init is 'gauss':
                global_centroid = centr(M.values())
                self.centroids = dict((p, [random.gauss(m, 2) for m in global_centroid]) for p in xrange(clusters_nr))
            elif m_init is 'samples':
                self.centroids = dict(enumerate(random.sample(M.itervalues(), clusters_nr)))
            else:
                print 'Bad centroids distribution initializer!'
                raise
        else: centroids = start_centroids

        c = dict((j , []) for j in centroids)
        self.nr_iteration = 0
        self.centroids_path = [centroids]

        ######################################################################
        # Iterating Lloyd's algorithm and initializing the clusters centroids
        ######################################################################
        while (not(all([one_centr in centroids.itervalues() for one_centr in c.itervalues()])) \
          or not(all([one_centr in c.itervalues() for one_centr in centroids.itervalues()]))) \
          and self.nr_iteration < max_nr_iter:
            
            self.nr_iteration += 1
            ###################################################################
            # The dict samples_for_cluster ties each neuron to a cluster index
            ###################################################################
            c = copy.copy(centroids)
            self.samples_for_cluster = {}
            for i, j in M.iteritems():
                dist_tuples = []
                for p, q in c.iteritems(): heapq.heappush(dist_tuples, (distance(j, q), p))
                self.samples_for_cluster.setdefault(heapq.heappop(dist_tuples)[1], []).append(i)
            clusters_data = dict((p, [M[i] for i in q]) for p, q in self.samples_for_cluster.iteritems())
            
            #####################################################################################
            # Updating clusters centroids.
            # If scaling down is permitted, centroid are calculated only for non empty clusters,
            # as those given by self.clusters_data, so the ordinal ID of empty clusters will 
            # not enter as key of self.centroids dict.
            # If scaling down is permitted, centroid of empty clusters are calculated ex-novo.
            # based on m_init inizialization options.
            #####################################################################################
            centroids = dict((p, centr(q)) for p, q in clusters_data.iteritems())
            if not scaling_down: 
                for k in xrange(clusters_nr):
                    if k not in self.samples_for_cluster:
                        if m_init is 'uniform':   centroids[k] = [random.uniform(m[0], m[1]) for m in mins_maxs]
                        elif m_init is 'gauss':   centroids[k] = [random.gauss(m, 2) for m in global_centroid]
                        elif m_init is 'samples': centroids[k] = random.sample(M.values(), 1)
            
            self.centroids_path.append(centroids)
        
        self.centroids = centroids
        if scatter_calc:
            #################################################################################
            # self.wilks_lambda contains the scattering values for the clustering performed.
            # Exceptions happens when e.g. a cluster contains only one data sample.
            #################################################################################
            try: self.wilks_lambda = scatter_method(clusters_data.values())
            except: raise

    def __str__(self): return 'K-means'
            
    def __repr__(self): return self.__str__() 


#######################################################################################
##########                                                                   ##########
##########                          K-MEANS REPLICA                          ##########
##########                                                                   ##########
#######################################################################################
class kmeans_replica(object):
    """Here iterating the kmeans performing for replicates_nr times and with the
       same clusters nr., as used also in Matlab, in order to achieve an optimal
       solution and to overcome local minima/maxima problems.
       Covariance matrix is singular when   nr.samples <= nr.variables.
       
       NB: - When: nr.clusters <= nr.variables
             or: nr.samples - nr.clusters < nr.variables
             the between-cluster covariance matrix is singular. 
           - If leaving unchanged 'Mahalanobis_threshold', better to set algorithm=kmeans,
             because fuzzy_kmeans shall change the clusters indexes, given by keys of 
             dicts self.clusters_data or self.centroids, converging very difficultly
             to a final solution!
           - fuzzy_kmeans change the final solution at every iterations,
             in general different from clusters centroids
             (calculated after 'Mahalanobis_threshold' based selection), 
             so Outliers Detection doesn't converge to a final unique solution. 
           - If setting algorithm=fuzzy_kmeans, then use an high 'Mahalanobis_threshold'".
       
       cluster_sorted_by_samples_nr = =1,0 or False,-1 for decrescent, none, crescent
       sorting of clusters based on the samples nr.
    """

    def __init__(self, D, data_name = None, \
                 min_clusters_nr = 2, \
                 max_clusters_nr = 6, \
                 normalization = False, \
                 replicates_nr = 100, \
                 outliers_detect = False, \
                 Mahalanobis_threshold = 2.0, \
                 move_poor_cluster_2_Outliers = False, \
                 cluster_sorted_by_samples_nr = 1, \
                 verbose = True,
                 max_nr_iter = 10**3,
                 **kwargs):
        
        ##################################################################################
        # Transforming D in dictionary M and applying data normalization if required
        ##################################################################################
        M = dict(enumerate(D)) if not(isinstance(D, dict)) else D
        self.data_name = filter(bool, (data_name, map(repr, xrange(len(D)))))[0]
        if normalization: M = dict(zip(M.keys(), statfuncs.normalize(M.values())))
        var_nr = len(D[0])
        if kwargs['scaling_down']: print 'Scaling_down parameter is on.' 

        #########################################################
        # Inspecting the goodness-of-fit for various clusters nr
        #########################################################
        bestclustering = kmeans(M, min_clusters_nr, **kwargs)
        scattering = 1. / filter(bool, (bestclustering.wilks_lambda, 10^-6))[0]
        for k in xrange(min_clusters_nr, max_clusters_nr + 1):
            ###################################################################
            # Inspecting the goodness-of-fit for various clustering iterations
            ###################################################################
            all_wilks = []
            for r in xrange(replicates_nr):
                clust_temp = kmeans(M, k, **kwargs)
                scattering_temp = 1. / filter(bool, (clust_temp.wilks_lambda, 10^-6))[0]
                all_wilks.append((clust_temp.wilks_lambda, len(clust_temp.samples_for_cluster)))
                ####################################################################
                # The 2nd condition aims to avoid degenarate cases 
                ####################################################################
                if scattering_temp > scattering and len(clust_temp.samples_for_cluster) != 1: 
                    bestclustering, scattering = clust_temp, scattering_temp
            ### print "Best goodness-of-fit for %d clusters is %f with %d clusters" \
            ###        % (k, min(all_wilks)[0], min(all_wilks)[1])
        
        self.centroids_path_total = bestclustering.centroids_path
        self.clusters_nr = len(bestclustering.samples_for_cluster)
        #if verbose:
        #    print "%s algorithm found %d clusters after %d iterations." \
        #           % (bestclustering, self.clusters_nr, bestclustering.nr_iteration)
        
        ################################################
        # Mahalanobis distance based Outliers Detection 
        ################################################
        if outliers_detect:
            #if verbose: print 'Detecting Outliers...'
            self.all_outliers_x_cluster = {}
            self.Mahalanobis_dist = {}
            mean_Mahalanobis_dist = {}
            std_Mahalanobis_dist = {}
            centr = statfuncs.centroid
            self.nr_iteration_outliers = 0
            self.centroids_no_Outliers = dict((p, [random.random() for m in D[0]]) for p in bestclustering.centroids) 
            
            ##############################################################################
            #####   N.B.: The following while exit condition holds for k-means!!!    #####
            ##############################################################################
            nr_iter = 0
            while set(statfuncs.tuplize(bestclustering.centroids.values())).symmetric_difference(set(statfuncs.tuplize(self.centroids_no_Outliers.values()))) \
                  and nr_iter < max_nr_iter:
            ###########################################################################################
            #####   N.B.: The following while exit condition is suitable for fuzzy k-means!!!     #####
            ###########################################################################################
            #while not(self.nr_iteration_outliers) or \
            # set(statfuncs.tuplize(self.centroids_path_total[-1].values())).symmetric_difference(set(statfuncs.tuplize(self.centroids_path_total[-2].values()))):
                
                nr_iter += 1
                for cluster, cluster_neurons in bestclustering.samples_for_cluster.iteritems():

                    cluster_centroid = bestclustering.centroids[cluster]
                    
                    ##################################################################
                    # Calculating Mahalanobis distances and clusters without Outliers
                    ##################################################################
                    try:
                        #################################
                        # If data have only one variable 
                        #################################
                        if var_nr == 1:
                            cluster_inv_cov = 1./numpy.var([M[i] for i in cluster_neurons])
                            self.Mahalanobis_dist[cluster] = dict((i, (M[i][0] - cluster_centroid[0])**2 / cluster_inv_cov) \
                                                                   for i in cluster_neurons)
                            #cluster_cov = numpy.var([M[i] for i in cluster_neurons])
                            #self.Mahalanobis_dist[cluster] = dict((i, cluster_cov*(M[i][0]-cluster_centroid[0])**2) for i in cluster_neurons)
                        ######################################
                        # If data have more than one variable 
                        ######################################
                        else:
                            cluster_inv_cov_matrix = numpy.linalg.inv(numpy.cov([M[i] for i in cluster_neurons], rowvar = 0, bias=1))
                            self.Mahalanobis_dist[cluster] = dict([(i, \
                                                numpy.dot(numpy.dot(numpy.array(M[i]) - numpy.array(cluster_centroid), \
                                                                    cluster_inv_cov_matrix), 
                                                          numpy.transpose(numpy.array(M[i]) - numpy.array(cluster_centroid)))) \
                                                for i in cluster_neurons])
                        mean_Mahalanobis_dist[cluster] = statfuncs.meanval(self.Mahalanobis_dist[cluster].values())
                        std_Mahalanobis_dist[cluster]  = statfuncs.standdev(self.Mahalanobis_dist[cluster].values())
                        self.all_outliers_x_cluster.setdefault(cluster, []).extend(filter(lambda i: 
                                                    math.fabs(self.Mahalanobis_dist[cluster][i] - mean_Mahalanobis_dist[cluster]) > \
                                                    Mahalanobis_threshold * std_Mahalanobis_dist[cluster], cluster_neurons))
                    ###################################################
                    # When there is only one sample in the cluster
                    # take it if move_poor_cluster_2_Outliers is False 
                    ###################################################
                    #except IndexError:
                    except:
                        if move_poor_cluster_2_Outliers: self.all_outliers_x_cluster[cluster] += cluster_neurons
                        continue
                    ############################################################################
                    # When covariance matrix is singular, e.g because samples nr.<variables nr.
                    ############################################################################
                    #except numpy.linalg.linalg.LinAlgError:
                    #    if move_poor_cluster_2_Outliers: self.all_outliers_x_cluster[cluster] += cluster_neurons
                    #    pass

                self.all_Outliers = numpy.asarray(statfuncs.deep_chain_new(self.all_outliers_x_cluster.values()))
                clusters_data_no_Outliers   = dict((c, [M[j] for j in d if j not in self.all_Outliers]) \
                                                    for c, d in bestclustering.samples_for_cluster.items())
                self.centroids_no_Outliers  = dict((p, centr(q)) for p, q in clusters_data_no_Outliers.items() if q)
                self.centroids_path_total.append(self.centroids_no_Outliers)
                M = dict((p,q) for p,q in M.iteritems() if p not in self.all_Outliers)
        
                clusters_nr_no_Outliers = len(self.centroids_no_Outliers)
                if verbose and clusters_nr_no_Outliers < self.clusters_nr:
                    print 'Scaling-down clusters number from %s to %s' % (self.clusters_nr, len(self.centroids_no_Outliers))
                self.clusters_nr = clusters_nr_no_Outliers

                #print 'Iter nr:\t', self.nr_iteration_outliers
                #print 'Outliers nr:\t', len(self.all_Outliers)
                #print 'Outliers:\t', self.all_outliers_x_cluster
                #print 'Clust no Out centr:\t', self.centroids_no_Outliers 
                
                ###################################################################
                # Inspecting the goodness-of-fit for various clustering iterations
                # kwargs_Outliers equal kwargs, but now when finding Outliers 
                # 'scaling_down' is set True by default in dict kwargs_Outliers.
                ###################################################################
                scattering_index = 0.0
                kwargs_Outliers = kwargs
                kwargs_Outliers['scaling_down'] = True
                for r in xrange(replicates_nr):
                    clust_temp = kmeans(M, self.clusters_nr, \
                                        start_centroids = self.centroids_no_Outliers, \
                                        **kwargs_Outliers)
                    scattering_index_temp = 1. / filter(bool, (clust_temp.wilks_lambda, 10^-6))[0]
                    
                    ######################################################################
                    # The following condition 'all_wilks[-1][0] is not 1' wants to avoid 
                    # the degenerate case of 1 cluster found via Wilks_lambda index.
                    # The condition 'len(clust_temp.samples_for_cluster.keys()) is not 1' 
                    # works when 'scaling_down' is True, because it is possible to obtain
                    # only one when searching two or more clusters. 
                    ######################################################################
                    if scattering_index_temp > scattering_index and \
                       len(clust_temp.samples_for_cluster) is not 1:
                        bestclustering, scattering_index = clust_temp, scattering_index_temp
                
                self.nr_iteration_outliers += 1

        #########################################################################
        # Storing in the following: 
        # - in clusters_names the names of neurons divided for each cluster,
        # - in clusters_centroids the centroids of each cluster,
        # - in cluster_nr the number k of clusters giving the lesser scattering,
        # - self.centroid_plus_sem is a dictionary with cluster index as key 
        #   cluster centroids with SEM as related value
        #########################################################################   
        self.clusters_nr        = len(bestclustering.samples_for_cluster)
        self.clusters_centroids = self.centroids = bestclustering.centroids
        self.clustering         = bestclustering
        self.centroids_plus_sem = dict((p, statfuncs.centroid_plus_sem([M.get(i, []) for i in q])) \
                                       for p, q in bestclustering.samples_for_cluster.iteritems())
        if verbose:
            print "Outliers Detections alg. found %d outliers in %d clusters after %d iterations." \
                   % (len(getattr(self, "all_Outliers", ())), self.clusters_nr, \
                      getattr(self, "nr_iteration_outliers", 0))
        
    def get_clusters_names(self):

        return dict((p, [self.data_name[i] for i in q]) \
                     for p, q in self.clustering.samples_for_cluster.iteritems() if q)
    
    def get_samples_for_cluster(self): return self.clustering.samples_for_cluster
   
    def get_nr_samples_for_cluster(self):
        
        return tuple(map(len, self.clustering.samples_for_cluster.itervalues()))
    
    def get_all_outliers_names(self):
        
        if not(hasattr(self, "all_outliers_x_cluster")): return {}
        return dict((p, [self.data_name[i] for i in q]) \
                    for p, q in self.all_outliers_x_cluster.iteritems())

    clusters_names          = property(get_clusters_names)
    samples_for_cluster     = property(get_samples_for_cluster)
    nr_samples_for_cluster  = property(get_nr_samples_for_cluster)
    all_outliers_names      = property(get_all_outliers_names)


"""
#####################################################################
##########                                                 ##########
##########          Testing clustering algorithms          ##########
##########                                                 ##########
#####################################################################
if __name__ == "kmeans":

    import pylab
    import pca
    import time

    pdf = random.gauss #expovariate
    
    M = [[pdf(1, 1.5) for j in range(2)] for i in range(100)] + \
        [[pdf(15, 1.5) for i in range(2)] for i in range(100)] 
        [[pdf(1.5, 1.5), pdf(15, 1.5)] for i in range(100)] + \
        [[pdf(14.5, 1.5), pdf(1, 1.5)] for i in range(100)] + \
        [[pdf(7, 1), pdf(7, 1)] for i in range(100)] + \
        [[pdf(20, 1), pdf(7, 1)] for i in range(100)] + \
        [[pdf(-7, 1), pdf(7.5, 1)] for i in range(100)] + \
        [[pdf(7, 1.5), pdf(20, 1.5)] for i in range(100)] + \
        [[pdf(7, 1), pdf(-7, 1)] for i in range(100)] + \
        [[random.uniform(-8, 21) for j in range(2)] for i in range(50)]
    kmeans_dict = {'min_clusters_nr':       8, \
                   'max_clusters_nr':       12, \
                   'replicates_nr':         50, \
                   'normalization':         False, \
                   'algorithm':             kmeans, \
                   'scatter_method':        scatterwithin_no_samplesnr, \
                   'outliers_detect':       True, \
                   'Mahalanobis_threshold': 1.5, \
                   'scaling_down':          True}
    
    start = time.clock()
    clustering = kmeans_replica(M, **kmeans_dict)
    print '>>> Time required to', str(kmeans_dict['algorithm']), 'algorithm:', time.clock() - start, 's.'

    colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', \
              '#880000', '#000088', '#008800', \
              '#AA0000', '#00AA00', '#0000AA']
    pylab.figure(1)
    pylab.axis('equal')
    pylab.grid(True)

    for K in clustering.clusters_data:
        if clustering.clusters_data[K]:
            _color = colors[clustering.clusters_data.keys().index(K)]
            S = zip(*clustering.clusters_data[K])
            #pylab.subplot(1, 2, 1)
            pylab.plot(S[0], S[1], color = _color, marker = '.', linestyle = 'None')
        
            one_centr_coord = zip(*[cp[K] for cp in clustering.centroids_path_total])
            pylab.plot(one_centr_coord[0], one_centr_coord[1], color = _color, marker = 'o', markersize = 3, linestyle = ':')
        
            if kmeans_dict['outliers_detect'] and clustering.all_outliers_x_cluster[K]:
                Out = zip(*[M[i] for i in clustering.all_outliers_x_cluster[K]])
                pylab.plot(Out[0], Out[1], color = _color, marker = 'x', linestyle = 'None')
        else: print 'No data in cluster', K


    print "Clustering on Principal Components.."
    P = pca.pca(M)
    scores = P.SCORE
    
    # Only if using the first Principal Component
    #scores = map(list, zip(map(operator.itemgetter(0), P.SCORE)))

    start = time.clock()
    clustering = kmeans_replica(scores, **kmeans_dict)
    print '>>> Time required to', str(kmeans_dict['algorithm']), 'algorithm:', time.clock() - start, 's.'
    
    pylab.figure(2)
    pylab.axis('equal')
    pylab.grid(True)
    
    for K in clustering.clusters_data:
        if clustering.clusters_data[K]:
            _color = colors[clustering.clusters_data.keys().index(K)]
            S = zip(*clustering.clusters_data[K])
            #pylab.subplot(1, 2, 2)
            pylab.plot(S[0], S[1], color = _color, marker = '.', linestyle = 'None')
            
            # Only if using the first Principal Component
            #pylab.plot(S[0], color = _color, marker = '.', linestyle = 'None')
                
            one_centr_coord = zip(*[cp[K] for cp in clustering.centroids_path_total])
            pylab.plot(one_centr_coord[0], one_centr_coord[1], color = _color, marker = 'o', markersize = 3, linestyle = ':')
            
            # Only if using the first Principal Component
            #pylab.plot(one_centr_coord[0], color = _color, marker = 'o', markersize = 3, linestyle = ':')
        
            if kmeans_dict['outliers_detect'] and clustering.all_outliers_x_cluster[K]:
                Out = zip(*[scores[i] for i in clustering.all_outliers_x_cluster[K]])
                pylab.plot(Out[0], Out[1], color = _color, marker = 'x', linestyle = 'None')
                
                # Only if using the first Principal Component
                #pylab.plot(Out[0], color = _color, marker = 'x', linestyle = 'None')
        else: print 'No data in cluster', K
    
    pylab.show()
    
    #print 'Testing performances of different kmeans implementations and algorithms...'
    #kmeans_dict = {'min_clusters_nr': 2, 'max_clusters_nr': 8, \
    #               'replicates_nr': 100, 'normalization': False, 'algorithm': None}
    #for alg in (kmeans, fuzzy_kmeans):
    #    kmeans_dict['algorithm'] = alg
    #    start = time.clock()
    #    clustering = kmeans_replica(M, **kmeans_dict)
    #    print '>>> Time required to', str(alg), 'algorithm:', time.clock() - start
"""
"""
    M = [[0, 0.1], [0.1, 0], [10.1, 10], [10.2, 10.4], [5, 5.05]]
    for nn in xrange(2, 7):
        C = kmeans(M, nn, scaling_down = False)
        for (cl, dl) in C.clusters_data.items(): print '- cluster', cl, ':', dl
    
    kmeans_dict = {'min_clusters_nr':       8,
                   'max_clusters_nr':       12,
                   'replicates_nr':         50,
                   'normalization':         False,
                   'algorithm':             kmeans,
                   'scatter_method':        scatterwithin_no_samplesnr,
                   'outliers_detect':       True,
                   'Mahalanobis_threshold': 1.5,
                   'scaling_down':          True}
"""

"""
    M = [[pdf(0, 1.5), pdf(0, 1)] for i in range(50)] + \
        [[pdf(4.5, 1), pdf(9, 1.5)] for i in range(50)] + \
        [[pdf(9, 1), pdf(0, 1.5)] for i in range(50)]
    #[[random.uniform(6, 8), random.uniform(-2, 2)] for i in range(100)]
    
    kmeans_dict = { 'min_clusters_nr':       3, \
                    'max_clusters_nr':       3, \
                    'replicates_nr':         25, \
                    'normalization':         False, \
                    'algorithm':             fuzzy_kmeans, \
                    'scatter_method':        scatterwithin_no_samplesnr, \
                    'outliers_detect':       True, \
                    'Mahalanobis_threshold': 1.5, \
                    'scaling_down':          True}
"""
