/***************************************
	Auteur : Pierre Aubert
	Mail : aubertp7@gmail.com
	Licence : CeCILL-C
****************************************/

#include "OptionParser.h"
#include "PTensor.h"

#include "Naive/naive_propagation.h"
#include "MatrixHdf5.h"

///Create the OptionParser of this program
/**	@return OptionParser of this program
*/
OptionParser createOptionParser(){
	OptionParser parser(true, __PROGRAM_VERSION__);
	parser.setExampleLongOption("naive_gray_scott --killrate=0.062 --feedrate=0.03 --nbimage=100 --nbrow=1080 --nbcol=1920 --output=outputFile.hdf5");
	parser.setExampleShortOption("naive_gray_scott -k 0.062 -f 0.03 -n 100 -r 1080 -c 1920 -o outputFile.hdf5");
	
	float killRate(0.054f), feedRate(0.014f);
	size_t nbImage(100lu), nbRow(1080lu), nbCol(1920lu);
	parser.addOption("killrate", "k", killRate, "rate of the process which converts V into P");
	parser.addOption("feedrate", "f", feedRate, "rate of the process which feeds U and drains U, V and P");
	parser.addOption("nbimage", "n", nbImage, "number of images to be created");
	size_t nbExtraStep(1lu);
	parser.addOption("nbextrastep", "e", nbExtraStep, "number of extra steps to be computed between images");
	
	parser.addOption("nbrow", "r", nbRow, "number of rows of the images to be created");
	parser.addOption("nbcol", "c", nbCol, "number of columns of the images to be created");
	float dt(1.0f);
	parser.addOption("deltat", "t", dt, "time interval between two computation");
	
	std::string defaultOutputFile("output.h5");
	parser.addOption("output", "o", defaultOutputFile, "Output file to be created with results");
	
// 	0.014,0.054,
// 	0,#000000,0.2,#00FF00,0.21,#FFFF00,0.4,#FF0000,0.6,#FFFFFF
	
	return parser;
}

///Swap two values
/**	@param[out] a : value will be b
 * 	@param[out] b : value will be a
*/
template<typename T>
void swapValue(T & a, T & b){
	T c(a);
	a = b;
	b = c;
}

///Simulate the images
/**	@param nbRow : number of rows of the images to be created
 * 	@param nbCol : number of columns of the images to be created
 * 	@param nbImage : number of images to be created
 * 	@param nbExtraStep : number of extra steps to be computed between images
 * 	@param killRate : rate of the process which converts V into P
 * 	@param feedRate : rate of the process which feeds U and drains U, V and P
 * 	@param dt : time interval between two computation
 * 	@param outputFile : name of the file to be created
 * 	@return true on succsess, false otherwise
*/
bool simulateImage(size_t nbRow, size_t nbCol, size_t nbImage, size_t nbExtraStep, float killRate, float feedRate, float dt, const std::string & outputFile){
	std::cout << "simulateImage : nbRow = " << nbRow << ", nbCol = " << nbCol << std::endl;
	
	MatrixHdf5 fullMat;
	fullMat.setAllDim(nbCol, nbRow);
	fullMat.resize(nbImage);
	
	PTensor<float> tmpInU(AllocMode::ALIGNED, nbRow, nbCol);
	PTensor<float> tmpInV(AllocMode::ALIGNED, nbRow, nbCol);
	PTensor<float> tmpOutU(AllocMode::ALIGNED, nbRow, nbCol);
	PTensor<float> tmpOutV(AllocMode::ALIGNED, nbRow, nbCol);
	tmpInU.fill(1.0f);
	tmpOutU.fill(1.0f);
	tmpInV.fill(0.0f);
	tmpOutV.fill(0.0f);
	
	for(size_t i((2lu*nbRow)/5lu); i < (3lu*nbRow)/5lu; ++i){
		for(size_t j((2lu*nbCol)/5lu); j < (3lu*nbCol)/5lu; ++j){
			tmpInU.setValue(i, j, 0.0f);
			tmpInV.setValue(i, j, 1.0f);
		}
	}
	
	float *tmpU1 = tmpInU.getData(), *tmpU2 = tmpOutU.getData();
	float *tmpV1 = tmpInV.getData(), *tmpV2 = tmpOutV.getData();
	
	long nbStencilRow(3l), nbStencilCol(3l);
	
	float diffudionRateU(0.01f), diffusionRateV(0.5f);
// 	float matDeltaSquare[] = 	{0.05f, 0.2f, 0.05f,
// 					0.2f, -1.0f, 0.2f,
// 					0.05f, 0.2f, 0.05f};
	float matDeltaSquare[] = 	{1.0f, 1.0f, 1.0f,
					1.0f, 1.0f, 1.0f,
					1.0f, 1.0f, 1.0f};
					
	for(size_t i(0lu); i < nbImage; ++i){
		std::cout << "simulateImage n°" << i << "..." << std::endl;
		for(size_t j(0lu); j < nbExtraStep; ++j){
			float * matInputU = tmpU2;
			float * matInputV = tmpV2;
			float * matOutputU = tmpU1;
			float * matOutputV = tmpV1;
	// 		if(i%2lu == 0lu){
				matInputU = tmpU1;
				matInputV = tmpV1;
				matOutputU = tmpU2;
				matOutputV = tmpV2;
	// 		}
			naive_propagation(matOutputU, matOutputV, matInputU, matInputV, nbRow, nbCol,
					matDeltaSquare, nbStencilRow, nbStencilCol,
					diffudionRateU, diffusionRateV, feedRate, killRate, dt);
			
			fullMat.setRow(i, matOutputV);
			
			///Let's swap the pointer
	// 		swapValue(tmpU1, tmpU2);
	// 		swapValue(tmpV1, tmpV2);
			
			tmpInU = tmpOutU;
			tmpInV = tmpOutV;
		}
	}
	std::cerr << "Done" << std::endl;
	//Let's save the output file
	fullMat.write(outputFile);
	return true;
}


int main(int argc, char** argv){
	OptionParser parser = createOptionParser();
	parser.parseArgument(argc, argv);
	
	const OptionMode & defaultMode = parser.getDefaultMode();
	float killRate(0.062f), feedRate(0.03f), dt(1.0f);
	size_t nbImage(100lu), nbRow(1080lu), nbCol(1920lu), nbExtraStep(1lu);
	defaultMode.getValue(killRate, "killrate");
	defaultMode.getValue(feedRate, "feedrate");
	defaultMode.getValue(nbImage, "nbimage");
	defaultMode.getValue(nbExtraStep, "nbextrastep");
	defaultMode.getValue(nbRow, "nbrow");
	defaultMode.getValue(nbCol, "nbcol");
	defaultMode.getValue(dt, "deltat");
	
	
	std::string outputFile("");
	defaultMode.getValue(outputFile, "output");
	
	bool b(simulateImage(nbRow, nbCol, nbImage, nbExtraStep, killRate, feedRate, dt, outputFile));
	return b - 1;
}