#include "include/interp_values.hpp"
//debug
#include <iostream>

namespace bvp{

  // ************************ interp_values stuff *****************
  badArgument interp_values::different_rbfs(int line, string file) const
  {
    badArgument exc;
    exc.reason = "Can't operate on interpolated values from "
      "incompatible interpolators";
    exc.line = line;
    exc.file = file;
    return exc;
  }

  interp_values 
  interp_values::operator+(const interp_values& w) const
  {
    if(rbfs_hash != w.rbfs_hash)
      throw different_rbfs(__LINE__, __FILE__);
    
    return interp_values(v + w.v,rbfs_hash,n,m);
  }

  interp_values 
  interp_values::operator-(const interp_values& w) const
  {
    if(rbfs_hash != w.rbfs_hash)
      throw different_rbfs(__LINE__, __FILE__);
    
    return interp_values(v - w.v,rbfs_hash,n,m);
  }

  interp_values 
  interp_values::operator*(const interp_values& w) const
  {
    if(rbfs_hash != w.rbfs_hash)
      throw different_rbfs(__LINE__, __FILE__);
    return interp_values(v * w.v,rbfs_hash,n,m);
  }

  interp_values 
  interp_values::operator/(const interp_values& w) const
  {
    if(rbfs_hash != w.rbfs_hash)
      throw different_rbfs(__LINE__, __FILE__);
    
    return interp_values(v / w.v,rbfs_hash,n,m);
  }

  interp_values
  interp_values::operator+(const double a) const
  {
    vector w(v.size(),a);
    return interp_values(v+w,rbfs_hash,n,m);
  }
  
  interp_values
  interp_values::operator-(const double a) const
  {
    vector w(v.size(),a);
    return interp_values(v-w, rbfs_hash,n,m);
  }

  interp_values
  interp_values::operator*(const double a) const
  {
    return interp_values(v*a,rbfs_hash,n,m);
  }
  
  interp_values
  interp_values::operator/(const double a) const
  {
    return interp_values(v/a, rbfs_hash,n,m);
  }

  //*************** bdry_vals stuff ************************************

  bdry_values::bdry_values(const interp_values& v_in) 
  : interp_values(v_in)
  {
    if(v.size() == m+n){       
      //Lower part of the vector contains the boundary
      linalg::slice s(n+1,v.size());
      v = v(s);
    }
    else if(v.size() != m){ 
      throw badArgument("Cannot initialise bdry_values from given "
			"interp_values", __FILE__,__LINE__);
    }
    // Else, we're initting with interp_values that are the right size
    // and there's no need to do anything, since they're bona fide
    // boundary values (parent copy ctor does the work).

      //debug
      using namespace std; 
      cout << "The boundary vector: " << endl
	   << v << endl;
  }

  //*************** normals stuff *************************************
  
  normals::normals(size_t rbfs_hash_in, 
	       const map<point, linalg::vector>& normals_in, 
	       size_t n_in) : rbfs_hash(rbfs_hash_in), n(n_in)
  {
    map<point, linalg::vector>::const_iterator I;
    size_t rows = normals_in.size();
    size_t cols = normals_in.begin() -> second.size();
    slice s(1,cols);
    size_t i = 1;
    matrix N_temp(rows,cols);
    for(I = normals_in.begin(); I != normals_in.end(); I++){
      N_temp(i,s) = (I->second)/(norm(I->second));
      i++;
    }
    N = N_temp;
  }

  bdry_values normals::operator()(size_t k) const
  {
    slice s(1,N.rows());
    return bdry_values(N(s,k),rbfs_hash,n,m);
  }

  bool normals::operator==(const normals& M) const
  {
    return rbfs_hash == M.rbfs_hash;
  }

  bool normals::operator!=(const normals& M) const
  {
    return rbfs_hash != M.rbfs_hash;
  }

  //*************** Non-member global friend functions ******************

  interp_values operator+(double a, const interp_values&  v)
  { 
    return v+a;
  } 

  interp_values operator-(double a, const interp_values&  v)
  { 
    linalg::vector w(v.v.size(),a);
    return interp_values(w-v.v, v.rbfs_hash,v.n,v.m);
  } 

  interp_values operator*(double a, const interp_values&  v)
  { 
    return v*a;
  } 

  interp_values operator/(double a, const interp_values&  v)
  { 
    linalg::vector w(v.v.size(),a);
    return interp_values(w/v.v, v.rbfs_hash,v.n,v.m);
  } 
}