using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Diagnostics;
using System.Drawing;
using System.Linq;
using System.Runtime.InteropServices;
using System.Security.Cryptography;
using System.Text;
using System.Threading;
using System.Threading.Tasks;
using System.Transactions;

namespace MulLinearRegression
{
    
    
    public enum NormType
    {
    
    
        L0,
        L1,
        L2,
        LP,
        P_INF,
        M_INF,
        INF,
        F
    }

    public class Matrix
    {
    
    

        ///   <summary> 
        ///   矩阵的转置 
        ///   </summary> 
        ///   <param   name= "iMatrix "> </param> 
        public static double[,] Transpose(double[,] iMatrix)
        {
    
    
            int row = iMatrix.GetLength(0);
            int column = iMatrix.GetLength(1);
            //double[,] iMatrix = new double[column, row];
            double[,] TempMatrix = new double[row, column];
            double[,] iMatrixT = new double[column, row];
            for (int i = 0; i < row; i++)
            {
    
    
                for (int j = 0; j < column; j++)
                {
    
    
                    TempMatrix[i, j] = iMatrix[i, j];
                }
            }
            for (int i = 0; i < column; i++)
            {
    
    
                for (int j = 0; j < row; j++)
                {
    
    
                    iMatrixT[i, j] = TempMatrix[j, i];
                }
            }
            return iMatrixT;

        }

        ///   <summary> 
        ///   矩阵的逆矩阵 
        ///   </summary> 
        ///   <param   name= "iMatrix "> </param> 
        public static double[,] Athwart(double[,] iMatrix)
        {
    
    
            int i = 0;
            int row = iMatrix.GetLength(0);
            double[,] MatrixZwei = new double[row, row * 2];
            double[,] iMatrixInv = new double[row, row];
            for (i = 0; i < row; i++)
            {
    
    
                for (int j = 0; j < row; j++)
                {
    
    
                    MatrixZwei[i, j] = iMatrix[i, j];
                }
            }
            for (i = 0; i < row; i++)
            {
    
    
                for (int j = row; j < row * 2; j++)
                {
    
    
                    MatrixZwei[i, j] = 0;
                    if (i + row == j)
                        MatrixZwei[i, j] = 1;
                }
            }

            for (i = 0; i < row; i++)
            {
    
    
                if (MatrixZwei[i, i] != 0)
                {
    
    
                    double intTemp = MatrixZwei[i, i];
                    for (int j = 0; j < row * 2; j++)
                    {
    
    
                        MatrixZwei[i, j] = MatrixZwei[i, j] / intTemp;
                    }
                }
                for (int j = 0; j < row; j++)
                {
    
    
                    if (j == i)
                        continue;
                    double intTemp = MatrixZwei[j, i];
                    for (int k = 0; k < row * 2; k++)
                    {
    
    
                        MatrixZwei[j, k] = MatrixZwei[j, k] - MatrixZwei[i, k] * intTemp;
                    }
                }
            }

            for (i = 0; i < row; i++)
            {
    
    
                for (int j = 0; j < row; j++)
                {
    
    
                    iMatrixInv[i, j] = MatrixZwei[i, j + row];
                }
            }
            return iMatrixInv;
        }

        ///   <summary> 
        ///   矩阵加法 
        ///   </summary> 
        ///   <param   name= "MatrixEin "> </param> 
        ///   <param   name= "MatrixZwei "> </param> 
        public static double[,] AddMatrix(double[,] MatrixEin, double[,] MatrixZwei)
        {
    
    
            double[,] MatrixResult = new double[MatrixEin.GetLength(0), MatrixZwei.GetLength(1)];
            for (int i = 0; i < MatrixEin.GetLength(0); i++)
                for (int j = 0; j < MatrixZwei.GetLength(1); j++)
                    MatrixResult[i, j] = MatrixEin[i, j] + MatrixZwei[i, j];
            return MatrixResult;
        }

        ///   <summary> 
        ///   矩阵乘法 
        ///   </summary> 
        ///   <param   name= "MatrixEin "> </param> 
        ///   <param   name= "MatrixZwei "> </param> 
        public static double[,] MultiplyMatrix(double[,] MatrixEin, double[,] MatrixZwei)
        {
    
    
            double[,] MatrixResult = new double[MatrixEin.GetLength(0), MatrixZwei.GetLength(1)];
            for (int i = 0; i < MatrixEin.GetLength(0); i++)
            {
    
    
                for (int j = 0; j < MatrixZwei.GetLength(1); j++)
                {
    
    
                    for (int k = 0; k < MatrixEin.GetLength(1); k++)
                    {
    
    
                        MatrixResult[i, j] += MatrixEin[i, k] * MatrixZwei[k, j];
                    }
                }
            }
            return MatrixResult;
        }

