You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
nmWTAI-Platform/Src/nmNum/nmData/nmDataSensitive.cpp

927 lines
37 KiB
C++

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

#include "nmDataSensitive.h"
#include <cmath>
#include <limits>
#include <random>
#include <QDateTime>
#include <QtCore/qnumeric.h>
bool nmDataSensitive::BatchResult::statisticValue(const RunResult& run, int variable, double& value) const
{
// 表格、直方图和散点图共用有效性判断,拟合参数不依赖压力误差是否可用。
if(!run.m_bSucceeded || run.m_bCancelled) return false;
if(variable == -1) {
value = run.m_dPressureRmse;
return qIsFinite(value) && value >= 0.0;
}
if(variable < 0 || variable >= m_fitParameters.size() || variable >= run.m_vecFittedParameters.size()) return false;
value = run.m_vecFittedParameters[variable];
return qIsFinite(value);
}
bool nmDataSensitive::BatchResult::filteredParameterBounds(int index, double& minimum, double& maximum, bool fitted) const
{
// 采样和拟合参数共用误差严格小于阈值的运行,不受曲线可见性影响。
bool found = false;
for(size_t i = 0; i < m_runs.size(); ++i) {
const RunResult& run = m_runs[i];
const QVector<double>& parameters = fitted ? run.m_vecFittedParameters : run.m_vecParameters;
if(!run.m_bSucceeded || run.m_bCancelled || !qIsFinite(run.m_dPressureRmse) ||
run.m_dPressureRmse < 0.0 || !(run.m_dPressureRmse < m_dFilterThreshold) ||
index < 0 || index >= parameters.size()) continue;
const double value = parameters[index];
if(!qIsFinite(value)) continue;
minimum = found ? qMin(minimum, value) : value;
maximum = found ? qMax(maximum, value) : value;
found = true;
}
return found;
}
namespace
{
double normalCdf(double z)
{
// 标准正态累积分布近似,最大绝对误差约 7.5e-8,兼容 VS2010。
const double x = std::fabs(z), t = 1.0 / (1.0 + 0.2316419 * x);
const double tail = 0.3989422804014327 * std::exp(-0.5 * x * x) * t *
(0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429))));
return z < 0.0 ? tail : 1.0 - tail;
}
double normalQuantile(double probability, double low = -10.0, double high = 10.0)
{
// 在累积分布上二分求逆,避免引入新的数学库。
for(int i = 0; i < 60; ++i) {
const double mid = (low + high) * 0.5;
if(normalCdf(mid) < probability) low = mid;
else high = mid;
}
return (low + high) * 0.5;
}
rapidjson::Value sensitiveNumbersToJson(const QVector<double>& values,
rapidjson::Document::AllocatorType& allocator)
{
// 使用数值数组保存曲线精度,不经过界面显示文本转换。
rapidjson::Value array(rapidjson::kArrayType);
for(int i = 0; i < values.size(); ++i) {
// 求解曲线中的无效点保存为空值,避免 NaN 或无穷值使整个项目 JSON 无法写入。
if(qIsFinite(values[i])) array.PushBack(values[i], allocator);
else array.PushBack(rapidjson::Value().Move(), allocator);
}
return array;
}
rapidjson::Value sensitiveCurvesToJson(const QVector<QVector<double> >& curves,
rapidjson::Document::AllocatorType& allocator)
{
// 曲线列顺序与求解器返回值一致。
rapidjson::Value array(rapidjson::kArrayType);
for(int i = 0; i < curves.size(); ++i)
array.PushBack(sensitiveNumbersToJson(curves[i], allocator), allocator);
return array;
}
bool sensitiveNumbersFromJson(const rapidjson::Value& array, QVector<double>& values)
{
// 损坏的数值列整体丢弃,避免跳过某个元素后时间与压力错位。
values.clear();
if(!array.IsArray()) return false;
for(rapidjson::SizeType i = 0; i < array.Size(); ++i) {
// 空值保留原位置,绘图时过滤无效点,不改变其他时间和值的对应关系。
if(array[i].IsNull()) {
values.append(std::numeric_limits<double>::quiet_NaN());
continue;
}
if(!array[i].IsNumber() || !qIsFinite(array[i].GetDouble())) {
values.clear();
return false;
}
values.append(array[i].GetDouble());
}
return true;
}
QVector<QVector<double> > sensitiveCurvesFromJson(const rapidjson::Value& value, const char* key)
{
// 兼容旧成果中不存在批量曲线的情况。
QVector<QVector<double> > curves;
if(!value.HasMember(key) || !value[key].IsArray()) return curves;
const rapidjson::Value& array = value[key];
for(rapidjson::SizeType i = 0; i < array.Size(); ++i) {
QVector<double> column;
if(!sensitiveNumbersFromJson(array[i], column)) return QVector<QVector<double> >();
curves.append(column);
}
return curves;
}
QString sensitiveTextFromJson(const rapidjson::Value& value, const char* key)
{
// 文本统一使用 UTF-8,参数名称和失败原因可包含中文。
return value.HasMember(key) && value[key].IsString()
? QString::fromUtf8(value[key].GetString()) : QString();
}
double sensitiveNumberFromJson(const rapidjson::Value& value, const char* key, double fallback)
{
// 新增数值字段缺失或无效时沿用默认值,兼容之前保存的成果。
return value.HasMember(key) && value[key].IsNumber() && qIsFinite(value[key].GetDouble())
? value[key].GetDouble() : fallback;
}
}
rapidjson::Value nmDataSensitive::BatchResult::ToJsonValue(
rapidjson::Document::AllocatorType& allocator) const
{
// 批次保存自己的配置及目标井,不引用之后可能被修改的当前模型。
rapidjson::Value value(rapidjson::kObjectType);
value.AddMember("CalculationType", int(m_calcType), allocator);
value.AddMember("WellCode", rapidjson::Value(m_sWellCode.toUtf8().constData(), allocator).Move(), allocator);
value.AddMember("WellName", rapidjson::Value(m_sWellName.toUtf8().constData(), allocator).Move(), allocator);
value.AddMember("AxisXName", rapidjson::Value(m_sAxisXName.toUtf8().constData(), allocator).Move(), allocator);
value.AddMember("AxisXUnit", rapidjson::Value(m_sAxisXUnit.toUtf8().constData(), allocator).Move(), allocator);
value.AddMember("AxisYName", rapidjson::Value(m_sAxisYName.toUtf8().constData(), allocator).Move(), allocator);
value.AddMember("AxisYUnit", rapidjson::Value(m_sAxisYUnit.toUtf8().constData(), allocator).Move(), allocator);
value.AddMember("Completed", m_bCompleted, allocator);
value.AddMember("Stopped", m_bStopped, allocator);
value.AddMember("FailureReason", rapidjson::Value(m_sFailureReason.toUtf8().constData(), allocator).Move(), allocator);
value.AddMember("HistoryLogLog", sensitiveCurvesToJson(m_vecHistoryLogLog, allocator), allocator);
value.AddMember("BaseLogLog", sensitiveCurvesToJson(m_vecBaseLogLog, allocator), allocator);
// 统计图和结果参数页直接使用已保存的误差、拟合值及筛选设置。
value.AddMember("FilterThreshold", qIsFinite(m_dFilterThreshold) ? m_dFilterThreshold : -1.0, allocator);
value.AddMember("HistogramVariable", m_nHistogramVariable, allocator);
value.AddMember("HistogramBins", m_nHistogramBins, allocator);
rapidjson::Value fitting(rapidjson::kArrayType);
for(int i = 0; i < m_fitParameters.size(); ++i) {
rapidjson::Value parameter(rapidjson::kObjectType);
parameter.AddMember("Attribute", m_fitParameters[i].ToJsonValue(allocator), allocator);
parameter.AddMember("Code", rapidjson::Value(m_fitParameterCodes.value(i).toUtf8().constData(), allocator).Move(), allocator);
fitting.PushBack(parameter, allocator);
}
value.AddMember("FitParameters", fitting, allocator);
rapidjson::Value variables(rapidjson::kArrayType);
for(int i = 0; i < m_variables.size(); ++i)
variables.PushBack(m_variables[i].ToJsonValue(allocator), allocator);
value.AddMember("Variables", variables, allocator);
rapidjson::Value runs(rapidjson::kArrayType);
for(size_t i = 0; i < m_runs.size(); ++i) {
const RunResult& run = m_runs[i];
rapidjson::Value item(rapidjson::kObjectType);
item.AddMember("RunIndex", run.m_nRunIndex, allocator);
item.AddMember("Parameters", sensitiveNumbersToJson(run.m_vecParameters, allocator), allocator);
item.AddMember("FittedParameters", sensitiveNumbersToJson(run.m_vecFittedParameters, allocator), allocator);
item.AddMember("PressureRmse", qIsFinite(run.m_dPressureRmse) ? run.m_dPressureRmse : -1.0, allocator);
item.AddMember("Pressure", sensitiveCurvesToJson(run.m_vecPressure, allocator), allocator);
item.AddMember("LogLog", sensitiveCurvesToJson(run.m_vecLogLog, allocator), allocator);
item.AddMember("SemiLog", sensitiveCurvesToJson(run.m_vecSemiLog, allocator), allocator);
item.AddMember("Succeeded", run.m_bSucceeded, allocator);
item.AddMember("Cancelled", run.m_bCancelled, allocator);
item.AddMember("Visible", run.m_bVisible, allocator);
item.AddMember("FailureReason", rapidjson::Value(run.m_sFailureReason.toUtf8().constData(), allocator).Move(), allocator);
runs.PushBack(item, allocator);
}
value.AddMember("Runs", runs, allocator);
return value;
}
void nmDataSensitive::BatchResult::FromJsonValue(const rapidjson::Value& value)
{
// 旧成果无需结果字段;新成果恢复原运行编号、状态和显示选择。
*this = BatchResult();
if(!value.IsObject()) return;
if(value.HasMember("CalculationType") && value["CalculationType"].IsInt() &&
value["CalculationType"].GetInt() >= CALC_DETERMINISTIC && value["CalculationType"].GetInt() <= CALC_MONTE_CARLO)
m_calcType = static_cast<CalculationType>(value["CalculationType"].GetInt());
m_sWellCode = sensitiveTextFromJson(value, "WellCode");
m_sWellName = sensitiveTextFromJson(value, "WellName");
m_sAxisXName = sensitiveTextFromJson(value, "AxisXName");
m_sAxisXUnit = sensitiveTextFromJson(value, "AxisXUnit");
m_sAxisYName = sensitiveTextFromJson(value, "AxisYName");
m_sAxisYUnit = sensitiveTextFromJson(value, "AxisYUnit");
m_sFailureReason = sensitiveTextFromJson(value, "FailureReason");
m_bCompleted = value.HasMember("Completed") && value["Completed"].IsBool() && value["Completed"].GetBool();
m_bStopped = value.HasMember("Stopped") && value["Stopped"].IsBool() && value["Stopped"].GetBool();
m_vecHistoryLogLog = sensitiveCurvesFromJson(value, "HistoryLogLog");
m_vecBaseLogLog = sensitiveCurvesFromJson(value, "BaseLogLog");
m_dFilterThreshold = sensitiveNumberFromJson(value, "FilterThreshold", -1.0);
if(value.HasMember("FitParameters") && value["FitParameters"].IsArray()) {
const rapidjson::Value& fitting = value["FitParameters"];
for(rapidjson::SizeType i = 0; i < fitting.Size(); ++i) {
// 保留参数位置,确保每组拟合值与名称、单位和 XML 编码一一对应。
nmDataAttribute parameter;
QString code;
if(fitting[i].IsObject()) {
if(fitting[i].HasMember("Attribute") && fitting[i]["Attribute"].IsObject())
parameter.FromJsonValue(fitting[i]["Attribute"]);
code = sensitiveTextFromJson(fitting[i], "Code");
}
m_fitParameters.append(parameter);
m_fitParameterCodes.append(code);
}
}
if(value.HasMember("HistogramVariable") && value["HistogramVariable"].IsInt()) {
const int index = value["HistogramVariable"].GetInt();
if(index >= -1 && index < m_fitParameters.size()) m_nHistogramVariable = index;
}
if(value.HasMember("HistogramBins") && value["HistogramBins"].IsInt())
m_nHistogramBins = qBound(1, value["HistogramBins"].GetInt(), 1000);
if(value.HasMember("Variables") && value["Variables"].IsArray()) {
const rapidjson::Value& variables = value["Variables"];
for(rapidjson::SizeType i = 0; i < variables.Size(); ++i) {
if(!variables[i].IsObject()) continue;
VariableSampling variable;
variable.FromJsonValue(variables[i]);
m_variables.append(variable);
}
}
if(!value.HasMember("Runs") || !value["Runs"].IsArray()) return;
const rapidjson::Value& runs = value["Runs"];
for(rapidjson::SizeType i = 0; i < runs.Size(); ++i) {
const rapidjson::Value& item = runs[i];
if(!item.IsObject() || !item.HasMember("RunIndex") || !item["RunIndex"].IsUint64()) continue;
RunResult run;
run.m_nRunIndex = item["RunIndex"].GetUint64();
if(item.HasMember("Parameters")) sensitiveNumbersFromJson(item["Parameters"], run.m_vecParameters);
if(item.HasMember("FittedParameters")) sensitiveNumbersFromJson(item["FittedParameters"], run.m_vecFittedParameters);
run.m_dPressureRmse = sensitiveNumberFromJson(item, "PressureRmse", -1.0);
run.m_vecPressure = sensitiveCurvesFromJson(item, "Pressure");
run.m_vecLogLog = sensitiveCurvesFromJson(item, "LogLog");
run.m_vecSemiLog = sensitiveCurvesFromJson(item, "SemiLog");
run.m_bSucceeded = item.HasMember("Succeeded") && item["Succeeded"].IsBool() && item["Succeeded"].GetBool();
run.m_bCancelled = item.HasMember("Cancelled") && item["Cancelled"].IsBool() && item["Cancelled"].GetBool();
if(item.HasMember("Visible") && item["Visible"].IsBool()) run.m_bVisible = item["Visible"].GetBool();
run.m_sFailureReason = sensitiveTextFromJson(item, "FailureReason");
m_runs.push_back(run);
}
}
nmDataSensitive::VariableSampling::VariableSampling()
: m_varGroup("")
, m_varName("")
, m_enabled(false)
, m_mode(MODE_AUTOMATIC)
, m_log(false)
, m_number(0)
, m_distribution(NORMAL), m_distributionInitialized(false), m_mean(0.0), m_stdDev(1.0), m_peak(0.0)
, m_distributionMin(0.0), m_distributionMax(1.0)
{
m_modelValue = nmDataAttribute("Model value", "", "");
m_minValue = nmDataAttribute("Minimum", "", "");
m_maxValue = nmDataAttribute("Maximum", "", "");
m_values.clear();
}
// 序列化 VariableSampling 为 RapidJSON Value
rapidjson::Value nmDataSensitive::VariableSampling::ToJsonValue(
rapidjson::Document::AllocatorType& allocator) const
{
rapidjson::Value obj(rapidjson::kObjectType);
// 字符串字段
obj.AddMember("VarGroup",
rapidjson::Value(m_varGroup.toUtf8().constData(), allocator).Move(),
allocator);
obj.AddMember("VarName",
rapidjson::Value(m_varName.toUtf8().constData(), allocator).Move(),
allocator);
// 简单类型
obj.AddMember("Enabled", m_enabled, allocator);
obj.AddMember("Mode", static_cast<int>(m_mode), allocator);
obj.AddMember("Log", m_log, allocator);
obj.AddMember("Number", m_number, allocator);
// 六种分布使用同一组字段,恢复时不重新生成用户设置的范围。
obj.AddMember("Distribution", int(m_distribution), allocator);
obj.AddMember("DistributionInitialized", m_distributionInitialized, allocator);
obj.AddMember("Mean", m_mean, allocator);
obj.AddMember("StdDev", m_stdDev, allocator);
obj.AddMember("Peak", m_peak, allocator);
obj.AddMember("DistributionMinimum", m_distributionMin, allocator);
obj.AddMember("DistributionMaximum", m_distributionMax, allocator);
// nmDataAttribute 成员
obj.AddMember("ModelValue", m_modelValue.ToJsonValue(allocator), allocator);
obj.AddMember("Minimum", m_minValue.ToJsonValue(allocator), allocator);
obj.AddMember("Maximum", m_maxValue.ToJsonValue(allocator), allocator);
// values 数组
rapidjson::Value arr(rapidjson::kArrayType);
for (int i = 0; i < m_values.size(); ++i) {
arr.PushBack(m_values.at(i), allocator);
}
obj.AddMember("Values", arr, allocator);
return obj;
}
// 从 RapidJSON Value 反序列化到 VariableSampling
void nmDataSensitive::VariableSampling::FromJsonValue(const rapidjson::Value& jsonValue)
{
// 旧文件未保存分布时保留默认状态,由分析窗口按模型值初始化。
*this = VariableSampling();
if(jsonValue.HasMember("Distribution") && jsonValue["Distribution"].IsInt()) {
const int distribution = jsonValue["Distribution"].GetInt();
if(distribution >= NORMAL && distribution <= TRIANGULAR)
m_distribution = static_cast<Distribution>(distribution);
}
m_distributionInitialized = jsonValue.HasMember("DistributionInitialized") &&
jsonValue["DistributionInitialized"].IsBool() && jsonValue["DistributionInitialized"].GetBool();
m_mean = sensitiveNumberFromJson(jsonValue, "Mean", m_mean);
m_stdDev = sensitiveNumberFromJson(jsonValue, "StdDev", m_stdDev);
m_peak = sensitiveNumberFromJson(jsonValue, "Peak", m_peak);
m_distributionMin = sensitiveNumberFromJson(jsonValue, "DistributionMinimum", m_distributionMin);
m_distributionMax = sensitiveNumberFromJson(jsonValue, "DistributionMaximum", m_distributionMax);
if (jsonValue.HasMember("VarGroup") && jsonValue["VarGroup"].IsString()) {
m_varGroup = QString::fromUtf8(jsonValue["VarGroup"].GetString());
}
if (jsonValue.HasMember("VarName") && jsonValue["VarName"].IsString()) {
m_varName = QString::fromUtf8(jsonValue["VarName"].GetString());
}
if (jsonValue.HasMember("Enabled") && jsonValue["Enabled"].IsBool()) {
m_enabled = jsonValue["Enabled"].GetBool();
}
if (jsonValue.HasMember("Mode") && jsonValue["Mode"].IsInt()) {
int m = jsonValue["Mode"].GetInt();
m_mode = (m == MODE_MANUAL) ? MODE_MANUAL : MODE_AUTOMATIC;
}
if (jsonValue.HasMember("Log") && jsonValue["Log"].IsBool()) {
m_log = jsonValue["Log"].GetBool();
}
if (jsonValue.HasMember("Number") && jsonValue["Number"].IsInt()) {
m_number = jsonValue["Number"].GetInt();
}
if (jsonValue.HasMember("ModelValue") && jsonValue["ModelValue"].IsObject()) {
m_modelValue.FromJsonValue(jsonValue["ModelValue"]);
}
if (jsonValue.HasMember("Minimum") && jsonValue["Minimum"].IsObject()) {
m_minValue.FromJsonValue(jsonValue["Minimum"]);
}
if (jsonValue.HasMember("Maximum") && jsonValue["Maximum"].IsObject()) {
m_maxValue.FromJsonValue(jsonValue["Maximum"]);
}
m_values.clear();
if (jsonValue.HasMember("Values") && jsonValue["Values"].IsArray()) {
const rapidjson::Value& arr = jsonValue["Values"];
for (rapidjson::SizeType i = 0; i < arr.Size(); ++i) {
if (arr[i].IsNumber()) {
m_values.append(arr[i].GetDouble());
}
}
}
}
// varGroup
QString nmDataSensitive::VariableSampling::getVarGroup() const
{
return m_varGroup;
}
void nmDataSensitive::VariableSampling::setVarGroup(const QString& group)
{
m_varGroup = group;
}
// varName
QString nmDataSensitive::VariableSampling::getVarName() const
{
return m_varName;
}
void nmDataSensitive::VariableSampling::setVarName(const QString& name)
{
m_varName = name;
}
// enabled
bool nmDataSensitive::VariableSampling::getEnabled() const
{
return m_enabled;
}
void nmDataSensitive::VariableSampling::setEnabled(bool e)
{
m_enabled = e;
}
// mode
nmDataSensitive::VariableSampling::Mode
nmDataSensitive::VariableSampling::getMode() const
{
return m_mode;
}
void nmDataSensitive::VariableSampling::setMode(Mode m)
{
m_mode = m;
}
// log
bool nmDataSensitive::VariableSampling::getLog() const
{
return m_log;
}
void nmDataSensitive::VariableSampling::setLog(bool l)
{
m_log = l;
}
// number
int nmDataSensitive::VariableSampling::getNumber() const
{
return m_number;
}
void nmDataSensitive::VariableSampling::setNumber(int n)
{
m_number = n;
}
// modelValue
nmDataAttribute& nmDataSensitive::VariableSampling::getModelValue()
{
return m_modelValue;
}
void nmDataSensitive::VariableSampling::setModelValue(const nmDataAttribute& attr)
{
m_modelValue = attr;
}
// minValue
nmDataAttribute& nmDataSensitive::VariableSampling::getMinValue()
{
return m_minValue;
}
void nmDataSensitive::VariableSampling::setMinValue(const nmDataAttribute& attr)
{
m_minValue = attr;
}
// maxValue
nmDataAttribute& nmDataSensitive::VariableSampling::getMaxValue()
{
return m_maxValue;
}
void nmDataSensitive::VariableSampling::setMaxValue(const nmDataAttribute& attr)
{
m_maxValue = attr;
}
// values
QVector<double>& nmDataSensitive::VariableSampling::getValues()
{
return m_values;
}
void nmDataSensitive::VariableSampling::setValues(const QVector<double>& vals)
{
m_values = vals;
}
QVector<double> nmDataSensitive::VariableSampling::samplingValues() const
{
// 表格显示精度不参与采样,包含零或负数的对数范围按线性处理。
if (m_mode == MODE_MANUAL) return m_values;
QVector<double> values;
if (m_number <= 0) return values;
const double minimum = m_minValue.getValue().toDouble();
const double maximum = m_maxValue.getValue().toDouble();
const bool logarithmic = m_log && minimum > 0.0 && maximum > 0.0;
for (int i = 0; i < m_number; ++i) {
const double t = m_number == 1 ? 0.0 : double(i) / (m_number - 1);
double value = logarithmic
? std::exp((1.0 - t) * std::log(minimum) + t * std::log(maximum))
: (1.0 - t) * minimum + t * maximum;
if (m_number == 1) value = m_modelValue.getValue().toDouble();
else if (i == 0) value = minimum;
else if (i == m_number - 1) value = maximum;
values.append(value);
}
return values;
}
bool nmDataSensitive::VariableSampling::distributionQuantile(double probability, double& value, bool limitToPhysicalRange) const
{
// 六种分布都由分位数取值;对数正态的均值、标准差仍使用参数原单位。
if(!qIsFinite(probability) || probability <= 0.0 || probability >= 1.0) return false;
const bool specifiedBounds = m_distribution >= TRUNCATED_NORMAL;
const bool bounded = specifiedBounds || limitToPhysicalRange;
const double low = specifiedBounds ? m_distributionMin : m_modelValue.getMinValue().toDouble();
const double high = specifiedBounds ? m_distributionMax : m_modelValue.getMaxValue().toDouble();
if(bounded && (!qIsFinite(low) || !qIsFinite(high) || low >= high)) return false;
if(m_distribution == UNIFORM) value = low + probability * (high - low);
else if(m_distribution == TRIANGULAR) {
if(!qIsFinite(m_peak) || m_peak < low || m_peak > high) return false;
const double split = (m_peak - low) / (high - low);
value = probability < split
? low + std::sqrt(probability * (high - low) * (m_peak - low))
: high - std::sqrt((1.0 - probability) * (high - low) * (high - m_peak));
} else {
if(m_distribution < NORMAL || m_distribution > TRUNCATED_LOGNORMAL ||
!qIsFinite(m_mean) || !qIsFinite(m_stdDev) || m_stdDev <= 0.0) return false;
const bool logarithmic = m_distribution == LOGNORMAL || m_distribution == TRUNCATED_LOGNORMAL;
if(logarithmic && (m_mean <= 0.0 || (bounded && low < 0.0))) return false;
double mean = m_mean, sigma = m_stdDev;
if(logarithmic) {
const double ratio = m_stdDev / m_mean;
const double variance = std::log(1.0 + ratio * ratio);
mean = std::log(m_mean) - variance * 0.5;
sigma = std::sqrt(variance);
}
if(!qIsFinite(mean) || !qIsFinite(sigma) || sigma <= 0.0) return false;
double quantile = 0.0;
if(bounded) {
// 在允许范围内重新分配概率后取样,保持拉丁超立方分层,不裁剪或重抽越界样本。
double zLow = logarithmic && low == 0.0 ? -std::numeric_limits<double>::infinity()
: ((logarithmic ? std::log(low) : low) - mean) / sigma;
double zHigh = ((logarithmic ? std::log(high) : high) - mean) / sigma;
// 右尾利用正态分布对称性转到左尾,避免两个接近 1 的累积概率相减。
const bool reflected = zLow > 0.0;
if(reflected) { const double previousLow = zLow; zLow = -zHigh; zHigh = -previousLow; }
const double a = normalCdf(zLow), b = normalCdf(zHigh);
if(!(b > a)) return false;
probability = a + (reflected ? 1.0 - probability : probability) * (b - a);
// 超过 40 个标准差时双精度尾概率已下溢;缩小二分区间,同时保留实际边界。
quantile = normalQuantile(probability, qMax(zLow, -40.0), qMin(zHigh, 40.0));
if(reflected) quantile = -quantile;
} else quantile = normalQuantile(probability);
value = mean + sigma * quantile;
if(logarithmic) value = std::exp(value);
}
return qIsFinite(value);
}
nmDataSensitive::nmDataSensitive()
: m_calcType(CALC_DETERMINISTIC)
, m_improveEnabled(false)
, m_totalModelCount(0)
, m_usesBaseUnits(false)
, m_variables()
, m_monteCarloCount(50)
{
}
// 拷贝构造
nmDataSensitive::nmDataSensitive(const nmDataSensitive& other)
{
*this = other;
}
// 赋值运算符
nmDataSensitive& nmDataSensitive::operator=(const nmDataSensitive& other)
{
if (this != &other) {
m_calcType = other.m_calcType;
m_improveEnabled = other.m_improveEnabled;
m_improveSettings = other.m_improveSettings;
m_totalModelCount = other.m_totalModelCount;
m_usesBaseUnits = other.m_usesBaseUnits;
m_variables = other.m_variables;
m_monteCarloCount = other.m_monteCarloCount;
m_monteCarloValues = other.m_monteCarloValues;
m_results = other.m_results;
}
return *this;
}
nmDataSensitive::~nmDataSensitive()
{
}
// 序列化 nmDataSensitive 为 RapidJSON Value
rapidjson::Value nmDataSensitive::ToJsonValue(
rapidjson::Document::AllocatorType& allocator) const
{
rapidjson::Value root(rapidjson::kObjectType);
// 旧配置中的采样值是显示单位,新配置统一保存业务基准单位。
root.AddMember("UsesBaseUnits", m_usesBaseUnits, allocator);
// 顶部 Calculation type
root.AddMember("CalculationType",
static_cast<int>(m_calcType),
allocator);
// 与曲线一起保存采样和拟合配置,重开项目后仍可查看同一批预览。
root.AddMember("ImproveEnabled", m_improveEnabled, allocator);
root.AddMember("MonteCarloCount", m_monteCarloCount, allocator);
root.AddMember("MonteCarloValues", sensitiveCurvesToJson(m_monteCarloValues, allocator), allocator);
if(!m_improveSettings.isNull()) {
const nmSensitivityFitSettings& settings = *m_improveSettings;
rapidjson::Value fitting(rapidjson::kObjectType);
fitting.AddMember("Confirmed", settings.confirmed, allocator);
fitting.AddMember("Parameters", settings.parameters.ToJsonValue(allocator), allocator);
fitting.AddMember("TargetWell", rapidjson::Value(settings.targetWell.toUtf8().constData(), allocator).Move(), allocator);
fitting.AddMember("InitialValues", sensitiveNumbersToJson(settings.initialValues, allocator), allocator);
root.AddMember("ImproveSettings", fitting, allocator);
}
// 总模型数
root.AddMember("TotalModelCount",
m_totalModelCount,
allocator);
// 变量数组
rapidjson::Value arr(rapidjson::kArrayType);
for (int i = 0; i < m_variables.size(); ++i) {
arr.PushBack(m_variables.at(i).ToJsonValue(allocator), allocator);
}
root.AddMember("Variables", arr, allocator);
// 批量结果随现有 Sensitive 节点保存,不增加独立成果文件。
if(!m_results.isNull()) root.AddMember("Results", m_results->ToJsonValue(allocator), allocator);
return root;
}
// 从 RapidJSON Value 反序列化数据到 nmDataSensitive
void nmDataSensitive::FromJsonValue(const rapidjson::Value& jsonValue)
{
// 先恢复默认值,旧成果缺失新增字段时不沿用上一项目的数据。
*this = nmDataSensitive();
if(!jsonValue.IsObject()) return;
m_improveEnabled = jsonValue.HasMember("ImproveEnabled") &&
jsonValue["ImproveEnabled"].IsBool() && jsonValue["ImproveEnabled"].GetBool();
if(jsonValue.HasMember("MonteCarloCount") && jsonValue["MonteCarloCount"].IsInt() &&
jsonValue["MonteCarloCount"].GetInt() > 0)
m_monteCarloCount = jsonValue["MonteCarloCount"].GetInt();
// 旧项目没有 Results 时保持空结果。
m_results.clear();
if(jsonValue.HasMember("Results") && jsonValue["Results"].IsObject()) {
m_results = QSharedPointer<BatchResult>(new BatchResult());
m_results->FromJsonValue(jsonValue["Results"]);
}
m_usesBaseUnits = jsonValue.HasMember("UsesBaseUnits") &&
jsonValue["UsesBaseUnits"].IsBool() && jsonValue["UsesBaseUnits"].GetBool();
if (jsonValue.HasMember("CalculationType") && jsonValue["CalculationType"].IsInt()) {
int t = jsonValue["CalculationType"].GetInt();
if (t >= CALC_DETERMINISTIC && t <= CALC_MONTE_CARLO)
m_calcType = static_cast<CalculationType>(t);
}
if (jsonValue.HasMember("TotalModelCount") && jsonValue["TotalModelCount"].IsUint64()) {
m_totalModelCount = jsonValue["TotalModelCount"].GetUint64();
}
m_variables.clear();
if (jsonValue.HasMember("Variables") && jsonValue["Variables"].IsArray()) {
const rapidjson::Value& arr = jsonValue["Variables"];
for (rapidjson::SizeType i = 0; i < arr.Size(); ++i) {
if (!arr[i].IsObject())
continue;
VariableSampling vs;
vs.FromJsonValue(arr[i]);
m_variables.append(vs);
}
}
// 恢复完整样本,列数必须与参数表对应,避免预览或生成时错用参数。
m_monteCarloValues = sensitiveCurvesFromJson(jsonValue, "MonteCarloValues");
bool validSamples = m_monteCarloValues.size() == m_monteCarloCount;
for(int row = 0; validSamples && row < m_monteCarloValues.size(); ++row) {
const QVector<double>& values = m_monteCarloValues[row];
validSamples = values.size() == m_variables.size();
for(int column = 0; validSamples && column < values.size(); ++column)
validSamples = qIsFinite(values[column]);
}
if(!validSamples) clearMonteCarloSamples();
if(jsonValue.HasMember("ImproveSettings") && jsonValue["ImproveSettings"].IsObject()) {
const rapidjson::Value& fitting = jsonValue["ImproveSettings"];
QSharedPointer<nmSensitivityFitSettings> settings(new nmSensitivityFitSettings());
settings->targetWell = sensitiveTextFromJson(fitting, "TargetWell");
const bool hasParameters = fitting.HasMember("Parameters") && fitting["Parameters"].IsObject();
if(hasParameters) settings->parameters.FromJsonValue(fitting["Parameters"]);
if(fitting.HasMember("InitialValues")) sensitiveNumbersFromJson(fitting["InitialValues"], settings->initialValues);
bool valid = hasParameters && settings->initialValues.size() == 6 && !settings->targetWell.isEmpty();
for(int i = 0; valid && i < settings->initialValues.size(); ++i)
valid = qIsFinite(settings->initialValues[i]);
settings->confirmed = valid && fitting.HasMember("Confirmed") &&
fitting["Confirmed"].IsBool() && fitting["Confirmed"].GetBool();
if(settings->confirmed) m_improveSettings = settings;
}
}
// === Calculation type ===
nmDataSensitive::CalculationType nmDataSensitive::getCalculationType() const
{
return m_calcType;
}
void nmDataSensitive::setCalculationType(CalculationType type)
{
// 切换计算方法后不再沿用上一种方法的预览样本。
if(type != m_calcType) clearMonteCarloSamples();
m_calcType = type;
}
bool nmDataSensitive::isImproveEnabled() const
{
// 开关独立于拟合设置,切回第三种后仍可保留上次确认的设置。
return m_improveEnabled;
}
void nmDataSensitive::setImproveEnabled(bool enabled)
{
// 只记录是否启用改善,不改变 Monte-Carlo 的采样方式。
m_improveEnabled = enabled;
}
QSharedPointer<const nmSensitivityFitSettings> nmDataSensitive::getImproveSettings() const
{
// 设置按只读引用共享,窗口重新打开时恢复,编辑时另建副本。
return m_improveSettings;
}
void nmDataSensitive::setImproveSettings(const QSharedPointer<const nmSensitivityFitSettings>& settings)
{
// 只接收拟合窗口确认后的完整配置。
m_improveSettings = settings;
}
// === 总模型数 ===
quint64 nmDataSensitive::getTotalModelCount() const
{
return m_totalModelCount;
}
void nmDataSensitive::setTotalModelCount(quint64 count)
{
m_totalModelCount = count;
}
quint64 nmDataSensitive::calculateModelCount(bool* overflow) const
{
// 未选择变量时为零;两种确定性方法分别累计和、积。
if (overflow) *overflow = false;
if(m_calcType == CALC_MONTE_CARLO) {
for(int i = 0; i < m_variables.size(); ++i)
if(m_variables[i].m_enabled) return qMax(0, m_monteCarloCount);
return 0;
}
quint64 count = 0;
bool selected = false;
for (int i = 0; i < m_variables.size(); ++i) {
const VariableSampling& variable = m_variables[i];
if (!variable.m_enabled) continue;
const quint64 number = variable.m_number > 0 ? variable.m_number : 0;
if (!number) return 0;
const quint64 limit = (std::numeric_limits<quint64>::max)();
if (m_calcType == CALC_DETERMINISTIC) {
if (count > limit - number) { if (overflow) *overflow = true; return 0; }
count += number;
} else {
if (!selected) count = 1;
if (count > limit / number) { if (overflow) *overflow = true; return 0; }
count *= number;
}
selected = true;
}
return count;
}
bool nmDataSensitive::getRunValues(quint64 runIndex, QVector<double>& values) const
{
// 按编号解码组合,不预先展开笛卡尔积,也不删除重复的基准组合。
bool overflow = false;
const quint64 count = calculateModelCount(&overflow);
if (overflow || runIndex >= count) return false;
if(m_calcType == CALC_MONTE_CARLO) {
if(m_monteCarloValues.size() != int(count)) return false;
values = m_monteCarloValues[int(runIndex)];
return true;
}
values.clear();
for (int i = 0; i < m_variables.size(); ++i)
values.append(m_variables[i].m_modelValue.getValue().toDouble());
if (m_calcType == CALC_DETERMINISTIC) {
for (int i = 0; i < m_variables.size(); ++i) {
const VariableSampling& variable = m_variables[i];
if (!variable.m_enabled) continue;
const QVector<double> samples = variable.samplingValues();
if (samples.size() != variable.m_number) return false;
if (runIndex < quint64(samples.size())) {
values[i] = samples[int(runIndex)];
return true;
}
runIndex -= samples.size();
}
} else {
// Saphir 按第一个选中参数最快变化的顺序列出多变量组合。
for (int i = 0; i < m_variables.size(); ++i) {
const VariableSampling& variable = m_variables[i];
if (!variable.m_enabled) continue;
const QVector<double> samples = variable.samplingValues();
if (samples.isEmpty() || samples.size() != variable.m_number) return false;
values[i] = samples[int(runIndex % samples.size())];
runIndex /= samples.size();
}
return true;
}
return false;
}
int nmDataSensitive::getMonteCarloCount() const
{
// Monte-Carlo 总次数由用户指定,不使用单参数次数的和或积。
return m_monteCarloCount;
}
void nmDataSensitive::setMonteCarloCount(int count)
{
// 次数变化后旧样本失效,下次预览或生成时重新抽样。
if(count != m_monteCarloCount) clearMonteCarloSamples();
m_monteCarloCount = count;
}
void nmDataSensitive::clearMonteCarloSamples()
{
// 只丢弃待计算样本,已完成批次的结果仍由结果窗口管理。
m_monteCarloValues.clear();
}
bool nmDataSensitive::resampleMonteCarlo(QString* error)
{
// 各参数的累积概率均分为 N 层,每层抽一点,再独立打乱,组成 N 组参数。
if(error) error->clear();
m_monteCarloValues.clear();
const int count = int(calculateModelCount());
if(m_calcType != CALC_MONTE_CARLO || count <= 0) return false;
std::mt19937 random(uint(QDateTime::currentMSecsSinceEpoch()) ^ uint(qrand()));
std::uniform_real_distribution<double> uniform(0.0, 1.0);
QVector<QVector<double> > runs(count, QVector<double>(m_variables.size()));
for(int column = 0; column < m_variables.size(); ++column) {
const VariableSampling& variable = m_variables[column];
QVector<double> samples(count, variable.m_modelValue.getValue().toDouble());
if(variable.m_enabled) {
for(int row = 0; row < count; ++row) {
const double jitter = qBound(1e-12, uniform(random), 1.0 - 1e-12);
if(!variable.distributionQuantile((row + jitter) / count, samples[row])) {
if(error) *error = tr("%1: invalid distribution parameters.").arg(variable.m_varName);
return false;
}
}
for(int row = count - 1; row > 0; --row) {
std::uniform_int_distribution<int> index(0, row);
qSwap(samples[row], samples[index(random)]);
}
}
for(int row = 0; row < count; ++row) runs[row][column] = samples[row];
}
m_monteCarloValues.swap(runs);
return true;
}
bool nmDataSensitive::usesBaseUnits() const
{
return m_usesBaseUnits;
}
QSharedPointer<nmDataSensitive::BatchResult> nmDataSensitive::getResults() const
{
// 曲线容器按共享引用交给界面,重复打开不复制整批结果。
return m_results;
}
void nmDataSensitive::setResults(const QSharedPointer<BatchResult>& results)
{
// 结果仅在线程结束后发布,后续参数编辑不会重写冻结的批次配置。
m_results = results;
}
void nmDataSensitive::setUsesBaseUnits(bool enabled)
{
// 单位标记随配置保存,避免再次打开时重复转换。
m_usesBaseUnits = enabled;
}
// === 变量列表 ===
QList<nmDataSensitive::VariableSampling>& nmDataSensitive::getVariables()
{
return m_variables;
}
const QList<nmDataSensitive::VariableSampling>& nmDataSensitive::getVariables() const
{
// 求解任务只读取采样配置,不取得可写参数列表。
return m_variables;
}
void nmDataSensitive::setVariables(const QList<VariableSampling>& vars)
{
m_variables = vars;
}
int nmDataSensitive::getVariableCount() const
{
return m_variables.size();
}
nmDataSensitive::VariableSampling& nmDataSensitive::getVariable(int index)
{
return m_variables[index];
}