From d5e911d01b68c2b044ebc576ff8b1c095f1e3302 Mon Sep 17 00:00:00 2001 From: Magne Sjaastad Date: Mon, 18 Aug 2014 12:03:21 +0200 Subject: [PATCH] Added RigStatisticsDataCache and RigStatisticsCalculator Moved cache from RigCaseCellResultsData to RigStatisticsDataCache Moved statistics computations from RigCaseCellResultsData to RigStatisticsCalculator --- .../ReservoirDataModel/CMakeLists_files.cmake | 4 + .../RigCaseCellResultsData.cpp | 417 ++---------------- .../RigCaseCellResultsData.h | 19 +- .../RigStatisticsCalculator.cpp | 240 ++++++++++ .../RigStatisticsCalculator.h | 87 ++++ .../RigStatisticsDataCache.cpp | 205 +++++++++ .../RigStatisticsDataCache.h | 65 +++ .../RiaPropertyDataCommands.cpp | 4 +- 8 files changed, 659 insertions(+), 382 deletions(-) create mode 100644 ApplicationCode/ReservoirDataModel/RigStatisticsCalculator.cpp create mode 100644 ApplicationCode/ReservoirDataModel/RigStatisticsCalculator.h create mode 100644 ApplicationCode/ReservoirDataModel/RigStatisticsDataCache.cpp create mode 100644 ApplicationCode/ReservoirDataModel/RigStatisticsDataCache.h diff --git a/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake b/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake index 25f9be656d..78e6bbc52a 100644 --- a/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake +++ b/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake @@ -32,6 +32,8 @@ ${CEE_CURRENT_LIST_DIR}cvfGeometryTools.inl ${CEE_CURRENT_LIST_DIR}RigPipeInCellEvaluator.h ${CEE_CURRENT_LIST_DIR}RigResultAccessor2d.h ${CEE_CURRENT_LIST_DIR}RigTernaryResultAccessor2d.h +${CEE_CURRENT_LIST_DIR}RigStatisticsDataCache.h +${CEE_CURRENT_LIST_DIR}RigStatisticsCalculator.h ) set (SOURCE_GROUP_SOURCE_FILES @@ -58,6 +60,8 @@ ${CEE_CURRENT_LIST_DIR}RigFault.cpp ${CEE_CURRENT_LIST_DIR}RigNNCData.cpp ${CEE_CURRENT_LIST_DIR}cvfGeometryTools.cpp ${CEE_CURRENT_LIST_DIR}RigTernaryResultAccessor2d.cpp +${CEE_CURRENT_LIST_DIR}RigStatisticsDataCache.cpp +${CEE_CURRENT_LIST_DIR}RigStatisticsCalculator.cpp ) list(APPEND CODE_HEADER_FILES diff --git a/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.cpp b/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.cpp index ea168688b6..9438237e6e 100644 --- a/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.cpp +++ b/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.cpp @@ -17,10 +17,11 @@ ///////////////////////////////////////////////////////////////////////////////// #include "RigCaseCellResultsData.h" -#include "RifReaderInterface.h" -#include "RigMainGrid.h" +#include "RigMainGrid.h" +#include "RigStatisticsDataCache.h" #include "RigStatisticsMath.h" +#include "RigStatisticsCalculator.h" #include #include @@ -33,11 +34,8 @@ RigCaseCellResultsData::RigCaseCellResultsData(RigMainGrid* ownerGrid) { CVF_ASSERT(ownerGrid != NULL); m_ownerMainGrid = ownerGrid; - - m_combinedTransmissibilityResultIndex = cvf::UNDEFINED_SIZE_T; } - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -49,39 +47,9 @@ void RigCaseCellResultsData::setMainGrid(RigMainGrid* ownerGrid) //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RigCaseCellResultsData::minMaxCellScalarValues( size_t scalarResultIndex, double& min, double& max ) +void RigCaseCellResultsData::minMaxCellScalarValues(size_t scalarResultIndex, double& min, double& max) { - min = HUGE_VAL; - max = -HUGE_VAL; - - CVF_ASSERT(scalarResultIndex < resultCount()); - - // Extend array and cache vars - - if (scalarResultIndex >= m_maxMinValues.size() ) - { - m_maxMinValues.resize(scalarResultIndex+1, std::make_pair(HUGE_VAL, -HUGE_VAL)); - } - - if (m_maxMinValues[scalarResultIndex].first != HUGE_VAL) - { - min = m_maxMinValues[scalarResultIndex].first; - max = m_maxMinValues[scalarResultIndex].second; - - return; - } - - size_t i; - for (i = 0; i < timeStepCount(scalarResultIndex); i++) - { - double tsmin, tsmax; - minMaxCellScalarValues(scalarResultIndex, i, tsmin, tsmax); - if (tsmin < min) min = tsmin; - if (tsmax > max) max = tsmax; - } - - m_maxMinValues[scalarResultIndex].first = min; - m_maxMinValues[scalarResultIndex].second= max; + m_statisticsDataCache[scalarResultIndex]->minMaxCellScalarValues(min, max); } //-------------------------------------------------------------------------------------------------- @@ -89,80 +57,23 @@ void RigCaseCellResultsData::minMaxCellScalarValues( size_t scalarResultIndex, d //-------------------------------------------------------------------------------------------------- void RigCaseCellResultsData::minMaxCellScalarValues(size_t scalarResultIndex, size_t timeStepIndex, double& min, double& max) { - min = HUGE_VAL; - max = -HUGE_VAL; + m_statisticsDataCache[scalarResultIndex]->minMaxCellScalarValues(timeStepIndex, min, max); +} - CVF_ASSERT(scalarResultIndex < resultCount()); +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigCaseCellResultsData::posNegClosestToZero(size_t scalarResultIndex, double& pos, double& neg) +{ + m_statisticsDataCache[scalarResultIndex]->posNegClosestToZero(pos, neg); +} - if (timeStepIndex >= m_cellScalarResults[scalarResultIndex].size()) - { - return; - } - - if (scalarResultIndex >= m_maxMinValuesPrTs.size()) - { - m_maxMinValuesPrTs.resize(scalarResultIndex+1); - } - - if (timeStepIndex >= m_maxMinValuesPrTs[scalarResultIndex].size()) - { - m_maxMinValuesPrTs[scalarResultIndex].resize(timeStepIndex+1, std::make_pair(HUGE_VAL, -HUGE_VAL)); - } - - if (m_maxMinValuesPrTs[scalarResultIndex][timeStepIndex].first != HUGE_VAL) - { - min = m_maxMinValuesPrTs[scalarResultIndex][timeStepIndex].first; - max = m_maxMinValuesPrTs[scalarResultIndex][timeStepIndex].second; - - return; - } - - if (scalarResultIndex == m_combinedTransmissibilityResultIndex) - { - size_t tranX, tranY, tranZ; - if (!findTransmissibilityResults(tranX, tranY, tranZ)) return; - - double tranMin; - double tranMax; - - minMaxCellScalarValues(tranX, timeStepIndex, tranMin, tranMax); - min = CVF_MIN(tranMin, min); - max = CVF_MAX(tranMax, max); - - minMaxCellScalarValues(tranY, timeStepIndex, tranMin, tranMax); - min = CVF_MIN(tranMin, min); - max = CVF_MAX(tranMax, max); - - minMaxCellScalarValues(tranZ, timeStepIndex, tranMin, tranMax); - min = CVF_MIN(tranMin, min); - max = CVF_MAX(tranMax, max); - - return; - } - - std::vector& values = m_cellScalarResults[scalarResultIndex][timeStepIndex]; - - size_t i; - for (i = 0; i < values.size(); i++) - { - if (values[i] == HUGE_VAL) - { - continue; - } - - if (values[i] < min) - { - min = values[i]; - } - - if (values[i] > max) - { - max = values[i]; - } - } - - m_maxMinValuesPrTs[scalarResultIndex][timeStepIndex].first = min; - m_maxMinValuesPrTs[scalarResultIndex][timeStepIndex].second= max; +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigCaseCellResultsData::posNegClosestToZero(size_t scalarResultIndex, size_t timeStepIndex, double& pos, double& neg) +{ + m_statisticsDataCache[scalarResultIndex]->posNegClosestToZero(timeStepIndex, pos, neg); } //-------------------------------------------------------------------------------------------------- @@ -170,68 +81,15 @@ void RigCaseCellResultsData::minMaxCellScalarValues(size_t scalarResultIndex, si //-------------------------------------------------------------------------------------------------- const std::vector& RigCaseCellResultsData::cellScalarValuesHistogram(size_t scalarResultIndex) { - CVF_ASSERT(scalarResultIndex < resultCount()); - - // Extend array and cache vars - - if (scalarResultIndex >= m_histograms.size() ) - { - m_histograms.resize(resultCount()); - m_p10p90.resize(resultCount(), std::make_pair(HUGE_VAL, HUGE_VAL)); - } - - if (m_histograms[scalarResultIndex].size()) - { - return m_histograms[scalarResultIndex]; - - } - - double min; - double max; - size_t nBins = 100; - this->minMaxCellScalarValues( scalarResultIndex, min, max ); - RigHistogramCalculator histCalc(min, max, nBins, &m_histograms[scalarResultIndex]); - - if (scalarResultIndex == m_combinedTransmissibilityResultIndex) - { - size_t tranX, tranY, tranZ; - if (findTransmissibilityResults(tranX, tranY, tranZ)) - { - for (size_t tsIdx = 0; tsIdx < this->timeStepCount(scalarResultIndex); tsIdx++) - { - histCalc.addData(m_cellScalarResults[tranX][tsIdx]); - histCalc.addData(m_cellScalarResults[tranY][tsIdx]); - histCalc.addData(m_cellScalarResults[tranZ][tsIdx]); - } - } - } - else - { - for (size_t tsIdx = 0; tsIdx < this->timeStepCount(scalarResultIndex); tsIdx++) - { - std::vector& values = m_cellScalarResults[scalarResultIndex][tsIdx]; - - histCalc.addData(values); - } - } - - m_p10p90[scalarResultIndex].first = histCalc.calculatePercentil(0.1); - m_p10p90[scalarResultIndex].second = histCalc.calculatePercentil(0.9); - - return m_histograms[scalarResultIndex]; + return m_statisticsDataCache[scalarResultIndex]->cellScalarValuesHistogram(); } - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RigCaseCellResultsData::p10p90CellScalarValues(size_t scalarResultIndex, double& p10, double& p90) { - // First make sure they are calculated - const std::vector& histogr = cellScalarValuesHistogram( scalarResultIndex); - // Then return them - p10 = m_p10p90[scalarResultIndex].first; - p90 = m_p10p90[scalarResultIndex].second; + m_statisticsDataCache[scalarResultIndex]->p10p90CellScalarValues(p10, p90); } //-------------------------------------------------------------------------------------------------- @@ -239,73 +97,7 @@ void RigCaseCellResultsData::p10p90CellScalarValues(size_t scalarResultIndex, do //-------------------------------------------------------------------------------------------------- void RigCaseCellResultsData::meanCellScalarValues(size_t scalarResultIndex, double& meanValue) { - CVF_ASSERT(scalarResultIndex < resultCount()); - - // Extend array and cache vars - - if (scalarResultIndex >= m_meanValues.size() ) - { - m_meanValues.resize(scalarResultIndex+1, HUGE_VAL); - } - - if (m_meanValues[scalarResultIndex] != HUGE_VAL) - { - meanValue = m_meanValues[scalarResultIndex]; - return; - } - - double valueSum = 0.0; - size_t count = 0; - - if (scalarResultIndex == m_combinedTransmissibilityResultIndex) - { - size_t tranX, tranY, tranZ; - if (findTransmissibilityResults(tranX, tranY, tranZ)) - { - for (size_t tIdx = 0; tIdx < timeStepCount(tranX); tIdx++) - { - { - std::vector& values = m_cellScalarResults[tranX][tIdx]; - for (size_t cIdx = 0; cIdx < values.size(); ++cIdx) - { - valueSum += values[cIdx]; - } - count += values.size(); - } - { - std::vector& values = m_cellScalarResults[tranY][tIdx]; - for (size_t cIdx = 0; cIdx < values.size(); ++cIdx) - { - valueSum += values[cIdx]; - } - count += values.size(); - } - { - std::vector& values = m_cellScalarResults[tranZ][tIdx]; - for (size_t cIdx = 0; cIdx < values.size(); ++cIdx) - { - valueSum += values[cIdx]; - } - count += values.size(); - } - } - } - } - else - { - for (size_t tIdx = 0; tIdx < timeStepCount(scalarResultIndex); tIdx++) - { - std::vector& values = m_cellScalarResults[scalarResultIndex][tIdx]; - for (size_t cIdx = 0; cIdx < values.size(); ++cIdx) - { - valueSum += values[cIdx]; - } - count += values.size(); - } - } - - m_meanValues[scalarResultIndex] = valueSum/count; - meanValue = m_meanValues[scalarResultIndex]; + m_statisticsDataCache[scalarResultIndex]->meanCellScalarValues(meanValue); } //-------------------------------------------------------------------------------------------------- @@ -435,6 +227,29 @@ size_t RigCaseCellResultsData::addEmptyScalarResult(RimDefines::ResultCatType ty m_cellScalarResults.push_back(std::vector >()); ResultInfo resInfo(type, needsToBeStored, false, resultName, scalarResultIndex); m_resultInfos.push_back(resInfo); + + // Create statistics calculator and add cache object + if (resultName == RimDefines::combinedTransmissibilityResultName()) + { + size_t tranX = findScalarResultIndex(RimDefines::STATIC_NATIVE, "TRANX"); + size_t tranY = findScalarResultIndex(RimDefines::STATIC_NATIVE, "TRANY"); + size_t tranZ = findScalarResultIndex(RimDefines::STATIC_NATIVE, "TRANZ"); + + cvf::ref calc = new RigMultipleDatasetStatCalc(); + calc->addStatisticsCalculator(new RigNativeStatCalc(this, tranX)); + calc->addStatisticsCalculator(new RigNativeStatCalc(this, tranY)); + calc->addStatisticsCalculator(new RigNativeStatCalc(this, tranZ)); + + cvf::ref dataCache = new RigStatisticsDataCache(calc.p()); + m_statisticsDataCache.push_back(dataCache.p()); + } + else + { + cvf::ref calc = new RigNativeStatCalc(this, scalarResultIndex); + + cvf::ref dataCache = new RigStatisticsDataCache(calc.p()); + m_statisticsDataCache.push_back(dataCache.p()); + } } return scalarResultIndex; @@ -460,20 +275,9 @@ QStringList RigCaseCellResultsData::resultNames(RimDefines::ResultCatType resTyp //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RigCaseCellResultsData::recalculateMinMax(size_t scalarResultIndex) +void RigCaseCellResultsData::recalculateStatistics(size_t scalarResultIndex) { - // Make sure cached max min values are recalculated next time asked for, since - // the data could be changed. - - if (scalarResultIndex < m_maxMinValues.size()) - { - m_maxMinValues[scalarResultIndex] = std::make_pair(HUGE_VAL, -HUGE_VAL); - } - - if (scalarResultIndex < m_maxMinValuesPrTs.size()) - { - m_maxMinValuesPrTs[scalarResultIndex].clear(); - } + m_statisticsDataCache[scalarResultIndex]->clearAllStatistics(); } //-------------------------------------------------------------------------------------------------- @@ -592,15 +396,9 @@ void RigCaseCellResultsData::removeResult(const QString& resultName) void RigCaseCellResultsData::clearAllResults() { m_cellScalarResults.clear(); - m_maxMinValues.clear(); - m_histograms.clear(); - m_p10p90.clear(); - m_meanValues.clear(); - m_maxMinValuesPrTs.clear(); m_resultInfos.clear(); } - //-------------------------------------------------------------------------------------------------- /// Removes all the actual numbers put into this object, and frees up the memory. /// Does not touch the metadata in any way @@ -673,122 +471,6 @@ void RigCaseCellResultsData::setMustBeCalculated(size_t scalarResultIndex) } } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RigCaseCellResultsData::posNegClosestToZero(size_t scalarResultIndex, double& pos, double& neg) -{ - pos = HUGE_VAL; - neg = -HUGE_VAL; - - CVF_ASSERT(scalarResultIndex < resultCount()); - - // Extend array and cache vars - - if (scalarResultIndex >= m_posNegClosestToZero.size() ) - { - m_posNegClosestToZero.resize(scalarResultIndex+1, std::make_pair(HUGE_VAL, -HUGE_VAL)); - } - - if (m_posNegClosestToZero[scalarResultIndex].first != HUGE_VAL) - { - pos = m_posNegClosestToZero[scalarResultIndex].first; - neg = m_posNegClosestToZero[scalarResultIndex].second; - - return; - } - - size_t i; - for (i = 0; i < timeStepCount(scalarResultIndex); i++) - { - double tsNeg, tsPos; - posNegClosestToZero(scalarResultIndex, i, tsPos, tsNeg); - if (tsNeg > neg && tsNeg < 0) neg = tsNeg; - if (tsPos < pos && tsPos > 0) pos = tsPos; - } - - m_posNegClosestToZero[scalarResultIndex].first = pos; - m_posNegClosestToZero[scalarResultIndex].second= neg; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RigCaseCellResultsData::posNegClosestToZero(size_t scalarResultIndex, size_t timeStepIndex, double& pos, double& neg) -{ - pos = HUGE_VAL; - neg = -HUGE_VAL; - - CVF_ASSERT(scalarResultIndex < resultCount()); - - if (timeStepIndex >= m_cellScalarResults[scalarResultIndex].size()) - { - return; - } - - if (scalarResultIndex >= m_posNegClosestToZeroPrTs.size()) - { - m_posNegClosestToZeroPrTs.resize(scalarResultIndex+1); - } - - if (timeStepIndex >= m_posNegClosestToZeroPrTs[scalarResultIndex].size()) - { - m_posNegClosestToZeroPrTs[scalarResultIndex].resize(timeStepIndex+1, std::make_pair(HUGE_VAL, -HUGE_VAL)); - } - - if (m_posNegClosestToZeroPrTs[scalarResultIndex][timeStepIndex].first != HUGE_VAL) - { - pos = m_posNegClosestToZeroPrTs[scalarResultIndex][timeStepIndex].first; - neg = m_posNegClosestToZeroPrTs[scalarResultIndex][timeStepIndex].second; - - return; - } - - if (scalarResultIndex == m_combinedTransmissibilityResultIndex) - { - size_t tranX, tranY, tranZ; - if (findTransmissibilityResults(tranX, tranY, tranZ)) - { - double traPos, traNeg; - posNegClosestToZero(tranX, timeStepIndex, traPos, traNeg); - if ( 0 < traPos && traPos < pos ) pos = traPos; - if ( neg < traNeg && traNeg < 0 ) neg = traNeg; - posNegClosestToZero(tranY, timeStepIndex, traPos, traNeg); - if ( 0 < traPos && traPos < pos ) pos = traPos; - if ( neg < traNeg && traNeg < 0 ) neg = traNeg; - posNegClosestToZero(tranZ, timeStepIndex, traPos, traNeg); - if ( 0 < traPos && traPos < pos ) pos = traPos; - if ( neg < traNeg && traNeg < 0 ) neg = traNeg; - } - - return; - } - - std::vector& values = m_cellScalarResults[scalarResultIndex][timeStepIndex]; - - size_t i; - for (i = 0; i < values.size(); i++) - { - if (values[i] == HUGE_VAL) - { - continue; - } - - if (values[i] < pos && values[i] > 0) - { - pos = values[i]; - } - - if (values[i] > neg && values[i] < 0) - { - neg = values[i]; - } - } - - m_posNegClosestToZeroPrTs[scalarResultIndex][timeStepIndex].first = pos; - m_posNegClosestToZeroPrTs[scalarResultIndex][timeStepIndex].second= neg; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -800,8 +482,7 @@ void RigCaseCellResultsData::createCombinedTransmissibilityResult() size_t tranX, tranY, tranZ; if (!findTransmissibilityResults(tranX, tranY, tranZ)) return; - - m_combinedTransmissibilityResultIndex = addStaticScalarResult(RimDefines::STATIC_NATIVE, RimDefines::combinedTransmissibilityResultName(), false, 0); + addStaticScalarResult(RimDefines::STATIC_NATIVE, RimDefines::combinedTransmissibilityResultName(), false, 0); } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.h b/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.h index ec88739efa..603fb41b4f 100644 --- a/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.h +++ b/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.h @@ -18,14 +18,18 @@ #pragma once +#include "RifReaderInterface.h" + #include "RimDefines.h" + #include + #include #include -#include "RifReaderInterface.h" class RifReaderInterface; class RigMainGrid; +class RigStatisticsDataCache; //================================================================================================== /// Class containing the results for the complete number of active cells. Both main grid and LGR's @@ -38,7 +42,7 @@ public: void setMainGrid(RigMainGrid* ownerGrid); // Max and min values of the results - void recalculateMinMax(size_t scalarResultIndex); + void recalculateStatistics(size_t scalarResultIndex); void minMaxCellScalarValues(size_t scalarResultIndex, double& min, double& max); void minMaxCellScalarValues(size_t scalarResultIndex, size_t timeStepIndex, double& min, double& max); void posNegClosestToZero(size_t scalarResultIndex, double& pos, double& neg); @@ -113,16 +117,7 @@ public: private: std::vector< std::vector< std::vector > > m_cellScalarResults; ///< Scalar results on the complete reservoir for each Result index (ResultVariable) and timestep - std::vector< std::pair > m_maxMinValues; ///< Max min values for each Result index - std::vector< std::pair > m_posNegClosestToZero; - std::vector< std::vector > m_histograms; ///< Histogram for each Result Index - std::vector< std::pair > m_p10p90; ///< P10 and p90 values for each Result Index - std::vector< double > m_meanValues; ///< Mean value for each Result Index - - std::vector< std::vector< std::pair > > m_maxMinValuesPrTs; ///< Max min values for each Result index and timestep - std::vector< std::vector< std::pair > > m_posNegClosestToZeroPrTs; - - size_t m_combinedTransmissibilityResultIndex; + cvf::Collection m_statisticsDataCache; private: std::vector m_resultInfos; diff --git a/ApplicationCode/ReservoirDataModel/RigStatisticsCalculator.cpp b/ApplicationCode/ReservoirDataModel/RigStatisticsCalculator.cpp new file mode 100644 index 0000000000..48d98e035c --- /dev/null +++ b/ApplicationCode/ReservoirDataModel/RigStatisticsCalculator.cpp @@ -0,0 +1,240 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2011-2012 Statoil ASA, Ceetron AS +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "RigStatisticsCalculator.h" + +#include "RigStatisticsMath.h" +#include "RigCaseCellResultsData.h" + +#include // Needed for HUGE_VAL on Linux + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigStatisticsCalculator::meanCellScalarValue(double& meanValue) +{ + double valueSum = 0.0; + size_t sampleCount = 0; + + this->valueSumAndSampleCount(valueSum, sampleCount); + + if (sampleCount == 0) + { + meanValue = HUGE_VAL; + } + else + { + meanValue = valueSum / sampleCount; + } +} + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigNativeStatCalc::RigNativeStatCalc(RigCaseCellResultsData* cellResultsData, size_t scalarResultIndex) + : m_resultsData(cellResultsData), + m_scalarResultIndex(scalarResultIndex) +{ + +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigNativeStatCalc::minMaxCellScalarValues(size_t timeStepIndex, double& min, double& max) +{ + std::vector& values = m_resultsData->cellScalarResults(m_scalarResultIndex, timeStepIndex); + + size_t i; + for (i = 0; i < values.size(); i++) + { + if (values[i] == HUGE_VAL) + { + continue; + } + + if (values[i] < min) + { + min = values[i]; + } + + if (values[i] > max) + { + max = values[i]; + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigNativeStatCalc::posNegClosestToZero(size_t timeStepIndex, double& pos, double& neg) +{ + std::vector& values = m_resultsData->cellScalarResults(m_scalarResultIndex, timeStepIndex); + + size_t i; + for (i = 0; i < values.size(); i++) + { + if (values[i] == HUGE_VAL) + { + continue; + } + + if (values[i] < pos && values[i] > 0) + { + pos = values[i]; + } + + if (values[i] > neg && values[i] < 0) + { + neg = values[i]; + } + } +} + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigNativeStatCalc::addDataToHistogramCalculator(RigHistogramCalculator& histogramCalculator) +{ + for (size_t tIdx = 0; tIdx < m_resultsData->timeStepCount(m_scalarResultIndex); tIdx++) + { + std::vector& values = m_resultsData->cellScalarResults(m_scalarResultIndex, tIdx); + + histogramCalculator.addData(values); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigNativeStatCalc::valueSumAndSampleCount(double& valueSum, size_t& sampleCount) +{ + for (size_t tIdx = 0; tIdx < m_resultsData->timeStepCount(m_scalarResultIndex); tIdx++) + { + std::vector& values = m_resultsData->cellScalarResults(m_scalarResultIndex, tIdx); + for (size_t cIdx = 0; cIdx < values.size(); ++cIdx) + { + valueSum += values[cIdx]; + } + sampleCount += values.size(); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigNativeStatCalc::timeStepCount() +{ + return m_resultsData->timeStepCount(m_scalarResultIndex); +} + + + + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigMultipleDatasetStatCalc::RigMultipleDatasetStatCalc() +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigMultipleDatasetStatCalc::addStatisticsCalculator(RigStatisticsCalculator* statisticsCalculator) +{ + if (statisticsCalculator) + { + m_nativeStatisticsCalculators.push_back(statisticsCalculator); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigMultipleDatasetStatCalc::minMaxCellScalarValues(size_t timeStepIndex, double& min, double& max) +{ + for (size_t i = 0; i < m_nativeStatisticsCalculators.size(); i++) + { + if (m_nativeStatisticsCalculators.at(i)) + { + m_nativeStatisticsCalculators.at(i)->minMaxCellScalarValues(timeStepIndex, min, max); + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigMultipleDatasetStatCalc::posNegClosestToZero(size_t timeStepIndex, double& pos, double& neg) +{ + for (size_t i = 0; i < m_nativeStatisticsCalculators.size(); i++) + { + if (m_nativeStatisticsCalculators.at(i)) + { + m_nativeStatisticsCalculators.at(i)->posNegClosestToZero(timeStepIndex, pos, neg); + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigMultipleDatasetStatCalc::valueSumAndSampleCount(double& valueSum, size_t& sampleCount) +{ + for (size_t i = 0; i < m_nativeStatisticsCalculators.size(); i++) + { + if (m_nativeStatisticsCalculators.at(i)) + { + m_nativeStatisticsCalculators.at(i)->valueSumAndSampleCount(valueSum, sampleCount); + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigMultipleDatasetStatCalc::addDataToHistogramCalculator(RigHistogramCalculator& histogramCalculator) +{ + for (size_t i = 0; i < m_nativeStatisticsCalculators.size(); i++) + { + if (m_nativeStatisticsCalculators.at(i)) + { + m_nativeStatisticsCalculators.at(i)->addDataToHistogramCalculator(histogramCalculator); + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigMultipleDatasetStatCalc::timeStepCount() +{ + if (m_nativeStatisticsCalculators.size() > 0) + { + return m_nativeStatisticsCalculators[0]->timeStepCount(); + } + + return 0; +} + diff --git a/ApplicationCode/ReservoirDataModel/RigStatisticsCalculator.h b/ApplicationCode/ReservoirDataModel/RigStatisticsCalculator.h new file mode 100644 index 0000000000..fbd4a43270 --- /dev/null +++ b/ApplicationCode/ReservoirDataModel/RigStatisticsCalculator.h @@ -0,0 +1,87 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) Statoil ASA, Ceetron Solutions AS +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "cvfBase.h" +#include "cvfObject.h" +#include "cvfCollection.h" + +#include + +class RigHistogramCalculator; +class RigCaseCellResultsData; + +//================================================================================================== +/// +//================================================================================================== +class RigStatisticsCalculator : public cvf::Object +{ +public: + virtual void minMaxCellScalarValues(size_t timeStepIndex, double& min, double& max) = 0; + virtual void posNegClosestToZero(size_t timeStepIndex, double& pos, double& neg) = 0; + + void meanCellScalarValue(double& meanValue); + virtual void valueSumAndSampleCount(double& valueSum, size_t& sampleCount) = 0; + virtual void addDataToHistogramCalculator(RigHistogramCalculator& histogramCalculator) = 0; + + virtual size_t timeStepCount() = 0; +}; + + +//================================================================================================== +/// +//================================================================================================== +class RigNativeStatCalc : public RigStatisticsCalculator +{ +public: + RigNativeStatCalc(RigCaseCellResultsData* cellResultsData, size_t scalarResultIndex); + + virtual void minMaxCellScalarValues(size_t timeStepIndex, double& min, double& max); + virtual void posNegClosestToZero(size_t timeStepIndex, double& pos, double& neg); + virtual void valueSumAndSampleCount(double& valueSum, size_t& sampleCount); + + virtual void addDataToHistogramCalculator(RigHistogramCalculator& histogramCalculator); + virtual size_t timeStepCount(); + +private: + cvf::ref m_resultsData; + size_t m_scalarResultIndex; +}; + + +//================================================================================================== +/// +//================================================================================================== +class RigMultipleDatasetStatCalc : public RigStatisticsCalculator +{ +public: + RigMultipleDatasetStatCalc(); + void addStatisticsCalculator(RigStatisticsCalculator* statisticsCalculator); + + virtual void minMaxCellScalarValues(size_t timeStepIndex, double& min, double& max); + virtual void posNegClosestToZero(size_t timeStepIndex, double& pos, double& neg); + + virtual void valueSumAndSampleCount(double& valueSum, size_t& sampleCount); + virtual void addDataToHistogramCalculator(RigHistogramCalculator& histogramCalculator); + + virtual size_t timeStepCount(); + +private: + cvf::Collection m_nativeStatisticsCalculators; +}; diff --git a/ApplicationCode/ReservoirDataModel/RigStatisticsDataCache.cpp b/ApplicationCode/ReservoirDataModel/RigStatisticsDataCache.cpp new file mode 100644 index 0000000000..f0dc353920 --- /dev/null +++ b/ApplicationCode/ReservoirDataModel/RigStatisticsDataCache.cpp @@ -0,0 +1,205 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) Statoil ASA, Ceetron Solutions AS +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "RigStatisticsDataCache.h" + +#include "RigStatisticsCalculator.h" +#include "RigStatisticsMath.h" + +#include // Needed for HUGE_VAL on Linux + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigStatisticsDataCache::RigStatisticsDataCache(RigStatisticsCalculator* statisticsCalculator) + : m_statisticsCalculator(statisticsCalculator) +{ + clearAllStatistics(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigStatisticsDataCache::clearAllStatistics() +{ + m_minValue = HUGE_VAL; + m_maxValue = -HUGE_VAL; + m_posClosestToZero = -HUGE_VAL; + m_negClosestToZero = HUGE_VAL; + m_p10 = HUGE_VAL; + m_p90 = HUGE_VAL; + m_meanValue = HUGE_VAL; + + m_histogram.clear(); + m_maxMinValuesPrTs.clear(); + m_posNegClosestToZeroPrTs.clear(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigStatisticsDataCache::minMaxCellScalarValues(double& min, double& max) +{ + if (m_minValue == HUGE_VAL) + { + min = HUGE_VAL; + max = -HUGE_VAL; + + size_t i; + for (i = 0; i < m_statisticsCalculator->timeStepCount(); i++) + { + double tsmin, tsmax; + this->minMaxCellScalarValues(i, tsmin, tsmax); + if (tsmin < min) min = tsmin; + if (tsmax > max) max = tsmax; + } + + m_minValue = min; + m_maxValue = max; + } + + min = m_minValue; + max = m_maxValue; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigStatisticsDataCache::minMaxCellScalarValues(size_t timeStepIndex, double& min, double& max) +{ + if (timeStepIndex >= m_maxMinValuesPrTs.size()) + { + m_maxMinValuesPrTs.resize(timeStepIndex + 1, std::make_pair(HUGE_VAL, -HUGE_VAL)); + } + + if (m_maxMinValuesPrTs[timeStepIndex].first == HUGE_VAL) + { + min = HUGE_VAL; + max = -HUGE_VAL; + + m_statisticsCalculator->minMaxCellScalarValues(timeStepIndex, min, max); + + m_maxMinValuesPrTs[timeStepIndex].first = min; + m_maxMinValuesPrTs[timeStepIndex].second = max; + } + + min = m_maxMinValuesPrTs[timeStepIndex].first; + max = m_maxMinValuesPrTs[timeStepIndex].second; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigStatisticsDataCache::posNegClosestToZero(double& pos, double& neg) +{ + if (m_posClosestToZero == HUGE_VAL) + { + pos = HUGE_VAL; + neg = -HUGE_VAL; + + size_t i; + for (i = 0; i < m_statisticsCalculator->timeStepCount(); i++) + { + double tsNeg, tsPos; + this->posNegClosestToZero(i, tsPos, tsNeg); + if (tsNeg > neg && tsNeg < 0) neg = tsNeg; + if (tsPos < pos && tsPos > 0) pos = tsPos; + } + + m_posClosestToZero = pos; + m_negClosestToZero = neg; + } + + pos = m_posClosestToZero; + neg = m_negClosestToZero; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigStatisticsDataCache::posNegClosestToZero(size_t timeStepIndex, double& pos, double& neg) +{ + if (timeStepIndex >= m_posNegClosestToZeroPrTs.size()) + { + m_posNegClosestToZeroPrTs.resize(timeStepIndex + 1, std::make_pair(HUGE_VAL, -HUGE_VAL)); + } + + if (m_posNegClosestToZeroPrTs[timeStepIndex].first == HUGE_VAL) + { + pos = HUGE_VAL; + neg = -HUGE_VAL; + + m_statisticsCalculator->posNegClosestToZero(timeStepIndex, pos, neg); + + m_posNegClosestToZeroPrTs[timeStepIndex].first = pos; + m_posNegClosestToZeroPrTs[timeStepIndex].second = neg; + } + + pos = m_posNegClosestToZeroPrTs[timeStepIndex].first; + neg = m_posNegClosestToZeroPrTs[timeStepIndex].second; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::vector& RigStatisticsDataCache::cellScalarValuesHistogram() +{ + if (m_histogram.size() == 0) + { + double min; + double max; + size_t nBins = 100; + this->minMaxCellScalarValues(min, max); + + RigHistogramCalculator histCalc(min, max, nBins, &m_histogram); + + m_statisticsCalculator->addDataToHistogramCalculator(histCalc); + + m_p10 = histCalc.calculatePercentil(0.1); + m_p90 = histCalc.calculatePercentil(0.9); + } + + return m_histogram; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigStatisticsDataCache::p10p90CellScalarValues(double& p10, double& p90) +{ + // First make sure they are calculated + const std::vector& histogr = this->cellScalarValuesHistogram(); + + p10 = m_p10; + p90 = m_p90; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigStatisticsDataCache::meanCellScalarValues(double& meanValue) +{ + if (m_meanValue == HUGE_VAL) + { + m_statisticsCalculator->meanCellScalarValue(m_meanValue); + } + + meanValue = m_meanValue; +} + diff --git a/ApplicationCode/ReservoirDataModel/RigStatisticsDataCache.h b/ApplicationCode/ReservoirDataModel/RigStatisticsDataCache.h new file mode 100644 index 0000000000..b588a8f511 --- /dev/null +++ b/ApplicationCode/ReservoirDataModel/RigStatisticsDataCache.h @@ -0,0 +1,65 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) Statoil ASA, Ceetron Solutions AS +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "RigStatisticsCalculator.h" + +#include "cvfBase.h" +#include "cvfObject.h" + +#include + + +//================================================================================================== +/// +//================================================================================================== +class RigStatisticsDataCache : public cvf::Object +{ +public: + RigStatisticsDataCache(RigStatisticsCalculator* statisticsCalculator); + + void clearAllStatistics(); + + void minMaxCellScalarValues(double& min, double& max); + void minMaxCellScalarValues(size_t timeStepIndex, double& min, double& max); + void posNegClosestToZero(double& pos, double& neg); + void posNegClosestToZero(size_t timeStepIndex, double& pos, double& neg); + + void p10p90CellScalarValues(double& p10, double& p90); + void meanCellScalarValues(double& meanValue); + const std::vector& cellScalarValuesHistogram(); + +private: + double m_minValue; + double m_maxValue; + double m_posClosestToZero; + double m_negClosestToZero; + + double m_p10; + double m_p90; + double m_meanValue; + + std::vector m_histogram; + + std::vector > m_maxMinValuesPrTs; ///< Max min values for each time step + std::vector > m_posNegClosestToZeroPrTs; ///< PosNeg values for each time step + + cvf::ref m_statisticsCalculator; +}; + diff --git a/ApplicationCode/SocketInterface/RiaPropertyDataCommands.cpp b/ApplicationCode/SocketInterface/RiaPropertyDataCommands.cpp index a8257f834b..7a7e83bf8e 100644 --- a/ApplicationCode/SocketInterface/RiaPropertyDataCommands.cpp +++ b/ApplicationCode/SocketInterface/RiaPropertyDataCommands.cpp @@ -640,7 +640,7 @@ public: m_currentReservoir->reservoirData() && m_currentReservoir->reservoirData()->results(m_porosityModelEnum) ) { - m_currentReservoir->reservoirData()->results(m_porosityModelEnum)->recalculateMinMax(m_currentScalarIndex); + m_currentReservoir->reservoirData()->results(m_porosityModelEnum)->recalculateStatistics(m_currentScalarIndex); } for (size_t i = 0; i < m_currentReservoir->reservoirViews.size(); ++i) @@ -974,7 +974,7 @@ public: m_currentReservoir->reservoirData() && m_currentReservoir->reservoirData()->results(m_porosityModelEnum) ) { - m_currentReservoir->reservoirData()->results(m_porosityModelEnum)->recalculateMinMax(m_currentScalarIndex); + m_currentReservoir->reservoirData()->results(m_porosityModelEnum)->recalculateStatistics(m_currentScalarIndex); } for (size_t i = 0; i < m_currentReservoir->reservoirViews.size(); ++i)