using System;
using System.Collections.Generic;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace ErrorAnalysis.Service
{
    public class Utility
    {
        /// <summary> 
        /// 单线性插值法
        /// </summary>
        /// <param name="value">插值</param>
        /// <param name="x0">起始值</param>
        /// <param name="y0">起始值对应值</param>
        /// <param name="x1">结束值</param>
        /// <param name="y1">结束值对应值</param>
        /// <returns>插值对应结果</returns>
        public static double Interpolate(double value, double x0, double y0, double x1, double y1)
        {
            // Calculate the interpolated value using linear interpolation formula
            double y = y0 + (value - x0) * (y1 - y0) / (x1 - x0);
            return y;
        }

        public static (double slope, double intercept) FitLine(PointF[] points)
        {
            // 计算各求和项
            int n = points.Length;
            double sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;

            foreach (var p in points)
            {
                sumX += p.X;
                sumY += p.Y;
                sumXY += p.X * p.Y;
                sumX2 += p.X * p.X;
            }

            // 计算斜率和截距
            double denominator = n * sumX2 - sumX * sumX;
            if (Math.Abs(denominator) < 1e-10) // 避免除以0
                throw new InvalidOperationException("点集过于集中，无法拟合直线");

            double slope = (n * sumXY - sumX * sumY) / denominator;
            double intercept = (sumY - slope * sumX) / n;

            return (slope, intercept);
        }

        public static double GetLineValue(double x, double slope, double intercept) => slope * x + intercept;

        public static (double a, double b, double c) FitParabola(PointF[] points)
        {
            if (points.Length < 3)
                throw new ArgumentException("至少需要3个点进行抛物线拟合");

            // 1. 计算各项求和值[1,6](@ref)
            int n = points.Length;
            double sx = 0, sy = 0, sx2 = 0, sx3 = 0, sx4 = 0, sxy = 0, sx2y = 0;

            foreach (var p in points)
            {
                double x = p.X;
                double y = p.Y;
                double x2 = x * x;
                double x3 = x2 * x;
                double x4 = x2 * x2;

                sx += x;
                sy += y;
                sx2 += x2;
                sx3 += x3;
                sx4 += x4;
                sxy += x * y;
                sx2y += x2 * y;
            }

            // 2. 构建正规方程矩阵[1,6](@ref)
            double[,] matrix = {
            { sx4, sx3, sx2, sx2y },
            { sx3, sx2, sx,  sxy  },
            { sx2, sx,  n,   sy   }
        };

            // 3. 使用高斯消元法求解
            double[] coefficients = GaussElimination(matrix);

            return (coefficients[0], coefficients[1], coefficients[2]);
        }

        // 高斯消元法实现（带部分主元选择）
        private static double[] GaussElimination(double[,] matrix)
        {
            int n = matrix.GetLength(0);
            double[] result = new double[n];

            // 前向消元
            for (int i = 0; i < n - 1; i++)
            {
                // 部分主元选择
                int maxRow = i;
                for (int k = i + 1; k < n; k++)
                {
                    if (Math.Abs(matrix[k, i]) > Math.Abs(matrix[maxRow, i]))
                        maxRow = k;
                }

                // 行交换
                if (maxRow != i)
                {
                    for (int j = 0; j <= n; j++)
                    {
                        double temp = matrix[i, j];
                        matrix[i, j] = matrix[maxRow, j];
                        matrix[maxRow, j] = temp;
                    }
                }

                // 消元过程
                for (int k = i + 1; k < n; k++)
                {
                    double factor = matrix[k, i] / matrix[i, i];
                    for (int j = i; j <= n; j++)
                    {
                        matrix[k, j] -= factor * matrix[i, j];
                    }
                }
            }

            // 回代求解
            for (int i = n - 1; i >= 0; i--)
            {
                result[i] = matrix[i, n];
                for (int j = i + 1; j < n; j++)
                {
                    result[i] -= matrix[i, j] * result[j];
                }
                result[i] /= matrix[i, i];
            }

            return result;
        }

        public static (double a, double b, double c) quadraticLine(double[] arrX, double[] arrY)
        {
            double[] coef = new double[3];
            double[,] AA = new double[3, 3];
            double[,] TAA = new double[3, 3];
            double[,] BAA = new double[3, 1];
            double[] B = new double[3];
            int m = arrX.Length;

            double[] arrXp = new double[m];
            double[] arrXc = new double[m];
            double[] arrXv = new double[m];
            double[] y = new double[m];

            for (int i = 0; i < m; i++)
            {
                arrXp[i] = arrX[i] * arrX[i];
            }
            for (int i = 0; i < m; i++)
            {
                arrXc[i] = arrX[i] * arrX[i] * arrX[i];
            }
            for (int i = 0; i < m; i++)
            {
                arrXv[i] = arrX[i] * arrX[i] * arrX[i] * arrX[i];
            }
            double a1 = mysum(arrX);
            double a2 = mysum(arrXp);
            double a3 = mysum(arrXc);
            double a4 = mysum(arrXv);

            AA[0, 0] = m;
            AA[0, 1] = a1;
            AA[0, 2] = a2;
            AA[1, 0] = a1;
            AA[1, 1] = a2;
            AA[1, 2] = a3;
            AA[2, 0] = a2;
            AA[2, 1] = a3;
            AA[2, 2] = a4;
            double[] xy = new double[m];
            double[] xxy = new double[m];
            for (int i = 0; i < m; i++)
            {
                xy[i] = arrX[i] * arrY[i];
            }
            for (int i = 0; i < m; i++)
            {
                xxy[i] = arrX[i] * arrX[i] * arrY[i];
            }
            double b1 = mysum(arrY);
            double b2 = mysum(xy);
            double b3 = mysum(xxy);
            B[0] = b1;
            B[1] = b2;
            B[2] = b3;
            TAA = ReverseMatrix(AA, 3);
            for (int i = 0; i < 3; i++)
            {
                for (int j = 0; j < 3; j++)
                {
                    coef[i] = coef[i] + TAA[i, j] * B[j];
                }
            }
            return (coef[2], coef[1], coef[0]);

        }

        public static double[,] ReverseMatrix(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;
                    else
                    {
                        // 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 static 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;
                    else
                    {
                        // 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;
        }

        public static double mysum(double[] XX)
        {
            double reout = 0;
            int mm = XX.Length;
            for (int i = 0; i < mm; i++)
            {
                reout = reout + XX[i];
            }
            return reout;
        }

        public static double GetParabolaValue(double x, double a, double b, double c)
        {
            return a * x * x + b * x + c;
        }
    }
}