// Code reads SWAN Hs output and converts to fort.63 format. (The SWAN INPUT file is included as a refrence)
// Also the parametric .63 files are compared to the SWAN results
// in terms of Mean Absolute Error (MAE) and Normalized Mean
// Absolute Error (NMAE), with outliers excluded. 
// The final output is a result file with the statistics as well
// as an XLS .csv format for these results.
//
// MAKE SURE TO ADJUST NDSETSE TO MATCH NUMBER OF TIME STEPS
//
// By: Samuel Carter Boyd

#include <cmath>  
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <stdlib.h>
#include <cstdlib>
#include <math.h> 
#include <vector>
#include <sstream>
#include <algorithm>
#include <ctime>
#include <bits/stdc++.h> 

using namespace std;

int line=0;
int x=0;
int x1=0;
int x2=0;
int timestep=0;
char content[10],temp1[10],temp2[10],content1[10],content2[10];
string s,s1,s2; 
int node_tot,dir_tot,node,k,k1,k2,direction,num_depths,headers,num_shared;

double pi = 4*atan(1);
bool check;

double g=9.7918;	//gravity at meanlat (28.1091 deg)
double u,Hw;
int di;
int NDSETSE,NP,DTDP,NSPOOLGE,IRTYPE,TIME,IT;

double residual_sum,RMSD_avg,RMSD_avg_sum,RMSD_avg_TOT,SI,SI_sum,SI_sum_TOT;
double depth_cur,depth_prev,denominator,HwINTERP,DS,HsSWAN,H_PARAM,H_SWAN;
double HsSWANavg,HsPARAMavg,DSavg;
double DS_sum,num_skipped,HsPARAM_sum,HsSWAN_sum,MAE_sum,MAE_TOT,NMAE_sum,NMAE_TOT=0;
double HsSWANavg_sum,HsPARAMavg_sum,DSavg_sum,HsPARAMavg_TOT,HsSWANavg_TOT,DSavg_TOT;
double MAE[36],NMAE[36];
int time_skipped=0;

double Savg[36];
double Pavg[36],MAEavg[36],NMAEavg[36];


