#5470 Add a new median calculator and make rkbdiff default to median well path rkb

ApplicationCode/Application/Tools/CMakeLists_files.cmake

@@ -32,6 +32,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RiaTimeHistoryCurveResampler.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaStatisticsTools.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaOffshoreSphericalCoords.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedMeanCalculator.h
+${CMAKE_CURRENT_LIST_DIR}/RiaMedianCalculator.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedMeanCalculator.inl
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedGeometricMeanCalculator.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedHarmonicMeanCalculator.h

ApplicationCode/Application/Tools/RiaMedianCalculator.h

@@ -0,0 +1,77 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2020-     Equinor ASA
+//
+//  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 <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+#pragma once
+
+#include "cafAssert.h"
+
+#include <algorithm>
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+template <class T>
+class RiaMedianCalculator
+{
+public:
+    RiaMedianCalculator() = default;
+
+    void   add( T value );
+    double median() const;
+    bool   valid() const;
+
+private:
+    std::vector<T> m_values;
+};
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+template <class T>
+void RiaMedianCalculator<T>::add( T value )
+{
+    m_values.insert( std::lower_bound( m_values.begin(), m_values.end(), value ), value );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+template <class T>
+double RiaMedianCalculator<T>::median() const
+{
+    CAF_ASSERT( valid() && "You must ensure the result is valid before calling this" );
+
+    auto count = m_values.size();
+    if ( count == 1u )
+        return m_values.front();
+    else if ( count % 2 == 0 )
+    {
+        return ( m_values[count / 2 - 1] + m_values[count / 2] ) * 0.5;
+    }
+    else
+    {
+        return m_values[count / 2];
+    }
+}
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+template <class T>
+bool RiaMedianCalculator<T>::valid() const
+{
+    return !m_values.empty();
+}

ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp

@@ -20,6 +20,10 @@
 
 #include "RimGeoMechResultDefinition.h"
 
+#include "RiaApplication.h"
+#include "RiaDefines.h"
+#include "RiaMedianCalculator.h"
+
 #include "RifGeoMechReaderInterface.h"
 
 #include "RigFemPartCollection.h"
@@ -29,8 +33,7 @@
 #include "RigFormationNames.h"
 #include "RigGeoMechCaseData.h"
 #include "RigWbsParameter.h"
-
+#include "RigWellPath.h"
-#include "RiaDefines.h"
 
 #include "Rim3dWellLogCurve.h"
 #include "RimGeoMechCase.h"
@@ -39,11 +42,16 @@
 #include "RimGeoMechView.h"
 #include "RimIntersectionCollection.h"
 #include "RimPlotCurve.h"
+#include "RimProject.h"
 #include "RimRegularLegendConfig.h"
 #include "RimViewLinker.h"
+#include "RimWellPath.h"
 
+#include "cafPdmUiDoubleValueEditor.h"
 #include "cafPdmUiListEditor.h"
 
+#include <QDoubleValidator>
+
 namespace caf
 {
 template <>
@@ -107,13 +115,15 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition( void )
                        "",
                        "" );
     CAF_PDM_InitField( &m_normalizationRkbDiff, "NormalizationRkbDiff", -1.0, "Air Gap", "", "", "" );
+    m_normalizationRkbDiff.uiCapability()->setUiEditorTypeName( caf::PdmUiDoubleValueEditor::uiEditorTypeName() );
 
     CAF_PDM_InitField( &m_compactionRefLayerUiField,
                        "CompactionRefLayerUi",
                        RigFemResultAddress::noCompactionValue(),
                        "Compaction Ref Layer",
                        "",
-                       "The compaction is calculated with reference to this layer. Default layer is the topmost layer "
+                       "The compaction is calculated with reference to this layer. Default layer is the topmost "
+                       "layer "
                        "with POR",
                        "" );
     m_compactionRefLayerUiField.xmlCapability()->disableIO();
@@ -151,7 +161,10 @@ void RimGeoMechResultDefinition::defineUiOrdering( QString uiConfigName, caf::Pd
     {
         caf::PdmUiGroup* normalizationGroup = uiOrdering.addNewGroup( "Result Normalization" );
         normalizationGroup->add( &m_normalizeByHydrostaticPressure );
-        normalizationGroup->add( &m_normalizationRkbDiff );
+        if ( m_normalizeByHydrostaticPressure )
+        {
+            normalizationGroup->add( &m_normalizationRkbDiff );
+        }
     }
 
     if ( m_resultPositionTypeUiField() != RIG_FORMATION_NAMES )
@@ -325,6 +338,29 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
             m_resultVariableUiField = "";
         }
     }
+    if ( &m_normalizeByHydrostaticPressure == changedField && m_normalizationRkbDiff < 0.0 )
+    {
+        RiaMedianCalculator<double> rkbDiffCalc;
+        for ( auto wellPath : RiaApplication::instance()->project()->allWellPaths() )
+        {
+            if ( wellPath->wellPathGeometry() )
+            {
+                double rkbDiff = wellPath->wellPathGeometry()->rkbDiff();
+                if ( rkbDiff > 0.0 )
+                {
+                    rkbDiffCalc.add( rkbDiff );
+                }
+            }
+        }
+        if ( rkbDiffCalc.valid() )
+        {
+            m_normalizationRkbDiff = rkbDiffCalc.median();
+        }
+        else
+        {
+            m_normalizationRkbDiff = 0.0;
+        }
+    }
 
     // Get the possible property filter owner
     RimGeoMechPropertyFilter* propFilter = dynamic_cast<RimGeoMechPropertyFilter*>( this->parentField()->ownerObject() );
