/*
    * 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.commons.math4.legacy.ml.clustering;
  
  import java.util.ArrayList;
  import java.util.Arrays;
  import java.util.Collection;
  import java.util.List;
  
  import org.apache.commons.rng.UniformRandomProvider;
  import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  import org.apache.commons.math4.legacy.ml.distance.DistanceMeasure;
  import org.apache.commons.math4.legacy.stat.descriptive.moment.VectorialMean;
  
  /**
   * Implementation of k-means++ algorithm.
   * It is based on
   * <blockquote>
   *  Elkan, Charles.
   *  "Using the triangle inequality to accelerate k-means."
   *  ICML. Vol. 3. 2003.
   * </blockquote>
   *
   * <p>
   * Algorithm uses triangle inequality to speed up computation, by reducing
   * the amount of distances calculations.  Towards the last iterations of
   * the algorithm, points which already assigned to some cluster are unlikely
   * to move to a new cluster; updates of cluster centers are also usually
   * relatively small.
   * Triangle inequality is thus used to determine the cases where distance
   * computation could be skipped since center move only a little, without
   * affecting points partitioning.
   *
   * <p>
   * For initial centers seeding, we apply the algorithm described in
   * <blockquote>
   *  Arthur, David, and Sergei Vassilvitskii.
   *  "k-means++: The advantages of careful seeding."
   *  Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms.
   *  Society for Industrial and Applied Mathematics, 2007.
   * </blockquote>
   *
   * @param <T> Type of the points to cluster.
   */
  public class ElkanKMeansPlusPlusClusterer<T extends Clusterable>
      extends KMeansPlusPlusClusterer<T> {
  
      /**
       * @param k Clustering parameter.
       */
      public ElkanKMeansPlusPlusClusterer(int k) {
          super(k);
      }
  
      /**
       * @param k Clustering parameter.
       * @param maxIterations Allowed number of iterations.
       * @param measure Distance measure.
       * @param random Random generator.
       */
      public ElkanKMeansPlusPlusClusterer(int k,
                                          int maxIterations,
                                          DistanceMeasure measure,
                                          UniformRandomProvider random) {
          super(k, maxIterations, measure, random);
      }
  
      /**
       * @param k Clustering parameter.
       * @param maxIterations Allowed number of iterations.
       * @param measure Distance measure.
       * @param random Random generator.
       * @param emptyStrategy Strategy for handling empty clusters that
       * may appear during algorithm progress.
       */
      public ElkanKMeansPlusPlusClusterer(int k,
                                          int maxIterations,
                                          DistanceMeasure measure,
                                          UniformRandomProvider random,
                                          EmptyClusterStrategy emptyStrategy) {
          super(k, maxIterations, measure, random, emptyStrategy);
      }
  
      /** {@inheritDoc} */
     @Override
     public List<CentroidCluster<T>> cluster(final Collection<T> points) {
         final int k = getNumberOfClusters();
 
         // Number of clusters has to be smaller or equal the number of data points.
         if (points.size() < k) {
             throw new NumberIsTooSmallException(points.size(), k, false);
         }
 
         final List<T> pointsList = new ArrayList<>(points);
         final int n = points.size();
         final int dim = pointsList.get(0).getPoint().length;
 
         // Keep minimum intra cluster distance, e.g. for given cluster c s[c] is
         // the distance to the closest cluster c' or s[c] = 1/2 * min_{c'!=c} dist(c', c)
         final double[] s = new double[k];
         Arrays.fill(s, Double.MAX_VALUE);
         // Store the matrix of distances between all cluster centers, e.g. dcc[c1][c2] = distance(c1, c2)
         final double[][] dcc = new double[k][k];
 
         // For each point keeps the upper bound distance to the cluster center.
         final double[] u = new double[n];
         Arrays.fill(u, Double.MAX_VALUE);
 
         // For each point and for each cluster keeps the lower bound for the distance between the point and cluster
         final double[][] l = new double[n][k];
 
         // Seed initial set of cluster centers.
         final double[][] centers = seed(pointsList);
 
         // Points partitioning induced by cluster centers, e.g. for point xi the value of partitions[xi] indicates
         // the cluster or index of the cluster center which is closest to xi. partitions[xi] = min_{c} distance(xi, c).
         final int[] partitions = partitionPoints(pointsList, centers, u, l);
 
         final double[] deltas = new double[k];
         VectorialMean[] means = new VectorialMean[k];
         for (int it = 0, max = getMaxIterations();
              it < max;
              it++) {
             int changes = 0;
             // Step I.
             // Compute inter-cluster distances.
             updateIntraCentersDistances(centers, dcc, s);
 
             for (int xi = 0; xi < n; xi++) {
                 boolean r = true;
 
                 // Step II.
                 if (u[xi] <= s[partitions[xi]]) {
                     continue;
                 }
 
                 for (int c = 0; c < k; c++) {
                     // Check condition III.
                     if (isSkipNext(partitions, u, l, dcc, xi, c)) {
                         continue;
                     }
 
                     final double[] x = pointsList.get(xi).getPoint();
 
                     // III(a)
                     if (r) {
                         u[xi] = distance(x, centers[partitions[xi]]);
                         l[xi][partitions[xi]] = u[xi];
                         r = false;
                     }
                     // III(b)
                     if (u[xi] > l[xi][c] || u[xi] > dcc[partitions[xi]][c]) {
                         l[xi][c] = distance(x, centers[c]);
                         if (l[xi][c] < u[xi]) {
                             partitions[xi] = c;
                             u[xi] = l[xi][c];
                             ++changes;
                         }
                     }
                 }
             }
 
             // Stopping criterion.
             if (changes == 0 &&
                 it != 0) { // First iteration needed (to update bounds).
                 break;
             }
 
             // Step IV.
             Arrays.fill(means, null);
             for (int i = 0; i < n; i++) {
                 if (means[partitions[i]] == null) {
                     means[partitions[i]] = new VectorialMean(dim);
                 }
                 means[partitions[i]].increment(pointsList.get(i).getPoint());
             }
 
             for (int i = 0; i < k; i++) {
                 deltas[i] = distance(centers[i], means[i].getResult());
                 centers[i] = means[i].getResult();
             }
 
             updateBounds(partitions, u, l, deltas);
         }
 
         return buildResults(pointsList, partitions, centers);
     }
 
     /**
      * kmeans++ seeding which provides guarantee of resulting with log(k) approximation
      * for final clustering results
      * <p>
      * Arthur, David, and Sergei Vassilvitskii. "k-means++: The advantages of careful seeding."
      * Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms.
      * Society for Industrial and Applied Mathematics, 2007.
      *
      * @param points input data points
      * @return an array of initial clusters centers
      *
      */
     private double[][] seed(final List<T> points) {
         final int k = getNumberOfClusters();
         final UniformRandomProvider random = getRandomGenerator();
 
         final double[][] result = new double[k][];
         final int n = points.size();
         final int pointIndex = random.nextInt(n);
 
         final double[] minDistances = new double[n];
 
         int idx = 0;
         result[idx] = points.get(pointIndex).getPoint();
 
         double sumSqDist = 0;
 
         for (int i = 0; i < n; i++) {
             final double d = distance(result[idx], points.get(i).getPoint());
             minDistances[i] = d * d;
             sumSqDist += minDistances[i];
         }
 
         while (++idx < k) {
             final double p = sumSqDist * random.nextDouble();
             int next = 0;
             for (double cdf = 0; cdf < p; next++) {
                 cdf += minDistances[next];
             }
 
             result[idx] = points.get(next - 1).getPoint();
             for (int i = 0; i < n; i++) {
                 final double d = distance(result[idx], points.get(i).getPoint());
                 sumSqDist -= minDistances[i];
                 minDistances[i] = Math.min(minDistances[i], d * d);
                 sumSqDist += minDistances[i];
             }
         }
 
         return result;
     }
 
 
     /**
      * Once initial centers are chosen, we can actually go through data points and assign points to the
      * cluster based on the distance between initial centers and points.
      *
      * @param pointsList data points list
      * @param centers current clusters centers
      * @param u points upper bounds
      * @param l lower bounds for points to clusters centers
      *
      * @return initial assignment of points into clusters
      */
     private int[] partitionPoints(List<T> pointsList,
                                   double[][] centers,
                                   double[] u,
                                   double[][] l) {
         final int k = getNumberOfClusters();
         final int n = pointsList.size();
         // Points assignments vector.
         final int[] assignments = new int[n];
         Arrays.fill(assignments, -1);
         // Need to assign points to the clusters for the first time and intitialize the lower bound l(x, c)
         for (int i = 0; i < n; i++) {
             final double[] x = pointsList.get(i).getPoint();
             for (int j = 0; j < k; j++) {
                 l[i][j] = distance(x, centers[j]); // l(x, c) = d(x, c)
                 if (u[i] > l[i][j]) {
                     u[i] = l[i][j]; // u(x) = min_c d(x, c)
                     assignments[i] = j; // c(x) = argmin_c d(x, c)
                 }
             }
         }
         return assignments;
     }
 
     /**
      * Updated distances between clusters centers and for each cluster
      * pick the closest neighbour and keep distance to it.
      *
      * @param centers cluster centers
      * @param dcc matrix of distance between clusters centers, e.g.
      * {@code dcc[i][j] = distance(centers[i], centers[j])}
      * @param s For a given cluster, {@code s[si]} holds distance value
      * to the closest cluster center.
      */
     private void updateIntraCentersDistances(double[][] centers,
                                              double[][] dcc,
                                              double[] s) {
         final int k = getNumberOfClusters();
         for (int i = 0; i < k; i++) {
             // Since distance(xi, xj) == distance(xj, xi), we need to update
             // only upper or lower triangle of the distances matrix and mirror
             // to the lower of upper triangle accordingly, trace has to be all
             // zeros, since distance(xi, xi) == 0.
             for (int j = i + 1; j < k; j++) {
                 dcc[i][j] = 0.5 * distance(centers[i], centers[j]);
                 dcc[j][i] = dcc[i][j];
                 if (dcc[i][j] < s[i]) {
                     s[i] = dcc[i][j];
                 }
                 if (dcc[j][i] < s[j]) {
                     s[j] = dcc[j][i];
                 }
             }
         }
     }
 
     /**
      * For given points and and cluster, check condition (3) of Elkan algorithm.
      *
      * <ul>
      *  <li>c is not the cluster xi assigned to</li>
      *  <li>{@code u[xi] > l[xi][x]} upper bound for point xi is greater than
      *   lower bound between xi and some cluster c</li>
      *  <li>{@code u[xi] > 1/2 * d(c(xi), c)} upper bound is greater than
      *   distance between center of xi's cluster and c</li>
      * </ul>
      *
      * @param partitions current partition of points into clusters
      * @param u upper bounds for points
      * @param l lower bounds for distance between cluster centers and points
      * @param dcc matrix of distance between clusters centers
      * @param xi index of the point
      * @param c index of the cluster
      * @return true if conditions above satisfied false otherwise
      */
     private static boolean isSkipNext(int[] partitions,
                                       double[] u,
                                       double[][] l,
                                       double[][] dcc,
                                       int xi,
                                       int c) {
         return c == partitions[xi] ||
                u[xi] <= l[xi][c] ||
                u[xi] <= dcc[partitions[xi]][c];
     }
 
     /**
      * Once kmeans iterations have been converged and no more movements, we can build up the final
      * resulted list of cluster centroids ({@link CentroidCluster}) and assign input points based
      * on the converged partitioning.
      *
      * @param pointsList list of data points
      * @param partitions current partition of points into clusters
      * @param centers cluster centers
      * @return cluster partitioning
      */
     private List<CentroidCluster<T>> buildResults(List<T> pointsList,
                                                   int[] partitions,
                                                   double[][] centers) {
         final int k = getNumberOfClusters();
         final List<CentroidCluster<T>> result = new ArrayList<>();
         for (int i = 0; i < k; i++) {
             final CentroidCluster<T> cluster = new CentroidCluster<>(new DoublePoint(centers[i]));
             result.add(cluster);
         }
         for (int i = 0; i < pointsList.size(); i++) {
             result.get(partitions[i]).addPoint(pointsList.get(i));
         }
         return result;
     }
 
     /**
      * Based on the distance that cluster center has moved we need to update our upper and lower bound.
      * Worst case assumption, the center of the assigned to given cluster moves away from the point, while
      * centers of over clusters become closer.
      *
      * @param partitions current points assiments to the clusters
      * @param u points upper bounds
      * @param l lower bounds for distances between point and corresponding cluster
      * @param deltas the movement delta for each cluster center
      */
     private void updateBounds(int[] partitions,
                               double[] u,
                               double[][] l,
                               double[] deltas) {
         final int k = getNumberOfClusters();
         for (int i = 0; i < partitions.length; i++) {
             u[i] += deltas[partitions[i]];
             for (int j = 0; j < k; j++) {
                 l[i][j] = Math.max(0, l[i][j] - deltas[j]);
             }
         }
     }
 
     /**
      * @param a Coordinates.
      * @param b Coordinates.
      * @return the distance between {@code a} and {@code b}.
      */
     private double distance(final double[] a,
                             final double[] b) {
         return getDistanceMeasure().compute(a, b);
     }
 }