YXS
/
nmWTAI-Platform
8
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

自动拟合损失函数改用普通对数 RMSE

Browse Source
feature/MultiWellAutoFit-20260805
lvjunjie 6 days ago
parent e8563506b9
commit 50606d3564
2 changed files with 42 additions and 93 deletions
Show Stats Download Patch File Download Diff File

6
Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h
Unescape Escape View File

@ -30,11 +30,11 @@ struct AutoFitObjectiveBreakdown {
double total; double total;
double pressureLoss; double pressureLoss;
double derivativeLoss; double derivativeLoss;
// 固定目标网格上的 Huber 等效残差。非代理搜索使用它建立完整 Jacobian， // 固定目标网格上的普通对数残差。非代理搜索使用它建立完整 Jacobian，
// 向量平方和与 total 的平方一致。 // 向量平方和与 total 的平方一致。
QVector<double> residualVector; QVector<double> residualVector;
// 上下偏差使用压力和导数残差共享的 Huber 稳健中心。 // 上下偏差使用压力和导数残差共享的算术平均中心。
// verticalCommonBias 为正表示模拟曲线整体偏高，为负表示整体偏低； // verticalCommonBias 为正表示模拟曲线整体偏高，为负表示整体偏低；
// verticalReliable=false 时仍保留数值，但不能据此确定参数调整方向。 // verticalReliable=false 时仍保留数值，但不能据此确定参数调整方向。
double verticalCommonBias; double verticalCommonBias;
@ -50,7 +50,7 @@ struct AutoFitObjectiveBreakdown {
// 诊断量选参，但去除公共中心后的 shapeLoss 仍可用于局部选参。 // 诊断量选参，但去除公共中心后的 shapeLoss 仍可用于局部选参。
bool registrationAmbiguous; bool registrationAmbiguous;
// 去除稳健公共中心和可信左右偏差后剩余的整体形状误差；verticalReliable // 去除公共均值中心和可信左右偏差后剩余的整体形状误差；verticalReliable
// 只控制能否把公共中心解释为上下参数方向，不改变 shape 的中心化公式。 // 只控制能否把公共中心解释为上下参数方向，不改变 shape 的中心化公式。
double shapeLoss; double shapeLoss;

129
Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp
Unescape Escape View File

@ -104,8 +104,8 @@ static inline bool isFiniteNumber(double value)
#endif #endif
} }
// 两个 Huber RMS 的差不能直接解释为被消除的独立误差。RMS 的平方才对应 // 两个 RMSE 的差不能直接解释为被消除的独立误差。RMSE 的平方才对应
// 稳健能量，因此先计算 reduced^2-full^2，再开方恢复原量纲。这里用于分别 // 均方能量，因此先计算 reduced^2-full^2，再开方恢复原量纲。这里用于分别
// 提取“消除公共上下偏差”和“消除水平位移”实际减少的误差贡献。 // 提取“消除公共上下偏差”和“消除水平位移”实际减少的误差贡献。
static double nestedRmsContribution(double reducedModelLoss, static double nestedRmsContribution(double reducedModelLoss,
double fullModelLoss) double fullModelLoss)
@ -3900,7 +3900,7 @@ struct TrustRegionEvaluation
{} {}
}; };
// LM 只使用固定长度、全部有限的稳健残差。代理路径不会进入本搜索器。 // LM 只使用固定长度、全部有限的普通残差。代理路径不会进入本搜索器。
static bool trustRegionResidualsValid( static bool trustRegionResidualsValid(
const AutoFitObjectiveBreakdown& breakdown) const AutoFitObjectiveBreakdown& breakdown)
{ {
@ -4216,7 +4216,7 @@ StopReasonPSO nmCalculationAutoFitPSO::runTrustRegionFitting()
bool modelRebuiltAtMinimumRadius = false; bool modelRebuiltAtMinimumRadius = false;
StopReasonPSO stopReason = PSO_MAX_ITERATIONS; StopReasonPSO stopReason = PSO_MAX_ITERATIONS;
// jacobian 的行对应固定 100 维稳健残差，列对应用户勾选的参数。 // jacobian 的行对应固定 100 维残差，列对应用户勾选的参数。
// 三个 gradient 单独描述诊断分量对参数的局部变化，只用于本轮选参。 // 三个 gradient 单独描述诊断分量对参数的局部变化，只用于本轮选参。
QVector<QVector<double> > jacobian; QVector<QVector<double> > jacobian;
QVector<double> verticalGradient(dimensions, 0.0); QVector<double> verticalGradient(dimensions, 0.0);
@ -4803,7 +4803,7 @@ StopReasonPSO nmCalculationAutoFitPSO::runTrustRegionFitting()
} }
stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep)); stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
// 用线性模型 r_new ~= r_current + J*step 预测稳健残差，再用平方能量 // 用线性模型 r_new ~= r_current + J*step 预测残差，再用平方能量
// 的下降量与真实候选下降量比较，作为调整阻尼和半径的依据。 // 的下降量与真实候选下降量比较，作为调整阻尼和半径的依据。
QVector<double> predictedResidual = QVector<double> predictedResidual =
current.breakdown.residualVector; current.breakdown.residualVector;
@ -6462,9 +6462,6 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
// 位于 log(time)-log(value) 坐标，因此得到的是相对尺度偏差而非原始压力量纲。 // 位于 log(time)-log(value) 坐标，因此得到的是相对尺度偏差而非原始压力量纲。
const double invalidLoss = 1.0e10; const double invalidLoss = 1.0e10;
const double valueFloor = 1.0e-12; const double valueFloor = 1.0e-12;
// Huber 转折点 qLn(1.2) 对应约 20% 的倍率偏差；小偏差保持平方惩罚，
// 更大的局部尖峰改为近似线性惩罚，避免单点支配整条曲线。
const double huberDelta = qLn(1.2);
const double minimumCoverage = 0.95; const double minimumCoverage = 0.95;
const int numPoints = 50; const int numPoints = 50;
m_lastObjectiveBreakdown = AutoFitObjectiveBreakdown(); m_lastObjectiveBreakdown = AutoFitObjectiveBreakdown();
@ -6704,7 +6701,7 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
} }
// 通过门槛后最多只缺少首尾少量目标点。按模拟曲线端点趋势补齐后， // 通过门槛后最多只缺少首尾少量目标点。按模拟曲线端点趋势补齐后，
// 每个候选仍在固定 50 点上计算 Huber 均值，不能靠少算难拟合端点获益。 // 每个候选仍在固定 50 点上计算均方根误差，不能靠少算难拟合端点获益。
for(int i = 0; i < numPoints; ++i) { for(int i = 0; i < numPoints; ++i) {
if(isFiniteNumber(pressureResidual[i]) && if(isFiniteNumber(pressureResidual[i]) &&
isFiniteNumber(derivativeResidual[i])) { isFiniteNumber(derivativeResidual[i])) {
@ -6726,8 +6723,8 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
resultLogDerivative - targetLogDerivative[i]; resultLogDerivative - targetLogDerivative[i];
} }
// Huber RMS 在小残差处保持平方损失，在异常点处转为线性增长。 // 在指定中心附近计算普通均方根误差。
auto huberRmsAround = [huberDelta]( auto rmseAround = [](
const QVector<double>& values, const QVector<double>& values,
int begin, int begin,
int end, int end,
@ -6744,12 +6741,7 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
} }
double difference = values[i] - center; double difference = values[i] - center;
double absoluteValue = qAbs(difference); sum += difference * difference;
double rho = absoluteValue <= huberDelta
? difference * difference
: 2.0 * huberDelta * absoluteValue -
huberDelta * huberDelta;
sum += rho;
++count; ++count;
} }
@ -6758,27 +6750,15 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
: std::numeric_limits<double>::quiet_NaN(); : std::numeric_limits<double>::quiet_NaN();
}; };
auto huberRms = [&huberRmsAround]( auto rmse = [&rmseAround](
const QVector<double>& values, const QVector<double>& values,
int begin, int begin,
int end) -> double { int end) -> double {
return huberRmsAround(values, begin, end, 0.0); return rmseAround(values, begin, end, 0.0);
}; };
// 将 Huber 能量转换成带符号的等效残差，使向量平方和与稳健损失一致。 // 普通算术平均中心保留上下偏差的符号。
// 非代理 LM 直接对该固定长度向量建立 Jacobian。 auto meanCenterRange = [](
auto huberEquivalentResidual = [huberDelta](double value) -> double {
double absoluteValue = qAbs(value);
double rho = absoluteValue <= huberDelta
? value * value
: 2.0 * huberDelta * absoluteValue -
huberDelta * huberDelta;
double magnitude = qSqrt(qMax(0.0, rho));
return value < 0.0 ? -magnitude : magnitude;
};
// Huber 加权中心保留上下偏差的符号，并降低局部尖峰的影响。
auto huberCenterRange = [huberDelta](
const QVector<double>& values, const QVector<double>& values,
int begin, int begin,
int end) -> double { int end) -> double {
@ -6797,42 +6777,12 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
if(count == 0) { if(count == 0) {
return std::numeric_limits<double>::quiet_NaN(); return std::numeric_limits<double>::quiet_NaN();
} }
center /= count; return center / count;
for(int iteration = 0; iteration < 8; ++iteration) {
double weightedSum = 0.0;
double weightTotal = 0.0;
for(int i = validBegin; i < validEnd; ++i) {
if(!isFiniteNumber(values[i])) {
continue;
}
double distance = qAbs(values[i] - center);
double weight =
distance <= huberDelta ||
distance < 1.0e-12
? 1.0
: huberDelta / distance;
weightedSum += weight * values[i];
weightTotal += weight;
}
if(weightTotal <= 1.0e-12) {
break;
}
double nextCenter = weightedSum / weightTotal;
if(qAbs(nextCenter - center) <= 1.0e-12) {
center = nextCenter;
break;
}
center = nextCenter;
}
return center;
}; };
// 压力和导数合并后只求一个公共中心，表示两条曲线共同的上下位移。 // 压力和导数合并后只求一个公共中心，表示两条曲线共同的上下位移。
// 分别去中心会把压力与导数之间真实的相对形状差异一并消除。 // 分别去中心会把压力与导数之间真实的相对形状差异一并消除。
auto huberCommonCenterRange = [&huberCenterRange]( auto commonMeanCenterRange = [&meanCenterRange](
const QVector<double>& pressureValues, const QVector<double>& pressureValues,
const QVector<double>& derivativeValues, const QVector<double>& derivativeValues,
int begin, int begin,
@ -6853,19 +6803,19 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
combined.append(derivativeValues[i]); combined.append(derivativeValues[i]);
} }
} }
return huberCenterRange(combined, 0, combined.size()); return meanCenterRange(combined, 0, combined.size());
}; };
// 两个通道按能量等权合并，返回值与单通道 Huber RMS 保持同一量纲。 // 两个通道按能量等权合并，返回值与单通道 RMSE 保持同一量纲。
auto jointHuberRmsAround = [&huberRmsAround]( auto jointRmseAround = [&rmseAround](
const QVector<double>& pressureValues, const QVector<double>& pressureValues,
const QVector<double>& derivativeValues, const QVector<double>& derivativeValues,
int begin, int begin,
int end, int end,
double center) -> double { double center) -> double {
double pressureLoss = huberRmsAround( double pressureLoss = rmseAround(
pressureValues, begin, end, center); pressureValues, begin, end, center);
double derivativeLoss = huberRmsAround( double derivativeLoss = rmseAround(
derivativeValues, begin, end, center); derivativeValues, begin, end, center);
if(!isFiniteNumber(pressureLoss) || if(!isFiniteNumber(pressureLoss) ||
!isFiniteNumber(derivativeLoss)) { !isFiniteNumber(derivativeLoss)) {
@ -6878,9 +6828,9 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
// 主目标始终使用未做上下或左右校正的完整曲线误差。 // 主目标始终使用未做上下或左右校正的完整曲线误差。
breakdown.pressureLoss = breakdown.pressureLoss =
huberRms(pressureResidual, 0, numPoints); rmse(pressureResidual, 0, numPoints);
breakdown.derivativeLoss = breakdown.derivativeLoss =
huberRms(derivativeResidual, 0, numPoints); rmse(derivativeResidual, 0, numPoints);
// 压力和导数各占一半权重。缩放后 residualVector 的二范数就是 // 压力和导数各占一半权重。缩放后 residualVector 的二范数就是
// sqrt(0.5 * pressureLoss^2 + 0.5 * derivativeLoss^2)。 // sqrt(0.5 * pressureLoss^2 + 0.5 * derivativeLoss^2)。
@ -6889,12 +6839,12 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
for(int i = 0; i < numPoints; ++i) { for(int i = 0; i < numPoints; ++i) {
breakdown.residualVector.append( breakdown.residualVector.append(
residualScale * residualScale *
huberEquivalentResidual(pressureResidual[i])); pressureResidual[i]);
} }
for(int i = 0; i < numPoints; ++i) { for(int i = 0; i < numPoints; ++i) {
breakdown.residualVector.append( breakdown.residualVector.append(
residualScale * residualScale *
huberEquivalentResidual(derivativeResidual[i])); derivativeResidual[i]);
} }
const double logGridStep = const double logGridStep =
@ -7049,31 +6999,31 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
continue; continue;
} }
double commonBias = huberCommonCenterRange( double commonBias = commonMeanCenterRange(
shiftedPressureResidual, shiftedPressureResidual,
shiftedDerivativeResidual, shiftedDerivativeResidual,
registrationBegin, registrationBegin,
registrationEnd); registrationEnd);
double centeredLoss = jointHuberRmsAround( double centeredLoss = jointRmseAround(
shiftedPressureResidual, shiftedPressureResidual,
shiftedDerivativeResidual, shiftedDerivativeResidual,
registrationBegin, registrationBegin,
registrationEnd, registrationEnd,
commonBias); commonBias);
double pressureBias = huberCenterRange( double pressureBias = meanCenterRange(
shiftedPressureResidual, shiftedPressureResidual,
registrationBegin, registrationBegin,
registrationEnd); registrationEnd);
double derivativeBias = huberCenterRange( double derivativeBias = meanCenterRange(
shiftedDerivativeResidual, shiftedDerivativeResidual,
registrationBegin, registrationBegin,
registrationEnd); registrationEnd);
double pressureLoss = huberRmsAround( double pressureLoss = rmseAround(
shiftedPressureResidual, shiftedPressureResidual,
registrationBegin, registrationBegin,
registrationEnd, registrationEnd,
pressureBias); pressureBias);
double derivativeLoss = huberRmsAround( double derivativeLoss = rmseAround(
shiftedDerivativeResidual, shiftedDerivativeResidual,
registrationBegin, registrationBegin,
registrationEnd, registrationEnd,
@ -7124,7 +7074,7 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
return invalidLoss; return invalidLoss;
} }
// horizontalGain 是“允许水平位移”相对“固定零位移”减少的稳健能量。 // horizontalGain 是“允许水平位移”相对“固定零位移”减少的均方能量。
// 只有改善足够明显且最优点不是边界，才把位移解释为可靠左右偏差。 // 只有改善足够明显且最优点不是边界，才把位移解释为可靠左右偏差。
double horizontalGain = nestedRmsContribution( double horizontalGain = nestedRmsContribution(
zeroShiftCenteredLoss, bestCenteredLoss); zeroShiftCenteredLoss, bestCenteredLoss);
@ -7183,8 +7133,7 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
isFiniteNumber(maximumNearOptimalBias) isFiniteNumber(maximumNearOptimalBias)
? maximumNearOptimalBias - minimumNearOptimalBias ? maximumNearOptimalBias - minimumNearOptimalBias
: std::numeric_limits<double>::infinity(); : std::numeric_limits<double>::infinity();
bool commonBiasStable = nearOptimalBiasSpread <= bool commonBiasStable = nearOptimalBiasSpread <= 1.0e-2;
qMax(1.0e-4, huberDelta * 0.05);
bool horizontalAtBoundary = bool horizontalAtBoundary =
halfShiftStepCount > 0 && halfShiftStepCount > 0 &&
qAbs(bestShiftStep) == halfShiftStepCount; qAbs(bestShiftStep) == halfShiftStepCount;
@ -7249,18 +7198,18 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
return invalidLoss; return invalidLoss;
} }
double commonBias = huberCommonCenterRange( double commonBias = commonMeanCenterRange(
shiftedPressureResidual, shiftedPressureResidual,
shiftedDerivativeResidual, shiftedDerivativeResidual,
diagnosticBegin, diagnosticBegin,
diagnosticEnd); diagnosticEnd);
double rawAlignedLoss = jointHuberRmsAround( double rawAlignedLoss = jointRmseAround(
shiftedPressureResidual, shiftedPressureResidual,
shiftedDerivativeResidual, shiftedDerivativeResidual,
diagnosticBegin, diagnosticBegin,
diagnosticEnd, diagnosticEnd,
0.0); 0.0);
double centeredAlignedLoss = jointHuberRmsAround( double centeredAlignedLoss = jointRmseAround(
shiftedPressureResidual, shiftedPressureResidual,
shiftedDerivativeResidual, shiftedDerivativeResidual,
diagnosticBegin, diagnosticBegin,
@ -7297,10 +7246,10 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
} }
} }
// shapeLoss 是去除可信左右位移和稳健公共中心后的剩余误差。 // shapeLoss 是去除可信左右位移和公共均值中心后的剩余误差。
double shapePressureLoss = huberRms( double shapePressureLoss = rmse(
shapePressure, diagnosticBegin, diagnosticEnd); shapePressure, diagnosticBegin, diagnosticEnd);
double shapeDerivativeLoss = huberRms( double shapeDerivativeLoss = rmse(
shapeDerivative, diagnosticBegin, diagnosticEnd); shapeDerivative, diagnosticBegin, diagnosticEnd);
if(!isFiniteNumber(shapePressureLoss) || if(!isFiniteNumber(shapePressureLoss) ||
!isFiniteNumber(shapeDerivativeLoss)) { !isFiniteNumber(shapeDerivativeLoss)) {
@ -7333,7 +7282,7 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
0.5 * breakdown.pressureLoss + 0.5 * breakdown.pressureLoss +
0.5 * breakdown.derivativeLoss; 0.5 * breakdown.derivativeLoss;
} else { } else {
// 非代理总目标等于固定 100 维 Huber 等效残差的二范数；压力和 // 非代理总目标等于固定 100 维普通残差的二范数；压力和
// 导数各占一半能量。上下、左右和形状分量不参与候选排序与接受。 // 导数各占一半能量。上下、左右和形状分量不参与候选排序与接受。
breakdown.total = qSqrt( breakdown.total = qSqrt(
0.5 * breakdown.pressureLoss * breakdown.pressureLoss + 0.5 * breakdown.pressureLoss * breakdown.pressureLoss +

Write Preview
Loading…
Cancel
Save