using ConeCalorimeter.Models;
using Serilog;

namespace ConeCalorimeter.Services;

public sealed class HeatReleaseCorrectionService : IHeatReleaseCorrectionService
{
    private const double TotalHeatReleaseDivisor = 1000;
    private const double NumericalTolerance = 1e-9;
    private const double MaximumCorrectedHeatReleaseRate = 5000;

    private readonly JsonHeatReleaseResponseCalibrationStore _calibrationStore;

    public HeatReleaseCorrectionService(JsonHeatReleaseResponseCalibrationStore calibrationStore)
    {
        _calibrationStore = calibrationStore;
    }

    public HeatReleaseCorrectionResult Correct(IReadOnlyList<RealtimeDataRecord> rawRecords)
    {
        var orderedRecords = rawRecords
            .Where(record => record.TestSeconds >= 0)
            .OrderBy(record => record.TestSeconds)
            .ThenBy(record => record.Timestamp)
            .ToList();

        if (orderedRecords.Count < 3)
        {
            return NotApplied(orderedRecords, "有效 HRR 记录不足，无法生成正式校正结果。");
        }

        var calibration = _calibrationStore.Load();
        if (!TryValidateCalibration(calibration, out var calibrationError))
        {
            return NotApplied(orderedRecords, calibrationError);
        }

        var rawHrr = orderedRecords
            .Select(record => double.IsFinite(record.HeatReleaseRate) ? Math.Max(0, record.HeatReleaseRate) : 0)
            .ToArray();
        var weights = BuildIntegrationWeights(orderedRecords);
        if (!TryValidateSampleInterval(weights, calibration!.SampleIntervalSeconds, out var intervalError))
        {
            return NotApplied(orderedRecords, intervalError);
        }

        var rawFinalThr = orderedRecords
            .Where(record => double.IsFinite(record.TotalHeatRelease))
            .Select(record => record.TotalHeatRelease)
            .DefaultIfEmpty(double.NaN)
            .Last();
        var targetArea = double.IsFinite(rawFinalThr)
            ? rawFinalThr * TotalHeatReleaseDivisor
            : Dot(weights, rawHrr);
        if (targetArea <= NumericalTolerance)
        {
            return NotApplied(orderedRecords, "原始 HRR 积分为零，无法生成正式校正结果。");
        }

        var kernel = NormalizeKernel(calibration.ResponseKernel);
        double[] correctedHrr;
        try
        {
            correctedHrr = SolveConstrainedDeconvolution(rawHrr, weights, targetArea, kernel, calibration);
        }
        catch (InvalidOperationException ex)
        {
            return NotApplied(orderedRecords, ex.Message);
        }

        if (correctedHrr.Any(value => !double.IsFinite(value) || value < 0 || value > MaximumCorrectedHeatReleaseRate))
        {
            return NotApplied(
                orderedRecords,
                $"HRR 校正结果超出允许范围 0-{MaximumCorrectedHeatReleaseRate:F0} kW/m²。");
        }

        var correctedRecords = BuildCorrectedRecords(orderedRecords, correctedHrr, weights, calibration.Version);
        var correctedFinalThr = correctedRecords[^1].TotalHeatRelease;
        var areaErrorPercent = Math.Abs(correctedFinalThr - rawFinalThr) / rawFinalThr * 100;

        if (areaErrorPercent > calibration.MaximumAreaErrorPercent)
        {
            return NotApplied(
                orderedRecords,
                $"HRR 校正面积守恒误差 {areaErrorPercent:F6}% 超过允许值 {calibration.MaximumAreaErrorPercent:F6}%。");
        }

        var rawPeak = rawHrr.Max();
        var correctedPeak = correctedHrr.Max();
        Log.Information(
            "Heat release correction applied. CalibrationVersion={CalibrationVersion}, Records={RecordCount}, RawPeak={RawPeak}, CorrectedPeak={CorrectedPeak}, RawFinalTHR={RawFinalTHR}, CorrectedFinalTHR={CorrectedFinalTHR}, AreaErrorPercent={AreaErrorPercent}.",
            calibration.Version,
            correctedRecords.Count,
            rawPeak,
            correctedPeak,
            rawFinalThr,
            correctedFinalThr,
            areaErrorPercent);

        return new HeatReleaseCorrectionResult(
            true,
            $"正式 HRR 校正已完成，标定版本：{calibration.Version}",
            calibration.Version,
            correctedRecords,
            rawPeak,
            correctedPeak,
            rawFinalThr,
            correctedFinalThr,
            areaErrorPercent);
    }

