#11954 Improve the curve merger algorithm

* Add flag to indicate identical X values when curves are merged
* Detect if x values are monotonically increasing
This commit is contained in:
Magne Sjaastad 2024-12-03 16:53:11 +01:00 committed by GitHub
parent d90dd91e29
commit 1ce2f9651c
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
3 changed files with 103 additions and 14 deletions

View File

@ -68,6 +68,7 @@ public:
private: private:
void computeUnionOfXValues( bool includeValuesFromPartialCurves ); void computeUnionOfXValues( bool includeValuesFromPartialCurves );
static bool isMonotonicallyIncreasing( const std::vector<XValueType>& curveXValues );
private: private:
std::vector<std::pair<std::vector<XValueType>, std::vector<double>>> m_originalValues; std::vector<std::pair<std::vector<XValueType>, std::vector<double>>> m_originalValues;
@ -76,6 +77,9 @@ private:
std::vector<XValueType> m_allXValues; std::vector<XValueType> m_allXValues;
std::vector<std::vector<double>> m_interpolatedValuesForAllCurves; std::vector<std::vector<double>> m_interpolatedValuesForAllCurves;
bool m_isXValuesSharedBetweenCurves;
bool m_isXValuesMonotonicallyIncreasing;
}; };
using RiaTimeHistoryCurveMerger = RiaCurveMerger<time_t>; using RiaTimeHistoryCurveMerger = RiaCurveMerger<time_t>;

View File

