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

using namespace std;

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


//implements a multilinear interpolator on rectangular D-dimensional grids
//interpolation is not very efficient (linear in sum_i grid_size_in_dimension_i), 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 GridMultilinearInterpolator{
private:
	vector<vector<double> > grid_points; // grid_points[i][j] is the value of the j-th grid point along the i-th dimension
	vector<double> z_values; 			// grid in 'row-major' format
	
	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(const vector<double> &location, double default_z_value, bool flatExtrapolate) const;
	
	static long numberOfVertices(long D);
	
public:
	
	//setup D-dimensional grid from a string
	//string should be of format <x1 values>,<x2 values>,..,<x_D values>,<z values>
	//where each of the N+1 parts is a space-separated list of real numbers (of length N1, N2, .., ND and prod_i N_i, respectively)
	//z-values are in row-major format
	bool load(long D, const string &s);
	
	//append a (D-1)-dimensional layer to the existing grid, with the 1st coordinate being given
	//z_values must be of size prod_{i=2}^D N_i
	//Assumes that setLayerPoints(..) has been called already
	bool appendLayer(double location, const vector<double> &new_z_values);
	
	//define grid points except for 1st dimension
	//this will delete the entire pre-existing grid and start a new empty grid
	//required if you are going to build the grid layer-by-layer using appendLayer(..)
	void setLayerPoints(const vector<vector<double> > &_layer_points);
	void setLayerPoints(long D, ...); // missing D-1 arguments must be (const char*), each defining the grid points along one of the dimensions 2,..,D
	
	//generate a string representation of the grid, in the same format as used by load(..)
	string getGridAsString() const;
		
	//print grid in scan format
	void printGridScann(ostream &ss) const;

	//returns closest grid value if point is not within grid
	double interpolate(const vector<double> &location) const;
	double interpolate(int dummy, ...) const; // missing D arguments must be of type double
		
	//returns default_z_value if point is not within grid
	double interpolateWithDefault(const vector<double> &location, double default_z_value) const; // location[] must be of size D
	double interpolateWithDefault(double default_z_value, ...) const; // missing D arguments must be of type double
	
	//returns a one-line summary of the grid setup
	string getShortDescription() const;
	
	bool isWithinGrid(const vector<double> &location) const;
	
	long getGridSize() const;
	long get_count_along_dimension(long d) const{ return grid_points[d].size(); }
	double get_grid_point_along_dimension(long d, long i) const{ return grid_points[d][i]; }
	double get_grid_z_value(const vector<long> &point) const;


	// 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