diff --git a/cmake/unit_testing.cmake b/cmake/unit_testing.cmake index c6d8c81..114621a 100644 --- a/cmake/unit_testing.cmake +++ b/cmake/unit_testing.cmake @@ -60,7 +60,6 @@ add_subdirectory( src/scion/interpolation/LinearLinear/test ) add_subdirectory( src/scion/interpolation/LinearLogarithmic/test ) add_subdirectory( src/scion/interpolation/LogarithmicLinear/test ) add_subdirectory( src/scion/interpolation/LogarithmicLogarithmic/test ) -add_subdirectory( src/scion/interpolation/Table/test ) add_subdirectory( src/scion/linearisation/ToleranceConvergence/test ) add_subdirectory( src/scion/linearisation/MidpointSplit/test ) @@ -76,10 +75,15 @@ add_subdirectory( src/scion/math/newton/test ) add_subdirectory( src/scion/math/IntervalDomain/test ) 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 ) +add_subdirectory( src/scion/math/LinearLogTableFunction/test ) add_subdirectory( src/scion/math/LogLogTable/test ) +add_subdirectory( src/scion/math/LogLogTableFunction/test ) add_subdirectory( src/scion/math/InterpolationTable/test ) add_subdirectory( src/scion/math/ChebyshevSeries/test ) add_subdirectory( src/scion/math/ChebyshevApproximation/test ) diff --git a/python/src/definitions.hpp b/python/src/definitions.hpp index 2567236..8853fa5 100644 --- a/python/src/definitions.hpp +++ b/python/src/definitions.hpp @@ -8,7 +8,6 @@ // other includes #include "scion/linearisation/ToleranceConvergence.hpp" -#include "scion/math/FunctionBase.hpp" // namespace aliases namespace python = pybind11; @@ -73,7 +72,7 @@ template < typename Component, typename X, typename Y, typename PythonClass > void addStandardFunctionDefinitions( PythonClass& component ) { using ToleranceConvergence = njoy::scion::linearisation::ToleranceConvergence< X, Y >; - using DomainVariant = typename njoy::scion::math::FunctionBase< X, Y >::DomainVariant; + using DomainVariant = typename Component::DomainVariant; // note: for is_same_domain to bind properly, all possible variant members // must have a default constructor diff --git a/src/scion/interpolation.hpp b/src/scion/interpolation.hpp index d234f45..8fab800 100644 --- a/src/scion/interpolation.hpp +++ b/src/scion/interpolation.hpp @@ -5,4 +5,4 @@ #include "scion/interpolation/LinearLogarithmic.hpp" #include "scion/interpolation/LogarithmicLinear.hpp" #include "scion/interpolation/LogarithmicLogarithmic.hpp" -#include "scion/interpolation/Table.hpp" + diff --git a/src/scion/interpolation/Histogram/test/Histogram.test.cpp b/src/scion/interpolation/Histogram/test/Histogram.test.cpp index 2c7f95e..8efbaa4 100644 --- a/src/scion/interpolation/Histogram/test/Histogram.test.cpp +++ b/src/scion/interpolation/Histogram/test/Histogram.test.cpp @@ -15,7 +15,7 @@ SCENARIO( "Histogram" ) { GIVEN( "Histogram interpolation object" ) { - WHEN( "interpolating an interval" ) { + WHEN( "interpolating an x,y interval" ) { interpolation::Histogram interpolator{}; @@ -31,6 +31,29 @@ SCENARIO( "Histogram" ) { CHECK_THAT( 1., WithinRel( interpolator( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,f(y) interval" ) { + + interpolation::Histogram interpolator{}; + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return - y; }; + + CHECK_THAT( 0.0, WithinRel( interpolator( 1.0, 0.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 0.0, WithinRel( interpolator( 1.5, 0.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 0.0, WithinRel( interpolator( 2.0, 0.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.0, WithinRel( interpolator( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.0, WithinRel( interpolator( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.0, WithinRel( interpolator( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( -1.0, WithinRel( interpolator( 1.0, -1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( -1.0, WithinRel( interpolator( 1.5, -1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( -1.0, WithinRel( interpolator( 2.0, -1.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN GIVEN( "histogram interpolation function" ) { @@ -49,5 +72,26 @@ SCENARIO( "Histogram" ) { CHECK_THAT( 1., WithinRel( interpolation::histogram( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,y interval" ) { + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return - y; }; + + CHECK_THAT( 0.0, WithinRel( interpolation::histogram( 1.0, 0.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 0.0, WithinRel( interpolation::histogram( 1.5, 0.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 0.0, WithinRel( interpolation::histogram( 2.0, 0.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.0, WithinRel( interpolation::histogram( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.0, WithinRel( interpolation::histogram( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.0, WithinRel( interpolation::histogram( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( -1.0, WithinRel( interpolation::histogram( 1.0, -1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( -1.0, WithinRel( interpolation::histogram( 1.5, -1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( -1.0, WithinRel( interpolation::histogram( 2.0, -1.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN } // SCENARIO diff --git a/src/scion/interpolation/InterpolatorBase.hpp b/src/scion/interpolation/InterpolatorBase.hpp index 97c4df8..64a89ed 100644 --- a/src/scion/interpolation/InterpolatorBase.hpp +++ b/src/scion/interpolation/InterpolatorBase.hpp @@ -23,7 +23,7 @@ namespace interpolation { /* methods */ /** - * @brief Perform interpolation + * @brief Perform interpolation on x,y tabulated data * * @param[in] x the value of x * @param[in] xLeft the left value on the x interval @@ -38,6 +38,43 @@ namespace interpolation { return static_cast< const Derived* >( this )->interpolate( x, xLeft, xRight, yLeft, yRight ); } + + /** + * @brief Perform interpolation on x,f(y) tabulated data + * + * This interpolation works by first evaluating the functions z = f(y) for the value + * of y, and then performing the normal interpolation on the interval defined by + * xl,zl and xr,zr as illustrated here: + * + * fl(y) + * xl - - - # - - - - > left function of y + * | + * x - - - * + * | + * xr - - - # - - - - > right function of y + * | + * fr(y) + * + * As such, this does not do unit-base interpolation. + * + * Since we can reuse the operator() on this base class so the derived class + * does not need to implement this particular interpolation operation. + * + * @param[in] x the value of x + * @param[in] y the value of y + * @param[in] xLeft the left value on the x interval + * @param[in] xRight the right value on the x interval + * @param[in] fLeft the left function on the y axis + * @param[in] fRight the right function on the y axis + */ + template < typename X, typename Y, typename FLeft, typename FRight = FLeft, + typename Z = decltype( std::declval< FLeft >()( std::declval< Y >() ) ) > + Z operator()( const X& x, const Y& y, + const X& xLeft, const X& xRight, + const FLeft& fLeft, const FRight& fRight ) const noexcept { + + return this->operator()( x, xLeft, xRight, fLeft( y ), fRight( y ) ); + } }; } // interpolation namespace diff --git a/src/scion/interpolation/LinearLinear/test/LinearLinear.test.cpp b/src/scion/interpolation/LinearLinear/test/LinearLinear.test.cpp index 87c9ffc..49ea9b0 100644 --- a/src/scion/interpolation/LinearLinear/test/LinearLinear.test.cpp +++ b/src/scion/interpolation/LinearLinear/test/LinearLinear.test.cpp @@ -15,7 +15,7 @@ SCENARIO( "LinearLinear" ) { GIVEN( "LinearLinear interpolation object" ) { - WHEN( "interpolating an interval" ) { + WHEN( "interpolating an x,y interval" ) { interpolation::LinearLinear interpolator{}; @@ -31,11 +31,34 @@ SCENARIO( "LinearLinear" ) { CHECK_THAT( 4.0, WithinRel( interpolator( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,f(y) interval" ) { + + interpolation::LinearLinear interpolator{}; + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return 2 * y; }; + + CHECK_THAT( 1.0 , WithinRel( interpolator( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.5 , WithinRel( interpolator( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.0 , WithinRel( interpolator( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.5 , WithinRel( interpolator( 1.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 3.75, WithinRel( interpolator( 1.5, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.0 , WithinRel( interpolator( 2.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 4.0 , WithinRel( interpolator( 1.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 6.0 , WithinRel( interpolator( 1.5, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 8.0 , WithinRel( interpolator( 2.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN GIVEN( "linlin interpolation function" ) { - WHEN( "interpolating an interval" ) { + WHEN( "interpolating an x,y interval" ) { THEN( "the interpolation is performed correctly" ) { @@ -49,5 +72,26 @@ SCENARIO( "LinearLinear" ) { CHECK_THAT( 4.0, WithinRel( interpolation::linlin( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,y interval" ) { + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return 2 * y; }; + + CHECK_THAT( 1.0 , WithinRel( interpolation::linlin( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.5 , WithinRel( interpolation::linlin( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.0 , WithinRel( interpolation::linlin( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.5 , WithinRel( interpolation::linlin( 1.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 3.75, WithinRel( interpolation::linlin( 1.5, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.0 , WithinRel( interpolation::linlin( 2.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 4.0 , WithinRel( interpolation::linlin( 1.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 6.0 , WithinRel( interpolation::linlin( 1.5, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 8.0 , WithinRel( interpolation::linlin( 2.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN } // SCENARIO diff --git a/src/scion/interpolation/LinearLogarithmic/test/LinearLogarithmic.test.cpp b/src/scion/interpolation/LinearLogarithmic/test/LinearLogarithmic.test.cpp index 896052c..8f1bb89 100644 --- a/src/scion/interpolation/LinearLogarithmic/test/LinearLogarithmic.test.cpp +++ b/src/scion/interpolation/LinearLogarithmic/test/LinearLogarithmic.test.cpp @@ -7,7 +7,8 @@ using Catch::Matchers::WithinRel; #include "scion/interpolation/LinearLogarithmic.hpp" // other includes - +#include +#include // convenience typedefs using namespace njoy::scion; @@ -15,7 +16,7 @@ SCENARIO( "LinearLogarithmic" ) { GIVEN( "LinearLogarithmic interpolation object" ) { - WHEN( "interpolating an interval" ) { + WHEN( "interpolating an x,y interval" ) { interpolation::LinearLogarithmic interpolator{}; @@ -31,6 +32,29 @@ SCENARIO( "LinearLogarithmic" ) { CHECK_THAT( 4.0 , WithinRel( interpolator( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,f(y) interval" ) { + + interpolation::LinearLogarithmic interpolator{}; + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return 2 * y; }; + + CHECK_THAT( 1.0 , WithinRel( interpolator( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.58496250072116, WithinRel( interpolator( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.0 , WithinRel( interpolator( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.5 , WithinRel( interpolator( 1.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 3.96240625180289, WithinRel( interpolator( 1.5, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.0 , WithinRel( interpolator( 2.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 4.0 , WithinRel( interpolator( 1.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 6.33985000288462, WithinRel( interpolator( 1.5, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 8.0 , WithinRel( interpolator( 2.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN GIVEN( "linlog interpolation function" ) { @@ -49,5 +73,26 @@ SCENARIO( "LinearLogarithmic" ) { CHECK_THAT( 4.0 , WithinRel( interpolation::linlog( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,y interval" ) { + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return 2 * y; }; + + CHECK_THAT( 1.0 , WithinRel( interpolation::linlog( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.58496250072116, WithinRel( interpolation::linlog( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.0 , WithinRel( interpolation::linlog( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.5 , WithinRel( interpolation::linlog( 1.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 3.96240625180289, WithinRel( interpolation::linlog( 1.5, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.0 , WithinRel( interpolation::linlog( 2.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 4.0 , WithinRel( interpolation::linlog( 1.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 6.33985000288462, WithinRel( interpolation::linlog( 1.5, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 8.0 , WithinRel( interpolation::linlog( 2.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN } // SCENARIO diff --git a/src/scion/interpolation/LogarithmicLinear/test/LogarithmicLinear.test.cpp b/src/scion/interpolation/LogarithmicLinear/test/LogarithmicLinear.test.cpp index 05160f1..ef5fb90 100644 --- a/src/scion/interpolation/LogarithmicLinear/test/LogarithmicLinear.test.cpp +++ b/src/scion/interpolation/LogarithmicLinear/test/LogarithmicLinear.test.cpp @@ -15,7 +15,7 @@ SCENARIO( "LogarithmicLinear" ) { GIVEN( "LogarithmicLinear interpolation object" ) { - WHEN( "interpolating an interval" ) { + WHEN( "interpolating an x,y interval" ) { interpolation::LogarithmicLinear interpolator{}; @@ -31,6 +31,29 @@ SCENARIO( "LogarithmicLinear" ) { CHECK_THAT( 4., WithinRel( interpolator( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,f(y) interval" ) { + + interpolation::LogarithmicLinear interpolator{}; + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return 2 * y; }; + + CHECK_THAT( 1.0 , WithinRel( interpolator( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.41421356237309, WithinRel( interpolator( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.0 , WithinRel( interpolator( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.5 , WithinRel( interpolator( 1.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 3.53553390593274, WithinRel( interpolator( 1.5, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.0 , WithinRel( interpolator( 2.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 4.0 , WithinRel( interpolator( 1.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.65685424949238, WithinRel( interpolator( 1.5, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 8.0 , WithinRel( interpolator( 2.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN GIVEN( "loglin interpolation function" ) { @@ -49,5 +72,26 @@ SCENARIO( "LogarithmicLinear" ) { CHECK_THAT( 4., WithinRel( interpolation::loglin( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,y interval" ) { + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return 2 * y; }; + + CHECK_THAT( 1.0 , WithinRel( interpolation::loglin( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.41421356237309, WithinRel( interpolation::loglin( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.0 , WithinRel( interpolation::loglin( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.5 , WithinRel( interpolation::loglin( 1.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 3.53553390593274, WithinRel( interpolation::loglin( 1.5, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.0 , WithinRel( interpolation::loglin( 2.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 4.0 , WithinRel( interpolation::loglin( 1.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.65685424949238, WithinRel( interpolation::loglin( 1.5, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 8.0 , WithinRel( interpolation::loglin( 2.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN } // SCENARIO diff --git a/src/scion/interpolation/LogarithmicLogarithmic/test/LogarithmicLogarithmic.test.cpp b/src/scion/interpolation/LogarithmicLogarithmic/test/LogarithmicLogarithmic.test.cpp index 5d71db7..15b93b5 100644 --- a/src/scion/interpolation/LogarithmicLogarithmic/test/LogarithmicLogarithmic.test.cpp +++ b/src/scion/interpolation/LogarithmicLogarithmic/test/LogarithmicLogarithmic.test.cpp @@ -15,7 +15,7 @@ SCENARIO( "LogarithmicLogarithmic" ) { GIVEN( "LogarithmicLogarithmic interpolation object" ) { - WHEN( "interpolating an interval" ) { + WHEN( "interpolating an x,y interval" ) { interpolation::LogarithmicLogarithmic interpolator{}; @@ -31,6 +31,29 @@ SCENARIO( "LogarithmicLogarithmic" ) { CHECK_THAT( 4.0 , WithinRel( interpolator( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,f(y) interval" ) { + + interpolation::LogarithmicLogarithmic interpolator{}; + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return 2 * y; }; + + CHECK_THAT( 1.0 , WithinRel( interpolator( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.5 , WithinRel( interpolator( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.0 , WithinRel( interpolator( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.5 , WithinRel( interpolator( 1.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 3.75, WithinRel( interpolator( 1.5, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.0 , WithinRel( interpolator( 2.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 4.0 , WithinRel( interpolator( 1.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 6.0 , WithinRel( interpolator( 1.5, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 8.0 , WithinRel( interpolator( 2.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN GIVEN( "loglog interpolation function" ) { @@ -49,5 +72,26 @@ SCENARIO( "LogarithmicLogarithmic" ) { CHECK_THAT( 4.0 , WithinRel( interpolation::loglog( 2.0, xLeft, xRight, yLeft, yRight ) ) ); } // THEN } // WHEN + + WHEN( "interpolating an x,y interval" ) { + + THEN( "the interpolation is performed correctly" ) { + + double xLeft = 1.0; + double xRight = 2.0; + auto fLeft = [] ( double y ) { return y; }; + auto fRight = [] ( double y ) { return 2 * y; }; + + CHECK_THAT( 1.0 , WithinRel( interpolation::loglog( 1.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 1.5 , WithinRel( interpolation::loglog( 1.5, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.0 , WithinRel( interpolation::loglog( 2.0, 1.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 2.5 , WithinRel( interpolation::loglog( 1.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 3.75, WithinRel( interpolation::loglog( 1.5, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 5.0 , WithinRel( interpolation::loglog( 2.0, 2.5, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 4.0 , WithinRel( interpolation::loglog( 1.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 6.0 , WithinRel( interpolation::loglog( 1.5, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + CHECK_THAT( 8.0 , WithinRel( interpolation::loglog( 2.0, 4.0, xLeft, xRight, fLeft, fRight ) ) ); + } // THEN + } // WHEN } // GIVEN } // SCENARIO diff --git a/src/scion/interpolation/Table.hpp b/src/scion/interpolation/Table.hpp deleted file mode 100644 index 9b83996..0000000 --- a/src/scion/interpolation/Table.hpp +++ /dev/null @@ -1,73 +0,0 @@ -#ifndef NJOY_SCION_INTERPOLATION_TABLE -#define NJOY_SCION_INTERPOLATION_TABLE - -// system includes -#include -#include - -// other includes -#include "tools/Log.hpp" -#include "scion/verification/ranges.hpp" - -namespace njoy { -namespace scion { -namespace interpolation { - - /** - * @class - * @brief Generic table for tabulated data - */ - template < typename Interpolator, typename XContainer, typename YContainer > - class Table { - - /* type aliases */ - using X = std::decay_t< decltype( *( std::declval< XContainer >().begin() ) ) >; - using Y = std::decay_t< decltype( *( std::declval< YContainer >().begin() ) ) >; - - /* fields */ - Interpolator interpolator_; - XContainer x_; - YContainer y_; - - /* auxiliary function */ - #include "scion/interpolation/Table/src/verifyTable.hpp" - - public: - - /* constructor */ - #include "scion/interpolation/Table/src/ctor.hpp" - - /* methods */ - - /** - * @brief Return the panel interpolator - */ - const Interpolator& interpolator() const noexcept { - - return this->interpolator_; - } - - /** - * @brief Return the x values of the table - */ - const XContainer& x() const noexcept { - - return this->x_; - } - - /** - * @brief Return the y values of the table - */ - const YContainer& y() const noexcept { - - return this->y_; - } - - #include "scion/interpolation/Table/src/evaluate.hpp" - }; - -} // math namespace -} // scion namespace -} // njoy namespace - -#endif diff --git a/src/scion/interpolation/Table/src/ctor.hpp b/src/scion/interpolation/Table/src/ctor.hpp deleted file mode 100644 index 64f272a..0000000 --- a/src/scion/interpolation/Table/src/ctor.hpp +++ /dev/null @@ -1,11 +0,0 @@ -/** - * @brief Constructor - * - * @param x the x values of the tabulated data - * @param y the y values of the tabulated data - */ -Table( XContainer x, YContainer y ) : - interpolator_(), x_( std::move( x ) ), y_( std::move( y ) ) { - - verifyTable( this->x(), this->y() ); -} diff --git a/src/scion/interpolation/Table/src/verifyTable.hpp b/src/scion/interpolation/Table/src/verifyTable.hpp deleted file mode 100644 index 36bae26..0000000 --- a/src/scion/interpolation/Table/src/verifyTable.hpp +++ /dev/null @@ -1,27 +0,0 @@ -static void verifyTable( const XContainer& x, - const YContainer& y ) { - - if ( ( ! verification::isAtLeastOfSize( x, 2 ) ) || - ( ! verification::isAtLeastOfSize( y, 2 ) ) ) { - - Log::error( "Insufficient x or y values defined for tabulated data " - "(at least 2 points are required)" ); - Log::info( "x.size(): {}", x.size() ); - Log::info( "y.size(): {}", y.size() ); - throw std::exception(); - } - - if ( ! verification::isSameSize( x, y ) ) { - - Log::error( "Inconsistent number of x and y values for tabulated data" ); - Log::info( "x.size(): {}", x.size() ); - Log::info( "y.size(): {}", y.size() ); - throw std::exception(); - } - - if ( ! verification::isSorted( x ) ) { - - Log::error( "The x values do not appear to be in ascending order" ); - throw std::exception(); - } -} diff --git a/src/scion/interpolation/Table/test/CMakeLists.txt b/src/scion/interpolation/Table/test/CMakeLists.txt deleted file mode 100644 index ecd5dbf..0000000 --- a/src/scion/interpolation/Table/test/CMakeLists.txt +++ /dev/null @@ -1 +0,0 @@ -add_cpp_test( interpolation.Table Table.test.cpp ) diff --git a/src/scion/interpolation/Table/test/Table.test.cpp b/src/scion/interpolation/Table/test/Table.test.cpp deleted file mode 100644 index 461ba2e..0000000 --- a/src/scion/interpolation/Table/test/Table.test.cpp +++ /dev/null @@ -1,102 +0,0 @@ -// include Catch2 -#include -#include -using Catch::Matchers::WithinRel; - -// what we are testing -#include "scion/interpolation/Table.hpp" - -// other includes -#include -#include "scion/interpolation/LinearLinear.hpp" - -// convenience typedefs -using namespace njoy::scion; - -SCENARIO( "Table" ) { - - GIVEN( "valid data for a Table object" ) { - - WHEN( "the data is given explicitly" ) { - - std::vector< double > x = { 1., 2., 3., 4. }; - std::vector< double > y = { 4., 3., 2., 1. }; - - interpolation::Table< interpolation::LinearLinear, - std::vector< double >, - std::vector< double > > chunk( std::move( x ), std::move( y ) ); - - THEN( "a Table can be constructed and members can be tested" ) { - - CHECK( 4 == chunk.x().size() ); - CHECK( 4 == chunk.y().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_THAT( 4., WithinRel( chunk.y()[0] ) ); - CHECK_THAT( 3., WithinRel( chunk.y()[1] ) ); - CHECK_THAT( 2., WithinRel( chunk.y()[2] ) ); - CHECK_THAT( 1., WithinRel( chunk.y()[3] ) ); - } // THEN - - THEN( "a Table can be evaluated" ) { - - // values of x in the x grid - CHECK_THAT( 4., WithinRel( chunk.evaluate( 1. ) ) ); - CHECK_THAT( 3., WithinRel( chunk.evaluate( 2. ) ) ); - CHECK_THAT( 2., WithinRel( chunk.evaluate( 3. ) ) ); - CHECK_THAT( 1., WithinRel( chunk.evaluate( 4. ) ) ); - - // values of x outside the x grid - CHECK_THAT( 0., WithinRel( chunk.evaluate( 0. ) ) ); - CHECK_THAT( 0., WithinRel( chunk.evaluate( 5. ) ) ); - - // values of x inside the x grid - CHECK_THAT( 3.5, WithinRel( chunk.evaluate( 1.5 ) ) ); - CHECK_THAT( 1.5, WithinRel( chunk.evaluate( 3.5 ) ) ); - } // THEN - } // WHEN - } // GIVEN - - GIVEN( "invalid data for a Table object" ) { - - using LinLinTable = interpolation::Table< interpolation::LinearLinear, - std::vector< double >, - std::vector< double > >; - - WHEN( "there are not enough values in the x or y grid" ) { - - std::vector< double > empty = {}; - std::vector< double > one = { 1. }; - - THEN( "an exception is thrown" ) { - - CHECK_THROWS( LinLinTable( empty, empty ) ); - CHECK_THROWS( LinLinTable( one, one ) ); - } // 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< double > y = { 4., 3., 2. }; - - THEN( "an exception is thrown" ) { - - CHECK_THROWS( LinLinTable( std::move( x ), std::move( y ) ) ); - } // THEN - } // WHEN - - WHEN( "the x grid is not sorted" ) { - - std::vector< double > x = { 1., 3., 2., 4. }; - std::vector< double > y = { 4., 3., 2., 1. }; - - THEN( "an exception is thrown" ) { - - CHECK_THROWS( LinLinTable( std::move( x ), std::move( y ) ) ); - } // THEN - } // WHEN - } // GIVEN -} // SCENARIO diff --git a/src/scion/math/ChebyshevApproximation.hpp b/src/scion/math/ChebyshevApproximation.hpp index 5f992ea..c04528f 100644 --- a/src/scion/math/ChebyshevApproximation.hpp +++ b/src/scion/math/ChebyshevApproximation.hpp @@ -5,7 +5,7 @@ #include // other includes -#include "scion/math/FunctionBase.hpp" +#include "scion/math/OneDimensionalFunctionBase.hpp" #include "scion/math/ChebyshevSeries.hpp" namespace njoy { @@ -40,10 +40,11 @@ namespace math { * equals 0. */ template < typename X, typename Y = X > - class ChebyshevApproximation : public FunctionBase< ChebyshevApproximation< X, Y >, X, Y > { + class ChebyshevApproximation : + public OneDimensionalFunctionBase< ChebyshevApproximation< X, Y >, X, Y > { /* type aliases */ - using Parent = FunctionBase< ChebyshevApproximation< X, Y >, X, Y >; + using Parent = OneDimensionalFunctionBase< ChebyshevApproximation< X, Y >, X, Y >; /* fields */ X lower_; diff --git a/src/scion/math/ChebyshevSeries.hpp b/src/scion/math/ChebyshevSeries.hpp index 7de002a..d82d6b6 100644 --- a/src/scion/math/ChebyshevSeries.hpp +++ b/src/scion/math/ChebyshevSeries.hpp @@ -6,6 +6,7 @@ // other includes #include "scion/math/clenshaw.hpp" #include "scion/math/matrix.hpp" +#include "scion/math/OneDimensionalFunctionBase.hpp" #include "scion/math/SeriesBase.hpp" namespace njoy { @@ -48,7 +49,7 @@ namespace math { /* friend declarations */ friend class SeriesBase< ChebyshevSeries< X, Y >, X, Y >; - friend class FunctionBase< ChebyshevSeries< X, Y >, X, Y >; + friend class OneDimensionalFunctionBase< ChebyshevSeries< X, Y >, X, Y >; /* type aliases */ using Parent = SeriesBase< ChebyshevSeries< X, Y >, X, Y >; diff --git a/src/scion/math/FunctionBase/src/ctor.hpp b/src/scion/math/FunctionBase/src/ctor.hpp deleted file mode 100644 index c9be75a..0000000 --- a/src/scion/math/FunctionBase/src/ctor.hpp +++ /dev/null @@ -1,7 +0,0 @@ -/** - * @brief Constructor - * - * @param domain the domain of the function - */ -FunctionBase( DomainVariant domain ) : - domain_( std::move( domain ) ) {} diff --git a/src/scion/math/HistogramTable.hpp b/src/scion/math/HistogramTable.hpp index f8a5133..0ded725 100644 --- a/src/scion/math/HistogramTable.hpp +++ b/src/scion/math/HistogramTable.hpp @@ -16,13 +16,13 @@ namespace math { /** * @class - * @brief Tabulated data with histogram interpolation (y is constant in x) + * @brief Tabulated x,y data with histogram interpolation (y is constant in x) * * The HistogramTable is templatised on the container type used for the * x and y values in addition to the actual x and y types. This allows us to * use something like utility::IteratorView instead of std::vector. */ - template < typename X, typename Y = X, + template < typename X, typename Y, typename XContainer = std::vector< X >, typename YContainer = std::vector< Y > > class HistogramTable : diff --git a/src/scion/math/HistogramTableFunction.hpp b/src/scion/math/HistogramTableFunction.hpp new file mode 100644 index 0000000..263e18e --- /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 histogram interpolation (f(y) is constant 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/InterpolationTable.hpp b/src/scion/math/InterpolationTable.hpp index baaadd8..265df99 100644 --- a/src/scion/math/InterpolationTable.hpp +++ b/src/scion/math/InterpolationTable.hpp @@ -11,7 +11,7 @@ #include "scion/linearisation/ToleranceConvergence.hpp" #include "scion/unionisation/unionise.hpp" #include "scion/math/newton.hpp" -#include "scion/math/FunctionBase.hpp" +#include "scion/math/OneDimensionalFunctionBase.hpp" #include "scion/math/HistogramTable.hpp" #include "scion/math/LinearLinearTable.hpp" #include "scion/math/LinearLogTable.hpp" @@ -29,13 +29,14 @@ namespace math { * @brief Tabulated data with one or more interpolation types */ template < typename X, typename Y = X > - class InterpolationTable : public FunctionBase< InterpolationTable< X, Y >, X, Y > { + class InterpolationTable : + public OneDimensionalFunctionBase< InterpolationTable< X, Y >, X, Y > { /* friend declarations */ - friend class FunctionBase< InterpolationTable< X, Y >, X, Y >; + friend class OneDimensionalFunctionBase< InterpolationTable< X, Y >, X, Y >; /* type aliases */ - using Parent = FunctionBase< InterpolationTable< X, Y >, X, Y >; + using Parent = OneDimensionalFunctionBase< InterpolationTable< X, Y >, X, Y >; using XIterator = typename std::vector< X >::const_iterator; using YIterator = typename std::vector< Y >::const_iterator; using XContainer = njoy::utility::IteratorView< XIterator >; diff --git a/src/scion/math/LegendreSeries.hpp b/src/scion/math/LegendreSeries.hpp index c6bd8a8..481e9c8 100644 --- a/src/scion/math/LegendreSeries.hpp +++ b/src/scion/math/LegendreSeries.hpp @@ -7,6 +7,7 @@ // other includes #include "scion/math/clenshaw.hpp" #include "scion/math/matrix.hpp" +#include "scion/math/OneDimensionalFunctionBase.hpp" #include "scion/math/SeriesBase.hpp" namespace njoy { @@ -54,7 +55,7 @@ namespace math { /* friend declarations */ friend class SeriesBase< LegendreSeries< X, Y >, X, Y >; - friend class FunctionBase< LegendreSeries< X, Y >, X, Y >; + friend class OneDimensionalFunctionBase< LegendreSeries< X, Y >, X, Y >; /* type aliases */ using Parent = SeriesBase< LegendreSeries< X, Y >, X, Y >; diff --git a/src/scion/math/LinearLinearTable.hpp b/src/scion/math/LinearLinearTable.hpp index f593ed4..4cd54ec 100644 --- a/src/scion/math/LinearLinearTable.hpp +++ b/src/scion/math/LinearLinearTable.hpp @@ -16,13 +16,13 @@ namespace math { /** * @class - * @brief Tabulated data with linear-linear interpolation (y is linear in x) + * @brief Tabulated x,y data with linear-linear interpolation (y is linear in x) * * The LinearLinearTable is templatised on the container type used for the * x and y values in addition to the actual x and y types. This allows us to * use something like utility::IteratorView instead of std::vector. */ - template < typename X, typename Y = X, + template < typename X, typename Y, typename XContainer = std::vector< X >, typename YContainer = std::vector< Y > > class LinearLinearTable : diff --git a/src/scion/math/LinearLinearTableFunction.hpp b/src/scion/math/LinearLinearTableFunction.hpp new file mode 100644 index 0000000..e98fb19 --- /dev/null +++ b/src/scion/math/LinearLinearTableFunction.hpp @@ -0,0 +1,75 @@ +#ifndef NJOY_SCION_MATH_LINEARLINEARTABLEFUNCTION +#define NJOY_SCION_MATH_LINEARLINEARTABLEFUNCTION + +// system includes +#include + +// other includes +#include "scion/interpolation/InterpolationType.hpp" +#include "scion/interpolation/LinearLinear.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 LinearLinearTableFunction 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 LinearLinearTableFunction : + public SingleTableFunctionBase< LinearLinearTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LinearLinear, X, Y, Z, F, + XContainer, FContainer > { + + /* friend declarations */ + friend class SingleTableFunctionBase< LinearLinearTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LinearLinear, X, Y, Z, F, + XContainer, FContainer >; + + /* type aliases */ + using Parent = SingleTableFunctionBase< LinearLinearTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LinearLinear, 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::LinearLinear; + } + + public: + + /* constructor */ + #include "scion/math/LinearLinearTableFunction/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/LinearLinearTableFunction/src/ctor.hpp b/src/scion/math/LinearLinearTableFunction/src/ctor.hpp new file mode 100644 index 0000000..2f5f532 --- /dev/null +++ b/src/scion/math/LinearLinearTableFunction/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 + */ +LinearLinearTableFunction( XContainer x, FContainer f ) : + Parent( std::move( x ), std::move( f ) ) {} diff --git a/src/scion/math/LinearLinearTableFunction/test/CMakeLists.txt b/src/scion/math/LinearLinearTableFunction/test/CMakeLists.txt new file mode 100644 index 0000000..c872af3 --- /dev/null +++ b/src/scion/math/LinearLinearTableFunction/test/CMakeLists.txt @@ -0,0 +1 @@ +add_cpp_test( math.LinearLinearTableFunction LinearLinearTableFunction.test.cpp ) diff --git a/src/scion/math/LinearLinearTableFunction/test/LinearLinearTableFunction.test.cpp b/src/scion/math/LinearLinearTableFunction/test/LinearLinearTableFunction.test.cpp new file mode 100644 index 0000000..2890ebb --- /dev/null +++ b/src/scion/math/LinearLinearTableFunction/test/LinearLinearTableFunction.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/LinearLinearTableFunction.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 LinearLinearTableFunction = math::LinearLinearTableFunction< X, Y, Z, F, XContainer, FContainer >; +using InterpolationType = interpolation::InterpolationType; + +SCENARIO( "LinearLinearTableFunction" ) { + + 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 } } + }; + + LinearLinearTableFunction< double > chunk( std::move( x ), std::move( f ) ); + + THEN( "a LinearLinearTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::LinearLinear == 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 LinearLinearTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.3 , 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.4975, WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.4725, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.375 , 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 ); + + LinearLinearTableFunction< double, double, double, InterpolationTable< double >, + XView, FView > chunk( std::move( x ), std::move( f ) ); + + THEN( "a LinearLinearTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::LinearLinear == 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 LinearLinearTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.3 , 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.4975, WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.4725, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.375 , WithinRel( chunk( 3.5, -0.5 ) ) ); + } // THEN + } // WHEN + } // GIVEN + + GIVEN( "invalid data for a LinearLinearTable 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( LinearLinearTableFunction< double >( xempty, fempty ) ); + CHECK_THROWS( LinearLinearTableFunction< double >( xone, fone ) ); + CHECK_THROWS( LinearLinearTableFunction< double >( xempty, fone ) ); + CHECK_THROWS( LinearLinearTableFunction< 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( LinearLinearTableFunction< 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( LinearLinearTableFunction< 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( LinearLinearTableFunction< double >( std::move( x ), std::move( f ) ) ); + } // THEN + } // WHEN + } // GIVEN +} // SCENARIO diff --git a/src/scion/math/LinearLogTable.hpp b/src/scion/math/LinearLogTable.hpp index c19fbd0..8ad6fe0 100644 --- a/src/scion/math/LinearLogTable.hpp +++ b/src/scion/math/LinearLogTable.hpp @@ -15,13 +15,13 @@ namespace math { /** * @class - * @brief Tabulated data with linear-log interpolation (y is linear in ln(x)) + * @brief Tabulated x,y data with linear-log interpolation (y is linear in ln(x)) * * The LinearLogTable is templatised on the container type used for the * x and y values in addition to the actual x and y types. This allows us to * use something like utility::IteratorView instead of std::vector. */ - template < typename X, typename Y = X, + template < typename X, typename Y, typename XContainer = std::vector< X >, typename YContainer = std::vector< Y > > class LinearLogTable : diff --git a/src/scion/math/LinearLogTableFunction.hpp b/src/scion/math/LinearLogTableFunction.hpp new file mode 100644 index 0000000..1da89ea --- /dev/null +++ b/src/scion/math/LinearLogTableFunction.hpp @@ -0,0 +1,75 @@ +#ifndef NJOY_SCION_MATH_LINEARLOGTABLEFUNCTION +#define NJOY_SCION_MATH_LINEARLOGTABLEFUNCTION + +// system includes +#include + +// other includes +#include "scion/interpolation/InterpolationType.hpp" +#include "scion/interpolation/LinearLogarithmic.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-log interpolation (f(y) is linear in ln(x)) + * + * The LinearLogTableFunction 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 LinearLogTableFunction : + public SingleTableFunctionBase< LinearLogTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LinearLogarithmic, X, Y, Z, F, + XContainer, FContainer > { + + /* friend declarations */ + friend class SingleTableFunctionBase< LinearLogTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LinearLogarithmic, X, Y, Z, F, + XContainer, FContainer >; + + /* type aliases */ + using Parent = SingleTableFunctionBase< LinearLogTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LinearLogarithmic, 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::LinearLog; + } + + public: + + /* constructor */ + #include "scion/math/LinearLogTableFunction/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/LinearLogTableFunction/src/ctor.hpp b/src/scion/math/LinearLogTableFunction/src/ctor.hpp new file mode 100644 index 0000000..9886030 --- /dev/null +++ b/src/scion/math/LinearLogTableFunction/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 + */ +LinearLogTableFunction( XContainer x, FContainer f ) : + Parent( std::move( x ), std::move( f ) ) {} diff --git a/src/scion/math/LinearLogTableFunction/test/CMakeLists.txt b/src/scion/math/LinearLogTableFunction/test/CMakeLists.txt new file mode 100644 index 0000000..70289a1 --- /dev/null +++ b/src/scion/math/LinearLogTableFunction/test/CMakeLists.txt @@ -0,0 +1 @@ +add_cpp_test( math.LinearLogTableFunction LinearLogTableFunction.test.cpp ) diff --git a/src/scion/math/LinearLogTableFunction/test/LinearLogTableFunction.test.cpp b/src/scion/math/LinearLogTableFunction/test/LinearLogTableFunction.test.cpp new file mode 100644 index 0000000..ae687a0 --- /dev/null +++ b/src/scion/math/LinearLogTableFunction/test/LinearLogTableFunction.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/LinearLogTableFunction.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 LinearLogTableFunction = math::LinearLogTableFunction< X, Y, Z, F, XContainer, FContainer >; +using InterpolationType = interpolation::InterpolationType; + +SCENARIO( "LinearLogTableFunction" ) { + + 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 } } + }; + + LinearLogTableFunction< double > chunk( std::move( x ), std::move( f ) ); + + THEN( "a LinearLogTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::LinearLog == 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 LinearLogTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.3 , 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.49707518749639, WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.47023471290541, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.36962445981765, 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 ); + + LinearLogTableFunction< double, double, double, InterpolationTable< double >, + XView, FView > chunk( std::move( x ), std::move( f ) ); + + THEN( "a LinearLogTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::LinearLog == 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 LinearLogTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.3 , 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.49707518749639, WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.47023471290541, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.36962445981765, WithinRel( chunk( 3.5, -0.5 ) ) ); + } // THEN + } // WHEN + } // GIVEN + + GIVEN( "invalid data for a LinearLogTable 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( LinearLogTableFunction< double >( xempty, fempty ) ); + CHECK_THROWS( LinearLogTableFunction< double >( xone, fone ) ); + CHECK_THROWS( LinearLogTableFunction< double >( xempty, fone ) ); + CHECK_THROWS( LinearLogTableFunction< 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( LinearLogTableFunction< 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( LinearLogTableFunction< 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( LinearLogTableFunction< double >( std::move( x ), std::move( f ) ) ); + } // THEN + } // WHEN + } // GIVEN +} // SCENARIO diff --git a/src/scion/math/LogLinearTable.hpp b/src/scion/math/LogLinearTable.hpp index 11ed320..8ce910c 100644 --- a/src/scion/math/LogLinearTable.hpp +++ b/src/scion/math/LogLinearTable.hpp @@ -15,13 +15,13 @@ namespace math { /** * @class - * @brief Tabulated data with log-linear interpolation (ln(y) is linear in x) + * @brief Tabulated x,y data with log-linear interpolation (ln(y) is linear in x) * * The LogLinearTable is templatised on the container type used for the * x and y values in addition to the actual x and y types. This allows us to * use something like utility::IteratorView instead of std::vector. */ - template < typename X, typename Y = X, + template < typename X, typename Y, typename XContainer = std::vector< X >, typename YContainer = std::vector< Y > > class LogLinearTable : diff --git a/src/scion/math/LogLinearTableFunction.hpp b/src/scion/math/LogLinearTableFunction.hpp new file mode 100644 index 0000000..ba07858 --- /dev/null +++ b/src/scion/math/LogLinearTableFunction.hpp @@ -0,0 +1,75 @@ +#ifndef NJOY_SCION_MATH_LOGLINEARTABLEFUNCTION +#define NJOY_SCION_MATH_LOGLINEARTABLEFUNCTION + +// system includes +#include + +// other includes +#include "scion/interpolation/InterpolationType.hpp" +#include "scion/interpolation/LogarithmicLinear.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 log-linear interpolation (ln(f(y)) is linear in x) + * + * The LogLinearTableFunction 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 LogLinearTableFunction : + public SingleTableFunctionBase< LogLinearTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LogarithmicLinear, X, Y, Z, F, + XContainer, FContainer > { + + /* friend declarations */ + friend class SingleTableFunctionBase< LogLinearTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LogarithmicLinear, X, Y, Z, F, + XContainer, FContainer >; + + /* type aliases */ + using Parent = SingleTableFunctionBase< LogLinearTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LogarithmicLinear, 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::LogLinear; + } + + public: + + /* constructor */ + #include "scion/math/LogLinearTableFunction/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/LogLinearTableFunction/src/ctor.hpp b/src/scion/math/LogLinearTableFunction/src/ctor.hpp new file mode 100644 index 0000000..24fb4a1 --- /dev/null +++ b/src/scion/math/LogLinearTableFunction/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 + */ +LogLinearTableFunction( XContainer x, FContainer f ) : + Parent( std::move( x ), std::move( f ) ) {} diff --git a/src/scion/math/LogLinearTableFunction/test/CMakeLists.txt b/src/scion/math/LogLinearTableFunction/test/CMakeLists.txt new file mode 100644 index 0000000..eb93008 --- /dev/null +++ b/src/scion/math/LogLinearTableFunction/test/CMakeLists.txt @@ -0,0 +1 @@ +add_cpp_test( math.LogLinearTableFunction LogLinearTableFunction.test.cpp ) diff --git a/src/scion/math/LogLinearTableFunction/test/LogLinearTableFunction.test.cpp b/src/scion/math/LogLinearTableFunction/test/LogLinearTableFunction.test.cpp new file mode 100644 index 0000000..6484d8a --- /dev/null +++ b/src/scion/math/LogLinearTableFunction/test/LogLinearTableFunction.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/LogLinearTableFunction.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 LogLinearTableFunction = math::LogLinearTableFunction< X, Y, Z, F, XContainer, FContainer >; +using InterpolationType = interpolation::InterpolationType; + +SCENARIO( "LogLinearTableFunction" ) { + + 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 } } + }; + + LogLinearTableFunction< double > chunk( std::move( x ), std::move( f ) ); + + THEN( "a LogLinearTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::LogLinear == 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 LogLinearTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.3 , 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.49749371855331, WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.47196398167657, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.36742346141748, 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 ); + + LogLinearTableFunction< double, double, double, InterpolationTable< double >, + XView, FView > chunk( std::move( x ), std::move( f ) ); + + THEN( "a LogLinearTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::LogLinear == 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 LogLinearTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.3 , 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.49749371855331, WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.47196398167657, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.36742346141748, WithinRel( chunk( 3.5, -0.5 ) ) ); + } // THEN + } // WHEN + } // GIVEN + + GIVEN( "invalid data for a LogLinearTable 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( LogLinearTableFunction< double >( xempty, fempty ) ); + CHECK_THROWS( LogLinearTableFunction< double >( xone, fone ) ); + CHECK_THROWS( LogLinearTableFunction< double >( xempty, fone ) ); + CHECK_THROWS( LogLinearTableFunction< 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( LogLinearTableFunction< 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( LogLinearTableFunction< 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( LogLinearTableFunction< double >( std::move( x ), std::move( f ) ) ); + } // THEN + } // WHEN + } // GIVEN +} // SCENARIO diff --git a/src/scion/math/LogLogTable.hpp b/src/scion/math/LogLogTable.hpp index 294bddf..c5fdf67 100644 --- a/src/scion/math/LogLogTable.hpp +++ b/src/scion/math/LogLogTable.hpp @@ -15,13 +15,13 @@ namespace math { /** * @class - * @brief Tabulated data with log-log interpolation (ln(y) is linear in ln(x)) + * @brief Tabulated x,y data with log-log interpolation (ln(y) is linear in ln(x)) * * The LogLogTable is templatised on the container type used for the * x and y values in addition to the actual x and y types. This allows us to * use something like utility::IteratorView instead of std::vector. */ - template < typename X, typename Y = X, + template < typename X, typename Y, typename XContainer = std::vector< X >, typename YContainer = std::vector< Y > > class LogLogTable : diff --git a/src/scion/math/LogLogTableFunction.hpp b/src/scion/math/LogLogTableFunction.hpp new file mode 100644 index 0000000..37e47b8 --- /dev/null +++ b/src/scion/math/LogLogTableFunction.hpp @@ -0,0 +1,75 @@ +#ifndef NJOY_SCION_MATH_LOGLOGTABLEFUNCTION +#define NJOY_SCION_MATH_LOGLOGTABLEFUNCTION + +// system includes +#include + +// other includes +#include "scion/interpolation/InterpolationType.hpp" +#include "scion/interpolation/LogarithmicLogarithmic.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 log-log interpolation (ln(f(y)) is linear in ln(x)) + * + * The LogLogTableFunction 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 LogLogTableFunction : + public SingleTableFunctionBase< LogLogTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LogarithmicLogarithmic, X, Y, Z, F, + XContainer, FContainer > { + + /* friend declarations */ + friend class SingleTableFunctionBase< LogLogTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LogarithmicLogarithmic, X, Y, Z, F, + XContainer, FContainer >; + + /* type aliases */ + using Parent = SingleTableFunctionBase< LogLogTableFunction< X, Y, Z, F, XContainer, FContainer >, + interpolation::LogarithmicLogarithmic, 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::LogLog; + } + + public: + + /* constructor */ + #include "scion/math/LogLogTableFunction/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/LogLogTableFunction/src/ctor.hpp b/src/scion/math/LogLogTableFunction/src/ctor.hpp new file mode 100644 index 0000000..ba02abb --- /dev/null +++ b/src/scion/math/LogLogTableFunction/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 + */ +LogLogTableFunction( XContainer x, FContainer f ) : + Parent( std::move( x ), std::move( f ) ) {} diff --git a/src/scion/math/LogLogTableFunction/test/CMakeLists.txt b/src/scion/math/LogLogTableFunction/test/CMakeLists.txt new file mode 100644 index 0000000..c05080a --- /dev/null +++ b/src/scion/math/LogLogTableFunction/test/CMakeLists.txt @@ -0,0 +1 @@ +add_cpp_test( math.LogLogTableFunction LogLogTableFunction.test.cpp ) diff --git a/src/scion/math/LogLogTableFunction/test/LogLogTableFunction.test.cpp b/src/scion/math/LogLogTableFunction/test/LogLogTableFunction.test.cpp new file mode 100644 index 0000000..e3f30ae --- /dev/null +++ b/src/scion/math/LogLogTableFunction/test/LogLogTableFunction.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/LogLogTableFunction.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 LogLogTableFunction = math::LogLogTableFunction< X, Y, Z, F, XContainer, FContainer >; +using InterpolationType = interpolation::InterpolationType; + +SCENARIO( "LogLogTableFunction" ) { + + 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 } } + }; + + LogLogTableFunction< double > chunk( std::move( x ), std::move( f ) ); + + THEN( "a LogLogTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::LogLog == 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 LogLogTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.3 , 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.49706908915929, WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.46970497533108, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.36212316985370, 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 ); + + LogLogTableFunction< double, double, double, InterpolationTable< double >, + XView, FView > chunk( std::move( x ), std::move( f ) ); + + THEN( "a LogLogTable can be constructed and members can be tested" ) { + + CHECK( InterpolationType::LogLog == 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 LogLogTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 0.5 , WithinRel( chunk( 1., -0.5 ) ) ); + CHECK_THAT( 0.495, WithinRel( chunk( 2., -0.5 ) ) ); + CHECK_THAT( 0.45 , WithinRel( chunk( 3., -0.5 ) ) ); + CHECK_THAT( 0.3 , 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.49706908915929, WithinRel( chunk( 1.5, -0.5 ) ) ); + CHECK_THAT( 0.46970497533108, WithinRel( chunk( 2.5, -0.5 ) ) ); + CHECK_THAT( 0.36212316985370, WithinRel( chunk( 3.5, -0.5 ) ) ); + } // THEN + } // WHEN + } // GIVEN + + GIVEN( "invalid data for a LogLogTable 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( LogLogTableFunction< double >( xempty, fempty ) ); + CHECK_THROWS( LogLogTableFunction< double >( xone, fone ) ); + CHECK_THROWS( LogLogTableFunction< double >( xempty, fone ) ); + CHECK_THROWS( LogLogTableFunction< 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( LogLogTableFunction< 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( LogLogTableFunction< 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( LogLogTableFunction< double >( std::move( x ), std::move( f ) ) ); + } // THEN + } // WHEN + } // GIVEN +} // SCENARIO diff --git a/src/scion/math/FunctionBase.hpp b/src/scion/math/OneDimensionalFunctionBase.hpp similarity index 77% rename from src/scion/math/FunctionBase.hpp rename to src/scion/math/OneDimensionalFunctionBase.hpp index 2f289d0..ea8e0c3 100644 --- a/src/scion/math/FunctionBase.hpp +++ b/src/scion/math/OneDimensionalFunctionBase.hpp @@ -1,5 +1,5 @@ -#ifndef NJOY_SCION_MATH_FUNCTIONBASE -#define NJOY_SCION_MATH_FUNCTIONBASE +#ifndef NJOY_SCION_MATH_ONEDIMENSIONALFUNCTIONBASE +#define NJOY_SCION_MATH_ONEDIMENSIONALFUNCTIONBASE // system includes #include @@ -14,13 +14,13 @@ namespace math { /** * @class - * @brief Base class for function objects + * @brief Base class for one dimensional function objects modelling y = f(x) * - * This base class provides the common interface for all function objects. - * This includes domain testing and function evaluation. + * This base class provides the common interface for all one dimensional + * function objects. This includes domain testing and function evaluation. */ - template < typename Derived, typename X, typename Y = X > - class FunctionBase { + template < typename Derived, typename X, typename Y > + class OneDimensionalFunctionBase { public: @@ -37,7 +37,14 @@ namespace math { protected: /* constructor */ - #include "scion/math/FunctionBase/src/ctor.hpp" + + /** + * @brief Constructor + * + * @param domain the domain of the function + */ + OneDimensionalFunctionBase( DomainVariant domain ) : + domain_( std::move( domain ) ) {} public: diff --git a/src/scion/math/PolynomialSeries.hpp b/src/scion/math/PolynomialSeries.hpp index e0f5b78..60f02fc 100644 --- a/src/scion/math/PolynomialSeries.hpp +++ b/src/scion/math/PolynomialSeries.hpp @@ -7,6 +7,7 @@ // other includes #include "scion/math/horner.hpp" #include "scion/math/matrix.hpp" +#include "scion/math/OneDimensionalFunctionBase.hpp" #include "scion/math/SeriesBase.hpp" namespace njoy { @@ -40,7 +41,7 @@ namespace math { /* friend declarations */ friend class SeriesBase< PolynomialSeries< X, Y >, X, Y >; - friend class FunctionBase< PolynomialSeries< X, Y >, X, Y >; + friend class OneDimensionalFunctionBase< PolynomialSeries< X, Y >, X, Y >; /* type aliases */ using Parent = SeriesBase< PolynomialSeries< X, Y >, X, Y >; diff --git a/src/scion/math/SeriesBase.hpp b/src/scion/math/SeriesBase.hpp index f273106..a19fa35 100644 --- a/src/scion/math/SeriesBase.hpp +++ b/src/scion/math/SeriesBase.hpp @@ -10,7 +10,7 @@ #include "scion/linearisation/MidpointSplit.hpp" #include "scion/linearisation/Lineariser.hpp" #include "scion/math/InterpolationTable.hpp" -#include "scion/math/FunctionBase.hpp" +#include "scion/math/OneDimensionalFunctionBase.hpp" #include "scion/math/compare.hpp" #include "scion/verification/ranges.hpp" @@ -26,10 +26,10 @@ namespace math { * as the polynomial series, Legendre series and Chebyshev series. */ template < typename Derived, typename X, typename Y = X > - class SeriesBase : public FunctionBase< Derived, X, Y > { + class SeriesBase : public OneDimensionalFunctionBase< Derived, X, Y > { /* type aliases */ - using Parent = FunctionBase< Derived, X, Y >; + using Parent = OneDimensionalFunctionBase< Derived, X, Y >; public: diff --git a/src/scion/math/SingleTableBase.hpp b/src/scion/math/SingleTableBase.hpp index 69aa85c..8ef0989 100644 --- a/src/scion/math/SingleTableBase.hpp +++ b/src/scion/math/SingleTableBase.hpp @@ -6,11 +6,10 @@ // other includes #include "tools/Log.hpp" #include "scion/interpolation/InterpolationType.hpp" -#include "scion/interpolation/Table.hpp" #include "scion/linearisation/ToleranceConvergence.hpp" #include "scion/linearisation/MidpointSplit.hpp" #include "scion/linearisation/Lineariser.hpp" -#include "scion/math/FunctionBase.hpp" +#include "scion/math/OneDimensionalFunctionBase.hpp" #include "scion/verification/ranges.hpp" namespace njoy { @@ -19,21 +18,19 @@ namespace math { /** * @class - * @brief Base class for tabulated data using a single interpolation region + * @brief Base class for x,y tabulated data using a single interpolation region * * This base class provides the common interface for single region * interpolation data such as the LinearLinearTable, LogLogTable, etc. */ template < typename Derived, typename Interpolator, - typename X, typename Y = X, + typename X, typename Y, typename XContainer = std::vector< X >, typename YContainer = std::vector< Y > > - class SingleTableBase : public FunctionBase< Derived, X, Y > { + class SingleTableBase : public OneDimensionalFunctionBase< Derived, X, Y > { /* type aliases */ - using Parent = FunctionBase< Derived, X, Y >; - using Table = interpolation::Table< Interpolator, - XContainer, YContainer >; + using Parent = OneDimensionalFunctionBase< Derived, X, Y >; public: @@ -43,11 +40,21 @@ namespace math { private: /* fields */ - Table table_; + Interpolator interpolator_; + XContainer x_; + YContainer y_; /* auxiliary function */ #include "scion/math/SingleTableBase/src/verifyTable.hpp" + /** + * @brief Return the panel interpolator + */ + const Interpolator& interpolator() const noexcept { + + return this->interpolator_; + } + protected: /* constructor */ @@ -70,7 +77,7 @@ namespace math { */ const XContainer& x() const noexcept { - return this->table_.x(); + return this->x_; } /** @@ -78,7 +85,7 @@ namespace math { */ const YContainer& y() const noexcept { - return this->table_.y(); + return this->y_; } /** @@ -89,16 +96,7 @@ namespace math { return this->x().size(); } - /** - * @brief Evaluate the function for a value of x - * - * @param x the value to be evaluated - */ - Y evaluate( const X& x ) const { - - return this->table_.evaluate( x ); - } - + #include "scion/math/SingleTableBase/src/evaluate.hpp" #include "scion/math/SingleTableBase/src/linearise.hpp" using Parent::domain; diff --git a/src/scion/math/SingleTableBase/src/ctor.hpp b/src/scion/math/SingleTableBase/src/ctor.hpp index 19af06b..2aa659d 100644 --- a/src/scion/math/SingleTableBase/src/ctor.hpp +++ b/src/scion/math/SingleTableBase/src/ctor.hpp @@ -1,14 +1,11 @@ private: /** - * @brief Private intermediate constructor + * @brief Private constructor */ -SingleTableBase( Table&& table ) : - Parent( IntervalDomain( table.x().front(), table.x().back() ) ), - table_( std::move( table ) ) { - - verifyTable( this->table_.x(), this->table_.y() ); -} +SingleTableBase( IntervalDomain< X >&& domain, XContainer&& x, YContainer&& y ) : + Parent( std::move( domain ) ), interpolator_(), + x_( std::move( x ) ), y_( std::move( y ) ) {} protected: @@ -19,4 +16,6 @@ SingleTableBase( Table&& table ) : * @param y the y values of the tabulated data */ SingleTableBase( XContainer x, YContainer y ) : - SingleTableBase( Table( std::move( x ), std::move( y ) ) ) {} + SingleTableBase( verifyTableAndRetrieveDomain( x, y ), + std::move( x ), std::move( y ) ) {} + diff --git a/src/scion/interpolation/Table/src/evaluate.hpp b/src/scion/math/SingleTableBase/src/evaluate.hpp similarity index 100% rename from src/scion/interpolation/Table/src/evaluate.hpp rename to src/scion/math/SingleTableBase/src/evaluate.hpp diff --git a/src/scion/math/SingleTableBase/src/verifyTable.hpp b/src/scion/math/SingleTableBase/src/verifyTable.hpp index 0764bb6..aa695ac 100644 --- a/src/scion/math/SingleTableBase/src/verifyTable.hpp +++ b/src/scion/math/SingleTableBase/src/verifyTable.hpp @@ -1,8 +1,36 @@ -static void verifyTable( const XContainer& x, const YContainer& y ) { +static IntervalDomain< X > +verifyTableAndRetrieveDomain( const XContainer& x, const YContainer& y ) { + + if ( ( ! verification::isAtLeastOfSize( x, 2 ) ) || + ( ! verification::isAtLeastOfSize( y, 2 ) ) ) { + + Log::error( "Insufficient x or y values defined for x,y tabulated data " + "with a single interpolation type (at least 2 points are " + "required)" ); + Log::info( "x.size(): {}", x.size() ); + Log::info( "y.size(): {}", y.size() ); + throw std::exception(); + } + + if ( ! verification::isSameSize( x, y ) ) { + + Log::error( "Inconsistent number of x and y values for x,y tabulated data " + "with a single interpolation type" ); + Log::info( "x.size(): {}", x.size() ); + Log::info( "y.size(): {}", y.size() ); + throw std::exception(); + } + + if ( ! verification::isSorted( x ) ) { + + Log::error( "The x values do not appear to be in ascending order for " + "x,y tabulated data with a single interpolation type" ); + throw std::exception(); + } if ( ! verification::isUnique( x ) ) { - Log::error( "The x values do not appear to be unique for tabulated values " + Log::error( "The x values do not appear to be unique for x,y tabulated values " "with a single interpolation type" ); auto iter = std::adjacent_find( x.begin(), x.end() ); @@ -13,4 +41,6 @@ static void verifyTable( const XContainer& x, const YContainer& y ) { } throw std::exception(); } + + return IntervalDomain( x.front(), x.back() ); } diff --git a/src/scion/math/SingleTableFunctionBase.hpp b/src/scion/math/SingleTableFunctionBase.hpp new file mode 100644 index 0000000..7bf687c --- /dev/null +++ b/src/scion/math/SingleTableFunctionBase.hpp @@ -0,0 +1,95 @@ +#ifndef NJOY_SCION_MATH_SINGLETABLEFUNCTIONBASE +#define NJOY_SCION_MATH_SINGLETABLEFUNCTIONBASE + +// system includes + +// other includes +#include "tools/Log.hpp" +#include "scion/interpolation/InterpolationType.hpp" +#include "scion/math/TwoDimensionalFunctionBase.hpp" +#include "scion/verification/ranges.hpp" + +namespace njoy { +namespace scion { +namespace math { + + /** + * @class + * @brief Base class for x,f(y) tabulated data using a single interpolation region + */ + template < typename Derived, typename Interpolator, + typename X, typename Y, typename Z, typename F, + typename XContainer = std::vector< X >, + typename FContainer = std::vector< F > > + class SingleTableFunctionBase : public TwoDimensionalFunctionBase< Derived, X, Y, Z > { + + /* type aliases */ + using Parent = TwoDimensionalFunctionBase< Derived, X, Y, Z >; + + /* fields */ + Interpolator interpolator_; + XContainer x_; + FContainer f_; + + /* auxiliary function */ + #include "scion/math/SingleTableFunctionBase/src/verifyTable.hpp" + + /** + * @brief Return the panel interpolator + */ + const Interpolator& interpolator() const noexcept { + + return this->interpolator_; + } + + protected: + + /* constructor */ + #include "scion/math/SingleTableFunctionBase/src/ctor.hpp" + + public: + + /* methods */ + + /** + * @brief Return the interpolation type + */ + constexpr interpolation::InterpolationType interpolation() const noexcept { + + return static_cast< const Derived* >( this )->type(); + } + + /** + * @brief Return the x values of the table + */ + const XContainer& x() const noexcept { + + return this->x_; + } + + /** + * @brief Return the f(y) functions of the table + */ + const FContainer& f() const noexcept { + + return this->f_; + } + + /** + * @brief Return the number of points in the table + */ + std::size_t numberPoints() const noexcept { + + return this->x().size(); + } + + #include "scion/math/SingleTableFunctionBase/src/evaluate.hpp" + + using Parent::operator(); + }; + +} // math namespace +} // scion namespace +} // njoy namespace + +#endif diff --git a/src/scion/math/SingleTableFunctionBase/src/ctor.hpp b/src/scion/math/SingleTableFunctionBase/src/ctor.hpp new file mode 100644 index 0000000..43e19ac --- /dev/null +++ b/src/scion/math/SingleTableFunctionBase/src/ctor.hpp @@ -0,0 +1,11 @@ +/** + * @brief Constructor + * + * @param x the x values of the tabulated data + * @param f the f(y) functions of the tabulated data + */ +SingleTableFunctionBase( XContainer x, FContainer f ) : + interpolator_(), x_( std::move( x ) ), f_( std::move( f ) ) { + + verifyTable( this->x(), this->f() ); +} diff --git a/src/scion/math/SingleTableFunctionBase/src/evaluate.hpp b/src/scion/math/SingleTableFunctionBase/src/evaluate.hpp new file mode 100644 index 0000000..05df693 --- /dev/null +++ b/src/scion/math/SingleTableFunctionBase/src/evaluate.hpp @@ -0,0 +1,31 @@ +/** + * @brief Evaluate the function for a value of x and y + * + * @param x the x value to be evaluated + * @param y the y value to be evaluated + */ +Z evaluate( const X& x, const Y& y ) const { + + auto xIter = std::next( this->x().begin() ); + if ( x != *this->x().begin() ) { + + xIter = std::lower_bound( this->x().begin(), this->x().end(), x ); + } + + if ( this->x().end() == xIter ) { + + return Z( 0. ); + } + else if ( this->x().begin() == xIter ) { + + return Z( 0. ); + } + else { + + auto fIter = this->f().begin(); + std::advance( fIter, std::distance( this->x().begin(), xIter ) ); + + return this->interpolator()( x, y, *std::prev( xIter ), *xIter, + *std::prev( fIter ), *fIter ); + } +} diff --git a/src/scion/math/SingleTableFunctionBase/src/verifyTable.hpp b/src/scion/math/SingleTableFunctionBase/src/verifyTable.hpp new file mode 100644 index 0000000..ee0c02a --- /dev/null +++ b/src/scion/math/SingleTableFunctionBase/src/verifyTable.hpp @@ -0,0 +1,43 @@ +static void verifyTable( const XContainer& x, const FContainer& f ) { + + if ( ( ! verification::isAtLeastOfSize( x, 2 ) ) || + ( ! verification::isAtLeastOfSize( f, 2 ) ) ) { + + Log::error( "Insufficient x or f values defined for x,f(y) tabulated data " + "with a single interpolation type (at least 2 points are " + "required)" ); + Log::info( "x.size(): {}", x.size() ); + Log::info( "f.size(): {}", f.size() ); + throw std::exception(); + } + + if ( ! verification::isSameSize( x, f ) ) { + + Log::error( "Inconsistent number of x and f values for x,f(y) tabulated data " + "with a single interpolation type" ); + Log::info( "x.size(): {}", x.size() ); + Log::info( "f.size(): {}", f.size() ); + throw std::exception(); + } + + if ( ! verification::isSorted( x ) ) { + + Log::error( "The x values do not appear to be in ascending order for " + "x,f(y) tabulated data with a single interpolation type" ); + throw std::exception(); + } + + if ( ! verification::isUnique( x ) ) { + + Log::error( "The x values do not appear to be unique for x,f(y) tabulated values " + "with a single interpolation type" ); + + auto iter = std::adjacent_find( x.begin(), x.end() ); + while ( iter != x.end() ) { + + Log::info( "Duplicate x value found: {}", *iter ); + iter = std::adjacent_find( ++iter, x.end() ); + } + throw std::exception(); + } +} diff --git a/src/scion/math/TwoDimensionalFunctionBase.hpp b/src/scion/math/TwoDimensionalFunctionBase.hpp new file mode 100644 index 0000000..3479173 --- /dev/null +++ b/src/scion/math/TwoDimensionalFunctionBase.hpp @@ -0,0 +1,41 @@ +#ifndef NJOY_SCION_MATH_TWODIMENSIONALFUNCTIONBASE +#define NJOY_SCION_MATH_TWODIMENSIONALFUNCTIONBASE + +// system includes + +// other includes + +namespace njoy { +namespace scion { +namespace math { + + /** + * @class + * @brief Base class for two dimensional function objects modelling y = f(x,y) + * + * This base class provides the common interface for all two dimensional + * function objects. This includes function evaluation. + */ + template < typename Derived, typename X, typename Y, typename Z > + class TwoDimensionalFunctionBase { + + public: + + /* methods */ + + /** + * @brief Evaluate the function for a value of x and y + * + * @param x the value to be evaluated + */ + Z operator()( const X& x, const Y& y ) const { + + return static_cast< const Derived* >( this )->evaluate( x, y ); + } + }; + +} // math namespace +} // scion namespace +} // njoy namespace + +#endif