nmWTAI-Platform/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp

#include "nmCalculationAutoFitLM.h"
#include "nmCalculationDllPebiSolverTask.h"
#include "nmCalculationUtils.h"
#include "nmDataAnalyzeManager.h"
#include "nmDataWellBase.h"
#include "nmDataVerticalWell.h"
#include "nmDataVerticalFracturedWell.h"
#include "nmDataHorizontalFracturedWell.h"
#include "nmDataReservoir.h"
#include "nmDataAutomaticFitting.h"

#include <QApplication>
#include <QDebug>
#include <QTime>
#include <QDir>
#include <QTextStream>
#include <QFileInfo>
#include <QDateTime>
#include <QtCore/qmath.h>
#include <cmath>
#include <algorithm>
#include <limits>

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

#ifdef QT_NO_DEBUG
static const bool kAutoFitDiagnosticTraceEnabled = false;
#else
static const bool kAutoFitDiagnosticTraceEnabled = true;
#endif

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

// 从嵌套 RMSE 中提取被消除的独立误差贡献。
static double nestedRmsContribution(double reducedModelLoss,
                                    double fullModelLoss)
{
    return qSqrt(qMax(0.0,
                      reducedModelLoss * reducedModelLoss -
                      fullModelLoss * fullModelLoss));
}

static inline void msleep(int ms)
{
#ifdef Q_OS_WIN
    Sleep(ms);
#else
    Q_UNUSED(ms);
#endif
}

static QString csvEscape(const QString& text)
{
    QString escaped = text;
    escaped.replace("\"", "\"\"");
    return QString("\"%1\"").arg(escaped);
}

static QString traceNumber(double value)
{
    return isFiniteNumber(value) ? QString::number(value, 'g', 17) : QString();
}

static QString traceParamAt(const QVector<double>& params, int index)
{
    return (index >= 0 && index < params.size())
           ? traceNumber(params[index])
           : QString();
}

static QString jsonEscape(const QString& text)
{
    QString escaped = text;
    escaped.replace("\\", "\\\\");
    escaped.replace("\"", "\\\"");
    escaped.replace("\b", "\\b");
    escaped.replace("\f", "\\f");
    escaped.replace("\n", "\\n");
    escaped.replace("\r", "\\r");
    escaped.replace("\t", "\\t");
    return QString("\"%1\"").arg(escaped);
}

static QString jsonNumber(double value)
{
    return isFiniteNumber(value) ? QString::number(value, 'g', 17) : QString("null");
}

static QString jsonDoubleArray(const QVector<double>& values)
{
    QStringList items;
    for(int i = 0; i < values.size(); ++i) {
        items << jsonNumber(values[i]);
    }
    return QString("[%1]").arg(items.join(","));
}

static QString jsonIntArray(const QVector<int>& values)
{
    QStringList items;
    for(int i = 0; i < values.size(); ++i) {
        items << QString::number(values[i]);
    }
    return QString("[%1]").arg(items.join(","));
}

static QString jsonBoolArray(const QVector<bool>& values)
{
    QStringList items;
    for(int i = 0; i < values.size(); ++i) {
        items << (values[i] ? "true" : "false");
    }
    return QString("[%1]").arg(items.join(","));
}

static QString jsonStringArray(const QStringList& values)
{
    QStringList items;
    for(int i = 0; i < values.size(); ++i) {
        items << jsonEscape(values[i]);
    }
    return QString("[%1]").arg(items.join(","));
}

// 信赖域搜索统一在 [0, 1] 内部坐标工作。正值参数使用对数坐标，使内部相同步长
// 表示近似相同的相对变化，避免 k、C、Ct、Cf 等跨数量级参数被线性尺度支配；
// skin 可为负数、Swi 的物理意义是线性比例，因此二者保持有界线性坐标。
static bool useTrustRegionLogScale(int parameterIndex, double lower, double upper)
{
    return parameterIndex != 1 && parameterIndex != 7 &&
           lower > 0.0 && upper > lower;
}

static double toTrustRegionCoordinate(double value,
                                      int parameterIndex,
                                      double lower,
                                      double upper)
{
    // 所有进入优化器的物理值先投影到用户上下界，再转换成无量纲坐标。
    // 这样有限差分步长、信赖半径和参数间相关性可以在统一尺度上比较。
    value = qMax(lower, qMin(upper, value));

    if(useTrustRegionLogScale(parameterIndex, lower, upper)) {
        return (qLn(value) - qLn(lower)) / (qLn(upper) - qLn(lower));
    }

    return upper > lower ? (value - lower) / (upper - lower) : 0.0;
}

static double fromTrustRegionCoordinate(double coordinate,
                                        int parameterIndex,
                                        double lower,
                                        double upper)
{
    // 候选内部坐标先限制在 [0,1]，再执行上述映射的逆变换，保证写回
    // DataManager 的参数始终位于用户设置的物理范围内。
    coordinate = qMax(0.0, qMin(1.0, coordinate));

    if(useTrustRegionLogScale(parameterIndex, lower, upper)) {
        return qExp(qLn(lower) + coordinate * (qLn(upper) - qLn(lower)));
    }

    return lower + coordinate * (upper - lower);
}

enum TrustRegionErrorComponent
{
    TRUST_REGION_VERTICAL_COMPONENT = 0,
    TRUST_REGION_HORIZONTAL_COMPONENT,
    TRUST_REGION_SHAPE_COMPONENT,
    TRUST_REGION_TOTAL_COMPONENT
};

// 一次真实求解的完整快照。除了参数和总误差，还保存内部坐标、诊断分量和
// 双对数曲线，因此拒绝候选后可以完整恢复上一个已接受工作点。
struct TrustRegionEvaluation
{
    QVector<double> parameters;
    QVector<double> coordinates;
    AutoFitObjectiveBreakdownLM breakdown;
    QVector<QVector<double> > curve;
    double fitness;
    int elapsedMs;
    bool valid;

    TrustRegionEvaluation()
        : fitness(1.0e10)
        , elapsedMs(-1)
        , valid(false)
    {}
};

// LM 只使用固定长度、全部有限的普通残差。
static bool trustRegionResidualsValid(
        const AutoFitObjectiveBreakdownLM& breakdown)
{
    // 损失函数固定使用 80 个压力点和 80 个导数点。严格校验长度，避免
    // Jacobian 沿用旧维度后访问另一候选的短残差向量。
    if(!breakdown.valid || breakdown.residualVector.size() != 160) {
        return false;
    }

    for(int i = 0; i < breakdown.residualVector.size(); ++i) {
        if(!isFiniteNumber(breakdown.residualVector[i])) {
            return false;
        }
    }
    return true;
}

// 计算向量二范数的平方，避免在只比较能量或计算正规方程时反复开方。
static double trustRegionSquaredNorm(const QVector<double>& values)
{
    double sum = 0.0;
    for(int i = 0; i < values.size(); ++i) {
        sum += values[i] * values[i];
    }
    return sum;
}

// 计算同维向量内积；维度不一致表示局部模型无效，返回零让调用方放弃修正。
static double trustRegionDotProduct(const QVector<double>& left,
                                    const QVector<double>& right)
{
    if(left.size() != right.size()) {
        return 0.0;
    }

    double sum = 0.0;
    for(int i = 0; i < left.size(); ++i) {
        sum += left[i] * right[i];
    }
    return sum;
}

// trace 和运行日志使用稳定的英文标识，便于现有离线脚本继续按字段筛选。
static QString trustRegionComponentName(int component)
{
    if(component == TRUST_REGION_VERTICAL_COMPONENT) {
        return "vertical";
    }
    if(component == TRUST_REGION_HORIZONTAL_COMPONENT) {
        return "horizontal";
    }
    if(component == TRUST_REGION_SHAPE_COMPONENT) {
        return "shape";
    }
    return "total";
}

// 三类损失量纲一致，直接选择当前最大的可靠分量；都很小时退回总残差梯度。
static int trustRegionDominantComponent(
        const AutoFitObjectiveBreakdownLM& breakdown,
        double diagnosisThreshold)
{
    int component = TRUST_REGION_TOTAL_COMPONENT;
    double largestLoss = diagnosisThreshold;

    if(breakdown.verticalReliable &&
            isFiniteNumber(breakdown.verticalLoss) &&
            breakdown.verticalLoss > largestLoss) {
        component = TRUST_REGION_VERTICAL_COMPONENT;
        largestLoss = breakdown.verticalLoss;
    }
    if(breakdown.horizontalReliable &&
            !breakdown.registrationAmbiguous &&
            isFiniteNumber(breakdown.horizontalLoss) &&
            breakdown.horizontalLoss > largestLoss) {
        component = TRUST_REGION_HORIZONTAL_COMPONENT;
        largestLoss = breakdown.horizontalLoss;
    }
    if(isFiniteNumber(breakdown.shapeLoss) &&
            breakdown.shapeLoss > largestLoss) {
        component = TRUST_REGION_SHAPE_COMPONENT;
    }

    return component;
}

// 求解选中参数对应的阻尼正规方程。上下和左右诊断量保留方向；形状没有
// 天然正负，因此使用 shapeLoss 对参数的局部导数。参数最多八维，使用带
// 部分主元的高斯消元即可处理该小矩阵，并在主元退化时明确返回失败。
static bool solveTrustRegionLinearSystem(
        QVector<QVector<double> > matrix,
        QVector<double> rightHandSide,
        QVector<double>* solution)
{
    if(!solution || matrix.isEmpty() ||
            matrix.size() != rightHandSide.size()) {
        return false;
    }

    const int size = matrix.size();
    for(int i = 0; i < size; ++i) {
        if(matrix[i].size() != size) {
            return false;
        }
    }

    for(int column = 0; column < size; ++column) {
        int pivotRow = column;
        double pivotMagnitude = qAbs(matrix[column][column]);
        for(int row = column + 1; row < size; ++row) {
            double magnitude = qAbs(matrix[row][column]);
            if(magnitude > pivotMagnitude) {
                pivotMagnitude = magnitude;
                pivotRow = row;
            }
        }
        if(pivotMagnitude <= 1.0e-14) {
            return false;
        }

        if(pivotRow != column) {
            qSwap(matrix[pivotRow], matrix[column]);
            qSwap(rightHandSide[pivotRow], rightHandSide[column]);
        }

        for(int row = column + 1; row < size; ++row) {
            double factor = matrix[row][column] /
                            matrix[column][column];
            matrix[row][column] = 0.0;
            for(int nextColumn = column + 1;
                    nextColumn < size; ++nextColumn) {
                matrix[row][nextColumn] -=
                    factor * matrix[column][nextColumn];
            }
            rightHandSide[row] -= factor * rightHandSide[column];
        }
    }

    solution->fill(0.0, size);
    for(int row = size - 1; row >= 0; --row) {
        double value = rightHandSide[row];
        for(int column = row + 1; column < size; ++column) {
            value -= matrix[row][column] * (*solution)[column];
        }
        double pivot = matrix[row][row];
        if(qAbs(pivot) <= 1.0e-14) {
            return false;
        }
        (*solution)[row] = value / pivot;
        if(!isFiniteNumber((*solution)[row])) {
            return false;
        }
    }
    return true;
}

// 计算两个 Jacobian 列向量的绝对余弦相似度。接近 1 表示两个参数在当前
// 工作点对曲线的影响几乎相同，联合调整容易产生不可辨识方向。
static double trustRegionJacobianColumnCorrelation(
        const QVector<QVector<double> >& jacobian,
        int leftColumn,
        int rightColumn)
{
    double product = 0.0;
    double leftNorm = 0.0;
    double rightNorm = 0.0;
    for(int row = 0; row < jacobian.size(); ++row) {
        if(leftColumn >= jacobian[row].size() ||
                rightColumn >= jacobian[row].size()) {
            return 1.0;
        }
        double left = jacobian[row][leftColumn];
        double right = jacobian[row][rightColumn];
        product += left * right;
        leftNorm += left * left;
        rightNorm += right * right;
    }

    if(leftNorm <= 1.0e-20 || rightNorm <= 1.0e-20) {
        return 0.0;
    }
    return qAbs(product) / qSqrt(leftNorm * rightNorm);
}

// 每次接受一个真实候选后，使用满足最新割线条件的秩一修正更新完整残差
// Jacobian。这样模型吸收了刚得到的真实变化，又不必立即逐参数重新试算。
static void updateTrustRegionJacobian(
        QVector<QVector<double> >* jacobian,
        const QVector<double>& oldResidual,
        const QVector<double>& newResidual,
        const QVector<double>& coordinateStep)
{
    if(!jacobian || jacobian->size() != oldResidual.size() ||
            oldResidual.size() != newResidual.size()) {
        return;
    }

    double denominator = trustRegionSquaredNorm(coordinateStep);
    if(denominator <= 1.0e-12) {
        return;
    }

    for(int row = 0; row < jacobian->size(); ++row) {
        if((*jacobian)[row].size() != coordinateStep.size()) {
            return;
        }

        double predictedChange = 0.0;
        for(int column = 0; column < coordinateStep.size(); ++column) {
            predictedChange +=
                (*jacobian)[row][column] * coordinateStep[column];
        }
        double correction =
            (newResidual[row] - oldResidual[row] - predictedChange) /
            denominator;
        for(int column = 0; column < coordinateStep.size(); ++column) {
            (*jacobian)[row][column] +=
                correction * coordinateStep[column];
        }
    }
}

// 对上下偏差、左右偏差和形状损失的梯度执行同样的割线秩一修正，使诊断
// 选参模型与完整残差 Jacobian 保持在同一个已接受工作点。
static void updateTrustRegionScalarGradient(
        QVector<double>* gradient,
        double oldValue,
        double newValue,
        const QVector<double>& coordinateStep)
{
    if(!gradient || gradient->size() != coordinateStep.size() ||
            !isFiniteNumber(oldValue) || !isFiniteNumber(newValue)) {
        return;
    }

    double denominator = trustRegionSquaredNorm(coordinateStep);
    if(denominator <= 1.0e-12) {
        return;
    }

    double predictedChange = trustRegionDotProduct(
        *gradient, coordinateStep);
    double correction =
        (newValue - oldValue - predictedChange) / denominator;
    for(int i = 0; i < gradient->size(); ++i) {
        (*gradient)[i] += correction * coordinateStep[i];
    }
}

static QStringList traceParameterNames()
{
    QStringList names;
    names << "k"
          << "skin"
          << "wellboreC"
          << "phi"
          << "h"
          << "Ct"
          << "Cf"
          << "Swi"
          << "Dfc"
          << "fractureHalfLength";
    return names;
}

nmCalculationAutoFitLM::nmCalculationAutoFitLM(QObject* parent)
    : QObject(parent)
    , m_isRunning(false)
    , m_shouldStop(false)
    , m_currentIteration(0)
    , m_globalBestFitness(1e10)
    , m_maxIterations(100)
    , m_targetError(0.001)
    , m_totalEvaluations(0)
    , m_successfulEvaluations(0)
    , m_evaluationInProgress(0)
    , m_consecutiveFailures(0)
    , m_userInitialFitness(1e10)
    , m_maxConsecutiveFailures(3)
    , m_hasValidUserSolution(false)
    , m_targetWellName("")
    , m_traceRunId("")
    , m_traceFilePath("")
    , m_traceMetaFilePath("")
{
    // 单次运行目录在输入校验通过后创建，避免只打开界面也产生临时文件。
    DEBUG_OUT("LM automatic fitting calculator initialized");
}

// 析构函数：停止仍在进行的拟合、关闭 trace 文件并清理临时目录。
// 自动拟合可能在 UI 线程中被窗口关闭打断，因此析构时要尽量温和地等待当前评价结束；
// 如果等待超时，再强制清除运行标志，避免对象销毁后还有信号回调访问成员变量。
nmCalculationAutoFitLM::~nmCalculationAutoFitLM()
{
    if(m_isRunning) {
        m_shouldStop = true;
        int waitCount = 0;
        while(m_isRunning && waitCount < 100) {
            QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, 50);
            msleep(50);
            waitCount++;
        }
        if(m_isRunning) {
            DEBUG_OUT("Force stopping LM fitting after timeout");
            m_isRunning = false;
        }
    }

    closeTraceFile();
    cleanupTemporaryDirectory();
    disconnect(this, nullptr, nullptr, nullptr);
}

// ===== 临时目录工具 =====
//
// 真实求解器和井曲线 CSV 共用一个受控的单次运行目录。创建失败时不允许
// 回退到程序目录，否则后续递归清理可能删除安装文件。
bool nmCalculationAutoFitLM::initializeTemporaryDirectory()
{
    nmCalculationUtils::cleanupStaleAutoFitTemporaryDirectories();
    m_tempDirectory = nmCalculationUtils::createAutoFitTemporaryDirectory();
    if(m_tempDirectory.isEmpty()) {
        DEBUG_OUT("Failed to create LM auto-fit temporary directory");
        return false;
    }

    m_traceRunId = QFileInfo(m_tempDirectory).fileName();
    DEBUG_OUT(QString("Initialized temp directory: %1").arg(m_tempDirectory));
    return true;
}

void nmCalculationAutoFitLM::cleanupTemporaryDirectory()
{
    if(m_tempDirectory.isEmpty()) {
        return;
    }
    const QString directoryPath = m_tempDirectory;
    m_tempDirectory.clear();
    if(nmCalculationUtils::removeAutoFitTemporaryDirectory(directoryPath)) {
        DEBUG_OUT("Temp directory cleaned up successfully");
    } else {
        DEBUG_OUT(QString("Failed to clean auto-fit temporary directory: %1")
                  .arg(directoryPath));
    }
}


// ==================== 公共接口方法 ====================

void nmCalculationAutoFitLM::setTargetLogLogData(const QVector<QVector<double> >& targetData)
{
    // 目标曲线由界面层从目标井 history log-log 传入。
    // 约定 targetData[0]=time，targetData[1]=pressure，targetData[2]=pressure derivative。
    m_targetLogLogData = targetData;
    DEBUG_OUT(QString("Target LogLog data set: %1 arrays").arg(targetData.size()));

    if(targetData.size() >= 3) {
        DEBUG_OUT(QString("LogLog data points: X=%1, Y1=%2, Y2=%3")
                  .arg(targetData[0].size())
                  .arg(targetData[1].size())
                  .arg(targetData[2].size()));
    }
}


void nmCalculationAutoFitLM::stopFitting()
{
    // 用户点击停止时只设置请求标志，让 LM 主循环和求解器等待逻辑自然退出。
    if(m_isRunning) {
        emit logMessageGenerated(tr("=== User Stop Request Received ==="));
        emit logMessageGenerated(tr("Gracefully stopping LM automatic fitting..."));
        m_shouldStop = true;

        // 给当前评价一个短暂的自然退出时间。若仍在运行，
        // runSolverDll() 会在下一个等待周期检查 m_shouldStop 并结束任务。
        int waitCount = 0;

        while(m_evaluationInProgress > 0 && waitCount < 30) {
            QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, 50);
            msleep(50);
            waitCount++;
        }

        if(m_evaluationInProgress > 0) {
            emit logMessageGenerated(tr("Waiting for current solver evaluation to stop..."));
        }

        emit logMessageGenerated(tr("LM automatic fitting stop request processed"));
    } else {
        emit logMessageGenerated(tr("Stop request received but optimization is not running"));
    }
}

bool nmCalculationAutoFitLM::isRunning() const
{
    // UI 查询当前是否处于自动拟合运行状态。
    return m_isRunning;
}

int nmCalculationAutoFitLM::getCurrentIteration() const
{
    // UI 进度条和日志展示用的当前迭代序号。
    return m_currentIteration;
}

int nmCalculationAutoFitLM::getTotalEvaluations() const
{
    // 返回本轮拟合实际调用真实求解器的总次数，供完成日志展示。
    return m_totalEvaluations;
}

QVector<double> nmCalculationAutoFitLM::getBestSolution() const
{
    // 返回紧凑的“启用参数向量”，顺序与 m_enabledParamIndices 一致。
    return m_globalBestPosition;
}

double nmCalculationAutoFitLM::getBestFitness() const
{
    // 当前全局最优真实误差。越小越好，1e10 附近通常表示尚无有效解。
    return m_globalBestFitness;
}

AutoFitObjectiveBreakdownLM nmCalculationAutoFitLM::getLastObjectiveBreakdown() const
{
    // 返回最近一次损失评价的误差分解，供界面或后续优化逻辑读取。
    return m_lastObjectiveBreakdown;
}

QString nmCalculationAutoFitLM::getLastError() const
{
    // 上一次失败的人类可读错误信息，主要给 UI 层弹窗或日志使用。
    return m_lastError;
}

void nmCalculationAutoFitLM::resetOptimizer()
{
    // 清空一次运行产生的状态，但不销毁对象。
    // 配置字段会在 startAutoFitting() 中重新从 DataManager 读取；
    // trace 文件先关闭，避免新一轮 run 继续写到旧 CSV。
    closeTraceFile();
    m_globalBestPosition.clear();
    m_globalBestFitness = 1e10;
    m_globalBestObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
    m_lastEvaluatedLogLogData.clear();
    m_globalBestLogLogData.clear();
    m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
    m_userInitialLogLogData.clear();
    m_userInitialObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
    m_currentIteration = 0;
    m_totalEvaluations = 0;
    m_successfulEvaluations = 0;
    m_lastError.clear();
    m_initialValues.clear();
    m_userInitialSolution.clear();
    m_userInitialFitness = 1e10;
    m_hasValidUserSolution = false;
    m_traceMetaFilePath.clear();

    DEBUG_OUT("LM optimizer reset");
}

void nmCalculationAutoFitLM::setTargetWellName(const QString& wellName)
{
    // 目标井名是贯穿拟合流程的关键索引：
    // 读目标曲线、写 skin/wellboreC、求解后取 resultLogLog 都依赖这个名字。
    m_targetWellName = wellName;
}

void nmCalculationAutoFitLM::initializeTraceFile()
{
    closeTraceFile();
    m_traceFilePath.clear();
    m_traceMetaFilePath.clear();
    if(!kAutoFitDiagnosticTraceEnabled) {
        return;
    }

    // 诊断轨迹仅在 Debug 构建中启用，并与源码、安装目录完全分离。
    if(m_traceRunId.isEmpty()) {
        m_traceRunId = QString("%1-%2")
            .arg(QString::number(QCoreApplication::applicationPid()))
            .arg(QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss_zzz"));
    }
    QDir traceDir(QDir(QDir::tempPath()).absoluteFilePath("WTAI/AutoFitTrace/LM"));

    if(!traceDir.exists() && !QDir().mkpath(traceDir.absolutePath())) {
        DEBUG_OUT(QString("Failed to create LM trace directory: %1").arg(traceDir.absolutePath()));
        m_traceFilePath.clear();
        m_traceMetaFilePath.clear();
        return;
    }

    m_traceFilePath = traceDir.absoluteFilePath(
        QString("lm_trust_region_trace_%1.csv").arg(m_traceRunId));
    m_traceMetaFilePath = traceDir.absoluteFilePath(
        QString("lm_trust_region_trace_%1.meta.json").arg(m_traceRunId));
    m_traceFile.setFileName(m_traceFilePath);

    if(!m_traceFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
        DEBUG_OUT(QString("Failed to open LM trace file: %1").arg(m_traceFilePath));
        m_traceFilePath.clear();
        m_traceMetaFilePath.clear();
        return;
    }

    writeTraceHeader();
    writeTraceMetaFile();
    emit logMessageGenerated(tr("LM fitting trace: %1").arg(m_traceFilePath));
    if(!m_traceMetaFilePath.isEmpty()) {
        emit logMessageGenerated(
            tr("LM fitting trace metadata: %1").arg(m_traceMetaFilePath));
    }
}

void nmCalculationAutoFitLM::closeTraceFile()
{
    if(m_traceFile.isOpen()) {
        m_traceFile.flush();
        m_traceFile.close();
    }
}

void nmCalculationAutoFitLM::writeTraceHeader()
{
    if(!m_traceFile.isOpen()) {
        return;
    }

    QStringList cols;
    cols << "run_id"
         << "iteration"
         << "parameter_index"
         << "phase"
         << "k"
         << "skin"
         << "wellboreC"
         << "phi"
         << "h"
         << "Ct"
         << "Cf"
         << "Swi"
         << "Dfc"
         << "fractureHalfLength"
         << "solver_objective"
         << "solver_success"
         << "elapsed_ms"
         << "decision"
         << "enabled_param_indices"
         << "pressure_loss"
         << "derivative_loss"
         << "vertical_common_bias"
         << "vertical_loss"
         << "vertical_reliable"
         << "horizontal_physical_shift"
         << "horizontal_loss"
         << "horizontal_reliable"
         << "shape_loss"
         << "late_trend_loss"
         << "late_slope_bias"
         << "late_trend_reliable"
         << "registration_ambiguous";

    QTextStream out(&m_traceFile);
    out << cols.join(",") << "\n";
}

void nmCalculationAutoFitLM::writeTraceMetaFile()
{
    if(m_traceMetaFilePath.isEmpty()) {
        return;
    }

    QFile metaFile(m_traceMetaFilePath);
    if(!metaFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
        DEBUG_OUT(QString("Failed to open LM trace meta file: %1").arg(m_traceMetaFilePath));
        m_traceMetaFilePath.clear();
        return;
    }

    QStringList parameterNames = traceParameterNames();
    QStringList enabledNames;
    for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
        int paramIndex = m_enabledParamIndices[i];
        enabledNames << ((paramIndex >= 0 && paramIndex < parameterNames.size())
                         ? parameterNames[paramIndex]
                         : QString::number(paramIndex));
    }

    QVector<double> initialFullParams = buildTraceParameterVector(m_initialValues);
    QVector<double> targetTime = m_targetLogLogData.size() > 0
                                 ? m_targetLogLogData[0] : QVector<double>();
    QVector<double> targetPressure = m_targetLogLogData.size() > 1
                                     ? m_targetLogLogData[1] : QVector<double>();
    QVector<double> targetDerivative = m_targetLogLogData.size() > 2
                                       ? m_targetLogLogData[2] : QVector<double>();

    QTextStream out(&metaFile);
    out << "{\n";
    out << "  \"schema_version\": 1,\n";
    out << "  \"trace_type\": \"finite_difference_lm_trust_region\",\n";
    out << "  \"run_id\": " << jsonEscape(m_traceRunId) << ",\n";
    out << "  \"created_at\": "
        << jsonEscape(QDateTime::currentDateTime().toString(Qt::ISODate)) << ",\n";
    out << "  \"trace_csv\": " << jsonEscape(QFileInfo(m_traceFilePath).fileName()) << ",\n";
    out << "  \"target\": {\n";
    out << "    \"well_name\": " << jsonEscape(m_targetWellName) << ",\n";
    out << "    \"time\": " << jsonDoubleArray(targetTime) << ",\n";
    out << "    \"pressure\": " << jsonDoubleArray(targetPressure) << ",\n";
    out << "    \"derivative\": " << jsonDoubleArray(targetDerivative) << "\n";
    out << "  },\n";
    out << "  \"lm\": {\n";
    out << "    \"max_iterations\": " << m_maxIterations << ",\n";
    out << "    \"target_error\": " << jsonNumber(m_targetError) << "\n";
    out << "  },\n";
    out << "  \"parameters\": {\n";
    out << "    \"names\": " << jsonStringArray(parameterNames) << ",\n";
    out << "    \"enabled_indices\": " << jsonIntArray(m_enabledParamIndices) << ",\n";
    out << "    \"enabled_names\": " << jsonStringArray(enabledNames) << ",\n";
    out << "    \"selected_flags\": " << jsonBoolArray(m_parameterSelected) << ",\n";
    out << "    \"lower\": " << jsonDoubleArray(m_parameterLower) << ",\n";
    out << "    \"upper\": " << jsonDoubleArray(m_parameterUpper) << ",\n";
    out << "    \"initial_selected\": " << jsonDoubleArray(m_initialValues) << ",\n";
    out << "    \"initial_full\": " << jsonDoubleArray(initialFullParams) << "\n";
    out << "  }\n";
    out << "}\n";
    metaFile.flush();
    metaFile.close();
}

void nmCalculationAutoFitLM::writeTraceRow(
        int iteration,
        int parameterIndex,
        const QString& phase,
        const QVector<double>& parameters,
        double solverObjective,
        bool solverSuccess,
        int elapsedMs,
        const QString& decision,
        const AutoFitObjectiveBreakdownLM* objectiveBreakdown)
{
    if(!m_traceFile.isOpen()) {
        return;
    }

    QVector<double> fullParams = buildTraceParameterVector(parameters);
    QStringList enabledIndices;
    for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
        enabledIndices << QString::number(m_enabledParamIndices[i]);
    }

    QStringList cols;
    cols << csvEscape(m_traceRunId)
         << QString::number(iteration)
         << QString::number(parameterIndex)
         << csvEscape(phase);
    for(int i = 0; i < 10; ++i) {
        cols << traceParamAt(fullParams, i);
    }
    cols << traceNumber(solverObjective)
         << QString::number(solverSuccess ? 1 : 0)
         << QString::number(elapsedMs)
         << csvEscape(decision)
         << csvEscape(enabledIndices.join(";"));

    if(objectiveBreakdown && objectiveBreakdown->valid) {
        cols << traceNumber(objectiveBreakdown->pressureLoss)
             << traceNumber(objectiveBreakdown->derivativeLoss)
             << traceNumber(objectiveBreakdown->verticalCommonBias)
             << traceNumber(objectiveBreakdown->verticalLoss)
             << QString::number(objectiveBreakdown->verticalReliable ? 1 : 0)
             << traceNumber(objectiveBreakdown->horizontalPhysicalShift)
             << traceNumber(objectiveBreakdown->horizontalLoss)
             << QString::number(objectiveBreakdown->horizontalReliable ? 1 : 0)
             << traceNumber(objectiveBreakdown->shapeLoss)
             << traceNumber(objectiveBreakdown->lateDerivativeTrendLoss)
             << traceNumber(objectiveBreakdown->lateDerivativeSlopeBias)
             << QString::number(objectiveBreakdown->lateDerivativeTrendReliable ? 1 : 0)
             << QString::number(objectiveBreakdown->registrationAmbiguous ? 1 : 0);
    } else {
        for(int i = 0; i < 13; ++i) {
            cols << QString();
        }
    }

    QTextStream out(&m_traceFile);
    out << cols.join(",") << "\n";
    m_traceFile.flush();
}

void nmCalculationAutoFitLM::emitRunSummary(bool success, StopReasonLM finalReason)
{
    // 汇总仅描述 LM 迭代和真实求解器评价。
    emit logMessageGenerated(tr("=== LM Run Summary ==="));
    emit logMessageGenerated(tr("Stop reason: %1").arg(getStopReasonDescription(finalReason)));
    emit logMessageGenerated(
        tr("Result: %1, final error=%2, iterations=%3, evaluations=%4 (successful=%5, failed=%6)")
        .arg(success ? tr("SUCCESS") : tr("FAILED"))
        .arg(m_globalBestFitness, 0, 'e', 4)
        .arg(m_currentIteration + 1)
        .arg(m_totalEvaluations)
        .arg(m_successfulEvaluations)
        .arg(m_totalEvaluations - m_successfulEvaluations));

    if(!m_traceFilePath.isEmpty()) {
        emit logMessageGenerated(tr("Artifacts: trace=%1").arg(m_traceFilePath));
    }
    if(!m_traceMetaFilePath.isEmpty()) {
        emit logMessageGenerated(tr("Artifacts: trace_meta=%1").arg(m_traceMetaFilePath));
    }
}

QVector<double> nmCalculationAutoFitLM::buildTraceParameterVector(const QVector<double>& selectedParameters) const
{
    // 将 LM 内部使用的“启用参数向量”还原成完整 10 维参数向量。
    // 未启用的参数从当前 DataManager 读取，启用的参数用 selectedParameters 覆盖。
    // trace CSV 和 meta 使用该完整向量记录一次候选评价。
    QVector<double> fullParams(10, 0.0);

    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();

    if(dataManager) {
        nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
        fullParams[0] = reservoirData.getPermeability().getValue().toDouble();
        fullParams[3] = reservoirData.getPorosity().getValue().toDouble();
        fullParams[4] = reservoirData.getThickness().getValue().toDouble();
        fullParams[5] = reservoirData.getCt().getValue().toDouble();
        fullParams[6] = reservoirData.getCf().getValue().toDouble();
        fullParams[7] = reservoirData.getSwi().getValue().toDouble();

        nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);

        if(pTargetWell) {
            nmDataPerforation* perforation = pTargetWell->getPerforation(0);
            if(perforation) {
                fullParams[1] = perforation->getSkin().getValue().toDouble();
            }
            fullParams[2] = pTargetWell->getWellboreStorage().getValue().toDouble();

            // Dfc 只存在于两类压裂井，普通井在完整向量中保持为 0。
            if(pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
                nmDataVerticalFracturedWell* fracturedWell =
                    dynamic_cast<nmDataVerticalFracturedWell*>(pTargetWell);
                if(fracturedWell) {
                    fullParams[8] = fracturedWell->getDfc().getValue().toDouble();
                    fullParams[9] = fracturedWell->getFractureHalfLength().getValue().toDouble();
                }
            } else if(pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
                nmDataHorizontalFracturedWell* fracturedWell =
                    dynamic_cast<nmDataHorizontalFracturedWell*>(pTargetWell);
                if(fracturedWell) {
                    fullParams[8] = fracturedWell->getDfc().getValue().toDouble();
                    fullParams[9] = fracturedWell->getFractureHalfLength().getValue().toDouble();
                }
            }
        }
    }

    for(int i = 0; i < selectedParameters.size() && i < m_enabledParamIndices.size(); ++i) {
        int paramIndex = m_enabledParamIndices[i];

        if(paramIndex >= 0 && paramIndex < fullParams.size()) {
            fullParams[paramIndex] = selectedParameters[i];
        }
    }

    return fullParams;
}

// ==================== 数据加载方法 ====================

bool nmCalculationAutoFitLM::loadAllConfigFromDataManager()
{
    // 统一从 DataManager 加载本次运行所需配置。
    // UI 层只负责把用户选择保存到 nmDataAutomaticFitting，本类从这里开始完全数据驱动。
    try {
        loadOptimizationConfig();
        loadParameterBounds();
        extractUserInitialValues(); // 直接提取初始值，无需条件判断
        return true;
    } catch(...) {
        m_lastError = "Failed to load configuration from data manager";
        return false;
    }
}

void nmCalculationAutoFitLM::loadOptimizationConfig()
{
    // LM 只读取迭代次数和目标误差，其他数值控制保留在现有算法实现中。
    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
    nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();

    m_maxIterations = fittingData.getIterationCount().getValue().toInt();
    m_targetError = fittingData.getErrorTolerance().getValue().toDouble();
    DEBUG_OUT(QString("Loaded LM config: iterations=%1, error=%2")
              .arg(m_maxIterations).arg(m_targetError));
}

void nmCalculationAutoFitLM::loadParameterBounds()
{
    // 读取用户勾选的拟合参数及上下界。
    //
    // 这里构建三个核心数组：
    // - m_parameterSelected[10]：完整参数体系中每个参数是否参与拟合；
    // - m_parameterLower/Upper[10]：完整参数体系的搜索上下界；
    // - m_enabledParamIndices：把粒子内部紧凑向量映射回完整参数索引。
    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
    nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();

    // 获取参数选择状态
    m_parameterSelected.resize(10);
    m_parameterSelected[0] = fittingData.getPermeabilitySelected();
    m_parameterSelected[1] = fittingData.getSkinSelected();
    m_parameterSelected[2] = fittingData.getWellboreStorageSelected();
    m_parameterSelected[3] = fittingData.getPorositySelected();
    m_parameterSelected[4] = fittingData.getThicknessSelected();
    m_parameterSelected[5] = fittingData.getCtSelected();
    m_parameterSelected[6] = fittingData.getCfSelected();
    m_parameterSelected[7] = fittingData.getSwiSelected();
    m_parameterSelected[8] = fittingData.getFractureConductivitySelected();
    m_parameterSelected[9] = fittingData.getFractureHalfLengthSelected();

    // 获取参数边界
    m_parameterLower.resize(10);
    m_parameterUpper.resize(10);

    m_parameterLower[0] = fittingData.getPermeabilityMin().getValue().toDouble();
    m_parameterUpper[0] = fittingData.getPermeabilityMax().getValue().toDouble();

    m_parameterLower[1] = fittingData.getSkinMin().getValue().toDouble();
    m_parameterUpper[1] = fittingData.getSkinMax().getValue().toDouble();

    m_parameterLower[2] = fittingData.getWellboreStorageMin().getValue().toDouble();
    m_parameterUpper[2] = fittingData.getWellboreStorageMax().getValue().toDouble();

    m_parameterLower[3] = fittingData.getPorosityMin().getValue().toDouble();
    m_parameterUpper[3] = fittingData.getPorosityMax().getValue().toDouble();

    m_parameterLower[4] = fittingData.getThicknessMin().getValue().toDouble();
    m_parameterUpper[4] = fittingData.getThicknessMax().getValue().toDouble();

    m_parameterLower[5] = fittingData.getCtMin().getValue().toDouble();
    m_parameterUpper[5] = fittingData.getCtMax().getValue().toDouble();

    m_parameterLower[6] = fittingData.getCfMin().getValue().toDouble();
    m_parameterUpper[6] = fittingData.getCfMax().getValue().toDouble();

    m_parameterLower[7] = fittingData.getSwiMin().getValue().toDouble();
    m_parameterUpper[7] = fittingData.getSwiMax().getValue().toDouble();

    m_parameterLower[8] = fittingData.getFractureConductivityMin().getValue().toDouble();
    m_parameterUpper[8] = fittingData.getFractureConductivityMax().getValue().toDouble();

    m_parameterLower[9] = fittingData.getFractureHalfLengthMin().getValue().toDouble();
    m_parameterUpper[9] = fittingData.getFractureHalfLengthMax().getValue().toDouble();

    // 更新启用参数索引
    m_enabledParamIndices.clear();

    for(int i = 0; i < m_parameterSelected.size(); ++i) {
        if(m_parameterSelected[i]) {
            m_enabledParamIndices.append(i);
        }
    }

    DEBUG_OUT(QString("Loaded parameter bounds: %1 enabled parameters")
              .arg(m_enabledParamIndices.size()));
}