    private static bool TryValidateCalibration(
        HeatReleaseResponseCalibration? calibration,
        out string error)
    {
        if (calibration is null)
        {
            error = "缺少热释放响应标定文件，无法生成正式校正结果。";
            return false;
        }

        if (string.IsNullOrWhiteSpace(calibration.Version)
            || string.IsNullOrWhiteSpace(calibration.Method)
            || calibration.CalibratedAtUtc == default)
        {
            error = "热释放响应标定缺少版本、方法或标定日期。";
            return false;
        }

        if (calibration.ResponseKernel.Count < 2
            || calibration.ResponseKernel.Any(value => !double.IsFinite(value) || value < 0)
            || calibration.ResponseKernel.Sum() <= NumericalTolerance)
        {
            error = "热释放响应标定核无效。";
            return false;
        }

        if (!double.IsFinite(calibration.SampleIntervalSeconds) || calibration.SampleIntervalSeconds <= 0
            || !double.IsFinite(calibration.RegularizationLambda) || calibration.RegularizationLambda < 0
            || !double.IsFinite(calibration.GradientStep) || calibration.GradientStep <= 0 || calibration.GradientStep > 1
            || calibration.Iterations < 1 || calibration.Iterations > 100000
            || calibration.ResponseKernel.Count > 600
            || !double.IsFinite(calibration.ConvergenceTolerance) || calibration.ConvergenceTolerance <= 0
            || !double.IsFinite(calibration.MaximumAreaErrorPercent) || calibration.MaximumAreaErrorPercent < 0
            || !calibration.ValidationErrorPercent.HasValue
            || !double.IsFinite(calibration.ValidationErrorPercent.Value)
            || calibration.ValidationErrorPercent.Value < 0)
        {
            error = "热释放响应标定参数无效。";
            return false;
        }

        error = string.Empty;
        return true;
    }

    private static double[] BuildIntegrationWeights(IReadOnlyList<RealtimeDataRecord> records)
    {
        var weights = new double[records.Count];
        for (var i = 1; i < records.Count; i++)
        {
            weights[i] = records[i].TestSeconds - records[i - 1].TestSeconds;
        }

        return weights;
    }

    private static bool TryValidateSampleInterval(
        IReadOnlyList<double> weights,
        double expectedInterval,
        out string error)
    {
        var validIntervals = weights.Skip(1).Where(value => value > 0).ToArray();
        if (validIntervals.Length != weights.Count - 1)
        {
            error = "HRR 记录时间不连续或存在重复，无法生成正式校正结果。";
            return false;
        }

        var tolerance = Math.Max(0.05, expectedInterval * 0.1);
        if (validIntervals.Any(value => Math.Abs(value - expectedInterval) > tolerance))
        {
            error = $"HRR 采样周期与标定周期 {expectedInterval:F3}s 不一致。";
            return false;
        }

        error = string.Empty;
        return true;
    }

    private static double[] NormalizeKernel(IReadOnlyList<double> kernel)
    {
        var sum = kernel.Sum();
        return kernel.Select(value => value / sum).ToArray();
    }

    private static double[] SolveConstrainedDeconvolution(
        IReadOnlyList<double> raw,
        IReadOnlyList<double> weights,
        double targetArea,
        IReadOnlyList<double> kernel,
        HeatReleaseResponseCalibration calibration)
    {
        var current = raw.ToArray();
        for (var iteration = 0; iteration < calibration.Iterations; iteration++)
        {
            var predicted = Convolve(current, kernel);
            var residual = predicted.Zip(raw, (left, right) => left - right).ToArray();
            var gradient = TransposedConvolve(residual, kernel);
            AddSmoothnessGradient(current, gradient, calibration.RegularizationLambda);

            var candidate = current
                .Zip(gradient, (value, derivative) => value - calibration.GradientStep * derivative)
                .ToArray();
            if (candidate.Any(value => !double.IsFinite(value)))
            {
                throw new InvalidOperationException("HRR 校正数值不稳定，请检查响应核、步长和正则参数。");
            }

            var projected = ProjectNonNegativeWeightedArea(candidate, weights, targetArea);
            var maximumChange = projected
                .Zip(current, (left, right) => Math.Abs(left - right))
                .Max();
            var scale = Math.Max(1, projected.Max());
            current = projected;
            if (maximumChange / scale <= calibration.ConvergenceTolerance)
            {
                break;
            }
        }

        return current;
    }

    private static double[] Convolve(IReadOnlyList<double> values, IReadOnlyList<double> kernel)
    {
        var result = new double[values.Count];
        for (var i = 0; i < result.Length; i++)
        {
            for (var k = 0; k < kernel.Count && k <= i; k++)
            {
                result[i] += kernel[k] * values[i - k];
            }
        }

        return result;
    }

