using Accord.Math.Optimization;
using MathNet.Numerics;
using MathNet.Numerics.Interpolation;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearRegression;
using MathNet.Numerics.Statistics;

namespace IStation.Test
{

    /// <summary>
    /// 泵曲线数据融合校正器（用于合并样条曲线和实测数据并进行优化）
    /// </summary>
    public class PumpCurveDataFusionCorrectorHelper
    {
        # region 配置参数
        // 异常值检测阈值（基于标准正态分布，2.5对应98.76%置信区间）
        public const double ZScoreThreshold = 2.5;
        // 过渡区域宽度（单位与输入数据相同，控制曲线衔接平滑度）
        public const double TransitionWidth = 500;
        // 多项式拟合阶数（三次平衡拟合能力与过拟合风险）
        public const int PolynomialDegree = 3;
        // 优化器最大迭代次数（Cobyla算法收敛阈值）
        private const int OptimizationMaxIterations = 1000;
        # endregion

        /// <summary>
        /// 主校正方法（入口函数）
        /// </summary>
        /// <param name="splineX">样条曲线X坐标</param>
        /// <param name="splineY">样条曲线Y坐标</param>
        /// <param name="measuredXAll">原始实测X坐标</param>
        /// <param name="measuredYAll">原始实测Y坐标</param>
        /// <returns>
        /// mergedX: 融合后的X坐标
        /// mergedY: 初步融合的Y坐标
        /// optimizedX: 优化后的X坐标采样
        /// optimizedY: 优化后的Y坐标结果
        /// </returns>
        public (double[] mergedX, double[] mergedY, double[] optimizedX, double[] optimizedY) Corrent(
            double[] splineX, double[] splineY,
            double[] measuredXAll, double[] measuredYAll)
        {
            # region 步骤1：稳健回归与异常值过滤（最小二乘法）
            // 执行简单线性回归获取基准线
            (double intercept, double slope) = SimpleRegression.Fit(measuredXAll, measuredYAll);

            // 计算预测值和残差
            double[] predictedY = measuredXAll.Select(x => slope * x + intercept).ToArray();
            double[] residuals = measuredYAll.Zip(predictedY, (m, p) => m - p).ToArray();

            // 基于Z-Score的异常值过滤
            var zScoreCalculator = new ZScore(residuals);
            bool[] validMask = zScoreCalculator.Scores
                .Select(score => Math.Abs(score) < ZScoreThreshold).ToArray();

            // 过滤得到有效实测数据
            var validData = measuredXAll.Zip(measuredYAll, validMask)
                .Where(t => t.Third)
                .Select(t => (t.First, t.Second)).ToList();
            double[] measuredXValid = validData.Select(t => t.First).ToArray();
            double[] measuredYValid = validData.Select(t => t.Second).ToArray();
            # endregion

            if (measuredXValid.Length < 5)
            {
                Console.WriteLine("异常值过滤后数据过小");
                return default;
            }

            return ProcessData(splineX, splineY, measuredXValid, measuredYValid,
                             PolynomialDegree, TransitionWidth);
        }

