#ifndef GRID_2D_INTERPOLATOR_DEF
#define GRID_2D_INTERPOLATOR_DEF
#include "Grid2DInterpolator.h"


void Grid2DInterpolator::string2numvector(const string &s, vector<double> &v){
	istringstream ss(s);
	v.clear();
	double d;
	do{
		ss >> d;
		if(ss.fail()) break;
		v.push_back(d);
	}while(true);
}


string Grid2DInterpolator::vector2string(const vector<double> &v, string separator){
	ostringstream ss;
	for(long n=0; n<v.size(); ++n){
		ss << (n==0 ? "" : separator) << v[n];
	}
	return ss.str();
}


void Grid2DInterpolator::explodeString(vector<string> &parts, const string &haystack, const string &separator, int maxPartsCount){
	parts.clear();
	int previousPos = 0, nextPos, separatorLen = separator.length();
	while(((nextPos = haystack.find(separator, previousPos)) != string::npos) && ((maxPartsCount>1) || (maxPartsCount<0))){
		parts.push_back(haystack.substr(previousPos, nextPos - previousPos));
		previousPos = nextPos + separatorLen;
		-- maxPartsCount;
	}
	parts.push_back(haystack.substr(previousPos));
}


bool Grid2DInterpolator::load(const string &s){
	vector<string> parts;
	explodeString(parts, s, ",", 3);
	if(parts.size()<3) return false;
	string2numvector(parts[0], r_values);
	string2numvector(parts[1], c_values);
	string2numvector(parts[2], z_values);
	if(z_values.size() != r_values.size()*c_values.size()) return false;
	return true;
}




string Grid2DInterpolator::getGridAsString() const{
	return vector2string(r_values," ") + "," + vector2string(c_values," ") + "," + vector2string(z_values," ");
}


bool Grid2DInterpolator::isWithinGrid(double r_value, double c_value) const{
	if((r_values.size()<2) || (c_values.size()<2)) return false;
	if(r_value+abs(r_value)*EPSILON<r_values.front() || r_value-abs(r_value)*EPSILON>r_values.back()) return false;
	if(c_value+abs(c_value)*EPSILON<c_values.front() || c_value-abs(c_value)*EPSILON>c_values.back()) return false;
	return true;
}

double Grid2DInterpolator::interpolate(double r_value, double c_value, double default_z_value) const{
	//find grid box
	long r1,r2,c1,c2;
	const long NC = c_values.size();
	double rv1, rv2, cv1, cv2, z11, z12, z21, z22, tz1, tz2, z;
	if(!isWithinGrid(r_value, c_value)) return default_z_value;
	const double r_value_du = r_value + abs(r_value)*EPSILON;
	const double r_value_dd = r_value - abs(r_value)*EPSILON;
	const double c_value_du = c_value + abs(c_value)*EPSILON;
	const double c_value_dd = c_value - abs(c_value)*EPSILON;
	
	//find appropriate grid box
	for(long r=1; r<r_values.size(); ++r){
		if(r_values[r-1]<=r_value_du && r_values[r]>=r_value_dd){ r1=r-1; r2=r; break; }
	}
	for(long c=1; c<c_values.size(); ++c){
		if(c_values[c-1]<=c_value_du && c_values[c]>=c_value_dd){ c1=c-1; c2=c; break; }
	}
	rv1 = r_values[r1];
	rv2 = r_values[r2];
	cv1 = c_values[c1];
	cv2 = c_values[c2];
	z11 = z_values[r1*NC + c1];
	z12 = z_values[r1*NC + c2];
	z21 = z_values[r2*NC + c1];
	z22 = z_values[r2*NC + c2];
		
	//bilinear interpolation: http://en.wikipedia.org/wiki/Bilinear_interpolation
	tz1 = z11*(rv2-r_value)/(rv2-rv1) + z21*(r_value-rv1)/(rv2-rv1);
	tz2 = z12*(rv2-r_value)/(rv2-rv1) + z22*(r_value-rv1)/(rv2-rv1);
	z 	= tz1*(cv2-c_value)/(cv2-cv1) + tz2*(c_value-cv1)/(cv2-cv1);
	return z;
}


double Grid2DInterpolator::interpolate(double r_value, double c_value) const{
	return interpolate(r_value, c_value, numeric_limits<double>::quiet_NaN());
}


void Grid2DInterpolator::setCvalues(const vector<double> &_c_values){
	c_values = _c_values;
	r_values.clear();
	z_values.clear();
}


bool Grid2DInterpolator::appendRow(double r_value, const vector<double> &new_z_values){
	if(new_z_values.size() != c_values.size()) return false;
	r_values.push_back(r_value);
	z_values.insert(z_values.end(), new_z_values.begin(), new_z_values.end());
	return true;
}


void Grid2DInterpolator::printGridMatrix(ostream &ss) const{
	long r,c;
	long NC = c_values.size();
	for(r=0; r<r_values.size(); ++r){
		for(c=0; c<NC; ++c){
			ss << (c==0 ? "" : " ") << z_values[r*NC + c];
		}
		ss << "\n";
	}
}


void Grid2DInterpolator::printGridScann(ostream &ss) const{
	long r,c;
	long NC = c_values.size();
	for(r=0; r<r_values.size(); ++r){
		for(c=0; c<NC; ++c){
			ss << r_values[r] << " " << c_values[c] << " " << z_values[r*NC + c] << "\n";
		}
		ss << "\n";
	}
}



//returns a one-line summary of the grid setup
string Grid2DInterpolator::getShortDescription() const{
	ostringstream ss;
	ss << "2D grid interpolator: " << r_values.size() << " x " << c_values.size() << " regular grid";
	return ss.str();
}

#endif