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

#include "nmCalculationPebiGrid.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

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 TrustRegionTimeWindow
{
    TRUST_REGION_WELLBORE_STORAGE_WINDOW = 0,
    TRUST_REGION_TRANSITION_WINDOW,
    TRUST_REGION_IARF_WINDOW,
    TRUST_REGION_LATE_WINDOW,
    TRUST_REGION_TIME_WINDOW_COUNT,
    TRUST_REGION_TOTAL_WINDOW
};

// 一次真实求解的完整快照。除了参数和总误差，还保存内部坐标、诊断分量和
// 双对数曲线，因此拒绝候选后可以完整恢复上一个已接受工作点。
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 ||
            breakdown.timeWindowBegin.size() !=
                TRUST_REGION_TIME_WINDOW_COUNT ||
            breakdown.timeWindowEnd.size() !=
                TRUST_REGION_TIME_WINDOW_COUNT ||
            breakdown.timeWindowLoss.size() !=
                TRUST_REGION_TIME_WINDOW_COUNT) {
        return false;
    }

    for(int i = 0; i < breakdown.residualVector.size(); ++i) {
        if(!isFiniteNumber(breakdown.residualVector[i])) {
            return false;
        }
    }
    for(int window = 0;
            window < TRUST_REGION_TIME_WINDOW_COUNT;
            ++window) {
        if(breakdown.timeWindowBegin[window] < 0 ||
                breakdown.timeWindowEnd[window] <=
                    breakdown.timeWindowBegin[window] ||
                breakdown.timeWindowEnd[window] > 80 ||
                !isFiniteNumber(breakdown.timeWindowLoss[window])) {
            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;
}

// trace 和运行日志使用稳定的英文标识，便于离线过程分析按窗口筛选。
static QString trustRegionWindowName(int window)
{
    if(window == TRUST_REGION_WELLBORE_STORAGE_WINDOW) {
        return "w0_storage";
    }
    if(window == TRUST_REGION_TRANSITION_WINDOW) {
        return "w1_transition";
    }
    if(window == TRUST_REGION_IARF_WINDOW) {
        return "w2_iarf";
    }
    if(window == TRUST_REGION_LATE_WINDOW) {
        return "w3_late";
    }
    return "total";
}

// 四个窗口在目标曲线的固定对数时间网格上识别。井储段依据早期单位斜率，
// IARF 依据压力导数平台；特征不清楚时按固定比例分段，保证每次评价维度一致。
static void buildTrustRegionTimeWindows(
        const QVector<double>& time,
        const QVector<double>& logPressure,
        const QVector<double>& logDerivative,
        QVector<int>* windowBegin,
        QVector<int>* windowEnd)
{
    if(!windowBegin || !windowEnd) {
        return;
    }

    const int pointCount = qMin(
        time.size(), qMin(logPressure.size(), logDerivative.size()));
    windowBegin->clear();
    windowEnd->clear();
    if(pointCount < 16) {
        return;
    }

    QVector<double> logTime(pointCount, 0.0);
    QVector<double> derivativeSlope(
        pointCount, std::numeric_limits<double>::quiet_NaN());
    for(int i = 0; i < pointCount; ++i) {
        logTime[i] = qLn(time[i]);
    }

    // 五点局部回归比相邻两点差分更不容易把目标曲线的小幅波动误判成流态边界。
    for(int center = 0; center < pointCount; ++center) {
        int begin = qMax(0, center - 2);
        int end = qMin(pointCount, center + 3);
        double meanX = 0.0;
        double meanY = 0.0;
        for(int i = begin; i < end; ++i) {
            meanX += logTime[i];
            meanY += logDerivative[i];
        }
        const int count = end - begin;
        meanX /= count;
        meanY /= count;
        double covariance = 0.0;
        double variance = 0.0;
        for(int i = begin; i < end; ++i) {
            double offsetX = logTime[i] - meanX;
            covariance += offsetX * (logDerivative[i] - meanY);
            variance += offsetX * offsetX;
        }
        if(variance > 1.0e-14) {
            derivativeSlope[center] = covariance / variance;
        }
    }

    const int minimumWindowPoints = qMax(5, pointCount / 16);
    int storageEnd = qBound(
        minimumWindowPoints,
        static_cast<int>(pointCount * 0.25 + 0.5),
        pointCount - 3 * minimumWindowPoints);
    bool storageIdentified = false;
    int lastStoragePoint = -1;
    int nonStorageCount = 0;
    const int maximumStorageSearch = qMin(
        pointCount - 3 * minimumWindowPoints,
        static_cast<int>(pointCount * 0.45));
    for(int i = 0; i < maximumStorageSearch; ++i) {
        bool storageLike = isFiniteNumber(derivativeSlope[i]) &&
                           derivativeSlope[i] >= 0.65 &&
                           derivativeSlope[i] <= 1.35 &&
                           qAbs(logPressure[i] - logDerivative[i]) <= 1.0;
        if(storageLike) {
            storageIdentified = true;
            lastStoragePoint = i;
            nonStorageCount = 0;
        } else if(storageIdentified) {
            ++nonStorageCount;
            if(nonStorageCount >= 3) {
                break;
            }
        } else if(i >= minimumWindowPoints) {
            break;
        }
    }
    if(storageIdentified && lastStoragePoint + 1 >= minimumWindowPoints) {
        storageEnd = qBound(
            minimumWindowPoints,
            lastStoragePoint + 1,
            pointCount - 3 * minimumWindowPoints);
    }

    int platformBegin = -1;
    int platformEnd = -1;
    int runBegin = -1;
    int longestRun = 0;
    const int platformSearchBegin = storageEnd + minimumWindowPoints;
    const int platformSearchEnd = pointCount - minimumWindowPoints;
    for(int i = platformSearchBegin; i <= platformSearchEnd; ++i) {
        bool platformLike = i < platformSearchEnd &&
                            isFiniteNumber(derivativeSlope[i]) &&
                            qAbs(derivativeSlope[i]) <= 0.15;
        if(platformLike && runBegin < 0) {
            runBegin = i;
        }
        if((!platformLike || i == platformSearchEnd) && runBegin >= 0) {
            int runEnd = i;
            if(runEnd - runBegin > longestRun) {
                longestRun = runEnd - runBegin;
                platformBegin = runBegin;
                platformEnd = runEnd;
            }
            runBegin = -1;
        }
    }

    if(longestRun < minimumWindowPoints) {
        // 没有清晰平台时使用固定对数时间比例，不对目标曲线强行作流态解释。
        storageEnd = qBound(
            minimumWindowPoints,
            static_cast<int>(pointCount * 0.25 + 0.5),
            pointCount - 3 * minimumWindowPoints);
        platformBegin = qBound(
            storageEnd + minimumWindowPoints,
            static_cast<int>(pointCount * 0.45 + 0.5),
            pointCount - 2 * minimumWindowPoints);
        platformEnd = qBound(
            platformBegin + minimumWindowPoints,
            static_cast<int>(pointCount * 0.75 + 0.5),
            pointCount - minimumWindowPoints);
    }

    // 边界两侧各重叠三个采样点，降低流态边界轻微识别偏差造成的选参跳变。
    const int overlapPoints = 3;
    *windowBegin << 0
                 << qMax(0, storageEnd - overlapPoints)
                 << qMax(0, platformBegin - overlapPoints)
                 << qMax(0, platformEnd - overlapPoints);
    *windowEnd << qMin(pointCount, storageEnd + overlapPoints)
               << qMin(pointCount, platformBegin + overlapPoints)
               << qMin(pointCount, platformEnd + overlapPoints)
               << pointCount;
}

// 求解选中参数对应的阻尼正规方程。参数最多八维，使用带部分主元的
// 高斯消元即可处理该小矩阵，并在主元退化时明确返回失败。
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];
        }
    }
}

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("")
{
    // LM 对象只初始化信赖域运行状态和求解器临时目录。
    initializeTemporaryDirectory();
    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);
}

