/**
 * Copyright (c) 2011, The University of Southampton and the individual contributors.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 *
 *   *  Redistributions of source code must retain the above copyright notice,
 *      this list of conditions and the following disclaimer.
 *
 *   *  Redistributions in binary form must reproduce the above copyright notice,
 *      this list of conditions and the following disclaimer in the documentation
 *      and/or other materials provided with the distribution.
 *
 *   *  Neither the name of the University of Southampton nor the names of its
 *      contributors may be used to endorse or promote products derived from this
 *      software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
package org.openimaj.ml.clustering.kmeans;

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

import org.openimaj.data.DataSource;
import org.openimaj.data.DoubleArrayBackedDataSource;
import org.openimaj.ml.clustering.IndexClusters;
import org.openimaj.ml.clustering.SpatialClusterer;
import org.openimaj.util.function.Operation;
import org.openimaj.util.parallel.Parallel;
import org.openimaj.util.parallel.Parallel.IntRange;

/**
 * Multithreaded (optionally) damped spherical k-means with support for
 * bigger-than-memory data.
 * <p>
 * Spherical K-Means uses the inner product as the similarity metric, and is
 * constrained to finding centroids that lie on the surface of the unit
 * hypersphere (i.e. their length is 1). More formally, it solves:
 * <p>
 * min_{D,s}(sum_i(||Ds^(i) - x^(i)||_2^2))
 * <p>
 * s.t. ||s^(i)||_0 <= 1, for all i and ||D^(j)||_2 = 1, for all i
 * <p>
 * where D is a dictionary of centroids (with unit length) and s is an indicator
 * vector that is all zeros, except for a non-zero value in the position
 * corresponding top the assigned centroid.
 * <p>
 * The optional damping operation includes the previous centroid position in the
 * update computation, ensuring smoother convergence.
 * <p>
 * This implementation performs initialisation by randomly sampling centroids
 * from a Gaussian distribution, and then normalising to unit length. Any
 * centroids that become empty during the iterations are replaced by a new
 * random centroid generated in the same manner.
 * <p>
 * This implementation is able to deal with larger-than-memory datasets by
 * streaming the samples from disk using an appropriate {@link DataSource}. The
 * only requirement is that there is enough memory to hold all the centroids
 * plus working memory for the batches of samples being assigned.
 *
 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
 *
 */
public class SphericalKMeans implements SpatialClusterer<SphericalKMeansResult, double[]> {
        /**
         * Object storing the result of the previous iteration of spherical kmeans.
         * The object should be considered to be immutable, and read only. The
         * kmeans implementation will reuse the same instance for each iteration.
         *
         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
         */
        public static class IterationResult {
                /**
                 * The iteration number, starting from 0
                 */
                public int iteration;
                /**
                 * The change in fitness from the previous iteration
                 */
                public double delta;
                /**
                 * The current results
                 */
                public SphericalKMeansResult result;
        }

        protected final Random rng = new Random();
        protected final boolean damped;
        protected final int maxIters;
        protected final int k;
        protected final double terminationEps = 0.1;
        protected List<Operation<IterationResult>> iterationListeners = new ArrayList<Operation<IterationResult>>(0);

        /**
         * Construct with the given parameters. Uses damped updates and terminates
         * if the change in fit between iterations is less than 0.1.
         *
         * @param k
         *            number of clusters
         * @param maxIters
         *            maximum number of iterations
         */
        public SphericalKMeans(int k, int maxIters) {
                this(k, maxIters, true);
        }

        /**
         * Construct with the given parameters. Terminates if the change in fit
         * between iterations is less than 0.1.
         *
         * @param k
         *            number of clusters
         * @param maxIters
         *            maximum number of iterations
         * @param damped
         *            use damped updates
         */
        public SphericalKMeans(int k, int maxIters, boolean damped) {
                this.k = k;
                this.maxIters = maxIters;
                this.damped = damped;
        }

        /**
         * Construct with the given parameters. Uses damped updates and terminates
         * if the change in fit between iterations is less than 0.1 or 10 iterations
         * is reached.
         *
         * @param k
         *            number of clusters
         */
        public SphericalKMeans(int k) {
                this(k, 10);
        }

