source: trunk/yat/regression/GSLInterpolation.h @ 1724

#ifndef _theplu_yat_regression_gslinterpolation_
#define _theplu_yat_regression_gslinterpolation_

// $Id: GSLInterpolation.h 1724 2009-01-15 14:32:36Z jari $

/*
  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009 Jari Häkkinen

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <gsl/gsl_interp.h>

namespace theplu {
namespace yat {
namespace utility {
  class VectorBase;
}
namespace regression {

  /**
     \brief Base class for interfacing GSL interpolation.

     The GSL interpolation is described in the
     <A HREF=
     "http://www.gnu.org/software/gsl/manual/html_node/Interpolation.html">
     GSL Manual.</A>
     The GSL library provides a variety of interpolation methods,
     including Cubic splines and Akima splines. Interpolations can be
     defined for both normal and periodic boundary
     conditions. Additional functions are available for computing
     derivatives and integrals of interpolating functions.

     Given a set of data points \f$ (x_1, y_1) \dots (x_n, y_n) \f$
     the sub classes compute a continuous interpolating function \f$
     y(x) \f$ such that \f$ y(x_i) = y_i \f$. The interpolation is
     piecewise smooth, and its behavior at the end-points is
     determined by the type of interpolation used.

     When underlying GSL functions return GSL error an GSL_error is
     thrown. Refer to the gsl_errno.h for the error code listing.

     \since New in yat 0.5
  */
  class GSLInterpolation
  {

  public:
    /**
       \brief The default constructor

       Initialization of the interpolation object for the data \f$ (x,
       y) \f$ where \a x and \a y are vector like objects of the same
       size. The content of \a x and \a y are copied for internal
       storage. \a x is always assumed to be strictly ordered, with
       increasing \a x values; the behavior for other arrangements is
       not defined.

       \throw GSL_error if some underlying GSL functions return GSL
       error. Refer to gsl_errno.h for the error code listing.
    */
    GSLInterpolation(const gsl_interp_type*, const utility::VectorBase& x,
                     const utility::VectorBase& y);

    /**
       \brief The destructor
    */
    virtual ~GSLInterpolation(void)=0;

    /**
       \brief Calculate the interpolated value for \a x.

       \return The interpolated value of \f$ y \f$ for a given point
       \a x.

       \throw GSL_error if evaluation is requested outside the range
       defined by the interpolation algorithm.
    */
    double evaluate(double x);

    /**
       \brief Calculate the derivative of the interpolated function at
       \a x.

       \return The derivative.

       \throw GSL_error if evaluation is requested outside the range
       defined by the interpolation algorithm.
    */
    double evaluate_derivative(double x);

    /**
       \brief Calculate the 2nd derivative of the interpolated
       function at \a x.

       \return The 2nd derivative.

       \throw GSL_error if evaluation is requested outside the range
       defined by the interpolation algorithm.
    */
    double evaluate_derivative2(double x);

    /**
       \brief Calculate the numerical integral of the interpolated
       function over the range \f$ [a,b] \f$.

       \return The integral.

       \throw GSL_error if evaluation is requested outside the range
       defined by the interpolation algorithm.
    */
    double evaluate_integral(double a, double b);

    /**
       \brief The result of the latest evaluaion function call is
       stored and can be retrieved with this function.

       \return The latest evaluated value.
     */
    double evaluation(void) const;

    /**
       \brief This function returns the minimum number of points
       required by the interpolation type.

       For example, Akima spline interpolation requires a minimum of 5
       points.

       \return The minimum number of points required.
    */
    unsigned int min_size(void) const;

  private:
    /**
       \brief Copy Constructor. (not implemented)
    */
    GSLInterpolation(const GSLInterpolation&);
    // no assignment
    GSLInterpolation& operator=(const GSLInterpolation&);

    gsl_interp_accel* accelerator_;
    double evaluation_;
    gsl_interp* interpolator_;
    double* x_;
    double* y_;
  };

}}} // of namespaces regression, yat, and theplu

#endif