K-Means initialization using k-means||

    This is algorithm 2 in Scalable K-Means++ (2012).
5027
def init_scalable(
    X, n_clusters, random_state=None, max_iter=None, oversampling_factor=2
):
    

    logger.info("Initializing with k-means||")
    # Step 1: Initialize Centers
    idx = 0
    centers = da.compute(X[idx, np.newaxis])[0]
    c_idx = {idx}

    # Step 2: Initialize cost
    cost, = compute(evaluate_cost(X, centers))

    if cost == 0:
        n_iter = 0
    else:
        n_iter = int(np.round(np.log(cost)))

    if max_iter is not None:
        n_iter = min(max_iter, n_iter)

    # Steps 3 - 6: update candidate Centers
    for i in range(n_iter):
        with _timer(
            "init iteration %2d/%2d , %2d centers" % (i + 1, n_iter, len(c_idx)),
            _logger=logger,
        ):
            new_idxs = _sample_points(X, centers, oversampling_factor, random_state)
            new_idxs = set(*compute(new_idxs))
            c_idx |= new_idxs

        # Sort before slicing, for better performance / memory
        # usage with the scheduler.
        # See https://github.com/dask/dask-ml/issues/39
        centers = X[sorted(c_idx)].compute()

    # XXX: scikit-learn doesn't have weighted k-means.
    # The paper weights each center by the number of points closest to it.
    # https://stackoverflow.com/a/37198799/1889400 claims you can scale the
    # features before clustering, but that doesn't seem right.
    # I think that replicating the *points*, proportional to the number of
    # original points closest to the candidate centers, would be a better way
    # to do that.

    if len(centers) < n_clusters:
        logger.warning("Found fewer than %d clusters in init.", n_clusters)
        # supplement with random
        need = n_clusters - len(centers)
        locs = sorted(
            random_state.choice(
                np.arange(0, len(X)), size=need, replace=False, chunks=len(X)
            )
        )
        extra = X[locs].compute()
        return np.vstack([centers, extra])
    else:
        # Step 7, 8 without weights
        # dask RandomState objects aren't valid for scikit-learn
        rng2 = (
            random_state.randint(0, np.iinfo("i4").max - 1, chunks=())
            .compute(scheduler="single-threaded")
            .item()
        )
        km = sk_k_means.KMeans(n_clusters, random_state=rng2)
        km.fit(centers)

    return km.cluster_centers_
350379
def K(self, X, X2=None):        
    # model :   -a d^2y/dx^2  + b dy/dt + c * y = U
    # kernel Kyy rbf spatiol temporal
    # vyt Y temporal variance  vyx Y spatiol variance   lyt Y temporal lengthscale   lyx Y spatiol lengthscale
    # kernel Kuu doper( doper(Kyy))
    # a   b    c    lyt   lyx    vyx*vyt
        """"""        
        X,slices = X[:,:-1],index_to_slices(X[:,-1])
        if X2 is None:
            X2,slices2 = X,slices
            K = np.zeros((X.shape[0], X.shape[0]))
        else:
            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
            K = np.zeros((X.shape[0], X2.shape[0]))


        tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
        xdist = (X[:,1][:,None] - X2[:,1][None,:])**2

        ttdist = (X[:,0][:,None] - X2[:,0][None,:])
        #rdist = [tdist,xdist]
        #dist = np.abs(X - X2.T)
        vyt = self.variance_Yt
        vyx = self.variance_Yx
        
        lyt=1/(2*self.lengthscale_Yt)
        lyx=1/(2*self.lengthscale_Yx)

        a = self.a ## -a is used in the model, negtive diffusion
        b = self.b
        c = self.c

        kyy = lambda tdist,xdist: np.exp(-lyt*(tdist) -lyx*(xdist))

        k1 = lambda tdist: (2*lyt - 4*lyt**2 * (tdist) )

        k2 = lambda xdist: ( 4*lyx**2 * (xdist)  - 2*lyx )

        k3 = lambda xdist: ( 3*4*lyx**2 - 6*8*xdist*lyx**3 + 16*xdist**2*lyx**4 )

        k4 = lambda ttdist: 2*lyt*(ttdist)

        for i, s1 in enumerate(slices):
            for j, s2 in enumerate(slices2):
                for ss1 in s1:
                    for ss2 in s2:
                        if i==0 and j==0:
                            K[ss1,ss2] = vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                        elif i==0 and j==1:
                            K[ss1,ss2] = (-a*k2(xdist[ss1,ss2]) + b*k4(ttdist[ss1,ss2]) + c)*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            #K[ss1,ss2]=  np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[ss1,ss2]) )   )
                            #K[ss1,ss2]=  np.where(  rdist[ss1,ss2]>0 , kuyp(rdist[ss1,ss2]), kuyn(rdist[ss1,ss2] )   )
                        elif i==1 and j==1:
                            K[ss1,ss2] = ( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 )* vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                        else:
                            K[ss1,ss2] = (-a*k2(xdist[ss1,ss2]) - b*k4(ttdist[ss1,ss2]) + c)*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            #K[ss1,ss2]= np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[ss1,ss2]) )   )
                            #K[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(rdist[ss1,ss2]), kyun(rdist[ss1,ss2] )   )
        
        #stop
        return K
