From 1ce2f9651c31df3a4f2b1ce659a00006de167b67 Mon Sep 17 00:00:00 2001 From: Magne Sjaastad Date: Tue, 3 Dec 2024 16:53:11 +0100 Subject: [PATCH] #11954 Improve the curve merger algorithm * Add flag to indicate identical X values when curves are merged * Detect if x values are monotonically increasing --- .../Application/Tools/RiaCurveMerger.h | 6 +- .../Application/Tools/RiaCurveMerger.inl | 87 ++++++++++++++++--- .../RigTimeCurveHistoryMerger-Test.cpp | 24 +++++ 3 files changed, 103 insertions(+), 14 deletions(-) diff --git a/ApplicationLibCode/Application/Tools/RiaCurveMerger.h b/ApplicationLibCode/Application/Tools/RiaCurveMerger.h index b6a4965f17..f8bf950288 100644 --- a/ApplicationLibCode/Application/Tools/RiaCurveMerger.h +++ b/ApplicationLibCode/Application/Tools/RiaCurveMerger.h @@ -67,7 +67,8 @@ public: interpolatedYValue( const XValueType& xValue, const std::vector& curveXValues, const std::vector& curveYValues ); private: - void computeUnionOfXValues( bool includeValuesFromPartialCurves ); + void computeUnionOfXValues( bool includeValuesFromPartialCurves ); + static bool isMonotonicallyIncreasing( const std::vector& curveXValues ); private: std::vector, std::vector>> m_originalValues; @@ -76,6 +77,9 @@ private: std::vector m_allXValues; std::vector> m_interpolatedValuesForAllCurves; + + bool m_isXValuesSharedBetweenCurves; + bool m_isXValuesMonotonicallyIncreasing; }; using RiaTimeHistoryCurveMerger = RiaCurveMerger; diff --git a/ApplicationLibCode/Application/Tools/RiaCurveMerger.inl b/ApplicationLibCode/Application/Tools/RiaCurveMerger.inl index d0e9fe06fd..b56590738c 100644 --- a/ApplicationLibCode/Application/Tools/RiaCurveMerger.inl +++ b/ApplicationLibCode/Application/Tools/RiaCurveMerger.inl @@ -55,6 +55,8 @@ bool XValueComparator::equals( const XValueType& lhs, const XValueTy //-------------------------------------------------------------------------------------------------- template RiaCurveMerger::RiaCurveMerger() + : m_isXValuesSharedBetweenCurves( false ) + , m_isXValuesMonotonicallyIncreasing( true ) { } @@ -68,6 +70,24 @@ void RiaCurveMerger::addCurveData( const std::vector& xV 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 ) ); } } @@ -131,14 +151,14 @@ void RiaCurveMerger::computeInterpolatedValues( bool includeValuesFr m_allXValues.clear(); m_interpolatedValuesForAllCurves.clear(); - computeUnionOfXValues( includeValuesFromPartialCurves ); - const size_t curveCount = m_originalValues.size(); if ( curveCount == 0 ) { return; } + computeUnionOfXValues( includeValuesFromPartialCurves ); + const size_t dataValueCount = m_allXValues.size(); if ( dataValueCount == 0 ) { @@ -157,8 +177,18 @@ void RiaCurveMerger::computeInterpolatedValues( bool includeValuesFr #pragma omp parallel for for ( int valueIndex = 0; valueIndex < static_cast( dataValueCount ); valueIndex++ ) { - double interpolValue = - interpolatedYValue( m_allXValues[valueIndex], m_originalValues[curveIdx].first, m_originalValues[curveIdx].second ); + 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 ); + } + if ( !RiaCurveDataTools::isValidValue( interpolValue, false ) ) { #pragma omp critical @@ -180,6 +210,13 @@ void RiaCurveMerger::computeUnionOfXValues( bool includeValuesForPar { 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 unionOfXValues; std::vector> originalXBounds; @@ -240,17 +277,27 @@ double RiaCurveMerger::interpolatedYValue( const XValueType& if ( xValues.empty() ) return HUGE_VAL; if ( yValues.size() != xValues.size() ) return HUGE_VAL; + size_t startIndex = 0; + + if ( isMonotonicallyIncreasing( xValues ) ) + { + // 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(); + auto it = std::lower_bound( xValues.begin(), xValues.end(), interpolationXValue - threshold ); + if ( it == xValues.end() ) return HUGE_VAL; + + startIndex = it - xValues.begin(); + if ( startIndex > 0 ) startIndex--; + } + const bool removeInterpolatedValues = false; - // 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 - XValueType threshold = 1.0e-6 * xValues.back(); - auto it = std::lower_bound( xValues.begin(), xValues.end(), interpolationXValue - threshold ); - if ( it == xValues.end() ) return HUGE_VAL; - - size_t startIndex = it - xValues.begin(); - if ( startIndex > 0 ) startIndex--; - for ( size_t firstI = startIndex; firstI < xValues.size(); firstI++ ) { if ( XComparator::equals( xValues.at( firstI ), interpolationXValue ) ) @@ -307,3 +354,17 @@ double RiaCurveMerger::interpolatedYValue( const XValueType& return HUGE_VAL; } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +template +bool RiaCurveMerger::isMonotonicallyIncreasing( const std::vector& 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(); +} diff --git a/ApplicationLibCode/UnitTests/RigTimeCurveHistoryMerger-Test.cpp b/ApplicationLibCode/UnitTests/RigTimeCurveHistoryMerger-Test.cpp index 45426a5762..b4d485a3c8 100644 --- a/ApplicationLibCode/UnitTests/RigTimeCurveHistoryMerger-Test.cpp +++ b/ApplicationLibCode/UnitTests/RigTimeCurveHistoryMerger-Test.cpp @@ -208,3 +208,27 @@ TEST( RiaTimeHistoryCurveMergerTest, NoTimeStepOverlap ) EXPECT_EQ( 0, static_cast( curveMerger.validIntervalsForAllXValues().size() ) ); } } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RiaTimeHistoryCurveMergerTest, SharedXValues ) +{ + std::vector valuesA{ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 }; + std::vector valuesB{ 10, 20, 30, 40, 50, 60, 70 }; + std::vector 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() ) ); +}