// ==================== 自动拟合核心方法 ====================
bool nmCalculationAutoFitLM::startAutoFitting()
{
    // 总入口只负责准备数据、调用有限差分 + LM/信赖域，并写回最终结果。
    StopReasonLM finalReason = LM_CONTINUE_OPTIMIZATION;

    if(m_isRunning) {
        m_lastError = "Auto fitting is already running";
        return false;
    }

    try {
        if(!loadAllConfigFromDataManager()) {
            emit logMessageGenerated(tr("ERROR: Failed to load configuration from data manager"));
            return false;
        }

        emit logMessageGenerated(tr("Algorithm: LM"));
        const int enabledParams = getEnabledParameterCount();
        emit logMessageGenerated(tr("Enabled parameters count: %1").arg(enabledParams));

        if(enabledParams == 0) {
            m_lastError = "No parameters enabled for optimization";
            emit logMessageGenerated(tr("ERROR: No parameters enabled for optimization"));
            return false;
        }

        if(m_targetLogLogData.size() < 3) {
            m_lastError = "Target LogLog data is empty or insufficient";
            emit logMessageGenerated(tr("ERROR: Target LogLog data is empty or insufficient"));
            return false;
        }

        if(m_targetLogLogData[0].size() != m_targetLogLogData[1].size() ||
                m_targetLogLogData[0].size() != m_targetLogLogData[2].size()) {
            m_lastError = "Target LogLog data arrays have inconsistent sizes";
            emit logMessageGenerated(tr("ERROR: Target LogLog data arrays have inconsistent sizes"));
            return false;
        }

        if(m_targetWellName.isEmpty()) {
            m_lastError = "Target well name is empty";
            emit logMessageGenerated(tr("ERROR: Target well name is empty"));
            return false;
        }

        emit logMessageGenerated(
            tr("Target data validation passed (%1 data points)")
            .arg(m_targetLogLogData[0].size()));

        // resetOptimizer() 会清空运行状态，因此先保存从 DataManager 提取的初始值。
        QVector<double> savedInitialValues = m_initialValues;
        resetOptimizer();
        // 输入校验通过后，为本次拟合创建新的临时计算目录。
        cleanupTemporaryDirectory();
        if(!initializeTemporaryDirectory()) {
            m_lastError = tr("Cannot create the automatic fitting temporary directory");
            emit logMessageGenerated(tr("ERROR: %1").arg(m_lastError));
            return false;
        }
        m_isRunning = true;
        m_shouldStop = false;
        m_currentIteration = 0;
        m_consecutiveFailures = 0;
        m_initialValues = savedInitialValues;
        initializeTraceFile();

        // 先用真实求解器评价用户当前模型，供最终精英保护使用。
        if(!savedInitialValues.isEmpty()) {
            m_userInitialSolution = savedInitialValues;
            emit logMessageGenerated(tr("=== Evaluating Initial Solution (Elite Protection) ==="));

            QString paramStr = tr("Initial parameters: ");
            for(int i = 0; i < m_userInitialSolution.size(); ++i) {
                paramStr += QString("[%1]=%2 ")
                            .arg(i).arg(m_userInitialSolution[i], 0, 'f', 6);
            }
            emit logMessageGenerated(paramStr);

            try {
                emit logMessageGenerated(
                    tr("Starting initial solution evaluation..."));
                QTime initialEvalTimer;
                initialEvalTimer.start();
                m_totalEvaluations++;
                m_userInitialFitness = evaluateFitness(m_userInitialSolution);
                const int initialEvalElapsedMs = initialEvalTimer.elapsed();

                if(m_userInitialFitness < 1e9) {
                    m_successfulEvaluations++;
                    m_hasValidUserSolution = true;
                    m_globalBestFitness = m_userInitialFitness;
                    m_globalBestPosition = m_userInitialSolution;
                    m_userInitialLogLogData = m_lastEvaluatedLogLogData;
                    m_globalBestLogLogData = m_userInitialLogLogData;
                    m_userInitialObjectiveBreakdown = m_lastObjectiveBreakdown;
                    m_globalBestObjectiveBreakdown = m_userInitialObjectiveBreakdown;
                    emit logMessageGenerated(tr("Initial solution evaluation successful"));
                    emit logMessageGenerated(
                        tr("Initial Error: %1").arg(m_userInitialFitness, 0, 'e', 4));
                    emit bestCurveUpdated(m_targetLogLogData,
                                          m_globalBestLogLogData,
                                          0,
                                          m_globalBestFitness);
                } else {
                    m_hasValidUserSolution = false;
                    emit logMessageGenerated(tr("Initial solution evaluation failed"));
                }

                writeTraceRow(-1,
                              -1,
                              "initial_solution",
                              m_userInitialSolution,
                              m_userInitialFitness,
                              m_userInitialFitness < 1e9,
                              initialEvalElapsedMs,
                              m_hasValidUserSolution ? "valid" : "invalid",
                              m_hasValidUserSolution
                              ? &m_userInitialObjectiveBreakdown : nullptr);
            } catch(...) {
                m_hasValidUserSolution = false;
                emit logMessageGenerated(tr("Exception during initial solution evaluation"));
            }

            m_initialValues = savedInitialValues;
        }

        finalReason = runTrustRegionFitting();
        validateAndProtectFinalResult();

        if(!m_globalBestPosition.isEmpty() && m_globalBestObjectiveBreakdown.valid) {
            // 精英保护之后记录最终行，保证轨迹与实际写回参数一致。
            writeTraceRow(m_currentIteration,
                          -1,
                          "trust_region_final",
                          m_globalBestPosition,
                          m_globalBestFitness,
                          m_globalBestFitness < 1.0e9,
                          -1,
                          "final_result",
                          &m_globalBestObjectiveBreakdown);
        }

        if(m_globalBestFitness < m_targetError) {
            finalReason = LM_TARGET_ACHIEVED;
        }
    } catch(const std::exception& e) {
        m_lastError = QString(tr("Critical exception in automatic fitting: %1")).arg(e.what());
        emit logMessageGenerated(tr("CRITICAL ERROR: %1").arg(e.what()));
        closeTraceFile();
        cleanupTemporaryDirectory();
        m_isRunning = false;
        emit fittingFinished(false, m_lastError);
        return false;
    } catch(...) {
        m_lastError = tr("Unknown critical exception in automatic fitting");
        emit logMessageGenerated(tr("CRITICAL ERROR: Unknown exception in automatic fitting"));
        closeTraceFile();
        cleanupTemporaryDirectory();
        m_isRunning = false;
        emit fittingFinished(false, m_lastError);
        return false;
    }

    bool finalFullSolverSucceeded = true;
    bool finalFullSolverExecuted = false;

    if(!m_globalBestPosition.isEmpty()) {
        try {
            emit logMessageGenerated(tr("Applying optimized parameters to model..."));
            applyParametersToDataManager(m_globalBestPosition);

            // 裂缝参数会改变网格输入；标记失效后，最终求解任务会基于新快照重建网格。
            const bool fractureGridParameterSelected =
                (m_parameterSelected.size() > 8 && m_parameterSelected[8]) ||
                (m_parameterSelected.size() > 9 && m_parameterSelected[9]);
            if(fractureGridParameterSelected) {
                nmDataAnalyzeManager* dataManager =
                        nmDataAnalyzeManager::getCurrentInstance();
                if(!dataManager) {
                    throw std::runtime_error("Data manager is unavailable");
                }
                dataManager->invalidatePebiGrid();
            }

            if(m_shouldStop) {
                emit logMessageGenerated(
                    tr("Final full-field calculation skipped after user stop"));
            } else {
                emit logMessageGenerated(
                    tr("Running final full-field calculation with optimized parameters..."));
                finalFullSolverExecuted = true;
                finalFullSolverSucceeded = runFinalFullSolver();

                if(finalFullSolverSucceeded) {
                    emit logMessageGenerated(
                        tr("Final full-field calculation completed successfully"));
                } else if(m_shouldStop) {
                    finalFullSolverExecuted = false;
                    finalFullSolverSucceeded = true;
                    emit logMessageGenerated(
                        tr("Final full-field calculation stopped by user"));
                } else {
                    m_lastError =
                        tr("Optimized parameters were found, but the final full-field calculation failed");
                    emit logMessageGenerated(
                        tr("ERROR: Final full-field calculation failed"));
                }
            }

            saveOptimizationResult();
            emit logMessageGenerated(tr("=== Optimization Results ==="));
            emit logMessageGenerated(
                tr("Final error: %1").arg(m_globalBestFitness, 0, 'e', 4));
            emit logMessageGenerated(
                tr("Total iterations: %1").arg(m_currentIteration + 1));
            emit logMessageGenerated(
                tr("Total evaluations: %1 (successful: %2)")
                .arg(m_totalEvaluations).arg(m_successfulEvaluations));

            QString finalParams = tr("Optimized parameters: ");
            for(int i = 0; i < m_globalBestPosition.size(); ++i) {
                finalParams += QString("[%1]=%2 ")
                               .arg(i).arg(m_globalBestPosition[i], 0, 'f', 6);
            }
            emit logMessageGenerated(finalParams);

            if(finalFullSolverExecuted && finalFullSolverSucceeded) {
                emit logMessageGenerated(
                    tr("Parameters and full-field results applied successfully to data manager"));
            } else if(!finalFullSolverExecuted) {
                emit logMessageGenerated(
                    tr("Optimized parameters applied to data manager"));
            }
        } catch(const std::exception& e) {
            finalFullSolverSucceeded = false;
            m_lastError = QString("Failed to apply final parameters: %1").arg(e.what());
            emit logMessageGenerated(
                tr("ERROR: Failed to apply final parameters: %1").arg(e.what()));
        } catch(...) {
            finalFullSolverSucceeded = false;
            m_lastError = "Failed to apply final parameters due to unknown error";
            emit logMessageGenerated(
                tr("ERROR: Unknown error applying final parameters"));
        }
    }

    m_isRunning = false;
    bool success = false;
    QString message;

    if(finalReason == LM_TARGET_ACHIEVED) {
        success = true;
        message = QString(tr("Target achieved. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING SUCCESSFUL ==="));
    } else if(finalReason == LM_TRUE_CONVERGENCE) {
        success = true;
        message = QString(
            tr("Automatic fitting converged to a stable solution. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING CONVERGED ==="));
    } else if(finalReason == LM_LOCAL_OPTIMUM) {
        success = true;
        message = QString(
            tr("Automatic fitting reached a local optimum. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING - LOCAL OPTIMUM ==="));
    } else if(finalReason == LM_MAX_ITERATIONS) {
        success = true;
        message = QString(tr("Max iterations reached. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING - MAX ITERATIONS ==="));
    } else if(finalReason == LM_USER_STOPPED) {
        success = true;
        message = QString(tr("Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING STOPPED BY USER ==="));
    } else if(finalReason == LM_CONSECUTIVE_FAILURES) {
        message = QString(
            tr("Automatic fitting failed due to consecutive failures. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING FAILED ==="));
    } else {
        message = QString(
            tr("Automatic fitting ended unexpectedly. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING - UNKNOWN END ==="));
    }

    if(!finalFullSolverSucceeded) {
        success = false;
        message = m_lastError;
    }

    emitRunSummary(success, finalReason);
    emit progressUpdated(m_maxIterations, m_globalBestFitness);
    QApplication::processEvents();
    msleep(200);
    QApplication::processEvents();
    closeTraceFile();
    cleanupTemporaryDirectory();
    emit fittingFinished(success, message);
    return success;
}

void nmCalculationAutoFitLM::extractUserInitialValues()
{
    // 从当前项目模型读取用户已有初始参数。
    // 只提取用户勾选的参数，并按 m_enabledParamIndices 的顺序写入 m_initialValues。
    // 这些值用于初始解真实评价、LM 起点和最终精英保护。
    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
    nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
    //QVector<nmDataWellBase*> wells = dataManager->getWellDataList();

    nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);

    m_initialValues.clear();

    // 按照启用参数的顺序提取初始值。井参数来自目标井，储层参数来自 reservoirData。
    for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
        int paramIndex = m_enabledParamIndices[i];
        double initialValue = 0.0;

        switch(paramIndex) {
            case 0: // 渗透率
                initialValue = reservoirData.getPermeability().getValue().toDouble();
                break;

            case 1: // 表皮系数
                if(pTargetWell) {
                    initialValue = pTargetWell->getPerforation(0)->getSkin().getValue().toDouble();
                }

                break;

            case 2: // 井筒储集系数
                if(pTargetWell) {
                    initialValue = pTargetWell->getWellboreStorage().getValue().toDouble();
                }

                break;

            case 3: // 孔隙度
                initialValue = reservoirData.getPorosity().getValue().toDouble();
                break;

            case 4: // 储层厚度
                initialValue = reservoirData.getThickness().getValue().toDouble();
                break;

            case 5: // 综合压缩系数
                initialValue = reservoirData.getCt().getValue().toDouble();
                break;

            case 6: // 岩石压缩系数
                initialValue = reservoirData.getCf().getValue().toDouble();
                break;

            case 7: // 初始含水饱和度
                initialValue = reservoirData.getSwi().getValue().toDouble();
                break;

            case 8: // 裂缝导流能力
                if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
                    nmDataVerticalFracturedWell* fracturedWell =
                        dynamic_cast<nmDataVerticalFracturedWell*>(pTargetWell);
                    if(fracturedWell) {
                        initialValue = fracturedWell->getDfc().getValue().toDouble();
                    }
                } else if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
                    nmDataHorizontalFracturedWell* fracturedWell =
                        dynamic_cast<nmDataHorizontalFracturedWell*>(pTargetWell);
                    if(fracturedWell) {
                        initialValue = fracturedWell->getDfc().getValue().toDouble();
                    }
                }
                break;

            case 9: // 裂缝半长
                if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
                    nmDataVerticalFracturedWell* fracturedWell =
                        dynamic_cast<nmDataVerticalFracturedWell*>(pTargetWell);
                    if(fracturedWell) {
                        initialValue = fracturedWell->getFractureHalfLength().getValue().toDouble();
                    }
                } else if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
                    nmDataHorizontalFracturedWell* fracturedWell =
                        dynamic_cast<nmDataHorizontalFracturedWell*>(pTargetWell);
                    if(fracturedWell) {
                        initialValue = fracturedWell->getFractureHalfLength().getValue().toDouble();
                    }
                }
                break;
        }

        m_initialValues.append(initialValue);
    }

    DEBUG_OUT(QString("Extracted %1 user initial values").arg(m_initialValues.size()));

    for(int i = 0; i < m_initialValues.size(); ++i) {
        DEBUG_OUT(QString("  Initial[%1] = %2").arg(i).arg(m_initialValues[i], 0, 'e', 3));
    }

    // 验证初始值
    if(!validateInitialValues()) {
        DEBUG_OUT("Warning: Some initial values are outside parameter bounds");
    }
}

bool nmCalculationAutoFitLM::evaluateTrustRegionPoint(
        const QVector<double>& parameters,
        double* fitness,
        AutoFitObjectiveBreakdownLM* breakdown,
        QVector<QVector<double> >* curve,
        int* elapsedMs)
{
    if(!fitness || !breakdown || !curve || !elapsedMs || m_shouldStop) {
        return false;
    }

    // evaluateFitness() 会写入 DataManager 并调用真实求解器。这里统一统计
    // 真实评价次数和耗时，同时严格要求固定残差、诊断结构和结果曲线均有效。
    QTime timer;
    timer.start();
    *fitness = evaluateFitness(parameters);
    *elapsedMs = timer.elapsed();
    *breakdown = m_lastObjectiveBreakdown;
    *curve = m_lastEvaluatedLogLogData;
    ++m_totalEvaluations;

    bool valid = isFiniteNumber(*fitness) && *fitness < 1.0e9 &&
                 breakdown->valid &&
                 trustRegionResidualsValid(*breakdown) &&
                 !curve->isEmpty();
    if(valid) {
        ++m_successfulEvaluations;
    }

    return valid;
}

StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
{
    const int dimensions = getEnabledParameterCount();
    if(dimensions <= 0 || m_enabledParamIndices.size() != dimensions) {
        m_lastError = tr("No valid parameters are available for trust-region fitting");
        return LM_OPTIMIZATION_FAILED;
    }

    // 真实求解次数比“外层迭代次数”更能反映耗时。预算至少允许完成一次全参数
    // 灵敏度和两次候选评价，同时避免连续重建 Jacobian 导致运行时间失控。
    const int maximumEvaluations = qMax(
        m_totalEvaluations + dimensions + 2,
        qMax(20, m_maxIterations * 3));
    // 下列步长均位于归一化内部坐标：0.04 表示参数范围的 4%，信赖半径
    // 限制一次联合移动的二范数，相关性门槛用于排除响应近乎共线的参数。
    const double sensitivityStep = 0.04;
    const double minimumCoordinateStep = 1.0e-5;
    const double minimumTrustRadius = 2.0e-3;
    const double maximumTrustRadius = 0.30;
    const double columnCorrelationLimit = 0.995;
    const double diagnosisThreshold = 1.0e-5;
    // 误差下降至少达到绝对 1e-5 且相对当前有效基准 0.2% 才算有效改善。
    // 更小的下降仍保留为最佳解，但不能反复清除停滞状态、延长拟合时间。
    const double effectiveRelativeImprovement = 2.0e-3;
    const double effectiveAbsoluteImprovement = 1.0e-5;
    const int maximumIneffectiveSteps = 3;

    // damping 是 LM 阻尼；拒绝或预测失准时增大，真实下降与预测一致时减小。
    // 两组累计量控制 Jacobian 重建，避免长期使用已偏离当前工作点的局部模型。
    double trustRadius = 0.12;
    double damping = 1.0e-2;
    int consecutiveRejectedSteps = 0;
    int consecutiveSolverFailures = 0;
    int acceptedSinceRebuild = 0;
    int consecutiveIneffectiveSteps = 0;
    double movementSinceRebuild = 0.0;
    bool rebuildRequested = true;
    bool modelRebuiltAtMinimumRadius = false;
    bool stagnationConfirmationRequested = false;
    StopReasonLM stopReason = LM_MAX_ITERATIONS;

    // jacobian 的行对应固定 160 维残差，列对应用户勾选的参数。
    // 三个 gradient 单独描述诊断分量对参数的局部变化，只用于本轮选参。
    QVector<QVector<double> > jacobian;
    QVector<double> verticalGradient(dimensions, 0.0);
    QVector<double> horizontalGradient(dimensions, 0.0);
    QVector<double> shapeGradient(dimensions, 0.0);
    QVector<bool> jacobianColumnValid(dimensions, false);

    // 参数向量的顺序始终与 m_enabledParamIndices 一致，不能按完整参数索引
    // 直接访问；下面两个转换函数集中维护这层映射关系。
    auto coordinatesFromParameters = [&](const QVector<double>& parameters)
            -> QVector<double> {
        QVector<double> coordinates(dimensions, 0.0);
        for(int i = 0; i < dimensions; ++i) {
            int parameterIndex = m_enabledParamIndices[i];
            coordinates[i] = toTrustRegionCoordinate(
                parameters[i], parameterIndex,
                m_parameterLower[parameterIndex],
                m_parameterUpper[parameterIndex]);
        }
        return coordinates;
    };

    auto parametersFromCoordinates = [&](const QVector<double>& coordinates)
            -> QVector<double> {
        QVector<double> parameters(dimensions, 0.0);
        for(int i = 0; i < dimensions; ++i) {
            int parameterIndex = m_enabledParamIndices[i];
            parameters[i] = fromTrustRegionCoordinate(
                coordinates[i], parameterIndex,
                m_parameterLower[parameterIndex],
                m_parameterUpper[parameterIndex]);
        }
        return parameters;
    };

    auto restoreEvaluationState = [&](const TrustRegionEvaluation& evaluation) {
        // evaluateFitness() 会把试算参数写入 DataManager。无论候选是否接受，
        // 下一次计算前都恢复到唯一的已接受工作点，防止失败试算污染后续求解。
        applyParametersToDataManager(evaluation.parameters);
        m_lastObjectiveBreakdown = evaluation.breakdown;
        m_lastEvaluatedLogLogData = evaluation.curve;
    };

    // 只有真实总误差更小的工作点才能发布为全局最优；曲线和诊断快照必须
    // 与参数同步更新，防止界面显示或最终精英保护使用错配的数据。
    auto publishAcceptedPoint = [&](const TrustRegionEvaluation& evaluation) {
        m_globalBestPosition = evaluation.parameters;
        m_globalBestFitness = evaluation.fitness;
        m_globalBestObjectiveBreakdown = evaluation.breakdown;
        m_globalBestLogLogData = evaluation.curve;
        emit bestCurveUpdated(m_targetLogLogData,
                              m_globalBestLogLogData,
                              m_currentIteration + 1,
                              m_globalBestFitness);
    };

    auto processPauseAndStop = [&]() -> bool {
        QApplication::processEvents();
        return !m_shouldStop;
    };

    // current 始终代表唯一已接受工作点。优先复用启动阶段已经真实验证的
    // 用户初始解，避免在信赖域入口重复调用一次昂贵求解器。
    TrustRegionEvaluation current;
    if(m_hasValidUserSolution &&
            m_globalBestPosition.size() == dimensions &&
            trustRegionResidualsValid(m_globalBestObjectiveBreakdown) &&
            !m_globalBestLogLogData.isEmpty()) {
        current.parameters = m_globalBestPosition;
        current.coordinates = coordinatesFromParameters(current.parameters);
        current.breakdown = m_globalBestObjectiveBreakdown;
        current.curve = m_globalBestLogLogData;
        current.fitness = m_globalBestFitness;
        current.elapsedMs = 0;
        current.valid = true;
    } else {
        // 用户初始解无效时只做一次确定性的范围中点回退；所有正值参数在对数
        // 坐标取中点，避免线性中点过分偏向跨数量级范围的上界。
        current.coordinates.fill(0.5, dimensions);
        current.parameters = parametersFromCoordinates(current.coordinates);
        current.valid = evaluateTrustRegionPoint(
            current.parameters,
            &current.fitness,
            &current.breakdown,
            &current.curve,
            &current.elapsedMs);
        writeTraceRow(-1, -1,
                      "trust_region_midpoint",
                      current.parameters,
                      current.fitness,
                      current.valid,
                      current.elapsedMs,
                      current.valid ? "midpoint_valid" : "midpoint_invalid",
                      current.valid ? &current.breakdown : nullptr);
        if(!current.valid) {
            m_lastError = tr("The initial solution and parameter-range midpoint are both invalid");
            return m_shouldStop
                   ? LM_USER_STOPPED
                   : LM_OPTIMIZATION_FAILED;
        }
        publishAcceptedPoint(current);
    }

    restoreEvaluationState(current);
    emit logMessageGenerated(tr("=== Starting LM Main Loop ==="));
    emit logMessageGenerated(
        tr("LM starting point error: %1; evaluation budget: %2")
        .arg(current.fitness, 0, 'e', 4)
        .arg(maximumEvaluations));

    // 有效改善始终相对“上一次有效改善后的误差”累计判断，避免一连串微小
    // 下降每次都清零计数；累计达到门槛后才开始新的有效改善基准。
    double effectiveImprovementBaseline = current.fitness;
    auto registerEffectiveImprovement = [&](double fitness) -> bool {
        const double requiredImprovement = qMax(
            effectiveAbsoluteImprovement,
            qAbs(effectiveImprovementBaseline) *
            effectiveRelativeImprovement);
        const double improvement = effectiveImprovementBaseline - fitness;
        if(improvement < requiredImprovement) {
            return false;
        }

        effectiveImprovementBaseline = fitness;
        consecutiveIneffectiveSteps = 0;
        stagnationConfirmationRequested = false;
        return true;
    };

    // 连续三次没有有效改善时只请求一次灵敏度重建。重建完成后由主循环
    // 直接检查累计改善，仍达不到门槛就判定局部收敛，不再继续微小试探。
    auto recordIneffectiveStep = [&]() -> bool {
        ++consecutiveIneffectiveSteps;
        if(consecutiveIneffectiveSteps < maximumIneffectiveSteps) {
            return false;
        }

        if(stagnationConfirmationRequested) {
            return true;
        }

        consecutiveIneffectiveSteps = 0;
        stagnationConfirmationRequested = true;
        rebuildRequested = true;
        emit logMessageGenerated(
            tr("No effective improvement for %1 consecutive steps; "
               "rebuilding sensitivity model for confirmation")
            .arg(maximumIneffectiveSteps));
        return false;
    };

    emit logMessageGenerated(
        tr("Effective improvement threshold: max(%1, %2% of baseline error); "
           "%3 consecutive ineffective steps trigger convergence confirmation")
        .arg(effectiveAbsoluteImprovement, 0, 'e', 2)
        .arg(effectiveRelativeImprovement * 100.0, 0, 'f', 2)
        .arg(maximumIneffectiveSteps));

    if(current.fitness < m_targetError) {
        return LM_TARGET_ACHIEVED;
    }

    // 在同一个真实工作点逐参数做单边差分。首选可用空间更大的方向；只有该方向
    // 求解失败时才补算反方向，因此初次建模通常每个参数只增加一次真实求解。
    auto rebuildSensitivity = [&]() -> bool {
        const TrustRegionEvaluation base = current;
        const int residualCount = base.breakdown.residualVector.size();
        if(residualCount <= 0) {
            return false;
        }

        jacobian = QVector<QVector<double> >(
            residualCount, QVector<double>(dimensions, 0.0));
        verticalGradient.fill(0.0, dimensions);
        horizontalGradient.fill(0.0, dimensions);
        shapeGradient.fill(0.0, dimensions);
        jacobianColumnValid.fill(false, dimensions);

        TrustRegionEvaluation bestProbe;
        int bestProbeColumn = -1;
        double bestProbeDelta = 0.0;
        // 差分步长不超过参数范围的 4%，信赖域收缩后同步减小，但保留 0.5%
        // 下限，避免步长太小使求解器数值噪声淹没真实灵敏度。
        const double finiteDifferenceStep = qMin(
            sensitivityStep,
            qMax(5.0e-3, trustRadius * 0.5));

        for(int column = 0;
                column < dimensions &&
                m_totalEvaluations < maximumEvaluations &&
                processPauseAndStop();
                ++column) {
            // 单边差分优先选择离边界空间更大的方向；首方向求解无效时才反向
            // 补算，因此正常情况下每个参数只消耗一次真实求解。
            double positiveRoom = 1.0 - base.coordinates[column];
            double negativeRoom = base.coordinates[column];
            double preferredSign = positiveRoom >= negativeRoom ? 1.0 : -1.0;
            bool columnBuilt = false;

            for(int directionAttempt = 0;
                    directionAttempt < 2 &&
                    !columnBuilt &&
                    m_totalEvaluations < maximumEvaluations;
                    ++directionAttempt) {
                double direction = directionAttempt == 0
                                   ? preferredSign : -preferredSign;
                double availableRoom = direction > 0.0
                                       ? positiveRoom : negativeRoom;
                double deltaMagnitude = qMin(
                    finiteDifferenceStep, availableRoom);
                if(deltaMagnitude < minimumCoordinateStep) {
                    continue;
                }

                TrustRegionEvaluation probe;
                probe.coordinates = base.coordinates;
                probe.coordinates[column] += direction * deltaMagnitude;
                probe.parameters = parametersFromCoordinates(probe.coordinates);
                probe.valid = evaluateTrustRegionPoint(
                    probe.parameters,
                    &probe.fitness,
                    &probe.breakdown,
                    &probe.curve,
                    &probe.elapsedMs);

                QString decision = probe.valid
                                   ? "sensitivity_valid"
                                   : (directionAttempt == 0
                                      ? "sensitivity_retry_opposite"
                                      : "sensitivity_invalid");
                writeTraceRow(m_currentIteration,
                              column,
                              "trust_region_sensitivity",
                              probe.parameters,
                              probe.fitness,
                              probe.valid,
                              probe.elapsedMs,
                              decision,
                              probe.valid ? &probe.breakdown : nullptr);

                if(!probe.valid) {
                    restoreEvaluationState(base);
                    continue;
                }

                double delta = probe.coordinates[column] -
                               base.coordinates[column];
                if(qAbs(delta) < minimumCoordinateStep ||
                        probe.breakdown.residualVector.size() != residualCount) {
                    restoreEvaluationState(base);
                    continue;
                }

                // 第 column 列是固定残差向量相对内部参数坐标的有限差分：
                // J[:,column] = (r_probe-r_base)/delta。
                for(int row = 0; row < residualCount; ++row) {
                    jacobian[row][column] =
                        (probe.breakdown.residualVector[row] -
                         base.breakdown.residualVector[row]) / delta;
                }

                // 有符号诊断量只有在基点和试算点都可靠时才能计算方向梯度；
                // shapeLoss 无方向可靠性标志，始终记录其局部变化率。
                if(base.breakdown.verticalReliable &&
                        probe.breakdown.verticalReliable &&
                        !base.breakdown.registrationAmbiguous &&
                        !probe.breakdown.registrationAmbiguous) {
                    verticalGradient[column] =
                        (probe.breakdown.verticalCommonBias -
                         base.breakdown.verticalCommonBias) / delta;
                }
                if(base.breakdown.horizontalReliable &&
                        probe.breakdown.horizontalReliable &&
                        !base.breakdown.registrationAmbiguous &&
                        !probe.breakdown.registrationAmbiguous) {
                    horizontalGradient[column] =
                        (probe.breakdown.horizontalPhysicalShift -
                         base.breakdown.horizontalPhysicalShift) / delta;
                }
                shapeGradient[column] =
                    (probe.breakdown.shapeLoss -
                     base.breakdown.shapeLoss) / delta;
                jacobianColumnValid[column] = true;
                columnBuilt = true;

                if(probe.fitness < base.fitness &&
                        (!bestProbe.valid ||
                         probe.fitness < bestProbe.fitness)) {
                    bestProbe = probe;
                    bestProbeColumn = column;
                    bestProbeDelta = delta;
                }
                restoreEvaluationState(base);
            }
        }

        int validColumnCount = 0;
        for(int i = 0; i < jacobianColumnValid.size(); ++i) {
            if(jacobianColumnValid[i]) {
                ++validColumnCount;
            }
        }
        if(validColumnCount == 0 || m_shouldStop) {
            restoreEvaluationState(base);
            return false;
        }

        // 灵敏度试算本身若找到更优真实解也应保留。所有列先基于同一个 base
        // 建完，再用该已知割线把 Jacobian 平移到新工作点，避免边算边移动基点。
        if(bestProbe.valid && bestProbeColumn >= 0) {
            QVector<double> acceptedStep(dimensions, 0.0);
            acceptedStep[bestProbeColumn] = bestProbeDelta;
            updateTrustRegionJacobian(
                &jacobian,
                base.breakdown.residualVector,
                bestProbe.breakdown.residualVector,
                acceptedStep);
            if(base.breakdown.verticalReliable &&
                    bestProbe.breakdown.verticalReliable) {
                updateTrustRegionScalarGradient(
                    &verticalGradient,
                    base.breakdown.verticalCommonBias,
                    bestProbe.breakdown.verticalCommonBias,
                    acceptedStep);
            }
            if(base.breakdown.horizontalReliable &&
                    bestProbe.breakdown.horizontalReliable) {
                updateTrustRegionScalarGradient(
                    &horizontalGradient,
                    base.breakdown.horizontalPhysicalShift,
                    bestProbe.breakdown.horizontalPhysicalShift,
                    acceptedStep);
            }
            updateTrustRegionScalarGradient(
                &shapeGradient,
                base.breakdown.shapeLoss,
                bestProbe.breakdown.shapeLoss,
                acceptedStep);

            current = bestProbe;
            publishAcceptedPoint(current);
            restoreEvaluationState(current);
            writeTraceRow(m_currentIteration,
                          bestProbeColumn,
                          "trust_region_sensitivity_accept",
                          current.parameters,
                          current.fitness,
                          true,
                          0,
                          "accepted_cached_probe",
                          &current.breakdown);
            emit logMessageGenerated(
                tr("Sensitivity probe accepted: error reduced to %1")
                .arg(current.fitness, 0, 'e', 4));
        } else {
            restoreEvaluationState(current);
        }

        acceptedSinceRebuild = 0;
        movementSinceRebuild = 0.0;
        consecutiveRejectedSteps = 0;
        rebuildRequested = false;
        // 若重建过程中接受了试算点，当前模型已通过割线平移而不是在新点完整
        // 重算；再遇到最小半径停滞时仍允许做一次真正的新点重建。
        modelRebuiltAtMinimumRadius =
            trustRadius <= minimumTrustRadius * 1.01 &&
            !bestProbe.valid;
        emit logMessageGenerated(
            tr("Sensitivity model rebuilt: %1/%2 parameter columns valid")
            .arg(validColumnCount)
            .arg(dimensions));
        return true;
    };

    int completedIterations = 0;
    for(int iteration = 0;
            iteration < m_maxIterations &&
            m_totalEvaluations < maximumEvaluations &&
            !m_shouldStop;
            ++iteration) {
        m_currentIteration = iteration;
        completedIterations = iteration + 1;

        if(!processPauseAndStop()) {
            break;
        }
        if(rebuildRequested) {
            const bool confirmingStagnation =
                stagnationConfirmationRequested;
            if(!rebuildSensitivity()) {
                stopReason = m_shouldStop
                             ? LM_USER_STOPPED
                             : LM_LOCAL_OPTIMUM;
                break;
            }
            if(current.fitness < m_targetError) {
                stopReason = LM_TARGET_ACHIEVED;
                break;
            }
            if(m_totalEvaluations >= maximumEvaluations) {
                stopReason = LM_MAX_ITERATIONS;
                break;
            }
            const bool rebuildEffective =
                registerEffectiveImprovement(current.fitness);
            if(confirmingStagnation && !rebuildEffective) {
                emit logMessageGenerated(
                    tr("Sensitivity rebuild produced no effective improvement; "
                       "local convergence detected"));
                stopReason = LM_LOCAL_OPTIMUM;
                break;
            }
        }

        // 先确定当前最突出的可靠诊断误差，用其梯度回答“哪些参数最能改善
        // 当前问题”；实际 LM 方向仍由完整残差梯度和 Jacobian 共同计算。
        int dominantComponent = trustRegionDominantComponent(
            current.breakdown, diagnosisThreshold);
        const QVector<double>* componentGradient = nullptr;
        if(dominantComponent == TRUST_REGION_VERTICAL_COMPONENT) {
            componentGradient = &verticalGradient;
        } else if(dominantComponent == TRUST_REGION_HORIZONTAL_COMPONENT) {
            componentGradient = &horizontalGradient;
        } else if(dominantComponent == TRUST_REGION_SHAPE_COMPONENT) {
            componentGradient = &shapeGradient;
        }

        // 主目标采用 0.5*||r||^2，其对参数的梯度为 J^T*r。这里不再叠加
        // vertical/horizontal/shape，保证诊断分量不会改变真实接受目标。
        QVector<double> totalGradient(dimensions, 0.0);
        for(int column = 0; column < dimensions; ++column) {
            if(!jacobianColumnValid[column]) {
                continue;
            }
            for(int row = 0; row < jacobian.size(); ++row) {
                totalGradient[column] +=
                    jacobian[row][column] *
                    current.breakdown.residualVector[row];
            }
        }

        // 每轮最多联合调整三个灵敏参数。按当前诊断梯度绝对值由大到小选取，
        // 并剔除 Jacobian 响应过度共线的列，降低弱可辨识参数互相补偿的风险。
        QVector<int> selectedColumns;
        QVector<bool> alreadyConsidered(dimensions, false);
        for(int selection = 0; selection < qMin(3, dimensions); ++selection) {
            int bestColumn = -1;
            double bestScore = 0.0;
            for(int column = 0; column < dimensions; ++column) {
                if(alreadyConsidered[column] ||
                        !jacobianColumnValid[column]) {
                    continue;
                }

                double score = componentGradient
                               ? qAbs((*componentGradient)[column])
                               : qAbs(totalGradient[column]);
                if(!isFiniteNumber(score) || score <= bestScore) {
                    continue;
                }

                bool excessivelyCorrelated = false;
                for(int selectedIndex = 0;
                        selectedIndex < selectedColumns.size();
                        ++selectedIndex) {
                    if(trustRegionJacobianColumnCorrelation(
                                jacobian,
                                column,
                                selectedColumns[selectedIndex]) >
                            columnCorrelationLimit) {
                        excessivelyCorrelated = true;
                        break;
                    }
                }
                if(!excessivelyCorrelated) {
                    bestColumn = column;
                    bestScore = score;
                }
            }
            if(bestColumn < 0 || bestScore <= 1.0e-12) {
                break;
            }
            selectedColumns.append(bestColumn);
            alreadyConsidered[bestColumn] = true;
        }

        // 诊断梯度接近零时，说明该分量在当前局部无法可靠选参，退回完整残差
        // 梯度，但接受标准仍然只有真实 total，诊断值不会重复计入目标函数。
        if(selectedColumns.isEmpty() && componentGradient) {
            dominantComponent = TRUST_REGION_TOTAL_COMPONENT;
            componentGradient = nullptr;
            alreadyConsidered.fill(false, dimensions);
            for(int selection = 0;
                    selection < qMin(3, dimensions);
                    ++selection) {
                int bestColumn = -1;
                double bestScore = 0.0;
                for(int column = 0; column < dimensions; ++column) {
                    if(alreadyConsidered[column] ||
                            !jacobianColumnValid[column]) {
                        continue;
                    }
                    double score = qAbs(totalGradient[column]);
                    if(score <= bestScore) {
                        continue;
                    }
                    bool excessivelyCorrelated = false;
                    for(int selectedIndex = 0;
                            selectedIndex < selectedColumns.size();
                            ++selectedIndex) {
                        if(trustRegionJacobianColumnCorrelation(
                                    jacobian,
                                    column,
                                    selectedColumns[selectedIndex]) >
                                columnCorrelationLimit) {
                            excessivelyCorrelated = true;
                            break;
                        }
                    }
                    if(!excessivelyCorrelated) {
                        bestColumn = column;
                        bestScore = score;
                    }
                }
                if(bestColumn < 0 || bestScore <= 1.0e-12) {
                    break;
                }
                selectedColumns.append(bestColumn);
                alreadyConsidered[bestColumn] = true;
            }
        }

        // 当前局部没有可用方向时先缩小半径并重建灵敏度；只有已经在最小
        // 半径完整重建后仍无方向，才把它判定为局部最优。
        if(selectedColumns.isEmpty()) {
            if(trustRadius <= minimumTrustRadius * 1.01 &&
                    modelRebuiltAtMinimumRadius) {
                stopReason = LM_LOCAL_OPTIMUM;
                break;
            }
            trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
            damping = qMin(1.0e8, damping * 4.0);
            rebuildRequested = true;
            if(recordIneffectiveStep()) {
                stopReason = LM_LOCAL_OPTIMUM;
                break;
            }
            continue;
        }

        // 在选中参数子空间构造 LM 正规方程：
        // (J^T*J + damping*diag(J^T*J))*step = -J^T*r。
        // 对角缩放使不同参数列的灵敏度量级差异不会直接改变阻尼强弱。
        const int selectedCount = selectedColumns.size();
        QVector<QVector<double> > normalMatrix(
            selectedCount, QVector<double>(selectedCount, 0.0));
        QVector<double> rightHandSide(selectedCount, 0.0);
        for(int left = 0; left < selectedCount; ++left) {
            int leftColumn = selectedColumns[left];
            rightHandSide[left] = -totalGradient[leftColumn];
            for(int right = 0; right < selectedCount; ++right) {
                int rightColumn = selectedColumns[right];
                for(int row = 0; row < jacobian.size(); ++row) {
                    normalMatrix[left][right] +=
                        jacobian[row][leftColumn] *
                        jacobian[row][rightColumn];
                }
            }
            double diagonalScale = qMax(
                1.0e-10, normalMatrix[left][left]);
            normalMatrix[left][left] += damping * diagonalScale;
        }

        QVector<double> selectedStep;
        bool solved = solveTrustRegionLinearSystem(
            normalMatrix, rightHandSide, &selectedStep);
        QVector<double> coordinateStep(dimensions, 0.0);
        if(solved) {
            for(int i = 0; i < selectedCount; ++i) {
                coordinateStep[selectedColumns[i]] = selectedStep[i];
            }
        }

        double stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
        if(!solved || !isFiniteNumber(stepNorm) ||
                stepNorm < minimumCoordinateStep) {
            // 正规方程退化时使用投影最速下降方向，仍只移动本轮已选择的参数。
            coordinateStep.fill(0.0, dimensions);
            double gradientNormSquared = 0.0;
            for(int i = 0; i < selectedCount; ++i) {
                int column = selectedColumns[i];
                double stepDirection = -totalGradient[column];
                if((current.coordinates[column] <= minimumCoordinateStep &&
                        stepDirection < 0.0) ||
                        (current.coordinates[column] >=
                         1.0 - minimumCoordinateStep &&
                         stepDirection > 0.0)) {
                    stepDirection = 0.0;
                }
                coordinateStep[column] = stepDirection;
                gradientNormSquared += stepDirection * stepDirection;
            }
            double gradientNorm = qSqrt(gradientNormSquared);
            if(gradientNorm > minimumCoordinateStep) {
                double scale = trustRadius / gradientNorm;
                for(int i = 0; i < selectedCount; ++i) {
                    int column = selectedColumns[i];
                    coordinateStep[column] *= scale;
                }
            }
            stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
        }

        // LM 解只给出局部模型建议方向；若超出当前信赖半径，保持方向不变并
        // 等比例截短，避免一次试算离开 Jacobian 有效的局部区域。
        if(stepNorm > trustRadius && stepNorm > 0.0) {
            double scale = trustRadius / stepNorm;
            for(int i = 0; i < coordinateStep.size(); ++i) {
                coordinateStep[i] *= scale;
            }
        }

        // 将 LM 步长投影到用户给定的参数范围，实际用于预测下降的也是投影后步长。
        QVector<double> candidateCoordinates = current.coordinates;
        for(int i = 0; i < dimensions; ++i) {
            candidateCoordinates[i] = qBound(
                0.0,
                current.coordinates[i] + coordinateStep[i],
                1.0);
            coordinateStep[i] = candidateCoordinates[i] -
                                current.coordinates[i];
        }
        stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));

        // 用线性模型 r_new ~= r_current + J*step 预测残差，再用平方能量
        // 的下降量与真实候选下降量比较，作为调整阻尼和半径的依据。
        QVector<double> predictedResidual =
            current.breakdown.residualVector;
        for(int row = 0; row < jacobian.size(); ++row) {
            for(int column = 0; column < dimensions; ++column) {
                predictedResidual[row] +=
                    jacobian[row][column] * coordinateStep[column];
            }
        }
        double predictedReduction = 0.5 *
            (trustRegionSquaredNorm(current.breakdown.residualVector) -
             trustRegionSquaredNorm(predictedResidual));

        // 无实际移动或模型预测不下降时没有必要调用昂贵求解器。将它按一次
        // 拒绝处理，并在连续发生后重建灵敏度，防止继续沿失效模型试算。
        if(stepNorm < minimumCoordinateStep ||
                !isFiniteNumber(predictedReduction) ||
                predictedReduction <= 1.0e-14) {
            trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
            damping = qMin(1.0e8, damping * 4.0);
            ++consecutiveRejectedSteps;
            if(consecutiveRejectedSteps >= 2) {
                if(trustRadius <= minimumTrustRadius * 1.01 &&
                        modelRebuiltAtMinimumRadius) {
                    stopReason = LM_LOCAL_OPTIMUM;
                    break;
                }
                rebuildRequested = true;
            }
            if(recordIneffectiveStep()) {
                stopReason = LM_LOCAL_OPTIMUM;
                break;
            }
            continue;
        }

        TrustRegionEvaluation candidate;
        candidate.coordinates = candidateCoordinates;
        candidate.parameters = parametersFromCoordinates(candidate.coordinates);
        candidate.valid = evaluateTrustRegionPoint(
            candidate.parameters,
            &candidate.fitness,
            &candidate.breakdown,
            &candidate.curve,
            &candidate.elapsedMs);

        if(!candidate.valid) {
            // 求解失败的候选不能改变 current。先完整恢复上一个已接受参数和
            // 对应误差快照，再缩小信赖域；连续失败达到上限才终止整个拟合。
            ++consecutiveSolverFailures;
            ++consecutiveRejectedSteps;
            trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
            damping = qMin(1.0e8, damping * 4.0);
            writeTraceRow(m_currentIteration,
                          -1,
                          "trust_region_candidate",
                          candidate.parameters,
                          candidate.fitness,
                          false,
                          candidate.elapsedMs,
                          "solver_invalid",
                          nullptr);
            restoreEvaluationState(current);
            if(consecutiveRejectedSteps >= 2) {
                rebuildRequested = true;
            }
            if(consecutiveSolverFailures >= m_maxConsecutiveFailures) {
                stopReason = LM_CONSECUTIVE_FAILURES;
                break;
            }
            if(recordIneffectiveStep()) {
                stopReason = LM_LOCAL_OPTIMUM;
                break;
            }
            continue;
        }

        consecutiveSolverFailures = 0;
        // 有效候选即使最终被拒绝，也提供了一条真实割线，可用于修正下一轮
        // 局部模型；是否成为新工作点仍只由下面的 total 严格比较决定。
        const AutoFitObjectiveBreakdownLM oldBreakdown = current.breakdown;
        updateTrustRegionJacobian(
            &jacobian,
            oldBreakdown.residualVector,
            candidate.breakdown.residualVector,
            coordinateStep);
        if(oldBreakdown.verticalReliable &&
                candidate.breakdown.verticalReliable &&
                !oldBreakdown.registrationAmbiguous &&
                !candidate.breakdown.registrationAmbiguous) {
            updateTrustRegionScalarGradient(
                &verticalGradient,
                oldBreakdown.verticalCommonBias,
                candidate.breakdown.verticalCommonBias,
                coordinateStep);
        }
        if(oldBreakdown.horizontalReliable &&
                candidate.breakdown.horizontalReliable &&
                !oldBreakdown.registrationAmbiguous &&
                !candidate.breakdown.registrationAmbiguous) {
            updateTrustRegionScalarGradient(
                &horizontalGradient,
                oldBreakdown.horizontalPhysicalShift,
                candidate.breakdown.horizontalPhysicalShift,
                coordinateStep);
        }
        updateTrustRegionScalarGradient(
            &shapeGradient,
            oldBreakdown.shapeLoss,
            candidate.breakdown.shapeLoss,
            coordinateStep);

        // reductionRatio 衡量局部线性模型的可信度：接近 1 表示预测准确；
        // 值较小表示虽然可能下降，但模型低估了非线性，需要收紧下一步。
        double actualReduction = 0.5 *
            (current.fitness * current.fitness -
             candidate.fitness * candidate.fitness);
        double reductionRatio = actualReduction / predictedReduction;
        bool accepted = candidate.fitness < current.fitness;
        QString componentName = trustRegionComponentName(dominantComponent);

        if(accepted) {
            // 真实总误差下降后才正式替换 current，并同步发布参数、曲线和诊断。
            // 模型预测可靠时减小阻尼并可扩大半径，预测较差时保守收缩。
            current = candidate;
            publishAcceptedPoint(current);
            restoreEvaluationState(current);
            ++acceptedSinceRebuild;
            movementSinceRebuild += stepNorm;
            consecutiveRejectedSteps = 0;

            if(reductionRatio > 0.75) {
                damping = qMax(1.0e-8, damping * 0.5);
                if(stepNorm >= trustRadius * 0.8) {
                    trustRadius = qMin(
                        maximumTrustRadius, trustRadius * 1.6);
                }
            } else if(reductionRatio > 0.25) {
                damping = qMax(1.0e-8, damping * 0.8);
            } else {
                damping = qMin(1.0e8, damping * 2.0);
                trustRadius = qMax(
                    minimumTrustRadius, trustRadius * 0.75);
            }

            if(acceptedSinceRebuild >= 6 ||
                    movementSinceRebuild >= 0.30) {
                rebuildRequested = true;
            }
            modelRebuiltAtMinimumRadius = false;
        } else {
            // 拒绝时 candidate 只保留在 trace 中，DataManager 和内存状态都恢复
            // 到 current。连续拒绝说明割线模型可能失真，因此请求重新试算灵敏度。
            ++consecutiveRejectedSteps;
            damping = qMin(1.0e8, damping * 4.0);
            trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
            restoreEvaluationState(current);
            if(consecutiveRejectedSteps >= 2) {
                rebuildRequested = true;
            }
        }

        // 候选只要更优就继续作为 current 保存；是否足以解除停滞，则统一
        // 相对上一次有效改善基准判断。拒绝和微小改善都会累计无效次数。
        const bool effectiveImprovement =
            registerEffectiveImprovement(current.fitness);
        if(!effectiveImprovement && recordIneffectiveStep()) {
            stopReason = LM_LOCAL_OPTIMUM;
        }

        writeTraceRow(m_currentIteration,
                      -1,
                      "trust_region_candidate",
                      candidate.parameters,
                      candidate.fitness,
                      true,
                      candidate.elapsedMs,
                      accepted
                      ? QString("accepted_%1").arg(componentName)
                      : QString("rejected_%1").arg(componentName),
                      &candidate.breakdown);

        QString componentDisplayName = componentName;
        if(componentName == "vertical") {
            componentDisplayName = tr("vertical deviation");
        } else if(componentName == "horizontal") {
            componentDisplayName = tr("horizontal deviation");
        } else if(componentName == "shape") {
            componentDisplayName = tr("shape deviation");
        } else if(componentName == "total") {
            componentDisplayName = tr("total error");
        }

        emit logMessageGenerated(
            tr("Iteration %1: focus=%2, parameters=%3, error=%4, result=%5")
            .arg(iteration + 1)
            .arg(componentDisplayName)
            .arg(selectedColumns.size())
            .arg(candidate.fitness, 0, 'e', 4)
            .arg(accepted ? tr("accepted") : tr("rejected")));
        emit progressUpdated(iteration + 1, m_globalBestFitness);

        if(stopReason == LM_LOCAL_OPTIMUM) {
            break;
        }

        if(current.fitness < m_targetError) {
            stopReason = LM_TARGET_ACHIEVED;
            break;
        }
        if(trustRadius <= minimumTrustRadius * 1.01 &&
                consecutiveRejectedSteps >= 2) {
            if(modelRebuiltAtMinimumRadius) {
                stopReason = LM_LOCAL_OPTIMUM;
                break;
            }
            rebuildRequested = true;
        }
    }

    if(completedIterations > 0) {
        m_currentIteration = completedIterations - 1;
    }
    restoreEvaluationState(current);

    if(m_shouldStop) {
        return LM_USER_STOPPED;
    }
    if(current.fitness < m_targetError) {
        return LM_TARGET_ACHIEVED;
    }
    if(stopReason == LM_CONSECUTIVE_FAILURES ||
            stopReason == LM_LOCAL_OPTIMUM ||
            stopReason == LM_OPTIMIZATION_FAILED) {
        return stopReason;
    }
    return LM_MAX_ITERATIONS;
}