675658
public void initialState(int K) {
        int i;

        int M = documents.length;

        // initialise count variables.
        word_topic_matrix = new int[V][K];
        nd = new int[M][K];
        nwsum = new int[K];
        ndsum = new int[M];

        // The z_i are are initialised to values in [1,K] to determine the
        // initial state of the Markov chain.

        z = new int[M][];
        for (int m = 0; m < M; m++) {
            int N = documents[m].length;
            z[m] = new int[N];
            for (int n = 0; n < N; n++) {
                int topic = (int) (Math.random() * K);
                z[m][n] = topic;
                // number of instances of word i assigned to topic j
                word_topic_matrix[documents[m][n]][topic]++;
                // number of words in document i assigned to topic j.
                nd[m][topic]++;
                // total number of words assigned to topic j.
                nwsum[topic]++;
            }
            // total number of words in document i
            ndsum[m] = N;
        }
    }
1020947
@Override
	public void trainOnInstanceImpl(Instance inst) {
		double[] x = inst.toDoubleArray();
		if (this.numDimensions != x.length) {
			System.out.println("Line skipped because line dimension is "
					+ x.length + " instead of " + this.numDimensions);
			return;
		}
		// Starts with the buffer phase to calculate the starting threshold
		if (this.bufferPhase) {
			// Calculates the pairwise distance to all unequal point in the
			// buffer
			for (double[] point : this.buffer) {
				double d = Metric.distanceSquared(point, x);
				if (d > 0) {
					this.pairwiseDifferent++;
					if (d < minDistance) {
						minDistance = d;
					}
				}
			}
			this.buffer.add(x);

			// Checks if the buffer is large enough
			if (this.pairwiseDifferent >= this.maxNumClusterFeatures + 1) {
				// Calculates the starting threshold
				this.T = 16 * minDistance;
				this.root.setThreshold(calcRSquared(1));
				this.bufferPhase = false;
				// Adds all points to the ClusteringFeature tree
				for (double[] point : this.buffer) {
					bicoUpdate(point);
				}
				this.buffer.clear();
				this.buffer = null;
			}
		} else {
			// Adds the point directly to the ClusteringFeature tree
			bicoUpdate(x);
		}
	}
