nmWTAI-Platform/Src/nmNum/nmData/nmDataDiagnostic.cpp

#include "nmDataDiagnostic.h"
#include "nmDataWellBase.h"
#include "nmDataReservoir.h"
#include <QDebug>
#include "qmath.h"
#include <QVector>
#include <algorithm>

#ifdef Q_OS_WIN
#include <windows.h>
#define DEBUG_OUT(msg) OutputDebugStringA(QString("%1\n").arg(msg).toLocal8Bit().data())
#endif

nmDataDiagnostic::nmDataDiagnostic()
{
	m_diagnosticSkin = nmDataAttribute("Skin", 0.0, "");
	m_diagnosticWellboreStorage = nmDataAttribute("Wellbore storage", 0.01, "m^3/MPa", UNIT_TYPE_COMPRESSIBILITY, QStringList(), QStringList() << "bbl/psi" << "m^3/bar" << "m^3/kPa" << "m^3/Pa" << "m^3.cm^2/kg" << "m^2" << "m^3/MPa");
	m_diagnosticPermeability = nmDataAttribute("Permeability", 0.025, "md", UNIT_TYPE_PERMEABILITY, QStringList(), QStringList() << "md" << "Darcy" << "m^2" << "cm^2" << "um^2");
	m_diagnosticTransmissibility = nmDataAttribute("Transmissibility", 1000.0, "md.m", UNIT_TYPE_CONDUCTIVITY, QStringList(), QStringList()<< "md.ft" << "md.m" << "m^3");
}

nmDataDiagnostic::~nmDataDiagnostic() {

}

// 序列化 nmDataDiagnostic 为 RapidJSON Value
rapidjson::Value nmDataDiagnostic::ToJsonValue(rapidjson::Document::AllocatorType& allocator) const
{
	// 创建一个 RapidJSON 对象类型的值
	rapidjson::Value diagnosticObject(rapidjson::kObjectType);
	// 序列化 nmDataAttribute 类型的成员
	// 调用 nmDataAttribute 自身的 ToJsonValue 方法进行递归序列化
	diagnosticObject.AddMember("DiagnosticSkin", m_diagnosticSkin.ToJsonValue(allocator), allocator);
	diagnosticObject.AddMember("DiagnosticWellboreStorage", m_diagnosticWellboreStorage.ToJsonValue(allocator), allocator);
	diagnosticObject.AddMember("DiagnosticPermeability", m_diagnosticPermeability.ToJsonValue(allocator), allocator);
	diagnosticObject.AddMember("DiagnosticTransmissibility", m_diagnosticTransmissibility.ToJsonValue(allocator), allocator);
	return diagnosticObject; // 返回序列化后的 RapidJSON Value
}

// 从 RapidJSON Value 反序列化数据到 nmDataDiagnostic
void nmDataDiagnostic::FromJsonValue(const rapidjson::Value& jsonValue)
{
	// 反序列化 nmDataAttribute 类型的成员
	// 调用 nmDataAttribute 自身的 FromJsonValue 方法进行递归反序列化
	if (jsonValue.HasMember("DiagnosticSkin") && jsonValue["DiagnosticSkin"].IsObject()) {
		m_diagnosticSkin.FromJsonValue(jsonValue["DiagnosticSkin"]);
	}
	if (jsonValue.HasMember("DiagnosticWellboreStorage") && jsonValue["DiagnosticWellboreStorage"].IsObject()) {
		m_diagnosticWellboreStorage.FromJsonValue(jsonValue["DiagnosticWellboreStorage"]);
	}
	if (jsonValue.HasMember("DiagnosticPermeability") && jsonValue["DiagnosticPermeability"].IsObject()) {
		m_diagnosticPermeability.FromJsonValue(jsonValue["DiagnosticPermeability"]);
	}
	if (jsonValue.HasMember("DiagnosticTransmissibility") && jsonValue["DiagnosticTransmissibility"].IsObject()) {
		m_diagnosticTransmissibility.FromJsonValue(jsonValue["DiagnosticTransmissibility"]);
	}
}

// 获取诊断结果参数
nmDataAttribute& nmDataDiagnostic::getDiagnosticPermeability() { return m_diagnosticPermeability; }
nmDataAttribute& nmDataDiagnostic::getDiagnosticSkin() { return m_diagnosticSkin; }
nmDataAttribute& nmDataDiagnostic::getDiagnosticWellboreStorage() { return m_diagnosticWellboreStorage; }
nmDataAttribute& nmDataDiagnostic::getDiagnosticTransmissibility() { return m_diagnosticTransmissibility; }



void nmDataDiagnostic::resetFromDiagnostic(const QVector<QVector<double>>& rawData, 
	nmDataWellBase* pWell, 
	nmDataReservoir* pReservoir)
{
	if (!pWell || !pReservoir) {
		DEBUG_OUT("Invalid well or reservoir data");
		return;
	}

	if (rawData.size() < 2 || rawData[0].isEmpty() || rawData[1].isEmpty()) {
		DEBUG_OUT("Invalid rawData format");
		return;
	}

	DEBUG_OUT("=== Starting Intelligent Diagnostic Analysis ===");

	// 调用各个参数的诊断计算方法
	calculateWellboreStorage(rawData, pWell, pReservoir);
	calculateTransmissibility(rawData, pWell, pReservoir);  // 先计算 kh
	calculatePermeability(pWell, pReservoir);              // 再基于 kh 计算 k
	//calculateSkin(rawData, pWell, pReservoir);

	DEBUG_OUT("=== Diagnostic Analysis Complete ===");
}

