Updating operation function

njoy · Nov 30, 2023 · 7edf07e · 7edf07e
1 parent 92a85e1
commit 7edf07e
Show file tree

Hide file tree

Showing 3 changed files with 19 additions and 20 deletions.
diff --git 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
@@ -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
@@ -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 ) );
     }