math/function.h

#ifndef RTS_FUNCTION_H
#define RTS_FUNCTION_H

#include <string>

namespace rts{

//template class for a one-dimensional function
template <class Tx, class Ty>
class function
{
	std::vector<Tx> X;
	std::vector<Ty> Y;

	//comparison function for searching lambda
	static bool findCeiling(Tx ax, Tx bx){
	    return (ax > bx);
	}

	//process a string to extract a function (generally delimited by tabs, commas, spaces, etc.)
	void process_string(std::string s){
		std::stringstream ss(s);

		//std::string test;
		//std::getline(ss, test);
		//std::cout<<test;
		//exit(1);

		Tx x;
		Ty y;
		std::string line;

		while(!ss.eof()){
			
			std::getline(ss, line);
			if(line[0] == '#') continue;

			std::stringstream lstream(line);

			lstream>>x;			//read the x value
			lstream>>y;			//read the y value

			//std::cout<<x<<", "<<y<<std::endl;

			if(ss.eof()) break;
			insert(x, y);			//insert the read value into the function
			
		}	
	}


public:

	//linear interpolation
	Ty linear(Tx x) const
	{
		if(X.size() == 0)	return (Ty)0;	//return zero if the function is empty

		unsigned int N = X.size();			//number of sample points

		//declare an iterator
		typedef typename std::vector< Tx >::iterator f_iter;	//declare an iterator
		f_iter it;

		//find the first X-coordinate that is greater than x
		unsigned int i;
		for(i = 0; i<N; i++){
			if(X[i] > x)
				break;
		}
		//i currently holds the ceiling

		//if the wavelength is past the end of the list, return the last sample point
		if(i == N) return Y[N];
		//if the wavelength is before the beginning of the list, return the front
		else if(i == 0) return Y[0];
		//otherwise interpolate
		else{
		    Tx xMax = X[i];
		    Tx xMin = X[i-1];

		    Tx a = (x - xMin) / (xMax - xMin);
		    Ty riMin = Y[i - 1];
		    Ty riMax = Y[i];
		    Ty interp = riMax * a + riMin * (1 - a);
		    return interp;
		}
	}

	///add a data point to a function
	void insert(Tx x, Ty y)
	{
		unsigned int N = X.size();		//number of sample points

		if(N == 0 || X[N-1] < x){
			X.push_back(x);
			Y.push_back(y);
			return;
		}

		//declare an iterator and search for the x value
        typename std::vector< Tx >::iterator it;		
        it = search(X.begin(), X.end(), &x, &x + 1, &function<Tx, Ty>::findCeiling);

        //if the function value is past the end of the vector, add it to the back
        if(*it == N){
        	X.push_back(x);
        	Y.push_back(y);
        }
        //otherwise add the value at the iterator position
        else{
			X.insert(it, x);
			Y.insert(Y.begin() + *it, y);
		}

	}

	Tx getX(unsigned int i) const{
		return X[i];
	}

	Ty getY(unsigned int i) const{
		return Y[i];
	}

	///get the number of data points in the function
	unsigned int getN() const{
		return X.size();
	}

	//look up an indexed component
	Ty operator[](int i) const{
		if(i <= X.size()){
			std::cout<<"ERROR: accessing non-existing sample point in 'function'"<<std::endl;
			exit(1);
		}
		return Y[i];
	}

	///linear interpolation
	Tx operator()(Tx x) const{
		return linear(x);
	}

	//add a constant to the function
	function<Tx, Ty> operator+(Ty r) const{

		function<Tx, Ty> result;

		//if there are points in the function
		if(X.size() > 0){
			//add r to every point in f
			for(unsigned int i=0; i<X.size(); i++){
				result.X.push_back(X[i]);
				result.Y.push_back(Y[i] + r);
			}
		}
		else{
			result = r;
		}

		return result;
	}

	function<Tx, Ty> & operator= (const Ty & rhs){		
		X.clear();
		Y.clear();
		if(rhs != 0)			//if the RHS is zero, just clear, otherwise add one value of RHS
			insert(0, rhs);

		return *this;
	}


	std::string str(){
		stringstream ss;

		unsigned int N = X.size();		//number of sample points
		//output each function value
		for(unsigned int i = 0; i<N; i++){
			ss<<X[i]<<", "<<Y[i]<<std::endl;
		}

		return ss.str();

	}

	void load(std::string filename){
		std::ifstream t(filename.c_str());
		std::string str((std::istreambuf_iterator<char>(t)),
		std::istreambuf_iterator<char>());

		process_string(str);
	}


};

}	//end namespace rts


#endif