// MatData.cpp
#include <omp.h>
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <iomanip>
#include <math.h>
#include "pugixml.hpp"
#include "MatData.h"
#include "easylogging++.h"

#ifdef _MSC_VER // tonipat:20131112
    #define SSCANF sscanf_s
#else
    #define SSCANF sscanf
			
#endif

namespace math1d_cl
{
	// Config constructor
	MatData::MatData(std::string fileName)
	{
		const std::string measuredDischargeVolumesCsvFileName = "MeasuredDischargeVolumes.csv";
		const std::string precipitationsCsvFileName = "Precipitations.csv";

		const std::string qCsvFileName = "Q.csv";
		const std::string hCsvFileName = "H.csv";
		const std::string vCsvFileName = "V.csv";

		m_configFilePath = fileName;
		
		CLOG(INFO, "model") << "Config file: " << m_configFilePath;

		pugi::xml_document doc;
		pugi::xml_parse_result result = doc.load_file(m_configFilePath.c_str());
		 
		if(result)
		{
			pugi::xml_node model = doc.child("conf").child("math1D");
			//int logTarget = model.child("logTarget").text().as_int(); TODO: handle log targe for easyloggingcpp
			std::string logFilePath = model.child("logFilePath").text().as_string();
			// m_logger = new Logger(logTarget, logFilePath);

			std::string resourcesPath = doc.child("conf").child("resourcesPath").text().as_string();

			std::string matDataXmlFileName = model.child("matDataXmlFileName").text().as_string();

			m_matDataXmlFilePath = resourcesPath + "/" + matDataXmlFileName;
			m_measuredDischargeVolumesCsvFilePath = resourcesPath + "/" + measuredDischargeVolumesCsvFileName;
			m_precipitationsCsvFilePath = resourcesPath + "/" + precipitationsCsvFileName;
			
			m_qCsvFilePath = resourcesPath + "/" + qCsvFileName;
			m_hCsvFilePath = resourcesPath + "/" + hCsvFileName;
			m_vCsvFilePath = resourcesPath + "/" + vCsvFileName;
		}
		else
		{
			CLOG(FATAL,"model") << "Config file " << m_configFilePath << " not loaded!";
			std::exit(-1);
		}

		CLOG(DEBUG, "model") << "Measured discharge file: " << m_measuredDischargeVolumesCsvFilePath;
		CLOG(DEBUG, "model") << "Precipitations file: " << m_precipitationsCsvFilePath;
		CLOG(DEBUG, "model") << "Q output file: " << m_qCsvFilePath;
		CLOG(DEBUG, "model") << "H output file: " << m_hCsvFilePath;
		CLOG(DEBUG, "model") << "V output file: " << m_vCsvFilePath;
	}//MatData

	// Default constructor
	/*MatData::MatData()
	{
		#ifdef _MSC_VER // tonipat:20131112
			m_configFilePath = "..\\Resources\\Config.xml";
			////getchar();
		#else
			m_configFilePath = "..//Resources//Config.linux.xml";	
			////getchar();
		#endif

		CLOG(INFO, "model") << "Config file: " << m_configFilePath;

		pugi::xml_document doc;
		 
		pugi::xml_parse_result result = doc.load_file(m_configFilePath.c_str());
		 
		if(result)
		{
			int logTarget = doc.child("math1D").child("logTarget").text().as_int();
			std::string logFilePath = doc.child("math1D").child("logFilePath").text().as_string();
			m_logger = new Logger(logTarget, logFilePath);

			m_matDataXmlFilePath = doc.child("math1D").child("matDataXmlFilePath").text().as_string();
			m_measuredDischargeVolumesCsvFilePath = doc.child("math1D").child("measuredDischargeVolumesCsvFilePath").text().as_string();
			m_precipitationsCsvFilePath =  doc.child("math1D").child("precipitationsCsvFilePath").text().as_string();
			m_qCsvFilePath =  doc.child("math1D").child("qCsvFilePath").text().as_string();
			m_hCsvFilePath =  doc.child("math1D").child("hCsvFilePath").text().as_string();
			m_vCsvFilePath =  doc.child("math1D").child("vCsvFilePath").text().as_string();
		}
		else
		{
			CLOG(ERROR,"model") << "Config file " << m_configFilePath << " not loaded!";
			std::exit(-1);
	}*/

