diff --git a/python/test/math/Test_scion_math_InterpolationTable.py b/python/test/math/Test_scion_math_InterpolationTable.py index aaac85d..7ec6498 100644 --- a/python/test/math/Test_scion_math_InterpolationTable.py +++ b/python/test/math/Test_scion_math_InterpolationTable.py @@ -94,7 +94,7 @@ def verify_chunk1( self, chunk ) : # verify arithmetic operators same = InterpolationTable( [ 1., 4. ], [ 0., 3. ] ) threshold = InterpolationTable( [ 2., 4. ], [ 0., 2. ] ) - nonzerothreshold = InterpolationTable( [ 2., 4. ], [ 1., 2. ] ) + nonzerothreshold = InterpolationTable( [ 2., 4. ], [ 1., 3. ] ) small = InterpolationTable( [ 1., 3. ], [ 0., 2. ] ) result = -chunk @@ -313,7 +313,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, result.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) - chunk += same; + chunk += same self.assertEqual( 4, chunk.number_points ) self.assertEqual( 1, chunk.number_regions ) self.assertEqual( 4, len( chunk.x ) ) @@ -331,7 +331,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, chunk.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) - chunk -= same; + chunk -= same self.assertEqual( 4, chunk.number_points ) self.assertEqual( 1, chunk.number_regions ) self.assertEqual( 4, len( chunk.x ) ) @@ -349,7 +349,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, chunk.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) - result = chunk + same; + result = chunk + same self.assertEqual( 4, result.number_points ) self.assertEqual( 1, result.number_regions ) self.assertEqual( 4, len( result.x ) ) @@ -367,7 +367,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, result.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) - result = chunk - same; + result = chunk - same self.assertEqual( 4, result.number_points ) self.assertEqual( 1, result.number_regions ) self.assertEqual( 4, len( result.x ) ) @@ -385,7 +385,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, result.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) - chunk += threshold; + chunk += threshold self.assertEqual( 4, chunk.number_points ) self.assertEqual( 1, chunk.number_regions ) self.assertEqual( 4, len( chunk.x ) ) @@ -403,7 +403,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, chunk.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) - chunk -= threshold; + chunk -= threshold self.assertEqual( 4, chunk.number_points ) self.assertEqual( 1, chunk.number_regions ) self.assertEqual( 4, len( chunk.x ) ) @@ -421,7 +421,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, chunk.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) - result = chunk + threshold; + result = chunk + threshold self.assertEqual( 4, result.number_points ) self.assertEqual( 1, result.number_regions ) self.assertEqual( 4, len( result.x ) ) @@ -439,7 +439,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, result.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) - result = chunk - threshold; + result = chunk - threshold self.assertEqual( 4, result.number_points ) self.assertEqual( 1, result.number_regions ) self.assertEqual( 4, len( result.x ) ) @@ -457,8 +457,92 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 3, result.boundaries[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) + chunk += nonzerothreshold + self.assertEqual( 5, chunk.number_points ) + self.assertEqual( 2, chunk.number_regions ) + self.assertEqual( 5, len( chunk.x ) ) + self.assertEqual( 5, len( chunk.y ) ) + self.assertEqual( 2, len( chunk.boundaries ) ) + self.assertEqual( 2, len( chunk.interpolants ) ) + self.assertAlmostEqual( 1. , chunk.x[0] ) + self.assertAlmostEqual( 2. , chunk.x[1] ) + self.assertAlmostEqual( 2. , chunk.x[2] ) + self.assertAlmostEqual( 3. , chunk.x[3] ) + self.assertAlmostEqual( 4. , chunk.x[4] ) + self.assertAlmostEqual( 4.0, chunk.y[0] ) + self.assertAlmostEqual( 3.0, chunk.y[1] ) + self.assertAlmostEqual( 4.0, chunk.y[2] ) + self.assertAlmostEqual( 4.0, chunk.y[3] ) + self.assertAlmostEqual( 4.0, chunk.y[4] ) + self.assertEqual( 1, chunk.boundaries[0] ) + self.assertEqual( 4, chunk.boundaries[1] ) + self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) + self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[1] ) + + chunk -= nonzerothreshold + self.assertEqual( 4, chunk.number_points ) + self.assertEqual( 1, chunk.number_regions ) + self.assertEqual( 4, len( chunk.x ) ) + self.assertEqual( 4, len( chunk.y ) ) + self.assertEqual( 1, len( chunk.boundaries ) ) + self.assertEqual( 1, len( chunk.interpolants ) ) + self.assertAlmostEqual( 1. , chunk.x[0] ) + self.assertAlmostEqual( 2. , chunk.x[1] ) + self.assertAlmostEqual( 3. , chunk.x[2] ) + self.assertAlmostEqual( 4. , chunk.x[3] ) + self.assertAlmostEqual( 4.0, chunk.y[0] ) + self.assertAlmostEqual( 3.0, chunk.y[1] ) + self.assertAlmostEqual( 2.0, chunk.y[2] ) + self.assertAlmostEqual( 1.0, chunk.y[3] ) + self.assertEqual( 3, chunk.boundaries[0] ) + self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) + + result = chunk + nonzerothreshold + self.assertEqual( 5, result.number_points ) + self.assertEqual( 2, result.number_regions ) + self.assertEqual( 5, len( result.x ) ) + self.assertEqual( 5, len( result.y ) ) + self.assertEqual( 2, len( result.boundaries ) ) + self.assertEqual( 2, len( result.interpolants ) ) + self.assertAlmostEqual( 1. , result.x[0] ) + self.assertAlmostEqual( 2. , result.x[1] ) + self.assertAlmostEqual( 2. , result.x[2] ) + self.assertAlmostEqual( 3. , result.x[3] ) + self.assertAlmostEqual( 4. , result.x[4] ) + self.assertAlmostEqual( 4.0, result.y[0] ) + self.assertAlmostEqual( 3.0, result.y[1] ) + self.assertAlmostEqual( 4.0, result.y[2] ) + self.assertAlmostEqual( 4.0, result.y[3] ) + self.assertAlmostEqual( 4.0, result.y[4] ) + self.assertEqual( 1, result.boundaries[0] ) + self.assertEqual( 4, result.boundaries[1] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) + + result = chunk - nonzerothreshold + self.assertEqual( 5, result.number_points ) + self.assertEqual( 2, result.number_regions ) + self.assertEqual( 5, len( result.x ) ) + self.assertEqual( 5, len( result.y ) ) + self.assertEqual( 2, len( result.boundaries ) ) + self.assertEqual( 2, len( result.interpolants ) ) + self.assertAlmostEqual( 1. , result.x[0] ) + self.assertAlmostEqual( 2. , result.x[1] ) + self.assertAlmostEqual( 2. , result.x[2] ) + self.assertAlmostEqual( 3. , result.x[3] ) + self.assertAlmostEqual( 4. , result.x[4] ) + self.assertAlmostEqual( 4., result.y[0] ) + self.assertAlmostEqual( 3., result.y[1] ) + self.assertAlmostEqual( 2., result.y[2] ) + self.assertAlmostEqual( 0., result.y[3] ) + self.assertAlmostEqual( -2., result.y[4] ) + self.assertEqual( 1, result.boundaries[0] ) + self.assertEqual( 4, result.boundaries[1] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) + # this will add a second point at the lower end point - result = chunk + small; + result = chunk + small self.assertEqual( 5, result.number_points ) self.assertEqual( 2, result.number_regions ) self.assertEqual( 5, len( result.x ) ) @@ -480,7 +564,7 @@ def verify_chunk1( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) # this will add a second point at the lower end point - result = chunk - small; + result = chunk - small self.assertEqual( 5, result.number_points ) self.assertEqual( 2, result.number_regions ) self.assertEqual( 5, len( result.x ) ) @@ -501,12 +585,6 @@ def verify_chunk1( self, chunk ) : self.assertEqual( 4, result.boundaries[1] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) - # the threshold table starts with a non-zero value - with self.assertRaises( Exception ) : chunk += nonzerothreshold - with self.assertRaises( Exception ) : chunk -= nonzerothreshold - with self.assertRaises( Exception ) : result = chunk + nonzerothreshold - with self.assertRaises( Exception ) : result = chunk - nonzerothreshold - def verify_chunk2( self, chunk ) : # verify content @@ -592,7 +670,7 @@ def verify_chunk2( self, chunk ) : # verify arithmetic operators same = InterpolationTable( [ 1., 4. ], [ 0., 3. ] ) threshold = InterpolationTable( [ 2., 4. ], [ 0., 2. ] ) - nonzerothreshold = InterpolationTable( [ 2., 4. ], [ 1., 2. ] ) + nonzerothreshold = InterpolationTable( [ 3., 4. ], [ 1., 2. ] ) small = InterpolationTable( [ 1., 3. ], [ 0., 2. ] ) result = -chunk @@ -835,7 +913,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) - chunk += same; + chunk += same self.assertEqual( 5, len( chunk.x ) ) self.assertEqual( 5, len( chunk.y ) ) self.assertEqual( 2, len( chunk.boundaries ) ) @@ -855,7 +933,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[1] ) - chunk -= same; + chunk -= same self.assertEqual( 5, len( chunk.x ) ) self.assertEqual( 5, len( chunk.y ) ) self.assertEqual( 2, len( chunk.boundaries ) ) @@ -875,7 +953,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[1] ) - result = chunk + same; + result = chunk + same self.assertEqual( 5, len( result.x ) ) self.assertEqual( 5, len( result.y ) ) self.assertEqual( 2, len( result.boundaries ) ) @@ -895,7 +973,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) - result = chunk - same; + result = chunk - same self.assertEqual( 5, len( result.x ) ) self.assertEqual( 5, len( result.y ) ) self.assertEqual( 2, len( result.boundaries ) ) @@ -915,7 +993,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) - chunk += threshold; + chunk += threshold self.assertEqual( 5, len( chunk.x ) ) self.assertEqual( 5, len( chunk.y ) ) self.assertEqual( 2, len( chunk.boundaries ) ) @@ -935,7 +1013,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[1] ) - chunk -= threshold; + chunk -= threshold self.assertEqual( 5, len( chunk.x ) ) self.assertEqual( 5, len( chunk.y ) ) self.assertEqual( 2, len( chunk.boundaries ) ) @@ -955,7 +1033,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[1] ) - result = chunk + threshold; + result = chunk + threshold self.assertEqual( 5, len( result.x ) ) self.assertEqual( 5, len( result.y ) ) self.assertEqual( 2, len( result.boundaries ) ) @@ -975,7 +1053,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) - result = chunk - threshold; + result = chunk - threshold self.assertEqual( 5, len( result.x ) ) self.assertEqual( 5, len( result.y ) ) self.assertEqual( 2, len( result.boundaries ) ) @@ -995,8 +1073,100 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) + chunk += nonzerothreshold + self.assertEqual( 6, len( chunk.x ) ) + self.assertEqual( 6, len( chunk.y ) ) + self.assertEqual( 3, len( chunk.boundaries ) ) + self.assertEqual( 3, len( chunk.interpolants ) ) + self.assertAlmostEqual( 1., chunk.x[0] ) + self.assertAlmostEqual( 2., chunk.x[1] ) + self.assertAlmostEqual( 2., chunk.x[2] ) + self.assertAlmostEqual( 3., chunk.x[3] ) + self.assertAlmostEqual( 3., chunk.x[4] ) + self.assertAlmostEqual( 4., chunk.x[5] ) + self.assertAlmostEqual( 4., chunk.y[0] ) + self.assertAlmostEqual( 3., chunk.y[1] ) + self.assertAlmostEqual( 4., chunk.y[2] ) + self.assertAlmostEqual( 3., chunk.y[3] ) + self.assertAlmostEqual( 4., chunk.y[4] ) + self.assertAlmostEqual( 4., chunk.y[5] ) + self.assertEqual( 1, chunk.boundaries[0] ) + self.assertEqual( 3, chunk.boundaries[1] ) + self.assertEqual( 5, chunk.boundaries[2] ) + self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) + self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[1] ) + self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[2] ) + + chunk -= nonzerothreshold + self.assertEqual( 5, len( chunk.x ) ) + self.assertEqual( 5, len( chunk.y ) ) + self.assertEqual( 2, len( chunk.boundaries ) ) + self.assertEqual( 2, len( chunk.interpolants ) ) + self.assertAlmostEqual( 1., chunk.x[0] ) + self.assertAlmostEqual( 2., chunk.x[1] ) + self.assertAlmostEqual( 2., chunk.x[2] ) + self.assertAlmostEqual( 3., chunk.x[3] ) + self.assertAlmostEqual( 4., chunk.x[4] ) + self.assertAlmostEqual( 4., chunk.y[0] ) + self.assertAlmostEqual( 3., chunk.y[1] ) + self.assertAlmostEqual( 4., chunk.y[2] ) + self.assertAlmostEqual( 3., chunk.y[3] ) + self.assertAlmostEqual( 2., chunk.y[4] ) + self.assertEqual( 1, chunk.boundaries[0] ) + self.assertEqual( 4, chunk.boundaries[1] ) + self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[0] ) + self.assertEqual( InterpolationType.LinearLinear, chunk.interpolants[1] ) + + result = chunk + nonzerothreshold + self.assertEqual( 6, len( result.x ) ) + self.assertEqual( 6, len( result.y ) ) + self.assertEqual( 3, len( result.boundaries ) ) + self.assertEqual( 3, len( result.interpolants ) ) + self.assertAlmostEqual( 1., result.x[0] ) + self.assertAlmostEqual( 2., result.x[1] ) + self.assertAlmostEqual( 2., result.x[2] ) + self.assertAlmostEqual( 3., result.x[3] ) + self.assertAlmostEqual( 3., result.x[4] ) + self.assertAlmostEqual( 4., result.x[5] ) + self.assertAlmostEqual( 4., result.y[0] ) + self.assertAlmostEqual( 3., result.y[1] ) + self.assertAlmostEqual( 4., result.y[2] ) + self.assertAlmostEqual( 3., result.y[3] ) + self.assertAlmostEqual( 4., result.y[4] ) + self.assertAlmostEqual( 4., result.y[5] ) + self.assertEqual( 1, result.boundaries[0] ) + self.assertEqual( 3, result.boundaries[1] ) + self.assertEqual( 5, result.boundaries[2] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[2] ) + + result = chunk - nonzerothreshold + self.assertEqual( 6, len( result.x ) ) + self.assertEqual( 6, len( result.y ) ) + self.assertEqual( 3, len( result.boundaries ) ) + self.assertEqual( 3, len( result.interpolants ) ) + self.assertAlmostEqual( 1., result.x[0] ) + self.assertAlmostEqual( 2., result.x[1] ) + self.assertAlmostEqual( 2., result.x[2] ) + self.assertAlmostEqual( 3., result.x[3] ) + self.assertAlmostEqual( 3., result.x[4] ) + self.assertAlmostEqual( 4., result.x[5] ) + self.assertAlmostEqual( 4., result.y[0] ) + self.assertAlmostEqual( 3., result.y[1] ) + self.assertAlmostEqual( 4., result.y[2] ) + self.assertAlmostEqual( 3., result.y[3] ) + self.assertAlmostEqual( 2., result.y[4] ) + self.assertAlmostEqual( 0., result.y[5] ) + self.assertEqual( 1, result.boundaries[0] ) + self.assertEqual( 3, result.boundaries[1] ) + self.assertEqual( 5, result.boundaries[2] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[0] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) + self.assertEqual( InterpolationType.LinearLinear, result.interpolants[2] ) + # this will add a second point at the lower end point - result = chunk + small; + result = chunk + small self.assertEqual( 6, len( result.x ) ) self.assertEqual( 6, len( result.y ) ) self.assertEqual( 3, len( result.boundaries ) ) @@ -1021,7 +1191,7 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, result.interpolants[2] ) # this will add a second point at the lower end point - result = chunk - small; + result = chunk - small self.assertEqual( 6, len( result.x ) ) self.assertEqual( 6, len( result.y ) ) self.assertEqual( 3, len( result.boundaries ) ) @@ -1045,12 +1215,6 @@ def verify_chunk2( self, chunk ) : self.assertEqual( InterpolationType.LinearLinear, result.interpolants[1] ) self.assertEqual( InterpolationType.LinearLinear, result.interpolants[2] ) - # the threshold table starts with a non-zero value - with self.assertRaises( Exception ) : chunk += nonzerothreshold - with self.assertRaises( Exception ) : chunk -= nonzerothreshold - with self.assertRaises( Exception ) : result = chunk + nonzerothreshold - with self.assertRaises( Exception ) : result = chunk - nonzerothreshold - def verify_chunk3( self, chunk ) : # verify content @@ -1366,6 +1530,12 @@ def test_failures( self ) : chunk = InterpolationTable( x = [ 1., 2., 2., 2., 3., 4. ], y = [ 4., 3., 3., 3., 2., 1. ] ) + # the x grid has a jump at the beginning + with self.assertRaises( Exception ) : + + chunk = InterpolationTable( x = [ 1., 1., 3., 4. ], + y = [ 4., 3., 1., 4. ] ) + # the x grid has a jump at the end with self.assertRaises( Exception ) : diff --git a/src/scion/math/InterpolationTable/src/evaluateOnGrid.hpp b/src/scion/math/InterpolationTable/src/evaluateOnGrid.hpp index 9a3238b..87d9c35 100644 --- a/src/scion/math/InterpolationTable/src/evaluateOnGrid.hpp +++ b/src/scion/math/InterpolationTable/src/evaluateOnGrid.hpp @@ -3,6 +3,10 @@ std::vector< Y > evaluateOnGrid( const std::vector< X >& x ) const { std::vector< Y > y( x.size(), Y( 0. ) ); auto xIter = std::lower_bound( x.begin(), x.end(), this->x().front() ); + if ( *std::next( xIter ) == this->x().front() ) { + + ++xIter; + } auto yIter = std::next( y.begin(), std::distance( x.begin(), xIter ) ); auto xTable = this->x().begin(); diff --git a/src/scion/math/InterpolationTable/src/generateTables.hpp b/src/scion/math/InterpolationTable/src/generateTables.hpp index ddd6b64..ba63dec 100644 --- a/src/scion/math/InterpolationTable/src/generateTables.hpp +++ b/src/scion/math/InterpolationTable/src/generateTables.hpp @@ -4,8 +4,6 @@ void generateTables() { auto xStart = this->x().begin(); auto yStart = this->y().begin(); - auto xEnd = xStart; - auto yEnd = yStart; std::size_t nr = this->boundaries().size(); bool linearised = true; for ( std::size_t i = 0; i < nr; ++i ) { diff --git a/src/scion/math/InterpolationTable/src/operation.hpp b/src/scion/math/InterpolationTable/src/operation.hpp index f16ed54..22c949b 100644 --- a/src/scion/math/InterpolationTable/src/operation.hpp +++ b/src/scion/math/InterpolationTable/src/operation.hpp @@ -33,30 +33,27 @@ InterpolationTable& operation( const InterpolationTable& right, } else { - // check for threshold tables - if ( this->x().front() != right.x().front() ) { - - Y ystart = this->x().front() < right.x().front() ? right.y().front() - : this->y().front(); - - if ( Y( 0. ) != ystart ) { - - X xstart = this->x().front() < right.x().front() ? right.x().front() - : this->x().front(); - - Log::error( "The threshold table's first y value is not zero" ); - Log::info( "Found x = {}", xstart ); - Log::info( "Found y = {}", ystart ); - throw std::exception(); - } - } - // unionise and evaluate on the new grid std::vector< X > x = unionisation::unionise( this->x(), right.x() ); std::vector< Y > y = this->evaluateOnGrid( x ); std::vector< Y > temp = right.evaluateOnGrid( x ); std::transform( y.begin(), y.end(), temp.begin(), y.begin(), operation ); + // check for threshold jump with the same y value + if ( this->x().front() != right.x().front() ) { + + X value = this->x().front() < right.x().front() ? right.x().front() + : this->x().front(); + auto xIter = std::lower_bound( x.begin(), x.end(), value ); + auto yIter = std::next( y.begin(), std::distance( x.begin(), xIter ) ); + if ( *std::next( yIter ) == *yIter ) { + + x.erase( xIter ); + y.erase( yIter ); + } + } + + // replace this with a new table *this = InterpolationTable( std::move( x ), std::move( y ) ); } diff --git a/src/scion/math/InterpolationTable/src/processBoundaries.hpp b/src/scion/math/InterpolationTable/src/processBoundaries.hpp index 2017ade..0812193 100644 --- a/src/scion/math/InterpolationTable/src/processBoundaries.hpp +++ b/src/scion/math/InterpolationTable/src/processBoundaries.hpp @@ -45,6 +45,12 @@ processBoundaries( const std::vector< X >& x, const std::vector< Y >& y, } auto xIter = std::adjacent_find( x.begin(), x.end() ); + if ( xIter == x.begin() ) { + + Log::error( "A jump in the x grid cannot occur at the beginning of the x grid" ); + throw std::exception(); + } + auto bIter = boundaries.begin(); auto iIter = interpolants.begin(); while ( xIter != x.end() ) { diff --git a/src/scion/math/InterpolationTable/test/InterpolationTable.test.cpp b/src/scion/math/InterpolationTable/test/InterpolationTable.test.cpp index 122aaf9..b878054 100644 --- a/src/scion/math/InterpolationTable/test/InterpolationTable.test.cpp +++ b/src/scion/math/InterpolationTable/test/InterpolationTable.test.cpp @@ -89,7 +89,7 @@ SCENARIO( "InterpolationTable" ) { InterpolationTable< double > result( { 1., 4. }, { 0., 0. } ); InterpolationTable< double > same( { 1., 4. }, { 0., 3. } ); InterpolationTable< double > threshold( { 2., 4. }, { 0., 2. } ); - InterpolationTable< double > nonzerothreshold( { 2., 4. }, { 1., 2. } ); + InterpolationTable< double > nonzerothreshold( { 2., 4. }, { 1., 3. } ); InterpolationTable< double > small( { 1., 3. }, { 0., 2. } ); chunk += 2.; @@ -492,6 +492,98 @@ SCENARIO( "InterpolationTable" ) { CHECK( InterpolationType::LinearLinear == result.interpolants()[0] ); CHECK( true == std::holds_alternative< IntervalDomain< double > >( result.domain() ) ); + chunk += nonzerothreshold; + + CHECK( 5 == chunk.numberPoints() ); + CHECK( 2 == chunk.numberRegions() ); + 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( 4., WithinRel( chunk.y()[3] ) ); + CHECK_THAT( 4., WithinRel( chunk.y()[4] ) ); + CHECK( 1 == chunk.boundaries()[0] ); + CHECK( 4 == chunk.boundaries()[1] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[0] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[1] ); + CHECK( true == std::holds_alternative< IntervalDomain< double > >( chunk.domain() ) ); + + chunk -= nonzerothreshold; + + CHECK( 4 == chunk.numberPoints() ); + CHECK( 1 == chunk.numberRegions() ); + CHECK( 4 == chunk.x().size() ); + CHECK( 4 == chunk.y().size() ); + CHECK( 1 == chunk.boundaries().size() ); + CHECK( 1 == chunk.interpolants().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.0, WithinRel( chunk.y()[0] ) ); + CHECK_THAT( 3.0, WithinRel( chunk.y()[1] ) ); + CHECK_THAT( 2.0, WithinRel( chunk.y()[2] ) ); + CHECK_THAT( 1.0, WithinRel( chunk.y()[3] ) ); + CHECK( 3 == chunk.boundaries()[0] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[0] ); + CHECK( true == std::holds_alternative< IntervalDomain< double > >( chunk.domain() ) ); + + result = chunk + nonzerothreshold; + + CHECK( 5 == result.numberPoints() ); + CHECK( 2 == result.numberRegions() ); + CHECK( 5 == result.x().size() ); + CHECK( 5 == result.y().size() ); + CHECK( 2 == result.boundaries().size() ); + CHECK( 2 == result.interpolants().size() ); + CHECK_THAT( 1., WithinRel( result.x()[0] ) ); + CHECK_THAT( 2., WithinRel( result.x()[1] ) ); + CHECK_THAT( 2., WithinRel( result.x()[2] ) ); + CHECK_THAT( 3., WithinRel( result.x()[3] ) ); + CHECK_THAT( 4., WithinRel( result.x()[4] ) ); + CHECK_THAT( 4., WithinRel( result.y()[0] ) ); + CHECK_THAT( 3., WithinRel( result.y()[1] ) ); + CHECK_THAT( 4., WithinRel( result.y()[2] ) ); + CHECK_THAT( 4., WithinRel( result.y()[3] ) ); + CHECK_THAT( 4., WithinRel( result.y()[4] ) ); + CHECK( 1 == result.boundaries()[0] ); + CHECK( 4 == result.boundaries()[1] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[0] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[1] ); + CHECK( true == std::holds_alternative< IntervalDomain< double > >( result.domain() ) ); + + result = chunk - nonzerothreshold; + + CHECK( 5 == result.numberPoints() ); + CHECK( 2 == result.numberRegions() ); + CHECK( 5 == result.x().size() ); + CHECK( 5 == result.y().size() ); + CHECK( 2 == result.boundaries().size() ); + CHECK( 2 == result.interpolants().size() ); + CHECK_THAT( 1., WithinRel( result.x()[0] ) ); + CHECK_THAT( 2., WithinRel( result.x()[1] ) ); + CHECK_THAT( 2., WithinRel( result.x()[2] ) ); + CHECK_THAT( 3., WithinRel( result.x()[3] ) ); + CHECK_THAT( 4., WithinRel( result.x()[4] ) ); + CHECK_THAT( 4., WithinRel( result.y()[0] ) ); + CHECK_THAT( 3., WithinRel( result.y()[1] ) ); + CHECK_THAT( 2., WithinRel( result.y()[2] ) ); + CHECK_THAT( 0., WithinRel( result.y()[3] ) ); + CHECK_THAT( -2., WithinRel( result.y()[4] ) ); + CHECK( 1 == result.boundaries()[0] ); + CHECK( 4 == result.boundaries()[1] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[0] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[1] ); + CHECK( true == std::holds_alternative< IntervalDomain< double > >( result.domain() ) ); + // this will add a second point at the lower end point result = chunk + small; @@ -539,12 +631,6 @@ SCENARIO( "InterpolationTable" ) { CHECK( 4 == result.boundaries()[1] ); CHECK( InterpolationType::LinearLinear == result.interpolants()[0] ); CHECK( true == std::holds_alternative< IntervalDomain< double > >( result.domain() ) ); - - // the threshold table starts with a non-zero value - CHECK_THROWS( chunk += nonzerothreshold ); - CHECK_THROWS( chunk -= nonzerothreshold ); - CHECK_THROWS( result = chunk + nonzerothreshold ); - CHECK_THROWS( result = chunk - nonzerothreshold ); } // THEN THEN( "an InterpolationTable can be linearised" ) { @@ -687,7 +773,7 @@ SCENARIO( "InterpolationTable" ) { InterpolationTable< double > result( { 1., 4. }, { 0., 0. } ); InterpolationTable< double > same( { 1., 4. }, { 0., 3. } ); InterpolationTable< double > threshold( { 2., 4. }, { 0., 2. } ); - InterpolationTable< double > nonzerothreshold( { 2., 4. }, { 1., 2. } ); + InterpolationTable< double > nonzerothreshold( { 3., 4. }, { 1., 2. } ); InterpolationTable< double > small( { 1., 3. }, { 0., 2. } ); chunk += 2.; @@ -1130,6 +1216,106 @@ SCENARIO( "InterpolationTable" ) { CHECK( InterpolationType::LinearLinear == result.interpolants()[1] ); CHECK( true == std::holds_alternative< IntervalDomain< double > >( result.domain() ) ); + chunk += nonzerothreshold; + + CHECK( 6 == chunk.x().size() ); + CHECK( 6 == chunk.y().size() ); + CHECK( 3 == chunk.boundaries().size() ); + CHECK( 3 == 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( 3., WithinRel( chunk.x()[4] ) ); + CHECK_THAT( 4., WithinRel( chunk.x()[5] ) ); + 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( 4., WithinRel( chunk.y()[4] ) ); + CHECK_THAT( 4., WithinRel( chunk.y()[5] ) ); + CHECK( 1 == chunk.boundaries()[0] ); + CHECK( 3 == chunk.boundaries()[1] ); + CHECK( 5 == chunk.boundaries()[2] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[0] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[1] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[2] ); + CHECK( true == std::holds_alternative< IntervalDomain< double > >( chunk.domain() ) ); + + chunk -= nonzerothreshold; + + 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] ); + CHECK( 4 == chunk.boundaries()[1] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[0] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[1] ); + CHECK( true == std::holds_alternative< IntervalDomain< double > >( chunk.domain() ) ); + + result = chunk + nonzerothreshold; + + CHECK( 6 == result.x().size() ); + CHECK( 6 == result.y().size() ); + CHECK( 3 == result.boundaries().size() ); + CHECK( 3 == result.interpolants().size() ); + CHECK_THAT( 1., WithinRel( result.x()[0] ) ); + CHECK_THAT( 2., WithinRel( result.x()[1] ) ); + CHECK_THAT( 2., WithinRel( result.x()[2] ) ); + CHECK_THAT( 3., WithinRel( result.x()[3] ) ); + CHECK_THAT( 3., WithinRel( result.x()[4] ) ); + CHECK_THAT( 4., WithinRel( result.x()[5] ) ); + CHECK_THAT( 4., WithinRel( result.y()[0] ) ); + CHECK_THAT( 3., WithinRel( result.y()[1] ) ); + CHECK_THAT( 4., WithinRel( result.y()[2] ) ); + CHECK_THAT( 3., WithinRel( result.y()[3] ) ); + CHECK_THAT( 4., WithinRel( result.y()[4] ) ); + CHECK_THAT( 4., WithinRel( result.y()[5] ) ); + CHECK( 1 == result.boundaries()[0] ); + CHECK( 3 == result.boundaries()[1] ); + CHECK( 5 == result.boundaries()[2] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[0] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[1] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[2] ); + CHECK( true == std::holds_alternative< IntervalDomain< double > >( result.domain() ) ); + + result = chunk - nonzerothreshold; + + CHECK( 6 == result.x().size() ); + CHECK( 6 == result.y().size() ); + CHECK( 3 == result.boundaries().size() ); + CHECK( 3 == result.interpolants().size() ); + CHECK_THAT( 1., WithinRel( result.x()[0] ) ); + CHECK_THAT( 2., WithinRel( result.x()[1] ) ); + CHECK_THAT( 2., WithinRel( result.x()[2] ) ); + CHECK_THAT( 3., WithinRel( result.x()[3] ) ); + CHECK_THAT( 3., WithinRel( result.x()[4] ) ); + CHECK_THAT( 4., WithinRel( result.x()[5] ) ); + CHECK_THAT( 4., WithinRel( result.y()[0] ) ); + CHECK_THAT( 3., WithinRel( result.y()[1] ) ); + CHECK_THAT( 4., WithinRel( result.y()[2] ) ); + CHECK_THAT( 3., WithinRel( result.y()[3] ) ); + CHECK_THAT( 2., WithinRel( result.y()[4] ) ); + CHECK_THAT( 0., WithinRel( result.y()[5] ) ); + CHECK( 1 == result.boundaries()[0] ); + CHECK( 3 == result.boundaries()[1] ); + CHECK( 5 == result.boundaries()[2] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[0] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[1] ); + CHECK( InterpolationType::LinearLinear == result.interpolants()[2] ); + CHECK( true == std::holds_alternative< IntervalDomain< double > >( result.domain() ) ); + // this will add a second point at the lower end point result = chunk + small; @@ -1183,12 +1369,6 @@ SCENARIO( "InterpolationTable" ) { CHECK( InterpolationType::LinearLinear == result.interpolants()[1] ); CHECK( InterpolationType::LinearLinear == result.interpolants()[2] ); CHECK( true == std::holds_alternative< IntervalDomain< double > >( result.domain() ) ); - - // the threshold table starts with a non-zero value - CHECK_THROWS( chunk += nonzerothreshold ); - CHECK_THROWS( chunk -= nonzerothreshold ); - CHECK_THROWS( result = chunk + nonzerothreshold ); - CHECK_THROWS( result = chunk - nonzerothreshold ); } // THEN } // WHEN } // GIVEN @@ -1569,6 +1749,17 @@ SCENARIO( "InterpolationTable" ) { } // THEN } // WHEN + WHEN( "the x grid has a jump at the beginning" ) { + + std::vector< double > x = { 1., 1., 3., 4. }; + std::vector< double > y = { 4., 3., 1., 4. }; + + THEN( "an exception is thrown" ) { + + CHECK_THROWS( InterpolationTable< double >( std::move( x ), std::move( y ) ) ); + } // THEN + } // WHEN + WHEN( "the x grid has a jump at the end" ) { std::vector< double > x = { 1., 2., 4., 4. }; diff --git a/src/scion/unionisation/unionise.hpp b/src/scion/unionisation/unionise.hpp index 99def99..39f9308 100644 --- a/src/scion/unionisation/unionise.hpp +++ b/src/scion/unionisation/unionise.hpp @@ -13,13 +13,9 @@ namespace unionisation { /** * @brief Unionise two grids and preserve duplicate points that appear in each * - * If the grids do not have the same end point, a duplicate point is inserted - * into the grid corresponding to the lowest end point (unless that is - * already a duplicate point). - * - * No special treatment is performed when the grids do not have the same - * starting point. We assume that the tabulated value at the higher starting - * point will be zero (we can enforce this behaviour when reading the data). + * If the grids do not have the same begin and/or end point, a duplicate point + * is inserted into the grid corresponding to the highest begining and/or + * lowest end point (unless those is already a duplicate point). * * @param first the first grid (assumed to be sorted) * @param second the second grid (assumed to be sorted) @@ -35,6 +31,17 @@ namespace unionisation { grid.begin() ); grid.erase( end, grid.end() ); + // special case: the begin points are not the same + if ( first.front() != second.front() ) { + + X x = first.front() < second.front() ? second.front() : first.front(); + auto iter = std::lower_bound( grid.begin(), grid.end(), x ); + if ( *std::next( iter ) != x ) { + + grid.insert( iter, x ); + } + } + // special case: the end points are not the same if ( first.back() != second.back() ) {