#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("WellCode", rapidjson::Value(m_wellCode.toUtf8().constData(), allocator).Move(), allocator);

	// 字符串字段
	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();
    m_wellCode = jsonValue.HasMember("WellCode") && jsonValue["WellCode"].IsString()
        ? QString::fromUtf8(jsonValue["WellCode"].GetString()) : QString();
    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::displayName() const
{
    // 预览与曲线标注同时显示井名，避免不同井的同名参数混淆。
    return m_wellCode.isEmpty() ? m_varName : m_varGroup + " / " + m_varName;
}

QString nmDataSensitive::VariableSampling::parameterKey() const
{
    // 参数互斥按所属井判断，共用储层参数仍只使用参数名。
    return m_wellCode.isEmpty() ? m_varName : m_wellCode + "/" + m_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];
}