	// Copy constructor
	MatData::MatData(const MatData& origin)
	{
		size_t i = 0;
		// Source stations
		for ( i = 0; i < origin.m_sourceStations.size(); i++)
		{
			this->m_sourceStations.push_back(*(new std::shared_ptr<Station>(new Station(*origin.m_sourceStations[i]))));
		}
		// River stations
		for (  i = 0; i < origin.m_riverStations.size(); i++)
		{
			this->m_riverStations.push_back(*(new std::shared_ptr<Station>(new Station(*origin.m_riverStations[i]))));
		}
		// Basin ID
		this->m_basinId = origin.m_basinId;
		// Weather stations
		for ( i = 0; i < origin.m_weatherStations.size(); i++)
		{
			this->m_weatherStations.push_back(*(new std::shared_ptr<Station>(new Station(*origin.m_weatherStations[i]))));
		}
		// Measure stations
		for (  i = 0; i < origin.m_measureStations.size(); i++)
		{
			this->m_measureStations.push_back(*(new std::shared_ptr<Station>(new Station(*origin.m_measureStations[i]))));
		}
		// Channels
		for ( i = 0; i < origin.m_channels.size(); i++)
		{
			this->m_channels.push_back(*(new std::shared_ptr<Channel>(new Channel(*origin.m_channels[i]))));
		}
		// Subbasins
		for (  i = 0; i < origin.m_subbasins.size(); i++)
		{
			this->m_subbasins.push_back(*(new std::shared_ptr<Subbasin>(new Subbasin(*origin.m_subbasins[i]))));
		}
		// Measured discharge volume
		this->m_measuredDischargeVolumes = origin.m_measuredDischargeVolumes;
		// Precipitations
		this->m_precipitations = origin.m_precipitations;
		// Measured hydrographs Q
		this->m_measuredHydrographsQ = origin.m_measuredHydrographsQ;
		// Options
		this->m_options = origin.m_options;
		// Logger
		//Logger logger = *origin.m_logger;
		// TODO: Preserve log target 
		//this->m_logger = new Logger(*(origin.m_logger));

		this->m_configFilePath = origin.m_configFilePath;
		this->m_matDataXmlFilePath = origin.m_matDataXmlFilePath;
		this->m_measuredDischargeVolumesCsvFilePath = origin.m_measuredDischargeVolumesCsvFilePath;
		this->m_precipitationsCsvFilePath = origin.m_precipitationsCsvFilePath;
		this->m_qCsvFilePath = origin.m_qCsvFilePath;
		this->m_hCsvFilePath = origin.m_hCsvFilePath;
		this->m_vCsvFilePath = origin.m_vCsvFilePath;
		this->m_nTimeSteps = origin.m_nTimeSteps;
		this->m_minuteStep = origin.m_minuteStep;
		this->m_precStationIds = origin.m_precStationIds;
		this->m_mdvStationIds = origin.m_mdvStationIds;
	}
	
	int MatData::runRR()
	{
		
		//TIMED_FUNC(rrtimer);

		// Set default options
		setOptions();

		// Load schematization and input data
		collectMatDataCsv();
		//PERFORMANCE_CHECKPOINT_WITH_ID(rrtimer, "loading csv data");

		CLOG(INFO,"model") << "Solving Rainfall-Runoff simulation...";
		CLOG(INFO,"model") << "Start time: " << printDateTime(m_options.getStartDate());
		CLOG(INFO,"model") << "End time: " << printDateTime(m_options.getEndDate());
		CLOG(INFO,"model") << "Minute step: " << intToString(m_options.getMinuteStep());
		CLOG(INFO,"model") << "Scheme name: " << m_options.getRrSchemeName();
		CLOG(INFO,"model") << "Meteorologic model: " << intToString(m_options.getMeteoModelId());
		CLOG(INFO,"model") << "Creation time: " << printDateTime(m_options.getCreationDate());
		
		// compute model
		rainfallRunoffModel();
		//PERFORMANCE_CHECKPOINT_WITH_ID(rrtimer, "RR simulation");

		// save results
		collectRRResultCsv();
		//PERFORMANCE_CHECKPOINT_WITH_ID(rrtimer, "saving results");

		return m_options.getSimulationId();
	}

