Skip to content
Snippets Groups Projects
KernelDensity.cpp 4.86 KiB
Newer Older
  • Learn to ignore specific revisions
  • Radim Vavřík's avatar
    Radim Vavřík committed
    #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;
    	}
    
    
    	
    
    }