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