From 41999d84637eae73ed891881a59014c364fecf84 Mon Sep 17 00:00:00 2001 From: Wim Haeck Date: Mon, 29 Jul 2024 16:50:05 -0600 Subject: [PATCH] Adding single zone histogram tabulated function --- cmake/unit_testing.cmake | 1 + src/scion/math/HistogramTableFunction.hpp | 75 ++++++ .../math/HistogramTableFunction/src/ctor.hpp | 8 + .../test/CMakeLists.txt | 1 + .../test/HistogramTableFunction.test.cpp | 249 ++++++++++++++++++ src/scion/math/LogLinearTableFunction.hpp | 4 +- 6 files changed, 336 insertions(+), 2 deletions(-) create mode 100644 src/scion/math/HistogramTableFunction.hpp create mode 100644 src/scion/math/HistogramTableFunction/src/ctor.hpp create mode 100644 src/scion/math/HistogramTableFunction/test/CMakeLists.txt create mode 100644 src/scion/math/HistogramTableFunction/test/HistogramTableFunction.test.cpp diff --git a/cmake/unit_testing.cmake b/cmake/unit_testing.cmake index eef9374..114621a 100644 --- a/cmake/unit_testing.cmake +++ b/cmake/unit_testing.cmake @@ -77,6 +77,7 @@ add_subdirectory( src/scion/math/OpenDomain/test ) add_subdirectory( src/scion/math/LinearLinearTable/test ) add_subdirectory( src/scion/math/LinearLinearTableFunction/test ) add_subdirectory( src/scion/math/HistogramTable/test ) +add_subdirectory( src/scion/math/HistogramTableFunction/test ) add_subdirectory( src/scion/math/LogLinearTable/test ) add_subdirectory( src/scion/math/LogLinearTableFunction/test ) add_subdirectory( src/scion/math/LinearLogTable/test ) diff --git a/src/scion/math/HistogramTableFunction.hpp b/src/scion/math/HistogramTableFunction.hpp new file mode 100644 index 0000000..a04dcae --- /dev/null +++ b/src/scion/math/HistogramTableFunction.hpp @@ -0,0 +1,75 @@ +#ifndef NJOY_SCION_MATH_HISTOGRAMTABLEFUNCTION +#define NJOY_SCION_MATH_HISTOGRAMTABLEFUNCTION + +// system includes +#include + +// other includes +#include "scion/interpolation/InterpolationType.hpp" +#include "scion/interpolation/Histogram.hpp" +#include "scion/linearisation/ToleranceConvergence.hpp" +#include "scion/math/SingleTableFunctionBase.hpp" + +namespace njoy { +namespace scion { +namespace math { + + /** + * @class + * @brief Tabulated x,f(y) data with linear-linear interpolation (f(y) is linear in x) + * + * The HistogramTableFunction is templatised on the actual x, y and z types, the function + * type F and the container type used for the x values and the functions. This allows us to + * use something like utility::IteratorView instead of std::vector. + */ + template < typename X, typename Y, typename Z, typename F, + typename XContainer = std::vector< X >, + typename FContainer = std::vector< F > > + class HistogramTableFunction : + public SingleTableFunctionBase< HistogramTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::Histogram, X, Y, Z, F, + XContainer, FContainer > { + + /* friend declarations */ + friend class SingleTableFunctionBase< HistogramTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::Histogram, X, Y, Z, F, + XContainer, FContainer >; + + /* type aliases */ + using Parent = SingleTableFunctionBase< HistogramTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::Histogram, X, Y, Z, F, + XContainer, FContainer >; + + /* fields */ + + /* auxiliary function */ + + /* interface implementation functions */ + + /** + * @brief Return the interpolation type + */ + static constexpr interpolation::InterpolationType type() noexcept { + + return interpolation::InterpolationType::Histogram; + } + + public: + + /* constructor */ + #include "scion/math/HistogramTableFunction/src/ctor.hpp" + + /* methods */ + + using Parent::interpolation; + using Parent::x; + using Parent::f; + using Parent::numberPoints; + using Parent::operator(); + }; + +} // math namespace +} // scion namespace +} // njoy namespace + +#endif diff --git a/src/scion/math/HistogramTableFunction/src/ctor.hpp b/src/scion/math/HistogramTableFunction/src/ctor.hpp new file mode 100644 index 0000000..14b73fc --- /dev/null +++ b/src/scion/math/HistogramTableFunction/src/ctor.hpp @@ -0,0 +1,8 @@ +/** + * @brief Constructor + * + * @param x the x values of the tabulated data + * @param f the f(y) functions of the tabulated data + */ +HistogramTableFunction( XContainer x, FContainer f ) : + Parent( std::move( x ), std::move( f ) ) {} diff --git a/src/scion/math/HistogramTableFunction/test/CMakeLists.txt b/src/scion/math/HistogramTableFunction/test/CMakeLists.txt new file mode 100644 index 0000000..445364e --- /dev/null +++ b/src/scion/math/HistogramTableFunction/test/CMakeLists.txt @@ -0,0 +1 @@ +add_cpp_test( math.HistogramTableFunction HistogramTableFunction.test.cpp ) diff --git a/src/scion/math/HistogramTableFunction/test/HistogramTableFunction.test.cpp b/src/scion/math/HistogramTableFunction/test/HistogramTableFunction.test.cpp new file mode 100644 index 0000000..b9fcf42 --- /dev/null +++ b/src/scion/math/HistogramTableFunction/test/HistogramTableFunction.test.cpp @@ -0,0 +1,249 @@ +// include Catch2 +#include +#include +using Catch::Matchers::WithinRel; + +// what we are testing +#include "scion/math/InterpolationTable.hpp" +#include "scion/math/HistogramTableFunction.hpp" + +// other includes +#include "utility/IteratorView.hpp" + +// convenience typedefs +using namespace njoy::scion; +template < typename X, typename Y = X > +using InterpolationTable = math::InterpolationTable< X, Y >; +template < typename X, typename Y = X, typename Z = X, + typename F = InterpolationTable< X >, + typename XContainer = std::vector< X >, + typename FContainer = std::vector< F > > +using HistogramTableFunction = math::HistogramTableFunction< X, Y, Z, F, XContainer, FContainer >; +using InterpolationType = interpolation::InterpolationType; + +SCENARIO( "HistogramTableFunction" ) { + + GIVEN( "tabulated x,f(y) data" ) { + + WHEN( "the data is given explicitly using vectors" ) { + + const std::vector< double > x = { 1., 2., 3., 4. }; + const std::vector< InterpolationTable< double > > f = { + + { { -1., +1. }, { 0.5, 0.5 } }, + { { -1., 0., +1. }, { 0.49, 0.5, 0.51 } }, + { { -1., 0., +1. }, { 0.4, 0.5, 0.6 } }, + { { -1., +1. }, { 0.1, 0.9 } } + }; + + HistogramTableFunction< double > chunk( std::move( x ), std::move( f ) ); + + THEN( "a HistogramTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::Histogram == chunk.interpolation() ); + CHECK( 4 == chunk.numberPoints() ); + CHECK( 4 == chunk.x().size() ); + CHECK( 4 == chunk.f().size() ); + CHECK_THAT( 1., WithinRel( chunk.x()[0] ) ); + CHECK_THAT( 2., WithinRel( chunk.x()[1] ) ); + CHECK_THAT( 3., WithinRel( chunk.x()[2] ) ); + CHECK_THAT( 4., WithinRel( chunk.x()[3] ) ); + CHECK( 2 == chunk.f()[0].x().size() ); + CHECK( 2 == chunk.f()[0].y().size() ); + CHECK( 3 == chunk.f()[1].x().size() ); + CHECK( 3 == chunk.f()[1].y().size() ); + CHECK( 3 == chunk.f()[2].x().size() ); + CHECK( 3 == chunk.f()[2].y().size() ); + CHECK( 2 == chunk.f()[3].x().size() ); + CHECK( 2 == chunk.f()[3].y().size() ); + CHECK_THAT( -1. , WithinRel( chunk.f()[0].x()[0] ) ); + CHECK_THAT( 1. , WithinRel( chunk.f()[0].x()[1] ) ); + CHECK_THAT( 0.5 , WithinRel( chunk.f()[0].y()[0] ) ); + CHECK_THAT( 0.5 , WithinRel( chunk.f()[0].y()[1] ) ); + CHECK_THAT( -1. , WithinRel( chunk.f()[1].x()[0] ) ); + CHECK_THAT( 0. , WithinRel( chunk.f()[1].x()[1] ) ); + CHECK_THAT( 1. , WithinRel( chunk.f()[1].x()[2] ) ); + CHECK_THAT( 0.49, WithinRel( chunk.f()[1].y()[0] ) ); + CHECK_THAT( 0.5 , WithinRel( chunk.f()[1].y()[1] ) ); + CHECK_THAT( 0.51, WithinRel( chunk.f()[1].y()[2] ) ); + CHECK_THAT( -1. , WithinRel( chunk.f()[2].x()[0] ) ); + CHECK_THAT( 0. , WithinRel( chunk.f()[2].x()[1] ) ); + CHECK_THAT( 1. , WithinRel( chunk.f()[2].x()[2] ) ); + CHECK_THAT( 0.4 , WithinRel( chunk.f()[2].y()[0] ) ); + CHECK_THAT( 0.5 , WithinRel( chunk.f()[2].y()[1] ) ); + CHECK_THAT( 0.6 , WithinRel( chunk.f()[2].y()[2] ) ); + CHECK_THAT( -1. , WithinRel( chunk.f()[3].x()[0] ) ); + CHECK_THAT( 1. , WithinRel( chunk.f()[3].x()[1] ) ); + CHECK_THAT( 0.1, WithinRel( chunk.f()[3].y()[0] ) ); + CHECK_THAT( 0.9, WithinRel( chunk.f()[3].y()[1] ) ); + } // THEN + + THEN( "a HistogramTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.5 , WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 4., -0.5 ) ) ); + + // values of x outside the x grid + CHECK_THAT( 0., WithinRel( chunk( 0., -0.5 ) ) ); + CHECK_THAT( 0., WithinRel( chunk( 5., -0.5 ) ) ); + + // values of x inside the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3.5, -0.5 ) ) ); + } // THEN + } // WHEN + + WHEN( "the data is given explicitly using iterator views" ) { + + using XView = njoy::utility::IteratorView< std::vector< double >::const_iterator >; + using FView = njoy::utility::IteratorView< std::vector< InterpolationTable< double > >::const_iterator >; + + const std::vector< double > xvalues = { 1., 2., 3., 4. }; + const std::vector< InterpolationTable< double > > fvalues = { + + { { -1., +1. }, { 0.5, 0.5 } }, + { { -1., 0., +1. }, { 0.49, 0.5, 0.51 } }, + { { -1., 0., +1. }, { 0.4, 0.5, 0.6 } }, + { { -1., +1. }, { 0.1, 0.9 } } + }; + + + XView x = njoy::utility::make_view( xvalues ); + FView f = njoy::utility::make_view( fvalues ); + + HistogramTableFunction< double, double, double, InterpolationTable< double >, + XView, FView > chunk( std::move( x ), std::move( f ) ); + + THEN( "a HistogramTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::Histogram == chunk.interpolation() ); + CHECK( 4 == chunk.numberPoints() ); + CHECK( 4 == chunk.x().size() ); + CHECK( 4 == chunk.f().size() ); + CHECK_THAT( 1., WithinRel( chunk.x()[0] ) ); + CHECK_THAT( 2., WithinRel( chunk.x()[1] ) ); + CHECK_THAT( 3., WithinRel( chunk.x()[2] ) ); + CHECK_THAT( 4., WithinRel( chunk.x()[3] ) ); + CHECK( 2 == chunk.f()[0].x().size() ); + CHECK( 2 == chunk.f()[0].y().size() ); + CHECK( 3 == chunk.f()[1].x().size() ); + CHECK( 3 == chunk.f()[1].y().size() ); + CHECK( 3 == chunk.f()[2].x().size() ); + CHECK( 3 == chunk.f()[2].y().size() ); + CHECK( 2 == chunk.f()[3].x().size() ); + CHECK( 2 == chunk.f()[3].y().size() ); + CHECK_THAT( -1. , WithinRel( chunk.f()[0].x()[0] ) ); + CHECK_THAT( 1. , WithinRel( chunk.f()[0].x()[1] ) ); + CHECK_THAT( 0.5 , WithinRel( chunk.f()[0].y()[0] ) ); + CHECK_THAT( 0.5 , WithinRel( chunk.f()[0].y()[1] ) ); + CHECK_THAT( -1. , WithinRel( chunk.f()[1].x()[0] ) ); + CHECK_THAT( 0. , WithinRel( chunk.f()[1].x()[1] ) ); + CHECK_THAT( 1. , WithinRel( chunk.f()[1].x()[2] ) ); + CHECK_THAT( 0.49, WithinRel( chunk.f()[1].y()[0] ) ); + CHECK_THAT( 0.5 , WithinRel( chunk.f()[1].y()[1] ) ); + CHECK_THAT( 0.51, WithinRel( chunk.f()[1].y()[2] ) ); + CHECK_THAT( -1. , WithinRel( chunk.f()[2].x()[0] ) ); + CHECK_THAT( 0. , WithinRel( chunk.f()[2].x()[1] ) ); + CHECK_THAT( 1. , WithinRel( chunk.f()[2].x()[2] ) ); + CHECK_THAT( 0.4 , WithinRel( chunk.f()[2].y()[0] ) ); + CHECK_THAT( 0.5 , WithinRel( chunk.f()[2].y()[1] ) ); + CHECK_THAT( 0.6 , WithinRel( chunk.f()[2].y()[2] ) ); + CHECK_THAT( -1. , WithinRel( chunk.f()[3].x()[0] ) ); + CHECK_THAT( 1. , WithinRel( chunk.f()[3].x()[1] ) ); + CHECK_THAT( 0.1, WithinRel( chunk.f()[3].y()[0] ) ); + CHECK_THAT( 0.9, WithinRel( chunk.f()[3].y()[1] ) ); + } // THEN + + THEN( "a HistogramTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.5 , WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 4., -0.5 ) ) ); + + // values of x outside the x grid + CHECK_THAT( 0., WithinRel( chunk( 0., -0.5 ) ) ); + CHECK_THAT( 0., WithinRel( chunk( 5., -0.5 ) ) ); + + // values of x inside the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3.5, -0.5 ) ) ); + } // THEN + } // WHEN + } // GIVEN + + GIVEN( "invalid data for a HistogramTable object" ) { + + WHEN( "there are not enough values in the x or y grid" ) { + + std::vector< double > xempty = {}; + std::vector< double > xone = { 1. }; + std::vector< InterpolationTable< double > > fempty = {}; + std::vector< InterpolationTable< double > > fone = { { { 1., 2. }, { 3., 4. } } }; + + THEN( "an exception is thrown" ) { + + CHECK_THROWS( HistogramTableFunction< double >( xempty, fempty ) ); + CHECK_THROWS( HistogramTableFunction< double >( xone, fone ) ); + CHECK_THROWS( HistogramTableFunction< double >( xempty, fone ) ); + CHECK_THROWS( HistogramTableFunction< double >( xone, fempty ) ); + } // THEN + } // WHEN + + WHEN( "the x and y grid do not have the same number of points" ) { + + std::vector< double > x = { 1., 2., 3., 4. }; + std::vector< InterpolationTable< double > > f = { + + { { -1., +1. }, { 0.5, 0.5 } }, + { { -1., 0., +1. }, { 0.49, 0.5, 0.51 } }, + { { -1., 0., +1. }, { 0.4, 0.5, 0.6 } } + }; + + THEN( "an exception is thrown" ) { + + CHECK_THROWS( HistogramTableFunction< double >( std::move( x ), std::move( f ) ) ); + } // THEN + } // WHEN + + WHEN( "the x grid is not sorted" ) { + + std::vector< double > x = { 1., 3., 2., 4. }; + std::vector< InterpolationTable< double > > f = { + + { { -1., +1. }, { 0.5, 0.5 } }, + { { -1., 0., +1. }, { 0.49, 0.5, 0.51 } }, + { { -1., 0., +1. }, { 0.4, 0.5, 0.6 } }, + { { -1., +1. }, { 0.1, 0.9 } } + }; + + THEN( "an exception is thrown" ) { + + CHECK_THROWS( HistogramTableFunction< double >( std::move( x ), std::move( f ) ) ); + } // THEN + } // WHEN + + WHEN( "the x grid has a duplicate point" ) { + + std::vector< double > x = { 1., 2., 2., 4. }; + std::vector< InterpolationTable< double > > f = { + + { { -1., +1. }, { 0.5, 0.5 } }, + { { -1., 0., +1. }, { 0.49, 0.5, 0.51 } }, + { { -1., 0., +1. }, { 0.4, 0.5, 0.6 } }, + { { -1., +1. }, { 0.1, 0.9 } } + }; + + THEN( "an exception is thrown" ) { + + CHECK_THROWS( HistogramTableFunction< double >( std::move( x ), std::move( f ) ) ); + } // THEN + } // WHEN + } // GIVEN +} // SCENARIO diff --git a/src/scion/math/LogLinearTableFunction.hpp b/src/scion/math/LogLinearTableFunction.hpp index 85a0620..c5c3e13 100644 --- a/src/scion/math/LogLinearTableFunction.hpp +++ b/src/scion/math/LogLinearTableFunction.hpp @@ -1,5 +1,5 @@ -#ifndef NJOY_SCION_MATH_LogLinearTABLEFUNCTION -#define NJOY_SCION_MATH_LogLinearTABLEFUNCTION +#ifndef NJOY_SCION_MATH_LOGLINEARTABLEFUNCTION +#define NJOY_SCION_MATH_LOGLINEARTABLEFUNCTION // system includes #include