	void MatData::scsMethod(std::vector<double>& q, size_t& index)
	{
		double f = m_subbasins[index]->getArea();
		double cn = m_subbasins[index]->getCn();
		double q0 = m_subbasins[index]->getBaseflow();
		int weatherStationIndex = m_subbasins[index]->getWeatherStationIndex();
		double lost = 0;

		// weather station index must be index in MatData.xml not id!
		// value 144 just for testing reasons
		// int wsId = 144;
		int wsId = m_weatherStations[weatherStationIndex]->getId();
		size_t order = -1;

		// Not sure if some weather station id can appear more than once in m_precStationIds!
		for(size_t i = 0; i < m_precStationIds.size(); i++)
		{
			if(m_precStationIds[i] == wsId)
			{
				order = i;
			}
		}

		std::vector<double> p(m_precipitations.size());
		// Cumulative sum of p
		std::vector<double> cumSumP(m_precipitations.size());
	    int lastNonZeroIndex = -1;

		double firstValue = m_precipitations[0].second[order];
		p[0] = firstValue;
		cumSumP[0] = firstValue;
		if(firstValue <= lost)
		{
			lastNonZeroIndex = 0;
		}

		for(size_t i = 1; i < m_precipitations.size(); i++)
		{
			p[i] = m_precipitations[i].second[order];
			cumSumP[i] = cumSumP[i - 1] + p[i];
			if(cumSumP[i] <= lost)
			{
				lastNonZeroIndex = i;
			}
		}

		if(lastNonZeroIndex != -1)
		{
			double sumP = 0;
			for(size_t i = 0; i < (size_t) lastNonZeroIndex + 1; i++)
			{
				if(lastNonZeroIndex < ((int)p.size() - 1))
				{
					sumP += p[i];
				}
				p[i] = 0;
			}
			if(lastNonZeroIndex < ((int)p.size() - 1))
			{
				p[lastNonZeroIndex + 1] -= lost - sumP;
			}
		}


		/* Unit hydrograph */
		double sumQ1 = 0.7 + 0.96 + 1 + 0.75 + 0.55 + 0.35;
		std::vector<double> q1 = {0, 0.7, 0.96, 1, 0.75, 0.55, 0.35};
		/*------------------*/

		/* Hypodermic flow */
		for(int i = 0; i < 73; i++)
		{
			double value = 0.18 * exp(-0.05*i);
			q1.push_back(value);
			sumQ1 += value;
		}
		/*------------------*/

		for(int i = 0; i < 80; i++)
		{
			// Take subbasin area into account
			// Suspected magic ??? 3.6 ???
			q1[i] = q1[i] / sumQ1 * f / 3.6;
		}

		std::vector<double> CN(m_precipitations.size(), cn);
		std::vector<double> p24(m_precipitations.size(), 0);
		size_t m = m_precipitations.size();
		size_t n = q1.size();
		std::vector<double> Q(m + n - 1); // Suspicious dimension ??? 

		for(size_t i = 0; i < m_precipitations.size(); i++)
		{
			//int iL = 0 < i - 120 ? i - 120 : 0;
			// Compute which values are predicted
			int iL = 0 < (int )(i - 120 + 1) ? (int )(i - 120 + 1 ): 0;
			double sumPTmp = 0;

			for(size_t j = iL; j < i + 1; j++)
			{
				sumPTmp += p[j];
			}

			if(sumPTmp > 53) // CHMI Official value = 40
			{
				// When cumulative precipitations for previous 5 days are
				// bigger than 53 <- ??? HOW, WHY ???, change CN, just a little bit...
				CN[i] = cn * exp(0.00673 * (100 - cn));
			}

			// Cumulative precipitations for last 24 hours
			iL = 0 < (int)(i - 24 + 1) ? (int ) (i - 24 + 1 ): 0;
			sumPTmp = 0;
			for(size_t j = iL; j < i + 1; j++)
			{
				sumPTmp += p[j];
			}
			p24[i] = sumPTmp;


			// Retention potential (S)
			std::vector<double> rp(m_precipitations.size());

			// Initial abstraction (Ia)
			std::vector<double> r1(m_precipitations.size());

			// Intermediate result
			std::vector<double> h24ef(m_precipitations.size());

			// Percentual infulence
			std::vector<double> proc(m_precipitations.size());

			// Effective height of overland flow
			std::vector<double> pef(m_precipitations.size());

			double h24 = 50; // <- ???? Educated guess
			for(size_t j = 0; j < (size_t) m; j++)
			{
				rp[j] = 25.4 * (1000 / CN[j] - 10);
				r1[j] = 0.2 * rp[j];
				//h24ef[j] = (h24 > r1[j]) * (pow((h24 - r1[j]), 2) / (h24 + rp[j] - r1[j]));
				h24ef[j] = (h24 > r1[j]) * ((h24 - r1[j]) * (h24 - r1[j])) / (h24 + rp[j] - r1[j]);
				proc[j] = h24ef[j] / h24;
				pef[j] = proc[j] * p[j];
			}

			/*
			testing data
			m = 5;
			n = 3;
			double u[] = {1,2,3,4,5};
			double v[] = {1,2,3};
			std::vector<double> d(m + n - 1, 0);*/

			// Convulation
			
			//for(int j = 0; j < m + n + 1; j++)
			//{
			//	int k = 0 > (j - n + 1) ? 0 : (j - n + 1);
			//	int kTo = (j < (m)) ? j : m;
			//	for(k; k < kTo + 1; k++)
			//	{
			//		//Q[j] += pef[k] * q1[j - k - 1];
			//		d[j] += u[k] * v[j - k];
			//	}
			//}

		
			/* Convolution of precipitations and unit hydrograph*/			
			int j=0, j1=0, k=0;
			//#pragma omp parallel for schedule(auto)
			for(j = 0; j < ((int)(m + n) - 1); j++)
			{
			    j1 = j;
				Q[j]=0;	
				
				for( k = 0; k < (int)n; k++)
				{
					if(j1 >= 0 && j1 < (int)m)
					{
					Q[j]  += (pef[j1]*q1[k]);
					
					}//if
					j1--;

				}//fork
			}//forj
		}//?

		/* Add base flow to computed runoff */
		for(size_t i = 0; i < (size_t ) m; i++)
		{
			q[i] = Q[i] + q0;
		}
	}

