#ifndef __INTERP_VALUES_HPP__
#define __INTERP_VALUES_HPP__

#include "linalg.hpp"
#include <map>

namespace bvp{
  using std::map;
  using namespace linalg;

  class bdry_values;

  //FIXME: Give this a better name
  /// A unit normals class associated to an interpolator
  class normals{
  public:
    /** \brief Return the kth coordinate of all normals along the
     * boundary as interp values
     */
    bdry_values operator()(size_t k) const;
    
    /*! \name comparisons 
     *
     * Normals are the same if they correspond to the same domain
     *(i.e. if their hash coincides)
     */
    //@{
    bool operator!=(const normals& M) const;
    bool operator==(const normals& M) const;
    //@}
  private:
    template<typename RBF> friend class interpolator;
    ///Only constructor available for interpolator class
    normals(size_t rbfs_hash_in, const map<point, linalg::vector>&
            normals_in, size_t n_in);
    ///Does nothing
    normals(){};
    ///Exception builder
    badArgument not_initted() const;
    
    ///A hash used to return bdry_values objects
    size_t rbfs_hash;

    ///Number of interior points in the matrix, which interp_values need
    size_t n;
    ///Number of boundary points
    size_t m;

    ///Store normals in a matrix
    linalg::matrix N;
  };




  /// 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;
    interp_values operator*(const double a) const;
    interp_values operator/(const double a) const;
    //@}

    /*! \name Non-member arithmetic operators
     *
     * For comfortable syntax.
     */
    //@{
    friend interp_values operator+(double a, const interp_values& v); 
    friend interp_values operator-(double a, const interp_values& v); 
    friend interp_values operator*(double a, const interp_values& v); 
    friend interp_values operator/(double a, const interp_values& v); 
    //@}

  protected:
    /// Interpolator creates interp_values with this
    interp_values(const linalg::vector& v_in, size_t rbfs_hash_in,
                  size_t n_in, size_t m_in) 
      : v(v_in), rbfs_hash(rbfs_hash_in), n(n_in), m(m_in) {};

    /// No construction!
    interp_values();

  protected:
    ///Exception builder
    badArgument different_rbfs(int line, string file) const;

    ///Interpolated values go here
    linalg::vector v;
    ///Two interp_values can be added iff these two values matches.
    size_t rbfs_hash;
    ///Number of interior points
    size_t n;
    ///Number of boundary points
    size_t m;

    template<typename RBF> friend class interpolator;
  };
  


  /// For referring to values purely on the boundary of an interpolator.
  class bdry_values : public interp_values{
  public:
    ///Empty initialisation
    bdry_values(){};
    ///Initialise the boundary values with interp_values
    bdry_values(const interp_values& v_in);
  private:
    friend class normals;
    ///Same as for parent class
    bdry_values(const linalg::vector& v_in, size_t rbfs_hash_in,
                size_t n_in, size_t m_in) 
      : interp_values(v_in, rbfs_hash_in, n_in, m_in) {};
  };

  /// For referring to values purely on the interior of an interpolator.
  class intr_values : public interp_values{
  public:
    ///Empty initialisation
    intr_values(){};
    ///Initialise the interior values with interp_values
    intr_values(const interp_values& v_in);
  private:
    ///Same as for parent class
    intr_values(const linalg::vector& v_in, size_t rbfs_hash_in,
                size_t n_in, size_t m_in) 
      : interp_values(v_in, rbfs_hash_in, n_in, m_in) {};
  };





  /*! \name Arithmetic operators
   * For comfortable syntax.
   */
  //@{
  interp_values operator+(double a, const interp_values& v); 
  interp_values operator-(double a, const interp_values& v); 
  interp_values operator*(double a, const interp_values& v); 
  interp_values operator/(double a, const interp_values& v); 
  //@}
}

#endif //__INTERP_VALUES_HPP__