double nmCalculationAutoFitLM::evaluateFitness(const QVector<double>& parameters)
{
    // LM 候选评价函数，也是自动拟合最核心的闭环：
    // 1. 校验候选参数是否在用户设置的上下界和基本物理范围内；
    // 2. 将参数写入 DataManager 的储层/目标井对象；
    // 3. 调用真实数值求解器，生成模拟结果；
    // 4. 从本次求解任务读取目标井 result log-log 曲线；
    // 5. 与目标 history log-log 曲线计算误差，误差越小代表拟合越好。
    //
    // 返回 1e10 表示候选评价失败或结果不可用。
    const QString funcName = QString("evaluateError[%1]").arg(m_currentIteration);
    static int callCount = 0;
    callCount++;
    m_lastEvaluatedLogLogData.clear();
    m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();

    try {
        DEBUG_OUT(QString("%1: Call #%2 - Starting evaluation with %3 parameters")
                  .arg(funcName).arg(callCount).arg(parameters.size()));

        // 打印参数值
        QString paramStr = "Parameters: ";

        for(int i = 0; i < parameters.size(); ++i) {
            paramStr += QString("[%1]=%2 ").arg(i).arg(parameters[i], 0, 'f', 6);
        }

        DEBUG_OUT(QString("%1: %2").arg(funcName).arg(paramStr));

        // 1. 参数有效性检查。这里先拦截明显非法的候选，
        // 避免把非有限数、越界值或极端危险值传给求解器。
        if(!validateParameters(parameters)) {
            DEBUG_OUT(QString("%1: Call #%2 - VALIDATION FAILED").arg(funcName).arg(callCount));

            // 详细检查每个参数
            for(int i = 0; i < parameters.size(); ++i) {
                if(!isFiniteNumber(parameters[i])) {
                    DEBUG_OUT(QString("  -> Param[%1] is NOT finite: %2").arg(i).arg(parameters[i]));
                }

                if(i < m_enabledParamIndices.size()) {
                    int paramIndex = m_enabledParamIndices[i];

                    if(paramIndex >= 0 && paramIndex < m_parameterLower.size()) {
                        double lower = m_parameterLower[paramIndex];
                        double upper = m_parameterUpper[paramIndex];

                        if(parameters[i] < lower) {
                            DEBUG_OUT(QString("  -> Param[%1]=%2 < lower bound %3")
                                      .arg(i).arg(parameters[i]).arg(lower));
                        }

                        if(parameters[i] > upper) {
                            DEBUG_OUT(QString("  -> Param[%1]=%2 > upper bound %3")
                                      .arg(i).arg(parameters[i]).arg(upper));
                        }

                        // 检查危险值。这些条件不是严格物理模型定义，
                        // 而是工程保护：避免求解器在明显异常输入下崩溃或返回无意义曲线。
                        switch(paramIndex) {
                            case 0: // 渗透率
                                if(parameters[i] <= 1e-6) {
                                    DEBUG_OUT(QString("  -> REJECTED: Permeability too small: %1").arg(parameters[i]));
                                }

                                break;

                            case 2: // 井筒储集系数
                                if(parameters[i] <= 1e-8) {
                                    DEBUG_OUT(QString("  -> REJECTED: Wellbore storage too small: %1").arg(parameters[i]));
                                }

                                break;

                            case 3: // 孔隙度
                                if(parameters[i] <= 1e-4 || parameters[i] >= 0.95) {
                                    DEBUG_OUT(QString("  -> REJECTED: Unrealistic porosity: %1").arg(parameters[i]));
                                }

                                break;
                        }
                    }
                }
            }

            return 1e10;
        }

        DEBUG_OUT(QString("%1: Call #%2 - Parameters validated OK").arg(funcName).arg(callCount));

        // 2. 数据管理器检查。后续参数写回和求解器组装都依赖当前 DataManager。
        nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();

        if(!dataManager) {
            DEBUG_OUT(QString("%1: Call #%2 - DataManager is NULL").arg(funcName).arg(callCount));
            return 1e10;
        }

        DEBUG_OUT(QString("%1: Call #%2 - DataManager OK").arg(funcName).arg(callCount));

        // 3. 应用参数。parameters 的顺序与 m_enabledParamIndices 对齐，
        // applyParametersToDataManager() 会把它们拆分写入储层参数和目标井参数。
        try {
            DEBUG_OUT(QString("%1: Call #%2 - Applying parameters...").arg(funcName).arg(callCount));
            applyParametersToDataManager(parameters);
            DEBUG_OUT(QString("%1: Call #%2 - Parameters applied successfully").arg(funcName).arg(callCount));
        } catch(const std::exception& e) {
            DEBUG_OUT(QString("%1: Call #%2 - FAILED to apply parameters: %3")
                      .arg(funcName).arg(callCount).arg(e.what()));
            return 1e10;
        }

        // Dfc 和裂缝半长属于网格输入。标记网格失效，使下一次任务基于当前参数快照重建。
        const bool fractureGridParameterSelected =
            (m_parameterSelected.size() > 8 && m_parameterSelected[8]) ||
            (m_parameterSelected.size() > 9 && m_parameterSelected[9]);
        if(fractureGridParameterSelected) {
            dataManager->invalidatePebiGrid();
        }

        // 4. 运行求解器。真实求解器偶发失败时允许重试，避免一次 DLL 调用异常
        // 直接让整个粒子评价失败。
        QVector<QVector<double>> solverResult;
        const int maxRetries = 2;
        bool solverSuccess = false;

        for(int retry = 0; retry <= maxRetries; ++retry) {
            if(m_shouldStop) {
                DEBUG_OUT(QString("%1: Call #%2 - User stop requested").arg(funcName).arg(callCount));
                return 1e10;
            }

            try {
                DEBUG_OUT(QString("%1: Call #%2 - Solver attempt %3/%4")
                          .arg(funcName).arg(callCount).arg(retry + 1).arg(maxRetries + 1));

                if(retry > 0) {
                    DEBUG_OUT(QString("%1: Call #%2 - Retry delay...").arg(funcName).arg(callCount));
                    msleep(1000);
                }

                solverResult = runSolver();

                // 详细检查求解器结果
                if(solverResult.isEmpty()) {
                    DEBUG_OUT(QString("%1: Call #%2 - Solver returned EMPTY result").arg(funcName).arg(callCount));
                } else {
                    DEBUG_OUT(QString("%1: Call #%2 - Solver returned %3 arrays")
                              .arg(funcName).arg(callCount).arg(solverResult.size()));

                    for(int i = 0; i < solverResult.size(); ++i) {
                        DEBUG_OUT(QString("  -> Array[%1] size: %2").arg(i).arg(solverResult[i].size()));
                    }

                    if(validateSolverResult(solverResult)) {
                        DEBUG_OUT(QString("%1: Call #%2 - Solver result VALIDATED on attempt %3")
                                  .arg(funcName).arg(callCount).arg(retry + 1));
                        solverSuccess = true;
                        break;
                    } else {
                        DEBUG_OUT(QString("%1: Call #%2 - Solver result VALIDATION FAILED on attempt %3")
                                  .arg(funcName).arg(callCount).arg(retry + 1));
                    }
                }

            } catch(const std::exception& e) {
                DEBUG_OUT(QString("%1: Call #%2 - Solver EXCEPTION on attempt %3: %4")
                          .arg(funcName).arg(callCount).arg(retry + 1).arg(e.what()));
            }
        }

        if(!solverSuccess) {
            DEBUG_OUT(QString("%1: Call #%2 - ALL SOLVER ATTEMPTS FAILED").arg(funcName).arg(callCount));
            return 1e10;
        }

        // 5. 获取 LogLog 数据。runSolverDll() 直接从求解任务复制目标井曲线，
        // 不再依赖 DataManager 中可能被其它井或上一粒子改写的共享结果。
        QVector<QVector<double>> resultLogLogData = m_lastEvaluatedLogLogData;

        try {
            if(!validateLogLogData(resultLogLogData)) {
                DEBUG_OUT(QString("%1: Call #%2 - LogLog data VALIDATION FAILED")
                          .arg(funcName).arg(callCount));

                // 详细输出LogLog数据问题
                if(resultLogLogData.size() < 3) {
                    DEBUG_OUT(QString("  -> LogLog arrays count: %1 (need 3)")
                              .arg(resultLogLogData.size()));
                } else {
                    DEBUG_OUT(QString("  -> LogLog array sizes: X=%1, Y1=%2, Y2=%3")
                              .arg(resultLogLogData[0].size())
                              .arg(resultLogLogData[1].size())
                              .arg(resultLogLogData[2].size()));

                    // 检查数据有效性
                    for(int i = 0; i < qMin(5, resultLogLogData[0].size()); ++i) {
                        if(!isFiniteNumber(resultLogLogData[0][i]) ||
                                !isFiniteNumber(resultLogLogData[1][i]) ||
                                !isFiniteNumber(resultLogLogData[2][i])) {
                            DEBUG_OUT(QString("  -> Invalid data at index %1: X=%2, Y1=%3, Y2=%4")
                                      .arg(i)
                                      .arg(resultLogLogData[0][i])
                                      .arg(resultLogLogData[1][i])
                                      .arg(resultLogLogData[2][i]));
                        }
                    }
                }

                return 1e10;
            }

            DEBUG_OUT(QString("%1: Call #%2 - LogLog data validated, size: %3")
                      .arg(funcName).arg(callCount).arg(resultLogLogData[0].size()));

        } catch(const std::exception& e) {
            DEBUG_OUT(QString("%1: Call #%2 - Error getting LogLog data: %3")
                      .arg(funcName).arg(callCount).arg(e.what()));
            return 1e10;
        }

        // 6. 计算误差。这里比较的是目标井 history log-log 与当前模拟 result log-log。
        double error;

        try {
            DEBUG_OUT(QString("%1: Call #%2 - Calculating error...")
                      .arg(funcName).arg(callCount));

            error = calculateLogLogCurveError(m_targetLogLogData, resultLogLogData);

            if(!isFiniteNumber(error) || error < 0) {
                DEBUG_OUT(QString("%1: Call #%2 - INVALID error value: %3")
                          .arg(funcName).arg(callCount).arg(error));
                return 1e10;
            }

            DEBUG_OUT(QString("%1: Call #%2 - SUCCESS! Error = %3")
                      .arg(funcName).arg(callCount).arg(error, 0, 'e', 6));
            // 保存最后一次有效曲线，供 LM 候选评价和精英保护复用。
            m_lastEvaluatedLogLogData = resultLogLogData;

        } catch(const std::exception& e) {
            DEBUG_OUT(QString("%1: Call #%2 - Error calculation FAILED: %3")
                      .arg(funcName).arg(callCount).arg(e.what()));
            return 1e10;
        }

        return error;

    } catch(const std::exception& e) {
        DEBUG_OUT(QString("%1: Call #%2 - TOP-LEVEL EXCEPTION: %3")
                  .arg(funcName).arg(callCount).arg(e.what()));
        return 1e10;
    } catch(...) {
        DEBUG_OUT(QString("%1: Call #%2 - UNKNOWN TOP-LEVEL EXCEPTION")
                  .arg(funcName).arg(callCount));
        return 1e10;
    }
}

