#ifndef _H_MATRIX4D
#define _H_MATRIX4D

#include <vector>
using namespace std;

typedef unsigned int uint;
class M4DRange
{
public:
	uint LowBound, HighBound;
	M4DRange(): Placeholder_All(true) {};
	M4DRange(uint low): LowBound(low), HighBound(low+1), Placeholder_All(false){};
	M4DRange(uint low, uint high): LowBound(low), HighBound(high), Placeholder_All(false){};
	void SetNMax(uint NumOfElements) 
	{		
		if (Placeholder_All) {
			LowBound=0; 
			HighBound=NumOfElements;
		};
		if (HighBound>NumOfElements) {
			HighBound=NumOfElements;
		}
		if (LowBound>HighBound) {
			LowBound=HighBound;
		}
	};
	uint GetDiff()
	{	
		return HighBound-LowBound;
	}
private:
	bool Placeholder_All;
};

class Matrix4D
{
private:
	uint DimA, DimB, DimC, DimD;
	vector <vector <vector <vector <float> > > > Data;
public:
	Matrix4D(){};
	Matrix4D(uint a, uint b, uint c, uint d): 
		DimA(a), DimB(b), DimC(c), DimD(d) 
	{
		Resize();
	};
	void GetDimensions(uint &a, uint &b, uint &c, uint &d) const
	{
		a=DimA;
		b=DimB;
		c=DimC;
		d=DimD;
	}
	
	void Resize(uint a, uint b, uint c, uint d)
	{
		DimA=a;
		DimB=b;
		DimC=c;
		DimD=d;
		Resize();
	}
	
	void Resize()
	{
		Data.resize(DimA);
		for(uint ia=0;ia<DimA;++ia) {
			Data[ia].resize(DimB);
			for (uint ib=0;ib<DimB;++ib) {
				Data[ia][ib].resize(DimC);
				for (uint ic=0;ic<DimC;++ic) {
					Data[ia][ib][ic].resize(DimD);
				}
			}
		}
	}
	
	
	void SetValue(float val)
	{
		for(uint ia=0;ia<DimA;++ia) {
			for (uint ib=0;ib<DimB;++ib) {
				for (uint ic=0;ic<DimC;++ic) {
					for (uint id=0;id<DimD;++id) {
						Data[ia][ib][ic][id]=val;
					}
				}
			}
		}
	}
	
	void SetValue(float val, uint a, uint b, uint c, uint d)
	{
		// if (a<DimA && b<DimB && c<DimC && d<DimD) 
		{
			Data[a][b][c][d]=val;
		}
	}
	
	/** Set Data in Matrix
	 * 
	 * @param InputData is a input vector. \n
	 * 				pay attention to data order!\n
	 * 				Data[a][b][c][d] == InputData[x] 
	 * 				Data[a][b][c][d+k] == InputData[x+k]
	 * 				Data[a][b][c+k][d] == InputData[x+k*DimD]
	 * @param RangeA 
	 * @param RangeB 
	 * @param RangeC 
	 * @param RangeD 
	 */
	void SetData(const vector<float>& InputData, 
	        M4DRange RangeA, M4DRange RangeB, 
	        M4DRange RangeC, M4DRange RangeD) 
	{
		RangeA.SetNMax(DimA);
		RangeB.SetNMax(DimB);
		RangeC.SetNMax(DimC);
		RangeD.SetNMax(DimD);
		size_t NElements = InputData.size();
		size_t count=0;
		for (uint a=RangeA.LowBound;a<RangeA.HighBound;++a) {
			for (uint b=RangeB.LowBound;b<RangeB.HighBound;++b) {
				for (uint c=RangeC.LowBound;c<RangeC.HighBound;++c) {
					for (uint d=RangeD.LowBound;d<RangeD.HighBound;++d) {
						if (count < NElements) {
							Data[a][b][c][d]=InputData[count];
							++count;
						}
					}
				}
			}
		}
	};
	
	/** Get Data from Matrix
	 * 
	 * @param OutputData is a return vector. \n
	 * 				pay attention to data order!\n
	 * 				Data[a][b][c][d] == OutputData[x] 
	 * 				Data[a][b][c][d+k] == OutputData[x+k]
	 * 				Data[a][b][c+k][d] == OutputData[x+k*DimD]
	 * @param RangeA 
	 * @param RangeB 
	 * @param RangeC 
	 * @param RangeD 
	 */
	void GetData(vector<float>& OutputData, 
	             M4DRange RangeA, M4DRange RangeB, 
	             M4DRange RangeC, M4DRange RangeD) const
	{
		RangeA.SetNMax(DimA);
		RangeB.SetNMax(DimB);
		RangeC.SetNMax(DimC);
		RangeD.SetNMax(DimD);
		size_t NElements = RangeA.GetDiff() 
			* RangeB.GetDiff() 
			* RangeC.GetDiff() 
			* RangeD.GetDiff();
		OutputData.resize(NElements);
		size_t count=0;
		for (uint a=RangeA.LowBound;a<RangeA.HighBound;++a) {
			for (uint b=RangeB.LowBound;b<RangeB.HighBound;++b) {
				for (uint c=RangeC.LowBound;c<RangeC.HighBound;++c) {
					for (uint d=RangeD.LowBound;d<RangeD.HighBound;++d) {
						if (count < NElements) {
							OutputData[count] = Data[a][b][c][d];
							++count;
						}
					}
				}
			}
		}
	};
	
	
};


#endif // _H_MATRIX4D