ManningUncertainity.cpp

#include "ManningUncertainity.h"
#include "Channel.h"
#include "Subbasin.h"
#include "easylogging++.h"
#include <iostream>
#include <mpi.h>

namespace math1d_cl
{
	ManningUncertainity::ManningUncertainity(std::shared_ptr<AbstractRandom> random, UncertainityOptions options,std::shared_ptr<MatData> matData) : AbstractParam(random, options, matData)
	{
		m_channelValues.reserve(m_options.simulationCount);
		m_subbasinValues.reserve(m_options.simulationCount);

		/* Get original values */
		std::vector<double> channelOrig, subbasinOrig;

		// Get channel's N
		//double chN = 0.0, subN = 0.0;
		for(std::vector<std::shared_ptr<Channel>>::const_iterator it = matData->getChannels().begin(); it != matData->getChannels().end(); ++it)
		{
			m_channelOriginalValues.push_back((*it)->getN());
			//m_channelOriginalValues.push_back(chN++);
		}

		// Get subbasin's N
		for(std::vector<std::shared_ptr<Subbasin>>::const_iterator it = matData->getSubbasins().begin(); it != matData->getSubbasins().end(); ++it)
		{
			m_subbasinOriginalValues.push_back((*it)->getN());
			//m_subbasinOriginalValues.push_back(subN++);
		}
	}

	void ManningUncertainity::setParam(MatData &matData)
	{
		/*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);*/

		if((m_currentIterationNumber > m_channelValues.size()) || (m_currentIterationNumber > m_subbasinValues.size()))
		{
			std::cerr << "ERROR: Manning uncertainity: Current iteration number bigger than generated values.";
			std::exit(-1);
		} else 
		{
			for(size_t i = 0; i < matData.getChannels().size(); ++i)
			{
				//CLOG(INFO,"montecarlo") << "[Rank " << rank <<"]N Channel - Orig: " << m_channelOriginalValues[i] << " New: " << m_channelValues[m_currentIterationNumber][i];
				matData.getChannels()[i]->setN(m_channelValues[m_currentIterationNumber][i]);
			}

			for(size_t i = 0; i < matData.getSubbasins().size(); ++i)
			{
				//CLOG(INFO,"montecarlo") << "[Rank " << rank <<"]N Subbasin - Orig: " << m_subbasinOriginalValues[i] << " New: " << m_subbasinValues[m_currentIterationNumber][i];
				matData.getSubbasins()[i]->setN(m_subbasinValues[m_currentIterationNumber][i]);
			}
		}	

		/* Do not forget to increment iteration number after each setParam call !*/
		m_currentIterationNumber++;
	}
	
	void ManningUncertainity::generateValues()
	{
		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);

    	const size_t chunkSize = m_options.simulationCount / numProc; // Monte Carlo iterations per chunk
    	const size_t valuesPerChunk = (m_subbasinOriginalValues.size() + m_channelOriginalValues.size()) * chunkSize; // Values per chunk

		// One chunk of values for one process
		// Allocate array for one chunk
		double *randomizedValuesChunk = new double[valuesPerChunk]; 

		// Allocate array for each process
		//m_subbasinValuesChunk = new double[m_subbasinOriginalValues.size() * chunkSize];

		int slaveRank = 1;