        public static double[,] HadamardProduct(double[,] MatrixEin, double[,] MatrixZwei)
        {
    
    
            int Row = MatrixEin.GetLength(0);
            int Col=MatrixEin.GetLength(1);
            double[,] outMatrix = new double[Row, Col];
            for(int i = 0; i < Row; i++)
            {
    
    
                for(int j = 0; j < Col; j++)
                {
    
    
                    outMatrix[i, j] = MatrixEin[i, j] * MatrixZwei[i, j];
                }
            }
            return outMatrix;
        }

        ///   <summary> 
        ///   3阶矩阵对应行列式 的值 
        ///   </summary> 
        ///   <param   name= "MatrixEin "> </param> 
        ///   <returns> </returns> 
        public static double ResultDeterminant(double[,] MatrixEin)
        {
    
    
            return MatrixEin[0, 0] * MatrixEin[1, 1] * MatrixEin[2, 2] + MatrixEin[0, 1] * MatrixEin[1, 2] * MatrixEin[2, 0] + MatrixEin[0, 2] * MatrixEin[1, 0] * MatrixEin[2, 1]
            - MatrixEin[0, 2] * MatrixEin[1, 1] * MatrixEin[2, 0] - MatrixEin[0, 1] * MatrixEin[1, 0] * MatrixEin[2, 2] - MatrixEin[0, 0] * MatrixEin[1, 2] * MatrixEin[2, 1];

        }

        /// <summary>
        /// 生成单位矩阵
        /// </summary>
        /// <param name="Row"></param>
        /// <returns></returns>
        public static double[,] UnitMatrix(int Row)
        {
    
    
            double[,] outMatrix = new double[Row, Row];

            for (int i = 0; i < Row; i++)
                outMatrix[i, i] = 1;
            return outMatrix;
        }
        
        /// <summary>
        ///向量转对角矩阵 
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <returns></returns>
        public static double[,] VectorToDiagonalMatrix(double[,] MatrixEin)
        {
    
    
            int Row = MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            double[] data = SubDatas(MatrixEin, 0, RowGet: Col > Row);
            Row = Row > Col ? Row : Col;
            double[,] outMatrix = new double[Row, Row];
            for(int i = 0; i < Row; i++)
            {
    
    
                outMatrix[i, i] = data[i];
            }

            return outMatrix;
        }

        /// <summary>
        /// 矩阵合并
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="MatrixZwei"></param>
        /// <param name="RowCombine">是否上下合并</param>
        /// <returns></returns>
        public static double[,] MatrixCombine(double[,] MatrixEin, double[,] MatrixZwei, bool RowCombine = false)
        {
    
    
            int rowE = MatrixEin.GetLength(0);
            int colE = MatrixEin.GetLength(1);
            int rowZ = MatrixZwei.GetLength(0);
            int colZ = MatrixZwei.GetLength(1);
            int dstRow = RowCombine ? rowE + rowZ : rowE;
            int dstCol = RowCombine ? colE : colE + colZ;
            double[,] outMatrix = new double[dstRow, dstCol];
            if (rowE *  colE == 0)
            {
    
    
                return (double[,])MatrixZwei.Clone();
            }
            if (rowZ * colZ == 0)
            {
    
    
                return (double[,])MatrixEin.Clone();
            }
            if (RowCombine)
            {
    
    
                Array.Copy(MatrixEin, 0, outMatrix, 0, rowE * colE);
                Array.Copy(MatrixZwei, 0, outMatrix, rowE * colE, rowZ * colZ);
            }
            else
            {
    
    
                outMatrix = new double[dstCol, dstRow];
                MatrixEin = Transpose(MatrixEin);
                MatrixZwei = Transpose(MatrixZwei);
                Array.Copy(MatrixEin, 0, outMatrix, 0, rowE * colE);
                Array.Copy(MatrixZwei, 0, outMatrix, rowE * colE, rowZ * colZ);
                outMatrix = Transpose(outMatrix);
            }
            return outMatrix;
        }