// ==================== 参数应用方法 ====================

void nmCalculationAutoFitLM::applyParametersToDataManager(const QVector<double>& parameters)
{
    // 将粒子的“启用参数向量”写回项目数据。
    // parameters 的维度必须等于用户勾选的参数数量，顺序由 m_enabledParamIndices 决定。
    // 这里不直接跑求解器，只负责把 DataManager 调整到该粒子对应的模型状态。
    if(parameters.size() != getEnabledParameterCount()) {
        return;
    }

    updateReservoirParameters(parameters);
    updateWellParameters(parameters);
}

void nmCalculationAutoFitLM::updateReservoirParameters(const QVector<double>& parameters)
{
    // 更新储层级参数。井级参数 skin/wellboreC 不在这里改，由 updateWellParameters() 负责。
    // 这里先取 DataManager 中 reservoir 的副本，修改后再整体写回 DataManager。
    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
    nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();

    // paramIndex 是粒子 position 中的索引；i 是完整 10 个参数体系中的索引。
    // 只有 m_parameterSelected[i] 为 true 时，才从 parameters 中消费一个值。
    int paramIndex = 0;

    for(int i = 0; i < m_parameterSelected.size(); ++i) {
        if(m_parameterSelected[i] && paramIndex < parameters.size()) {
            double value = parameters[paramIndex];

            switch(i) {
                case 0: // 渗透率
                    reservoirData.getPermeability().setValue(value);
                    break;

                case 3: // 孔隙度
                    reservoirData.getPorosity().setValue(value);
                    break;

                case 4: // 储层厚度
                    reservoirData.getThickness().setValue(value);
                    break;

                case 5: // 综合压缩系数
                    reservoirData.getCt().setValue(value);
                    break;

                case 6: // 岩石压缩系数
                    reservoirData.getCf().setValue(value);
                    break;

                case 7: // 初始含水饱和度
                    reservoirData.getSwi().setValue(value);
                    break;
            }

            paramIndex++;
        }
    }

    // 更新数据管理器
    dataManager->updateReservoirData(reservoirData);
}

void nmCalculationAutoFitLM::updateWellParameters(const QVector<double>& parameters)
{
    // 更新目标井上的拟合参数。目前井级可拟合参数主要是：
    // - skin：写入第一个 perforation；
    // - wellboreC：写入井筒储集系数；
    // - Dfc/裂缝半长：只写入垂直压裂井或多段压裂水平井。
    // 如果目标井不存在或没有射孔数据，这里只记录 debug，不抛异常。
    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();

    // 只更新当前目标井的井参数，避免多井项目中误改其他井。
    //QVector<nmDataWellBase*> wells = dataManager->getWellDataList();
    nmDataWellBase* pWell = dataManager->findWellByName(m_targetWellName);

    if(!pWell) return;

    //nmDataWellBase* pWell = wells[0]; // 使用第一口井

    int paramIndex = 0;

    for(int i = 0; i < m_parameterSelected.size(); ++i) {
        if(m_parameterSelected[i] && paramIndex < parameters.size()) {
            double value = parameters[paramIndex];

            switch(i) {
                case 1: { // 表皮系数
                        nmDataPerforation* perf = pWell->getPerforation(0);

                        if(perf) {
                            nmDataAttribute skinAttr = perf->getSkin();
                            skinAttr.setValue(value);
                            perf->setSkin(skinAttr);
                        }
                    }
                    break;

                case 2: { // 井筒储集系数
                        nmDataAttribute wellboreAttr = pWell->getWellboreStorage();
                        wellboreAttr.setValue(value);
                        pWell->setWellboreStorage(wellboreAttr);
                    }
                    break;

                case 8: { // 裂缝导流能力
                        if(pWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
                            nmDataVerticalFracturedWell* fracturedWell =
                                dynamic_cast<nmDataVerticalFracturedWell*>(pWell);
                            if(fracturedWell) {
                                nmDataAttribute dfc = fracturedWell->getDfc();
                                dfc.setValue(value);
                                fracturedWell->setDfc(dfc);
                            }
                        } else if(pWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
                            nmDataHorizontalFracturedWell* fracturedWell =
                                dynamic_cast<nmDataHorizontalFracturedWell*>(pWell);
                            if(fracturedWell) {
                                nmDataAttribute dfc = fracturedWell->getDfc();
                                dfc.setValue(value);
                                fracturedWell->setDfc(dfc);
                            }
                        }
                    }
                    break;

                case 9: { // 裂缝半长
                        // 直接修改井对象中的属性，复用已有信号重算裂缝端点。
                        if(pWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
                            nmDataVerticalFracturedWell* fracturedWell =
                                dynamic_cast<nmDataVerticalFracturedWell*>(pWell);
                            if(fracturedWell) {
                                fracturedWell->getFractureHalfLength().setValue(value);
                            }
                        } else if(pWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
                            nmDataHorizontalFracturedWell* fracturedWell =
                                dynamic_cast<nmDataHorizontalFracturedWell*>(pWell);
                            if(fracturedWell) {
                                fracturedWell->getFractureHalfLength().setValue(value);
                            }
                        }
                    }
                    break;
            }

            paramIndex++;
        }
    }

    // 根据井类型更新到数据管理器
    updateWellToDataManager(pWell);
}

void nmCalculationAutoFitLM::updateWellToDataManager(nmDataWellBase* pWell)
{
    if(!pWell) return;

    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
    NM_WELL_MODEL wellType = pWell->getWellType();

    // DataManager 内部按井型维护不同容器。修改基类指针后，需要根据实际井型
    // 调用对应 update 接口，才能让后续求解器组装读到最新 skin / wellboreC。
    switch(wellType) {
        case NM_WELL_MODEL::Vertical_Well: {
                nmDataVerticalWell* pVerticalWell = dynamic_cast<nmDataVerticalWell*>(pWell);

                if(pVerticalWell) {
                    QVector<nmDataVerticalWell> wells;
                    wells.append(*pVerticalWell);
                    dataManager->updateVerticalWells(wells);
                }

                break;
            }

        case NM_WELL_MODEL::Vertical_Fractured_Well: {
                nmDataVerticalFracturedWell* pVFracturedWell = dynamic_cast<nmDataVerticalFracturedWell*>(pWell);

                if(pVFracturedWell) {
                    QVector<nmDataVerticalFracturedWell> wells;
                    wells.append(*pVFracturedWell);
                    dataManager->updateVerticalFracturedWells(wells);
                }

                break;
            }

        case NM_WELL_MODEL::Horizontal_Fractured_Well: {
                nmDataHorizontalFracturedWell* pHFracturedWell = dynamic_cast<nmDataHorizontalFracturedWell*>(pWell);

                if(pHFracturedWell) {
                    QVector<nmDataHorizontalFracturedWell> wells;
                    wells.append(*pHFracturedWell);
                    dataManager->updateHorizontalFracturedWells(wells);
                }

                break;
            }

        default:
            break;
    }
}

// ==================== 求解器相关方法 ====================

QVector<QVector<double>> nmCalculationAutoFitLM::runSolver()
{
    // 真实求解器统一走 DLL 方式，返回值由 evaluateFitness() 继续校验。
    return runSolverDll();
}

// ==================== 数据处理方法 ====================
// ==================== 算法辅助方法 ====================

void nmCalculationAutoFitLM::saveOptimizationResult()
{
    // 当前函数只做日志记录。真正把最优参数写回项目数据的是
    // startAutoFitting() 结束阶段的 applyParametersToDataManager(m_globalBestPosition)。
    DEBUG_OUT(QString("Optimization result: error=%1, evaluations=%2/%3")
              .arg(m_globalBestFitness, 0, 'e', 4)
              .arg(m_successfulEvaluations)
              .arg(m_totalEvaluations));
}

void nmCalculationAutoFitLM::validateAndProtectFinalResult()
{
    // 最终精英保护只阻止无效结果或真正变差的结果。任何真实误差下降都应保留，
    // 不能再用固定百分比门槛把已经找到的更优解恢复成初始值。
    if(!m_hasValidUserSolution) {
        emit logMessageGenerated(tr("No initial solution for elite protection"));
        return;
    }

    emit logMessageGenerated(tr("=== Final Result Validation (Elite Protection) ==="));

    // 使用已有的评估结果
    double finalFitness = m_globalBestFitness;
    double initialFitness = m_userInitialFitness;

    emit logMessageGenerated(tr("Comparing results: Initial=%1, Final=%2")
                             .arg(initialFitness, 0, 'e', 4).arg(finalFitness, 0, 'e', 4));

    bool finalValid = isFiniteNumber(finalFitness) &&
                      finalFitness < 1.0e9 &&
                      m_globalBestPosition.size() ==
                      m_userInitialSolution.size() &&
                      !m_globalBestLogLogData.isEmpty() &&
                      m_globalBestObjectiveBreakdown.valid;
    if(finalValid) {
        double improvement = initialFitness - finalFitness;
        double relativeImprovement =
            improvement / qMax(1.0e-10, qAbs(initialFitness));
        emit logMessageGenerated(tr("Improvement: %1 (%2%)")
                                 .arg(improvement, 0, 'e', 4)
                                 .arg(relativeImprovement * 100, 0, 'f', 2));
    }

    if(!finalValid || finalFitness > initialFitness) {
        emit logMessageGenerated(
            tr("Elite protection triggered: final result is invalid or worse than initial"));
        emit logMessageGenerated(tr("Restoring initial solution as final result"));

        m_globalBestFitness = initialFitness;
        m_globalBestPosition = m_userInitialSolution;
        m_globalBestLogLogData = m_userInitialLogLogData;
        m_globalBestObjectiveBreakdown = m_userInitialObjectiveBreakdown;
        emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, m_globalBestFitness);

        emit logMessageGenerated(tr("Initial solution restored successfully"));
    } else {
        emit logMessageGenerated(
            tr("Final result validated - solution is not worse than initial"));
    }
}

// ==================== 工具方法 ====================

int nmCalculationAutoFitLM::getEnabledParameterCount() const
{
    // 返回粒子维度，即用户勾选参与拟合的参数数量。
    int count = 0;

    for(int i = 0; i < m_parameterSelected.size(); ++i) {
        if(m_parameterSelected[i]) count++;
    }

    return count;
}

// ==================== 验证和处理方法 ====================

bool nmCalculationAutoFitLM::validateParameters(const QVector<double>& parameters) const
{
    // 参数物理范围已由拟合窗口统一校验；候选评价只检查维度、有限数和
    // 用户设置的上下界，避免另一套硬编码阈值与实际搜索范围冲突。
    if(parameters.size() != getEnabledParameterCount()) {
        return false;
    }

    for(int i = 0; i < parameters.size(); ++i) {
        if(!isFiniteNumber(parameters[i])) {
            return false;
        }

        // 检查参数范围
        if(i < m_enabledParamIndices.size()) {
            int paramIndex = m_enabledParamIndices[i];

            if(paramIndex >= 0 && paramIndex < m_parameterLower.size()) {
                if(parameters[i] < m_parameterLower[paramIndex] ||
                        parameters[i] > m_parameterUpper[paramIndex]) {
                    return false;
                }
            }
        }
    }

    return true;
}

bool nmCalculationAutoFitLM::validateLogLogData(const QVector<QVector<double>>& logLogData) const
{
    // 校验双对数曲线结构。约定：
    // logLogData[0]=time，logLogData[1]=pressure，logLogData[2]=pressure derivative。
    // 三列必须长度一致，且至少有足够点数用于插值和误差计算。
    if(logLogData.size() < 3) {
        DEBUG_OUT("LogLog data has less than 3 arrays");
        return false;
    }

    // 检查数组大小一致性
    int size = logLogData[0].size();

    if(size == 0) {
        DEBUG_OUT("Empty LogLog data");
        return false;
    }

    if(logLogData[1].size() != size || logLogData[2].size() != size) {
        DEBUG_OUT(QString("LogLog data size mismatch: X=%1, Y1=%2, Y2=%3")
                  .arg(logLogData[0].size())
                  .arg(logLogData[1].size())
                  .arg(logLogData[2].size()));
        return false;
    }

    // 检查最小数据点数
    if(size < 5) {
        DEBUG_OUT(QString("Too few LogLog data points: %1").arg(size));
        return false;
    }

    // 数据有效性检查
    for(int i = 0; i < size; ++i) {
        if(!isFiniteNumber(logLogData[0][i]) ||
                !isFiniteNumber(logLogData[1][i]) ||
                !isFiniteNumber(logLogData[2][i])) {
            DEBUG_OUT(QString("Invalid LogLog data at index %1").arg(i));
            return false;
        }
    }

    return true;
}

bool nmCalculationAutoFitLM::validateInitialValues() const
{
    // 检查当前模型读取出的初始参数是否和用户勾选维度一致，并且在上下界内。
    // 如果初始值越界，算法仍可继续，但日志会提示，因为精英保护可能不可用或效果变差。
    if(m_initialValues.size() != m_enabledParamIndices.size()) {
        DEBUG_OUT("Initial values count mismatch with enabled parameters");
        return false;
    }

    bool allValid = true;

    for(int i = 0; i < m_initialValues.size(); ++i) {
        int paramIndex = m_enabledParamIndices[i];
        double value = m_initialValues[i];

        if(!isFiniteNumber(value)) {
            DEBUG_OUT(QString("Initial value[%1] is not finite: %2").arg(i).arg(value));
            allValid = false;
            continue;
        }

        if(paramIndex < m_parameterLower.size() && paramIndex < m_parameterUpper.size()) {
            double minVal = m_parameterLower[paramIndex];
            double maxVal = m_parameterUpper[paramIndex];

            if(value < minVal || value > maxVal) {
                DEBUG_OUT(QString("Initial value[%1] = %2 is outside bounds [%3, %4]")
                          .arg(i).arg(value).arg(minVal).arg(maxVal));
                allValid = false;
            }
        }
    }

    return allValid;
}

bool nmCalculationAutoFitLM::validateSolverResult(const QVector<QVector<double>>& result) const
{
    // 校验求解器压力结果。这里检查的是 pressure result，至少需要 time 和 pressure 两列。
    // result log-log 的结构会在 validateLogLogData() 中另行检查。
    if(result.size() < 2) {
        DEBUG_OUT("Solver result has less than 2 arrays");
        return false;
    }

    if(result[0].size() != result[1].size()) {
        DEBUG_OUT(QString("Size mismatch: X=%1, Y=%2").arg(result[0].size()).arg(result[1].size()));
        return false;
    }

    if(result[0].size() == 0) {
        DEBUG_OUT("Empty solver result");
        return false;
    }

    // 检查最小数据点数
    if(result[0].size() < 10) {
        DEBUG_OUT(QString("Too few data points: %1").arg(result[0].size()));
        return false;
    }

    // 数据有效性检查
    for(int i = 0; i < result[0].size(); ++i) {
        if(!isFiniteNumber(result[0][i]) || !isFiniteNumber(result[1][i])) {
            DEBUG_OUT(QString("Invalid data at index %1: X=%2, Y=%3")
                      .arg(i).arg(result[0][i]).arg(result[1][i]));
            return false;
        }
    }

    // 检查X值单调性
    bool isMonotonic = true;

    for(int i = 1; i < result[0].size(); ++i) {
        if(result[0][i] <= result[0][i - 1]) {
            isMonotonic = false;
            break;
        }
    }

    if(!isMonotonic) {
        DEBUG_OUT("X values are not monotonically increasing");
    }

    return true;
}

