Skip to content
Snippets Groups Projects
Select Git revision
  • develop
  • master default protected
  • feature/triangulation-qhull
  • jst
  • ti_lab_build
  • features/splines_and_piecewise_polynomials
  • ma_2018/erraji
  • fabian
  • ITABase_v2024a
  • VA_v2023b
  • VA_v2023a
  • VA_v2022a
  • before_cmake_rework
  • v2021.a
  • v2020.a
  • v2019.a
  • v2018.b
  • v2018.a
  • v2017.d
  • v2017.c
  • v2017.b
  • v2017.a
  • v2016.a
23 results

ITABaseUnitTestPiesewisePolynomials.cpp

Blame
  • 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 );
    		}
    	}
    }