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

#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 <QDebug>
#include <QTime>
#include <QDir>
#include <QTextStream>
#include <QFileInfo>
#include <QDateTime>
#include <QtCore/qmath.h>
#include <cmath>
#include <algorithm>
#include <limits>

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

static const bool kAutoFitDiagnosticTraceEnabled = true;
static const int kAutoFitTimeWindowCount = 4;
static const double kAutoFitTimeWindowOverlapRatio = 0.20;

static double autoFitTimeWindowWeight(double coordinate, int windowIndex)
{
    // 重叠区以基础边界为中心，总宽度为窗口宽度的 20%；相邻权重互补。
    const double width = 1.0 / kAutoFitTimeWindowCount;
    const double overlap = width * kAutoFitTimeWindowOverlapRatio;
    const double halfOverlap = overlap * 0.5;
    const double pi = 3.14159265358979323846;
    double weight = 1.0;
    if(windowIndex > 0) {
        double u = qBound(0.0,
                         (coordinate - windowIndex * width + halfOverlap) / overlap,
                         1.0);
        weight *= 0.5 * (1.0 - qCos(pi * u));
    }
    if(windowIndex + 1 < kAutoFitTimeWindowCount) {
        double u = qBound(0.0,
                         (coordinate - (windowIndex + 1) * width + halfOverlap) / overlap,
                         1.0);
        weight *= 0.5 * (1.0 + qCos(pi * u));
    }
    return weight;
}

// 目标时间点构成嵌套子集。窗口中心、重叠区两端及边界使用原始点两侧作锚点。
static QVector<int> autoFitSamplingIndices(const QVector<double>& coordinates, int stride)
{
    QVector<bool> selected(coordinates.size(), false);
    for(int i = 0; i < selected.size(); i += stride) {
        selected[i] = true;
    }
    selected[0] = true;
    selected[selected.size() - 1] = true;
    for(int k = 0; k < kAutoFitTimeWindowCount; ++k) {
        const double width = 1.0 / kAutoFitTimeWindowCount;
        const double anchors[] = {(k + 0.5) * width,
            k * width - width * kAutoFitTimeWindowOverlapRatio * 0.5,
            k * width, k * width + width * kAutoFitTimeWindowOverlapRatio * 0.5};
        for(int a = 0; a < 4; ++a) {
            const int right = static_cast<int>(std::lower_bound(
                coordinates.begin(), coordinates.end(), anchors[a]) - coordinates.begin());
            if(right < selected.size()) selected[right] = true;
            if(right > 0) selected[right - 1] = true;
        }
    }
    QVector<int> indices;
    for(int i = 0; i < selected.size(); ++i) {
        if(selected[i]) indices.append(i);
    }
    return indices;
}

// 在归一化 log(time) 上使用梯形积分权重，避免原始数据密集段被重复放大。
static QVector<double> autoFitLogTimeWeights(const QVector<double>& coordinates)
{
    QVector<double> weights(coordinates.size(), 0.0);
    for(int i = 1; i < coordinates.size(); ++i) {
        const double halfWidth = 0.5 * (coordinates[i] - coordinates[i - 1]);
        weights[i - 1] += halfWidth;
        weights[i] += halfWidth;
    }
    return weights;
}

static QVector<AutoFitTimeWindowLM> calculateAutoFitTimeWindows(
        const QVector<double>& residualVector, double timeMin, double timeMax,
        const QVector<double>& coordinates = QVector<double>(),
        const QVector<double>& timeWeights = QVector<double>())
{
    // 残差前后两半分别为压力和导数，已包含各占一半及时间采样权重的归一化。
    const int pointCount = residualVector.size() / 2;
    const double logMin = qLn(timeMin);
    const double logSpan = qLn(timeMax) - logMin;
    QVector<AutoFitTimeWindowLM> windows(kAutoFitTimeWindowCount);
    for(int k = 0; k < windows.size(); ++k) {
        AutoFitTimeWindowLM& window = windows[k];
        window.timeMin = k == 0 ? timeMin
                : qExp(logMin + logSpan * k / kAutoFitTimeWindowCount);
        window.timeMax = k + 1 == windows.size() ? timeMax
                : qExp(logMin + logSpan * (k + 1) / kAutoFitTimeWindowCount);
        for(int i = 0; i < pointCount; ++i) {
            const double weight = autoFitTimeWindowWeight(
                        coordinates.isEmpty() ? static_cast<double>(i) / (pointCount - 1)
                                              : coordinates[i], k);
            window.weightSum += weight * (timeWeights.isEmpty() ? 1.0 : timeWeights[i]);
            window.energy += weight *
                    (residualVector[i] * residualVector[i] +
                     residualVector[pointCount + i] * residualVector[pointCount + i]);
        }
        window.rmsError = window.weightSum > 0.0
                ? qSqrt(window.energy * (timeWeights.isEmpty() ? pointCount : 1.0)
                        / window.weightSum) : 0.0;
    }
    return windows;
}

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 等跨数量级参数被线性尺度支配；
// skin 可为负数、Swi 的物理意义是线性比例，因此二者保持有界线性坐标。
static bool useTrustRegionLogScale(int parameterIndex, double lower, double upper)
{
    return parameterIndex != 1 && parameterIndex != 4 &&
           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)
{
    // 固定模式仍要求 160 维；分层模式按当前时间点数校验，切层后重建 J。
    const int pointCount = breakdown.sampleCoordinates.isEmpty()
                         ? 80 : breakdown.sampleCoordinates.size();
    if(!breakdown.valid || pointCount < 3 ||
            breakdown.residualVector.size() != 2 * pointCount) {
        return false;
    }

    for(int i = 0; i < breakdown.residualVector.size(); ++i) {
        if(!isFiniteNumber(breakdown.residualVector[i])) {
            return false;
        }
    }
    if(breakdown.shapeResiduals.size() != 146 || !isFiniteNumber(breakdown.shapeLoss)) return false;
    for(int i = 0; i < breakdown.shapeResiduals.size(); ++i) {
        if(!isFiniteNumber(breakdown.shapeResiduals[i])) return false;
    }
    if(breakdown.earlyValueResiduals.size() != 162 || breakdown.earlyParallelResiduals.size() != 146 ||
            !isFiniteNumber(breakdown.earlyValueLoss) || !isFiniteNumber(breakdown.earlyParallelLoss) ||
            !isFiniteNumber(breakdown.earlyParallelBias)) return false;
    for(int i = 0; i < breakdown.earlyValueResiduals.size(); ++i) {
        if(!isFiniteNumber(breakdown.earlyValueResiduals[i])) return false;
    }
    for(int i = 0; i < breakdown.earlyParallelResiduals.size(); ++i) {
        if(!isFiniteNumber(breakdown.earlyParallelResiduals[i])) return false;
    }
    return true;
}

// 数值、整体形状和前期形状残差共用一次求解，联合缓存供各子阶段复用。
static QVector<double> trustRegionFullResidual(const AutoFitObjectiveBreakdownLM& objective)
{
    QVector<double> residual = objective.residualVector;
    residual += objective.shapeResiduals;
    residual += objective.earlyValueResiduals;
    residual += objective.earlyParallelResiduals;
    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;
}

// 数值残差已含 sqrt(0.5 * 窗口权重)，相减并换底后得到 log10 间距的加权 RMS。
// 不减去平均间距，保留“两条曲线整体离得太远”的信息；共同纵移会自然抵消。
static double trustRegionEarlyGapLoss(const AutoFitObjectiveBreakdownLM& objective)
{
    const int count = objective.earlyValueResiduals.size() / 2;
    double energy = 0.0;
    for(int i = 0; i < count; ++i) {
        const double difference = objective.earlyValueResiduals[i] - objective.earlyValueResiduals[count + i];
        energy += 2.0 * difference * difference;
    }
    return qSqrt(energy) / qLn(10.0);
}

// 前期带符号的平均间距差：正值表示模拟比目标分得更开，负值表示更靠近。
// 复用现有加权数值残差，恢复第一窗口的归一化权重，不用 RMS 推断偏差符号。
static double trustRegionEarlyGapBias(const AutoFitObjectiveBreakdownLM& objective)
{
    const int count = objective.earlyValueResiduals.size() / 2;
    QVector<double> weights(count, 0.0);
    double weightSum = 0.0;
    for(int i = 0; i < count; ++i) {
        weights[i] = autoFitTimeWindowWeight(static_cast<double>(i) / (count - 1), 0) *
            ((i == 0 || i == count - 1) ? 0.5 : 1.0);
        weightSum += weights[i];
    }
    double bias = 0.0;
    for(int i = 0; i < count; ++i) {
        bias += qSqrt(2.0 * weights[i] / weightSum) *
            (objective.earlyValueResiduals[i] - objective.earlyValueResiduals[count + i]);
    }
    return bias / qLn(10.0);
}

// 按用户指定的形态规则定向：初始模拟间距大则增大井储/表皮，小则减小。
static double trustRegionEarlyGapDirection(double initialGapBias)
{
    if(!isFiniteNumber(initialGapBias) || qAbs(initialGapBias) <= 1.0e-12) return 0.0;
    return initialGapBias > 0.0 ? 1.0 : -1.0;
}

// 第一窗口的形状能量与局部 Fisher 使用相同的中心时间和重叠权重，
// 压力、导数同时参与；不重新插值或调用求解器。
static double trustRegionEarlyShapeEnergy(const QVector<double>& residual)
{
    const int pointCount = residual.size() / 2;
    double energy = 0.0;
    for(int row = 0; row < residual.size(); ++row) {
        const double coordinate = (row % pointCount + 4.0) / 80.0;
        energy += autoFitTimeWindowWeight(coordinate, 0) * residual[row] * residual[row];
    }
    return energy;
}