    private static double[] TransposedConvolve(IReadOnlyList<double> values, IReadOnlyList<double> kernel)
    {
        var result = new double[values.Count];
        for (var i = 0; i < result.Length; i++)
        {
            for (var k = 0; k < kernel.Count && i + k < values.Count; k++)
            {
                result[i] += kernel[k] * values[i + k];
            }
        }

        return result;
    }

    private static void AddSmoothnessGradient(
        IReadOnlyList<double> values,
        double[] gradient,
        double lambda)
    {
        if (lambda <= 0)
        {
            return;
        }

        for (var i = 1; i < values.Count - 1; i++)
        {
            var secondDifference = values[i - 1] - (2 * values[i]) + values[i + 1];
            gradient[i - 1] += lambda * secondDifference;
            gradient[i] -= 2 * lambda * secondDifference;
            gradient[i + 1] += lambda * secondDifference;
        }
    }

    private static double[] ProjectNonNegativeWeightedArea(
        IReadOnlyList<double> values,
        IReadOnlyList<double> weights,
        double targetArea)
    {
        var lower = -1d;
        var upper = 1d;
        while (WeightedProjectedArea(values, weights, lower) < targetArea)
        {
            lower *= 2;
        }

        while (WeightedProjectedArea(values, weights, upper) > targetArea)
        {
            upper *= 2;
        }

        for (var iteration = 0; iteration < 100; iteration++)
        {
            var midpoint = (lower + upper) / 2;
            if (WeightedProjectedArea(values, weights, midpoint) > targetArea)
            {
                lower = midpoint;
            }
            else
            {
                upper = midpoint;
            }
        }

        var threshold = (lower + upper) / 2;
        return values
            .Select((value, index) => Math.Max(0, value - threshold * weights[index]))
            .ToArray();
    }

    private static double WeightedProjectedArea(
        IReadOnlyList<double> values,
        IReadOnlyList<double> weights,
        double threshold)
    {
        var area = 0d;
        for (var i = 0; i < values.Count; i++)
        {
            area += weights[i] * Math.Max(0, values[i] - threshold * weights[i]);
        }

        return area;
    }

    private static List<RealtimeDataRecord> BuildCorrectedRecords(
        IReadOnlyList<RealtimeDataRecord> rawRecords,
        IReadOnlyList<double> correctedHrr,
        IReadOnlyList<double> weights,
        string calibrationVersion)
    {
        var result = new List<RealtimeDataRecord>(rawRecords.Count);
        var accumulatedThr = 0d;
        var correctedPeak = correctedHrr.Max();
        for (var i = 0; i < rawRecords.Count; i++)
        {
            accumulatedThr += correctedHrr[i] * weights[i] / TotalHeatReleaseDivisor;
            var raw = rawRecords[i];
            result.Add(raw with
            {
                HeatReleaseRate = correctedHrr[i],
                PeakHeatReleaseRate = correctedPeak,
                TotalHeatRelease = accumulatedThr,
                RawHeatReleaseRate = double.IsFinite(raw.RawHeatReleaseRate)
                    ? raw.RawHeatReleaseRate
                    : raw.HeatReleaseRate,
                RawTotalHeatRelease = double.IsFinite(raw.RawTotalHeatRelease)
                    ? raw.RawTotalHeatRelease
                    : raw.TotalHeatRelease,
                IsHeatReleaseCorrected = true,
                HeatReleaseCorrectionVersion = calibrationVersion
            });
        }

        return result;
    }

    private static HeatReleaseCorrectionResult NotApplied(
        IReadOnlyList<RealtimeDataRecord> rawRecords,
        string status)
    {
        var rawPeak = rawRecords
            .Where(record => double.IsFinite(record.HeatReleaseRate))
            .Select(record => record.HeatReleaseRate)
            .DefaultIfEmpty(double.NaN)
            .Max();
        var rawFinalThr = rawRecords
            .Where(record => double.IsFinite(record.TotalHeatRelease))
            .Select(record => record.TotalHeatRelease)
            .DefaultIfEmpty(double.NaN)
            .Last();

        Log.Warning("Heat release correction not applied. Reason={Reason}", status);
        return new HeatReleaseCorrectionResult(
            false,
            status,
            string.Empty,
            rawRecords.ToList(),
            rawPeak,
            double.NaN,
            rawFinalThr,
            double.NaN,
            double.NaN);
    }

    private static double Dot(IReadOnlyList<double> left, IReadOnlyList<double> right)
    {
        var result = 0d;
        for (var i = 0; i < left.Count; i++)
        {
            result += left[i] * right[i];
        }

        return result;
    }
}