@ -55,6 +55,8 @@ bool XValueComparator<XValueType>::equals( const XValueType& lhs, const XValueTy
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
template <typename XValueType> template <typename XValueType>
RiaCurveMerger<XValueType>::RiaCurveMerger() RiaCurveMerger<XValueType>::RiaCurveMerger()
: m_isXValuesSharedBetweenCurves( false )
, m_isXValuesMonotonicallyIncreasing( true )
{ {
} }
@ -68,6 +70,24 @@ void RiaCurveMerger<XValueType>::addCurveData( const std::vector<XValueType>& xV
if ( !xValues.empty() ) if ( !xValues.empty() )
{ {
if ( m_originalValues.empty() )
{
m_isXValuesSharedBetweenCurves = true;
m_isXValuesMonotonicallyIncreasing = isMonotonicallyIncreasing( xValues );
}
else
{
if ( m_isXValuesSharedBetweenCurves )
{
const auto& firstXValues = m_originalValues.front().first;
m_isXValuesSharedBetweenCurves = std::equal( firstXValues.begin(), firstXValues.end(), xValues.begin() );
}
if ( m_isXValuesMonotonicallyIncreasing )
{
m_isXValuesMonotonicallyIncreasing = isMonotonicallyIncreasing( xValues );
}
}
m_originalValues.push_back( std::make_pair( xValues, yValues ) ); m_originalValues.push_back( std::make_pair( xValues, yValues ) );
} }
} }
@ -131,14 +151,14 @@ void RiaCurveMerger<XValueType>::computeInterpolatedValues( bool includeValuesFr
m_allXValues.clear(); m_allXValues.clear();
m_interpolatedValuesForAllCurves.clear(); m_interpolatedValuesForAllCurves.clear();
computeUnionOfXValues( includeValuesFromPartialCurves );
const size_t curveCount = m_originalValues.size(); const size_t curveCount = m_originalValues.size();
if ( curveCount == 0 ) if ( curveCount == 0 )
{ {
return; return;
} }
computeUnionOfXValues( includeValuesFromPartialCurves );
const size_t dataValueCount = m_allXValues.size(); const size_t dataValueCount = m_allXValues.size();
if ( dataValueCount == 0 ) if ( dataValueCount == 0 )
{ {
@ -157,8 +177,18 @@ void RiaCurveMerger<XValueType>::computeInterpolatedValues( bool includeValuesFr
#pragma omp parallel for #pragma omp parallel for
for ( int valueIndex = 0; valueIndex < static_cast<int>( dataValueCount ); valueIndex++ ) for ( int valueIndex = 0; valueIndex < static_cast<int>( dataValueCount ); valueIndex++ )
{ {
double interpolValue = double interpolValue = 0.0;
if ( m_isXValuesSharedBetweenCurves )
{
interpolValue = m_originalValues[curveIdx].second[valueIndex];
}
else
{
interpolValue =
interpolatedYValue( m_allXValues[valueIndex], m_originalValues[curveIdx].first, m_originalValues[curveIdx].second ); interpolatedYValue( m_allXValues[valueIndex], m_originalValues[curveIdx].first, m_originalValues[curveIdx].second );
}
if ( !RiaCurveDataTools::isValidValue( interpolValue, false ) ) if ( !RiaCurveDataTools::isValidValue( interpolValue, false ) )
{ {
#pragma omp critical #pragma omp critical
@ -180,6 +210,13 @@ void RiaCurveMerger<XValueType>::computeUnionOfXValues( bool includeValuesForPar
{ {
m_allXValues.clear(); m_allXValues.clear();
if ( m_isXValuesSharedBetweenCurves && !m_originalValues.empty() )
{
// If all curves have the same X values, use the X values from the first curve.
m_allXValues = m_originalValues.front().first;
return;
}
std::set<XValueType, XComparator> unionOfXValues; std::set<XValueType, XComparator> unionOfXValues;
std::vector<std::pair<XValueType, XValueType>> originalXBounds; std::vector<std::pair<XValueType, XValueType>> originalXBounds;
@ -240,16 +277,26 @@ double RiaCurveMerger<XValueType>::interpolatedYValue( const XValueType&
if ( xValues.empty() ) return HUGE_VAL; if ( xValues.empty() ) return HUGE_VAL;
if ( yValues.size() != xValues.size() ) return HUGE_VAL; if ( yValues.size() != xValues.size() ) return HUGE_VAL;
const bool removeInterpolatedValues = false; size_t startIndex = 0;
// Use lower_bound to find the first element that is not less than the interpolation value using a threshold that is larger than the if ( isMonotonicallyIncreasing( xValues ) )
// threshold used in XComparator::equals {
// Use lower_bound to find the first element that is not less than the interpolation value using a threshold that is larger than
// the threshold used in XComparator::equals
//
// Using this method will improve the performance significantly for large datasets, as std::lower_bound is much faster than the
// search in the loop below. One relevant use case is computation of delta summary values for large datasets. Here the time
// steps are specified in increasing order
//
XValueType threshold = 1.0e-6 * xValues.back(); XValueType threshold = 1.0e-6 * xValues.back();
auto it = std::lower_bound( xValues.begin(), xValues.end(), interpolationXValue - threshold ); auto it = std::lower_bound( xValues.begin(), xValues.end(), interpolationXValue - threshold );
if ( it == xValues.end() ) return HUGE_VAL; if ( it == xValues.end() ) return HUGE_VAL;
size_t startIndex = it - xValues.begin(); startIndex = it - xValues.begin();
if ( startIndex > 0 ) startIndex--; if ( startIndex > 0 ) startIndex--;
}
const bool removeInterpolatedValues = false;
for ( size_t firstI = startIndex; firstI < xValues.size(); firstI++ ) for ( size_t firstI = startIndex; firstI < xValues.size(); firstI++ )
{ {
@ -307,3 +354,17 @@ double RiaCurveMerger<XValueType>::interpolatedYValue( const XValueType&
return HUGE_VAL; return HUGE_VAL;
} }
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
template <typename XValueType>
bool RiaCurveMerger<XValueType>::isMonotonicallyIncreasing( const std::vector<XValueType>& container )
{
return std::adjacent_find( container.begin(),
container.end(),
[]( const auto& a, const auto& b )
{
return b < a; // Returns true if not monotonically increasing
} ) == container.end();
}

View File

@ -208,3 +208,27 @@ TEST( RiaTimeHistoryCurveMergerTest, NoTimeStepOverlap )
EXPECT_EQ( 0, static_cast<int>( curveMerger.validIntervalsForAllXValues().size() ) ); EXPECT_EQ( 0, static_cast<int>( curveMerger.validIntervalsForAllXValues().size() ) );
} }
} }
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST( RiaTimeHistoryCurveMergerTest, SharedXValues )
{
std::vector<double> valuesA{ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
std::vector<double> valuesB{ 10, 20, 30, 40, 50, 60, 70 };
std::vector<time_t> timeSteps{ 1, 2, 3, 4, 5, 6, 7 };
RiaTimeHistoryCurveMerger interpolate;
interpolate.addCurveData( timeSteps, valuesA );
interpolate.addCurveData( timeSteps, valuesB );
interpolate.computeInterpolatedValues( true );
auto interpolatedTimeSteps = interpolate.allXValues();
EXPECT_TRUE( std::equal( timeSteps.begin(), timeSteps.end(), interpolatedTimeSteps.begin() ) );
auto generatedYValuesA = interpolate.interpolatedYValuesForAllXValues( 0 );
EXPECT_TRUE( std::equal( valuesA.begin(), valuesA.end(), generatedYValuesA.begin() ) );
auto generatedYValuesB = interpolate.interpolatedYValuesForAllXValues( 1 );
EXPECT_TRUE( std::equal( valuesB.begin(), valuesB.end(), generatedYValuesB.begin() ) );
}