#ifndef __INTERPOLATOR_HPP__
#define __INTERPOLATOR_HPP__

#include <vector>
#include <map>
#include <boost/shared_ptr.hpp>
#include "bvp.hpp"
#include "linalg.hpp"
#include "func.hpp"
#include "diff_op.hpp"

namespace bvp{
  using std::map;
  using boost::shared_ptr;

  class interp_values;

  template<typename RBF>
  class interpolator : public realfunc{
  public:

    /** @name Constructors
     * 
     * Constructors that take interpolation data perform the
     * interpolation as part of the initialisation, which includes
     * factoring a matrix
     */
    //@{
    ///Does not initialise the interpolator.
    interpolator();
    
    ///Interpolate given a BVP
    interpolator(shared_ptr<linear_BVP2> bvp);

    ///Interpolate given some data points and the value at those points
    interpolator(const map<point, double>& Xi);

    /** \short Interpolate given a domain, data points and the value
     *         at those points
     *
     *  Must provide domain information. The values of Xi must match
     *  points on the given domain Omega.
     */
    interpolator(shared_ptr<domain> Omega, const map<point, double>& Xi);
    //@}

    /** @name Interpolation
     *
     * Interpolate again either given new data or a different BVP.
     */
    //@{
    void interpolate(const map<point, double>& Xi); 
    void interpolate(shared_ptr<linear_BVP2> bvp);
    
    /** \short Interpolate given values on the domain
     * 
     * The given values must match the values returned by other
     * interpolators on compatible domains.
     */
    void interpolate(const interp_values& values); 
    //@}
    

    /** @name Pointwise evaluation and derivatives
     */
    //@{
    ///Point evaluation
    double operator()(const point& p) const;
    ///Point evaluation
    double at(const point& p) const;
    /// Pointwise first derivative
    double d(const point& p, size_t k) const;
    /// Pointwise second derivatives
    double d2(const point &p, size_t k1, size_t k2) const;
    //@}
    
    /** @name Full domain evaluation and derivatives
     */
    //@{
    ///Full domain evaluation
    interp_values operator()() const; 
    ///Full domain evaluation
    interp_values at() const; 
    /// Full domain first erivative
    interp_values d(size_t k) const; 
    /// Full domain second derivatives
    interp_values d2(size_t k1, size_t k2) const; 
    //@}


    /** @name Precomputation of RBFs
     *
     * The following functions will precompute upon request all the
     * interpolator's RBFs at all points of the domain. This may speed
     * up some computations at the expense of increasing memory
     * usage.
     */
    //@{
    ///Precompute all the function evaluation of the RBFs.
    void precompute_ev();
    ///Precompute all relevant first derivatives.
    void precompute_d1() ;
    ///Precompute all relevant second derivatives.
    void precompute_d2();
    //@}

    /** @name Partial redefinitions
     *
     * These functions allow for partial redefinition of the BVP as
     * equired for the additive Schwartz domain decomposition method,
     * and for other methods. They do not factor a matrix again.
     */
    //@{
    void set_f(const realfunc &f);
    void set_g(const realfunc &g);
    void set_f(const map<point, double>& f);
    void set_g(const map<point, double>& g);
    //@}

    /** @name Linear arithmetic operators
     *
     *  These functions return a new interpolator. They are pointwise
     *  linear operations and work by manipulating interpolator
     *  coefficients.
     */
    //@{
    /// Needs two operators on the same domain.
    interpolator<RBF> operator+(const interpolator<RBF>& u) const;
    /// Needs two operators on the same domain.
    interpolator<RBF> operator-(const interpolator<RBF>& u) const;
    /// Returns a scaled interpolator.
    interpolator<RBF> operator*(double a) const;
    /// Returns a scaled interpolator.
    interpolator<RBF> operator/(double a) const;
    //@}
    
  private:
    ///Once the matrix is defined, this function inverts it.
    void computecoeffs(); 
	 
    ///Perform the actual interpolation.
    void init(shared_ptr<linear_BVP2> bvp);

    shared_ptr<linear_BVP2> thebvp;

    ///Number of interior points. 
    size_t n; 
    ///Number of boundary points.
    size_t m; 

    ///The matrix to invert.
    matrix M; 

    ///Is the interpolator ready for use?
    bool initted;

    /** @name Exception throwers
     */
    //@{
    void not_initted(int line, string file) const; 
    void not_precomputed(int line, string file) const;
    //@}
    
    ///Coefficients of the RBFs
    linalg::vector coeffs;
    ///The RBFS
    std::vector<RBF> rbfs;
    ///Their hash
    size_t rbfs_hash;    

    ///An RBF hasher
    size_t hash_value(const std::vector<RBF>& rbfs_in);
       
    /** \short Precomputed RBFs 
     * 
     * For all points on the domain, this stores all the RBFs
     * evaluated at those points, as well as the derivatives on where
     * the RBFs have been evaluated. This is to speed up successive
     * evaluations of the interpolator when the interpolator's domain
     * doesn't change.
     *
     * The keys in this map are vectors that represent the multindex
     * for the derivative, where missing trailing entries in the
     * vector represents zeros. Thus, an empty vector represents
     * evaluation instead of derivatives.
     */
    mutable map<std::vector<size_t>, matrix> precomp_rbfs;
    
    /** @name  Evaluation flags
     *Have the RBFs or their derivatives been precomputed?
     */
    //@{
    bool ev_precomp;
    bool d1_precomp;
    bool d2_precomp;
    //@}

    /// Precomputed values using precomp_rbfs
    mutable map<std::vector<size_t>, map<point, double> >
    precomp_values;
    /// Same, but storing the vectors.
    mutable map<std::vector<size_t>, linalg::vector >
    precomp_values_vec;
    /// Convert vector of values to a map
    map<point, double> valsvec2map(const linalg::vector& values) const;
  };

  /// Interpolated values manipulated by the interpolator class.
  class interp_values{
  public:
    /** @name Arithmetic with interpolated values
     *
     * These operators allow pointwise arithmetic with interpolated
     * values returned by interpolators. They must work with
     * interpolated values from compatible interpolators
     * (i.e. interpolators with the same RBFs).
     */
    //@{
    interp_values operator+(const interp_values& w) const;
    interp_values operator-(const interp_values& w) const;
    interp_values operator*(const interp_values& w) const;
    interp_values operator/(const interp_values& w) const;
    interp_values operator*(const double a) const;
    interp_values operator/(const double a) const;
    //@}
  private:
    
    /// No construction!
    interp_values();
    /// Interpolator creates interp_values with this
    interp_values(const linalg::vector& v_in, size_t rbfs_hash_in) 
      : v(v_in), rbfs_hash(rbfs_hash_in) {};

    ///Interpolated values go here
    linalg::vector v;
    ///Two interp_values can be added iff these two values matches.
    size_t rbfs_hash;
    ///Exception thrower
    void different_rbfs(int line, string file) const;

    template<typename RBF> friend class interpolator;
    friend interp_values operator*(const interp_values& v, const
				   interp_values& w); 
  };

  ///For comfortable syntax
  template <typename RBF>
  interpolator<RBF> operator*(double a, const interpolator<RBF>& u);
}

#endif