#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;

  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);
    //@}
    
    /** @name Evaluations and derivatives
     */
    //@{
    ///Evaluation
    double operator()(const point& p) const;
    ///Evaluation
    double at(const point& p) const;

    /// First derivative
    double d(const point& p, size_t k) const;
    /// Second derivatives
    double d2(const point &p, 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_at();
    ///Precompute all relevant first derivatives.
    void precompute_d() ;
    ///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;
    void not_initted(int line, string file) const; //Exception thrower
    
    ///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 vector 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;
    
    /// 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;
    /// Invert phis\coeffs and return result vector in a map.
    map<point, double> precompute_values(const matrix& phis) const;
  };

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

#endif