// 前期复用整体形状的双曲线斜率残差，只按第一窗口重新加权并归一化。
// 单独平移任一曲线不改变此指标；前期数值误差仅作诊断，不参与井储、表皮验收。
static void populateEarlyWellboreMetrics(AutoFitObjectiveBreakdownLM* objective,
        const QVector<double> targetLogs[2], const QVector<double> resultLogs[2], double logTimeSpan)
{
    const int count = 81;
    QVector<double> weights(count, 0.0);
    double weightSum = 0.0;
    for(int i = 0; i < count; ++i) {
        weights[i] = autoFitTimeWindowWeight(i / 80.0, 0) * ((i == 0 || i == count - 1) ? 0.5 : 1.0);
        weightSum += weights[i];
    }
    objective->earlyValueResiduals.fill(0.0, 2 * count);
    for(int i = 0; i < count; ++i) {
        const double weight = weights[i] / weightSum;
        for(int component = 0; component < 2; ++component) {
            const int row = component * count + i;
            const double scale = qSqrt(0.5 * weight);
            objective->earlyValueResiduals[row] = scale * (resultLogs[component][i] - targetLogs[component][i]);
        }
    }
    objective->earlyValueLoss = qSqrt(trustRegionSquaredNorm(objective->earlyValueResiduals));

    // 复用形状指标的 8 点跨度，窗口权重取斜率区间中心，减少相邻点噪声。
    const int lag = 8;
    const int slopeCount = count - lag;
    const double logTimeStep = logTimeSpan * lag / (count - 1);
    QVector<double> parallelWeights(slopeCount, 0.0);
    double parallelWeightSum = 0.0;
    for(int i = 0; i < slopeCount; ++i) {
        parallelWeights[i] = autoFitTimeWindowWeight((i + lag * 0.5) / (count - 1), 0);
        parallelWeightSum += parallelWeights[i];
    }
    objective->earlyParallelResiduals.fill(0.0, 2 * slopeCount);
    objective->earlyParallelBias = 0.0;
    for(int i = 0; i < slopeCount; ++i) {
        const double resultSlopeDifference = ((resultLogs[0][i + lag] - resultLogs[0][i]) -
            (resultLogs[1][i + lag] - resultLogs[1][i])) / logTimeStep;
        const double targetSlopeDifference = ((targetLogs[0][i + lag] - targetLogs[0][i]) -
            (targetLogs[1][i + lag] - targetLogs[1][i])) / logTimeStep;
        const double relativeSlopeError = resultSlopeDifference - targetSlopeDifference;
        const double weight = parallelWeights[i] / parallelWeightSum;
        objective->earlyParallelBias += weight * relativeSlopeError;
        // 整体残差已含 sqrt(0.5 / slopeCount)，换成第一窗口权重后仍保持双曲线等权。
        // 压力、导数分别匹配目标；相对斜率偏差只保留为诊断，不参与前期验收。
        const double scale = qSqrt(slopeCount * weight);
        objective->earlyParallelResiduals[i] = scale * objective->shapeResiduals[i];
        objective->earlyParallelResiduals[slopeCount + i] = scale * objective->shapeResiduals[slopeCount + i];
    }
    objective->earlyParallelLoss = qSqrt(trustRegionSquaredNorm(objective->earlyParallelResiduals));
}

// 计算同维向量内积；维度不一致表示局部模型无效，返回零让调用方放弃修正。
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;
}

// Fisher 仅作为当前归一化坐标下的局部信息矩阵，不用于统计置信区间。
struct TrustRegionFisher
{
    QVector<QVector<double> > matrix;
    QVector<double> gradient;

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

static QVector<TrustRegionFisher> buildTrustRegionFisher(
        const QVector<QVector<double> >& jacobian,
        const QVector<double>& residual,
        const QVector<bool>& columnValid,
        const QVector<double>& coordinates = QVector<double>())
{
    // 残差和 J 已包含压力/导数及时间采样权重，只再乘一次窗口权重。
    // 前半段是压力，后半段是导数；同一时间点的两行使用相同权重。
    const int dimensions = columnValid.size();
    const int pointCount = residual.size() / 2;
    QVector<TrustRegionFisher> information(
        kAutoFitTimeWindowCount + 1, TrustRegionFisher(dimensions));
    for(int k = 0; k < kAutoFitTimeWindowCount; ++k) {
        TrustRegionFisher& local = information[k];
        for(int row = 0; row < residual.size(); ++row) {
            const double weight = autoFitTimeWindowWeight(
                coordinates.isEmpty() ? static_cast<double>(row % pointCount) / (pointCount - 1)
                                      : coordinates[row % pointCount], k);
            for(int p = 0; p < dimensions; ++p) {
                if(!columnValid[p]) {
                    continue;
                }
                local.gradient[p] += weight * jacobian[row][p] * residual[row];
                for(int q = 0; q <= p; ++q) {
                    if(columnValid[q]) {
                        local.matrix[p][q] +=
                            weight * jacobian[row][p] * jacobian[row][q];
                    }
                }
            }
        }
        // 互补窗口之和恢复全局 J^T*J 和 J^T*r，不能对各窗口再单独平均。
        TrustRegionFisher& global = information[kAutoFitTimeWindowCount];
        for(int p = 0; p < dimensions; ++p) {
            global.gradient[p] += local.gradient[p];
            for(int q = 0; q <= p; ++q) {
                local.matrix[q][p] = local.matrix[p][q];
                global.matrix[p][q] += local.matrix[p][q];
                global.matrix[q][p] = global.matrix[p][q];
            }
        }
    }
    return information;
}

static bool buildTrustRegionFisherStep(
        const TrustRegionFisher& 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();
}

// 为每个可调参数独立求解 LM 步，按边界投影后的预测下降量选择一个参数。
static bool buildBestSingleParameterStep(const TrustRegionFisher& information,
        const QVector<int>& available, const QVector<double>& coordinates,
        double damping, double trustRadius, double minimumStep,
        QVector<int>* selected, QVector<double>* step, double* predictedReduction)
{
    selected->clear();
    step->fill(0.0, coordinates.size());
    *predictedReduction = 0.0;
    for(int i = 0; i < available.size(); ++i) {
        QVector<int> singleColumn(1, available[i]);
        QVector<double> trialStep;
        double trialPrediction = 0.0;
        if(buildTrustRegionFisherStep(information, singleColumn, coordinates,
                damping, trustRadius, minimumStep, &trialStep, &trialPrediction) &&
                trialPrediction > *predictedReduction) {
            *selected = singleColumn;
            *step = trialStep;
            *predictedReduction = trialPrediction;
        }
    }
    return !selected->isEmpty();
}

// 前期每次只调一个参数：间距锁定符号，形状 LM 提供幅度和后续验收依据。
static bool buildEarlyGapGuidedStep(const TrustRegionFisher& shapeInformation,
        const QVector<int>& selected, const QVector<double>& coordinates, double direction,
        double damping, double trustRadius, double minimumStep,
        QVector<double>* step, double* predictedReduction)
{
    if(selected.size() != 1 || direction == 0.0) return false;
    const int column = selected[0];
    TrustRegionFisher guided = shapeInformation;
    guided.gradient[column] = -direction * qAbs(shapeInformation.gradient[column]);
    if(!buildTrustRegionFisherStep(guided, selected, coordinates, damping,
            trustRadius, minimumStep, step, predictedReduction)) return false;
    // 改方向后的预测必须用原形状模型重算，不能把上坡伪装成预测下降。
    const double delta = (*step)[column];
    *predictedReduction = -shapeInformation.gradient[column] * delta -
        0.5 * shapeInformation.matrix[column][column] * delta * delta;
    return isFiniteNumber(*predictedReduction);
}

// 放大当前方向，仍受传入的最大半径和参数边界限制。
// 前期初始探路可跨过预测上坡区，最终仍按真实形状验收。
static bool buildExpandedTrustRegionStep(const TrustRegionFisher& information,
        const QVector<double>& coordinates, const QVector<double>& originalStep,
        double maximumRadius, QVector<double>* expandedStep, double* prediction,
        bool requirePredictedDescent = true)
{
    const double norm = qSqrt(trustRegionSquaredNorm(originalStep));
    if(norm <= 1.0e-12) return false;
    const double scale = qMin(2.0, maximumRadius / norm);
    if(scale <= 1.01) return false;
    expandedStep->resize(originalStep.size());
    double difference = 0.0;
    for(int i = 0; i < originalStep.size(); ++i) {
        (*expandedStep)[i] = qBound(0.0, coordinates[i] + scale * originalStep[i], 1.0) - coordinates[i];
        difference += qAbs((*expandedStep)[i] - originalStep[i]);
    }
    if(difference <= 1.0e-10) return false;
    *prediction = -trustRegionDotProduct(information.gradient, *expandedStep);
    for(int i = 0; i < expandedStep->size(); ++i)
        for(int j = 0; j < expandedStep->size(); ++j)
            *prediction -= 0.5 * (*expandedStep)[i] * information.matrix[i][j] * (*expandedStep)[j];
    return isFiniteNumber(*prediction) && (!requirePredictedDescent || *prediction > 1.0e-14);
}

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

    DEBUG_OUT("LM optimizer reset");
}