void nmDataDiagnostic::calculateWellboreStorage(const QVector<QVector<double>>& rawData,
	nmDataWellBase* pWell, nmDataReservoir* pReservoir)
{
	DEBUG_OUT("=== Intelligent Wellbore Storage Analysis ===");

	const QVector<double>& timeData = rawData[0];      // 时间数据 (hr)
	const QVector<double>& pressureDropData = rawData[1];  // 压力差数据 (MPa)

	// 确保有足够的数据点
	int availablePoints = qMin(timeData.size(), pressureDropData.size());

	if (availablePoints < 3) {
		DEBUG_OUT("Not enough data points for wellbore storage analysis (need at least 3 points)");
		return;
	}

	// 只分析早期数据点（前20%或前30个点，取较小值）
	int earlyPointsToUse = qMin(qMax(5, static_cast<int>(availablePoints * 0.2)), 30);
	earlyPointsToUse = qMin(earlyPointsToUse, availablePoints);

	DEBUG_OUT(QString("Total points: %1, analyzing early %2 points (%3%) for flow regime identification")
		.arg(availablePoints).arg(earlyPointsToUse).arg(100.0 * earlyPointsToUse / availablePoints, 0, 'f', 1));

	// 获取储层和井参数
	double viscosity = pReservoir->getMiuo().getValue().toDouble();  // 粘度
	double bo = pReservoir->getBo().getValue().toDouble();  // 体积系数


	DEBUG_OUT(QString("Reservoir parameters: viscosity=%1, Bo=%2").arg(viscosity).arg(bo));

	// 智能识别流动阶段：分析wellbore storage和径向流两种机制
	FlowRegimeAnalysis wellboreStorageAnalysis = analyzeFlowRegime(timeData, pressureDropData, 
		earlyPointsToUse, 1.0, 
		"Wellbore Storage", 
		pWell, pReservoir);

	FlowRegimeAnalysis radialFlowAnalysis = analyzeFlowRegime(timeData, pressureDropData, 
		earlyPointsToUse, 0.5, 
		"Radial Flow", 
		pWell, pReservoir);

	// 智能选择最佳结果 - 改进的选择逻辑
	FlowRegimeAnalysis* selectedAnalysis = nullptr;
	QString selectionReason = "";

	if (wellboreStorageAnalysis.isValid && radialFlowAnalysis.isValid) {
		// 检查每个斜率更接近哪个理论值
		double wbSlope = wellboreStorageAnalysis.slope;
		double rfSlope = radialFlowAnalysis.slope;

		// 计算斜率到两个理论值的距离
		double wbDistanceTo1 = qAbs(wbSlope - 1.0);
		double wbDistanceTo05 = qAbs(wbSlope - 0.5);
		double rfDistanceTo1 = qAbs(rfSlope - 1.0);
		double rfDistanceTo05 = qAbs(rfSlope - 0.5);

		// 判断wellbore storage斜率更接近哪个理论值
		bool wbCloserTo1 = wbDistanceTo1 < wbDistanceTo05;
		// 判断radial flow斜率更接近哪个理论值  
		bool rfCloserTo05 = rfDistanceTo05 < rfDistanceTo1;

		DEBUG_OUT(QString("=== Selection Analysis ==="));
		DEBUG_OUT(QString("WB slope %1: distance to 1.0=%2, to 0.5=%3, closer to %4")
			.arg(wbSlope, 0, 'f', 4).arg(wbDistanceTo1, 0, 'f', 4).arg(wbDistanceTo05, 0, 'f', 4)
			.arg(wbCloserTo1 ? "1.0" : "0.5"));
		DEBUG_OUT(QString("RF slope %1: distance to 1.0=%2, to 0.5=%3, closer to %4")
			.arg(rfSlope, 0, 'f', 4).arg(rfDistanceTo1, 0, 'f', 4).arg(rfDistanceTo05, 0, 'f', 4)
			.arg(rfCloserTo05 ? "0.5" : "1.0"));

		// 决策逻辑：
		if (wbCloserTo1 && wellboreStorageAnalysis.avgTime < 0.05) {
			// 如果wellbore storage斜率确实更接近1.0且时间很早，优先选择wellbore storage
			selectedAnalysis = &wellboreStorageAnalysis;
			selectionReason = QString("Wellbore storage slope (%1) closer to 1.0 and early time (%2hr)")
				.arg(wbSlope, 0, 'f', 4).arg(wellboreStorageAnalysis.avgTime, 0, 'f', 6);
		} else if (rfCloserTo05) {
			// 如果radial flow斜率确实更接近0.5，选择radial flow
			selectedAnalysis = &radialFlowAnalysis;
			selectionReason = QString("Radial flow slope (%1) closer to 0.5, difference=%2")
				.arg(rfSlope, 0, 'f', 4).arg(radialFlowAnalysis.slopeDifference, 0, 'f', 4);
		} else if (wbCloserTo1) {
			// 如果wellbore storage斜率更接近1.0，选择wellbore storage
			selectedAnalysis = &wellboreStorageAnalysis;
			selectionReason = QString("Wellbore storage slope (%1) closer to 1.0, difference=%2")
				.arg(wbSlope, 0, 'f', 4).arg(wellboreStorageAnalysis.slopeDifference, 0, 'f', 4);
		} else {
			// 兜底：选择斜率差异更小的
			if (wellboreStorageAnalysis.slopeDifference < radialFlowAnalysis.slopeDifference) {
				selectedAnalysis = &wellboreStorageAnalysis;
				selectionReason = QString("Fallback: WB slope difference (%1) < RF difference (%2)")
					.arg(wellboreStorageAnalysis.slopeDifference, 0, 'f', 4)
					.arg(radialFlowAnalysis.slopeDifference, 0, 'f', 4);
			} else {
				selectedAnalysis = &radialFlowAnalysis;
				selectionReason = QString("Fallback: RF slope difference (%1) < WB difference (%2)")
					.arg(radialFlowAnalysis.slopeDifference, 0, 'f', 4)
					.arg(wellboreStorageAnalysis.slopeDifference, 0, 'f', 4);
			}
		}
	} else if (wellboreStorageAnalysis.isValid) {
		selectedAnalysis = &wellboreStorageAnalysis;
		selectionReason = "Only wellbore storage analysis is valid";
	} else if (radialFlowAnalysis.isValid) {
		selectedAnalysis = &radialFlowAnalysis;
		selectionReason = "Only radial flow analysis is valid";
	}

	// 设置诊断结果
	if (selectedAnalysis && selectedAnalysis->result > 0) {
		m_diagnosticWellboreStorage.setValue(selectedAnalysis->result);

		DEBUG_OUT(QString("=== Final Result ==="));
		DEBUG_OUT(QString("Selected method: %1").arg(selectedAnalysis->description));
		DEBUG_OUT(QString("Selection reason: %1").arg(selectionReason));
		DEBUG_OUT(QString("Final wellbore storage: %1").arg(selectedAnalysis->result));
		DEBUG_OUT(QString("Based on points [%1,%2] at times %3hr and %4hr")
			.arg(selectedAnalysis->pointPair.first).arg(selectedAnalysis->pointPair.second)
			.arg(timeData[selectedAnalysis->pointPair.first])
			.arg(timeData[selectedAnalysis->pointPair.second]));
	} else {
		DEBUG_OUT("Unable to calculate wellbore storage - no valid flow regime identified");
	}
}

