#include "gr4j_sl_v1.h"
#include "model_tools.h"
#include <fstream>

Cgr4j_v1::Cgr4j_v1()
{
	//-----------------------------------------------------
	// 1. Info...
	Set_Name("GR4J Version 1.0");

	Set_Author("copyrights (c) 2007 Stefan Liersch");

	Set_Description("The Rainfall-Runoff Model GR4J \n \n \n"
		"Reference: \n \n"
		"Perrin, C., Michel, C., Andreassian, V. (2003). \n"
		"Improvement of a parsimonious model for streamflow simulation \n"
		"Journal of Hydrology, 279, p. 275-289. \n"
	);

	///////////////////////////////////////////////////////////////////
	// FIRST MODULE DIALOG
	_CreateDialog1();
	///////////////////////////////////////////////////////////////////
}

Cgr4j_v1::~Cgr4j_v1()
{
}

bool Cgr4j_v1::On_Execute()
{
	//---------------------------------------------------------------
	//				PARAMETERS
	//---------------------------------------------------------------
	CSG_Table		*pTable, *pTable_out, *pTable_Settings;
	// table field numbers
	int				dateField;
	int				dischargeField;
	int				pcpField;
	int				pet_tmpField;
	// used time series (0 = ETp, 1 = temperature)
	int				Choice_ETp_TMP;

	CSG_String		date1, date2;
	int				first, last;
	double			Latitude;
	double			area;

	// Assign Parameters from Module Dialog
	//---------------------------------------------------------
	pTable			= Parameters("TABLE")			->asTable();
	// Field numbers
	dateField		= Parameters("DATE_Field")		->asInt();
	dischargeField	= Parameters("DISCHARGE_Field")	->asInt();
	pcpField		= Parameters("PCP_Field")		->asInt();
	pet_tmpField		= Parameters("ETP_TMP_Field")		->asInt();
	// used time series
	Choice_ETp_TMP	= Parameters("ETP_TMP_TS")		->asInt();
	// Parameters
	X.X1			= Parameters("X1")				->asDouble();
	X.X2			= Parameters("X2")				->asDouble();
	X.X3			= Parameters("X3")				->asDouble();
	X.X4			= Parameters("X4")				->asDouble();
	//---------------------------------------------------------

	//---------------------------------------------------------
	// Create Second Dialog and assign linear module parameters
	_CreateDialog2(pTable, dateField, date1, date2, Choice_ETp_TMP, Latitude, area);
	//---------------------------------------------------------

	//---------------------------------------------------------
	// searching the first and the last record of the time range
	gr4j.AssignFirstLastRec(*pTable, first, last, date1, date2, dateField);
	sizeAll = last - first + 1;
	//---------------------------------------------------------

	//---------------------------------------------------------
	// initialize struct Input
	Input.date.resize(sizeAll);
	Input.Q_obs = new double[sizeAll];
	Input.precipitation = new double[sizeAll];
	Input.pET = new double[sizeAll];
	Input.size = sizeAll;
	if (Choice_ETp_TMP == 1) temperature = new double[sizeAll];
	//---------------------------------------------------------
	// assign data to struct Input
	for (int i = 0; i < sizeAll; i++)
	{
		Input.date[i].append(pTable->Get_Record(i)->asString(dateField));
		// convert Q_obs from m3/s to mm/day
		Input.Q_obs[i] = (pTable->Get_Record(i)->asDouble(dischargeField) * 86.4 / area);
		Input.precipitation[i] = (pTable->Get_Record(i)->asDouble(pcpField));
		if (Choice_ETp_TMP == 0) // if ETp data was selected
			Input.pET[i] = (pTable->Get_Record(i)->asDouble(pet_tmpField));
		else
			temperature[i] = (pTable->Get_Record(i)->asDouble(pet_tmpField));
	}
	// if no ETp data available the temperature time series is used
	// to calculate ETp after Hamon (1961)
	if (Choice_ETp_TMP == 1) // if temperature data was selected
	{
		int first_year, month, day;
		_Get_Date_Int(Input.date[0], first_year, month, day);
		// calculate ETp
		ETp = new CpET(temperature, sizeAll, Latitude, first_year);
		// assign ETp to Input.pET
		for (int i = 0; i < sizeAll; i++)
			Input.pET[i] = ETp->Get_pET_Hamon(i);
	}
	//---------------------------------------------------------

	//---------------------------------------------------------
	//		SIMULATION
	gr4j.Simulation(Input, X);
	//---------------------------------------------------------
	
	double NSE = model_tools::CalcEfficiency(Input.Q_obs, gr4j.Q_sim, sizeAll);

	//---------------------------------------------------------
	// write output in table
	pTable_out = SG_Create_Table();
	_CreateTableSim(pTable_out, Input, gr4j.Q_sim, area);
	pTable_out->Set_Name("GR4J_Sim");
	Parameters("TABLEout")->Set_Value(pTable_out);

	pTable_Settings = SG_Create_Table();
	_CreateTableSettings(pTable_Settings, X, NSE, area, Latitude);
	pTable_out->Set_Name("GR4J_Settings");
	Parameters("TABLEsettings")->Set_Value(pTable_Settings);

/*
	ofstream file("GR4J_out.txt");
	file << "Date, pcp, pET, Q_obs, Q_sim" << endl;
	for (int i = 0; i < sizeAll; i++)
	{
		file << Input.date[i].c_str() << ", ";
		file << Input.precipitation[i] << ", ";
		file << Input.pET[i] << ", ";
		file << Input.Q_obs[i] * area / 86.4 << ", ";
		file << gr4j.Q_sim[i] * area / 86.4 << endl;
	}
	file.close();
*/

	Input.date.resize(0);
	delete[] Input.pET;
	delete[] Input.precipitation;
	delete[] Input.Q_obs;
	
	if (Choice_ETp_TMP == 1)
	{
		delete[] temperature;
		delete ETp;
	}

	
	return(true);
}

