Table Integration

.gitignore

@@ -1 +1,2 @@
 bin
+build

CMakeLists.txt

@@ -62,7 +62,6 @@ target_link_libraries( interpolation
     INTERFACE dimwits
     )
 
-
 #######################################################################
 # Top-level Only
 #######################################################################

README.md

@@ -43,6 +43,7 @@ int main(){
   std::cout << myTable( 2.5 ) << std::endl;
   std::cout << myTable.domainMin() << ' ' << myTable.domainMax() << std::endl;
   std::cout << myTable.tableMin() << ' ' << myTable.tableMax() << std::endl;
+  std::cout << myTable.integrate(1., 3.) << std::endl;
 }
 ```
 

src/interpolation/Interpolant/Histogram.hpp

@@ -5,6 +5,11 @@ struct Histogram : public Interpolant {
   template< typename Yarg, typename X, typename Y >
   static auto invert( Yarg&&, X&& xLeft, X&&, Y&&, Y&& ){ return xLeft; }
 
+  template< typename Xarg, typename X, typename Y >
+  static auto integrate(Xarg&& xLow, Xarg&& xHi, X&&, X&&, Y&& yLeft, Y&&) {
+    return yLeft * (xHi - xLow);
+  }
+
   template< typename Range >
   static void verifyXGridAssumptions( Range&& range ){
     auto it = ranges::adjacent_find( range );

src/interpolation/Interpolant/Histogram/test/Histogram.test.cpp

@@ -22,6 +22,16 @@ SCENARIO("The Histogram interpolant computes the correct inversion",
   REQUIRE( xLeft == Histogram::invert( 15.0, xLeft, xRight, yLeft, yRight ) );
 }
 
+SCENARIO("The Histogram interpolant computes the correct integral",
+         "[interpolant], [Histogram]"){
+  double xLeft = 0.0, xRight = 1.0, yLeft = 10.0, yRight = 20.0;
+  double xLow = 0., xHi = 0.5;
+  REQUIRE( 0.5*yLeft == Histogram::integrate( xLow, xHi, xLeft, xRight, yLeft, yRight ) );
+  xLow = 0.25;
+  xHi = xRight;
+  REQUIRE( 0.75*yLeft == Histogram::integrate( xLow, xHi, xLeft, xRight, yLeft, yRight ) );
+}
+
 SCENARIO("The Histogram interpolant is compatible with units",
          "[interpolant], [Histogram]"){
   auto xLeft = 0.0 * electronVolts;
@@ -32,6 +42,10 @@ SCENARIO("The Histogram interpolant is compatible with units",
 	   Histogram::apply( 0.5 * electronVolts, xLeft, xRight, yLeft, yRight ) );
   REQUIRE( xLeft ==
 	   Histogram::invert( 15.0 * barns, xLeft, xRight, yLeft, yRight ) );  
+  auto xLow = 0. * electronVolts;
+  auto xHi = 0.5 * electronVolts;
+  REQUIRE( 5. * electronVolts * barn ==
+           Histogram::integrate(xLow, xHi, xLeft, xRight, yLeft, yRight));
 }
 
 SCENARIO("The Histogram rejects grids with redundant x-values",

src/interpolation/Interpolant/LinearLinear.hpp

@@ -8,4 +8,10 @@ struct LinearLinear : public Interpolant {
   static auto invert( Yarg&& y, X&& xLeft, X&& xRight, Y&& yLeft, Y&& yRight ){
     return xLeft + ( xRight - xLeft ) * ( y - yLeft ) / ( yRight - yLeft );
   }
+
+  template< typename Xarg, typename X, typename Y >
+  static auto integrate(Xarg&& xLow, Xarg&& xHi, X&& xLeft, X&& xRight, Y&& yLeft, Y&& yRight) {
+    const auto m = (yRight - yLeft) / (xRight - xLeft);
+    return 0.5*((xHi*xHi) - (xLow*xLow))*m + (xHi - xLow)*(yLeft - m*xLeft);
+  }
 };

src/interpolation/Interpolant/LinearLinear/test/LinearLinear.test.cpp

@@ -30,6 +30,18 @@ SCENARIO("The LinearLinear interpolant computes the correct inversion",
   REQUIRE( xRight == LinearLinear::invert( yRight, xLeft, xRight, yLeft, yRight ) );
 }
 
+SCENARIO("The LinearLinear interpolant computes the correct integral",
+         "[interpolant], [LinearLinear]"){
+  double xLeft = 0.0, xRight = 1.0;
+  auto f = []( double x ){ return 10.0 * x + 10; };
+  double yLeft = f( xLeft ), yRight = f( xRight );
+  auto F = [](double x){ return 5.*(x*x) + (10.*x); };
+
+  REQUIRE( F(0.5) == LinearLinear::integrate( 0., 0.5, xLeft, xRight, yLeft, yRight ) );
+  REQUIRE( F(0.75) == LinearLinear::integrate( 0., 0.75, xLeft, xRight, yLeft, yRight ) );
+  REQUIRE( F(0.75) - F(0.25) == LinearLinear::integrate( 0.25, 0.75, xLeft, xRight, yLeft, yRight ) );
+}
+
 SCENARIO("The LinearLinear interpolant is compatible with units",
          "[interpolant], [Histogram]"){
   auto xLeft = 0.0 * electronVolts;
@@ -40,4 +52,8 @@ SCENARIO("The LinearLinear interpolant is compatible with units",
 	   LinearLinear::apply( 0.5 * electronVolts, xLeft, xRight, yLeft, yRight ) );
   REQUIRE( 0.5 * electronVolts ==
 	   LinearLinear::invert( 15.0 * barns, xLeft, xRight, yLeft, yRight ) );  
+  auto xLow = 0. * electronVolts;
+  auto xHi  = 1. * electronVolts;
+  REQUIRE(15.*electronVolts*barns ==
+          LinearLinear::integrate(xLow, xHi, xLeft, xRight, yLeft, yRight));
 }

src/interpolation/Interpolant/LinearLogarithmic.hpp

@@ -16,6 +16,17 @@ struct LinearLogarithmic : public Interpolant {
                         ( ( safe(y) - yLeft ) / ( yRight - yLeft ) ) ) : xLeft;
   }
 
+  template< typename Xarg, typename X, typename Y >
+  static auto integrate( Xarg&& xLow, Xarg&& xHi, X&& xLeft, X&& xRight, Y&& yLeft, Y&& yRight ){
+    using safe = std::decay_t<X>;
+    const auto numerator_hi = xHi * ((yRight - yLeft) * std::log(safe(xHi) / xLeft) +
+                              yLeft * std::log(xRight / xLeft) + yLeft - yRight);
+    const auto numerator_low = xLow * ((yRight - yLeft) * std::log(safe(xLow) / xLeft) +
+                               yLeft * std::log(xRight / xLeft) + yLeft - yRight);
+    const auto denominator = std::log(xRight / xLeft);
+    return (numerator_hi / denominator) - (numerator_low / denominator);
+  }
+
   template< typename Range >
   static void verifyXGridAssumptions( Range&& range ){
     const auto zero = 0 * range.front();

src/interpolation/Interpolant/LinearLogarithmic/test/CMakeLists.txt

@@ -3,6 +3,7 @@ add_executable( interpolation.Interpolant.LinearLogarithmic.test
                 verifyXGridAssumptions.test.cpp
                 invert.test.cpp
                 LinearLogarithmic.test.cpp
+                integration.test.cpp
                 interpolate.test.cpp )
 target_compile_options( interpolation.Interpolant.LinearLogarithmic.test PRIVATE ${${PREFIX}_common_flags}
 $<$<BOOL:${strict}>:${${PREFIX}_strict_flags}>$<$<CONFIG:DEBUG>:

src/interpolation/Interpolant/LinearLogarithmic/test/integration.test.cpp

@@ -0,0 +1,91 @@
+#include "catch.hpp"
+#include "interpolation.hpp"
+
+#include "dimwits.hpp"
+
+#include "range/v3/view/take.hpp"
+
+using namespace njoy::interpolation;
+using namespace dimwits;
+
+SCENARIO("The LinearLogarithmic interpolant computes the correct integral",
+         "[interpolant], [LinearLogarithmic]"){
+  std::vector<double> xValues{ 1E-2, 1, 10 };
+  auto f1 = []( double x ){ return 13.2 * std::log(3.2*x); };
+  auto f2 = []( double x ){ return -5.6 * std::log(2.4*x); };
+  auto F1 = [](double x){return 13.2*x*(std::log(3.2*x)-1.);};
+  auto F2 = [](double x){return -5.6*x*(std::log(2.4*x)-1.);};
+
+  auto avg =
+    []( double xLeft, double xRight ){ return 0.5 * ( xLeft + xRight ); };
+
+  auto excessiveError =
+    []( double reference, double trial ){
+    auto error = std::abs( (trial - reference)
+			   / ( ( reference != 0 ) ? reference : 1.0 ) );
+    return error > 1E-15; };
+
+  auto iterator = xValues.begin();
+  auto last = std::prev( xValues.end() );
+  do {
+    double xLeft = *iterator, xRight = *( ++iterator );
+    double y1Left = f1( xLeft ), y1Right = f1( xRight );
+    double y2Left = f2( xLeft ), y2Right = f2( xRight );
+    double xBar = avg( xLeft, xRight );
+    double xLow = xLeft, xHi = xRight;
+
+    REQUIRE( not excessiveError( F1(xBar) - F1(xLow),
+                                 LinearLogarithmic::integrate
+                                 ( xLow, xBar, xLeft, xRight, y1Left, y1Right ) ) );
+    REQUIRE( not excessiveError( F2(xBar) - F2(xLow),
+                                 LinearLogarithmic::integrate
+                                 ( xLow, xBar, xLeft, xRight, y2Left, y2Right ) ) );
+
+    REQUIRE( not excessiveError( F1(xHi) - F1(xBar),
+                                 LinearLogarithmic::integrate
+                                 ( xBar, xHi, xLeft, xRight, y1Left, y1Right ) ) );
+    REQUIRE( not excessiveError( F2(xHi) - F2(xBar),
+                                 LinearLogarithmic::integrate
+                                 ( xBar, xHi, xLeft, xRight, y2Left, y2Right ) ) );
+  } while( iterator != last );
+}
+
+SCENARIO("The LinearLogarithmic interpolant integration is compatible with units",
+         "[interpolant], [LinearLogarithmic]"){
+  std::vector< double > xValues{ 1E-2, 1, 10 };
+  auto f1 = []( auto x ){ return 13.2 * std::log(3.2 * x.value) * barn; };
+  auto f2 = []( auto x ){ return -5.6 * std::log(2.4 * x.value) * barn; };
+  auto F1 = []( auto x ){return 13.2*x.value*(std::log(3.2*x.value)-1.) * barn * electronVolts;};
+  auto F2 = []( auto x ){return -5.6*x.value*(std::log(2.4*x.value)-1.) * barn * electronVolts;};
+
+  auto avg = []( auto xLeft, auto xRight ) { return 0.5 * ( xLeft + xRight ); };
+
+  auto excessiveError =
+    []( auto reference, auto trial ){
+    auto error = std::abs( (trial - reference).value
+			   / ( ( reference.value != 0 ) ?
+			       reference.value : 1.0 ) );
+    return error > 1E-15; };
+
+  auto units = xValues |
+    ranges::view::take( xValues.size() - 1 ) |
+    ranges::view::transform( []( auto arg ){ return arg * electronVolts; } );
+
+  auto iterator = units.begin();
+  auto last = units.end();
+  do {
+    auto xLeft = *iterator, xRight = *( ++iterator );
+    auto y1Left = f1( xLeft ), y1Right = f1( xRight );
+    auto y2Left = f2( xLeft ), y2Right = f2( xRight );
+    auto xBar = avg( xLeft, xRight );
+    auto xLow = xLeft, xHi = xRight;
+
+    REQUIRE( not excessiveError( F2(xBar) - F2(xLow),
+                                 LinearLogarithmic::integrate
+                                 ( xLow, xBar, xLeft, xRight, y2Left, y2Right ) ) );
+
+    REQUIRE( not excessiveError( F1(xHi) - F1(xBar),
+                                 LinearLogarithmic::integrate
+                                 ( xBar, xHi, xLeft, xRight, y1Left, y1Right ) ) );
+  } while( iterator != last );
+}

src/interpolation/Interpolant/LinearLogarithmic/test/interpolate.test.cpp

@@ -52,7 +52,9 @@ SCENARIO("The LinearLogarithmic interpolant computes is compatible with units",
   std::vector< double > xValues{ 1E-2, 1, 10 };
   auto f1 = []( auto x ){ return 13.2 * std::log(3.2 * x.value) * barn; };
   auto f2 = []( auto x ){ return -5.6 * std::log(2.4 * x.value) * barn; };
-
+  auto F1 = []( auto x ){return 13.2*x.value*(std::log(3.2*x.value)-1.) * barn * electronVolts;};
+  auto F2 = []( auto x ){return -5.6*x.value*(std::log(2.4*x.value)-1.) * barn * electronVolts;};
+
   auto avg = []( auto xLeft, auto xRight ) { return 0.5 * ( xLeft + xRight ); };
 
   auto excessiveError =
@@ -73,6 +75,7 @@ SCENARIO("The LinearLogarithmic interpolant computes is compatible with units",
     auto y1Left = f1( xLeft ), y1Right = f1( xRight );
     auto y2Left = f2( xLeft ), y2Right = f2( xRight );
     auto xBar = avg( xLeft, xRight );
+    auto xLow = xLeft, xHi = xRight;
 
     REQUIRE( not excessiveError( y1Left,
                                  LinearLogarithmic::apply
@@ -87,5 +90,6 @@ SCENARIO("The LinearLogarithmic interpolant computes is compatible with units",
     REQUIRE( not excessiveError( f2( xBar ),
                                  LinearLogarithmic::apply
                                  ( xBar, xLeft, xRight, y2Left, y2Right ) ) );
+
   } while( iterator != last );
 }

src/interpolation/Interpolant/LogarithmicLinear.hpp

@@ -17,6 +17,18 @@ struct LogarithmicLinear : public Interpolant {
       xLeft;
   }
 
+  template< typename Xarg, typename X, typename Y >
+  static auto integrate( Xarg&& xLow, Xarg&& xHi, X&& xLeft, X&& xRight, Y&& yLeft, Y&& yRight ){
+    using safe = std::decay_t<X>;
+    const auto base = yRight / yLeft;
+    const auto denominator = std::log(base);
+    const auto coefficient = yLeft * (xRight - xLeft);
+    const auto exponent_hi = (xLeft - safe(xHi)) / (xLeft - xRight);
+    const auto exponent_low = (xLeft - safe(xLow)) / (xLeft - xRight);
+    return (coefficient / denominator) * 
+           (std::pow(base, exponent_hi) - std::pow(base, exponent_low));
+  }
+
   template< typename Range >
   static void verifyYGridAssumptions( Range&& range ){
     auto it = Interpolant::findChangeOfSign( range );

src/interpolation/Interpolant/LogarithmicLinear/test/CMakeLists.txt

@@ -2,6 +2,7 @@
 add_executable( interpolation.Interpolant.LogarithmicLinear.test
                 apply.test.cpp
                 LogarithmicLinear.test.cpp
+                integration.test.cpp
                 invert.test.cpp )
 target_compile_options( interpolation.Interpolant.LogarithmicLinear.test PRIVATE ${${PREFIX}_common_flags}
 $<$<BOOL:${strict}>:${${PREFIX}_strict_flags}>$<$<CONFIG:DEBUG>:

src/interpolation/Interpolant/LogarithmicLinear/test/integration.test.cpp

@@ -0,0 +1,91 @@
+#include "catch.hpp"
+#include "interpolation.hpp"
+
+#include "dimwits.hpp"
+
+#include "range/v3/view/take.hpp"
+
+using namespace njoy::interpolation;
+
+SCENARIO("The LogarithmicLinear interpolant computes the correct integral",
+         "[interpolant], [LogarithmicLinear]"){
+  std::vector<double> xValues {10, 1, 1E-1};
+  auto f1 = []( double x ){ return 13.2 * pow(3.14, 0.002*x); };
+  auto f2 = []( double x ){ return -13.2 * pow(3.14, 0.002*x); };
+  auto F1 = []( double xl, double xh ){ return (13.2/(0.002*std::log(3.14)))*(pow(3.14, 0.002*xh) - pow(3.14, 0.002*xl));};
+  auto F2 = []( double xl, double xh ){ return (-13.2/(0.002*std::log(3.14)))*(pow(3.14, 0.002*xh) - pow(3.14, 0.002*xl));};
+  //auto F2 = []( double x ){ return (-13.2*pow(3.14, 0.002*x))/(0.002*std::log(3.14));};
+  auto avg = []( double xLeft, double xRight ){ return 0.5 * ( xLeft + xRight ); };
+
+  auto excessiveError =
+    []( double reference, double trial ){
+    auto error = std::abs( (trial - reference)
+                           / ( ( reference != 0 ) ? reference : 1.0 ) );
+    return error > 1E-13; };
+
+  auto iterator = xValues.begin();
+  auto last = std::prev( xValues.end() );
+  do {
+    double xLeft = *iterator, xRight = *( ++iterator );
+    double y1Left = f1( xLeft ), y1Right = f1( xRight );
+    double y2Left = f2( xLeft ), y2Right = f2( xRight );
+    double xBar = avg( xLeft, xRight );
+    double xLow = xLeft, xHi = xRight;
+
+    REQUIRE( not excessiveError( F1(xLow, xBar),
+                                 LogarithmicLinear::integrate
+                                 ( xLow, xBar, xLeft, xRight, y1Left, y1Right ) ) );
+    REQUIRE( not excessiveError( F2(xLow, xBar),
+                                 LogarithmicLinear::integrate
+                                 ( xLow, xBar, xLeft, xRight, y2Left, y2Right ) ) );
+
+    REQUIRE( not excessiveError( F1(xBar, xHi),
+                                 LogarithmicLinear::integrate
+                                 ( xBar, xHi, xLeft, xRight, y1Left, y1Right ) ) );
+    REQUIRE( not excessiveError( F2(xBar, xHi),
+                                 LogarithmicLinear::integrate
+                                 ( xBar, xHi, xLeft, xRight, y2Left, y2Right ) ) );
+  } while( iterator != last );
+}
+
+using namespace dimwits;
+
+SCENARIO("The LogarithmicLinear integration is compatible with units",
+         "[interpolant], [LinearLogarithmic]"){
+  std::vector< double > xValues{ 1E-2, 1, 10 };
+  auto f1 = []( auto x ){ return 13.2 * pow(3.14, 0.002*x.value) * barn; };
+  auto f2 = []( auto x ){ return -13.2 * pow(3.14, 0.002*x.value) * barn; };
+  auto F1 = []( auto x ){ return (13.2*pow(3.14, 0.002*x.value))/(0.002*std::log(3.14))*barn*electronVolts;};
+  auto F2 = []( auto x ){ return (-13.2*pow(3.14, 0.002*x.value))/(0.002*std::log(3.14))*barn*electronVolts;};
+  auto avg = []( auto xLeft, auto xRight ) { return 0.5 * ( xLeft + xRight ); };
+
+  auto excessiveError =
+    []( auto reference, auto trial ){
+    auto error = std::abs( (reference - trial).value
+			   / ( ( reference.value != 0 ) ?
+			       reference.value : 1.0 ) );
+    return error > 1E-13; };
+
+  auto units = xValues |
+    ranges::view::take( xValues.size() - 1 ) |
+    ranges::view::transform( []( auto arg ){ return arg * electronVolts; } );
+
+  auto iterator = units.begin();
+  auto last = units.end();
+  do {
+    auto xLeft = *iterator, xRight = *( ++iterator );
+    auto y1Left = f1( xLeft ), y1Right = f1( xRight );
+    auto y2Left = f2( xLeft ), y2Right = f2( xRight );
+    auto xBar = avg( xLeft, xRight );
+    auto xLow = xLeft, xHi = xRight;
+
+    REQUIRE( not excessiveError( F2(xBar) - F2(xLow),
+                                 LogarithmicLinear::integrate
+                                 ( xLow, xBar, xLeft, xRight, y2Left, y2Right ) ) );
+
+    REQUIRE( not excessiveError( F1(xHi) - F1(xBar),
+                                 LogarithmicLinear::integrate
+                                 ( xBar, xHi, xLeft, xRight, y1Left, y1Right ) ) );
+
+  } while( iterator != last );
+}

src/interpolation/Interpolant/LogarithmicLogarithmic.hpp

@@ -17,6 +17,20 @@ struct LogarithmicLogarithmic : public Interpolant {
                           / std::log( yRight * inverseYLeft );
     return xLeft * std::pow(xRight / xLeft, logRatio);
   }
+
+  template< typename Xarg, typename X, typename Y >
+  static auto integrate( Xarg&& xLow, Xarg&& xHi, X&& xLeft, X&& xRight, Y&& yLeft, Y&& yRight ){
+    using safe = std::decay_t<X>;
+    const auto yRatio = yRight / yLeft;
+    const auto xRatio = xRight / xLeft;
+    const auto log_yRatio = std::log(yRatio);
+    const auto log_xRatio = std::log(xRatio);
+    const auto exponent = log_yRatio / log_xRatio;
+    const auto denominator = exponent + 1.0;
+    if(std::abs(denominator) <= 1.E-12) {return (yLeft*xLeft) * std::log(safe(xHi)/safe(xLow));}
+    return (yLeft / denominator) * (safe(xHi) * std::pow(safe(xHi) / xLeft, exponent) -
+                                safe(xLow) * std::pow(safe(xLow) / xLeft, exponent));
+  }
 
   template< typename Range >
   static void verifyXGridAssumptions( Range&& range ){

src/interpolation/Interpolant/LogarithmicLogarithmic/test/CMakeLists.txt

@@ -1,6 +1,7 @@
 
 add_executable( interpolation.Interpolant.LogarithmicLogarithmic.test
                 verifyXGridAssumptions.test.cpp
+                integration.test.cpp
                 invert.test.cpp
                 LogarithmicLogarithmic.test.cpp
                 verifyYGridAssumptions.test.cpp