void nmDataDiagnostic::calculateTransmissibility(const QVector<QVector<double>>& rawData,
    nmDataWellBase* pWell, nmDataReservoir* pReservoir)
{
    DEBUG_OUT("=== Transmissibility (kh) Analysis ===");

    // 检查数据格式：需要3个向量 [时间, 压力, 导数]
    if (rawData.size() < 3 || rawData[0].isEmpty() || rawData[1].isEmpty() || rawData[2].isEmpty()) {
        DEBUG_OUT("Invalid rawData format: need [time, pressure, derivative] vectors");
        return;
    }

    const QVector<double>& timeData = rawData[0];          // 时间数据 (hr)
    const QVector<double>& derivativeData = rawData[2];    // 压力导数数据 (MPa)

    // 确保有足够的数据点
    int availablePoints = timeData.size();
    if (availablePoints < 10) {
        DEBUG_OUT("Not enough data points for transmissibility analysis (need at least 10 points)");
        return;
    }

    DEBUG_OUT("Trying multiple analysis methods to find best transmissibility estimate...");

    // 方法1：标准径向流平台分析
    FlowRegimeAnalysis standardRadialFlow = analyzeStandardRadialFlow(
        timeData, derivativeData, availablePoints, pWell, pReservoir);

    // 方法2：晚期分析（使用极小导数值）
    FlowRegimeAnalysis saphirStyleAnalysis = analyzeSaphirStyleLateTime(
        timeData, derivativeData, availablePoints, pWell, pReservoir);

    // 方法3：早中期稳定段分析
    FlowRegimeAnalysis earlyStableAnalysis = analyzeEarlyStableFlow(
        timeData, derivativeData, availablePoints, pWell, pReservoir);

    // 智能选择最佳方法
    FlowRegimeAnalysis* selectedAnalysis = selectBestTransmissibilityAnalysis(
        standardRadialFlow, saphirStyleAnalysis, earlyStableAnalysis);

    if (selectedAnalysis && selectedAnalysis->isValid && selectedAnalysis->result > 0) {
        double khValue = selectedAnalysis->result;
        m_diagnosticTransmissibility.setValue(khValue);

        DEBUG_OUT(QString("=== Transmissibility Analysis Result ==="));
        DEBUG_OUT(QString("Selected method: %1").arg(selectedAnalysis->description));
        DEBUG_OUT(QString("Transmissibility (kh): %1 md*m").arg(khValue, 0, 'f', 6));
        DEBUG_OUT(QString("Analysis time: %1 hr").arg(selectedAnalysis->avgTime, 0, 'f', 4));
        DEBUG_OUT(QString("Derivative value used: %1 MPa").arg(selectedAnalysis->avgPressure, 0, 'f', 6));
    } else {
        DEBUG_OUT("All transmissibility analysis methods failed");
    }
}