        /// <summary>
        /// 子阵获取
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Start"></param>
        /// <param name="len"></param>
        /// <param name="RowGet"></param>
        /// <returns></returns>
        public static double[,] SubMatrix(double[,] MatrixEin,int Start,int len=1,bool RowGet = false)
        {
    
    
            double[,] temp = RowGet? MatrixEin:Transpose(MatrixEin);
            int Row=temp.GetLength(0);
            int Col=temp.GetLength(1);
            double[,] outMatrix = new double[len, Col];
            Buffer.BlockCopy(temp,sizeof(double)*Start*Col, outMatrix, 0, sizeof(double) * Col * len);
            return RowGet? outMatrix:Transpose(outMatrix);
        }

        public static double[,] SubMatrix(double[,] MatrixEin,Point LeftTop,int Width,int Height)
        {
    
    
            double[,] outMatrix;
            int uRow = LeftTop.Y;
            int lCol = LeftTop.X;
            int dRow = uRow + Height;
            int rCol = lCol + Width;
            outMatrix = SubMatrix(MatrixEin, uRow, Height, RowGet: true);
            outMatrix = SubMatrix(outMatrix, lCol, Width, RowGet: false);
            return outMatrix;
        }

        /// <summary>
        /// 子数组（一维）
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Start"></param>
        /// <param name="len"></param>
        /// <param name="RowGet"></param>
        /// <returns></returns>
        public static double[] SubDatas(double[,] MatrixEin,int Start,int len=1,bool RowGet = false)
        {
    
    
            double[,] temp = RowGet ? (double[,])MatrixEin.Clone() : Transpose(MatrixEin);
            int Row = temp.GetLength(0);
            int Col = temp.GetLength(1);
            double[] data = new double[Col*len];
            Buffer.BlockCopy(temp, Start * Col * sizeof(double), data, 0, Col * len * sizeof(double));
            return data;
        }

        /// <summary>
        /// 二维数据向量生成
        /// </summary>
        /// <param name="Data">填充到向量中的数据</param>
        /// <param name="RowVec">是否行向量</param>
        /// <returns></returns>
        public static double[,] VectorGenerate(double[] Data, bool RowVec = false)
        {
    
    
            double[,] outVector = new double[1, Data.Length];
            Buffer.BlockCopy(Data, 0, outVector, 0, sizeof(double) * Data.Length);
            outVector = RowVec ? outVector : Transpose(outVector);
            return outVector;
        }
        public static double[,] VectorGenerate(int len, double value = 1.0, bool RowVec = false)
        {
    
    
            double[,] array = new double[1, len];
            for (int i = 0; i < len; i++)
            {
    
    
                array[0, i] = value;
            }
            if (RowVec)
            {
    
    
                return array;
            }
            else
            {
    
    
                return Transpose(array);
            }
        }

        /// <summary>
        /// 删除矩阵行列向量
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Id">向量索引</param>
        /// <param name="RowRmv">是否行向量</param>
        /// <returns></returns>
        public static double[,] MatrixRemoveVector(double[,] MatrixEin, int Id, bool RowRmv = false)
        {
    
    
            double[,] uMatrix = MultiplyConst(MatrixEin, 1);
            uMatrix = RowRmv ? uMatrix : Transpose(uMatrix);
            int eRow = uMatrix.GetLength(0);
            int eCol = uMatrix.GetLength(1);
            int dRow = eRow - 1;
            double[,] outMatrix = new double[dRow, eCol];
            Array.Copy(uMatrix, 0, outMatrix, 0, Id * eCol);
            if (Id + 1 < eRow)
            {
    
    
                Array.Copy(uMatrix, (Id + 1) * eCol, outMatrix, Id * eCol,
                    (eRow - Id - 1) * eCol);
            }
            outMatrix = RowRmv ? outMatrix : Transpose(outMatrix);
            return outMatrix;
        }

