#include "calibrate.h"
#include <math.h>
#include <iostream>
#include <fstream>
#include <string>
#include <iostream>
#include <sstream>

//CRF are CRiteria Functions of estimating the goodness of fit for observed and modeled hydrographs, explained in Wang, J. PhD Dissertation at SUNY ESF (2002)
//CRF1 = for emphasizing peak flow, based on Nash and Sutcliffe (1970) Eq 3 and Perrin et al (2001) Eq 4; From 1970 Eq 3 R^2 = (Fo^2 - F^2)/Fo^2 = 1 - F^2/Fo^2, where F^2 = (qmod - qobs)^2, Fo^2 = (qobs - qobs_mean)^2
//CRF2 = for emphasizing baseflow, based on Ye et al. (1997) Eq E, and derived by Perrin et al (2001) Eq 8; From 1997 Eq E = 1 - mu_e / sigma_q, where mu_e = 2nd moment of residuals ei = (qobs - qmod), sigma_q = variance of qobs
//CRF3 = for balancing peak and baseflow, based on Chiew and McMahon (1994) Eq 1, and derived by Perrin et al (2001) Eq 5; From 1994 OBJ = (sqrt(qsim) - sqrt(qobs))^2

//Nash, J. E., & Sutcliffe, J. V. (1970). River flow forecasting through conceptual models. Part I - a discussion of principles. Journal of Hydrology, 27(3), 282-290. 
//Ye, W., Bates, B. C., Viney, N. R., Silvapan, M., & Jakeman, A. J. (1997). Performance of conceptual rainfall-runoff models in low-yielding ephemeral catchments. Water Resources Research, 33(1), 153-166. 
//Chiew, F. H. S., & McMahon, T. (1994). Application of the daily rainfall-runoff model MODHYDROLOG to 28 Australian catchments. Journal of Hydrology, 153, 383-416. 
//Perrin, C., Michel, C., & Andreassian, V. (2001).Does a large number of parameters enhance model performance ? Comparative assessment of common catchment model structures on 429 catchments.Journal of Hydrology, 242(3 - 4), 275 - 301. doi : 10.1016 / s0022 - 1694(00)00393 - 0