///////////////////////////////////////////////////////////////////////
//
//                          CREATE OUTPUT TABLES
//
///////////////////////////////////////////////////////////////////////

void Cgr4j_v1::_CreateTableSim(CSG_Table *pTable,
							   InputData& Input, double* Q_sim,
							   double area) // convert mm to m3/s
{
	int i = 0; // used in function Get_Record(i)
	CSG_Table_Record	*pRecord;

	// creating the column titles
	pTable->Add_Field("Date",		TABLE_FIELDTYPE_String);
	pTable->Add_Field("Flow OBS",	TABLE_FIELDTYPE_Double);
	pTable->Add_Field("Flow SIM",	TABLE_FIELDTYPE_Double);

	for (int j = 0; j < Input.size; j++)
	{
		pTable->Add_Record();
		pRecord = pTable->Get_Record(i);

		// writing the data into the rows
		pRecord->Set_Value(0,Input.date[j].c_str());
		pRecord->Set_Value(1,Input.Q_obs[j] * area / 86.4);
		pRecord->Set_Value(2,Q_sim[j] * area / 86.4);
		i++;
	}
}
//---------------------------------------------------------------------

void Cgr4j_v1::_CreateTableSettings(CSG_Table *pTable_Settings, GR4J_Parameters X, 
									double NSE, double area, double latitude)
{
	int i = 0; // used in function Get_Record(i)
	CSG_Table_Record	*pRecord;

	// creating the column titles
	pTable_Settings->Add_Field("Parameters",	TABLE_FIELDTYPE_String);
	pTable_Settings->Add_Field("Settings",	TABLE_FIELDTYPE_Double);

	pTable_Settings->Add_Record();
	pRecord = pTable_Settings->Get_Record(0);
	pRecord->Set_Value(0,"Capacity of the production store [mm]");
	pRecord->Set_Value(1,X.X1);	

	pTable_Settings->Add_Record();
	pRecord = pTable_Settings->Get_Record(1);
	pRecord->Set_Value(0,"Water exchange coefficient [mm]");
	pRecord->Set_Value(1,X.X2);	

	pTable_Settings->Add_Record();
	pRecord = pTable_Settings->Get_Record(2);
	pRecord->Set_Value(0,"Capacity of the routing store [mm]");
	pRecord->Set_Value(1,X.X3);	

	pTable_Settings->Add_Record();
	pRecord = pTable_Settings->Get_Record(3);
	pRecord->Set_Value(0,"Time base of the unit hydrograph [days]");
	pRecord->Set_Value(1,X.X4);	

	pTable_Settings->Add_Record();
	pRecord = pTable_Settings->Get_Record(4);
	pRecord->Set_Value(0,"Area");
	pRecord->Set_Value(1,area);

	pTable_Settings->Add_Record();
	pRecord = pTable_Settings->Get_Record(5);
	pRecord->Set_Value(0,"Latitude");
	pRecord->Set_Value(1,latitude);

	pTable_Settings->Add_Record();
	pRecord = pTable_Settings->Get_Record(6);
	pRecord->Set_Value(0,"Nash-Sutcliffe");
	pRecord->Set_Value(1,NSE);


}
//---------------------------------------------------------------------
// Converting a date string 'YYYYMMDD' to three separate integer values

void Cgr4j_v1::_Get_Date_Int(string date, int &year, int &month, int &day)
{
	// divide date string to year, month, and day string
	string y(date.c_str(),0,4);
	string m(date.c_str(),4,2);
	string d(date.c_str(),6,2);

	// convert to integer
	sscanf(y.c_str(),"%d",&year);
	sscanf(m.c_str(),"%d",&month);
	sscanf(d.c_str(),"%d",&day);
}
//---------------------------------------------------------------------

///////////////////////////////////////////////////////////////////////
//
//		                          DIALOGS
//
///////////////////////////////////////////////////////////////////////

