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

#include <QApplication>
#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 const int kAutoFitIntervalsPerDecade = 20;

static int autoFitShapeLag(int pointCount)
{
    // 保持斜率跨度约为完整对数时间范围的 10%，四舍五入到间隔数且至少跨一个间隔。
    return qMax(1, qRound((pointCount - 1) * 0.10));
}

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、Dfc 等跨数量级参数被线性尺度支配；
// 表皮可为负数，使用线性坐标。
static bool useTrustRegionLogScale(int parameterIndex, double lower, double upper)
{
    return parameterIndex != 1 &&
           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)
{
    coordinate = qMax(0.0, qMin(1.0, coordinate));

    // 端点直接返回原始边界，避免对数或线性逆变换的舍入误差造成越界。
    if(coordinate <= 0.0) {
        return lower;
    }
    if(coordinate >= 1.0) {
        return upper;
    }

    double value;
    if(useTrustRegionLogScale(parameterIndex, lower, upper)) {
        value = qExp(qLn(lower) + coordinate * (qLn(upper) - qLn(lower)));
    } else {
        value = lower + coordinate * (upper - lower);
    }

    // 内部点转换后也限制到实际上下界，保证候选参数通过严格的边界检查。
    return qMax(lower, qMin(upper, value));
}

// 一次真实求解的完整快照。除了参数和总误差，还保存内部坐标、诊断分量和
// 双对数曲线，因此拒绝候选后可以完整恢复上一个已接受工作点。
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)
{
    // 压力和导数共用采样时间，两个通道点数相同且各至少三个点。
    if(!breakdown.valid || breakdown.residualVector.size() < 6 ||
            breakdown.residualVector.size() % 2 != 0) {
        return false;
    }

    for(int i = 0; i < breakdown.residualVector.size(); ++i) {
        if(!isFiniteNumber(breakdown.residualVector[i])) {
            return false;
        }
    }
    const int pointCount = breakdown.residualVector.size() / 2;
    const int slopeCount = pointCount - autoFitShapeLag(pointCount);
    if(breakdown.shapeResiduals.size() != 2 * slopeCount || !isFiniteNumber(breakdown.shapeLoss)) return false;
    for(int i = 0; i < breakdown.shapeResiduals.size(); ++i) {
        if(!isFiniteNumber(breakdown.shapeResiduals[i])) return false;
    }
    return true;
}

// 数值和整体形状残差共用一次求解，联合缓存供两个 LM 目标复用。
static QVector<double> trustRegionFullResidual(const AutoFitObjectiveBreakdownLM& objective)
{
    QVector<double> residual = objective.residualVector;
    residual += objective.shapeResiduals;
    return residual;
}

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

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

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

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

// LM 正规方程所需的 J^T*J 与 J^T*r，使用当前归一化参数坐标。
struct TrustRegionNormalEquations
{
    QVector<QVector<double> > matrix;
    QVector<double> gradient;

    explicit TrustRegionNormalEquations(int dimensions = 0)
        : matrix(dimensions, QVector<double>(dimensions, 0.0))
        , gradient(dimensions, 0.0)
    {}
};

static TrustRegionNormalEquations buildTrustRegionNormalEquations(
        const QVector<QVector<double> >& jacobian,
        const QVector<double>& residual,
        const QVector<bool>& columnValid)
{
    // 残差和 J 已包含压力、导数及采样权重，直接累加整体矩阵，避免重复计权。
    const int dimensions = columnValid.size();
    TrustRegionNormalEquations equations(dimensions);
    for(int row = 0; row < residual.size(); ++row) {
        for(int p = 0; p < dimensions; ++p) {
            if(!columnValid[p]) continue;
            equations.gradient[p] += jacobian[row][p] * residual[row];
            for(int q = 0; q <= p; ++q) {
                if(columnValid[q])
                    equations.matrix[p][q] += jacobian[row][p] * jacobian[row][q];
            }
        }
    }
    for(int p = 0; p < dimensions; ++p) {
        for(int q = 0; q < p; ++q)
            equations.matrix[q][p] = equations.matrix[p][q];
    }
    return equations;
}

static bool buildTrustRegionStep(
        const TrustRegionNormalEquations& global,
        const QVector<int>& selected,
        const QVector<double>& coordinates,
        double damping, double trustRadius, double minimumStep,
        QVector<double>* step, double* predictedReduction)
{
    // 复用原有全局 LM 方程、阻尼缩放、退化回退、信赖半径和边界投影。
    const int dimensions = coordinates.size();
    const int count = selected.size();
    QVector<QVector<double> > normal(count, QVector<double>(count, 0.0));
    QVector<double> rhs(count, 0.0);
    for(int i = 0; i < count; ++i) {
        const int p = selected[i];
        rhs[i] = -global.gradient[p];
        for(int j = 0; j < count; ++j) {
            normal[i][j] = global.matrix[p][selected[j]];
        }
        normal[i][i] += damping * qMax(1.0e-10, normal[i][i]);
    }
    QVector<double> solution;
    bool solved = solveTrustRegionLinearSystem(normal, rhs, &solution);
    step->fill(0.0, dimensions);
    if(solved) {
        for(int i = 0; i < count; ++i) {
            (*step)[selected[i]] = solution[i];
        }
    }
    auto projectAndPredict = [&]() -> bool {
        const double norm = qSqrt(trustRegionSquaredNorm(*step));
        if(!isFiniteNumber(norm) || norm < minimumStep) return false;
        if(norm > trustRadius) {
            for(int p = 0; p < dimensions; ++p) (*step)[p] *= trustRadius / norm;
        }
        for(int p = 0; p < dimensions; ++p) {
            (*step)[p] = qBound(0.0, coordinates[p] + (*step)[p], 1.0) - coordinates[p];
        }
        // 只比较投影后可执行步长的下降，阻尼项不属于真实拟合目标。
        *predictedReduction = -trustRegionDotProduct(global.gradient, *step);
        for(int p = 0; p < dimensions; ++p) {
            for(int q = 0; q < dimensions; ++q) {
                *predictedReduction -= 0.5 * (*step)[p] * global.matrix[p][q] * (*step)[q];
            }
        }
        return isFiniteNumber(*predictedReduction) && *predictedReduction > 1.0e-14 &&
            qSqrt(trustRegionSquaredNorm(*step)) >= minimumStep;
    };
    if(solved && projectAndPredict()) return true;

    // 联合解可能被边界投影破坏；此时尝试可行梯度方向，不直接判为没有下降方向。
    step->fill(0.0, dimensions);
    for(int i = 0; i < count; ++i) {
        const int p = selected[i];
        double direction = -global.gradient[p];
        if((coordinates[p] <= minimumStep && direction < 0.0) ||
                (coordinates[p] >= 1.0 - minimumStep && direction > 0.0)) direction = 0.0;
        (*step)[p] = direction;
    }
    const double norm = qSqrt(trustRegionSquaredNorm(*step));
    if(!isFiniteNumber(norm) || norm < minimumStep) return false;
    for(int p = 0; p < dimensions; ++p) {
        (*step)[p] = qBound(0.0, coordinates[p] + (*step)[p] * trustRadius / norm, 1.0) - coordinates[p];
    }
    const double linearReduction = -trustRegionDotProduct(global.gradient, *step);
    double curvature = 0.0;
    for(int p = 0; p < dimensions; ++p) {
        for(int q = 0; q < dimensions; ++q) {
            curvature += (*step)[p] * global.matrix[p][q] * (*step)[q];
        }
    }
    if(!isFiniteNumber(linearReduction) || !isFiniteNumber(curvature) || linearReduction <= 0.0) return false;
    // 沿投影梯度最小化局部二次模型，只缩步，不突破信赖域或参数边界。
    const double scale = curvature > 0.0 ? qMin(1.0, linearReduction / curvature) : 1.0;
    for(int p = 0; p < dimensions; ++p) (*step)[p] *= scale;
    return projectAndPredict();
}

// 每次得到有效真实候选后，使用满足最新割线条件的秩一修正更新完整残差
// 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"
          << "Dfc"
          << "fractureHalfLength";
    return names;
}

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

// 请求停止并限时等待，然后关闭轨迹文件、清理临时目录和信号连接。
// 等待超时后仅清除运行标志，不代表外层同步调用已经退出。
nmCalculationAutoFitLM::~nmCalculationAutoFitLM()
{
    if(m_isRunning) {
        m_shouldStop = true;
        int waitCount = 0;
        while(m_isRunning && waitCount < 100) {
            QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, 50);
            msleep(50);
            waitCount++;
        }
        if(m_isRunning) {
            DEBUG_OUT("Force stopping LM fitting after timeout");
            m_isRunning = false;
        }
    }

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

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

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

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


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

