/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.mahout.clustering.streaming.cluster;

import java.util.ArrayList;
import java.util.Collections;
import java.util.Iterator;
import java.util.List;
import java.util.Random;

import com.google.common.base.Function;
import com.google.common.base.Preconditions;
import com.google.common.collect.Iterables;
import com.google.common.collect.Iterators;
import org.apache.mahout.clustering.ClusteringUtils;
import org.apache.mahout.common.Pair;
import org.apache.mahout.common.RandomUtils;
import org.apache.mahout.common.distance.DistanceMeasure;
import org.apache.mahout.math.Centroid;
import org.apache.mahout.math.Vector;
import org.apache.mahout.math.WeightedVector;
import org.apache.mahout.math.neighborhood.UpdatableSearcher;
import org.apache.mahout.math.random.Multinomial;
import org.apache.mahout.math.random.WeightedThing;

/**
 * Implements a ball k-means algorithm for weighted vectors with probabilistic seeding similar to k-means++.
 * The idea is that k-means++ gives good starting clusters and ball k-means can tune up the final result very nicely
 * in only a few passes (or even in a single iteration for well-clusterable data).
 *
 * A good reference for this class of algorithms is "The Effectiveness of Lloyd-Type Methods for the k-Means Problem"
 * by Rafail Ostrovsky, Yuval Rabani, Leonard J. Schulman and Chaitanya Swamy.  The code here uses the seeding strategy
 * as described in section 4.1.1 of that paper and the ball k-means step as described in section 4.2.  We support
 * multiple iterations in contrast to the algorithm described in the paper.
 */
public class BallKMeans implements Iterable<Centroid> {
  /**
   * The searcher containing the centroids.
   */
  private final UpdatableSearcher centroids;

  /**
   * The number of clusters to cluster the data into.
   */
  private final int numClusters;

  /**
   * The maximum number of iterations of the algorithm to run waiting for the cluster assignments
   * to stabilize. If there are no changes in cluster assignment earlier, we can finish early.
   */
  private final int maxNumIterations;

  /**
   * When deciding which points to include in the new centroid calculation,
   * it's preferable to exclude outliers since it increases the rate of convergence.
   * So, we calculate the distance from each cluster to its closest neighboring cluster. When
   * evaluating the points assigned to a cluster, we compare the distance between the centroid to
   * the point with the distance between the centroid and its closest centroid neighbor
   * multiplied by this trimFraction. If the distance between the centroid and the point is
   * greater, we consider it an outlier and we don't use it.
   */
  private final double trimFraction;

  /**
   * Selecting the initial centroids is the most important part of the ball k-means clustering. Poor choices, like two
   * centroids in the same actual cluster result in a low-quality final result.
   * k-means++ initialization yields good quality clusters, especially when using BallKMeans after StreamingKMeans as
   * the points have weights.
   * Simple, random selection of the points based on their weights is faster but sometimes fails to produce the
   * desired number of clusters.
   * This field is true if the initialization should be done with k-means++.
   */
  private final boolean kMeansPlusPlusInit;

  /**
   * When using trimFraction, the weight of each centroid will not be the sum of the weights of
   * the vectors assigned to that cluster because outliers are not used to compute the updated
   * centroid.
   * So, the total weight is probably wrong. This can be fixed by doing another pass over the
   * data points and adjusting the weights of each centroid. This doesn't update the coordinates
   * of the centroids, but is useful if the weights matter.
   */
  private final boolean correctWeights;

  /**
   * When running multiple ball k-means passes to get the one with the smallest total cost, can compute the
   * overall cost, using all the points for clustering, or reserve a fraction of them, testProbability in a test set.
   * The cost is the sum of the distances between each point and its corresponding centroid.
   * We then use this set of points to compute the total cost on. We're therefore trying to select the clustering
   * that best describes the underlying distribution of the clusters.
   * This field is the probability of assigning a given point to the test set. If this is 0, the cost will be computed
   * on the entire set of points.
   */
  private final double testProbability;

