#ifndef GRID_3D_INTERPOLATOR_DECL
#define GRID_3D_INTERPOLATOR_DECL
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#include <limits>
#include <sstream>

using namespace std;

#ifndef EPSILON
#define EPSILON 1e-15
#endif


//implements a bilinear interpolator on rectangular 2D grids
//interpolation is not very efficient (linear in row_counts+col_counts+sheet_counts), so this is meant for cases where only a few interpolations are needed
//extrapolation of outliers is flat, i.e. by choosing the closest grid value
//alternatively, a default value can be specified at call-time for outliers
class Grid3DInterpolator{
private:
	vector<double> s_values, r_values, c_values; // grid coordinates
	vector<double> z_values; //qube in sheet-row-major format of size NS x NR x NC
	
	static void string2numvector(const string &s, vector<double> &v);
	static string vector2string(const vector<double> &v, string separator);
		
	void explodeString(vector<string> &parts, const string &haystack, const string &separator, int maxPartsCount);

	double interpolate(double s_value, double r_value, double c_value, double default_z_value, bool flatExtrapolate) const;
	
public:
	
	//setup grid from a string
	//string should be of format <s values>,<r values>,<c values>,<z values>
	//where each of the 4 parts is a space-separated list of real numbers (of length NS, NR, NC and NS*NR*NC, respectively)
	//z-values are in sheet-row-major format
	bool load(const string &s);
	
	//append a sheet to the existing grid
	//z_values must be of size get_r_count()*get_c_count()
	//Assumes that setRCvalues(..) has been called already
	bool appendSheet(double s_value, const vector<double> &new_z_values);
	
	//define sheet r_values x c_values
	//this will delete the entire pre-existing grid and start a new empty grid
	//required if you are going to build the grid sheet-by-sheet using appendSheet(..)
	void setRCvalues(const vector<double> &_r_values, const vector<double> &_c_values);
	void setRCvalues(const string &_r_values, const string &_c_values);
	
	//generate a string representation of the grid, in the same format as used by load(..)
	string getGridAsString() const;
	
	//print grid z_values in matrix format
	void printGridMatrix(ostream &ss) const;
	
	//print grid s,r,c and z_values in scan format
	void printGridScann(ostream &ss) const;

	//returns closest grid value if point is not within grid
	double interpolate(double s_value, double r_value, double c_value) const;
		
	//returns default_z_value if point is not within grid
	double interpolate(double s_value, double r_value, double c_value, double default_z_value) const;
	
	//returns a one-line summary of the grid setup
	string getShortDescription() const;
	
	bool isWithinGrid(double s_value, double r_value, double c_value) const;
	
	long get_s_count() const{ return s_values.size(); }
	long get_r_count() const{ return r_values.size(); }
	long get_c_count() const{ return c_values.size(); }
	double get_s_value(long s) const{ return s_values[s]; }
	double get_r_value(long r) const{ return r_values[r]; }
	double get_c_value(long c) const{ return c_values[c]; }
	double get_z_value(long s, long r, long c) const{ return z_values[s*r_values.size()*c_values.size() + r*c_values.size()+c]; }


	// returns left-index of grid box (0,..,lastIndex-1) or -1 if left-outlier or lastIndex+1 if right-outlier
	static long findOn1DGrid(const vector<double> &gridValues, double value);

};



#endif