/home/mhx/git/github/libmoost/include/moost/algorithm/spline_interpolation.hpp

Go to the documentation of this file.

/* vim:set ts=3 sw=3 sts=3 et: */
#ifndef MOOST_ALGORITHM_SPLINE_INTERPOLATION_HPP
#define MOOST_ALGORITHM_SPLINE_INTERPOLATION_HPP

#include <cassert>
#include <vector>
#include <algorithm>

#include <gsl/gsl_spline.h> // gnu scientific library

#include <boost/noncopyable.hpp>

#include "../utils/scope_exit.hpp"


namespace moost { namespace algorithm {

   class spline_interpolation : private boost::noncopyable
   {
      private:
         typedef moost::utils::scope_exit::type<
            gsl_interp *>::call_free_function_with_val interp_t;

         typedef moost::utils::scope_exit::type<
            gsl_interp_accel *>::call_free_function_with_val accel_t;

      public:


         spline_interpolation(
               std::vector<double> const & x,
               std::vector<double> const & y)
            : interp_(gsl_interp_alloc(gsl_interp_cspline, x.size()), &gsl_interp_free)
              , accel_(gsl_interp_accel_alloc(), &gsl_interp_accel_free)
              , x_(x)
              , y_(y)

         {
            validate_construction(); // throws on fail
            gsl_interp_init(interp_->get(), &x_[0], &y_[0], x_.size());
         }


         bool operator () (double const x, double & y, double const z = 0) const
         {
            bool ok =  0 == gsl_interp_eval_e(
                  interp_->get(),
                  &x_[0], &y_[0], x,
                  accel_->get(), &y
                  );

            if(!ok) { y = z;}

            return ok;
         }

         double operator () (double const x) const
         {
            double y = 0.0;

            if(!(*this)(x, y))
            {
               throw std::range_error("out of range");
            }

            return y;
         }

      private:

         void validate_construction() const
         {
            // x and y must be the same size
            assert(x_.size() == y_.size());
            if(x_.size() != y_.size())
            {
               throw std::runtime_error("the size of x and y must be the same");
            }

            // x must be sorted by-asc (not my rule, gsl demands it)
            bool const sorted = (std::adjacent_find(x_.begin(), x_.end(), std::greater_equal<double>()) == x_.end());
            assert(sorted);
            if(!sorted)
            {
               throw std::runtime_error("x must be sorted in assending order");
            }
         }

      private:
         interp_t interp_;
         accel_t  accel_;
         std::vector<double> const & x_;
         std::vector<double> const & y_;
   };


}}

#endif // MOOST_ALGORITHM_SPLINE_INTERPOLATION_HPP