void Cgr4j_v1::_CreateDialog1()
{
	CSG_Parameter	*pNode;
	CSG_String s;

	//---------------------------------------------------------------
	// Input Table
	pNode = Parameters.Add_Table(
		NULL	, "TABLE"	, _TL("Table"),
		"",
		PARAMETER_INPUT
	);

	Parameters.Add_Table_Field(
		pNode	, "DATE_Field"	, _TL("Date Column"),
		"Select the column containing the Date"
	);

	Parameters.Add_Table_Field(
		pNode	, "DISCHARGE_Field"	, _TL("Streamflow Column [m3/s]"),
		"Select the Column containing Discharge Values"
	);

	Parameters.Add_Table_Field(
		pNode	, "PCP_Field"	, _TL("Precipitation Column [mm]"),
		"Select the Column containing precipitation Values"
	);

	Parameters.Add_Table_Field(
		pNode	, "ETP_TMP_Field"	, _TL("ETp or Temperature Column"),
		"Select the Column containing Temperature Values"
	);

	Parameters.Add_Choice(
		pNode	, "ETP_TMP_TS"		, "Specify which time series you are using!",
		"",
		CSG_String::Format("%s|%s|",
			_TL("ETp [mm]"),			//  0
			_TL("mean air temperature [°C]")	//  1 
		)
	);
	//---------------------------------------------------------------

	//---------------------------------------------------------------
	// Input / Calibration Parameters
	s.Printf("Node2", 2);
		pNode = Parameters.Add_Node(NULL,s,"Free Parameters","");

	Parameters.Add_Value(
		pNode,	"X1",	_TL("Capacity of the production store [mm]"),
		_TL("(x1) Capacity of the production store [mm]"),
		PARAMETER_TYPE_Double,
		0.0, 0.0, true, 50000.0, true
	);
	Parameters.Add_Value(
		pNode,	"X2",	_TL("Water exchange coefficient [mm]"),
		_TL("(x2) Water exchange coefficient [mm]"),
		PARAMETER_TYPE_Double,
		0.0, -100.0, true, 100.0, true
	);
	Parameters.Add_Value(
		pNode,	"X3",	_TL("Capacity of the routing store [mm]"),
		_TL("(x3) Capacity of the routing store [mm]"),
		PARAMETER_TYPE_Double,
		0.0, 0.0, true, 50000.0, true
	);
	Parameters.Add_Value(
		pNode,	"X4",	_TL("Time base of the unit hydrograph [days]"),
		_TL("(x4) Time base of the unit hydrograph [days]"),
		PARAMETER_TYPE_Double,
		0.0, 0.0, true, 50000.0, true
	);
	//---------------------------------------------------------------

	// SAGA output tables
	Parameters.Add_Table_Output(
		NULL	, "TABLEout"	, _TL("Table"),
		""
	);

	Parameters.Add_Table_Output(
		NULL	, "TABLEsettings"	, _TL("Table"),
		""
	);

	//---------------------------------------------------------------
}

//-------------------------------------------------------------------
//						DIALOG 2
//-------------------------------------------------------------------

void Cgr4j_v1::_CreateDialog2(CSG_Table* pTable, int dateField,
							  CSG_String& date1, CSG_String& date2,
							  int Choice_ETp_TMP,
							  double& Latitude, double& area)
{
	CSG_String		s;
	CSG_Parameters	P;  // used to add Parameters in the second dialog
	CSG_Parameter*	pNode;

	P.Set_Name(_TL("Choose Time Range"));

	s.Printf("Node1", 1);
	pNode = P.Add_Node(NULL,s,"Time Range","");

	s.Printf("FDAY" , 1-1);
	P.Add_String(pNode,s,_TL("First Day"),_TL(""),
				 pTable->Get_Record(0)->asString(dateField));

	s.Printf("LDAY" , 1-2);
	P.Add_String(pNode,s,_TL("Last Day"),_TL(""),
				 pTable->Get_Record(pTable->Get_Record_Count()-1)->asString(dateField));

	s.Printf("Node2", 2);
	pNode = P.Add_Node(NULL,s,"Catchment Area","");

	P.Add_Value(
		pNode,	"AREA",	_TL("Catchment Area in [km2]"),
		_TL("Catchment area in square kilometers, used to convert m3/s to mm"),
		PARAMETER_TYPE_Double,
		1.0, -0.0, true, 10000.0, true
	);

	if (Choice_ETp_TMP ==1)
	{
		s.Printf("Node2-1", 2-1);
		pNode = P.Add_Node(NULL,s,"Latitude","");

		P.Add_Value(
			pNode,	"LATITUDE",	_TL("Latitude of temperature station in [degrees]"),
			_TL("Latitude of temperature station in [degrees]"),
			PARAMETER_TYPE_Double,
			50.0, -90.0, true, 90.0, true
		);
	}

	if( SG_UI_Dlg_Parameters(&P, _TL("Choose Time Range")) )
	{
		///////////////////////////////////////////////////////////////
		//
		//                ASSIGN DATA FROM SECOND DIALOG
		//
		///////////////////////////////////////////////////////////////
		date1		= P(CSG_String::Format("FDAY",dateField).c_str())
									->asString();
		date2		= P(CSG_String::Format("LDAY",dateField).c_str())
									->asString();

		area		= P("AREA")			->asDouble();

		if (Choice_ETp_TMP ==1)
			Latitude = P("LATITUDE")	->asDouble();

	}

}