void nmCalculationAutoFitLM::setLayeredSamplingEnabled(bool enabled)
{
    // 拟合运行期间不允许改变残差定义；新一轮由 resetOptimizer 初始化层级。
    if(!m_isRunning) {
        m_layeredSampling = enabled;
    }
}

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

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

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

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

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

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

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

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

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

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

    for(int k = 0; k < kAutoFitTimeWindowCount; ++k) {
        const QString prefix = QString("window_%1_").arg(k + 1);
        cols << prefix + "time_min" << prefix + "time_max"
             << prefix + "weight_sum" << prefix + "rms_error" << prefix + "energy";
    }

    cols << "sampling_mode" << "sampling_stride" << "sampling_points"
         << "full_target_points" << "layer_objective"
         << "pressure_vertical_bias" << "derivative_vertical_bias" << "first_window_shape_loss"
         << "early_value_loss" << "early_parallel_loss" << "early_parallel_bias" << "early_gap_loss" << "early_gap_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\": 39,\n";
    out << "  \"strategy\": \"permeability_height_then_shape_then_joint_lm\",\n";
    out << "  \"height_acceptance\": \"reliable_vertical_loss_decrease; no_shape_or_total_loss_constraint\",\n";
    out << "  \"shape_priority\": \"storage_then_skin_then_shape_without_wellbore_then_optional_wellbore_recheck\",\n";
    out << "  \"shape_metric\": \"pressure_and_derivative_log_slopes_81_points_lag_8\",\n";
    out << "  \"early_parallel_metric\": \"pressure_and_derivative_log_slopes_first_window_81_points_lag_8\",\n";
    out << "  \"early_gap_metric\": \"weighted_rms_log10_pressure_derivative_ratio_error_first_window_81_points\",\n";
    out << "  \"early_gap_bias_metric\": \"weighted_mean_signed_log10_gap_simulation_minus_target_first_window\",\n";
    out << "  \"early_direction_policy\": \"initial_gap_bias_positive_increase_negative_decrease; direction_shared_and_locked_for_storage_and_skin; recompute_only_at_wellbore_recheck_entry\",\n";
    out << "  \"early_initial_step_policy\": \"after_gap_direction_selection; storage_ratio_cap_1.05; skin_local_scale_fraction_0.02; no_cached_probe_acceptance_before_first_shape_improvement\",\n";
    out << "  \"early_initial_expansion_policy\": \"before_first_accepted_shape_improvement; if_shape_worsens_double_coordinate_step_from_unchanged_base_up_to_parameter_bound; retain_base_on_failure\",\n";
    out << "  \"early_wellbore_acceptance\": \"locked_gap_direction_and_early_shape_decrease\",\n";
    out << "  \"early_wellbore_total_constraint\": false,\n";
    out << "  \"shape_stage_total_constraint\": false,\n";
    out << "  \"shape_wellbore_parameters_frozen\": true,\n";
    out << "  \"shape_sensitivity_refresh\": \"non_wellbore_columns_once_at_sweep_entry; secant updates during sweep; full Jacobian rebuilt at total-stage entry\",\n";
    out << "  \"shape_expansion_policy\": \"lock parameter and direction; after each strict shape decrease double next coordinate step from accepted point; otherwise retain point and halve step\",\n";
    out << "  \"shape_radius_policy\": \"initial LM step capped at 0.24; subsequent coordinate steps capped at 0.60; no prediction-ratio gate\",\n";
    out << "  \"shape_global_parameter_selection\": \"largest predicted reduction among unvisited non-wellbore parameters; finish current parameter before selecting next; each parameter visited once\",\n";
    out << "  \"stage2_budget\": \"one parameter sweep; each parameter stops after 3 consecutive non-improving trials, a bound, or minimum step; existing solver failure limit retained\",\n";
    out << "  \"shape_step_improvement_threshold\": \"any strict shape loss decrease; no cumulative improvement threshold\",\n";
    out << "  \"iteration_count_scope\": \"max_iterations applies to total-stage LM; trace iteration remains cumulative\",\n";
    out << "  \"wellbore_recheck_trigger\": \"early loss exceeds initial early-stage exit by max(0.0001, 10%)\",\n";
    out << "  \"wellbore_recheck_acceptance\": \"early decrease with fixed entry shape limit max(0.0001, 5%); no_total_loss_constraint\",\n";
    out << "  \"wellbore_recheck_budget\": \"once when early shape degrades; each wellbore parameter stops after 2 ineffective steps, independent of total-stage budget\",\n";
    out << "  \"total_stage_shape_constraint\": false,\n";
    out << "  \"total_parameter_selection\": \"all_valid_free_columns\",\n";
    out << "  \"skin_difference_policy\": \"local_scale_in_all_stages\",\n";
    out << "  \"stage2_difference_policy\": \"coordinate_step_cap_0.10_with_full_trust_radius; log_parameter_ratio_cap_1.50; skin_local_scale_fraction_0.10; total_stage_unchanged\",\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(m_layeredSampling ? "layered_target" : "fixed_80") << ",\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 << "  \"time_windows\": {\n";
    out << "    \"count\": " << kAutoFitTimeWindowCount << ",\n";
    out << "    \"overlap_ratio\": " << jsonNumber(kAutoFitTimeWindowOverlapRatio) << ",\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 < 7; ++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();
        }
    }

    // 无效候选也补齐窗口列，保证轨迹每行结构一致。
    for(int k = 0; k < kAutoFitTimeWindowCount; ++k) {
        if(objectiveBreakdown && objectiveBreakdown->valid &&
                k < objectiveBreakdown->timeWindows.size()) {
            const AutoFitTimeWindowLM& window = objectiveBreakdown->timeWindows[k];
            cols << traceNumber(window.timeMin) << traceNumber(window.timeMax)
                 << traceNumber(window.weightSum) << traceNumber(window.rmsError)
                 << traceNumber(window.energy);
        } else {
            for(int i = 0; i < 5; ++i) {
                cols << QString();
            }
        }
    }

    cols << (m_layeredSampling ? "layered_target" : "fixed_80")
         << QString::number(objectiveBreakdown ? objectiveBreakdown->samplingStride : m_samplingStride)
         << (objectiveBreakdown ? QString::number(objectiveBreakdown->residualVector.size() / 2) : QString())
         << (objectiveBreakdown ? QString::number(objectiveBreakdown->fullPointCount) : QString())
         << (objectiveBreakdown ? traceNumber(objectiveBreakdown->layerError) : QString())
         << (objectiveBreakdown ? traceNumber(objectiveBreakdown->pressureVerticalBias) : QString())
         << (objectiveBreakdown ? traceNumber(objectiveBreakdown->derivativeVerticalBias) : QString())
         << (objectiveBreakdown && objectiveBreakdown->valid
             ? traceNumber(qSqrt(trustRegionEarlyShapeEnergy(objectiveBreakdown->shapeResiduals))) : QString())
         << (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(objectiveBreakdown->earlyValueLoss) : QString())
         << (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(objectiveBreakdown->earlyParallelLoss) : QString())
         << (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(objectiveBreakdown->earlyParallelBias) : QString())
         << (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(trustRegionEarlyGapLoss(*objectiveBreakdown)) : QString())
         << (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(trustRegionEarlyGapBias(*objectiveBreakdown)) : QString());
    QTextStream out(&m_traceFile);
    out << cols.join(",") << "\n";
    m_traceFile.flush();
}

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

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

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

    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();

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

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

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

            // Dfc 只存在于两类压裂井，普通井在完整向量中保持为 0。
            if(pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
                nmDataVerticalFracturedWell* fracturedWell =
                    dynamic_cast<nmDataVerticalFracturedWell*>(pTargetWell);
                if(fracturedWell) {
                    fullParams[5] = fracturedWell->getDfc().getValue().toDouble();
                    fullParams[6] = fracturedWell->getFractureHalfLength().getValue().toDouble();
                }
            } else if(pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
                nmDataHorizontalFracturedWell* fracturedWell =
                    dynamic_cast<nmDataHorizontalFracturedWell*>(pTargetWell);
                if(fracturedWell) {
                    fullParams[5] = fracturedWell->getDfc().getValue().toDouble();
                    fullParams[6] = 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 只读取迭代次数和目标误差，其他数值控制保留在现有算法实现中。
    nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
    nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        if(enabledParams == 0) {
            m_lastError = 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;
        }

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

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

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

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

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

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

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

            m_initialValues = savedInitialValues;
        }

        finalReason = runTrustRegionFitting();
        validateAndProtectFinalResult();

        // 两种模式共用全目标点比较指标，仅重算已有曲线，不增加真实求解次数。
        if(m_comparisonTimeMin > 0.0 && !m_globalBestLogLogData.isEmpty()) {
            const bool savedMode = m_layeredSampling;
            const int savedStride = m_samplingStride;
            const AutoFitObjectiveBreakdownLM savedBreakdown = m_lastObjectiveBreakdown;
            m_layeredSampling = true;
            m_samplingStride = 1;
            const double comparisonFinal = calculateLogLogCurveError(m_targetLogLogData, m_globalBestLogLogData);
            const double comparisonInitial = m_userInitialLogLogData.isEmpty() ? 1.0e10
                : calculateLogLogCurveError(m_targetLogLogData, m_userInitialLogLogData);
            m_layeredSampling = savedMode;
            m_samplingStride = savedStride;
            m_lastObjectiveBreakdown = savedBreakdown;
            if(comparisonFinal < 1.0e9) {
                emit logMessageGenerated(tr("Full-target comparison error: initial=%1; final=%2")
                    .arg(comparisonInitial < 1.0e9 ? QString::number(comparisonInitial, 'e', 6) : tr("Unavailable"))
                    .arg(comparisonFinal, 0, 'e', 6));
            }
            if(savedMode && savedStride > 1 && finalReason == LM_MAX_ITERATIONS) {
                emit logMessageGenerated(tr("Budget exhausted before the full sampling layer; convergence is not confirmed."));
            }
        }

        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_layeredSampling || m_samplingStride == 1) && 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;
        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;
        emit fittingFinished(false, m_lastError);
        return false;
    }

    // 只有最终完整求解实际执行并提交快照后才能置为成功。
    bool finalFullSolverSucceeded = false;
    bool finalFullSolverExecuted = false;

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

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

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

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

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

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

            if(finalFullSolverExecuted && finalFullSolverSucceeded) {
                emit logMessageGenerated(
                    tr("Parameters and full-field results applied successfully to data manager"));
            } else if(!finalFullSolverExecuted) {
                emit logMessageGenerated(
                    tr("Optimized parameters applied to data manager"));
            }
        } catch(const std::exception& e) {
            finalFullSolverSucceeded = false;
            m_lastError = tr("Failed to apply final parameters: %1").arg(e.what());
            emit logMessageGenerated(
                tr("ERROR: Failed to apply final parameters: %1").arg(e.what()));
        } catch(...) {
            finalFullSolverSucceeded = false;
            m_lastError = tr("Failed to apply final parameters due to unknown error");
            emit logMessageGenerated(
                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);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING SUCCESSFUL ==="));
    } else if(finalReason == LM_TRUE_CONVERGENCE) {
        success = true;
        message = QString(
            tr("Automatic fitting converged to a stable solution. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING CONVERGED ==="));
    } else if(finalReason == LM_LOCAL_OPTIMUM) {
        success = true;
        message = QString(
            tr("Automatic fitting reached a local optimum. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING - LOCAL OPTIMUM ==="));
    } else if(finalReason == LM_MAX_ITERATIONS) {
        success = true;
        message = QString(tr("Total-stage budget reached. Best error: %1, Cumulative iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING - MAX ITERATIONS ==="));
    } else if(finalReason == LM_USER_STOPPED) {
        success = true;
        message = QString(tr("Stopped by user. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING STOPPED BY USER ==="));
    } else if(finalReason == LM_CONSECUTIVE_FAILURES) {
        message = QString(
            tr("Automatic fitting failed due to consecutive failures. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING FAILED ==="));
    } else {
        message = QString(
            tr("Automatic fitting ended unexpectedly. Best error: %1, Iterations: %2"))
                  .arg(m_globalBestFitness, 0, 'e', 4)
                  .arg(m_currentIteration + 1);
        emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING - UNKNOWN END ==="));
    }

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

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

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

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

    m_initialValues.clear();

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

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

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

                break;

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

                break;

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

            case 4: // 初始含水饱和度
                initialValue = reservoirData.getSwi().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->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 6: // 裂缝半长
                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()
{
    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 totalEvaluationBudget = qMax(dimensions + 2, qMax(20, m_maxIterations * 3));
    int maximumEvaluations = m_totalEvaluations + totalEvaluationBudget;
    int totalIterations = 0;
    // 下列步长均位于归一化内部坐标：0.04 表示参数范围的 4%，信赖半径
    // 限制一次联合移动的二范数，相关性门槛用于排除响应近乎共线的参数。
    const double sensitivityStep = 0.04;
    const double minimumCoordinateStep = 1.0e-5;
    const double minimumTrustRadius = 2.0e-3;
    const double maximumTrustRadius = 0.30;
    // 第二阶段整体形状允许更大步幅；第三阶段继续使用原来的半径和预测验收规则。
    const double initialShapeTrustRadius = 0.24;
    const double maximumShapeTrustRadius = 0.60;
    // 误差下降至少达到绝对 1e-5 且相对当前有效基准 0.2% 才算有效改善。
    // 更小的下降仍保留为最佳解，但不能反复清除停滞状态、延长拟合时间。
    const double effectiveRelativeImprovement = 2.0e-3;
    const double effectiveAbsoluteImprovement = 1.0e-5;
    const int maximumIneffectiveSteps = 3;

    // damping 是 LM 阻尼；拒绝或预测失准时增大，真实下降与预测一致时减小。
    // 按连续预测失准刷新 Jacobian；接受步数只作兜底，不因累计移动距离强制重建。
    double trustRadius = 0.12;
    double damping = 1.0e-2;
    int consecutiveRejectedSteps = 0;
    int consecutiveSolverFailures = 0;
    int acceptedSinceRebuild = 0;
    int consecutiveIneffectiveSteps = 0;
    int consecutivePoorPredictions = 0;
    bool rebuildRequested = true;
    // 跟踪文件保留稳定英文原因，界面输出时再翻译，避免受本地编码影响。
    const char* rebuildReason = QT_TR_NOOP("initial sensitivity model");
    bool modelRebuiltAtMinimumRadius = false;
    bool stagnationConfirmationRequested = false;
    bool globalFallbackAttempted = false;
    StopReasonLM stopReason = LM_MAX_ITERATIONS;

    // jacobian 的行对应当前采样层的残差，列对应用户勾选的参数。
    // Fisher 按当前阶段取数值或形状残差行；冻结列不参与形状差分，整体阶段重建全部列。
    QVector<QVector<double> > jacobian;
    QVector<bool> jacobianColumnValid(dimensions, false);

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

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

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

    // 各阶段只发布通过当前目标和约束检查的工作点；曲线和诊断快照必须
    // 与参数同步更新，防止界面显示或最终精英保护使用错配的数据。
    auto publishAcceptedPoint = [&](const TrustRegionEvaluation& evaluation) {
        m_globalBestPosition = evaluation.parameters;
        m_globalBestFitness = evaluation.fitness;
        m_globalBestObjectiveBreakdown = evaluation.breakdown;
        m_globalBestLogLogData = evaluation.curve;
        emit bestCurveUpdated(m_targetLogLogData,
                              m_globalBestLogLogData,
                              m_currentIteration + 1,
                              m_globalBestFitness);
    };

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

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

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

    // 高度预调整只改变用户勾选的渗透率。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;
    if(permeabilityColumn >= 0) {
        emit logMessageGenerated(tr("LM stage 1: align curve height using permeability only; no fixed evaluation limit, effective improvement threshold 10%."));
        // 不限制试算次数，由高度对齐、改善停滞及可行步长决定何时结束。
        while(processPauseAndStop()) {
            const double bias = current.breakdown.verticalCommonBias;
            if(!current.breakdown.verticalReliable || qAbs(bias) <= 0.01 || 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);
            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;
            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(-1, 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;
            emit logMessageGenerated(tr("Permeability alignment: k=%1, height error=%2, shape error=%3, result=%4")
                .arg(candidate.parameters[permeabilityColumn], 0, 'g', 6)
                .arg(candidate.breakdown.verticalLoss, 0, 'e', 4).arg(candidate.breakdown.shapeLoss, 0, 'e', 4)
                .arg(accepted ? tr("accepted") : tr("rejected")));
        }
        QString heightReason = tr("no feasible permeability step or step too small");
        if(m_shouldStop) heightReason = tr("stopped by user");
        else if(!current.breakdown.verticalReliable) heightReason = tr("pressure and derivative height directions conflict");
        else if(current.breakdown.verticalLoss <= 0.01) 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("Permeability alignment ended: %1").arg(heightReason));
    }

    if(m_shouldStop) return LM_USER_STOPPED;

    // 第二阶段先测灵敏度，再逐个调完非井筒参数；每个参数只访问一次。
    // 普通形状调整冻结井储、表皮，只按形状验收；总误差达标只在第三阶段判断。
    bool shapeStage = true;
    bool singleShapeStarted = false;
    int shapeStepCount = 0;
    QVector<bool> shapeParameterFinished(dimensions, false);
    int activeShapeColumn = -1;
    double shapeCoordinateStep = 0.0;
    int consecutiveShapeRejections = 0;
    const int maximumShapeRejections = 3;
    const int storageColumn = m_enabledParamIndices.indexOf(2);
    const int skinColumn = m_enabledParamIndices.indexOf(1);
    bool hasRemainingShapeParameters = false;
    for(int column = 0; column < dimensions; ++column) {
        if(m_enabledParamIndices[column] != 1 && m_enabledParamIndices[column] != 2)
            hasRemainingShapeParameters = true;
    }
    const int wellboreParameterCount = (storageColumn >= 0 ? 1 : 0) + (skinColumn >= 0 ? 1 : 0);
    bool wellboreRecheck = false;
    bool wellboreRecheckDone = false;
    double earlyWellboreBaseline = current.breakdown.earlyParallelLoss;
    double recheckShapeLimit = 0.0;
    // 先井储、再表皮；回检沿用相同顺序和前期指标，不回到形状阶段反复循环。
    bool earlyShapeStage = shapeStage && wellboreParameterCount > 0;
    int earlyShapeParameter = storageColumn >= 0 ? 2 : 1;
    double initialEarlyGapBias = trustRegionEarlyGapBias(current.breakdown);
    double earlyGapDirection = trustRegionEarlyGapDirection(initialEarlyGapBias);
    bool earlyShapeHasImproved = false;
    auto announceEarlyGapDirection = [&]() {
        // 两个参数共用本轮入口的方向；日志保留初始偏差，避免误读为探针响应方向。
        const QString directionName = earlyGapDirection > 0.0 ? "increase" :
            (earlyGapDirection < 0.0 ? "decrease" : "matched");
        writeTraceRow(m_currentIteration, m_enabledParamIndices.indexOf(earlyShapeParameter),
            "early_gap_direction", current.parameters, current.fitness, true, 0,
            QString("%1_initial_gap_bias_%2").arg(directionName).arg(initialEarlyGapBias, 0, 'g', 12),
            &current.breakdown);
        emit logMessageGenerated(tr("Early gap guidance: %1, initial signed gap=%2, direction=%3; fixed for storage and skin.")
            .arg(earlyShapeParameter == 2 ? tr("wellbore storage") : tr("skin"))
            .arg(initialEarlyGapBias, 0, 'g', 6)
            .arg(earlyGapDirection > 0.0 ? tr("increase") :
                (earlyGapDirection < 0.0 ? tr("decrease") : tr("matched"))));
    };
    auto parameterAllowedInStage = [&](int column) -> bool {
        const int parameterIndex = m_enabledParamIndices[column];
        // 候选和缓存灵敏度探针共用此限制，防止探针绕过形状阶段的冻结。
        if(!shapeStage) return true;
        if(earlyShapeStage) return parameterIndex == earlyShapeParameter;
        return parameterIndex != 1 && parameterIndex != 2;
    };
    auto stageError = [&](const TrustRegionEvaluation& point) -> double {
        if(earlyShapeStage) return point.breakdown.earlyParallelLoss;
        return shapeStage ? point.breakdown.shapeLoss : point.fitness;
    };
    auto acceptable = [&](const TrustRegionEvaluation& point, const TrustRegionEvaluation& base) -> bool {
        if(!point.valid) return false;
        if(earlyShapeStage) {
            // 缓存探针也只能沿间距确定的方向接受；停止仍看原来的前期形状改善。
            const int column = m_enabledParamIndices.indexOf(earlyShapeParameter);
            return earlyGapDirection * (point.coordinates[column] - base.coordinates[column]) > 0.0 &&
                point.breakdown.earlyParallelLoss < base.breakdown.earlyParallelLoss &&
                (!wellboreRecheck || point.breakdown.shapeLoss <= recheckShapeLimit);
        }
        return shapeStage
            ? point.breakdown.shapeLoss < base.breakdown.shapeLoss
            // 预调整结束后恢复原 LM：有效候选只按整体误差下降接受。
            : point.fitness < base.fitness;
    };
    auto acceptPoint = [&](const TrustRegionEvaluation& point) {
        // 缓存探针和正式候选共用此标记；首次接受形状改善后，不再启用初始扩步探路。
        if(earlyShapeStage) earlyShapeHasImproved = true;
        current = point;
        publishAcceptedPoint(current);
        restoreEvaluationState(current);
    };
    emit logMessageGenerated(shapeStage
        ? tr("LM stage 2: adjust each shape parameter once; double the step after improvement, halve it after rejection, and switch after 3 consecutive rejections.")
        : tr("LM stage 3: original LM fitting; accept by total error only."));
    if(earlyShapeStage) {
        emit logMessageGenerated(earlyShapeParameter == 2
            ? tr("Early adjustment: storage then skin; pressure-derivative gap selects direction, early shape error controls acceptance and stopping.")
            : tr("Early adjustment: skin; pressure-derivative gap selects direction, early shape error controls acceptance and stopping."));
        announceEarlyGapDirection();
    }

    // 有效改善始终相对“上一次有效改善后的误差”累计判断，避免一连串微小
    // 下降每次都清零计数；累计达到门槛后才开始新的有效改善基准。
    double effectiveImprovementBaseline = stageError(current);
    int fullDataRejections = 0;
    bool samplingRefreshFailed = false;
    auto promoteSampling = [&](bool complete) -> bool {
        if(!m_layeredSampling || m_samplingStride == 1 || m_shouldStop) {
            return false;
        }
        const int previousStride = m_samplingStride;
        const AutoFitObjectiveBreakdownLM previousBreakdown = current.breakdown;
        m_samplingStride = complete ? 1 : m_samplingStride / 2;
        const double refreshedFitness = calculateLogLogCurveError(m_targetLogLogData, current.curve);
        if(refreshedFitness >= 1.0e9 || !trustRegionResidualsValid(m_lastObjectiveBreakdown)) {
            m_samplingStride = previousStride;
            m_lastObjectiveBreakdown = previousBreakdown;
            m_lastError = tr("Failed to refresh the sampling layer from the current curve.");
            samplingRefreshFailed = true;
            return false;
        }
        current.fitness = refreshedFitness;
        current.breakdown = m_lastObjectiveBreakdown;
        publishAcceptedPoint(current);
        restoreEvaluationState(current);
        jacobian.clear();
        jacobianColumnValid.fill(false);
        rebuildRequested = true;
        rebuildReason = QT_TR_NOOP("sampling layer changed");
        trustRadius = 0.12;
        damping = 1.0e-2;
        consecutiveRejectedSteps = 0;
        consecutiveIneffectiveSteps = 0;
        acceptedSinceRebuild = 0;
        consecutivePoorPredictions = 0;
        modelRebuiltAtMinimumRadius = false;
        stagnationConfirmationRequested = false;
        globalFallbackAttempted = false;
        fullDataRejections = 0;
        effectiveImprovementBaseline = stageError(current);
        emit logMessageGenerated(tr("LM sampling refined: %1 / %2 target points; full-target error: %3")
            .arg(current.breakdown.residualVector.size() / 2)
            .arg(current.breakdown.fullPointCount).arg(current.fitness, 0, 'e', 4));
        writeTraceRow(m_currentIteration, -1, "sampling_refinement", current.parameters,
                      current.fitness, true, 0, "rebuild_required", &current.breakdown);
        return true;
    };

    emit logMessageGenerated(m_layeredSampling
        ? tr("LM sampling: layered target points (%1 / %2); total-stage acceptance uses all valid target points.")
            .arg(current.breakdown.residualVector.size() / 2).arg(current.breakdown.fullPointCount)
        : tr("LM sampling: fixed 80 points (original mode)."));

    auto reuseSensitivityForNextStage = [&]() {
        // 更换目标时复用完整 J，清除旧目标下的拒绝和预测失准计数。
        rebuildRequested = jacobian.isEmpty() || consecutiveSolverFailures > 0 || acceptedSinceRebuild >= 10;
        if(jacobian.isEmpty()) rebuildReason = QT_TR_NOOP("no valid sensitivity model");
        else if(consecutiveSolverFailures > 0) rebuildReason = QT_TR_NOOP("solver failure before stage switch");
        else if(acceptedSinceRebuild >= 10) rebuildReason = QT_TR_NOOP("10 accepted steps since last rebuild");
        consecutivePoorPredictions = 0;
        modelRebuiltAtMinimumRadius = false;
        consecutiveIneffectiveSteps = 0;
        consecutiveRejectedSteps = 0;
        consecutiveSolverFailures = 0;
        stagnationConfirmationRequested = false;
        globalFallbackAttempted = false;
        trustRadius = shapeStage && !earlyShapeStage ? initialShapeTrustRadius : 0.12;
        damping = 0.01;
    };
    auto enterTotalStage = [&](const QString& reason) {
        const bool finishedRecheck = wellboreRecheck;
        totalIterations = 0;
        maximumEvaluations = m_totalEvaluations + totalEvaluationBudget;
        shapeStage = false;
        earlyShapeStage = false;
        wellboreRecheck = false;
        effectiveImprovementBaseline = current.fitness;
        reuseSensitivityForNextStage();
        // 第二阶段只测各子阶段需要的列；所有出口都在当前点重建完整 J，不能把缺列模型带入整体 LM。
        jacobian.clear();
        rebuildRequested = true;
        rebuildReason = QT_TR_NOOP("full sensitivity at total-stage entry");
        emit progressUpdated(0, current.fitness);
        emit logMessageGenerated(tr("Total-stage budget starts now: %1 iterations, %2 evaluations; pre-adjustment is counted separately.")
            .arg(m_maxIterations).arg(totalEvaluationBudget));
        emit logMessageGenerated(tr("Shape stage ended: %1").arg(reason));
        emit logMessageGenerated(tr("LM stage 3: original LM fitting; accept by total error only."));
        writeTraceRow(m_currentIteration, -1, "stage_switch", current.parameters, current.fitness,
                      true, 0, finishedRecheck ? "wellbore_recheck_to_total" : "shape_to_total", &current.breakdown);
    };

    auto finishShapeStage = [&](const QString& reason) {
        if(m_shouldStop) return;
        // 只检查一次；用绝对加相对容差区分有意义的退化和接近零时的比例放大。
        if(wellboreRecheckDone || wellboreParameterCount == 0 || !hasRemainingShapeParameters) {
            enterTotalStage(reason);
            return;
        }
        wellboreRecheckDone = true;
        const double earlyTolerance = qMax(1.0e-4, 0.10 * earlyWellboreBaseline);
        if(current.breakdown.earlyParallelLoss <= earlyWellboreBaseline + earlyTolerance) {
            writeTraceRow(m_currentIteration, -1, "stage_switch", current.parameters, current.fitness,
                          true, 0, "wellbore_recheck_skipped_no_degradation", &current.breakdown);
            enterTotalStage(reason);
            return;
        }
        wellboreRecheck = true;
        earlyShapeStage = true;
        earlyShapeParameter = storageColumn >= 0 ? 2 : 1;
        initialEarlyGapBias = trustRegionEarlyGapBias(current.breakdown);
        earlyGapDirection = trustRegionEarlyGapDirection(initialEarlyGapBias);
        earlyShapeHasImproved = false;
        // 回检只约束形状：允许总误差上升，但保护其他参数已调好的整体形状。
        recheckShapeLimit = current.breakdown.shapeLoss + qMax(1.0e-4, 0.05 * current.breakdown.shapeLoss);
        effectiveImprovementBaseline = current.breakdown.earlyParallelLoss;
        // 其他参数已改变，回检入口重测形状灵敏度，方向由本轮初始间距重新确定。
        jacobian.clear();
        reuseSensitivityForNextStage();
        rebuildRequested = true;
        rebuildReason = QT_TR_NOOP("new sensitivity for wellbore recheck");
        emit logMessageGenerated(tr("Wellbore recheck: early shape error=%1, shape limit=%2.")
            .arg(current.breakdown.earlyParallelLoss, 0, 'g', 6)
            .arg(recheckShapeLimit, 0, 'g', 6));
        writeTraceRow(m_currentIteration, -1, "stage_switch", current.parameters, current.fitness,
                      true, 0, "shape_to_wellbore_recheck", &current.breakdown);
        announceEarlyGapDirection();
    };

    auto finishEarlyShapeStage = [&](const QString& reason) {
        if(earlyShapeParameter == 2 && skinColumn >= 0) {
            earlyShapeParameter = 1;
            earlyShapeHasImproved = false;
            effectiveImprovementBaseline = current.breakdown.earlyParallelLoss;
            reuseSensitivityForNextStage();
            // 井储已变，重测表皮的形状灵敏度和步幅，方向仍保持本轮入口的判断。
            rebuildRequested = true;
            rebuildReason = QT_TR_NOOP("refresh shape sensitivity at skin adjustment entry");
            emit logMessageGenerated(tr("Wellbore storage adjustment ended: %1; now adjusting skin.").arg(reason));
            writeTraceRow(m_currentIteration, -1, "stage_switch", current.parameters, current.fitness,
                          true, 0, wellboreRecheck ? "recheck_storage_to_skin" : "storage_to_skin", &current.breakdown);
            announceEarlyGapDirection();
            return;
        }
        if(wellboreRecheck) {
            enterTotalStage(reason);
            return;
        }
        earlyShapeStage = false;
        earlyWellboreBaseline = current.breakdown.earlyParallelLoss;
        emit logMessageGenerated(tr("Early wellbore adjustment ended: %1").arg(reason));
        if(!hasRemainingShapeParameters) {
            enterTotalStage(tr("no remaining shape parameters"));
            return;
        }
        effectiveImprovementBaseline = current.breakdown.shapeLoss;
        reuseSensitivityForNextStage();
        emit logMessageGenerated(tr("Shape fitting continues with storage and skin fixed; a conditional wellbore recheck follows."));
        writeTraceRow(m_currentIteration, -1, "stage_switch", current.parameters, current.fitness,
                      true, 0, "early_wellbore_to_shape", &current.breakdown);
    };

    auto finishShapeParameter = [&](const QString& reason) {
        // 已完成的参数不再参与后续选优，保证所有可调参数只扫一遍。
        if(activeShapeColumn >= 0) {
            shapeParameterFinished[activeShapeColumn] = true;
            writeTraceRow(m_currentIteration, activeShapeColumn, "shape_parameter_end",
                current.parameters, current.fitness, true, 0, reason, &current.breakdown);
        }
        activeShapeColumn = -1;
        shapeCoordinateStep = 0.0;
        consecutiveShapeRejections = 0;
        for(int column = 0; column < dimensions; ++column) {
            if(parameterAllowedInStage(column) && !shapeParameterFinished[column]) return;
        }
        finishShapeStage(tr("all shape parameters visited once"));
    };
    auto completeShapeSearchStep = [&](bool accepted) {
        if(!shapeStage || earlyShapeStage || m_shouldStop) return;
        // 只按本次真实下降控制实际步长；不再使用累计 10% 门槛提前结束整轮。
        ++shapeStepCount;
        consecutiveShapeRejections = accepted ? 0 : consecutiveShapeRejections + 1;
        shapeCoordinateStep = qBound(-maximumShapeTrustRadius,
            shapeCoordinateStep * (accepted ? 2.0 : 0.5), maximumShapeTrustRadius);
        writeTraceRow(m_currentIteration, activeShapeColumn, "shape_step_end", current.parameters,
            current.fitness, true, 0,
            QString("%1_step_%2_rejections_%3_next_delta_%4")
                .arg(accepted ? "improved" : "rejected").arg(shapeStepCount)
                .arg(consecutiveShapeRejections).arg(shapeCoordinateStep, 0, 'g', 12), &current.breakdown);
        if(consecutiveShapeRejections >= maximumShapeRejections)
            finishShapeParameter("3_consecutive_non_improving_trials");
        else if(qAbs(shapeCoordinateStep) < minimumCoordinateStep)
            finishShapeParameter("minimum_coordinate_step");
    };

    auto registerEffectiveImprovement = [&](double fitness) -> bool {
        if(shapeStage && !earlyShapeStage) return false;
        const double requiredImprovement = qMax(
            shapeStage ? 1.0e-4 : effectiveAbsoluteImprovement,
            qAbs(effectiveImprovementBaseline) *
            (shapeStage ? 0.01 : effectiveRelativeImprovement));
        const double improvement = effectiveImprovementBaseline - fitness;
        if(improvement < requiredImprovement) {
            return false;
        }

        // 有效改善后更新累计基准，并清除停滞计数。
        globalFallbackAttempted = false;
        effectiveImprovementBaseline = fitness;
        consecutiveIneffectiveSteps = 0;
        stagnationConfirmationRequested = false;
        return true;
    };

    // 普通形状的失败试算缩步并累计当前参数的拒绝次数；
    // 前期井储、表皮与整体 LM 沿用各自的无效次数处理。
    auto recordIneffectiveStep = [&]() -> bool {
        if(shapeStage && !earlyShapeStage) {
            completeShapeSearchStep(false);
            return false;
        }
        ++consecutiveIneffectiveSteps;
        if(earlyShapeStage) {
            // 一次拒绝只缩步，再给当前参数一次真实尝试；连续无改善才交给下一段。
            if(consecutiveIneffectiveSteps >= 2)
                finishEarlyShapeStage(tr("2 consecutive steps without effective early improvement"));
            return false;
        }
        if(!globalFallbackAttempted) {
            return false;
        }
        globalFallbackAttempted = false;
        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;
    };

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

    if(!shapeStage && current.fitness < m_targetError) {
        promoteSampling(true);
        return samplingRefreshFailed ? LM_OPTIMIZATION_FAILED : LM_TARGET_ACHIEVED;
    }

    // 在同一个真实工作点逐参数做单边差分。失败先在原方向缩步，再反向尝试；
    // 正常情况下每列仍只需一次真实求解，重试也计入总评价预算。
    auto rebuildSensitivity = [&]() -> bool {
        // 第二阶段差分不受第三阶段额度限制；每列仍只有有限的缩步和反向重试。
        const int evaluationLimit = shapeStage ? (std::numeric_limits<int>::max)() : maximumEvaluations;
        const int sensitivityColumnCount = shapeStage
            ? (earlyShapeStage ? wellboreParameterCount : dimensions - wellboreParameterCount) : dimensions;
        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;
        // 第二阶段允许内部坐标范围的 10% 及完整信赖半径，避免较大比例试探被旧上限截小。
        // 第三阶段仍为 4% 及半个信赖半径；均保留 0.5% 下限以减少数值噪声影响。
        const double finiteDifferenceStep = qMin(
            shapeStage ? 0.10 : sensitivityStep,
            qMax(5.0e-3, trustRadius * (shapeStage ? 1.0 : 0.5)));

        for(int column = 0;
                column < dimensions &&
                m_totalEvaluations < evaluationLimit &&
                processPauseAndStop();
                ++column) {
            const int parameterIndex = m_enabledParamIndices[column];
            if(earlyShapeStage && parameterIndex != 1 && parameterIndex != 2) continue;
            // 普通形状阶段跳过冻结的井储、表皮，避免无用正演；回检和整体阶段会重新差分。
            if(shapeStage && !earlyShapeStage && !parameterAllowedInStage(column)) continue;
            const double lower = m_parameterLower[parameterIndex];
            const double upper = m_parameterUpper[parameterIndex];
            if(upper <= lower) continue;
            // 第二阶段扩大形状试探范围；第三阶段仍沿用原来的局部差分尺度。
            // 表皮包含零和负值，按物理尺度扰动；正值参数按对数比例限制幅度。
            double localStep = finiteDifferenceStep;
            if(parameterIndex == 1) {
                const double skinStepFraction = shapeStage ? 0.10 : 0.02;
                localStep = qMin(localStep, skinStepFraction * qMax(0.1, qAbs(base.parameters[column])) / (upper - lower));
            } else if(shapeStage && useTrustRegionLogScale(parameterIndex, lower, upper)) {
                localStep = qMin(localStep, qLn(1.50) / (qLn(upper) - qLn(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,
                                  (wellboreRecheck ? "wellbore_recheck_" : (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;

                    // 后半段形状探针仅用于灵敏度，正式调整必须经过单参数 LM 预测选优。
                    // 前期同时测井储和表皮，但缓存探针不能绕过当前子阶段的选参限制。
                    // 前期首次正式调整必须从小步开始，测形状灵敏度的大探针不能提前成为工作点。
                    if((!shapeStage || earlyShapeStage) &&
                            (!earlyShapeStage || earlyShapeHasImproved) &&
                            parameterAllowedInStage(column) && acceptable(probe, base) &&
                            (!bestProbe.valid || stageError(probe) < stageError(bestProbe))) {
                        bestProbe = probe;
                        bestProbeColumn = column;
                        bestProbeDelta = delta;
                    }
                    restoreEvaluationState(base);
                }
            }
        }

        int validColumnCount = 0;
        bool hasActiveSensitivity = false;
        for(int i = 0; i < jacobianColumnValid.size(); ++i) {
            if(jacobianColumnValid[i]) {
                ++validColumnCount;
                if(parameterAllowedInStage(i)) hasActiveSensitivity = true;
            }
        }
        // 冻结列成功不能掩盖所有可调形状列的求解失败；有效零梯度交给逐参数选步处理。
        if(validColumnCount == 0 || (shapeStage && !earlyShapeStage && !hasActiveSensitivity) || m_shouldStop) {
            restoreEvaluationState(base);
            return false;
        }

        // 前期井储/表皮和整体阶段仍可接受更优缓存探针；后半段形状只由 LM 选步。
        // 所有列先基于同一个 base 建完，再用已知割线平移 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,
                          wellboreRecheck ? "wellbore_recheck_accepted_cached_probe" :
                              (earlyShapeStage ? "early_shape_accepted_cached_probe" :
                              (shapeStage ? "shape_accepted_cached_probe" : "total_accepted_cached_probe")),
                          &current.breakdown);
            emit logMessageGenerated(
                (earlyShapeStage ? tr("Sensitivity probe accepted: early shape error=%1")
                                 : (shapeStage ? tr("Sensitivity probe accepted: shape error=%1")
                                               : tr("Sensitivity probe accepted: total error=%1")))
                .arg(stageError(current), 0, 'e', 4));
        } else {
            restoreEvaluationState(current);
        }

        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(sensitivityColumnCount));
        return true;
    };

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

        if(!processPauseAndStop()) {
            break;
        }
        if(wellboreRecheck && current.breakdown.earlyParallelLoss <=
                earlyWellboreBaseline + qMax(1.0e-4, 0.10 * earlyWellboreBaseline))
            enterTotalStage(tr("early shape error restored within tolerance"));
        if(shapeStage && !earlyShapeStage && !singleShapeStarted) {
            singleShapeStarted = true;
            // 初调已改变井储/表皮，普通形状必须在新的工作点重测其余可调参数的响应。
            jacobian.clear();
            reuseSensitivityForNextStage();
            rebuildRequested = true;
            rebuildReason = QT_TR_NOOP("fresh sensitivity at single-parameter shape entry");
        }
        if(m_shouldStop) break;
        if(!shapeStage && (totalIterations >= m_maxIterations || m_totalEvaluations >= maximumEvaluations)) break;
        if(!shapeStage && m_layeredSampling && m_samplingStride > 1) {
            const int remainingEvaluations = maximumEvaluations - m_totalEvaluations;
            // 已准备重建时，提前计入本轮差分及候选的开销，避免紧接着因预算再加密。
            const bool reserveFinalBudget = remainingEvaluations <= 2 * (dimensions + 1) ||
                (rebuildRequested && remainingEvaluations <= 3 * (dimensions + 1));
            const int layerDeadline = qMax(1, m_maxIterations * (m_samplingStride == 4 ? 1 : 2) / 3);
            const bool mergeRefinement = rebuildRequested && totalIterations < layerDeadline &&
                (totalIterations + 1 >= layerDeadline || reserveFinalBudget);
            if(reserveFinalBudget || totalIterations >= layerDeadline || mergeRefinement) {
                const bool promoted = promoteSampling(reserveFinalBudget);
                if(samplingRefreshFailed) break;
                if(promoted && mergeRefinement) {
                    rebuildReason = QT_TR_NOOP("sampling refinement merged with pending rebuild");
                }
            }
        }
        if(rebuildRequested) {
            emit logMessageGenerated(tr("Rebuilding sensitivity model: %1").arg(tr(rebuildReason)));
            writeTraceRow(m_currentIteration, -1, "sensitivity_rebuild", current.parameters,
                          current.fitness, true, 0, rebuildReason, &current.breakdown);
            if(!rebuildSensitivity()) {
                if(shapeStage && !m_shouldStop) {
                    if(earlyShapeStage) {
                        finishEarlyShapeStage(tr("no valid early sensitivity model"));
                        rebuildRequested = true;
                        rebuildReason = QT_TR_NOOP("no valid sensitivity model");
                        continue;
                    }
                    // 真实差分全失败不能冒充停滞收敛。
                    m_lastError = tr("Unable to build a valid shape sensitivity model.");
                    stopReason = LM_OPTIMIZATION_FAILED;
                    break;
                }
                if(promoteSampling(false)) continue;
                stopReason = m_shouldStop
                             ? LM_USER_STOPPED
                             : LM_LOCAL_OPTIMUM;
                break;
            }
            if(!shapeStage && current.fitness < m_targetError) {
                promoteSampling(true);
                stopReason = LM_TARGET_ACHIEVED;
                break;
            }
            if(!shapeStage && m_totalEvaluations >= maximumEvaluations) {
                stopReason = LM_MAX_ITERATIONS;
                break;
            }
            if(wellboreRecheck && current.breakdown.earlyParallelLoss <=
                    earlyWellboreBaseline + qMax(1.0e-4, 0.10 * earlyWellboreBaseline)) {
                enterTotalStage(tr("wellbore recheck completed during sensitivity evaluation"));
                --iteration;
                continue;
            }
            // 差分探针不等于 LM 候选；重建后继续按当前阶段求步并真实评价，再确认停滞。
            registerEffectiveImprovement(stageError(current));
        }

        // 第三阶段独立计数；方向不可行也消耗一次局部尝试，不能无限缩步循环。
        if(!shapeStage) ++totalIterations;
        // 每次从最新 J 和当前阶段残差重算全局信息，供单参数或联合 LM 求步。
        const QVector<double> objectiveResidual = earlyShapeStage
            ? current.breakdown.earlyParallelResiduals
            : (shapeStage ? current.breakdown.shapeResiduals : current.breakdown.residualVector);
        const int rowOffset = earlyShapeStage
            ? current.breakdown.residualVector.size() + current.breakdown.shapeResiduals.size() +
                current.breakdown.earlyValueResiduals.size()
            : (shapeStage ? current.breakdown.residualVector.size() : 0);
        const QVector<QVector<double> > objectiveJacobian = jacobian.mid(rowOffset, objectiveResidual.size());
        QVector<double> objectiveCoordinates = current.breakdown.sampleCoordinates;
        if(shapeStage) {
            objectiveCoordinates.clear();
            for(int i = 0; i < 73; ++i) objectiveCoordinates.append((i + 4.0) / 80.0);
        }
        const QVector<TrustRegionFisher> information = buildTrustRegionFisher(
            objectiveJacobian, objectiveResidual, jacobianColumnValid, objectiveCoordinates);
        const TrustRegionFisher& global = information[kAutoFitTimeWindowCount];
        // 前期残差已含第一窗口权重，使用全局和，不能再次乘局部窗口权重。
        const TrustRegionFisher& stepInformation = global;
        QVector<bool> stageColumnValid = jacobianColumnValid;
        for(int column = 0; column < dimensions; ++column) {
            if(!parameterAllowedInStage(column)) stageColumnValid[column] = false;
        }
        QVector<int> selectedColumns;
        QVector<double> coordinateStep;
        double predictedReduction = 0.0;
        const int selectedWindow = earlyShapeStage ? 0 : -1;

        // 前期只开放当前井储或表皮；后半段形状从有效自由列中选一个，整体 LM 仍联合求解。
        // 保留弱敏感、相关及边界列供预测比较；窗口只保留作误差诊断。
        for(int column = 0; column < dimensions; ++column) {
            if(stageColumnValid[column] && (earlyShapeStage || global.matrix[column][column] > 0.0))
                selectedColumns.append(column);
        }
        if(earlyShapeStage) {
            double earlyStepRadius = trustRadius;
            if(!earlyShapeHasImproved && selectedColumns.size() == 1) {
                // 首次调整沿用小步起调，形状变差时由后面的同向扩步循环处理。
                const int column = selectedColumns[0];
                const int parameterIndex = m_enabledParamIndices[column];
                const double lower = m_parameterLower[parameterIndex];
                const double upper = m_parameterUpper[parameterIndex];
                const double smallStep = useTrustRegionLogScale(parameterIndex, lower, upper)
                    ? qLn(1.05) / (qLn(upper) - qLn(lower))
                    : (parameterIndex == 1 ? 0.02 * qMax(0.1, qAbs(current.parameters[column]))
                                          : 0.05 * qMax(1.0e-8, qAbs(current.parameters[column]))) / (upper - lower);
                earlyStepRadius = qMin(trustRadius, smallStep);
            }
            if(selectedColumns.isEmpty() || !buildEarlyGapGuidedStep(
                    stepInformation, selectedColumns, current.coordinates, earlyGapDirection,
                    damping, earlyStepRadius, minimumCoordinateStep, &coordinateStep, &predictedReduction)) {
                selectedColumns.clear();
            }
        } else if(shapeStage) {
            // 首步由灵敏度选出尚未调过的参数；后续锁定方向，直接使用扩缩后的步长。
            if(activeShapeColumn < 0) {
                QVector<int> availableColumns;
                for(int i = 0; i < selectedColumns.size(); ++i) {
                    if(!shapeParameterFinished[selectedColumns[i]]) availableColumns.append(selectedColumns[i]);
                }
                buildBestSingleParameterStep(global, availableColumns, current.coordinates,
                    0.01, initialShapeTrustRadius, minimumCoordinateStep,
                    &selectedColumns, &coordinateStep, &predictedReduction);
                if(!selectedColumns.isEmpty()) {
                    activeShapeColumn = selectedColumns[0];
                    shapeCoordinateStep = coordinateStep[activeShapeColumn];
                    consecutiveShapeRejections = 0;
                    writeTraceRow(m_currentIteration, activeShapeColumn, "shape_parameter_start",
                        current.parameters, current.fitness, true, 0,
                        "largest_predicted_reduction_unvisited", &current.breakdown);
                } else {
                    // 剩余参数均无可行预测下降步时逐项记录，不能重复选择已完成参数。
                    for(int column = 0; column < dimensions; ++column) {
                        if(!parameterAllowedInStage(column) || shapeParameterFinished[column]) continue;
                        activeShapeColumn = column;
                        finishShapeParameter(jacobianColumnValid[column]
                            ? "no_feasible_descent_step" : "invalid_sensitivity");
                        if(!shapeStage || earlyShapeStage) break;
                    }
                    --iteration;
                    continue;
                }
            } else {
                selectedColumns.clear();
                selectedColumns.append(activeShapeColumn);
                coordinateStep.fill(0.0, dimensions);
                shapeCoordinateStep = qBound(0.0,
                    current.coordinates[activeShapeColumn] + shapeCoordinateStep, 1.0)
                    - current.coordinates[activeShapeColumn];
                if(qAbs(shapeCoordinateStep) < minimumCoordinateStep) {
                    finishShapeParameter("parameter_bound_or_minimum_step");
                    --iteration;
                    continue;
                }
                coordinateStep[activeShapeColumn] = shapeCoordinateStep;
                predictedReduction = -global.gradient[activeShapeColumn] * shapeCoordinateStep
                    - 0.5 * global.matrix[activeShapeColumn][activeShapeColumn]
                        * shapeCoordinateStep * shapeCoordinateStep;
            }
        } else {
            // 第三阶段保持全部有效自由参数的联合 LM 调整。
            globalFallbackAttempted = true;
            if(selectedColumns.isEmpty() || !buildTrustRegionFisherStep(
                    global, selectedColumns, current.coordinates, damping,
                    trustRadius, minimumCoordinateStep, &coordinateStep, &predictedReduction)) {
                selectedColumns.clear();
            }
        }

        // 当前阶段没有可行下降步时，收缩半径并进入原有重建/收敛处理。
        if(selectedColumns.isEmpty()) {
            if(earlyShapeStage) {
                finishEarlyShapeStage(tr("no feasible early adjustment direction or parameter at bound"));
                // 尚未求解候选，原迭代留给其余参数，不额外占用迭代或求解预算。
                --iteration;
                continue;
            }
            if(trustRadius <= minimumTrustRadius * 1.01 &&
                    modelRebuiltAtMinimumRadius) {
                if(promoteSampling(false)) continue;
                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()) {
                if(promoteSampling(false)) continue;
                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 = (wellboreRecheck ? QString("recheck_") : QString()) + QString(earlyShapeStage
            ? (earlyShapeParameter == 2 ? "storage_parallel_" : "skin_parallel_") : (shapeStage ? "shape_" : "total_")) + (selectedWindow >= 0
            ? QString("window_%1").arg(selectedWindow + 1) : QString("global")) +
            "_params_" + selectedParameterIndices.join("_");

        TrustRegionEvaluation candidate;
        candidate.coordinates = candidateCoordinates;
        candidate.parameters = parametersFromCoordinates(candidate.coordinates);
        if(earlyShapeStage && earlyShapeParameter == 2) {
            emit logMessageGenerated(tr("Wellbore storage trial: early shape error=%1, C=%2 -> %3")
                .arg(current.breakdown.earlyParallelLoss, 0, 'g', 5)
                .arg(current.parameters[storageColumn], 0, 'g', 6)
                .arg(candidate.parameters[storageColumn], 0, 'g', 6));
        }
        candidate.valid = evaluateTrustRegionPoint(
            candidate.parameters,
            &candidate.fitness,
            &candidate.breakdown,
            &candidate.curve,
            &candidate.elapsedMs);

        // 初始形状变差时先沿锁定方向跨大步探路，不立即缩步或累计无效次数。
        // 始终以未移动的 current 为基点；每次翻倍并裁到边界，差解不发布、不接受。
        bool earlyInitialExpanded = false;
        bool earlyInitialExpansionAtBound = false;
        while(earlyShapeStage && !earlyShapeHasImproved && candidate.valid &&
                candidate.breakdown.earlyParallelLoss > current.breakdown.earlyParallelLoss &&
                processPauseAndStop()) {
            QVector<double> expandedStep;
            double expandedPrediction = 0.0;
            if(!buildExpandedTrustRegionStep(stepInformation, current.coordinates, coordinateStep,
                    1.0, &expandedStep, &expandedPrediction, false)) {
                earlyInitialExpansionAtBound = true;
                break;
            }
            writeTraceRow(m_currentIteration, selectedColumns[0], "early_initial_expansion",
                candidate.parameters, candidate.fitness, true, candidate.elapsedMs,
                "shape_worse_expand_same_direction", &candidate.breakdown);
            restoreEvaluationState(current);
            TrustRegionEvaluation expanded;
            expanded.coordinates = current.coordinates;
            for(int column = 0; column < dimensions; ++column)
                expanded.coordinates[column] += expandedStep[column];
            expanded.parameters = parametersFromCoordinates(expanded.coordinates);
            const int column = selectedColumns[0];
            emit logMessageGenerated(tr("Early shape initially worsened: %1, trial=%2 -> %3; expanding in the same direction.")
                .arg(earlyShapeParameter == 2 ? tr("wellbore storage") : tr("skin"))
                .arg(candidate.parameters[column], 0, 'g', 6)
                .arg(expanded.parameters[column], 0, 'g', 6));
            expanded.valid = evaluateTrustRegionPoint(expanded.parameters, &expanded.fitness,
                &expanded.breakdown, &expanded.curve, &expanded.elapsedMs);
            candidate = expanded;
            coordinateStep = expandedStep;
            predictedReduction = expandedPrediction;
            stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
            earlyInitialExpanded = true;
        }
        if(earlyInitialExpanded) selectionName += "_initial_expanded";
        if(m_shouldStop) {
            restoreEvaluationState(current);
            stopReason = LM_USER_STOPPED;
            break;
        }

        if(!candidate.valid) {
            // 求解失败的候选不能改变 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 && (!shapeStage || earlyShapeStage)) {
                rebuildRequested = true;
                rebuildReason = QT_TR_NOOP("2 consecutive solver failures");
            }
            if(consecutiveSolverFailures >= m_maxConsecutiveFailures) {
                stopReason = LM_CONSECUTIVE_FAILURES;
                break;
            }
            if(recordIneffectiveStep()) {
                if(promoteSampling(false)) continue;
                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 = earlyShapeStage
            ? candidate.breakdown.earlyParallelResiduals
            : (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 && (!shapeStage || earlyShapeStage)) {
            rebuildRequested = true;
            rebuildReason = QT_TR_NOOP("2 consecutive steps with actual improvement below 25% of prediction");
        }
        bool accepted = acceptable(candidate, current);
        // 粗层认为下降而完整数据不认可时累计，连续两次就提前加密。
        if(!shapeStage && m_layeredSampling && m_samplingStride > 1 && !accepted && actualReduction > 0.0) {
            ++fullDataRejections;
        } else {
            fullDataRejections = 0;
        }
        QString componentName = earlyShapeStage
            ? (earlyShapeParameter == 2 ? "storage_parallel" : "skin_parallel") : (shapeStage ? "shape" : "total");

        if(accepted) {
            // 当前阶段接受候选后同步发布参数、曲线和诊断。
            // 模型预测可靠时减小阻尼并可扩大半径，预测较差时保守收缩。
            acceptPoint(candidate);
            ++acceptedSinceRebuild;
            consecutiveRejectedSteps = 0;

            // 井筒调整及第三阶段保留原 LM 控制；逐参数粗调在步末直接扩缩实际步长。
            if(!shapeStage || earlyShapeStage) {
                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 && (!shapeStage || earlyShapeStage)) {
                rebuildRequested = true;
                rebuildReason = QT_TR_NOOP("10 accepted steps since last rebuild");
            }
            modelRebuiltAtMinimumRadius = false;
        } else {
            // 拒绝时 candidate 只保留在 trace 中，DataManager 和内存状态都恢复
            // 到 current。形状上限或全目标点验收也可能拒绝候选，不能仅凭拒绝
            // 次数认定模型失准；重建由上面的预测质量判断，约束冲突先缩步。
            ++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);

        QString componentDisplayName = componentName;
        if(componentName == "vertical") {
            componentDisplayName = tr("vertical deviation");
        } else if(componentName == "horizontal") {
            componentDisplayName = tr("horizontal deviation");
        } else if(componentName == "storage_parallel") {
            componentDisplayName = tr("wellbore storage early shape");
        } else if(componentName == "skin_parallel") {
            componentDisplayName = tr("skin early shape");
        } else if(componentName == "shape") {
            componentDisplayName = tr("shape deviation");
        } else if(componentName == "total") {
            componentDisplayName = tr("total error");
        }

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

        // 先记录当前试调结果，再发布阶段切换，避免日志显示为未试调就结束。
        if(earlyInitialExpansionAtBound && !accepted) {
            finishEarlyShapeStage(tr("initial same-direction expansion reached the parameter bound without early shape improvement"));
            continue;
        }
        if(shapeStage && !earlyShapeStage) {
            completeShapeSearchStep(accepted);
            continue;
        }
        const bool effectiveImprovement = registerEffectiveImprovement(stageError(current));
        if(!effectiveImprovement && recordIneffectiveStep()) stopReason = LM_LOCAL_OPTIMUM;

        if(stopReason == LM_LOCAL_OPTIMUM || fullDataRejections >= 2) {
            if(promoteSampling(false)) {
                stopReason = LM_MAX_ITERATIONS;
                continue;
            }
            if(stopReason == LM_LOCAL_OPTIMUM || samplingRefreshFailed) break;
        }

        if(!shapeStage && current.fitness < m_targetError) {
            promoteSampling(true);
            stopReason = LM_TARGET_ACHIEVED;
            break;
        }
        if(!shapeStage && selectedWindow < 0 && trustRadius <= minimumTrustRadius * 1.01 &&
                consecutiveRejectedSteps >= 2 && consecutivePoorPredictions >= 2) {
            if(modelRebuiltAtMinimumRadius) {
                if(promoteSampling(false)) continue;
                stopReason = LM_LOCAL_OPTIMUM;
                break;
            }
            rebuildRequested = true;
            rebuildReason = QT_TR_NOOP("inaccurate model at minimum trust radius");
        }
    }

    if(completedIterations > 0) {
        m_currentIteration = completedIterations - 1;
    }
    restoreEvaluationState(current);
    emit logMessageGenerated(tr("Adaptive fitting counts: %1 shape attempts, %2 total-stage iterations, %3 total evaluations.")
        .arg(shapeStepCount).arg(totalIterations).arg(m_totalEvaluations));

    if(m_shouldStop) {
        return LM_USER_STOPPED;
    }
    if(samplingRefreshFailed) return LM_OPTIMIZATION_FAILED;
    if(!shapeStage && current.fitness < m_targetError) {
        promoteSampling(true);
        return samplingRefreshFailed ? LM_OPTIMIZATION_FAILED : LM_TARGET_ACHIEVED;
    }
    if(stopReason == LM_CONSECUTIVE_FAILURES ||
            stopReason == LM_LOCAL_OPTIMUM ||
            stopReason == LM_OPTIMIZATION_FAILED) {
        return stopReason;
    }
    return LM_MAX_ITERATIONS;
}

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

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

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

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

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

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

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

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

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

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

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

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

                                break;

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

                                break;

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

                                break;
                        }
                    }
                }
            }

            return 1e10;
        }

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

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

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

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

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

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

        // 4. 运行求解器。真实求解器偶发失败时允许重试，避免一次 DLL 调用异常
        // 直接让初始评价失败；迭代候选不重复相同参数，交给信赖域缩步。
        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 = runSolver();

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

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

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

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

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

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

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

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

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

                return 1e10;
            }

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

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

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

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

            error = calculateLogLogCurveError(m_targetLogLogData, resultLogLogData);

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

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

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

        return error;

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

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

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

    updateReservoirParameters(parameters);
    updateWellParameters(parameters);
}

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

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

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

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

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

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

            paramIndex++;
        }
    }

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

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

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

    if(!pWell) return;

    // 先在参数副本中组装本次候选值，全部解析完成后再写回现有井对象。
    // 这样不会触发整井赋值，也不会删除并重建井内已有的射孔对象。
    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 5: { // 裂缝导流能力
                        if(pVerticalFracturedWell || pHorizontalFracturedWell) {
                            dfcAttr.setValue(value);
                            updateDfc = true;
                        }
                    }
                    break;

                case 6: { // 裂缝半长
                        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());
        }
    }
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    return count;
}

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

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

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

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

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

    return true;
}

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

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

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

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

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

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

    return true;
}

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

    bool allValid = true;

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

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

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

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

    return allValid;
}

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

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

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

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

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

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

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

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

    return true;
}

