using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading;

namespace SMAT_CE.Matlab
{
    /// <summary>
    /// 矩阵运算
    /// </summary>
    public static class MatrixCompute
    {
        #region inv 高斯-约当法求逆矩阵
        #region Float
        /// <summary>
        /// 高斯-约当法矩阵求逆
        /// </summary>
        /// <remarks>
        /// 实矩阵求逆的全选主元高斯-约当法
        /// 算法原理：
        ///     对于k从0到n-1做如下操作：
        ///     （1）从第k行，第k列开始的右下角子阵选取绝对值最大的元素，并标记此元素所在的行号和列号，再通过行交换与列交换将它交换到主元素位置上（这一步称为全选主元）
        ///     （2）1/akk => akk
        ///     （3）akjakk => akj, j=0,1,...,n-1;j!=k
        ///     （4）aij - aikakj => aij, i,j=0,1,...,n-1;i,j!=k
        ///     （5）-aikakk => aik, i=0,1,...,n-1; i!=k
        ///     根据在全选主元过程户所记录的行、列交换信息进行恢复，恢复的原则如下：
        ///     在全选主元过程中，先交换的行、列后进行恢复；原来的行(列)交换用列(行)交换来恢复。
        /// </remarks>
        /// <param name="matrix">矩阵</param>
        public static bool Inv(Matrix<float> matrix, out Matrix<float> result)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new Exception("Matrix must be square.");
            }
            result = matrix.Clone();
            int i, j, k;
            int cols = result.ColumnCount;
            float d, p;
            // 分配内存
            int[] pnRow = new int[cols];
            int[] pnCol = new int[cols];
            // 消元
            for (k = 0; k < cols; k++)
            {
                d = 0.0f;
                for (i = k; i < cols; i++)
                {
                    for (j = k; j < cols; j++)
                    {
                        p = Math.Abs(result.At(i, j));
                        if (p > d)
                        {
                            d = p;
                            pnRow[k] = i;
                            pnCol[k] = j;
                        }
                    }
                }
                // 失败
                if (d == 0.0f || d + 1.0f == 1.0f)
                {
                    pnRow = new int[0];
                    pnCol = new int[0];
                    return false;
                }

                if (pnRow[k] != k)
                {
                    for (j = 0; j < cols; j++)
                    {
                        p = result.At(k, j);
                        result.At(k, j, result.At(pnRow[k], j));
                        result.At(pnRow[k], j, p);
                    }
                }

                if (pnCol[k] != k)
                {
                    for (i = 0; i < cols; i++)
                    {
                        p = result.At(i, k);
                        result.At(i, k, result.At(i, pnCol[k]));
                        result.At(i, pnCol[k], p);
                    }
                }

                result.At(k, k, 1.0f / result.At(k, k));
                for (j = 0; j < cols; j++)
                {
                    if (j != k)
                    {
                        result.At(k, j, result.At(k, j) * result.At(k, k));
                    }
                }

                for (i = 0; i < cols; i++)
                {
                    if (i != k)
                    {
                        for (j = 0; j < cols; j++)
                        {
                            if (j != k)
                            {
                                result.At(i, j, result.At(i, j) - result.At(i, k) * result.At(k, j));
                            }
                        }
                    }
                }

                for (i = 0; i < cols; i++)
                {
                    if (i != k)
                    {
                        result.At(i, k, -result.At(i, k) * result.At(k, k));
                    }
                }
            }

            // 调整恢复行列次序
            for (k = cols - 1; k >= 0; k--)
            {
                if (pnCol[k] != k)
                {
                    for (j = 0; j < cols; j++)
                    {
                        p = result.At(k, j);
                        result.At(k, j, result.At(pnCol[k], j));
                        result.At(pnCol[k], j, p);
                    }
                }

                if (pnRow[k] != k)
                {
                    for (i = 0; i < cols; i++)
                    {
                        p = result.At(i, k);
                        result.At(i, k, result.At(i, pnRow[k]));
                        result.At(i, pnRow[k], p);
                    }
                }
            }