double nmCalculationAutoFitLM::calculateLogLogCurveError(
    const QVector<QVector<double> >& target,
    const QVector<QVector<double> >& result) const
{
    // 主目标在目标与模拟曲线的公共时间范围内比较压力和导数残差；上下、左右
    // 和形状只负责诊断误差来源和选择参数，避免同一残差在 total 中被重复计算。
    // 整个计算过程均位于 log(time)-log(value) 坐标。
    const double invalidLoss = 1.0e10;
    const double valueFloor = 1.0e-12;
    const int numPoints = 80;
    m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();

    if(!validateLogLogData(target) || !validateLogLogData(result)) {
        return invalidLoss;
    }

    try {
        // 无法比较的采样行先跳过；有限但非正的导数无法进入双对数空间，
        // 当前数据又没有逐点有效掩码，因此遇到这种导数时判本次评价无效。
        auto prepareCurve = [valueFloor](const QVector<QVector<double> >& data,
                                         int firstIndex,
                                         QVector<QPointF>* pressure,
                                         QVector<QPointF>* derivative) -> bool {
            if(!pressure || !derivative || data.size() < 3 ||
                    data[0].size() != data[1].size() ||
                    data[0].size() != data[2].size() ||
                    firstIndex < 0 || firstIndex >= data[0].size()) {
                return false;
            }

            for(int i = firstIndex; i < data[0].size(); ++i) {
                if(!isFiniteNumber(data[0][i]) ||
                        !isFiniteNumber(data[1][i]) ||
                        !isFiniteNumber(data[2][i]) ||
                        data[0][i] <= 0.0 ||
                        data[1][i] <= 0.0) {
                    continue;
                }
                if(data[2][i] <= 0.0) {
                    return false;
                }

                pressure->append(QPointF(data[0][i], data[1][i]));
                derivative->append(
                    QPointF(data[0][i], qMax(data[2][i], valueFloor)));
            }

            // 求解器输出可能不是严格升序，且同一时刻可能出现重复记录。
            // 插值前统一排序并让后出现的记录覆盖同时间旧值，保证横坐标严格递增。
            auto sortAndUnique = [](QVector<QPointF>* curve) {
                std::stable_sort(
                    curve->begin(), curve->end(),
                    [](const QPointF& left, const QPointF& right) {
                        return left.x() < right.x();
                    });

                QVector<QPointF> unique;
                unique.reserve(curve->size());
                for(int i = 0; i < curve->size(); ++i) {
                    if(unique.isEmpty() ||
                            curve->at(i).x() > unique.last().x()) {
                        unique.append(curve->at(i));
                    } else {
                        unique[unique.size() - 1] = curve->at(i);
                    }
                }
                *curve = unique;
            };

            sortAndUnique(pressure);
            sortAndUnique(derivative);
            return pressure->size() >= 3 && derivative->size() >= 3;
        };

        QVector<QPointF> targetPressure;
        QVector<QPointF> targetDerivative;
        QVector<QPointF> resultPressure;
        QVector<QPointF> resultDerivative;
        // 模拟结果已跳过 DLL 首点，因此误差计算同步忽略目标曲线首点。
        if(!prepareCurve(target, 1, &targetPressure, &targetDerivative) ||
                !prepareCurve(result, 0, &resultPressure, &resultDerivative)) {
            return invalidLoss;
        }

        // 在双对数坐标中插值。二分定位用于后面的多次水平配准试算。
        auto interpolateLogValue = [valueFloor](
                const QVector<QPointF>& curve,
                double x,
                double* value) -> bool {
            if(!value || curve.size() < 2 || x <= 0.0 ||
                    x < curve.first().x() || x > curve.last().x()) {
                return false;
            }

            int low = 0;
            int high = curve.size() - 1;
            while(low < high) {
                int middle = low + (high - low) / 2;
                if(curve[middle].x() < x) {
                    low = middle + 1;
                } else {
                    high = middle;
                }
            }

            int right = qBound(1, low, curve.size() - 1);
            int left = right - 1;
            double leftLogX = qLn(curve[left].x());
            double rightLogX = qLn(curve[right].x());
            double denominator = rightLogX - leftLogX;
            double leftLogY =
                qLn(qMax(qAbs(curve[left].y()), valueFloor));
            double rightLogY =
                qLn(qMax(qAbs(curve[right].y()), valueFloor));

            if(qAbs(denominator) <= 1.0e-12) {
                *value = leftLogY;
            } else {
                double ratio = (qLn(x) - leftLogX) / denominator;
                *value = leftLogY +
                         ratio * (rightLogY - leftLogY);
            }
            return isFiniteNumber(*value);
        };

        const double targetMinX = targetPressure.first().x();
        const double targetMaxX = targetPressure.last().x();
        const double resultMinX = resultPressure.first().x();
        const double resultMaxX = resultPressure.last().x();
        if(targetMinX <= 0.0 || targetMaxX <= targetMinX ||
                resultMinX <= 0.0 || resultMaxX <= resultMinX) {
            return invalidLoss;
        }

        // 与 PSO 保持一致：只在目标与模拟曲线的时间交集内比较，不再设置
        // 覆盖率门槛，也不对交集之外的首尾数据做外推。
        const double overlapMinX = qMax(targetMinX, resultMinX);
        const double overlapMaxX = qMin(targetMaxX, resultMaxX);
        if(overlapMinX >= overlapMaxX) {
            return invalidLoss;
        }

        QVector<double> commonX(numPoints);
        QVector<double> commonLogX(numPoints);
        QVector<double> targetLogPressure(numPoints);
        QVector<double> targetLogDerivative(numPoints);
        const double comparisonLogMinX = qLn(overlapMinX);
        const double comparisonLogMaxX = qLn(overlapMaxX);

        // 在公共时间范围内生成固定维度的 log-time 网格，保持 LM 残差向量为 160 维。
        for(int i = 0; i < numPoints; ++i) {
            double logX = comparisonLogMinX +
                          static_cast<double>(i) *
                          (comparisonLogMaxX - comparisonLogMinX) /
                          (numPoints - 1);
            commonLogX[i] = logX;
			// 首尾直接使用原始端点，避免 exp(log(t)) 的舍入误差越过严格插值边界。
			if(i == 0) {
				commonX[i] = overlapMinX;
			} else if(i == numPoints - 1) {
				commonX[i] = overlapMaxX;
			} else {
				commonX[i] = qExp(logX);
			}

            if(!interpolateLogValue(
                        targetPressure, commonX[i],
                        &targetLogPressure[i]) ||
                    !interpolateLogValue(
                        targetDerivative, commonX[i],
                        &targetLogDerivative[i])) {
                return invalidLoss;
            }
        }

        QVector<double> pressureResidual(
            numPoints, std::numeric_limits<double>::quiet_NaN());
        QVector<double> derivativeResidual(
            numPoints, std::numeric_limits<double>::quiet_NaN());

        // 残差定义为“模拟减目标”：正值表示模拟曲线偏高，负值表示偏低。
        for(int i = 0; i < numPoints; ++i) {
            double resultLogPressure = 0.0;
            double resultLogDerivative = 0.0;
            if(!interpolateLogValue(
                        resultPressure, commonX[i],
                        &resultLogPressure) ||
                    !interpolateLogValue(
                        resultDerivative, commonX[i],
                        &resultLogDerivative)) {
                // 已位于结果时间范围内却无法插值说明数据存在内部断点，不能补线。
                return invalidLoss;
            }

            pressureResidual[i] =
                resultLogPressure - targetLogPressure[i];
            derivativeResidual[i] =
                resultLogDerivative - targetLogDerivative[i];
        }

        AutoFitObjectiveBreakdownLM breakdown;

        // 在指定中心附近计算普通均方根误差。
        auto rmseAround = [](
                const QVector<double>& values,
                int begin,
                int end,
                double center) -> double {
            double sum = 0.0;
            int count = 0;
            int validBegin = qMax(0, begin);
            int validEnd =
                qMin(end, static_cast<int>(values.size()));

            for(int i = validBegin; i < validEnd; ++i) {
                if(!isFiniteNumber(values[i])) {
                    continue;
                }

                double difference = values[i] - center;
                sum += difference * difference;
                ++count;
            }

            return count > 0
                   ? qSqrt(sum / count)
                   : std::numeric_limits<double>::quiet_NaN();
        };

        auto rmse = [&rmseAround](
                const QVector<double>& values,
                int begin,
                int end) -> double {
            return rmseAround(values, begin, end, 0.0);
        };

        // 普通算术平均中心保留上下偏差的符号。
        auto meanCenterRange = [](
                const QVector<double>& values,
                int begin,
                int end) -> double {
            int validBegin = qMax(0, begin);
            int validEnd =
                qMin(end, static_cast<int>(values.size()));
            double center = 0.0;
            int count = 0;

            for(int i = validBegin; i < validEnd; ++i) {
                if(isFiniteNumber(values[i])) {
                    center += values[i];
                    ++count;
                }
            }
            if(count == 0) {
                return std::numeric_limits<double>::quiet_NaN();
            }
            return center / count;
        };

        // 压力和导数合并后只求一个公共中心，表示两条曲线共同的上下位移。
        // 分别去中心会把压力与导数之间真实的相对形状差异一并消除。
        auto commonMeanCenterRange = [&meanCenterRange](
                const QVector<double>& pressureValues,
                const QVector<double>& derivativeValues,
                int begin,
                int end) -> double {
            QVector<double> combined;
            int validBegin = qMax(0, begin);
            int validEnd = qMin(
                end,
                qMin(static_cast<int>(pressureValues.size()),
                     static_cast<int>(derivativeValues.size())));
            combined.reserve(2 * qMax(0, validEnd - validBegin));

            for(int i = validBegin; i < validEnd; ++i) {
                if(isFiniteNumber(pressureValues[i])) {
                    combined.append(pressureValues[i]);
                }
                if(isFiniteNumber(derivativeValues[i])) {
                    combined.append(derivativeValues[i]);
                }
            }
            return meanCenterRange(combined, 0, combined.size());
        };

        // 两个通道按能量等权合并，返回值与单通道 RMSE 保持同一量纲。
        auto jointRmseAround = [&rmseAround](
                const QVector<double>& pressureValues,
                const QVector<double>& derivativeValues,
                int begin,
                int end,
                double center) -> double {
            double pressureLoss = rmseAround(
                pressureValues, begin, end, center);
            double derivativeLoss = rmseAround(
                derivativeValues, begin, end, center);
            if(!isFiniteNumber(pressureLoss) ||
                    !isFiniteNumber(derivativeLoss)) {
                return std::numeric_limits<double>::quiet_NaN();
            }
            return qSqrt(0.5 *
                         (pressureLoss * pressureLoss +
                          derivativeLoss * derivativeLoss));
        };

        // 主目标始终使用未做上下或左右校正的完整曲线误差。
        breakdown.pressureLoss =
            rmse(pressureResidual, 0, numPoints);
        breakdown.derivativeLoss =
            rmse(derivativeResidual, 0, numPoints);

        // 压力和导数各占一半权重。缩放后 residualVector 的二范数就是
        // sqrt(0.5 * pressureLoss^2 + 0.5 * derivativeLoss^2)。
        const double residualScale = qSqrt(0.5 / numPoints);
        breakdown.residualVector.reserve(2 * numPoints);
        for(int i = 0; i < numPoints; ++i) {
            breakdown.residualVector.append(
                residualScale *
                pressureResidual[i]);
        }
        for(int i = 0; i < numPoints; ++i) {
            breakdown.residualVector.append(
                residualScale *
                derivativeResidual[i]);
        }

        const double logGridStep =
            (comparisonLogMaxX - comparisonLogMinX) /
            (numPoints - 1);
        const double resultLogMinX = qLn(resultMinX);
        const double resultLogMaxX = qLn(resultMaxX);
        // 左右配准只在所有候选位移都共同覆盖的固定区间比较，至少保留 80%
        // 目标点；每个 log-time 网格间隔再细分为 8 份，提高位移分辨率。
        const int minimumRegistrationPoints =
            (numPoints * 4) / 5;
        const int shiftSubdivisions = 8;
        int maximumShiftIntervals = 4;
        int registrationBegin = 0;
        int registrationEnd = numPoints;
        double maximumPhysicalShift =
            maximumShiftIntervals * logGridStep;

        // 所有位移候选使用同一组目标点。结果范围不足时逐步缩小最大位移，
        // 但用于配准的固定公共区间不得少于目标网格的 80%。
        auto updateRegistrationRange = [&](double maximumShift) {
            registrationBegin = 0;
            registrationEnd = numPoints;
            while(registrationBegin < registrationEnd &&
                    commonLogX[registrationBegin] - maximumShift <
                    resultLogMinX - 1.0e-12) {
                ++registrationBegin;
            }
            while(registrationEnd > registrationBegin &&
                    commonLogX[registrationEnd - 1] + maximumShift >
                    resultLogMaxX + 1.0e-12) {
                --registrationEnd;
            }
        };

        updateRegistrationRange(maximumPhysicalShift);
        while(maximumShiftIntervals > 0 &&
                registrationEnd - registrationBegin <
                minimumRegistrationPoints) {
            --maximumShiftIntervals;
            maximumPhysicalShift =
                maximumShiftIntervals * logGridStep;
            updateRegistrationRange(maximumPhysicalShift);
        }
        if(registrationEnd - registrationBegin <
                minimumRegistrationPoints) {
            return invalidLoss;
        }

        // physicalShift 为正表示模拟曲线偏右；对齐时在目标时刻右侧读取模拟值。
        auto buildShiftResidual = [&](
                double physicalShift,
                int compareBegin,
                int compareEnd,
                QVector<double>* shiftedPressure,
                QVector<double>* shiftedDerivative) -> bool {
            if(!shiftedPressure || !shiftedDerivative) {
                return false;
            }

            shiftedPressure->fill(
                std::numeric_limits<double>::quiet_NaN(),
                numPoints);
            shiftedDerivative->fill(
                std::numeric_limits<double>::quiet_NaN(),
                numPoints);

            int validBegin = qMax(0, compareBegin);
            int validEnd = qMin(numPoints, compareEnd);
            for(int i = validBegin; i < validEnd; ++i) {
                double shiftedLogX =
                    commonLogX[i] + physicalShift;
                if(shiftedLogX < resultLogMinX - 1.0e-12 ||
                        shiftedLogX >
                        resultLogMaxX + 1.0e-12) {
                    return false;
                }

                // 对数时间还原后再次限制到原始端点，避免 exp(log(t)) 的
                // 舍入误差越过严格插值边界。
                double shiftedX = qBound(
                    resultMinX,
                    qExp(qBound(resultLogMinX,
                                shiftedLogX,
                                resultLogMaxX)),
                    resultMaxX);
                double resultLogPressure = 0.0;
                double resultLogDerivative = 0.0;
                if(!interpolateLogValue(
                            resultPressure, shiftedX,
                            &resultLogPressure) ||
                        !interpolateLogValue(
                            resultDerivative, shiftedX,
                            &resultLogDerivative)) {
                    return false;
                }

                (*shiftedPressure)[i] =
                    resultLogPressure - targetLogPressure[i];
                (*shiftedDerivative)[i] =
                    resultLogDerivative - targetLogDerivative[i];
            }
            return true;
        };

        // 损失相同时优先选择绝对位移更小的候选，避免平坦 profile 在数值噪声
        // 下无故偏向搜索边界。
        auto isBetterProfileValue = [](
                double loss,
                double shift,
                double bestLoss,
                double bestShift) -> bool {
            const double tolerance = 1.0e-12;
            return loss < bestLoss - tolerance ||
                   (qAbs(loss - bestLoss) <= tolerance &&
                    qAbs(shift) < qAbs(bestShift));
        };

        const int halfShiftStepCount =
            maximumShiftIntervals * shiftSubdivisions;
        const double physicalShiftStep =
            logGridStep / shiftSubdivisions;
        double zeroShiftCenteredLoss =
            std::numeric_limits<double>::quiet_NaN();
        double bestCenteredLoss =
            std::numeric_limits<double>::infinity();
        double bestPhysicalShift = 0.0;
        int bestShiftStep = 0;
        double bestPressureLoss =
            std::numeric_limits<double>::infinity();
        double bestPressureShift = 0.0;
        double bestDerivativeLoss =
            std::numeric_limits<double>::infinity();
        double bestDerivativeShift = 0.0;
        QVector<double> profileLosses(
            2 * halfShiftStepCount + 1,
            std::numeric_limits<double>::quiet_NaN());
        QVector<double> profileCommonBiases(
            2 * halfShiftStepCount + 1,
            std::numeric_limits<double>::quiet_NaN());
        QVector<double> shiftedPressureResidual;
        QVector<double> shiftedDerivativeResidual;

        // 位移和公共上下偏移联合求解，避免“先扣上下还是先扣左右”的顺序依赖。
        for(int shiftStep = -halfShiftStepCount;
                shiftStep <= halfShiftStepCount;
                ++shiftStep) {
            double physicalShift =
                shiftStep * physicalShiftStep;
            if(!buildShiftResidual(
                        physicalShift,
                        registrationBegin,
                        registrationEnd,
                        &shiftedPressureResidual,
                        &shiftedDerivativeResidual)) {
                continue;
            }

            double commonBias = commonMeanCenterRange(
                shiftedPressureResidual,
                shiftedDerivativeResidual,
                registrationBegin,
                registrationEnd);
            double centeredLoss = jointRmseAround(
                shiftedPressureResidual,
                shiftedDerivativeResidual,
                registrationBegin,
                registrationEnd,
                commonBias);
            double pressureBias = meanCenterRange(
                shiftedPressureResidual,
                registrationBegin,
                registrationEnd);
            double derivativeBias = meanCenterRange(
                shiftedDerivativeResidual,
                registrationBegin,
                registrationEnd);
            double pressureLoss = rmseAround(
                shiftedPressureResidual,
                registrationBegin,
                registrationEnd,
                pressureBias);
            double derivativeLoss = rmseAround(
                shiftedDerivativeResidual,
                registrationBegin,
                registrationEnd,
                derivativeBias);

            if(!isFiniteNumber(centeredLoss) ||
                    !isFiniteNumber(pressureLoss) ||
                    !isFiniteNumber(derivativeLoss)) {
                continue;
            }
            int profileIndex = shiftStep + halfShiftStepCount;
            profileLosses[profileIndex] = centeredLoss;
            profileCommonBiases[profileIndex] = commonBias;
            if(shiftStep == 0) {
                zeroShiftCenteredLoss = centeredLoss;
            }
            if(isBetterProfileValue(
                        centeredLoss,
                        physicalShift,
                        bestCenteredLoss,
                        bestPhysicalShift)) {
                bestCenteredLoss = centeredLoss;
                bestPhysicalShift = physicalShift;
                bestShiftStep = shiftStep;
            }
            if(isBetterProfileValue(
                        pressureLoss,
                        physicalShift,
                        bestPressureLoss,
                        bestPressureShift)) {
                bestPressureLoss = pressureLoss;
                bestPressureShift = physicalShift;
            }
            if(isBetterProfileValue(
                        derivativeLoss,
                        physicalShift,
                        bestDerivativeLoss,
                        bestDerivativeShift)) {
                bestDerivativeLoss = derivativeLoss;
                bestDerivativeShift = physicalShift;
            }
        }

        if(!isFiniteNumber(zeroShiftCenteredLoss) ||
                !isFiniteNumber(bestCenteredLoss) ||
                !isFiniteNumber(bestPressureLoss) ||
                !isFiniteNumber(bestDerivativeLoss)) {
            return invalidLoss;
        }

        // horizontalGain 是“允许水平位移”相对“固定零位移”减少的均方能量。
        // 只有改善足够明显且最优点不是边界，才把位移解释为可靠左右偏差。
        double horizontalGain = nestedRmsContribution(
            zeroShiftCenteredLoss, bestCenteredLoss);
        double horizontalSignalThreshold =
            qMax(1.0e-5, zeroShiftCenteredLoss * 0.02);
        int bestProfileIndex = bestShiftStep + halfShiftStepCount;
        double nearbyProfileLoss =
            std::numeric_limits<double>::infinity();
        int leftProfileIndex =
            bestProfileIndex - shiftSubdivisions;
        int rightProfileIndex =
            bestProfileIndex + shiftSubdivisions;
        if(leftProfileIndex >= 0 &&
                leftProfileIndex < profileLosses.size() &&
                isFiniteNumber(profileLosses[leftProfileIndex])) {
            nearbyProfileLoss = qMin(
                nearbyProfileLoss,
                profileLosses[leftProfileIndex]);
        }
        if(rightProfileIndex >= 0 &&
                rightProfileIndex < profileLosses.size() &&
                isFiniteNumber(profileLosses[rightProfileIndex])) {
            nearbyProfileLoss = qMin(
                nearbyProfileLoss,
                profileLosses[rightProfileIndex]);
        }
        double profileContrast = isFiniteNumber(nearbyProfileLoss)
                                 ? nestedRmsContribution(
                                       nearbyProfileLoss,
                                       bestCenteredLoss)
                                 : 0.0;
        bool flatRegistrationProfile =
            profileContrast <= horizontalSignalThreshold;

        // 平台曲线的 profile 也可能很平，但公共 bias 在各个位移下保持稳定，
        // 此时仍能可靠判断上下。只有近优位移会明显改变 bias 才说明上下/左右不可辨识。
        double minimumNearOptimalBias =
            std::numeric_limits<double>::infinity();
        double maximumNearOptimalBias =
            -std::numeric_limits<double>::infinity();
        for(int i = 0; i < profileLosses.size(); ++i) {
            if(isFiniteNumber(profileLosses[i]) &&
                    isFiniteNumber(profileCommonBiases[i]) &&
                    profileLosses[i] <=
                    bestCenteredLoss + horizontalSignalThreshold) {
                minimumNearOptimalBias = qMin(
                    minimumNearOptimalBias,
                    profileCommonBiases[i]);
                maximumNearOptimalBias = qMax(
                    maximumNearOptimalBias,
                    profileCommonBiases[i]);
            }
        }
        double nearOptimalBiasSpread =
            isFiniteNumber(minimumNearOptimalBias) &&
            isFiniteNumber(maximumNearOptimalBias)
            ? maximumNearOptimalBias - minimumNearOptimalBias
            : std::numeric_limits<double>::infinity();
        bool commonBiasStable = nearOptimalBiasSpread <= 1.0e-2;
        bool horizontalAtBoundary =
            halfShiftStepCount > 0 &&
            qAbs(bestShiftStep) == halfShiftStepCount;
        // 压力和导数通道分别求出的最佳位移若方向相反或相差过大，说明一个
        // 单一水平平移无法解释两条曲线，此时标记配准歧义并禁用左右引导。
        bool pressureShiftDetected =
            qAbs(bestPressureShift) >=
            0.5 * physicalShiftStep;
        bool derivativeShiftDetected =
            qAbs(bestDerivativeShift) >=
            0.5 * physicalShiftStep;
        bool channelShiftConflict =
            pressureShiftDetected &&
            derivativeShiftDetected &&
            (bestPressureShift * bestDerivativeShift < 0.0 ||
             qAbs(bestPressureShift - bestDerivativeShift) >
             2.0 * logGridStep);

        breakdown.horizontalLoss = horizontalGain;
        breakdown.horizontalReliable =
            maximumShiftIntervals > 0 &&
            !horizontalAtBoundary &&
            !channelShiftConflict &&
            !flatRegistrationProfile &&
            horizontalGain > horizontalSignalThreshold &&
            qAbs(bestPhysicalShift) >=
            0.5 * physicalShiftStep;
        breakdown.registrationAmbiguous =
            channelShiftConflict ||
            (qAbs(bestPhysicalShift) >=
             0.5 * physicalShiftStep &&
             !breakdown.horizontalReliable) ||
            (flatRegistrationProfile && !commonBiasStable);
        breakdown.horizontalPhysicalShift =
            breakdown.horizontalReliable
            ? bestPhysicalShift
            : 0.0;

        // 可信水平位移确定后，在该位移实际覆盖的最大区间重新计算上下和形状。
        int diagnosticBegin = 0;
        int diagnosticEnd = numPoints;
        while(diagnosticBegin < diagnosticEnd &&
                commonLogX[diagnosticBegin] +
                breakdown.horizontalPhysicalShift <
                resultLogMinX - 1.0e-12) {
            ++diagnosticBegin;
        }
        while(diagnosticEnd > diagnosticBegin &&
                commonLogX[diagnosticEnd - 1] +
                breakdown.horizontalPhysicalShift >
                resultLogMaxX + 1.0e-12) {
            --diagnosticEnd;
        }
        if(diagnosticEnd - diagnosticBegin <
                minimumRegistrationPoints ||
                !buildShiftResidual(
                    breakdown.horizontalPhysicalShift,
                    diagnosticBegin,
                    diagnosticEnd,
                    &shiftedPressureResidual,
                    &shiftedDerivativeResidual)) {
            return invalidLoss;
        }

        double commonBias = commonMeanCenterRange(
            shiftedPressureResidual,
            shiftedDerivativeResidual,
            diagnosticBegin,
            diagnosticEnd);
        double rawAlignedLoss = jointRmseAround(
            shiftedPressureResidual,
            shiftedDerivativeResidual,
            diagnosticBegin,
            diagnosticEnd,
            0.0);
        double centeredAlignedLoss = jointRmseAround(
            shiftedPressureResidual,
            shiftedDerivativeResidual,
            diagnosticBegin,
            diagnosticEnd,
            commonBias);
        if(!isFiniteNumber(commonBias) ||
                !isFiniteNumber(rawAlignedLoss) ||
                !isFiniteNumber(centeredAlignedLoss)) {
            return invalidLoss;
        }

        // 原始对齐误差减去公共中心后的能量差定义为上下误差贡献。只有它相对
        // 当前对齐误差足够明显，且配准无歧义时，公共 bias 才可用于有符号选参。
        breakdown.verticalCommonBias = commonBias;
        breakdown.verticalLoss = nestedRmsContribution(
            rawAlignedLoss, centeredAlignedLoss);
        breakdown.verticalReliable =
            !breakdown.registrationAmbiguous &&
            breakdown.verticalLoss >
            qMax(1.0e-5, rawAlignedLoss * 0.02);

        QVector<double> shapePressure(
            numPoints, std::numeric_limits<double>::quiet_NaN());
        QVector<double> shapeDerivative(
            numPoints, std::numeric_limits<double>::quiet_NaN());
        for(int i = diagnosticBegin; i < diagnosticEnd; ++i) {
            if(isFiniteNumber(shiftedPressureResidual[i])) {
                shapePressure[i] =
                    shiftedPressureResidual[i] - commonBias;
            }
            if(isFiniteNumber(shiftedDerivativeResidual[i])) {
                shapeDerivative[i] =
                    shiftedDerivativeResidual[i] - commonBias;
            }
        }

        // shapeLoss 是去除可信左右位移和公共均值中心后的剩余误差。
        double shapePressureLoss = rmse(
            shapePressure, diagnosticBegin, diagnosticEnd);
        double shapeDerivativeLoss = rmse(
            shapeDerivative, diagnosticBegin, diagnosticEnd);
        if(!isFiniteNumber(shapePressureLoss) ||
                !isFiniteNumber(shapeDerivativeLoss)) {
            return invalidLoss;
        }
        breakdown.shapeLoss = qSqrt(
            0.5 *
            (shapePressureLoss * shapePressureLoss +
             shapeDerivativeLoss * shapeDerivativeLoss));

        // 现阶段不识别或单独调度晚期流动段；保留字段只为了维持现有 trace 列。
        breakdown.lateDerivativeSlopeBias = 0.0;
        breakdown.lateDerivativeTrendLoss = 0.0;
        breakdown.lateDerivativeTrendReliable = false;

        if(!isFiniteNumber(breakdown.pressureLoss) ||
                !isFiniteNumber(breakdown.derivativeLoss) ||
                !isFiniteNumber(breakdown.verticalCommonBias) ||
                !isFiniteNumber(breakdown.verticalLoss) ||
                !isFiniteNumber(breakdown.horizontalPhysicalShift) ||
                !isFiniteNumber(breakdown.horizontalLoss) ||
                !isFiniteNumber(breakdown.shapeLoss) ||
                breakdown.residualVector.size() != 2 * numPoints) {
            return invalidLoss;
        }

        // LM 总目标等于固定残差向量的二范数；压力和导数各占一半能量。
        // 上下、左右和形状分量不参与候选排序与接受。
        breakdown.total = qSqrt(
            0.5 * breakdown.pressureLoss * breakdown.pressureLoss +
            0.5 * breakdown.derivativeLoss * breakdown.derivativeLoss);
        breakdown.valid =
            isFiniteNumber(breakdown.total) &&
            breakdown.total >= 0.0;
        m_lastObjectiveBreakdown = breakdown;

        DEBUG_OUT(
            QString("LogLog objective: pressure=%1, derivative=%2, vertical=%3, horizontal=%4, shape=%5, ambiguous=%6, shift=%7, total=%8")
            .arg(breakdown.pressureLoss, 0, 'e', 4)
            .arg(breakdown.derivativeLoss, 0, 'e', 4)
            .arg(breakdown.verticalLoss, 0, 'e', 4)
            .arg(breakdown.horizontalLoss, 0, 'e', 4)
            .arg(breakdown.shapeLoss, 0, 'e', 4)
            .arg(breakdown.registrationAmbiguous)
            .arg(breakdown.horizontalPhysicalShift, 0, 'e', 4)
            .arg(breakdown.total, 0, 'e', 4));

        return breakdown.valid
               ? qMin(1.0e9, breakdown.total)
               : invalidLoss;
    } catch(const std::exception& e) {
        DEBUG_OUT(
            QString("Exception in LogLog error calculation: %1")
            .arg(e.what()));
        return invalidLoss;
    } catch(...) {
        DEBUG_OUT("Unknown exception in LogLog error calculation");
        return invalidLoss;
    }
}