        /// <summary>
        /// 核心数据处理流程（包含数据融合、过渡处理和优化）
        /// </summary>
        private (double[] mergedX, double[] mergedY, double[] optimizedX, double[] optimizedY) ProcessData(
            double[] splineX, double[] splineY,
            double[] measuredXValid, double[] measuredYValid,
            int polyDegree, double transitionWidth)
        {
            # region 步骤2：多项式拟合与初步融合
            // 对有效实测数据进行多项式拟合
            var polyFunc = Polynomial.Fit(measuredXValid, measuredYValid, polyDegree);
            double[] polyCoeffs = polyFunc.Coefficients.Reverse().ToArray(); // 注意系数顺序转换

            // 生成融合后的X坐标序列（合并样条和实测范围）
            double[] mergedX = GenerateMergedX(splineX, measuredXValid, 200);

            // 创建样条插值器并生成初始融合曲线
            var splineInterpolator = LinearSpline.InterpolateSorted(splineX, splineY);
            //double[] mergedY = mergedX.Select(x =>
            //    x < splineX[0] ? splineY[0] :         // 左外推
            //    x > splineX.Last() ? splineY.Last() :  // 右外推
            //    splineInterpolator.Interpolate(x)      // 插值区域
            //).ToArray();

            double[] mergedY = mergedX.Select(x =>
             //x < splineX[0] ? splineY[0] :         // 左外推
             //x > splineX.Last() ? splineY.Last() :  // 右外推
             splineInterpolator.Interpolate(x)      // 插值区域
         ).ToArray();
            # endregion

            # region 步骤3：核心区域修正
            // 在实测数据范围内应用多项式拟合结果
            double minX = measuredXValid.Min(), maxX = measuredXValid.Max();
            bool[] coreMask = mergedX.Select(x => x >= minX && x <= maxX).ToArray();
            for (int i = 0; i < mergedX.Length; i++)
            {
                if (coreMask[i])
                    mergedY[i] = EvaluatePolynomial(polyCoeffs, mergedX[i]);
            }
            # endregion

            # region 步骤4：过渡区域处理
            // 对左右过渡区域进行平滑处理
            var splinePoly = Polynomial.Fit(splineX, splineY, polyDegree);
            double[] splineCoeffs = splinePoly.Coefficients.Reverse().ToArray();

            // 左过渡区处理（样条->多项式）
            ApplyTransition(
                (minX - transitionWidth, minX),
                splineCoeffs,
                polyCoeffs,
                mergedX,
                ref mergedY
            );

            // 右过渡区处理（多项式->样条）
            ApplyTransition(
                (maxX, maxX + transitionWidth),
                polyCoeffs,
                splineCoeffs,
                mergedX,
                ref mergedY
            );
            # endregion

            # region 步骤5：带约束优化
            // 数据标准化处理
            double xMean = mergedX.Average();
            double xStd = ArrayStatistics.PopulationStandardDeviation(mergedX);
            double[] xNorm = mergedX.Select(x => (x - xMean) / xStd).ToArray();

            // 构建三次多项式优化问题
            var X = Matrix<double>.Build.DenseOfRowArrays(
                xNorm.Select(x => new[] { Math.Pow(x, 3), x * x, x, 1.0 }).ToArray()
            );
            var initialParams = MultipleRegression.QR(X, Vector<double>.Build.DenseOfArray(mergedY));

            // 设置单调性约束（导数非负）
            var constraints = GenerateLinspace(xNorm.Min(),  xNorm.Max(), 50)
                .Select(x => new NonlinearConstraint(
                    numberOfVariables: 4,
                    function: p => -Derivative(p, x), // 导数>=0
                    shouldBe: ConstraintType.GreaterThanOrEqualTo,
                    value: 0
                )).ToList();

            // 执行COBYLA优化
            var optimizer = new Cobyla(
                function: BuildObjectiveFunction(xNorm, mergedY),
                constraints: constraints.ToArray()
            )
            { MaxIterations = OptimizationMaxIterations };

            bool success = optimizer.Minimize(initialParams.ToArray());
            double[] optimizedParams = optimizer.Solution;
            # endregion

            // 生成最终优化曲线
            double[] optimizedX = GenerateLinspace(
                Math.Min(splineX.Min(), minX),
                Math.Max(splineX.Max(), maxX),
                300
            );
            double[] optimizedY = optimizedX
                .Select(x => CubicModel(optimizedParams, (x - xMean) / xStd))
                .ToArray();

            return (mergedX, mergedY, optimizedX, optimizedY);
        }

        # region 工具方法

        /// <summary>
        /// 三次多项式模型（降幂排列）
        /// </summary>
        private static double CubicModel(double[] param, double x) =>
            param[0] * x * x * x + param[1] * x * x + param[2] * x + param[3];

        /// <summary>
        /// 多项式导数（用于单调性约束）
        /// </summary>
        private static double Derivative(double[] param, double x) =>
            3 * param[0] * x * x + 2 * param[1] * x + param[2];

        /// <summary>
        /// 构建目标函数（均方误差）
        /// </summary>
        private NonlinearObjectiveFunction BuildObjectiveFunction(double[] xNorm, double[] y) =>
            new NonlinearObjectiveFunction(4, p =>
                xNorm.Select((x, i) => Math.Pow(CubicModel(p, x) - y[i], 2)).Average()
            );

        /// <summary>
        /// 多项式求值（霍纳法则）
        /// </summary>
        private static double EvaluatePolynomial(double[] coeffs, double x) =>
            coeffs.Aggregate((sum, c) => sum * x + c);

        /// <summary>
        /// 生成等间距采样点
        /// </summary>
        private static double[] GenerateLinspace(double start, double end, int numPoints) =>
            Enumerable.Range(0, numPoints)
                .Select(i => start + i * (end - start) / (numPoints - 1))
                .ToArray();
        # endregion

        # region 过渡处理
        /// <summary>
        /// 应用过渡混合（在两个多项式之间渐变）
        /// </summary>
        private static void ApplyTransition(
            (double start, double end) range,
            double[] baseCoeffs,
            double[] targetCoeffs,
            double[] xValues,
            ref double[] yValues)
        {
            for (int i = 0; i < xValues.Length; i++)
            {
                double x = xValues[i];
                if (x < range.start || x > range.end) continue;

                double weight = (x - range.start) / (range.end - range.start);
                double baseVal = EvaluatePolynomial(baseCoeffs, x);
                double targetVal = EvaluatePolynomial(targetCoeffs, x);
                yValues[i] = (1 - weight) * baseVal + weight * targetVal;
            }
        }

        /// <summary>
        /// 生成合并的 x 值
        /// </summary> 
        public static double[] GenerateMergedX(double[] splineX, double[] measuredXValid, int numPoints)
        {
            // 生成等间距的点
            double[] linspace = GenerateLinspace(measuredXValid.Min(), measuredXValid.Max(), numPoints);
            // 拼接并去重排序
            var concat = splineX.Concat(linspace).ToArray();
            var merged = concat.Distinct().OrderBy(x => x).ToArray();
            return merged;
        }

        # endregion

        # region 辅助类

        /// <summary>
        /// Z分数计算器（用于异常值检测）
        /// </summary>
        private class ZScore
        {
            public readonly double[] Scores;

            public ZScore(double[] data)
            {
                double mean = data.Average();
                double std = data.StandardDeviation();
                Scores = data.Select(x => (x - mean) / std).ToArray();
            }
        }

        # endregion
    }
}