#2659 Ensemble statistics. First version

This commit is contained in:
Bjørn Erik Jensen
2018-06-07 10:02:54 +02:00
parent b73c8e0565
commit 69accdc4c6
10 changed files with 632 additions and 7 deletions

View File

@@ -25,10 +25,14 @@
#include "RifReaderEclipseSummary.h"
#include "RigStatisticsMath.h"
#include "RigTimeHistoryCurveMerger.h"
#include "RimEnsembleCurveFilter.h"
#include "RimEnsembleCurveFilterCollection.h"
#include "RimEnsembleCurveSetCollection.h"
#include "RimEnsembleCurveSetColorManager.h"
#include "RimEnsembleStatistics.h"
#include "RimProject.h"
#include "RimRegularLegendConfig.h"
#include "RimSummaryAddress.h"
@@ -137,6 +141,10 @@ RimEnsembleCurveSet::RimEnsembleCurveSet()
CAF_PDM_InitFieldNoDefault(&m_curveFilters, "CurveFilters", "Curve Filters", "", "", "");
m_curveFilters = new RimEnsembleCurveFilterCollection();
CAF_PDM_InitFieldNoDefault(&m_statistics, "Statistics", "Statistics", "", "", "");
m_statistics = new RimEnsembleStatistics();
m_statistics.uiCapability()->setUiTreeHidden(true);
CAF_PDM_InitField(&m_userDefinedName, "UserDefinedName", QString("Ensemble Curve Set"), "Curve Set Name", "", "", "");
CAF_PDM_InitFieldNoDefault(&m_autoGeneratedName, "AutoGeneratedName", "Curve Set Name", "", "", "");
@@ -153,6 +161,9 @@ RimEnsembleCurveSet::RimEnsembleCurveSet()
m_qwtPlotCurveForLegendText = new QwtPlotCurve;
m_qwtPlotCurveForLegendText->setLegendAttribute(QwtPlotCurve::LegendShowSymbol, true);
m_ensembleStatCase.reset(new RimEnsembleStatisticsCase(this));
m_ensembleStatCase->createSummaryReaderInterface();
}
//--------------------------------------------------------------------------------------------------
@@ -200,6 +211,7 @@ void RimEnsembleCurveSet::loadDataAndUpdate(bool updateParentPlot)
m_yValuesUiFilterResultSelection = m_yValuesCurveVariable->address();
updateAllCurves();
updateStatisticsCurves(true);
if (updateParentPlot)
{
@@ -288,9 +300,45 @@ std::vector<RimSummaryCurve*> RimEnsembleCurveSet::curves() const
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimEnsembleCurveSet::deleteAllCurves()
void RimEnsembleCurveSet::deleteEnsembleCurves()
{
m_curves.deleteAllChildObjects();
std::vector<int> curvesIndexesToDelete;
for (int c = 0; c < m_curves.size(); c++)
{
RimSummaryCurve* curve = m_curves[c];
if (curve->summaryAddressY().category() != RifEclipseSummaryAddress::SUMMARY_ENSEMBLE_STATISTICS)
curvesIndexesToDelete.push_back(c);
}
while (curvesIndexesToDelete.size() > 0)
{
int currIndex = curvesIndexesToDelete.back();
delete m_curves[currIndex];
m_curves.erase(currIndex);
curvesIndexesToDelete.pop_back();
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimEnsembleCurveSet::deleteStatisticsCurves()
{
std::vector<int> curvesIndexesToDelete;
for (int c = 0; c < m_curves.size(); c++)
{
RimSummaryCurve* curve = m_curves[c];
if (curve->summaryAddressY().category() == RifEclipseSummaryAddress::SUMMARY_ENSEMBLE_STATISTICS)
curvesIndexesToDelete.push_back(c);
}
while (curvesIndexesToDelete.size() > 0)
{
int currIndex = curvesIndexesToDelete.back();
delete m_curves[currIndex];
m_curves.erase(currIndex);
curvesIndexesToDelete.pop_back();
}
}
//--------------------------------------------------------------------------------------------------
@@ -380,6 +428,7 @@ void RimEnsembleCurveSet::fieldChangedByUi(const caf::PdmFieldHandle* changedFie
m_yValuesCurveVariable->setAddress(m_yValuesUiFilterResultSelection());
updateAllCurves();
updateStatisticsCurves(true);
updateTextInPlot = true;
}
@@ -388,7 +437,8 @@ void RimEnsembleCurveSet::fieldChangedByUi(const caf::PdmFieldHandle* changedFie
// Empty address cache
m_allAddressesCache.clear();
updateAllCurves();
updateStatisticsCurves(true);
updateTextInPlot = true;
}
else if (changedField == &m_color)
@@ -524,6 +574,9 @@ void RimEnsembleCurveSet::defineUiOrdering(QString uiConfigName, caf::PdmUiOrder
}
}
caf::PdmUiGroup* statGroup = uiOrdering.addNewGroup("Statistics");
m_statistics->defineUiOrdering(uiConfigName, *statGroup);
uiOrdering.skipRemainingFields(true);
}
@@ -818,7 +871,7 @@ void RimEnsembleCurveSet::updateAllCurves()
firstAncestorOrThisOfType(plot);
CVF_ASSERT(plot);
deleteAllCurves();
deleteEnsembleCurves();
m_qwtPlotCurveForLegendText->detach();
RimSummaryCaseCollection* group = m_yValuesSummaryGroup();
@@ -862,6 +915,57 @@ void RimEnsembleCurveSet::updateAllCurves()
updateCurveColors();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimEnsembleCurveSet::updateStatisticsCurves(bool calculate = true)
{
using SAddr = RifEclipseSummaryAddress;
if (calculate)
{
// Calculate
m_ensembleStatCase->calculate(m_yValuesSummaryGroup, m_yValuesCurveVariable->address());
}
RimSummaryPlot* plot = nullptr;
firstAncestorOrThisOfType(plot);
CVF_ASSERT(plot);
std::vector<RifEclipseSummaryAddress> addresses;
if (m_statistics->isActive())
{
if (m_statistics->showP10Curve()) addresses.push_back(SAddr::ensembleStatisticsAddress(ENSEMBLE_STAT_P10_QUANTITY_NAME));
if (m_statistics->showP50Curve()) addresses.push_back(SAddr::ensembleStatisticsAddress(ENSEMBLE_STAT_P50_QUANTITY_NAME));
if (m_statistics->showP90Curve()) addresses.push_back(SAddr::ensembleStatisticsAddress(ENSEMBLE_STAT_P90_QUANTITY_NAME));
if (m_statistics->showMeanCurve()) addresses.push_back(SAddr::ensembleStatisticsAddress(ENSEMBLE_STAT_MEAN_QUANTITY_NAME));
}
deleteStatisticsCurves();
for (auto address : addresses)
{
auto curve = new RimSummaryCurve();
curve->setParentQwtPlotNoReplot(plot->qwtPlot());
m_curves.push_back(curve);
curve->setSymbol(RimPlotCurve::SYMBOL_ELLIPSE);
curve->setSymbolSkipDinstance(50);
curve->setLineStyle(RimPlotCurve::STYLE_SOLID);
curve->setSummaryCaseY(m_ensembleStatCase.get());
curve->setSummaryAddressY(address);
curve->setZOrder(1000);
curve->updateCurveVisibility(false);
curve->loadDataAndUpdate(false);
}
if (plot->qwtPlot())
{
plot->qwtPlot()->updateLegend();
plot->qwtPlot()->replot();
plot->updateAxes();
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@@ -1035,3 +1139,244 @@ void RimEnsembleCurveSet::updateLegendMappingMode()
break;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RifEnsembleStatisticsReader::RifEnsembleStatisticsReader(RimEnsembleStatisticsCase* ensStatCase)
{
CVF_ASSERT(ensStatCase);
m_ensembleStatCase = ensStatCase;
m_allResultAddresses = std::vector<RifEclipseSummaryAddress>(
{
RifEclipseSummaryAddress::ensembleStatisticsAddress(ENSEMBLE_STAT_P10_QUANTITY_NAME),
RifEclipseSummaryAddress::ensembleStatisticsAddress(ENSEMBLE_STAT_P50_QUANTITY_NAME),
RifEclipseSummaryAddress::ensembleStatisticsAddress(ENSEMBLE_STAT_P90_QUANTITY_NAME),
RifEclipseSummaryAddress::ensembleStatisticsAddress(ENSEMBLE_STAT_MEAN_QUANTITY_NAME)
});
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<time_t>& RifEnsembleStatisticsReader::timeSteps(const RifEclipseSummaryAddress& resultAddress) const
{
return m_ensembleStatCase->timeSteps();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RifEnsembleStatisticsReader::values(const RifEclipseSummaryAddress& resultAddress, std::vector<double>* values) const
{
if (!validateAddress(resultAddress)) return false;
const std::vector<double>* sourceData = nullptr;
if (resultAddress.quantityName() == ENSEMBLE_STAT_P10_QUANTITY_NAME) sourceData = &m_ensembleStatCase->p10();
else if (resultAddress.quantityName() == ENSEMBLE_STAT_P50_QUANTITY_NAME) sourceData = &m_ensembleStatCase->p50();
else if (resultAddress.quantityName() == ENSEMBLE_STAT_P90_QUANTITY_NAME) sourceData = &m_ensembleStatCase->p90();
else if (resultAddress.quantityName() == ENSEMBLE_STAT_MEAN_QUANTITY_NAME) sourceData = &m_ensembleStatCase->mean();
if (!sourceData) return false;
values->clear();
values->reserve(sourceData->size());
for (auto val : *sourceData) values->push_back(val);
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::string RifEnsembleStatisticsReader::unitName(const RifEclipseSummaryAddress& resultAddress) const
{
return "(RifEnsembleStatisticsReader::unitName)";
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RifEnsembleStatisticsReader::validateAddress(const RifEclipseSummaryAddress& address) const
{
return address.category() == RifEclipseSummaryAddress::SUMMARY_ENSEMBLE_STATISTICS &&
!address.quantityName().empty();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimEnsembleStatisticsCase::RimEnsembleStatisticsCase(RimEnsembleCurveSet* curveSet)
{
m_curveSet = curveSet;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<time_t>& RimEnsembleStatisticsCase::timeSteps() const
{
return m_timeSteps;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RimEnsembleStatisticsCase::p10() const
{
return m_p10Data;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RimEnsembleStatisticsCase::p50() const
{
return m_p50Data;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RimEnsembleStatisticsCase::p90() const
{
return m_p90Data;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RimEnsembleStatisticsCase::mean() const
{
return m_meanData;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimEnsembleStatisticsCase::caseName()
{
return "Ensemble Statistics";
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimEnsembleStatisticsCase::createSummaryReaderInterface()
{
m_statisticsReader.reset(new RifEnsembleStatisticsReader(this));
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RifSummaryReaderInterface* RimEnsembleStatisticsCase::summaryReader()
{
return m_statisticsReader.get();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const RimEnsembleCurveSet* RimEnsembleStatisticsCase::curveSet() const
{
return m_curveSet;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimEnsembleStatisticsCase::calculate(const RimSummaryCaseCollection* ensemble, const RifEclipseSummaryAddress& inputAddress)
{
RigTimeHistoryCurveMerger curveMerger;
int caseCount = (int)ensemble->allSummaryCases().size();
{
for (const auto& sumCase : ensemble->allSummaryCases())
{
const auto& reader = sumCase->summaryReader();
if (reader)
{
std::vector<time_t> timeSteps = reader->timeSteps(inputAddress);
std::vector<double> values;
reader->values(inputAddress, &values);
curveMerger.addCurveData(values, timeSteps);
}
}
curveMerger.computeInterpolatedValues();
}
const std::vector<time_t>& allTimeSteps = curveMerger.allTimeSteps();
std::vector<std::vector<double>> allValues;
{
for (int c = 0; c < caseCount; c++)
{
allValues.push_back(curveMerger.interpolatedCurveValuesForAllTimeSteps(c));
}
}
clearData();
m_timeSteps = allTimeSteps;
for (int t = 0; t < (int)allTimeSteps.size(); t++)
{
std::vector<double> valuesForTimeSteps;
valuesForTimeSteps.reserve(caseCount);
for (int c = 0; c < caseCount; c++)
{
valuesForTimeSteps.push_back(allValues[c][t]);
}
{
double min, max, range, mean, stdev;
RigStatisticsMath::calculateBasicStatistics(valuesForTimeSteps, &min, &max, nullptr, &range, &mean, &stdev);
std::vector<size_t> histogram;
RigHistogramCalculator histCalc(min, max, 100, &histogram);
histCalc.addData(valuesForTimeSteps);
double p10, p50, p90;
p10 = histCalc.calculatePercentil(0.1);
p50 = histCalc.calculatePercentil(0.5);
p90 = histCalc.calculatePercentil(0.9);
m_p10Data.push_back(p10);
m_p50Data.push_back(p50);
m_p90Data.push_back(p90);
m_meanData.push_back(mean);
}
}
m_addressUsedInLastCalculation = inputAddress;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimEnsembleStatisticsCase::calculate()
{
auto inputAddress = m_curveSet->summaryAddress();
auto ensemble = m_curveSet->summaryCaseCollection();
if (m_statisticsReader && ensemble && inputAddress.isValid())
{
calculate(ensemble, inputAddress);
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimEnsembleStatisticsCase::clearData()
{
m_timeSteps.clear();
m_p10Data.clear();
m_p50Data.clear();
m_p90Data.clear();
m_meanData.clear();
m_addressUsedInLastCalculation = RifEclipseSummaryAddress();
}