void nmCalculationAutoFitLM::setTargetLogLogData(const QVector<QVector<double> >& targetData)
{
    // 目标曲线由界面层从目标井 history log-log 传入。
    // 约定 targetData[0]=time，targetData[1]=pressure，targetData[2]=pressure derivative。
    m_targetLogLogData = targetData;
    m_comparisonTimeMin = 0.0;
    m_comparisonTimeMax = 0.0;
    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()
{
    if(m_isFinalizing) {
        DEBUG_OUT(
                    tr("The final full-field calculation is already running"));
        return;
    }

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

        // 此槽由求解等待循环派发，必须立即返回，外层循环才能转发取消请求。
        if(m_evaluationInProgress > 0) {
            emit logMessageGenerated(tr("Waiting for current solver evaluation to stop..."));
        }

        DEBUG_OUT(tr("LM automatic fitting stop request processed"));
    } else {
        DEBUG_OUT(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;
}

int nmCalculationAutoFitLM::getSolverCallCount() const
{
    // 包含失败调用、重试和最终完整求解，不包含未进入求解器的评价。
    return m_solverCallCount;
}

void nmCalculationAutoFitLM::recordSolverCall()
{
    // 信号来自求解线程，使用原子累加避免跨线程读写竞争。
    m_solverCallCount.fetchAndAddOrdered(1);
}

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

double nmCalculationAutoFitLM::getBestFitness() const
{
    // 当前已接受工作点的总误差；高度和形状阶段不保证总误差单调下降。
    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_comparisonTimeMin = 0.0;
    m_comparisonTimeMax = 0.0;
    m_currentIteration = -1; // 尚未执行迭代；对外日志使用索引加一。
    m_totalEvaluations = 0;
    m_solverCallCount = 0;
    m_successfulEvaluations = 0;
    m_lastError.clear();
    m_initialValues.clear();
    m_userInitialSolution.clear();
    m_userInitialFitness = 1e10;
    m_hasValidUserSolution = false;
    m_traceMetaFilePath.clear();
    m_isFinalizing = false;

    DEBUG_OUT("LM optimizer reset");
}

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

void nmCalculationAutoFitLM::initializeTraceFile()
{
    closeTraceFile();
    m_traceFilePath.clear();
    m_traceMetaFilePath.clear();

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

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

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

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

    writeTraceHeader();
    writeTraceMetaFile();
    DEBUG_OUT(tr("LM fitting trace: %1").arg(m_traceFilePath));
    if(!m_traceMetaFilePath.isEmpty()) {
        DEBUG_OUT(
            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"
         << "Dfc"
         << "fractureHalfLength"
         << "solver_objective"
         << "solver_success"
         << "elapsed_ms"
         << "decision"
         << "enabled_param_indices"
         << "pressure_loss"
         << "derivative_loss"
         << "vertical_common_bias"
         << "vertical_loss"
         << "vertical_reliable"
         << "horizontal_physical_shift"
         << "horizontal_loss"
         << "horizontal_reliable"
         << "shape_loss"
         << "late_trend_loss"
         << "late_slope_bias"
         << "late_trend_reliable"
         << "registration_ambiguous";

    cols << "sampling_mode" << "sampling_points"
         << "pressure_vertical_bias" << "derivative_vertical_bias";
    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\": 56,\n";
    out << "  \"strategy\": \"permeability_height_then_joint_lm_total_for_vertical_well_else_shape_then_total\",\n";
    out << "  \"height_acceptance\": \"reliable_vertical_loss_decrease; no_shape_or_total_loss_constraint\",\n";
    out << "  \"shape_metric\": \"pressure_and_derivative_log_slopes_shared_sampling_grid_lag_round_10_percent_min_1\",\n";
    out << "  \"shape_parameter_selection\": \"all_valid_free_columns_joint_LM; no_single_parameter_sweep_or_wellbore_recheck\",\n";
    out << "  \"shape_acceptance\": \"strict_full_shape_loss_decrease; no_total_loss_constraint\",\n";
    out << "  \"lm_step_policy\": \"same_joint_LM_damping_trust_radius_secant_updates_and_sensitivity_rebuilds_in_both_phases\",\n";
    out << "  \"lm_phase_switch\": \"shape_target_reached_or_confirmed_stagnation_or_phase_evaluation_budget_exhausted_then_total_if_iterations_remain; shared_iteration_limit_or_user_stop_or_consecutive_solver_failure_abort\",\n";
    out << "  \"lm_phase_budget\": \"shared_max_iterations_including_height; each_LM_phase_retains_max(dimensions+2,20,3*max_iterations)_evaluations\",\n";
    out << "  \"lm_effective_improvement\": \"max(0.01,0.10*baseline_stage_error); after_sensitivity_rebuild_reset_baseline_to_accepted_probe_error_if_improved_else_rebuild_start_error; reset_ineffective_count_preserving_stagnation_confirmation; 3_ineffective_steps_trigger_stagnation_confirmation\",\n";
    out << "  \"iteration_count_scope\": \"height_shape_and_total_share_max_iterations; rejected_and_failed_candidates_count; sensitivity_evaluations_are_not_iterations; trace_iteration_is_zero_based_cumulative\",\n";
    out << "  \"total_stage_shape_constraint\": false,\n";
    out << "  \"total_parameter_selection\": \"all_valid_free_columns\",\n";
    out << "  \"skin_difference_policy\": \"local_scale_fraction_0.02_in_both_LM_phases\",\n";
    out << "  \"difference_policy\": \"coordinate_step_cap_0.04_with_half_trust_radius_in_both_LM_phases\",\n";
    out << "  \"difference_failure_policy\": \"halve_before_opposite_direction\",\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 << "  \"sampling_mode\": "
        << jsonEscape("log_uniform_per_decade") << ",\n";
    out << "  \"sampling_intervals_per_decade\": " << kAutoFitIntervalsPerDecade << ",\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";
    // 区间尚未冻结时写 null；首次有效评价后记录固定的误差比较范围。
    out << "  \"comparison_time_range\": {\n";
    out << "    \"comparison_time_min\": "
        << (m_comparisonTimeMin > 0.0 ? jsonNumber(m_comparisonTimeMin) : QString("null")) << ",\n";
    out << "    \"comparison_time_max\": "
        << (m_comparisonTimeMax > 0.0 ? jsonNumber(m_comparisonTimeMax) : QString("null")) << "\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 < 6; ++i) {
        cols << traceParamAt(fullParams, i);
    }
    cols << traceNumber(solverObjective)
         << QString::number(solverSuccess ? 1 : 0)
         << QString::number(elapsedMs)
         << csvEscape(decision)
         << csvEscape(enabledIndices.join(";"));

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

    cols << "log_uniform_per_decade"
         << (objectiveBreakdown ? QString::number(objectiveBreakdown->residualVector.size() / 2) : QString())
         << (objectiveBreakdown ? traceNumber(objectiveBreakdown->pressureVerticalBias) : QString())
         << (objectiveBreakdown ? traceNumber(objectiveBreakdown->derivativeVerticalBias) : QString());
    QTextStream out(&m_traceFile);
    out << cols.join(",") << "\n";
    m_traceFile.flush();
}

// 拟合界面统一显示普通小数；失败评价不展示惩罚值或非有限数。
static QString lmLogError(double value)
{
    return isFiniteNumber(value) && value < 1.0e9
        ? QString::number(value, 'f', 6) : QString("--");
}

// 改善率统一相对上次已接受误差计算；起点为零时不作除法。
static QString lmLogImprovement(double beforeError, double afterError)
{
    return beforeError > 0.0
        ? QString::number((beforeError - afterError) / beforeError * 100.0, 'f', 2)
        : QString("--");
}

void nmCalculationAutoFitLM::emitRunSummary(bool success, StopReasonLM finalReason)
{
    // 最终只汇总总误差，耗时由进度窗口在收到结束信号后追加。
    emit logMessageGenerated(success ? tr("[Fitting completed]") : tr("[Fitting ended]"));
    emit logMessageGenerated(tr("Stop reason: %1").arg(success
        ? getStopReasonDescription(finalReason)
        : (m_lastError.isEmpty() ? getStopReasonDescription(finalReason) : m_lastError)));
    emit logMessageGenerated(tr("Total error: %1 -> %2")
        .arg(lmLogError(m_userInitialFitness)).arg(lmLogError(m_globalBestFitness)));
    emit logMessageGenerated(tr("Solver calls: %1").arg(getSolverCallCount()));
}

void nmCalculationAutoFitLM::emitTrialLog(const QString& title,
    const AutoFitObjectiveBreakdownLM& before, const AutoFitObjectiveBreakdownLM& after,
    bool valid, bool accepted, bool shapeStage)
{
    // 只显示当前阶段目标；改善率相对上次已接受的误差计算，变差时为负值。
    if(!valid) {
        emit logMessageGenerated(tr("%1 | Solver failed\nNo valid error was returned; keep the last accepted result and reduce the step.").arg(title));
        return;
    }
    const QString status = accepted ? tr("accepted") : tr("rejected");
    const double beforeError = shapeStage ? before.shapeLoss : before.total;
    const double afterError = shapeStage ? after.shapeLoss : after.total;
    const QString objective = (shapeStage ? tr("Shape error: %1 -> %2")
        : tr("Total error: %1 -> %2"))
        .arg(lmLogError(beforeError)).arg(lmLogError(afterError));
    QString message = tr("%1 | %2 | %3 | Improvement: %4%")
        .arg(title).arg(status).arg(objective).arg(lmLogImprovement(beforeError, afterError));
    if(!accepted) message += tr("\n    Objective error did not decrease; keep the last accepted result and reduce the step.");
    emit logMessageGenerated(message);
}

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

    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();

    if(dataManager) {
        nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
        fullParams[0] = reservoirData.getPermeability().getValue().toDouble();
        fullParams[3] = reservoirData.getPorosity().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[4] = fracturedWell->getDfc().getValue().toDouble();
                    fullParams[5] = fracturedWell->getFractureHalfLength().getValue().toDouble();
                }
            } else if(pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
                nmDataHorizontalFracturedWell* fracturedWell =
                    dynamic_cast<nmDataHorizontalFracturedWell*>(pTargetWell);
                if(fracturedWell) {
                    fullParams[4] = fracturedWell->getDfc().getValue().toDouble();
                    fullParams[5] = 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 = tr("Failed to load configuration from data manager");
        return false;
    }
}

void nmCalculationAutoFitLM::loadOptimizationConfig()
{
    // LM 只读取迭代次数和目标误差，其他数值控制保留在现有算法实现中。
    // 敏感性运行只读取确认过的设置，普通拟合保留原来的配置来源。
    nmDataAutomaticFitting fittingData = m_sensitivitySettings
        ? m_sensitivitySettings->parameters
        : nmDataAnalyzeManager::getCurrentInstance()->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[6]：完整参数体系中每个参数是否参与拟合；
    // - m_parameterLower/Upper[6]：完整参数体系的搜索上下界；
    // - m_enabledParamIndices：将勾选参数向量映射回完整参数索引。
    // 敏感性运行只读取确认过的设置，普通拟合保留原来的配置来源。
    nmDataAutomaticFitting fittingData = m_sensitivitySettings
        ? m_sensitivitySettings->parameters
        : nmDataAnalyzeManager::getCurrentInstance()->getAutomaticFittingDataCopy();

    m_parameterSelected.resize(6);
    m_parameterSelected[0] = fittingData.getPermeabilitySelected();
    m_parameterSelected[1] = fittingData.getSkinSelected();
    m_parameterSelected[2] = fittingData.getWellboreStorageSelected();
    m_parameterSelected[3] = fittingData.getPorositySelected();
    m_parameterSelected[4] = fittingData.getFractureConductivitySelected();
    m_parameterSelected[5] = fittingData.getFractureHalfLengthSelected();

    m_parameterLower.resize(6);
    m_parameterUpper.resize(6);

    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.getFractureConductivityMin().getValue().toDouble();
    m_parameterUpper[4] = fittingData.getFractureConductivityMax().getValue().toDouble();

    m_parameterLower[5] = fittingData.getFractureHalfLengthMin().getValue().toDouble();
    m_parameterUpper[5] = 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::startSensitivityFitting(
        const nmSensitivityFitSettings& settings,
        const std::function<QVector<QVector<double> >(const QVector<double>&)>& evaluate,
        const std::function<bool()>& cancelled, bool verticalWell)
{
    // 同步复用原算法；返回前释放回调，避免持有任务局部对象。
    m_sensitivitySettings = &settings;
    m_sensitivityEvaluate = evaluate;
    m_sensitivityCancelled = cancelled;
    m_sensitivityVerticalWell = verticalWell;
    m_targetWellName = settings.targetWell;
    const bool succeeded = startAutoFitting();
    m_sensitivitySettings = nullptr;
    m_sensitivityEvaluate = std::function<QVector<QVector<double> >(const QVector<double>&)>();
    m_sensitivityCancelled = std::function<bool()>();
    return succeeded;
}

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

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

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

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

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

        if(m_targetLogLogData.size() < 3) {
            m_lastError = tr("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 = tr("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 = tr("Target well name is empty");
            emit logMessageGenerated(tr("ERROR: Target well name is empty"));
            return false;
        }

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

        // resetOptimizer() 会清空运行状态，因此先保存从 DataManager 提取的初始值。
        QVector<double> savedInitialValues = m_initialValues;
        resetOptimizer();
        // 输入校验通过后，为本次拟合创建新的临时计算目录。
        cleanupTemporaryDirectory();
        if(!m_sensitivitySettings && !initializeTemporaryDirectory()) {
            m_lastError = tr("Cannot create the automatic fitting temporary directory");
            emit logMessageGenerated(tr("ERROR: %1").arg(m_lastError));
            return false;
        }
        m_isRunning = true;
        m_shouldStop = false;
        m_isFinalizing = false;
        m_currentIteration = -1;
        m_initialValues = savedInitialValues;
        if(!m_sensitivitySettings) initializeTraceFile();
        // 开头只展示一次初始参数和总误差，不展示形状误差。
        emit logMessageGenerated(tr("========== LM fitting started =========="));
        emit logMessageGenerated(tr("Target well: %1").arg(m_targetWellName));
        emit logMessageGenerated(tr("Target error: %1").arg(lmLogError(m_targetError)));
        emit fittingStageChanged(0);

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

            const char* names[] = { QT_TR_NOOP("Permeability"), QT_TR_NOOP("Skin"),
                QT_TR_NOOP("Wellbore storage"), QT_TR_NOOP("Porosity"),
                QT_TR_NOOP("Fracture conductivity"), QT_TR_NOOP("Fracture half length") };
            const char* units[] = { "mD", "", "m^3/MPa", "", "mD.m", "m" };
            emit logMessageGenerated(tr("[Initial parameters]"));
            for(int i = 0; i < m_userInitialSolution.size(); ++i) {
                const int index = m_enabledParamIndices[i];
                emit logMessageGenerated(QString("%1: %2 %3").arg(tr(names[index]))
                    .arg(m_userInitialSolution[i], 0, 'g', 10).arg(tr(units[index])));
            }

            try {
                DEBUG_OUT(
                    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;
                    DEBUG_OUT(tr("Initial solution evaluation successful"));
                    emit logMessageGenerated(tr("Initial total error: %1").arg(lmLogError(m_userInitialFitness)));
                    emit bestCurveUpdated(m_targetLogLogData,
                                          m_globalBestLogLogData,
                                          0,
                                          m_globalBestFitness);
                } else {
                    m_hasValidUserSolution = false;
                    if(!m_shouldStop) emit logMessageGenerated(tr("Initial solution failed; trying the parameter-range midpoint."));
                }

                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;
                if(!m_shouldStop) emit logMessageGenerated(tr("Initial solution failed; trying the parameter-range midpoint."));
            }

            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(finalReason != LM_USER_STOPPED && finalReason != LM_OPTIMIZATION_FAILED &&
                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_isFinalizing = false;
        m_isRunning = false;
        emitRunSummary(false, LM_OPTIMIZATION_FAILED);
        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_isFinalizing = false;
        m_isRunning = false;
        emitRunSummary(false, LM_OPTIMIZATION_FAILED);
        emit fittingFinished(false, m_lastError);
        return false;
    }

    // 敏感性任务自己用最优参数收尾，不能进入普通拟合的模型回写流程。
    if(m_sensitivitySettings) {
        m_isRunning = false;
        return !m_shouldStop && !m_sensitivityCancelled() &&
               !m_globalBestPosition.isEmpty() && m_globalBestFitness < 1.0e9;
    }

    // 普通拟合只有完整求解并提交快照后才算成功。
    bool finalFullSolverSucceeded = false;
    bool finalFullSolverExecuted = false;

    if(!m_globalBestPosition.isEmpty()) {
        try {
            // Stop 只结束优化迭代；从这里开始必须用当前最优参数生成并发布正式快照。
            m_isFinalizing = true;
            emit finalizingStarted();
            DEBUG_OUT(tr("Applying optimized parameters to model..."));
            applyParametersToDataManager(m_globalBestPosition);

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

            DEBUG_OUT(
                tr("Running final full-field calculation with optimized parameters..."));
            finalFullSolverExecuted = true;
            finalFullSolverSucceeded = runFinalFullSolver();

            if(finalFullSolverSucceeded) {
                DEBUG_OUT(
                    tr("Final full-field calculation completed successfully"));
            } else {
                m_lastError =
                    tr("Optimized parameters were found, but the final full-field calculation failed");
                DEBUG_OUT(
                    tr("ERROR: Final full-field calculation failed"));
            }

            DEBUG_OUT(tr("=== Optimization Results ==="));
            DEBUG_OUT(
                tr("Final error: %1").arg(m_globalBestFitness, 0, 'e', 4));
            DEBUG_OUT(
                tr("Total iterations: %1").arg(m_currentIteration + 1));
            DEBUG_OUT(
                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);
            }
            DEBUG_OUT(finalParams);

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

    m_isFinalizing = false;
    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);
        DEBUG_OUT(tr("=== LM AUTOMATIC FITTING SUCCESSFUL ==="));
    } 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);
        DEBUG_OUT(tr("=== LM AUTOMATIC FITTING - LOCAL OPTIMUM ==="));
    } else if(finalReason == LM_MAX_ITERATIONS) {
        success = true;
        message = QString(tr("Iteration or evaluation budget reached. Best error: %1, Cumulative iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        DEBUG_OUT(tr("=== LM AUTOMATIC FITTING - MAX ITERATIONS ==="));
    } else if(finalReason == LM_USER_STOPPED) {
        success = true;
        message = QString(tr("Stopped by user. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        DEBUG_OUT(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);
        DEBUG_OUT(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);
        DEBUG_OUT(tr("=== LM AUTOMATIC FITTING - UNKNOWN END ==="));
    }

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

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

void nmCalculationAutoFitLM::extractUserInitialValues()
{
    // 按勾选顺序提取快照初值，并沿用窗口确认的上下限。
    if(m_sensitivitySettings) {
        m_initialValues.clear();
        for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
            const int index = m_enabledParamIndices[i];
            m_initialValues.append(qBound(m_parameterLower[index],
                m_sensitivitySettings->initialValues[index], m_parameterUpper[index]));
        }
        return;
    }

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

    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: // 裂缝导流能力
                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 5: // 裂缝半长
                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, false);
    *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()
{
    // 本次拟合开始时按目标井的实际井型确定流程，普通直井跳过形状调整。
    nmDataAnalyzeManager* dataManager = m_sensitivitySettings ? nullptr : nmDataAnalyzeManager::getCurrentInstance();
    nmDataWellBase* targetWell = dataManager ? dataManager->findWellByName(m_targetWellName) : nullptr;
    const bool skipShapeStage = m_sensitivitySettings ? m_sensitivityVerticalWell
        : targetWell && targetWell->getWellType() == NM_WELL_MODEL::Vertical_Well;
    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;
    }

    // 高度、形状、总误差共用迭代上限；两个 LM 目标分别设置求解评价额度。
    const int phaseEvaluationBudget = qMax(dimensions + 2, qMax(20, m_maxIterations * 3));
    int maximumEvaluations = 0;
    int phaseIterations = 0;
    int shapeIterations = 0;
    int totalIterations = 0;
    int completedIterations = 0;
    // 步长使用归一化内部坐标：0.04 表示内部范围的 4%，并非物理值变化 4%。
    // 信赖半径限制所有参数一次联合移动的二范数。
    const double sensitivityStep = 0.04;
    const double minimumCoordinateStep = 1.0e-5;
    const double minimumTrustRadius = 2.0e-3;
    const double maximumTrustRadius = 0.30;
    // 误差下降至少达到绝对 0.01 且相对当前有效基准 10% 才算有效改善。
    // 更小的下降仍保留为最佳解，但不能反复清除停滞状态、延长拟合时间。
    const double effectiveRelativeImprovement = 1.0e-1;
    const double effectiveAbsoluteImprovement = 1.0e-2;
    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;
    int consecutivePoorPredictions = 0;
    bool rebuildRequested = true;
    // 跟踪文件保留稳定英文原因，界面输出时再翻译，避免受本地编码影响。
    const char* rebuildReason = QT_TR_NOOP("initial sensitivity model");
    bool modelRebuiltAtMinimumRadius = false;
    bool stagnationConfirmationRequested = false;
    StopReasonLM stopReason = LM_MAX_ITERATIONS;

    // jacobian 的行对应固定采样的残差，列对应用户勾选的参数。
    // LM 按当前目标取数值或形状残差行，两段均使用所有勾选参数。
    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。无论候选是否接受，
        // 下一次计算前都恢复到唯一的已接受工作点，防止失败试算污染后续求解。
        if(!m_sensitivitySettings) 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 {
        if(m_sensitivitySettings) m_shouldStop = m_sensitivityCancelled();
        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 progressUpdated(0, current.fitness);
    DEBUG_OUT(tr("=== Starting LM Main Loop ==="));
    DEBUG_OUT(
        tr("LM starting point error: %1; evaluation budget per LM phase: %2")
        .arg(current.fitness, 0, 'e', 4)
        .arg(phaseEvaluationBudget));

    // 高度预调整只改变用户勾选的渗透率。ln(k) 的初始变化由有符号高度差
    // 给出，真实求解后用割线估计修正；拒绝时缩步，不让其他参数补偿高度。
    const int permeabilityColumn = m_enabledParamIndices.indexOf(0);
    double heightBaseline = current.breakdown.verticalLoss;
    int heightStagnation = 0;
    double heightScale = 1.0;
    double heightSlope = -1.0;
    int heightTrials = 0;
    if(permeabilityColumn >= 0) {
        emit fittingStageChanged(1);
        emit logMessageGenerated(tr("[Stage 1: Height alignment]"));
        // 高度调整计入总迭代数；对齐、停滞或无可行步长时结束。
        while(completedIterations < m_maxIterations && processPauseAndStop()) {
            const double bias = current.breakdown.verticalCommonBias;
            if(!current.breakdown.verticalReliable || qAbs(bias) <= 0.1 || heightStagnation >= 2) break;
            const double k = current.parameters[permeabilityColumn];
            if(k <= 0.0 || m_parameterLower[0] <= 0.0 || m_parameterUpper[0] <= m_parameterLower[0]) break;
            const double logRange = qLn(m_parameterUpper[0]) - qLn(m_parameterLower[0]);
            const double change = qBound(-0.30 * logRange, -bias / heightSlope * heightScale, 0.30 * logRange);
            const AutoFitObjectiveBreakdownLM previous = current.breakdown;
            TrustRegionEvaluation candidate;
            candidate.parameters = current.parameters;
            candidate.parameters[permeabilityColumn] = qBound(m_parameterLower[0], k * qExp(change), m_parameterUpper[0]);
            candidate.coordinates = coordinatesFromParameters(candidate.parameters);
            const double actualStep = qLn(candidate.parameters[permeabilityColumn] / k);
            if(qAbs(actualStep) < 1.0e-5) break;
            m_currentIteration = completedIterations++;
            ++heightTrials;
            emit progressUpdated(completedIterations, current.fitness);
            candidate.valid = evaluateTrustRegionPoint(candidate.parameters, &candidate.fitness,
                &candidate.breakdown, &candidate.curve, &candidate.elapsedMs);
            // 第一阶段只按上下偏移改善验收，不限制形状；仍防止两条曲线一高一低互相抵消。
            const bool accepted = candidate.valid && candidate.breakdown.verticalReliable &&
                candidate.breakdown.verticalLoss < current.breakdown.verticalLoss;
            if(candidate.valid) {
                const double slope = (candidate.breakdown.verticalCommonBias - bias) / actualStep;
                if(slope < -0.05 && isFiniteNumber(slope)) heightSlope = slope;
            }
            writeTraceRow(m_currentIteration, 0, "permeability_height", candidate.parameters, candidate.fitness,
                candidate.valid, candidate.elapsedMs, accepted ? "accepted_height" : "rejected_height",
                candidate.valid ? &candidate.breakdown : nullptr);
            if(accepted) {
                current = candidate;
                publishAcceptedPoint(current);
            } else {
                heightScale *= 0.5;
            }
            restoreEvaluationState(current);
            if(heightBaseline - current.breakdown.verticalLoss >= qMax(1.0e-4, 0.10 * heightBaseline)) {
                heightBaseline = current.breakdown.verticalLoss;
                heightStagnation = 0;
            } else ++heightStagnation;
            if(candidate.valid) {
                const QString objective = tr("Height error: %1 -> %2")
                    .arg(lmLogError(previous.verticalLoss)).arg(lmLogError(candidate.breakdown.verticalLoss));
                emit logMessageGenerated(tr("%1 | %2 | %3 | Improvement: %4%")
                    .arg(tr("Iteration %1").arg(completedIterations))
                    .arg(accepted ? tr("accepted") : tr("rejected")).arg(objective)
                    .arg(lmLogImprovement(previous.verticalLoss, candidate.breakdown.verticalLoss)));
                if(!accepted) emit logMessageGenerated(candidate.breakdown.verticalReliable
                    ? tr("Height error did not decrease; keep the last accepted result and reduce the step.")
                    : tr("Pressure and derivative height directions conflict; keep the last accepted result and reduce the step."));
            } else {
                emit logMessageGenerated(tr("%1 | Solver failed\nNo valid error was returned; keep the last accepted result and reduce the step.")
                    .arg(tr("Iteration %1").arg(completedIterations)));
            }
            emit progressUpdated(completedIterations, current.fitness);
        }
        QString heightReason = tr("no feasible permeability step or step too small");
        if(m_shouldStop) heightReason = tr("stopped by user");
        else if(completedIterations >= m_maxIterations) heightReason = getStopReasonDescription(LM_MAX_ITERATIONS);
        else if(!current.breakdown.verticalReliable) heightReason = tr("pressure and derivative height directions conflict");
        else if(current.breakdown.verticalLoss <= 0.1) heightReason = tr("curve height is approximately aligned");
        else if(heightStagnation >= 2) heightReason = tr("2 consecutive steps without effective height improvement");
        else if(current.parameters[permeabilityColumn] <= m_parameterLower[0] ||
                current.parameters[permeabilityColumn] >= m_parameterUpper[0]) heightReason = tr("permeability reached its bound");
        emit logMessageGenerated(tr("Height alignment ended: %1").arg(heightReason));
    }

    if(permeabilityColumn < 0) emit logMessageGenerated(tr("Height alignment skipped: permeability is not selected."));
    if(m_shouldStop) return LM_USER_STOPPED;
    if(completedIterations >= m_maxIterations) return LM_MAX_ITERATIONS;

    // 直井直接匹配总误差，其他井型仍先匹配整体形状，再匹配总误差。
    bool shapeStage = !skipShapeStage;
    auto stageError = [&](const TrustRegionEvaluation& point) -> double {
        return shapeStage ? point.breakdown.shapeLoss : point.fitness;
    };
    auto acceptable = [&](const TrustRegionEvaluation& point, const TrustRegionEvaluation& base) -> bool {
        return point.valid && stageError(point) < stageError(base);
    };
    auto acceptPoint = [&](const TrustRegionEvaluation& point) {
        // 同步发布当前目标下已接受的参数、曲线和诊断。
        current = point;
        publishAcceptedPoint(current);
        restoreEvaluationState(current);
    };
    double effectiveImprovementBaseline = stageError(current);
    DEBUG_OUT(tr("LM sampling: %1 points (%2 intervals per log-time decade).")
        .arg(current.breakdown.residualVector.size() / 2).arg(kAutoFitIntervalsPerDecade));

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

    // 两个 LM 目标共用有效改善门槛和重建后的停滞确认规则。
    auto recordIneffectiveStep = [&]() -> bool {
        ++consecutiveIneffectiveSteps;
        if(consecutiveIneffectiveSteps < maximumIneffectiveSteps) {
            return false;
        }

        if(stagnationConfirmationRequested) {
            return true;
        }

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

    DEBUG_OUT(
        tr("LM 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));

    // 在同一个真实工作点逐参数做单边差分。失败先在原方向缩步，再反向尝试；
    // 正常情况下每列仍只需一次真实求解，重试也计入总评价预算。
    auto rebuildSensitivity = [&]() -> bool {
        // 灵敏度试算计入当前目标的评价额度，形状和总误差使用相同差分规则。
        const int evaluationLimit = maximumEvaluations;
        const TrustRegionEvaluation base = current;
        const QVector<double> baseResidual = trustRegionFullResidual(base.breakdown);
        const int residualCount = baseResidual.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;
        // 两个目标共用局部差分尺度，信赖半径缩小时同步减小扰动。
        const double finiteDifferenceStep = qMin(
            sensitivityStep, qMax(5.0e-3, trustRadius * 0.5));

        for(int column = 0;
                column < dimensions &&
                m_totalEvaluations < evaluationLimit &&
                processPauseAndStop();
                ++column) {
            const int parameterIndex = m_enabledParamIndices[column];

            const double lower = m_parameterLower[parameterIndex];
            const double upper = m_parameterUpper[parameterIndex];
            if(upper <= lower) continue;
            // 表皮包含零和负值，按当前物理值的局部尺度确定扰动。
            double localStep = finiteDifferenceStep;
            if(parameterIndex == 1) {
                localStep = qMin(localStep, 0.02 * qMax(0.1, qAbs(base.parameters[column])) / (upper - lower));
            }
            double positiveRoom = 1.0 - base.coordinates[column];
            double negativeRoom = base.coordinates[column];
            double preferredSign = positiveRoom >= negativeRoom ? 1.0 : -1.0;
            bool columnBuilt = false;
            double probeStep = localStep;

            for(int directionAttempt = 0;
                    directionAttempt < 2 &&
                    !columnBuilt &&
                    m_totalEvaluations < evaluationLimit && processPauseAndStop();
                    ++directionAttempt) {
                double direction = directionAttempt == 0
                                   ? preferredSign : -preferredSign;
                double availableRoom = direction > 0.0
                                       ? positiveRoom : negativeRoom;
                double deltaMagnitude = qMin(probeStep, availableRoom);
                for(int shrinkAttempt = 0; shrinkAttempt < 3 && !columnBuilt &&
                        m_totalEvaluations < evaluationLimit && processPauseAndStop(); ++shrinkAttempt) {
                    if(deltaMagnitude < minimumCoordinateStep) {
                        break;
                    }
                    const bool canShrink = shrinkAttempt < 2 && deltaMagnitude * 0.5 >= minimumCoordinateStep;

                    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"
                                        : (canShrink ? "sensitivity_retry_smaller" : (directionAttempt == 0
                                           ? "sensitivity_retry_opposite"
                                           : "sensitivity_invalid"));
                    writeTraceRow(m_currentIteration,
                                  column,
                                  "trust_region_sensitivity",
                                  probe.parameters,
                                  probe.fitness,
                                  probe.valid,
                                  probe.elapsedMs,
                                  (shapeStage ? "shape_" : "total_") + decision,
                                  probe.valid ? &probe.breakdown : nullptr);

                    if(!probe.valid) {
                        restoreEvaluationState(base);
                        if(!canShrink) break;
                        deltaMagnitude *= 0.5;
                        probeStep = deltaMagnitude;
                        continue;
                    }

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

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

                    jacobianColumnValid[column] = true;
                    columnBuilt = true;

                    // 缓存当前目标下更优的探针，全部灵敏度列建完后再接受。
                    if(acceptable(probe, base) &&
                            (!bestProbe.valid || stageError(probe) < stageError(bestProbe))) {
                        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;
        }

        // 所有列基于同一工作点建立，再接受更优探针并用割线更新 Jacobian。
        if(bestProbe.valid && bestProbeColumn >= 0) {
            QVector<double> acceptedStep(dimensions, 0.0);
            acceptedStep[bestProbeColumn] = bestProbeDelta;
            updateTrustRegionJacobian(
                &jacobian,
                trustRegionFullResidual(base.breakdown),
                trustRegionFullResidual(bestProbe.breakdown),
                acceptedStep);
            acceptPoint(bestProbe);
            writeTraceRow(m_currentIteration,
                          bestProbeColumn,
                          "trust_region_sensitivity_accept",
                          current.parameters,
                          current.fitness,
                          true,
                          0,
                          shapeStage ? "shape_accepted_cached_probe" : "total_accepted_cached_probe",
                          &current.breakdown);
            emitTrialLog(tr("Sensitivity trial"), base.breakdown, current.breakdown, true, true, shapeStage);
        } else {
            restoreEvaluationState(current);
        }

        // 有更优探针时以接受后的误差为新基准，否则仍为重建起点误差。
        // 重新计数但保留停滞确认状态，确认阶段连续三次改善不足仍会结束。
        effectiveImprovementBaseline = stageError(current);
        consecutiveIneffectiveSteps = 0;
        acceptedSinceRebuild = 0;
        consecutivePoorPredictions = 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;
    };

    for(int phase = skipShapeStage ? 1 : 0;
            phase < 2 && !m_shouldStop && completedIterations < m_maxIterations; ++phase) {
        shapeStage = phase == 0;
        phaseIterations = 0;
        maximumEvaluations = m_totalEvaluations + phaseEvaluationBudget;
        stopReason = LM_MAX_ITERATIONS;
        effectiveImprovementBaseline = stageError(current);
        trustRadius = 0.12;
        damping = 0.01;
        consecutiveRejectedSteps = 0;
        consecutiveSolverFailures = 0;
        acceptedSinceRebuild = 0;
        consecutiveIneffectiveSteps = 0;
        consecutivePoorPredictions = 0;
        modelRebuiltAtMinimumRadius = false;
        stagnationConfirmationRequested = false;
        jacobian.clear();
        rebuildRequested = true;
        rebuildReason = shapeStage ? QT_TR_NOOP("full sensitivity at shape LM entry")
                                   : QT_TR_NOOP("full sensitivity at total-stage entry");
        emit fittingStageChanged(shapeStage ? 2 : 3);
        emit progressUpdated(completedIterations, current.fitness);
        emit logMessageGenerated(shapeStage
            ? tr("[Stage 2: Joint LM - Shape]")
            : tr("[Stage 2: Joint LM - Total error]"));
        DEBUG_OUT(tr("LM phase budget: %1 remaining shared iterations, %2 evaluations.")
            .arg(m_maxIterations - completedIterations).arg(phaseEvaluationBudget));
        writeTraceRow(completedIterations, -1, "stage_switch", current.parameters, current.fitness,
            true, 0, shapeStage ? "height_to_shape_lm" : (skipShapeStage ? "height_to_total" : "shape_to_total"),
            &current.breakdown);

        while(!m_shouldStop && completedIterations < m_maxIterations && m_totalEvaluations < maximumEvaluations) {
            m_currentIteration = completedIterations;
            if(!processPauseAndStop()) break;
            if(stageError(current) < m_targetError) {
                stopReason = LM_TARGET_ACHIEVED;
                break;
            }
            if(rebuildRequested) {
                emit logMessageGenerated(tr("Calculating sensitivity..."));
                writeTraceRow(m_currentIteration, -1, "sensitivity_rebuild", current.parameters,
                              current.fitness, true, 0, rebuildReason, &current.breakdown);
                if(!rebuildSensitivity()) {
                    stopReason = m_shouldStop
                                 ? LM_USER_STOPPED
                                 : LM_LOCAL_OPTIMUM;
                    break;
                }
                if(stageError(current) < m_targetError) {
                    stopReason = LM_TARGET_ACHIEVED;
                    break;
                }
                if(m_totalEvaluations >= maximumEvaluations) {
                    stopReason = LM_MAX_ITERATIONS;
                    break;
                }
                // 重建已设置新基准；后续 LM 候选从该基准累计改善并确认停滞。
            }

            // 全程累计迭代；无可行方向、拒绝和求解失败也占用一次额度。
            ++phaseIterations;
            ++completedIterations;
            emit progressUpdated(completedIterations, current.fitness);
            const QVector<double> objectiveResidual = shapeStage
                ? current.breakdown.shapeResiduals : current.breakdown.residualVector;
            const int rowOffset = shapeStage ? current.breakdown.residualVector.size() : 0;
            const QVector<QVector<double> > objectiveJacobian = jacobian.mid(rowOffset, objectiveResidual.size());
            const TrustRegionNormalEquations global = buildTrustRegionNormalEquations(
                objectiveJacobian, objectiveResidual, jacobianColumnValid);
            QVector<int> selectedColumns;
            QVector<double> coordinateStep;
            double predictedReduction = 0.0;

            // 两段都联合调整全部有效自由参数，只切换残差及其对应的 Jacobian 行。
            for(int column = 0; column < dimensions; ++column) {
                if(jacobianColumnValid[column] && global.matrix[column][column] > 0.0)
                    selectedColumns.append(column);
            }
            if(selectedColumns.isEmpty() || !buildTrustRegionStep(
                    global, selectedColumns, current.coordinates, damping,
                    trustRadius, minimumCoordinateStep, &coordinateStep, &predictedReduction)) {
                selectedColumns.clear();
            }

            // 没有可行下降步时收缩半径，按当前停滞状态重建灵敏度或结束。
            if(selectedColumns.isEmpty()) {
                emit logMessageGenerated(tr("Iteration %1 | No feasible descent step; keep the last accepted result.").arg(completedIterations));
                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;
                rebuildReason = QT_TR_NOOP("no feasible descent step");
                if(recordIneffectiveStep()) {
                    stopReason = LM_LOCAL_OPTIMUM;
                    break;
                }
                continue;
            }

            double stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
            QVector<double> candidateCoordinates = current.coordinates;
            for(int column = 0; column < dimensions; ++column) {
                candidateCoordinates[column] += coordinateStep[column];
            }
            // 参数索引保留在跟踪文件中，便于检查每次联合调整的参数。
            QStringList selectedParameterIndices;
            for(int i = 0; i < selectedColumns.size(); ++i) {
                selectedParameterIndices << QString::number(m_enabledParamIndices[selectedColumns[i]]);
            }
            QString selectionName = QString(shapeStage ? "shape_global_params_" : "total_global_params_") +
                selectedParameterIndices.join("_");

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

            if(m_shouldStop) {
                restoreEvaluationState(current);
                stopReason = LM_USER_STOPPED;
                break;
            }

            if(!candidate.valid) {
                emitTrialLog(tr("Iteration %1").arg(completedIterations), current.breakdown, candidate.breakdown,
                    false, false, shapeStage);
                // 求解失败的候选不能改变 current。先完整恢复上一个已接受参数和
                // 对应误差快照，再缩小信赖域；连续失败达到上限才终止整个拟合。
                ++consecutiveSolverFailures;
                ++consecutiveRejectedSteps;
                consecutivePoorPredictions = 0;
                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_" + selectionName,
                              nullptr);
                restoreEvaluationState(current);
                if(consecutiveSolverFailures >= 2) {
                    rebuildRequested = true;
                    rebuildReason = QT_TR_NOOP("2 consecutive solver failures");
                }
                if(consecutiveSolverFailures >= m_maxConsecutiveFailures) {
                    stopReason = LM_CONSECUTIVE_FAILURES;
                    break;
                }
                if(recordIneffectiveStep()) {
                    stopReason = LM_LOCAL_OPTIMUM;
                    break;
                }
                continue;
            }

            if(m_shouldStop) {
                restoreEvaluationState(current);
                stopReason = LM_USER_STOPPED;
                break;
            }
            consecutiveSolverFailures = 0;
            // 有效候选即使最终被拒绝，也提供了一条真实割线，可用于修正下一轮
            // 局部模型；是否成为新工作点由当前阶段的目标和约束共同决定。
            const AutoFitObjectiveBreakdownLM oldBreakdown = current.breakdown;
            updateTrustRegionJacobian(
                &jacobian,
                trustRegionFullResidual(oldBreakdown),
                trustRegionFullResidual(candidate.breakdown),
                coordinateStep);
            // reductionRatio 衡量局部线性模型的可信度：接近 1 表示预测准确；
            // 值较小表示虽然可能下降，但模型低估了非线性，需要收紧下一步。
            const double objectiveEnergy = 0.5 * trustRegionSquaredNorm(objectiveResidual);
            const QVector<double> candidateResidual = shapeStage
                ? candidate.breakdown.shapeResiduals : candidate.breakdown.residualVector;
            const double candidateEnergy = 0.5 * trustRegionSquaredNorm(candidateResidual);
            double actualReduction = objectiveEnergy - candidateEnergy;
            double reductionRatio = predictedReduction > 1.0e-14 ? actualReduction / predictedReduction : 0.0;
            // 使用固定采样和同一阶段目标比较预测与实际改善。接近收敛时的微小
            // 预测量交给停滞逻辑处理，避免比例数值波动反复触发昂贵的全参数重建。
            const double predictionFloor = qMax(1.0e-14, objectiveEnergy * 1.0e-8);
            const bool poorPrediction = predictedReduction > predictionFloor &&
                (!isFiniteNumber(reductionRatio) || reductionRatio < 0.25);
            consecutivePoorPredictions = poorPrediction ? consecutivePoorPredictions + 1 : 0;
            if(consecutivePoorPredictions >= 2) {
                rebuildRequested = true;
                rebuildReason = QT_TR_NOOP("2 consecutive steps with actual improvement below 25% of prediction");
            }
            bool accepted = acceptable(candidate, current);

            if(accepted) {
                // 当前阶段接受候选后同步发布参数、曲线和诊断。
                // 模型预测可靠时减小阻尼并可扩大半径，预测较差时保守收缩。
                acceptPoint(candidate);
                ++acceptedSinceRebuild;
                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 >= 10 && !rebuildRequested) {
                    rebuildRequested = true;
                    rebuildReason = QT_TR_NOOP("10 accepted steps since last rebuild");
                }
                modelRebuiltAtMinimumRadius = false;
            } else {
                // 拒绝后恢复已接受参数、增大阻尼并缩小半径；按预测质量判断是否重建。
                ++consecutiveRejectedSteps;
                damping = qMin(1.0e8, damping * 4.0);
                trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
                restoreEvaluationState(current);
            }

            writeTraceRow(m_currentIteration,
                          -1,
                          "trust_region_candidate",
                          candidate.parameters,
                          candidate.fitness,
                          true,
                          candidate.elapsedMs,
                          accepted
                          ? "accepted_" + selectionName
                          : "rejected_" + selectionName,
                          &candidate.breakdown);

            emitTrialLog(tr("Iteration %1").arg(completedIterations), oldBreakdown, candidate.breakdown,
                true, accepted, shapeStage);
            emit progressUpdated(completedIterations, m_globalBestFitness);

            const bool effectiveImprovement = registerEffectiveImprovement(stageError(current));
            if(!effectiveImprovement && recordIneffectiveStep()) stopReason = LM_LOCAL_OPTIMUM;

            if(stopReason == LM_LOCAL_OPTIMUM) break;

            if(stageError(current) < m_targetError) {
                stopReason = LM_TARGET_ACHIEVED;
                break;
            }
            if(trustRadius <= minimumTrustRadius * 1.01 &&
                    consecutiveRejectedSteps >= 2 && consecutivePoorPredictions >= 2) {
                if(modelRebuiltAtMinimumRadius) {
                    stopReason = LM_LOCAL_OPTIMUM;
                    break;
                }
                rebuildRequested = true;
                rebuildReason = QT_TR_NOOP("inaccurate model at minimum trust radius");
            }
        }
        if(shapeStage) shapeIterations = phaseIterations;
        else totalIterations = phaseIterations;
        if(m_shouldStop || stopReason == LM_CONSECUTIVE_FAILURES || stopReason == LM_OPTIMIZATION_FAILED)
            break;
        // 形状部分结束后仅在总迭代额度有剩余时进入总误差部分。
        if(shapeStage) {
            emit logMessageGenerated(tr("Shape phase ended: %1").arg(getStopReasonDescription(stopReason)));
            writeTraceRow(m_currentIteration, -1, "lm_phase_end", current.parameters, current.fitness,
                true, 0, QString("shape_stop_reason_%1").arg(static_cast<int>(stopReason)), &current.breakdown);
        } else {
            emit logMessageGenerated(tr("Total-error phase ended: %1").arg(getStopReasonDescription(stopReason)));
        }
    }

    // 灵敏度建模或初始误差已达标时可能没有候选迭代，不多计一次。
    m_currentIteration = completedIterations - 1;
    restoreEvaluationState(current);
    DEBUG_OUT(tr("Adaptive fitting counts: %1 height iterations, %2 shape LM iterations, %3 total LM iterations, %4 total evaluations.")
        .arg(heightTrials).arg(shapeIterations).arg(totalIterations).arg(m_totalEvaluations));

    if(m_shouldStop) {
        return LM_USER_STOPPED;
    }
    if(!shapeStage && current.fitness < m_targetError) {
        return LM_TARGET_ACHIEVED;
    }
    if(completedIterations >= m_maxIterations && stopReason != LM_CONSECUTIVE_FAILURES && stopReason != LM_OPTIMIZATION_FAILED) {
        return LM_MAX_ITERATIONS;
    }
    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, bool retrySolver)
{
    // 校验并应用候选参数，调用求解器后比较目标井曲线；失败返回惩罚值 1e10。
    const QString funcName = QString("evaluateError[%1]").arg(m_currentIteration);
    static int callCount = 0;
    callCount++;
    m_lastEvaluatedLogLogData.clear();
    m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();

    // 候选求解在冻结输入副本上执行，误差口径仍使用本算法现有实现。
    if(m_sensitivitySettings) {
        m_shouldStop = m_sensitivityCancelled();
        if(m_shouldStop || !validateParameters(parameters)) return 1.0e10;
        const QVector<QVector<double> > curve = m_sensitivityEvaluate(parameters);
        m_shouldStop = m_sensitivityCancelled();
        if(m_shouldStop || !validateLogLogData(curve)) return 1.0e10;
        const double fitness = calculateLogLogCurveError(m_targetLogLogData, curve);
        if(fitness < 1.0e9) m_lastEvaluatedLogLogData = curve;
        return fitness;
    }

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

        // 拦截维度不匹配、非有限数或超出用户范围的参数。
        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));
                        }
                    }
                }
            }

            return 1e10;
        }

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

        // 参数写回和求解任务依赖当前数据管理器。
        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));

        // 按勾选顺序将候选值写入储层和目标井。
        try {
            DEBUG_OUT(QString("%1: Call #%2 - Applying parameters...").arg(funcName).arg(callCount));
            applyParametersToDataManager(parameters);
            DEBUG_OUT(QString("%1: Call #%2 - Parameters applied successfully").arg(funcName).arg(callCount));
        } catch(const std::exception& e) {
            DEBUG_OUT(QString("%1: Call #%2 - FAILED to apply parameters: %3")
                      .arg(funcName).arg(callCount).arg(e.what()));
            return 1e10;
        }

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

        // 初始评价失败可重试；LM 候选失败则交给信赖域缩步，不重复相同参数。
        QVector<QVector<double>> solverResult;
        const int maxRetries = retrySolver ? 2 : 0;
        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 = runSolverDll();

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

        // 使用求解任务复制的目标井曲线，避免读取其他井或上一候选的共享结果。
        QVector<QVector<double>> resultLogLogData = m_lastEvaluatedLogLogData;

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

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

        // 比较目标井的实测与模拟双对数曲线。
        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)
{
    // 按 m_enabledParamIndices 对应的顺序，将勾选参数写入储层和目标井。
    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 是勾选参数向量中的索引；i 是完整 6 个参数体系中的索引。
    // 只有 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;

            }

            paramIndex++;
        }
    }

    dataManager->updateReservoirData(reservoirData);
}

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

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

    if(!pWell) return;

    // 先在参数副本中组装本次候选值，全部解析完成后再写回现有井对象。
    // 这样不会触发整井赋值，也不会删除并重建井内已有的射孔对象。
    nmDataPerforation* pPerforation = pWell->getPerforation(0);
    nmDataAttribute skinAttr;
    bool updateSkin = false;
    if(pPerforation) {
        skinAttr = pPerforation->getSkin();
    }

    nmDataAttribute wellboreAttr = pWell->getWellboreStorage();
    bool updateWellboreStorage = false;

    nmDataVerticalFracturedWell* pVerticalFracturedWell =
        dynamic_cast<nmDataVerticalFracturedWell*>(pWell);
    nmDataHorizontalFracturedWell* pHorizontalFracturedWell =
        dynamic_cast<nmDataHorizontalFracturedWell*>(pWell);

    nmDataAttribute dfcAttr;
    nmDataAttribute fractureHalfLengthAttr;
    if(pVerticalFracturedWell) {
        dfcAttr = pVerticalFracturedWell->getDfc();
        fractureHalfLengthAttr = pVerticalFracturedWell->getFractureHalfLength();
    } else if(pHorizontalFracturedWell) {
        dfcAttr = pHorizontalFracturedWell->getDfc();
        fractureHalfLengthAttr = pHorizontalFracturedWell->getFractureHalfLength();
    }
    bool updateDfc = false;
    bool updateFractureHalfLength = false;

    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: { // 表皮系数
                        if(pPerforation) {
                            skinAttr.setValue(value);
                            updateSkin = true;
                        }
                    }
                    break;

                case 2: { // 井筒储集系数
                        wellboreAttr.setValue(value);
                        updateWellboreStorage = true;
                    }
                    break;

                case 4: { // 裂缝导流能力
                        if(pVerticalFracturedWell || pHorizontalFracturedWell) {
                            dfcAttr.setValue(value);
                            updateDfc = true;
                        }
                    }
                    break;

                case 5: { // 裂缝半长
                        if(pVerticalFracturedWell || pHorizontalFracturedWell) {
                            fractureHalfLengthAttr.setValue(value);
                            updateFractureHalfLength = true;
                        }
                    }
                    break;
            }

            paramIndex++;
        }
    }

    if(updateSkin) {
        pPerforation->getSkin().setValue(skinAttr.getValue());
    }
    if(updateWellboreStorage) {
        pWell->getWellboreStorage().setValue(wellboreAttr.getValue());
    }
    if(pVerticalFracturedWell) {
        if(updateDfc) {
            pVerticalFracturedWell->getDfc().setValue(dfcAttr.getValue());
        }
        if(updateFractureHalfLength) {
            pVerticalFracturedWell->getFractureHalfLength().setValue(
                fractureHalfLengthAttr.getValue());
        }
    } else if(pHorizontalFracturedWell) {
        if(updateDfc) {
            pHorizontalFracturedWell->getDfc().setValue(dfcAttr.getValue());
        }
        if(updateFractureHalfLength) {
            pHorizontalFracturedWell->getFractureHalfLength().setValue(
                fractureHalfLengthAttr.getValue());
        }
    }
}

