#ifndef GRID_3D_INTERPOLATOR_DEF
#define GRID_3D_INTERPOLATOR_DEF
#include "Grid3DInterpolator.h"


void Grid3DInterpolator::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 Grid3DInterpolator::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 Grid3DInterpolator::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 Grid3DInterpolator::load(const string &s){
	vector<string> parts;
	explodeString(parts, s, ",", 4);
	if(parts.size()<4) return false;
	string2numvector(parts[0], s_values);
	string2numvector(parts[1], r_values);
	string2numvector(parts[2], c_values);
	string2numvector(parts[3], z_values);
	if(z_values.size() != s_values.size()*r_values.size()*c_values.size()) return false;
	return true;
}




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


bool Grid3DInterpolator::isWithinGrid(double s_value, double r_value, double c_value) const{
	if((s_values.size()<2) || (r_values.size()<2) || (c_values.size()<2)) return false;
	if(s_value+abs(s_value)*EPSILON<s_values.front() || s_value-abs(s_value)*EPSILON>s_values.back()) 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;
}



long Grid3DInterpolator::findOn1DGrid(const vector<double> &gridValues, double value){
	if(value<gridValues.front()) return -1;
	if(value>gridValues.back()) return gridValues.size();
	const double value_du = value + abs(value)*EPSILON;
	const double value_dd = value - abs(value)*EPSILON;
	for(long i=1; i<gridValues.size(); ++i){
		if(gridValues[i-1]<=value_du && gridValues[i]>=value_dd){ return i-1; }
	}
	return 0; // this should not happen
}


double Grid3DInterpolator::interpolate(	double s_value, 
										double r_value, 
										double c_value, 
										double default_z_value, 
										bool flatExtrapolate) const{
	//find grid box
	long s,r,c;
	double sv1, sv2, rv1, rv2, cv1, cv2, tz1, tz2, z;
	const long NS = s_values.size();
	const long NR = r_values.size();
	const long NC = c_values.size();
	double V[2][2][2];
	
	//find appropriate grid box (linear search)
	s = findOn1DGrid(s_values, s_value);
	if(s<0 || s>=s_values.size()){
		// not within s-grid
		if(!flatExtrapolate) return default_z_value;
	}else{
		// within s-grid
		sv1 = s_values[s];
		sv2 = s_values[s+1];
	}
	r = findOn1DGrid(r_values, r_value);
	if(r<0 || r>=NR){
		// not within s-grid
		if(!flatExtrapolate) return default_z_value;
	}else{
		// within s-grid
		rv1 = r_values[r];
		rv2 = r_values[r+1];
	}
	c = findOn1DGrid(c_values, c_value);
	if(c<0 || c>=NC){
		// not within s-grid
		if(!flatExtrapolate) return default_z_value;
	}else{
		// within s-grid
		cv1 = c_values[c];
		cv2 = c_values[c+1];
	}
	int i,j,k;
	for(i=0; i<=1; ++i){
		for(j=0; j<=1; ++j){
			for(k=0; k<=1; ++k){
				V[i][j][k] = z_values[(s<0 ? 0 : (s>=NS ? NS-1 : (i==0 ? s : s+1)))*NR*NC + (r<0 ? 0 : (r>=NR ? NR-1 : (j==0 ? r : r+1)))*NC + (c<0 ? 0 : (c>=NC ? NC-1 : (k==0 ? c : c+1)))];
			}
		}
	}
		
	//trilinear interpolation:  http://en.wikipedia.org/wiki/Trilinear_interpolation#Method
	double ds = (s<0 || s>=NS ? 0 : (s_value - sv1)/(sv2-sv1));
	double dr = (r<0 || r>=NR ? 0 : (r_value - rv1)/(rv2-rv1));
	double dc = (c<0 || c>=NC ? 0 : (c_value - cv1)/(cv2-cv1));
	double tzv[2][2];
	for(r=0; r<=1; ++r){
		for(c=0; c<=1; ++c){
			tzv[r][c] = V[0][r][c] * (1-ds) + V[1][r][c] * ds;
		}
	}
	return (tzv[0][0]*(1-dr)+tzv[1][0]*dr) * (1-dc)  +  (tzv[0][1]*(1-dr)+tzv[1][1]*dr)*dc;
}



double Grid3DInterpolator::interpolate(double s_value, double r_value, double c_value, double default_z_value) const{
	return interpolate(	s_value, 
						r_value, 
						c_value, 
						default_z_value, 
						false);
}


double Grid3DInterpolator::interpolate(double s_value, double r_value, double c_value) const{
	return interpolate(	s_value, 
						r_value, 
						c_value, 
						0, 
						true);
}



void Grid3DInterpolator::setRCvalues(const vector<double> &_r_values, const vector<double> &_c_values){
	r_values = _r_values;
	c_values = _c_values;
	s_values.clear();
	z_values.clear();
}



void Grid3DInterpolator::setRCvalues(const string &_r_values, const string &_c_values){
	string2numvector(_r_values, r_values);
	string2numvector(_c_values, c_values);
	s_values.clear();
	z_values.clear();
}



bool Grid3DInterpolator::appendSheet(double s_value, const vector<double> &new_z_values){
	if(new_z_values.size() != r_values.size()*c_values.size()) return false;
	s_values.push_back(s_value);
	z_values.insert(z_values.end(), new_z_values.begin(), new_z_values.end());
	return true;
}


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


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



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

#endif