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