void nmDataDiagnostic::calculatePermeability(nmDataWellBase* pWell, nmDataReservoir* pReservoir)
{
	DEBUG_OUT("=== Permeability (k) Analysis ===");

	// 获取已计算的 kh 值
	double khValue = m_diagnosticTransmissibility.getValue().toDouble();

	if (khValue <= 0) {
		DEBUG_OUT("Invalid transmissibility value - cannot calculate permeability");
		DEBUG_OUT("Please ensure transmissibility is calculated first");
		return;
	}

	// 获取地层厚度
	double thickness = pReservoir->getThickness().getValue().toDouble();    // 地层厚度 (m)

	if (thickness <= 0) {
		DEBUG_OUT("Invalid reservoir thickness for permeability calculation");
		return;
	}

	// 计算渗透率 k = kh / h
	double kValue = khValue / thickness;
	m_diagnosticPermeability.setValue(kValue);

	DEBUG_OUT(QString("=== Permeability Analysis Result ==="));
	DEBUG_OUT(QString("Transmissibility (kh): %1 md*m").arg(khValue, 0, 'f', 6));
	DEBUG_OUT(QString("Reservoir thickness (h): %1 m").arg(thickness, 0, 'f', 2));
	DEBUG_OUT(QString("Permeability (k): %1 md").arg(kValue, 0, 'f', 6));
}

//void nmDataDiagnostic::calculateSkin(const QVector<QVector<double>>& rawData,
//	nmDataWellBase* pWell, nmDataReservoir* pReservoir)
//{
//
//}

nmDataDiagnostic::FlowRegimeAnalysis nmDataDiagnostic::analyzeFlowRegime(
	const QVector<double>& timeData,
	const QVector<double>& pressureDropData,
	int pointCount, double targetSlope,
	const QString& regimeName,
	nmDataWellBase* pWell, nmDataReservoir* pReservoir)
{
	FlowRegimeAnalysis analysis;
	analysis.pointPair = findClosestSlopeToValue(timeData, pressureDropData, pointCount, targetSlope);

	if (analysis.pointPair.first >= 0 && analysis.pointPair.second >= 0) {
		int i = analysis.pointPair.first;
		int j = analysis.pointPair.second;

		double t1 = timeData[i];
		double t2 = timeData[j];
		double deltaP1 = pressureDropData[i];
		double deltaP2 = pressureDropData[j];

		analysis.slope = calculateSlopeBetweenPoints(t1, deltaP1, t2, deltaP2);
		analysis.slopeDifference = qAbs(analysis.slope - targetSlope);
		analysis.avgTime = (t1 + t2) / 2.0;
		analysis.avgPressure = (deltaP1 + deltaP2) / 2.0;

		// 获取流量
		double flowRate = pWell->getFlowRateAtTime(analysis.avgTime);

		if (flowRate > 0) {
			// 获取储层参数
			double viscosity = pReservoir->getMiuo().getValue().toDouble();
			double bo = pReservoir->getBo().getValue().toDouble();

			// 使用公式计算wellbore storage: c = qb*μo*Bo*t/(24*ΔP)
			analysis.result = (flowRate * viscosity * bo * analysis.avgTime) / (24.0 * analysis.avgPressure);
			analysis.isValid = true;

			analysis.description = QString("%1 (slope=%2, target=%3, diff=%4)")
				.arg(regimeName).arg(analysis.slope, 0, 'f', 4)
				.arg(targetSlope, 0, 'f', 1).arg(analysis.slopeDifference, 0, 'f', 4);

			DEBUG_OUT(QString("=== %1 Analysis ===").arg(regimeName));
			DEBUG_OUT(QString("Best points: [%1,%2] at times %3hr - %4hr")
				.arg(i).arg(j).arg(t1, 0, 'f', 6).arg(t2, 0, 'f', 6));
			DEBUG_OUT(QString("Slope: %1 (target: %2, difference: %3)")
				.arg(analysis.slope, 0, 'f', 4).arg(targetSlope, 0, 'f', 1)
				.arg(analysis.slopeDifference, 0, 'f', 4));
			DEBUG_OUT(QString("Average time: %1hr, Average derTaP: %2").arg(analysis.avgTime, 0, 'f', 6).arg(analysis.avgPressure, 0, 'f', 6));
			DEBUG_OUT(QString("Flow rate: %1, Result: %2").arg(flowRate).arg(analysis.result, 0, 'f', 6));
		}
	}

	return analysis;
}

QPair<int, int> nmDataDiagnostic::findClosestSlopeToValue(const QVector<double>& timeData, 
	const QVector<double>& pressureDropData, 
	int pointCount, double targetSlope)
{
	double bestSlope = -1.0;
	double minDifference = 999999.0;
	QPair<int, int> bestPair(-1, -1);

	// 遍历所有可能的点对，寻找斜率最接近目标值的
	for (int i = 0; i < pointCount - 1; ++i) {
		for (int j = i + 1; j < pointCount; ++j) {
			double t1 = timeData[i];
			double t2 = timeData[j];
			double deltaP1 = pressureDropData[i];
			double deltaP2 = pressureDropData[j];

			// 确保数据有效
			if (t1 <= 0 || t2 <= 0 || deltaP1 <= 0 || deltaP2 <= 0 || t2 <= t1) {
				continue;
			}

			// 计算斜率
			double slope = calculateSlopeBetweenPoints(t1, deltaP1, t2, deltaP2);

			// 计算与目标值的差异
			double difference = qAbs(slope - targetSlope);

			// 如果这个斜率更接近目标值，更新最佳结果
			if (difference < minDifference) {
				minDifference = difference;
				bestSlope = slope;
				bestPair.first = i;
				bestPair.second = j;
			}
		}
	}

	return bestPair;
}