        /// <summary>
        /// 范数
        /// 向量：L0,L1,L2,LP=>[Σ(|x[i]|^p)]^(1/p);P_INF=>max|x[i],|M_INF=>min|x[i]
        /// 矩阵:L1=>max,colΣ(|x[i,j]|);INF=>max,rowΣ(|x[i,j]|);F=>[Σ[(x[i,j])^2]]^0.5;L2=>[λmax(X'X)]^0.5
        /// </summary>
        /// <param name="Datas"></param>
        /// <param name="Type">范数类型</param>
        /// <param name="P">|x[i,j]|^P</param>
        /// <returns></returns>
        public static double Norm(double[,] MatrixEin, NormType type, double P = 0)
        {
    
    
            double norm = double.MinValue;
            int Row = MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            if (Row * Col == Row || Row * Col == Col)
            {
    
    
                double[] data = SubDatas(MatrixEin, 0, RowGet: Col > Row);
                switch (type)
                {
    
    
                    case NormType.L0:
                        {
    
    
                            norm = data.ToList().Count(x => x != 0);
                            break;
                        }
                    case NormType.L1:
                        {
    
    
                            norm = data.Sum(x => Math.Abs(x));
                            break;
                        }
                    case NormType.L2:
                        {
    
    
                            norm = Math.Sqrt(data.Sum(x => x * x));
                            break;
                        }
                    case NormType.LP:
                        {
    
    
                            norm = Math.Pow(data.Sum(x => Math.Pow(Math.Abs(x), P)), 1.0 / (double)P);
                            break;
                        }
                    case NormType.P_INF:
                        {
    
    
                            norm = data.Max(x => Math.Abs(x));
                            break;
                        }
                    case NormType.M_INF:
                        {
    
    
                            norm = data.Min(x => Math.Abs(x));
                            break;
                        }
                }
            }
            else
            {
    
    
                switch (type)
                {
    
    
                    case NormType.L1:
                        //列元素绝对值和的最大值
                        {
    
    
                            for (int i = 0; i < Col; i++)
                            {
    
    
                                norm = Math.Max(norm, SubDatas(MatrixEin, i, RowGet: false).Sum(x => Math.Abs(x)));
                            }
                            break;
                        }
                    case NormType.INF:
                        //行元素绝对值和的最大值
                        {
    
    
                            for (int i = 0; i < Row; i++)
                            {
    
    
                                norm = Math.Max(norm, SubDatas(MatrixEin, i, RowGet: true).Sum(x => Math.Abs(x)));
                            }
                            break;
                        }
                    case NormType.F:
                        //[Σ(x[i,j]^2)]^0.5
                        {
    
    
                            norm = 0;
                            for (int i = 0; i < Row; i++)
                            {
    
    
                                norm += SubDatas(MatrixEin, i, RowGet: true).Sum(x => x * x);
                            }
                            norm = Math.Sqrt(norm);
                            break;
                        }
                    case NormType.L2:
                        //λmax(ATA)^0.5
                        {
    
    
                            double[,] nnd = Transpose(MatrixEin);
                            nnd = MultiplyMatrix(nnd, MatrixEin);
                            double lambda = MatrixEigenvalue(nnd, ref nnd, Type: 0);
                            norm = Math.Sqrt(lambda);
                            break;
                        }
                }

            }

            return norm;
        }

        /// <summary>
        /// 数据判定
        /// </summary>
        /// <param name="Datas"></param>
        /// <returns>0:标量，1：向量，2：矩阵</returns>
        public static int MatrixTypeJudge(double[,] Datas)
        {
    
    
            int Row = Datas.GetLength(0);
            int Col = Datas.GetLength(1);

            if (Row == 1 && Col == 1)
            {
    
    
                return 0;
            }
            else if (Row * Col == Row || Row * Col == Col)
            {
    
    
                return 1;
            }
            else
            {
    
    
                return 2;
            }
        }

        /// <summary>
        /// 矩阵最大最小值
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Min"></param>
        /// <param name="Max"></param>
        /// <param name="Abs">是否通过绝对值判定</param>
        public static void MinMax(double[,] MatrixEin, out double Min, out double Max)
        {
    
    
            Max = double.MinValue;
            Min = double.MaxValue;
            int Row = MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            for (int i = 0; i < Row; i++)
            {
    
    
                Max = Math.Max(Max, SubDatas(MatrixEin, i, RowGet: true).Max());
                Min = Math.Min(Min, SubDatas(MatrixEin, i, RowGet: true).Min());
            }
        }


