#ifndef TIRA_TENSORFIELD
#define TIRA_TENSORFIELD

#include <valarray>
#include <string>
#include <fstream>
#include <initializer_list>

enum tira_floattype {unused0, unused1, float16, unused3, float32, unused4, unused5, unused6, unused7, float64};
enum tira_error {TIRA_SUCCESS, TIRA_UNDEFINED, TIRA_FILE_OPEN_ERROR, TIRA_UNKNOWN_FILE_TYPE};
enum tira_filetype {unused, spectrum, tensor2};
struct tira_header {
	std::valarray<size_t> shape;				//stores the resolution of the tensor field along each dimension
	tira_filetype filetype;						//stores the file type
	tira_floattype datatype;					//stores the data type of the tensor field
	float version;								//file version
};

tira_error tira_load_header(std::ifstream &infile, tira_header &header) {

	//see if the file represents a TIRA structure
	unsigned char TIRA[5];
	TIRA[4] = 0;
	infile.read((char*)TIRA, 4);
	if (strcmp((char*)TIRA, "TIRA") != 0)
		return TIRA_UNKNOWN_FILE_TYPE;

	//get the file type
	infile.read((char*)&(header.filetype), sizeof(tira_filetype));

	//get the version
	infile.read((char*)&(header.version), sizeof(float));

	//get the number of dimensions in the field
	size_t ndim;
	infile.read((char*)&ndim, sizeof(size_t));
	header.shape.resize(ndim);

	//read each dimension
	for (int d = 0; d < ndim; d++)
		infile.read((char*)&(header.shape[d]), sizeof(size_t));

	//read the data type
	infile.read((char*)&(header.datatype), sizeof(tira_floattype));

	return TIRA_SUCCESS;
}

template <typename T>
class tensorfield {
public:

	std::valarray<size_t> shape;				//number of samples in each dimension
	std::valarray<size_t> cum_offsets;			//cumulative factors for fast array accesses
	tira_filetype filetype;						//stores the file type
	tira_floattype datatype;					//stores the data type of the tensor field
	T* ptr;									//pointer to the raw tensor data

	void init() {
		size_t ndims = shape.size();			//get the number of dimensions
		cum_offsets.resize(ndims);				//allocate space for the cumulative offsets
		size_t cum = 1;
		for (int i = 0; i < ndims; i++) {
			cum_offsets[ndims - 1 - i] = cum;
			cum *= shape[ndims - 1 - i];
		}
	}

	size_t bytes() {
		size_t product = 1;
		for (size_t i = 0; i < shape.size(); i++)
			product *= shape[i];
		return product * (size_t)datatype;
	}

	// Function turns a series of array indices into a 1D index to the corresponding element
	inline size_t idx(std::initializer_list<size_t> coords) {
		size_t index = 0;
		size_t i = 0;
		for (auto c : coords) {
			index += c * cum_offsets[i];
			i++;
		}
		return index;
	}

	tira_error load_tira(std::string filename) {
		
		//open the file for binary reading
		std::ifstream infile(filename, std::ios::binary);
		if (!infile) return TIRA_FILE_OPEN_ERROR;

		tira_header header;
		tira_load_header(infile, header);
		shape = header.shape;
		filetype = header.filetype;
		datatype = header.datatype;
		
		//allocate space in memory for the tensor field
		size_t size_bytes = bytes();
		ptr = (T*)malloc(size_bytes);

		//copy data from the file to memory
		infile.read((char*)ptr, bytes());

		//close the file
		infile.close();

		//initialize all internal variables
		init();

		return TIRA_SUCCESS;
	}

	T& operator()(size_t x, size_t y, size_t z, size_t r, size_t c) {
		size_t i = idx({ z, y, x, r, c });
		return ptr[i];
	}

};

#endif