#include "nmCalculationAutoFitPSO.h"

#include <QtCore/qmath.h>
#include <algorithm>
#include <cmath>

#ifdef Q_OS_WIN
#include <float.h>
#endif

namespace {

static bool autoFitMetricIsFinite(double value)
{
#ifdef Q_OS_WIN
	return _finite(value) != 0;
#else
	return std::isfinite(value);
#endif
}

static void autoFitAppendMetricReason(QString* reasons,
	const QString& reason)
{
	if(!reasons || reason.isEmpty() || reasons->contains(reason)) {
		return;
	}
	if(!reasons->isEmpty()) {
		reasons->append(';');
	}
	reasons->append(reason);
}

struct AutoFitPressurePoint
{
	double timeHr;
	double pressureMpa;
};

static bool autoFitPressurePointLess(const AutoFitPressurePoint& left,
	const AutoFitPressurePoint& right)
{
	return left.timeHr < right.timeHr;
}

// 清理压力曲线并按时间升序排列。相同时间只保留最后一个值，避免插值区间为零。
static QVector<AutoFitPressurePoint> autoFitNormalizePressureCurve(
	const QVector<QVector<double> >& pressureData)
{
	QVector<AutoFitPressurePoint> points;
	if(pressureData.size() < 2) {
		return points;
	}

	const int count = qMin(pressureData[0].size(), pressureData[1].size());
	points.reserve(count);
	for(int i = 0; i < count; ++i) {
		const double timeHr = pressureData[0][i];
		const double pressureMpa = pressureData[1][i];
		if(!autoFitMetricIsFinite(timeHr) || timeHr <= 0.0 ||
		   !autoFitMetricIsFinite(pressureMpa)) {
			continue;
		}
		AutoFitPressurePoint point;
		point.timeHr = timeHr;
		point.pressureMpa = pressureMpa;
		points.append(point);
	}

	std::sort(points.begin(), points.end(), autoFitPressurePointLess);
	QVector<AutoFitPressurePoint> uniquePoints;
	uniquePoints.reserve(points.size());
	for(int i = 0; i < points.size(); ++i) {
		if(!uniquePoints.isEmpty() &&
		   qAbs(uniquePoints.last().timeHr - points[i].timeHr) <=
		   qMax(1.0e-14, points[i].timeHr * 1.0e-12)) {
			uniquePoints.last() = points[i];
		} else {
			uniquePoints.append(points[i]);
		}
	}
	return uniquePoints;
}

// 压力在 log10(t) 坐标中线性插值，与七算例统一误差协议保持一致。
static bool autoFitInterpolatePressure(
	const QVector<AutoFitPressurePoint>& points,
	double logTime,
	double* pressureMpa)
{
	if(!pressureMpa || points.size() < 2) {
		return false;
	}

	const double timeHr = qPow(10.0, logTime);
	if(timeHr < points.first().timeHr || timeHr > points.last().timeHr) {
		return false;
	}

	int left = 0;
	int right = points.size() - 1;
	while(right - left > 1) {
		const int middle = left + (right - left) / 2;
		if(points[middle].timeHr <= timeHr) {
			left = middle;
		} else {
			right = middle;
		}
	}

	if(qAbs(timeHr - points[left].timeHr) <=
	   qMax(1.0e-14, timeHr * 1.0e-12)) {
		*pressureMpa = points[left].pressureMpa;
		return true;
	}
	if(qAbs(timeHr - points[right].timeHr) <=
	   qMax(1.0e-14, timeHr * 1.0e-12)) {
		*pressureMpa = points[right].pressureMpa;
		return true;
	}

	const double leftLogTime = qLn(points[left].timeHr) / qLn(10.0);
	const double rightLogTime = qLn(points[right].timeHr) / qLn(10.0);
	const double denominator = rightLogTime - leftLogTime;
	if(denominator <= 0.0) {
		return false;
	}
	const double ratio = (logTime - leftLogTime) / denominator;
	*pressureMpa = points[left].pressureMpa +
		(points[right].pressureMpa - points[left].pressureMpa) * ratio;
	return autoFitMetricIsFinite(*pressureMpa);
}

// 在统一时间网格上按 Bourdet 三点公式计算 d(DeltaP)/d(ln t)。首尾点
// 没有完整邻点，保持 NaN 并且不参与导数 RMSE。
static bool autoFitCalculateBourdetDerivative(
	const QVector<double>& timeHr,
	const QVector<double>& deltaPMpa,
	QVector<double>* derivativeMpa,
	QString* invalidReason,
	bool* hasInvalidDerivative)
{
	if(hasInvalidDerivative) {
		*hasInvalidDerivative = false;
	}
	if(!derivativeMpa || timeHr.size() != deltaPMpa.size() ||
	   timeHr.size() < 3) {
		if(invalidReason) {
			*invalidReason = "INSUFFICIENT_DERIVATIVE_POINTS";
		}
		return false;
	}

	const double invalidValue = std::numeric_limits<double>::quiet_NaN();
	derivativeMpa->fill(invalidValue, timeHr.size());
	for(int i = 1; i < timeHr.size() - 1; ++i) {
		const double leftInterval = qLn(timeHr[i] / timeHr[i - 1]);
		const double rightInterval = qLn(timeHr[i + 1] / timeHr[i]);
		const double totalInterval = leftInterval + rightInterval;
		if(!autoFitMetricIsFinite(leftInterval) || leftInterval <= 0.0 ||
		   !autoFitMetricIsFinite(rightInterval) || rightInterval <= 0.0 ||
		   !autoFitMetricIsFinite(totalInterval) || totalInterval <= 0.0) {
			if(invalidReason) {
				*invalidReason = "INVALID_TIME_INTERVAL";
			}
			return false;
		}

		const double leftSlope =
			(deltaPMpa[i] - deltaPMpa[i - 1]) / leftInterval;
		const double rightSlope =
			(deltaPMpa[i + 1] - deltaPMpa[i]) / rightInterval;
		const double derivative =
			leftSlope * rightInterval / totalInterval +
			rightSlope * leftInterval / totalInterval;
		if(!autoFitMetricIsFinite(derivative) || derivative <= 0.0) {
			// 单点异常不破坏其余 Bourdet 点。该位置保持 NaN，由调用方从
			// 双对数导数 RMSE 中剔除，同时保留整条曲线的数据质量失败标志。
			if(hasInvalidDerivative) {
				*hasInvalidDerivative = true;
			}
			continue;
		}
		(*derivativeMpa)[i] = derivative;
	}
	return true;
}

} // namespace

