using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra.Complex;
using MathNet.Numerics.LinearRegression;

namespace GeologyToolkit
{
    internal class KrigingMethod
    {
        /// <summary>
        ///     求已知点的距离矩阵
        /// </summary>
        /// <param name="points"></param>
        /// <returns></returns>
        public DenseMatrix DistanceDenseMatrix(List<XYPoint> points)
        {
            int n = points.Count;

            var matrix = new DenseMatrix(n, n);

            Complex[,] dis = new Complex[n, n];
            for (int i = 0; i < n; i++) //行
            {
                for (int j = 0; j < n; j++) //列
                {
                    Complex distance = points[i].DistanceTo(points[j]);

                    dis.SetValue(distance, i, j);
                }
            }

            matrix = DenseMatrix.OfArray(dis);

            return matrix;
        }

        /// <summary>
        ///     求属性值的半方差矩阵
        /// </summary>
        /// <param name="points"></param>
        public DenseMatrix ElevationVarD(List<XYPoint> points)
        {
            int n = points.Count;
            Complex[,] dis = new Complex[points.Count, points.Count];
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    //属性半方差
                    Complex var_D = Math.Pow(points[i].Elevation - points[j].Elevation, 2) / 2;
                    dis.SetValue(var_D, i, j);
                }
            }