  /**
   * Whether or not testProbability > 0, i.e., there exists a non-empty 'test' set.
   */
  private final boolean splitTrainTest;

  /**
   * How many k-means runs to have. If there's more than one run, we compute the cost of each clustering as described
   * above and select the clustering that minimizes the cost.
   * Multiple runs are a lot more useful when using the random initialization. With kmeans++, 1-2 runs are enough and
   * more runs are not likely to help quality much.
   */
  private final int numRuns;

  /**
   * Random object to sample values from.
   */
  private final Random random;

  public BallKMeans(UpdatableSearcher searcher, int numClusters, int maxNumIterations) {
    // By default, the trimFraction is 0.9, k-means++ is used, the weights will be corrected at the end,
    // there will be 0 points in the test set and 1 run.
    this(searcher, numClusters, maxNumIterations, 0.9, true, true, 0.0, 1);
  }

  public BallKMeans(UpdatableSearcher searcher, int numClusters, int maxNumIterations,
                    boolean kMeansPlusPlusInit, int numRuns) {
    // By default, the trimFraction is 0.9, k-means++ is used, the weights will be corrected at the end,
    // there will be 10% points of in the test set.
    this(searcher, numClusters, maxNumIterations, 0.9, kMeansPlusPlusInit, true, 0.1, numRuns);
  }

  public BallKMeans(UpdatableSearcher searcher, int numClusters, int maxNumIterations,
                    double trimFraction, boolean kMeansPlusPlusInit, boolean correctWeights,
                    double testProbability, int numRuns) {
    Preconditions.checkArgument(searcher.size() == 0, "Searcher must be empty initially to populate with centroids");
    Preconditions.checkArgument(numClusters > 0, "The requested number of clusters must be positive");
    Preconditions.checkArgument(maxNumIterations > 0, "The maximum number of iterations must be positive");
    Preconditions.checkArgument(trimFraction > 0, "The trim fraction must be positive");
    Preconditions.checkArgument(testProbability >= 0 && testProbability < 1, "The testProbability must be in [0, 1)");
    Preconditions.checkArgument(numRuns > 0, "There has to be at least one run");

    this.centroids = searcher;
    this.numClusters = numClusters;
    this.maxNumIterations = maxNumIterations;

    this.trimFraction = trimFraction;
    this.kMeansPlusPlusInit = kMeansPlusPlusInit;
    this.correctWeights = correctWeights;

    this.testProbability = testProbability;
    this.splitTrainTest = testProbability > 0;
    this.numRuns = numRuns;

    this.random = RandomUtils.getRandom();
  }

  public Pair<List<? extends WeightedVector>, List<? extends WeightedVector>> splitTrainTest(
      List<? extends WeightedVector> datapoints) {
    // If there will be no points assigned to the test set, return now.
    if (testProbability == 0) {
      return new Pair<List<? extends WeightedVector>, List<? extends WeightedVector>>(datapoints,
          new ArrayList<WeightedVector>());
    }

    int numTest = (int) (testProbability * datapoints.size());
    Preconditions.checkArgument(numTest > 0 && numTest < datapoints.size(),
        "Must have nonzero number of training and test vectors. Asked for %.1f %% of %d vectors for test",
        testProbability * 100, datapoints.size());

    Collections.shuffle(datapoints);
    return new Pair<List<? extends WeightedVector>, List<? extends WeightedVector>>(
        datapoints.subList(numTest, datapoints.size()), datapoints.subList(0, numTest));
  }

