/*!\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
 */
#pragma once

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

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

/** \defgroup Interpolator */
/// \{
class vtkplot;

//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< kwantix::point, double>). Once an interpolator is
 * initialised, it becomes a bona fide kwantix::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 intr_values &f);
  void set_g(const bdry_values &g);
  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;
    
  /** \brief This will put into an ostream a kwantix::matrix where
   *   the first \f$n-1\f$ columns will be all the points where the
   *   interpolator is defined and the \f$n^{th}\f$ column is the
   *   values that the interpolator takes on that point.
   */
  void into_os(std::ostream& os) const;

  ///Needed in order to plot the data
  friend class vtkplot;
    
private:
  ///Once the matrix is defined, this function inverts it.
  void computecoeffs(bool rhs_defined = false); 
	 
  ///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.
  shared_ptr<matrix> M; 

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


  /** @name Exception builders
   */
  //@{
  kwantix::badArgument not_initted(int line, string file) const; 
  kwantix::badArgument not_precomputed(int line, string file) const;
  //@}
    
  ///Coefficients of the RBFs
  kwantix::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>, 
              shared_ptr<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>, kwantix::vector >
  precomp_values_vec;
  /// Convert vector of values to a map
  map<point, double> valsvec2map(const  kwantix::vector& values) const;
};
/// \}


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

/// Stream insertion
template <typename RBF>
std::ostream& operator<< (std::ostream& os, const interpolator<RBF>& u)
{
  u.into_os(os);
  return os;
}

}//namespace kwantix