	void MatData::rainfallRunoffModel()
	{
		// Clear hydrograms
		clearResults();
		
		int rectangleProfile = 0;
		if(m_options.getRrProfileShape() == "rectangle")
		{
			rectangleProfile = 1;
		}
		
		CLOG(DEBUG, "model") << "Runoff computations...";
		
		// Iterate through channels
		for(size_t i = 0; i < m_channels.size(); i++)
		{
			CLOG(DEBUG, "model") << "Computing channel "  << i;
			//m_logger->log("Computing channel " + intToString((int )i));
			
			std::vector<double> qIn(m_nTimeSteps, 0);
			std::vector<double> qOut(m_nTimeSteps, 0);
			std::vector<double> hIn(m_nTimeSteps, 0);

			// Initial channels
			if(m_channels[i]->getUpstreams().size() == 0)
			{
				// Subbasin contribution
				if(m_options.getRrType() == SCS_CN)
				{
					scsMethod(qOut, i);
					m_channels[i]->getHydrograph().setQOut(qOut);
				}
			}
			else // Ordinary channels
			{
				m_channels[i]->getHydrograph().setQOut(qOut);

				// Add upstreams contributions to QIn
				for(size_t upstream = 0; upstream < m_channels[i]->getUpstreams().size(); upstream++)
				{
					int index = m_channels[i]->getUpstreams()[upstream];

					if(index >= 0)
					{
						for(int j = 0; j < m_nTimeSteps; j++ )
						{
							qIn[j] += m_channels[index]->getHydrograph().getQOut()[j];
						}
					}
				}

				m_channels[i]->getHydrograph().setQIn(qIn);

				// Compute corresponding Hin
				m_channels[i]->getHydrograph().setHIn(qToH(qIn, m_channels[i], rectangleProfile));

				// Compute channel contributions
				switch(m_options.getRrHdType())
				{
				case SV1D:
					break;
				case KWA_Comsol:
					break;
				case KWA_FV:
					break;
				case Vel:
					{
						double v = velocity(m_channels[i]);
						double it = floor(m_channels[i]->getLength() / (3600 * v) + 0.5);
						for(int j = 0; j < it; j++)
						{
							qOut[j] = qIn[0];
						}
						for(int j = (int )it; j < m_nTimeSteps; j++)
						{
							qOut[j] = qIn[j - (int )it];
						}
						m_channels[i]->getHydrograph().setQOut(qOut);
					}
					break;
				case FDM:
					break;
				}

				// Add subbasin contributions
				std::vector<double> qSub(m_nTimeSteps, 0);

				switch(m_options.getRrType())
				{
				case SCS_CN:
					size_t subbasin = m_channels[i]->getSubbasinIndex();
					scsMethod(qSub, subbasin);
					break;
				}

				for(int j = 0; j < m_nTimeSteps; j++)
				{
					m_channels[i]->getHydrograph().getQOut()[j] += qSub[j];
				}
			}

			// Fitting
			if(m_options.getFitting())
			{
				for(size_t j = 0; j < m_measureStations.size(); j++)
				{
					if(m_measureStations[j]->getChannelIndex() == m_channels[i]->getSubbasinIndex())
					{
						// remove all zero elements over last nonzero element
						size_t size = 0;
						for(size_t k = m_measuredHydrographsQ[j].size(); k-- > 0;)
						{
							while(m_measuredHydrographsQ[j][k] == 0)
								continue;
							size = k + 1;
							break;
						}

						std::vector<double> measQ(size, 0);
						std::vector<double> compQ(size, 0);
						std::vector<double> wageQ(size, 1);

						std::vector<double> added = m_channels[i]->getHydrograph().getQOut();

						for(size_t k = 0; k < size; k++)
						{
							measQ[k] = m_measuredHydrographsQ[j][k];
							compQ[k] = added[k];
						}

						double alpha, beta;
						fit(alpha, beta, measQ, compQ, wageQ);

						for(int k = 0; k < m_nTimeSteps; k++)
						{
							added[k] = alpha * added[k] + beta;
							if(added[k] <= 0)
							{
								added[k] = 0.1;
							}
						}

						m_channels[i]->getHydrograph().setQOut(added);
					}
				}
			}
		}

		// Iterate through channels
		for(size_t i = 0; i < m_channels.size(); i++)
		{
			if(m_channels[i]->getUpstreams().size() == 1 && m_channels[i]->getUpstreams()[0] == 0)
			{
				int index = m_channels[i]->getDownstream();
				if(index < 0)
					continue;
				m_channels[i]->getHydrograph().setHOut(qToH(m_channels[i]->getHydrograph().getQOut(),
					m_channels[index], rectangleProfile));
			}
			else
			{
				m_channels[i]->getHydrograph().setHOut(qToH(m_channels[i]->getHydrograph().getQOut(),
					m_channels[i],	rectangleProfile));
			}
		}

		m_simulationDone = true;
	}

	void MatData::collectRRResultCsv()
	{
		CLOG(INFO, "model") << "Saving results...";
		int nChannels =(int ) m_channels.size();

		std::ofstream qFile(m_qCsvFilePath.c_str());
		// std::ofstream vFile(m_vCsvFilePath.c_str()); // only zeros
		std::ofstream hFile(m_hCsvFilePath.c_str());

		std::string firstLine = "time\\id;";

		for(size_t i = 0; i < (size_t ) nChannels; i++)
		{
			firstLine += intToString(m_channels[i]->getStationId());
			if(i <(unsigned int ) (nChannels - 1))
			{
				firstLine += ";";
			}
			else
			{
				firstLine += "\n";
			}
		}

		qFile << firstLine;
		//vFile << firstLine;
		hFile << firstLine;

		for(int i = 0; i < m_nTimeSteps; i++)
		{	
			qFile << printDateTime(m_precipitations[i].first) << ";";
			//vFile << printDateTime(m_precipitations[i].first);
			hFile << printDateTime(m_precipitations[i].first) << ";";
			for(size_t j = 0; j < (size_t) nChannels; j++)
			{
				qFile << std::fixed << std::setprecision(6) << m_channels[j]->getHydrograph().getQOut()[i];
				
				//if(m_channels[j]->getHydrograph().getHOut().size() > i)
				hFile << std::fixed << std::setprecision(6) << m_channels[j]->getHydrograph().getHOut()[i];
				//else
				//	hFile << ";";
				if(j < (unsigned int )(nChannels - 1))
				{
					qFile << ";  ";
					hFile << ";  ";
				}
			}
			qFile << "\n";
			//vFile << "\n";
			hFile << "\n";
		}

		qFile.close();
		//vFile.close();
		hFile.close();
	}

	void MatData::fit(double& alpha, double& beta, std::vector<double>& x, std::vector<double>& y, std::vector<double>& wage)
	{
		size_t size = x.size();
		std::vector<double> yTest(size, 0);
		double tmp = 0;

		for(size_t i = 0; i < size; i++)
		{
			yTest[i] = y[i] - y[0];
			tmp += yTest[i] * yTest[i];
		}

		// 2-norm of vector yTest == 0
		if(sqrt(tmp) == 0)
		{
			alpha = 1;
			double diff = 0;
			
			for(size_t i = 0; i < size; i++)
			{
				diff += x[i] - y[i];
			}
			beta = diff / size;
		}
		else
		{
			double wageX = 0;
			double wageY = 0;
			double a = 0;
			double b = 0;
			std::vector<double> yMod(size, 0);
			std::vector<double> xMod(size, 0);
			std::vector<double> yyMod(size, 0);
			std::vector<double> yxMod(size, 0);

			for(size_t i = 0; i < size; i++)
			{
				wageX += x[i];
				wageY += y[i];
			}

			for(size_t i = 0; i < size; i++)
			{
				xMod[i] = x[i] - (wageX / size);
				yMod[i] = y[i] * size - wageY;
				yxMod[i] = y[i] * xMod[i];
				yyMod[i] = y[i] * yMod[i];
				a += yxMod[i];
				b += yyMod[i] / size;
			}

			alpha = a / b;
			beta = 0;

			if(alpha < 0)
			{
				alpha = 1;
				for(size_t i = 0; i < size; i++)
				{
					beta += (x[i] - y[i]) / size;
				}
			}
			else
			{
				std::vector<double> alphaY(size);
				for(size_t i = 0; i < size; i++)
				{
					alphaY[i] = alpha * y[i];
					beta += (x[i] - alphaY[i]) / size;
				}
			}
		}
	}

	double MatData::velocity(std::shared_ptr<Channel> channel)
	{
		double n = channel->getN(); // Manning coefficient
		double iS = channel->getSlope(); // Channel slope
		double w = channel->getWidth(); // Channel width
		double iB = channel->getBankSlope(); // Channel bank slope
		double q = channel->getHydrograph().getQIn()[0];

		if(iS == 0)
		{
			iS += 0.001;
		}

		// Version 2 (rectangle profile)
		double x = q * n / (w * sqrt(iS));
		double h = pow(x, 3.0/5.0);
		double s = h * (w + h / iB);
		return q/s;
	}

	std::vector<double> MatData::qToH(std::vector<double>& q, std::shared_ptr<Channel> channel, int& rectangleProfile)
	{
		std::vector<double> h(m_nTimeSteps, 0);
		double n = channel->getN(); // Manning coefficient
		double iS = channel->getSlope(); // Channel slope
		double w = channel->getWidth(); // Channel width
		//double iB = channel->getBankSlope(); // Channel bank slope

		if(rectangleProfile == 1) // Version 1 (rectangle profile)
		{
			for(int i = 0; i < m_nTimeSteps; i++)
			{
				h[i] = pow(n / (w * sqrt(iS)) * q[i], 3/5);
			}
		}
		else // Version 1 (rectangle profile)
		{
			// not implemented yet
		}

		return h;
	}

	void MatData::setOptions()
	{
		m_options.setUri("http://release.floreon.vsb.cz:8088/WS_Database/M9_M11.asmx?WSDL");

		tm startDate;
		startDate.tm_mday = 10;
		startDate.tm_mon = 8;
		startDate.tm_year = 107;
		startDate.tm_hour = 0;
		startDate.tm_min = 0;
		startDate.tm_sec = 0;
		time_t start = mktime(&startDate);
		m_options.setStartDate(start);

		tm endDate;
		endDate.tm_mday = 15;
		endDate.tm_mon = 8;
		endDate.tm_year = 107;
		endDate.tm_hour = 10;
		endDate.tm_min = 0;
		endDate.tm_sec = 0;
		time_t end = mktime(&endDate);
		m_options.setEndDate(end);

		m_options.setMinuteStep(60);

		m_options.setMeteoModelId(2);
		m_options.setMeteoModelName("Measured");
		m_options.setRr(true);
		m_options.setHd(false);
		   
		m_options.setRrSchemeId(1);
		m_options.setRrSchemeName("RRschematizace");
		m_options.setRrModelId(3);
		m_options.setRrModelName("Math_1D");
		m_options.setRrType(SCS_CN);
		m_options.setRrHdType(Vel);
		m_options.setRrProfileShape("rectangle");
		m_options.setHdModelId(0);
		m_options.setHdType("KWA_Comsol");
		m_options.setHdProfileShape("rectangle");
		m_options.setLog(true);
			   
		m_options.setSimulationId(0);
		m_options.setFitting(true);
		m_options.setCalibtemp(true);

		time_t now = time(NULL);
		m_options.setCreationDate(now);
	}
	
	void MatData::collectMatDataCsv()
	{
		CLOG(INFO, "model") << "Loading MatData and meteodata...";
		
		// Loads MatData schematization
		MatData::loadSchematization();

		// Measured Discharge Volumes
		loadMeasuredDischargeVolumesFromCsv();
		
		// Precipitations
		loadPrecipitationsFromCsv();

		m_nTimeSteps =(int ) m_precipitations.size();
		m_minuteStep = 60; // Hardcoded
		m_options.setStartDate(m_precipitations[0].first);
		m_options.setEndDate(m_precipitations[m_precipitations.size() - 1].first);


		// Prefill vectors with default values
		for(size_t i = 0; i < m_channels.size(); i++)
		{
			for(size_t j = 0; j < m_precipitations.size(); j++)
			{
				Hydrograph h = m_channels[i]->getHydrograph();
				h.getHIn().push_back(0);
				h.getQIn().push_back(0);
				h.getHOut().push_back(0);
				h.getQOut().push_back(0);
			}
		}

		for(size_t i = 0; i < m_precipitations.size(); i++)
		{
			for(size_t j = 0; j < m_precStationIds.size(); j++)
			{
				if(m_precipitations[i].second[j] == -999)
				{
					m_precipitations[i].second[j] = -1;
				}
			}
		}

		// Creates Hydrographs
		getMeasuredHydrographs();
	}

