Select Git revision
ITABaseUnitTestPiesewisePolynomials.cpp
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
ITABaseUnitTestPiesewisePolynomials.cpp 9.21 KiB
/*
* ----------------------------------------------------------------
*
* ITA Base
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*
*
*/
#include "ITABaseTestUtils.h"
#include <ITABase/Math/PiecewisePolynomial.h>
#include <random>
#include <limits>
using namespace ITABase::Math;
const double X0 = -3.0;
const double X1 = -1.0;
const double X2 = 1.0;
const double X3 = 3.0;
//----------------------------------------//
// ---- POLYNOM COMPARISON FUNCTIONS ---- //
//----------------------------------------//
/// Piecewise linear function (x+2 | 2 | -x+2)
/// not continous
void LinearLinearLinear( const double& x, double& y, double& dY )
{
//NOTE: Extrapolation is implicitely included for linear functions
if( x < X1 )
{
y = x + 2;
dY = 1;
}
else if( x < X2 )
{
y = 2;
dY = 0;
}
else
{
y = -x + 2;
dY = -1;
}
}
void LinearLinearLinear( const std::vector<double>& vdX, std::vector<double>& vdY, std::vector<double>& vdDerivedY )
{
vdY.resize( vdX.size( ) );
vdDerivedY.resize( vdX.size( ) );
for( int idx = 0; idx < vdX.size( ); idx++ )
LinearLinearLinear( vdX[idx], vdY[idx], vdDerivedY[idx] );
}
/// Piecewise linear|quadratic|linear function (-2x+3 | x^2 | 2x-1)
/// C1 continous
void LinearQuadraticLinear( const double& x, double& y, double& dY )
{
// NOTE: Extrapolation is implicitely included for linear functions
if( x < X1 )
{
y = -2 * x + 3;
dY = -2;
}
else if( x < X2 )
{
y = x * x;
dY = 2 * x;
}
else
{
y = 2 * x - 1;
dY = 2;
}
}
void LinearQuadraticLinear( const std::vector<double>& vdX, std::vector<double>& vdY, std::vector<double>& vdDerivedY )
{
vdY.resize( vdX.size( ) );
vdDerivedY.resize( vdX.size( ) );
for( int idx = 0; idx < vdX.size( ); idx++ )
LinearQuadraticLinear( vdX[idx], vdY[idx], vdDerivedY[idx] );
}
/// Piecewise linear|cubic|quadratic function (x | x^3 | x^2)
/// C0 continous
void LinearCubicQuadratic( const double& x, double& y, double& dY )
{
if( x < X1 ) // Left-sided extrapolation implicitely included
{
y = x;
dY = 1;
}
else if( x < X2 )
{
y = x * x * x;
dY = 3 * x * x;
}
else if( x <= X3 )
{
y = x * x;
dY = 2 * x;
}
else // Right-sided extrapolation
{
y = X3 * X3 + ( 2 * X3 ) * ( x - X3 ); // f(X3) + f'(X3) * Delta-X
dY = 2 * X3 + 2 * ( x - X3 );
}
}
void LinearCubicQuadratic( const std::vector<double>& vdX, std::vector<double>& vdY, std::vector<double>& vdDerivedY )
{
vdY.resize( vdX.size( ) );
vdDerivedY.resize( vdX.size( ) );
for( int idx = 0; idx < vdX.size( ); idx++ )
LinearCubicQuadratic( vdX[idx], vdY[idx], vdDerivedY[idx] );
}
std::vector<double> GetRandomVector( const int iNumValues, const double dMin, const double dMax )
{
return GENERATE_COPY( chunk( iNumValues, take( iNumValues, random( dMin, dMax ) ) ) );
}
TEST_CASE( "ITABase::PiecewisePolynomial/General", "[ITABase][PiecewisePolynomial][Recommended]" )
{
const int iOrder = 3;
const std::vector<double> vdBreakPoints = { 0, 1, 2 };
const int iNumCoeffsPerPoly = iOrder + 1;
const int iNumIntervals = vdBreakPoints.size( ) - 1;
const int iNumCoeffsTotal = iNumIntervals * iNumCoeffsPerPoly;
const std::vector<double> vdCoeffs = GetRandomVector( iNumCoeffsTotal, -1000.0, 1000.0 );
const auto myPoly = CPiecewisePolynomial( vdBreakPoints, vdCoeffs, iOrder );
SECTION( "Exceptions" )
{
SECTION( "Constructor" )
{
REQUIRE_THROWS( CPiecewisePolynomial( vdBreakPoints, vdCoeffs, -1 ) );
REQUIRE_THROWS( CPiecewisePolynomial( std::vector<double>( 1 ), vdCoeffs, iOrder ) );
REQUIRE_THROWS( CPiecewisePolynomial( vdBreakPoints, std::vector<double>( vdCoeffs.size( ) + 1 ), iOrder ) );
}
SECTION( "Derivation( )" )
{
REQUIRE_NOTHROW( myPoly.Derivation( ) );
}
SECTION( "Evaluate( )" )
{
REQUIRE_NOTHROW( myPoly.Evaluate( 0.5 ) );
REQUIRE_NOTHROW( myPoly.Evaluate( vdBreakPoints.front( ) ) );
REQUIRE_NOTHROW( myPoly.Evaluate( vdBreakPoints.back( ) ) );
REQUIRE_NOTHROW( myPoly.Evaluate( vdBreakPoints.front( ) - 1.0, true ) );
REQUIRE_NOTHROW( myPoly.Evaluate( vdBreakPoints.back( ) + 1.0, true ) );
REQUIRE_THROWS( myPoly.Evaluate( vdBreakPoints.front( )-1.0 ) );
REQUIRE_THROWS( myPoly.Evaluate( vdBreakPoints.back( ) + 1.0 ) );
}
}
SECTION( "Get Methods" )
{
REQUIRE( myPoly.Order( ) == iOrder );
REQUIRE( myPoly.BreakPoints( ).size( ) == vdBreakPoints.size( ) );
REQUIRE( myPoly.Coefficients( ).size( ) == iNumCoeffsTotal );
REQUIRE( myPoly.NumCoefficients( ) == iNumCoeffsPerPoly );
REQUIRE( myPoly.NumIntervals( ) == iNumIntervals );
}
}
/// Creates a piecewise polynomial and its derivation and compares this to a set of expected values
void ComparePPEvaluationWithExpectedData( const std::vector<double>& vdBreakPoints, const std::vector<double>& vdCoeffs, const int iOrder, const std::vector<double>& vdEvalPoints,
const std::vector<double>& vdExpectedY, const std::vector<double>& vdExpectedDerivedY, bool bExtrapolate = false )
{
CPiecewisePolynomial myPoly = CPiecewisePolynomial( vdBreakPoints, vdCoeffs, iOrder );
CPiecewisePolynomial myDerivedPoly = myPoly.Derivation( );
const std::vector<double> vdActualY = myPoly.Evaluate( vdEvalPoints, bExtrapolate );
const std::vector<double> vdActualDerivedY = myDerivedPoly.Evaluate( vdEvalPoints, bExtrapolate );
constexpr double D_EPSILON = std::numeric_limits<double>::epsilon( ) * 10.0;
for( int idx = 0; idx < vdExpectedY.size( ); idx++ )
{
REQUIRE_THAT( vdActualY[idx], Catch::Matchers::WithinAbs( vdExpectedY[idx], D_EPSILON ) );
REQUIRE_THAT( vdActualDerivedY[idx], Catch::Matchers::WithinAbs( vdExpectedDerivedY[idx], D_EPSILON ) );
}
}
TEST_CASE( "ITABase::PiecewisePolynomial/Polynom-Evaluation", "[ITABase][PiecewisePolynomial][Required]" )
{
const int iNumRandEvalPoints = 50;
const std::vector<double> vdBreakPoints = { X0, X1, X2, X3 };
const int iNumIntervals = vdBreakPoints.size( ) - 1;
std::vector<double> vdCoeffs;
std::vector<double> vdRandEvalPoints = GetRandomVector( iNumRandEvalPoints, X0, X3);
std::vector<double> vdExpectedY;
std::vector<double> vdExpectedDerivedY;
int iOrder = 1;
SECTION( "LinearLinearLinear" )
{
iOrder = GENERATE( 1, 2, 3 );
if( iOrder == 3 )
vdCoeffs = { 0, 0, 1, -1, 0, 0, 0, 2, 0, 0, -1, 1 };
else if( iOrder == 2 )
vdCoeffs = { 0, 1, -1, 0, 0, 2, 0, -1, 1 };
else if( iOrder == 1 )
vdCoeffs = { 1, -1, 0, 2, -1, 1 };
SECTION( "Random evaluation points" )
{
LinearLinearLinear( vdRandEvalPoints, vdExpectedY, vdExpectedDerivedY );
ComparePPEvaluationWithExpectedData( vdBreakPoints, vdCoeffs, iOrder, vdRandEvalPoints, vdExpectedY, vdExpectedDerivedY );
}
SECTION( "Evaluation at interval limits" )
{
LinearLinearLinear( vdBreakPoints, vdExpectedY, vdExpectedDerivedY );
ComparePPEvaluationWithExpectedData( vdBreakPoints, vdCoeffs, iOrder, vdBreakPoints, vdExpectedY, vdExpectedDerivedY );
}
}
SECTION( "LinearQuadraticLinear" )
{
iOrder = GENERATE( 2, 3 );
if( iOrder == 3 )
vdCoeffs = { 0, 0, -2, 9, 0, 1, -2, 1, 0, 0, 2, 1 };
else if( iOrder == 2 )
vdCoeffs = { 0, -2, 9, 1, -2, 1, 0, 2, 1 };
SECTION( "Random evaluation points" )
{
LinearQuadraticLinear( vdRandEvalPoints, vdExpectedY, vdExpectedDerivedY );
ComparePPEvaluationWithExpectedData( vdBreakPoints, vdCoeffs, iOrder, vdRandEvalPoints, vdExpectedY, vdExpectedDerivedY );
}
SECTION( "Evaluation at interval limits" )
{
LinearQuadraticLinear( vdBreakPoints, vdExpectedY, vdExpectedDerivedY );
ComparePPEvaluationWithExpectedData( vdBreakPoints, vdCoeffs, iOrder, vdBreakPoints, vdExpectedY, vdExpectedDerivedY );
}
}
SECTION( "LinearCubicQuadratic" )
{
iOrder = 3;
vdCoeffs = { 0, 0, 1, -3, 1, -3, 3, -1, 0, 1, 2, 1 };
SECTION( "Random evaluation points" )
{
LinearCubicQuadratic( vdRandEvalPoints, vdExpectedY, vdExpectedDerivedY );
ComparePPEvaluationWithExpectedData( vdBreakPoints, vdCoeffs, iOrder, vdRandEvalPoints, vdExpectedY, vdExpectedDerivedY );
}
SECTION( "Evaluation at interval limits" )
{
LinearCubicQuadratic( vdBreakPoints, vdExpectedY, vdExpectedDerivedY );
ComparePPEvaluationWithExpectedData( vdBreakPoints, vdCoeffs, iOrder, vdBreakPoints, vdExpectedY, vdExpectedDerivedY );
}
SECTION( "Left-sided extrapolation" )
{
std::vector<double> vdExtrapolPoints = GetRandomVector( iNumRandEvalPoints, X0 - 100.0, X0 - 0.1 );
LinearCubicQuadratic( vdExtrapolPoints, vdExpectedY, vdExpectedDerivedY );
ComparePPEvaluationWithExpectedData( vdBreakPoints, vdCoeffs, iOrder, vdExtrapolPoints, vdExpectedY, vdExpectedDerivedY, true );
}
SECTION( "Right-sided extrapolation" )
{
std::vector<double> vdExtrapolPoints = GetRandomVector( iNumRandEvalPoints, X3 + 0.1, X3 + 100.0 );
LinearCubicQuadratic( vdExtrapolPoints, vdExpectedY, vdExpectedDerivedY );
ComparePPEvaluationWithExpectedData( vdBreakPoints, vdCoeffs, iOrder, vdExtrapolPoints, vdExpectedY, vdExpectedDerivedY, true );
}
}
}