  /**
   * Clusters the datapoints in the list doing either random seeding of the centroids or k-means++.
   *
   * @param datapoints the points to be clustered.
   * @return an UpdatableSearcher with the resulting clusters.
   */
  public UpdatableSearcher cluster(List<? extends WeightedVector> datapoints) {
    Pair<List<? extends WeightedVector>, List<? extends WeightedVector>> trainTestSplit = splitTrainTest(datapoints);
    List<Vector> bestCentroids = new ArrayList<>();
    double cost = Double.POSITIVE_INFINITY;
    double bestCost = Double.POSITIVE_INFINITY;
    for (int i = 0; i < numRuns; ++i) {
      centroids.clear();
      if (kMeansPlusPlusInit) {
        // Use k-means++ to set initial centroids.
        initializeSeedsKMeansPlusPlus(trainTestSplit.getFirst());
      } else {
        // Randomly select the initial centroids.
        initializeSeedsRandomly(trainTestSplit.getFirst());
      }
      // Do k-means iterations with trimmed mean computation (aka ball k-means).
      if (numRuns > 1) {
        // If the clustering is successful (there are no zero-weight centroids).
        iterativeAssignment(trainTestSplit.getFirst());
        // Compute the cost of the clustering and possibly save the centroids.
        cost = ClusteringUtils.totalClusterCost(
            splitTrainTest ? datapoints : trainTestSplit.getSecond(), centroids);
        if (cost < bestCost) {
          bestCost = cost;
          bestCentroids.clear();
          Iterables.addAll(bestCentroids, centroids);
        }
      } else {
        // If there is only going to be one run, the cost doesn't need to be computed, so we just return the clustering.
        iterativeAssignment(datapoints);
        return centroids;
      }
    }
    if (bestCost == Double.POSITIVE_INFINITY) {
      throw new RuntimeException("No valid clustering was found");
    }
    if (cost != bestCost) {
      centroids.clear();
      centroids.addAll(bestCentroids);
    }
    if (correctWeights) {
      for (WeightedVector testDatapoint : trainTestSplit.getSecond()) {
        WeightedVector closest = (WeightedVector) centroids.searchFirst(testDatapoint, false).getValue();
        closest.setWeight(closest.getWeight() + testDatapoint.getWeight());
      }
    }
    return centroids;
  }

  /**
   * Selects some of the original points randomly with probability proportional to their weights. This is much
   * less sophisticated than the kmeans++ approach, however it is faster and coupled with
   *
   * The side effect of this method is to fill the centroids structure itself.
   *
   * @param datapoints The datapoints to select from.  These datapoints should be WeightedVectors of some kind.
   */
  private void initializeSeedsRandomly(List<? extends WeightedVector> datapoints) {
    int numDatapoints = datapoints.size();
    double totalWeight = 0;
    for (WeightedVector datapoint : datapoints) {
      totalWeight += datapoint.getWeight();
    }
    Multinomial<Integer> seedSelector = new Multinomial<>();
    for (int i = 0; i < numDatapoints; ++i) {
      seedSelector.add(i, datapoints.get(i).getWeight() / totalWeight);
    }
    for (int i = 0; i < numClusters; ++i) {
      int sample = seedSelector.sample();
      seedSelector.delete(sample);
      Centroid centroid = new Centroid(datapoints.get(sample));
      centroid.setIndex(i);
      centroids.add(centroid);
    }
  }