	void MatData::loadSchematization()
	{
		CLOG(DEBUG, "model") << "Loading schematization..";
		pugi::xml_document doc;
		pugi::xml_parse_result result = doc.load_file(m_matDataXmlFilePath.c_str());
		if(result.status != pugi::status_ok)
		{
			CLOG(FATAL, "model") << "Schematization load result: " << result.description();
			std::exit(EXIT_FAILURE);
		}

		// Source stations
		pugi::xml_node sourceStations = doc.child("MatData").child("SourceStations");
		loadStationsFromXml(m_sourceStations, sourceStations);

		// River stations
		pugi::xml_node riverStations = doc.child("MatData").child("RiverStations");
		loadStationsFromXml(m_riverStations, riverStations);

		// Weather stations
		// Weather station contains only Id, Name, Code and Location
		pugi::xml_node weatherStations = doc.child("MatData").child("WeatherStations");
		loadStationsFromXml(m_weatherStations, weatherStations);

		// Channels
		pugi::xml_node channels = doc.child("MatData").child("Channels");
		loadChannelsFromXml(m_channels, channels);

		// Subbasins
		pugi::xml_node subbasins = doc.child("MatData").child("Subbasins");
		loadSubbasinsFromXml(m_subbasins, subbasins);
	}

	void MatData::loadStationsFromXml(std::vector<std::shared_ptr<Station>>& stations, pugi::xml_node& xml_stations)
	{
		for (pugi::xml_node xml_station = xml_stations.first_child(); xml_station; xml_station = xml_station.next_sibling())
		{
			std::shared_ptr<Station> station(new Station());

			station->setId(xml_station.child("Id").text().as_int());
			station->setDescription(xml_station.child_value("Description"));
			station->setName(xml_station.child_value("Name"));
			Location location;
			location.setX(xml_station.child("Location").child("X").text().as_float());
			location.setY(xml_station.child("Location").child("Y").text().as_float());
			location.setZ(xml_station.child("Location").child("Z").text().as_float());
			station->setLocation(location);
			station->setClockStep(xml_station.child("ClockStep").text().as_double());
			station->setCode(xml_station.child_value("Code"));
			station->setPrecipitationFlag(xml_station.child("PrecipitationFlag").text().as_bool());
			station->setTemperatureFlag(xml_station.child("TemperatureFlag").text().as_bool());
			station->setSnowFlag(xml_station.child("SnowFlag").text().as_bool());
			station->setWindVelocityFlag(xml_station.child("WindVelocityFlag").text().as_bool());
			station->setRadarFlag(xml_station.child("RadarFlag").text().as_bool());
			station->setIsVirtual(xml_station.child("Virtual").text().as_bool());
			station->setSpa1(xml_station.child("Spa1").text().as_double());
			station->setSpa2(xml_station.child("Spa2").text().as_double());
			station->setSpa3(xml_station.child("Spa3").text().as_double());
			station->setQ(xml_station.child("Q").text().as_double());
			//station->setOwner(xml_station.child_value("Owner")); // Not known what is owner yet.
			std::string wgcString = xml_station.child_value("WaterGaugingCategory");
			WaterGaugingCategory wgcResult;
			if(wgcString == "A") wgcResult = A;
			else if(wgcString == "B") wgcResult = B;
			else if(wgcString == "N") wgcResult = N;
			station->setWaterGaugingCategory(wgcResult);
			station->setHSpa1(xml_station.child("HSpa1").text().as_double());
			station->setHSpa2(xml_station.child("HSpa2").text().as_double());
			station->setHSpa3(xml_station.child("HSpa3").text().as_double());
			station->setH(xml_station.child("H").text().as_int());
			//station->setChmuProfileId(xml_station.child("ChmuProfileId").text().as_int());
			
			//station->setChannelIndex(xml_station.child("Channel").child("Id").text().as_int());
			station->setChannelIndex(xml_station.child("ChannelIndex").text().as_int() - 1);
			
			stations.push_back(station);
		}
	}

	void MatData::loadChannelsFromXml(std::vector<std::shared_ptr<Channel>>& channels, pugi::xml_node& xml_channels)
	{
		for (pugi::xml_node xml_channel = xml_channels.first_child(); xml_channel; xml_channel = xml_channel.next_sibling())
		{
			std::shared_ptr<Channel> channel(new Channel());

			channel->setId(xml_channel.child("Id").text().as_int());
			channel->setName(xml_channel.child_value("Name"));
			channel->setSourceStationId(xml_channel.child("SourceStationId").text().as_int());
			channel->setSourceStationIndex(xml_channel.child("SourceStationIndex").text().as_int() - 1);
			channel->setStationId(xml_channel.child("StationId").text().as_int());
			channel->setStationIndex(xml_channel.child("StationIndex").text().as_int() - 1);
			channel->setH(xml_channel.child("H").text().as_double());
			channel->setD(xml_channel.child("D").text().as_double());
			channel->setLength(xml_channel.child("Length").text().as_double());
			channel->setSlope(xml_channel.child("Slope").text().as_double());
			channel->setBankSlope(xml_channel.child("BankSlope").text().as_double());
			channel->setDepth(xml_channel.child("Depth").text().as_double());
			channel->setWidth(xml_channel.child("Width").text().as_double());
			channel->setN(xml_channel.child("N").text().as_double());
			channel->setSubbasinId(xml_channel.child("SubbasinId").text().as_int());
			channel->setSubbasinIndex(xml_channel.child("SubbasinIndex").text().as_int() - 1);

			pugi::xml_node xml_upstreams = xml_channel.child("Upstreams");
			for (pugi::xml_node xml_upstream = xml_upstreams.first_child(); xml_upstream; xml_upstream = xml_upstream.next_sibling())
			{
				channel->getUpstreams().push_back(xml_upstream.child("Index").text().as_int() - 1);
			}
			//channel->setDownstream(xml_channel.child("Downstreams").child("Downstream").child("Index").text().as_int());
			channel->setDownstream(xml_channel.child("DownstreamIndex").text().as_int() - 1);
			channels.push_back(channel);
		}
	}