double nmDataDiagnostic::calculateSlopeBetweenPoints(double t1, double deltaP1, double t2, double deltaP2)
{
	if (t1 <= 0 || t2 <= 0 || deltaP1 <= 0 || deltaP2 <= 0 || t2 <= t1) {
		return 0.0;
	}

	// 计算双对数斜率: slope = log(ΔP2/ΔP1) / log(t2/t1)
	double logTimeRatio = log10(t2 / t1);
	double logPressureRatio = log10(deltaP2 / deltaP1);

	if (qAbs(logTimeRatio) < 1e-10) {  // 避免除以零
		return 0.0;
	}

	double slope = logPressureRatio / logTimeRatio;

	return slope;
}

double nmDataDiagnostic::interpolatePressureAtTime(const QVector<double>& timeData, 
	const QVector<double>& pressureDropData, 
	double targetTime)
{
	for (int i = 0; i < timeData.size() - 1; ++i) {
		if (targetTime >= timeData[i] && targetTime <= timeData[i + 1]) {
			// 线性插值
			double ratio = (targetTime - timeData[i]) / (timeData[i + 1] - timeData[i]);
			return pressureDropData[i] + ratio * (pressureDropData[i + 1] - pressureDropData[i]);
		}
	}
	return -1.0; // 未找到合适的区间
}

nmDataDiagnostic::FlowRegimeAnalysis nmDataDiagnostic::analyzeRadialFlowForKh(
	const QVector<double>& timeData,
	const QVector<double>& derivativeData,
	int startPoint, int pointCount,
	nmDataWellBase* pWell, nmDataReservoir* pReservoir)
{
	FlowRegimeAnalysis analysis;
	analysis.isValid = false;

	DEBUG_OUT(QString("Starting radial flow analysis from point %1, analyzing %2 points")
		.arg(startPoint).arg(pointCount));

	// 寻找导数平台段
	QPair<int, int> bestPlatform = findDerivativePlatform(timeData, derivativeData, 
		startPoint, pointCount);

	if (bestPlatform.first >= 0 && bestPlatform.second >= 0) {
		int platformLength = bestPlatform.second - bestPlatform.first + 1;

		// 计算整个平台段的平均值，而不是只使用两个端点
		double totalDeriv = 0.0;
		double totalTime = 0.0;
		for (int i = bestPlatform.first; i <= bestPlatform.second; ++i) {
			totalDeriv += derivativeData[i];
			totalTime += timeData[i];
		}

		analysis.avgPressure = totalDeriv / platformLength;  // 平台段平均值
		analysis.avgTime = totalTime / platformLength;       // 平台段时间平均值
		analysis.pointPair = bestPlatform;

		// 获取对应时间的流量
		double flowRate = pWell->getFlowRateAtTime(analysis.avgTime);
		double bo = pReservoir->getBo().getValue().toDouble();
		double miuo = pReservoir->getMiuo().getValue().toDouble();

		DEBUG_OUT(QString("Radial flow platform found: points [%1-%2], length=%3 points")
			.arg(bestPlatform.first).arg(bestPlatform.second).arg(platformLength));
		DEBUG_OUT(QString("Time range: %1hr to %2hr")
			.arg(timeData[bestPlatform.first], 0, 'f', 4)
			.arg(timeData[bestPlatform.second], 0, 'f', 4));
		DEBUG_OUT(QString("Platform average derivative: %1 (from %2 points)")
			.arg(analysis.avgPressure, 0, 'f', 6).arg(platformLength));
		DEBUG_OUT(QString("Analysis time: %1 hr, flow rate: %2")
			.arg(analysis.avgTime, 0, 'f', 4).arg(flowRate, 0, 'f', 4));

		if (flowRate > 0 && analysis.avgPressure > 0) {
			// 公式: kh = kh = (70.6 * q * B * μ) / (ΔP'_{platform})
			double kh = flowRate * bo * miuo * 70.6 / analysis.avgPressure;

			analysis.result = kh;
			analysis.isValid = true;

			// 计算平台段的整体斜率（使用首尾点）
			double t1 = timeData[bestPlatform.first];
			double t2 = timeData[bestPlatform.second];
			double deriv1 = derivativeData[bestPlatform.first];
			double deriv2 = derivativeData[bestPlatform.second];

			analysis.slope = calculateSlopeBetweenPoints(t1, deriv1, t2, deriv2);
			analysis.slopeDifference = qAbs(analysis.slope - 0.0);

			analysis.description = QString("Radial Flow Platform (length=%1, DertaP'=%2, kh=%3)")
				.arg(platformLength)
				.arg(analysis.avgPressure, 0, 'f', 6)
				.arg(kh, 0, 'f', 6);

			DEBUG_OUT(QString("Platform overall slope: %1 (over %2 hours)")
				.arg(analysis.slope, 0, 'f', 4)
				.arg(t2 - t1, 0, 'f', 2));
			DEBUG_OUT(QString("Calculated kh: %1").arg(kh, 0, 'f', 6));
		}
	}

	return analysis;
}

