Commit 05c3fe81 authored by Philipp Schäfer's avatar Philipp Schäfer
Browse files

PiecewisePolynomials and Spline

- cleaned up and optimized alot of code
parent ba70cb25
......@@ -43,8 +43,10 @@ namespace ITABase
class ITA_BASE_API CPiecewisePolynomial
{
public:
//! Main constructor using double vectors of break points and coefficients for the intervals between breakpoints to define a piecewise polynomial
CPiecewisePolynomial(const std::vector<double>& vdBreakPoints, const std::vector<double>& vdCoeffs, const int iOrder);
CPiecewisePolynomial(const std::vector<float>& vdBreakPnts, const std::vector<float>& vdCoeffs, const int iOrder);
//! Constructor using float vectors, which are internally converted to double vectors
CPiecewisePolynomial(const std::vector<float>& vdBreakPoints, const std::vector<float>& vdCoeffs, const int iOrder);
//! Returns the polynomial order
inline int Order() const { return iOrder; };
//! Returns the number of intervals (pieces) between break points
......@@ -58,13 +60,11 @@ namespace ITABase
//! Evaluate the piecewise polynomial at a given x value. If desired, xValues from outsite the breakpoints can be extrapolated linearly.
double Value(const double& xValue, bool bExtrapolate = false) const;
//! Evaluate the piecewise polynomial at a given x value. If desired, xValues from outsite the breakpoints can be extrapolated linearly.
double Value(const float& xValue, bool bExtrapolate = false) const;
//! Evaluates the piecewise polynomial at all given x values. If desired, xValues from outsite the breakpoints can be extrapolated linearly.
std::vector<double> Value(const std::vector<double>& vdXValues, bool bExtrapolate = false) const;
//! Evaluates the piecewise polynomial at all given x values. If desired, xValues from outsite the breakpoints can be extrapolated linearly.
std::vector<double> Value(const std::vector<float>& vdXValues, bool bExtrapolate = false) const;
std::vector<double> Value(const std::vector<float>& vfXValues, bool bExtrapolate = false) const;
//! Finds the index of Interval at a given x value. If extrapolation is desired, it will return -1 for xValue outsite the breakpoints
......
......@@ -27,26 +27,27 @@ namespace ITABase
{
namespace Math
{
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs vdSupportingPoints and vdDataPoints. vdSupportingPoints need to be sorted in ascending order
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints);
//All other Permutations of data Types
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given supporting and data points. The supporting points need to be sorted in ascending order.
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints);
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs vdSupportingPoints and vdDataPoints. vdSupportingPoints need to be sorted in ascending order
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints);
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs vdSupportingPoints and vdDataPoints. vdSupportingPoints need to be sorted in ascending order
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<double>& vdDataPoints);
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs vdSupportingPoints and vdDataPoints. vdSupportingPoints need to be sorted in ascending order
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<float>& vdDataPoints);
//All other permutations of data types
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given supporting and data points. The supporting points need to be sorted in ascending order. BEWARE of the internal conversion to from float double!
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<float>& vfSupportingPoints, const std::vector<float>& vfDataPoints);
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given supporting and data points. The supporting points need to be sorted in ascending order. BEWARE of the internal conversion to from float double!
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<float>& vfSupportingPoints, const std::vector<double>& vdDataPoints);
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given supporting and data points. The supporting points need to be sorted in ascending order. BEWARE of the internal conversion to from float double!
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<float>& vfDataPoints);
//! Uses cubic splines to interpolate the function defined by the given data pairs vdSupportingPoints and vdDataPoints to a vector of points
ITA_BASE_API std::vector<double> CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints, const std::vector<double>& vdEvalPoints);
//! Uses cubic splines to interpolate the function defined by the given data pairs vdSupportingPoints and vdDataPoints to a vector of points
ITA_BASE_API std::vector<double> CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, const std::vector<float>& vdEvalPoints);
//! Uses cubic splines to interpolate the function defined by the given supporting and data points to a vector of evaluation points
ITA_BASE_API std::vector<double> CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints, const std::vector<double>& vdEvalPoints);
//! Uses cubic splines to interpolate the function defined by the given supporting and data points to a vector of evaluation points. BEWARE of the internal conversion to from float double!
ITA_BASE_API std::vector<double> CubicSpline(const std::vector<float>& vfSupportingPoints, const std::vector<float>& vfDataPoints, const std::vector<float>& vfEvalPoints);
//! Uses cubic splines to interpolate the function defined by the given data pairs vdSupportingPoints and vdDataPoints to a single point
ITA_BASE_API double CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints, double dEvalPoint);
//! Uses cubic splines to interpolate the function defined by the given data pairs vdSupportingPoints and vdDataPoints to a single point. BEWARE of the internal conversion to type double!
ITA_BASE_API float CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, float dEvalPoint);
//! Uses cubic splines to interpolate the function defined by the given supporting and data points to a single evaluation points
ITA_BASE_API double CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints, double dEvalPoint);
//! Uses cubic splines to interpolate the function defined by the given supporting and data points to a single evaluation points. BEWARE of the internal conversion to from float double!
ITA_BASE_API float CubicSpline(const std::vector<float>& vfSupportingPoints, const std::vector<float>& vfDataPoints, float dEvalPoint);
}
}
......
......@@ -21,35 +21,17 @@ CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<double>& vdBreakPoi
if ( vdBreakPoints.size() < iMinimumNBreakPoints )
ITA_EXCEPT_INVALID_PARAMETER("Number of break points not sufficient for polynom order: Expected " + std::to_string(iMinimumNBreakPoints) + " but got " + std::to_string(vdBreakPoints.size()) + ".");
const int nCoeffs = NumIntervals() * (iOrder + 1);
if (vdCoeffs.size() != nCoeffs)
ITA_EXCEPT_INVALID_PARAMETER("Invalid number of coefficients: Expected " + std::to_string(nCoeffs) + " but got " + std::to_string(vdCoeffs.size()) + ".");
const int iNumCoeffs = NumIntervals() * (iOrder + 1);
if (vdCoeffs.size() != iNumCoeffs)
ITA_EXCEPT_INVALID_PARAMETER("Invalid number of coefficients: Expected " + std::to_string(iNumCoeffs) + " but got " + std::to_string(vdCoeffs.size()) + ".");
}
CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<float>& vdBreakPnts, const std::vector<float>& vdCoeffs, const int iOrder)
: iOrder(iOrder)
CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<float>& vfBreakPoints, const std::vector<float>& vfCoeffs, const int iOrder)
: CPiecewisePolynomial( std::vector<double>( vfBreakPoints.begin(), vfBreakPoints.end() ), std::vector<double>( vfCoeffs.begin(), vfCoeffs.end() ), iOrder )
{
/*
auto doubles_BreakPoints = std::vector<double>(vdBreakPnts.begin(), vdBreakPnts.end());
auto double_vdCoeffs = std::vector<double>(vdCoeffs.begin(), vdCoeffs.end());
CPiecewisePolynomial(doubles_BreakPoints, double_vdCoeffs, iOrder);
*/
auto doubles_BreakPoints = std::vector<double>(vdBreakPnts.begin(), vdBreakPnts.end());
auto double_vdCoeffs = std::vector<double>(vdCoeffs.begin(), vdCoeffs.end());
vdCoefficients = double_vdCoeffs;
vdBreakPoints = doubles_BreakPoints;
if (iOrder < 0)
ITA_EXCEPT_INVALID_PARAMETER("Polynom order must be greater or equal zero.");
const int iMinimumNBreakPoints = 2; // At least one interval
if (vdBreakPnts.size() < iMinimumNBreakPoints)
ITA_EXCEPT_INVALID_PARAMETER("Number of break points not sufficient for polynom order: Expected " + std::to_string(iMinimumNBreakPoints) + " but got " + std::to_string(vdBreakPoints.size()) + ".");
}
const int nCoeffs = NumIntervals() * (iOrder + 1);
if (vdCoeffs.size() != nCoeffs)
ITA_EXCEPT_INVALID_PARAMETER("Invalid number of coefficients: Expected " + std::to_string(nCoeffs) + " but got " + std::to_string(vdCoeffs.size()) + ".");
}
double CPiecewisePolynomial::Value(const double & xValue, bool bExtrapolate) const
{
//a(x-x0)^(n) + b(x-x0)^(n-1) + c(x-x0)^(n-2) +...
......@@ -64,6 +46,7 @@ double CPiecewisePolynomial::Value(const double & xValue, bool bExtrapolate) con
double fX0 = vdBreakPoints[iInterval];
return ppDerivative.Value(fX0) * ( xValue - fX0) + this->Value(fX0); //y=m(x-x0)+d
}
//inside bounds
std::vector<double>::const_iterator itCoef = CoefficientsOfInterval(iInterval);
const double x0 = vdBreakPoints[iInterval];
......@@ -73,10 +56,6 @@ double CPiecewisePolynomial::Value(const double & xValue, bool bExtrapolate) con
}
return fResult;
}
double CPiecewisePolynomial::Value(const float& xValue, bool bExtrapolate) const
{
return Value(double(xValue), bExtrapolate);
}
std::vector<double> CPiecewisePolynomial::Value(const std::vector<double>& vdXValues, bool bExtrapolate) const
{
std::vector<double> result(vdXValues.size());
......@@ -84,11 +63,15 @@ std::vector<double> CPiecewisePolynomial::Value(const std::vector<double>& vdXVa
result[i] = Value(vdXValues[i], bExtrapolate);
return result;
}
std::vector<double> CPiecewisePolynomial::Value(const std::vector<float>& vdXValues, bool bExtrapolate) const
std::vector<double> CPiecewisePolynomial::Value(const std::vector<float>& vfXValues, bool bExtrapolate) const
{
auto doubles_xValues = std::vector<double>(vdXValues.begin(), vdXValues.end());
return Value(doubles_xValues, bExtrapolate);
std::vector<double> result(vfXValues.size());
for (int i = 0; i < vfXValues.size(); i++)
result[i] = Value(vfXValues[i], bExtrapolate);
return result;
}
int CPiecewisePolynomial::IntervalIndex(const double xValue, bool bExtrapolate) const
{
if (xValue < vdBreakPoints[0] || xValue > vdBreakPoints[NumIntervals()])
......@@ -100,12 +83,12 @@ int CPiecewisePolynomial::IntervalIndex(const double xValue, bool bExtrapolate)
}
return NumIntervals()-1; //last possibility is xValue == last BreakPoint --> xValue in last interval
}
int CPiecewisePolynomial::IntervalIndex(const float xValue, bool bExtrapolate) const
{
return IntervalIndex(double(xValue), bExtrapolate);
}
CPiecewisePolynomial CPiecewisePolynomial::Derivation() const
{
const int iNewOrder = std::max(0, iOrder - 1); // Minimum order is zero
......
......@@ -6,7 +6,11 @@
using namespace ITABase;
using namespace ITABase::Math;
//basic Constructors
//---FUNCTIONS CREATING POLYNOMIALS---
//------------------------------------
CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints)
{
/*setUp the bandmatrix M with the the h values. The vector with the coefficent b is called b.
......@@ -71,7 +75,6 @@ CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<double>& vdSup
//sort for Constructor CPiecewisePolynomial
int iCoefficients = (iPoints - 1) * (iPolynomialOrder + 1);
std::vector<double> vdCoefficients(iCoefficients);
int k = 0;
for (int k = 0; k < iPoints - 1; k++) {
vdCoefficients[4 * k] = vdA[k];
vdCoefficients[4 * k + 1] = vdB[k];
......@@ -81,49 +84,48 @@ CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<double>& vdSup
return CPiecewisePolynomial(vdSupportingPoints, vdCoefficients, iPolynomialOrder);
}
CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints)
CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<float>& vfSupportingPoints, const std::vector<float>& vfDataPoints)
{
//Typecasting upwards to doubles
auto doubles_SupportingPoints = std::vector<double>(vdSupportingPoints.begin(), vdSupportingPoints.end());
auto doubles_DataPoints = std::vector<double>(vdDataPoints.begin(), vdDataPoints.end());
return CubicSpline(doubles_SupportingPoints, doubles_DataPoints);//Constructor working with doubles
auto vdSupportingPoints = std::vector<double>(vfSupportingPoints.begin(), vfSupportingPoints.end());
auto vdDataPoints = std::vector<double>(vfDataPoints.begin(), vfDataPoints.end());
return CubicSpline(vdSupportingPoints, vdDataPoints); //Calling main function working with doubles
}
CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<double>& vdDataPoints)
CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<float>& vfSupportingPoints, const std::vector<double>& vdDataPoints)
{
//Typecasting upwards to doubles
auto doubles_SupportingPoints = std::vector<double>(vdSupportingPoints.begin(), vdSupportingPoints.end());
return CubicSpline(doubles_SupportingPoints, vdDataPoints);//Constructor working with doubles
auto vdSupportingPoints = std::vector<double>(vfSupportingPoints.begin(), vfSupportingPoints.end());
return CubicSpline(vdSupportingPoints, vdDataPoints); //Calling main function working with doubles
}
CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<float>& vdDataPoints)
CPiecewisePolynomial ITABase::Math::CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<float>& vfDataPoints)
{
//Typecasting upwards to doubles
auto doubles_DataPoints = std::vector<double>(vdDataPoints.begin(), vdDataPoints.end());
return CubicSpline(vdSupportingPoints, doubles_DataPoints);//Constructor working with doubles
auto vdDataPoints = std::vector<double>(vfDataPoints.begin(), vfDataPoints.end());
return CubicSpline(vdSupportingPoints, vdDataPoints); //Calling main function working with doubles
}
//Constructors returning evaluated Points
//---FUNCTIONS EVALUATING POLYNOMIALS---
//--------------------------------------
std::vector<double> ITABase::Math::CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints, const std::vector<double>& vdEvalPoints)
{
CPiecewisePolynomial myPoly = CubicSpline(vdSupportingPoints, vdDataPoints);
return myPoly.Value(vdEvalPoints);
return myPoly.Value( vdEvalPoints );
}
std::vector<double> ITABase::Math::CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, const std::vector<float>& vdEvalPoints)
std::vector<double> ITABase::Math::CubicSpline(const std::vector<float>& vfSupportingPoints, const std::vector<float>& vfDataPoints, const std::vector<float>& vfEvalPoints)
{
CPiecewisePolynomial myPoly = CubicSpline(vdSupportingPoints, vdDataPoints);
return myPoly.Value(vdEvalPoints);
CPiecewisePolynomial myPoly = CubicSpline(vfSupportingPoints, vfDataPoints);
return myPoly.Value( vfEvalPoints );
}
double ITABase::Math::CubicSpline(const std::vector<double>& vdSupportingPoints, const std::vector<double>& vdDataPoints, double dEvalPoint)
{
CPiecewisePolynomial myPoly = CubicSpline(vdSupportingPoints, vdDataPoints);
return myPoly.Value(dEvalPoint);
return myPoly.Value( dEvalPoint );
}
float ITABase::Math::CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, float dEvalPoint)
float ITABase::Math::CubicSpline(const std::vector<float>& vfSupportingPoints, const std::vector<float>& vfDataPoints, float fEvalPoint)
{
CPiecewisePolynomial myPoly = CubicSpline(vdSupportingPoints, vdDataPoints);
return float(myPoly.Value(dEvalPoint));
CPiecewisePolynomial myPoly = CubicSpline(vfSupportingPoints, vfDataPoints);
return float( myPoly.Value(fEvalPoint) );
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment