#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; } }