Approximates a 1D function using Chebyshev points. More...

#include <ChebyshevInterpolatingFunction.hh>

Inheritance diagram for marley::ChebyshevInterpolatingFunction:

Public Member Functions
	ChebyshevInterpolatingFunction (const std::function< double(double)> &func, double x_min, double x_max, size_t N=0)

ChebyshevInterpolatingFunction	cdf () const

const std::vector< double > &	chebyshev_coeffs () const

double	evaluate (double x) const override
	Approximates the represented function using the barycentric formula.

const std::vector< double > &	Fs () const

double	integral () const

int	N () const

double	x_max () const

double	x_min () const

const std::vector< double > &	Xs () const

Public Member Functions inherited from marley::InterpolatingFunction
	InterpolatingFunction (double x_min, double x_max)

Protected Member Functions
	ChebyshevInterpolatingFunction ()
	Default constructor used by cdf()

double	chebyshev_point (size_t j) const
	For a given N, returns the x position of the jth Chebyshev point (of the second kind)

void	compute_integral ()

Protected Member Functions inherited from marley::InterpolatingFunction
	InterpolatingFunction ()
	Allows construction by derived classes without immediately defining the bounds of the interpolation range.

Protected Attributes
std::vector< double >	chebyshev_coeffs_
	Coefficients of the Chebyshev expansion of this function.

std::vector< double >	Fs_
	Function values at the grid points.

std::vector< int >	ifac_

double	integral_
	Approximate integral of the function on [x_min_, x_max_].

size_t	N_
	Grid size parameter (N + 1 total points)

std::vector< double >	wsave_

std::vector< double >	Xs_
	Chebyshev points at which the function was evaluated.

Protected Attributes inherited from marley::InterpolatingFunction
double	x_max_
	Upper bound of the interpolation range.

double	x_min_
	Lower bound of the interpolation range.

Static Protected Attributes
static constexpr size_t	N_MAX_ = 65536u
	Maximum allowed value of the grid size parameter.

Detailed Description

Approximates a 1D function using Chebyshev points.

Definition at line 29 of file ChebyshevInterpolatingFunction.hh.

Constructor & Destructor Documentation

◆ ChebyshevInterpolatingFunction() [1/2]

marley::ChebyshevInterpolatingFunction::ChebyshevInterpolatingFunction	(	const std::function< double(double)> &	func,
		double	x_min,
		double	x_max,
		size_t	N = 0 )

Todo: PRINT WARNING MESSAGE

Definition at line 34 of file ChebyshevInterpolatingFunction.cc.

  : marley::InterpolatingFunction( x_min, x_max )
{
  bool ok;
 
  if ( N != 0 ) {
    ok = true;
    N_ = N;
  }
  else {
    ok = false;
    N_ = 1;
  }
 
  // Adaptively find a good grid size
  // @todo add more explanation
  do {
 
    if ( !ok ) N_ *= 2;
 
    //std::cout << "N = " << N_ << '\n';
 
    // @todo optimize this
    Xs_.clear();
    Fs_.clear();
 
    for ( size_t j = 0; j <= N_; ++j ) {
      double x = chebyshev_point( j );
      double f = func(x);
      Xs_.push_back( x );
      Fs_.push_back( f );
    }
 
    chebyshev_coeffs_ = Fs_;
    wsave_ = std::vector<double>(3*Fs_.size() + 15, 0.);
    ifac_ = std::vector<int>(Fs_.size() / 2, 0);
 
    int my_size = Fs_.size();
    costi( &my_size, wsave_.data(), ifac_.data() );
    cost( &my_size, chebyshev_coeffs_.data(), wsave_.data(), ifac_.data() );
 
    double biggest_coeff = *std::max_element(chebyshev_coeffs_.cbegin(),
      chebyshev_coeffs_.cend(), [](double left, double right) -> double
      { return std::abs(left) < std::abs(right); });
 
    double upper_limit = 2. * std::abs( biggest_coeff ) * MY_EPSILON;
 
    double last_coeff_mag = std::abs( chebyshev_coeffs_.back() );
    double next_to_last_coeff_mag = std::abs(
      chebyshev_coeffs_.at( chebyshev_coeffs_.size() - 2 ));
    if ( last_coeff_mag < upper_limit && next_to_last_coeff_mag < upper_limit )
    {
      ok = true;
    }
 
    if ( N_ >= N_MAX_ ) {
      ok = true;
    }
 
  } while ( !ok );
 
  // Normalize the Chebyshev coefficients by dividing by N
  for (double& d : chebyshev_coeffs_) d /= N_;
 
  compute_integral();
}

References chebyshev_coeffs_, chebyshev_point(), Fs_, N_, N_MAX_, and Xs_.

◆ ChebyshevInterpolatingFunction() [2/2]

marley::ChebyshevInterpolatingFunction::ChebyshevInterpolatingFunction ( )

inlineprotected

Default constructor used by cdf()

Definition at line 85 of file ChebyshevInterpolatingFunction.hh.

85{};

Member Function Documentation

◆ cdf()

marley::ChebyshevInterpolatingFunction marley::ChebyshevInterpolatingFunction::cdf ( ) const

Definition at line 122 of file ChebyshevInterpolatingFunction.cc.