	void MatData::loadSubbasinsFromXml(std::vector<std::shared_ptr<Subbasin>>& subbasins, pugi::xml_node& xml_subbasins)
	{
		for (pugi::xml_node xml_subbasin = xml_subbasins.first_child(); xml_subbasin; xml_subbasin = xml_subbasin.next_sibling())
		{
			std::shared_ptr<Subbasin> subbasin(new Subbasin());

			subbasin->setId(xml_subbasin.child("Id").text().as_int());
			subbasin->setName(xml_subbasin.child_value("Name"));
			subbasin->setArea(xml_subbasin.child("Area").text().as_double());
			subbasin->setH(xml_subbasin.child("H").text().as_double());
			subbasin->setD(xml_subbasin.child("D").text().as_double());
			subbasin->setLength(xml_subbasin.child("Length").text().as_double());
			subbasin->setSlope(xml_subbasin.child("Slope").text().as_double());
			subbasin->setBaseflow(xml_subbasin.child("BaseFlow").text().as_double());
			subbasin->setCn(xml_subbasin.child("Cn").text().as_double());
			subbasin->setN(xml_subbasin.child("N").text().as_double());
			subbasin->setLai(xml_subbasin.child("Lai").text().as_double());
			subbasin->setInitAbstraction(xml_subbasin.child("InitAbstraction").text().as_double());
			subbasin->setTimeConcentration(xml_subbasin.child("TimeConcentration").text().as_double());
			subbasin->setStorageCoeff(xml_subbasin.child("StorageCoeff").text().as_double());
			subbasin->setChannelIndex((int )(xml_subbasin.child("ChannelIndex").text().as_double() - 1));
			subbasin->setWeatherStationIndex(xml_subbasin.child("WeatherStationIndex").text().as_int() - 1); // Change to Index!

			subbasins.push_back(subbasin);
		}
	}

	void MatData::loadMeasuredDischargeVolumesFromCsv()
	{
		std::ifstream measuredDischargeVolumesFile(m_measuredDischargeVolumesCsvFilePath.c_str());
		std::string value, firstLine, mdv;
		int nColumns = 0;
		

		// Gets number of columns
		// First line is header consisting of string "Time" and Station Ids
		if(getline(measuredDischargeVolumesFile, firstLine))
		{
			std::istringstream iss(firstLine);
			while(getline(iss, value, ';'))
			{
				if(value != "Time")
				{
					m_mdvStationIds.push_back(atoi(value.c_str()));
				}
				nColumns++;
			}
		}

		// Parse data
		while(getline(measuredDischargeVolumesFile, value))
		{
			std::istringstream iss(value);
			std::vector<double> tmp;
			for(int i=0; i<nColumns; i++)
			{
				getline(iss, mdv, ';');
				if(i>0) // First token is Time, not station id
				{
					tmp.push_back(atof(mdv.c_str()));
				}
			}
			m_measuredDischargeVolumes.push_back(tmp);
		}
		measuredDischargeVolumesFile.close();
	}

	void MatData::loadPrecipitationsFromCsv()
	{
		std::ifstream precipitationsFile(m_precipitationsCsvFilePath.c_str());
		std::string value, firstLine, precipitation;
		int nColumns = 0;

		// Gets number of columns
		// First line is header consisting of string "Time" and Station Ids
		if(getline(precipitationsFile, firstLine))
		{
			std::istringstream iss(firstLine);
			while(getline(iss, value, ';'))
			{
				if(value != "Time")
				{
					m_precStationIds.push_back(atoi(value.c_str()));
				}
				nColumns++;
			}
		}

		// Parse data
		while(getline(precipitationsFile, value))
		{
			std::istringstream iss(value);
			std::pair<time_t ,std::vector<double>> tmp;
			for(int i=0; i<nColumns; i++)
			{
				getline(iss, precipitation, ';');
				if(i==0) // First token is Time, not station id
				{
					tm dateTime;
					int day, month, year, hour, minute;
					//14.5.2010 6:00
					SSCANF(precipitation.c_str(), "%d-%d-%d %d:%d", &year, &month, &day, &hour, &minute);

					dateTime.tm_year = year - 1900;
					dateTime.tm_mon = month - 1;
					dateTime.tm_mday = day;
					dateTime.tm_hour = hour;
					dateTime.tm_min = minute;
					dateTime.tm_sec = 0;

					tmp.first = mktime(&dateTime);
				}
				else
				{
					tmp.second.push_back(std::atof(precipitation.c_str()));
				}
			}
			m_precipitations.push_back(tmp);
		}
		precipitationsFile.close();
	}

