OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-10 23:46:00 -06:00
Code Issues Packages Projects Releases Wiki Activity

#6035 Add method for interpolating missing values (with extrapolation)

Browse Source
This commit is contained in:
Kristian Bendiksen 2020-06-08 19:48:16 +02:00
parent 56753f688d
commit 3085fe22e8
3 changed files with 195 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

125
ApplicationCode/Application/Tools/RiaInterpolationTools.cpp
View File

@ -64,3 +64,128 @@ double RiaInterpolationTools::linear( const std::vector<double>& x, const std::v
return lowerY + ( ( value - lowerX ) / deltaX ) * deltaY;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RiaInterpolationTools::extrapolate( const std::vector<double>& x, const std::vector<double>& y, double value )
{
return y[0] + ( value - x[0] ) / ( x[1] - x[0] ) * ( y[1] - y[0] );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RiaInterpolationTools::findNextDataPoint( const std::vector<double>& values, int index )
{
for ( size_t i = index; i < values.size(); i++ )
{
if ( values[i] != std::numeric_limits<double>::infinity() ) return static_cast<int>( i );
}
return -1;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RiaInterpolationTools::findPreviousDataPoint( const std::vector<double>& values, int index )
{
assert( index >= 0 );
for ( int i = index; i >= 0; i-- )
{
if ( values[i] != std::numeric_limits<double>::infinity() ) return static_cast<int>( i );
}
return -1;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RiaInterpolationTools::extrapolateRange( int start,
int end,
int firstPoint,
int lastPoint,
const std::vector<double>& x,
std::vector<double>& y )
{
std::vector<double> xs = {x[firstPoint], x[lastPoint]};
std::vector<double> ys = {y[firstPoint], y[lastPoint]};
for ( int index = start; index < end; index++ )
{
y[index] = extrapolate( xs, ys, x[index] );
}
return end;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RiaInterpolationTools::interpolateRange( int start,
int end,
int firstPoint,
int lastPoint,
const std::vector<double>& x,
std::vector<double>& y )
{
assert( start <= end );
std::vector<double> xs = {x[firstPoint], x[lastPoint]};
std::vector<double> ys = {y[firstPoint], y[lastPoint]};
for ( int index = start; index < end; index++ )
{
y[index] = RiaInterpolationTools::linear( xs, ys, x[index] );
}
return end;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RiaInterpolationTools::interpolateMissingValues( const std::vector<double>& x, std::vector<double>& y )
{
assert( x.size() == y.size() );
int index = 0;
// Previous index which is not inf
int prevSetIndex = -1;
while ( index < static_cast<int>( y.size() ) )
{
// Missing values are inf in the input data
if ( y[index] == std::numeric_limits<double>::infinity() )
{
// Find the next index with a value
int nextSetIndex = findNextDataPoint( y, index + 1 );
if ( prevSetIndex == -1 )
{
// The first value is inf: need to find next two valid points and extrapolate
int nextSetIndex2 = findNextDataPoint( y, nextSetIndex + 1 );
index = extrapolateRange( index, nextSetIndex, nextSetIndex, nextSetIndex2, x, y );
}
else if ( nextSetIndex == -1 )
{
// The last value is inf: extrapolate from two last data points
int prevSetIndex2 = findPreviousDataPoint( y, prevSetIndex - 1 );
index = extrapolateRange( index, y.size(), prevSetIndex2, prevSetIndex, x, y );
}
else if ( nextSetIndex != static_cast<int>( y.size() ) )
{
// The missing values somewhere between non-inf data: interpolate all the values
index = interpolateRange( index, nextSetIndex, prevSetIndex, nextSetIndex, x, y );
}
}
else
{
// Nothing to do for the values which are not missing
prevSetIndex = index;
++index;
}
}
}

20
ApplicationCode/Application/Tools/RiaInterpolationTools.h
View File

@ -27,4 +27,24 @@ class RiaInterpolationTools
{
public:
static double linear( const std::vector<double>& x, const std::vector<double>& y, double value );
// Interpolate/extrapolate away inf values in y vector.
static void interpolateMissingValues( const std::vector<double>& x, std::vector<double>& y );
private:
static int interpolateRange( int start,
int end,
int firstPoint,
int lastPoint,
const std::vector<double>& x,
std::vector<double>& y );
static int extrapolateRange( int start,
int end,
int firstPoint,
int lastPoint,
const std::vector<double>& x,
std::vector<double>& y );
static int findNextDataPoint( const std::vector<double>& values, int index );
static int findPreviousDataPoint( const std::vector<double>& values, int index );
static double extrapolate( const std::vector<double>& x, const std::vector<double>& y, double value );
};

50
ApplicationCode/UnitTests/RiaInterpolationTools-Test.cpp
View File

@ -77,3 +77,53 @@ TEST( RiaInterpolationToolsTest, ValidIntervalValueTooHigh )
double res = RiaInterpolationTools::linear( x, y, 100.0 );
EXPECT_DOUBLE_EQ( std::numeric_limits<double>::infinity(), res );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST( RiaInterpolationToolsTest, InterpolateMisssingValuesStraightLine )
{
double inf = std::numeric_limits<double>::infinity();
std::vector<double> x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
std::vector<double> y = {0.0, 1.0, inf, inf, inf, 5.0};
RiaInterpolationTools::interpolateMissingValues( x, y );
EXPECT_DOUBLE_EQ( y[2], 2.0 );
EXPECT_DOUBLE_EQ( y[3], 3.0 );
EXPECT_DOUBLE_EQ( y[4], 4.0 );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST( RiaInterpolationToolsTest, InterpolateMissingValuesStraightLineExtrapolateStart )
{
double inf = std::numeric_limits<double>::infinity();
std::vector<double> x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
std::vector<double> y = {inf, inf, 2.0, inf, 4.0, 5.0};
RiaInterpolationTools::interpolateMissingValues( x, y );
EXPECT_DOUBLE_EQ( y[0], 0.0 );
EXPECT_DOUBLE_EQ( y[1], 1.0 );
EXPECT_DOUBLE_EQ( y[2], 2.0 );
EXPECT_DOUBLE_EQ( y[3], 3.0 );
EXPECT_DOUBLE_EQ( y[4], 4.0 );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST( RiaInterpolationToolsTest, InterpolateMissingValuesStraightLineExtrapolateEnd )
{
double inf = std::numeric_limits<double>::infinity();
std::vector<double> x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
std::vector<double> y = {0.0, inf, 2.0, inf, 4.0, inf};
RiaInterpolationTools::interpolateMissingValues( x, y );
EXPECT_DOUBLE_EQ( y[0], 0.0 );
EXPECT_DOUBLE_EQ( y[1], 1.0 );
EXPECT_DOUBLE_EQ( y[2], 2.0 );
EXPECT_DOUBLE_EQ( y[3], 3.0 );
EXPECT_DOUBLE_EQ( y[4], 4.0 );
EXPECT_DOUBLE_EQ( y[5], 5.0 );
}