Commit 2c1d051d authored by pia.stammer's avatar pia.stammer
Browse files

Added documentation for interpolation class

parent 34c6676b
Pipeline #119880 canceled with stages
/*!
* @file interpolation.h
* @brief Class for interpolation of tabulated values (used for database interpolation)
*/
#ifndef INTERPOLATION_H
#define INTERPOLATION_H
......@@ -16,18 +20,73 @@ class Interpolation
Interpolation() = delete;
/*!
* @brief Check if vector lengths match and values are sorted ascendingly
*/
void Setup();
/*!
* @return index of element closest to 'value' in 'v'
*/
unsigned IndexOfClosestValue( double value, const Vector& v ) const;
/*!
* @return value of third degree/cubic polynomial with parameters 'param' at 'x'
*/
inline double EvalCubic1DSpline( double param[4], double x ) const;
public:
//In constructor tabulated values are initialised
/*!
* @brief constructor linear interpolation for std::vector input
* @param[in] x - table values for x
* @param[in] y - table values for y
* @param[in] type - linear or cubic interpolation
*/
Interpolation( const std::vector<double>& x, const std::vector<double>& y, TYPE type = linear );
/*!
* @brief constructor linear interpolation for Vector input
* @param[in] x - table values for x
* @param[in] y - table values for y
* @param[in] type - linear or cubic interpolation
*/
Interpolation( const Vector& x, const Vector& y, TYPE type = linear );
/*!
* @brief constructor cubic interpolation
* @param[in] x - table values for x
* @param[in] y - table values for y
* @param[in] type - linear or cubic interpolation
*/
Interpolation( const Vector& x, const Vector& y, const Matrix& data, TYPE type = cubic );
//Here the interpolation between the previously defined table values is performed
/*!
* @brief defines one dimensional interpolation at x
* @param[in] x - value at which to interpolate
* @param[out] y - corresponding interpolated y value
*/
double operator()( double x ) const;
/*!
* @brief defines interpolation for a Vector of values
* @param[in] v - values at which to interpolate
* @param[out] y - corresponding interpolated values
*/
Vector operator()( Vector v ) const;
/*!
* @brief defines interpolation for a std::vector of values
* @param[in] v - values at which to interpolate
* @param[out] y - corresponding interpolated values
*/
std::vector<double> operator()( std::vector<double> v ) const;
/*!
* @brief defines 2D interpolation at x and y
* @param[in] x - value at which to interpolate
* @param[in] y - value at which to interpolate
* @param[out] data - corresponding interpolated value
*/
double operator()( double x, double y ) const;
};
......
......@@ -2,31 +2,43 @@
#include "toolboxes/errormessages.h"
#include <blaze/math/lapack/posv.h>
// Change Vector type to blaze for std input
Interpolation::Interpolation( const std::vector<double>& x, const std::vector<double>& y, TYPE type ) : _dim( 1u ), _type( type ) {
_x = Vector( x.size(), x.data() );
_y = Vector( y.size(), y.data() );
Setup();
}
// 1D Constructor: Initialise values and call Set-up
Interpolation::Interpolation( const Vector& x, const Vector& y, TYPE type ) : _dim( 1u ), _x( x ), _y( y ), _type( type ) { Setup(); }
// 2D Constructor: Initialise values and call Set-up
Interpolation::Interpolation( const Vector& x, const Vector& y, const Matrix& data, TYPE type )
: _dim( 2u ), _x( x ), _y( y ), _data( data ), _type( type ) {
Setup();
}
// Set-up: check validity of input
void Interpolation::Setup() {
// 1D
if( _dim == 1u ) {
// Length of tables must be equal
if( _x.size() != _y.size() ) ErrorMessages::Error( "Vectors are of unequal length!", CURRENT_FUNCTION );
if( _type == loglinear ) _y = blaze::log( _y );
// Values must be sorted ascendingly (property is used later when searching tables)
for( unsigned i = 0; i < _x.size() - 1u; i++ ) {
if( !( _x[i] < _x[i + 1] ) ) ErrorMessages::Error( "x is not sorted ascendingly!", CURRENT_FUNCTION );
if( !( _x[i] <_x[i + 1] ) ) ErrorMessages::Error( "x is not sorted ascendingly!", CURRENT_FUNCTION );
}
}
// 2D
else if( _dim == 2u ) {
// Length of tables must be equal
if( _x.size() != _data.rows() ) ErrorMessages::Error( "x and data are of unequal length!", CURRENT_FUNCTION );
if( _y.size() != _data.columns() ) ErrorMessages::Error( "y and data are of unequal length!", CURRENT_FUNCTION );
// Values must be sorted ascendingly (property is used later when searching tables)
for( unsigned i = 0; i < _x.size() - 1u; i++ ) {
if( !( _x[i] < _x[i + 1] ) ) ErrorMessages::Error( "x is not sorted ascendingly!", CURRENT_FUNCTION );
}
......@@ -36,15 +48,22 @@ void Interpolation::Setup() {
}
}
// 1D interpolation
double Interpolation::operator()( double x ) const {
// Check whether 1D
if( _dim != 1u ) ErrorMessages::Error( "Invalid data dimension for operator(x)!", CURRENT_FUNCTION );
// x must be between min and max of table values
if( x < _x[0] || x > _x[_x.size() - 1u] ) ErrorMessages::Error( "Extrapolation is not supported!", CURRENT_FUNCTION );
// points to first value in _x that is not smaller than x (upper bound)
Vector::ConstIterator it = std::lower_bound( _x.begin(), _x.end(), x );
// index of the lower bound to x in _x
unsigned idx = static_cast<unsigned>( std::max( int( it - _x.begin() ) - 1, 0 ) );
// linear interpolation
if( _type == linear || _type == loglinear ) {
// interpolate between lower bound and upper bound of x in _x
double res = _y[idx] + ( _y[idx + 1] - _y[idx] ) / ( _x[idx + 1] - _x[idx] ) * ( x - _x[idx] );
if( _type == loglinear )
......@@ -52,6 +71,8 @@ double Interpolation::operator()( double x ) const {
else
return res;
}
// cubic interpolation
else if( _type == cubic ) {
double param[4] = { _y[idx - 1], _y[idx], _y[idx + 1], _y[idx + 2] };
double t = ( x - _x[idx] ) / ( _x[idx + 1] - _x[idx] );
......@@ -63,14 +84,17 @@ double Interpolation::operator()( double x ) const {
}
}
// 2D interpolation
double Interpolation::operator()( double x, double y ) const {
// Check whether 2D
if( _dim != 2u ) ErrorMessages::Error( "Invalid data dimension for operator(x,y)!", CURRENT_FUNCTION );
if( _type == cubic ) {
// find closest values to x and y in table (lower bounds)
unsigned xId = IndexOfClosestValue( x, _x );
unsigned yId = IndexOfClosestValue( y, _y );
// store all 16 interpolation points needed
// store all 16 interpolation points needed
double points[4][4];
for( int i = -1; i < 3; ++i ) {
unsigned idx_y;
......@@ -103,6 +127,7 @@ double Interpolation::operator()( double x, double y ) const {
}
}
// extension of 1D interpolation to Vector input
Vector Interpolation::operator()( Vector v ) const {
Vector res( v.size() );
for( unsigned i = 0; i < v.size(); ++i ) {
......@@ -111,6 +136,7 @@ Vector Interpolation::operator()( Vector v ) const {
return res;
}
// extension of 1D interpolation to Vector input for std::vectors
std::vector<double> Interpolation::operator()( std::vector<double> v ) const {
std::vector<double> res( v.size() );
for( unsigned i = 0; i < v.size(); ++i ) {
......@@ -119,6 +145,7 @@ std::vector<double> Interpolation::operator()( std::vector<double> v ) const {
return res;
}
// implementation of third degree polynomial
inline double Interpolation::EvalCubic1DSpline( double param[4], double x ) const {
return ( param[1] + 0.5 * x *
( param[2] - param[0] +
......@@ -126,6 +153,8 @@ inline double Interpolation::EvalCubic1DSpline( double param[4], double x ) cons
x * ( 3.0 * ( param[1] - param[2] ) + param[3] - param[0] ) ) ) );
}
//find index of closes value to 'value' in vector 'v'
unsigned Interpolation::IndexOfClosestValue( double value, const Vector& v ) const {
auto i = std::min_element( begin( v ), end( v ), [=]( double x, double y ) { return abs( x - value ) < abs( y - value ); } );
return std::distance( begin( v ), i );
......
Supports Markdown
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