	void MatData::getMeasuredHydrographs()
	{
		int nH = 0;
		int rows = (int )m_measuredDischargeVolumes.size();

		if(rows > 0)
		{
			// Get number of columns
			nH = (int )m_measuredDischargeVolumes[0].size();
		}

		for(int i = 0; i < nH; i++)
		{
			bool noValue = true;

			for(size_t j = 0; j < m_measuredDischargeVolumes[i].size(); j++ )
			{
				if(m_measuredDischargeVolumes[j][i] > -999)
				{
					noValue = false;
					break;
				}
			}

			if(noValue)	continue;
			
			std::vector<double> tmp;
			for(int j = 0; j < rows; j++)
			{
				if(m_measuredDischargeVolumes[j][i] > -999)
				{
					// int k = j * 60 / m_minuteStep + 1; // Don't understand
					
					tmp.push_back(m_measuredDischargeVolumes[j][i]);
				}
				else
				{
					tmp.push_back(0);
				}
			}	
			m_measuredHydrographsQ.push_back(tmp);

			m_measureStations.push_back(findStation(m_riverStations, m_mdvStationIds[i]));
		}
	}

	std::vector<std::shared_ptr<Station>>& MatData::getSourceStations()
	{
		return MatData::m_sourceStations;
	}

	void MatData::setSourceStations(const std::vector<std::shared_ptr<Station>>& sourceStations)
	{
		m_sourceStations = sourceStations;
	}

	const int& MatData::getBasinId()
	{
		return m_basinId;
	}
	
	std::string MatData::intToString(const int& number)
	{
		std::ostringstream os;
		os << number;
		return os.str();
	}

	std::string MatData::printDateTime(const time_t& dateTime)
	{
		std::ostringstream str;
		struct tm * timeInfo;
		char buffer[6];
		timeInfo = localtime(&dateTime);
		strftime(buffer, 6, "%H:%M", timeInfo);
		std::string month = intToString(timeInfo->tm_mon + 1);
		if(month.length() == 1)
		{
			month = "0" + month;
		}

		str << timeInfo->tm_year + 1900 << "-" << month << "-" << timeInfo->tm_mday << " " 
			<< buffer;

		return str.str();
	}

	std::shared_ptr<Station> MatData::findStation(std::vector<std::shared_ptr<Station>>& stations, int& id)
	{
		std::shared_ptr<Station> station;
		for(size_t i = 0; i < stations.size(); i++)
		{
			if(stations[i].get()->getId() == id)
				station = stations[i];
		}

		return station;
	}

	Options& MatData::getOptions()
	{
		return this->m_options;
	}

	void MatData::setOptions(Options options)
	{
		this->m_options = options;
	}

	std::vector<std::vector<double>> MatData::getMeasuredHydrographsQ()
	{
		return this->m_measuredHydrographsQ;
	}

	std::vector<std::shared_ptr<Station>>& MatData::getMeasureStations()
	{
		return this->m_measureStations;
	}

	std::vector<std::shared_ptr<Station>>& MatData::getRiverStations()
	{
		return this->m_riverStations;
	}

	std::vector<std::shared_ptr<Channel>>& MatData::getChannels()
	{
		return this->m_channels;
	}

	std::shared_ptr<Channel> MatData::getChannelById(int id)
	{
		std::shared_ptr<Channel> res;
		for (size_t i = 0; i < this->m_channels.size(); i++)
		{
			if(m_channels[i]->getId() == id)
			{
				res = m_channels[i];
				break;
			}
		}
		return res;
	}

	std::shared_ptr<Subbasin> MatData::getSubbasinById(int id)
	{
		std::shared_ptr<Subbasin> res;
		for (size_t i = 0; i < this->m_subbasins.size(); i++)
		{
			if(m_subbasins[i]->getId() == id)
			{
				res = m_subbasins[i];
				break;
			}
		}
		return res;
	}

	std::vector<std::shared_ptr<Subbasin>>& MatData::getSubbasins()
	{
		return m_subbasins;
	}
	
	void MatData::clearResults()
	{
		for (size_t i = 0; i < m_channels.size(); i++)
		{
			m_channels[i]->getHydrograph().clear();
		}
	}

	precipitationsVector MatData::getPrecipitations()
	{
		return m_precipitations;
	}

	void MatData::setPrecipitations(const precipitationsVector precipitations)
	{
		m_precipitations = precipitations;
	}

	int MatData::checkFWL()
	{
		int FWL = -1;
		if(!m_simulationDone)
		{
			CLOG(FATAL,"model") << "Cannot check FWL exceeding, simulation is not done!";	
		}
		for (size_t i = 0; i < m_measureStations.size(); i++)
		{
			std::vector<double> qOut = m_channels[m_measureStations[i]->getChannelIndex()]->getHydrograph().getQOut();
			double spa1 = m_measureStations[i]->getSpa1();
			double spa2 = m_measureStations[i]->getSpa2();
			double spa3 = m_measureStations[i]->getSpa3();
			int basinFWL = -1;

			for (size_t j = 0; j < qOut.size(); j++)
			{
				int channelFWL = -1;
				
				if(qOut[j] < spa1)
					channelFWL = 0;
				else if(qOut[j] < spa2)
					channelFWL = 1;
				else if(qOut[j] < spa3)
					channelFWL = 2;
				else
					channelFWL = 3;

				if(channelFWL > basinFWL)
					basinFWL = channelFWL;
			}

			if(basinFWL > FWL)
				FWL = basinFWL;
		}

		return FWL;
	}
}