@@ -492,6 +528,24 @@ void RimGeoMechResultDefinition::initAfterRead()
     m_timeLapseBaseTimestep.uiCapability()->setUiReadOnly( resultPositionType() == RIG_WELLPATH_DERIVED );
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimGeoMechResultDefinition::defineEditorAttribute( const caf::PdmFieldHandle* field,
+                                                        QString                    uiConfigName,
+                                                        caf::PdmUiEditorAttribute* attribute )
+{
+    if ( field == &m_normalizationRkbDiff )
+    {
+        auto attr = dynamic_cast<caf::PdmUiDoubleValueEditorAttribute*>( attribute );
+        if ( attr )
+        {
+            attr->m_decimals  = 2;
+            attr->m_validator = new QDoubleValidator( 0.0, std::numeric_limits<double>::max(), 2 );
+        }
+    }
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h

@@ -95,7 +95,9 @@ private:
                                                     const QVariant&            oldValue,
                                                     const QVariant&            newValue ) override;
     void                          initAfterRead() override;
-
+    void                          defineEditorAttribute( const caf::PdmFieldHandle* field,
+                                                         QString                    uiConfigName,
+                                                         caf::PdmUiEditorAttribute* attribute ) override;
     // Metadata and option build tools
 
     std::map<std::string, std::vector<std::string>> getResultMetaDataForUIFieldSetting();

ApplicationCode/UnitTests/CMakeLists_files.cmake

@@ -49,6 +49,7 @@ ${CMAKE_CURRENT_LIST_DIR}/SolveSpaceSolver-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaPolyArcLineSampler-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RifTextDataTableFormatter-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedMean-Test.cpp
+${CMAKE_CURRENT_LIST_DIR}/RiaMedianCalculator-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedGeometricMeanCalculator-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedHarmonicMeanCalculator-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaCellDividingTools-Test.cpp

ApplicationCode/UnitTests/RiaMedianCalculator-Test.cpp

@@ -0,0 +1,127 @@
+#include "gtest/gtest.h"
+
+#include "RiaMedianCalculator.h"
+
+#include <algorithm>
+#include <random>
+#include <vector>
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+TEST( RiaMedianCalculator, EmptyDataSet )
+{
+    RiaMedianCalculator<double> calc;
+
+    EXPECT_FALSE( calc.valid() );
+}
+
+TEST( RiaMedianCalculator, SingleValue )
+{
+    RiaMedianCalculator<double> calc;
+    calc.add( 4.0 );
+    EXPECT_DOUBLE_EQ( 4.0, calc.median() );
+}
+
+TEST( RiaMedianCalculator, TwoValues )
+{
+    RiaMedianCalculator<double> calc;
+
+    std::vector<double> values{3.0, 6.0};
+
+    for ( size_t i = 0; i < values.size(); i++ )
+    {
+        calc.add( values[i] );
+    }
+    EXPECT_TRUE( calc.valid() );
+    EXPECT_DOUBLE_EQ( 4.5, calc.median() );
+}
+
+TEST( RiaMedianCalculator, ThreeValues )
+{
+    RiaMedianCalculator<double> calc;
+
+    std::vector<double> values{3.0, 6.0, 2.0};
+
+    for ( size_t i = 0; i < values.size(); i++ )
+    {
+        calc.add( values[i] );
+    }
+    EXPECT_TRUE( calc.valid() );
+    EXPECT_DOUBLE_EQ( 3.0, calc.median() );
+}
+
+TEST( RiaMedianCalculator, SameValues )
+{
+    RiaMedianCalculator<double> calc;
+    const double                value = 113.0;
+    for ( size_t i = 0; i < 1000; ++i )
+    {
+        calc.add( value );
+    }
+    EXPECT_TRUE( calc.valid() );
+    EXPECT_DOUBLE_EQ( value, calc.median() );
+}
+
+TEST( RiaMedianCalculator, OrderedRangeOdd )
+{
+    RiaMedianCalculator<double> calc;
+    for ( size_t i = 0; i < 101; ++i )
+    {
+        calc.add( (double)i );
+    }
+    EXPECT_TRUE( calc.valid() );
+    EXPECT_DOUBLE_EQ( 50.0, calc.median() );
+}
+
+TEST( RiaMedianCalculator, OrderedRangeEven )
+{
+    RiaMedianCalculator<double> calc;
+    for ( size_t i = 0; i < 100; ++i )
+    {
+        calc.add( (double)i );
+    }
+    EXPECT_TRUE( calc.valid() );
+    EXPECT_DOUBLE_EQ( 49.5, calc.median() );
+}
+
+TEST( RiaMedianCalculator, ShuffledRangeOdd )
+{
+    RiaMedianCalculator<double> calc;
+
+    std::vector<int> values;
+    for ( int i = 0; i < 101; ++i )
+    {
+        values.push_back( i );
+    }
+    std::random_device rd;
+    std::mt19937       g( rd() );
+    std::shuffle( values.begin(), values.end(), g );
+    for ( double value : values )
+    {
+        calc.add( value );
+    }
+    EXPECT_TRUE( calc.valid() );
+    EXPECT_DOUBLE_EQ( 50.0, calc.median() );
+}
+
+TEST( RiaMedianCalculator, ShuffledRangeEven )
+{
+    RiaMedianCalculator<double> calc;
+
+    std::vector<int> values;
+    for ( int i = 0; i < 200; ++i )
+    {
+        values.push_back( i );
+    }
+    std::random_device rd;
+    std::mt19937       g( rd() );
+    std::shuffle( values.begin(), values.end(), g );
+    for ( double value : values )
+    {
+        calc.add( value );
+    }
+
+    EXPECT_TRUE( calc.valid() );
+    EXPECT_DOUBLE_EQ( 99.5, calc.median() );
+}