using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using ClusterAnalyze;

namespace _160240
{
    class Program
    {
        static void Main(string[] args)
        {
            uint rows, cols;
            rows = InputUInt("Введите количество строк матрицы:");
            cols = InputUInt("Введите количество столбцов матрицы:");
            double[,] matrix;
            try
            {
                matrix = new double[rows, cols];
                Random rndGen = new Random();
                for (uint i = 0; i < rows; i++)
                {
                    for (uint j = 0; j < cols; j++)
                    {
                        matrix[i, j] = rndGen.NextDouble() * 100 - 50;
                    }
                }
                Console.WriteLine("Исходная матрица:");
                for (uint i = 0; i < rows; i++)
                {
                    for (uint j = 0; j < cols; j++)
                    {
                        Console.Write("{0,7:F3} ", matrix[i, j]);
                    }
                    Console.WriteLine();
                }
                double deviation = InputDouble("Введите желаемую погрешность:");
                try
                {
                    ClusterAnalyze.ClusterResult result = KMeans.Solve(matrix, rows, deviation);
                    uint clusterNo = 0;
                    foreach (ClusterAnalyze.ListOfRecords cluster in result)
                    {
                        Console.WriteLine("Кластер {0}", ++clusterNo);
                        foreach (uint line in cluster)
                        {
                            for (uint i = 0; i < cols; i++)
                            {
                                Console.Write("{0,7:F3} ", matrix[line, i]);
                            }
                            Console.WriteLine();
                        }
                    }
                    Console.WriteLine("Среднее отклонение:{0}", result.AverageDeviation);
                    Console.WriteLine("Максимальное отклонение:{0}", result.MaxDeviation);
                }
                catch (ClusterAnalyze.ClusterAnalyzeException ex)
                {
                    ErrorMsg(ex.Message);
                }
            }
            catch (OutOfMemoryException)
            {
                ErrorMsg("Не могу выделить память для матрицы");
            }
            Console.WriteLine("Нажмите любую клавишу для выхода");
            Console.ReadKey();
        }
        /// <summary>
        /// Ввод числа типа double
        /// </summary>
        /// <param name="msg">Сообщение</param>
        /// <returns>Число</returns>
        static double InputDouble(string msg)
        {
            while (true)
            {
                Console.Write(msg);
                try
                {
                    double val = Convert.ToDouble(Console.ReadLine());
                    return val;
                }
                catch (FormatException)
                { }
                catch (OverflowException)
                { }
                ErrorMsg("Неверный ввод: ожидается действительное число");
            }
        }
        /// <summary>
        /// Ввод числа типа uint
        /// </summary>
        /// <param name="msg">Сообщение</param>
        /// <returns>Число</returns>
        static uint InputUInt(string msg)
        {
            while (true)
            {
                Console.Write(msg);
                try
                {
                    uint val = Convert.ToUInt32(Console.ReadLine());
                    return val;
                }
                catch (FormatException)
                { }
                catch (OverflowException)
                { }
                ErrorMsg("Неверный ввод: ожидается целое положительное число");
            }
        }
        /// <summary>
        /// Ввод числа типа uint в заданном диапазоне
        /// </summary>
        /// <param name="msg">Сообщение</param>
        /// <param name="lo">Нижняя граница</param>
        /// <param name="high">Верхняя граница</param>
        /// <returns>Число</returns>
        static uint InputUInt(string msg, uint lo, uint high)
        {
            while (true)
            {
                uint val = InputUInt(msg);
                if (val >= lo && val <= high) return val;
                ErrorMsg(String.Format("Неверный ввод: ожидается значение в интервале {0}..{1}", lo, high));
            }
        }
        /// <summary>
        /// Выводит сообщение об ошибке
        /// </summary>
        /// <param name="msg">Сообщение</param>
        static void ErrorMsg(string msg)
        {
            ConsoleColor color = Console.ForegroundColor;
            Console.ForegroundColor = ConsoleColor.Red;
            Console.WriteLine(msg);
            Console.ForegroundColor = color;
        }
    }
}

