diff --git a/.gitignore b/.gitignore index ba077a4..c745919 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ bin +build diff --git a/CMakeLists.txt b/CMakeLists.txt index fc8abc7..30a1781 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -62,7 +62,6 @@ target_link_libraries( interpolation INTERFACE dimwits ) - ####################################################################### # Top-level Only ####################################################################### diff --git a/README.md b/README.md index 5b70a10..f9f5f7c 100644 --- a/README.md +++ b/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; } ``` diff --git a/src/interpolation/Interpolant/Histogram.hpp b/src/interpolation/Interpolant/Histogram.hpp index f6f3ee2..7077d91 100644 --- a/src/interpolation/Interpolant/Histogram.hpp +++ b/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 ); diff --git a/src/interpolation/Interpolant/Histogram/test/Histogram.test.cpp b/src/interpolation/Interpolant/Histogram/test/Histogram.test.cpp index 97734e4..c3812d4 100644 --- a/src/interpolation/Interpolant/Histogram/test/Histogram.test.cpp +++ b/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", diff --git a/src/interpolation/Interpolant/LinearLinear.hpp b/src/interpolation/Interpolant/LinearLinear.hpp index 2dbf325..781ac60 100644 --- a/src/interpolation/Interpolant/LinearLinear.hpp +++ b/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); + } }; diff --git a/src/interpolation/Interpolant/LinearLinear/test/LinearLinear.test.cpp b/src/interpolation/Interpolant/LinearLinear/test/LinearLinear.test.cpp index bfda716..f92564b 100644 --- a/src/interpolation/Interpolant/LinearLinear/test/LinearLinear.test.cpp +++ b/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)); } diff --git a/src/interpolation/Interpolant/LinearLogarithmic.hpp b/src/interpolation/Interpolant/LinearLogarithmic.hpp index 0c473ba..61a3951 100644 --- a/src/interpolation/Interpolant/LinearLogarithmic.hpp +++ b/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; + 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(); diff --git a/src/interpolation/Interpolant/LinearLogarithmic/test/CMakeLists.txt b/src/interpolation/Interpolant/LinearLogarithmic/test/CMakeLists.txt index 6f7dd47..9393fd7 100644 --- a/src/interpolation/Interpolant/LinearLogarithmic/test/CMakeLists.txt +++ b/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} $<$:${${PREFIX}_strict_flags}>$<$: diff --git a/src/interpolation/Interpolant/LinearLogarithmic/test/integration.test.cpp b/src/interpolation/Interpolant/LinearLogarithmic/test/integration.test.cpp new file mode 100644 index 0000000..482b9b4 --- /dev/null +++ b/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 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 ); +} \ No newline at end of file diff --git a/src/interpolation/Interpolant/LinearLogarithmic/test/interpolate.test.cpp b/src/interpolation/Interpolant/LinearLogarithmic/test/interpolate.test.cpp index ac1a3b2..8d67d33 100644 --- a/src/interpolation/Interpolant/LinearLogarithmic/test/interpolate.test.cpp +++ b/src/interpolation/Interpolant/LinearLogarithmic/test/interpolate.test.cpp @@ -52,7 +52,7 @@ 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 avg = []( auto xLeft, auto xRight ) { return 0.5 * ( xLeft + xRight ); }; auto excessiveError = @@ -73,6 +73,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 +88,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 ); } diff --git a/src/interpolation/Interpolant/LogarithmicLinear.hpp b/src/interpolation/Interpolant/LogarithmicLinear.hpp index f031ee9..8884f70 100644 --- a/src/interpolation/Interpolant/LogarithmicLinear.hpp +++ b/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; + 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 ); diff --git a/src/interpolation/Interpolant/LogarithmicLinear/test/CMakeLists.txt b/src/interpolation/Interpolant/LogarithmicLinear/test/CMakeLists.txt index 47079fd..df434db 100644 --- a/src/interpolation/Interpolant/LogarithmicLinear/test/CMakeLists.txt +++ b/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} $<$:${${PREFIX}_strict_flags}>$<$: diff --git a/src/interpolation/Interpolant/LogarithmicLinear/test/integration.test.cpp b/src/interpolation/Interpolant/LogarithmicLinear/test/integration.test.cpp new file mode 100644 index 0000000..d96f29e --- /dev/null +++ b/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 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 ); +} \ No newline at end of file diff --git a/src/interpolation/Interpolant/LogarithmicLogarithmic.hpp b/src/interpolation/Interpolant/LogarithmicLogarithmic.hpp index bbb1c67..2ffe987 100644 --- a/src/interpolation/Interpolant/LogarithmicLogarithmic.hpp +++ b/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; + 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 ){ diff --git a/src/interpolation/Interpolant/LogarithmicLogarithmic/test/CMakeLists.txt b/src/interpolation/Interpolant/LogarithmicLogarithmic/test/CMakeLists.txt index aaec36f..1561070 100644 --- a/src/interpolation/Interpolant/LogarithmicLogarithmic/test/CMakeLists.txt +++ b/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 diff --git a/src/interpolation/Interpolant/LogarithmicLogarithmic/test/integration.test.cpp b/src/interpolation/Interpolant/LogarithmicLogarithmic/test/integration.test.cpp new file mode 100644 index 0000000..eb53098 --- /dev/null +++ b/src/interpolation/Interpolant/LogarithmicLogarithmic/test/integration.test.cpp @@ -0,0 +1,164 @@ +#include "catch.hpp" +#include "interpolation.hpp" + +#include "dimwits.hpp" + +#include "range/v3/view/take.hpp" + +using namespace njoy::interpolation; + +SCENARIO("The LogarithmicLogarithmic interpolant computes the correct integral", + "[interpolant], [LogarithmicLogarithmic]"){ + std::vector xValues {100, 10, 1, 1E-1, 1E-2}; + + // Function Definitions + auto f1 = []( double x ){ return 10.0 * x; }; + auto f2 = []( double x ){ return 10.0 * std::pow(x,2); }; + auto f3 = []( double x ){ return 5.0 * std::pow(x,2); }; + auto f4 = []( double x ){ return 5.0 * std::pow(x,-1.0); }; + auto f5 = []( double x ){ return 5.0 * std::pow(x,-2.0); }; + + // Integral Definitions + auto F1 = [](double xLow, double xHi){return 5. * (std::pow(xHi,2.) - std::pow(xLow,2.));}; + auto F2 = [](double xLow, double xHi){return (10./3.)*(std::pow(xHi,3.) - std::pow(xLow,3.));}; + auto F3 = [](double xLow, double xHi){return (5./3.)*(std::pow(xHi,3.) - std::pow(xLow,3.));}; + auto F4 = [](double xLow, double xHi){return 5.*(std::log(xHi) - std::log(xLow));}; + auto F5 = [](double xLow, double xHi){return 5.*((1./xLow) - (1./xHi));}; + + 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 y3Left = f3( xLeft ), y3Right = f3( xRight ); + double y4Left = f4( xLeft ), y4Right = f4( xRight ); + double y5Left = f5( xLeft ), y5Right = f5( xRight ); + double xBar = avg( xLeft, xRight ); + double xLow = xLeft, xHi = xRight; + + REQUIRE( not excessiveError( F1(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y1Left, y1Right ) ) ); + REQUIRE( not excessiveError( F2(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y2Left, y2Right ) ) ); + REQUIRE( not excessiveError( F3(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y3Left, y3Right ) ) ); + REQUIRE( not excessiveError( F4(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y4Left, y4Right ) ) ); + REQUIRE( not excessiveError( F5(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y5Left, y5Right ) ) ); + + + REQUIRE( not excessiveError( F1(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y1Left, y1Right ) ) ); + REQUIRE( not excessiveError( F2(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y2Left, y2Right ) ) ); + REQUIRE( not excessiveError( F3(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y3Left, y3Right ) ) ); + REQUIRE( not excessiveError( F4(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y4Left, y4Right ) ) ); + REQUIRE( not excessiveError( F5(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y5Left, y5Right ) ) ); + } while( iterator != last ); +} + +using namespace dimwits; + +SCENARIO("The LogarithmicLogarithmic interpolant computes the correct integral with units", + "[interpolant], [LogarithmicLogarithmic]"){ + std::vector xValues {100, 10, 1, 1E-1, 1E-2}; + + // Function Definitions + auto f1 = []( auto x ){ return 10.0 * x.value * barns; }; + auto f2 = []( auto x ){ return 10.0 * std::pow(x.value,2) * barns; }; + auto f3 = []( auto x ){ return 5.0 * std::pow(x.value,2) * barns; }; + auto f4 = []( auto x ){ return 5.0 * std::pow(x.value,-1.0) * barns; }; + auto f5 = []( auto x ){ return 5.0 * std::pow(x.value,-2.0) * barns; }; + + // Integral Definitions + auto F1 = [](auto xLow, auto xHi){return 5. * (std::pow(xHi.value,2.) - std::pow(xLow.value,2.)) * barns * electronVolts;}; + auto F2 = [](auto xLow, auto xHi){return (10./3.)*(std::pow(xHi.value,3.) - std::pow(xLow.value,3.)) * barns * electronVolts;}; + auto F3 = [](auto xLow, auto xHi){return (5./3.)*(std::pow(xHi.value,3.) - std::pow(xLow.value,3.)) * barns * electronVolts;}; + auto F4 = [](auto xLow, auto xHi){return 5.*(std::log(xHi.value) - std::log(xLow.value)) * barns * electronVolts;}; + auto F5 = [](auto xLow, auto xHi){return 5.*((1./xLow.value) - (1./xHi.value)) * barns * 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 y3Left = f3( xLeft ), y3Right = f3( xRight ); + auto y4Left = f4( xLeft ), y4Right = f4( xRight ); + auto y5Left = f5( xLeft ), y5Right = f5( xRight ); + auto xBar = avg( xLeft, xRight ); + auto xLow = xLeft, xHi = xRight; + + REQUIRE( not excessiveError( F1(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y1Left, y1Right ) ) ); + REQUIRE( not excessiveError( F2(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y2Left, y2Right ) ) ); + REQUIRE( not excessiveError( F3(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y3Left, y3Right ) ) ); + REQUIRE( not excessiveError( F4(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y4Left, y4Right ) ) ); + REQUIRE( not excessiveError( F5(xLow, xBar), + LogarithmicLogarithmic::integrate + ( xLow, xBar, xLeft, xRight, y5Left, y5Right ) ) ); + + + REQUIRE( not excessiveError( F1(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y1Left, y1Right ) ) ); + REQUIRE( not excessiveError( F2(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y2Left, y2Right ) ) ); + REQUIRE( not excessiveError( F3(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y3Left, y3Right ) ) ); + REQUIRE( not excessiveError( F4(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y4Left, y4Right ) ) ); + REQUIRE( not excessiveError( F5(xBar, xHi), + LogarithmicLogarithmic::integrate + ( xBar, xHi, xLeft, xRight, y5Left, y5Right ) ) ); + } while( iterator != last ); +} \ No newline at end of file diff --git a/src/interpolation/table/Type.hpp b/src/interpolation/table/Type.hpp index 3bb9ca5..f2f1867 100644 --- a/src/interpolation/table/Type.hpp +++ b/src/interpolation/table/Type.hpp @@ -40,6 +40,7 @@ class Type { static constexpr Xdata domainMax(){ return infinity(); } #include "interpolation/table/Type/src/evaluate.hpp" + #include "interpolation/table/Type/src/integrate.hpp" }; template< typename InterpolationAlgorithm, diff --git a/src/interpolation/table/Type/Table.hpp b/src/interpolation/table/Type/Table.hpp index a510294..85dc12f 100644 --- a/src/interpolation/table/Type/Table.hpp +++ b/src/interpolation/table/Type/Table.hpp @@ -41,6 +41,17 @@ class Table< table::Type< Args... >, Decorators... > : return Parent:: template evaluate( x, std::forward(args)... ); } + + template< typename Xdata > + auto integrate(const Xdata& xLow, const Xdata& xHi) const { + return Parent:: template integrate( xLow, xHi, Parent::search() ); + } + + template< typename Xdata, typename... CallArgs > + auto integrate(const Xdata& xLow, const Xdata& xHi, CallArgs&&... args) const { + return Parent:: template integrate( xLow, xHi, + std::forward(args)...); + } }; diff --git a/src/interpolation/table/Type/src/integrate.hpp b/src/interpolation/table/Type/src/integrate.hpp new file mode 100644 index 0000000..e19c2cc --- /dev/null +++ b/src/interpolation/table/Type/src/integrate.hpp @@ -0,0 +1,78 @@ +template< typename Table, typename Algorithm > +auto integrate( const Xdata& xL, const Xdata& xH, + Algorithm&& searchAlgorithm ) const { + auto xLow = xL; + auto xHi = xH; + + const bool inverted = xLow > xHi; + if(inverted) { + const auto xLowTmp = xLow; + xLow = xHi; + xHi = xLowTmp; + } + + // Integration may only be carried out over the function's valid domain + if (xLow < this->tableMin()) { + xLow = this->tableMin(); + } else if (xLow > this->tableMax()) { + xLow = this->tableMax(); + } + + if (xHi > this->tableMax()) { + xHi = this->tableMax(); + } else if (xHi < this->tableMin()) { + xHi = this->tableMin(); + } + + // Initialize to zero by doing integral over range of 0 width + auto integral = (xLow - xLow) * this->y().front(); + if (xLow == xHi) { + return integral; + } + + // Get iterator for lower bound of first interval + auto lowIt = searchAlgorithm.apply(this->x(), xLow); + if (*lowIt > xLow) { + lowIt = ranges::prev(lowIt); + } + + auto xLowLim = xLow; + auto xUppLim = xHi; + bool integrating = true; + while (integrating) { + auto hiIt = lowIt; + hiIt = ranges::next(hiIt); + + auto i = ranges::distance(this->x().begin(), lowIt); + + auto x1 = *lowIt; + auto x2 = *hiIt; + auto y1 = *ranges::next(this->y().begin(), i); + auto y2 = *ranges::next(this->y().begin(), i+1); + + if (xLowLim < x1) { + xLowLim = x1; + } + if (xUppLim > x2) { + xUppLim = x2; + } + + if(x1 != x2) { + integral += this->interpolant().integrate(xLowLim, xUppLim, x1, x2, y1, y2); + } + + if (xUppLim == xHi) { + integrating = false; + } else { + xLowLim = xUppLim; + xUppLim = xHi; + lowIt = ranges::next(lowIt); + } + } + + if(inverted) { + integral *= -1.; + } + + return integral; +} \ No newline at end of file diff --git a/src/interpolation/table/Type/test/Type.test.cpp b/src/interpolation/table/Type/test/Type.test.cpp index 636e9bf..8fe61e6 100644 --- a/src/interpolation/table/Type/test/Type.test.cpp +++ b/src/interpolation/table/Type/test/Type.test.cpp @@ -46,6 +46,20 @@ SCENARIO("An interpolation table can be constructed"){ REQUIRE( 6.5 == myTable( 2.5, myTable.search() ) ); } + THEN("the table can be integrated"){ + REQUIRE(4. == myTable.integrate(1., 2.)); + REQUIRE(6.5 == myTable.integrate(2., 3.)); + REQUIRE(10.5 == myTable.integrate(1., 3.)); + // When doing integral over wider range, integration limits get + // truncated to the table range by default. + REQUIRE(10.5 == myTable.integrate(0., 4.)); + + REQUIRE(-4. == myTable.integrate(2., 1.)); + REQUIRE(-6.5 == myTable.integrate(3., 2.)); + REQUIRE(-10.5 == myTable.integrate(3., 1.)); + REQUIRE(-10.5 == myTable.integrate(4., 0.)); + } + THEN("the table handles discontinuities"){ REQUIRE( 5.0 == myTable( 2.0, myTable.search() ) ); } @@ -87,6 +101,13 @@ SCENARIO("An interpolation table can be constructed"){ REQUIRE( 6.5 == myTable( 2.5, myTable.search() ) ); } + THEN("the table can be integrated"){ + REQUIRE(4. == myTable.integrate(1., 2.)); + REQUIRE(6.5 == myTable.integrate(2., 3.)); + REQUIRE(10.5 == myTable.integrate(1., 3.)); + REQUIRE(10.5 == myTable.integrate(0., 4.)); + } + THEN("the table can provide x-values"){ REQUIRE( ranges::equal( myTable.x(), xGrid ) ); } @@ -111,6 +132,13 @@ SCENARIO("An interpolation table can be constructed"){ REQUIRE( 6.5 == myTable( 2.5, myTable.search() ) ); } + THEN("the table can be integrated"){ + REQUIRE(4. == myTable.integrate(1., 2.)); + REQUIRE(6.5 == myTable.integrate(2., 3.)); + REQUIRE(10.5 == myTable.integrate(1., 3.)); + REQUIRE(10.5 == myTable.integrate(0., 4.)); + } + THEN("the table handles discontinuities"){ REQUIRE( 5.0 == myTable( 2.0, myTable.search() ) ); } @@ -135,6 +163,18 @@ SCENARIO("An interpolation table can be constructed"){ REQUIRE( 6.5 == myTable( 2.5, search ) ); } + THEN("the table can be integrated"){ + REQUIRE(4. == myTable.integrate(1., 2.)); + REQUIRE(6.5 == myTable.integrate(2., 3.)); + REQUIRE(10.5 == myTable.integrate(1., 3.)); + REQUIRE(10.5 == myTable.integrate(0., 4.)); + + REQUIRE(-4. == myTable.integrate(2., 1.)); + REQUIRE(-6.5 == myTable.integrate(3., 2.)); + REQUIRE(-10.5 == myTable.integrate(3., 1.)); + REQUIRE(-10.5 == myTable.integrate(4., 0.)); + } + THEN("the table handles discontinuities"){ REQUIRE( 5.0 == myTable( 2.0, search ) ); REQUIRE( 5.0 == myTable( 2.0, search ) ); @@ -156,6 +196,18 @@ SCENARIO("An interpolation table can be constructed"){ REQUIRE( 6.5 * seconds == myTable( 2.5 * meter, myTable.search() ) ); } + THEN("the table can be integrated"){ + REQUIRE(4. * meter * seconds == myTable.integrate(1. * meter, 2. * meter)); + REQUIRE(6.5 * meter * seconds == myTable.integrate(2. * meter, 3. * meter)); + REQUIRE(10.5 * meter * seconds == myTable.integrate(1. * meter, 3. * meter)); + REQUIRE(10.5 * meter * seconds == myTable.integrate(0. * meter, 4. * meter)); + + REQUIRE(-4. * meter * seconds == myTable.integrate(2. * meter, 1. * meter)); + REQUIRE(-6.5 * meter * seconds == myTable.integrate(3. * meter, 2. * meter)); + REQUIRE(-10.5 * meter * seconds == myTable.integrate(3. * meter, 1. * meter)); + REQUIRE(-10.5 * meter * seconds == myTable.integrate(4. * meter, 0. * meter)); + } + THEN("the table can provide x-values"){ REQUIRE( ranges::equal( myTable.x(), x ) ); } diff --git a/src/interpolation/table/Variant.hpp b/src/interpolation/table/Variant.hpp index 67f0620..1f3b3d8 100644 --- a/src/interpolation/table/Variant.hpp +++ b/src/interpolation/table/Variant.hpp @@ -121,6 +121,17 @@ class Variant { }, this->core ); } + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xLow, Arg&& xHi, Args&&... args ) const { + return std::visit( [&]( auto&& core ){ return core.integrate( + std::forward(xLow), + std::forward(xHi), + std::forward(args)... ); + }, + this->core + ); + } + public: template< typename Arg > Variant( Arg&& arg ) : core( std::forward< Arg >( arg ) ){} diff --git a/src/interpolation/table/Variant/Table.hpp b/src/interpolation/table/Variant/Table.hpp index 932e104..c727a6f 100644 --- a/src/interpolation/table/Variant/Table.hpp +++ b/src/interpolation/table/Variant/Table.hpp @@ -32,6 +32,12 @@ class Table< table::Variant< Args... >, Decorators... > : return Parent:: template evaluate( std::forward(x), std::forward(args)... ); } + + template< typename Arg, typename... CallArgs > + auto integrate(Arg&& xLow, Arg&& xHi, CallArgs&&... args) const { + return Parent:: template integrate( std::forward(xLow), + std::forward(xHi), std::forward(args)... ); + } }; namespace table { diff --git a/src/interpolation/table/Variant/test/Variant.test.cpp b/src/interpolation/table/Variant/test/Variant.test.cpp index 376651e..95782ae 100644 --- a/src/interpolation/table/Variant/test/Variant.test.cpp +++ b/src/interpolation/table/Variant/test/Variant.test.cpp @@ -39,6 +39,8 @@ SCENARIO("An variant interpolation table can be constructed"){ REQUIRE( correct ); REQUIRE( 6.0 == myTable(2.5) ); + REQUIRE( 10. == myTable.integrate(1., 3.) ); + REQUIRE( -10. == myTable.integrate(3., 1.) ); REQUIRE( -infinity() == myTable.domainMin() ); REQUIRE( infinity() == myTable.domainMax() ); REQUIRE( xGrid.front() == myTable.tableMin() ); diff --git a/src/interpolation/table/Vector.hpp b/src/interpolation/table/Vector.hpp index 473ad63..2166918 100644 --- a/src/interpolation/table/Vector.hpp +++ b/src/interpolation/table/Vector.hpp @@ -59,6 +59,70 @@ struct Vector< TableType, Adjoining > { return (*iterator)( std::forward(x), std::forward(args)... ); } + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xL, Arg&& xH, Args&&... args ) const { + auto xLow = xL; + auto xHi = xH; + const bool reversed = xLow > xHi; + + if ( reversed ) { + auto xLowTmp = xLow; + xLow = xHi; + xHi = xLowTmp; + } + + // Integration may only be carried out over the function's valid domain + if (xLow < tableMin()) { + xLow = tableMin(); + } else if (xLow > tableMax()) { + xLow = tableMax(); + } + + if (xHi > tableMax()) { + xHi = tableMax(); + } else if (xHi < tableMin()) { + xHi = tableMin(); + } + + // Get region which contains xLow + auto region = ranges::lower_bound(this->core, xLow, ranges::less(), + []( auto&& table ){ return table.tableMax(); } ); + + if (region->tableMin() > xLow) { + region = ranges::prev(region); + } + + auto integral = (xLow - xLow) * region->y().front(); + double xLowLim = xLow; + double xUppLim = xHi; + bool integrating = true; + while (integrating) { + if (xLowLim < region->tableMin()) { + xLowLim = region->tableMin(); + } + + if (xUppLim > region->tableMax()) { + xUppLim = region->tableMax(); + } + + integral += region->integrate(xLowLim, xUppLim); + + if (xUppLim == xHi) { + integrating = false; + } else { + xLowLim = xUppLim; + xUppLim = xHi; + region = ranges::next(region); + } + } + + if (reversed) { + integral *= -1.; + } + + return integral; + } + auto cachedSearch() const { return this->core.front().cachedSearch(); } diff --git a/src/interpolation/table/Vector/Table.hpp b/src/interpolation/table/Vector/Table.hpp index d85472f..5a84371 100644 --- a/src/interpolation/table/Vector/Table.hpp +++ b/src/interpolation/table/Vector/Table.hpp @@ -36,6 +36,13 @@ class Table< table::Vector< Args... >, Decorators... > : return Parent:: template evaluate( std::forward(x), std::forward(args)... ); } + + template< typename Xdata, typename... CallArgs > + auto integrate(const Xdata& xLow, const Xdata& xHi, CallArgs&&... args) const { + return Parent:: template integrate( xLow, xHi, + std::forward(args)...); + } + }; namespace table { diff --git a/src/interpolation/table/Vector/test/Vector.test.cpp b/src/interpolation/table/Vector/test/Vector.test.cpp index 8798046..9a84cd0 100644 --- a/src/interpolation/table/Vector/test/Vector.test.cpp +++ b/src/interpolation/table/Vector/test/Vector.test.cpp @@ -64,6 +64,9 @@ SCENARIO("An variant interpolation table can be constructed"){ REQUIRE( core[0](2.5) == vc(2.5) ); REQUIRE( core[1](6.5) == vc(6.5) ); + + REQUIRE( core[0].integrate(1.,3.) + core[1].integrate(3.,8.) == vc.integrate(1., 8.) ); + REQUIRE( core[0].integrate(1.,3.) + core[1].integrate(3.,5.) == vc.integrate(1., 5.) ); } SECTION("exceptional behavior"){ diff --git a/src/interpolation/table/domain/Max.hpp b/src/interpolation/table/domain/Max.hpp index a544dc0..e70faf4 100644 --- a/src/interpolation/table/domain/Max.hpp +++ b/src/interpolation/table/domain/Max.hpp @@ -15,6 +15,29 @@ class Max : public Parent { std::forward(args)... ); } + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xLow, Arg&& xHi, Args&&... args ) const { + // Captures child leftLimit method + const auto& table = static_cast< const Table& >( *this ); + if ( xLow > table.domainMax() ){ + Log::error( "Lower integration limit is greater than right domain limit." ); + Log::info( "Lower integration limit: {}", xLow ); + Log::info( "Right domain limit: {}", table.domainMax() ); + throw std::out_of_range("Lower integration limit is outside table domain."); + } + + if ( xHi > table.domainMax() ){ + Log::error( "Upper integration limit is greater than right domain limit." ); + Log::info( "Uppper integration limit: {}", xHi ); + Log::info( "Right domain limit: {}", table.domainMax() ); + throw std::out_of_range("Upper integration limit is outside table domain."); + } + + return Parent:: template integrate( std::forward(xLow), + std::forward(xHi), + std::forward(args)... ); + } + public: static constexpr auto specifiesDomainMax(){ return true; } diff --git a/src/interpolation/table/domain/Min.hpp b/src/interpolation/table/domain/Min.hpp index 9182476..8aaa580 100644 --- a/src/interpolation/table/domain/Min.hpp +++ b/src/interpolation/table/domain/Min.hpp @@ -15,6 +15,29 @@ class Min : public Parent { std::forward(args)... ); } + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xLow, Arg&& xHi, Args&&... args ) const { + // Captures child leftLimit method + const auto& table = static_cast< const Table& >( *this ); + if ( xLow < table.domainMin() ){ + Log::error( "Lower integration limit is less than left domain limit." ); + Log::info( "Lower integration limit: {}", xLow ); + Log::info( "Left domain limit: {}", table.domainMin() ); + throw std::out_of_range("Lower integration limit is outside table domain."); + } + + if ( xHi < table.domainMin() ){ + Log::error( "Upper integration limit is less than left domain limit." ); + Log::info( "Uppper integration limit: {}", xHi ); + Log::info( "Left domain limit: {}", table.domainMin() ); + throw std::out_of_range("Upper integration limit is outside table domain."); + } + + return Parent:: template integrate( std::forward(xLow), + std::forward(xHi), + std::forward(args)... ); + } + public: static constexpr auto specifiesDomainMin(){ return true; } diff --git a/src/interpolation/table/domain/max/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp b/src/interpolation/table/domain/max/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp index 1829662..1e57b2d 100644 --- a/src/interpolation/table/domain/max/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp +++ b/src/interpolation/table/domain/max/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp @@ -36,5 +36,13 @@ SCENARIO("An interpolation table can be constructed with a domain maximum"){ THEN("the table will throw for values grater than the domain max"){ REQUIRE_THROWS( myTable( 5.0 ) ); } + + THEN("the table can be integrated in the center interval") { + REQUIRE( 10. == myTable.integrate(1., 3.) ); + } + + THEN("the table will throw for integration limits greater than the domain max") { + REQUIRE_THROWS( myTable.integrate(1., 5.) ); + } } } diff --git a/src/interpolation/table/domain/max/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp b/src/interpolation/table/domain/max/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp index c378046..aa6f353 100644 --- a/src/interpolation/table/domain/max/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp +++ b/src/interpolation/table/domain/max/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp @@ -33,5 +33,13 @@ SCENARIO("An interpolation table can be constructed with a domain maximum"){ THEN("the table will throw for values greater than the domain max"){ REQUIRE_THROWS( myTable( 5.0 ) ); } + + THEN("the table can be integrated in the center interval") { + REQUIRE( 10. == myTable.integrate(1., 3.) ); + } + + THEN("the table will throw for integration limits greater than the domain max") { + REQUIRE_THROWS( myTable.integrate(1., 5.) ); + } } } diff --git a/src/interpolation/table/domain/min/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp b/src/interpolation/table/domain/min/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp index d5880b2..1bb839f 100644 --- a/src/interpolation/table/domain/min/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp +++ b/src/interpolation/table/domain/min/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp @@ -35,5 +35,15 @@ SCENARIO("An interpolation table can be constructed" REQUIRE_THROWS( myTable( -1.0 ) ); REQUIRE_THROWS( myTable( -1.0, myTable.search() ) ); } + + THEN("the table can be integrated in the center interval"){ + REQUIRE( 10. == myTable.integrate(1., 3.) ); + REQUIRE( 10. == myTable.integrate(1., 3., myTable.search()) ); + } + + THEN("the table will throw for integration limits less than the domain min"){ + REQUIRE_THROWS( myTable.integrate(-1., 3.) ); + REQUIRE_THROWS( myTable.integrate(-1., 3., myTable.search()) ); + } } } diff --git a/src/interpolation/table/domain/min/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp b/src/interpolation/table/domain/min/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp index 28a1e3b..3710ed3 100644 --- a/src/interpolation/table/domain/min/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp +++ b/src/interpolation/table/domain/min/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp @@ -32,5 +32,15 @@ SCENARIO("An interpolation table can be constructed with a domain minimum"){ REQUIRE_THROWS( myTable( -1.0 ) ); REQUIRE_THROWS( myTable( -1.0, myTable.search() ) ); } + + THEN("the table can be integrated in the center interval"){ + REQUIRE( 10. == myTable.integrate(1., 3.) ); + REQUIRE( 10. == myTable.integrate(1., 3., myTable.search()) ); + } + + THEN("the table will throw for integration limits less than the domain min"){ + REQUIRE_THROWS( myTable.integrate(-1., 3.) ); + REQUIRE_THROWS( myTable.integrate(-1., 3., myTable.search()) ); + } } } diff --git a/src/interpolation/table/left/interval.hpp b/src/interpolation/table/left/interval.hpp index 7b9e45a..d463065 100644 --- a/src/interpolation/table/left/interval.hpp +++ b/src/interpolation/table/left/interval.hpp @@ -14,6 +14,35 @@ class Interval : public Parent { Parent::template evaluate( std::forward(x), std::forward(args)... ); } + + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xL, Arg&& xH, Args&&... args ) const { + // Captures child leftLimit method + const auto& table = static_cast< const Table& >( *this ); + + const auto xLow = xL; + const auto xHi = xH; + + // Integral automatically cuts off to range to table + auto integral = Parent:: template integrate( std::forward(xL), + std::forward(xH), + std::forward(args)... ); + + // Must how add portions which occur outside of table + const bool reversed = xLow > xHi; + if(xLow < table.tableMin() && xHi < table.tableMin()) { + integral += (xHi - xLow) * table.leftIntervalValue(); + } else if(reversed) { + // If didn't select first case, and we are reversed, it must be that xLow is + // within the interval, but xHi isn't, so we only need to count tableMin + // to xHi, in negative sense. + integral += (xHi - table.tableMin()) * table.leftIntervalValue(); + } else if(xLow < table.tableMin()){ + integral += (table.tableMin() - xLow) * table.leftIntervalValue(); + } + + return integral; + } public: static constexpr auto specifiesLeftInterval(){ return true; } diff --git a/src/interpolation/table/left/interval/IsAsymptotic/test/IsAsymptotic.test.cpp b/src/interpolation/table/left/interval/IsAsymptotic/test/IsAsymptotic.test.cpp index 99cfb5f..66c5f08 100644 --- a/src/interpolation/table/left/interval/IsAsymptotic/test/IsAsymptotic.test.cpp +++ b/src/interpolation/table/left/interval/IsAsymptotic/test/IsAsymptotic.test.cpp @@ -25,6 +25,11 @@ SCENARIO("An asymptotic right interval can be applied to a table"){ REQUIRE( myTable(3.5) == 8.0 ); REQUIRE( myTable(0.0) == 3.0 ); + + // Integration tests for left + REQUIRE( 4. == myTable.integrate(1., 2.) ); + REQUIRE( 7. == myTable.integrate(0., 2.) ); + REQUIRE( -7. == myTable.integrate(2., 0.) ); } WHEN("constructing a Table without the search method"){ @@ -63,6 +68,11 @@ SCENARIO("An asymptotic right interval can be applied to a table"){ REQUIRE( myTable(3.5) == 8.0 ); REQUIRE( myTable(0.0) == 3.0 ); + + // Integration tests for left + REQUIRE( 4. == myTable.integrate(1., 2.) ); + REQUIRE( 7. == myTable.integrate(0., 2.) ); + REQUIRE( -7. == myTable.integrate(2., 0.) ); } } } diff --git a/src/interpolation/table/left/interval/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp b/src/interpolation/table/left/interval/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp index 48cede9..d562b8e 100644 --- a/src/interpolation/table/left/interval/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp +++ b/src/interpolation/table/left/interval/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp @@ -33,5 +33,13 @@ SCENARIO("An interpolation table can be constructed" THEN("the table will throw in the left interval"){ REQUIRE( 10.0 == myTable( 0.0, myTable.search() ) ); } + + THEN("the table can be integrated in the center interval") { + REQUIRE( 10. == myTable.integrate(1., 3.) ); + } + + THEN("the table can be integrated left of the interval") { + REQUIRE( 20. == myTable.integrate(0., 3.) ); + } } } diff --git a/src/interpolation/table/left/interval/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp b/src/interpolation/table/left/interval/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp index 26e3764..1125417 100644 --- a/src/interpolation/table/left/interval/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp +++ b/src/interpolation/table/left/interval/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp @@ -34,5 +34,15 @@ SCENARIO("An interpolation table can be constructed" REQUIRE( 10.0 == myTable( 0.0 ) ); REQUIRE( 10.0 == myTable( 0.0, myTable.search() ) ); } + + THEN("the table can be integrated in the center interval") { + REQUIRE( 10. == myTable.integrate(1., 3.) ); + REQUIRE( 10. == myTable.integrate(1., 3., myTable.search()) ); + } + + THEN("the table can be integrated left of the interval") { + REQUIRE( 20. == myTable.integrate(0., 3.) ); + REQUIRE( 20. == myTable.integrate(0., 3., myTable.search()) ); + } } } diff --git a/src/interpolation/table/left/interval/Throws.hpp b/src/interpolation/table/left/interval/Throws.hpp index 4496fc9..c92ca5a 100644 --- a/src/interpolation/table/left/interval/Throws.hpp +++ b/src/interpolation/table/left/interval/Throws.hpp @@ -16,6 +16,33 @@ struct Throws { return Parent::template evaluate( std::forward(x), std::forward(args)... ); } + + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xLow, Arg&& xHi, Args&&... args ) const { + // Captures child leftLimit method + const auto& table = static_cast< const Table& >( *this ); + if ( xLow < table.tableMin() ){ + Log::error("Lower integration limit in interval with undefined behavior."); + Log::info("Lower integration limit: {}", xLow); + Log::info("The function is undefined for values less than: {}", + table.tableMin() ); + throw std::out_of_range( + "Lower integration limit in interval with undefined behavior." ); + } + + if ( xHi < table.tableMin() ){ + Log::error("Upper integration limit in interval with undefined behavior."); + Log::info("Upper integration limit: {}", xHi); + Log::info("The function is undefined for values less than: {}", + table.tableMin() ); + throw std::out_of_range( + "Upper integration limit in interval with undefined behavior." ); + } + + return Parent:: template integrate( std::forward(xLow), + std::forward(xHi), + std::forward(args)... ); + } public: template< typename... Args > diff --git a/src/interpolation/table/left/interval/Throws/test/Throws.test.cpp b/src/interpolation/table/left/interval/Throws/test/Throws.test.cpp index c0121d1..d034eaf 100644 --- a/src/interpolation/table/left/interval/Throws/test/Throws.test.cpp +++ b/src/interpolation/table/left/interval/Throws/test/Throws.test.cpp @@ -32,5 +32,13 @@ SCENARIO("An interpolation table can be constructed" THEN("the table will throw in the left interval"){ REQUIRE_THROWS( myTable( 0.0, myTable.search() ) ); } + + THEN("the table can be integrated in the center interval") { + REQUIRE( 10. == myTable.integrate(1., 3.) ); + } + + THEN("the table will throw when attempting to integrate left of the interval") { + REQUIRE_THROWS( myTable.integrate(0., 3.) ); + } } } diff --git a/src/interpolation/table/right/interval.hpp b/src/interpolation/table/right/interval.hpp index 8ec7170..17c12c2 100644 --- a/src/interpolation/table/right/interval.hpp +++ b/src/interpolation/table/right/interval.hpp @@ -14,6 +14,35 @@ class Interval : public Parent { Parent::template evaluate( std::forward(x), std::forward(args)... ); } + + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xL, Arg&& xH, Args&&... args ) const { + // Captures child leftLimit method + const auto& table = static_cast< const Table& >( *this ); + + const auto xLow = xL; + const auto xHi = xH; + + // Integral automatically cuts off to range to table + auto integral = Parent:: template integrate( std::forward(xL), + std::forward(xH), + std::forward(args)... ); + + // Must how add portions which occur outside of table + const bool reversed = xLow > xHi; + if(xLow > table.tableMax() && xHi > table.tableMax()) { + integral += (xHi - xLow) * table.rightIntervalValue(); + } else if(reversed) { + // If didn't select first case, and we are reversed, it must be that xHi is + // within the interval, but xLow isn't, so we only need to count tableMax + // to xLow, in negative sense. + integral += (table.tableMax() - xLow) * table.rightIntervalValue(); + } else if(xHi > table.tableMax()){ + integral += (xHi - table.tableMax()) * table.rightIntervalValue(); + } + + return integral; + } public: static constexpr auto specifiesRightInterval(){ return true; } diff --git a/src/interpolation/table/right/interval/IsAsymptotic/test/IsAsymptotic.test.cpp b/src/interpolation/table/right/interval/IsAsymptotic/test/IsAsymptotic.test.cpp index 4886f0e..5d5609f 100644 --- a/src/interpolation/table/right/interval/IsAsymptotic/test/IsAsymptotic.test.cpp +++ b/src/interpolation/table/right/interval/IsAsymptotic/test/IsAsymptotic.test.cpp @@ -25,6 +25,9 @@ SCENARIO("An asymptotic right interval can be applied to a table"){ REQUIRE( myTable(3.5) == 8.0 ); REQUIRE( myTable(9.0) == 17.0 ); + + REQUIRE( myTable.integrate(1., 8.) == 70.); + REQUIRE( myTable.integrate(1., 9.) == 87.); } WHEN("constructing a Table without the search method"){ @@ -63,6 +66,9 @@ SCENARIO("An asymptotic right interval can be applied to a table"){ REQUIRE( myTable(3.5) == 8.0 ); REQUIRE( myTable(9.0) == 17.0 ); + + REQUIRE( myTable.integrate(1., 8.) == 70.); + REQUIRE( myTable.integrate(1., 9.) == 87.); } } } diff --git a/src/interpolation/table/right/interval/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp b/src/interpolation/table/right/interval/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp index fd4fabe..56f2a80 100644 --- a/src/interpolation/table/right/interval/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp +++ b/src/interpolation/table/right/interval/IsCompiletimeConstant/test/IsCompiletimeConstant.test.cpp @@ -33,5 +33,13 @@ SCENARIO("An interpolation table can be constructed" THEN("the table will throw in the right interval"){ REQUIRE( 10.0 == myTable( 4.0, myTable.search() ) ); } + + THEN("the table can be integrated in the center interval") { + REQUIRE( 10. == myTable.integrate(1., 3.) ); + } + + THEN("the table can be integrated left of the interval") { + REQUIRE( 20. == myTable.integrate(1., 4.) ); + } } } diff --git a/src/interpolation/table/right/interval/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp b/src/interpolation/table/right/interval/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp index 2203d87..0e13770 100644 --- a/src/interpolation/table/right/interval/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp +++ b/src/interpolation/table/right/interval/IsRuntimeConstant/test/IsRuntimeConstant.test.cpp @@ -31,5 +31,13 @@ SCENARIO("An interpolation table can be constructed" THEN("the table will throw in the right interval"){ REQUIRE( 10.0 == myTable( 4.0, myTable.search() ) ); } + + THEN("the table can be integrated in the center interval") { + REQUIRE( 10. == myTable.integrate(1., 3., myTable.search()) ); + } + + THEN("the table can be integrated left of the interval") { + REQUIRE( 20. == myTable.integrate(1., 4., myTable.search()) ); + } } } diff --git a/src/interpolation/table/right/interval/Throws.hpp b/src/interpolation/table/right/interval/Throws.hpp index d403f3e..d1a12ea 100644 --- a/src/interpolation/table/right/interval/Throws.hpp +++ b/src/interpolation/table/right/interval/Throws.hpp @@ -16,6 +16,33 @@ struct Throws { std::forward(args)... ); } + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xLow, Arg&& xHi, Args&&... args ) const { + // Captures child leftLimit method + const auto& table = static_cast< const Table& >( *this ); + if ( xLow > table.tableMax() ){ + Log::error("Lower integration limit in interval with undefined behavior."); + Log::info( "Lower integration limit: {}", xLow ); + Log::info("The function is undefined for values greater than: {}", + table.tableMax() ); + throw std::out_of_range( + "Lower integration limit in interval with undefined behavior." ); + } + + if ( xHi > table.tableMax() ){ + Log::error("Upper integration limit in interval with undefined behavior."); + Log::info( "Upper integration limit: {}", xHi ); + Log::info("The function is undefined for values greater than: {}", + table.tableMax() ); + throw std::out_of_range( + "Upper integration limit in interval with undefined behavior." ); + } + + return Parent:: template integrate( std::forward(xLow), + std::forward(xHi), + std::forward(args)... ); + } + public: template< typename... Args > Type( Args&&... args ) : Interval( std::forward< Args >( args )... ){} diff --git a/src/interpolation/table/right/interval/Throws/test/Throws.test.cpp b/src/interpolation/table/right/interval/Throws/test/Throws.test.cpp index 78f1656..4363779 100644 --- a/src/interpolation/table/right/interval/Throws/test/Throws.test.cpp +++ b/src/interpolation/table/right/interval/Throws/test/Throws.test.cpp @@ -32,5 +32,13 @@ SCENARIO("An interpolation table can be constructed" THEN("the table will throw in the right interval"){ REQUIRE_THROWS( myTable( 4.0, myTable.search() ) ); } + + THEN("the table can be integrated in the center interval") { + REQUIRE( 10. == myTable.integrate(1., 3., myTable.search()) ); + } + + THEN("the table will throw when attempting to integrate right of the interval") { + REQUIRE_THROWS( myTable.integrate(1., 4., myTable.search()) ); + } } } diff --git a/src/interpolation/table/transform.hpp b/src/interpolation/table/transform.hpp index 3e72391..6e19538 100644 --- a/src/interpolation/table/transform.hpp +++ b/src/interpolation/table/transform.hpp @@ -12,6 +12,19 @@ struct X { ( Transform::invert( std::forward(x) ), std::forward(args)... ); } + + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xLow, Arg&& xHi, Args&&... args ) const { + // Throw and error immediately + Log::error("Integration of a table with X Transform is undefined."); + throw std::logic_error( + "Integration of a table with X Transform is undefined." ); + + // This will NEVER get called in theory, but is there just for continuity + return Parent:: template integrate( std::forward(xLow), + std::forward(xHi), + std::forward(args)... ); + } public: using Xdata = @@ -46,6 +59,19 @@ struct Y { ( Parent::template evaluate( std::forward(x), std::forward(args)... ) ); } + + template< typename Table, typename Arg, typename... Args > + auto integrate( Arg&& xLow, Arg&& xHi, Args&&... args ) const { + // Throw and error immediately + Log::error("Integration of a table with Y Transform is undefined."); + throw std::logic_error( + "Integration of a table with Y Transform is undefined." ); + + // This will NEVER get called in theory, but is there just for continuity + return Parent:: template integrate( std::forward(xLow), + std::forward(xHi), + std::forward(args)... ); + } public: using Ydata = diff --git a/src/interpolation/table/transform/test/transform.test.cpp b/src/interpolation/table/transform/test/transform.test.cpp index 3b3db6e..a37b9a8 100644 --- a/src/interpolation/table/transform/test/transform.test.cpp +++ b/src/interpolation/table/transform/test/transform.test.cpp @@ -77,6 +77,10 @@ SCENARIO("An interpolation table can be constructed with transforms"){ REQUIRE( std::exp(4.0) == trial( std::exp(1.5), trial.search() ) ) ; } + THEN("tables with a transform will throw upon attempted integration") { + REQUIRE_THROWS(trial.integrate(1., 3.)); + } + THEN("lifetime of iterator is not limited to the "){ auto persistance = []( auto&& arg ){ return arg.begin(); }; auto maybe = persistance( trial.x() ); diff --git a/src/interpolation/test/Example1.cpp b/src/interpolation/test/Example1.cpp index 11fb27c..a5bfedf 100644 --- a/src/interpolation/test/Example1.cpp +++ b/src/interpolation/test/Example1.cpp @@ -26,6 +26,7 @@ TEST_CASE("Example1"){ // Interpolation tables are function objects and are callable REQUIRE( myTable( 2.5 ) == 6.0 ); + REQUIRE( myTable.integrate(1., 2.5) == 6.75 ); REQUIRE( myTable.domainMin() == -infinity() ); REQUIRE( myTable.domainMax() == infinity() ); REQUIRE( myTable.tableMin() == 1.0 ); diff --git a/src/interpolation/test/Example5.cpp b/src/interpolation/test/Example5.cpp index a51850f..cecd527 100644 --- a/src/interpolation/test/Example5.cpp +++ b/src/interpolation/test/Example5.cpp @@ -23,4 +23,5 @@ TEST_CASE("Example5"){ REQUIRE( myTable( 2.5 * electronVolts ) == 6.0 * barns ); REQUIRE( myTable( 2.5E-3 * kilo(electronVolts) ) == 6.0 * barns ); + REQUIRE( myTable.integrate(1. * electronVolts, 3. * electronVolts) == 10. * electronVolts * barns ); } diff --git a/src/interpolation/test/Example8.cpp b/src/interpolation/test/Example8.cpp index d14a2cf..f67d77d 100644 --- a/src/interpolation/test/Example8.cpp +++ b/src/interpolation/test/Example8.cpp @@ -52,4 +52,5 @@ TEST_CASE("Example 8"){ Tab1 myTable( std::move(core) ); REQUIRE( myTable( 2.5 ) == 6.0 ); + REQUIRE( myTable.integrate(1., 3.) == 10. ); }