#include "KernelDensity.h"
#include "CSVReader.h"
#include "easylogging++.h"
#include <string>
#include <iostream>
#include <random>
#include <omp.h>
#include <mpi.h>

namespace math1d_cl {

	KernelDensity::KernelDensity(UncertainityOptions options) : AbstractRandom(options)
	{	
		
		int rank = -1;
		MPI_Comm_rank(MPI_COMM_WORLD, &rank);
		if(rank == 0) // Only root needs to load CDF
		{
			// Load probability density curves for each limit
			CLOG(INFO,"montecarlo") << "[KernelDensity] Loading CDF...";
			m_kernelEst.resize(m_options.limits.size());

			for (size_t lim = 0; lim < m_options.limits.size(); ++lim)
			{
				// Open CDF CSV file
				CLOG(INFO,"montecarlo") << "[KernelDensity] Limit " << m_options.limits[lim].upperLimit << " File: " << m_options.limits[lim].pathCDF;
				math1d_cl::CSVReader reader('.' /*decimal*/,','/*separator*/);
				reader.open(m_options.limits[lim].pathCDF);

				// Read file
				bool hasRows = false;
				while(!reader.eof())
				{
					std::shared_ptr<std::vector<std::string>> row = reader.getRow();
					if(row->size() == 5) // Check for expected row size
					{
						m_kernelEst[lim].insert( std::pair<double, double>(std::stod(row->at(2).c_str()),std::stod(row->at(1).c_str())));
						hasRows = true;
					}
					row.reset();
				}
				reader.close();
				if(!hasRows)
				{
					std::cerr << "[KernelDensity] ERROR: No points loaded for limit " << m_options.limits[lim].upperLimit;
					std::exit(-1);
				}

			}
			CLOG(INFO,"montecarlo") << "[KernelDensity] Loading CDF Done.";
		}
	}

	void KernelDensity::fillRandom(math1d_cl::precipitationsVector &prec, double deviation)
	{
		int rank = -1;
		MPI_Comm_rank(MPI_COMM_WORLD, &rank);
		if(rank != 0)
		{
			std::cerr << "[KernelDensity][MPI Rank "<< rank << "] Only root can generate values !";
			std::exit(-1);
		}

		// Compute start of Alladin forecast
		// Holds number of offsets from the beginning of prediction
		//std::normal_distribution<double> normal(0.48,1.3);
		//std::normal_distribution<double> normal(0.197,1.0);
		size_t i = 0, j = 0;
		#pragma omp parallel for schedule(dynamic) private(i,j)
		for (i = (m_options.daysMeasured == 0) ? 0 : (m_options.daysMeasured * m_options.valuesPerDay) + 1; i < prec.size(); ++i) // TODO: To config.xml as Measured Days, Predictions per Day
		{
			//int offset = (int)floor(cnt/6);

			for (j = 0; j < prec[i].second.size() ; ++j)
			{
				if(prec[i].second[j] >= 0.1)
				{
					//printf("[Time %d Channel %d] Orig: %.4f -- ",i,j,prec[i].second[j]);
					// Compute bounds given by deviation
					double upper = prec[i].second[j] + (prec[i].second[j]*deviation);
					double lower = prec[i].second[j] - (prec[i].second[j]*deviation);

					// Generate random precipitation within bounds
					double randomized = getRandDouble(lower, upper);
					
					//double error = normal(m_uRand);
		 			//double error = 0;

		 			// Get index for limit
		 			// If the value is out of limits, error is zero
		 			int limitIdx = getLimitIdx(prec[i].second[j]);
		 			double error = (limitIdx < 0) ? 0 : getRandom(limitIdx);

		 			//if(error > 0.3) error = 0.3;

	 				prec[i].second[j] = randomized + error;
	 				//prec[i].second[j] = 6.197;
					
					// Correct if precipitation is negative
					prec[i].second[j] = (prec[i].second[j] < 0) ? 0 : prec[i].second[j];
		 			//printf("New: %.4f Randomized: %.4f Error: %.4f Low: %.4f Up: %.4f\n",prec[i].second[j],randomized,error, lower, upper);
		 			//printf("%.4f\n",prec[i].second[j]);

	 				//printf("New: %.4f -- Error: %.4f Low: %.4f Up: %.4f\n",prec[i].second[j],error, lower, upper);
					
				}
			} // Precipitations
		} // Times
	}

	std::vector<double> KernelDensity::test(int num, int limitIdx)
	{
		std::vector<double> result(num);
		
		#pragma omp parallel for schedule(dynamic)
		for(int i = 0; i < num; ++i)
		{
			//CLOG(INFO,"montecarlo") << "Thread: " << omp_get_thread_num();
			result[i] = getRandom(limitIdx);
			//printf("[prob: %f delta: %f] [y0: %f] [x0: %f] [y1: %f] [x1: %f]\n",y,x,y0,x0,y1,x1);
		}
	
		return result;
	}

	int KernelDensity::getLimitIdx(double precipValue)
	{
		int idx = -1;

		for(size_t lim = 0; lim < m_options.limits.size(); ++lim) {
			if(precipValue < m_options.limits[lim].upperLimit)
			{
				idx = (int)lim;
				break;
			}
		}
		return idx;
	}

	double KernelDensity::getRandom(int limitIdx) {
		// Generate random probability
		//double y = m_randProb(m_uRand);
		double y = std::generate_canonical<double,std::numeric_limits<double>::digits>(m_uRand);
	
		// Interpolate corresponding error
		// TODO: Add offsets, limits, also here...
		std::map<double, double>::iterator upper = m_kernelEst[limitIdx].upper_bound(y);
		std::map<double, double>::iterator lower = (--m_kernelEst[limitIdx].lower_bound(y));

		double x1 = upper->second;
		double y1 = upper->first;
		double x0 = lower->second;
		double y0 = lower->first;

		// Return the linear interpolation 
		return ( ((y-y0)*(x1-x0))/(y1-y0))+x0;
	}


	

}