673899
private static void clusterIteration(Matrix matrix,
                                         int numClusters,
                                         KMeansSeed seedType,
                                         CriterionFunction criterion) {
        DoubleVector[] centers = seedType.chooseSeeds(numClusters, matrix);

        // Compute the initial set of assignments for each data point based on
        // the initial assignments.
        int[] initialAssignments = new int[matrix.rows()];

        // If there is to be only one cluster, then everything will be auto
        // assigned to the first cluster.  This is just a special case that only
        // comes up when comparing other solutions to the non-solution.
        if (numClusters != 1) {
            int nc = 0;
            for (int i = 0; i < matrix.rows(); ++i) {

                DoubleVector vector = matrix.getRowVector(i);
                double bestSimilarity = 0;
                for (int c = 0; c < numClusters; ++c) {
                    double similarity = Similarity.cosineSimilarity(
                            centers[c], vector);
                    nc++;
                    if (similarity >= bestSimilarity) {
                        bestSimilarity = similarity;
                        initialAssignments[i] = c;
                    }
                }
            }
        }

        // Setup the criterion function with it's meta data.
        criterion.setup(matrix, initialAssignments, numClusters);

        // Iteratively swap each data point to a better cluster if such an
        // assignment exists.
        List<Integer> indices = new ArrayList<Integer>(matrix.rows());
        for (int i = 0; i < matrix.rows(); ++i)
            indices.add(i);

        // Iterate through each data point in random order.  For each data
        // point, try to assign it to a new cluster.  If no data point is moved
        // in an iteration, end the iterations.
        boolean changed = true;
        while (changed) {
            changed = false;
            Collections.shuffle(indices);
            for (int index : indices)
                changed |= criterion.update(index);
        }
    }
112257
def fit_cluster(data, assignments, k, P_init, R_init, means):
    
    for c in range(k):
        if data[:,assignments==c].shape[1] == 0:
            _, assignments = kmeans_pp(data, k)
    genes, cells = data.shape
    nb_gene_indices = np.array([True for i in range(genes)])
    for c in range(k):
        c_data = data[:,assignments==c]
        nb_gene_indices = nb_gene_indices & find_nb_genes(c_data)
    for c in range(k):
        c_data = data[:,assignments==c]
        nb_genes = c_data[nb_gene_indices,:]
        poisson_genes = c_data[~nb_gene_indices, :]
        P_init[nb_gene_indices, c], R_init[nb_gene_indices, c] = nb_fit(nb_genes)
        means[~nb_gene_indices, c] = poisson_genes.mean(1)
    return nb_gene_indices
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]
1277210
private void initialClustering() {
		//System.out.println("INITIAL CLUSTERING CALLED");
		//printDStreamState();
		// 1. Update the density of all grids in grid_list

		updateGridListDensity();
		//printGridList();
				
		// 2. Assign each dense grid to a distinct cluster
		// and
		// 3. Label all other grids as NO_CLASS	
		Iterator<Map.Entry<DensityGrid, CharacteristicVector>> glIter = this.grid_list.entrySet().iterator();
		HashMap<DensityGrid, CharacteristicVector> newGL = new HashMap<DensityGrid, CharacteristicVector>();
		
		while(glIter.hasNext())
		{
			Map.Entry<DensityGrid, CharacteristicVector> grid = glIter.next();
			DensityGrid dg = grid.getKey();
			CharacteristicVector cvOfG = grid.getValue();

			//System.out.print(dg.toString());
			if(cvOfG.getAttribute() == DENSE)
			{
				int gridClass = this.cluster_list.size();
				cvOfG.setLabel(gridClass);
				GridCluster gc = new GridCluster ((CFCluster)dg, new ArrayList<CFCluster>(), gridClass);
				gc.addGrid(dg);
				this.cluster_list.add(gc);
				//System.out.print(" was dense (class "+gridClass+")");
			}
			else
				cvOfG.setLabel(NO_CLASS);

			//System.out.println();
			newGL.put(dg, cvOfG);
		}

		this.grid_list = newGL;
		//printGridClusters();
		
		// 4. Make changes to grid labels by doing:
		//    a. For each cluster c
		//    b. For each outside grid g of c
		//    c. For each neighbouring grid h of g
		//    d. If h belongs to c', label c and c' with 
		//       the label of the largest cluster
		//    e. Else if h is transitional, assign it to c
		//    f. While changes can be made

		boolean changesMade;

		do{	
			changesMade = adjustLabels();
		}while(changesMade);	// while changes are being made
		
		//printGridList();
		//printGridClusters();
	}