  /**
   * Selects some of the original points according to the k-means++ algorithm.  The basic idea is that
   * points are selected with probability proportional to their distance from any selected point.  In
   * this version, points have weights which multiply their likelihood of being selected.  This is the
   * same as if there were as many copies of the same point as indicated by the weight.
   *
   * This is pretty expensive, but it vastly improves the quality and convergences of the k-means algorithm.
   * The basic idea can be made much faster by only processing a random subset of the original points.
   * In the context of streaming k-means, the total number of possible seeds will be about k log n so this
   * selection will cost O(k^2 (log n)^2) which isn't much worse than the random sampling idea.  At
   * n = 10^9, the cost of this initialization will be about 10x worse than a reasonable random sampling
   * implementation.
   *
   * The side effect of this method is to fill the centroids structure itself.
   *
   * @param datapoints The datapoints to select from.  These datapoints should be WeightedVectors of some kind.
   */
  private void initializeSeedsKMeansPlusPlus(List<? extends WeightedVector> datapoints) {
    Preconditions.checkArgument(datapoints.size() > 1, "Must have at least two datapoints points to cluster " +
        "sensibly");
    Preconditions.checkArgument(datapoints.size() >= numClusters,
        String.format("Must have more datapoints [%d] than clusters [%d]", datapoints.size(), numClusters));
    // Compute the centroid of all of the datapoints.  This is then used to compute the squared radius of the datapoints.
    Centroid center = new Centroid(datapoints.iterator().next());
    for (WeightedVector row : Iterables.skip(datapoints, 1)) {
      center.update(row);
    }

    // Given the centroid, we can compute \Delta_1^2(X), the total squared distance for the datapoints
    // this accelerates seed selection.
    double deltaX = 0;
    DistanceMeasure distanceMeasure = centroids.getDistanceMeasure();
    for (WeightedVector row : datapoints) {
      deltaX += distanceMeasure.distance(row, center);
    }

    // Find the first seed c_1 (and conceptually the second, c_2) as might be done in the 2-means clustering so that
    // the probability of selecting c_1 and c_2 is proportional to || c_1 - c_2 ||^2.  This is done
    // by first selecting c_1 with probability:
    //
    // p(c_1) = sum_{c_1} || c_1 - c_2 ||^2 \over sum_{c_1, c_2} || c_1 - c_2 ||^2
    //
    // This can be simplified to:
    //
    // p(c_1) = \Delta_1^2(X) + n || c_1 - c ||^2 / (2 n \Delta_1^2(X))
    //
    // where c = \sum x / n and \Delta_1^2(X) = sum || x - c ||^2
    //
    // All subsequent seeds c_i (including c_2) can then be selected from the remaining points with probability
    // proportional to Pr(c_i == x_j) = min_{m < i} || c_m - x_j ||^2.

    // Multinomial distribution of vector indices for the selection seeds. These correspond to
    // the indices of the vectors in the original datapoints list.
    Multinomial<Integer> seedSelector = new Multinomial<>();
    for (int i = 0; i < datapoints.size(); ++i) {
      double selectionProbability =
          deltaX + datapoints.size() * distanceMeasure.distance(datapoints.get(i), center);
      seedSelector.add(i, selectionProbability);
    }

    int selected = random.nextInt(datapoints.size());
    Centroid c_1 = new Centroid(datapoints.get(selected).clone());
    c_1.setIndex(0);
    // Construct a set of weighted things which can be used for random selection.  Initial weights are
    // set to the squared distance from c_1
    for (int i = 0; i < datapoints.size(); ++i) {
      WeightedVector row = datapoints.get(i);
      double w = distanceMeasure.distance(c_1, row) * 2 * Math.log(1 + row.getWeight());
      seedSelector.set(i, w);
    }

    // From here, seeds are selected with probability proportional to:
    //
    // r_i = min_{c_j} || x_i - c_j ||^2
    //
    // when we only have c_1, we have already set these distances and as we select each new
    // seed, we update the minimum distances.
    centroids.add(c_1);
    int clusterIndex = 1;
    while (centroids.size() < numClusters) {
      // Select according to weights.
      int seedIndex = seedSelector.sample();
      Centroid nextSeed = new Centroid(datapoints.get(seedIndex));
      nextSeed.setIndex(clusterIndex++);
      centroids.add(nextSeed);
      // Don't select this one again.
      seedSelector.delete(seedIndex);
      // Re-weight everything according to the minimum distance to a seed.
      for (int currSeedIndex : seedSelector) {
        WeightedVector curr = datapoints.get(currSeedIndex);
        double newWeight = nextSeed.getWeight() * distanceMeasure.distance(nextSeed, curr);
        if (newWeight < seedSelector.getWeight(currSeedIndex)) {
          seedSelector.set(currSeedIndex, newWeight);
        }
      }
    }
  }