// ===== 临时目录工具 =====
//
// 真实求解器 DLL 和自动拟合过程会产生中间文件，因此每次创建独立的
// autofit_temp_<pid>_<timestamp> 目录。退出时只删除本类创建的目录，启动时顺便清理
// 旧进程遗留的 autofit_temp_*，避免长期调试后应用目录被临时文件堆满。

void nmCalculationAutoFitLM::initializeTemporaryDirectory()
{
    // 先清理历史遗留目录，再为本次对象创建唯一目录。
    // 目录名包含进程 ID 和毫秒时间戳，通常足够唯一；counter 是极端重名时的兜底。
    cleanupOldTemporaryDirectories();

    QString timestamp = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss_zzz");
    QString processId = QString::number(QCoreApplication::applicationPid());

    m_tempDirectory = QApplication::applicationDirPath() +
                      "/autofit_temp_" + processId + "_" + timestamp;

    int counter = 0;
    QString originalPath = m_tempDirectory;

    while(QDir(m_tempDirectory).exists() && counter < 100) {
        m_tempDirectory = originalPath + "_" + QString::number(counter);
        counter++;
    }

    if(QDir().mkpath(m_tempDirectory)) {
        DEBUG_OUT(QString("Initialized temp directory: %1").arg(m_tempDirectory));
    } else {
        DEBUG_OUT(QString("Warning: Failed to create temp directory: %1").arg(m_tempDirectory));
        m_tempDirectory = QApplication::applicationDirPath();
    }
}