980034
@Reference(authors = "D. Pelleg, A. Moore", //
      title = "X-means: Extending K-means with Efficient Estimation on the Number of Clusters", //
      booktitle = "Proc. 17th Int. Conf. on Machine Learning (ICML 2000)", //
      url = "http://www.pelleg.org/shared/hp/download/xmeans.ps", //
      bibkey = "DBLP:conf/icml/PellegM00")
  public static <V extends NumberVector> double logLikelihood(Relation<V> relation, Clustering<? extends MeanModel> clustering, NumberVectorDistanceFunction<? super V> distanceFunction) {
    List<? extends Cluster<? extends MeanModel>> clusters = clustering.getAllClusters();
    // number of clusters
    final int m = clusters.size();

    // number of objects in the clustering
    int n = 0;
    // cluster sizes
    int[] n_i = new int[m];
    // total variance
    double d = 0.;
    // variances
    double[] d_i = new double[m];

    // Iterate over clusters:
    Iterator<? extends Cluster<? extends MeanModel>> it = clusters.iterator();
    for(int i = 0; it.hasNext(); ++i) {
      Cluster<? extends MeanModel> cluster = it.next();
      n += n_i[i] = cluster.size();
      d += d_i[i] = varianceOfCluster(cluster, distanceFunction, relation);
    }

    // No remaining variance, if every point is on its own:
    if(n <= m) {
      return Double.NEGATIVE_INFINITY;
    }

    // Total variance (corrected for bias)
    final double logv = FastMath.log(d / (n - m));

    final int dim = RelationUtil.dimensionality(relation);
    // log likelihood of this clustering
    double logLikelihood = 0.;

    // Aggregate
    for(int i = 0; i < m; i++) {
      logLikelihood += n_i[i] * FastMath.log(n_i[i]) // Post. entropy Rn log Rn
          - n_i[i] * .5 * MathUtil.LOGTWOPI // Rn/2 log2pi
          - n_i[i] * dim * .5 * logv // Rn M/2 log sigma^2
          - (d_i[i] - m) * .5; // (Rn-K)/2
    }
    logLikelihood -= n * FastMath.log(n); // Prior entropy, sum_i Rn log R
    return logLikelihood;
  }
473196
public IModel KMeans(int numCentroids, int maxIters, int runs, String initMode)
      throws DDFException {
    return this.train("kmeans", numCentroids, maxIters, runs, initMode);
  }
448157
public static KMeansClustering setup(int clusterCount, int maxIterationCount, Distance distanceFunction,
                    boolean inverse) {
        ClusteringStrategy clusteringStrategy =
                        FixedClusterCountStrategy.setup(clusterCount, distanceFunction, inverse);
        clusteringStrategy.endWhenIterationCountEquals(maxIterationCount);
        return new KMeansClustering(clusteringStrategy);
    }