QPair<int, int> nmDataDiagnostic::findDerivativePlatform(
	const QVector<double>& timeData,
	const QVector<double>& derivativeData, 
	int startPoint, int pointCount)
{
	QPair<int, int> bestSegment(-1, -1);
	if (startPoint < 0 || pointCount <= 0) {
		DEBUG_OUT("Not enough points for platform identification");
		return bestSegment;
	}

	int endPoint = qMin(startPoint + pointCount, qMin(timeData.size(), derivativeData.size()));

	if (endPoint - startPoint < 3) {
		DEBUG_OUT("Not enough points for platform identification");
		return bestSegment;
	}

	DEBUG_OUT(QString("Searching for longest stable platform in points [%1-%2]").arg(startPoint).arg(endPoint-1));

	// 首先计算整个分析区间的导数平均值，作为平台参考值
	double totalDeriv = 0.0;
	int count = 0;
	for (int i = startPoint; i < endPoint; ++i) {
		totalDeriv += derivativeData[i];
		count++;
	}
	if (count <= 0) {
		return bestSegment;
	}
	double averageDeriv = totalDeriv / count;

	DEBUG_OUT(QString("Overall derivative average in analysis range: %1").arg(averageDeriv, 0, 'f', 6));

	// 定义平台稳定性阈值（相对于平均值的允许偏差）
	const double stabilityThreshold = 0.01; // 1%的偏差
	double maxAllowedDeviation = averageDeriv * stabilityThreshold;

	DEBUG_OUT(QString("Stability threshold: +-%1").arg(maxAllowedDeviation, 0, 'f', 6));

	// 寻找最长的连续稳定段
	int maxLength = 0;
	int currentStart = -1;
	int currentLength = 0;

	for (int i = startPoint; i < endPoint; ++i) {
		double deviation = qAbs(derivativeData[i] - averageDeriv);

		if (deviation <= maxAllowedDeviation) {
			// 点在稳定范围内
			if (currentStart == -1) {
				currentStart = i;
			}
			currentLength++;
		} else {
			// 点超出稳定范围，结束当前段
			if (currentLength > maxLength) {
				maxLength = currentLength;
				bestSegment = qMakePair(currentStart, currentStart + currentLength - 1);
			}
			currentStart = -1;
			currentLength = 0;
		}
	}

	// 检查最后一段
	if (currentLength > maxLength) {
		maxLength = currentLength;
		bestSegment = qMakePair(currentStart, currentStart + currentLength - 1);
	}

	if (bestSegment.first >= 0) {
		DEBUG_OUT(QString("Longest stable platform found: points [%1-%2], length=%3 points")
			.arg(bestSegment.first).arg(bestSegment.second).arg(maxLength));
		DEBUG_OUT(QString("Time range: %1hr to %2hr")
			.arg(timeData[bestSegment.first], 0, 'f', 4)
			.arg(timeData[bestSegment.second], 0, 'f', 4));

		// 计算该平台段的统计信息
		double platformAvg = 0.0;
		for (int i = bestSegment.first; i <= bestSegment.second; ++i) {
			platformAvg += derivativeData[i];
		}
		platformAvg /= maxLength;

		DEBUG_OUT(QString("Platform average derivative: %1, overall average: %2")
			.arg(platformAvg, 0, 'f', 6).arg(averageDeriv, 0, 'f', 6));
	} else {
		DEBUG_OUT("No suitable platform segment found using stability criteria");
	}

	return bestSegment;
}

// 方法1：标准径向流平台分析（保持原有逻辑）
nmDataDiagnostic::FlowRegimeAnalysis nmDataDiagnostic::analyzeStandardRadialFlow(
	const QVector<double>& timeData,
	const QVector<double>& derivativeData,
	int availablePoints,
	nmDataWellBase* pWell, nmDataReservoir* pReservoir)
{
	DEBUG_OUT("=== Method 1: Standard Radial Flow Platform Analysis ===");

	// 使用中后期数据寻找平台
	int startPoint = qMin(static_cast<int>(availablePoints * 0.3), availablePoints - 20);
	int pointsToAnalyze = availablePoints - startPoint;

	if (pointsToAnalyze < 5) {
		FlowRegimeAnalysis failedAnalysis;
		failedAnalysis.isValid = false;
		failedAnalysis.description = "Standard Radial Flow (failed - insufficient points)";
		return failedAnalysis;
	}

	return analyzeRadialFlowForKh(timeData, derivativeData, startPoint, pointsToAnalyze, pWell, pReservoir);
}

