Making a correction to InterpolationTable

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() ) {

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" ) {