QVector<QVector<double>> nmCalculationAutoFitLM::runSolverDll()
{
    // DLL 求解器路径。
    // 这个函数负责把当前 DataManager 中的项目状态交给底层数值求解器，
    // 数值求解仍包含全部计算井，以保留井间干扰；后处理只提取目标井曲线。
    //
    // 如果这里失败，通常需要优先检查：HX_NWTM.dll、license、网格/井数据是否完整、
    // 目标井是否存在，以及 DataManager 中刚写入的参数是否导致求解器异常。
    DEBUG_OUT("SOLVER DLL START");

    // 创建任务时绑定当前分析，并在当前线程冻结本次求解所需的全部输入。
    nmDataAnalyzeManager* dataManager =
            nmDataAnalyzeManager::getCurrentInstance();
    if(!dataManager || m_targetWellName.isEmpty()) {
        DEBUG_OUT("Data manager or target well is unavailable");
        return QVector<QVector<double> >();
    }

    if(m_evaluationInProgress > 0) {
        DEBUG_OUT("DLL Solver already running, skipping");
        return QVector<QVector<double>>();
    }

    ++m_evaluationInProgress;
    m_lastEvaluatedLogLogData.clear();
    QVector<QVector<double>> result;
    nmCalculationDllPebiSolverTask* dllTask = nullptr;

    try {
        DEBUG_OUT("Creating DLL solver task");
        dllTask = new nmCalculationDllPebiSolverTask(
                    m_tempDirectory,
                    dataManager,
                    m_targetWellName);

        if(m_shouldStop) {
            DEBUG_OUT("Should stop - cleaning up and returning empty result");
            delete dllTask;
            --m_evaluationInProgress;
            return result;
        }

        DEBUG_OUT("Starting DLL solver execution...");

        // 异步执行 DLL 任务。循环等待期间持续 processEvents，保证界面不会完全卡死。
        dllTask->start();

        // 等待完成。最大等待 1 小时，适配大模型慢算；用户停止时会 terminate。
        const int maxWait = 3600000;  // 1h超时
        const int checkInterval = 50;
        QTime waitTimer;
        waitTimer.start();

        while(waitTimer.elapsed() < maxWait) {
            // wait(timeout) 会在线程一完成时立即返回，避免原来固定 msleep(500)
            // 带来的每次 0～500ms 额外等待；50ms 间隔仍可及时处理停止请求和界面事件。
            if(dllTask->wait(checkInterval)) {
                DEBUG_OUT("DLL solver task completed");
                break;
            }

            QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, checkInterval);

            if(m_shouldStop) {
                DEBUG_OUT("DLL solver task terminated by user");
                dllTask->terminate();
                break;
            }
        }

        // 超时处理
        if(dllTask->isRunning()) {
            DEBUG_OUT("DLL solver task timeout, terminating...");
            dllTask->terminate();
            dllTask->wait(2000);

            delete dllTask;
            dllTask = nullptr;
            --m_evaluationInProgress;
            m_consecutiveFailures++;
            return result;
        }

        // 线程结束后检查真实执行结果，防止失败时复用上一粒子的旧曲线。
        dllTask->wait();
        if(!dllTask->wasSuccessful()) {
            DEBUG_OUT("DLL solver task reported failure");
            delete dllTask;
            dllTask = nullptr;
            --m_evaluationInProgress;
            m_consecutiveFailures++;
            return result;
        }

        // 任务结束后复制其局部结果，删除任务前不再持有任务内部引用。
        QVector<QVector<double>> pressureResult = dllTask->getAutoFitResultPressure();
        QVector<QVector<double>> logLogResult = dllTask->getAutoFitResultLogLog();

        DEBUG_OUT(QString("DLL result verification - Pressure arrays: %1, LogLog arrays: %2")
                  .arg(pressureResult.size()).arg(logLogResult.size()));

        if(pressureResult.size() >= 2) {
            DEBUG_OUT(QString("Pressure result - Time points: %1, Pressure points: %2")
                      .arg(pressureResult[0].size()).arg(pressureResult[1].size()));

            if(pressureResult[0].size() > 0) {
                DEBUG_OUT(QString("Sample pressure data - Time[0]: %1, Time[last]: %2, P[0]: %3, P[last]: %4")
                          .arg(pressureResult[0][0])
                          .arg(pressureResult[0][pressureResult[0].size() - 1])
                          .arg(pressureResult[1][0])
                          .arg(pressureResult[1][pressureResult[1].size() - 1]));
            }
        }

        // 数据有效性检查
        if(pressureResult.size() >= 2
                && pressureResult[0].size() > 0
                && pressureResult[1].size() > 0
                && validateLogLogData(logLogResult)) {
            result = pressureResult;
            m_lastEvaluatedLogLogData = logLogResult;
            DEBUG_OUT(QString("Got DLL solver result: %1 points").arg(result[0].size()));
            m_consecutiveFailures = 0;

            // 检查结果是否与之前不同。若连续粒子得到完全相同的压力曲线，
            // 可能说明参数没有正确写入 DataManager，或求解器缓存/状态没有刷新。
            static QVector<double> lastPressureResult;
            bool isDifferentFromLast = false;

            if(lastPressureResult.isEmpty() || lastPressureResult.size() != pressureResult[1].size()) {
                isDifferentFromLast = true;
            } else {
                for(int i = 0; i < qMin(5, pressureResult[1].size()); ++i) {
                    if(qAbs(lastPressureResult[i] - pressureResult[1][i]) > 1e-12) {
                        isDifferentFromLast = true;
                        break;
                    }
                }
            }

            if(isDifferentFromLast) {
                DEBUG_OUT("RESULT VERIFICATION: Got NEW result data from DLL");
                lastPressureResult = pressureResult[1];
            } else {
                DEBUG_OUT("!!! WARNING: Result data appears to be identical to previous run !!!");
            }

        } else {
            DEBUG_OUT("DLL solver result is empty or invalid");
            DEBUG_OUT(QString("Pressure result size: %1, Array sizes: %2, %3")
                      .arg(pressureResult.size())
                      .arg(pressureResult.size() > 0 ? pressureResult[0].size() : 0)
                      .arg(pressureResult.size() > 1 ? pressureResult[1].size() : 0));
            m_consecutiveFailures++;
        }

    } catch(const std::bad_alloc& e) {
        DEBUG_OUT(QString("Memory allocation failed in DLL solver: %1").arg(e.what()));
        m_consecutiveFailures++;
    } catch(const std::exception& e) {
        DEBUG_OUT(QString("Exception in DLL solver: %1").arg(e.what()));
        m_consecutiveFailures++;
    } catch(...) {
        DEBUG_OUT("Unknown exception in DLL solver");
        m_consecutiveFailures++;
    }

    // 清理DLL任务
    if(dllTask) {
        if(dllTask->isRunning()) {
            dllTask->terminate();
            dllTask->wait(3000);
        }

        DEBUG_OUT("Cleaning up DLL solver task...");
        delete dllTask;
        dllTask = nullptr;
    }

    QApplication::processEvents(QEventLoop::AllEvents, 100);
    --m_evaluationInProgress;

    DEBUG_OUT(QString("SOLVER DLL END - ResultPoints: %1")
              .arg(result.isEmpty() ? 0 : result[0].size()));

    return result;
}

bool nmCalculationAutoFitLM::runFinalFullSolver()
{
    // 不设置目标井名，任务按原完整模式保存全部井和网格结果。
    nmDataAnalyzeManager* dataManager =
            nmDataAnalyzeManager::getCurrentInstance();
    if(!dataManager) {
        DEBUG_OUT("Cannot start final full-field solver without a data manager");
        return false;
    }

    if(m_evaluationInProgress > 0) {
        DEBUG_OUT("Cannot start final full-field solver while another evaluation is running");
        return false;
    }

    ++m_evaluationInProgress;
    nmCalculationDllPebiSolverTask dllTask(m_tempDirectory, dataManager);
    dllTask.start();

    const int maxWait = 3600000;
    const int checkInterval = 50;
    QTime waitTimer;
    waitTimer.start();

    while(waitTimer.elapsed() < maxWait) {
        if(dllTask.wait(checkInterval)) {
            break;
        }

        // 最终完整计算可能持续较长时间，此处需处理停止按钮事件。
        QApplication::processEvents(QEventLoop::AllEvents, checkInterval);

        if(!dllTask.isRunning()) {
            break;
        }

        if(m_shouldStop) {
            DEBUG_OUT("Final full-field solver terminated by user");
            dllTask.terminate();
            dllTask.wait(2000);
            --m_evaluationInProgress;
            return false;
        }
    }

    if(dllTask.isRunning()) {
        DEBUG_OUT("Final full-field solver timeout, terminating task");
        dllTask.terminate();
        dllTask.wait(2000);
        --m_evaluationInProgress;
        return false;
    }

    dllTask.wait();
    // 后台只生成局部结果快照；确认求解和输入版本均有效后再一次性写回当前分析。
    const bool succeeded =
            dllTask.wasSuccessful() && dllTask.commitResult(dataManager);
    --m_evaluationInProgress;
    return succeeded;
}

QString nmCalculationAutoFitLM::getStopReasonDescription(StopReasonLM reason) const
{
    // 将停止枚举转成人类可读文本，用于日志和运行摘要。
    switch(reason) {
        case LM_TARGET_ACHIEVED:
            return tr("Target error achieved");

        case LM_TRUE_CONVERGENCE:
            return tr("Algorithm converged to stable solution");

        case LM_LOCAL_OPTIMUM:
            return tr("Local optimum detected");

        case LM_MAX_ITERATIONS:
            return tr("Maximum iterations reached");

        case LM_USER_STOPPED:
            return tr("Stopped by user request");

        case LM_CONSECUTIVE_FAILURES:
            return tr("Too many consecutive failures");

        case LM_OPTIMIZATION_FAILED:
            return tr("Optimization failed");

        default:
            return tr("Unknown reason");
    }
}