namespace ClusterAnalyze
{
    /// <summary>
    /// Класс для кластерного анализа
    /// </summary>
    static class KMeans
    {
        /// <summary>
        /// Решает задачу кластеризации
        /// </summary>
        /// <param name="matrix">Исходная матрица</param>
        /// <param name="numCluster">Количество кластеров</param>
        /// <returns>Результат анализа</returns>
        public static ClusterResult Solve(double[,] matrix, uint numCluster)
        {
            return Clustering(Matrix2JadMatrix(matrix), numCluster);
        }
        /// <summary>
        /// Решает задачу кластеризации
        /// </summary>
        /// <param name="matrix">Исходная матрица</param>
        /// <param name="maxNumCluster">Ограничение на количество кластеров сверху</param>
        /// <param name="deviation">Допустимое отклонение записи от центроида</param>
        /// <returns>Результат анализа</returns>
        public static ClusterResult Solve(double[,] matrix, uint maxNumCluster, double deviation)
        {
            if (deviation == 0) deviation = 1e-10;
            else deviation = deviation < 0 ? -deviation : deviation;
            ClusterResult res;
            double[][] jadMatrix = Matrix2JadMatrix(matrix);
            uint clusters = 0;
            do
            {
                res = Clustering(jadMatrix, ++clusters);
            } while ((clusters < maxNumCluster) && (res.MaxDeviation > deviation));
            GC.Collect();
            return res;
        }
        /// <summary>
        /// Собственно решение здесь
        /// </summary>
        /// <param name="matrix">Матрица</param>
        /// <param name="numCluster">Количество кластеров</param>
        /// <returns>Результат анализа</returns>
        private static ClusterResult Clustering(double[][] matrix, uint numCluster)
        {
            uint rows = (uint)matrix.Length;
            // Количество кластеров должно быть не менее 1х и не более чем записей
            if (numCluster > rows || numCluster < 1)
                throw new ClusterAnalyzeException("ClusterResult Clustering(double[][] matrix, uint numCluster):\n" +
                    "Количество кластеров меньше 1 или больше чем число строк в матрице");
            uint cols = (uint)matrix[0].Length;
            // Задаем центроиды для начала
            double[][] centroids = InitialCentroids(matrix, numCluster);
            bool isCont;
            // Словарь ключями которого являются номера центроидов, а значениями списки записей в матрице
            Dictionary<uint, ListOfRecords> clustersDict;
            double averageDeviation, maxDeviation;
            do
            {
                clustersDict = new Dictionary<uint, ListOfRecords>();
                // Отклонения от центроидов
                averageDeviation = 0;
                maxDeviation = 0;
                // Для заданных центроидов определяем к какому кластеру относятся записи
                for (uint i = 0; i < rows; i++)
                {
                    // Индекс ближайшего
                    uint nearest = 0;
                    // Минимальная дистанция
                    double minDistance = Range(centroids[0], matrix[i]);
                    // Перебор центроидов
                    for (uint j = 1; j < numCluster; j++)
                    {
                        double distance = Range(centroids[j], matrix[i]);
                        if (distance < minDistance)
                        {
                            minDistance = distance;
                            nearest = j;
                        }
                    }
                    averageDeviation += minDistance / rows;
                    maxDeviation = maxDeviation < minDistance ? minDistance : maxDeviation;
                    // Добавляем в словарь
                    if (!clustersDict.ContainsKey(nearest))
                    {
                        clustersDict.Add(nearest, new ListOfRecords());
                    }
                    clustersDict[nearest].Add(i);
                }
                // Если есть не рабочие центроиды то избавимся от них
                if (clustersDict.Count != numCluster)
                {
                    numCluster = (uint)clustersDict.Count;
                    double[][] tmpCentroids = new double[numCluster][];
                    Dictionary<uint, ListOfRecords> tmpClusterDict = new Dictionary<uint, ListOfRecords>();
                    uint newNo = 0;
                    foreach (uint oldNo in clustersDict.Keys)
                    {
                        tmpCentroids[newNo] = centroids[oldNo];
                        tmpClusterDict.Add(newNo, clustersDict[oldNo]);
                        newNo++;
                    }
                    centroids = tmpCentroids;
                    clustersDict = tmpClusterDict;
                }
                // Флаг продолжения итераций
                isCont = false;
                // Считаем новые центроиды
                double[][] newCentroids = new double[numCluster][];
                foreach (uint i in clustersDict.Keys)
                {
                    newCentroids[i] = new double[cols];
                    uint itemsInCluster = (uint)clustersDict[i].Count;
                    foreach (uint line in clustersDict[i])
                    {
                        for (uint j = 0; j < cols; j++)
                        {
                            newCentroids[i][j] += matrix[line][j] / itemsInCluster;
                        }
                    }
                    // Сравниваем со старым центроидом и копируем в него, если есть различие
                    for (uint j = 0; j < cols; j++)
                    {
                        if (centroids[i][j] != newCentroids[i][j])
                        {
                            centroids[i][j] = newCentroids[i][j];
                            isCont = true;
                        }
                    }
                }
            } while (isCont);
            return new ClusterResult(clustersDict, averageDeviation, maxDeviation);
        }
        /// <summary>
        /// Возвращает начальные значения для центроидов
        /// </summary>
        /// <param name="matrix">Матрица</param>
        /// <param name="numCluster">Количество кластеров</param>
        /// <returns>Двухмерный массив центроидов</returns>
        private static double[][] InitialCentroids(double[][] matrix, uint numCluster)
        {
            uint rows = (uint)matrix.Length,
                 cols = (uint)matrix[0].Length;
            double[][] centroids = new double[numCluster][];
            for (uint i = 0; i < numCluster; i++)
            {
                centroids[i] = new double[cols];
                uint line = i * rows / numCluster;
                for (uint j = 0; j < cols; j++)
                {
                    centroids[i][j] = matrix[line][j];
                }
            }
            return centroids;
        }
        /// <summary>
        /// Вычисляет среднеквадратичное отклонение между двумя векторами
        /// </summary>
        /// <param name="array1">1й вектор</param>
        /// <param name="array2">2й вектор</param>
        /// <returns>Отклонение</returns>
        private static double Range(double[] array1, double[] array2)
        {
            uint len;
            if ((len = (uint)array1.Length) != array2.Length)
                throw new ClusterAnalyzeException("double Range(double[] array1, double[] array2):\n" +
                    "Массивы имеют разную длину");
            double sum = 0;
            for (uint i = 0; i < len; i++) sum += Math.Pow(array1[i] - array2[i], 2);
            return Math.Sqrt(sum / len);
        }
        /// <summary>
        /// Преобразует массив вида [,] к [][]
        /// </summary>
        /// <param name="matrix"></param>
        /// <returns></returns>
        private static double[][] Matrix2JadMatrix(double[,] matrix)
        {
            uint rows = (uint)matrix.GetLength(0),
                 cols = (uint)matrix.GetLength(1);
            double[][] jadMatrix = new double[rows][];
            for (uint i = 0; i < rows; i++)
            {
                jadMatrix[i] = new double[cols];
                for (uint j = 0; j < cols; j++)
                {
                    jadMatrix[i][j] = matrix[i, j];
                }
            }
            return jadMatrix;
        }
    }
    /// <summary>
    /// Результат анализа возвращается в виде этого класса
    /// </summary>
    public class ClusterResult : IEnumerable
    {
        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="dict"></param>
        /// <param name="averageDeviation">Среднее отклонение</param>
        /// <param name="maxDeviation">Максимальное отклонение</param>
        public ClusterResult(Dictionary<uint, ListOfRecords> dict, double averageDeviation, double maxDeviation)
        {
            AverageDeviation = averageDeviation;
            MaxDeviation = maxDeviation;
            data = dict.Values.ToArray();
        }
        /// <summary>
        /// Список индексов записей в кластере
        /// </summary>
        /// <param name="Index">Номер кластера</param>
        /// <returns>Список индексов записей</returns>
        public ListOfRecords this[uint Index]
        {
            get
            {
                return data[Index];
            }
        }
        /// <summary>
        /// Количество кластеров
        /// </summary>
        public uint Size
        {
            get
            {
                return (uint)data.Length;
            }
        }
        /// <summary>
        /// Возвращает тип перечислителя
        /// </summary>
        /// <returns>Перечислитель</returns>
        public IEnumerator GetEnumerator()
        {
            return data.GetEnumerator();
        }
        /// <summary>
        /// Средне отклонение записей от центроидов кластеров
        /// </summary>
        public readonly double AverageDeviation;
        /// <summary>
        /// Максимальное отклонение от центроидов кластеров
        /// </summary>
        public readonly double MaxDeviation;
        private readonly ListOfRecords[] data;
    }
    /// <summary>
    /// Класс-исключение для перехвата ошибок возникающих в классе ClusterAnalyze
    /// </summary>
    public class ClusterAnalyzeException : ApplicationException
    {
        /// <summary>
        /// Конструктор
        /// </summary>
        /// <param name="msg">Сообщение об ошибке</param>
        public ClusterAnalyzeException(string msg)
            : base(msg)
        { }
    }
    /// <summary>
    /// Класс список номеров записей в файле. Служит исключительно как синоним базового класса
    /// </summary>
    public class ListOfRecords : List<uint>
    { }
}