From 40f58f74abc6737ce07a7792531aa64fea078ba0 Mon Sep 17 00:00:00 2001 From: Wim Haeck Date: Thu, 6 Jun 2024 11:11:54 -0600 Subject: [PATCH] Making a correction to InterpolationTable --- .../src/processBoundaries.hpp | 39 ++++- .../test/InterpolationTable.test.cpp | 157 ++++++++++++++++++ 2 files changed, 192 insertions(+), 4 deletions(-) diff --git a/src/scion/math/InterpolationTable/src/processBoundaries.hpp b/src/scion/math/InterpolationTable/src/processBoundaries.hpp index 0812193..72def51 100644 --- a/src/scion/math/InterpolationTable/src/processBoundaries.hpp +++ b/src/scion/math/InterpolationTable/src/processBoundaries.hpp @@ -1,3 +1,27 @@ +/** + * @brief Verify and correct boundaries and interpolants + * + * This function does a lot, so here's an overview of what it does. First of all, it verifies + * the following things: + * - There are at least 2 values in the x and y grid + * - The x and y grid have the same size + * - The number of boundaries and interpolants are the same + * - The last boundary index is equal to the index of the last x value + * - The x grid is sorted + * - There is no jump at the beginning or end of the x grid + * - The x values appear only a maximum of two times in the grid + * + * Next, this function will look for every jump in the x grid and check if the jump corresponds + * to a change in interpolation region (meaning that the index of the first x value in the jump + * is in the boundaries). If that is not the case, an additional interpolation region wil be + * inserted. This ensures that none of the interpolation regions will contain a jump in their + * local x grid. + * + * In some cases, the boundary values can point to the second point of a jump. While this is not + * an error (we will never interpolate on a jump), we need the boundaries to point to the first + * point in the jump instead of the second one. When this is encountered, the boundary value is + * adjusted. + */ static std::tuple< std::vector< std::size_t >, std::vector< interpolation::InterpolationType > > processBoundaries( const std::vector< X >& x, const std::vector< Y >& y, @@ -59,10 +83,17 @@ processBoundaries( const std::vector< X >& x, const std::vector< Y >& y, bIter = std::lower_bound( bIter, boundaries.end(), index ); if ( *bIter != index ) { - iIter = std::next( interpolants.begin(), - std::distance( boundaries.begin(), bIter ) ); - bIter = boundaries.insert( bIter, index ); - iIter = interpolants.insert( iIter, *iIter ); + if ( *bIter == index + 1 ) { + + *bIter -= 1; + } + else { + + iIter = std::next( interpolants.begin(), + std::distance( boundaries.begin(), bIter ) ); + bIter = boundaries.insert( bIter, index ); + iIter = interpolants.insert( iIter, *iIter ); + } } ++xIter; if ( std::next( xIter ) == x.end() ) { diff --git a/src/scion/math/InterpolationTable/test/InterpolationTable.test.cpp b/src/scion/math/InterpolationTable/test/InterpolationTable.test.cpp index b878054..8ad3aa1 100644 --- a/src/scion/math/InterpolationTable/test/InterpolationTable.test.cpp +++ b/src/scion/math/InterpolationTable/test/InterpolationTable.test.cpp @@ -1687,6 +1687,163 @@ SCENARIO( "InterpolationTable" ) { } // WHEN } // GIVEN + GIVEN( "non-linearised data with multiple regions with a jump and boundaries " + "that point to the second x value in the jump" ) { + + // note: at construction time, the boundary value will be set to the first point in + // the jump. As a result, the final data contained in this InterpolationTable is the + // same as the previous test. + + WHEN( "the data is given explicitly" ) { + + const std::vector< double > x = { 1., 2., 2., 3., 4. }; + const std::vector< double > y = { 4., 3., 4., 3., 2. }; + const std::vector< std::size_t > boundaries = { 2, 4 }; + const std::vector< InterpolationType > interpolants = { + + InterpolationType::LinearLinear, + InterpolationType::LinearLog + }; + + InterpolationTable< double > chunk( std::move( x ), std::move( y ), + std::move( boundaries ), + std::move( interpolants ) ); + + THEN( "a InterpolationTable can be constructed and members can be tested" ) { + + CHECK( 5 == chunk.x().size() ); + CHECK( 5 == chunk.y().size() ); + CHECK( 2 == chunk.boundaries().size() ); + CHECK( 2 == chunk.interpolants().size() ); + CHECK_THAT( 1., WithinRel( chunk.x()[0] ) ); + CHECK_THAT( 2., WithinRel( chunk.x()[1] ) ); + CHECK_THAT( 2., WithinRel( chunk.x()[2] ) ); + CHECK_THAT( 3., WithinRel( chunk.x()[3] ) ); + CHECK_THAT( 4., WithinRel( chunk.x()[4] ) ); + CHECK_THAT( 4., WithinRel( chunk.y()[0] ) ); + CHECK_THAT( 3., WithinRel( chunk.y()[1] ) ); + CHECK_THAT( 4., WithinRel( chunk.y()[2] ) ); + CHECK_THAT( 3., WithinRel( chunk.y()[3] ) ); + CHECK_THAT( 2., WithinRel( chunk.y()[4] ) ); + CHECK( 1 == chunk.boundaries()[0] ); // <-- this is changed from 2 to 1 + CHECK( 4 == chunk.boundaries()[1] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[0] ); + CHECK( InterpolationType::LinearLog == chunk.interpolants()[1] ); + CHECK( false == chunk.isLinearised() ); + + CHECK( true == std::holds_alternative< IntervalDomain< double > >( chunk.domain() ) ); + } // THEN + + THEN( "a InterpolationTable can be evaluated" ) { + + // values of x in the x grid + CHECK_THAT( 4., WithinRel( chunk( 1. ) ) ); + CHECK_THAT( 4., WithinRel( chunk( 2. ) ) ); + CHECK_THAT( 3., WithinRel( chunk( 3. ) ) ); + CHECK_THAT( 2., WithinRel( chunk( 4. ) ) ); + + // values of x outside the x grid + CHECK_THAT( 0., WithinRel( chunk( 0. ) ) ); + CHECK_THAT( 0., WithinRel( chunk( 5. ) ) ); + + // values of x inside the x grid (lin-lin piece) + CHECK_THAT( 3.5, WithinRel( chunk( 1.5 ) ) ); + + // values of x inside the x grid (lin-log piece) + CHECK_THAT( 3.44966028678679, WithinRel( chunk( 2.5 ) ) ); + CHECK_THAT( 2.46416306545103, WithinRel( chunk( 3.5 ) ) ); + } // THEN + + THEN( "the domain can be tested" ) { + + CHECK( true == chunk.isInside( 1.0 ) ); + CHECK( true == chunk.isInside( 2.5 ) ); + CHECK( true == chunk.isInside( 4.0 ) ); + + CHECK( false == chunk.isContained( 1.0 ) ); + CHECK( true == chunk.isContained( 2.5 ) ); + CHECK( false == chunk.isContained( 4.0 ) ); + + CHECK( true == chunk.isSameDomain( IntervalDomain< double >( 1., 4. ) ) ); + CHECK( false == chunk.isSameDomain( IntervalDomain< double >( 0., 4. ) ) ); + CHECK( false == chunk.isSameDomain( OpenDomain< double >() ) ); + } // THEN + + THEN( "an InterpolationTable can be linearised" ) { + + InterpolationTable< double > linear = chunk.linearise(); + + CHECK( 12 == linear.numberPoints() ); + CHECK( 2 == linear.numberRegions() ); + + CHECK( 12 == linear.x().size() ); + CHECK( 12 == linear.y().size() ); + CHECK( 2 == linear.boundaries().size() ); + CHECK( 2 == linear.interpolants().size() ); + + CHECK( 1 == linear.boundaries()[0] ); + CHECK( 11 == linear.boundaries()[1] ); + + CHECK( InterpolationType::LinearLinear == linear.interpolants()[0] ); + CHECK( InterpolationType::LinearLinear == linear.interpolants()[1] ); + + CHECK_THAT( 1. , WithinRel( linear.x()[0] ) ); + CHECK_THAT( 2. , WithinRel( linear.x()[1] ) ); + CHECK_THAT( 2. , WithinRel( linear.x()[2] ) ); + CHECK_THAT( 2.125, WithinRel( linear.x()[3] ) ); + CHECK_THAT( 2.25 , WithinRel( linear.x()[4] ) ); + CHECK_THAT( 2.5 , WithinRel( linear.x()[5] ) ); + CHECK_THAT( 2.75 , WithinRel( linear.x()[6] ) ); + CHECK_THAT( 3. , WithinRel( linear.x()[7] ) ); + CHECK_THAT( 3.25 , WithinRel( linear.x()[8] ) ); + CHECK_THAT( 3.5 , WithinRel( linear.x()[9] ) ); + CHECK_THAT( 3.75 , WithinRel( linear.x()[10] ) ); + CHECK_THAT( 4. , WithinRel( linear.x()[11] ) ); + + CHECK_THAT( 4. , WithinRel( linear.y()[0] ) ); + CHECK_THAT( 3. , WithinRel( linear.y()[1] ) ); + CHECK_THAT( 4. , WithinRel( linear.y()[2] ) ); + CHECK_THAT( 3.85048128530886, WithinRel( linear.y()[3] ) ); + CHECK_THAT( 3.70951129135145, WithinRel( linear.y()[4] ) ); + CHECK_THAT( 3.44966028678679, WithinRel( linear.y()[5] ) ); + CHECK_THAT( 3.21459646033567, WithinRel( linear.y()[6] ) ); + CHECK_THAT( 3. , WithinRel( linear.y()[7] ) ); + CHECK_THAT( 2.72176678584324, WithinRel( linear.y()[8] ) ); + CHECK_THAT( 2.46416306545103, WithinRel( linear.y()[9] ) ); + CHECK_THAT( 2.22433973930853, WithinRel( linear.y()[10] ) ); + CHECK_THAT( 2. , WithinRel( linear.y()[11] ) ); + + CHECK( true == linear.isLinearised() ); + } // THEN + + THEN( "arithmetic operations cannot be performed" ) { + + InterpolationTable< double > result( { 1., 4. }, { 0., 0. } ); + InterpolationTable< double > right( { 1., 4. }, { 0., 0. } ); + + // scalar operations + CHECK_THROWS( chunk += 2. ); + CHECK_THROWS( chunk -= 2. ); + CHECK_THROWS( chunk *= 2. ); + CHECK_THROWS( chunk /= 2. ); + CHECK_THROWS( result = -chunk ); + CHECK_THROWS( result = chunk + 2. ); + CHECK_THROWS( result = chunk - 2. ); + CHECK_THROWS( result = chunk * 2. ); + CHECK_THROWS( result = chunk / 2. ); + CHECK_THROWS( result = 2. + chunk ); + CHECK_THROWS( result = 2. - chunk ); + CHECK_THROWS( result = 2. * chunk ); + + // table operations + CHECK_THROWS( chunk += right ); + CHECK_THROWS( chunk -= right ); + CHECK_THROWS( result = chunk + right ); + CHECK_THROWS( result = chunk - right ); + } // THEN + } // WHEN + } // GIVEN + GIVEN( "invalid data for an InterpolationTable object" ) { WHEN( "there are not enough values in the x or y grid" ) {