Skip to content
Snippets Groups Projects
PrecipitationUncertainity.cpp 4.04 KiB
Newer Older
  • Learn to ignore specific revisions
  • Radim Vavřík's avatar
    Radim Vavřík committed
    #include "PrecipitationUncertainity.h"
    #include "UncertainityOptions.h"
    #include "easylogging++.h"
    #include <iostream>
    #include <mpi.h>
    
    
    
    namespace math1d_cl
    {
    	PrecipitationUncertainity::PrecipitationUncertainity(std::shared_ptr<AbstractRandom> random, UncertainityOptions options, std::shared_ptr<MatData> matData) : AbstractParam(random, options, matData)
    	{
    		m_randomizedValues.reserve(m_options.simulationCount); // Allocate space for randomized values
    		m_originalValues = matData->getPrecipitations(); // Get original precipitations
    	}
    
    	void PrecipitationUncertainity::generateValues()
    	{
    		size_t timeSteps = m_originalValues.size(); // Values per station
    		size_t stationsCount = m_originalValues.front().second.size(); // Weather stations count
    
    		int rank = -1, numProc = -1; // Holds rank of current process and total number of processes
     		MPI_Comm_size(MPI_COMM_WORLD, &numProc);
        	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
        	
        	// Determine chunks of precipitaion values
       	 	std::vector<int> precipChunkSizes; // Holds numbers of values in each chunk
       	 	std::vector<int> displValues; // Holds displacement values
        	for (size_t p = 0; p < m_options.chunkSizes.size(); ++p)
        	{
        		precipChunkSizes.push_back(m_options.chunkSizes[p] * timeSteps * stationsCount);
        		displValues.push_back(p * m_options.chunkSizes[p] * timeSteps * stationsCount);
        	}
    
    		std::vector<double> allRandomizedValues; // Holds precipitation values for all iterations
    		std::vector<double> recvRandomizedValues; // Holds precipitation values for received chunk
    
       		/* Root rank generate chunks of values */
    		if(rank == 0)
    		{
    			/* Generate randomized values */
    			CLOG(INFO,"montecarlo") << "Randomizing precipitations...";
    			
    			// Randomize them for each monte carlo run
    			allRandomizedValues.reserve(m_options.simulationCount * timeSteps * stationsCount); // Reserve space for values
    
    			for(size_t mc = 0; mc < (size_t)m_options.simulationCount; ++mc)
    			{
    				precipitationsVector mcPrec = m_originalValues; // Copy precipitations
    				m_random->fillRandom(mcPrec, m_options.precipDeviation); // Randomize precipitations
    				
    				// For each timestep
    				// Append precipitations to buffer 
    				for (size_t i = 0; i < mcPrec.size(); ++i)
    				{
    					allRandomizedValues.insert(allRandomizedValues.end(), mcPrec[i].second.begin(), mcPrec[i].second.end());
    				}
    				mcPrec.clear();
    			}
    			CLOG(INFO,"montecarlo") << "DONE.";
    		}
    
    
    		/* MPI Collective scatter */
    		recvRandomizedValues.resize(precipChunkSizes[rank]); // Reserve space for received chunk
    		MPI_Scatterv(&allRandomizedValues.front(),&precipChunkSizes.front(), &displValues.front(), MPI_DOUBLE, 
    			&recvRandomizedValues.front(), precipChunkSizes[rank], MPI_DOUBLE, 0, MPI_COMM_WORLD);
    
    		// All randomized values can be now freed
    		if(rank == 0) allRandomizedValues.clear();
    
    		CLOG(INFO,"montecarlo") << "[MPI Rank " << rank <<  "]: Received " << recvRandomizedValues.size() << " values";
    		
    		
    
    		// Copy the received values to member vector
    		for (int i = 0; i < m_options.chunkSizes[rank]; ++i)
    		{
    			precipitationsVector iterationBuffer;
    			for (size_t tim = 0; tim < timeSteps; ++tim)
    			{
    			 	std::pair<time_t, std::vector<double>> timestepValues;
    			 	// Copy time value
    				timestepValues.first = m_originalValues[tim].first;
    				// Copy precipitaions for each station
    				std::vector<double>::iterator start = recvRandomizedValues.begin() + (i * timeSteps * stationsCount) + (tim * stationsCount);
    
    				timestepValues.second.insert(timestepValues.second.begin(), start, start + stationsCount);
    				iterationBuffer.push_back(timestepValues);
    			}
    			m_randomizedValues.push_back(iterationBuffer);
    			iterationBuffer.clear();
    		}
    		recvRandomizedValues.clear();
    		if(rank == 0) CLOG(INFO,"montecarlo") << "Done.";	
    	}
    
    	void PrecipitationUncertainity::setParam(MatData &matData)
    	{
    		if(m_currentIterationNumber >= (size_t)m_options.simulationCount) 
    		{
    			std::cerr << "ERROR: Iteration count bigger than monte carlo loop count!";
    			std::exit(-1);
    		}
    		matData.setPrecipitations(m_randomizedValues[m_currentIterationNumber++]);
    	}
    }