void nmCalculationAutoFitLM::cleanupTemporaryDirectory()
{
    // 析构或用户停止时调用。删除失败一般是文件仍被 DLL/系统占用，
    // 这里只记录 debug 信息，不让清理失败影响 UI 退出。
    if(QDir(m_tempDirectory).exists()) {
        if(removeDirectoryRecursively(m_tempDirectory)) {
            DEBUG_OUT("Temp directory cleaned up successfully");
        } else {
            DEBUG_OUT("Warning: Failed to clean up temp directory completely");
        }
    }
}

bool nmCalculationAutoFitLM::removeDirectoryRecursively(const QString& path)
{
    // Qt 旧版本没有统一可用的 removeRecursively 行为时，用本函数递归删除。
    // 调用方传入的是本类创建的临时目录或旧 autofit_temp_* 目录。
    QDir dir(path);

    if(!dir.exists()) {
        return true;
    }

    // 递归删除子目录和文件。包含 Hidden，避免隐藏中间文件阻塞目录删除。
    QFileInfoList entries = dir.entryInfoList(QDir::NoDotAndDotDot | QDir::AllEntries | QDir::Hidden);
    bool allRemoved = true;

    for(int i = 0; i < entries.size(); ++i) {
        const QFileInfo& entry = entries[i];

        if(entry.isDir()) {
            if(!removeDirectoryRecursively(entry.absoluteFilePath())) {
                allRemoved = false;
            }
        } else {
            QFile file(entry.absoluteFilePath());

            // 处理只读文件。某些求解器输出可能带只读属性，删除前先补写权限。
            if(!file.permissions().testFlag(QFile::WriteUser)) {
                file.setPermissions(file.permissions() | QFile::WriteUser);
            }

            if(!file.remove()) {
                DEBUG_OUT(QString("Failed to remove file: %1").arg(entry.absoluteFilePath()));
                allRemoved = false;
            }
        }
    }

    // 删除目录本身
    if(allRemoved) {
        return dir.rmdir(path);
    }

    return false;
}