		if(rank == 0)
		{
			/* Generate randomized values */
			/*
					Data are origanised as follows:

							CH0: N
							CH1: N 
							.
							.
							.
							CH48: N
							CH49: N
				Iteration 0	-------
							SUB0: N
							SUB1: N
							.
							.
							.
							SUB48: N
							SUB49: N
							------------------------
							CH0: N
							CH1: N 
							.
							.
							.
							CH48: N
							CH49: N
				Iteration 1	-------
							SUB0: N
							SUB1: N
							.
							.
							.
							SUB48: N
							SUB49: N
					etc...
			*/
			CLOG(INFO,"montecarlo") << "Randomizing Manning's coefficients...";
			CLOG(INFO,"montecarlo") << "[MPI]Total chunk count: " << m_options.simulationCount/chunkSize;
			CLOG(INFO,"montecarlo") << "[MPI]Chunksize: " << chunkSize;
		
			for(size_t mc = 0; mc < (size_t)m_options.simulationCount; ++mc)
			{
				double min = 0.0, max = 0.0;
				std::vector<double> randomChannel, randomSubbasin;

				// Channels
				for(size_t i = 0; i < m_channelOriginalValues.size(); ++i)
				{
					int chunkNum = i + ((mc % chunkSize) *(m_subbasinOriginalValues.size() + m_channelOriginalValues.size())); // Compute index
					//CLOG(INFO,"montecarlo") << "[MPI]ChunkNum: " << chunkNum;


					if(m_channelOriginalValues[i] >= 0 && m_options.nDeviation != 0) // Skip not defined values (-999.999)
					{
						switch(m_options.nRandType) {
							case RandomType::KERNEL_DENSITY :	
								std::cerr << "ERROR: N cannot be modeled by kernel density.";
								std::exit(-1);
								break;
							case RandomType::UNIFORM :
								// Lower and upper boundary of selected deviation
								// Mind limits given in options
								min = keepLimits(m_channelOriginalValues[i] - (m_channelOriginalValues[i] * m_options.nDeviation), m_options.nLimits);
								max = keepLimits(m_channelOriginalValues[i] + (m_channelOriginalValues[i] * m_options.nDeviation), m_options.nLimits);

								// Sample random value within computed boundaries
								//randomchannel.push_back(m_random->getRandDouble(min, max));
								randomizedValuesChunk[chunkNum] = m_random->getRandDouble(min, max);
								break;
							case RandomType::NORMAL :
								{
									double mean = 0.0;
									double stddev = (m_channelOriginalValues[i] * 0.02) / 3; 
									randomizedValuesChunk[chunkNum] = keepLimits(m_channelOriginalValues[i] + m_random->getRandDouble(mean, stddev), m_options.nLimits);
								}
								break;
							default:
								std::cerr << "ERROR: Bad random function, check random function options.";
								std::exit(-1);
								break;
						}
					
					} else
					{
						randomizedValuesChunk[chunkNum] = m_channelOriginalValues[i];
					}	
				}
				//m_channelValues.push_back(randomChannel);
				//randomChannel.clear();

				// Subbasins
				for(size_t i = 0; i < m_subbasinOriginalValues.size(); ++i)
				{
					int chunkNum = i + m_channelOriginalValues.size() + ((mc % chunkSize) * (m_subbasinOriginalValues.size() + m_channelOriginalValues.size())); // Compute index
					
					if(m_subbasinOriginalValues[i] >= 0 && m_options.nDeviation != 0) // Skip not defined values (-999.999)
					{
						switch(m_options.nRandType) {
							case RandomType::KERNEL_DENSITY :	
								std::cerr << "ERROR: N cannot be modeled by kernel density.";
								std::exit(-1);
								break;
							case RandomType::UNIFORM :
								// Lower and upper boundary of selected deviation
								// Mind limits given in options
								min = keepLimits(m_subbasinOriginalValues[i] - (m_subbasinOriginalValues[i] * m_options.nDeviation), m_options.nLimits);
								max = keepLimits(m_subbasinOriginalValues[i] + (m_subbasinOriginalValues[i] * m_options.nDeviation), m_options.nLimits);

								// Sample random value within computed boundaries
								//randomsubbasin.push_back(m_random->getRandDouble(min, max));
								randomizedValuesChunk[chunkNum] = m_random->getRandDouble(min, max);
								break;
							case RandomType::NORMAL :
								{
									double mean = 0.0;
									double stddev = (m_subbasinOriginalValues[i] * 1.5) / 3; 
									randomizedValuesChunk[chunkNum] = keepLimits(m_subbasinOriginalValues[i] + m_random->getRandDouble(mean, stddev), m_options.nLimits);
								}
								break;
							default:
								std::cerr << "ERROR: Bad random function, check random function options.";
								std::exit(-1);
								break;
						}
					} else
					{
						randomizedValuesChunk[chunkNum] = m_subbasinOriginalValues[i];
						//randomSubbasin.push_back(m_subbasinOriginalValues[i]);
					}	

				}
				//m_subbasinValues.push_back(randomSubbasin);
				//randomSubbasin.clear();
			

				if(((mc+1) % chunkSize) == 0) // One chunk generated, send it to its process
				{
					if((mc+1) == chunkSize) // Keep first chunk to root process
					{
						//std::copy(&randomizedValuesChunk[0], &randomizedValuesChunk[(chunkSize-1)], std::back_inserter(m_subbasinValues));
						//memcpy(randomizedValuesChunk, m_subbasinValuesChunk, chunkSize);

						for(size_t iter = 0; iter < chunkSize; ++iter)
						{
							std::vector<double> iterationBuffer;
							for(size_t ch = 0; ch < m_channelOriginalValues.size(); ++ch)
							{
								iterationBuffer.push_back(randomizedValuesChunk[ch + (iter * (m_subbasinOriginalValues.size() + m_channelOriginalValues.size()))]);	
							}
							m_channelValues.push_back(iterationBuffer);
							iterationBuffer.clear();

							for(size_t sub = 0; sub < m_subbasinOriginalValues.size(); ++sub)
							{
								iterationBuffer.push_back(randomizedValuesChunk[sub + m_channelOriginalValues.size() + (iter * (m_subbasinOriginalValues.size() + m_channelOriginalValues.size()))]);	
							}
							m_subbasinValues.push_back(iterationBuffer);
							iterationBuffer.clear();
						}

						CLOG(INFO,"montecarlo") << "[MPI]Root rank: Got " << m_subbasinValues.size() << " iterations in chunk";

					} else { 
						// Send chunk to its respective process
						//int slaveRank = (int) ((mc+1) / chunkSize); // Determine destination chunk
						CLOG(INFO,"montecarlo") << "[MPI]Sent chunk no." << mc/chunkSize << " to slave rank:  " << slaveRank;
						MPI_Send(randomizedValuesChunk, valuesPerChunk, MPI_DOUBLE, slaveRank, m_tag, MPI_COMM_WORLD);
						
						// Increment slave rank to next process
						slaveRank++;
					}
				}
			}
			CLOG(INFO,"montecarlo") << "DONE.";
		}

