From 04a8759b648a9c05adee9a50da7158e1a9c7b92e Mon Sep 17 00:00:00 2001
From: lvjunjie <lvjunjie@abc.com>
Date: Thu, 15 Oct 2026 11:39:30 +0800
Subject: [PATCH] =?UTF-8?q?feat(nmNum):=20=E5=A2=9E=E5=8A=A0=20Monte-Carlo?=
 =?UTF-8?q?=20=E5=8E=8B=E5=8A=9B=E8=AF=AF=E5=B7=AE=E8=AE=A1=E7=AE=97?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- 使用求解快照中的实测压力，插值计算每组 RMSE。
- 失败、取消或无有效数据时标记误差无效。
- 结果暂存内存，供后续绘图和筛选使用。
---
 Include/nmNum/nmData/nmDataSensitive.h        |  3 +-
 .../nmCalculationDllPebiSolverTask.cpp        | 48 +++++++++++++++++++
 2 files changed, 50 insertions(+), 1 deletion(-)

diff --git a/Include/nmNum/nmData/nmDataSensitive.h b/Include/nmNum/nmData/nmDataSensitive.h
index 327133af..c63bc469 100644
--- a/Include/nmNum/nmData/nmDataSensitive.h
+++ b/Include/nmNum/nmData/nmDataSensitive.h
@@ -103,12 +103,13 @@ public:
     // 运行结果属于分析数据，求解器与界面共用同一结构，避免重复复制整批曲线。
     struct RunResult
     {
-        RunResult() : m_nRunIndex(0), m_bSucceeded(false), m_bCancelled(false), m_bVisible(true) {}
+        RunResult() : m_nRunIndex(0), m_dPressureRmse(-1.0), m_bSucceeded(false), m_bCancelled(false), m_bVisible(true) {}
         quint64 m_nRunIndex;
         QVector<double> m_vecParameters;
         QVector<QVector<double> > m_vecPressure;
         QVector<QVector<double> > m_vecLogLog;
         QVector<QVector<double> > m_vecSemiLog;
+        double m_dPressureRmse; // Monte-Carlo 压力均方根误差（MPa），-1 表示尚无有效误差。
         bool m_bSucceeded;
         bool m_bCancelled;
         bool m_bVisible;
diff --git a/Src/nmNum/nmCalculation/nmCalculationDllPebiSolverTask.cpp b/Src/nmNum/nmCalculation/nmCalculationDllPebiSolverTask.cpp
index ac9be4d4..7b63eb35 100644
--- a/Src/nmNum/nmCalculation/nmCalculationDllPebiSolverTask.cpp
+++ b/Src/nmNum/nmCalculation/nmCalculationDllPebiSolverTask.cpp
@@ -211,6 +211,46 @@ bool isFiniteSolverNumber(double value)
 #endif
 }
 
+double calculateSensitivityPressureRmse(const QVector<QVector<double> >& measured,
+                                        const QVector<QVector<double> >& simulated,
+                                        const QAtomicInt* cancelRequested)
+{
+    if(measured.size() < 2 || simulated.size() < 2 || simulated[0].isEmpty() ||
+       measured[0].size() != measured[1].size() ||
+       simulated[0].size() != simulated[1].size()) return -1.0;
+    const QVector<double>& times = simulated[0];
+    const QVector<double>& pressures = simulated[1];
+    for(int i = 0; i < times.size(); ++i) {
+        if((i % 256) == 0 && isTaskCancellationRequested(cancelRequested)) return -1.0;
+        if(!isFiniteSolverNumber(times[i]) || !isFiniteSolverNumber(pressures[i]) ||
+           (i > 0 && times[i] <= times[i - 1])) return -1.0;
+    }
+
+    // 在实测时间上线性插值原始压力，只统计覆盖范围内的有效点，不外推。
+    double squaredError = 0.0;
+    int count = 0;
+    for(int i = 0; i < measured[0].size(); ++i) {
+        if((i % 256) == 0 && isTaskCancellationRequested(cancelRequested)) return -1.0;
+        const double time = measured[0][i];
+        const double pressure = measured[1][i];
+        if(!isFiniteSolverNumber(time) || !isFiniteSolverNumber(pressure) ||
+           time < times.first() || time > times.last()) continue;
+        const int right = static_cast<int>(
+                    std::lower_bound(times.constBegin(), times.constEnd(), time) - times.constBegin());
+        double predicted = pressures[right];
+        if(times[right] != time) {
+            const double fraction = (time - times[right - 1]) / (times[right] - times[right - 1]);
+            predicted = pressures[right - 1] + fraction * (pressures[right] - pressures[right - 1]);
+        }
+        const double residual = predicted - pressure;
+        squaredError += residual * residual;
+        ++count;
+    }
+    // 普通 Monte-Carlo 的各点权重均为 1；无有效点不能当作完全拟合。
+    return count > 0 && isFiniteSolverNumber(squaredError)
+            ? std::sqrt(squaredError / count) : -1.0;
+}
+
 bool savedPseudoPressure(const nmDataBinaryTools::NM_PEBI_SCENE& scene,
                          double pressure, double& result)
 {
@@ -1012,6 +1052,13 @@ bool nmCalculationDllPebiSolverTask::executeSensitivityBatch()
                     result.m_vecPressure = m_autoFitResultPressure;
                     result.m_vecLogLog = m_autoFitResultLogLog;
                     result.m_vecSemiLog = m_autoFitResultSemiLog;
+                    if(m_pSensitivityResult->m_calcType == nmDataSensitive::CALC_MONTE_CARLO) {
+                        // 使用求解时冻结的目标井实测压力，不读取界面后来修改的数据。
+                        const int target = m_pInputSnapshot->m_oSensitivityBaseline.m_nTarget;
+                        result.m_dPressureRmse = calculateSensitivityPressureRmse(
+                                    m_pInputSnapshot->m_oGridInput.m_vecWellInputs[target].m_vecHistoryPressure,
+                                    result.m_vecPressure, &m_nCancelRequested);
+                    }
                 } else {
                     result.m_sFailureReason = m_sFailureReason;
                 }
@@ -1023,6 +1070,7 @@ bool nmCalculationDllPebiSolverTask::executeSensitivityBatch()
                 result.m_vecPressure.clear();
                 result.m_vecLogLog.clear();
                 result.m_vecSemiLog.clear();
+                result.m_dPressureRmse = -1.0;
             } else if(!result.m_bSucceeded) {
                 ++failed;
             }