//calibrate function will read in Runoff_Observed_for_Calibration.csv and use for comparison with HydroPlus generated Q values
calibrate::calibrate(Inputs* input, string Directory_Input_CLArg, double startDate, double endDate)
{
	CRF1 = 0.0;
	CRF2 = 0.0;
	CRF3 = 0.0;
	SE = 0.0;
	VE = 0.0;
	string temp;
	string HrMinSec;
	double tempQ = 0.0;
	int size = 0;
	double date = -1.0;
	double tempQobs;
	double date_InputDataRow_HH;
	double date_InputDataRow_GDH;

	//Runoff_Observed_vec_m (m) vector defined for local operations, using input->Runoff_Observed_m (m) vector
	Runoff_Observed_vec_m = input->Runoff_Observed_m;

	//timeStep_Calibration_str set to TimeStep_CompareToObservedRunoff_str
	string timeStep_Calibration_str = input->SimulationStringParams["TimeStep_CompareToObservedRunoff_str"];
	//If timeStep_Calibration_str is Hour or hour then calculate day, week, and month time step versions of runoff
	if (timeStep_Calibration_str == "Hour" || timeStep_Calibration_str == "hour") {
		calcDailyQobs(input);
		calcWeeklyQobs(input);
		calcMonthlyQobs(input);
	}
}
//calibrate::objFunRunoff function receives simulated runoff vector, Qi_Simulated_m, which is vector of Qi_Routed_m (m) or unrouted Qi_m runoff
//Note: HydroPlusConfig.xmls element Flag_CompareToObservedRunoff can equal 0, 1 (routed runoff), or 2 (unrouted runoff)
void calibrate::objFunRunoff(Inputs* input, vector<double>& Qi_Simulated_m) {
	
	double SUMQobs = 0.0;
	double SSQ = 0.0;

	double VARQ1 = 0.0;
    double VARQ2 = 0.0;
    double VARQ3 = 0.0;

	double aveQobs =0.0;
	double sqrtAveQobs =0.0;
    double sqrtSUMQobs =0.0;
	int size = 1;
	double sF1 = 0.0;
	double sF2 = 0.0;
	double sF3 = 0.0;
	
	CRF1 =0.0;
	CRF2 = 0.0;
	CRF3 =0.0;
	SE = 0.0;
	VE = 0.0;

	int k =0;
	int j=0;
	
	//size defined with Ternary operator ?: If condition at left is true, then first item after ? taken, otherwise item after : is taken
	//Note: Ternary operator explained: For (x_var<=y_var?x_var:y_var), if x_var<=y_var is true, then x_var is taken, otherwise y_var taken.
	size = (Qi_Simulated_m.size() > Runoff_Observed_vec_m.size()) ?  Runoff_Observed_vec_m.size() : Qi_Simulated_m.size();
	
	//calculate sum of observed flow
	for(int i=0; i < size; i++) {
		SUMQobs += Runoff_Observed_vec_m[i];
		sqrtSUMQobs += sqrt(Runoff_Observed_vec_m[i]);
	}

	//CRF1 = Eq 4 in Perrin et al. (2001) is Nash Sutcliffe (1970) Efficiency for peak flow
	//CRF2 = Eq 8 in Perrin et al. (2001) is Ye et al. (1997) Efficiency for baseflow
	//CRF3 = Eq 5 in Perrin et al. (2001) is Chiew and McMahon (1994) Efficiency for balancing peak and baseflow

	//calculate average observed flow
	aveQobs = SUMQobs/size;         
    sqrtAveQobs = sqrtSUMQobs/size;
	for(int i=0; i<size; i++) {

		sF1 += (Qi_Simulated_m[i] - Runoff_Observed_vec_m[i])*(Qi_Simulated_m[i] - Runoff_Observed_vec_m[i]);
		sF2 += fabs(Runoff_Observed_vec_m[i] - Qi_Simulated_m[i]);
        sF3 += (sqrt(Qi_Simulated_m[i]) - sqrt(Runoff_Observed_vec_m[i]))*(sqrt(Qi_Simulated_m[i]) - sqrt(Runoff_Observed_vec_m[i]));
		
		VARQ1 += (Runoff_Observed_vec_m[i] - aveQobs)*(Runoff_Observed_vec_m[i] - aveQobs);
	    VARQ2 += fabs(Runoff_Observed_vec_m[i] - aveQobs);	
		VARQ3 += (sqrt(Runoff_Observed_vec_m[i]) - sqrtAveQobs)*(sqrt(Runoff_Observed_vec_m[i]) - sqrtAveQobs);
	}
	//If size equals 1 then VARQ1 VARQ2 or VARQ3 were -nan, due to aveQobs = 1st item in vector, and equation above 
	if (size == 1) {
		VARQ1 = Runoff_Observed_vec_m[0];
		VARQ2 = Runoff_Observed_vec_m[0];
		VARQ3 = Runoff_Observed_vec_m[0];
	}

	// mean squared error
	SE = sF1 / size;
	// 1-SUM(Qobs-Qsim)^2/SUM(Qobs-aveQobs)^2
	CRF1 = 1 - sF1 / VARQ1;
	// 1-SUM|Qobs-Qsim|/SUM|Qobs-aveQobs|
	CRF2 = 1 - sF2 / VARQ2;
	// 1-SUM(sqrt(Qobs)-sqrt(Qsim))^2/SUM(sqrt(Qobs)-sqrt(aveQobs))^2
	CRF3 = 1 - sF3 / VARQ3;
	// 1 - SUM(abs(Qsim - Qobs)) / SUM (Qobs)
	VE = 1 - sF2 / SUMQobs;
}

