|
|
|
|
@ -53,10 +53,52 @@ static double autoFitTimeWindowWeight(double coordinate, int windowIndex)
|
|
|
|
|
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>& 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;
|
|
|
|
|
@ -69,13 +111,16 @@ static QVector<AutoFitTimeWindowLM> calculateAutoFitTimeWindows(
|
|
|
|
|
: qExp(logMin + logSpan * (k + 1) / kAutoFitTimeWindowCount);
|
|
|
|
|
for(int i = 0; i < pointCount; ++i) {
|
|
|
|
|
const double weight = autoFitTimeWindowWeight(
|
|
|
|
|
static_cast<double>(i) / (pointCount - 1), k);
|
|
|
|
|
window.weightSum += weight;
|
|
|
|
|
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 = 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<TrustRegionFisher> buildTrustRegionFisher(
|
|
|
|
|
const QVector<QVector<double> >& jacobian,
|
|
|
|
|
const QVector<double>& residual,
|
|
|
|
|
const QVector<bool>& columnValid)
|
|
|
|
|
const QVector<bool>& columnValid,
|
|
|
|
|
const QVector<double>& coordinates = QVector<double>())
|
|
|
|
|
{
|
|
|
|
|
// 残差和 J 已包含压力/导数及点数归一化,只再乘一次窗口权重。
|
|
|
|
|
// 残差和 J 已包含压力/导数及时间采样权重,只再乘一次窗口权重。
|
|
|
|
|
// 前半段是压力,后半段是导数;同一时间点的两行使用相同权重。
|
|
|
|
|
const int dimensions = columnValid.size();
|
|
|
|
|
const int pointCount = residual.size() / 2;
|
|
|
|
|
@ -441,7 +489,8 @@ static QVector<TrustRegionFisher> buildTrustRegionFisher(
|
|
|
|
|
TrustRegionFisher& local = information[k];
|
|
|
|
|
for(int row = 0; row < residual.size(); ++row) {
|
|
|
|
|
const double weight = autoFitTimeWindowWeight(
|
|
|
|
|
static_cast<double>(row % pointCount) / (pointCount - 1), k);
|
|
|
|
|
coordinates.isEmpty() ? static_cast<double>(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<QVector<double> > jacobian;
|
|
|
|
|
QVector<bool> 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<TrustRegionFisher> 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<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(!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);
|
|
|
|
|
@ -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;
|
|
|
|
|
|