        /// <summary>
        /// 矩阵归一化
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="type">
        /// 0:=>X[i,j]^2/Σ(X[i,j]^2);
        /// 1:=>x[i,col]^2/Σ(x[i,col]^2)
        /// 2:=>(x[i,col]-x[col]')/[(Σ((x[i,col]-x[i,col]')^2))^0.5]-standard
        /// 3:=>x[i,col]/max|x[i,col]|
        /// </param>
        /// <returns></returns>
        public static double[,] Normalization(double[,] MatrixEin,int type = 0)
        {
    
    
            int Row=MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            double[,] outMatrix = new double[Row, Col];
            switch (type)
            {
    
    
                case 0: 
                    {
    
    
                        double squSum = Math.Pow(Norm(MatrixEin, NormType.F), 2);
                        for (int i = 0; i < Row; i++)
                        {
    
    
                            for (int j = 0; j < Col; j++)
                            {
    
    
                                outMatrix[i, j] = MatrixEin[i, j] * MatrixEin[i, j] / squSum;
                            }
                        }
                        break;
                    }
                case 1:
                    {
    
    
                        for (int c = 0; c < Col; c++)
                        {
    
    
                            double colSquSum = Math.Pow(Norm(SubMatrix(MatrixEin, c), NormType.L2), 2);
                            for (int r = 0; r < Row; r++)
                            {
    
    
                                outMatrix[r, c] = MatrixEin[r, c] * MatrixEin[r, c] / colSquSum;
                            }
                        }
                        break;
                    }
                case 2:
                    {
    
    
                        double[] ljj = new double[Col];
                        double[] avgCol = new double[Col];
                        for (int c = 0; c < Col; c++)
                        {
    
    
                            avgCol[c] = SubDatas(MatrixEin, c).Average();
                            for (int r = 0; r < Row; r++)
                            {
    
    
                                ljj[c] += Math.Pow(MatrixEin[r, c] - avgCol[c], 2);
                            }
                        }
                        for(int c = 0; c < Col; c++)
                        {
    
    
                            for(int r = 0; r < Row; r++)
                            {
    
    
                                outMatrix[r, c] = (MatrixEin[r, c] - avgCol[c]) / Math.Sqrt(ljj[c]);
                            }
                        }
                        break;
                    }
                case 3:
                    {
    
    
                        for(int c = 0; c < Col; c++)
                        {
    
    
                            double colAbsMax = Norm(SubMatrix(MatrixEin, c), NormType.P_INF);
                            for(int r = 0; r < Row; r++)
                            {
    
    
                                outMatrix[r, c] = MatrixEin[r, c] / colAbsMax;
                            }
                        }
                        break;
                    }
            }

            return outMatrix;
        }


        /// <summary>
        /// 矩阵最大、最小特征值(power-method，invpower-method)
        /// </summary>
        /// <param name="MatrixEin">非奇异矩阵</param>
        /// <param name="Eigenvector">输出特征向量</param>
        /// <param name="epsilon">迭代精度</param>
        /// <param name="Type">0:最大特征，1:最小特征</param>
        /// <returns>特征值</returns>
        public static double MatrixEigenvalue(double[,] MatrixEin, ref double[,] Eigenvector, double epsilon = 1e-6, int Type = 0)
        {
    
    
            MatrixEin = Type == 0 ? MatrixEin : Athwart(MatrixEin);
            double Eigenvalue = double.MaxValue;
            int Row = MatrixEin.GetLength(0);
            double[,] cVec = VectorGenerate(Row, 1, RowVec: false);
            double[,] uVec = VectorGenerate(Row, 1, RowVec: false);
            double diff = double.MaxValue;
            while (diff > epsilon)
            {
    
    
                cVec = MultiplyMatrix(MatrixEin, uVec);
                uVec = Normalization(cVec, 3);
                MinMax(cVec, out double min, out double max);
                diff = Math.Abs(max - Eigenvalue);
                Eigenvalue = max;
                Eigenvector = uVec;
            }
            return Eigenvalue;
        }

