using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra.Double;
using NPOI.SS.Formula.Functions;
using System;
using System.Collections.Generic;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Hydro.CommonBase
{
    public class CurveFitHelper
    {


        public DenseVector coefficients;
        int degree;
        double xMax;
        double xMin;
        List<PointF> data;
        bool isFitted = false;

        public CurveFitHelper(List<PointF> data, int degree = 2)
        {
            this.degree = degree;
            this.data = data;

            double[] xlist = data.Select(p => (double)p.X).ToArray();
            coefficients = Fit.Polynomial(xlist, data.Select(p => (double)p.Y).ToArray(), 2);
            isFitted = true;
            var list = xlist.ToList();
            xMax = list.Max();
            xMin = list.Min();
        }

        public double Evaluate(double x)
        {
            while (!isFitted) ;
            double y = 0;
            for (int i = 0; i < coefficients.Count; i++)
            {
                y += (float)(coefficients[i] * Math.Pow(x, i));
            }
            return y;
        }


        public List<PointF> GetFitCurve(int count)
        {
            List<PointF> curve = new List<PointF>();
            for (double x = xMin; x < xMax; x += (xMax - xMin) / count)
            {
                double y = Evaluate(x);
                curve.Add(new PointF((float)x, (float)y));
            }
            return curve;

        }

        public List<List<PointF>> GetConfidenceCurve(int count)
        {
            List<PointF> up_list = new List<PointF>();
            List<PointF> down_list = new List<PointF>();
            for (double x = xMin; x < xMax; x += (xMax - xMin) / count)
            {
                var list = data.FindAll(p => x <= p.X && p.X < x + (xMax - xMin) / count).Select(p => (double)p.Y).ToList();
                if (list.Count < 3) continue;
                var Bounds = ConfidenceRangeHelper.GetConfidenceRange(list, 95);
                var x_avg = x;
                if (Bounds[0] == double.NaN || Bounds[1] == double.NaN) continue;
                down_list.Add(new PointF((float)x_avg, (float)Bounds[0]));
                up_list.Add(new PointF((float)x_avg, (float)Bounds[1]));
            }

            //CurveFitHelper h0 = new CurveFitHelper(down_list);
            //var curve0 = h0.GetFitCurve(count);
            //CurveFitHelper h1 = new CurveFitHelper(up_list);
            //var curve1 = h1.GetFitCurve(count);
            return new List<List<PointF>> { down_list, up_list };
        }
    }

    public static class ConfidenceRangeHelper
    {
        public static double[] GetConfidenceRange(List<double> data, double confidenceLevel = 95)
        {
            // 使用循环计算数据的最小值和最大值
            double minValue = double.MaxValue;
            double maxValue = double.MinValue;
            foreach (double value in data)
            {
                if (value < minValue)
                {
                    minValue = value;
                }
                if (value > maxValue)
                {
                    maxValue = value;
                }
            }

            // 使用循环计算样本均值
            double sum = 0;
            foreach (double value in data)
            {
                sum += value;
            }
            double mean = sum / data.Count;

            // 使用循环计算标准差
            double sumOfSquares = 0;
            foreach (double value in data)
            {
                sumOfSquares += Math.Pow(value - mean, 2);
            }
            double stdDev = Math.Sqrt(sumOfSquares / (data.Count - 1));

            double zValue = 0;
            double alpha = 1 - confidenceLevel / 100.0;

            double zValueOneTail = GetZValueOneTail(alpha / 2);
            double marginOfError = zValueOneTail * stdDev;
            double lowerBound = mean - marginOfError;
            double upperBound = mean + marginOfError;

            return new double[] { lowerBound, upperBound };
        }
        static double[,] zTable = {
      {0.0000, 0.0039, 0.0078, 0.0117, 0.0156, 0.0195, 0.0234, 0.0274, 0.0314, 0.0353},
      {0.0392, 0.0432, 0.0471, 0.0510, 0.0550, 0.0589, 0.0628, 0.0668, 0.0708, 0.0749},
      {0.0789, 0.0830, 0.0871, 0.0912, 0.0953, 0.0994, 0.1035, 0.1076, 0.1117, 0.1158},
      {0.1190, 0.1231, 0.1272, 0.1314, 0.1357, 0.1398, 0.1441, 0.1483, 0.1525, 0.1568},
      {0.1611, 0.1654, 0.1697, 0.1740, 0.1783, 0.1826, 0.1869, 0.1912, 0.1955, 0.1999},
      {0.2042, 0.2085, 0.2129, 0.2172, 0.2216, 0.2260, 0.2304, 0.2348, 0.2392, 0.2436},
      {0.2480, 0.2524, 0.2568, 0.2612, 0.2658, 0.2703, 0.2747, 0.2793, 0.2838, 0.2883},
      {0.2929, 0.2974, 0.3020, 0.3066, 0.3112, 0.3159, 0.3205, 0.3251, 0.3299, 0.3346},
      {0.3393, 0.3441, 0.3487, 0.3535, 0.3583, 0.3632, 0.3679, 0.3727, 0.3776, 0.3823},
      {0.3872, 0.3920, 0.3968, 0.4017, 0.4066, 0.4115, 0.4164, 0.4214, 0.4264, 0.4313}
    };
        static double GetZValueOneTail(double alpha)
        {
            if (alpha == 0.5)
            {
                return 0;
            }
            else if (alpha > 0.5)
            {
                // Use symmetry property of the normal distribution
                return -GetZValueOneTail(1 - alpha);
            }
            double[] alphaValues = { 0.1, 0.05, 0.025, 0.01, 0.005, 0.001 };
            double[] zValues = { 1.28, 1.645, 1.96, 2.33, 2.575, 3.09 };
            int n = alphaValues.Length;
            int i = 0;
            while (i < n && alpha < alphaValues[i])
            {
                i++;
            }

            if (i == 0)
            {
                return 0;
            }
            else if (i == n)
            {
                return zValues[n - 1];
            }

            double x0 = alphaValues[i - 1];
            double x1 = alphaValues[i];
            double y0 = zValues[i - 1];
            double y1 = zValues[i];

            double zValue = ((y1 - y0) / (x1 - x0)) * (alpha - x0) + y0;

            return zValue;
        }

    }
}