            return DenseMatrix.OfArray(dis);
        }

        /// <summary>
        ///     分组后，每组距离的平均值-半方差值平均值（用于拟合半方差函数）
        /// </summary>
        /// <param name="pairs"></param>
        /// <param name="groupCount"></param>
        /// <returns></returns>
        public List<Pair> GetAveragePair(List<Pair> pairs, int groupCount)
        {
            var orderedData = pairs.OrderBy(x => x.Distance.Real);
            var max = orderedData.Last().Distance;
            var min = orderedData.First().Distance;
            double step = (max.Real - min.Real + 0.0001) / groupCount;
            //分组
            List<List<Pair>> ds = new List<List<Pair>>();
            for (int i = 0; i < groupCount; i++)
            {
                List<Pair> d = new List<Pair>();
                foreach (var pair in pairs)
                {
                    if (pair.Distance.Real >= step * i + min.Real && pair.Distance.Real <= (i + 1) * step + min.Real)
                    {
                        d.Add(pair);
                    }
                }

                ds.Add(d);
                //IEnumerable<Pair> d =
                //    orderedData.TakeWhile(n =>
                //        n.Distance.Real > step * i + min.Real && n.Distance.Real <= (i + 1) * step + min.Real);
                //ds.Add(d);
            }

            List<Pair> averagePairs = new List<Pair>();
            foreach (var d in ds)
            {
                Complex totalDis = 0.0;
                Complex totalD = 0.0;
                foreach (var complex in d)
                {
                    totalDis += complex.Distance;
                    totalD += complex.VarD;
                }

                var pair = new Pair
                {
                    Distance = totalDis / d.Count(),
                    VarD = totalD / d.Count()
                };
                averagePairs.Add(pair);
            }

            return averagePairs;
        }

        public List<Pair> GroupyDistances(DenseMatrix distances, DenseMatrix elevationD)
        {
            List<Pair> groups = new List<Pair>();
            for (int i = 0; i < distances.RowCount; i++)
            {
                for (int j = 0; j < elevationD.RowCount; j++)
                {
                    var r = new Pair
                    {
                        Distance = distances[i, j],
                        VarD = elevationD[i, j]
                    };
                    groups.Add(r);
                }
            }

            return groups;
        }

        public double MainMethod(List<XYPoint> points, XYPoint pointToInvestigate, int groupCount = 10)
        {
            var distancesMatrix = DistanceDenseMatrix(points);
            var elevationsVarD = ElevationVarD(points); //高程属性值的半方差
            List<Pair> pairs = new List<Pair>();
            for (int i = 0; i < distancesMatrix.RowCount; i++)
            {
                for (int j = 0; j < elevationsVarD.RowCount; j++)
                {
                    for (int k = 0; k < distancesMatrix.ColumnCount; k++)
                    {
                        Pair pair = new Pair
                        {
                            Distance = distancesMatrix.At(i, k),
                            VarD = elevationsVarD.At(j, k)
                        };
                        pairs.Add(pair);
                    }
                }
            }

            var fitData = GetAveragePair(pairs, groupCount);
            //通过平均距离-半方差，得到拟合函数（距离-半方差关系函数）
            double[] x = new double[fitData.Count];
            double[] y = new double[fitData.Count];
            for (int i = 0; i < x.Length; i++)
            {
                x[i] = fitData[i].Distance.Real;
                y[i] = fitData[i].VarD.Real;
            }

            //指数模型
            var formula_P = Fit.ExponentialFunc(x, y);
            //球形模型
            var formula_S = Fit.PolynomialFunc(x, y, 3, DirectRegressionMethod.NormalEquations);


            //求未知点的距离列矩阵
            var unknownDistancesMatrix = UnknownPointDistances(pointToInvestigate, points);
            //代入拟合函数，求未知点的半方差矩阵
            Complex[] rows = new Complex[unknownDistancesMatrix.Values.Length];
            for (var i = 0; i < unknownDistancesMatrix.Values.Length; i++)
            {
                var d = unknownDistancesMatrix.Values[i];
                var distance = d.Real;
                //半方差值
                var ps = formula_S(distance);
                rows[i] = ps;
            }

            var unknownVarDMatrix = DenseMatrix.OfColumnArrays(rows);
            //得到权重（矩阵的逆乘以未知点的距离半方差矩阵)
            var inverse = elevationsVarD.Inverse();
            var wi = inverse * unknownVarDMatrix;
            double n = default;
            for (int i = 0; i < wi.RowCount; i++)
            {
                double na = wi[i, 0].Real;
                n += na;
            }

            //权重乘以对应属性，加权得到
            double prop = default;
            for (int i = 0; i < points.Count; i++)
            {
                XYPoint point = points[i];
                prop += point.Elevation * wi[i, 0].Real;
            }

            return prop;
            //var data = matrix.Values;
            //var orderedData = data.OrderBy(x => x);
            //var max = data.Max();
            //double step = (max.Real + 10) / groupCount;
            //List<IEnumerable<Complex>> ds = new List<IEnumerable<Complex>>();
            //for (int i = 0; i < 10; i++)
            //{
            //    IEnumerable<Complex> d = orderedData.TakeWhile(n => n.Real > step * i && n.Real < (i + 1) * step);
            //    ds.Add(d);
            //}
            //foreach (var d in ds)
            //{
            //    Complex total = 0.0;
            //    foreach (var complex in d)
            //    {
            //        total += complex;
            //    }
            //    var average = total / d.Count();
            //}
            //matrix.Inverse();
        }

        /// <summary>
        ///     求未知点距离列矩阵(代入拟合半方差函数，可求得半方差)
        /// </summary>
        /// <param name="pointToInvestigate"></param>
        /// <param name="points"></param>
        public DenseMatrix UnknownPointDistances(XYPoint pointToInvestigate, List<XYPoint> points)
        {
            //DenseMatrix matrix = new DenseMatrix(points.Count, 1);
            Complex[] rows = new Complex[points.Count];
            for (int i = 0; i < points.Count; i++)
            {
                Complex complex = pointToInvestigate.DistanceTo(points[i]);
                rows[i] = complex;
            }

            return DenseMatrix.OfColumnArrays(rows);
        }

        private double[,] Deinv(double[,] dMatrix, int Level)
        {
            double dMatrixValue = MatrixValue(dMatrix, Level);
            if (dMatrixValue == 0) return null;

            double[,] dReverseMatrix = new double[Level, 2 * Level];
            double x, c;
            // Init Reverse matrix
            for (int i = 0; i < Level; i++)
            {
                for (int j = 0; j < 2 * Level; j++)
                {
                    if (j < Level)
                        dReverseMatrix[i, j] = dMatrix[i, j];
                    else
                        dReverseMatrix[i, j] = 0;
                }

                dReverseMatrix[i, Level + i] = 1;
            }

            for (int i = 0, j = 0; i < Level && j < Level; i++, j++)
            {
                if (dReverseMatrix[i, j] == 0)
                {
                    int m = i;
                    for (; dMatrix[m, j] == 0; m++) ;
                    if (m == Level)
                        return null;
                    // Add i-row with m-row
                    for (int n = j; n < 2 * Level; n++)
                        dReverseMatrix[i, n] += dReverseMatrix[m, n];
                }

                // Format the i-row with "1" start
                x = dReverseMatrix[i, j];
                if (x != 1)
                {
                    for (int n = j; n < 2 * Level; n++)
                        if (dReverseMatrix[i, n] != 0)
                            dReverseMatrix[i, n] /= x;
                }

                // Set 0 to the current column in the rows after current row
                for (int s = Level - 1; s > i; s--)
                {
                    x = dReverseMatrix[s, j];
                    for (int t = j; t < 2 * Level; t++)
                        dReverseMatrix[s, t] -= dReverseMatrix[i, t] * x;
                }
            }

            // Format the first matrix into unit-matrix
            for (int i = Level - 2; i >= 0; i--)
            {
                for (int j = i + 1; j < Level; j++)
                    if (dReverseMatrix[i, j] != 0)
                    {
                        c = dReverseMatrix[i, j];
                        for (int n = j; n < 2 * Level; n++)
                            dReverseMatrix[i, n] -= c * dReverseMatrix[j, n];
                    }
            }

            double[,] dReturn = new double[Level, Level];
            for (int i = 0; i < Level; i++)
            for (int j = 0; j < Level; j++)
                dReturn[i, j] = dReverseMatrix[i, j + Level];
            return dReturn;
        }
        //private List<double> pointToInvestigate(XYPoint pointToInvestigate, List<XYPoint> points)
        //{
        //    List<double> dis = new List<double>();
        //    int n = points.Count;
        //    for (int i = 0; i < n; i++)
        //    {

        //        var distance = pointToInvestigate.DistanceTo(points[i]);
        //        dis.Add(distance);
        //    }
        //    return dis;
        //}
        /// <summary>
        ///     求已知点的距离矩阵数组（无依赖）
        /// </summary>
        /// <param name="points"></param>
        /// <returns></returns>
        private double[,] DistanceMatrix(List<XYPoint> points)
        {
            double[,] dis = new double[points.Count + 1, points.Count];
            int n = points.Count;
            for (int i = 0; i <= n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    double distance = points[i].DistanceTo(points[j]);

                    dis.SetValue(distance, i, j);
                }
            }

            return dis;
        }

        private void FitExponential(double[] x, double[] y)
        {
            var f = Fit.ExponentialFunc(x, y);
            double d = 15.1;
            var height = f.Invoke(d);
        }

        /// <summary>
        ///     高斯模型
        /// </summary>
        /// <param name="h">采样点距离</param>
        /// <param name="a">变程，与起始采样点距离，</param>
        private void FitMethod(double h, double a)
        {
            if (h == 0)
            {
            }
            else if (h > 0)
            {
            }
        }

        private double MatrixValue(double[,] MatrixList, int Level)
        {
            double[,] dMatrix = new double[Level, Level];
            for (int i = 0; i < Level; i++)
            for (int j = 0; j < Level; j++)
                dMatrix[i, j] = MatrixList[i, j];
            double c, x;
            int k = 1;
            for (int i = 0, j = 0; i < Level && j < Level; i++, j++)
            {
                if (dMatrix[i, j] == 0)
                {
                    int m = i;
                    for (; dMatrix[m, j] == 0; m++) ;
                    if (m == Level)
                        return 0;
                    // Row change between i-row and m-row
                    for (int n = j; n < Level; n++)
                    {
                        c = dMatrix[i, n];
                        dMatrix[i, n] = dMatrix[m, n];
                        dMatrix[m, n] = c;
                    }

                    // Change value pre-value
                    k *= -1;
                }

                // Set 0 to the current column in the rows after current row
                for (int s = Level - 1; s > i; s--)
                {
                    x = dMatrix[s, j];
                    for (int t = j; t < Level; t++)
                        dMatrix[s, t] -= dMatrix[i, t] * (x / dMatrix[i, j]);
                }
            }

            double sn = 1;
            for (int i = 0; i < Level; i++)
            {
                if (dMatrix[i, i] != 0)
                    sn *= dMatrix[i, i];
                else
                    return 0;
            }

            return k * sn;
        }

        /// <summary>
        ///     普通克里金
        /// </summary>
        /// <param name="dpx">数据点坐标x</param>
        /// <param name="dpy">数据点坐标y</param>
        /// <param name="dpz">数据点坐标z</param>
        /// <param name="ndp">数据点个数</param>
        /// <param name="ipx">待插入点的坐标x</param>
        /// <param name="ipy">待插入点的坐标y</param>
        /// <param name="ipz">待插入点的坐标z</param>
        /// <param name="nip">待插入点的个数</param>
        private void OrdinaryKriging(double[] dpx, double[] dpy, double[] dpz, int ndp, double[] ipx, double[] ipy, double[] ipz, int nip)
        {
            int LM = 10; //距离分组时的组数
            double DR1, DZ;
            double[] a = new double[3]; //存储球状模型三个系数
            double[] p, A, B, C; //矩阵
            int[] LD = new int[LM];
            double[] AR = new double[LM];
            double[] AZ = new double[LM];
            double[] DR = new double[LM + 1];

            p = new double[(ndp - 1) * ndp / 2]; //距离矩阵
            A = new double[Convert.ToInt16(Math.Pow(ndp + 1, 2))]; //系数矩阵
            for (int i = 0; i < Math.Pow(ndp + 1, 2); i++)
            {
                A[i] = 0;
            }

            int numOfDis = 0;
            for (int i = 0; i < ndp - 1; i++)
            {
                for (int j = 0; j < ndp; j++)
                {
                    p[numOfDis] = Math.Sqrt(Math.Pow(dpx[i] - dpx[j], 2) + Math.Pow(dpx[i] - dpx[j], 2));
                    A[i * (ndp + 1) + j] = p[numOfDis]; //记录距离矩阵上半部分
                    A[j * (ndp + 1) + i] = p[numOfDis]; //记录距离矩阵下半部分
                    numOfDis++;
                }
            }

            //堆排序，由小到大
            p = p.OrderBy(x => x).ToArray();
            DR[0] = 0.0;
            for (int i = 1; i < LM + 1; i++)
            {
                var k = i * numOfDis / LM - 1;
                DR[i] = p[k];
            }

            for (int i = 1; i < LM; i++)
            {
                LD[i] = 0;
                AR[i] = 0;
                AZ[i] = 0;
            }

            for (int i = 0; i < ndp - 1; i++)
            {
                for (int j = 0; j < ndp; j++)
                {
                    DR1 = Math.Sqrt(Math.Pow(dpx[i] - dpx[j], 2) + Math.Pow(dpx[i] - dpx[j], 2)); //距离
                    DZ = Math.Pow(dpz[i] - dpx[j], 2); //平方增量
                    for (int k = 0; k <= LM; k++)
                    {
                        if (DR1 > DR[k - 1] && DR1 <= DR[k])
                        {
                            LD[k - 1] += 1; //点个数
                            AR[k - 1] += DR1; //每组距离和
                            AZ[k - 1] += DZ; //每组平方增量和
                        }
                    }
                }
            }

            //二元线性回归
            //L11*X1+L12*X2=L10
            //L21*X1+L22*X2=L20
            double X1 = 0.0;
            double X2 = 0.0;
            double Y = 0.0;
            double L11 = 0.0;
            double L12 = 0.0;
            double L21 = 0.0;
            double L22 = 0.0;
            double L10 = 0.0;
            double L20 = 0.0;
            for (int i = 0; i < LM; i++)
            {
                AR[i] = AR[i] / LD[i]; //距离平均值
                AZ[i] = AZ[i] / LD[i] / 2; //平方增量和
                X1 += AR[i] / LM;
                X2 += Math.Pow(AR[i], 3);
                Y += AZ[i] / LM;
            }

            for (int i = 0; i < LM; i++)
            {
                L10 += (AR[i] - X1) * (AZ[i] - Y);
                L11 += Math.Pow(AR[i] - X1, 2);
                L12 += (AR[i] - X1) * (Math.Pow(AR[i], 3) - X2);
                L12 += (Math.Pow(AR[i], 3) - X2) * (AZ[i] - Y);
                L21 += (Math.Pow(AR[i], 3) - X2) * (AR[i] - X1);
                L22 += (Math.Pow(AR[i], 3) - X2) * (Math.Pow(AR[i], 3) - X2);
                //球状模型系数
                double B1 = (L10 * L22 - L20 * L12) / (L11 * L22 - L12 * L12);
                double B2 = (L20 * L11 - L10 * L21) / (L11 * L22 - L12 * L12);
                double B0 = Y - B1 * X1 - B2 * X2;
                if (B0 > 0)
                {
                    a[0] = B0;
                    a[1] = Math.Sqrt(-B1 / 3 / B2);
                    a[2] = 2 * B1 / 3 * Math.Sqrt(-B1 / 3 / B2);
                }
            }

            for (int i = 0; i < ndp; i++)
            {
                for (int j = 0; j < ndp; j++)
                {
                    if (i == j)
                    {
                        a[i * (ndp + 1) + j] = 0.0;
                    }
                    else
                    {
                        if (A[i * (ndp + 1) + j] < a[1]) //小于变程
                        {
                            A[i * (ndp + 1) + j] = a[0] + a[2] * (1.5 * A[i * (ndp + 1) + j] / a[1]) - 0.5 * Math.Pow(A[i * (ndp + 1) + j] / a[1], 3);
                        }
                        else
                        {
                            A[i * (ndp + 1) + j] = a[0] + a[2]; //大于等于变程
                        }
                    }
                }
            }

            for (int i = 0; i < ndp; i++) //对每个数据点加入值为1的项
            {
                A[(i + 1) * (ndp + 1) - 1] = 1.0;
                A[ndp * (ndp + 1) + i] = 1.0;
            }

            //系数矩阵求逆
            //Deinv(A, ndp + 1);
            B = new double[ndp + 1]; //存放矩阵方程组等号右边的列向量
            C = new double[ndp + 1]; //存放结果
            //待插入点循环
            for (int i = 0; i < nip; i++)
            {
                for (int j = 0; j < ndp; j++)
                {
                    DR1 = Math.Sqrt(Math.Pow(ipx[i] - dpx[j], 2) + Math.Pow(ipy[i] - dpy[j], 2));
                    if (DR1 < 0.00000000001)
                    {
                        B[j] = 0.0;
                    }
                    else
                    {
                        if (DR1 < a[1])
                        {
                            B[j] = a[0] + a[2] * (1.5 * DR1 / a[1] - 0.5 * Math.Pow(DR1 / a[1], 3));
                        }
                        else
                        {
                            B[j] = a[0] + a[2];
                        }
                    }
                }

                B[ndp] = 1.0; //加入值为1的项
                for (int k = 0; k < ndp + 1; k++) //插值方程的解
                {
                    C[k] = 0.0;
                    for (int j = 0; j < ndp + 1; j++)
                    {
                        C[k] += A[k * (ndp + 1) + j] * B[j];
                    }
                }

                ipz[i] = 0.0;
                for (int j = 0; j < ndp; j++)
                {
                    ipz[i] += C[j] * dpz[j]; //普通克里金估值
                }
            }
        }
    }

    public class XYPoint
    {
        public XYPoint(double x, double y)
        {
            X = x;
            Y = y;
        }

        public XYPoint()
        {
            X = 0.0;
            Y = 0.0;
        }

        public double Elevation { get; set; }

        public double X { get; set; }

        public double Y { get; set; }

        /// <summary>
        ///     水平距离
        /// </summary>
        /// <param name="point"></param>
        /// <returns></returns>
        public double DistanceTo(XYPoint point)
        {
            double d = Math.Pow(X - point.X, 2) + Math.Pow(Y - point.Y, 2);
            return Math.Sqrt(d);
        }
    }

    public class Pair
    {
        /// <summary>
        ///     距离
        /// </summary>
        public Complex Distance { get; set; }

        /// <summary>
        ///     半方差
        /// </summary>
        public Complex VarD { get; set; }
    }
}