1081307
private Clustering<OPTICSModel> extractClusters(ClusterOrder clusterOrderResult, Relation<?> relation, double ixi, int minpts) {
    ArrayDBIDs clusterOrder = clusterOrderResult.ids;
    DoubleDataStore reach = clusterOrderResult.reachability;

    DBIDArrayIter tmp = clusterOrder.iter();
    DBIDVar tmp2 = DBIDUtil.newVar();
    double mib = 0.0;
    List<SteepArea> salist = keepsteep ? new ArrayList<SteepArea>() : null;
    List<SteepDownArea> sdaset = new ArrayList<>();
    final Clustering<OPTICSModel> clustering = new Clustering<>("OPTICS Xi-Clusters", "optics");
    HashSet<Cluster<OPTICSModel>> curclusters = new HashSet<>();
    HashSetModifiableDBIDs unclaimedids = DBIDUtil.newHashSet(relation.getDBIDs());

    FiniteProgress scanprog = LOG.isVerbose() ? new FiniteProgress("OPTICS Xi cluster extraction", clusterOrder.size(), LOG) : null;
    for(SteepScanPosition scan = new SteepScanPosition(clusterOrderResult); scan.hasNext();) {
      if(scanprog != null) {
        scanprog.setProcessed(scan.index, LOG);
      }
      // Update maximum-inbetween
      mib = MathUtil.max(mib, scan.getReachability());
      // The last point cannot be the start of a steep area.
      if(!scan.next.valid()) {
        break;
      }
      // Xi-steep down area
      if(scan.steepDown(ixi)) {
        // Update mib values with current mib and filter
        updateFilterSDASet(mib, sdaset, ixi);
        final double startval = scan.getReachability();
        mib = 0.;
        int startsteep = scan.index, endsteep = scan.index;
        for(scan.next(); scan.hasNext(); scan.next()) {
          // still steep - continue.
          if(scan.steepDown(ixi)) {
            endsteep = scan.index;
            continue;
          }
          // Not going downward at all - stop here.
          // Always stop looking after minpts "flat" steps.
          if(!scan.steepDown(1.0) || scan.index - endsteep > minpts) {
            break;
          }
        }
        final SteepDownArea sda = new SteepDownArea(startsteep, endsteep, startval, 0);
        if(LOG.isDebuggingFinest()) {
          LOG.debugFinest("New steep down area: " + sda.toString());
        }
        sdaset.add(sda);
        if(salist != null) {
          salist.add(sda);
        }
        continue;
      }
      // Xi-steep up area
      if(scan.steepUp(ixi)) {
        // Update mib values with current mib and filter
        updateFilterSDASet(mib, sdaset, ixi);
        final SteepUpArea sua;
        // Compute steep-up area
        {
          int startsteep = scan.index, endsteep = scan.index;
          mib = scan.getReachability();
          double esuccr = scan.getNextReachability();
          // Find end of steep-up-area, eventually updating mib again
          while(!Double.isInfinite(esuccr) && scan.hasNext()) {
            scan.next();
            // still steep - continue.
            if(scan.steepUp(ixi)) {
              endsteep = scan.index;
              mib = scan.getReachability();
              esuccr = scan.getNextReachability();
              continue;
            }
            // Not going upward - stop here.
            // Stop looking after minpts non-up steps.
            if(!scan.steepUp(1.0) || scan.index - endsteep > minpts) {
              break;
            }
          }
          if(Double.isInfinite(esuccr)) {
            scan.next();
          }
          sua = new SteepUpArea(startsteep, endsteep, esuccr);
          if(LOG.isDebuggingFinest()) {
            LOG.debugFinest("New steep up area: " + sua.toString());
          }
          if(salist != null) {
            salist.add(sua);
          }
        }
        // Validate and computer clusters
        // LOG.debug("SDA size:"+sdaset.size()+" "+sdaset);
        ListIterator<SteepDownArea> sdaiter = sdaset.listIterator(sdaset.size());
        // Iterate backwards for correct hierarchy generation.
        while(sdaiter.hasPrevious()) {
          SteepDownArea sda = sdaiter.previous();
          if(LOG.isDebuggingFinest()) {
            LOG.debugFinest("Comparing: eU=" + mib + " SDA: " + sda.toString());
          }
          // Condition 3b: end-of-steep-up > maximum-in-between lower
          if(mib * ixi < sda.getMib()) {
            if(LOG.isDebuggingFinest()) {
              LOG.debugFinest("mib * ixi = " + mib * ixi + " >= sda.getMib() = " + sda.getMib());
            }
            continue;
          }
          // By default, clusters cover both the steep up and steep down area
          int cstart = sda.getStartIndex(),
              cend = MathUtil.min(sua.getEndIndex(), clusterOrder.size() - 1);
          // However, we sometimes have to adjust this (Condition 4):
          {
            // Case b)
            if(sda.getMaximum() * ixi >= sua.getMaximum()) {
              while(cstart < cend && //
                  reach.doubleValue(tmp.seek(cstart + 1)) > sua.getMaximum()) {
                cstart++;
              }
            }
            // Case c)
            else if(sua.getMaximum() * ixi >= sda.getMaximum()) {
              while(cend > cstart && //
                  reach.doubleValue(tmp.seek(cend - 1)) > sda.getMaximum()) {
                cend--;
              }
            }
            // Case a) is the default
          }
          // MST-based filtering technique of Schubert:
          // ensure that the predecessor is in the current cluster.
          // This filter removes common artifacts from the Xi method
          if(!nocorrect) {
            double startval = clusterOrderResult.reachability.doubleValue(tmp.seek(cstart));
            simplify: while(cend > cstart) {
              tmp.seek(cend);
              // If the reachability is less than the first points', then there
              // must be some point inbetween that "reached" the current point.
              if(clusterOrderResult.reachability.doubleValue(tmp) < startval) {
                break simplify;
              }
              clusterOrderResult.predecessor.assignVar(tmp, tmp2);
              // "Slow" search for predecessor:
              for(int i = cstart; i < cend; i++) {
                if(DBIDUtil.equal(tmp2, tmp.seek(i))) {
                  break simplify;
                }
              }
              // Not found, prune the last point
              if(LOG.isDebuggingFinest()) {
                LOG.debugFinest("Pruned one point by predecessor rule.");
              }
              --cend;
            }
          }
          // Condition 3a: obey minpts
          if(cend - cstart + 1 < minpts) {
            if(LOG.isDebuggingFinest()) {
              LOG.debugFinest("MinPts not satisfied.");
            }
            continue;
          }
          // Build the cluster
          ModifiableDBIDs dbids = DBIDUtil.newArray();
          for(int idx = cstart; idx <= cend; idx++) {
            tmp.seek(idx);
            // Collect only unclaimed IDs.
            if(unclaimedids.remove(tmp)) {
              dbids.add(tmp);
            }
          }
          if(LOG.isDebuggingFine()) {
            LOG.debugFine("Found cluster with " + dbids.size() + " new objects, length " + (cend - cstart + 1));
          }
          OPTICSModel model = new OPTICSModel(cstart, cend);
          Cluster<OPTICSModel> cluster = new Cluster<>("Cluster_" + cstart + "_" + cend, dbids, model);
          // Build the hierarchy
          for(Iterator<Cluster<OPTICSModel>> iter = curclusters.iterator(); iter.hasNext();) {
            Cluster<OPTICSModel> clus = iter.next();
            OPTICSModel omodel = clus.getModel();
            if(model.getStartIndex() <= omodel.getStartIndex() && omodel.getEndIndex() <= model.getEndIndex()) {
              clustering.addChildCluster(cluster, clus);
              iter.remove(); // Not a top-level cluster
            }
          }
          curclusters.add(cluster);
        }
        continue;
      }
      // Flat - advance anyway.
      scan.next();
    }
    if(scanprog != null) {
      scanprog.setProcessed(clusterOrder.size(), LOG);
    }
    if(!unclaimedids.isEmpty()) {
      boolean noise = reach.doubleValue(tmp.seek(clusterOrder.size() - 1)) >= Double.POSITIVE_INFINITY;
      Cluster<OPTICSModel> allcluster = new Cluster<>(noise ? "Noise" : "Cluster", unclaimedids, noise, new OPTICSModel(0, clusterOrder.size() - 1));
      for(Cluster<OPTICSModel> cluster : curclusters) {
        clustering.addChildCluster(allcluster, cluster);
      }
      clustering.addToplevelCluster(allcluster);
    }
    else {
      for(Cluster<OPTICSModel> cluster : curclusters) {
        clustering.addToplevelCluster(cluster);
      }
    }
    clustering.addChildResult(clusterOrderResult);
    if(salist != null) {
      clusterOrderResult.addChildResult(new SteepAreaResult(salist));
    }
    return clustering;
  }