void nmCalculationAutoFitLM::cleanupOldTemporaryDirectories()
{
    // 应用启动或新建自动拟合对象时清理旧目录。
    // 不删除当前进程 ID 对应目录，防止同进程内多个拟合对象或正在运行的求解器被误删。
    QString appDir = QApplication::applicationDirPath();
    QDir dir(appDir);

    // 获取当前进程ID，避免删除当前进程可能使用的目录
    QString currentProcessId = QString::number(QCoreApplication::applicationPid());

    // 查找所有以 "autofit_temp_" 开头的目录
    QStringList filters;
    filters << "autofit_temp_*";
    QFileInfoList tempDirs = dir.entryInfoList(filters, QDir::Dirs | QDir::NoDotAndDotDot);

    if(tempDirs.isEmpty()) {
        DEBUG_OUT("No old temporary directories found");
        return;
    }

    DEBUG_OUT(QString("Found %1 potential old temporary directories to clean up").arg(tempDirs.size()));

    int cleanedCount = 0;
    int failedCount = 0;
    int skippedCount = 0;

    for(int i = 0; i < tempDirs.size(); ++i) {
        const QFileInfo& dirInfo = tempDirs[i];
        QString dirPath = dirInfo.absoluteFilePath();
        QString dirName = dirInfo.fileName();

        // 检查是否是当前进程的目录（虽然理论上不应该存在，但为了安全起见）
        if(dirName.contains("_" + currentProcessId + "_")) {
            DEBUG_OUT(QString("Skipping current process directory: %1").arg(dirName));
            skippedCount++;
            continue;
        }

        // 尝试删除目录
        DEBUG_OUT(QString("Attempting to remove old temp directory: %1").arg(dirName));

        if(removeDirectoryRecursively(dirPath)) {
            DEBUG_OUT(QString("Successfully cleaned up: %1").arg(dirName));
            cleanedCount++;
        } else {
            DEBUG_OUT(QString("Failed to clean up: %1 (may be in use by another process)").arg(dirName));
            failedCount++;
        }
    }

    // 输出清理统计信息
    DEBUG_OUT(QString("Old temp directories cleanup summary: %1 removed, %2 failed, %3 skipped")
              .arg(cleanedCount).arg(failedCount).arg(skippedCount));
}


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

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()
{
    // LM 轨迹独立保存在应用目录，不依赖机器学习目录或外部评分进程。
    closeTraceFile();
    m_traceRunId = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss_zzz");
    QDir traceDir(QDir(QApplication::applicationDirPath()).absoluteFilePath("autofit_lm_trace"));

    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"
          << "w0_begin_time"
          << "w0_end_time"
          << "w0_pressure_loss"
          << "w0_derivative_loss"
          << "w0_loss"
          << "w1_begin_time"
          << "w1_end_time"
          << "w1_pressure_loss"
          << "w1_derivative_loss"
          << "w1_loss"
          << "w2_begin_time"
          << "w2_end_time"
          << "w2_pressure_loss"
          << "w2_derivative_loss"
          << "w2_loss"
          << "w3_begin_time"
          << "w3_end_time"
          << "w3_pressure_loss"
          << "w3_derivative_loss"
          << "w3_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"
         << "coverage";

    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\": 2,\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);
        for(int window = 0;
                window < TRUST_REGION_TIME_WINDOW_COUNT;
                ++window) {
            cols << (window < objectiveBreakdown->timeWindowStart.size()
                     ? traceNumber(objectiveBreakdown->timeWindowStart[window])
                     : QString())
                 << (window < objectiveBreakdown->timeWindowEndTime.size()
                     ? traceNumber(objectiveBreakdown->timeWindowEndTime[window])
                     : QString())
                 << (window < objectiveBreakdown->timeWindowPressureLoss.size()
                     ? traceNumber(
                         objectiveBreakdown->timeWindowPressureLoss[window])
                     : QString())
                 << (window < objectiveBreakdown->timeWindowDerivativeLoss.size()
                     ? traceNumber(
                         objectiveBreakdown->timeWindowDerivativeLoss[window])
                     : QString())
                 << (window < objectiveBreakdown->timeWindowLoss.size()
                     ? traceNumber(objectiveBreakdown->timeWindowLoss[window])
                     : QString());
        }
        cols << 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)
             << traceNumber(objectiveBreakdown->coverage);
    } else {
        for(int i = 0; i < 34; ++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();
        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);

            // 裂缝参数改变时恢复最终已接受参数对应的 PEBI 缓存。
            const bool fractureGridParameterSelected =
                (m_parameterSelected.size() > 8 && m_parameterSelected[8]) ||
                (m_parameterSelected.size() > 9 && m_parameterSelected[9]);
            if(fractureGridParameterSelected) {
                nmCalculationPebiGrid* pebiGrid = nmCalculationPebiGrid::getInstance();
                if(!pebiGrid || !pebiGrid->generateOutputPara()) {
                    throw std::runtime_error("Failed to refresh final fracture parameters");
                }
            }

            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();
    emit fittingFinished(success, message);
    cleanupTemporaryDirectory();
    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;
    // 误差下降至少达到绝对 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 维残差，列对应用户勾选的参数。窗口选参直接
    // 使用对应残差行和 Jacobian，不再维护额外的上下、左右或形状梯度。
    QVector<QVector<double> > jacobian;
    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));
        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;
                }

                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);

            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;
            }
        }

        // 主目标采用 0.5*||r||^2，其对参数的梯度为 J^T*r。这里不再叠加
        // 窗口损失，保证分段只改变选参，不改变真实接受目标和 LM 方向。
        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];
            }
        }

        // 每轮最多联合调整三个灵敏参数，并继续剔除全曲线响应过度共线的列。
        auto selectColumnsByScore = [&](const QVector<double>& scores)
                -> QVector<int> {
            QVector<int> selected;
            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] ||
                            column >= scores.size() ||
                            !isFiniteNumber(scores[column]) ||
                            scores[column] <= bestScore) {
                        continue;
                    }
                    bool excessivelyCorrelated = false;
                    for(int i = 0; i < selected.size(); ++i) {
                        if(trustRegionJacobianColumnCorrelation(
                                    jacobian, column, selected[i]) >
                                columnCorrelationLimit) {
                            excessivelyCorrelated = true;
                            break;
                        }
                    }
                    if(!excessivelyCorrelated) {
                        bestColumn = column;
                        bestScore = scores[column];
                    }
                }
                if(bestColumn < 0 || bestScore <= 1.0e-12) {
                    break;
                }
                selected.append(bestColumn);
                alreadyConsidered[bestColumn] = true;
            }
            return selected;
        };

        QVector<int> windowOrder;
        for(int window = 0;
                window < TRUST_REGION_TIME_WINDOW_COUNT;
                ++window) {
            windowOrder.append(window);
        }
        std::sort(windowOrder.begin(), windowOrder.end(),
                  [&](int left, int right) -> bool {
            double leftLoss = left < current.breakdown.timeWindowLoss.size()
                              ? current.breakdown.timeWindowLoss[left] : -1.0;
            double rightLoss = right < current.breakdown.timeWindowLoss.size()
                               ? current.breakdown.timeWindowLoss[right] : -1.0;
            return leftLoss > rightLoss;
        });

        QVector<int> selectedColumns;
        int dominantWindow = TRUST_REGION_TOTAL_WINDOW;
        const int residualPointCount = jacobian.size() / 2;
        bool windowsValid = jacobian.size() % 2 == 0 &&
            current.breakdown.timeWindowBegin.size() ==
                TRUST_REGION_TIME_WINDOW_COUNT &&
            current.breakdown.timeWindowEnd.size() ==
                TRUST_REGION_TIME_WINDOW_COUNT &&
            current.breakdown.timeWindowLoss.size() ==
                TRUST_REGION_TIME_WINDOW_COUNT;

        // score 是该参数在窗口内按一维最小二乘预计可消除的残差能量。
        // 它只用于回答“调整哪些参数”，正负方向仍由后面的完整 LM 方程决定。
        if(windowsValid) {
            for(int order = 0;
                    order < windowOrder.size() && selectedColumns.isEmpty();
                    ++order) {
                int window = windowOrder[order];
                int begin = qBound(
                    0, current.breakdown.timeWindowBegin[window],
                    residualPointCount);
                int end = qBound(
                    begin, current.breakdown.timeWindowEnd[window],
                    residualPointCount);
                int pointCount = end - begin;
                if(pointCount <= 0 ||
                        !isFiniteNumber(
                            current.breakdown.timeWindowLoss[window])) {
                    continue;
                }

                QVector<double> windowScores(dimensions, 0.0);
                for(int column = 0; column < dimensions; ++column) {
                    if(!jacobianColumnValid[column]) {
                        continue;
                    }
                    double gradient = 0.0;
                    double curvature = 0.0;
                    for(int point = begin; point < end; ++point) {
                        int pressureRow = point;
                        int derivativeRow = residualPointCount + point;
                        gradient +=
                            jacobian[pressureRow][column] *
                            current.breakdown.residualVector[pressureRow] +
                            jacobian[derivativeRow][column] *
                            current.breakdown.residualVector[derivativeRow];
                        curvature +=
                            jacobian[pressureRow][column] *
                            jacobian[pressureRow][column] +
                            jacobian[derivativeRow][column] *
                            jacobian[derivativeRow][column];
                    }
                    windowScores[column] =
                        gradient * gradient /
                        (curvature + 1.0e-20) / pointCount;
                }
                selectedColumns = selectColumnsByScore(windowScores);
                if(!selectedColumns.isEmpty()) {
                    dominantWindow = window;
                }
            }
        }

        // 所有窗口都缺少可用灵敏方向时退回完整残差梯度，避免分段识别异常
        // 阻断细搜；候选仍然只按全曲线 total 是否下降来接受。
        if(selectedColumns.isEmpty()) {
            QVector<double> totalScores(dimensions, 0.0);
            for(int column = 0; column < dimensions; ++column) {
                totalScores[column] = qAbs(totalGradient[column]);
            }
            selectedColumns = selectColumnsByScore(totalScores);
        }

        // 当前局部没有可用方向时先缩小半径并重建灵敏度；只有已经在最小
        // 半径完整重建后仍无方向，才把它判定为局部最优。
        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 严格比较决定。
        updateTrustRegionJacobian(
            &jacobian,
            current.breakdown.residualVector,
            candidate.breakdown.residualVector,
            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 windowName = trustRegionWindowName(dominantWindow);

        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(windowName)
                       : QString("rejected_%1").arg(windowName),
                       &candidate.breakdown);

        QString windowDisplayName = windowName;
        if(dominantWindow == TRUST_REGION_WELLBORE_STORAGE_WINDOW) {
            windowDisplayName = tr("wellbore-storage interval");
        } else if(dominantWindow == TRUST_REGION_TRANSITION_WINDOW) {
            windowDisplayName = tr("early transition interval");
        } else if(dominantWindow == TRUST_REGION_IARF_WINDOW) {
            windowDisplayName = tr("IARF interval");
        } else if(dominantWindow == TRUST_REGION_LATE_WINDOW) {
            windowDisplayName = tr("late-time interval");
        } else {
            windowDisplayName = tr("total error");
        }

        emit logMessageGenerated(
            tr("Iteration %1: focus=%2, parameters=%3, error=%4, result=%5")
            .arg(iteration + 1)
            .arg(windowDisplayName)
            .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 和裂缝半长都通过 PEBI 裂缝数组传入，不属于每次求解都会重新组装的 Base/CS 参数。
		// 勾选任一裂缝参数时刷新网格输出，保证本次真实试算使用新的导流能力和端点坐标。
		const bool fractureGridParameterSelected =
			(m_parameterSelected.size() > 8 && m_parameterSelected[8]) ||
			(m_parameterSelected.size() > 9 && m_parameterSelected[9]);
		if(fractureGridParameterSelected) {
			nmCalculationPebiGrid* pebiGrid = nmCalculationPebiGrid::getInstance();
			if(!pebiGrid || !pebiGrid->generateOutputPara()) {
				DEBUG_OUT(QString("%1: Call #%2 - Failed to refresh PEBI fracture parameters")
				          .arg(funcName).arg(callCount));
				return 1e10;
			}
		}

        // 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 double minimumCoverage = 0.95;
    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);
        };

        // 覆盖率通过后若只缺少首尾少量点，用模拟曲线自身的端点斜率作短距离
        // 双对数外推。该外推只用于主损失的固定网格，不参与水平配准搜索。
        auto extrapolateEndpointLogValue = [valueFloor](
                const QVector<QPointF>& curve,
                double x,
                double* value) -> bool {
            if(!value || curve.size() < 2 || x <= 0.0) {
                return false;
            }
            int left = x < curve.first().x()
                       ? 0
                       : curve.size() - 2;
            int right = left + 1;
            double leftLogX = qLn(curve[left].x());
            double rightLogX = qLn(curve[right].x());
            double denominator = rightLogX - leftLogX;
            if(qAbs(denominator) <= 1.0e-12) {
                return false;
            }
            double leftLogY =
                qLn(qMax(qAbs(curve[left].y()), valueFloor));
            double rightLogY =
                qLn(qMax(qAbs(curve[right].y()), valueFloor));
            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;
        }

        QVector<double> commonX(numPoints);
        QVector<double> commonLogX(numPoints);
        QVector<double> targetLogPressure(numPoints);
        QVector<double> targetLogDerivative(numPoints);
        const double targetLogMinX = qLn(targetMinX);
        const double targetLogMaxX = qLn(targetMaxX);

        // 固定使用目标曲线的完整 log-time 网格，候选之间不会因采样点不同而失去可比性。
        for(int i = 0; i < numPoints; ++i) {
            double logX = targetLogMinX +
                          static_cast<double>(i) *
                          (targetLogMaxX - targetLogMinX) /
                          (numPoints - 1);
            commonLogX[i] = logX;
			// 首尾直接使用原始端点，避免 exp(log(t)) 的舍入误差越过严格插值边界。
			if(i == 0) {
				commonX[i] = targetMinX;
			} else if(i == numPoints - 1) {
				commonX[i] = targetMaxX;
			} 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());
        int firstSupported = -1;
        int lastSupported = -1;
        int supportedCount = 0;

        // 残差定义为“模拟减目标”：正值表示模拟曲线偏高，负值表示偏低。
        for(int i = 0; i < numPoints; ++i) {
            if(commonX[i] < resultMinX || commonX[i] > resultMaxX) {
                continue;
            }

            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];
            ++supportedCount;
            if(firstSupported < 0) {
                firstSupported = i;
            }
            lastSupported = i;
        }

        // 覆盖率同时约束“有效点数量”和“连续时间跨度”。取两者较小值可避免
        // 点数很多但只集中在局部时段的候选被误认为覆盖充分。
        AutoFitObjectiveBreakdownLM breakdown;
        breakdown.coverage = supportedCount > 0
                             ? static_cast<double>(supportedCount) /
                               numPoints
                             : 0.0;
        if(firstSupported >= 0 && lastSupported >= firstSupported) {
            double targetSpan =
                qMax(1.0e-12, targetLogMaxX - targetLogMinX);
            double coveredSpan =
                commonLogX[lastSupported] - commonLogX[firstSupported];
            breakdown.coverage = qMin(
                breakdown.coverage,
                qMax(0.0, coveredSpan / targetSpan));
        }

        // 覆盖率只判断候选是否有效，不再加入固定惩罚，避免所有误差被整体抬高。
        if(breakdown.coverage < minimumCoverage) {
            breakdown.total = invalidLoss;
            m_lastObjectiveBreakdown = breakdown;
            return invalidLoss;
        }

        // 通过门槛后最多只缺少首尾少量目标点。按模拟曲线端点趋势补齐后，
        // 每个候选仍在固定 80 点上计算均方根误差，不能靠少算难拟合端点获益。
        for(int i = 0; i < numPoints; ++i) {
            if(isFiniteNumber(pressureResidual[i]) &&
                    isFiniteNumber(derivativeResidual[i])) {
                continue;
            }
            double resultLogPressure = 0.0;
            double resultLogDerivative = 0.0;
            if(!extrapolateEndpointLogValue(
                        resultPressure, commonX[i],
                        &resultLogPressure) ||
                    !extrapolateEndpointLogValue(
                        resultDerivative, commonX[i],
                        &resultLogDerivative)) {
                return invalidLoss;
            }
            pressureResidual[i] =
                resultLogPressure - targetLogPressure[i];
            derivativeResidual[i] =
                resultLogDerivative - targetLogDerivative[i];
        }

        // 在指定中心附近计算普通均方根误差。
        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]);
        }

        // 窗口只由目标曲线确定，因此所有候选使用完全相同的时间分段。每个窗口
        // 同时统计压力、导数和二者等权联合 RMS，供细搜选择当前关注区间。
        buildTrustRegionTimeWindows(
            commonX,
            targetLogPressure,
            targetLogDerivative,
            &breakdown.timeWindowBegin,
            &breakdown.timeWindowEnd);
        if(breakdown.timeWindowBegin.size() !=
                TRUST_REGION_TIME_WINDOW_COUNT ||
                breakdown.timeWindowEnd.size() !=
                TRUST_REGION_TIME_WINDOW_COUNT) {
            return invalidLoss;
        }
        for(int window = 0;
                window < TRUST_REGION_TIME_WINDOW_COUNT;
                ++window) {
            int begin = breakdown.timeWindowBegin[window];
            int end = breakdown.timeWindowEnd[window];
            if(begin < 0 || end <= begin || end > numPoints) {
                return invalidLoss;
            }
            double windowPressureLoss = rmse(
                pressureResidual, begin, end);
            double windowDerivativeLoss = rmse(
                derivativeResidual, begin, end);
            if(!isFiniteNumber(windowPressureLoss) ||
                    !isFiniteNumber(windowDerivativeLoss)) {
                return invalidLoss;
            }
            breakdown.timeWindowStart.append(commonX[begin]);
            breakdown.timeWindowEndTime.append(commonX[end - 1]);
            breakdown.timeWindowPressureLoss.append(windowPressureLoss);
            breakdown.timeWindowDerivativeLoss.append(windowDerivativeLoss);
            breakdown.timeWindowLoss.append(qSqrt(
                0.5 * windowPressureLoss * windowPressureLoss +
                0.5 * windowDerivativeLoss * windowDerivativeLoss));
        }

        const double logGridStep =
            (targetLogMaxX - targetLogMinX) /
            (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, coverage=%8, total=%9")
            .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.coverage, 0, 'f', 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");

    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);

        if(m_targetWellName.isEmpty()) {
            DEBUG_OUT("Target well name is empty - target-only solver cannot start");
            delete dllTask;
            --m_evaluationInProgress;
            return result;
        }

        dllTask->setAutoFitTargetWell(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()
{
    // 不设置目标井名，任务按原完整模式保存全部井和网格结果。
    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);
    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();
    --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");
    }
}