AutoFitCurveMetrics::AutoFitCurveMetrics()
	: valid(false)
	, passed(false)
	, coverage(std::numeric_limits<double>::quiet_NaN())
	, pressureRmseMpa(std::numeric_limits<double>::quiet_NaN())
	, pressureMaxAbsErrorMpa(std::numeric_limits<double>::quiet_NaN())
	, logDeltaPRmseDecade(std::numeric_limits<double>::quiet_NaN())
	, logDerivativeRmseDecade(std::numeric_limits<double>::quiet_NaN())
	, unifiedCurveError(std::numeric_limits<double>::quiet_NaN())
	, sampleCount(0)
	, validDerivativeCount(0)
{
}

AutoFitParameterResult::AutoFitParameterResult()
	: initialValue(std::numeric_limits<double>::quiet_NaN())
	, lowerBound(std::numeric_limits<double>::quiet_NaN())
	, upperBound(std::numeric_limits<double>::quiet_NaN())
	, finalValue(std::numeric_limits<double>::quiet_NaN())
{
}

AutoFitRunResult::AutoFitRunResult()
	: ompThreads(-1)
	, iluReuseSteps(-1)
	, optimizationWallTimeMs(-1)
	, workflowWallTimeMs(-1)
	, solverTimeSumMs(0)
	, finalSolverTimeMs(-1)
	, iterationCount(0)
	, parameterEvaluationCount(0)
	, modelSolverCallCount(0)
	, finalSolverCallCount(0)
	, solverSuccessCount(0)
	, solverFailureCount(0)
	, solverTimeoutCount(0)
	, optimizationPebiCount(-1)
	, finalPebiCount(-1)
	, pebiCount(-1)
	, finalSolverStatus("NOT_RUN")
	, initialPressureMpa(std::numeric_limits<double>::quiet_NaN())
	, initialInternalError(std::numeric_limits<double>::quiet_NaN())
	, finalInternalError(std::numeric_limits<double>::quiet_NaN())
{
}

void AutoFitRunResult::recordOptimizationSolverCall(bool success,
	bool timeout,
	qint64 solveTimeMs,
	int pebiCount)
{
	++modelSolverCallCount;
	if(timeout) {
		++solverTimeoutCount;
	} else if(success) {
		++solverSuccessCount;
	} else {
		++solverFailureCount;
	}
	if(solveTimeMs >= 0) {
		solverTimeSumMs += solveTimeMs;
	}
	if(pebiCount >= 0) {
		optimizationPebiCount = pebiCount;
		this->pebiCount = pebiCount;
	}
}

void AutoFitRunResult::recordFinalSolverCall(bool success,
	bool timeout,
	qint64 solveTimeMs,
	int pebiCount)
{
	++finalSolverCallCount;
	finalSolverTimeMs = solveTimeMs >= 0 ? solveTimeMs : -1;
	finalPebiCount = pebiCount >= 0 ? pebiCount : -1;
	if(pebiCount >= 0) {
		this->pebiCount = pebiCount;
	}
	if(timeout) {
		finalSolverStatus = "TIMEOUT";
	} else {
		finalSolverStatus = success ? "SUCCESS" : "FAILED";
	}
}

