diff --git a/Bin/Config/Lang/cn/nmNum_cn.qm b/Bin/Config/Lang/cn/nmNum_cn.qm index 51dbc19b..409f506c 100644 Binary files a/Bin/Config/Lang/cn/nmNum_cn.qm and b/Bin/Config/Lang/cn/nmNum_cn.qm differ diff --git a/Bin/Config/Lang/cn/nmNum_cn.ts b/Bin/Config/Lang/cn/nmNum_cn.ts index eb5e59ba..92d6829e 100644 --- a/Bin/Config/Lang/cn/nmNum_cn.ts +++ b/Bin/Config/Lang/cn/nmNum_cn.ts @@ -664,6 +664,34 @@ Reason: %1 nmCalculationAutoFitLM + + Failed to refresh the sampling layer from the current curve. + 无法根据当前曲线刷新采样层级。 + + + LM sampling refined: %1 / %2 target points; full-target error: %3 + LM 采样加密：%1 / %2 个目标点；全目标点误差：%3 + + + LM sampling: layered target points (%1 / %2); acceptance uses all valid target points. + LM 采样方式：目标点分层采样（%1 / %2 点），使用全部有效目标点验收。 + + + LM sampling: fixed 80 points (original mode). + LM 采样方式：固定 80 点（原模式）。 + + + Full-target comparison error: initial=%1; final=%2 + 全目标点对比误差：初始=%1；最终=%2 + + + Unavailable + 不可用 + + + Budget exhausted before the full sampling layer; convergence is not confirmed. + 预算已耗尽，尚未进入完整采样层，未确认收敛。 + === User Stop Request Received === === 用户停止请求已接收 === @@ -3794,6 +3822,22 @@ Supported types: Vertical, Vertical Fractured, and Horizontal Multi-Fractured We nmWxAutomaticFitting + + LM sampling: + LM 采样方式： + + + Fixed sampling (80 points) + 固定采样（80 点） + + + Layered target sampling + 目标点分层采样 + + + Layered sampling uses target-point subsets for directions and all valid target points for acceptance. + 逐层增加目标曲线采样点，并始终使用全部有效目标点判断拟合是否改善。 + Automatic fitting 自动拟合 diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h index bea69df0..b2048714 100644 --- a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h +++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h @@ -12,7 +12,7 @@ #include "nmCalculation_global.h" // 固定对数时间窗口的误差诊断。时间边界为不含重叠区的基础边界； -// rmsError 是局部加权均方根，energy 是对 total 平方的贡献。 +// rmsError 是局部加权均方根，energy 是对当前层残差能量的贡献。 struct AutoFitTimeWindowLM { double timeMin; double timeMax; @@ -25,7 +25,7 @@ struct AutoFitTimeWindowLM { {} }; -// 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据， +// 双对数曲线误差分解。total 是候选接受和排序的统一依据；分层模式使用全部目标点， // 时间窗口和残差用于 Fisher 选参，其余诊断量用于解释曲线失配。 struct AutoFitObjectiveBreakdownLM { bool valid; @@ -33,6 +33,11 @@ struct AutoFitObjectiveBreakdownLM { double pressureLoss; double derivativeLoss; QVector residualVector; + // 分层模式的残差属于当前层，total 属于全部目标点；固定模式坐标数组为空。 + QVector sampleCoordinates; + int samplingStride; + int fullPointCount; + double layerError; QVector timeWindows; double verticalCommonBias; double verticalLoss; @@ -51,6 +56,9 @@ struct AutoFitObjectiveBreakdownLM { , total(1.0e10) , pressureLoss(std::numeric_limits::quiet_NaN()) , derivativeLoss(std::numeric_limits::quiet_NaN()) + , samplingStride(1) + , fullPointCount(80) + , layerError(1.0e10) , verticalCommonBias(std::numeric_limits::quiet_NaN()) , verticalLoss(std::numeric_limits::quiet_NaN()) , verticalReliable(false) @@ -97,6 +105,7 @@ public: int getTotalEvaluations() const; void resetOptimizer(); void setTargetWellName(const QString& wellName); + void setLayeredSamplingEnabled(bool enabled); signals: void progressUpdated(int iteration, double bestFitness); @@ -191,6 +200,9 @@ private: double m_comparisonTimeMin; double m_comparisonTimeMax; QString m_targetWellName; + // 此开关由拟合界面传入；每轮重置层级，固定模式仍使用原来的 80 点目标。 + bool m_layeredSampling; + int m_samplingStride; int m_maxIterations; double m_targetError; diff --git a/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h b/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h index 3f526c4d..4c6c6947 100644 --- a/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h +++ b/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h @@ -82,6 +82,8 @@ private: QLineEdit* m_errorLimitEdit; QComboBox* m_targetWellCombo; QComboBox* m_algorithmCombo; + QLabel* m_samplingLabel; + QComboBox* m_samplingCombo; QLabel* m_surrogateLabel; QComboBox* m_surrogateCombo; diff --git a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp index 8e4ae45b..fece2678 100644 --- a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp +++ b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp @@ -53,10 +53,52 @@ static double autoFitTimeWindowWeight(double coordinate, int windowIndex) return weight; } +// 目标时间点构成嵌套子集。窗口中心、重叠区两端及边界使用原始点两侧作锚点。 +static QVector autoFitSamplingIndices(const QVector& coordinates, int stride) +{ + QVector 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(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 indices; + for(int i = 0; i < selected.size(); ++i) { + if(selected[i]) indices.append(i); + } + return indices; +} + +// 在归一化 log(time) 上使用梯形积分权重，避免原始数据密集段被重复放大。 +static QVector autoFitLogTimeWeights(const QVector& coordinates) +{ + QVector 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 calculateAutoFitTimeWindows( - const QVector& residualVector, double timeMin, double timeMax) + const QVector& residualVector, double timeMin, double timeMax, + const QVector& coordinates = QVector(), + const QVector& timeWeights = QVector()) { - // 残差前后两半分别为压力和导数，已包含各占一半及采样点数的归一化。 + // 残差前后两半分别为压力和导数，已包含各占一半及时间采样权重的归一化。 const int pointCount = residualVector.size() / 2; const double logMin = qLn(timeMin); const double logSpan = qLn(timeMax) - logMin; @@ -69,13 +111,16 @@ static QVector calculateAutoFitTimeWindows( : qExp(logMin + logSpan * (k + 1) / kAutoFitTimeWindowCount); for(int i = 0; i < pointCount; ++i) { const double weight = autoFitTimeWindowWeight( - static_cast(i) / (pointCount - 1), k); - window.weightSum += weight; + coordinates.isEmpty() ? static_cast(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 = qSqrt(window.energy * pointCount / window.weightSum); + window.rmsError = window.weightSum > 0.0 + ? qSqrt(window.energy * (timeWeights.isEmpty() ? pointCount : 1.0) + / window.weightSum) : 0.0; } return windows; } @@ -262,9 +307,11 @@ struct TrustRegionEvaluation static bool trustRegionResidualsValid( const AutoFitObjectiveBreakdownLM& breakdown) { - // 损失函数固定使用 80 个压力点和 80 个导数点。严格校验长度，避免 - // Jacobian 沿用旧维度后访问另一候选的短残差向量。 - if(!breakdown.valid || breakdown.residualVector.size() != 160) { + // 固定模式仍要求 160 维；分层模式按当前时间点数校验，切层后重建 J。 + const int pointCount = breakdown.sampleCoordinates.isEmpty() + ? 80 : breakdown.sampleCoordinates.size(); + if(!breakdown.valid || pointCount < 3 || + breakdown.residualVector.size() != 2 * pointCount) { return false; } @@ -429,9 +476,10 @@ struct TrustRegionFisher static QVector buildTrustRegionFisher( const QVector >& jacobian, const QVector& residual, - const QVector& columnValid) + const QVector& columnValid, + const QVector& coordinates = QVector()) { - // 残差和 J 已包含压力/导数及点数归一化，只再乘一次窗口权重。 + // 残差和 J 已包含压力/导数及时间采样权重，只再乘一次窗口权重。 // 前半段是压力，后半段是导数；同一时间点的两行使用相同权重。 const int dimensions = columnValid.size(); const int pointCount = residual.size() / 2; @@ -441,7 +489,8 @@ static QVector buildTrustRegionFisher( TrustRegionFisher& local = information[k]; for(int row = 0; row < residual.size(); ++row) { const double weight = autoFitTimeWindowWeight( - static_cast(row % pointCount) / (pointCount - 1), k); + coordinates.isEmpty() ? static_cast(row % pointCount) / (pointCount - 1) + : coordinates[row % pointCount], k); for(int p = 0; p < dimensions; ++p) { if(!columnValid[p]) { continue; @@ -685,6 +734,8 @@ nmCalculationAutoFitLM::nmCalculationAutoFitLM(QObject* parent) , 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) @@ -864,6 +915,7 @@ void nmCalculationAutoFitLM::resetOptimizer() 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; @@ -878,6 +930,14 @@ void nmCalculationAutoFitLM::resetOptimizer() DEBUG_OUT("LM optimizer reset"); } +void nmCalculationAutoFitLM::setLayeredSamplingEnabled(bool enabled) +{ + // 拟合运行期间不允许改变残差定义；新一轮由 resetOptimizer 初始化层级。 + if(!m_isRunning) { + m_layeredSampling = enabled; + } +} + void nmCalculationAutoFitLM::setTargetWellName(const QString& wellName) { // 目标井名是贯穿拟合流程的关键索引： @@ -982,6 +1042,8 @@ void nmCalculationAutoFitLM::writeTraceHeader() << prefix + "weight_sum" << prefix + "rms_error" << prefix + "energy"; } + cols << "sampling_mode" << "sampling_stride" << "sampling_points" + << "full_target_points" << "layer_objective"; QTextStream out(&m_traceFile); out << cols.join(",") << "\n"; } @@ -1018,12 +1080,14 @@ void nmCalculationAutoFitLM::writeTraceMetaFile() QTextStream out(&metaFile); out << "{\n"; - out << " \"schema_version\": 3,\n"; + out << " \"schema_version\": 4,\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"; @@ -1128,6 +1192,11 @@ void nmCalculationAutoFitLM::writeTraceRow( } } + 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()); QTextStream out(&m_traceFile); out << cols.join(",") << "\n"; m_traceFile.flush(); @@ -1428,6 +1497,29 @@ bool nmCalculationAutoFitLM::startAutoFitting() 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, @@ -1441,8 +1533,8 @@ bool nmCalculationAutoFitLM::startAutoFitting() &m_globalBestObjectiveBreakdown); } - if(finalReason != LM_USER_STOPPED && - m_globalBestFitness < m_targetError) { + 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) { @@ -1714,7 +1806,7 @@ bool nmCalculationAutoFitLM::evaluateTrustRegionPoint( } // evaluateFitness() 会写入 DataManager 并调用真实求解器。这里统一统计 - // 真实评价次数和耗时，同时严格要求固定残差、诊断结构和结果曲线均有效。 + // 真实评价次数和耗时，同时要求当前层残差、误差结构和结果曲线均有效。 QTime timer; timer.start(); *fitness = evaluateFitness(parameters); @@ -1776,7 +1868,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() bool globalFallbackAttempted = false; StopReasonLM stopReason = LM_MAX_ITERATIONS; - // jacobian 的行对应固定 160 维残差，列对应用户勾选的参数。 + // jacobian 的行对应当前采样层的残差，列对应用户勾选的参数。 // Fisher 直接复用残差 Jacobian，上下/左右/形状诊断仅保留用于结果说明。 QVector > jacobian; QVector jacobianColumnValid(dimensions, false); @@ -1887,6 +1979,55 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() // 有效改善始终相对“上一次有效改善后的误差”累计判断，避免一连串微小 // 下降每次都清零计数；累计达到门槛后才开始新的有效改善基准。 double effectiveImprovementBaseline = current.fitness; + 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; + trustRadius = 0.12; + damping = 1.0e-2; + consecutiveRejectedSteps = 0; + consecutiveIneffectiveSteps = 0; + acceptedSinceRebuild = 0; + movementSinceRebuild = 0.0; + modelRebuiltAtMinimumRadius = false; + stagnationConfirmationRequested = false; + attemptedWindows.fill(false); + globalFallbackAttempted = false; + fullDataRejections = 0; + effectiveImprovementBaseline = current.fitness; + 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", ¤t.breakdown); + return true; + }; + + emit logMessageGenerated(m_layeredSampling + ? tr("LM sampling: layered target points (%1 / %2); 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 registerEffectiveImprovement = [&](double fitness) -> bool { const double requiredImprovement = qMax( effectiveAbsoluteImprovement, @@ -1941,7 +2082,8 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() .arg(maximumIneffectiveSteps)); if(current.fitness < m_targetError) { - return LM_TARGET_ACHIEVED; + promoteSampling(true); + return samplingRefreshFailed ? LM_OPTIMIZATION_FAILED : LM_TARGET_ACHIEVED; } // 在同一个真实工作点逐参数做单边差分。首选可用空间更大的方向；只有该方向 @@ -2122,16 +2264,26 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() if(!processPauseAndStop()) { break; } + if(m_layeredSampling && m_samplingStride > 1) { + const bool reserveFinalBudget = maximumEvaluations - m_totalEvaluations <= 2 * (dimensions + 1); + const int layerDeadline = qMax(1, m_maxIterations * (m_samplingStride == 4 ? 1 : 2) / 3); + if(reserveFinalBudget || iteration >= layerDeadline) { + promoteSampling(reserveFinalBudget); + if(samplingRefreshFailed) break; + } + } if(rebuildRequested) { const bool confirmingStagnation = stagnationConfirmationRequested; if(!rebuildSensitivity()) { + if(promoteSampling(false)) continue; stopReason = m_shouldStop ? LM_USER_STOPPED : LM_LOCAL_OPTIMUM; break; } if(current.fitness < m_targetError) { + promoteSampling(true); stopReason = LM_TARGET_ACHIEVED; break; } @@ -2142,6 +2294,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() const bool rebuildEffective = registerEffectiveImprovement(current.fitness); if(confirmingStagnation && !rebuildEffective) { + if(promoteSampling(false)) continue; emit logMessageGenerated( tr("Sensitivity rebuild produced no effective improvement; " "local convergence detected")); @@ -2152,7 +2305,8 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() // 每轮从最新 J 和当前残差重算窗口 Fisher，包含有限差分与割线更新的变化。 const QVector information = buildTrustRegionFisher( - jacobian, current.breakdown.residualVector, jacobianColumnValid); + jacobian, current.breakdown.residualVector, jacobianColumnValid, + current.breakdown.sampleCoordinates); const TrustRegionFisher& global = information[kAutoFitTimeWindowCount]; const int dominantComponent = trustRegionDominantComponent( current.breakdown, diagnosisThreshold); // 仅用于现有诊断日志。 @@ -2211,6 +2365,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() if(selectedColumns.isEmpty()) { if(trustRadius <= minimumTrustRadius * 1.01 && modelRebuiltAtMinimumRadius) { + if(promoteSampling(false)) continue; stopReason = LM_LOCAL_OPTIMUM; break; } @@ -2218,6 +2373,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() damping = qMin(1.0e8, damping * 4.0); rebuildRequested = true; if(recordIneffectiveStep()) { + if(promoteSampling(false)) continue; stopReason = LM_LOCAL_OPTIMUM; break; } @@ -2273,6 +2429,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() break; } if(recordIneffectiveStep()) { + if(promoteSampling(false)) continue; stopReason = LM_LOCAL_OPTIMUM; break; } @@ -2290,11 +2447,18 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() coordinateStep); // reductionRatio 衡量局部线性模型的可信度：接近 1 表示预测准确； // 值较小表示虽然可能下降，但模型低估了非线性，需要收紧下一步。 - double actualReduction = 0.5 * - (current.fitness * current.fitness - - candidate.fitness * candidate.fitness); + double actualReduction = m_layeredSampling + ? 0.5 * (trustRegionSquaredNorm(current.breakdown.residualVector) - + trustRegionSquaredNorm(candidate.breakdown.residualVector)) + : 0.5 * (current.fitness * current.fitness - candidate.fitness * candidate.fitness); double reductionRatio = actualReduction / predictedReduction; bool accepted = candidate.fitness < current.fitness; + // 粗层认为下降而完整数据不认可时累计，连续两次就提前加密。 + if(m_layeredSampling && m_samplingStride > 1 && !accepted && actualReduction > 0.0) { + ++fullDataRejections; + } else { + fullDataRejections = 0; + } QString componentName = trustRegionComponentName(dominantComponent); if(accepted) { @@ -2378,17 +2542,23 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() .arg(accepted ? tr("accepted") : tr("rejected"))); emit progressUpdated(iteration + 1, m_globalBestFitness); - if(stopReason == LM_LOCAL_OPTIMUM) { - break; + if(stopReason == LM_LOCAL_OPTIMUM || fullDataRejections >= 2) { + if(promoteSampling(false)) { + stopReason = LM_MAX_ITERATIONS; + continue; + } + if(stopReason == LM_LOCAL_OPTIMUM || samplingRefreshFailed) break; } if(current.fitness < m_targetError) { + promoteSampling(true); stopReason = LM_TARGET_ACHIEVED; break; } if(selectedWindow < 0 && trustRadius <= minimumTrustRadius * 1.01 && consecutiveRejectedSteps >= 2) { if(modelRebuiltAtMinimumRadius) { + if(promoteSampling(false)) continue; stopReason = LM_LOCAL_OPTIMUM; break; } @@ -2404,8 +2574,10 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() if(m_shouldStop) { return LM_USER_STOPPED; } + if(samplingRefreshFailed) return LM_OPTIMIZATION_FAILED; if(current.fitness < m_targetError) { - return LM_TARGET_ACHIEVED; + promoteSampling(true); + return samplingRefreshFailed ? LM_OPTIMIZATION_FAILED : LM_TARGET_ACHIEVED; } if(stopReason == LM_CONSECUTIVE_FAILURES || stopReason == LM_LOCAL_OPTIMUM || @@ -3221,9 +3393,9 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError( // 仅首次有效评价使用交集建立基准；失败试算不能冻结区间，后续候选 // 必须覆盖完整基准，不允许靠丢失首尾点缩小误差或改变窗口位置。 const bool comparisonRangeFixed = m_comparisonTimeMin > 0.0; - const double overlapMinX = comparisonRangeFixed + double overlapMinX = comparisonRangeFixed ? m_comparisonTimeMin : qMax(targetMinX, resultMinX); - const double overlapMaxX = comparisonRangeFixed + double overlapMaxX = comparisonRangeFixed ? m_comparisonTimeMax : qMin(targetMaxX, resultMaxX); if(overlapMinX >= overlapMaxX) { return invalidLoss; @@ -3234,6 +3406,86 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError( return invalidLoss; } + if(m_layeredSampling) { + // 完整基准只取公共范围内的目标原始时间点；模拟点数不改变评价标准。 + QVector times; + QVector pressureResidual; + QVector 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 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 indices = autoFitSamplingIndices(coordinates, stride); + const QVector fullWeights = autoFitLogTimeWeights(coordinates); + QVector layerCoordinates; + for(int i = 0; i < indices.size(); ++i) { + layerCoordinates.append(coordinates[indices[i]]); + } + const QVector 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(!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 commonX(numPoints); QVector commonLogX(numPoints); QVector targetLogPressure(numPoints); @@ -3856,6 +4108,7 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError( 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; diff --git a/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp b/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp index 24200afd..1f0e2f3c 100644 --- a/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp +++ b/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp @@ -724,6 +724,13 @@ void nmWxAutomaticFitting::setupControlPanel() m_surrogateCombo->setCurrentIndex(automaticFittingData.getSurrogateScreeningEnabled() ? 1 : 0); m_surrogateCombo->setMaximumWidth(160); m_surrogateCombo->setMinimumWidth(160); + // 分层采样用于 LM 对比试验，每次打开默认使用原来的固定 80 点模式。 + m_samplingLabel = new QLabel(tr("LM sampling:")); + m_samplingCombo = new QComboBox(); + m_samplingCombo->addItem(tr("Fixed sampling (80 points)")); + m_samplingCombo->addItem(tr("Layered target sampling")); + m_samplingCombo->setMinimumWidth(180); + m_samplingCombo->setToolTip(tr("Layered sampling uses target-point subsets for directions and all valid target points for acceptance.")); connect(m_algorithmCombo, SIGNAL(currentIndexChanged(int)), this, SLOT(onAlgorithmChanged(int))); onAlgorithmChanged(m_algorithmCombo->currentIndex()); @@ -781,6 +788,10 @@ void nmWxAutomaticFitting::setupControlPanel() m_controlLayout->addWidget(m_surrogateCombo); m_controlLayout->addSpacing(15); + m_controlLayout->addWidget(m_samplingLabel); + m_controlLayout->addWidget(m_samplingCombo); + m_controlLayout->addSpacing(15); + m_controlLayout->addWidget(iterationLabel); m_controlLayout->addWidget(m_iterationEdit); m_controlLayout->addSpacing(15); @@ -802,6 +813,8 @@ void nmWxAutomaticFitting::onAlgorithmChanged(int index) const bool usePSO = (index == 0); m_surrogateLabel->setEnabled(usePSO); m_surrogateCombo->setEnabled(usePSO); + m_samplingLabel->setEnabled(!usePSO); + m_samplingCombo->setEnabled(!usePSO); } void nmWxAutomaticFitting::setupButtons() @@ -1189,6 +1202,7 @@ void nmWxAutomaticFitting::startAutoFitting(const QVector>& targ m_autoFitterLM = new nmCalculationAutoFitLM(this); m_autoFitterLM->setTargetLogLogData(targetData); m_autoFitterLM->setTargetWellName(targetWellName); + m_autoFitterLM->setLayeredSamplingEnabled(m_samplingCombo->currentIndex() == 1); } else { DEBUG_UI("Creating PSO auto fitter"); m_autoFitterPSO = new nmCalculationAutoFitPSO(this);