|
|
#include "nmCalculationAutoFitLM.h"
|
|
|
#include "nmCalculationDllPebiSolverTask.h"
|
|
|
#include "nmCalculationUtils.h"
|
|
|
#include "nmDataAnalyzeManager.h"
|
|
|
#include "nmDataWellBase.h"
|
|
|
#include "nmDataVerticalFracturedWell.h"
|
|
|
#include "nmDataHorizontalFracturedWell.h"
|
|
|
#include "nmDataReservoir.h"
|
|
|
#include "nmDataAutomaticFitting.h"
|
|
|
|
|
|
#include <QApplication>
|
|
|
#include <QTime>
|
|
|
#include <QDir>
|
|
|
#include <QTextStream>
|
|
|
#include <QFileInfo>
|
|
|
#include <QDateTime>
|
|
|
#include <QtCore/qmath.h>
|
|
|
#include <cmath>
|
|
|
#include <algorithm>
|
|
|
#include <limits>
|
|
|
|
|
|
#ifdef Q_OS_WIN
|
|
|
#include <windows.h>
|
|
|
#include <float.h>
|
|
|
#define DEBUG_OUT(msg) OutputDebugStringA(QString("[AutoFitLM] %1\n").arg(msg).toLocal8Bit().data())
|
|
|
#endif
|
|
|
|
|
|
static const int kAutoFitIntervalsPerDecade = 20;
|
|
|
|
|
|
static int autoFitShapeLag(int pointCount)
|
|
|
{
|
|
|
// 保持斜率跨度约为完整对数时间范围的 10%，四舍五入到间隔数且至少跨一个间隔。
|
|
|
return qMax(1, qRound((pointCount - 1) * 0.10));
|
|
|
}
|
|
|
|
|
|
static inline bool isFiniteNumber(double value)
|
|
|
{
|
|
|
#ifdef Q_OS_WIN
|
|
|
return _finite(value) != 0;
|
|
|
#else
|
|
|
return std::isfinite(value);
|
|
|
#endif
|
|
|
}
|
|
|
|
|
|
// 从嵌套 RMSE 中提取被消除的独立误差贡献。
|
|
|
static double nestedRmsContribution(double reducedModelLoss,
|
|
|
double fullModelLoss)
|
|
|
{
|
|
|
return qSqrt(qMax(0.0,
|
|
|
reducedModelLoss * reducedModelLoss -
|
|
|
fullModelLoss * fullModelLoss));
|
|
|
}
|
|
|
|
|
|
static inline void msleep(int ms)
|
|
|
{
|
|
|
#ifdef Q_OS_WIN
|
|
|
Sleep(ms);
|
|
|
#else
|
|
|
Q_UNUSED(ms);
|
|
|
#endif
|
|
|
}
|
|
|
|
|
|
static QString csvEscape(const QString& text)
|
|
|
{
|
|
|
QString escaped = text;
|
|
|
escaped.replace("\"", "\"\"");
|
|
|
return QString("\"%1\"").arg(escaped);
|
|
|
}
|
|
|
|
|
|
static QString traceNumber(double value)
|
|
|
{
|
|
|
return isFiniteNumber(value) ? QString::number(value, 'g', 17) : QString();
|
|
|
}
|
|
|
|
|
|
static QString traceParamAt(const QVector<double>& params, int index)
|
|
|
{
|
|
|
return (index >= 0 && index < params.size())
|
|
|
? traceNumber(params[index])
|
|
|
: QString();
|
|
|
}
|
|
|
|
|
|
static QString jsonEscape(const QString& text)
|
|
|
{
|
|
|
QString escaped = text;
|
|
|
escaped.replace("\\", "\\\\");
|
|
|
escaped.replace("\"", "\\\"");
|
|
|
escaped.replace("\b", "\\b");
|
|
|
escaped.replace("\f", "\\f");
|
|
|
escaped.replace("\n", "\\n");
|
|
|
escaped.replace("\r", "\\r");
|
|
|
escaped.replace("\t", "\\t");
|
|
|
return QString("\"%1\"").arg(escaped);
|
|
|
}
|
|
|
|
|
|
static QString jsonNumber(double value)
|
|
|
{
|
|
|
return isFiniteNumber(value) ? QString::number(value, 'g', 17) : QString("null");
|
|
|
}
|
|
|
|
|
|
static QString jsonDoubleArray(const QVector<double>& values)
|
|
|
{
|
|
|
QStringList items;
|
|
|
for(int i = 0; i < values.size(); ++i) {
|
|
|
items << jsonNumber(values[i]);
|
|
|
}
|
|
|
return QString("[%1]").arg(items.join(","));
|
|
|
}
|
|
|
|
|
|
static QString jsonIntArray(const QVector<int>& values)
|
|
|
{
|
|
|
QStringList items;
|
|
|
for(int i = 0; i < values.size(); ++i) {
|
|
|
items << QString::number(values[i]);
|
|
|
}
|
|
|
return QString("[%1]").arg(items.join(","));
|
|
|
}
|
|
|
|
|
|
static QString jsonBoolArray(const QVector<bool>& values)
|
|
|
{
|
|
|
QStringList items;
|
|
|
for(int i = 0; i < values.size(); ++i) {
|
|
|
items << (values[i] ? "true" : "false");
|
|
|
}
|
|
|
return QString("[%1]").arg(items.join(","));
|
|
|
}
|
|
|
|
|
|
static QString jsonStringArray(const QStringList& values)
|
|
|
{
|
|
|
QStringList items;
|
|
|
for(int i = 0; i < values.size(); ++i) {
|
|
|
items << jsonEscape(values[i]);
|
|
|
}
|
|
|
return QString("[%1]").arg(items.join(","));
|
|
|
}
|
|
|
|
|
|
// 信赖域搜索统一在 [0, 1] 内部坐标工作。正值参数使用对数坐标，使内部相同步长
|
|
|
// 表示近似相同的相对变化，避免 k、C、Dfc 等跨数量级参数被线性尺度支配；
|
|
|
// 表皮可为负数，使用线性坐标。
|
|
|
static bool useTrustRegionLogScale(int parameterIndex, double lower, double upper)
|
|
|
{
|
|
|
return parameterIndex != 1 &&
|
|
|
lower > 0.0 && upper > lower;
|
|
|
}
|
|
|
|
|
|
static double toTrustRegionCoordinate(double value,
|
|
|
int parameterIndex,
|
|
|
double lower,
|
|
|
double upper)
|
|
|
{
|
|
|
// 所有进入优化器的物理值先投影到用户上下界，再转换成无量纲坐标。
|
|
|
// 有限差分步长和信赖半径因此不受参数物理单位的直接影响。
|
|
|
value = qMax(lower, qMin(upper, value));
|
|
|
|
|
|
if(useTrustRegionLogScale(parameterIndex, lower, upper)) {
|
|
|
return (qLn(value) - qLn(lower)) / (qLn(upper) - qLn(lower));
|
|
|
}
|
|
|
|
|
|
return upper > lower ? (value - lower) / (upper - lower) : 0.0;
|
|
|
}
|
|
|
|
|
|
static double fromTrustRegionCoordinate(double coordinate,
|
|
|
int parameterIndex,
|
|
|
double lower,
|
|
|
double upper)
|
|
|
{
|
|
|
coordinate = qMax(0.0, qMin(1.0, coordinate));
|
|
|
|
|
|
// 端点直接返回原始边界，避免对数或线性逆变换的舍入误差造成越界。
|
|
|
if(coordinate <= 0.0) {
|
|
|
return lower;
|
|
|
}
|
|
|
if(coordinate >= 1.0) {
|
|
|
return upper;
|
|
|
}
|
|
|
|
|
|
double value;
|
|
|
if(useTrustRegionLogScale(parameterIndex, lower, upper)) {
|
|
|
value = qExp(qLn(lower) + coordinate * (qLn(upper) - qLn(lower)));
|
|
|
} else {
|
|
|
value = lower + coordinate * (upper - lower);
|
|
|
}
|
|
|
|
|
|
// 内部点转换后也限制到实际上下界，保证候选参数通过严格的边界检查。
|
|
|
return qMax(lower, qMin(upper, value));
|
|
|
}
|
|
|
|
|
|
// 一次真实求解的完整快照。除了参数和总误差，还保存内部坐标、诊断分量和
|
|
|
// 双对数曲线，因此拒绝候选后可以完整恢复上一个已接受工作点。
|
|
|
struct TrustRegionEvaluation
|
|
|
{
|
|
|
QVector<double> parameters;
|
|
|
QVector<double> coordinates;
|
|
|
AutoFitObjectiveBreakdownLM breakdown;
|
|
|
QVector<QVector<double> > curve;
|
|
|
double fitness;
|
|
|
int elapsedMs;
|
|
|
bool valid;
|
|
|
|
|
|
TrustRegionEvaluation()
|
|
|
: fitness(1.0e10)
|
|
|
, elapsedMs(-1)
|
|
|
, valid(false)
|
|
|
{}
|
|
|
};
|
|
|
|
|
|
// LM 残差长度由本轮固定比较区间决定，所有残差必须有限。
|
|
|
static bool trustRegionResidualsValid(
|
|
|
const AutoFitObjectiveBreakdownLM& breakdown)
|
|
|
{
|
|
|
// 压力和导数共用采样时间，两个通道点数相同且各至少三个点。
|
|
|
if(!breakdown.valid || breakdown.residualVector.size() < 6 ||
|
|
|
breakdown.residualVector.size() % 2 != 0) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
for(int i = 0; i < breakdown.residualVector.size(); ++i) {
|
|
|
if(!isFiniteNumber(breakdown.residualVector[i])) {
|
|
|
return false;
|
|
|
}
|
|
|
}
|
|
|
const int pointCount = breakdown.residualVector.size() / 2;
|
|
|
const int slopeCount = pointCount - autoFitShapeLag(pointCount);
|
|
|
if(breakdown.shapeResiduals.size() != 2 * slopeCount || !isFiniteNumber(breakdown.shapeLoss)) return false;
|
|
|
for(int i = 0; i < breakdown.shapeResiduals.size(); ++i) {
|
|
|
if(!isFiniteNumber(breakdown.shapeResiduals[i])) return false;
|
|
|
}
|
|
|
return true;
|
|
|
}
|
|
|
|
|
|
// 数值和整体形状残差共用一次求解，联合缓存供两个 LM 目标复用。
|
|
|
static QVector<double> trustRegionFullResidual(const AutoFitObjectiveBreakdownLM& objective)
|
|
|
{
|
|
|
QVector<double> residual = objective.residualVector;
|
|
|
residual += objective.shapeResiduals;
|
|
|
return residual;
|
|
|
}
|
|
|
|
|
|
// 计算向量二范数的平方，避免在只比较能量或计算正规方程时反复开方。
|
|
|
static double trustRegionSquaredNorm(const QVector<double>& values)
|
|
|
{
|
|
|
double sum = 0.0;
|
|
|
for(int i = 0; i < values.size(); ++i) {
|
|
|
sum += values[i] * values[i];
|
|
|
}
|
|
|
return sum;
|
|
|
}
|
|
|
|
|
|
// 计算同维向量内积；维度不一致表示局部模型无效，返回零让调用方放弃修正。
|
|
|
static double trustRegionDotProduct(const QVector<double>& left,
|
|
|
const QVector<double>& right)
|
|
|
{
|
|
|
if(left.size() != right.size()) {
|
|
|
return 0.0;
|
|
|
}
|
|
|
|
|
|
double sum = 0.0;
|
|
|
for(int i = 0; i < left.size(); ++i) {
|
|
|
sum += left[i] * right[i];
|
|
|
}
|
|
|
return sum;
|
|
|
}
|
|
|
|
|
|
// 求解选中参数对应的阻尼正规方程。当前最多六个拟合参数，使用带部分主元的
|
|
|
// 高斯消元处理该小矩阵，并在主元退化时明确返回失败。
|
|
|
static bool solveTrustRegionLinearSystem(
|
|
|
QVector<QVector<double> > matrix,
|
|
|
QVector<double> rightHandSide,
|
|
|
QVector<double>* solution)
|
|
|
{
|
|
|
if(!solution || matrix.isEmpty() ||
|
|
|
matrix.size() != rightHandSide.size()) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
const int size = matrix.size();
|
|
|
for(int i = 0; i < size; ++i) {
|
|
|
if(matrix[i].size() != size) {
|
|
|
return false;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
for(int column = 0; column < size; ++column) {
|
|
|
int pivotRow = column;
|
|
|
double pivotMagnitude = qAbs(matrix[column][column]);
|
|
|
for(int row = column + 1; row < size; ++row) {
|
|
|
double magnitude = qAbs(matrix[row][column]);
|
|
|
if(magnitude > pivotMagnitude) {
|
|
|
pivotMagnitude = magnitude;
|
|
|
pivotRow = row;
|
|
|
}
|
|
|
}
|
|
|
if(pivotMagnitude <= 1.0e-14) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(pivotRow != column) {
|
|
|
qSwap(matrix[pivotRow], matrix[column]);
|
|
|
qSwap(rightHandSide[pivotRow], rightHandSide[column]);
|
|
|
}
|
|
|
|
|
|
for(int row = column + 1; row < size; ++row) {
|
|
|
double factor = matrix[row][column] /
|
|
|
matrix[column][column];
|
|
|
matrix[row][column] = 0.0;
|
|
|
for(int nextColumn = column + 1;
|
|
|
nextColumn < size; ++nextColumn) {
|
|
|
matrix[row][nextColumn] -=
|
|
|
factor * matrix[column][nextColumn];
|
|
|
}
|
|
|
rightHandSide[row] -= factor * rightHandSide[column];
|
|
|
}
|
|
|
}
|
|
|
|
|
|
solution->fill(0.0, size);
|
|
|
for(int row = size - 1; row >= 0; --row) {
|
|
|
double value = rightHandSide[row];
|
|
|
for(int column = row + 1; column < size; ++column) {
|
|
|
value -= matrix[row][column] * (*solution)[column];
|
|
|
}
|
|
|
double pivot = matrix[row][row];
|
|
|
if(qAbs(pivot) <= 1.0e-14) {
|
|
|
return false;
|
|
|
}
|
|
|
(*solution)[row] = value / pivot;
|
|
|
if(!isFiniteNumber((*solution)[row])) {
|
|
|
return false;
|
|
|
}
|
|
|
}
|
|
|
return true;
|
|
|
}
|
|
|
|
|
|
// LM 正规方程所需的 J^T*J 与 J^T*r，使用当前归一化参数坐标。
|
|
|
struct TrustRegionNormalEquations
|
|
|
{
|
|
|
QVector<QVector<double> > matrix;
|
|
|
QVector<double> gradient;
|
|
|
|
|
|
explicit TrustRegionNormalEquations(int dimensions = 0)
|
|
|
: matrix(dimensions, QVector<double>(dimensions, 0.0))
|
|
|
, gradient(dimensions, 0.0)
|
|
|
{}
|
|
|
};
|
|
|
|
|
|
static TrustRegionNormalEquations buildTrustRegionNormalEquations(
|
|
|
const QVector<QVector<double> >& jacobian,
|
|
|
const QVector<double>& residual,
|
|
|
const QVector<bool>& columnValid)
|
|
|
{
|
|
|
// 残差和 J 已包含压力、导数及采样权重，直接累加整体矩阵，避免重复计权。
|
|
|
const int dimensions = columnValid.size();
|
|
|
TrustRegionNormalEquations equations(dimensions);
|
|
|
for(int row = 0; row < residual.size(); ++row) {
|
|
|
for(int p = 0; p < dimensions; ++p) {
|
|
|
if(!columnValid[p]) continue;
|
|
|
equations.gradient[p] += jacobian[row][p] * residual[row];
|
|
|
for(int q = 0; q <= p; ++q) {
|
|
|
if(columnValid[q])
|
|
|
equations.matrix[p][q] += jacobian[row][p] * jacobian[row][q];
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
for(int p = 0; p < dimensions; ++p) {
|
|
|
for(int q = 0; q < p; ++q)
|
|
|
equations.matrix[q][p] = equations.matrix[p][q];
|
|
|
}
|
|
|
return equations;
|
|
|
}
|
|
|
|
|
|
static bool buildTrustRegionStep(
|
|
|
const TrustRegionNormalEquations& global,
|
|
|
const QVector<int>& selected,
|
|
|
const QVector<double>& coordinates,
|
|
|
double damping, double trustRadius, double minimumStep,
|
|
|
QVector<double>* step, double* predictedReduction)
|
|
|
{
|
|
|
// 复用原有全局 LM 方程、阻尼缩放、退化回退、信赖半径和边界投影。
|
|
|
const int dimensions = coordinates.size();
|
|
|
const int count = selected.size();
|
|
|
QVector<QVector<double> > normal(count, QVector<double>(count, 0.0));
|
|
|
QVector<double> rhs(count, 0.0);
|
|
|
for(int i = 0; i < count; ++i) {
|
|
|
const int p = selected[i];
|
|
|
rhs[i] = -global.gradient[p];
|
|
|
for(int j = 0; j < count; ++j) {
|
|
|
normal[i][j] = global.matrix[p][selected[j]];
|
|
|
}
|
|
|
normal[i][i] += damping * qMax(1.0e-10, normal[i][i]);
|
|
|
}
|
|
|
QVector<double> solution;
|
|
|
bool solved = solveTrustRegionLinearSystem(normal, rhs, &solution);
|
|
|
step->fill(0.0, dimensions);
|
|
|
if(solved) {
|
|
|
for(int i = 0; i < count; ++i) {
|
|
|
(*step)[selected[i]] = solution[i];
|
|
|
}
|
|
|
}
|
|
|
auto projectAndPredict = [&]() -> bool {
|
|
|
const double norm = qSqrt(trustRegionSquaredNorm(*step));
|
|
|
if(!isFiniteNumber(norm) || norm < minimumStep) return false;
|
|
|
if(norm > trustRadius) {
|
|
|
for(int p = 0; p < dimensions; ++p) (*step)[p] *= trustRadius / norm;
|
|
|
}
|
|
|
for(int p = 0; p < dimensions; ++p) {
|
|
|
(*step)[p] = qBound(0.0, coordinates[p] + (*step)[p], 1.0) - coordinates[p];
|
|
|
}
|
|
|
// 只比较投影后可执行步长的下降，阻尼项不属于真实拟合目标。
|
|
|
*predictedReduction = -trustRegionDotProduct(global.gradient, *step);
|
|
|
for(int p = 0; p < dimensions; ++p) {
|
|
|
for(int q = 0; q < dimensions; ++q) {
|
|
|
*predictedReduction -= 0.5 * (*step)[p] * global.matrix[p][q] * (*step)[q];
|
|
|
}
|
|
|
}
|
|
|
return isFiniteNumber(*predictedReduction) && *predictedReduction > 1.0e-14 &&
|
|
|
qSqrt(trustRegionSquaredNorm(*step)) >= minimumStep;
|
|
|
};
|
|
|
if(solved && projectAndPredict()) return true;
|
|
|
|
|
|
// 联合解可能被边界投影破坏；此时尝试可行梯度方向，不直接判为没有下降方向。
|
|
|
step->fill(0.0, dimensions);
|
|
|
for(int i = 0; i < count; ++i) {
|
|
|
const int p = selected[i];
|
|
|
double direction = -global.gradient[p];
|
|
|
if((coordinates[p] <= minimumStep && direction < 0.0) ||
|
|
|
(coordinates[p] >= 1.0 - minimumStep && direction > 0.0)) direction = 0.0;
|
|
|
(*step)[p] = direction;
|
|
|
}
|
|
|
const double norm = qSqrt(trustRegionSquaredNorm(*step));
|
|
|
if(!isFiniteNumber(norm) || norm < minimumStep) return false;
|
|
|
for(int p = 0; p < dimensions; ++p) {
|
|
|
(*step)[p] = qBound(0.0, coordinates[p] + (*step)[p] * trustRadius / norm, 1.0) - coordinates[p];
|
|
|
}
|
|
|
const double linearReduction = -trustRegionDotProduct(global.gradient, *step);
|
|
|
double curvature = 0.0;
|
|
|
for(int p = 0; p < dimensions; ++p) {
|
|
|
for(int q = 0; q < dimensions; ++q) {
|
|
|
curvature += (*step)[p] * global.matrix[p][q] * (*step)[q];
|
|
|
}
|
|
|
}
|
|
|
if(!isFiniteNumber(linearReduction) || !isFiniteNumber(curvature) || linearReduction <= 0.0) return false;
|
|
|
// 沿投影梯度最小化局部二次模型，只缩步，不突破信赖域或参数边界。
|
|
|
const double scale = curvature > 0.0 ? qMin(1.0, linearReduction / curvature) : 1.0;
|
|
|
for(int p = 0; p < dimensions; ++p) (*step)[p] *= scale;
|
|
|
return projectAndPredict();
|
|
|
}
|
|
|
|
|
|
// 每次得到有效真实候选后，使用满足最新割线条件的秩一修正更新完整残差
|
|
|
// Jacobian。这样模型吸收了刚得到的真实变化，又不必立即逐参数重新试算。
|
|
|
static void updateTrustRegionJacobian(
|
|
|
QVector<QVector<double> >* jacobian,
|
|
|
const QVector<double>& oldResidual,
|
|
|
const QVector<double>& newResidual,
|
|
|
const QVector<double>& coordinateStep)
|
|
|
{
|
|
|
if(!jacobian || jacobian->size() != oldResidual.size() ||
|
|
|
oldResidual.size() != newResidual.size()) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
double denominator = trustRegionSquaredNorm(coordinateStep);
|
|
|
if(denominator <= 1.0e-12) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
for(int row = 0; row < jacobian->size(); ++row) {
|
|
|
if((*jacobian)[row].size() != coordinateStep.size()) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
double predictedChange = 0.0;
|
|
|
for(int column = 0; column < coordinateStep.size(); ++column) {
|
|
|
predictedChange +=
|
|
|
(*jacobian)[row][column] * coordinateStep[column];
|
|
|
}
|
|
|
double correction =
|
|
|
(newResidual[row] - oldResidual[row] - predictedChange) /
|
|
|
denominator;
|
|
|
for(int column = 0; column < coordinateStep.size(); ++column) {
|
|
|
(*jacobian)[row][column] +=
|
|
|
correction * coordinateStep[column];
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
|
|
|
static QStringList traceParameterNames()
|
|
|
{
|
|
|
QStringList names;
|
|
|
names << "k"
|
|
|
<< "skin"
|
|
|
<< "wellboreC"
|
|
|
<< "phi"
|
|
|
<< "Dfc"
|
|
|
<< "fractureHalfLength";
|
|
|
return names;
|
|
|
}
|
|
|
|
|
|
nmCalculationAutoFitLM::nmCalculationAutoFitLM(QObject* parent)
|
|
|
: QObject(parent)
|
|
|
, m_isRunning(false)
|
|
|
, m_shouldStop(false)
|
|
|
, m_isFinalizing(false)
|
|
|
, m_currentIteration(-1)
|
|
|
, m_globalBestFitness(1e10)
|
|
|
, m_comparisonTimeMin(0.0)
|
|
|
, m_comparisonTimeMax(0.0)
|
|
|
, m_maxIterations(100)
|
|
|
, m_targetError(0.001)
|
|
|
, m_totalEvaluations(0)
|
|
|
, m_solverCallCount(0)
|
|
|
, m_successfulEvaluations(0)
|
|
|
, m_evaluationInProgress(0)
|
|
|
, m_userInitialFitness(1e10)
|
|
|
, m_maxConsecutiveFailures(3)
|
|
|
, m_hasValidUserSolution(false)
|
|
|
, m_targetWellName("")
|
|
|
, m_traceRunId("")
|
|
|
, m_traceFilePath("")
|
|
|
, m_traceMetaFilePath("")
|
|
|
{
|
|
|
// 单次运行目录在输入校验通过后创建，避免只打开界面也产生临时文件。
|
|
|
DEBUG_OUT("LM automatic fitting calculator initialized");
|
|
|
}
|
|
|
|
|
|
// 请求停止并限时等待，然后关闭轨迹文件、清理临时目录和信号连接。
|
|
|
// 等待超时后仅清除运行标志，不代表外层同步调用已经退出。
|
|
|
nmCalculationAutoFitLM::~nmCalculationAutoFitLM()
|
|
|
{
|
|
|
if(m_isRunning) {
|
|
|
m_shouldStop = true;
|
|
|
int waitCount = 0;
|
|
|
while(m_isRunning && waitCount < 100) {
|
|
|
QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, 50);
|
|
|
msleep(50);
|
|
|
waitCount++;
|
|
|
}
|
|
|
if(m_isRunning) {
|
|
|
DEBUG_OUT("Force stopping LM fitting after timeout");
|
|
|
m_isRunning = false;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
closeTraceFile();
|
|
|
cleanupTemporaryDirectory();
|
|
|
disconnect(this, nullptr, nullptr, nullptr);
|
|
|
}
|
|
|
|
|
|
// ==================== 临时目录工具 ====================
|
|
|
// 真实求解器和井曲线 CSV 共用一个受控的单次运行目录。创建失败时不允许
|
|
|
// 回退到程序目录，否则后续递归清理可能删除安装文件。
|
|
|
bool nmCalculationAutoFitLM::initializeTemporaryDirectory()
|
|
|
{
|
|
|
nmCalculationUtils::cleanupStaleAutoFitTemporaryDirectories();
|
|
|
m_tempDirectory = nmCalculationUtils::createAutoFitTemporaryDirectory();
|
|
|
if(m_tempDirectory.isEmpty()) {
|
|
|
DEBUG_OUT("Failed to create LM auto-fit temporary directory");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
m_traceRunId = QFileInfo(m_tempDirectory).fileName();
|
|
|
DEBUG_OUT(QString("Initialized temp directory: %1").arg(m_tempDirectory));
|
|
|
return true;
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::cleanupTemporaryDirectory()
|
|
|
{
|
|
|
if(m_tempDirectory.isEmpty()) {
|
|
|
return;
|
|
|
}
|
|
|
const QString directoryPath = m_tempDirectory;
|
|
|
m_tempDirectory.clear();
|
|
|
if(nmCalculationUtils::removeAutoFitTemporaryDirectory(directoryPath)) {
|
|
|
DEBUG_OUT("Temp directory cleaned up successfully");
|
|
|
} else {
|
|
|
DEBUG_OUT(QString("Failed to clean auto-fit temporary directory: %1")
|
|
|
.arg(directoryPath));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
|
|
|
// ==================== 公共接口方法 ====================
|
|
|
|
|
|
void nmCalculationAutoFitLM::setTargetLogLogData(const QVector<QVector<double> >& targetData)
|
|
|
{
|
|
|
// 目标曲线由界面层从目标井 history log-log 传入。
|
|
|
// 约定 targetData[0]=time，targetData[1]=pressure，targetData[2]=pressure derivative。
|
|
|
m_targetLogLogData = targetData;
|
|
|
m_comparisonTimeMin = 0.0;
|
|
|
m_comparisonTimeMax = 0.0;
|
|
|
DEBUG_OUT(QString("Target LogLog data set: %1 arrays").arg(targetData.size()));
|
|
|
|
|
|
if(targetData.size() >= 3) {
|
|
|
DEBUG_OUT(QString("LogLog data points: X=%1, Y1=%2, Y2=%3")
|
|
|
.arg(targetData[0].size())
|
|
|
.arg(targetData[1].size())
|
|
|
.arg(targetData[2].size()));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
|
|
|
void nmCalculationAutoFitLM::stopFitting()
|
|
|
{
|
|
|
if(m_isFinalizing) {
|
|
|
DEBUG_OUT(
|
|
|
tr("The final full-field calculation is already running"));
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
// 用户点击停止时只设置请求标志，让 LM 主循环和求解器等待逻辑自然退出。
|
|
|
if(m_isRunning) {
|
|
|
DEBUG_OUT(tr("=== User Stop Request Received ==="));
|
|
|
DEBUG_OUT(tr("Gracefully stopping LM automatic fitting..."));
|
|
|
m_shouldStop = true;
|
|
|
|
|
|
// 此槽由求解等待循环派发，必须立即返回，外层循环才能转发取消请求。
|
|
|
if(m_evaluationInProgress > 0) {
|
|
|
emit logMessageGenerated(tr("Waiting for current solver evaluation to stop..."));
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(tr("LM automatic fitting stop request processed"));
|
|
|
} else {
|
|
|
DEBUG_OUT(tr("Stop request received but optimization is not running"));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
bool nmCalculationAutoFitLM::isRunning() const
|
|
|
{
|
|
|
// UI 查询当前是否处于自动拟合运行状态。
|
|
|
return m_isRunning;
|
|
|
}
|
|
|
|
|
|
int nmCalculationAutoFitLM::getCurrentIteration() const
|
|
|
{
|
|
|
// UI 进度条和日志展示用的当前迭代序号。
|
|
|
return m_currentIteration;
|
|
|
}
|
|
|
|
|
|
int nmCalculationAutoFitLM::getTotalEvaluations() const
|
|
|
{
|
|
|
// 返回评价次数，用于迭代预算；真实求解器调用次数另行统计。
|
|
|
return m_totalEvaluations;
|
|
|
}
|
|
|
|
|
|
int nmCalculationAutoFitLM::getSolverCallCount() const
|
|
|
{
|
|
|
// 包含失败调用、重试和最终完整求解，不包含未进入求解器的评价。
|
|
|
return m_solverCallCount;
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::recordSolverCall()
|
|
|
{
|
|
|
// 信号来自求解线程，使用原子累加避免跨线程读写竞争。
|
|
|
m_solverCallCount.fetchAndAddOrdered(1);
|
|
|
}
|
|
|
|
|
|
QVector<double> nmCalculationAutoFitLM::getBestSolution() const
|
|
|
{
|
|
|
// 返回紧凑的“启用参数向量”，顺序与 m_enabledParamIndices 一致。
|
|
|
return m_globalBestPosition;
|
|
|
}
|
|
|
|
|
|
double nmCalculationAutoFitLM::getBestFitness() const
|
|
|
{
|
|
|
// 当前已接受工作点的总误差；高度和形状阶段不保证总误差单调下降。
|
|
|
return m_globalBestFitness;
|
|
|
}
|
|
|
|
|
|
AutoFitObjectiveBreakdownLM nmCalculationAutoFitLM::getLastObjectiveBreakdown() const
|
|
|
{
|
|
|
// 返回最近一次损失评价的误差分解，供界面或后续优化逻辑读取。
|
|
|
return m_lastObjectiveBreakdown;
|
|
|
}
|
|
|
|
|
|
QString nmCalculationAutoFitLM::getLastError() const
|
|
|
{
|
|
|
// 上一次失败的人类可读错误信息，主要给 UI 层弹窗或日志使用。
|
|
|
return m_lastError;
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::resetOptimizer()
|
|
|
{
|
|
|
// 清空一次运行产生的状态，但不销毁对象。
|
|
|
// 配置字段会在 startAutoFitting() 中重新从 DataManager 读取；
|
|
|
// trace 文件先关闭，避免新一轮 run 继续写到旧 CSV。
|
|
|
closeTraceFile();
|
|
|
m_globalBestPosition.clear();
|
|
|
m_globalBestFitness = 1e10;
|
|
|
m_globalBestObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
|
|
|
m_lastEvaluatedLogLogData.clear();
|
|
|
m_globalBestLogLogData.clear();
|
|
|
m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
|
|
|
m_userInitialLogLogData.clear();
|
|
|
m_userInitialObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
|
|
|
m_comparisonTimeMin = 0.0;
|
|
|
m_comparisonTimeMax = 0.0;
|
|
|
m_currentIteration = -1; // 尚未执行迭代；对外日志使用索引加一。
|
|
|
m_totalEvaluations = 0;
|
|
|
m_solverCallCount = 0;
|
|
|
m_successfulEvaluations = 0;
|
|
|
m_lastError.clear();
|
|
|
m_initialValues.clear();
|
|
|
m_userInitialSolution.clear();
|
|
|
m_userInitialFitness = 1e10;
|
|
|
m_hasValidUserSolution = false;
|
|
|
m_traceMetaFilePath.clear();
|
|
|
m_isFinalizing = false;
|
|
|
|
|
|
DEBUG_OUT("LM optimizer reset");
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::setTargetWellName(const QString& wellName)
|
|
|
{
|
|
|
// 目标井名是贯穿拟合流程的关键索引：
|
|
|
// 读目标曲线、写 skin/wellboreC、求解后取 resultLogLog 都依赖这个名字。
|
|
|
m_targetWellName = wellName;
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::initializeTraceFile()
|
|
|
{
|
|
|
closeTraceFile();
|
|
|
m_traceFilePath.clear();
|
|
|
m_traceMetaFilePath.clear();
|
|
|
|
|
|
// 诊断轨迹在所有构建中启用，并与源码、安装目录完全分离。
|
|
|
if(m_traceRunId.isEmpty()) {
|
|
|
m_traceRunId = QString("%1-%2")
|
|
|
.arg(QString::number(QCoreApplication::applicationPid()))
|
|
|
.arg(QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss_zzz"));
|
|
|
}
|
|
|
QDir traceDir(QDir(QDir::tempPath()).absoluteFilePath("WTAI/AutoFitTrace/LM"));
|
|
|
|
|
|
if(!traceDir.exists() && !QDir().mkpath(traceDir.absolutePath())) {
|
|
|
DEBUG_OUT(QString("Failed to create LM trace directory: %1").arg(traceDir.absolutePath()));
|
|
|
m_traceFilePath.clear();
|
|
|
m_traceMetaFilePath.clear();
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
m_traceFilePath = traceDir.absoluteFilePath(
|
|
|
QString("lm_trust_region_trace_%1.csv").arg(m_traceRunId));
|
|
|
m_traceMetaFilePath = traceDir.absoluteFilePath(
|
|
|
QString("lm_trust_region_trace_%1.meta.json").arg(m_traceRunId));
|
|
|
m_traceFile.setFileName(m_traceFilePath);
|
|
|
|
|
|
if(!m_traceFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
|
|
|
DEBUG_OUT(QString("Failed to open LM trace file: %1").arg(m_traceFilePath));
|
|
|
m_traceFilePath.clear();
|
|
|
m_traceMetaFilePath.clear();
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
writeTraceHeader();
|
|
|
writeTraceMetaFile();
|
|
|
DEBUG_OUT(tr("LM fitting trace: %1").arg(m_traceFilePath));
|
|
|
if(!m_traceMetaFilePath.isEmpty()) {
|
|
|
DEBUG_OUT(
|
|
|
tr("LM fitting trace metadata: %1").arg(m_traceMetaFilePath));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::closeTraceFile()
|
|
|
{
|
|
|
if(m_traceFile.isOpen()) {
|
|
|
m_traceFile.flush();
|
|
|
m_traceFile.close();
|
|
|
}
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::writeTraceHeader()
|
|
|
{
|
|
|
if(!m_traceFile.isOpen()) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
QStringList cols;
|
|
|
cols << "run_id"
|
|
|
<< "iteration"
|
|
|
<< "parameter_index"
|
|
|
<< "phase"
|
|
|
<< "k"
|
|
|
<< "skin"
|
|
|
<< "wellboreC"
|
|
|
<< "phi"
|
|
|
<< "Dfc"
|
|
|
<< "fractureHalfLength"
|
|
|
<< "solver_objective"
|
|
|
<< "solver_success"
|
|
|
<< "elapsed_ms"
|
|
|
<< "decision"
|
|
|
<< "enabled_param_indices"
|
|
|
<< "pressure_loss"
|
|
|
<< "derivative_loss"
|
|
|
<< "vertical_common_bias"
|
|
|
<< "vertical_loss"
|
|
|
<< "vertical_reliable"
|
|
|
<< "horizontal_physical_shift"
|
|
|
<< "horizontal_loss"
|
|
|
<< "horizontal_reliable"
|
|
|
<< "shape_loss"
|
|
|
<< "late_trend_loss"
|
|
|
<< "late_slope_bias"
|
|
|
<< "late_trend_reliable"
|
|
|
<< "registration_ambiguous";
|
|
|
|
|
|
cols << "sampling_mode" << "sampling_points"
|
|
|
<< "pressure_vertical_bias" << "derivative_vertical_bias";
|
|
|
QTextStream out(&m_traceFile);
|
|
|
out << cols.join(",") << "\n";
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::writeTraceMetaFile()
|
|
|
{
|
|
|
if(m_traceMetaFilePath.isEmpty()) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
QFile metaFile(m_traceMetaFilePath);
|
|
|
if(!metaFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
|
|
|
DEBUG_OUT(QString("Failed to open LM trace meta file: %1").arg(m_traceMetaFilePath));
|
|
|
m_traceMetaFilePath.clear();
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
QStringList parameterNames = traceParameterNames();
|
|
|
QStringList enabledNames;
|
|
|
for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
|
|
|
int paramIndex = m_enabledParamIndices[i];
|
|
|
enabledNames << ((paramIndex >= 0 && paramIndex < parameterNames.size())
|
|
|
? parameterNames[paramIndex]
|
|
|
: QString::number(paramIndex));
|
|
|
}
|
|
|
|
|
|
QVector<double> initialFullParams = buildTraceParameterVector(m_initialValues);
|
|
|
QVector<double> targetTime = m_targetLogLogData.size() > 0
|
|
|
? m_targetLogLogData[0] : QVector<double>();
|
|
|
QVector<double> targetPressure = m_targetLogLogData.size() > 1
|
|
|
? m_targetLogLogData[1] : QVector<double>();
|
|
|
QVector<double> targetDerivative = m_targetLogLogData.size() > 2
|
|
|
? m_targetLogLogData[2] : QVector<double>();
|
|
|
|
|
|
QTextStream out(&metaFile);
|
|
|
out << "{\n";
|
|
|
out << " \"schema_version\": 56,\n";
|
|
|
out << " \"strategy\": \"permeability_height_then_joint_lm_total_for_vertical_well_else_shape_then_total\",\n";
|
|
|
out << " \"height_acceptance\": \"reliable_vertical_loss_decrease; no_shape_or_total_loss_constraint\",\n";
|
|
|
out << " \"shape_metric\": \"pressure_and_derivative_log_slopes_shared_sampling_grid_lag_round_10_percent_min_1\",\n";
|
|
|
out << " \"shape_parameter_selection\": \"all_valid_free_columns_joint_LM; no_single_parameter_sweep_or_wellbore_recheck\",\n";
|
|
|
out << " \"shape_acceptance\": \"strict_full_shape_loss_decrease; no_total_loss_constraint\",\n";
|
|
|
out << " \"lm_step_policy\": \"same_joint_LM_damping_trust_radius_secant_updates_and_sensitivity_rebuilds_in_both_phases\",\n";
|
|
|
out << " \"lm_phase_switch\": \"shape_target_reached_or_confirmed_stagnation_or_phase_evaluation_budget_exhausted_then_total_if_iterations_remain; shared_iteration_limit_or_user_stop_or_consecutive_solver_failure_abort\",\n";
|
|
|
out << " \"lm_phase_budget\": \"shared_max_iterations_including_height; each_LM_phase_retains_max(dimensions+2,20,3*max_iterations)_evaluations\",\n";
|
|
|
out << " \"lm_effective_improvement\": \"max(0.01,0.10*baseline_stage_error); after_sensitivity_rebuild_reset_baseline_to_accepted_probe_error_if_improved_else_rebuild_start_error; reset_ineffective_count_preserving_stagnation_confirmation; 3_ineffective_steps_trigger_stagnation_confirmation\",\n";
|
|
|
out << " \"iteration_count_scope\": \"height_shape_and_total_share_max_iterations; rejected_and_failed_candidates_count; sensitivity_evaluations_are_not_iterations; trace_iteration_is_zero_based_cumulative\",\n";
|
|
|
out << " \"total_stage_shape_constraint\": false,\n";
|
|
|
out << " \"total_parameter_selection\": \"all_valid_free_columns\",\n";
|
|
|
out << " \"skin_difference_policy\": \"local_scale_fraction_0.02_in_both_LM_phases\",\n";
|
|
|
out << " \"difference_policy\": \"coordinate_step_cap_0.04_with_half_trust_radius_in_both_LM_phases\",\n";
|
|
|
out << " \"difference_failure_policy\": \"halve_before_opposite_direction\",\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("log_uniform_per_decade") << ",\n";
|
|
|
out << " \"sampling_intervals_per_decade\": " << kAutoFitIntervalsPerDecade << ",\n";
|
|
|
out << " \"target\": {\n";
|
|
|
out << " \"well_name\": " << jsonEscape(m_targetWellName) << ",\n";
|
|
|
out << " \"time\": " << jsonDoubleArray(targetTime) << ",\n";
|
|
|
out << " \"pressure\": " << jsonDoubleArray(targetPressure) << ",\n";
|
|
|
out << " \"derivative\": " << jsonDoubleArray(targetDerivative) << "\n";
|
|
|
out << " },\n";
|
|
|
out << " \"lm\": {\n";
|
|
|
out << " \"max_iterations\": " << m_maxIterations << ",\n";
|
|
|
out << " \"target_error\": " << jsonNumber(m_targetError) << "\n";
|
|
|
out << " },\n";
|
|
|
// 区间尚未冻结时写 null；首次有效评价后记录固定的误差比较范围。
|
|
|
out << " \"comparison_time_range\": {\n";
|
|
|
out << " \"comparison_time_min\": "
|
|
|
<< (m_comparisonTimeMin > 0.0 ? jsonNumber(m_comparisonTimeMin) : QString("null")) << ",\n";
|
|
|
out << " \"comparison_time_max\": "
|
|
|
<< (m_comparisonTimeMax > 0.0 ? jsonNumber(m_comparisonTimeMax) : QString("null")) << "\n";
|
|
|
out << " },\n";
|
|
|
out << " \"parameters\": {\n";
|
|
|
out << " \"names\": " << jsonStringArray(parameterNames) << ",\n";
|
|
|
out << " \"enabled_indices\": " << jsonIntArray(m_enabledParamIndices) << ",\n";
|
|
|
out << " \"enabled_names\": " << jsonStringArray(enabledNames) << ",\n";
|
|
|
out << " \"selected_flags\": " << jsonBoolArray(m_parameterSelected) << ",\n";
|
|
|
out << " \"lower\": " << jsonDoubleArray(m_parameterLower) << ",\n";
|
|
|
out << " \"upper\": " << jsonDoubleArray(m_parameterUpper) << ",\n";
|
|
|
out << " \"initial_selected\": " << jsonDoubleArray(m_initialValues) << ",\n";
|
|
|
out << " \"initial_full\": " << jsonDoubleArray(initialFullParams) << "\n";
|
|
|
out << " }\n";
|
|
|
out << "}\n";
|
|
|
metaFile.flush();
|
|
|
metaFile.close();
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::writeTraceRow(
|
|
|
int iteration,
|
|
|
int parameterIndex,
|
|
|
const QString& phase,
|
|
|
const QVector<double>& parameters,
|
|
|
double solverObjective,
|
|
|
bool solverSuccess,
|
|
|
int elapsedMs,
|
|
|
const QString& decision,
|
|
|
const AutoFitObjectiveBreakdownLM* objectiveBreakdown)
|
|
|
{
|
|
|
if(!m_traceFile.isOpen()) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
QVector<double> fullParams = buildTraceParameterVector(parameters);
|
|
|
QStringList enabledIndices;
|
|
|
for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
|
|
|
enabledIndices << QString::number(m_enabledParamIndices[i]);
|
|
|
}
|
|
|
|
|
|
QStringList cols;
|
|
|
cols << csvEscape(m_traceRunId)
|
|
|
<< QString::number(iteration)
|
|
|
<< QString::number(parameterIndex)
|
|
|
<< csvEscape(phase);
|
|
|
for(int i = 0; i < 6; ++i) {
|
|
|
cols << traceParamAt(fullParams, i);
|
|
|
}
|
|
|
cols << traceNumber(solverObjective)
|
|
|
<< QString::number(solverSuccess ? 1 : 0)
|
|
|
<< QString::number(elapsedMs)
|
|
|
<< csvEscape(decision)
|
|
|
<< csvEscape(enabledIndices.join(";"));
|
|
|
|
|
|
if(objectiveBreakdown && objectiveBreakdown->valid) {
|
|
|
cols << traceNumber(objectiveBreakdown->pressureLoss)
|
|
|
<< traceNumber(objectiveBreakdown->derivativeLoss)
|
|
|
<< traceNumber(objectiveBreakdown->verticalCommonBias)
|
|
|
<< traceNumber(objectiveBreakdown->verticalLoss)
|
|
|
<< QString::number(objectiveBreakdown->verticalReliable ? 1 : 0)
|
|
|
<< traceNumber(objectiveBreakdown->horizontalPhysicalShift)
|
|
|
<< traceNumber(objectiveBreakdown->horizontalLoss)
|
|
|
<< QString::number(objectiveBreakdown->horizontalReliable ? 1 : 0)
|
|
|
<< traceNumber(objectiveBreakdown->shapeLoss)
|
|
|
<< traceNumber(objectiveBreakdown->lateDerivativeTrendLoss)
|
|
|
<< traceNumber(objectiveBreakdown->lateDerivativeSlopeBias)
|
|
|
<< QString::number(objectiveBreakdown->lateDerivativeTrendReliable ? 1 : 0)
|
|
|
<< QString::number(objectiveBreakdown->registrationAmbiguous ? 1 : 0);
|
|
|
} else {
|
|
|
for(int i = 0; i < 13; ++i) {
|
|
|
cols << QString();
|
|
|
}
|
|
|
}
|
|
|
|
|
|
cols << "log_uniform_per_decade"
|
|
|
<< (objectiveBreakdown ? QString::number(objectiveBreakdown->residualVector.size() / 2) : QString())
|
|
|
<< (objectiveBreakdown ? traceNumber(objectiveBreakdown->pressureVerticalBias) : QString())
|
|
|
<< (objectiveBreakdown ? traceNumber(objectiveBreakdown->derivativeVerticalBias) : QString());
|
|
|
QTextStream out(&m_traceFile);
|
|
|
out << cols.join(",") << "\n";
|
|
|
m_traceFile.flush();
|
|
|
}
|
|
|
|
|
|
// 拟合界面统一显示普通小数；失败评价不展示惩罚值或非有限数。
|
|
|
static QString lmLogError(double value)
|
|
|
{
|
|
|
return isFiniteNumber(value) && value < 1.0e9
|
|
|
? QString::number(value, 'f', 6) : QString("--");
|
|
|
}
|
|
|
|
|
|
// 改善率统一相对上次已接受误差计算；起点为零时不作除法。
|
|
|
static QString lmLogImprovement(double beforeError, double afterError)
|
|
|
{
|
|
|
return beforeError > 0.0
|
|
|
? QString::number((beforeError - afterError) / beforeError * 100.0, 'f', 2)
|
|
|
: QString("--");
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::emitRunSummary(bool success, StopReasonLM finalReason)
|
|
|
{
|
|
|
// 最终只汇总总误差，耗时由进度窗口在收到结束信号后追加。
|
|
|
emit logMessageGenerated(success ? tr("[Fitting completed]") : tr("[Fitting ended]"));
|
|
|
emit logMessageGenerated(tr("Stop reason: %1").arg(success
|
|
|
? getStopReasonDescription(finalReason)
|
|
|
: (m_lastError.isEmpty() ? getStopReasonDescription(finalReason) : m_lastError)));
|
|
|
emit logMessageGenerated(tr("Total error: %1 -> %2")
|
|
|
.arg(lmLogError(m_userInitialFitness)).arg(lmLogError(m_globalBestFitness)));
|
|
|
emit logMessageGenerated(tr("Solver calls: %1").arg(getSolverCallCount()));
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::emitTrialLog(const QString& title,
|
|
|
const AutoFitObjectiveBreakdownLM& before, const AutoFitObjectiveBreakdownLM& after,
|
|
|
bool valid, bool accepted, bool shapeStage)
|
|
|
{
|
|
|
// 只显示当前阶段目标；改善率相对上次已接受的误差计算，变差时为负值。
|
|
|
if(!valid) {
|
|
|
emit logMessageGenerated(tr("%1 | Solver failed\nNo valid error was returned; keep the last accepted result and reduce the step.").arg(title));
|
|
|
return;
|
|
|
}
|
|
|
const QString status = accepted ? tr("accepted") : tr("rejected");
|
|
|
const double beforeError = shapeStage ? before.shapeLoss : before.total;
|
|
|
const double afterError = shapeStage ? after.shapeLoss : after.total;
|
|
|
const QString objective = (shapeStage ? tr("Shape error: %1 -> %2")
|
|
|
: tr("Total error: %1 -> %2"))
|
|
|
.arg(lmLogError(beforeError)).arg(lmLogError(afterError));
|
|
|
QString message = tr("%1 | %2 | %3 | Improvement: %4%")
|
|
|
.arg(title).arg(status).arg(objective).arg(lmLogImprovement(beforeError, afterError));
|
|
|
if(!accepted) message += tr("\n Objective error did not decrease; keep the last accepted result and reduce the step.");
|
|
|
emit logMessageGenerated(message);
|
|
|
}
|
|
|
|
|
|
QVector<double> nmCalculationAutoFitLM::buildTraceParameterVector(const QVector<double>& selectedParameters) const
|
|
|
{
|
|
|
// 将 LM 内部使用的“启用参数向量”还原成完整 6 维参数向量。
|
|
|
// 未启用的参数从当前 DataManager 读取，启用的参数用 selectedParameters 覆盖。
|
|
|
// trace CSV 和 meta 使用该完整向量记录一次候选评价。
|
|
|
QVector<double> fullParams(6, 0.0);
|
|
|
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
|
|
|
if(dataManager) {
|
|
|
nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
|
|
|
fullParams[0] = reservoirData.getPermeability().getValue().toDouble();
|
|
|
fullParams[3] = reservoirData.getPorosity().getValue().toDouble();
|
|
|
|
|
|
nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
|
|
|
|
|
|
if(pTargetWell) {
|
|
|
nmDataPerforation* perforation = pTargetWell->getPerforation(0);
|
|
|
if(perforation) {
|
|
|
fullParams[1] = perforation->getSkin().getValue().toDouble();
|
|
|
}
|
|
|
fullParams[2] = pTargetWell->getWellboreStorage().getValue().toDouble();
|
|
|
|
|
|
// Dfc 只存在于两类压裂井，普通井在完整向量中保持为 0。
|
|
|
if(pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
|
|
|
nmDataVerticalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataVerticalFracturedWell*>(pTargetWell);
|
|
|
if(fracturedWell) {
|
|
|
fullParams[4] = fracturedWell->getDfc().getValue().toDouble();
|
|
|
fullParams[5] = fracturedWell->getFractureHalfLength().getValue().toDouble();
|
|
|
}
|
|
|
} else if(pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
|
|
|
nmDataHorizontalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataHorizontalFracturedWell*>(pTargetWell);
|
|
|
if(fracturedWell) {
|
|
|
fullParams[4] = fracturedWell->getDfc().getValue().toDouble();
|
|
|
fullParams[5] = fracturedWell->getFractureHalfLength().getValue().toDouble();
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
|
|
|
for(int i = 0; i < selectedParameters.size() && i < m_enabledParamIndices.size(); ++i) {
|
|
|
int paramIndex = m_enabledParamIndices[i];
|
|
|
|
|
|
if(paramIndex >= 0 && paramIndex < fullParams.size()) {
|
|
|
fullParams[paramIndex] = selectedParameters[i];
|
|
|
}
|
|
|
}
|
|
|
|
|
|
return fullParams;
|
|
|
}
|
|
|
|
|
|
// ==================== 数据加载方法 ====================
|
|
|
|
|
|
bool nmCalculationAutoFitLM::loadAllConfigFromDataManager()
|
|
|
{
|
|
|
// 统一从 DataManager 加载本次运行所需配置。
|
|
|
// UI 层只负责把用户选择保存到 nmDataAutomaticFitting，本类从这里开始完全数据驱动。
|
|
|
try {
|
|
|
loadOptimizationConfig();
|
|
|
loadParameterBounds();
|
|
|
extractUserInitialValues();
|
|
|
return true;
|
|
|
} catch(...) {
|
|
|
m_lastError = tr("Failed to load configuration from data manager");
|
|
|
return false;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::loadOptimizationConfig()
|
|
|
{
|
|
|
// LM 只读取迭代次数和目标误差，其他数值控制保留在现有算法实现中。
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();
|
|
|
|
|
|
m_maxIterations = fittingData.getIterationCount().getValue().toInt();
|
|
|
m_targetError = fittingData.getErrorTolerance().getValue().toDouble();
|
|
|
DEBUG_OUT(QString("Loaded LM config: iterations=%1, error=%2")
|
|
|
.arg(m_maxIterations).arg(m_targetError));
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::loadParameterBounds()
|
|
|
{
|
|
|
// 读取用户勾选的拟合参数及上下界。
|
|
|
//
|
|
|
// 这里构建三个核心数组：
|
|
|
// - m_parameterSelected[6]：完整参数体系中每个参数是否参与拟合；
|
|
|
// - m_parameterLower/Upper[6]：完整参数体系的搜索上下界；
|
|
|
// - m_enabledParamIndices：将勾选参数向量映射回完整参数索引。
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();
|
|
|
|
|
|
m_parameterSelected.resize(6);
|
|
|
m_parameterSelected[0] = fittingData.getPermeabilitySelected();
|
|
|
m_parameterSelected[1] = fittingData.getSkinSelected();
|
|
|
m_parameterSelected[2] = fittingData.getWellboreStorageSelected();
|
|
|
m_parameterSelected[3] = fittingData.getPorositySelected();
|
|
|
m_parameterSelected[4] = fittingData.getFractureConductivitySelected();
|
|
|
m_parameterSelected[5] = fittingData.getFractureHalfLengthSelected();
|
|
|
|
|
|
m_parameterLower.resize(6);
|
|
|
m_parameterUpper.resize(6);
|
|
|
|
|
|
m_parameterLower[0] = fittingData.getPermeabilityMin().getValue().toDouble();
|
|
|
m_parameterUpper[0] = fittingData.getPermeabilityMax().getValue().toDouble();
|
|
|
|
|
|
m_parameterLower[1] = fittingData.getSkinMin().getValue().toDouble();
|
|
|
m_parameterUpper[1] = fittingData.getSkinMax().getValue().toDouble();
|
|
|
|
|
|
m_parameterLower[2] = fittingData.getWellboreStorageMin().getValue().toDouble();
|
|
|
m_parameterUpper[2] = fittingData.getWellboreStorageMax().getValue().toDouble();
|
|
|
|
|
|
m_parameterLower[3] = fittingData.getPorosityMin().getValue().toDouble();
|
|
|
m_parameterUpper[3] = fittingData.getPorosityMax().getValue().toDouble();
|
|
|
|
|
|
m_parameterLower[4] = fittingData.getFractureConductivityMin().getValue().toDouble();
|
|
|
m_parameterUpper[4] = fittingData.getFractureConductivityMax().getValue().toDouble();
|
|
|
|
|
|
m_parameterLower[5] = fittingData.getFractureHalfLengthMin().getValue().toDouble();
|
|
|
m_parameterUpper[5] = fittingData.getFractureHalfLengthMax().getValue().toDouble();
|
|
|
|
|
|
m_enabledParamIndices.clear();
|
|
|
|
|
|
for(int i = 0; i < m_parameterSelected.size(); ++i) {
|
|
|
if(m_parameterSelected[i]) {
|
|
|
m_enabledParamIndices.append(i);
|
|
|
}
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(QString("Loaded parameter bounds: %1 enabled parameters")
|
|
|
.arg(m_enabledParamIndices.size()));
|
|
|
}
|
|
|
|
|
|
// ==================== 自动拟合核心方法 ====================
|
|
|
bool nmCalculationAutoFitLM::startAutoFitting()
|
|
|
{
|
|
|
// 总入口只负责准备数据、调用有限差分 + LM/信赖域，并写回最终结果。
|
|
|
StopReasonLM finalReason = LM_CONTINUE_OPTIMIZATION;
|
|
|
|
|
|
if(m_isRunning) {
|
|
|
m_lastError = tr("Auto fitting is already running");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
try {
|
|
|
if(!loadAllConfigFromDataManager()) {
|
|
|
emit logMessageGenerated(tr("ERROR: Failed to load configuration from data manager"));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(tr("Algorithm: LM"));
|
|
|
const int enabledParams = getEnabledParameterCount();
|
|
|
DEBUG_OUT(tr("Enabled parameters count: %1").arg(enabledParams));
|
|
|
|
|
|
if(enabledParams == 0) {
|
|
|
m_lastError = tr("No parameters enabled for optimization");
|
|
|
emit logMessageGenerated(tr("ERROR: No parameters enabled for optimization"));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(m_targetLogLogData.size() < 3) {
|
|
|
m_lastError = tr("Target LogLog data is empty or insufficient");
|
|
|
emit logMessageGenerated(tr("ERROR: Target LogLog data is empty or insufficient"));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(m_targetLogLogData[0].size() != m_targetLogLogData[1].size() ||
|
|
|
m_targetLogLogData[0].size() != m_targetLogLogData[2].size()) {
|
|
|
m_lastError = tr("Target LogLog data arrays have inconsistent sizes");
|
|
|
emit logMessageGenerated(tr("ERROR: Target LogLog data arrays have inconsistent sizes"));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(m_targetWellName.isEmpty()) {
|
|
|
m_lastError = tr("Target well name is empty");
|
|
|
emit logMessageGenerated(tr("ERROR: Target well name is empty"));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(
|
|
|
tr("Target data validation passed (%1 data points)")
|
|
|
.arg(m_targetLogLogData[0].size()));
|
|
|
|
|
|
// resetOptimizer() 会清空运行状态，因此先保存从 DataManager 提取的初始值。
|
|
|
QVector<double> savedInitialValues = m_initialValues;
|
|
|
resetOptimizer();
|
|
|
// 输入校验通过后，为本次拟合创建新的临时计算目录。
|
|
|
cleanupTemporaryDirectory();
|
|
|
if(!initializeTemporaryDirectory()) {
|
|
|
m_lastError = tr("Cannot create the automatic fitting temporary directory");
|
|
|
emit logMessageGenerated(tr("ERROR: %1").arg(m_lastError));
|
|
|
return false;
|
|
|
}
|
|
|
m_isRunning = true;
|
|
|
m_shouldStop = false;
|
|
|
m_isFinalizing = false;
|
|
|
m_currentIteration = -1;
|
|
|
m_initialValues = savedInitialValues;
|
|
|
initializeTraceFile();
|
|
|
// 开头只展示一次初始参数和总误差，不展示形状误差。
|
|
|
emit logMessageGenerated(tr("========== LM fitting started =========="));
|
|
|
emit logMessageGenerated(tr("Target well: %1").arg(m_targetWellName));
|
|
|
emit logMessageGenerated(tr("Target error: %1").arg(lmLogError(m_targetError)));
|
|
|
emit fittingStageChanged(0);
|
|
|
|
|
|
// 先用真实求解器评价用户当前模型，供最终精英保护使用。
|
|
|
if(!savedInitialValues.isEmpty()) {
|
|
|
m_userInitialSolution = savedInitialValues;
|
|
|
DEBUG_OUT(tr("=== Evaluating Initial Solution (Elite Protection) ==="));
|
|
|
|
|
|
const char* names[] = { QT_TR_NOOP("Permeability"), QT_TR_NOOP("Skin"),
|
|
|
QT_TR_NOOP("Wellbore storage"), QT_TR_NOOP("Porosity"),
|
|
|
QT_TR_NOOP("Fracture conductivity"), QT_TR_NOOP("Fracture half length") };
|
|
|
const char* units[] = { "mD", "", "m^3/MPa", "", "mD.m", "m" };
|
|
|
emit logMessageGenerated(tr("[Initial parameters]"));
|
|
|
for(int i = 0; i < m_userInitialSolution.size(); ++i) {
|
|
|
const int index = m_enabledParamIndices[i];
|
|
|
emit logMessageGenerated(QString("%1: %2 %3").arg(tr(names[index]))
|
|
|
.arg(m_userInitialSolution[i], 0, 'g', 10).arg(tr(units[index])));
|
|
|
}
|
|
|
|
|
|
try {
|
|
|
DEBUG_OUT(
|
|
|
tr("Starting initial solution evaluation..."));
|
|
|
QTime initialEvalTimer;
|
|
|
initialEvalTimer.start();
|
|
|
m_totalEvaluations++;
|
|
|
m_userInitialFitness = evaluateFitness(m_userInitialSolution);
|
|
|
const int initialEvalElapsedMs = initialEvalTimer.elapsed();
|
|
|
|
|
|
if(m_userInitialFitness < 1e9) {
|
|
|
m_successfulEvaluations++;
|
|
|
m_hasValidUserSolution = true;
|
|
|
m_globalBestFitness = m_userInitialFitness;
|
|
|
m_globalBestPosition = m_userInitialSolution;
|
|
|
m_userInitialLogLogData = m_lastEvaluatedLogLogData;
|
|
|
m_globalBestLogLogData = m_userInitialLogLogData;
|
|
|
m_userInitialObjectiveBreakdown = m_lastObjectiveBreakdown;
|
|
|
m_globalBestObjectiveBreakdown = m_userInitialObjectiveBreakdown;
|
|
|
DEBUG_OUT(tr("Initial solution evaluation successful"));
|
|
|
emit logMessageGenerated(tr("Initial total error: %1").arg(lmLogError(m_userInitialFitness)));
|
|
|
emit bestCurveUpdated(m_targetLogLogData,
|
|
|
m_globalBestLogLogData,
|
|
|
0,
|
|
|
m_globalBestFitness);
|
|
|
} else {
|
|
|
m_hasValidUserSolution = false;
|
|
|
if(!m_shouldStop) emit logMessageGenerated(tr("Initial solution failed; trying the parameter-range midpoint."));
|
|
|
}
|
|
|
|
|
|
writeTraceRow(-1,
|
|
|
-1,
|
|
|
"initial_solution",
|
|
|
m_userInitialSolution,
|
|
|
m_userInitialFitness,
|
|
|
m_userInitialFitness < 1e9,
|
|
|
initialEvalElapsedMs,
|
|
|
m_hasValidUserSolution ? "valid" : "invalid",
|
|
|
m_hasValidUserSolution
|
|
|
? &m_userInitialObjectiveBreakdown : nullptr);
|
|
|
} catch(...) {
|
|
|
m_hasValidUserSolution = false;
|
|
|
if(!m_shouldStop) emit logMessageGenerated(tr("Initial solution failed; trying the parameter-range midpoint."));
|
|
|
}
|
|
|
|
|
|
m_initialValues = savedInitialValues;
|
|
|
}
|
|
|
|
|
|
finalReason = runTrustRegionFitting();
|
|
|
validateAndProtectFinalResult();
|
|
|
|
|
|
if(!m_globalBestPosition.isEmpty() && m_globalBestObjectiveBreakdown.valid) {
|
|
|
// 精英保护之后记录最终行，保证轨迹与实际写回参数一致。
|
|
|
writeTraceRow(m_currentIteration,
|
|
|
-1,
|
|
|
"trust_region_final",
|
|
|
m_globalBestPosition,
|
|
|
m_globalBestFitness,
|
|
|
m_globalBestFitness < 1.0e9,
|
|
|
-1,
|
|
|
"final_result",
|
|
|
&m_globalBestObjectiveBreakdown);
|
|
|
}
|
|
|
|
|
|
if(finalReason != LM_USER_STOPPED && finalReason != LM_OPTIMIZATION_FAILED &&
|
|
|
m_globalBestFitness < m_targetError) {
|
|
|
finalReason = LM_TARGET_ACHIEVED;
|
|
|
}
|
|
|
} catch(const std::exception& e) {
|
|
|
m_lastError = QString(tr("Critical exception in automatic fitting: %1")).arg(e.what());
|
|
|
emit logMessageGenerated(tr("CRITICAL ERROR: %1").arg(e.what()));
|
|
|
closeTraceFile();
|
|
|
cleanupTemporaryDirectory();
|
|
|
m_isFinalizing = false;
|
|
|
m_isRunning = false;
|
|
|
emitRunSummary(false, LM_OPTIMIZATION_FAILED);
|
|
|
emit fittingFinished(false, m_lastError);
|
|
|
return false;
|
|
|
} catch(...) {
|
|
|
m_lastError = tr("Unknown critical exception in automatic fitting");
|
|
|
emit logMessageGenerated(tr("CRITICAL ERROR: Unknown exception in automatic fitting"));
|
|
|
closeTraceFile();
|
|
|
cleanupTemporaryDirectory();
|
|
|
m_isFinalizing = false;
|
|
|
m_isRunning = false;
|
|
|
emitRunSummary(false, LM_OPTIMIZATION_FAILED);
|
|
|
emit fittingFinished(false, m_lastError);
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
// 只有最终完整求解实际执行并提交快照后才能置为成功。
|
|
|
bool finalFullSolverSucceeded = false;
|
|
|
bool finalFullSolverExecuted = false;
|
|
|
|
|
|
if(!m_globalBestPosition.isEmpty()) {
|
|
|
try {
|
|
|
// Stop 只结束优化迭代；从这里开始必须用当前最优参数生成并发布正式快照。
|
|
|
m_isFinalizing = true;
|
|
|
emit finalizingStarted();
|
|
|
DEBUG_OUT(tr("Applying optimized parameters to model..."));
|
|
|
applyParametersToDataManager(m_globalBestPosition);
|
|
|
|
|
|
// 裂缝参数会改变网格输入；标记失效后，最终求解任务会基于新快照重建网格。
|
|
|
const bool fractureGridParameterSelected =
|
|
|
(m_parameterSelected.size() > 4 && m_parameterSelected[4]) ||
|
|
|
(m_parameterSelected.size() > 5 && m_parameterSelected[5]);
|
|
|
if(fractureGridParameterSelected) {
|
|
|
nmDataAnalyzeManager* dataManager =
|
|
|
nmDataAnalyzeManager::getCurrentInstance();
|
|
|
if(!dataManager) {
|
|
|
throw std::runtime_error("Data manager is unavailable");
|
|
|
}
|
|
|
dataManager->invalidatePebiGrid();
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(
|
|
|
tr("Running final full-field calculation with optimized parameters..."));
|
|
|
finalFullSolverExecuted = true;
|
|
|
finalFullSolverSucceeded = runFinalFullSolver();
|
|
|
|
|
|
if(finalFullSolverSucceeded) {
|
|
|
DEBUG_OUT(
|
|
|
tr("Final full-field calculation completed successfully"));
|
|
|
} else {
|
|
|
m_lastError =
|
|
|
tr("Optimized parameters were found, but the final full-field calculation failed");
|
|
|
DEBUG_OUT(
|
|
|
tr("ERROR: Final full-field calculation failed"));
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(tr("=== Optimization Results ==="));
|
|
|
DEBUG_OUT(
|
|
|
tr("Final error: %1").arg(m_globalBestFitness, 0, 'e', 4));
|
|
|
DEBUG_OUT(
|
|
|
tr("Total iterations: %1").arg(m_currentIteration + 1));
|
|
|
DEBUG_OUT(
|
|
|
tr("Total evaluations: %1 (successful: %2)")
|
|
|
.arg(m_totalEvaluations).arg(m_successfulEvaluations));
|
|
|
|
|
|
QString finalParams = tr("Optimized parameters: ");
|
|
|
for(int i = 0; i < m_globalBestPosition.size(); ++i) {
|
|
|
finalParams += QString("[%1]=%2 ")
|
|
|
.arg(i).arg(m_globalBestPosition[i], 0, 'f', 6);
|
|
|
}
|
|
|
DEBUG_OUT(finalParams);
|
|
|
|
|
|
if(finalFullSolverExecuted && finalFullSolverSucceeded) {
|
|
|
DEBUG_OUT(
|
|
|
tr("Parameters and full-field results applied successfully to data manager"));
|
|
|
} else if(!finalFullSolverExecuted) {
|
|
|
DEBUG_OUT(
|
|
|
tr("Optimized parameters applied to data manager"));
|
|
|
}
|
|
|
} catch(const std::exception& e) {
|
|
|
finalFullSolverSucceeded = false;
|
|
|
m_lastError = tr("Failed to apply final parameters: %1").arg(e.what());
|
|
|
DEBUG_OUT(
|
|
|
tr("ERROR: Failed to apply final parameters: %1").arg(e.what()));
|
|
|
} catch(...) {
|
|
|
finalFullSolverSucceeded = false;
|
|
|
m_lastError = tr("Failed to apply final parameters due to unknown error");
|
|
|
DEBUG_OUT(
|
|
|
tr("ERROR: Unknown error applying final parameters"));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
m_isFinalizing = false;
|
|
|
m_isRunning = false;
|
|
|
bool success = false;
|
|
|
QString message;
|
|
|
|
|
|
if(finalReason == LM_TARGET_ACHIEVED) {
|
|
|
success = true;
|
|
|
message = QString(tr("Target achieved. Best error: %1, Iterations: %2"))
|
|
|
.arg(m_globalBestFitness, 0, 'e', 4)
|
|
|
.arg(m_currentIteration + 1);
|
|
|
DEBUG_OUT(tr("=== LM AUTOMATIC FITTING SUCCESSFUL ==="));
|
|
|
} else if(finalReason == LM_LOCAL_OPTIMUM) {
|
|
|
success = true;
|
|
|
message = QString(
|
|
|
tr("Automatic fitting reached a local optimum. Best error: %1, Iterations: %2"))
|
|
|
.arg(m_globalBestFitness, 0, 'e', 4)
|
|
|
.arg(m_currentIteration + 1);
|
|
|
DEBUG_OUT(tr("=== LM AUTOMATIC FITTING - LOCAL OPTIMUM ==="));
|
|
|
} else if(finalReason == LM_MAX_ITERATIONS) {
|
|
|
success = true;
|
|
|
message = QString(tr("Iteration or evaluation budget reached. Best error: %1, Cumulative iterations: %2"))
|
|
|
.arg(m_globalBestFitness, 0, 'e', 4)
|
|
|
.arg(m_currentIteration + 1);
|
|
|
DEBUG_OUT(tr("=== LM AUTOMATIC FITTING - MAX ITERATIONS ==="));
|
|
|
} else if(finalReason == LM_USER_STOPPED) {
|
|
|
success = true;
|
|
|
message = QString(tr("Stopped by user. Best error: %1, Iterations: %2"))
|
|
|
.arg(m_globalBestFitness, 0, 'e', 4)
|
|
|
.arg(m_currentIteration + 1);
|
|
|
DEBUG_OUT(tr("=== LM AUTOMATIC FITTING STOPPED BY USER ==="));
|
|
|
} else if(finalReason == LM_CONSECUTIVE_FAILURES) {
|
|
|
message = QString(
|
|
|
tr("Automatic fitting failed due to consecutive failures. Best error: %1, Iterations: %2"))
|
|
|
.arg(m_globalBestFitness, 0, 'e', 4)
|
|
|
.arg(m_currentIteration + 1);
|
|
|
DEBUG_OUT(tr("=== LM AUTOMATIC FITTING FAILED ==="));
|
|
|
} else {
|
|
|
message = QString(
|
|
|
tr("Automatic fitting ended unexpectedly. Best error: %1, Iterations: %2"))
|
|
|
.arg(m_globalBestFitness, 0, 'e', 4)
|
|
|
.arg(m_currentIteration + 1);
|
|
|
DEBUG_OUT(tr("=== LM AUTOMATIC FITTING - UNKNOWN END ==="));
|
|
|
}
|
|
|
|
|
|
if(!finalFullSolverSucceeded) {
|
|
|
success = false;
|
|
|
message = m_lastError;
|
|
|
}
|
|
|
|
|
|
emitRunSummary(success, finalReason);
|
|
|
QApplication::processEvents();
|
|
|
msleep(200);
|
|
|
QApplication::processEvents();
|
|
|
closeTraceFile();
|
|
|
cleanupTemporaryDirectory();
|
|
|
emit fittingFinished(success, message);
|
|
|
return success;
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::extractUserInitialValues()
|
|
|
{
|
|
|
// 从当前项目模型读取用户已有初始参数。
|
|
|
// 只提取用户勾选的参数，并按 m_enabledParamIndices 的顺序写入 m_initialValues。
|
|
|
// 这些值用于初始解真实评价、LM 起点和最终精英保护。
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
|
|
|
|
|
|
nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
|
|
|
|
|
|
m_initialValues.clear();
|
|
|
|
|
|
// 按照启用参数的顺序提取初始值。井参数来自目标井，储层参数来自 reservoirData。
|
|
|
for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
|
|
|
int paramIndex = m_enabledParamIndices[i];
|
|
|
double initialValue = 0.0;
|
|
|
|
|
|
switch(paramIndex) {
|
|
|
case 0: // 渗透率
|
|
|
initialValue = reservoirData.getPermeability().getValue().toDouble();
|
|
|
break;
|
|
|
|
|
|
case 1: // 表皮系数
|
|
|
if(pTargetWell) {
|
|
|
initialValue = pTargetWell->getPerforation(0)->getSkin().getValue().toDouble();
|
|
|
}
|
|
|
|
|
|
break;
|
|
|
|
|
|
case 2: // 井筒储集系数
|
|
|
if(pTargetWell) {
|
|
|
initialValue = pTargetWell->getWellboreStorage().getValue().toDouble();
|
|
|
}
|
|
|
|
|
|
break;
|
|
|
|
|
|
case 3: // 孔隙度
|
|
|
initialValue = reservoirData.getPorosity().getValue().toDouble();
|
|
|
break;
|
|
|
|
|
|
case 4: // 裂缝导流能力
|
|
|
if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
|
|
|
nmDataVerticalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataVerticalFracturedWell*>(pTargetWell);
|
|
|
if(fracturedWell) {
|
|
|
initialValue = fracturedWell->getDfc().getValue().toDouble();
|
|
|
}
|
|
|
} else if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
|
|
|
nmDataHorizontalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataHorizontalFracturedWell*>(pTargetWell);
|
|
|
if(fracturedWell) {
|
|
|
initialValue = fracturedWell->getDfc().getValue().toDouble();
|
|
|
}
|
|
|
}
|
|
|
break;
|
|
|
|
|
|
case 5: // 裂缝半长
|
|
|
if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
|
|
|
nmDataVerticalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataVerticalFracturedWell*>(pTargetWell);
|
|
|
if(fracturedWell) {
|
|
|
initialValue = fracturedWell->getFractureHalfLength().getValue().toDouble();
|
|
|
}
|
|
|
} else if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
|
|
|
nmDataHorizontalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataHorizontalFracturedWell*>(pTargetWell);
|
|
|
if(fracturedWell) {
|
|
|
initialValue = fracturedWell->getFractureHalfLength().getValue().toDouble();
|
|
|
}
|
|
|
}
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
m_initialValues.append(initialValue);
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(QString("Extracted %1 user initial values").arg(m_initialValues.size()));
|
|
|
|
|
|
for(int i = 0; i < m_initialValues.size(); ++i) {
|
|
|
DEBUG_OUT(QString(" Initial[%1] = %2").arg(i).arg(m_initialValues[i], 0, 'e', 3));
|
|
|
}
|
|
|
|
|
|
if(!validateInitialValues()) {
|
|
|
DEBUG_OUT("Warning: Some initial values are outside parameter bounds");
|
|
|
}
|
|
|
}
|
|
|
|
|
|
bool nmCalculationAutoFitLM::evaluateTrustRegionPoint(
|
|
|
const QVector<double>& parameters,
|
|
|
double* fitness,
|
|
|
AutoFitObjectiveBreakdownLM* breakdown,
|
|
|
QVector<QVector<double> >* curve,
|
|
|
int* elapsedMs)
|
|
|
{
|
|
|
if(!fitness || !breakdown || !curve || !elapsedMs || m_shouldStop) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
// evaluateFitness() 会写入 DataManager 并调用真实求解器。这里统一统计
|
|
|
// 真实评价次数和耗时，同时要求固定采样残差、误差结构和结果曲线均有效。
|
|
|
QTime timer;
|
|
|
timer.start();
|
|
|
*fitness = evaluateFitness(parameters, false);
|
|
|
*elapsedMs = timer.elapsed();
|
|
|
*breakdown = m_lastObjectiveBreakdown;
|
|
|
*curve = m_lastEvaluatedLogLogData;
|
|
|
++m_totalEvaluations;
|
|
|
|
|
|
bool valid = isFiniteNumber(*fitness) && *fitness < 1.0e9 &&
|
|
|
breakdown->valid &&
|
|
|
trustRegionResidualsValid(*breakdown) &&
|
|
|
!curve->isEmpty();
|
|
|
if(valid) {
|
|
|
++m_successfulEvaluations;
|
|
|
}
|
|
|
|
|
|
return valid;
|
|
|
}
|
|
|
|
|
|
StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
|
|
|
{
|
|
|
// 本次拟合开始时按目标井的实际井型确定流程，普通直井跳过形状调整。
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
nmDataWellBase* targetWell = dataManager ? dataManager->findWellByName(m_targetWellName) : nullptr;
|
|
|
const bool skipShapeStage = targetWell && targetWell->getWellType() == NM_WELL_MODEL::Vertical_Well;
|
|
|
const int dimensions = getEnabledParameterCount();
|
|
|
if(dimensions <= 0 || m_enabledParamIndices.size() != dimensions) {
|
|
|
m_lastError = tr("No valid parameters are available for trust-region fitting");
|
|
|
return LM_OPTIMIZATION_FAILED;
|
|
|
}
|
|
|
|
|
|
// 高度、形状、总误差共用迭代上限；两个 LM 目标分别设置求解评价额度。
|
|
|
const int phaseEvaluationBudget = qMax(dimensions + 2, qMax(20, m_maxIterations * 3));
|
|
|
int maximumEvaluations = 0;
|
|
|
int phaseIterations = 0;
|
|
|
int shapeIterations = 0;
|
|
|
int totalIterations = 0;
|
|
|
int completedIterations = 0;
|
|
|
// 步长使用归一化内部坐标：0.04 表示内部范围的 4%，并非物理值变化 4%。
|
|
|
// 信赖半径限制所有参数一次联合移动的二范数。
|
|
|
const double sensitivityStep = 0.04;
|
|
|
const double minimumCoordinateStep = 1.0e-5;
|
|
|
const double minimumTrustRadius = 2.0e-3;
|
|
|
const double maximumTrustRadius = 0.30;
|
|
|
// 误差下降至少达到绝对 0.01 且相对当前有效基准 10% 才算有效改善。
|
|
|
// 更小的下降仍保留为最佳解，但不能反复清除停滞状态、延长拟合时间。
|
|
|
const double effectiveRelativeImprovement = 1.0e-1;
|
|
|
const double effectiveAbsoluteImprovement = 1.0e-2;
|
|
|
const int maximumIneffectiveSteps = 3;
|
|
|
|
|
|
// damping 是 LM 阻尼；拒绝或预测失准时增大，真实下降与预测一致时减小。
|
|
|
// 按连续预测失准刷新 Jacobian；接受步数只作兜底，不因累计移动距离强制重建。
|
|
|
double trustRadius = 0.12;
|
|
|
double damping = 1.0e-2;
|
|
|
int consecutiveRejectedSteps = 0;
|
|
|
int consecutiveSolverFailures = 0;
|
|
|
int acceptedSinceRebuild = 0;
|
|
|
int consecutiveIneffectiveSteps = 0;
|
|
|
int consecutivePoorPredictions = 0;
|
|
|
bool rebuildRequested = true;
|
|
|
// 跟踪文件保留稳定英文原因，界面输出时再翻译，避免受本地编码影响。
|
|
|
const char* rebuildReason = QT_TR_NOOP("initial sensitivity model");
|
|
|
bool modelRebuiltAtMinimumRadius = false;
|
|
|
bool stagnationConfirmationRequested = false;
|
|
|
StopReasonLM stopReason = LM_MAX_ITERATIONS;
|
|
|
|
|
|
// jacobian 的行对应固定采样的残差，列对应用户勾选的参数。
|
|
|
// LM 按当前目标取数值或形状残差行，两段均使用所有勾选参数。
|
|
|
QVector<QVector<double> > jacobian;
|
|
|
QVector<bool> jacobianColumnValid(dimensions, false);
|
|
|
|
|
|
// 参数向量的顺序始终与 m_enabledParamIndices 一致，不能按完整参数索引
|
|
|
// 直接访问；下面两个转换函数集中维护这层映射关系。
|
|
|
auto coordinatesFromParameters = [&](const QVector<double>& parameters)
|
|
|
-> QVector<double> {
|
|
|
QVector<double> coordinates(dimensions, 0.0);
|
|
|
for(int i = 0; i < dimensions; ++i) {
|
|
|
int parameterIndex = m_enabledParamIndices[i];
|
|
|
coordinates[i] = toTrustRegionCoordinate(
|
|
|
parameters[i], parameterIndex,
|
|
|
m_parameterLower[parameterIndex],
|
|
|
m_parameterUpper[parameterIndex]);
|
|
|
}
|
|
|
return coordinates;
|
|
|
};
|
|
|
|
|
|
auto parametersFromCoordinates = [&](const QVector<double>& coordinates)
|
|
|
-> QVector<double> {
|
|
|
QVector<double> parameters(dimensions, 0.0);
|
|
|
for(int i = 0; i < dimensions; ++i) {
|
|
|
int parameterIndex = m_enabledParamIndices[i];
|
|
|
parameters[i] = fromTrustRegionCoordinate(
|
|
|
coordinates[i], parameterIndex,
|
|
|
m_parameterLower[parameterIndex],
|
|
|
m_parameterUpper[parameterIndex]);
|
|
|
}
|
|
|
return parameters;
|
|
|
};
|
|
|
|
|
|
auto restoreEvaluationState = [&](const TrustRegionEvaluation& evaluation) {
|
|
|
// evaluateFitness() 会把试算参数写入 DataManager。无论候选是否接受，
|
|
|
// 下一次计算前都恢复到唯一的已接受工作点，防止失败试算污染后续求解。
|
|
|
applyParametersToDataManager(evaluation.parameters);
|
|
|
m_lastObjectiveBreakdown = evaluation.breakdown;
|
|
|
m_lastEvaluatedLogLogData = evaluation.curve;
|
|
|
};
|
|
|
|
|
|
// 各阶段只发布通过当前目标和约束检查的工作点；曲线和诊断快照必须
|
|
|
// 与参数同步更新，防止界面显示或最终精英保护使用错配的数据。
|
|
|
auto publishAcceptedPoint = [&](const TrustRegionEvaluation& evaluation) {
|
|
|
m_globalBestPosition = evaluation.parameters;
|
|
|
m_globalBestFitness = evaluation.fitness;
|
|
|
m_globalBestObjectiveBreakdown = evaluation.breakdown;
|
|
|
m_globalBestLogLogData = evaluation.curve;
|
|
|
emit bestCurveUpdated(m_targetLogLogData,
|
|
|
m_globalBestLogLogData,
|
|
|
m_currentIteration + 1,
|
|
|
m_globalBestFitness);
|
|
|
};
|
|
|
|
|
|
auto processPauseAndStop = [&]() -> bool {
|
|
|
QApplication::processEvents();
|
|
|
return !m_shouldStop;
|
|
|
};
|
|
|
|
|
|
// current 始终代表唯一已接受工作点。优先复用启动阶段已经真实验证的
|
|
|
// 用户初始解，避免在信赖域入口重复调用一次昂贵求解器。
|
|
|
TrustRegionEvaluation current;
|
|
|
if(m_hasValidUserSolution &&
|
|
|
m_globalBestPosition.size() == dimensions &&
|
|
|
trustRegionResidualsValid(m_globalBestObjectiveBreakdown) &&
|
|
|
!m_globalBestLogLogData.isEmpty()) {
|
|
|
current.parameters = m_globalBestPosition;
|
|
|
current.coordinates = coordinatesFromParameters(current.parameters);
|
|
|
current.breakdown = m_globalBestObjectiveBreakdown;
|
|
|
current.curve = m_globalBestLogLogData;
|
|
|
current.fitness = m_globalBestFitness;
|
|
|
current.elapsedMs = 0;
|
|
|
current.valid = true;
|
|
|
} else {
|
|
|
// 用户初始解无效时只做一次确定性的范围中点回退；所有正值参数在对数
|
|
|
// 坐标取中点，避免线性中点过分偏向跨数量级范围的上界。
|
|
|
current.coordinates.fill(0.5, dimensions);
|
|
|
current.parameters = parametersFromCoordinates(current.coordinates);
|
|
|
current.valid = evaluateTrustRegionPoint(
|
|
|
current.parameters,
|
|
|
¤t.fitness,
|
|
|
¤t.breakdown,
|
|
|
¤t.curve,
|
|
|
¤t.elapsedMs);
|
|
|
writeTraceRow(-1, -1,
|
|
|
"trust_region_midpoint",
|
|
|
current.parameters,
|
|
|
current.fitness,
|
|
|
current.valid,
|
|
|
current.elapsedMs,
|
|
|
current.valid ? "midpoint_valid" : "midpoint_invalid",
|
|
|
current.valid ? ¤t.breakdown : nullptr);
|
|
|
if(!current.valid) {
|
|
|
m_lastError = tr("The initial solution and parameter-range midpoint are both invalid");
|
|
|
return m_shouldStop
|
|
|
? LM_USER_STOPPED
|
|
|
: LM_OPTIMIZATION_FAILED;
|
|
|
}
|
|
|
publishAcceptedPoint(current);
|
|
|
}
|
|
|
|
|
|
restoreEvaluationState(current);
|
|
|
emit progressUpdated(0, current.fitness);
|
|
|
DEBUG_OUT(tr("=== Starting LM Main Loop ==="));
|
|
|
DEBUG_OUT(
|
|
|
tr("LM starting point error: %1; evaluation budget per LM phase: %2")
|
|
|
.arg(current.fitness, 0, 'e', 4)
|
|
|
.arg(phaseEvaluationBudget));
|
|
|
|
|
|
// 高度预调整只改变用户勾选的渗透率。ln(k) 的初始变化由有符号高度差
|
|
|
// 给出，真实求解后用割线估计修正；拒绝时缩步，不让其他参数补偿高度。
|
|
|
const int permeabilityColumn = m_enabledParamIndices.indexOf(0);
|
|
|
double heightBaseline = current.breakdown.verticalLoss;
|
|
|
int heightStagnation = 0;
|
|
|
double heightScale = 1.0;
|
|
|
double heightSlope = -1.0;
|
|
|
int heightTrials = 0;
|
|
|
if(permeabilityColumn >= 0) {
|
|
|
emit fittingStageChanged(1);
|
|
|
emit logMessageGenerated(tr("[Stage 1: Height alignment]"));
|
|
|
// 高度调整计入总迭代数；对齐、停滞或无可行步长时结束。
|
|
|
while(completedIterations < m_maxIterations && processPauseAndStop()) {
|
|
|
const double bias = current.breakdown.verticalCommonBias;
|
|
|
if(!current.breakdown.verticalReliable || qAbs(bias) <= 0.1 || heightStagnation >= 2) break;
|
|
|
const double k = current.parameters[permeabilityColumn];
|
|
|
if(k <= 0.0 || m_parameterLower[0] <= 0.0 || m_parameterUpper[0] <= m_parameterLower[0]) break;
|
|
|
const double logRange = qLn(m_parameterUpper[0]) - qLn(m_parameterLower[0]);
|
|
|
const double change = qBound(-0.30 * logRange, -bias / heightSlope * heightScale, 0.30 * logRange);
|
|
|
const AutoFitObjectiveBreakdownLM previous = current.breakdown;
|
|
|
TrustRegionEvaluation candidate;
|
|
|
candidate.parameters = current.parameters;
|
|
|
candidate.parameters[permeabilityColumn] = qBound(m_parameterLower[0], k * qExp(change), m_parameterUpper[0]);
|
|
|
candidate.coordinates = coordinatesFromParameters(candidate.parameters);
|
|
|
const double actualStep = qLn(candidate.parameters[permeabilityColumn] / k);
|
|
|
if(qAbs(actualStep) < 1.0e-5) break;
|
|
|
m_currentIteration = completedIterations++;
|
|
|
++heightTrials;
|
|
|
emit progressUpdated(completedIterations, current.fitness);
|
|
|
candidate.valid = evaluateTrustRegionPoint(candidate.parameters, &candidate.fitness,
|
|
|
&candidate.breakdown, &candidate.curve, &candidate.elapsedMs);
|
|
|
// 第一阶段只按上下偏移改善验收，不限制形状；仍防止两条曲线一高一低互相抵消。
|
|
|
const bool accepted = candidate.valid && candidate.breakdown.verticalReliable &&
|
|
|
candidate.breakdown.verticalLoss < current.breakdown.verticalLoss;
|
|
|
if(candidate.valid) {
|
|
|
const double slope = (candidate.breakdown.verticalCommonBias - bias) / actualStep;
|
|
|
if(slope < -0.05 && isFiniteNumber(slope)) heightSlope = slope;
|
|
|
}
|
|
|
writeTraceRow(m_currentIteration, 0, "permeability_height", candidate.parameters, candidate.fitness,
|
|
|
candidate.valid, candidate.elapsedMs, accepted ? "accepted_height" : "rejected_height",
|
|
|
candidate.valid ? &candidate.breakdown : nullptr);
|
|
|
if(accepted) {
|
|
|
current = candidate;
|
|
|
publishAcceptedPoint(current);
|
|
|
} else {
|
|
|
heightScale *= 0.5;
|
|
|
}
|
|
|
restoreEvaluationState(current);
|
|
|
if(heightBaseline - current.breakdown.verticalLoss >= qMax(1.0e-4, 0.10 * heightBaseline)) {
|
|
|
heightBaseline = current.breakdown.verticalLoss;
|
|
|
heightStagnation = 0;
|
|
|
} else ++heightStagnation;
|
|
|
if(candidate.valid) {
|
|
|
const QString objective = tr("Height error: %1 -> %2")
|
|
|
.arg(lmLogError(previous.verticalLoss)).arg(lmLogError(candidate.breakdown.verticalLoss));
|
|
|
emit logMessageGenerated(tr("%1 | %2 | %3 | Improvement: %4%")
|
|
|
.arg(tr("Iteration %1").arg(completedIterations))
|
|
|
.arg(accepted ? tr("accepted") : tr("rejected")).arg(objective)
|
|
|
.arg(lmLogImprovement(previous.verticalLoss, candidate.breakdown.verticalLoss)));
|
|
|
if(!accepted) emit logMessageGenerated(candidate.breakdown.verticalReliable
|
|
|
? tr("Height error did not decrease; keep the last accepted result and reduce the step.")
|
|
|
: tr("Pressure and derivative height directions conflict; keep the last accepted result and reduce the step."));
|
|
|
} else {
|
|
|
emit logMessageGenerated(tr("%1 | Solver failed\nNo valid error was returned; keep the last accepted result and reduce the step.")
|
|
|
.arg(tr("Iteration %1").arg(completedIterations)));
|
|
|
}
|
|
|
emit progressUpdated(completedIterations, current.fitness);
|
|
|
}
|
|
|
QString heightReason = tr("no feasible permeability step or step too small");
|
|
|
if(m_shouldStop) heightReason = tr("stopped by user");
|
|
|
else if(completedIterations >= m_maxIterations) heightReason = getStopReasonDescription(LM_MAX_ITERATIONS);
|
|
|
else if(!current.breakdown.verticalReliable) heightReason = tr("pressure and derivative height directions conflict");
|
|
|
else if(current.breakdown.verticalLoss <= 0.1) heightReason = tr("curve height is approximately aligned");
|
|
|
else if(heightStagnation >= 2) heightReason = tr("2 consecutive steps without effective height improvement");
|
|
|
else if(current.parameters[permeabilityColumn] <= m_parameterLower[0] ||
|
|
|
current.parameters[permeabilityColumn] >= m_parameterUpper[0]) heightReason = tr("permeability reached its bound");
|
|
|
emit logMessageGenerated(tr("Height alignment ended: %1").arg(heightReason));
|
|
|
}
|
|
|
|
|
|
if(permeabilityColumn < 0) emit logMessageGenerated(tr("Height alignment skipped: permeability is not selected."));
|
|
|
if(m_shouldStop) return LM_USER_STOPPED;
|
|
|
if(completedIterations >= m_maxIterations) return LM_MAX_ITERATIONS;
|
|
|
|
|
|
// 直井直接匹配总误差，其他井型仍先匹配整体形状，再匹配总误差。
|
|
|
bool shapeStage = !skipShapeStage;
|
|
|
auto stageError = [&](const TrustRegionEvaluation& point) -> double {
|
|
|
return shapeStage ? point.breakdown.shapeLoss : point.fitness;
|
|
|
};
|
|
|
auto acceptable = [&](const TrustRegionEvaluation& point, const TrustRegionEvaluation& base) -> bool {
|
|
|
return point.valid && stageError(point) < stageError(base);
|
|
|
};
|
|
|
auto acceptPoint = [&](const TrustRegionEvaluation& point) {
|
|
|
// 同步发布当前目标下已接受的参数、曲线和诊断。
|
|
|
current = point;
|
|
|
publishAcceptedPoint(current);
|
|
|
restoreEvaluationState(current);
|
|
|
};
|
|
|
double effectiveImprovementBaseline = stageError(current);
|
|
|
DEBUG_OUT(tr("LM sampling: %1 points (%2 intervals per log-time decade).")
|
|
|
.arg(current.breakdown.residualVector.size() / 2).arg(kAutoFitIntervalsPerDecade));
|
|
|
|
|
|
auto registerEffectiveImprovement = [&](double fitness) -> bool {
|
|
|
const double requiredImprovement = qMax(
|
|
|
effectiveAbsoluteImprovement,
|
|
|
qAbs(effectiveImprovementBaseline) * effectiveRelativeImprovement);
|
|
|
const double improvement = effectiveImprovementBaseline - fitness;
|
|
|
if(improvement < requiredImprovement) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
// 有效改善后更新累计基准，并清除停滞计数。
|
|
|
effectiveImprovementBaseline = fitness;
|
|
|
consecutiveIneffectiveSteps = 0;
|
|
|
stagnationConfirmationRequested = false;
|
|
|
return true;
|
|
|
};
|
|
|
|
|
|
// 两个 LM 目标共用有效改善门槛和重建后的停滞确认规则。
|
|
|
auto recordIneffectiveStep = [&]() -> bool {
|
|
|
++consecutiveIneffectiveSteps;
|
|
|
if(consecutiveIneffectiveSteps < maximumIneffectiveSteps) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(stagnationConfirmationRequested) {
|
|
|
return true;
|
|
|
}
|
|
|
|
|
|
stagnationConfirmationRequested = true;
|
|
|
rebuildRequested = true;
|
|
|
rebuildReason = QT_TR_NOOP("confirm stagnation after ineffective steps");
|
|
|
emit logMessageGenerated(
|
|
|
tr("No effective improvement for %1 consecutive steps; "
|
|
|
"rebuilding sensitivity model for confirmation")
|
|
|
.arg(consecutiveIneffectiveSteps));
|
|
|
consecutiveIneffectiveSteps = 0;
|
|
|
return false;
|
|
|
};
|
|
|
|
|
|
DEBUG_OUT(
|
|
|
tr("LM effective improvement threshold: max(%1, %2% of baseline error); "
|
|
|
"%3 consecutive ineffective steps trigger convergence confirmation")
|
|
|
.arg(effectiveAbsoluteImprovement, 0, 'e', 2)
|
|
|
.arg(effectiveRelativeImprovement * 100.0, 0, 'f', 2)
|
|
|
.arg(maximumIneffectiveSteps));
|
|
|
|
|
|
// 在同一个真实工作点逐参数做单边差分。失败先在原方向缩步，再反向尝试；
|
|
|
// 正常情况下每列仍只需一次真实求解，重试也计入总评价预算。
|
|
|
auto rebuildSensitivity = [&]() -> bool {
|
|
|
// 灵敏度试算计入当前目标的评价额度，形状和总误差使用相同差分规则。
|
|
|
const int evaluationLimit = maximumEvaluations;
|
|
|
const TrustRegionEvaluation base = current;
|
|
|
const QVector<double> baseResidual = trustRegionFullResidual(base.breakdown);
|
|
|
const int residualCount = baseResidual.size();
|
|
|
if(residualCount <= 0) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
jacobian = QVector<QVector<double> >(
|
|
|
residualCount, QVector<double>(dimensions, 0.0));
|
|
|
jacobianColumnValid.fill(false, dimensions);
|
|
|
|
|
|
TrustRegionEvaluation bestProbe;
|
|
|
int bestProbeColumn = -1;
|
|
|
double bestProbeDelta = 0.0;
|
|
|
// 两个目标共用局部差分尺度，信赖半径缩小时同步减小扰动。
|
|
|
const double finiteDifferenceStep = qMin(
|
|
|
sensitivityStep, qMax(5.0e-3, trustRadius * 0.5));
|
|
|
|
|
|
for(int column = 0;
|
|
|
column < dimensions &&
|
|
|
m_totalEvaluations < evaluationLimit &&
|
|
|
processPauseAndStop();
|
|
|
++column) {
|
|
|
const int parameterIndex = m_enabledParamIndices[column];
|
|
|
|
|
|
const double lower = m_parameterLower[parameterIndex];
|
|
|
const double upper = m_parameterUpper[parameterIndex];
|
|
|
if(upper <= lower) continue;
|
|
|
// 表皮包含零和负值，按当前物理值的局部尺度确定扰动。
|
|
|
double localStep = finiteDifferenceStep;
|
|
|
if(parameterIndex == 1) {
|
|
|
localStep = qMin(localStep, 0.02 * qMax(0.1, qAbs(base.parameters[column])) / (upper - lower));
|
|
|
}
|
|
|
double positiveRoom = 1.0 - base.coordinates[column];
|
|
|
double negativeRoom = base.coordinates[column];
|
|
|
double preferredSign = positiveRoom >= negativeRoom ? 1.0 : -1.0;
|
|
|
bool columnBuilt = false;
|
|
|
double probeStep = localStep;
|
|
|
|
|
|
for(int directionAttempt = 0;
|
|
|
directionAttempt < 2 &&
|
|
|
!columnBuilt &&
|
|
|
m_totalEvaluations < evaluationLimit && processPauseAndStop();
|
|
|
++directionAttempt) {
|
|
|
double direction = directionAttempt == 0
|
|
|
? preferredSign : -preferredSign;
|
|
|
double availableRoom = direction > 0.0
|
|
|
? positiveRoom : negativeRoom;
|
|
|
double deltaMagnitude = qMin(probeStep, availableRoom);
|
|
|
for(int shrinkAttempt = 0; shrinkAttempt < 3 && !columnBuilt &&
|
|
|
m_totalEvaluations < evaluationLimit && processPauseAndStop(); ++shrinkAttempt) {
|
|
|
if(deltaMagnitude < minimumCoordinateStep) {
|
|
|
break;
|
|
|
}
|
|
|
const bool canShrink = shrinkAttempt < 2 && deltaMagnitude * 0.5 >= minimumCoordinateStep;
|
|
|
|
|
|
TrustRegionEvaluation probe;
|
|
|
probe.coordinates = base.coordinates;
|
|
|
probe.coordinates[column] += direction * deltaMagnitude;
|
|
|
probe.parameters = parametersFromCoordinates(probe.coordinates);
|
|
|
probe.valid = evaluateTrustRegionPoint(
|
|
|
probe.parameters,
|
|
|
&probe.fitness,
|
|
|
&probe.breakdown,
|
|
|
&probe.curve,
|
|
|
&probe.elapsedMs);
|
|
|
|
|
|
QString decision = probe.valid
|
|
|
? "sensitivity_valid"
|
|
|
: (canShrink ? "sensitivity_retry_smaller" : (directionAttempt == 0
|
|
|
? "sensitivity_retry_opposite"
|
|
|
: "sensitivity_invalid"));
|
|
|
writeTraceRow(m_currentIteration,
|
|
|
column,
|
|
|
"trust_region_sensitivity",
|
|
|
probe.parameters,
|
|
|
probe.fitness,
|
|
|
probe.valid,
|
|
|
probe.elapsedMs,
|
|
|
(shapeStage ? "shape_" : "total_") + decision,
|
|
|
probe.valid ? &probe.breakdown : nullptr);
|
|
|
|
|
|
if(!probe.valid) {
|
|
|
restoreEvaluationState(base);
|
|
|
if(!canShrink) break;
|
|
|
deltaMagnitude *= 0.5;
|
|
|
probeStep = deltaMagnitude;
|
|
|
continue;
|
|
|
}
|
|
|
|
|
|
double delta = probe.coordinates[column] -
|
|
|
base.coordinates[column];
|
|
|
if(qAbs(delta) < minimumCoordinateStep ||
|
|
|
trustRegionFullResidual(probe.breakdown).size() != residualCount) {
|
|
|
restoreEvaluationState(base);
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
// 第 column 列是固定残差向量相对内部参数坐标的有限差分：
|
|
|
// J[:,column] = (r_probe-r_base)/delta。
|
|
|
const QVector<double> probeResidual = trustRegionFullResidual(probe.breakdown);
|
|
|
for(int row = 0; row < residualCount; ++row) {
|
|
|
jacobian[row][column] = (probeResidual[row] - baseResidual[row]) / delta;
|
|
|
}
|
|
|
|
|
|
jacobianColumnValid[column] = true;
|
|
|
columnBuilt = true;
|
|
|
|
|
|
// 缓存当前目标下更优的探针，全部灵敏度列建完后再接受。
|
|
|
if(acceptable(probe, base) &&
|
|
|
(!bestProbe.valid || stageError(probe) < stageError(bestProbe))) {
|
|
|
bestProbe = probe;
|
|
|
bestProbeColumn = column;
|
|
|
bestProbeDelta = delta;
|
|
|
}
|
|
|
restoreEvaluationState(base);
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
|
|
|
int validColumnCount = 0;
|
|
|
for(int i = 0; i < jacobianColumnValid.size(); ++i) {
|
|
|
if(jacobianColumnValid[i]) {
|
|
|
++validColumnCount;
|
|
|
}
|
|
|
}
|
|
|
// 至少需要一列有效灵敏度；零梯度由后续可行方向和停滞判断处理。
|
|
|
if(validColumnCount == 0 || m_shouldStop) {
|
|
|
restoreEvaluationState(base);
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
// 所有列基于同一工作点建立，再接受更优探针并用割线更新 Jacobian。
|
|
|
if(bestProbe.valid && bestProbeColumn >= 0) {
|
|
|
QVector<double> acceptedStep(dimensions, 0.0);
|
|
|
acceptedStep[bestProbeColumn] = bestProbeDelta;
|
|
|
updateTrustRegionJacobian(
|
|
|
&jacobian,
|
|
|
trustRegionFullResidual(base.breakdown),
|
|
|
trustRegionFullResidual(bestProbe.breakdown),
|
|
|
acceptedStep);
|
|
|
acceptPoint(bestProbe);
|
|
|
writeTraceRow(m_currentIteration,
|
|
|
bestProbeColumn,
|
|
|
"trust_region_sensitivity_accept",
|
|
|
current.parameters,
|
|
|
current.fitness,
|
|
|
true,
|
|
|
0,
|
|
|
shapeStage ? "shape_accepted_cached_probe" : "total_accepted_cached_probe",
|
|
|
¤t.breakdown);
|
|
|
emitTrialLog(tr("Sensitivity trial"), base.breakdown, current.breakdown, true, true, shapeStage);
|
|
|
} else {
|
|
|
restoreEvaluationState(current);
|
|
|
}
|
|
|
|
|
|
// 有更优探针时以接受后的误差为新基准，否则仍为重建起点误差。
|
|
|
// 重新计数但保留停滞确认状态，确认阶段连续三次改善不足仍会结束。
|
|
|
effectiveImprovementBaseline = stageError(current);
|
|
|
consecutiveIneffectiveSteps = 0;
|
|
|
acceptedSinceRebuild = 0;
|
|
|
consecutivePoorPredictions = 0;
|
|
|
consecutiveRejectedSteps = 0;
|
|
|
rebuildRequested = false;
|
|
|
// 若重建过程中接受了试算点，当前模型已通过割线平移而不是在新点完整
|
|
|
// 重算；再遇到最小半径停滞时仍允许做一次真正的新点重建。
|
|
|
modelRebuiltAtMinimumRadius =
|
|
|
trustRadius <= minimumTrustRadius * 1.01 &&
|
|
|
!bestProbe.valid;
|
|
|
emit logMessageGenerated(
|
|
|
tr("Sensitivity model rebuilt: %1/%2 parameter columns valid")
|
|
|
.arg(validColumnCount)
|
|
|
.arg(dimensions));
|
|
|
return true;
|
|
|
};
|
|
|
|
|
|
for(int phase = skipShapeStage ? 1 : 0;
|
|
|
phase < 2 && !m_shouldStop && completedIterations < m_maxIterations; ++phase) {
|
|
|
shapeStage = phase == 0;
|
|
|
phaseIterations = 0;
|
|
|
maximumEvaluations = m_totalEvaluations + phaseEvaluationBudget;
|
|
|
stopReason = LM_MAX_ITERATIONS;
|
|
|
effectiveImprovementBaseline = stageError(current);
|
|
|
trustRadius = 0.12;
|
|
|
damping = 0.01;
|
|
|
consecutiveRejectedSteps = 0;
|
|
|
consecutiveSolverFailures = 0;
|
|
|
acceptedSinceRebuild = 0;
|
|
|
consecutiveIneffectiveSteps = 0;
|
|
|
consecutivePoorPredictions = 0;
|
|
|
modelRebuiltAtMinimumRadius = false;
|
|
|
stagnationConfirmationRequested = false;
|
|
|
jacobian.clear();
|
|
|
rebuildRequested = true;
|
|
|
rebuildReason = shapeStage ? QT_TR_NOOP("full sensitivity at shape LM entry")
|
|
|
: QT_TR_NOOP("full sensitivity at total-stage entry");
|
|
|
emit fittingStageChanged(shapeStage ? 2 : 3);
|
|
|
emit progressUpdated(completedIterations, current.fitness);
|
|
|
emit logMessageGenerated(shapeStage
|
|
|
? tr("[Stage 2: Joint LM - Shape]")
|
|
|
: tr("[Stage 2: Joint LM - Total error]"));
|
|
|
DEBUG_OUT(tr("LM phase budget: %1 remaining shared iterations, %2 evaluations.")
|
|
|
.arg(m_maxIterations - completedIterations).arg(phaseEvaluationBudget));
|
|
|
writeTraceRow(completedIterations, -1, "stage_switch", current.parameters, current.fitness,
|
|
|
true, 0, shapeStage ? "height_to_shape_lm" : (skipShapeStage ? "height_to_total" : "shape_to_total"),
|
|
|
¤t.breakdown);
|
|
|
|
|
|
while(!m_shouldStop && completedIterations < m_maxIterations && m_totalEvaluations < maximumEvaluations) {
|
|
|
m_currentIteration = completedIterations;
|
|
|
if(!processPauseAndStop()) break;
|
|
|
if(stageError(current) < m_targetError) {
|
|
|
stopReason = LM_TARGET_ACHIEVED;
|
|
|
break;
|
|
|
}
|
|
|
if(rebuildRequested) {
|
|
|
emit logMessageGenerated(tr("Calculating sensitivity..."));
|
|
|
writeTraceRow(m_currentIteration, -1, "sensitivity_rebuild", current.parameters,
|
|
|
current.fitness, true, 0, rebuildReason, ¤t.breakdown);
|
|
|
if(!rebuildSensitivity()) {
|
|
|
stopReason = m_shouldStop
|
|
|
? LM_USER_STOPPED
|
|
|
: LM_LOCAL_OPTIMUM;
|
|
|
break;
|
|
|
}
|
|
|
if(stageError(current) < m_targetError) {
|
|
|
stopReason = LM_TARGET_ACHIEVED;
|
|
|
break;
|
|
|
}
|
|
|
if(m_totalEvaluations >= maximumEvaluations) {
|
|
|
stopReason = LM_MAX_ITERATIONS;
|
|
|
break;
|
|
|
}
|
|
|
// 重建已设置新基准；后续 LM 候选从该基准累计改善并确认停滞。
|
|
|
}
|
|
|
|
|
|
// 全程累计迭代；无可行方向、拒绝和求解失败也占用一次额度。
|
|
|
++phaseIterations;
|
|
|
++completedIterations;
|
|
|
emit progressUpdated(completedIterations, current.fitness);
|
|
|
const QVector<double> objectiveResidual = shapeStage
|
|
|
? current.breakdown.shapeResiduals : current.breakdown.residualVector;
|
|
|
const int rowOffset = shapeStage ? current.breakdown.residualVector.size() : 0;
|
|
|
const QVector<QVector<double> > objectiveJacobian = jacobian.mid(rowOffset, objectiveResidual.size());
|
|
|
const TrustRegionNormalEquations global = buildTrustRegionNormalEquations(
|
|
|
objectiveJacobian, objectiveResidual, jacobianColumnValid);
|
|
|
QVector<int> selectedColumns;
|
|
|
QVector<double> coordinateStep;
|
|
|
double predictedReduction = 0.0;
|
|
|
|
|
|
// 两段都联合调整全部有效自由参数，只切换残差及其对应的 Jacobian 行。
|
|
|
for(int column = 0; column < dimensions; ++column) {
|
|
|
if(jacobianColumnValid[column] && global.matrix[column][column] > 0.0)
|
|
|
selectedColumns.append(column);
|
|
|
}
|
|
|
if(selectedColumns.isEmpty() || !buildTrustRegionStep(
|
|
|
global, selectedColumns, current.coordinates, damping,
|
|
|
trustRadius, minimumCoordinateStep, &coordinateStep, &predictedReduction)) {
|
|
|
selectedColumns.clear();
|
|
|
}
|
|
|
|
|
|
// 没有可行下降步时收缩半径，按当前停滞状态重建灵敏度或结束。
|
|
|
if(selectedColumns.isEmpty()) {
|
|
|
emit logMessageGenerated(tr("Iteration %1 | No feasible descent step; keep the last accepted result.").arg(completedIterations));
|
|
|
if(trustRadius <= minimumTrustRadius * 1.01 &&
|
|
|
modelRebuiltAtMinimumRadius) {
|
|
|
stopReason = LM_LOCAL_OPTIMUM;
|
|
|
break;
|
|
|
}
|
|
|
trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
|
|
|
damping = qMin(1.0e8, damping * 4.0);
|
|
|
rebuildRequested = true;
|
|
|
rebuildReason = QT_TR_NOOP("no feasible descent step");
|
|
|
if(recordIneffectiveStep()) {
|
|
|
stopReason = LM_LOCAL_OPTIMUM;
|
|
|
break;
|
|
|
}
|
|
|
continue;
|
|
|
}
|
|
|
|
|
|
double stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
|
|
|
QVector<double> candidateCoordinates = current.coordinates;
|
|
|
for(int column = 0; column < dimensions; ++column) {
|
|
|
candidateCoordinates[column] += coordinateStep[column];
|
|
|
}
|
|
|
// 参数索引保留在跟踪文件中，便于检查每次联合调整的参数。
|
|
|
QStringList selectedParameterIndices;
|
|
|
for(int i = 0; i < selectedColumns.size(); ++i) {
|
|
|
selectedParameterIndices << QString::number(m_enabledParamIndices[selectedColumns[i]]);
|
|
|
}
|
|
|
QString selectionName = QString(shapeStage ? "shape_global_params_" : "total_global_params_") +
|
|
|
selectedParameterIndices.join("_");
|
|
|
|
|
|
TrustRegionEvaluation candidate;
|
|
|
candidate.coordinates = candidateCoordinates;
|
|
|
candidate.parameters = parametersFromCoordinates(candidate.coordinates);
|
|
|
candidate.valid = evaluateTrustRegionPoint(
|
|
|
candidate.parameters,
|
|
|
&candidate.fitness,
|
|
|
&candidate.breakdown,
|
|
|
&candidate.curve,
|
|
|
&candidate.elapsedMs);
|
|
|
|
|
|
if(m_shouldStop) {
|
|
|
restoreEvaluationState(current);
|
|
|
stopReason = LM_USER_STOPPED;
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
if(!candidate.valid) {
|
|
|
emitTrialLog(tr("Iteration %1").arg(completedIterations), current.breakdown, candidate.breakdown,
|
|
|
false, false, shapeStage);
|
|
|
// 求解失败的候选不能改变 current。先完整恢复上一个已接受参数和
|
|
|
// 对应误差快照，再缩小信赖域；连续失败达到上限才终止整个拟合。
|
|
|
++consecutiveSolverFailures;
|
|
|
++consecutiveRejectedSteps;
|
|
|
consecutivePoorPredictions = 0;
|
|
|
trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
|
|
|
damping = qMin(1.0e8, damping * 4.0);
|
|
|
writeTraceRow(m_currentIteration,
|
|
|
-1,
|
|
|
"trust_region_candidate",
|
|
|
candidate.parameters,
|
|
|
candidate.fitness,
|
|
|
false,
|
|
|
candidate.elapsedMs,
|
|
|
"solver_invalid_" + selectionName,
|
|
|
nullptr);
|
|
|
restoreEvaluationState(current);
|
|
|
if(consecutiveSolverFailures >= 2) {
|
|
|
rebuildRequested = true;
|
|
|
rebuildReason = QT_TR_NOOP("2 consecutive solver failures");
|
|
|
}
|
|
|
if(consecutiveSolverFailures >= m_maxConsecutiveFailures) {
|
|
|
stopReason = LM_CONSECUTIVE_FAILURES;
|
|
|
break;
|
|
|
}
|
|
|
if(recordIneffectiveStep()) {
|
|
|
stopReason = LM_LOCAL_OPTIMUM;
|
|
|
break;
|
|
|
}
|
|
|
continue;
|
|
|
}
|
|
|
|
|
|
if(m_shouldStop) {
|
|
|
restoreEvaluationState(current);
|
|
|
stopReason = LM_USER_STOPPED;
|
|
|
break;
|
|
|
}
|
|
|
consecutiveSolverFailures = 0;
|
|
|
// 有效候选即使最终被拒绝，也提供了一条真实割线，可用于修正下一轮
|
|
|
// 局部模型；是否成为新工作点由当前阶段的目标和约束共同决定。
|
|
|
const AutoFitObjectiveBreakdownLM oldBreakdown = current.breakdown;
|
|
|
updateTrustRegionJacobian(
|
|
|
&jacobian,
|
|
|
trustRegionFullResidual(oldBreakdown),
|
|
|
trustRegionFullResidual(candidate.breakdown),
|
|
|
coordinateStep);
|
|
|
// reductionRatio 衡量局部线性模型的可信度：接近 1 表示预测准确；
|
|
|
// 值较小表示虽然可能下降，但模型低估了非线性，需要收紧下一步。
|
|
|
const double objectiveEnergy = 0.5 * trustRegionSquaredNorm(objectiveResidual);
|
|
|
const QVector<double> candidateResidual = shapeStage
|
|
|
? candidate.breakdown.shapeResiduals : candidate.breakdown.residualVector;
|
|
|
const double candidateEnergy = 0.5 * trustRegionSquaredNorm(candidateResidual);
|
|
|
double actualReduction = objectiveEnergy - candidateEnergy;
|
|
|
double reductionRatio = predictedReduction > 1.0e-14 ? actualReduction / predictedReduction : 0.0;
|
|
|
// 使用固定采样和同一阶段目标比较预测与实际改善。接近收敛时的微小
|
|
|
// 预测量交给停滞逻辑处理，避免比例数值波动反复触发昂贵的全参数重建。
|
|
|
const double predictionFloor = qMax(1.0e-14, objectiveEnergy * 1.0e-8);
|
|
|
const bool poorPrediction = predictedReduction > predictionFloor &&
|
|
|
(!isFiniteNumber(reductionRatio) || reductionRatio < 0.25);
|
|
|
consecutivePoorPredictions = poorPrediction ? consecutivePoorPredictions + 1 : 0;
|
|
|
if(consecutivePoorPredictions >= 2) {
|
|
|
rebuildRequested = true;
|
|
|
rebuildReason = QT_TR_NOOP("2 consecutive steps with actual improvement below 25% of prediction");
|
|
|
}
|
|
|
bool accepted = acceptable(candidate, current);
|
|
|
|
|
|
if(accepted) {
|
|
|
// 当前阶段接受候选后同步发布参数、曲线和诊断。
|
|
|
// 模型预测可靠时减小阻尼并可扩大半径，预测较差时保守收缩。
|
|
|
acceptPoint(candidate);
|
|
|
++acceptedSinceRebuild;
|
|
|
consecutiveRejectedSteps = 0;
|
|
|
|
|
|
// 预测比只控制阻尼和信赖半径，接受候选仍以当前目标误差下降为准。
|
|
|
if(reductionRatio > 0.75) {
|
|
|
damping = qMax(1.0e-8, damping * 0.5);
|
|
|
if(stepNorm >= trustRadius * 0.8) {
|
|
|
trustRadius = qMin(
|
|
|
maximumTrustRadius, trustRadius * 1.6);
|
|
|
}
|
|
|
} else if(reductionRatio > 0.25) {
|
|
|
damping = qMax(1.0e-8, damping * 0.8);
|
|
|
} else {
|
|
|
damping = qMin(1.0e8, damping * 2.0);
|
|
|
trustRadius = qMax(
|
|
|
minimumTrustRadius, trustRadius * 0.75);
|
|
|
}
|
|
|
|
|
|
if(acceptedSinceRebuild >= 10 && !rebuildRequested) {
|
|
|
rebuildRequested = true;
|
|
|
rebuildReason = QT_TR_NOOP("10 accepted steps since last rebuild");
|
|
|
}
|
|
|
modelRebuiltAtMinimumRadius = false;
|
|
|
} else {
|
|
|
// 拒绝后恢复已接受参数、增大阻尼并缩小半径；按预测质量判断是否重建。
|
|
|
++consecutiveRejectedSteps;
|
|
|
damping = qMin(1.0e8, damping * 4.0);
|
|
|
trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
|
|
|
restoreEvaluationState(current);
|
|
|
}
|
|
|
|
|
|
writeTraceRow(m_currentIteration,
|
|
|
-1,
|
|
|
"trust_region_candidate",
|
|
|
candidate.parameters,
|
|
|
candidate.fitness,
|
|
|
true,
|
|
|
candidate.elapsedMs,
|
|
|
accepted
|
|
|
? "accepted_" + selectionName
|
|
|
: "rejected_" + selectionName,
|
|
|
&candidate.breakdown);
|
|
|
|
|
|
emitTrialLog(tr("Iteration %1").arg(completedIterations), oldBreakdown, candidate.breakdown,
|
|
|
true, accepted, shapeStage);
|
|
|
emit progressUpdated(completedIterations, m_globalBestFitness);
|
|
|
|
|
|
const bool effectiveImprovement = registerEffectiveImprovement(stageError(current));
|
|
|
if(!effectiveImprovement && recordIneffectiveStep()) stopReason = LM_LOCAL_OPTIMUM;
|
|
|
|
|
|
if(stopReason == LM_LOCAL_OPTIMUM) break;
|
|
|
|
|
|
if(stageError(current) < m_targetError) {
|
|
|
stopReason = LM_TARGET_ACHIEVED;
|
|
|
break;
|
|
|
}
|
|
|
if(trustRadius <= minimumTrustRadius * 1.01 &&
|
|
|
consecutiveRejectedSteps >= 2 && consecutivePoorPredictions >= 2) {
|
|
|
if(modelRebuiltAtMinimumRadius) {
|
|
|
stopReason = LM_LOCAL_OPTIMUM;
|
|
|
break;
|
|
|
}
|
|
|
rebuildRequested = true;
|
|
|
rebuildReason = QT_TR_NOOP("inaccurate model at minimum trust radius");
|
|
|
}
|
|
|
}
|
|
|
if(shapeStage) shapeIterations = phaseIterations;
|
|
|
else totalIterations = phaseIterations;
|
|
|
if(m_shouldStop || stopReason == LM_CONSECUTIVE_FAILURES || stopReason == LM_OPTIMIZATION_FAILED)
|
|
|
break;
|
|
|
// 形状部分结束后仅在总迭代额度有剩余时进入总误差部分。
|
|
|
if(shapeStage) {
|
|
|
emit logMessageGenerated(tr("Shape phase ended: %1").arg(getStopReasonDescription(stopReason)));
|
|
|
writeTraceRow(m_currentIteration, -1, "lm_phase_end", current.parameters, current.fitness,
|
|
|
true, 0, QString("shape_stop_reason_%1").arg(static_cast<int>(stopReason)), ¤t.breakdown);
|
|
|
} else {
|
|
|
emit logMessageGenerated(tr("Total-error phase ended: %1").arg(getStopReasonDescription(stopReason)));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// 灵敏度建模或初始误差已达标时可能没有候选迭代，不多计一次。
|
|
|
m_currentIteration = completedIterations - 1;
|
|
|
restoreEvaluationState(current);
|
|
|
DEBUG_OUT(tr("Adaptive fitting counts: %1 height iterations, %2 shape LM iterations, %3 total LM iterations, %4 total evaluations.")
|
|
|
.arg(heightTrials).arg(shapeIterations).arg(totalIterations).arg(m_totalEvaluations));
|
|
|
|
|
|
if(m_shouldStop) {
|
|
|
return LM_USER_STOPPED;
|
|
|
}
|
|
|
if(!shapeStage && current.fitness < m_targetError) {
|
|
|
return LM_TARGET_ACHIEVED;
|
|
|
}
|
|
|
if(completedIterations >= m_maxIterations && stopReason != LM_CONSECUTIVE_FAILURES && stopReason != LM_OPTIMIZATION_FAILED) {
|
|
|
return LM_MAX_ITERATIONS;
|
|
|
}
|
|
|
if(stopReason == LM_CONSECUTIVE_FAILURES ||
|
|
|
stopReason == LM_LOCAL_OPTIMUM ||
|
|
|
stopReason == LM_OPTIMIZATION_FAILED) {
|
|
|
return stopReason;
|
|
|
}
|
|
|
return LM_MAX_ITERATIONS;
|
|
|
}
|
|
|
|
|
|
double nmCalculationAutoFitLM::evaluateFitness(const QVector<double>& parameters, bool retrySolver)
|
|
|
{
|
|
|
// 校验并应用候选参数，调用求解器后比较目标井曲线；失败返回惩罚值 1e10。
|
|
|
const QString funcName = QString("evaluateError[%1]").arg(m_currentIteration);
|
|
|
static int callCount = 0;
|
|
|
callCount++;
|
|
|
m_lastEvaluatedLogLogData.clear();
|
|
|
m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
|
|
|
|
|
|
try {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Starting evaluation with %3 parameters")
|
|
|
.arg(funcName).arg(callCount).arg(parameters.size()));
|
|
|
|
|
|
QString paramStr = "Parameters: ";
|
|
|
|
|
|
for(int i = 0; i < parameters.size(); ++i) {
|
|
|
paramStr += QString("[%1]=%2 ").arg(i).arg(parameters[i], 0, 'f', 6);
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(QString("%1: %2").arg(funcName).arg(paramStr));
|
|
|
|
|
|
// 拦截维度不匹配、非有限数或超出用户范围的参数。
|
|
|
if(!validateParameters(parameters)) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - VALIDATION FAILED").arg(funcName).arg(callCount));
|
|
|
|
|
|
for(int i = 0; i < parameters.size(); ++i) {
|
|
|
if(!isFiniteNumber(parameters[i])) {
|
|
|
DEBUG_OUT(QString(" -> Param[%1] is NOT finite: %2").arg(i).arg(parameters[i]));
|
|
|
}
|
|
|
|
|
|
if(i < m_enabledParamIndices.size()) {
|
|
|
int paramIndex = m_enabledParamIndices[i];
|
|
|
|
|
|
if(paramIndex >= 0 && paramIndex < m_parameterLower.size()) {
|
|
|
double lower = m_parameterLower[paramIndex];
|
|
|
double upper = m_parameterUpper[paramIndex];
|
|
|
|
|
|
if(parameters[i] < lower) {
|
|
|
DEBUG_OUT(QString(" -> Param[%1]=%2 < lower bound %3")
|
|
|
.arg(i).arg(parameters[i]).arg(lower));
|
|
|
}
|
|
|
|
|
|
if(parameters[i] > upper) {
|
|
|
DEBUG_OUT(QString(" -> Param[%1]=%2 > upper bound %3")
|
|
|
.arg(i).arg(parameters[i]).arg(upper));
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Parameters validated OK").arg(funcName).arg(callCount));
|
|
|
|
|
|
// 参数写回和求解任务依赖当前数据管理器。
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
|
|
|
if(!dataManager) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - DataManager is NULL").arg(funcName).arg(callCount));
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - DataManager OK").arg(funcName).arg(callCount));
|
|
|
|
|
|
// 按勾选顺序将候选值写入储层和目标井。
|
|
|
try {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Applying parameters...").arg(funcName).arg(callCount));
|
|
|
applyParametersToDataManager(parameters);
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Parameters applied successfully").arg(funcName).arg(callCount));
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - FAILED to apply parameters: %3")
|
|
|
.arg(funcName).arg(callCount).arg(e.what()));
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
// Dfc 和裂缝半长属于网格输入。标记网格失效，使下一次任务基于当前参数快照重建。
|
|
|
const bool fractureGridParameterSelected =
|
|
|
(m_parameterSelected.size() > 4 && m_parameterSelected[4]) ||
|
|
|
(m_parameterSelected.size() > 5 && m_parameterSelected[5]);
|
|
|
if(fractureGridParameterSelected) {
|
|
|
dataManager->invalidatePebiGrid();
|
|
|
}
|
|
|
|
|
|
// 初始评价失败可重试；LM 候选失败则交给信赖域缩步，不重复相同参数。
|
|
|
QVector<QVector<double>> solverResult;
|
|
|
const int maxRetries = retrySolver ? 2 : 0;
|
|
|
bool solverSuccess = false;
|
|
|
|
|
|
for(int retry = 0; retry <= maxRetries; ++retry) {
|
|
|
if(m_shouldStop) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - User stop requested").arg(funcName).arg(callCount));
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
try {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Solver attempt %3/%4")
|
|
|
.arg(funcName).arg(callCount).arg(retry + 1).arg(maxRetries + 1));
|
|
|
|
|
|
if(retry > 0) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Retry delay...").arg(funcName).arg(callCount));
|
|
|
msleep(1000);
|
|
|
}
|
|
|
|
|
|
solverResult = runSolverDll();
|
|
|
|
|
|
if(solverResult.isEmpty()) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Solver returned EMPTY result").arg(funcName).arg(callCount));
|
|
|
} else {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Solver returned %3 arrays")
|
|
|
.arg(funcName).arg(callCount).arg(solverResult.size()));
|
|
|
|
|
|
for(int i = 0; i < solverResult.size(); ++i) {
|
|
|
DEBUG_OUT(QString(" -> Array[%1] size: %2").arg(i).arg(solverResult[i].size()));
|
|
|
}
|
|
|
|
|
|
if(validateSolverResult(solverResult)) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Solver result VALIDATED on attempt %3")
|
|
|
.arg(funcName).arg(callCount).arg(retry + 1));
|
|
|
solverSuccess = true;
|
|
|
break;
|
|
|
} else {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Solver result VALIDATION FAILED on attempt %3")
|
|
|
.arg(funcName).arg(callCount).arg(retry + 1));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Solver EXCEPTION on attempt %3: %4")
|
|
|
.arg(funcName).arg(callCount).arg(retry + 1).arg(e.what()));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(!solverSuccess) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - ALL SOLVER ATTEMPTS FAILED").arg(funcName).arg(callCount));
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
// 使用求解任务复制的目标井曲线，避免读取其他井或上一候选的共享结果。
|
|
|
QVector<QVector<double>> resultLogLogData = m_lastEvaluatedLogLogData;
|
|
|
|
|
|
try {
|
|
|
if(!validateLogLogData(resultLogLogData)) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - LogLog data VALIDATION FAILED")
|
|
|
.arg(funcName).arg(callCount));
|
|
|
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - LogLog data validated, size: %3")
|
|
|
.arg(funcName).arg(callCount).arg(resultLogLogData[0].size()));
|
|
|
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Error getting LogLog data: %3")
|
|
|
.arg(funcName).arg(callCount).arg(e.what()));
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
// 比较目标井的实测与模拟双对数曲线。
|
|
|
double error;
|
|
|
|
|
|
try {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Calculating error...")
|
|
|
.arg(funcName).arg(callCount));
|
|
|
|
|
|
error = calculateLogLogCurveError(m_targetLogLogData, resultLogLogData);
|
|
|
|
|
|
if(!isFiniteNumber(error) || error < 0) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - INVALID error value: %3")
|
|
|
.arg(funcName).arg(callCount).arg(error));
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - SUCCESS! Error = %3")
|
|
|
.arg(funcName).arg(callCount).arg(error, 0, 'e', 6));
|
|
|
// 保存最后一次有效曲线，供 LM 候选评价和精英保护复用。
|
|
|
m_lastEvaluatedLogLogData = resultLogLogData;
|
|
|
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - Error calculation FAILED: %3")
|
|
|
.arg(funcName).arg(callCount).arg(e.what()));
|
|
|
return 1e10;
|
|
|
}
|
|
|
|
|
|
return error;
|
|
|
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - TOP-LEVEL EXCEPTION: %3")
|
|
|
.arg(funcName).arg(callCount).arg(e.what()));
|
|
|
return 1e10;
|
|
|
} catch(...) {
|
|
|
DEBUG_OUT(QString("%1: Call #%2 - UNKNOWN TOP-LEVEL EXCEPTION")
|
|
|
.arg(funcName).arg(callCount));
|
|
|
return 1e10;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// ==================== 参数应用方法 ====================
|
|
|
|
|
|
void nmCalculationAutoFitLM::applyParametersToDataManager(const QVector<double>& parameters)
|
|
|
{
|
|
|
// 按 m_enabledParamIndices 对应的顺序，将勾选参数写入储层和目标井。
|
|
|
if(parameters.size() != getEnabledParameterCount()) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
updateReservoirParameters(parameters);
|
|
|
updateWellParameters(parameters);
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::updateReservoirParameters(const QVector<double>& parameters)
|
|
|
{
|
|
|
// 更新储层级参数。井级参数 skin/wellboreC 不在这里改，由 updateWellParameters() 负责。
|
|
|
// 这里先取 DataManager 中 reservoir 的副本，修改后再整体写回 DataManager。
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
|
|
|
|
|
|
// paramIndex 是勾选参数向量中的索引；i 是完整 6 个参数体系中的索引。
|
|
|
// 只有 m_parameterSelected[i] 为 true 时，才从 parameters 中消费一个值。
|
|
|
int paramIndex = 0;
|
|
|
|
|
|
for(int i = 0; i < m_parameterSelected.size(); ++i) {
|
|
|
if(m_parameterSelected[i] && paramIndex < parameters.size()) {
|
|
|
double value = parameters[paramIndex];
|
|
|
|
|
|
switch(i) {
|
|
|
case 0: // 渗透率
|
|
|
reservoirData.getPermeability().setValue(value);
|
|
|
break;
|
|
|
|
|
|
case 3: // 孔隙度
|
|
|
reservoirData.getPorosity().setValue(value);
|
|
|
break;
|
|
|
|
|
|
}
|
|
|
|
|
|
paramIndex++;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
dataManager->updateReservoirData(reservoirData);
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::updateWellParameters(const QVector<double>& parameters)
|
|
|
{
|
|
|
// 更新目标井上的拟合参数。目前井级可拟合参数主要是：
|
|
|
// - skin：写入第一个 perforation；
|
|
|
// - wellboreC：写入井筒储集系数；
|
|
|
// - Dfc/裂缝半长：只写入垂直压裂井或多段压裂水平井。
|
|
|
// 目标井不存在时返回；没有射孔数据时仅跳过表皮更新。
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
|
|
|
// 只更新当前目标井的井参数，避免多井项目中误改其他井。
|
|
|
nmDataWellBase* pWell = dataManager->findWellByName(m_targetWellName);
|
|
|
|
|
|
if(!pWell) return;
|
|
|
|
|
|
// 先在参数副本中组装本次候选值，全部解析完成后再写回现有井对象。
|
|
|
// 这样不会触发整井赋值，也不会删除并重建井内已有的射孔对象。
|
|
|
nmDataPerforation* pPerforation = pWell->getPerforation(0);
|
|
|
nmDataAttribute skinAttr;
|
|
|
bool updateSkin = false;
|
|
|
if(pPerforation) {
|
|
|
skinAttr = pPerforation->getSkin();
|
|
|
}
|
|
|
|
|
|
nmDataAttribute wellboreAttr = pWell->getWellboreStorage();
|
|
|
bool updateWellboreStorage = false;
|
|
|
|
|
|
nmDataVerticalFracturedWell* pVerticalFracturedWell =
|
|
|
dynamic_cast<nmDataVerticalFracturedWell*>(pWell);
|
|
|
nmDataHorizontalFracturedWell* pHorizontalFracturedWell =
|
|
|
dynamic_cast<nmDataHorizontalFracturedWell*>(pWell);
|
|
|
|
|
|
nmDataAttribute dfcAttr;
|
|
|
nmDataAttribute fractureHalfLengthAttr;
|
|
|
if(pVerticalFracturedWell) {
|
|
|
dfcAttr = pVerticalFracturedWell->getDfc();
|
|
|
fractureHalfLengthAttr = pVerticalFracturedWell->getFractureHalfLength();
|
|
|
} else if(pHorizontalFracturedWell) {
|
|
|
dfcAttr = pHorizontalFracturedWell->getDfc();
|
|
|
fractureHalfLengthAttr = pHorizontalFracturedWell->getFractureHalfLength();
|
|
|
}
|
|
|
bool updateDfc = false;
|
|
|
bool updateFractureHalfLength = false;
|
|
|
|
|
|
int paramIndex = 0;
|
|
|
|
|
|
for(int i = 0; i < m_parameterSelected.size(); ++i) {
|
|
|
if(m_parameterSelected[i] && paramIndex < parameters.size()) {
|
|
|
double value = parameters[paramIndex];
|
|
|
|
|
|
switch(i) {
|
|
|
case 1: { // 表皮系数
|
|
|
if(pPerforation) {
|
|
|
skinAttr.setValue(value);
|
|
|
updateSkin = true;
|
|
|
}
|
|
|
}
|
|
|
break;
|
|
|
|
|
|
case 2: { // 井筒储集系数
|
|
|
wellboreAttr.setValue(value);
|
|
|
updateWellboreStorage = true;
|
|
|
}
|
|
|
break;
|
|
|
|
|
|
case 4: { // 裂缝导流能力
|
|
|
if(pVerticalFracturedWell || pHorizontalFracturedWell) {
|
|
|
dfcAttr.setValue(value);
|
|
|
updateDfc = true;
|
|
|
}
|
|
|
}
|
|
|
break;
|
|
|
|
|
|
case 5: { // 裂缝半长
|
|
|
if(pVerticalFracturedWell || pHorizontalFracturedWell) {
|
|
|
fractureHalfLengthAttr.setValue(value);
|
|
|
updateFractureHalfLength = true;
|
|
|
}
|
|
|
}
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
paramIndex++;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(updateSkin) {
|
|
|
pPerforation->getSkin().setValue(skinAttr.getValue());
|
|
|
}
|
|
|
if(updateWellboreStorage) {
|
|
|
pWell->getWellboreStorage().setValue(wellboreAttr.getValue());
|
|
|
}
|
|
|
if(pVerticalFracturedWell) {
|
|
|
if(updateDfc) {
|
|
|
pVerticalFracturedWell->getDfc().setValue(dfcAttr.getValue());
|
|
|
}
|
|
|
if(updateFractureHalfLength) {
|
|
|
pVerticalFracturedWell->getFractureHalfLength().setValue(
|
|
|
fractureHalfLengthAttr.getValue());
|
|
|
}
|
|
|
} else if(pHorizontalFracturedWell) {
|
|
|
if(updateDfc) {
|
|
|
pHorizontalFracturedWell->getDfc().setValue(dfcAttr.getValue());
|
|
|
}
|
|
|
if(updateFractureHalfLength) {
|
|
|
pHorizontalFracturedWell->getFractureHalfLength().setValue(
|
|
|
fractureHalfLengthAttr.getValue());
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// ==================== 算法辅助方法 ====================
|
|
|
|
|
|
void nmCalculationAutoFitLM::validateAndProtectFinalResult()
|
|
|
{
|
|
|
// 最终精英保护只阻止无效结果或真正变差的结果。任何真实误差下降都应保留，
|
|
|
// 不能再用固定百分比门槛把已经找到的更优解恢复成初始值。
|
|
|
if(!m_hasValidUserSolution) {
|
|
|
DEBUG_OUT(tr("No initial solution for elite protection"));
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT(tr("=== Final Result Validation (Elite Protection) ==="));
|
|
|
|
|
|
double finalFitness = m_globalBestFitness;
|
|
|
double initialFitness = m_userInitialFitness;
|
|
|
|
|
|
DEBUG_OUT(tr("Comparing results: Initial=%1, Final=%2")
|
|
|
.arg(initialFitness, 0, 'e', 4).arg(finalFitness, 0, 'e', 4));
|
|
|
|
|
|
bool finalValid = isFiniteNumber(finalFitness) &&
|
|
|
finalFitness < 1.0e9 &&
|
|
|
m_globalBestPosition.size() ==
|
|
|
m_userInitialSolution.size() &&
|
|
|
!m_globalBestLogLogData.isEmpty() &&
|
|
|
m_globalBestObjectiveBreakdown.valid;
|
|
|
if(finalValid) {
|
|
|
double improvement = initialFitness - finalFitness;
|
|
|
double relativeImprovement =
|
|
|
improvement / qMax(1.0e-10, qAbs(initialFitness));
|
|
|
DEBUG_OUT(tr("Improvement: %1 (%2%)")
|
|
|
.arg(improvement, 0, 'e', 4)
|
|
|
.arg(relativeImprovement * 100, 0, 'f', 2));
|
|
|
}
|
|
|
|
|
|
if(!finalValid || finalFitness > initialFitness) {
|
|
|
DEBUG_OUT(
|
|
|
tr("Elite protection triggered: final result is invalid or worse than initial"));
|
|
|
emit logMessageGenerated(tr("Final total error did not improve; restoring the initial result."));
|
|
|
|
|
|
m_globalBestFitness = initialFitness;
|
|
|
m_globalBestPosition = m_userInitialSolution;
|
|
|
m_globalBestLogLogData = m_userInitialLogLogData;
|
|
|
m_globalBestObjectiveBreakdown = m_userInitialObjectiveBreakdown;
|
|
|
emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, m_globalBestFitness);
|
|
|
|
|
|
DEBUG_OUT(tr("Initial solution restored successfully"));
|
|
|
} else {
|
|
|
DEBUG_OUT(
|
|
|
tr("Final result validated - solution is not worse than initial"));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// ==================== 工具方法 ====================
|
|
|
|
|
|
int nmCalculationAutoFitLM::getEnabledParameterCount() const
|
|
|
{
|
|
|
// 返回用户勾选参与拟合的参数数量。
|
|
|
int count = 0;
|
|
|
|
|
|
for(int i = 0; i < m_parameterSelected.size(); ++i) {
|
|
|
if(m_parameterSelected[i]) count++;
|
|
|
}
|
|
|
|
|
|
return count;
|
|
|
}
|
|
|
|
|
|
// ==================== 验证和处理方法 ====================
|
|
|
|
|
|
bool nmCalculationAutoFitLM::validateParameters(const QVector<double>& parameters) const
|
|
|
{
|
|
|
// 参数物理范围已由拟合窗口统一校验；候选评价只检查维度、有限数和
|
|
|
// 用户设置的上下界，避免另一套硬编码阈值与实际搜索范围冲突。
|
|
|
if(parameters.size() != getEnabledParameterCount()) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
for(int i = 0; i < parameters.size(); ++i) {
|
|
|
if(!isFiniteNumber(parameters[i])) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(i < m_enabledParamIndices.size()) {
|
|
|
int paramIndex = m_enabledParamIndices[i];
|
|
|
|
|
|
if(paramIndex >= 0 && paramIndex < m_parameterLower.size()) {
|
|
|
if(parameters[i] < m_parameterLower[paramIndex] ||
|
|
|
parameters[i] > m_parameterUpper[paramIndex]) {
|
|
|
return false;
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
|
|
|
return true;
|
|
|
}
|
|
|
|
|
|
bool nmCalculationAutoFitLM::validateLogLogData(const QVector<QVector<double>>& logLogData) const
|
|
|
{
|
|
|
// 校验双对数曲线结构。约定：
|
|
|
// logLogData[0]=time，logLogData[1]=pressure，logLogData[2]=pressure derivative。
|
|
|
// 三列必须长度一致、至少包含五个点，且所有数值有限。
|
|
|
if(logLogData.size() < 3) {
|
|
|
DEBUG_OUT("LogLog data has less than 3 arrays");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
int size = logLogData[0].size();
|
|
|
|
|
|
if(size == 0) {
|
|
|
DEBUG_OUT("Empty LogLog data");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(logLogData[1].size() != size || logLogData[2].size() != size) {
|
|
|
DEBUG_OUT(QString("LogLog data size mismatch: X=%1, Y1=%2, Y2=%3")
|
|
|
.arg(logLogData[0].size())
|
|
|
.arg(logLogData[1].size())
|
|
|
.arg(logLogData[2].size()));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(size < 5) {
|
|
|
DEBUG_OUT(QString("Too few LogLog data points: %1").arg(size));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
for(int i = 0; i < size; ++i) {
|
|
|
if(!isFiniteNumber(logLogData[0][i]) ||
|
|
|
!isFiniteNumber(logLogData[1][i]) ||
|
|
|
!isFiniteNumber(logLogData[2][i])) {
|
|
|
DEBUG_OUT(QString("Invalid LogLog data at index %1").arg(i));
|
|
|
return false;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
return true;
|
|
|
}
|
|
|
|
|
|
bool nmCalculationAutoFitLM::validateInitialValues() const
|
|
|
{
|
|
|
// 检查当前模型读取出的初始参数是否和用户勾选维度一致，并且在上下界内。
|
|
|
// 如果初始值越界，算法仍可继续，但日志会提示，因为精英保护可能不可用或效果变差。
|
|
|
if(m_initialValues.size() != m_enabledParamIndices.size()) {
|
|
|
DEBUG_OUT("Initial values count mismatch with enabled parameters");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
bool allValid = true;
|
|
|
|
|
|
for(int i = 0; i < m_initialValues.size(); ++i) {
|
|
|
int paramIndex = m_enabledParamIndices[i];
|
|
|
double value = m_initialValues[i];
|
|
|
|
|
|
if(!isFiniteNumber(value)) {
|
|
|
DEBUG_OUT(QString("Initial value[%1] is not finite: %2").arg(i).arg(value));
|
|
|
allValid = false;
|
|
|
continue;
|
|
|
}
|
|
|
|
|
|
if(paramIndex < m_parameterLower.size() && paramIndex < m_parameterUpper.size()) {
|
|
|
double minVal = m_parameterLower[paramIndex];
|
|
|
double maxVal = m_parameterUpper[paramIndex];
|
|
|
|
|
|
if(value < minVal || value > maxVal) {
|
|
|
DEBUG_OUT(QString("Initial value[%1] = %2 is outside bounds [%3, %4]")
|
|
|
.arg(i).arg(value).arg(minVal).arg(maxVal));
|
|
|
allValid = false;
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
|
|
|
return allValid;
|
|
|
}
|
|
|
|
|
|
bool nmCalculationAutoFitLM::validateSolverResult(const QVector<QVector<double>>& result) const
|
|
|
{
|
|
|
// 校验求解器压力结果。这里检查的是 pressure result，至少需要 time 和 pressure 两列。
|
|
|
// result log-log 的结构会在 validateLogLogData() 中另行检查。
|
|
|
if(result.size() < 2) {
|
|
|
DEBUG_OUT("Solver result has less than 2 arrays");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(result[0].size() != result[1].size()) {
|
|
|
DEBUG_OUT(QString("Size mismatch: X=%1, Y=%2").arg(result[0].size()).arg(result[1].size()));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(result[0].size() == 0) {
|
|
|
DEBUG_OUT("Empty solver result");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(result[0].size() < 10) {
|
|
|
DEBUG_OUT(QString("Too few data points: %1").arg(result[0].size()));
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
for(int i = 0; i < result[0].size(); ++i) {
|
|
|
if(!isFiniteNumber(result[0][i]) || !isFiniteNumber(result[1][i])) {
|
|
|
DEBUG_OUT(QString("Invalid data at index %1: X=%2, Y=%3")
|
|
|
.arg(i).arg(result[0][i]).arg(result[1][i]));
|
|
|
return false;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// 单调性仅作诊断；误差计算会重新排序，不在此拒绝结果。
|
|
|
bool isMonotonic = true;
|
|
|
|
|
|
for(int i = 1; i < result[0].size(); ++i) {
|
|
|
if(result[0][i] <= result[0][i - 1]) {
|
|
|
isMonotonic = false;
|
|
|
break;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(!isMonotonic) {
|
|
|
DEBUG_OUT("X values are not monotonically increasing");
|
|
|
}
|
|
|
|
|
|
return true;
|
|
|
}
|
|
|
|
|
|
double nmCalculationAutoFitLM::calculateLogLogCurveError(
|
|
|
const QVector<QVector<double> >& target,
|
|
|
const QVector<QVector<double> >& result)
|
|
|
{
|
|
|
// 首次有效评价后固定公共时间范围，保证各候选在相同时间点比较。
|
|
|
// 总误差由完整压力和导数残差计算，形状阶段单独使用斜率残差。
|
|
|
// 整个计算过程均位于 log(time)-log(value) 坐标。
|
|
|
const double invalidLoss = 1.0e10;
|
|
|
const double valueFloor = 1.0e-12;
|
|
|
m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
|
|
|
|
|
|
if(!validateLogLogData(target) || !validateLogLogData(result)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
try {
|
|
|
// 非正导数无法进入双对数空间，跳过对应采样行，使用剩余有效点比较。
|
|
|
auto prepareCurve = [valueFloor](const QVector<QVector<double> >& data,
|
|
|
int firstIndex,
|
|
|
QVector<QPointF>* pressure,
|
|
|
QVector<QPointF>* derivative) -> bool {
|
|
|
if(!pressure || !derivative || data.size() < 3 ||
|
|
|
data[0].size() != data[1].size() ||
|
|
|
data[0].size() != data[2].size() ||
|
|
|
firstIndex < 0 || firstIndex >= data[0].size()) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
for(int i = firstIndex; i < data[0].size(); ++i) {
|
|
|
if(!isFiniteNumber(data[0][i]) ||
|
|
|
!isFiniteNumber(data[1][i]) ||
|
|
|
!isFiniteNumber(data[2][i]) ||
|
|
|
data[0][i] <= 0.0 ||
|
|
|
data[1][i] <= 0.0 ||
|
|
|
data[2][i] <= 0.0) {
|
|
|
continue;
|
|
|
}
|
|
|
|
|
|
pressure->append(QPointF(data[0][i], data[1][i]));
|
|
|
derivative->append(
|
|
|
QPointF(data[0][i], qMax(data[2][i], valueFloor)));
|
|
|
}
|
|
|
|
|
|
// 求解器输出可能不是严格升序，且同一时刻可能出现重复记录。
|
|
|
// 插值前统一排序并让后出现的记录覆盖同时间旧值，保证横坐标严格递增。
|
|
|
auto sortAndUnique = [](QVector<QPointF>* curve) {
|
|
|
std::stable_sort(
|
|
|
curve->begin(), curve->end(),
|
|
|
[](const QPointF& left, const QPointF& right) {
|
|
|
return left.x() < right.x();
|
|
|
});
|
|
|
|
|
|
QVector<QPointF> unique;
|
|
|
unique.reserve(curve->size());
|
|
|
for(int i = 0; i < curve->size(); ++i) {
|
|
|
if(unique.isEmpty() ||
|
|
|
curve->at(i).x() > unique.last().x()) {
|
|
|
unique.append(curve->at(i));
|
|
|
} else {
|
|
|
unique[unique.size() - 1] = curve->at(i);
|
|
|
}
|
|
|
}
|
|
|
*curve = unique;
|
|
|
};
|
|
|
|
|
|
sortAndUnique(pressure);
|
|
|
sortAndUnique(derivative);
|
|
|
return pressure->size() >= 3 && derivative->size() >= 3;
|
|
|
};
|
|
|
|
|
|
QVector<QPointF> targetPressure;
|
|
|
QVector<QPointF> targetDerivative;
|
|
|
QVector<QPointF> resultPressure;
|
|
|
QVector<QPointF> resultDerivative;
|
|
|
// 模拟结果已跳过 DLL 首点，因此误差计算同步忽略目标曲线首点。
|
|
|
if(!prepareCurve(target, 1, &targetPressure, &targetDerivative) ||
|
|
|
!prepareCurve(result, 0, &resultPressure, &resultDerivative)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// 在双对数坐标中插值。二分定位用于后面的多次水平配准试算。
|
|
|
auto interpolateLogValue = [valueFloor](
|
|
|
const QVector<QPointF>& curve,
|
|
|
double x,
|
|
|
double* value) -> bool {
|
|
|
if(!value || curve.size() < 2 || x <= 0.0 ||
|
|
|
x < curve.first().x() || x > curve.last().x()) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
int low = 0;
|
|
|
int high = curve.size() - 1;
|
|
|
while(low < high) {
|
|
|
int middle = low + (high - low) / 2;
|
|
|
if(curve[middle].x() < x) {
|
|
|
low = middle + 1;
|
|
|
} else {
|
|
|
high = middle;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
int right = qBound(1, low, curve.size() - 1);
|
|
|
int left = right - 1;
|
|
|
double leftLogX = qLn(curve[left].x());
|
|
|
double rightLogX = qLn(curve[right].x());
|
|
|
double denominator = rightLogX - leftLogX;
|
|
|
double leftLogY =
|
|
|
qLn(qMax(qAbs(curve[left].y()), valueFloor));
|
|
|
double rightLogY =
|
|
|
qLn(qMax(qAbs(curve[right].y()), valueFloor));
|
|
|
|
|
|
if(qAbs(denominator) <= 1.0e-12) {
|
|
|
*value = leftLogY;
|
|
|
} else {
|
|
|
double ratio = (qLn(x) - leftLogX) / denominator;
|
|
|
*value = leftLogY +
|
|
|
ratio * (rightLogY - leftLogY);
|
|
|
}
|
|
|
return isFiniteNumber(*value);
|
|
|
};
|
|
|
|
|
|
const double targetMinX = targetPressure.first().x();
|
|
|
const double targetMaxX = targetPressure.last().x();
|
|
|
const double resultMinX = resultPressure.first().x();
|
|
|
const double resultMaxX = resultPressure.last().x();
|
|
|
if(targetMinX <= 0.0 || targetMaxX <= targetMinX ||
|
|
|
resultMinX <= 0.0 || resultMaxX <= resultMinX) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// 仅首次有效评价使用交集建立基准；失败试算不能冻结区间，后续候选
|
|
|
// 必须覆盖完整基准，不允许靠丢失首尾点缩小误差。
|
|
|
const bool comparisonRangeFixed = m_comparisonTimeMin > 0.0;
|
|
|
double overlapMinX = comparisonRangeFixed
|
|
|
? m_comparisonTimeMin : qMax(targetMinX, resultMinX);
|
|
|
double overlapMaxX = comparisonRangeFixed
|
|
|
? m_comparisonTimeMax : qMin(targetMaxX, resultMaxX);
|
|
|
if(overlapMinX >= overlapMaxX) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
if(targetMinX > overlapMinX || targetMaxX < overlapMaxX ||
|
|
|
resultMinX > overlapMinX || resultMaxX < overlapMaxX) {
|
|
|
DEBUG_OUT("Candidate does not cover the fixed LM comparison time range");
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// 所有指标共用每个对数时间数量级 20 个间隔，向上取整并包含两端，至少三个点。
|
|
|
// 比较区间在本轮固定，因此点数也固定；微小容差避免整数量级因舍入多取一点。
|
|
|
const double timeDecades = (qLn(overlapMaxX) - qLn(overlapMinX)) / qLn(10.0);
|
|
|
const int numPoints = qMax(3, qCeil(timeDecades * kAutoFitIntervalsPerDecade - 1.0e-10) + 1);
|
|
|
QVector<double> commonX(numPoints);
|
|
|
QVector<double> commonLogX(numPoints);
|
|
|
QVector<double> targetLogPressure(numPoints);
|
|
|
QVector<double> targetLogDerivative(numPoints);
|
|
|
const double comparisonLogMinX = qLn(overlapMinX);
|
|
|
const double comparisonLogMaxX = qLn(overlapMaxX);
|
|
|
|
|
|
// 在公共时间范围内生成等距 log-time 网格，压力和导数各使用 numPoints 个残差。
|
|
|
for(int i = 0; i < numPoints; ++i) {
|
|
|
double logX = comparisonLogMinX +
|
|
|
static_cast<double>(i) *
|
|
|
(comparisonLogMaxX - comparisonLogMinX) /
|
|
|
(numPoints - 1);
|
|
|
commonLogX[i] = logX;
|
|
|
// 首尾直接使用原始端点，避免 exp(log(t)) 的舍入误差越过严格插值边界。
|
|
|
if(i == 0) {
|
|
|
commonX[i] = overlapMinX;
|
|
|
} else if(i == numPoints - 1) {
|
|
|
commonX[i] = overlapMaxX;
|
|
|
} else {
|
|
|
commonX[i] = qExp(logX);
|
|
|
}
|
|
|
|
|
|
if(!interpolateLogValue(
|
|
|
targetPressure, commonX[i],
|
|
|
&targetLogPressure[i]) ||
|
|
|
!interpolateLogValue(
|
|
|
targetDerivative, commonX[i],
|
|
|
&targetLogDerivative[i])) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
QVector<double> pressureResidual(
|
|
|
numPoints, std::numeric_limits<double>::quiet_NaN());
|
|
|
QVector<double> derivativeResidual(
|
|
|
numPoints, std::numeric_limits<double>::quiet_NaN());
|
|
|
|
|
|
// 残差定义为“模拟减目标”：正值表示模拟曲线偏高，负值表示偏低。
|
|
|
for(int i = 0; i < numPoints; ++i) {
|
|
|
double resultLogPressure = 0.0;
|
|
|
double resultLogDerivative = 0.0;
|
|
|
if(!interpolateLogValue(
|
|
|
resultPressure, commonX[i],
|
|
|
&resultLogPressure) ||
|
|
|
!interpolateLogValue(
|
|
|
resultDerivative, commonX[i],
|
|
|
&resultLogDerivative)) {
|
|
|
// 已位于结果时间范围内却无法插值说明数据存在内部断点，不能补线。
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
pressureResidual[i] =
|
|
|
resultLogPressure - targetLogPressure[i];
|
|
|
derivativeResidual[i] =
|
|
|
resultLogDerivative - targetLogDerivative[i];
|
|
|
}
|
|
|
|
|
|
auto populateStageMetrics = [&](AutoFitObjectiveBreakdownLM* objective) -> bool {
|
|
|
// 复用同一网格上已算出的原始残差，避免为高度和形状再次插值。
|
|
|
// 斜率跨度保持约为全对数时域的 10%，压力和导数仍等权。
|
|
|
const int count = numPoints;
|
|
|
const int lag = autoFitShapeLag(count);
|
|
|
const double span = qLn(overlapMaxX) - qLn(overlapMinX);
|
|
|
const QVector<double>* residuals[] = { &pressureResidual, &derivativeResidual };
|
|
|
double biases[2] = {0.0, 0.0};
|
|
|
for(int component = 0; component < 2; ++component) {
|
|
|
for(int i = 0; i < count; ++i) {
|
|
|
const double residual = (*residuals[component])[i];
|
|
|
// 使用对数时间梯形权重，避免密集段主导高度。
|
|
|
biases[component] += residual * ((i == 0 || i == count - 1) ? 0.5 : 1.0)
|
|
|
/ (count - 1);
|
|
|
}
|
|
|
}
|
|
|
objective->pressureVerticalBias = biases[0];
|
|
|
objective->derivativeVerticalBias = biases[1];
|
|
|
objective->verticalCommonBias = 0.5 * (biases[0] + biases[1]);
|
|
|
objective->verticalLoss = qAbs(objective->verticalCommonBias);
|
|
|
objective->verticalReliable = !(biases[0] * biases[1] < 0.0 &&
|
|
|
qMin(qAbs(biases[0]), qAbs(biases[1])) > 0.01);
|
|
|
objective->shapeResiduals.clear();
|
|
|
const double scale = qSqrt(0.5 / (count - lag));
|
|
|
for(int component = 0; component < 2; ++component) {
|
|
|
for(int i = 0; i < count - lag; ++i) {
|
|
|
objective->shapeResiduals.append(scale *
|
|
|
((*residuals[component])[i + lag] - (*residuals[component])[i]) /
|
|
|
(span * lag / (count - 1)));
|
|
|
}
|
|
|
}
|
|
|
objective->shapeLoss = qSqrt(trustRegionSquaredNorm(objective->shapeResiduals));
|
|
|
return isFiniteNumber(objective->shapeLoss);
|
|
|
};
|
|
|
|
|
|
AutoFitObjectiveBreakdownLM breakdown;
|
|
|
|
|
|
// 在指定中心附近计算普通均方根误差。
|
|
|
auto rmseAround = [](
|
|
|
const QVector<double>& values,
|
|
|
int begin,
|
|
|
int end,
|
|
|
double center) -> double {
|
|
|
double sum = 0.0;
|
|
|
int count = 0;
|
|
|
int validBegin = qMax(0, begin);
|
|
|
int validEnd =
|
|
|
qMin(end, static_cast<int>(values.size()));
|
|
|
|
|
|
for(int i = validBegin; i < validEnd; ++i) {
|
|
|
if(!isFiniteNumber(values[i])) {
|
|
|
continue;
|
|
|
}
|
|
|
|
|
|
double difference = values[i] - center;
|
|
|
sum += difference * difference;
|
|
|
++count;
|
|
|
}
|
|
|
|
|
|
return count > 0
|
|
|
? qSqrt(sum / count)
|
|
|
: std::numeric_limits<double>::quiet_NaN();
|
|
|
};
|
|
|
|
|
|
auto rmse = [&rmseAround](
|
|
|
const QVector<double>& values,
|
|
|
int begin,
|
|
|
int end) -> double {
|
|
|
return rmseAround(values, begin, end, 0.0);
|
|
|
};
|
|
|
|
|
|
// 普通算术平均中心保留上下偏差的符号。
|
|
|
auto meanCenterRange = [](
|
|
|
const QVector<double>& values,
|
|
|
int begin,
|
|
|
int end) -> double {
|
|
|
int validBegin = qMax(0, begin);
|
|
|
int validEnd =
|
|
|
qMin(end, static_cast<int>(values.size()));
|
|
|
double center = 0.0;
|
|
|
int count = 0;
|
|
|
|
|
|
for(int i = validBegin; i < validEnd; ++i) {
|
|
|
if(isFiniteNumber(values[i])) {
|
|
|
center += values[i];
|
|
|
++count;
|
|
|
}
|
|
|
}
|
|
|
if(count == 0) {
|
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
|
}
|
|
|
return center / count;
|
|
|
};
|
|
|
|
|
|
// 压力和导数合并后只求一个公共中心，表示两条曲线共同的上下位移。
|
|
|
// 分别去中心会把压力与导数之间真实的相对形状差异一并消除。
|
|
|
auto commonMeanCenterRange = [&meanCenterRange](
|
|
|
const QVector<double>& pressureValues,
|
|
|
const QVector<double>& derivativeValues,
|
|
|
int begin,
|
|
|
int end) -> double {
|
|
|
QVector<double> combined;
|
|
|
int validBegin = qMax(0, begin);
|
|
|
int validEnd = qMin(
|
|
|
end,
|
|
|
qMin(static_cast<int>(pressureValues.size()),
|
|
|
static_cast<int>(derivativeValues.size())));
|
|
|
combined.reserve(2 * qMax(0, validEnd - validBegin));
|
|
|
|
|
|
for(int i = validBegin; i < validEnd; ++i) {
|
|
|
if(isFiniteNumber(pressureValues[i])) {
|
|
|
combined.append(pressureValues[i]);
|
|
|
}
|
|
|
if(isFiniteNumber(derivativeValues[i])) {
|
|
|
combined.append(derivativeValues[i]);
|
|
|
}
|
|
|
}
|
|
|
return meanCenterRange(combined, 0, combined.size());
|
|
|
};
|
|
|
|
|
|
// 两个通道按能量等权合并，返回值与单通道 RMSE 保持同一量纲。
|
|
|
auto jointRmseAround = [&rmseAround](
|
|
|
const QVector<double>& pressureValues,
|
|
|
const QVector<double>& derivativeValues,
|
|
|
int begin,
|
|
|
int end,
|
|
|
double center) -> double {
|
|
|
double pressureLoss = rmseAround(
|
|
|
pressureValues, begin, end, center);
|
|
|
double derivativeLoss = rmseAround(
|
|
|
derivativeValues, begin, end, center);
|
|
|
if(!isFiniteNumber(pressureLoss) ||
|
|
|
!isFiniteNumber(derivativeLoss)) {
|
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
|
}
|
|
|
return qSqrt(0.5 *
|
|
|
(pressureLoss * pressureLoss +
|
|
|
derivativeLoss * derivativeLoss));
|
|
|
};
|
|
|
|
|
|
// 主目标始终使用未做上下或左右校正的完整曲线误差。
|
|
|
breakdown.pressureLoss =
|
|
|
rmse(pressureResidual, 0, numPoints);
|
|
|
breakdown.derivativeLoss =
|
|
|
rmse(derivativeResidual, 0, numPoints);
|
|
|
|
|
|
// 压力和导数各占一半权重。缩放后 residualVector 的二范数就是
|
|
|
// sqrt(0.5 * pressureLoss^2 + 0.5 * derivativeLoss^2)。
|
|
|
const double residualScale = qSqrt(0.5 / numPoints);
|
|
|
breakdown.residualVector.reserve(2 * numPoints);
|
|
|
for(int i = 0; i < numPoints; ++i) {
|
|
|
breakdown.residualVector.append(
|
|
|
residualScale *
|
|
|
pressureResidual[i]);
|
|
|
}
|
|
|
for(int i = 0; i < numPoints; ++i) {
|
|
|
breakdown.residualVector.append(
|
|
|
residualScale *
|
|
|
derivativeResidual[i]);
|
|
|
}
|
|
|
|
|
|
const double logGridStep =
|
|
|
(comparisonLogMaxX - comparisonLogMinX) /
|
|
|
(numPoints - 1);
|
|
|
const double resultLogMinX = qLn(resultMinX);
|
|
|
const double resultLogMaxX = qLn(resultMaxX);
|
|
|
// 左右配准只在所有候选位移都共同覆盖的固定区间比较，至少保留 80%
|
|
|
// 目标点；每个 log-time 网格间隔再细分为 8 份，提高位移分辨率。
|
|
|
const int minimumRegistrationPoints =
|
|
|
(numPoints * 4) / 5;
|
|
|
const int shiftSubdivisions = 8;
|
|
|
int maximumShiftIntervals = 4;
|
|
|
int registrationBegin = 0;
|
|
|
int registrationEnd = numPoints;
|
|
|
double maximumPhysicalShift =
|
|
|
maximumShiftIntervals * logGridStep;
|
|
|
|
|
|
// 所有位移候选使用同一组目标点。结果范围不足时逐步缩小最大位移，
|
|
|
// 但用于配准的固定公共区间不得少于目标网格的 80%。
|
|
|
auto updateRegistrationRange = [&](double maximumShift) {
|
|
|
registrationBegin = 0;
|
|
|
registrationEnd = numPoints;
|
|
|
while(registrationBegin < registrationEnd &&
|
|
|
commonLogX[registrationBegin] - maximumShift <
|
|
|
resultLogMinX - 1.0e-12) {
|
|
|
++registrationBegin;
|
|
|
}
|
|
|
while(registrationEnd > registrationBegin &&
|
|
|
commonLogX[registrationEnd - 1] + maximumShift >
|
|
|
resultLogMaxX + 1.0e-12) {
|
|
|
--registrationEnd;
|
|
|
}
|
|
|
};
|
|
|
|
|
|
updateRegistrationRange(maximumPhysicalShift);
|
|
|
while(maximumShiftIntervals > 0 &&
|
|
|
registrationEnd - registrationBegin <
|
|
|
minimumRegistrationPoints) {
|
|
|
--maximumShiftIntervals;
|
|
|
maximumPhysicalShift =
|
|
|
maximumShiftIntervals * logGridStep;
|
|
|
updateRegistrationRange(maximumPhysicalShift);
|
|
|
}
|
|
|
if(registrationEnd - registrationBegin <
|
|
|
minimumRegistrationPoints) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// physicalShift 为正表示模拟曲线偏右；对齐时在目标时刻右侧读取模拟值。
|
|
|
auto buildShiftResidual = [&](
|
|
|
double physicalShift,
|
|
|
int compareBegin,
|
|
|
int compareEnd,
|
|
|
QVector<double>* shiftedPressure,
|
|
|
QVector<double>* shiftedDerivative) -> bool {
|
|
|
if(!shiftedPressure || !shiftedDerivative) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
shiftedPressure->fill(
|
|
|
std::numeric_limits<double>::quiet_NaN(),
|
|
|
numPoints);
|
|
|
shiftedDerivative->fill(
|
|
|
std::numeric_limits<double>::quiet_NaN(),
|
|
|
numPoints);
|
|
|
|
|
|
int validBegin = qMax(0, compareBegin);
|
|
|
int validEnd = qMin(numPoints, compareEnd);
|
|
|
for(int i = validBegin; i < validEnd; ++i) {
|
|
|
double shiftedLogX =
|
|
|
commonLogX[i] + physicalShift;
|
|
|
if(shiftedLogX < resultLogMinX - 1.0e-12 ||
|
|
|
shiftedLogX >
|
|
|
resultLogMaxX + 1.0e-12) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
// 对数时间还原后再次限制到原始端点，避免 exp(log(t)) 的
|
|
|
// 舍入误差越过严格插值边界。
|
|
|
double shiftedX = qBound(
|
|
|
resultMinX,
|
|
|
qExp(qBound(resultLogMinX,
|
|
|
shiftedLogX,
|
|
|
resultLogMaxX)),
|
|
|
resultMaxX);
|
|
|
double resultLogPressure = 0.0;
|
|
|
double resultLogDerivative = 0.0;
|
|
|
if(!interpolateLogValue(
|
|
|
resultPressure, shiftedX,
|
|
|
&resultLogPressure) ||
|
|
|
!interpolateLogValue(
|
|
|
resultDerivative, shiftedX,
|
|
|
&resultLogDerivative)) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
(*shiftedPressure)[i] =
|
|
|
resultLogPressure - targetLogPressure[i];
|
|
|
(*shiftedDerivative)[i] =
|
|
|
resultLogDerivative - targetLogDerivative[i];
|
|
|
}
|
|
|
return true;
|
|
|
};
|
|
|
|
|
|
// 损失相同时优先选择绝对位移更小的候选，避免平坦 profile 在数值噪声
|
|
|
// 下无故偏向搜索边界。
|
|
|
auto isBetterProfileValue = [](
|
|
|
double loss,
|
|
|
double shift,
|
|
|
double bestLoss,
|
|
|
double bestShift) -> bool {
|
|
|
const double tolerance = 1.0e-12;
|
|
|
return loss < bestLoss - tolerance ||
|
|
|
(qAbs(loss - bestLoss) <= tolerance &&
|
|
|
qAbs(shift) < qAbs(bestShift));
|
|
|
};
|
|
|
|
|
|
const int halfShiftStepCount =
|
|
|
maximumShiftIntervals * shiftSubdivisions;
|
|
|
const double physicalShiftStep =
|
|
|
logGridStep / shiftSubdivisions;
|
|
|
double zeroShiftCenteredLoss =
|
|
|
std::numeric_limits<double>::quiet_NaN();
|
|
|
double bestCenteredLoss =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
double bestPhysicalShift = 0.0;
|
|
|
int bestShiftStep = 0;
|
|
|
double bestPressureLoss =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
double bestPressureShift = 0.0;
|
|
|
double bestDerivativeLoss =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
double bestDerivativeShift = 0.0;
|
|
|
QVector<double> profileLosses(
|
|
|
2 * halfShiftStepCount + 1,
|
|
|
std::numeric_limits<double>::quiet_NaN());
|
|
|
QVector<double> profileCommonBiases(
|
|
|
2 * halfShiftStepCount + 1,
|
|
|
std::numeric_limits<double>::quiet_NaN());
|
|
|
QVector<double> shiftedPressureResidual;
|
|
|
QVector<double> shiftedDerivativeResidual;
|
|
|
|
|
|
// 位移和公共上下偏移联合求解，避免“先扣上下还是先扣左右”的顺序依赖。
|
|
|
for(int shiftStep = -halfShiftStepCount;
|
|
|
shiftStep <= halfShiftStepCount;
|
|
|
++shiftStep) {
|
|
|
double physicalShift =
|
|
|
shiftStep * physicalShiftStep;
|
|
|
if(!buildShiftResidual(
|
|
|
physicalShift,
|
|
|
registrationBegin,
|
|
|
registrationEnd,
|
|
|
&shiftedPressureResidual,
|
|
|
&shiftedDerivativeResidual)) {
|
|
|
continue;
|
|
|
}
|
|
|
|
|
|
double commonBias = commonMeanCenterRange(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd);
|
|
|
double centeredLoss = jointRmseAround(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd,
|
|
|
commonBias);
|
|
|
double pressureBias = meanCenterRange(
|
|
|
shiftedPressureResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd);
|
|
|
double derivativeBias = meanCenterRange(
|
|
|
shiftedDerivativeResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd);
|
|
|
double pressureLoss = rmseAround(
|
|
|
shiftedPressureResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd,
|
|
|
pressureBias);
|
|
|
double derivativeLoss = rmseAround(
|
|
|
shiftedDerivativeResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd,
|
|
|
derivativeBias);
|
|
|
|
|
|
if(!isFiniteNumber(centeredLoss) ||
|
|
|
!isFiniteNumber(pressureLoss) ||
|
|
|
!isFiniteNumber(derivativeLoss)) {
|
|
|
continue;
|
|
|
}
|
|
|
int profileIndex = shiftStep + halfShiftStepCount;
|
|
|
profileLosses[profileIndex] = centeredLoss;
|
|
|
profileCommonBiases[profileIndex] = commonBias;
|
|
|
if(shiftStep == 0) {
|
|
|
zeroShiftCenteredLoss = centeredLoss;
|
|
|
}
|
|
|
if(isBetterProfileValue(
|
|
|
centeredLoss,
|
|
|
physicalShift,
|
|
|
bestCenteredLoss,
|
|
|
bestPhysicalShift)) {
|
|
|
bestCenteredLoss = centeredLoss;
|
|
|
bestPhysicalShift = physicalShift;
|
|
|
bestShiftStep = shiftStep;
|
|
|
}
|
|
|
if(isBetterProfileValue(
|
|
|
pressureLoss,
|
|
|
physicalShift,
|
|
|
bestPressureLoss,
|
|
|
bestPressureShift)) {
|
|
|
bestPressureLoss = pressureLoss;
|
|
|
bestPressureShift = physicalShift;
|
|
|
}
|
|
|
if(isBetterProfileValue(
|
|
|
derivativeLoss,
|
|
|
physicalShift,
|
|
|
bestDerivativeLoss,
|
|
|
bestDerivativeShift)) {
|
|
|
bestDerivativeLoss = derivativeLoss;
|
|
|
bestDerivativeShift = physicalShift;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(!isFiniteNumber(zeroShiftCenteredLoss) ||
|
|
|
!isFiniteNumber(bestCenteredLoss) ||
|
|
|
!isFiniteNumber(bestPressureLoss) ||
|
|
|
!isFiniteNumber(bestDerivativeLoss)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// horizontalGain 是允许水平位移后消除的均方误差贡献的平方根。
|
|
|
// 只有改善足够明显且最优点不是边界，才把位移解释为可靠左右偏差。
|
|
|
double horizontalGain = nestedRmsContribution(
|
|
|
zeroShiftCenteredLoss, bestCenteredLoss);
|
|
|
double horizontalSignalThreshold =
|
|
|
qMax(1.0e-5, zeroShiftCenteredLoss * 0.02);
|
|
|
int bestProfileIndex = bestShiftStep + halfShiftStepCount;
|
|
|
double nearbyProfileLoss =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
int leftProfileIndex =
|
|
|
bestProfileIndex - shiftSubdivisions;
|
|
|
int rightProfileIndex =
|
|
|
bestProfileIndex + shiftSubdivisions;
|
|
|
if(leftProfileIndex >= 0 &&
|
|
|
leftProfileIndex < profileLosses.size() &&
|
|
|
isFiniteNumber(profileLosses[leftProfileIndex])) {
|
|
|
nearbyProfileLoss = qMin(
|
|
|
nearbyProfileLoss,
|
|
|
profileLosses[leftProfileIndex]);
|
|
|
}
|
|
|
if(rightProfileIndex >= 0 &&
|
|
|
rightProfileIndex < profileLosses.size() &&
|
|
|
isFiniteNumber(profileLosses[rightProfileIndex])) {
|
|
|
nearbyProfileLoss = qMin(
|
|
|
nearbyProfileLoss,
|
|
|
profileLosses[rightProfileIndex]);
|
|
|
}
|
|
|
double profileContrast = isFiniteNumber(nearbyProfileLoss)
|
|
|
? nestedRmsContribution(
|
|
|
nearbyProfileLoss,
|
|
|
bestCenteredLoss)
|
|
|
: 0.0;
|
|
|
bool flatRegistrationProfile =
|
|
|
profileContrast <= horizontalSignalThreshold;
|
|
|
|
|
|
// 平坦的配准误差不一定表示歧义；只有近优位移显著改变公共偏差时，
|
|
|
// 才在旧配准诊断中标记高度和水平位移无法区分。
|
|
|
double minimumNearOptimalBias =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
double maximumNearOptimalBias =
|
|
|
-std::numeric_limits<double>::infinity();
|
|
|
for(int i = 0; i < profileLosses.size(); ++i) {
|
|
|
if(isFiniteNumber(profileLosses[i]) &&
|
|
|
isFiniteNumber(profileCommonBiases[i]) &&
|
|
|
profileLosses[i] <=
|
|
|
bestCenteredLoss + horizontalSignalThreshold) {
|
|
|
minimumNearOptimalBias = qMin(
|
|
|
minimumNearOptimalBias,
|
|
|
profileCommonBiases[i]);
|
|
|
maximumNearOptimalBias = qMax(
|
|
|
maximumNearOptimalBias,
|
|
|
profileCommonBiases[i]);
|
|
|
}
|
|
|
}
|
|
|
double nearOptimalBiasSpread =
|
|
|
isFiniteNumber(minimumNearOptimalBias) &&
|
|
|
isFiniteNumber(maximumNearOptimalBias)
|
|
|
? maximumNearOptimalBias - minimumNearOptimalBias
|
|
|
: std::numeric_limits<double>::infinity();
|
|
|
bool commonBiasStable = nearOptimalBiasSpread <= 1.0e-2;
|
|
|
bool horizontalAtBoundary =
|
|
|
halfShiftStepCount > 0 &&
|
|
|
qAbs(bestShiftStep) == halfShiftStepCount;
|
|
|
// 压力和导数通道分别求出的最佳位移若方向相反或相差过大，说明一个
|
|
|
// 单一水平平移无法解释两条曲线，此时将水平诊断标记为不可靠。
|
|
|
bool pressureShiftDetected =
|
|
|
qAbs(bestPressureShift) >=
|
|
|
0.5 * physicalShiftStep;
|
|
|
bool derivativeShiftDetected =
|
|
|
qAbs(bestDerivativeShift) >=
|
|
|
0.5 * physicalShiftStep;
|
|
|
bool channelShiftConflict =
|
|
|
pressureShiftDetected &&
|
|
|
derivativeShiftDetected &&
|
|
|
(bestPressureShift * bestDerivativeShift < 0.0 ||
|
|
|
qAbs(bestPressureShift - bestDerivativeShift) >
|
|
|
2.0 * logGridStep);
|
|
|
|
|
|
breakdown.horizontalLoss = horizontalGain;
|
|
|
breakdown.horizontalReliable =
|
|
|
maximumShiftIntervals > 0 &&
|
|
|
!horizontalAtBoundary &&
|
|
|
!channelShiftConflict &&
|
|
|
!flatRegistrationProfile &&
|
|
|
horizontalGain > horizontalSignalThreshold &&
|
|
|
qAbs(bestPhysicalShift) >=
|
|
|
0.5 * physicalShiftStep;
|
|
|
breakdown.registrationAmbiguous =
|
|
|
channelShiftConflict ||
|
|
|
(qAbs(bestPhysicalShift) >=
|
|
|
0.5 * physicalShiftStep &&
|
|
|
!breakdown.horizontalReliable) ||
|
|
|
(flatRegistrationProfile && !commonBiasStable);
|
|
|
breakdown.horizontalPhysicalShift =
|
|
|
breakdown.horizontalReliable
|
|
|
? bestPhysicalShift
|
|
|
: 0.0;
|
|
|
|
|
|
// 在可靠水平位移覆盖的区间计算旧配准诊断，保留其有效性检查。
|
|
|
int diagnosticBegin = 0;
|
|
|
int diagnosticEnd = numPoints;
|
|
|
while(diagnosticBegin < diagnosticEnd &&
|
|
|
commonLogX[diagnosticBegin] +
|
|
|
breakdown.horizontalPhysicalShift <
|
|
|
resultLogMinX - 1.0e-12) {
|
|
|
++diagnosticBegin;
|
|
|
}
|
|
|
while(diagnosticEnd > diagnosticBegin &&
|
|
|
commonLogX[diagnosticEnd - 1] +
|
|
|
breakdown.horizontalPhysicalShift >
|
|
|
resultLogMaxX + 1.0e-12) {
|
|
|
--diagnosticEnd;
|
|
|
}
|
|
|
if(diagnosticEnd - diagnosticBegin <
|
|
|
minimumRegistrationPoints ||
|
|
|
!buildShiftResidual(
|
|
|
breakdown.horizontalPhysicalShift,
|
|
|
diagnosticBegin,
|
|
|
diagnosticEnd,
|
|
|
&shiftedPressureResidual,
|
|
|
&shiftedDerivativeResidual)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
double commonBias = commonMeanCenterRange(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
diagnosticBegin,
|
|
|
diagnosticEnd);
|
|
|
double rawAlignedLoss = jointRmseAround(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
diagnosticBegin,
|
|
|
diagnosticEnd,
|
|
|
0.0);
|
|
|
double centeredAlignedLoss = jointRmseAround(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
diagnosticBegin,
|
|
|
diagnosticEnd,
|
|
|
commonBias);
|
|
|
if(!isFiniteNumber(commonBias) ||
|
|
|
!isFiniteNumber(rawAlignedLoss) ||
|
|
|
!isFiniteNumber(centeredAlignedLoss)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// 此处保留旧配准诊断的高度有效性检查。真正用于高度调整的指标，
|
|
|
// 在 populateStageMetrics() 中按未平移曲线的公共偏差重新计算。
|
|
|
breakdown.verticalCommonBias = commonBias;
|
|
|
breakdown.verticalLoss = nestedRmsContribution(
|
|
|
rawAlignedLoss, centeredAlignedLoss);
|
|
|
breakdown.verticalReliable =
|
|
|
!breakdown.registrationAmbiguous &&
|
|
|
breakdown.verticalLoss >
|
|
|
qMax(1.0e-5, rawAlignedLoss * 0.02);
|
|
|
|
|
|
QVector<double> shapePressure(
|
|
|
numPoints, std::numeric_limits<double>::quiet_NaN());
|
|
|
QVector<double> shapeDerivative(
|
|
|
numPoints, std::numeric_limits<double>::quiet_NaN());
|
|
|
for(int i = diagnosticBegin; i < diagnosticEnd; ++i) {
|
|
|
if(isFiniteNumber(shiftedPressureResidual[i])) {
|
|
|
shapePressure[i] =
|
|
|
shiftedPressureResidual[i] - commonBias;
|
|
|
}
|
|
|
if(isFiniteNumber(shiftedDerivativeResidual[i])) {
|
|
|
shapeDerivative[i] =
|
|
|
shiftedDerivativeResidual[i] - commonBias;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// 此处检查配准后的剩余误差；LM 使用的斜率形状误差在后面覆盖写入。
|
|
|
double shapePressureLoss = rmse(
|
|
|
shapePressure, diagnosticBegin, diagnosticEnd);
|
|
|
double shapeDerivativeLoss = rmse(
|
|
|
shapeDerivative, diagnosticBegin, diagnosticEnd);
|
|
|
if(!isFiniteNumber(shapePressureLoss) ||
|
|
|
!isFiniteNumber(shapeDerivativeLoss)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
breakdown.shapeLoss = qSqrt(
|
|
|
0.5 *
|
|
|
(shapePressureLoss * shapePressureLoss +
|
|
|
shapeDerivativeLoss * shapeDerivativeLoss));
|
|
|
|
|
|
// 现阶段不识别或单独调度晚期流动段；保留字段只为了维持现有 trace 列。
|
|
|
breakdown.lateDerivativeSlopeBias = 0.0;
|
|
|
breakdown.lateDerivativeTrendLoss = 0.0;
|
|
|
breakdown.lateDerivativeTrendReliable = false;
|
|
|
|
|
|
if(!isFiniteNumber(breakdown.pressureLoss) ||
|
|
|
!isFiniteNumber(breakdown.derivativeLoss) ||
|
|
|
!isFiniteNumber(breakdown.verticalCommonBias) ||
|
|
|
!isFiniteNumber(breakdown.verticalLoss) ||
|
|
|
!isFiniteNumber(breakdown.horizontalPhysicalShift) ||
|
|
|
!isFiniteNumber(breakdown.horizontalLoss) ||
|
|
|
!isFiniteNumber(breakdown.shapeLoss) ||
|
|
|
breakdown.residualVector.size() != 2 * numPoints) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// LM 总目标等于固定残差向量的二范数；压力和导数各占一半能量。
|
|
|
// 高度和斜率形状指标随后由同一固定网格的原始残差计算。
|
|
|
breakdown.total = qSqrt(
|
|
|
0.5 * breakdown.pressureLoss * breakdown.pressureLoss +
|
|
|
0.5 * breakdown.derivativeLoss * breakdown.derivativeLoss);
|
|
|
breakdown.valid =
|
|
|
isFiniteNumber(breakdown.total) &&
|
|
|
breakdown.total >= 0.0;
|
|
|
if(breakdown.valid && breakdown.total < 1.0e9) {
|
|
|
if(!comparisonRangeFixed) {
|
|
|
// 整个目标及诊断均有效后才提交基准，中点回退可重新建立区间。
|
|
|
m_comparisonTimeMin = overlapMinX;
|
|
|
m_comparisonTimeMax = overlapMaxX;
|
|
|
writeTraceMetaFile();
|
|
|
}
|
|
|
}
|
|
|
if(!populateStageMetrics(&breakdown)) return invalidLoss;
|
|
|
m_lastObjectiveBreakdown = breakdown;
|
|
|
|
|
|
DEBUG_OUT(
|
|
|
QString("LogLog objective: pressure=%1, derivative=%2, vertical=%3, horizontal=%4, shape=%5, ambiguous=%6, shift=%7, total=%8")
|
|
|
.arg(breakdown.pressureLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.derivativeLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.verticalLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.horizontalLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.shapeLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.registrationAmbiguous)
|
|
|
.arg(breakdown.horizontalPhysicalShift, 0, 'e', 4)
|
|
|
.arg(breakdown.total, 0, 'e', 4));
|
|
|
|
|
|
return breakdown.valid
|
|
|
? qMin(1.0e9, breakdown.total)
|
|
|
: invalidLoss;
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(
|
|
|
QString("Exception in LogLog error calculation: %1")
|
|
|
.arg(e.what()));
|
|
|
return invalidLoss;
|
|
|
} catch(...) {
|
|
|
DEBUG_OUT("Unknown exception in LogLog error calculation");
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
QVector<QVector<double>> nmCalculationAutoFitLM::runSolverDll()
|
|
|
{
|
|
|
// DLL 求解器路径。
|
|
|
// 这个函数负责把当前 DataManager 中的项目状态交给底层数值求解器，
|
|
|
// 数值求解仍包含全部计算井，以保留井间干扰；后处理只提取目标井曲线。
|
|
|
//
|
|
|
// 如果这里失败，通常需要优先检查：HX_NWTM.dll、license、网格/井数据是否完整、
|
|
|
// 目标井是否存在，以及 DataManager 中刚写入的参数是否导致求解器异常。
|
|
|
DEBUG_OUT("SOLVER DLL START");
|
|
|
|
|
|
// 创建任务时绑定当前分析，并在当前线程冻结本次求解所需的全部输入。
|
|
|
nmDataAnalyzeManager* dataManager =
|
|
|
nmDataAnalyzeManager::getCurrentInstance();
|
|
|
if(!dataManager || m_targetWellName.isEmpty()) {
|
|
|
DEBUG_OUT("Data manager or target well is unavailable");
|
|
|
return QVector<QVector<double> >();
|
|
|
}
|
|
|
|
|
|
if(m_evaluationInProgress > 0) {
|
|
|
DEBUG_OUT("DLL Solver already running, skipping");
|
|
|
return QVector<QVector<double>>();
|
|
|
}
|
|
|
|
|
|
++m_evaluationInProgress;
|
|
|
m_lastEvaluatedLogLogData.clear();
|
|
|
QVector<QVector<double>> result;
|
|
|
nmCalculationDllPebiSolverTask* dllTask = nullptr;
|
|
|
|
|
|
try {
|
|
|
DEBUG_OUT("Creating DLL solver task");
|
|
|
dllTask = new nmCalculationDllPebiSolverTask(
|
|
|
m_tempDirectory,
|
|
|
dataManager,
|
|
|
m_targetWellName);
|
|
|
|
|
|
if(m_shouldStop) {
|
|
|
DEBUG_OUT("Should stop - cleaning up and returning empty result");
|
|
|
delete dllTask;
|
|
|
--m_evaluationInProgress;
|
|
|
return result;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT("Starting DLL solver execution...");
|
|
|
|
|
|
// 异步执行 DLL 任务。循环等待期间持续 processEvents，保证界面不会完全卡死。
|
|
|
connect(dllTask, SIGNAL(solverInvoked()), this, SLOT(recordSolverCall()), Qt::DirectConnection);
|
|
|
dllTask->start();
|
|
|
|
|
|
// 等待完成。最大等待 1 小时；用户停止使用任务已有的协作取消接口。
|
|
|
const int maxWait = 3600000;
|
|
|
const int checkInterval = 50;
|
|
|
QTime waitTimer;
|
|
|
waitTimer.start();
|
|
|
|
|
|
while(waitTimer.elapsed() < maxWait) {
|
|
|
// 求解完成时立即返回；等待期间每 50 ms 处理停止请求和界面事件。
|
|
|
if(dllTask->wait(checkInterval)) {
|
|
|
DEBUG_OUT("DLL solver task completed");
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
// 等待期间也派发鼠标、键盘事件，让停止按钮能及时登记请求。
|
|
|
QApplication::processEvents(QEventLoop::AllEvents, checkInterval);
|
|
|
|
|
|
if(m_shouldStop) {
|
|
|
// DLL 没有中断接口，返回后丢弃结果，避免强制终止破坏其内部锁。
|
|
|
dllTask->requestCancel();
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(dllTask->isRunning()) {
|
|
|
DEBUG_OUT("DLL solver task timeout, terminating...");
|
|
|
dllTask->terminate();
|
|
|
dllTask->wait(2000);
|
|
|
|
|
|
delete dllTask;
|
|
|
dllTask = nullptr;
|
|
|
--m_evaluationInProgress;
|
|
|
return result;
|
|
|
}
|
|
|
|
|
|
// 线程结束后检查真实执行结果，防止失败时复用上一候选的旧曲线。
|
|
|
dllTask->wait();
|
|
|
if(m_shouldStop) {
|
|
|
DEBUG_OUT("DLL solver evaluation cancelled by user");
|
|
|
delete dllTask;
|
|
|
dllTask = nullptr;
|
|
|
--m_evaluationInProgress;
|
|
|
return result;
|
|
|
}
|
|
|
if(!dllTask->wasSuccessful()) {
|
|
|
DEBUG_OUT("DLL solver task reported failure");
|
|
|
delete dllTask;
|
|
|
dllTask = nullptr;
|
|
|
--m_evaluationInProgress;
|
|
|
return result;
|
|
|
}
|
|
|
|
|
|
// 任务结束后复制其局部结果，删除任务前不再持有任务内部引用。
|
|
|
QVector<QVector<double>> pressureResult = dllTask->getAutoFitResultPressure();
|
|
|
QVector<QVector<double>> logLogResult = dllTask->getAutoFitResultLogLog();
|
|
|
|
|
|
DEBUG_OUT(QString("DLL result verification - Pressure arrays: %1, LogLog arrays: %2")
|
|
|
.arg(pressureResult.size()).arg(logLogResult.size()));
|
|
|
|
|
|
if(pressureResult.size() >= 2) {
|
|
|
DEBUG_OUT(QString("Pressure result - Time points: %1, Pressure points: %2")
|
|
|
.arg(pressureResult[0].size()).arg(pressureResult[1].size()));
|
|
|
|
|
|
if(pressureResult[0].size() > 0) {
|
|
|
DEBUG_OUT(QString("Sample pressure data - Time[0]: %1, Time[last]: %2, P[0]: %3, P[last]: %4")
|
|
|
.arg(pressureResult[0][0])
|
|
|
.arg(pressureResult[0][pressureResult[0].size() - 1])
|
|
|
.arg(pressureResult[1][0])
|
|
|
.arg(pressureResult[1][pressureResult[1].size() - 1]));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(pressureResult.size() >= 2
|
|
|
&& pressureResult[0].size() > 0
|
|
|
&& pressureResult[1].size() > 0
|
|
|
&& validateLogLogData(logLogResult)) {
|
|
|
result = pressureResult;
|
|
|
m_lastEvaluatedLogLogData = logLogResult;
|
|
|
DEBUG_OUT(QString("Got DLL solver result: %1 points").arg(result[0].size()));
|
|
|
|
|
|
} else {
|
|
|
DEBUG_OUT("DLL solver result is empty or invalid");
|
|
|
DEBUG_OUT(QString("Pressure result size: %1, Array sizes: %2, %3")
|
|
|
.arg(pressureResult.size())
|
|
|
.arg(pressureResult.size() > 0 ? pressureResult[0].size() : 0)
|
|
|
.arg(pressureResult.size() > 1 ? pressureResult[1].size() : 0));
|
|
|
}
|
|
|
|
|
|
} catch(const std::bad_alloc& e) {
|
|
|
DEBUG_OUT(QString("Memory allocation failed in DLL solver: %1").arg(e.what()));
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(QString("Exception in DLL solver: %1").arg(e.what()));
|
|
|
} catch(...) {
|
|
|
DEBUG_OUT("Unknown exception in DLL solver");
|
|
|
}
|
|
|
|
|
|
// 释放求解任务，退出评价前处理待派发的界面事件。
|
|
|
if(dllTask) {
|
|
|
if(dllTask->isRunning()) {
|
|
|
dllTask->terminate();
|
|
|
dllTask->wait(3000);
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT("Cleaning up DLL solver task...");
|
|
|
delete dllTask;
|
|
|
dllTask = nullptr;
|
|
|
}
|
|
|
|
|
|
QApplication::processEvents(QEventLoop::AllEvents, 100);
|
|
|
--m_evaluationInProgress;
|
|
|
|
|
|
DEBUG_OUT(QString("SOLVER DLL END - ResultPoints: %1")
|
|
|
.arg(result.isEmpty() ? 0 : result[0].size()));
|
|
|
|
|
|
return result;
|
|
|
}
|
|
|
|
|
|
bool nmCalculationAutoFitLM::runFinalFullSolver()
|
|
|
{
|
|
|
// 不设置目标井名，任务按原完整模式保存全部井和网格结果。
|
|
|
nmDataAnalyzeManager* dataManager =
|
|
|
nmDataAnalyzeManager::getCurrentInstance();
|
|
|
if(!dataManager) {
|
|
|
DEBUG_OUT("Cannot start final full-field solver without a data manager");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(m_evaluationInProgress > 0) {
|
|
|
DEBUG_OUT("Cannot start final full-field solver while another evaluation is running");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
++m_evaluationInProgress;
|
|
|
nmCalculationDllPebiSolverTask dllTask(m_tempDirectory, dataManager);
|
|
|
connect(&dllTask, SIGNAL(solverInvoked()), this, SLOT(recordSolverCall()), Qt::DirectConnection);
|
|
|
dllTask.start();
|
|
|
|
|
|
const int maxWait = 3600000;
|
|
|
const int checkInterval = 50;
|
|
|
QTime waitTimer;
|
|
|
waitTimer.start();
|
|
|
|
|
|
while(waitTimer.elapsed() < maxWait) {
|
|
|
if(dllTask.wait(checkInterval)) {
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
// 最终完整计算期间处理界面事件，保持窗口响应。
|
|
|
QApplication::processEvents(QEventLoop::AllEvents, checkInterval);
|
|
|
|
|
|
if(!dllTask.isRunning()) {
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
}
|
|
|
|
|
|
if(dllTask.isRunning()) {
|
|
|
DEBUG_OUT("Final full-field solver timeout, terminating task");
|
|
|
dllTask.terminate();
|
|
|
dllTask.wait(2000);
|
|
|
--m_evaluationInProgress;
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
dllTask.wait();
|
|
|
// 后台只生成局部结果快照；确认求解和输入版本均有效后再一次性写回当前分析。
|
|
|
const bool succeeded =
|
|
|
dllTask.wasSuccessful() && dllTask.commitResult(dataManager);
|
|
|
--m_evaluationInProgress;
|
|
|
return succeeded;
|
|
|
}
|
|
|
|
|
|
QString nmCalculationAutoFitLM::getStopReasonDescription(StopReasonLM reason) const
|
|
|
{
|
|
|
// 将停止枚举转成人类可读文本，用于日志和运行摘要。
|
|
|
switch(reason) {
|
|
|
case LM_TARGET_ACHIEVED:
|
|
|
return tr("Target error achieved");
|
|
|
|
|
|
case LM_LOCAL_OPTIMUM:
|
|
|
return tr("Local optimum detected");
|
|
|
|
|
|
case LM_MAX_ITERATIONS:
|
|
|
return tr("Iteration or evaluation budget reached");
|
|
|
|
|
|
case LM_USER_STOPPED:
|
|
|
return tr("Stopped by user request");
|
|
|
|
|
|
case LM_CONSECUTIVE_FAILURES:
|
|
|
return tr("Too many consecutive failures");
|
|
|
|
|
|
case LM_OPTIMIZATION_FAILED:
|
|
|
return tr("Optimization failed");
|
|
|
|
|
|
default:
|
|
|
return tr("Unknown reason");
|
|
|
}
|
|
|
}
|