1105152
public void merge(ClusteringFeature x) {
		assert (this.sumPoints.length == x.sumPoints.length);
		this.numPoints += x.numPoints;
		super.setWeight(this.numPoints);
		for (int i = 0; i < this.sumPoints.length; i++) {
			this.sumPoints[i] += x.sumPoints[i];
		}
		this.sumSquaredLength += x.sumSquaredLength;
	}
90508
public static void main(String[] args) throws Exception {
		
		if(!parseParameters(args)) {
			return;
		}
	
		// set up execution environment
		ExecutionEnvironment env = ExecutionEnvironment.getExecutionEnvironment();
		
		// get input data
		DataSet<Point> points = getPointDataSet(env);
		DataSet<Centroid> centroids = getCentroidDataSet(env);
		
		// set number of bulk iterations for KMeans algorithm
		IterativeDataSet<Centroid> loop = centroids.iterate(numIterations);
		
		DataSet<Centroid> newCentroids = points
			// compute closest centroid for each point
			.map(new SelectNearestCenter()).withBroadcastSet(loop, "centroids")
			// count and sum point coordinates for each centroid
			.map(new CountAppender())
			.groupBy(0).reduce(new CentroidAccumulator())
			// compute new centroids from point counts and coordinate sums
			.map(new CentroidAverager());
		
		// feed new centroids back into next iteration
		DataSet<Centroid> finalCentroids = loop.closeWith(newCentroids);
		
		DataSet<Tuple2<Integer, Point>> clusteredPoints = points
				// assign points to final clusters
				.map(new SelectNearestCenter()).withBroadcastSet(finalCentroids, "centroids");
		
		// emit result
		if(fileOutput) {
			clusteredPoints.writeAsCsv(outputPath, "\n", " ");
		} else {
			clusteredPoints.print();
		}

		// execute program
		env.execute("KMeans Example");
		
	}
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
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)
997394
def hierarchical_clusterization_2d(no_clusters = 0, distance_method = 0,
                                       vec = @values, debug = false)
      clusters = []

      if vec.length == 1
        hash = { vec[0] => 1 }
        cluster = PairCluster.new(hash)
        clusters.push(cluster)
        clusters
      end

      # Thresholds
      # threshold_distance = (0.25 * (vec.max-vec.min))
      threshold_density = (0.5 * vec.length).to_i

      # make a histogram from the input vector
      histogram = Hash[vec.group_by { |a| a }.map { |k, vs| [k, vs.length] }]

      # clusters = array of clusters
      # initially each length belongs to a different cluster
      histogram.each do |e|
        warn "pair (#{e[0].x} #{e[0].y}) appears #{e[1]} times" if debug
        hash = { e[0] => e[1] }
        cluster = PairCluster.new(hash)
        clusters.push(cluster)
      end

      clusters.each(&:print) if debug

      return clusters if clusters.length == 1

      # each iteration merge the closest two adiacent cluster
      # the loop stops according to the stop conditions
      iteration = 0
      loop do
        # stop condition 1
        break if no_clusters != 0 && clusters.length == no_clusters

        iteration += iteration
        warn "\nIteration #{iteration}" if debug

        min_distance = 100_000_000
        cluster1     = 0
        cluster2     = 0
        density      = 0

        [*(0..(clusters.length - 2))].each do |i|
          [*((i + 1)..(clusters.length - 1))].each do |j|
            dist = clusters[i].distance(clusters[j], distance_method)
            warn "distance between clusters #{i} and #{j} is #{dist}" if debug
            current_density = clusters[i].density + clusters[j].density
            if dist < min_distance
              min_distance = dist
              cluster1     = i
              cluster2     = j
              density      = current_density
            elsif dist == min_distance && density < current_density
              cluster1 = i
              cluster2 = j
              density  = current_density
            end
          end
        end

        # merge clusters 'cluster1' and 'cluster2'
        warn "clusters to merge #{cluster1} and #{cluster2}" if debug

        clusters[cluster1].add(clusters[cluster2])
        clusters.delete_at(cluster2)

        if debug
          clusters.each_with_index do |elem, i|
            warn "cluster #{i}"
            elem.print
          end
        end

        # stop condition 3
        # the density of the biggest clusters exceeds the threshold
        if no_clusters == 0 && clusters[cluster].density > threshold_density
          break
        end
      end

      @clusters = clusters
    end