		/* Receive values from slave */
		if(rank > 0)
		{
			// Receive values from root
			MPI_Status mpiStat;
			MPI_Recv(randomizedValuesChunk, valuesPerChunk, MPI_DOUBLE , 0, m_tag, MPI_COMM_WORLD, &mpiStat);
			CLOG(INFO,"montecarlo") << "[MPI]Slave rank " << rank << ": Received chunk from " << mpiStat.MPI_SOURCE  << slaveRank;
			
			// Copy received values to memeber vector
			std::vector<double> iterationBuffer;
			
			// Channels
			for(size_t iter = 0; iter < chunkSize; ++iter)
			{
				for(size_t ch = 0; ch < m_channelOriginalValues.size(); ++ch)
				{
					iterationBuffer.push_back(randomizedValuesChunk[ch + (iter * (m_subbasinOriginalValues.size() + m_channelOriginalValues.size()))]);	
				}
				m_channelValues.push_back(iterationBuffer);
				iterationBuffer.clear();

				for(size_t sub = 0; sub < m_subbasinOriginalValues.size(); ++sub)
				{
					iterationBuffer.push_back(randomizedValuesChunk[sub + m_channelOriginalValues.size() + (iter * (m_subbasinOriginalValues.size() + m_channelOriginalValues.size()))]);	
				}
				m_subbasinValues.push_back(iterationBuffer);
				iterationBuffer.clear();
			}
			
			iterationBuffer.clear();

			// Print first element of each iteration in received chunk
			/*for(std::vector<std::vector<double>>::const_iterator it = m_subbasinValues.begin(); it != m_subbasinValues.end(); ++it)
			{
				CLOG(INFO,"montecarlo") << (*it).at(0);
			}*/

			delete[] randomizedValuesChunk;
		}
		


		
	}
}