//Note: Consider comments to better explain function
double calibrate::operator [](size_t i) const
{	
	static double safe;
	if(i<0 || Runoff_Observed_vec_m.size()<=i) {
		cout << "index out of bounds in calibrate Qobs[], i: "<<endl<<endl;
		return safe;
	}
    return Runoff_Observed_vec_m[i];
}	

//Note: Consider refactor so Runoff_Observed_vec_m is replaced by input->Runoff_Observed_m
//calibrate::calcDailyQobs function converts input->Runoff_Observed_m from time step hour to day
void calibrate::calcDailyQobs(Inputs* input){
	//dailySize defined as Runoff_Observed_vec_m.size() / 24
	int dailySize = Runoff_Observed_vec_m.size() / 24;
	//dailyQobs vector resized to dailySize
	dailyQobs.resize(dailySize);
	//For loop for i from 0 to < dailySize
	for (int i = 0; i<dailySize; i++) {
		//Variable_Day_sum reset to 0
		double Variable_Day_sum = 0.0;
		//For Loop for j from 0 to < 24
		for (int j = 0; j<24; j++) {
			//Variable_Day_sum increased by Runoff_Observed_vec_m[i * 24 + j], hour flow values summed for 1 day
			Variable_Day_sum = Variable_Day_sum + Runoff_Observed_vec_m[i * 24 + j];
		}
		//dailyQobs vector stores Variable_Day_sum
		dailyQobs[i] = Variable_Day_sum;
	}
}

//calibrate::calcWeeklyQobs function converts input->Runoff_Observed_m from time step hour to week
void calibrate::calcWeeklyQobs(Inputs* input){
	//weeklySize defined as Runoff_Observed_vec_m.size() / (24 * 7)
	int weeklySize = Runoff_Observed_vec_m.size() / (24 * 7);
	//weeklyQobs vector resized to weeklySize
	weeklyQobs.resize(weeklySize);
	//For loop for i from 0 to < weeklySize
	for (int i = 0; i<weeklySize; i++) {
		//Variable_Day_sum reset to 0
		double weeklySum = 0.0;
		//For Loop for j from 0 to < (24 * 7)
		for (int j = 0; j<(24 * 7); j++) {
			//Variable_Day_sum increased by Runoff_Observed_vec_m[i * 24 * 7 + j], hour flow values summed for 1 week
			weeklySum = weeklySum + Runoff_Observed_vec_m[i * 24 * 7 + j];
		}
		//weeklyQobs vector stores Variable_Day_sum
		weeklyQobs[i] = weeklySum;
	}
}

//calibrate::calcMonthlyQobs function converts input->Runoff_Observed_m from time step hour to month
void calibrate::calcMonthlyQobs(Inputs* input){
	//monthlySize defined as Runoff_Observed_vec_m.size() / (24 * 30)
	int monthlySize = Runoff_Observed_vec_m.size() / (24 * 30);
	//monthlyQobs vector resized to monthlySize
	monthlyQobs.resize(monthlySize);
	//For loop for i from 0 to < monthlySize
	for (int i = 0; i< monthlySize; i++) {
		//Variable_Day_sum reset to 0
		double monthlySum = 0.0;
		//For Loop for j from 0 to < (24 * 30)
		for (int j = 0; j<(24 * 30); j++) {
			//Variable_Day_sum increased by Runoff_Observed_vec_m[i * 24 * 30 + j], hour flow values summed for 1 week
			monthlySum = monthlySum + Runoff_Observed_vec_m[i * 24 * 30 + j];
		}
		//monthlyQobs vector stores Variable_Day_sum
		monthlyQobs[i] = monthlySum;
	}
}