  /**
   * Examines the datapoints and updates cluster centers to be the centroid of the nearest datapoints points.  To
   * compute a new center for cluster c_i, we average all points that are closer than d_i * trimFraction
   * where d_i is
   *
   * d_i = min_j \sqrt ||c_j - c_i||^2
   *
   * By ignoring distant points, the centroids converge more quickly to a good approximation of the
   * optimal k-means solution (given good starting points).
   *
   * @param datapoints the points to cluster.
   */
  private void iterativeAssignment(List<? extends WeightedVector> datapoints) {
    DistanceMeasure distanceMeasure = centroids.getDistanceMeasure();
    // closestClusterDistances.get(i) is the distance from the i'th cluster to its closest
    // neighboring cluster.
    List<Double> closestClusterDistances = new ArrayList<>(numClusters);
    // clusterAssignments[i] == j means that the i'th point is assigned to the j'th cluster. When
    // these don't change, we are done.
    // Each point is assigned to the invalid "-1" cluster initially.
    List<Integer> clusterAssignments = new ArrayList<>(Collections.nCopies(datapoints.size(), -1));

    boolean changed = true;
    for (int i = 0; changed && i < maxNumIterations; i++) {
      changed = false;
      // We compute what the distance between each cluster and its closest neighbor is to set a
      // proportional distance threshold for points that should be involved in calculating the
      // centroid.
      closestClusterDistances.clear();
      for (Vector center : centroids) {
        // If a centroid has no points assigned to it, the clustering failed.
        Vector closestOtherCluster = centroids.searchFirst(center, true).getValue();
        closestClusterDistances.add(distanceMeasure.distance(center, closestOtherCluster));
      }

      // Copies the current cluster centroids to newClusters and sets their weights to 0. This is
      // so we calculate the new centroids as we go through the datapoints.
      List<Centroid> newCentroids = new ArrayList<>();
      for (Vector centroid : centroids) {
        // need a deep copy because we will mutate these values
        Centroid newCentroid = (Centroid)centroid.clone();
        newCentroid.setWeight(0);
        newCentroids.add(newCentroid);
      }

      // Pass over the datapoints computing new centroids.
      for (int j = 0; j < datapoints.size(); ++j) {
        WeightedVector datapoint = datapoints.get(j);
        // Get the closest cluster this point belongs to.
        WeightedThing<Vector> closestPair = centroids.searchFirst(datapoint, false);
        int closestIndex = ((WeightedVector) closestPair.getValue()).getIndex();
        double closestDistance = closestPair.getWeight();
        // Update its cluster assignment if necessary.
        if (closestIndex != clusterAssignments.get(j)) {
          changed = true;
          clusterAssignments.set(j, closestIndex);
        }
        // Only update if the datapoints point is near enough. What this means is that the weight
        // of outliers is NOT taken into account and the final weights of the centroids will
        // reflect this (it will be less or equal to the initial sum of the weights).
        if (closestDistance < trimFraction * closestClusterDistances.get(closestIndex)) {
          newCentroids.get(closestIndex).update(datapoint);
        }
      }
      // Add the new centers back into searcher.
      centroids.clear();
      centroids.addAll(newCentroids);
    }

    if (correctWeights) {
      for (Vector v : centroids) {
        ((Centroid)v).setWeight(0);
      }
      for (WeightedVector datapoint : datapoints) {
        Centroid closestCentroid = (Centroid) centroids.searchFirst(datapoint, false).getValue();
        closestCentroid.setWeight(closestCentroid.getWeight() + datapoint.getWeight());
      }
    }
  }

  @Override
  public Iterator<Centroid> iterator() {
    return Iterators.transform(centroids.iterator(), new Function<Vector, Centroid>() {
      @Override
      public Centroid apply(Vector input) {
        Preconditions.checkArgument(input instanceof Centroid, "Non-centroid in centroids " +
            "searcher");
        //noinspection ConstantConditions
        return (Centroid)input;
      }
    });
  }
}