double nmCalculationAutoFitLM::calculateLogLogCurveError(
    const QVector<QVector<double> >& target,
    const QVector<QVector<double> >& result)
{
    // 首次有效评价后固定公共时间范围。窗口与上下、左右、形状只负责诊断，
    // 主目标仍由完整压力和导数残差计算，避免窗口重叠造成重复计权。
    // 整个计算过程均位于 log(time)-log(value) 坐标。
    const double invalidLoss = 1.0e10;
    const double valueFloor = 1.0e-12;
    const int numPoints = 80;
    m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();

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

    try {
        // 非正导数无法进入双对数空间，跳过对应采样行，使用剩余有效点比较。
        auto prepareCurve = [valueFloor](const QVector<QVector<double> >& data,
                                         int firstIndex,
                                         QVector<QPointF>* pressure,
                                         QVector<QPointF>* derivative) -> bool {
            if(!pressure || !derivative || data.size() < 3 ||
                    data[0].size() != data[1].size() ||
                    data[0].size() != data[2].size() ||
                    firstIndex < 0 || firstIndex >= data[0].size()) {
                return false;
            }

            for(int i = firstIndex; i < data[0].size(); ++i) {
                if(!isFiniteNumber(data[0][i]) ||
                        !isFiniteNumber(data[1][i]) ||
                        !isFiniteNumber(data[2][i]) ||
                        data[0][i] <= 0.0 ||
                        data[1][i] <= 0.0 ||
                        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;
        }

        auto populateStageMetrics = [&](AutoFitObjectiveBreakdownLM* objective) -> bool {
            // 81 个固定对数时间点只用于插值评价，不增加求解点数或 DLL 调用。
            // 斜率用跨度为全时域 10% 的差分，避免相邻点噪声；两条曲线等权。
            const int count = 81;
            const int lag = 8;
            const double span = qLn(overlapMaxX) - qLn(overlapMinX);
            QVector<double> residuals[2];
            QVector<double> targetLogs[2], resultLogs[2];
            double biases[2] = {0.0, 0.0};
            for(int component = 0; component < 2; ++component) {
                for(int i = 0; i < count; ++i) {
                    const double time = i == 0 ? overlapMinX : (i == count - 1 ? overlapMaxX
                        : qExp(qLn(overlapMinX) + span * i / (count - 1)));
                    double targetValue = 0.0, resultValue = 0.0;
                    if(!interpolateLogValue(component == 0 ? targetPressure : targetDerivative,
                                            time, &targetValue) ||
                       !interpolateLogValue(component == 0 ? resultPressure : resultDerivative,
                                            time, &resultValue)) return false;
                    const double residual = resultValue - targetValue;
                    residuals[component].append(residual);
                    targetLogs[component].append(targetValue);
                    resultLogs[component].append(resultValue);
                    // 对数时间梯形权重等价于互补窗口加权求和，避免密集段主导高度。
                    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));
            populateEarlyWellboreMetrics(objective, targetLogs, resultLogs, span);
            return isFiniteNumber(objective->shapeLoss);
        };

        if(m_layeredSampling) {
            // 完整基准只取公共范围内的目标原始时间点；模拟点数不改变评价标准。
            QVector<double> times;
            QVector<double> pressureResidual;
            QVector<double> derivativeResidual;
            for(int i = 0; i < targetPressure.size(); ++i) {
                const double time = targetPressure[i].x();
                if(time < overlapMinX || time > overlapMaxX) {
                    continue;
                }
                double pressure = 0.0;
                double derivative = 0.0;
                if(!interpolateLogValue(resultPressure, time, &pressure) ||
                        !interpolateLogValue(resultDerivative, time, &derivative)) {
                    return invalidLoss;
                }
                times.append(time);
                pressureResidual.append(pressure - qLn(qMax(targetPressure[i].y(), valueFloor)));
                derivativeResidual.append(derivative - qLn(targetDerivative[i].y()));
            }
            if(times.size() < 3) {
                return invalidLoss;
            }
            overlapMinX = times.first();
            overlapMaxX = times.last();
            QVector<double> coordinates;
            const double logMin = qLn(overlapMinX);
            const double logSpan = qLn(overlapMaxX) - logMin;
            for(int i = 0; i < times.size(); ++i) {
                coordinates.append((qLn(times[i]) - logMin) / logSpan);
            }

            // 数据较少时直接全量。补点只引用原始目标点，且各层共享同一组锚点。
            const int stride = times.size() <= 41 ? 1 : m_samplingStride;
            const QVector<int> indices = autoFitSamplingIndices(coordinates, stride);
            const QVector<double> fullWeights = autoFitLogTimeWeights(coordinates);
            QVector<double> layerCoordinates;
            for(int i = 0; i < indices.size(); ++i) {
                layerCoordinates.append(coordinates[indices[i]]);
            }
            const QVector<double> layerWeights = autoFitLogTimeWeights(layerCoordinates);
            AutoFitObjectiveBreakdownLM breakdown;
            breakdown.sampleCoordinates = layerCoordinates;
            breakdown.samplingStride = stride;
            breakdown.fullPointCount = times.size();
            double pressureEnergy = 0.0;
            double derivativeEnergy = 0.0;
            for(int i = 0; i < times.size(); ++i) {
                pressureEnergy += fullWeights[i] * pressureResidual[i] * pressureResidual[i];
                derivativeEnergy += fullWeights[i] * derivativeResidual[i] * derivativeResidual[i];
            }
            breakdown.pressureLoss = qSqrt(pressureEnergy);
            breakdown.derivativeLoss = qSqrt(derivativeEnergy);
            breakdown.total = qSqrt(0.5 * (pressureEnergy + derivativeEnergy));
            for(int component = 0; component < 2; ++component) {
                for(int i = 0; i < indices.size(); ++i) {
                    const double residual = component == 0
                            ? pressureResidual[indices[i]] : derivativeResidual[indices[i]];
                    breakdown.residualVector.append(qSqrt(0.5 * layerWeights[i]) * residual);
                }
            }
            breakdown.layerError = qSqrt(trustRegionSquaredNorm(breakdown.residualVector));
            breakdown.valid = isFiniteNumber(breakdown.total) && breakdown.total < 1.0e9;
            if(!populateStageMetrics(&breakdown)) return invalidLoss;
            if(!trustRegionResidualsValid(breakdown)) {
                return invalidLoss;
            }
            breakdown.timeWindows = calculateAutoFitTimeWindows(
                breakdown.residualVector, overlapMinX, overlapMaxX,
                layerCoordinates, layerWeights);
            m_lastObjectiveBreakdown = breakdown;
            // 首次完整评价有效后再冻结区间和实际层级，失败候选不能改变基准。
            m_samplingStride = stride;
            if(!comparisonRangeFixed) {
                m_comparisonTimeMin = overlapMinX;
                m_comparisonTimeMax = overlapMaxX;
                writeTraceMetaFile();
            }
            return breakdown.total;
        }

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

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

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

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

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

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

        AutoFitObjectiveBreakdownLM breakdown;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        // LM 总目标等于固定残差向量的二范数；压力和导数各占一半能量。
        // total 的定义不变；阶段形状和高度指标在下方以固定网格单独计算。
        breakdown.total = qSqrt(
            0.5 * breakdown.pressureLoss * breakdown.pressureLoss +
            0.5 * breakdown.derivativeLoss * breakdown.derivativeLoss);
        breakdown.layerError = breakdown.total;
        breakdown.valid =
            isFiniteNumber(breakdown.total) &&
            breakdown.total >= 0.0;
        if(breakdown.valid && breakdown.total < 1.0e9) {
            breakdown.timeWindows = calculateAutoFitTimeWindows(
                        breakdown.residualVector, overlapMinX, overlapMaxX);
            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，保证界面不会完全卡死。
        dllTask->start();

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

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

            // 等待期间也派发鼠标、键盘事件，让停止按钮能及时登记请求。
            QApplication::processEvents(QEventLoop::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;
            m_consecutiveFailures++;
            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;
            m_consecutiveFailures++;
            return result;
        }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    return result;
}

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

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

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

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

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

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

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

    }

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

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

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

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

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

        case LM_MAX_ITERATIONS:
            return tr("Total-stage 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");
    }
}