// 
//  NNeighbour.h
//  Qdos2 - Nearest Neighbour Algorithm Template Functions
//
//  Janet Barnes
//  Sep 2000
//

#ifndef QDOS_NNEIGHBOUR_H
#define QDOS_NNEIGHBOUR_H

#include "Algorithm.h"

namespace QDOS {

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ nneighbour
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------- nneighbour
template <class Tx, class Ty>
bool nneighbour( const NumericSequence<Tx>& x_data, \
    const NumericSequence<Ty>& y_data, const Tx& x_interp, Ty& y_interp)
{
// calculate the value of a function y=f(x) (as given by values of y, y_data,
// at values of x, x_data) at x_interp using nearest neighbour method
// floor_idx is the index of the greatest value of x such that 
// x[floor_idx] <= x_interp
  
  int floor_idx = floor_index( x_data, x_interp );
  
  const int max_idx = x_data.size()-1;
  if (floor_idx==-1)
      y_interp=y_data[0];
  else if (floor_idx==max_idx)
      y_interp=y_data[max_idx];
  else
  {
      double factor=( x_interp - x_data[floor_idx]) /  \
          (x_data[floor_idx+1]-x_data[floor_idx]);
      y_interp=(factor<0.5)?y_data[floor_idx]:y_data[floor_idx+1];
  }
  return true;
}

template <class Tx, class Ty>
bool nneighbour( const NumericSequence<Tx>& x_data, \
    const NumericSequence<Ty>& y_data, const Tx& x_interp, Ty& y_interp,\
    const Tx& gap_criterion)
{
// calculate the value of a function y=f(x) (as given by values of y, y_data,
// at values of x, x_data) at x_interp using nearest neighbour method
// floor_idx is the index of the greatest value of x such that 
// x[floor_idx] <= x_interp
  
  int floor_idx = floor_index( x_data, x_interp );
  if (!interp_at_pt(x_data, x_interp, floor_idx, gap_criterion))
      return false; 
  
  const int max_idx = x_data.size()-1;
  if (floor_idx==-1)
      y_interp=y_data[0];
  else if (floor_idx==max_idx)
      y_interp=y_data[max_idx];
  else
  {
      double factor=( x_interp - x_data[floor_idx]) /  \
          (x_data[floor_idx+1]-x_data[floor_idx]);
      y_interp=(factor<0.5)?y_data[floor_idx]:y_data[floor_idx+1];
  }
  return true;
}
// ++++++++++++++++++++++++++++++++++++++++++++++++++++ series_nneighbour
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ---------------------------------------------------- series_nneighbour
template<class Tx, class Ty>
vector<bool> series_nneighbour( const NumericSequence<Tx>& x_data, \
    const NumericSequence<Ty>& y_data, const NumericSequence<Tx>& x_interp, \
    NumericSequence<Ty>& y_interp )
{
// find value of x_data "nearest" x_interp and use its correspoinding value
// of y_data to estimate value of y=f(x), at x=x_interp
    data_is_ok(x_data,y_data);
    int sz=x_interp.size();
    vector<bool> mask(sz,false);
    y_interp.resize(sz);
    for(register int i=0; i< sz; i++)
        mask[i]=nneighbour( x_data, y_data, x_interp[i],y_interp[i] );
    return mask;
}

}  // end namespace QDOS
#endif