            // 成功返回
            pnRow = new int[0];
            pnCol = new int[0];
            return true;
        }
        #endregion

        #region Double
        /// <summary>
        /// 高斯-约当法矩阵求逆
        /// </summary>
        /// <remarks>
        /// 实矩阵求逆的全选主元高斯-约当法
        /// 算法原理：
        ///     对于k从0到n-1做如下操作：
        ///     （1）从第k行，第k列开始的右下角子阵选取绝对值最大的元素，并标记此元素所在的行号和列号，再通过行交换与列交换将它交换到主元素位置上（这一步称为全选主元）
        ///     （2）1/akk => akk
        ///     （3）akjakk => akj, j=0,1,...,n-1;j!=k
        ///     （4）aij - aikakj => aij, i,j=0,1,...,n-1;i,j!=k
        ///     （5）-aikakk => aik, i=0,1,...,n-1; i!=k
        ///     根据在全选主元过程户所记录的行、列交换信息进行恢复，恢复的原则如下：
        ///     在全选主元过程中，先交换的行、列后进行恢复；原来的行(列)交换用列(行)交换来恢复。
        /// </remarks>
        /// <param name="matrix">矩阵</param>
        public static bool Inv(Matrix<double> matrix, out Matrix<double> result)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new Exception("Matrix must be square.");
            }
            result = matrix.Clone();
            int i, j, k;
            int cols = result.ColumnCount;
            double d, p;
            // 分配内存
            int[] pnRow = new int[cols];
            int[] pnCol = new int[cols];
            // 消元
            for (k = 0; k < cols; k++)
            {
                d = 0.0;
                for (i = k; i < cols; i++)
                {
                    for (j = k; j < cols; j++)
                    {
                        p = Math.Abs(result.At(i, j));
                        if (p > d)
                        {
                            d = p;
                            pnRow[k] = i;
                            pnCol[k] = j;
                        }
                    }
                }
                // 失败
                if (d == 0.0 || d + 1.0 == 1.0)
                {
                    pnRow = new int[0];
                    pnCol = new int[0];
                    return false;
                }

                if (pnRow[k] != k)
                {
                    for (j = 0; j < cols; j++)
                    {
                        p = result.At(k, j);
                        result.At(k, j, result.At(pnRow[k], j));
                        result.At(pnRow[k], j, p);
                    }
                }

                if (pnCol[k] != k)
                {
                    for (i = 0; i < cols; i++)
                    {
                        p = result.At(i, k);
                        result.At(i, k, result.At(i, pnCol[k]));
                        result.At(i, pnCol[k], p);
                    }
                }

                result.At(k, k, 1.0 / result.At(k, k));
                for (j = 0; j < cols; j++)
                {
                    if (j != k)
                    {
                        result.At(k, j, result.At(k, j) * result.At(k, k));
                    }
                }

                for (i = 0; i < cols; i++)
                {
                    if (i != k)
                    {
                        for (j = 0; j < cols; j++)
                        {
                            if (j != k)
                            {
                                result.At(i, j, result.At(i, j) - result.At(i, k) * result.At(k, j));
                            }
                        }
                    }
                }

                for (i = 0; i < cols; i++)
                {
                    if (i != k)
                    {
                        result.At(i, k, -result.At(i, k) * result.At(k, k));
                    }
                }
            }

            // 调整恢复行列次序
            for (k = cols - 1; k >= 0; k--)
            {
                if (pnCol[k] != k)
                {
                    for (j = 0; j < cols; j++)
                    {
                        p = result.At(k, j);
                        result.At(k, j, result.At(pnCol[k], j));
                        result.At(pnCol[k], j, p);
                    }
                }

                if (pnRow[k] != k)
                {
                    for (i = 0; i < cols; i++)
                    {
                        p = result.At(i, k);
                        result.At(i, k, result.At(i, pnRow[k]));
                        result.At(i, pnRow[k], p);
                    }
                }
            }

            // 成功返回
            pnRow = new int[0];
            pnCol = new int[0];
            return true;
        }
        #endregion
        #endregion

        #region ismember 检查元素是否存在另一元素中
        #region Float
        /// <summary>
        /// 矩阵各个元素是否出现在常数中
        /// 返回矩阵按列排列的数组，索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="scalar">常数</param>
        public static Tuple<bool[], int[]> MatrixElementIsmemberScalar(Matrix<float> matrix, float scalar)
        {
            float[] data = matrix.ToColumnWiseArray();
            int len = data.Length;
            bool[] boolResult = new bool[len];
            int[] intResult = new int[len];
            int index = Array.IndexOf(data, scalar);
            if (index >= 0)
            {
                boolResult.SetValue(true, index);
                intResult.SetValue(1, index);
            }
            data = new float[0];
            return Tuple.Create<bool[], int[]>(boolResult, intResult);
        }
        /// <summary>
        /// 常数是否出现在矩阵各个元素中
        /// 返回的索引从1开始
        /// </summary>
        /// <param name="scalar">常数</param>
        /// <param name="matrix">矩阵</param>
        public static Tuple<bool, int> ScalarIsmemberMatrixElement(float scalar, Matrix<float> matrix)
        {
            float[] data = matrix.ToColumnWiseArray();
            int index = Array.IndexOf(data, scalar);
            data = new float[0];
            if (index >= 0)
            {
                return Tuple.Create<bool, int>(true, index + 1);
            }
            return Tuple.Create<bool, int>(false, 0);
        }
        /// <summary>
        /// 矩阵各个元素是否出现在另一矩阵中
        /// 返回矩阵按列排列的数组，索引从1开始
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        public static Tuple<bool[], int[]> MatrixElementIsmemberMatrixElement(Matrix<float> matrix1, Matrix<float> matrix2)
        {
            float[] data1 = matrix1.ToColumnWiseArray();
            float[] data2 = matrix2.ToColumnWiseArray();
            int len = data1.Length;
            bool[] boolResult = new bool[len];
            int[] intResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        int index = Array.IndexOf(data2, data1[i]);
            //        if (index >= 0)
            //        {
            //            boolResult.SetValue(true, i);
            //            intResult.SetValue(index + 1, i);
            //        }
            //    }
            //});
            for (int i = 0; i < len; i++)
            {
                int index = Array.IndexOf(data2, data1[i]);
                if (index >= 0)
                {
                    boolResult.SetValue(true, i);
                    intResult.SetValue(index + 1, i);
                }
            }
            data1 = new float[0]; data2 = new float[0];
            return Tuple.Create<bool[], int[]>(boolResult, intResult);
        }
        /// <summary>
        /// 矩阵各行是否出现在另一矩阵中
        /// 返回的行索引从1开始
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        public static Tuple<bool[], int[]> MatrixRowIsmemberMatrixRow(Matrix<float> matrix1, Matrix<float> matrix2)
        {
            if (matrix1.ColumnCount != matrix2.ColumnCount)
            {
                throw new Exception("Inputs matrix1 and matrix2 must be matrices with the same number of columns in the 'rows' case.");
            }
            int len = matrix1.RowCount;
            bool[] boolResult = new bool[len];
            int[] intResult = new int[len];
            List<SingleVector> data2 = new List<SingleVector>();
            float[] row;
            SingleVector vector;
            for (int i = 0; i < matrix2.RowCount; i++)
            {
                row = matrix2.Row(i).ToArray();
                vector = new SingleVector(row);
                data2.Add(vector);
            }
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        row = matrix1.Row(i).ToArray();
            //        vector = new SingleVector(row);
            //        int index = data2.IndexOf(vector);
            //        if (index >= 0)
            //        {
            //            boolResult.SetValue(true, i);
            //            intResult.SetValue(index + 1, i);
            //        }
            //    }
            //});
            for (int i = 0; i < len; i++)
            {
                row = matrix1.Row(i).ToArray();
                vector = new SingleVector(row);
                int index = data2.IndexOf(vector);
                if (index >= 0)
                {
                    boolResult.SetValue(true, i);
                    intResult.SetValue(index + 1, i);
                }
            }
            data2.Clear(); row = new float[0];
            return Tuple.Create<bool[], int[]>(boolResult, intResult);
        }
        #endregion

        #region Double
        /// <summary>
        /// 矩阵各个元素是否出现在常数中
        /// 返回矩阵按列排列的数组，索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="scalar">常数</param>
        public static Tuple<bool[], int[]> MatrixElementIsmemberScalar(Matrix<double> matrix, double scalar)
        {
            double[] data = matrix.ToColumnWiseArray();
            int len = data.Length;
            bool[] boolResult = new bool[len];
            int[] intResult = new int[len];
            int index = Array.IndexOf(data, scalar);
            if (index >= 0)
            {
                boolResult.SetValue(true, index);
                intResult.SetValue(1, index);
            }
            data = new double[0];
            return Tuple.Create<bool[], int[]>(boolResult, intResult);
        }
        /// <summary>
        /// 常数是否出现在矩阵各个元素中
        /// 返回的索引从1开始
        /// </summary>
        /// <param name="scalar">常数</param>
        /// <param name="matrix">矩阵</param>
        public static Tuple<bool, int> ScalarIsmemberMatrixElement(double scalar, Matrix<double> matrix)
        {
            double[] data = matrix.ToColumnWiseArray();
            int index = Array.IndexOf(data, scalar);
            data = new double[0];
            if (index >= 0)
            {
                return Tuple.Create<bool, int>(true, index + 1);
            }
            return Tuple.Create<bool, int>(false, 0);
        }
        /// <summary>
        /// 矩阵各个元素是否出现在另一矩阵中
        /// 返回矩阵按列排列的数组，索引从1开始
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        public static Tuple<bool[], int[]> MatrixElementIsmemberMatrixElement(Matrix<double> matrix1, Matrix<double> matrix2)
        {
            double[] data1 = matrix1.ToColumnWiseArray();
            double[] data2 = matrix2.ToColumnWiseArray();
            int len = data1.Length;
            bool[] boolResult = new bool[len];
            int[] intResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        int index = Array.IndexOf(data2, data1[i]);
            //        if (index >= 0)
            //        {
            //            boolResult.SetValue(true, i);
            //            intResult.SetValue(index + 1, i);
            //        }
            //    }
            //});
            for (int i = 0; i < len; i++)
            {
                int index = Array.IndexOf(data2, data1[i]);
                if (index >= 0)
                {
                    boolResult.SetValue(true, i);
                    intResult.SetValue(index + 1, i);
                }
            }
            data1 = new double[0]; data2 = new double[0];
            return Tuple.Create<bool[], int[]>(boolResult, intResult);
        }
        /// <summary>
        /// 矩阵各行是否出现在另一矩阵中
        /// 返回的行索引从1开始
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        public static Tuple<bool[], int[]> MatrixRowIsmemberMatrixRow(Matrix<double> matrix1, Matrix<double> matrix2)
        {
            if (matrix1.ColumnCount != matrix2.ColumnCount)
            {
                throw new Exception("Inputs matrix1 and matrix2 must be matrices with the same number of columns in the 'rows' case.");
            }
            int len = matrix1.RowCount;
            bool[] boolResult = new bool[len];
            int[] intResult = new int[len];
            List<DoubleVector> data2 = new List<DoubleVector>();
            double[] row;
            DoubleVector vector;
            for (int i = 0; i < matrix2.RowCount; i++)
            {
                row = matrix2.Row(i).ToArray();
                vector = new DoubleVector(row);
                data2.Add(vector);
            }
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        row = matrix1.Row(i).ToArray();
            //        vector = new DoubleVector(row);
            //        int index = data2.IndexOf(vector);
            //        if (index >= 0)
            //        {
            //            boolResult.SetValue(true, i);
            //            intResult.SetValue(index + 1, i);
            //        }
            //    }
            //});
            for (int i = 0; i < len; i++)
            {
                row = matrix1.Row(i).ToArray();
                vector = new DoubleVector(row);
                int index = data2.IndexOf(vector);
                if (index >= 0)
                {
                    boolResult.SetValue(true, i);
                    intResult.SetValue(index + 1, i);
                }
            }
            data2.Clear(); row = new double[0];
            return Tuple.Create<bool[], int[]>(boolResult, intResult);
        }
        #endregion
        #endregion

        #region max 行或列最大值
        #region Float
        /// <summary>
        /// 行方向的最大值和索引
        /// 索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Tuple<float[], int[]> MaxRowWise(Matrix<float> matrix)
        {
            int len = matrix.RowCount;
            float[] valueResult = new float[len];
            int[] indexResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        float[] data = matrix.Row(i).ToArray();
            //        float max = data.Max();
            //        int index = Array.IndexOf(data, max);
            //        valueResult[i] = max;
            //        indexResult[i] = index + 1;
            //    }
            //});
            float[] data;
            for (int i = 0; i < len; i++)
            {
                data = matrix.Row(i).ToArray();
                float max = data.Max();
                int index = Array.IndexOf(data, max);
                valueResult[i] = max;
                indexResult[i] = index + 1;
            }
            data = new float[0];
            return Tuple.Create<float[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 列方向的最大值和索引
        /// 索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Tuple<float[], int[]> MaxColumnWise(Matrix<float> matrix)
        {
            int len = matrix.ColumnCount;
            float[] valueResult = new float[len];
            int[] indexResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        float[] data = matrix.Column(i).ToArray();
            //        float max = data.Max();
            //        int index = Array.IndexOf(data, max);
            //        valueResult[i] = max;
            //        indexResult[i] = index + 1;
            //    }
            //});
            float[] data;
            for (int i = 0; i < len; i++)
            {
                data = matrix.Column(i).ToArray();
                float max = data.Max();
                int index = Array.IndexOf(data, max);
                valueResult[i] = max;
                indexResult[i] = index + 1;
            }
            data = new float[0];
            return Tuple.Create<float[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 与常量比较返回最大值矩阵
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="scalar">常量</param>
        public static Matrix<float> MaxCompareToScalar(Matrix<float> matrix, float scalar)
        {
            float[] data = matrix.ToColumnWiseArray();
            //CommonParallel.For(0, data.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        data[i] = Math.Max(data[i], scalar);
            //    }
            //});
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = Math.Max(data[i], scalar);
            }
            return Matrix<float>.Build.DenseOfColumnMajor(matrix.RowCount, matrix.ColumnCount, data);
        }
        /// <summary>
        /// 与矩阵比较返回最大值矩阵
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        public static Matrix<float> MaxCompareToMatrix(Matrix<float> matrix1, Matrix<float> matrix2)
        {
            if (matrix1.RowCount != matrix2.RowCount || matrix1.ColumnCount != matrix2.ColumnCount)
            {
                throw new Exception("Matrix dimensions must agree.");
            }
            float[] data1 = matrix1.ToColumnWiseArray();
            float[] data2 = matrix2.ToColumnWiseArray();
            //CommonParallel.For(0, data1.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        data1[i] = Math.Max(data1[i], data2[i]);
            //    }
            //});
            for (int i = 0; i < data1.Length; i++)
            {
                data1[i] = Math.Max(data1[i], data2[i]);
            }
            data2 = new float[0];
            return Matrix<float>.Build.DenseOfColumnMajor(matrix1.RowCount, matrix1.ColumnCount, data1);
        }
        #endregion

        #region Double
        /// <summary>
        /// 行方向的最大值和索引
        /// 索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Tuple<double[], int[]> MaxRowWise(Matrix<double> matrix)
        {
            int len = matrix.RowCount;
            double[] valueResult = new double[len];
            int[] indexResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        double[] data = matrix.Row(i).ToArray();
            //        double max = data.Max();
            //        int index = Array.IndexOf(data, max);
            //        valueResult[i] = max;
            //        indexResult[i] = index + 1;
            //    }
            //});
            double[] data;
            for (int i = 0; i < len; i++)
            {
                data = matrix.Row(i).ToArray();
                double max = data.Max();
                int index = Array.IndexOf(data, max);
                valueResult[i] = max;
                indexResult[i] = index + 1;
            }
            data = new double[0];
            return Tuple.Create<double[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 列方向的最大值和索引
        /// 索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Tuple<double[], int[]> MaxColumnWise(Matrix<double> matrix)
        {
            int len = matrix.ColumnCount;
            double[] valueResult = new double[len];
            int[] indexResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        double[] data = matrix.Column(i).ToArray();
            //        double max = data.Max();
            //        int index = Array.IndexOf(data, max);
            //        valueResult[i] = max;
            //        indexResult[i] = index + 1;
            //    }
            //});
            double[] data;
            for (int i = 0; i < len; i++)
            {
                data = matrix.Column(i).ToArray();
                double max = data.Max();
                int index = Array.IndexOf(data, max);
                valueResult[i] = max;
                indexResult[i] = index + 1;
            }
            data = new double[0];
            return Tuple.Create<double[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 与常量比较返回最大值矩阵
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="scalar">常量</param>
        public static Matrix<double> MaxCompareToScalar(Matrix<double> matrix, double scalar)
        {
            double[] data = matrix.ToColumnWiseArray();
            //CommonParallel.For(0, data.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        data[i] = Math.Max(data[i], scalar);
            //    }
            //});
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = Math.Max(data[i], scalar);
            }
            return Matrix<double>.Build.DenseOfColumnMajor(matrix.RowCount, matrix.ColumnCount, data);
        }
        /// <summary>
        /// 与矩阵比较返回最大值矩阵
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        public static Matrix<double> MaxCompareToMatrix(Matrix<double> matrix1, Matrix<double> matrix2)
        {
            if (matrix1.RowCount != matrix2.RowCount || matrix1.ColumnCount != matrix2.ColumnCount)
            {
                throw new Exception("Matrix dimensions must agree.");
            }
            double[] data1 = matrix1.ToColumnWiseArray();
            double[] data2 = matrix2.ToColumnWiseArray();
            //CommonParallel.For(0, data1.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        data1[i] = Math.Max(data1[i], data2[i]);
            //    }
            //});
            for (int i = 0; i < data1.Length; i++)
            {
                data1[i] = Math.Max(data1[i], data2[i]);
            }
            data2 = new double[0];
            return Matrix<double>.Build.DenseOfColumnMajor(matrix1.RowCount, matrix1.ColumnCount, data1);
        }
        #endregion
        #endregion

        #region min 行或列最小值
        #region Float
        /// <summary>
        /// 行方向的最小值和索引
        /// 索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Tuple<float[], int[]> MinRowWise(Matrix<float> matrix)
        {
            int len = matrix.RowCount;
            float[] valueResult = new float[len];
            int[] indexResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        float[] data = matrix.Row(i).ToArray();
            //        float min = data.Min();
            //        int index = Array.IndexOf(data, min);
            //        valueResult[i] = min;
            //        indexResult[i] = index + 1;
            //    }
            //});
            float[] data;
            for (int i = 0; i < len; i++)
            {
                data = matrix.Row(i).ToArray();
                float min = data.Min();
                int index = Array.IndexOf(data, min);
                valueResult[i] = min;
                indexResult[i] = index + 1;
            }
            data = new float[0];
            return Tuple.Create<float[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 列方向的最小值和索引
        /// 索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Tuple<float[], int[]> MinColumnWise(Matrix<float> matrix)
        {
            int len = matrix.ColumnCount;
            float[] valueResult = new float[len];
            int[] indexResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        float[] data = matrix.Column(i).ToArray();
            //        float min = data.Min();
            //        int index = Array.IndexOf(data, min);
            //        valueResult[i] = min;
            //        indexResult[i] = index + 1;
            //    }
            //});
            float[] data;
            for (int i = 0; i < len; i++)
            {
                data = matrix.Column(i).ToArray();
                float min = data.Min();
                int index = Array.IndexOf(data, min);
                valueResult[i] = min;
                indexResult[i] = index + 1;
            }
            data = new float[0];
            return Tuple.Create<float[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 与常量比较返回最小值矩阵
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="scalar">常量</param>
        public static Matrix<float> MinCompareToScalar(Matrix<float> matrix, float scalar)
        {
            float[] data = matrix.ToColumnWiseArray();
            //CommonParallel.For(0, data.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        data[i] = Math.Min(data[i], scalar);
            //    }
            //});
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = Math.Min(data[i], scalar);
            }
            return Matrix<float>.Build.DenseOfColumnMajor(matrix.RowCount, matrix.ColumnCount, data);
        }
        /// <summary>
        /// 与矩阵比较返回最小值矩阵
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        public static Matrix<float> MinCompareToMatrix(Matrix<float> matrix1, Matrix<float> matrix2)
        {
            if (matrix1.RowCount != matrix2.RowCount || matrix1.ColumnCount != matrix2.ColumnCount)
            {
                throw new Exception("Matrix dimensions must agree.");
            }
            float[] data1 = matrix1.ToColumnWiseArray();
            float[] data2 = matrix2.ToColumnWiseArray();
            //CommonParallel.For(0, data1.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        data1[i] = Math.Min(data1[i], data2[i]);
            //    }
            //});
            for (int i = 0; i < data1.Length; i++)
            {
                data1[i] = Math.Min(data1[i], data2[i]);
            }
            data2 = new float[0];
            return Matrix<float>.Build.DenseOfColumnMajor(matrix1.RowCount, matrix1.ColumnCount, data1);
        }
        #endregion

        #region Double
        /// <summary>
        /// 行方向的最小值和索引
        /// 索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Tuple<double[], int[]> MinRowWise(Matrix<double> matrix)
        {
            int len = matrix.RowCount;
            double[] valueResult = new double[len];
            int[] indexResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        double[] data = matrix.Row(i).ToArray();
            //        double min = data.Min();
            //        int index = Array.IndexOf(data, min);
            //        valueResult[i] = min;
            //        indexResult[i] = index + 1;
            //    }
            //});
            double[] data;
            for (int i = 0; i < len; i++)
            {
                data = matrix.Row(i).ToArray();
                double min = data.Min();
                int index = Array.IndexOf(data, min);
                valueResult[i] = min;
                indexResult[i] = index + 1;
            }
            data = new double[0];
            return Tuple.Create<double[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 列方向的最小值和索引
        /// 索引从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Tuple<double[], int[]> MinColumnWise(Matrix<double> matrix)
        {
            int len = matrix.ColumnCount;
            double[] valueResult = new double[len];
            int[] indexResult = new int[len];
            //CommonParallel.For(0, len, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        double[] data = matrix.Column(i).ToArray();
            //        double min = data.Min();
            //        int index = Array.IndexOf(data, min);
            //        valueResult[i] = min;
            //        indexResult[i] = index + 1;
            //    }
            //});
            double[] data;
            for (int i = 0; i < len; i++)
            {
                data = matrix.Column(i).ToArray();
                double min = data.Min();
                int index = Array.IndexOf(data, min);
                valueResult[i] = min;
                indexResult[i] = index + 1;
            }
            data = new double[0];
            return Tuple.Create<double[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 与常量比较返回最小值矩阵
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="scalar">常量</param>
        public static Matrix<double> MinCompareToScalar(Matrix<double> matrix, double scalar)
        {
            double[] data = matrix.ToColumnWiseArray();
            //CommonParallel.For(0, data.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        data[i] = Math.Min(data[i], scalar);
            //    }
            //});
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = Math.Min(data[i], scalar);
            }
            return Matrix<double>.Build.DenseOfColumnMajor(matrix.RowCount, matrix.ColumnCount, data);
        }
        /// <summary>
        /// 与矩阵比较返回最小值矩阵
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        public static Matrix<double> MinCompareToMatrix(Matrix<double> matrix1, Matrix<double> matrix2)
        {
            if (matrix1.RowCount != matrix2.RowCount || matrix1.ColumnCount != matrix2.ColumnCount)
            {
                throw new Exception("Matrix dimensions must agree.");
            }
            double[] data1 = matrix1.ToColumnWiseArray();
            double[] data2 = matrix2.ToColumnWiseArray();
            //CommonParallel.For(0, data1.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        data1[i] = Math.Min(data1[i], data2[i]);
            //    }
            //});
            for (int i = 0; i < data1.Length; i++)
            {
                data1[i] = Math.Min(data1[i], data2[i]);
            }
            data2 = new double[0];
            return Matrix<double>.Build.DenseOfColumnMajor(matrix1.RowCount, matrix1.ColumnCount, data1);
        }
        #endregion
        #endregion

        #region pdist2 两个矩阵的欧几里得距离
        #region Float
        /// <summary>
        /// 计算两个矩阵的欧几里得距离
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        /// <param name="threadCount">线程数</param>
        public static Matrix<float> Pdist2(Matrix<float> matrix1, Matrix<float> matrix2, int threadCount = 0)
        {
            if (matrix1.ColumnCount != matrix2.ColumnCount)
            {
                throw new Exception("matrix1 and matrix2 must have the same number of columns.");
            }
            #region 一般方式
            //Matrix<float> result = Matrix<float>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //for (int i = 0; i < result.RowCount; i++)
            //{
            //    float[] left = matrix1.Row(i).ToArray();
            //    for (int j = 0; j < result.ColumnCount; j++)
            //    {
            //        float[] right = matrix2.Row(j).ToArray();
            //        float d = MathNet.Numerics.Distance.Euclidean(left, right);
            //        result.At(i, j, d);
            //    }
            //}
            //return result;

            int len1 = matrix1.RowCount, len2 = matrix2.RowCount;
            float[] data = new float[len1 * len2], datarow, left, right;
            for (int j = 0; j < len2; j++)
            {
                datarow = new float[len1];
                right = matrix2.Row(j).ToArray();
                for (int i = 0; i < len1; i++)
                {
                    left = matrix1.Row(i).ToArray();
                    float d = MathNet.Numerics.Distance.Euclidean(left, right);
                    datarow[i] = d;
                }
                Array.Copy(datarow, 0, data, j * len1, len1);
            }
            datarow = new float[0]; left = new float[0]; right = new float[0];
            return Matrix<float>.Build.DenseOfColumnMajor(len1, len2, data);
            #endregion
            #region 并行方式
            //ConcurrentDictionary<int, float[]> dic = new ConcurrentDictionary<int, float[]>();
            //CommonParallel.For(0, matrix1.RowCount, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        float[] v = new double[matrix2.RowCount];
            //        float[] left = matrix1.Row(i).ToArray();
            //        for (int j = 0; j < matrix2.RowCount; j++)
            //        {
            //            float[] right = matrix2.Row(j).ToArray();
            //            float d = MathNet.Numerics.Distance.Euclidean(left, right);
            //            v[j] = d;
            //        }
            //        dic.TryAdd(i, v);
            //    }
            //});
            //Matrix<double> result = Matrix<double>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //// 要不要排序？
            //for (int i = 0; i < dic.Count; i++)
            //{
            //    float[] row;
            //    dic.TryGetValue(i, out row);
            //    result.SetRow(i, row);
            //}
            //return result;
            #endregion
            #region 多线程方式，线程返回值
            //int threadCount = MathNet.Numerics.Control.MaxDegreeOfParallelism;
            //int row_base = matrix1.RowCount / threadCount;
            //int row_remain = matrix1.RowCount % threadCount;
            //List<Tuple<int, int>> indexList = new List<Tuple<int, int>>();
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}
            //Task<Matrix<float>>[] tasks = new Task<Matrix<float>>[threadCount];
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    int to = indexList[i].Item2;
            //    tasks[i] = Task<Matrix<float>>.Factory.StartNew(() =>
            //    {
            //        Matrix<float> m = Matrix<float>.Build.Dense(to - from, matrix2.RowCount);
            //        for (int p = 0; p < m.RowCount; p++)
            //        {
            //            float[] left = matrix1.Row(p).ToArray();
            //            for (int q = 0; q < m.ColumnCount; q++)
            //            {
            //                float[] right = matrix2.Row(q).ToArray();
            //                float d = MathNet.Numerics.Distance.Euclidean(left, right);
            //                m.At(p, q, d);
            //            }
            //        }
            //        return m;
            //    });
            //}
            //Task.WaitAll(tasks);
            //Matrix<float> result = Matrix<float>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    var submatrix = tasks[i].Result;
            //    result.SetSubMatrix(from, 0, submatrix);
            //}
            //return result;

            //threadCount = (threadCount > MathNet.Numerics.Control.MaxDegreeOfParallelism || threadCount <= 0) ?
            //    MathNet.Numerics.Control.MaxDegreeOfParallelism :
            //    threadCount;
            //if (matrix1.RowCount <= threadCount)
            //    threadCount = 1;
            //int row_base = matrix1.RowCount / threadCount;
            //int row_remain = matrix1.RowCount % threadCount;

            //List<Tuple<int, int>> indexList = new List<Tuple<int, int>>();
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}

            //Task<Matrix<float>>[] tasks = new Task<Matrix<float>>[threadCount];
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    int to = indexList[i].Item2;
            //    tasks[i] = Task<Matrix<float>>.Factory.StartNew(() =>
            //    {
            //        int len1 = to - from, len2 = matrix2.RowCount;
            //        float[] resultrow = new float[len1 * len2];
            //        for (int j = 0; j < len2; j++)
            //        {
            //            int n = 0;
            //            float[] datarow = new float[len1];
            //            float[] right = matrix2.Row(j).ToArray();
            //            for (int k = from; k < to; k++)
            //            {
            //                float[] left = matrix1.Row(k).ToArray();
            //                float d = MathNet.Numerics.Distance.Euclidean(left, right);
            //                datarow[n] = d;
            //                n++;
            //            }
            //            Array.Copy(datarow, 0, resultrow, j * len1, len1);
            //        }
            //        return Matrix<float>.Build.DenseOfColumnMajor(len1, len2, resultrow);
            //    }, TaskCreationOptions.PreferFairness);
            //}
            //Task.WaitAll(tasks);

            //Matrix<float> result = Matrix<float>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    var submatrix = tasks[i].Result;
            //    result.SetSubMatrix(from, 0, submatrix);
            //}
            //return result;
            #endregion
            #region 多线程方式，任务工厂
            //// 确定线程数
            //threadCount = (threadCount > MathNet.Numerics.Control.MaxDegreeOfParallelism || threadCount <= 0) ?
            //    MathNet.Numerics.Control.MaxDegreeOfParallelism :
            //    threadCount;
            //if (matrix1.RowCount <= threadCount)
            //    threadCount = 1;
            //// 子任务分块
            //int row_base = matrix1.RowCount / threadCount;
            //int row_remain = matrix1.RowCount % threadCount;
            //List<Tuple<int, int>> indexList = new List<Tuple<int, int>>();
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}
            //// 创建一个父任务，由父任务再创建子任务，父任务返回每个子任务的计算结果
            //Task<Matrix<float>[]> parentTask = new Task<Matrix<float>[]>(() =>
            //{
            //    CancellationTokenSource cts = new CancellationTokenSource();
            //    // 任务工厂统一创建子任务
            //    TaskFactory<Matrix<float>> factory = new TaskFactory<Matrix<float>>(cts.Token, TaskCreationOptions.AttachedToParent, TaskContinuationOptions.PreferFairness, TaskScheduler.Default);
            //    Matrix<float>[] submatrixs = new Matrix<float>[threadCount];
            //    Task<Matrix<float>>[] childTasks = new Task<Matrix<float>>[threadCount];
            //    for (int i = 0; i < threadCount; i++)
            //    {
            //        int n = i;
            //        int from = indexList[n].Item1;
            //        int to = indexList[n].Item2;
            //        childTasks[n] = factory.StartNew(() => submatrixs[n] = GetDoubleSubMatrix(cts.Token, from, to, matrix1, matrix2));
            //        // 这里，任何子任务抛出异常就取消其余子任务
            //        childTasks[n].ContinueWith(t => cts.Cancel(), TaskContinuationOptions.OnlyOnFaulted);
            //    }
            //    return submatrixs;
            //});
            //// 拼装结果矩阵，使用ContinueWith，父任务完成时启动另一个任务去输出结果
            //// 这里，应该也可以用parentTask.Wait(); parentTask.Result; 得到结果后再输出
            //Matrix<float> result = Matrix<float>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //var wt = parentTask.ContinueWith(task =>
            //{
            //    for (int i = 0; i < threadCount; i++)
            //    {
            //        int from = indexList[i].Item1;
            //        var submatrix = task.Result[i];
            //        if (submatrix != null)
            //            result.SetSubMatrix(from, 0, submatrix);
            //    }
            //});
            //parentTask.Start();
            //wt.Wait();
            //return result;
            #endregion
        }
        /// <summary>
        /// 子任务事件，得到子矩阵
        /// </summary>
        [Obsolete("单线程计算欧氏距离用不到")]
        static Matrix<float> GetDoubleSubMatrix(CancellationToken ct, int from, int to, Matrix<float> matrix1, Matrix<float> matrix2)
        {
            try
            {
                int len1 = to - from, len2 = matrix2.RowCount;
                float[] resultrow = new float[len1 * len2];
                for (int j = 0; j < len2; j++)
                {
                    int n = 0;
                    float[] datarow = new float[len1];
                    float[] right = matrix2.Row(j).ToArray();
                    for (int k = from; k < to; k++)
                    {
                        float[] left = matrix1.Row(k).ToArray();
                        float d = MathNet.Numerics.Distance.Euclidean(left, right);
                        datarow[n] = d;
                        n++;
                    }
                    Array.Copy(datarow, 0, resultrow, j * len1, len1);
                }
                return Matrix<float>.Build.DenseOfColumnMajor(len1, len2, resultrow);
            }
            catch
            {
                ct.ThrowIfCancellationRequested();
            }
            return null;
        }
        #endregion

        #region Double
        /// <summary>
        /// 计算两个矩阵的欧几里得距离
        /// </summary>
        /// <param name="matrix1">矩阵</param>
        /// <param name="matrix2">矩阵</param>
        /// <param name="threadCount">线程数</param>
        public static Matrix<double> Pdist2(Matrix<double> matrix1, Matrix<double> matrix2, int threadCount = 0)
        {
            if (matrix1.ColumnCount != matrix2.ColumnCount)
            {
                throw new Exception("matrix1 and matrix2 must have the same number of columns.");
            }
            #region 一般方式
            //Matrix<double> result = Matrix<double>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //for (int i = 0; i < result.RowCount; i++)
            //{
            //    double[] left = matrix1.Row(i).ToArray();
            //    for (int j = 0; j < result.ColumnCount; j++)
            //    {
            //        double[] right = matrix2.Row(j).ToArray();
            //        double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //        result.At(i, j, d);
            //    }
            //}
            //return result;

            int len1 = matrix1.RowCount, len2 = matrix2.RowCount;
            double[] data = new double[len1 * len2], datarow, left, right;
            for (int j = 0; j < len2; j++)
            {
                datarow = new double[len1];
                right = matrix2.Row(j).ToArray();
                for (int i = 0; i < len1; i++)
                {
                    left = matrix1.Row(i).ToArray();
                    double d = MathNet.Numerics.Distance.Euclidean(left, right);
                    datarow[i] = d;
                }
                Array.Copy(datarow, 0, data, j * len1, len1);
            }
            datarow = new double[0]; left = new double[0]; right = new double[0];
            return Matrix<double>.Build.DenseOfColumnMajor(len1, len2, data);
            #endregion
            #region 并行方式
            //ConcurrentDictionary<int, double[]> dic = new ConcurrentDictionary<int, double[]>();
            //CommonParallel.For(0, matrix1.RowCount, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        double[] v = new double[matrix2.RowCount];
            //        double[] left = matrix1.Row(i).ToArray();
            //        for (int j = 0; j < matrix2.RowCount; j++)
            //        {
            //            double[] right = matrix2.Row(j).ToArray();
            //            double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //            v[j] = d;
            //        }
            //        dic.TryAdd(i, v);
            //    }
            //});
            //Matrix<double> result = Matrix<double>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //// 要不要排序？
            //for (int i = 0; i < dic.Count; i++)
            //{
            //    double[] row;
            //    dic.TryGetValue(i, out row);
            //    result.SetRow(i, row);
            //}
            //return result;
            #endregion
            #region 多线程方式，线程返回值
            //int threadCount = MathNet.Numerics.Control.MaxDegreeOfParallelism;
            //int row_base = matrix1.RowCount / threadCount;
            //int row_remain = matrix1.RowCount % threadCount;
            //List<Tuple<int, int>> indexList = new List<Tuple<int, int>>();
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}
            //Task<Matrix<double>>[] tasks = new Task<Matrix<double>>[threadCount];
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    int to = indexList[i].Item2;
            //    tasks[i] = Task<Matrix<double>>.Factory.StartNew(() =>
            //    {
            //        Matrix<double> m = Matrix<double>.Build.Dense(to - from, matrix2.RowCount);
            //        for (int p = 0; p < m.RowCount; p++)
            //        {
            //            double[] left = matrix1.Row(p).ToArray();
            //            for (int q = 0; q < m.ColumnCount; q++)
            //            {
            //                double[] right = matrix2.Row(q).ToArray();
            //                double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //                m.At(p, q, d);
            //            }
            //        }
            //        return m;
            //    });
            //}
            //Task.WaitAll(tasks);
            //Matrix<double> result = Matrix<double>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    var submatrix = tasks[i].Result;
            //    result.SetSubMatrix(from, 0, submatrix);
            //}
            //return result;

            //threadCount = (threadCount > MathNet.Numerics.Control.MaxDegreeOfParallelism || threadCount <= 0) ?
            //    MathNet.Numerics.Control.MaxDegreeOfParallelism :
            //    threadCount;
            //if (matrix1.RowCount <= threadCount)
            //    threadCount = 1;
            //int row_base = matrix1.RowCount / threadCount;
            //int row_remain = matrix1.RowCount % threadCount;

            //List<Tuple<int, int>> indexList = new List<Tuple<int, int>>();
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}

            //Task<Matrix<double>>[] tasks = new Task<Matrix<double>>[threadCount];
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    int to = indexList[i].Item2;
            //    tasks[i] = Task<Matrix<double>>.Factory.StartNew(() =>
            //    {
            //        int len1 = to - from, len2 = matrix2.RowCount;
            //        double[] resultrow = new double[len1 * len2];
            //        for (int j = 0; j < len2; j++)
            //        {
            //            int n = 0;
            //            double[] datarow = new double[len1];
            //            double[] right = matrix2.Row(j).ToArray();
            //            for (int k = from; k < to; k++)
            //            {
            //                double[] left = matrix1.Row(k).ToArray();
            //                double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //                datarow[n] = d;
            //                n++;
            //            }
            //            Array.Copy(datarow, 0, resultrow, j * len1, len1);
            //        }
            //        return Matrix<double>.Build.DenseOfColumnMajor(len1, len2, resultrow);
            //    });
            //}
            //Task.WaitAll(tasks);

            //Matrix<double> result = Matrix<double>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    var submatrix = tasks[i].Result;
            //    result.SetSubMatrix(from, 0, submatrix);
            //}
            //return result;
            #endregion
            #region 多线程方式，任务工厂
            //// 确定线程数
            //threadCount = (threadCount > MathNet.Numerics.Control.MaxDegreeOfParallelism || threadCount <= 0) ?
            //    MathNet.Numerics.Control.MaxDegreeOfParallelism :
            //    threadCount;
            //if (matrix1.RowCount <= threadCount)
            //    threadCount = 1;
            //// 子任务分块
            //int row_base = matrix1.RowCount / threadCount;
            //int row_remain = matrix1.RowCount % threadCount;
            //List<Tuple<int, int>> indexList = new List<Tuple<int, int>>();
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}
            //// 创建一个父任务，由父任务再创建子任务，父任务返回每个子任务的计算结果
            //Task<Matrix<double>[]> parentTask = new Task<Matrix<double>[]>(() =>
            //{
            //    CancellationTokenSource cts = new CancellationTokenSource();
            //    // 任务工厂统一创建子任务
            //    TaskFactory<Matrix<double>> factory = new TaskFactory<Matrix<double>>(cts.Token, TaskCreationOptions.AttachedToParent, TaskContinuationOptions.None, TaskScheduler.Default);
            //    Matrix<double>[] submatrixs = new Matrix<double>[threadCount];
            //    Task<Matrix<double>>[] childTasks = new Task<Matrix<double>>[threadCount];
            //    for (int i = 0; i < threadCount; i++)
            //    {
            //        int n = i;
            //        int from = indexList[n].Item1;
            //        int to = indexList[n].Item2;
            //        childTasks[n] = factory.StartNew(() => submatrixs[n] = GetDoubleSubMatrix(cts.Token, from, to, matrix1, matrix2));
            //        // 这里，任何子任务抛出异常就取消其余子任务
            //        childTasks[n].ContinueWith(t => cts.Cancel(), TaskContinuationOptions.OnlyOnFaulted);
            //    }
            //    return submatrixs;
            //});
            //// 拼装结果矩阵，使用ContinueWith，父任务完成时启动另一个任务去输出结果
            //// 这里，应该也可以用parentTask.Wait(); parentTask.Result; 得到结果后再输出
            //Matrix<double> result = Matrix<double>.Build.Dense(matrix1.RowCount, matrix2.RowCount);
            //var wt = parentTask.ContinueWith(task =>
            //{
            //    for (int i = 0; i < threadCount; i++)
            //    {
            //        int from = indexList[i].Item1;
            //        var submatrix = task.Result[i];
            //        if (submatrix != null)
            //            result.SetSubMatrix(from, 0, submatrix);
            //    }
            //});
            //parentTask.Start();
            //wt.Wait();
            //return result;
            #endregion
        }
        /// <summary>
        /// 子任务事件，得到子矩阵
        /// </summary>
        [Obsolete("单线程计算欧氏距离用不到")]
        static Matrix<double> GetDoubleSubMatrix(CancellationToken ct, int from, int to, Matrix<double> matrix1, Matrix<double> matrix2)
        {
            try
            {
                int len1 = to - from, len2 = matrix2.RowCount;
                double[] resultrow = new double[len1 * len2];
                for (int j = 0; j < len2; j++)
                {
                    int n = 0;
                    double[] datarow = new double[len1];
                    double[] right = matrix2.Row(j).ToArray();
                    for (int k = from; k < to; k++)
                    {
                        double[] left = matrix1.Row(k).ToArray();
                        double d = MathNet.Numerics.Distance.Euclidean(left, right);
                        datarow[n] = d;
                        n++;
                    }
                    Array.Copy(datarow, 0, resultrow, j * len1, len1);
                }
                return Matrix<double>.Build.DenseOfColumnMajor(len1, len2, resultrow);
            }
            catch
            {
                ct.ThrowIfCancellationRequested();
            }
            return null;
        }
        #endregion
        #endregion

        #region regstats 一元线性回归
        /// <summary>
        /// 暂时只拟合一元线性回归
        /// </summary>
        /// <param name="x">自变量</param>
        /// <param name="y">因变量</param>
        /// <returns>返回回归显著性检验</returns>
        public static LinearRegression Regstats(double[] x, double[] y)
        {
            if (x == null || y == null)
            {
                return null;
            }
            if (x.Length <= 0 || y.Length <= 0 || x.Length != y.Length)
            {
                return null;
            }
            LinearRegression result = new LinearRegression();
            int i = 0, l = x.Length;
            double Lxx = 0.0, Lxy = 0.0, xa = 0.0, ya = 0.0;

            xa = x.Average();// x的平均值
            ya = y.Average();// y的平均值

            for (i = 0; i < l; i++)
            {
                Lxx += (x[i] - xa) * (x[i] - xa);// Sum((X-xa)*(X-xa))
                Lxy += (x[i] - xa) * (y[i] - ya);// Sum((X-xa)*(Y-ya))
            }

            // 回归系数
            // Y = Item1 + Item2 * X1
            result.Coefficients.Item2 = Lxy / Lxx;
            result.Coefficients.Item1 = ya - result.Coefficients.Item2 * xa;

            // 方差分析
            for (i = 0; i < l; i++)
            {
                Lxy = result.Coefficients.Item1 + result.Coefficients.Item2 * x[i];
                result.RegressionSumOfSquares += (Lxy - ya) * (Lxy - ya);// 回归平方和
                result.ResidualSumOfSquares += (y[i] - Lxy) * (y[i] - Lxy);// 剩余平方和
            }

            result.SumOfSquaresOfDeviations = result.RegressionSumOfSquares + result.ResidualSumOfSquares;// 离差平方和
            result.RegressionVariance = result.RegressionSumOfSquares;// 回归方差
            result.ResidualVariance = result.ResidualSumOfSquares / (l - 2);// 剩余方差
            result.StandardError = Math.Sqrt(result.ResidualVariance);// 标准误差
            result.VerificationF = result.RegressionVariance / result.ResidualVariance;// F检验
            result.CorrelationCoefficient = Math.Sqrt(result.RegressionSumOfSquares / result.SumOfSquaresOfDeviations);// 相关系数

            return result;
        }
        #endregion

        #region repmat 矩阵以块的形式复制
        #region Float
        /// <summary>
        /// 矩阵以块的形式复制
        /// </summary>
        /// <param name="matrix">原始矩阵</param>
        /// <param name="m">行倍数</param>
        /// <param name="n">列倍数</param>
        public static Matrix<float> Repmat(Matrix<float> matrix, int m, int n)
        {
            int rowNum = matrix.RowCount,
                colNum = matrix.ColumnCount,
                newRowNum = rowNum * m,
                newColNum = colNum * n,
                i, j, p, q;
            Matrix<float> result = Matrix<float>.Build.Dense(newRowNum, newColNum);
            for (i = 0; i < newRowNum; i++)
            {
                p = i % rowNum;
                for (j = 0; j < newColNum; j++)
                {
                    q = j % colNum;
                    result.At(i, j, matrix.At(p, q));
                }
            }
            return result;
        }
        #endregion

        #region Double
        /// <summary>
        /// 矩阵以块的形式复制
        /// </summary>
        /// <param name="matrix">原始矩阵</param>
        /// <param name="m">行倍数</param>
        /// <param name="n">列倍数</param>
        public static Matrix<double> Repmat(Matrix<double> matrix, int m, int n)
        {
            int rowNum = matrix.RowCount,
                colNum = matrix.ColumnCount,
                newRowNum = rowNum * m,
                newColNum = colNum * n,
                i, j, p, q;
            Matrix<double> result = Matrix<double>.Build.Dense(newRowNum, newColNum);
            for (i = 0; i < newRowNum; i++)
            {
                p = i % rowNum;
                for (j = 0; j < newColNum; j++)
                {
                    q = j % colNum;
                    result.At(i, j, matrix.At(p, q));
                }
            }
            return result;
        }
        #endregion
        #endregion

        #region reshape 调整矩阵的行和列
        #region Float
        /// <summary>
        /// 调整矩阵的行和列
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="m">新行数</param>
        /// <param name="n">新列数</param>
        public static Matrix<float> Reshape(Matrix<float> matrix, int m, int n)
        {
            if (matrix.RowCount * matrix.ColumnCount != m * n)
            {
                throw new Exception("To RESHAPE the number of elements must not change.");
            }
            float[] data = matrix.ToColumnWiseArray();
            return Matrix<float>.Build.DenseOfColumnMajor(m, n, data);
        }
        #endregion

        #region Double
        /// <summary>
        /// 调整矩阵的行和列
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="m">新行数</param>
        /// <param name="n">新列数</param>
        public static Matrix<double> Reshape(Matrix<double> matrix, int m, int n)
        {
            if (matrix.RowCount * matrix.ColumnCount != m * n)
            {
                throw new Exception("To RESHAPE the number of elements must not change.");
            }
            double[] data = matrix.ToColumnWiseArray();
            return Matrix<double>.Build.DenseOfColumnMajor(m, n, data);
        }
        #endregion
        #endregion

        #region sort 排序
        #region Float
        /// <summary>
        /// 对数组排序
        /// </summary>
        /// <param name="array">数组</param>
        /// <param name="model">排序类型，默认升序</param>
        public static Tuple<float[], int[]> SortArray(float[] array, SortType model = SortType.Ascend)
        {
            int len = array.Length;
            float[] source = (float[])array.Clone();
            float[] valueResult = (float[])array.Clone();
            int[] indexResult = new int[len];
            Sorting.Sort<float>(valueResult, new SingleComparer(model));
            int index;
            for (int i = 0; i < len; i++)
            {
                index = Array.IndexOf(source, valueResult[i]);
                indexResult[i] = index + 1;
                source[index] = float.NaN;
            }
            source = new float[0];
            return Tuple.Create<float[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 对矩阵列排序
        /// 返回排序后的各列集合和各行索引，从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="model">排序类型，默认升序</param>
        public static Tuple<float[][], int[][]> SortColumnWise(Matrix<float> matrix, SortType model = SortType.Ascend)
        {
            int rowNum = matrix.RowCount;
            int colNum = matrix.ColumnCount;
            float[][] valueResult = new float[colNum][];
            int[][] indexResult = new int[colNum][];
            float[] valueArray;
            int[] indexArray;
            float[] columnData;
            int index;
            if (rowNum == 1)
            {
                for (int i = 0; i < colNum; i++)
                {
                    valueArray = new float[] { matrix.At(0, i) };
                    indexArray = new int[] { 1 };
                    valueResult[i] = valueArray;
                    indexResult[i] = indexArray;
                }
            }
            else if (colNum == 1)
            {
                columnData = matrix.Column(0).ToArray();
                valueArray = (float[])columnData.Clone();
                Sorting.Sort<float>(valueArray, new SingleComparer(model));
                valueResult[0] = valueArray;
                indexArray = new int[rowNum];
                for (int i = 0; i < rowNum; i++)
                {
                    index = Array.IndexOf(columnData, valueArray[i]);
                    indexArray[i] = index + 1;
                    columnData[index] = float.NaN;
                }
                indexResult[0] = indexArray;
            }
            else
            {
                for (int i = 0; i < colNum; i++)
                {
                    columnData = matrix.Column(i).ToArray();
                    valueArray = (float[])columnData.Clone();
                    Sorting.Sort<float>(valueArray, new SingleComparer(model));
                    valueResult[i] = valueArray;
                    indexArray = new int[rowNum];
                    for (int j = 0; j < rowNum; j++)
                    {
                        index = Array.IndexOf(columnData, valueArray[j]);
                        indexArray[j] = index + 1;
                        columnData[index] = float.NaN;
                    }
                    indexResult[i] = indexArray;
                }
            }
            valueArray = new float[0]; indexArray = new int[0]; columnData = new float[0];
            return Tuple.Create<float[][], int[][]>(valueResult, indexResult);
        }
        /// <summary>
        /// 对矩阵行排序
        /// 返回排序后的各行集合和各行索引，从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="model">排序类型，默认升序</param>
        public static Tuple<float[][], int[][]> SortRowWise(Matrix<float> matrix, SortType model = SortType.Ascend)
        {
            int rowNum = matrix.RowCount;
            int colNum = matrix.ColumnCount;
            float[][] valueResult = new float[rowNum][];
            int[][] indexResult = new int[rowNum][];
            float[] valueArray;
            int[] indexArray;
            float[] rowData;
            int index;
            if (rowNum == 1)
            {
                rowData = matrix.Row(0).ToArray();
                valueArray = (float[])rowData.Clone();
                Sorting.Sort<float>(valueArray, new SingleComparer(model));
                valueResult[0] = valueArray;
                indexArray = new int[colNum];
                for (int i = 0; i < colNum; i++)
                {
                    index = Array.IndexOf(rowData, valueArray[i]);
                    indexArray[i] = index + 1;
                    rowData[index] = float.NaN;
                }
                indexResult[0] = indexArray;
            }
            else if (colNum == 1)
            {
                for (int i = 0; i < rowNum; i++)
                {
                    valueArray = new float[] { matrix.At(i, 0) };
                    indexArray = new int[] { 1 };
                    valueResult[i] = valueArray;
                    indexResult[i] = indexArray;
                }
            }
            else
            {
                for (int i = 0; i < rowNum; i++)
                {
                    rowData = matrix.Row(i).ToArray();
                    valueArray = (float[])rowData.Clone();
                    Sorting.Sort<float>(valueArray, new SingleComparer(model));
                    valueResult[i] = valueArray;
                    indexArray = new int[colNum];
                    for (int j = 0; j < colNum; j++)
                    {
                        index = Array.IndexOf(rowData, valueArray[j]);
                        indexArray[j] = index + 1;
                        rowData[index] = float.NaN;
                    }
                    indexResult[i] = indexArray;
                }
            }
            valueArray = new float[0]; indexArray = new int[0]; rowData = new float[0];
            return Tuple.Create<float[][], int[][]>(valueResult, indexResult);
        }
        #endregion

        #region Double
        /// <summary>
        /// 对数组排序
        /// </summary>
        /// <param name="array">数组</param>
        /// <param name="model">排序类型，默认升序</param>
        public static Tuple<double[], int[]> SortArray(double[] array, SortType model = SortType.Ascend)
        {
            int len = array.Length;
            double[] source = (double[])array.Clone();
            double[] valueResult = (double[])array.Clone();
            int[] indexResult = new int[len];
            Sorting.Sort<double>(valueResult, new DoubleComparer(model));
            int index;
            for (int i = 0; i < len; i++)
            {
                index = Array.IndexOf(source, valueResult[i]);
                indexResult[i] = index + 1;
                source[index] = double.NaN;
            }
            source = new double[0];
            return Tuple.Create<double[], int[]>(valueResult, indexResult);
        }
        /// <summary>
        /// 对矩阵列排序
        /// 返回排序后的各列集合和各行索引，从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="model">排序类型，默认升序</param>
        public static Tuple<double[][], int[][]> SortColumnWise(Matrix<double> matrix, SortType model = SortType.Ascend)
        {
            int rowNum = matrix.RowCount;
            int colNum = matrix.ColumnCount;
            double[][] valueResult = new double[colNum][];
            int[][] indexResult = new int[colNum][];
            double[] valueArray;
            int[] indexArray;
            double[] columnData;
            int index;
            if (rowNum == 1)
            {
                for (int i = 0; i < colNum; i++)
                {
                    valueArray = new double[] { matrix.At(0, i) };
                    indexArray = new int[] { 1 };
                    valueResult[i] = valueArray;
                    indexResult[i] = indexArray;
                }
            }
            else if (colNum == 1)
            {
                columnData = matrix.Column(0).ToArray();
                valueArray = (double[])columnData.Clone();
                Sorting.Sort<double>(valueArray, new DoubleComparer(model));
                valueResult[0] = valueArray;
                indexArray = new int[rowNum];
                for (int i = 0; i < rowNum; i++)
                {
                    index = Array.IndexOf(columnData, valueArray[i]);
                    indexArray[i] = index + 1;
                    columnData[index] = double.NaN;
                }
                indexResult[0] = indexArray;
            }
            else
            {
                for (int i = 0; i < colNum; i++)
                {
                    columnData = matrix.Column(i).ToArray();
                    valueArray = (double[])columnData.Clone();
                    Sorting.Sort<double>(valueArray, new DoubleComparer(model));
                    valueResult[i] = valueArray;
                    indexArray = new int[rowNum];
                    for (int j = 0; j < rowNum; j++)
                    {
                        index = Array.IndexOf(columnData, valueArray[j]);
                        indexArray[j] = index + 1;
                        columnData[index] = double.NaN;
                    }
                    indexResult[i] = indexArray;
                }
            }
            valueArray = new double[0]; indexArray = new int[0]; columnData = new double[0];
            return Tuple.Create<double[][], int[][]>(valueResult, indexResult);
        }
        /// <summary>
        /// 对矩阵行排序
        /// 返回排序后的各行集合和各行索引，从1开始
        /// </summary>
        /// <param name="matrix">矩阵</param>
        /// <param name="model">排序类型，默认升序</param>
        public static Tuple<double[][], int[][]> SortRowWise(Matrix<double> matrix, SortType model = SortType.Ascend)
        {
            int rowNum = matrix.RowCount;
            int colNum = matrix.ColumnCount;
            double[][] valueResult = new double[rowNum][];
            int[][] indexResult = new int[rowNum][];
            double[] valueArray;
            int[] indexArray;
            double[] rowData;
            int index;
            if (rowNum == 1)
            {
                rowData = matrix.Row(0).ToArray();
                valueArray = (double[])rowData.Clone();
                Sorting.Sort<double>(valueArray, new DoubleComparer(model));
                valueResult[0] = valueArray;
                indexArray = new int[colNum];
                for (int i = 0; i < colNum; i++)
                {
                    index = Array.IndexOf(rowData, valueArray[i]);
                    indexArray[i] = index + 1;
                    rowData[index] = double.NaN;
                }
                indexResult[0] = indexArray;
            }
            else if (colNum == 1)
            {
                for (int i = 0; i < rowNum; i++)
                {
                    valueArray = new double[] { matrix.At(i, 0) };
                    indexArray = new int[] { 1 };
                    valueResult[i] = valueArray;
                    indexResult[i] = indexArray;
                }
            }
            else
            {
                for (int i = 0; i < rowNum; i++)
                {
                    rowData = matrix.Row(i).ToArray();
                    valueArray = (double[])rowData.Clone();
                    Sorting.Sort<double>(valueArray, new DoubleComparer(model));
                    valueResult[i] = valueArray;
                    indexArray = new int[colNum];
                    for (int j = 0; j < colNum; j++)
                    {
                        index = Array.IndexOf(rowData, valueArray[j]);
                        indexArray[j] = index + 1;
                        rowData[index] = double.NaN;
                    }
                    indexResult[i] = indexArray;
                }
            }
            valueArray = new double[0]; indexArray = new int[0]; rowData = new double[0];
            return Tuple.Create<double[][], int[][]>(valueResult, indexResult);
        }
        #endregion
        #endregion

        #region sum 求和
        #region Float
        /// <summary>
        /// 对矩阵列求和
        /// 返回数组
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static float[] SumColumnWise(Matrix<float> matrix)
        {
            if (matrix.RowCount == 1)
            {
                return matrix.ToColumnWiseArray();
            }
            else if (matrix.ColumnCount == 1)
            {
                return new float[1] { matrix.ToColumnWiseArray().Sum() };
            }
            else
            {
                int len = matrix.ColumnCount;
                float[] result = new float[len];
                //CommonParallel.For(0, len, delegate(int a, int b)
                //{
                //    for (int i = a; i < b; i++)
                //    {
                //        result[i] = matrix.Column(i).ToArray().Sum();
                //    }
                //});
                for (int i = 0; i < len; i++)
                {
                    result[i] = matrix.Column(i).ToArray().Sum();
                }
                return result;
            }
        }
        /// <summary>
        /// 对矩阵行求和
        /// 返回数组
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static float[] SumRowWise(Matrix<float> matrix)
        {
            if (matrix.ColumnCount == 1)
            {
                return matrix.ToColumnWiseArray();
            }
            else if (matrix.RowCount == 1)
            {
                return new float[1] { matrix.ToColumnWiseArray().Sum() };
            }
            else
            {
                int len = matrix.RowCount;
                float[] result = new float[len];
                //CommonParallel.For(0, len, delegate(int a, int b)
                //{
                //    for (int i = a; i < b; i++)
                //    {
                //        result[i] = matrix.Row(i).ToArray().Sum();
                //    }
                //});
                for (int i = 0; i < len; i++)
                {
                    result[i] = matrix.Row(i).ToArray().Sum();
                }
                return result;
            }
        }
        #endregion

        #region Double
        /// <summary>
        /// 对矩阵列求和
        /// 返回数组
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static double[] SumColumnWise(Matrix<double> matrix)
        {
            if (matrix.RowCount == 1)
            {
                return matrix.ToColumnWiseArray();
            }
            else if (matrix.ColumnCount == 1)
            {
                return new double[1] { matrix.ToColumnWiseArray().Sum() };
            }
            else
            {
                int len = matrix.ColumnCount;
                double[] result = new double[len];
                //CommonParallel.For(0, len, delegate(int a, int b)
                //{
                //    for (int i = a; i < b; i++)
                //    {
                //        result[i] = matrix.Column(i).ToArray().Sum();
                //    }
                //});
                for (int i = 0; i < len; i++)
                {
                    result[i] = matrix.Column(i).ToArray().Sum();
                }
                return result;
            }
        }
        /// <summary>
        /// 对矩阵行求和
        /// 返回数组
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static double[] SumRowWise(Matrix<double> matrix)
        {
            if (matrix.ColumnCount == 1)
            {
                return matrix.ToColumnWiseArray();
            }
            else if (matrix.RowCount == 1)
            {
                return new double[1] { matrix.ToColumnWiseArray().Sum() };
            }
            else
            {
                int len = matrix.RowCount;
                double[] result = new double[len];
                //CommonParallel.For(0, len, delegate(int a, int b)
                //{
                //    for (int i = a; i < b; i++)
                //    {
                //        result[i] = matrix.Row(i).ToArray().Sum();
                //    }
                //});
                for (int i = 0; i < len; i++)
                {
                    result[i] = matrix.Row(i).ToArray().Sum();
                }
                return result;
            }
        }
        #endregion
        #endregion

        #region tabulate 数组出现的次数和频率
        #region Int
        /// <summary>
        /// tabulate 函数
        /// 非连续
        /// </summary>
        /// <param name="x">数组</param>
        /// <returns>返回三元组。第一为唯一数值，第二为数值出现的次数，第三为出现频率</returns>
        public static Tuple<int[], int[], double[]> TabulateUnContinuous(int[] x)
        {
            if (x == null || x.Length <= 0)
            {
                return null;
            }
            int len = x.Length;
            Sorting.Sort<int>(x);
            List<int> lstVal = new List<int>();
            List<int> lstCount = new List<int>();
            List<double> lstPercent = new List<double>();
            Hashtable hshArray = new Hashtable();
            int num;
            int item;
            for (int i = 0; i < len; i++)
            {
                item = x[i];
                if (!hshArray.ContainsKey(item))
                {
                    hshArray.Add(item, item);
                    lstVal.Add(item);
                    num = x.Count((int a) => a == item);
                    lstCount.Add(num);
                    lstPercent.Add((double)num / len);
                }
            }
            hshArray.Clear();
            hshArray = null;
            return Tuple.Create<int[], int[], double[]>(lstVal.ToArray(), lstCount.ToArray(), lstPercent.ToArray());
        }
        /// <summary>
        /// tabulate 函数
        /// 连续
        /// </summary>
        /// <param name="x">数组</param>
        /// <returns>返回三元组。第一为唯一数值，第二为数值出现的次数，第三为出现频率</returns>
        public static Tuple<int[], int[], double[]> TabulateContinuous(int[] x)
        {
            if (x == null || x.Length <= 0)
            {
                return null;
            }
            int len = x.Length;
            Sorting.Sort<int>(x);

            List<int> lstVal = new List<int>();
            List<int> lstCount = new List<int>();
            List<double> lstPercent = new List<double>();
            int max = x[len - 1];
            int num;
            for (int i = 1; i <= max; i++)
            {
                lstVal.Add(i);
                num = x.Count((int a) => a == i);
                lstCount.Add(num);
                lstPercent.Add((double)num / len);
            }
            return Tuple.Create<int[], int[], double[]>(lstVal.ToArray(), lstCount.ToArray(), lstPercent.ToArray());
        }
        #endregion

        #region Float
        /// <summary>
        /// tabulate 函数
        /// 非连续
        /// </summary>
        /// <param name="x">数组</param>
        /// <returns>返回三元组。第一为唯一数值，第二为数值出现的次数，第三为出现频率</returns>
        public static Tuple<float[], int[], double[]> TabulateUnContinuous(float[] x)
        {
            if (x == null || x.Length <= 0)
            {
                return null;
            }
            int len = x.Length;
            Sorting.Sort<float>(x);
            List<float> lstVal = new List<float>();
            List<int> lstCount = new List<int>();
            List<double> lstPercent = new List<double>();
            Hashtable hshArray = new Hashtable();
            int num;
            float item;
            for (int i = 0; i < len; i++)
            {
                item = x[i];
                if (!hshArray.ContainsKey(item))
                {
                    hshArray.Add(item, item);
                    lstVal.Add(item);
                    num = x.Count((float a) => a == item);
                    lstCount.Add(num);
                    lstPercent.Add((double)num / len);
                }
            }
            hshArray.Clear();
            hshArray = null;
            return Tuple.Create<float[], int[], double[]>(lstVal.ToArray(), lstCount.ToArray(), lstPercent.ToArray());
        }
        /// <summary>
        /// tabulate 函数
        /// 连续
        /// </summary>
        /// <param name="x">数组</param>
        /// <returns>返回三元组。第一为唯一数值，第二为数值出现的次数，第三为出现频率</returns>
        public static Tuple<float[], int[], double[]> TabulateContinuous(float[] x)
        {
            if (x == null || x.Length <= 0)
            {
                return null;
            }
            int len = x.Length;
            Sorting.Sort<float>(x);
            List<float> lstVal = new List<float>();
            List<int> lstCount = new List<int>();
            List<double> lstPercent = new List<double>();
            float max = x[len - 1];
            int num;
            for (float i = 1.0f; i <= max; i += 1.0f)
            {
                lstVal.Add(i);
                num = x.Count((float a) => a == i);
                lstCount.Add(num);
                lstPercent.Add((double)num / len);
            }
            return Tuple.Create<float[], int[], double[]>(lstVal.ToArray(), lstCount.ToArray(), lstPercent.ToArray());
        }
        #endregion

        #region Double
        /// <summary>
        /// tabulate 函数
        /// 非连续
        /// </summary>
        /// <param name="x">数组</param>
        /// <returns>返回三元组。第一为唯一数值，第二为数值出现的次数，第三为出现频率</returns>
        public static Tuple<double[], int[], double[]> TabulateUnContinuous(double[] x)
        {
            if (x == null || x.Length <= 0)
            {
                return null;
            }
            int len = x.Length;
            Sorting.Sort<double>(x);
            List<double> lstVal = new List<double>();
            List<int> lstCount = new List<int>();
            List<double> lstPercent = new List<double>();
            Hashtable hshArray = new Hashtable();
            int num;
            double item;
            for (int i = 0; i < len; i++)
            {
                item = x[i];
                if (!hshArray.ContainsKey(item))
                {
                    hshArray.Add(item, item);
                    lstVal.Add(item);
                    num = x.Count((double a) => a == item);
                    lstCount.Add(num);
                    lstPercent.Add((double)num / len);
                }
            }
            hshArray.Clear();
            hshArray = null;
            return Tuple.Create<double[], int[], double[]>(lstVal.ToArray(), lstCount.ToArray(), lstPercent.ToArray());
        }
        /// <summary>
        /// tabulate 函数
        /// 连续
        /// </summary>
        /// <param name="x">数组</param>
        /// <returns>返回三元组。第一为唯一数值，第二为数值出现的次数，第三为出现频率</returns>
        public static Tuple<double[], int[], double[]> TabulateContinuous(double[] x)
        {
            if (x == null || x.Length <= 0)
            {
                return null;
            }
            int len = x.Length;
            Sorting.Sort<double>(x);
            List<double> lstVal = new List<double>();
            List<int> lstCount = new List<int>();
            List<double> lstPercent = new List<double>();
            double max = x[len - 1];
            int num;
            for (double i = 1.0; i <= max; i += 1.0)
            {
                lstVal.Add(i);
                num = x.Count((double a) => a == i);
                lstCount.Add(num);
                lstPercent.Add((double)num / len);
            }
            return Tuple.Create<double[], int[], double[]>(lstVal.ToArray(), lstCount.ToArray(), lstPercent.ToArray());
        }
        #endregion
        #endregion

        #region unique 唯一元素和唯一行
        #region Float
        /// <summary>
        /// 数组返回唯一元素
        /// </summary>
        /// <param name="array">数组</param>
        public static float[] UniqueArray(float[] array)
        {
            float[] result = array.Distinct().ToArray();
            Sorting.Sort<float>(result);
            return result;
        }
        /// <summary>
        /// 返回唯一元素
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static float[] UniqueElement(Matrix<float> matrix)
        {
            float[] data = matrix.ToColumnWiseArray();
            //Hashtable htSyn = Hashtable.Synchronized(new Hashtable());
            //CommonParallel.For(0, data.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        float key = data[i];
            //        if (!htSyn.ContainsKey(key))
            //        {
            //            htSyn.Add(key, key);
            //        }
            //    }
            //});
            Hashtable htSyn = new Hashtable();
            for (int i = 0; i < data.Length; i++)
            {
                float key = data[i];
                if (!htSyn.ContainsKey(key))
                {
                    htSyn.Add(key, key);
                }
            }
            float[] result = new float[htSyn.Count];
            htSyn.Values.CopyTo(result, 0);
            Sorting.Sort<float>(result);
            htSyn.Clear();
            htSyn = null;
            data = new float[0];
            return result;
        }
        /// <summary>
        /// 返回唯一行
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Matrix<float> UniqueRow(Matrix<float> matrix)
        {
            int colNum = matrix.ColumnCount;
            //Hashtable htSyn = Hashtable.Synchronized(new Hashtable());
            //CommonParallel.For(0, matrix.RowCount, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        float[] data = matrix.Row(i).ToArray();
            //        SingleVector vector = new SingleVector(data);
            //        string key = vector.ToString();
            //        if (!htSyn.ContainsKey(key))
            //        {
            //            htSyn.Add(key, vector);
            //        }
            //    }
            //});
            Hashtable htSyn = new Hashtable();
            float[] data;
            for (int i = 0; i < matrix.RowCount; i++)
            {
                data = matrix.Row(i).ToArray();
                SingleVector vector = new SingleVector(data);
                string key = vector.ToString();
                if (!htSyn.ContainsKey(key))
                {
                    htSyn.Add(key, vector);
                }
            }
            SingleVector[] udata = new SingleVector[htSyn.Count];
            htSyn.Values.CopyTo(udata, 0);
            Sorting.Sort<SingleVector>(udata);
            htSyn.Clear();
            htSyn = null;
            data = new float[0];
            int rowNum = udata.Length;
            var result = Matrix<float>.Build.Dense(rowNum, colNum);
            for (int i = 0; i < rowNum; i++)
            {
                result.SetRow(i, udata[i].Data);
            }
            udata = new SingleVector[0];
            return result;
        }
        #endregion

        #region Double
        /// <summary>
        /// 数组返回唯一元素
        /// </summary>
        /// <param name="array">数组</param>
        public static double[] UniqueArray(double[] array)
        {
            double[] result = array.Distinct().ToArray();
            Sorting.Sort<double>(result);
            return result;
        }
        /// <summary>
        /// 矩阵返回唯一元素
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static double[] UniqueElement(Matrix<double> matrix)
        {
            double[] data = matrix.ToColumnWiseArray();
            //Hashtable htSyn = Hashtable.Synchronized(new Hashtable());
            //CommonParallel.For(0, data.Length, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        double key = data[i];
            //        if (!htSyn.ContainsKey(key))
            //        {
            //            htSyn.Add(key, key);
            //        }
            //    }
            //});
            Hashtable htSyn = new Hashtable();
            for (int i = 0; i < data.Length; i++)
            {
                double key = data[i];
                if (!htSyn.ContainsKey(key))
                {
                    htSyn.Add(key, key);
                }
            }
            double[] result = new double[htSyn.Count];
            htSyn.Values.CopyTo(result, 0);
            Sorting.Sort<double>(result);
            htSyn.Clear();
            htSyn = null;
            data = new double[0];
            return result;
        }
        /// <summary>
        /// 矩阵返回唯一行
        /// </summary>
        /// <param name="matrix">矩阵</param>
        public static Matrix<double> UniqueRow(Matrix<double> matrix)
        {
            int colNum = matrix.ColumnCount;
            //Hashtable htSyn = Hashtable.Synchronized(new Hashtable());
            //CommonParallel.For(0, matrix.RowCount, delegate(int a, int b)
            //{
            //    for (int i = a; i < b; i++)
            //    {
            //        double[] data = matrix.Row(i).ToArray();
            //        DoubleVector vector = new DoubleVector(data);
            //        string key = vector.ToString();
            //        if (!htSyn.ContainsKey(key))
            //        {
            //            htSyn.Add(key, vector);
            //        }
            //    }
            //});
            Hashtable htSyn = new Hashtable();
            double[] data;
            for (int i = 0; i < matrix.RowCount; i++)
            {
                data = matrix.Row(i).ToArray();
                DoubleVector vector = new DoubleVector(data);
                string key = vector.ToString();
                if (!htSyn.ContainsKey(key))
                {
                    htSyn.Add(key, vector);
                }
            }
            DoubleVector[] udata = new DoubleVector[htSyn.Count];
            htSyn.Values.CopyTo(udata, 0);
            Sorting.Sort<DoubleVector>(udata);
            htSyn.Clear();
            htSyn = null;
            data = new double[0];
            int rowNum = udata.Length;
            var result = Matrix<double>.Build.Dense(rowNum, colNum);
            for (int i = 0; i < rowNum; i++)
            {
                result.SetRow(i, udata[i].Data);
            }
            udata = new DoubleVector[0];
            return result;
        }
        #endregion
        #endregion
    }
}