//calibrate::calibrate function defined in calibrate.h as calibrate(vector<double>& Qsim, vector<double>& Qobs);
//Note: calibrate should be renamed to be more specific calibrate_Runoff_Simulated_vs_Observed
calibrate::calibrate(Inputs* input, vector<double>& Qi_Simulated_m, vector<double>& Qobs)
{
	//Runoff_Observed_vec_m (m) defined second time by input->Runoff_Observed_m
	Runoff_Observed_vec_m = input->Runoff_Observed_m;

	CRF1 = 0.0;
	CRF2 = 0.0;
	CRF3 = 0.0;
	SE = 0.0;
	VE = 0.0;

	double tempQ = 0.0;
	int size = 0;

	double SUMQobs = 0.0;
	double SSQ = 0.0;

	double VARQ1 = 0.0;
	double VARQ2 = 0.0;
	double VARQ3 = 0.0;

	double aveQobs = 0.0;
	double sqrtAveQobs = 0.0;
	double sqrtSUMQobs = 0.0;

	double sF1 = 0.0;
	double sF2 = 0.0;
	double sF3 = 0.0;

	int k = 0;
	int j = 0;

	size = Qi_Simulated_m.size();

	//Note: Vector error thrown if Qi_Simulated_m.size() different than Weather.csv
	//calculate sum of observed flow
	for (int i = 0; i<size; i++) {
		SUMQobs += Runoff_Observed_vec_m[i];
		sqrtSUMQobs += sqrt(Runoff_Observed_vec_m[i]);
	}

	//CRF1 = Eq 4 in Perrin et al. (2001) is Nash Sutcliffe (1970) Efficiency for peak flow
	//CRF2 = Eq 8 in Perrin et al. (2001) is Ye et al. (1997) Efficiency for baseflow
	//CRF3 = Eq 5 in Perrin et al. (2001) is Chiew and McMahon (1994) Efficiency for balancing peak and baseflow

	aveQobs = SUMQobs / size;         //calculate average observed flow
	sqrtAveQobs = sqrtSUMQobs / size;
	for (int i = 0; i<size; i++) {
		sF1 += (Qi_Simulated_m[i] - Runoff_Observed_vec_m[i])*(Qi_Simulated_m[i] - Runoff_Observed_vec_m[i]);
		sF2 += fabs(Runoff_Observed_vec_m[i] - Qi_Simulated_m[i]);
		sF3 += (sqrt(Qi_Simulated_m[i]) - sqrt(Runoff_Observed_vec_m[i]))*(sqrt(Qi_Simulated_m[i]) - sqrt(Runoff_Observed_vec_m[i]));

		VARQ1 += (Runoff_Observed_vec_m[i] - aveQobs)*(Runoff_Observed_vec_m[i] - aveQobs);
		VARQ2 += fabs(Runoff_Observed_vec_m[i] - aveQobs);
		VARQ3 += (sqrt(Runoff_Observed_vec_m[i]) - sqrtAveQobs)*(sqrt(Runoff_Observed_vec_m[i]) - sqrtAveQobs);
	}

	//If size equals 1 then VARQ1 VARQ2 or VARQ3 were -nan, due to aveQobs = 1st item in vector, and equation above 
	if (size == 1) {
		VARQ1 = Runoff_Observed_vec_m[0];
		VARQ2 = Runoff_Observed_vec_m[0];
		VARQ3 = Runoff_Observed_vec_m[0];
	}

	// mean squared error
	SE = sF1 / size;
	// 1-SUM(Runoff_Observed_vec_m-Qsim)^2/SUM(Qobsi-aveQobs)^2
	CRF1 = 1 - sF1 / VARQ1;
	// 1-SUM|Runoff_Observed_vec_m-Qsim|/SUM|Qobsi-aveQobs|
	CRF2 = 1 - sF2 / VARQ2;
	// 1-SUM(sqrt(Runoff_Observed_vec_m)-sqrt(Qsim))^2/SUM(sqrt(Runoff_Observed_vec_m)-sqrt(aveQobs))^2
	CRF3 = 1 - sF3 / VARQ3;
	// 1 - SUM(abs(Qsim - Runoff_Observed_vec_m)) / SUM (Runoff_Observed_vec_m)
	VE = 1 - sF2 / SUMQobs;
}