        /// <summary>
        /// 雅可比计算对称矩阵全部特征值及相关特征向量(JacobiMethod)
        /// </summary>
        /// <param name="MatrixEin">input matrix</param>
        /// <param name="EigenVector">eigenvectors</param>
        /// <param name="epsilon">precision requirement</param>
        /// <param name="iterations">iterate max times</param>
        /// <returns>TotalEigenvalues</returns>
        public static double[] SymmetricMatricEigenvalud(double[,] MatrixEin, ref double[,] EigenVector, double epsilon = 1e-7, int iterations = 1000)
        {
    
    
            double[,] A = (double[,])MatrixEin.Clone();
            int Row = MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            double[,] Q = UnitMatrix(Row);
            double[,] R = UnitMatrix(Row);
            double max = double.MaxValue;
            int mx = -1, my = -1;
            int times = 0;
            while (max > epsilon && times < iterations)
            {
    
    
                max = double.MinValue;
                for (int i = 0; i < Row; i++)
                {
    
    
                    for (int j = i + 1; j < Col; j++)
                    {
    
    
                        if (Math.Abs(A[i, j]) > max)
                        {
    
    
                            max = Math.Abs(A[i, j]);
                            mx = i;
                            my = j;
                        }
                    }
                }
                times += 1;
                double d = (A[mx, mx] - A[my, my]) / (2.0 * A[mx, my]);
                double t = d > 0 ? -d + Math.Sqrt(d * d + 1) : d == 0 ? -1 : -d - Math.Sqrt(d * d + 1);
                double c = Math.Pow(1 + t * t, -0.5);
                double s = t * c;
                R = UnitMatrix(Row);
                R[mx, mx] = R[my, my] = c;
                R[mx, my] = s;
                R[my, mx] = -s;
                double[,] RT = Transpose(R);
                double[,] B = MultiplyMatrix(R, A);
                A = MultiplyMatrix(B, RT);
                Q = MultiplyMatrix(Q, RT);
            }
            EigenVector = Q;
            double[] Eigenvalue = new double[Col];
            for (int i = 0; i < Row; i++)
            {
    
    
                Eigenvalue[i] = A[i, i];
            }
            return Eigenvalue;
        }

        /// <summary>
        /// 矩阵洗牌
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="shuffled"></param>
        public static void MatrixShuffle(double[,] MatrixEin, out double[,] shuffled)
        {
    
    
            int Row = MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            shuffled = new double[Row, Col];
            int[] randIds = Enumerable.Range(0, Row).OrderBy(x => Guid.NewGuid()).Take(Row).ToArray();
            for (int i = 0; i < randIds.Length; i++)
            {
    
    
                Buffer.BlockCopy(MatrixEin, randIds[i] * Col * sizeof(double), shuffled, i * Col * sizeof(double), Col * sizeof(double));
            }
        }

        /// <summary>
        /// 矩阵数乘
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Coefficient"></param>
        /// <returns></returns>
        public static double[,] MultiplyConst(double[,] MatrixEin, double Coefficient)
        {
    
    
            int Row = MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            double[,] outMatrix = new double[Row, Col];
            for (int i = 0; i < outMatrix.GetLength(0); i++)
            {
    
    
                for (int j = 0; j < outMatrix.GetLength(1); j++)
                {
    
    
                    outMatrix[i, j] = MatrixEin[i, j] * Coefficient;
                }
            }
            return outMatrix;
        }

        /// <summary>
        /// 矩阵内部元素开方
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <returns></returns>
        public static double[,] MatrixElementSqrt(double[,] MatrixEin)
        {
    
    
            int Row = MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            double[,] outMatrix = new double[Row, Col];
            for (int i = 0; i < Row; i++)
            {
    
    
                for (int j = 0; j < Col; j++)
                {
    
    
                    outMatrix[i, j] = Math.Sqrt(MatrixEin[i, j]);
                }
            }
            return outMatrix;
        }

        /// <summary>
        /// 矩阵元素取倒数
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <returns></returns>
        public static double[,] MatrixElementReciprocal(double[,] MatrixEin)
        {
    
    
            int Row = MatrixEin.GetLength(0);
            int Col = MatrixEin.GetLength(1);
            double[,] outMatrix = new double[Row, Col];
            for (int i = 0; i < Row; i++)
            {
    
    
                for (int j = 0; j < Col; j++)
                {
    
    
                    outMatrix[i, j] = 1.0 / MatrixEin[i, j];
                }
            }
            return outMatrix;
        }