int main(){

	ifstream infile3;
	infile3.open("SWAN_const_U_30_Hs.txt");		//reading in SWAN Hs outputs

	ofstream outfile4("SWAN_Hs_const.63");	//creating fort.63 file for output

NDSETSE=36;	//number of outputs snapshots requested
NP=126772;
DTDP=1;
NSPOOLGE=3600;
IRTYPE=1;
TIME=0;
IT=1;


outfile4<<"SWAN Generated Hs  U="<<"VARYING  Dir= 0(deg N)"<<"  from Mesh.grd"<<endl;
outfile4<<NDSETSE<<"     "<<NP<<"     "<<DTDP*(NSPOOLGE)<<"     "<<NSPOOLGE<<"     "<<IRTYPE<<endl;
outfile4<<NSPOOLGE<<"  "<<NSPOOLGE<<endl;


	if(infile3.is_open()){

		cout<<"Reading Hs Data (SWAN)"<<endl;
		cout<<"Time Step: "<<IT<<" at "<<NSPOOLGE<<" seconds"<<endl;
		line=0;	

		while(getline(infile3,s)){

			if(line==126772){

				line=0;
				IT=IT+1;					//model time step number
				TIME=NSPOOLGE*IT; 				//model time (in seconds)
			outfile4<<TIME<<"  "<<TIME<<endl;//should be TIME<<"  "<<IT<<endl; if they are not equal
			cout<<"Time Step: "<<IT<<" at "<<TIME<<" seconds"<<endl;

			}

			line=line+1;

			stringstream ss;
			ss<<s;			//read the line as string (s)
			ss>>content;

			HsSWAN=atof(content);

			if(HsSWAN==-9){

				HsSWAN=-99999;

			}

			outfile4<<line<<"  "<<HsSWAN<<endl;

		}

	}
	outfile4.close();


	ifstream infile4;
	infile4.open("WEMO_Hcum_const_wind.63");			//PARAM .63 file

	ifstream infile5;
	infile5.open("SWAN_Hs_const.63");			//SWAN .63 file

	ofstream outfile5("Result_const_WEMO.txt");		//DS results

	ofstream outfile6("XLS_Result_const_WEMO.txt");	//results in excel format

	if(infile4.is_open() && infile5.is_open()){

		cout<<"Both Files Open"<<endl;

		x=0;
		line=0;
		timestep=1;

		while(getline(infile4,s1) && getline(infile5,s2)){

			line=line+1;

			if(timestep==1){

				headers=3;
			
			}


			if(timestep>1){

				headers=1;

			}

			if(line>headers && line<=126772+headers){	//ignoring header lines

				stringstream ss1;
				ss1<<s1;			//read the line as string (s1)
				ss1>>content1;

				while(!ss1.eof()){	//while values in the line (string ss1)

					ss1>>temp1;	//read each value

					H_PARAM=atof(temp1);

				}

				stringstream ss2;
				ss2<<s2;			//read the line as string (s2)
				ss2>>content2;

				while(!ss2.eof()){	//while values in the line (string ss2)

					ss2>>temp2;	//read each value

					node=atof(content2);
					H_SWAN=atof(temp2);

				}
				
				//cout<<"H_PARAM: "<<H_PARAM<<"  "<<"node: "<<node<<endl;
				//cout<<"H_SWAN : "<<H_SWAN<<"  "<<"node: "<<node<<endl;


				if(H_PARAM>0 && H_SWAN>0){ //if both models have H values NOT being -99999 

	//cout<<node<<" H_PARAM: "<<H_PARAM<<" H_SWAN: "<<H_SWAN<<endl;
					DS=abs(H_SWAN-H_PARAM);
					DS_sum=DS_sum+DS;
					HsPARAM_sum=HsPARAM_sum+H_PARAM;
					HsSWAN_sum=HsSWAN_sum+H_SWAN;
					residual_sum=residual_sum+pow((H_SWAN-H_PARAM),2);
					num_shared=num_shared+1;

				}


				else{	//otherwise at least 1 value is -99999 and therefore dont account for it

					//DS=0;
					num_skipped=num_skipped+1;

				}


			}	//close if lines < or > headers


			if(line==126772+headers){	//last value in a timestep


				HsSWANavg=HsSWAN_sum/num_shared;	//each timestep avg
				HsPARAMavg=HsPARAM_sum/num_shared;
				DSavg=DS_sum/num_shared;

				MAE[timestep]=DSavg;
				NMAE[timestep]=100*(MAE[timestep]/HsSWANavg);

				RMSD_avg=pow((residual_sum/num_shared),0.5);
				SI=100*RMSD_avg/HsSWANavg;

				if(num_shared==0){

					cout<<"NO NODES SHARED"<<endl;

					HsSWANavg=0;
					HsPARAMavg=0;
					DSavg=0;
					MAE[timestep]=0;
					NMAE[timestep]=0;
					RMSD_avg=0;
					SI=0;

					time_skipped=time_skipped+1;

				}


				if(time_skipped>NDSETSE){
					time_skipped=0;			//to prevent error
				}

				HsSWANavg_sum=HsSWANavg_sum+HsSWANavg;	//sum of AVGs for total simulation avg
				HsPARAMavg_sum=HsPARAMavg_sum+HsPARAMavg;
				DSavg_sum=DSavg_sum+DSavg;
			
				MAE_sum=MAE_sum+MAE[timestep];
				NMAE_sum=NMAE_sum+NMAE[timestep];

				RMSD_avg_sum=RMSD_avg_sum+RMSD_avg;
				SI_sum=SI_sum+SI;
									//total simulation averages

				Savg[timestep]=HsSWANavg;
				Pavg[timestep]=HsPARAMavg;


//cout<<"timestep: "<<timestep<<" MAE: "<<fixed<<setprecision(3)<<MAE[timestep]<<"(m)  NMAE: "<<fixed<<setprecision(1)<<NMAE[timestep]<<"(%)"<<endl;	//<<fixed<<setprecision(3)<<" RMSD: "<<RMSD_avg<<"(m)"<<" SI: "<<fixed<<setprecision(1)<<SI<<"(%)"<<endl;

	outfile5<<setw(20)<<"**********TIMESTEP: "<<timestep<<"**********"<<endl;
	outfile5<<setw(10)<<"          |   SWAN   | "<<setw(10)<<"   PARAM   "<<endl;
	outfile5<<setw(6)<<"Hs(avg)[m]| "<<fixed<<setprecision(6)<<HsSWANavg;
	outfile5<<" | "<<setw(6)<<fixed<<setprecision(6)<<HsPARAMavg<<endl;
	outfile5<<"   MAE    | "<<fixed<<setprecision(3)<<MAE[timestep]<<"(m)"<<endl;
	outfile5<<"  NMAE    | "<<fixed<<setprecision(1)<<NMAE[timestep]<<"(%)"<<endl;


	cout<<setw(20)<<"**********TIMESTEP: "<<timestep<<"**********"<<endl;
	cout<<setw(10)<<"          |   SWAN   | "<<setw(10)<<"   PARAM   "<<endl;
	cout<<setw(6)<<"Hs(avg)[m]| "<<fixed<<setprecision(6)<<HsSWANavg;
	cout<<" | "<<setw(6)<<fixed<<setprecision(6)<<HsPARAMavg<<endl;
	cout<<"   MAE    | "<<fixed<<setprecision(3)<<MAE[timestep]<<"(m)"<<endl;
	cout<<"  NMAE    | "<<fixed<<setprecision(1)<<NMAE[timestep]<<"(%)"<<endl;


				if(timestep==NDSETSE){ 	//total simulation avg
				
					if(time_skipped<0 || time_skipped>NDSETSE){	//to prevent error

						time_skipped=0;

					}

					HsPARAMavg_TOT=HsPARAMavg_sum/(timestep-time_skipped);
					HsSWANavg_TOT=HsPARAMavg_sum/(timestep-time_skipped);
					DSavg_TOT=DSavg_sum/(timestep-time_skipped);

					MAE_TOT=MAE_sum/(timestep-time_skipped);
					NMAE_TOT=NMAE_sum/(timestep-time_skipped);

					RMSD_avg_TOT=RMSD_avg_sum/(timestep-time_skipped);
					SI_sum_TOT=SI_sum/(timestep-time_skipped);
				

//sorting arrays for outlier investigation
//outliers are found if they are outside the
//inter-quartile range (for NMAE)	
					sort(MAE, MAE + (NDSETSE-time_skipped));
					sort(NMAE, NMAE + (NDSETSE-time_skipped));


					int median_index=floor((NDSETSE-time_skipped)/2);
					int Q1_index=floor(median_index/2);
					int Q3_index=floor(median_index+(median_index/2));
					double median=NMAE[median_index];
					double Q1=NMAE[Q1_index];
					double Q3=NMAE[Q3_index];
					double IQR=Q3-Q1;
					double NMAE_stat[NDSETSE];
					double MAE_stat[NDSETSE];
					double Savg_stat[NDSETSE];
					double Pavg_stat[NDSETSE];


					int j=0;	//used to count new array length without outliers
	
				cout<<endl;
			cout<<"(NMAE) median: "<<median<<"  Q1: "<<Q1<<"  Q3: "<<Q3<<"  IQR: "<<IQR<<endl;
				cout<<"outliers above NMAE: "<<Q3+(1.5*IQR)<<" and below "<<Q1-(1.5*IQR)<<endl; 
				cout<<endl;

				outfile5<<endl;
				outfile5<<"median: "<<median<<"  Q1: "<<Q1<<"  Q3: "<<Q3<<"  IQR: "<<IQR<<endl;
			   outfile5<<"outliers above NMAE: "<<Q3+(1.5*IQR)<<" and below "<<Q1-(1.5*IQR)<<endl; 
				outfile5<<endl;

			    		for (int i=1; i<=NDSETSE-time_skipped; i++){ 

//non-outliers
						if(Q1-(1.5*IQR)<NMAE[i] && NMAE[i]<Q3+(1.5*IQR)){

							j=j+1;
		
							NMAE_stat[j]=NMAE[i];
							MAE_stat[j]=MAE[i];
							Savg_stat[j]=Savg[i];
							Pavg_stat[j]=Pavg[i];

//cout<<"MAE_stat["<<j<<"]: "<<fixed<<setprecision(5)<<MAE_stat[j]<<"  NMAE_stat["<<j<<"]: "<<NMAE_stat[j]<<"  Savg_stat["<<j<<"]: "<<Savg_stat[j]<<"  Pavg_stat["<<j<<"]: "<<Pavg_stat[j]<<endl;

						}

						else{
						    cout<<"NMAE: "<<NMAE[i]<<" is an outlier"<<endl;
						    outfile5<<"NMAE: "<<NMAE[i]<<" is an outlier"<<endl;
						} 
					}
				
					double MAE_stat_sum=0;
					double MAE_stat_avg;
					double NMAE_stat_sum=0;
					double NMAE_stat_avg;
					double Savg_stat_sum=0;
					double Savg_stat_avg;
					double Pavg_stat_sum=0;
					double Pavg_stat_avg;

					for(int i=1; i<=j; i++){

						MAE_stat_sum=MAE_stat_sum+MAE_stat[i];
						NMAE_stat_sum=NMAE_stat_sum+NMAE_stat[i];
						Savg_stat_sum=Savg_stat_sum+Savg_stat[i];
						Pavg_stat_sum=Pavg_stat_sum+Pavg_stat[i];
				
					}

					MAE_stat_avg=MAE_stat_sum/j;
					NMAE_stat_avg=NMAE_stat_sum/j;
					Savg_stat_avg=Savg_stat_sum/j;
					Pavg_stat_avg=Pavg_stat_sum/j;



//cout<<"MAE_stat_sum: "<<MAE_stat_sum<<"  NMAE_stat_sum: "<<NMAE_stat_sum<<"  Savg_stat_sum: "<<Savg_stat_sum<<"  Pavg_stat_sum: "<<Pavg_stat_sum<<" j: "<<j<<endl;
//cout<<"MAE_stat_avg: "<<MAE_stat_avg<<"  NMAE_stat_avg: "<<NMAE_stat_avg<<"  Savg_stat_avg: "<<Savg_stat_avg<<"  Pavg_stat_avg: "<<Pavg_stat_avg<<" j: "<<j<<endl;

outfile5<<fixed<<setprecision(5)<<endl;
outfile5<<"Total simulation avgs (outliers included)"<<endl;
outfile5<<"HSWAN: "<<HsSWANavg_TOT<<"(m)  HPARAM: "<<HsPARAMavg_TOT<<"(m)  MAE: "<<MAE_TOT<"(m)";
outfile5<<"  N.MAE: "<<NMAE_TOT<<endl;
//"(%)  RMSD: "<<RMSD_avg_TOT<<"(m)"<<" SI: "<<SI_sum_TOT<<"(%)"<<endl;
outfile5<<endl;
outfile5<<"Total simulation avgs (without outliers)"<<endl;
outfile5<<"HSWAN: "<<Savg_stat_avg<<"(m)  HPARAM: "<<Pavg_stat_avg<<"(m)  MAE: "<<MAE_stat_avg<"(m)";
outfile5<<"  N.MAE: "<<NMAE_stat_avg<<endl;
//"(%)  RMSD: "<<RMSD_avg_TOT<<"(m)"<<" SI: "<<SI_sum_TOT<<"(%)"<<endl;
outfile5<<endl;

cout<<fixed<<setprecision(5)<<endl;
cout<<"Total simulation avgs (outliers included)"<<endl;
cout<<"HSWAN: "<<HsSWANavg_TOT<<"(m)  HPARAM: "<<HsPARAMavg_TOT<<"(m)  MAE: "<<MAE_TOT<<"(m)";
cout<<"  N.MAE: "<<NMAE_TOT<<"(%)"<<endl;
//"(%)  RMSD: "<<RMSD_avg_TOT<<"(m)"<<" SI: "<<SI_sum_TOT<<"(%)"<<endl;
cout<<endl;
cout<<"Total simulation avgs (without outliers)"<<endl;
cout<<"HSWAN: "<<Savg_stat_avg<<"(m)  HPARAM: "<<Pavg_stat_avg<<"(m)  MAE: "<<MAE_stat_avg<<"(m)";
cout<<"  N.MAE: "<<NMAE_stat_avg<<"(%)"<<endl;
//"(%)  RMSD: "<<RMSD_avg_TOT<<"(m)"<<" SI: "<<SI_sum_TOT<<"(%)"<<endl;
cout<<endl;


				break;

				}			//close if timestep==NDSETSE

				HsSWANavg=0;			//resetting values
				HsPARAMavg=0;
				DSavg=0;
				//NMAE=0;
				num_shared=0;
				num_skipped=0;
				HsPARAM_sum=0;
				HsSWAN_sum=0;
				residual_sum=0;
				DS_sum=0;
				node=0;
				H_PARAM=0;
				H_SWAN=0;
				line=0;
				SI=0;
				timestep=timestep+1;

			}	//close if line<126772+headers
				

		}	//close while both files are reading content

		infile4.close();
		infile5.close();
		outfile5.close();

	}	//close if both files open

//XLS output (all values including outliers)

	for(int i=1; i<=NDSETSE; i++){
		if(i==1){
			outfile6<<"SWAN,PARAM,MAE[m],NMAE[%]"<<endl;
		}
		outfile6<<Savg[i]<<","<<Pavg[i]<<","<<MAE[i]<<","<<NMAE[i]<<endl;
	}


outfile6.close();


return 0;
}