JPolfit.hh

Go to the documentation of this file.

#ifndef __JTOOLS__JPOLFIT__
#define __JTOOLS__JPOLFIT__
 
#include <cmath>
#include <iterator>
#include <algorithm>
#include <vector>
#include <utility>
#include <sstream>
 
#include "JLang/JException.hh"
#include "JLang/JAssert.hh"
#include "JLang/JStreamAvailable.hh"
#include "JMath/JMathSupportkit.hh"
#include "JMath/JLegendre.hh"
#include "JMath/JZero.hh"
#include "JTools/JFunctional.hh"
#include "JTools/JDistance.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JTOOLS {}
namespace JPP { using namespace JTOOLS; }
 
namespace JTOOLS {
 
  using JLANG::JException;
  using JLANG::JFunctionalException;
  using JLANG::JNumericalPrecision;
  using JLANG::JValueOutOfRange;
 
 
  /**
   * Functional collection with Legendre polynomial interpolation.
   */
  template<unsigned int N,      // number of points to interpolate
           unsigned int M,      // degree of polynomial
           class JElement_t,
           template<class, class> class JCollection_t,
           class JResult_t,
           class JDistance_t>
  class JPolfitFunction :
    public JCollection_t<JElement_t, JDistance_t>,
    public virtual JFunctional<typename JElement_t::abscissa_type,
                               typename JResultType<typename JElement_t::ordinate_type>::result_type>
  {
  public:
 
    typedef JCollection_t<JElement_t, JDistance_t>                               collection_type;
 
    typedef typename collection_type::abscissa_type                              abscissa_type;
    typedef typename collection_type::ordinate_type                              ordinate_type;
    typedef typename collection_type::value_type                                 value_type;
    typedef typename collection_type::distance_type                              distance_type;
 
    typedef typename collection_type::const_iterator                             const_iterator;
    typedef typename collection_type::const_reverse_iterator                     const_reverse_iterator;
    typedef typename collection_type::iterator                                   iterator;
    typedef typename collection_type::reverse_iterator                           reverse_iterator;
 
    typedef typename JResultType<ordinate_type>::result_type                     data_type;
    typedef JFunctional<abscissa_type, data_type>                                  function_type;           
 
    typedef typename function_type::argument_type                                argument_type;
    typedef typename function_type::result_type                                  result_type;
    typedef typename function_type::JExceptionHandler                            exceptionhandler_type;
 
 
    /**
     * Default constructor.                     
     */
    JPolfitFunction()
    {
      STATIC_CHECK((N > M));
    }
 
 
    /**
     * Recursive interpolation method implementation.
     *
     * \param  pX              pointer to abscissa values
     * \return                 function value
     */
    virtual result_type evaluate(const argument_type* pX) const override 
    {
      using namespace std;
      using namespace JPP;
 
      const argument_type x = *pX;
 
      if (this->size() <= 1u) {
 
        try {
          return this->getExceptionHandler().action();
        }
        catch (const JException& error) {
 
          std::ostringstream os;
 
          os << __FILE__ << ':' << __LINE__ << " not enough data " << STREAM("?") << x;
        
          throw JFunctionalException(os.str());
        }
      }
 
      const_iterator p = this->lower_bound(x);
 
      if ((p == this->begin() && this->getDistance(x, (p++)->getX()) > distance_type::precision) ||
          (p == this->end()   && this->getDistance((--p)->getX(), x) > distance_type::precision)) {
 
        try {
          return this->getExceptionHandler().action();
        }
        catch (const JException& error) {
 
          std::ostringstream os;
 
          os << __FILE__ << ':' << __LINE__ << " abscissa out of range " 
             << STREAM("?") << x                      << " <> " 
             << STREAM("?") << this->begin() ->getX() << ' '
             << STREAM("?") << this->rbegin()->getX();
          
          throw JValueOutOfRange(os.str());
        }
      }
 
      ++pX;  // next argument value
 
 
      int n = std::min((int) (N + 1), (int) this->size());               // number of points to interpolate      
 
      for (int i = n/2; i != 0 && p != this->end();   --i, ++p) {}       // move p to begin of data
      for (int i = n  ; i != 0 && p != this->begin(); --i, --p) {}
 
      buffer.clear();
 
      for (const_iterator q = p; n != 0; ++q, --n) {
        buffer.push_back(make_pair(this->getDistance(p->getX(), q->getX()), function_type::getValue(q->getY(), pX)));
      }
 
      const JLegendre<result_type, M> g1(buffer.begin(), buffer.end());
 
      return g1(this->getDistance(p->getX(), x));
    }
 
  private: 
    /**
     * Function compilation.
     */
    virtual void do_compile() override 
    {}
 
    mutable std::vector< std::pair<abscissa_type, result_type> > buffer;
  };
 
 
  /**
   * Template class for polynomial interpolation in 1D
   *
   * This class implements the JFunction1D interface.
   */
  template<unsigned int N,
           unsigned int M,
           class JElement_t,
           template<class, class> class JCollection_t,
           class JResult_t   = typename JElement_t::ordinate_type,
           class JDistance_t = JDistance<typename JElement_t::abscissa_type> >
  class JPolfitFunction1D :
    public JPolfitFunction<N, M, JElement_t, JCollection_t, JResult_t, JDistance_t>,
    public JFunction1D<typename JElement_t::abscissa_type, JResult_t> 
  {
  public:
 
    typedef JCollection_t<JElement_t, JDistance_t>                               collection_type;
 
    typedef typename collection_type::abscissa_type                              abscissa_type;
    typedef typename collection_type::ordinate_type                              ordinate_type;
    typedef typename collection_type::value_type                                 value_type;
    typedef typename collection_type::distance_type                              distance_type;
 
    typedef typename collection_type::const_iterator                             const_iterator;
    typedef typename collection_type::const_reverse_iterator                     const_reverse_iterator;
    typedef typename collection_type::iterator                                   iterator;
    typedef typename collection_type::reverse_iterator                           reverse_iterator;
 
    typedef JFunction1D<abscissa_type, JResult_t>                                function_type;           
 
    typedef typename function_type::argument_type                                argument_type;
    typedef typename function_type::result_type                                  result_type;
    typedef typename function_type::JExceptionHandler                            exceptionhandler_type;
 
 
    /**
     * Default contructor.
     */
    JPolfitFunction1D()
    {}
  };
}
 
#endif