// 
//  Algorithm.h
//  Qdos2 - Common Algorithm Routines
//
//  Janet Barnes
//  Nov 1999
//

#ifndef QDOS_ALGORITHM_H
#define QDOS_ALGORITHM_H

#include "Exception.h"
#include "Sequence.h"
#include <algorithm>
#include <iterator>

namespace QDOS {

class AlgorithmException: public Exception
{
public:
    AlgorithmException(ExceptionMsg m=""):\
        Exception(m<<"[AlgorithmException]"){;}
};

class BadData: public AlgorithmException
{
public:
    BadData(ExceptionMsg m=""):\
        AlgorithmException(m<<"[BadData]"){;}
};

class DataNotConformant: public AlgorithmException
{
public:
    DataNotConformant(ExceptionMsg m=""):\
        AlgorithmException(m<<"[DataNotConformant]"){;}
};

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ data_is_ok
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------- data_is_ok
template <class Tx, class Ty>
void data_is_ok( const NumericSequence<Tx>& x_data, \
    const NumericSequence<Ty>& y_data )
{
  if( y_data.size() != x_data.size() )
    throw DataNotConformant("ERROR: x and y sequences are different sizes!");
  if( y_data.size() == 0 ) 
    throw DataNotConformant("ERROR: y sequence has 0 elements!");
  if(!x_data.is_strictly_monotonic_increasing())
  {	  
    throw DataNotConformant("ERROR: x sequence is not strictly mononotic increaing (x[i] is NOT > x[i-1])!");
  }
}


// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ floor_index
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------- floor_index
template <class Tx>
int floor_index( const NumericSequence<Tx>& x_data, \
    const Tx& x_interp )
{
// Returns the index in x_data of the last element >=x_interp
// If x_interp <= x_data[0] returns -1
// If x_interp >= x_data[x.data.size()-1] returns x.data.size()-1
// upper_bound returns itr such that x[itr] > x_interp 
    typename vector<Tx>::const_iterator itr=\
	upper_bound(x_data.begin(),x_data.end(),\
        x_interp);         
    int dist=-1;
    dist+=distance(x_data.begin(), itr);
    return dist;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ interp_at_pt
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- interp_at_pt
template<class Tx>
bool interp_at_pt(const NumericSequence<Tx>& x_data, const Tx& x_interp,\
    int  floor_idx,const Tx& gap_criterion)
{
    const int max_idx=x_data.size()-1;       
    if(floor_idx==-1)
        return x_data[0]-x_interp<gap_criterion;
    if(floor_idx>=max_idx)
        return x_interp-x_data[max_idx]<gap_criterion;
    return x_data[floor_idx+1]-x_data[floor_idx] < gap_criterion;     
}              


}  // end namespace QDOS
#endif