{
  marley::ChebyshevInterpolatingFunction result;
  result.N_ = this->N_ + 1;
  result.x_min_ = this->x_min_;
  result.x_max_ = this->x_max_;
 
  double beta_0 = this->chebyshev_coeffs_.at(0);
  beta_0 += -0.5 * this->chebyshev_coeffs_.at(1);
  for ( size_t j = 2; j <= this->N_; ++j ) {
    beta_0 += 2.*this->chebyshev_coeffs_.at(j)
      * std::pow(-1.0, j + 1) / (j*j - 1.0);
  }
 
  result.chebyshev_coeffs_.push_back( beta_0 );
  for ( size_t k = 1; k < this->N_; ++k ) {
    double ak_minus_one = this->chebyshev_coeffs_.at( k - 1 );
    double ak_plus_one = this->chebyshev_coeffs_.at( k + 1 );
    result.chebyshev_coeffs_.push_back(
      (ak_minus_one - ak_plus_one) / (2 * k) );
  }
  result.chebyshev_coeffs_.push_back( this->chebyshev_coeffs_.at( N_ - 1 )
    / (2 * N_) );
  result.chebyshev_coeffs_.push_back( this->chebyshev_coeffs_.at( N_ )
    / (2 * (N_ + 1)) );
 
  for (double& d : result.chebyshev_coeffs_) d *= (x_max_ - x_min_) / 2.;
 
  result.compute_integral();
 
  for ( size_t j = 0; j <= result.N_; ++j ) {
    double x = result.chebyshev_point( j );
    result.Xs_.push_back( x );
  }
 
  result.wsave_ = std::vector<double>(
    3*result.chebyshev_coeffs_.size() + 15, 0.);
  result.ifac_ = std::vector<int>(result.chebyshev_coeffs_.size() / 2, 0);
 
  int my_size = result.chebyshev_coeffs_.size();
  result.Fs_ = result.chebyshev_coeffs_;
  costi( &my_size, result.wsave_.data(), result.ifac_.data() );
  cost( &my_size, result.Fs_.data(), result.wsave_.data(), result.ifac_.data() );
  for (double& d : result.Fs_) d *= 0.5;
 
  return result;
}

◆ chebyshev_coeffs()

const std::vector< double > & marley::ChebyshevInterpolatingFunction::chebyshev_coeffs ( ) const

inline

Definition at line 66 of file ChebyshevInterpolatingFunction.hh.

67 { return chebyshev_coeffs_; }

◆ chebyshev_point()

double marley::ChebyshevInterpolatingFunction::chebyshev_point ( size_t j ) const

inlineprotected

For a given N, returns the x position of the jth Chebyshev point (of the second kind)

Todo: If needed for speed, consider caching the std::cos evaluations used here

Definition at line 113 of file ChebyshevInterpolatingFunction.hh.

                                                    {
        // jth Chebyshev point (of the second kind) on [-1, 1]
        double x = std::cos( static_cast<double>(marley_utils::pi * j) / N_ );
        // Affine transformation to [x_min_, x_max_] (primed basis)
        double x_prime = (( x_max_ - x_min_ )*x + ( x_max_ + x_min_ )) / 2.;
        return x_prime;
      }

References N_, marley::InterpolatingFunction::x_max_, and marley::InterpolatingFunction::x_min_.

Referenced by ChebyshevInterpolatingFunction().

◆ compute_integral()

void marley::ChebyshevInterpolatingFunction::compute_integral ( )

protected

Definition at line 104 of file ChebyshevInterpolatingFunction.cc.

                                                            {
  // Compute the integral over [x_min_, x_max_] via Clenshaw-Curtis
  // quadrature
  integral_ = 0.;
  for (size_t k = 0; k <= N_; k += 2) {
    // Only even Chebyshev polynomials contribute to this integral
    double term = chebyshev_coeffs_.at( k ) * 2.0 / ( 1. - std::pow(k, 2) );
    if ( k == 0 || k == N_ ) term /= 2.;
    integral_ += term;
  }
 
  // Extra factor to account for integration over [x_min_, x_max_] instead
  // of the standard interval [-1, 1] for the Chebyshev polynomials
  integral_ *= (x_max_ - x_min_) / 2.;
}

◆ evaluate()

double marley::ChebyshevInterpolatingFunction::evaluate ( double x ) const

inlineoverridevirtual

Approximates the represented function using the barycentric formula.

Implements marley::InterpolatingFunction.

Definition at line 39 of file ChebyshevInterpolatingFunction.hh.

                                                        {
        double numer = 0.;
        double denom = 0.;
        double w_j = -1.;
        for ( size_t j = 0; j <= N_; ++j ) {
          w_j = -w_j; // w_j = (-1)^(j)
          double x_j = Xs_.at( j );
          double f_j = Fs_.at( j );
          // If the requested x value is exactly equal to one
          // of the grid points where we previously evaluated the
          // function, then just return that
          if ( x_j == x ) return f_j;
 
          double temp = w_j / ( x - x_j );
          if ( j == 0 || j == N_ ) temp /= 2.;
          denom += temp;
          numer += temp * f_j;
        }
 
        double px = numer / denom;
        return px;
      }