/*!\file interpolator.hpp
 *
 * \brief An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in \Xi}
 * \lambda_\xi \varphi_\xi(x) \f$, templated by RBFs, and its helper
 * classes are declared here.
 * 
 *  \file interpolator.cpp
 *  \brief Implementations and instantiations of the functions and
 *  classes declared in interpolator.hpp
 */
#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"
#include "interp_values.hpp"

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

  //FIXME: Break this up into an ansatz class and an interpolator
  //class. 
  /*! \brief An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in
   *  \Xi} \lambda_\xi \varphi_\xi(x) \f$, where the \f$ \varphi_\xi
   *  \f$ are RBFs.
   *
   * An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in \Xi}
   * \lambda_\xi \varphi_\xi(x) \f$, where the \f$\varphi_\xi\f$ are
   * RBFs centred at \f$\xi\f$, and \f$\Xi\f$ is a set of points on
   * the interior and boundary of some domain \f$\Omega\f$.
   *
   * The interpolator can be used in various ways, initialised by
   * either a BVP::linear_bvp2 or by interpolation data
   * (i.e. std::map< linalg::point, double>). Once an interpolator is
   * initialised, it becomes a bona fide bvp::realfunc, so that its
   * evaluations and derivatives become available. For certain
   * problems, it is convenient to be able to evaluate the
   * interpolator or its derivatives at all points of the domain. It
   * is also useful to be able to interpolate again on the same domain
   * or with a slightly modified BVP. Interfaces to both of these
   * functions are provided.
   */
  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. The reason why the domain
     *  might be relevant is in case it is important to specify which
     *  points are on the interior and which are on the boundary, so
     *  that it's possible to interpolate again by changing the values
     *  on either the interior or the boundary.
     */
    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 derivative
    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 Reinterpolate again 
     *
     * As above group of functions, but for purely interpolating
     * interpolators 
     */
    //@{
    void set_bdry_values(const bdry_values& b_new);
    //@}

    /** @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;
    //@}

    /** \brief Provide a unit normals object associated to the
     *  interpolator.
     */
    normals get_normals() const;
    
  private:
    ///Once the matrix is defined, this function inverts it.
    void computecoeffs(); 
	 
    ///Perform the actual interpolation.
    void init(shared_ptr<linear_BVP2> bvp);

    ///If the coefficients change, need to recompute precomputed values 
    void recompute_values_vec();
    
    ///BVP associated to this interpolator
    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 builders
     */
    //@{
    error_handling::badArgument not_initted(int line, string file) const; 
    error_handling::badArgument 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);

    ///The normals at the domain's boundary
    normals nrmls;
       
    /** \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;
  };


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

#endif