AutoFitRunResult nmCalculationAutoFitPSO::getLastRunResult() const
{
	return m_lastRunResult;
}

AutoFitCurveMetrics nmCalculationAutoFitPSO::calculateUnifiedCurveMetrics(
	const QVector<QVector<double> >& targetPressureData,
	const QVector<QVector<double> >& fittedPressureData,
	double initialPressureMpa,
	int sampleCount)
{
	AutoFitCurveMetrics metrics;
	if(!autoFitMetricIsFinite(initialPressureMpa)) {
		metrics.invalidReason = "INVALID_INITIAL_PRESSURE";
		return metrics;
	}
	if(sampleCount < 3) {
		metrics.invalidReason = "INVALID_SAMPLE_COUNT";
		return metrics;
	}

	const QVector<AutoFitPressurePoint> targetPoints =
		autoFitNormalizePressureCurve(targetPressureData);
	const QVector<AutoFitPressurePoint> fittedPoints =
		autoFitNormalizePressureCurve(fittedPressureData);
	if(targetPoints.size() < 2) {
		metrics.invalidReason = "MISSING_TARGET_PRESSURE_DATA";
		return metrics;
	}
	if(fittedPoints.size() < 2) {
		metrics.invalidReason = "MISSING_FITTED_PRESSURE_DATA";
		return metrics;
	}

	const double targetMinLogTime =
		qLn(targetPoints.first().timeHr) / qLn(10.0);
	const double targetMaxLogTime =
		qLn(targetPoints.last().timeHr) / qLn(10.0);
	const double fittedMinLogTime =
		qLn(fittedPoints.first().timeHr) / qLn(10.0);
	const double fittedMaxLogTime =
		qLn(fittedPoints.last().timeHr) / qLn(10.0);
	const double targetLogSpan = targetMaxLogTime - targetMinLogTime;
	if(!autoFitMetricIsFinite(targetLogSpan) || targetLogSpan <= 0.0) {
		metrics.invalidReason = "INVALID_TARGET_TIME_RANGE";
		return metrics;
	}

	const double commonMinLogTime = qMax(targetMinLogTime, fittedMinLogTime);
	const double commonMaxLogTime = qMin(targetMaxLogTime, fittedMaxLogTime);
	const double commonLogSpan = commonMaxLogTime - commonMinLogTime;
	metrics.coverage = qBound(0.0, commonLogSpan / targetLogSpan, 1.0);
	if(!autoFitMetricIsFinite(commonLogSpan) || commonLogSpan <= 0.0) {
		metrics.invalidReason = "NO_COMMON_TIME_RANGE";
		return metrics;
	}

	metrics.sampleCount = sampleCount;
	metrics.timeHr.reserve(sampleCount);
	metrics.targetPressureMpa.reserve(sampleCount);
	metrics.fittedPressureMpa.reserve(sampleCount);
	metrics.targetDeltaPMpa.reserve(sampleCount);
	metrics.fittedDeltaPMpa.reserve(sampleCount);

	double pressureSquaredSum = 0.0;
	double pressureMaxAbsError = 0.0;
	double logDeltaPSquaredSum = 0.0;
	int validLogDeltaPCount = 0;
	bool hasInvalidDeltaP = false;
	for(int i = 0; i < sampleCount; ++i) {
		const double ratio = sampleCount > 1
			? static_cast<double>(i) / static_cast<double>(sampleCount - 1)
			: 0.0;
		const double logTime = commonMinLogTime + commonLogSpan * ratio;
		double targetPressure = 0.0;
		double fittedPressure = 0.0;
		if(!autoFitInterpolatePressure(targetPoints, logTime, &targetPressure) ||
		   !autoFitInterpolatePressure(fittedPoints, logTime, &fittedPressure)) {
			metrics.invalidReason = "PRESSURE_INTERPOLATION_FAILED";
			return metrics;
		}

		const double timeHr = qPow(10.0, logTime);
		const double targetDeltaP = qAbs(initialPressureMpa - targetPressure);
		const double fittedDeltaP = qAbs(initialPressureMpa - fittedPressure);
		metrics.timeHr.append(timeHr);
		metrics.targetPressureMpa.append(targetPressure);
		metrics.fittedPressureMpa.append(fittedPressure);
		metrics.targetDeltaPMpa.append(targetDeltaP);
		metrics.fittedDeltaPMpa.append(fittedDeltaP);

		const double pressureResidual = fittedPressure - targetPressure;
		pressureSquaredSum += pressureResidual * pressureResidual;
		pressureMaxAbsError = qMax(pressureMaxAbsError, qAbs(pressureResidual));
		if(!autoFitMetricIsFinite(targetDeltaP) || targetDeltaP <= 0.0 ||
		   !autoFitMetricIsFinite(fittedDeltaP) || fittedDeltaP <= 0.0) {
			hasInvalidDeltaP = true;
			continue;
		}

		const double logResidual =
			qLn(fittedDeltaP / targetDeltaP) / qLn(10.0);
		if(autoFitMetricIsFinite(logResidual)) {
			logDeltaPSquaredSum += logResidual * logResidual;
			++validLogDeltaPCount;
		} else {
			hasInvalidDeltaP = true;
		}
	}

	metrics.pressureRmseMpa =
		qSqrt(pressureSquaredSum / static_cast<double>(sampleCount));
	metrics.pressureMaxAbsErrorMpa = pressureMaxAbsError;
	if(validLogDeltaPCount > 0) {
		metrics.logDeltaPRmseDecade = qSqrt(
			logDeltaPSquaredSum /
			static_cast<double>(validLogDeltaPCount));
	} else {
		autoFitAppendMetricReason(&metrics.invalidReason,
			"NO_VALID_LOG_DELTA_P_SAMPLES");
	}
	if(hasInvalidDeltaP) {
		autoFitAppendMetricReason(&metrics.invalidReason,
			"NON_POSITIVE_OR_INVALID_DELTA_P");
	}

	QString targetDerivativeReason;
	QString fittedDerivativeReason;
	bool targetHasInvalidDerivative = false;
	bool fittedHasInvalidDerivative = false;
	if(!autoFitCalculateBourdetDerivative(metrics.timeHr,
		metrics.targetDeltaPMpa,
		&metrics.targetDerivativeMpa,
		&targetDerivativeReason,
		&targetHasInvalidDerivative)) {
		autoFitAppendMetricReason(&metrics.invalidReason,
			targetDerivativeReason);
		return metrics;
	}
	if(!autoFitCalculateBourdetDerivative(metrics.timeHr,
		metrics.fittedDeltaPMpa,
		&metrics.fittedDerivativeMpa,
		&fittedDerivativeReason,
		&fittedHasInvalidDerivative)) {
		autoFitAppendMetricReason(&metrics.invalidReason,
			fittedDerivativeReason);
		return metrics;
	}

	double logDerivativeSquaredSum = 0.0;
	for(int i = 1; i < sampleCount - 1; ++i) {
		const double targetDerivative = metrics.targetDerivativeMpa[i];
		const double fittedDerivative = metrics.fittedDerivativeMpa[i];
		if(!autoFitMetricIsFinite(targetDerivative) ||
		   targetDerivative <= 0.0 ||
		   !autoFitMetricIsFinite(fittedDerivative) ||
		   fittedDerivative <= 0.0) {
			continue;
		}
		const double logResidual =
			qLn(fittedDerivative / targetDerivative) / qLn(10.0);
		if(!autoFitMetricIsFinite(logResidual)) {
			targetHasInvalidDerivative = true;
			continue;
		}
		logDerivativeSquaredSum += logResidual * logResidual;
		++metrics.validDerivativeCount;
	}

	if(metrics.validDerivativeCount > 0) {
		metrics.logDerivativeRmseDecade = qSqrt(
			logDerivativeSquaredSum /
			static_cast<double>(metrics.validDerivativeCount));
	} else {
		autoFitAppendMetricReason(&metrics.invalidReason,
			"NO_VALID_LOG_DERIVATIVE_SAMPLES");
	}
	if(targetHasInvalidDerivative || fittedHasInvalidDerivative) {
		autoFitAppendMetricReason(&metrics.invalidReason,
			"NON_POSITIVE_OR_INVALID_DERIVATIVE");
	}
	if(autoFitMetricIsFinite(metrics.logDeltaPRmseDecade) &&
	   autoFitMetricIsFinite(metrics.logDerivativeRmseDecade)) {
		metrics.unifiedCurveError = qSqrt(
			(metrics.logDeltaPRmseDecade * metrics.logDeltaPRmseDecade +
			 metrics.logDerivativeRmseDecade * metrics.logDerivativeRmseDecade) /
			2.0);
	}

	if(metrics.coverage < 0.95) {
		autoFitAppendMetricReason(&metrics.invalidReason,
			"TIME_COVERAGE_BELOW_0_95");
	}

	metrics.valid = metrics.invalidReason.isEmpty();
	metrics.passed = metrics.valid &&
		metrics.logDeltaPRmseDecade <= 0.02 &&
		metrics.logDerivativeRmseDecade <= 0.02;
	return metrics;
}