OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2026-09-03 20:53:13 -05:00
Code Issues Packages Projects Releases Wiki Activity

#14627 Grid statistics: Fix NaN deviation for near-constant results

Browse Source 
The single-pass sum-of-squares formula subtracts two large and nearly equal numbers. Where a result such as PORO is identical or near-identical across the realizations, the true variance is below the rounding noise floor, the radicand turns negative and sqrt() returns NaN.

Compute the deviation in two passes instead, centering on the exact mean. This also removes the loss of significance the old formula had for values with a large common offset.

Min, max, range, mean and deviation all report HUGE_VAL when the input holds no valid values, but the sum reported the zero it was accumulated from. A caller testing the result with isValidNumber would accept that zero as real data.
This commit is contained in:
Magne Sjaastad
2026-08-28 07:01:57 +02:00
committed by GitHub
parent ca8c1b8d04
commit 5e17f9986d
2 changed files with 105 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
ApplicationLibCode/ResultStatisticsCache/RigStatisticsMath.cpp
+35 -24
View File
@@ -26,7 +26,6 @@
#include <cassert>
#include <cmath>
#include <expected>
#include <numeric>
namespace
{
@@ -58,15 +57,12 @@ bool areValidPercentiles( const std::vector<double>& percentiles )
/// mean = sum(x) / n
///
/// Standard deviation (population):
/// stdev = sqrt((n * sum(x^2) - (sum(x))^2)) / n
///
/// Which is equivalent to: sqrt(sum((x - mean)^2) / n)
/// stdev = sqrt(sum((x - mean)^2) / n)
///
/// range = max - min
///
/// References:
/// Standard deviation: https://en.wikipedia.org/wiki/Standard_deviation
/// Rapid calculation method: https://en.wikipedia.org/wiki/Standard_deviation#Rapid_calculation_methods
//--------------------------------------------------------------------------------------------------
void RigStatisticsMath::calculateBasicStatistics( const std::vector<double>& values,
@@ -82,8 +78,7 @@ void RigStatisticsMath::calculateBasicStatistics( const std::vector<double>& val
double m_mean( HUGE_VAL );
double m_dev( HUGE_VAL );
double m_sum = 0.0;
double sumSquared = 0.0;
double m_sum = 0.0;
size_t validValueCount = 0;
@@ -98,26 +93,40 @@ void RigStatisticsMath::calculateBasicStatistics( const std::vector<double>& val
if ( val > m_max ) m_max = val;
m_sum += val;
sumSquared += ( val * val );
}
if ( validValueCount > 0 )
{
m_mean = m_sum / validValueCount;
// http://en.wikipedia.org/wiki/Standard_deviation#Rapid_calculation_methods
// Running standard deviation
// Two-pass computation of the deviation. The single-pass sum-of-squares formula
// n*sum(x*x) - sum(x)*sum(x) subtracts two large and nearly equal numbers, and for values with a
// spread below the rounding noise floor the result can become negative, making sqrt() return NaN.
//
// Welford's online algorithm avoids the same cancellation in a single pass, but is only required when the
// values cannot be revisited. All values are available here, and centering on the exact mean is more
// accurate than Welford's running estimate.
//
// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm
double s0 = static_cast<double>( validValueCount );
double s1 = m_sum;
double s2 = sumSquared;
double sumSquaredDeviations = 0.0;
m_dev = sqrt( ( s0 * s2 ) - ( s1 * s1 ) ) / s0;
for ( size_t i = 0; i < values.size(); i++ )
{
double val = values[i];
if ( RigStatisticsTools::isInvalidNumber<double>( val ) ) continue;
double deviation = val - m_mean;
sumSquaredDeviations += deviation * deviation;
}
m_dev = sqrt( sumSquaredDeviations / validValueCount );
}
if ( min ) *min = m_min;
if ( max ) *max = m_max;
if ( sum ) *sum = m_sum;
// Report the sum as undefined when there is nothing to sum, so it is not mistaken for a computed zero
if ( sum ) *sum = ( validValueCount > 0 ) ? m_sum : HUGE_VAL;
if ( range ) *range = m_max - m_min;
if ( mean ) *mean = m_mean;
@@ -433,19 +442,21 @@ std::expected<std::vector<double>, std::string>
//--------------------------------------------------------------------------------------------------
double RigStatisticsMath::calculateMean( const std::vector<double>& values )
{
std::vector<double> validValues = values;
validValues.erase( std::remove_if( validValues.begin(),
validValues.end(),
[]( double x ) { return !RigStatisticsTools::isValidNumber( x ); } ),
validValues.end() );
double valueSum = 0.0;
size_t validValueCount = 0;
if ( !validValues.empty() )
for ( size_t i = 0; i < values.size(); i++ )
{
double valueSum = std::accumulate( validValues.begin(), validValues.end(), 0.0 );
return valueSum / validValues.size();
double val = values[i];
if ( RigStatisticsTools::isInvalidNumber<double>( val ) ) continue;
valueSum += val;
validValueCount++;
}
return HUGE_VAL;
if ( validValueCount == 0 ) return HUGE_VAL;
return valueSum / validValueCount;
}
//--------------------------------------------------------------------------------------------------
ApplicationLibCode/UnitTests/RigStatisticsMath-Test.cpp
+70
View File
@@ -26,6 +26,7 @@
#include "QElapsedTimer"
#include <cmath>
#include <limits>
#include <numeric>
//--------------------------------------------------------------------------------------------------
@@ -63,6 +64,75 @@ TEST( RigStatisticsMath, BasicTest )
EXPECT_DOUBLE_EQ( 66104.391542887200, stdev );
}
//--------------------------------------------------------------------------------------------------
/// All outputs must report as undefined when the input holds no valid values. The sum in particular
/// must not report the zero it was accumulated from, as that reads as a computed value.
//--------------------------------------------------------------------------------------------------
TEST( RigStatisticsMath, NoValidValues )
{
std::vector<double> values = { HUGE_VAL, -HUGE_VAL, std::numeric_limits<double>::quiet_NaN() };
double min, max, sum, range, mean, stdev;
RigStatisticsMath::calculateBasicStatistics( values, &min, &max, &sum, &range, &mean, &stdev );
EXPECT_FALSE( RigStatisticsTools::isValidNumber( min ) );
EXPECT_FALSE( RigStatisticsTools::isValidNumber( max ) );
EXPECT_FALSE( RigStatisticsTools::isValidNumber( sum ) );
EXPECT_FALSE( RigStatisticsTools::isValidNumber( range ) );
EXPECT_FALSE( RigStatisticsTools::isValidNumber( mean ) );
EXPECT_FALSE( RigStatisticsTools::isValidNumber( stdev ) );
}
//--------------------------------------------------------------------------------------------------
/// A porosity value that is identical in all realizations has zero deviation. The sum-of-squares
/// formula computes n*sum(x*x) - sum(x)*sum(x), a difference of two large and nearly equal numbers,
/// so rounding can make the radicand negative and sqrt() return NaN.
//--------------------------------------------------------------------------------------------------
TEST( RigStatisticsMath, DeviationOfIdenticalValues )
{
const std::vector<double> porosities = { 0.2345, 0.1, 0.31415926, 0.27, 0.185, 0.3, 0.15 };
const std::vector<size_t> sourceCounts = { 3, 5, 10, 25, 50, 100 };
for ( double porosity : porosities )
{
for ( size_t sourceCount : sourceCounts )
{
SCOPED_TRACE( "porosity " + std::to_string( porosity ) + " in " + std::to_string( sourceCount ) + " cases" );
std::vector<double> values( sourceCount, porosity );
double min, max, sum, range, mean, stdev;
RigStatisticsMath::calculateBasicStatistics( values, &min, &max, &sum, &range, &mean, &stdev );
EXPECT_FALSE( std::isnan( stdev ) );
EXPECT_NEAR( 0.0, stdev, 1e-12 );
}
}
}
//--------------------------------------------------------------------------------------------------
/// Same cancellation, but for values that differ far below the rounding noise floor.
//--------------------------------------------------------------------------------------------------
TEST( RigStatisticsMath, DeviationOfNearlyIdenticalValues )
{
const double porosity = 0.2345;
const double relativeNoise = 1e-9;
std::vector<double> values;
for ( size_t i = 0; i < 100; i++ )
{
// Alternating tiny perturbation around the base value
double sign = ( i % 2 == 0 ) ? 1.0 : -1.0;
values.push_back( porosity + sign * porosity * relativeNoise );
}
double min, max, sum, range, mean, stdev;
RigStatisticsMath::calculateBasicStatistics( values, &min, &max, &sum, &range, &mean, &stdev );
EXPECT_FALSE( std::isnan( stdev ) );
EXPECT_NEAR( porosity * relativeNoise, stdev, 1e-12 );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------