        /// <summary>
        /// 向量拉伸矩阵
        /// </summary>
        /// <param name="Vector"></param>
        /// <param name="Size"></param>
        /// <returns></returns>
        public static double[,] VectorStretch(double[,] Vector,int Size = 1)
        {
    
    
            int Row=Vector.GetLength(0);
            int Col=Vector.GetLength(1);
            double[,] temp = Row > Col ? Transpose(Vector) : (double[,])Vector.Clone();
            double[,] outMatrix = new double[Size, Row > Col ? Row : Col];
            for(int i = 0; i < Size; i++)
            {
    
    
                Buffer.BlockCopy(temp, 0, outMatrix, i * sizeof(double) * (Row > Col ? Row : Col), sizeof(double) * (Row > Col ? Row : Col));
            }
            return Row > Col ? Transpose(outMatrix) : outMatrix;
        }
        /// <summary>
        /// 矩阵变形
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Row"></param>
        /// <returns></returns>
        public static double[,] MatrixReshape(double[,] MatrixEin,int Row)
        {
    
    
            int Width = MatrixEin.GetLength(0) * MatrixEin.GetLength(1) / Row;
            double[,] outMatrix = new double[Row,Width];
            for(int i = 0; i < Row; i++)
            {
    
    
                Buffer.BlockCopy(MatrixEin, i * Width * sizeof(double), outMatrix, i * Width * sizeof(double), sizeof(double) * Width);
            }
            return outMatrix;
        }
        /// <summary>
        /// 矩阵sigmoid
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Grad"></param>
        /// <returns></returns>
        public static double[,] MatrixSigmoid(double[,] MatrixEin,bool Grad=false)
        {
    
    
            double[,] outMatrix = new double[MatrixEin.GetLength(0), MatrixEin.GetLength(1)];
            if (!Grad)
            {
    
    
                for (int i = 0; i < outMatrix.GetLength(0); i++)
                    for (int j = 0; j < outMatrix.GetLength(1); j++)
                        outMatrix[i, j] = 1.0 / (1 + Math.Exp(-1 * MatrixEin[i, j]));
            }
            else
            {
    
    
                for (int i = 0; i < outMatrix.GetLength(0); i++)
                    for (int j = 0; j < outMatrix.GetLength(1); j++)
                        outMatrix[i, j] = MatrixEin[i, j] * (1 - MatrixEin[i, j]);
            }
            return outMatrix;
        }
        /// <summary>
        /// 矩阵双曲正切
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Grad"></param>
        /// <returns></returns>
        public static double[,] MatrixTanh(double[,] MatrixEin, bool Grad = false)
        {
    
    
            double[,] outMatrix = new double[MatrixEin.GetLength(0), MatrixEin.GetLength(1)];
            if (!Grad)
            {
    
    
                for (int i = 0; i < outMatrix.GetLength(0); i++)
                    for (int j = 0; j < outMatrix.GetLength(1); j++)
                        outMatrix[i, j] = (1.0 - Math.Exp(-2 * MatrixEin[i, j])) / (1.0 + Math.Exp(-2 * MatrixEin[i, j]));
            }
            else
            {
    
    
                for (int i = 0; i < outMatrix.GetLength(0); i++)
                    for (int j = 0; j < outMatrix.GetLength(1); j++)
                        outMatrix[i, j] = 1-(MatrixEin[i, j]* MatrixEin[i, j]);
            }
            return outMatrix;
        }
        /// <summary>
        /// 矩阵Relu
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="Grad"></param>
        /// <returns></returns>
        public static double[,] MatrixReLU(double[,] MatrixEin, bool Grad = false)
        {
    
    
            double[,] outMatrix = new double[MatrixEin.GetLength(0), MatrixEin.GetLength(1)];
            if (!Grad)
            {
    
    
                for (int i = 0; i < outMatrix.GetLength(0); i++)
                    for (int j = 0; j < outMatrix.GetLength(1); j++)
                        outMatrix[i, j] = Math.Max(0, MatrixEin[i, j]);
            }
            else
            {
    
    
                for (int i = 0; i < outMatrix.GetLength(0); i++)
                    for (int j = 0; j < outMatrix.GetLength(1); j++)
                        outMatrix[i, j] = MatrixEin[i, j] > 0 ? 1 : 0;
            }
            return outMatrix;
        }
        /// <summary>
        /// 矩阵取对数
        /// </summary>
        /// <param name="MatrixEin"></param>
        /// <param name="basic"></param>
        /// <returns></returns>
        public static double[,] MatrixLog(double[,] MatrixEin,double basic=2)
        {
    
    
            int row = MatrixEin.GetLength(0);
            int col = MatrixEin.GetLength(1);
            double[,] outMatrix = new double[row, col];
            for(int i = 0; i < row; i++)
            {
    
    
                for(int j = 0; j < col; j++)
                {
    
    
                    outMatrix[i, j] =
                        basic == 2 ? Math.Log2(MatrixEin[i, j]) :
                        basic == 10 ? Math.Log10(MatrixEin[i, j]) :
                        Math.Log(MatrixEin[i, j]);
                }
            }
            return outMatrix;
        }
    }
}