K-means algorithm for clustering

    Differences from scikit-learn:

    * init='k-means||'
    * oversampling_factor keyword
    * n_jobs=-1
14594
def k_means(
    X,
    n_clusters,
    init="k-means||",
    precompute_distances="auto",
    n_init=1,
    max_iter=300,
    verbose=False,
    tol=1e-4,
    random_state=None,
    copy_x=True,
    n_jobs=-1,
    algorithm="full",
    return_n_iter=False,
    oversampling_factor=2,
    init_max_iter=None,
):
    
    labels, inertia, centers, n_iter = _kmeans_single_lloyd(
        X,
        n_clusters,
        max_iter=max_iter,
        init=init,
        verbose=verbose,
        tol=tol,
        random_state=random_state,
        oversampling_factor=oversampling_factor,
        init_max_iter=init_max_iter,
    )
    if return_n_iter:
        return labels, centers, inertia, n_iter
    else:
        return labels, centers, inertia
513323
def cluster(args):
    

    args = _check_args(args)
    read_stats_file = op.join(args.dir_out, "read_stats.tsv")
    if file_exists(read_stats_file):
        os.remove(read_stats_file)

    bam_file, seq_obj = _clean_alignment(args)

    logger.info("Parsing matrix file")
    seqL, y, l = parse_ma_file(seq_obj, args.ffile)
    # y, l = _total_counts(seqL.keys(), seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'aligned'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    if len(seqL.keys()) < 10:
        logger.error("It seems you have low coverage. Please check your fastq files have enough sequences.")
        raise ValueError("So few sequences.")

    logger.info("Cleaning bam file")
    y, l = _total_counts(seqL.keys(), seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'cleaned'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    clusL = _create_clusters(seqL, bam_file, args)
    y, l = _total_counts(clusL.seq.keys(), clusL.seq, aligned=True)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'clusters'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    logger.info("Solving multi-mapping events in the network of clusters")
    clusLred = _cleaning(clusL, args.dir_out)
    y, l = _total_counts(clusLred.clus, seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'meta-cluster'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')
    logger.info("Clusters up to %s" % (len(clusLred.clus.keys())))

    if args.show:
        logger.info("Creating sequences alignment to precursor")
        clusLred = show_seq(clusLred, args.index)

    clusLred = peak_calling(clusLred)

    clusLred = _annotate(args, clusLred)

    logger.info("Creating json and count matrix")
    json_file = _create_json(clusLred, args)

    logger.info("Output file in: %s" % args.dir_out)

    if args.db:
        name = args.db + ".db"
        logger.info("Create database: database/" + name)
        data = load_data(json_file)
        out_dir = op.join(args.dir_out, "database")
        make_database(data, name, out_dir)
    logger.info("Finished")
268993
def _kmeans_run(X, n_clusters, max_iter, tol):
    
    membs = np.empty(shape=X.shape[0], dtype=int)
    centers = _kmeans_init(X, n_clusters)

    sse_last = 9999.9
    n_iter = 0
    for it in range(1,max_iter):
        membs = _assign_clusters(X, centers)
        centers,sse_arr = _update_centers(X, membs, n_clusters)
        sse_total = np.sum(sse_arr)
        if np.abs(sse_total - sse_last) < tol:
            n_iter = it
            break
        sse_last = sse_total

    return(centers, membs, sse_total, sse_arr, n_iter)
1125515
def fit(self, X):
        
        self.centers_, self.labels_, self.sse_arr_, self.n_iter_ = \
              _kmedoids(X, self.n_clusters, self.distance, self.max_iter, self.n_trials, self.tol, self.rng)
880466
def poisson_consensus_se(data, k, n_runs=10, **se_params):
    
    clusters = []
    for i in range(n_runs):
        assignments, means  = poisson_cluster(data, k)
        clusters.append(assignments)
    clusterings = np.vstack(clusters)
    consensus = CE.cluster_ensembles(clusterings, verbose=False, N_clusters_max=k)
    init_m, init_w = nmf_init(data, consensus, k, 'basic')
    M, W, ll = poisson_estimate_state(data, k, init_means=init_m, init_weights=init_w, **se_params)
    return M, W, ll
471086
def cluster(list_of_texts, num_clusters=3):
    
    pipeline = Pipeline([
        ("vect", CountVectorizer()),
        ("tfidf", TfidfTransformer()),
        ("clust", KMeans(n_clusters=num_clusters))
    ])

    try:
        clusters = pipeline.fit_predict(list_of_texts)
    except ValueError:
        clusters = list(range(len(list_of_texts)))

    return clusters
780374
function cluster(k) {
                var Centroids = this.createCentroids(k);
                var LastDistances = Matrix$d.Zero(this.Observations.rows(), this.Observations.cols());
                var Distances = this.distanceFrom(Centroids);
                var Groups;

                while(!(LastDistances.eql(Distances))) {
                Groups = Distances.minColumnIndexes();
                LastDistances = Distances;

                var newCentroids = [];

                for(var i = 1; i <= Centroids.rows(); i++) {
                    var centroid = [];

                    for(var j = 1; j <= Centroids.cols(); j++) {
                    var sum = 0;
                    var count = 0;

                    for(var l = 1; l <= this.Observations.rows(); l++) {
                        if(Groups.e(l) == i) {
                        count++;
                        sum += this.Observations.e(l, j);
                        }
                    }

                    centroid.push(sum / count);
                    }

                    newCentroids.push(centroid);
                }

                Centroids = $M(newCentroids);
                Distances = this.distanceFrom(Centroids);
                }

                return Groups;
            }
592952
def _regroup(self):
        
        # clean up old maps
        for j in range(self.nout):
            self.inv_map[j] = []

        # find smallest divergence between input component i
        # and output component j of the cluster mixture density
        for i in range(self.nin):
            self.min_kl[i] = np.inf
            j_min = None
            for j in range(self.nout):
                kl = kullback_leibler(self.f.components[i], self.g.components[j])
                if kl < self.min_kl[i]:
                    self.min_kl[i] = kl
                    j_min = j
            assert j_min is not None
            self.inv_map[j_min].append(i)
190371
def _k_prototypes_iter(Xnum, Xcat, centroids, cl_attr_sum, cl_memb_sum, cl_attr_freq,
                       membship, num_dissim, cat_dissim, gamma, random_state):
    """"""
    moves = 0
    for ipoint in range(Xnum.shape[0]):
        clust = np.argmin(
            num_dissim(centroids[0], Xnum[ipoint]) +
            gamma * cat_dissim(centroids[1], Xcat[ipoint], X=Xcat, membship=membship)
        )
        if membship[clust, ipoint]:
            # Point is already in its right place.
            continue

        # Move point, and update old/new cluster frequencies and centroids.
        moves += 1
        old_clust = np.argwhere(membship[:, ipoint])[0][0]

        # Note that membship gets updated by kmodes.move_point_cat.
        # move_point_num only updates things specific to the k-means part.
        cl_attr_sum, cl_memb_sum = move_point_num(
            Xnum[ipoint], clust, old_clust, cl_attr_sum, cl_memb_sum
        )
        cl_attr_freq, membship, centroids[1] = kmodes.move_point_cat(
            Xcat[ipoint], ipoint, clust, old_clust,
            cl_attr_freq, membship, centroids[1]
        )

        # Update old and new centroids for numerical attributes using
        # the means and sums of all values
        for iattr in range(len(Xnum[ipoint])):
            for curc in (clust, old_clust):
                if cl_memb_sum[curc]:
                    centroids[0][curc, iattr] = cl_attr_sum[curc, iattr] / cl_memb_sum[curc]
                else:
                    centroids[0][curc, iattr] = 0.

        # In case of an empty cluster, reinitialize with a random point
        # from largest cluster.
        if not cl_memb_sum[old_clust]:
            from_clust = membship.sum(axis=1).argmax()
            choices = [ii for ii, ch in enumerate(membship[from_clust, :]) if ch]
            rindx = random_state.choice(choices)

            cl_attr_sum, cl_memb_sum = move_point_num(
                Xnum[rindx], old_clust, from_clust, cl_attr_sum, cl_memb_sum
            )
            cl_attr_freq, membship, centroids[1] = kmodes.move_point_cat(
                Xcat[rindx], rindx, old_clust, from_clust,
                cl_attr_freq, membship, centroids[1]
            )

    return centroids, moves
1117715
def k_prototypes(X, categorical, n_clusters, max_iter, num_dissim, cat_dissim,
                 gamma, init, n_init, verbose, random_state, n_jobs):
    """"""
    random_state = check_random_state(random_state)
    if sparse.issparse(X):
        raise TypeError("k-prototypes does not support sparse data.")

    if categorical is None or not categorical:
        raise NotImplementedError(
            "No categorical data selected, effectively doing k-means. "
            "Present a list of categorical columns, or use scikit-learn's "
            "KMeans instead."
        )
    if isinstance(categorical, int):
        categorical = [categorical]
    assert len(categorical) != X.shape[1], \
        "All columns are categorical, use k-modes instead of k-prototypes."
    assert max(categorical) < X.shape[1], \
        "Categorical index larger than number of columns."

    ncatattrs = len(categorical)
    nnumattrs = X.shape[1] - ncatattrs
    n_points = X.shape[0]
    assert n_clusters <= n_points, "Cannot have more clusters ({}) " \
                                   "than data points ({}).".format(n_clusters, n_points)

    Xnum, Xcat = _split_num_cat(X, categorical)
    Xnum, Xcat = check_array(Xnum), check_array(Xcat, dtype=None)

    # Convert the categorical values in Xcat to integers for speed.
    # Based on the unique values in Xcat, we can make a mapping to achieve this.
    Xcat, enc_map = encode_features(Xcat)

    # Are there more n_clusters than unique rows? Then set the unique
    # rows as initial values and skip iteration.
    unique = get_unique_rows(X)
    n_unique = unique.shape[0]
    if n_unique <= n_clusters:
        max_iter = 0
        n_init = 1
        n_clusters = n_unique
        init = list(_split_num_cat(unique, categorical))
        init[1], _ = encode_features(init[1], enc_map)

    # Estimate a good value for gamma, which determines the weighing of
    # categorical values in clusters (see Huang [1997]).
    if gamma is None:
        gamma = 0.5 * Xnum.std()

    results = []
    seeds = random_state.randint(np.iinfo(np.int32).max, size=n_init)
    if n_jobs == 1:
        for init_no in range(n_init):
            results.append(k_prototypes_single(Xnum, Xcat, nnumattrs, ncatattrs,
                                               n_clusters, n_points, max_iter,
                                               num_dissim, cat_dissim, gamma,
                                               init, init_no, verbose, seeds[init_no]))
    else:
        results = Parallel(n_jobs=n_jobs, verbose=0)(
            delayed(k_prototypes_single)(Xnum, Xcat, nnumattrs, ncatattrs,
                                         n_clusters, n_points, max_iter,
                                         num_dissim, cat_dissim, gamma,
                                         init, init_no, verbose, seed)
            for init_no, seed in enumerate(seeds))
    all_centroids, all_labels, all_costs, all_n_iters = zip(*results)

    best = np.argmin(all_costs)
    if n_init > 1 and verbose:
        print("Best run was number {}".format(best + 1))

    # Note: return gamma in case it was automatically determined.
    return all_centroids[best], enc_map, all_labels[best], \
        all_costs[best], all_n_iters[best], gamma
1019134
def fit(self, X, y=None, **kwargs):
        
        X = pandas_to_numpy(X)

        random_state = check_random_state(self.random_state)
        self._enc_cluster_centroids, self._enc_map, self.labels_,\
            self.cost_, self.n_iter_ = k_modes(X,
                                               self.n_clusters,
                                               self.max_iter,
                                               self.cat_dissim,
                                               self.init,
                                               self.n_init,
                                               self.verbose,
                                               random_state,
                                               self.n_jobs)
        return self
1241492
def eeg_microstates_clustering(data, n_microstates=4, clustering_method="kmeans", n_jobs=1, n_init=25, occurence_rejection_treshold=0.05, max_refitting=5, verbose=True):
    
    # Create training set
    training_set = data.copy()

    if verbose is True:
        print("- Initializing the clustering algorithm...")
    if clustering_method == "kmeans":
        algorithm = sklearn.cluster.KMeans(init='k-means++', n_clusters=n_microstates, n_init=n_init, n_jobs=n_jobs)
    elif clustering_method == "spectral":
        algorithm = sklearn.cluster.SpectralClustering(n_clusters=n_microstates, n_init=n_init, n_jobs=n_jobs)
    elif clustering_method == "agglom":
        algorithm = sklearn.cluster.AgglomerativeClustering(n_clusters=n_microstates, linkage="complete")
    elif clustering_method == "dbscan":
        algorithm = sklearn.cluster.DBSCAN(min_samples=100)
    elif clustering_method == "affinity":
        algorithm = sklearn.cluster.AffinityPropagation(damping=0.5)
    else:
        print("NeuroKit Error: eeg_microstates(): clustering_method must be 'kmeans', 'spectral', 'dbscan', 'affinity' or 'agglom'")


    refitting = 0  # Initialize the number of refittings
    good_fit_achieved = False
    while good_fit_achieved is False:
        good_fit_achieved = True
        if verbose is True:
            print("- Fitting the classifier...")
        # Fit the algorithm
        algorithm.fit(training_set)

        if verbose is True:
            print("- Clustering back the initial data...")
        # Predict the more likely cluster for each observation
        predicted = algorithm.fit_predict(training_set)

        if verbose is True:
            print("- Check for abnormalities...")
        # Check for abnormalities and prune the training set until none found
        occurences = dict(collections.Counter(predicted))
        masks = [np.array([True]*len(training_set))]
        for microstate in occurences:
            # is the frequency of one microstate inferior to a treshold
            if occurences[microstate] < len(data)*occurence_rejection_treshold:
                good_fit_achieved = False
                refitting += 1  # Increment the refitting
                print("NeuroKit Warning: eeg_microstates(): detected some outliers: refitting the classifier (n=" + str(refitting) + ").")
                masks.append(predicted!=microstate)
        mask = np.all(masks, axis=0)
        training_set = training_set[mask]

    return(algorithm)
1332418
def calc_cat_clust_order(net, inst_rc):
  '''
  
  '''
  from .__init__ import Network
  from copy import deepcopy
  from . import calc_clust, run_filter

  inst_keys = list(net.dat['node_info'][inst_rc].keys())
  all_cats = [x for x in inst_keys if 'cat-' in x]

  if len(all_cats) > 0:

    for inst_name_cat in all_cats:

      tmp_name = 'dict_' + inst_name_cat.replace('-', '_')
      dict_cat = net.dat['node_info'][inst_rc][tmp_name]

      unordered_cats = dict_cat.keys()

      ordered_cats = order_categories(unordered_cats)

      # this is the ordering of the columns based on their category, not
      # including their clustering ordering within category
      all_cat_orders = []
      tmp_names_list = []
      for inst_cat in ordered_cats:

        inst_nodes = dict_cat[inst_cat]

        tmp_names_list.extend(inst_nodes)

      #   cat_net = deepcopy(Network())

      #   cat_net.dat['mat'] = deepcopy(net.dat['mat'])
      #   cat_net.dat['nodes'] = deepcopy(net.dat['nodes'])

      #   cat_df = cat_net.dat_to_df()

      #   sub_df = {}
      #   if inst_rc == 'col':
      #     sub_df['mat'] = cat_df['mat'][inst_nodes]
      #   elif inst_rc == 'row':
      #     # need to transpose df
      #     cat_df['mat'] = cat_df['mat'].transpose()
      #     sub_df['mat'] = cat_df['mat'][inst_nodes]
      #     sub_df['mat'] = sub_df['mat'].transpose()

      #   # filter matrix before clustering
      #   ###################################
      #   threshold = 0.0001
      #   sub_df = run_filter.df_filter_row_sum(sub_df, threshold)
      #   sub_df = run_filter.df_filter_col_sum(sub_df, threshold)

      #   # load back to dat
      #   cat_net.df_to_dat(sub_df)

      #   cat_mat_shape = cat_net.dat['mat'].shape

      #   print('***************')
      #   try:
      #     if cat_mat_shape[0]>1 and cat_mat_shape[1] > 1 and all_are_numbers == False:

      #       calc_clust.cluster_row_and_col(cat_net, 'cos')
      #       inst_cat_order = cat_net.dat['node_info'][inst_rc]['clust']
      #     else:
      #       inst_cat_order = list(range(len(cat_net.dat['nodes'][inst_rc])))

      #   except:
      #     inst_cat_order = list(range(len(cat_net.dat['nodes'][inst_rc])))


      #   prev_order_len = len(all_cat_orders)

      #   # add prev order length to the current order number
      #   inst_cat_order = [i + prev_order_len for i in inst_cat_order]
      #   all_cat_orders.extend(inst_cat_order)

      # # generate ordered list of row/col names, which will be used to
      # # assign the order to specific nodes
      # names_clust_list = [x for (y, x) in sorted(zip(all_cat_orders,
      #                     tmp_names_list))]

      names_clust_list = tmp_names_list

      # calc category-cluster order
      final_order = []

      for i in range(len(net.dat['nodes'][inst_rc])):

        inst_node_name = net.dat['nodes'][inst_rc][i]
        inst_node_num = names_clust_list.index(inst_node_name)

        final_order.append(inst_node_num)

      inst_index_cat = inst_name_cat.replace('-', '_') + '_index'

      net.dat['node_info'][inst_rc][inst_index_cat] = final_order
127092
def k_modes(X, n_clusters, max_iter, dissim, init, n_init, verbose, random_state, n_jobs):
    """"""
    random_state = check_random_state(random_state)
    if sparse.issparse(X):
        raise TypeError("k-modes does not support sparse data.")

    X = check_array(X, dtype=None)

    # Convert the categorical values in X to integers for speed.
    # Based on the unique values in X, we can make a mapping to achieve this.
    X, enc_map = encode_features(X)

    n_points, n_attrs = X.shape
    assert n_clusters <= n_points, "Cannot have more clusters ({}) " \
                                   "than data points ({}).".format(n_clusters, n_points)

    # Are there more n_clusters than unique rows? Then set the unique
    # rows as initial values and skip iteration.
    unique = get_unique_rows(X)
    n_unique = unique.shape[0]
    if n_unique <= n_clusters:
        max_iter = 0
        n_init = 1
        n_clusters = n_unique
        init = unique

    results = []
    seeds = random_state.randint(np.iinfo(np.int32).max, size=n_init)
    if n_jobs == 1:
        for init_no in range(n_init):
            results.append(k_modes_single(X, n_clusters, n_points, n_attrs, max_iter,
                                          dissim, init, init_no, verbose, seeds[init_no]))
    else:
        results = Parallel(n_jobs=n_jobs, verbose=0)(
            delayed(k_modes_single)(X, n_clusters, n_points, n_attrs, max_iter,
                                    dissim, init, init_no, verbose, seed)
            for init_no, seed in enumerate(seeds))
    all_centroids, all_labels, all_costs, all_n_iters = zip(*results)

    best = np.argmin(all_costs)
    if n_init > 1 and verbose:
        print("Best run was number {}".format(best + 1))

    return all_centroids[best], enc_map, all_labels[best], \
        all_costs[best], all_n_iters[best]
384844
def fit(self, X):
        
        self.centers_, self.labels_, self.sse_arr_, self.n_iter_ = \
              _kmeans(X, self.n_clusters, self.max_iter, self.n_trials, self.tol)
641275
def _kmedoids_run(X, n_clusters, distance, max_iter, tol, rng):
    
    membs = np.empty(shape=X.shape[0], dtype=int)
    centers = kmeans._kmeans_init(X, n_clusters, method='', rng=rng)

    sse_last = 9999.9
    n_iter = 0
    for it in range(1,max_iter):
        membs = kmeans._assign_clusters(X, centers)
        centers,sse_arr = _update_centers(X, membs, n_clusters, distance)
        sse_total = np.sum(sse_arr)
        if np.abs(sse_total - sse_last) < tol:
            n_iter = it
            break
        sse_last = sse_total

    return(centers, membs, sse_total, sse_arr, n_iter)
982125
def _make_clusters(self, matrix, num_clusters_per_roi, metric):
        """"""

        from scipy.cluster.hierarchy import fclusterdata

        # maxclust needed to ensure t is interpreted as # clusters in heirarchical clustering
        group_ids = fclusterdata(matrix, metric=metric, t=num_clusters_per_roi,
                                 criterion='maxclust')
        group_set = np.unique(group_ids)
        clusters = [
            self._summary_func(matrix[group_ids == group, :], axis=0, keepdims=True)
            for group in group_set]

        return np.vstack(clusters).squeeze()