// 方法2：晚期分析
nmDataDiagnostic::FlowRegimeAnalysis nmDataDiagnostic::analyzeSaphirStyleLateTime(
	const QVector<double>& timeData,
	const QVector<double>& derivativeData,
	int availablePoints,
	nmDataWellBase* pWell, nmDataReservoir* pReservoir)
{
	DEBUG_OUT("=== Method 2: Late-Time Analysis ===");

	FlowRegimeAnalysis analysis;
	analysis.isValid = false;
	analysis.description = "Late-Time Analysis";

	// 使用后50%的数据寻找极小导数值
	int startPoint = qMax(static_cast<int>(availablePoints * 0.5), availablePoints - 50);
	int endPoint = availablePoints - 1;

	if (endPoint - startPoint < 10) {
		DEBUG_OUT("Not enough late-time data for analysis");
		return analysis;
	}

	DEBUG_OUT(QString("Analyzing late-time points [%1-%2] for minimum derivative")
		.arg(startPoint).arg(endPoint));

	// 寻找晚期导数的最小值段
	QPair<int, int> minDerivativeSegment = findMinimumDerivativeSegment(
		timeData, derivativeData, startPoint, endPoint);

	if (minDerivativeSegment.first >= 0) {
		// 计算该段的平均值
		int segmentStart = minDerivativeSegment.first;
		int segmentEnd = minDerivativeSegment.second;
		int segmentLength = segmentEnd - segmentStart + 1;

		double totalDeriv = 0.0;
		double totalTime = 0.0;
		for (int i = segmentStart; i <= segmentEnd; ++i) {
			totalDeriv += derivativeData[i];
			totalTime += timeData[i];
		}

		analysis.avgPressure = totalDeriv / segmentLength;  // 最小导数段平均值
		analysis.avgTime = totalTime / segmentLength;       // 对应时间平均值
		analysis.pointPair = minDerivativeSegment;

		// 获取流量和储层参数
		double flowRate = pWell->getFlowRateAtTime(analysis.avgTime);
		double bo = pReservoir->getBo().getValue().toDouble();
		double miuo = pReservoir->getMiuo().getValue().toDouble();

		DEBUG_OUT(QString("Minimum derivative segment found: points [%1-%2], length=%3")
			.arg(segmentStart).arg(segmentEnd).arg(segmentLength));
		DEBUG_OUT(QString("Time range: %1hr to %2hr")
			.arg(timeData[segmentStart], 0, 'f', 2).arg(timeData[segmentEnd], 0, 'f', 2));
		DEBUG_OUT(QString("Minimum derivative average: %1 MPa")
			.arg(analysis.avgPressure, 0, 'f', 6));
		DEBUG_OUT(QString("Flow rate at analysis time: %1")
			.arg(flowRate, 0, 'f', 4));

		if (flowRate > 0 && analysis.avgPressure > 0) {
			// kh = (70.6 * q * B * μ) / (DertaP'_min)
			analysis.result = (70.6 * flowRate * bo * miuo) / analysis.avgPressure;
			analysis.isValid = true;

			// 计算质量指标
			analysis.slope = 0.0;  // 晚期分析不关注斜率
			analysis.slopeDifference = 0.0;

			analysis.description = QString("Late-Time (segment=%1pts, DertaP'=%2, kh=%3)")
				.arg(segmentLength)
				.arg(analysis.avgPressure, 0, 'f', 6)
				.arg(analysis.result, 0, 'f', 0);

			DEBUG_OUT(QString("Calculation successful:"));
			DEBUG_OUT(QString("Formula: kh = (70.6 * %1 * %2 * %3) / %4")
				.arg(flowRate, 0, 'f', 1).arg(bo, 0, 'f', 2)
				.arg(miuo, 0, 'f', 2).arg(analysis.avgPressure, 0, 'f', 6));
			DEBUG_OUT(QString("Result: kh = %1 md*m").arg(analysis.result, 0, 'f', 0));
		} else {
			DEBUG_OUT("Invalid flow rate or derivative value for analysis");
		}
	} else {
		DEBUG_OUT("Failed to find suitable minimum derivative segment");
	}

	return analysis;
}

// 方法3：早中期稳定段分析
nmDataDiagnostic::FlowRegimeAnalysis nmDataDiagnostic::analyzeEarlyStableFlow(
	const QVector<double>& timeData,
	const QVector<double>& derivativeData,
	int availablePoints,
	nmDataWellBase* pWell, nmDataReservoir* pReservoir)
{
	DEBUG_OUT("=== Method 3: Early-Middle Stable Flow Analysis ===");

	FlowRegimeAnalysis analysis;
	analysis.isValid = false;
	analysis.description = "Early-Middle Stable Flow";

	// 使用10%-60%的数据寻找相对稳定段
	int startPoint = qMax(static_cast<int>(availablePoints * 0.1), 5);
	int endPoint = qMin(static_cast<int>(availablePoints * 0.6), availablePoints - 5);

	if (endPoint - startPoint < 10) {
		DEBUG_OUT("Not enough data for early-middle stable flow analysis");
		return analysis;
	}

	DEBUG_OUT(QString("Searching for stable flow in points [%1-%2]")
		.arg(startPoint).arg(endPoint));

	// 寻找变异系数最小的连续段
	QPair<int, int> mostStableSegment = findMostStableDerivativeSegment(
		derivativeData, startPoint, endPoint);

	if (mostStableSegment.first >= 0) {
		int segmentStart = mostStableSegment.first;
		int segmentEnd = mostStableSegment.second;
		int segmentLength = segmentEnd - segmentStart + 1;

		// 计算该段的平均值
		double totalDeriv = 0.0;
		double totalTime = 0.0;
		for (int i = segmentStart; i <= segmentEnd; ++i) {
			totalDeriv += derivativeData[i];
			totalTime += timeData[i];
		}

		analysis.avgPressure = totalDeriv / segmentLength;
		analysis.avgTime = totalTime / segmentLength;
		analysis.pointPair = mostStableSegment;

		// 获取流量和储层参数
		double flowRate = pWell->getFlowRateAtTime(analysis.avgTime);
		double bo = pReservoir->getBo().getValue().toDouble();
		double miuo = pReservoir->getMiuo().getValue().toDouble();

		if (flowRate > 0 && analysis.avgPressure > 0) {
			analysis.result = (70.6 * flowRate * bo * miuo) / analysis.avgPressure;
			analysis.isValid = true;

			analysis.description = QString("Early-Middle Stable (segment=%1pts, DertaP'=%2, kh=%3)")
				.arg(segmentLength)
				.arg(analysis.avgPressure, 0, 'f', 6)
				.arg(analysis.result, 0, 'f', 0);

			DEBUG_OUT(QString("Early-middle stable flow found: points [%1-%2], kh=%3")
				.arg(segmentStart).arg(segmentEnd).arg(analysis.result, 0, 'f', 0));
		}
	}

	return analysis;
}