// ==================== 算法辅助方法 ====================

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

    DEBUG_OUT(tr("=== Final Result Validation (Elite Protection) ==="));

    double finalFitness = m_globalBestFitness;
    double initialFitness = m_userInitialFitness;

    DEBUG_OUT(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));
        DEBUG_OUT(tr("Improvement: %1 (%2%)")
                                 .arg(improvement, 0, 'e', 4)
                                 .arg(relativeImprovement * 100, 0, 'f', 2));
    }

    if(!finalValid || finalFitness > initialFitness) {
        DEBUG_OUT(
            tr("Elite protection triggered: final result is invalid or worse than initial"));
        emit logMessageGenerated(tr("Final total error did not improve; restoring the initial 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);

        DEBUG_OUT(tr("Initial solution restored successfully"));
    } else {
        DEBUG_OUT(
            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;
        }
    }

    // 单调性仅作诊断；误差计算会重新排序，不在此拒绝结果。
    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)
{
    // 首次有效评价后固定公共时间范围，保证各候选在相同时间点比较。
    // 总误差由完整压力和导数残差计算，形状阶段单独使用斜率残差。
    // 整个计算过程均位于 log(time)-log(value) 坐标。
    const double invalidLoss = 1.0e10;
    const double valueFloor = 1.0e-12;
    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 ||
                        data[2][i] <= 0.0) {
                    continue;
                }

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

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

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

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

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

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

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

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

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

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

        // 仅首次有效评价使用交集建立基准；失败试算不能冻结区间，后续候选
        // 必须覆盖完整基准，不允许靠丢失首尾点缩小误差。
        const bool comparisonRangeFixed = m_comparisonTimeMin > 0.0;
        double overlapMinX = comparisonRangeFixed
                ? m_comparisonTimeMin : qMax(targetMinX, resultMinX);
        double overlapMaxX = comparisonRangeFixed
                ? m_comparisonTimeMax : qMin(targetMaxX, resultMaxX);
        if(overlapMinX >= overlapMaxX) {
            return invalidLoss;
        }
        if(targetMinX > overlapMinX || targetMaxX < overlapMaxX ||
                resultMinX > overlapMinX || resultMaxX < overlapMaxX) {
            DEBUG_OUT("Candidate does not cover the fixed LM comparison time range");
            return invalidLoss;
        }

        // 所有指标共用每个对数时间数量级 20 个间隔，向上取整并包含两端，至少三个点。
        // 比较区间在本轮固定，因此点数也固定；微小容差避免整数量级因舍入多取一点。
        const double timeDecades = (qLn(overlapMaxX) - qLn(overlapMinX)) / qLn(10.0);
        const int numPoints = qMax(3, qCeil(timeDecades * kAutoFitIntervalsPerDecade - 1.0e-10) + 1);
        QVector<double> commonX(numPoints);
        QVector<double> commonLogX(numPoints);
        QVector<double> targetLogPressure(numPoints);
        QVector<double> targetLogDerivative(numPoints);
        const double comparisonLogMinX = qLn(overlapMinX);
        const double comparisonLogMaxX = qLn(overlapMaxX);

        // 在公共时间范围内生成等距 log-time 网格，压力和导数各使用 numPoints 个残差。
        for(int i = 0; i < numPoints; ++i) {
            double logX = comparisonLogMinX +
                          static_cast<double>(i) *
                          (comparisonLogMaxX - comparisonLogMinX) /
                          (numPoints - 1);
            commonLogX[i] = logX;
            // 首尾直接使用原始端点，避免 exp(log(t)) 的舍入误差越过严格插值边界。
            if(i == 0) {
                commonX[i] = overlapMinX;
            } else if(i == numPoints - 1) {
                commonX[i] = overlapMaxX;
            } else {
                commonX[i] = qExp(logX);
            }

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

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

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

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

        auto populateStageMetrics = [&](AutoFitObjectiveBreakdownLM* objective) -> bool {
            // 复用同一网格上已算出的原始残差，避免为高度和形状再次插值。
            // 斜率跨度保持约为全对数时域的 10%，压力和导数仍等权。
            const int count = numPoints;
            const int lag = autoFitShapeLag(count);
            const double span = qLn(overlapMaxX) - qLn(overlapMinX);
            const QVector<double>* residuals[] = { &pressureResidual, &derivativeResidual };
            double biases[2] = {0.0, 0.0};
            for(int component = 0; component < 2; ++component) {
                for(int i = 0; i < count; ++i) {
                    const double residual = (*residuals[component])[i];
                    // 使用对数时间梯形权重，避免密集段主导高度。
                    biases[component] += residual * ((i == 0 || i == count - 1) ? 0.5 : 1.0)
                                         / (count - 1);
                }
            }
            objective->pressureVerticalBias = biases[0];
            objective->derivativeVerticalBias = biases[1];
            objective->verticalCommonBias = 0.5 * (biases[0] + biases[1]);
            objective->verticalLoss = qAbs(objective->verticalCommonBias);
            objective->verticalReliable = !(biases[0] * biases[1] < 0.0 &&
                qMin(qAbs(biases[0]), qAbs(biases[1])) > 0.01);
            objective->shapeResiduals.clear();
            const double scale = qSqrt(0.5 / (count - lag));
            for(int component = 0; component < 2; ++component) {
                for(int i = 0; i < count - lag; ++i) {
                    objective->shapeResiduals.append(scale *
                        ((*residuals[component])[i + lag] - (*residuals[component])[i]) /
                        (span * lag / (count - 1)));
                }
            }
            objective->shapeLoss = qSqrt(trustRegionSquaredNorm(objective->shapeResiduals));
            return isFiniteNumber(objective->shapeLoss);
        };

        AutoFitObjectiveBreakdownLM breakdown;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        // 平坦的配准误差不一定表示歧义；只有近优位移显著改变公共偏差时，
        // 才在旧配准诊断中标记高度和水平位移无法区分。
        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;
        }

        // 此处保留旧配准诊断的高度有效性检查。真正用于高度调整的指标，
        // 在 populateStageMetrics() 中按未平移曲线的公共偏差重新计算。
        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;
            }
        }

        // 此处检查配准后的剩余误差；LM 使用的斜率形状误差在后面覆盖写入。
        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;
        if(breakdown.valid && breakdown.total < 1.0e9) {
            if(!comparisonRangeFixed) {
                // 整个目标及诊断均有效后才提交基准，中点回退可重新建立区间。
                m_comparisonTimeMin = overlapMinX;
                m_comparisonTimeMax = overlapMaxX;
                writeTraceMetaFile();
            }
        }
        if(!populateStageMetrics(&breakdown)) return invalidLoss;
        m_lastObjectiveBreakdown = breakdown;

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

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

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

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

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

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

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

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

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

        // 异步执行 DLL 任务。循环等待期间持续 processEvents，保证界面不会完全卡死。
        connect(dllTask, SIGNAL(solverInvoked()), this, SLOT(recordSolverCall()), Qt::DirectConnection);
        dllTask->start();

        // 等待完成。最大等待 1 小时；用户停止使用任务已有的协作取消接口。
        const int maxWait = 3600000;
        const int checkInterval = 50;
        QTime waitTimer;
        waitTimer.start();

        while(waitTimer.elapsed() < maxWait) {
            // 求解完成时立即返回；等待期间每 50 ms 处理停止请求和界面事件。
            if(dllTask->wait(checkInterval)) {
                DEBUG_OUT("DLL solver task completed");
                break;
            }

            // 等待期间也派发鼠标、键盘事件，让停止按钮能及时登记请求。
            QApplication::processEvents(QEventLoop::AllEvents, checkInterval);

            if(m_shouldStop) {
                // DLL 没有中断接口，返回后丢弃结果，避免强制终止破坏其内部锁。
                dllTask->requestCancel();
            }
        }

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

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

        // 线程结束后检查真实执行结果，防止失败时复用上一候选的旧曲线。
        dllTask->wait();
        if(m_shouldStop) {
            DEBUG_OUT("DLL solver evaluation cancelled by user");
            delete dllTask;
            dllTask = nullptr;
            --m_evaluationInProgress;
            return result;
        }
        if(!dllTask->wasSuccessful()) {
            DEBUG_OUT("DLL solver task reported failure");
            delete dllTask;
            dllTask = nullptr;
            --m_evaluationInProgress;
            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()));

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

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

    // 释放求解任务，退出评价前处理待派发的界面事件。
    if(dllTask) {
        if(dllTask->isRunning()) {
            dllTask->terminate();
            dllTask->wait(3000);
        }

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

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

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

    return result;
}

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

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

    ++m_evaluationInProgress;
    nmCalculationDllPebiSolverTask dllTask(m_tempDirectory, dataManager);
    connect(&dllTask, SIGNAL(solverInvoked()), this, SLOT(recordSolverCall()), Qt::DirectConnection);
    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(dllTask.isRunning()) {
        DEBUG_OUT("Final full-field solver timeout, terminating task");
        dllTask.terminate();
        dllTask.wait(2000);
        --m_evaluationInProgress;
        return false;
    }

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

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

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

        case LM_MAX_ITERATIONS:
            return tr("Iteration or evaluation budget 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");
    }
}