ChebyshevInterpolatingFunction.hh

//
// This file is part of MARLEY (Model of Argon Reaction Low Energy Yields)
//
// MARLEY is free software: you can redistribute it and/or modify it under the
// terms of version 3 of the GNU General Public License as published by the
// Free Software Foundation.
//
// For the full text of the license please see COPYING or
// visit http://opensource.org/licenses/GPL-3.0
//
// Please respect the MCnet academic usage guidelines. See GUIDELINES
// or visit https://www.montecarlonet.org/GUIDELINES for details.
 
#pragma once
#include <functional>
#include <vector>
 
#include "marley/InterpolatingFunction.hh"
 
namespace marley {
 
  // Default Chebyshev grid size to use (when not using adaptive sizing)
  constexpr size_t DEFAULT_N_CHEBYSHEV = 64u;

  class ChebyshevInterpolatingFunction : public InterpolatingFunction {
 
    public:
 
      // N = 0 triggers adaptive grid sizing
      ChebyshevInterpolatingFunction( const std::function<double(double)>& func,
        double x_min, double x_max, size_t N = 0 );

      inline double evaluate( double x ) const override {
        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;
      }
 
      // @brief Returns the integral of this function on the interval [x_min_,
      // x_max_]
      inline double integral() const { return integral_; };
 
      inline const std::vector<double>& chebyshev_coeffs() const
        { return chebyshev_coeffs_; }
 
      inline const std::vector<double>& Fs() const
        { return Fs_; }
 
      inline const std::vector<double>& Xs() const
        { return Xs_; }
 
      inline int N() const { return N_; }
 
      ChebyshevInterpolatingFunction cdf() const;
 
      inline double x_min() const { return x_min_; }
      inline double x_max() const { return x_max_; }
 
    protected:

      inline ChebyshevInterpolatingFunction() {};

      static constexpr size_t N_MAX_ = 65536u; // 16 adaptive iterations

      size_t N_;

      double integral_;

      std::vector<double> Xs_;

      std::vector<double> Fs_;

      std::vector<double> chebyshev_coeffs_;
 
      // Helper arrays for discrete cosine transforms calculated using FFTPACK4
      std::vector<double> wsave_;
      std::vector<int> ifac_;

      inline double chebyshev_point(size_t j) const {
        // 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;
      }
 
      void compute_integral();
  };
 
}