// 寻找最小导数段
QPair<int, int> nmDataDiagnostic::findMinimumDerivativeSegment(
	const QVector<double>& timeData,
	const QVector<double>& derivativeData,
	int startPoint, int endPoint)
{
	QPair<int, int> bestSegment(-1, -1);

	// 寻找连续的低值段（至少5个点）
	const int minSegmentLength = 5;
	double minAvgDerivative = 999999.0;

	for (int i = startPoint; i <= endPoint - minSegmentLength; ++i) {
		for (int len = minSegmentLength; len <= qMin(20, endPoint - i + 1); ++len) {
			int j = i + len - 1;

			// 计算这个段的平均导数
			double totalDeriv = 0.0;
			for (int k = i; k <= j; ++k) {
				totalDeriv += derivativeData[k];
			}
			double avgDeriv = totalDeriv / len;

			// 如果这个段的平均导数更小，更新最佳结果
			if (avgDeriv < minAvgDerivative) {
				minAvgDerivative = avgDeriv;
				bestSegment = qMakePair(i, j);
			}
		}
	}

	if (bestSegment.first >= 0) {
		DEBUG_OUT(QString("Found minimum derivative segment [%1-%2] with average DertaP'=%3")
			.arg(bestSegment.first).arg(bestSegment.second)
			.arg(minAvgDerivative, 0, 'f', 6));
	}

	return bestSegment;
}

// 寻找最稳定的导数段
QPair<int, int> nmDataDiagnostic::findMostStableDerivativeSegment(
	const QVector<double>& derivativeData,
	int startPoint, int endPoint)
{
	QPair<int, int> bestSegment(-1, -1);

	const int minSegmentLength = 8;
	double minVariationCoeff = 999999.0;

	for (int i = startPoint; i <= endPoint - minSegmentLength; ++i) {
		for (int len = minSegmentLength; len <= qMin(25, endPoint - i + 1); ++len) {
			int j = i + len - 1;

			// 计算变异系数 (标准差/均值)
			double variationCoeff = calculateVariationCoefficient(derivativeData, i, j);

			if (variationCoeff < minVariationCoeff) {
				minVariationCoeff = variationCoeff;
				bestSegment = qMakePair(i, j);
			}
		}
	}

	return bestSegment;
}

// 辅助函数：计算变异系数
double nmDataDiagnostic::calculateVariationCoefficient(
	const QVector<double>& data, int start, int end)
{
	if (end <= start) return 999999.0;

	int count = end - start + 1;

	// 计算均值
	double mean = 0.0;
	for (int i = start; i <= end; ++i) {
		mean += data[i];
	}
	mean /= count;

	if (mean <= 0) return 999999.0;

	// 计算标准差
	double variance = 0.0;
	for (int i = start; i <= end; ++i) {
		double diff = data[i] - mean;
		variance += diff * diff;
	}
	double stdDev = qSqrt(variance / count);

	return stdDev / mean;  // 变异系数
}

// 智能选择最佳分析方法
nmDataDiagnostic::FlowRegimeAnalysis* nmDataDiagnostic::selectBestTransmissibilityAnalysis(
	FlowRegimeAnalysis& standardRadial, FlowRegimeAnalysis& saphirStyle, 
	FlowRegimeAnalysis& earlyStable)
{
	QVector<FlowRegimeAnalysis*> validMethods;

	if (standardRadial.isValid) validMethods.append(&standardRadial);
	if (saphirStyle.isValid) validMethods.append(&saphirStyle);
	if (earlyStable.isValid) validMethods.append(&earlyStable);

	if (validMethods.isEmpty()) {
		DEBUG_OUT("No valid transmissibility calculation method found");
		return NULL;
	}

	DEBUG_OUT("=== Method Selection Analysis ===");
	for (int i = 0; i < validMethods.size(); ++i) {
		FlowRegimeAnalysis* method = validMethods[i];
		DEBUG_OUT(QString("Valid method: %1, kh=%2")
			.arg(method->description).arg(method->result, 0, 'f', 0));
	}

	// 选择策略：
	// 1. 优先选择标准径向流（如果存在且质量好）
	// 2. 如果没有标准径向流，选择模拟商业软件行为
	// 3. 最后选择早期稳定流

	if (standardRadial.isValid && standardRadial.slopeDifference < 0.3) {
		DEBUG_OUT("Selected: Standard radial flow (good platform quality)");
		return &standardRadial;
	}

	if (saphirStyle.isValid) {
		DEBUG_OUT("Selected: Saphir-style late-time analysis (following commercial software approach)");
		return &saphirStyle;
	}

	if (earlyStable.isValid) {
		DEBUG_OUT("Selected: Early stable flow (fallback method)");
		return &earlyStable;
	}

	// 兜底：返回第一个有效方法
	DEBUG_OUT(QString("Selected: First available method (%1)").arg(validMethods[0]->description));
	return validMethods[0];
}

double nmDataDiagnostic::calculateSlopeVariation(
	const QVector<double>& timeData,
	const QVector<double>& derivativeData,
	int start, int end)
{
	if (end <= start + 1) return 999999.0;

	double totalSlope = 0.0;
	int slopeCount = 0;

	// 计算段内所有相邻点之间的斜率绝对值
	for (int i = start; i < end; ++i) {
		if (i + 1 >= timeData.size()) break;

		double slope = calculateSlopeBetweenPoints(timeData[i], derivativeData[i],
			timeData[i+1], derivativeData[i+1]);
		totalSlope += qAbs(slope); // 理想平台斜率应该接近0
		slopeCount++;
	}

	return slopeCount > 0 ? totalSlope / slopeCount : 999999.0;
}