        private void makeRandomCentroid(double[] ds) {
                double sumsq = 0;
                for (int i = 0; i < ds.length; i++) {
                        ds[i] = rng.nextGaussian();
                        sumsq += (ds[i] * ds[i]);
                }

                sumsq = 1 / Math.sqrt(sumsq);
                for (int i = 0; i < ds.length; i++) {
                        ds[i] *= sumsq;
                }
        }

        double performIteration(final DataSource<double[]> data, final SphericalKMeansResult result) {
                final int[] clusterSizes = new int[result.centroids.length];
                final double[][] newCentroids = new double[result.centroids.length][result.centroids[0].length];

                final double[] delta = { 0 };

                // perform the assignments
                Parallel.forRange(0, data.size(), 1, new Operation<Parallel.IntRange>() {
                        @Override
                        public void perform(IntRange range) {
                                for (int i = range.start; i < range.stop; i++) {
                                        final double[] vector = data.getData(i);
                                        double assignmentWeight = Double.MIN_VALUE;
                                        for (int j = 0; j < result.centroids.length; j++) {
                                                final double[] centroid = result.centroids[j];

                                                double dp = 0;
                                                for (int k = 0; k < centroid.length; k++) {
                                                        dp += centroid[k] * vector[k];
                                                }

                                                if (dp > assignmentWeight) {
                                                        assignmentWeight = dp;
                                                        result.assignments[i] = j;
                                                }
                                        }

                                        // aggregate the assignments to the relevant cluster
                                        synchronized (newCentroids) {
                                                clusterSizes[result.assignments[i]]++;
                                                delta[0] += assignmentWeight;
                                                for (int k = 0; k < newCentroids[0].length; k++) {
                                                        newCentroids[result.assignments[i]][k] += vector[k];
                                                }
                                        }
                                }
                        }
                });

                // update the centroids
                Parallel.forRange(0, result.centroids.length, 1, new Operation<Parallel.IntRange>() {
                        @Override
                        public void perform(IntRange range) {
                                for (int j = range.start; j < range.stop; j++) {
                                        if (clusterSizes[j] == 0) {
                                                // reinit to random vector
                                                makeRandomCentroid(result.centroids[j]);
                                        } else {
                                                final double[] centroid = result.centroids[j];
                                                final double[] ncentroid = newCentroids[j];

                                                double norm = 0;
                                                if (damped) {
                                                        for (int k = 0; k < centroid.length; k++) {
                                                                centroid[k] += ncentroid[k];
                                                                norm += centroid[k] * centroid[k];
                                                        }
                                                } else {
                                                        for (int k = 0; k < centroid.length; k++) {
                                                                centroid[k] = ncentroid[k];
                                                                norm += centroid[k] * centroid[k];
                                                        }
                                                }
                                                norm = 1.0 / Math.sqrt(norm);

                                                for (int k = 0; k < ncentroid.length; k++) {
                                                        centroid[k] *= norm;
                                                }
                                        }
                                }
                        }
                });

                return delta[0];
        }

        @Override
        public int[][] performClustering(double[][] data) {
                return new IndexClusters(cluster(data).assignments).clusters();
        }

        @Override
        public SphericalKMeansResult cluster(double[][] data) {
                return cluster(new DoubleArrayBackedDataSource(data));
        }

        @Override
        public SphericalKMeansResult cluster(DataSource<double[]> data) {
                final IterationResult ir = new IterationResult();
                ir.result = new SphericalKMeansResult();
                ir.result.centroids = new double[k][data.numDimensions()];
                ir.result.assignments = new int[data.size()];

                for (int j = 0; j < ir.result.centroids.length; j++) {
                        makeRandomCentroid(ir.result.centroids[j]);
                }

                double last = 0;
                for (ir.iteration = 0; ir.iteration < maxIters; ir.iteration++) {
                        for (final Operation<IterationResult> l : iterationListeners)
                                l.perform(ir);

                        final double d = performIteration(data, ir.result);
                        ir.delta = d - last;

                        if (ir.delta < terminationEps)
                                break;
                        last = d;
                }

                return ir.result;
        }

        /**
         * Add a listener that will be called before every iteration.
         *
         * @param op
         *            the listener
         */
        public void addIterationListener(Operation<IterationResult> op) {
                iterationListeners.add(op);
        }
}