-
Reineking Bjoern authored818b3175
/*========================================================*/
/* Seeds Dispersal Models Class */
/* - severals ways to dispers are implemented - */
/*========================================================*/
#include "Disp.h"
#include <vector>
#include <numeric>
#include <cmath>
#include <math.h>
#include <cstdlib>
#include <ctime>
#include <chrono>
#include <boost/random.hpp>
#include <boost/random/uniform_int.hpp>
using namespace std;
// boost::mt19937 est un Mersenne twister generator, ou générateur de nombres pseudo-aléatoires
// boost::uniform_01<RandomGenerator> génère une distribution aléatoire uniforme
// boost::variate_generator<RandomGenerator&, Normal> est un bivariate generator (générateur MT19937+distribution standard uniforme)
typedef boost::mt19937 RandomGenerator;
typedef boost::uniform_01<RandomGenerator&> Uni01;
typedef boost::uniform_int<int> UniInt;
typedef boost::variate_generator<RandomGenerator&, UniInt> GeneratorUniInt;
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
/* Constructors */
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
Disp::Disp() : m_FGparams(), m_SeedMapIn(NULL), m_SeedMapOut(NULL),
m_FGdistCircle(0,vector< vector<int> >(0, vector<int>(0))),
m_prop_d1(0,vector<float>(0.0)), m_prop_d2(0,vector<float>(0.0)),
m_prob_d1(0,vector<float>(0.0)), m_prob_d2(0,vector<float>(0.0))
{
/* Nothing to do */
}
Disp::Disp(vecFGPtr FGparams, DoubleMapPtr seedMapIn, DoubleMapPtr seedMapOut, bool doDisp) :
m_FGparams(FGparams), m_SeedMapIn(seedMapIn), m_SeedMapOut(seedMapOut),
m_FGdistCircle(0,vector< vector<int> >(0, vector<int>(0))),
m_prop_d1(0,vector<float>(0.0)), m_prop_d2(0,vector<float>(0.0)),
m_prob_d1(0,vector<float>(0.0)), m_prob_d2(0,vector<float>(0.0))
{
/* Nothing to do */
}
Disp::Disp(vecFGPtr FGparams, DoubleMapPtr seedMapIn, DoubleMapPtr seedMapOut) :
m_FGparams(FGparams), m_SeedMapIn(seedMapIn), m_SeedMapOut(seedMapOut)
{
m_FGdistCircle.resize(m_FGparams->size(),vector< vector<int> >(0));
m_prop_d1.resize(m_FGparams->size(),vector<float>(0.0)); /* proportion of seeds into crown d50 */
m_prop_d2.resize(m_FGparams->size(),vector<float>(0.0)); /* proportion of seeds into crown d99 */
m_prob_d1.resize(m_FGparams->size(),vector<float>(0.0)); /* probability vector of receiving seeds into crown d50 */
m_prob_d2.resize(m_FGparams->size(),vector<float>(0.0)); /* probability vector of receiving seeds into crown d99 */
cout << "> GetDistancesXY & GetPropProb..." << endl;
GetDistancesXY();
GetPropProb();
}
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
/* Destructor */
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
Disp::~Disp()
{
/* Nothing to do */
}
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
/* Getters & Setters */
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
const vector< vector< vector<int> > >& Disp::getFGdistCircle() const{ return m_FGdistCircle; }
const vector< vector<int> >& Disp::getFGdistCircle(int fg) const{ return m_FGdistCircle[fg]; }
const vector<int>& Disp::getFGdistCircle(int fg, int d_xy) const{ return m_FGdistCircle[fg][d_xy]; }
const vector< vector<float> >& Disp::getprop_d1() const{ return m_prop_d1; }
const vector<float>& Disp::getprop_d1(int fg) const{ return m_prop_d1[fg]; }
const vector< vector<float> >& Disp::getprop_d2() const{ return m_prop_d2; }
const vector<float>& Disp::getprop_d2(int fg) const{ return m_prop_d2[fg]; }
const vector< vector<float> >& Disp::getprob_d1() const{ return m_prob_d1; }
const vector<float>& Disp::getprob_d1(int fg) const{ return m_prob_d1[fg]; }
const vector< vector<float> >& Disp::getprob_d2() const{ return m_prob_d2; }
const vector<float>& Disp::getprob_d2(int fg) const{ return m_prob_d2[fg]; }
vecFGPtr Disp::getFGparams_() { return m_FGparams; }
DoubleMapPtr Disp::getSeedMapIn_() { return m_SeedMapIn; }
DoubleMapPtr Disp::getSeedMapOut_() { return m_SeedMapOut; }
void Disp::setFGdistCircle(const vector< vector< vector<int> > >& FGdistCircle){ m_FGdistCircle = FGdistCircle; }
void Disp::setFGdistCircle(const int& fg, const vector< vector<int> >& FGdistCircle){ m_FGdistCircle[fg] = FGdistCircle; }
void Disp::setFGdistCircle(const int& fg, const int& d_xy, const vector<int>& FGdistCircle){ m_FGdistCircle[fg][d_xy] = FGdistCircle; }
void Disp::setprop_d1(const vector< vector<float> >& prop_d1){ m_prop_d1 = prop_d1; }
void Disp::setprop_d1(const int& fg, const vector<float>& prop_d1){ m_prop_d1[fg] = prop_d1; }
void Disp::setprop_d2(const vector< vector<float> >& prop_d2){ m_prop_d2 = prop_d2; }
void Disp::setprop_d2(const int& fg, const vector<float>& prop_d2){ m_prop_d2[fg] = prop_d2; }
void Disp::setprob_d1(const vector< vector<float> >& prob_d1){ m_prob_d1 = prob_d1; }
void Disp::setprob_d1(const int& fg, const vector<float>& prob_d1){ m_prob_d1[fg] = prob_d1; }
void Disp::setprob_d2(const vector< vector<float> >& prob_d2){ m_prob_d2 = prob_d2; }
void Disp::setprob_d2(const int& fg, const vector<float>& prob_d2){ m_prob_d2[fg] = prob_d2; }
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
/* Other functions */
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
/* Exponential law density function */
float Disp::ExpDensity(float x, float lambda)
{
return lambda * exp(-lambda * x);
}
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
float Disp::ContinuousDecrProba(int d, int d99)
{
float res = 1.0 - d / (float) d99;
if (res!=res) { return 0.0; // nan case
} else { return res; }
}
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
/* Several Dispersal functions */
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
int gaussCircleProblem(int radius)
{
int allPoints=0; //holds the sum of points
double y=0; //will hold the precise y coordinate of a point on the circle edge for a given x coordinate.
long inscribedSquare=(long) sqrt(radius*radius/2); //the length of the side of an inscribed square in the upper right quarter of the circle
int x=(int)inscribedSquare; //will hold x coordinate - starts on the edge of the inscribed square
while(x<=radius){
allPoints+=(long) y; //returns floor of y, which is initially 0
x++; //because we need to start behind the inscribed square and move outwards from there
y=sqrt(radius*radius-x*x); // Pythagorean equation - returns how many points there are vertically between the X axis and the edge of the circle for given x
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
}
allPoints*=8; //because we were counting points in the right half of the upper right corner of that circle, so we had just one-eightth
allPoints+=(4*inscribedSquare*inscribedSquare); //how many points there are in the inscribed square
allPoints+=(4*radius+1); //the loop and the inscribed square calculations did not touch the points on the axis and in the center
return allPoints;
}
void Disp::GetDistancesXY()
{
for(unsigned fg=0; fg<m_FGparams->size(); fg++)
{
/* conversion of dispersal distances in meters into pixels */
int d1 = m_FGparams->at(fg).getDisp50() / m_SeedMapIn->getXres(); // disp50 into pixels
int d2 = m_FGparams->at(fg).getDisp99() / m_SeedMapIn->getXres(); // disp99 into pixels
int dld = m_FGparams->at(fg).getDispLD() / m_SeedMapIn->getXres(); // dispLD into pixels
/* Determination of indices of cells relative to seed pool for each crown */
vector<int> v1x, v2x, vldx, v1y, v2y, vldy;
v1x.reserve(gaussCircleProblem(d1)); v2x.reserve(gaussCircleProblem(d2)); vldx.reserve(gaussCircleProblem(dld));
v1y.reserve(gaussCircleProblem(d1)); v2y.reserve(gaussCircleProblem(d2)); vldy.reserve(gaussCircleProblem(dld));
for (int x= - (max(max(d1,d2),dld)); x<=(max(max(d1,d2),dld)); x++)
{
for (int y= - (max(max(d1,d2),dld)); y<=max(max(d1,d2),dld); y++)
{
if(x*x + y*y - d1 * d1 <= 0) /* into the r=d1 disk */
{
v1x.emplace_back(x);
v1y.emplace_back(y);
} else if(x*x + y*y - d2 * d2 <= 0) /* into the r1=d1, r2=d2 crown */
{
v2x.emplace_back(x);
v2y.emplace_back(y);
} else if(x*x + y*y - dld * dld <= 0) /* into the r1=d2, r2=dld crown */
{
vldx.emplace_back(x);
vldy.emplace_back(y);
}
}
}
v1x.shrink_to_fit(); v2x.shrink_to_fit(); vldx.shrink_to_fit();
v1y.shrink_to_fit(); v2y.shrink_to_fit(); vldy.shrink_to_fit();
m_FGdistCircle[fg].clear();
m_FGdistCircle[fg].reserve(6);
m_FGdistCircle[fg].emplace_back(v1x); m_FGdistCircle[fg].emplace_back(v2x); m_FGdistCircle[fg].emplace_back(vldx);
m_FGdistCircle[fg].emplace_back(v1y); m_FGdistCircle[fg].emplace_back(v2y); m_FGdistCircle[fg].emplace_back(vldy);
}
}
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
void Disp::GetPropProb()
{
for(unsigned fg=0; fg<m_FGparams->size(); fg++)
{
/* conversion of dispersal distances in meters into pixels */
int d1 = m_FGparams->at(fg).getDisp50() / m_SeedMapIn->getXres(); // disp50 into pixels
int d2 = m_FGparams->at(fg).getDisp99() / m_SeedMapIn->getXres(); // disp99 into pixels
/* 50 % of seeds are under 0 and ln(2) according to exponential density law (lambda=1) */
/* 0-ln(2) interval is divide into d1+1 parts */
m_prop_d1[fg].clear();
m_prob_d1[fg].clear();
for (int i=0; i<=d1; i++)
{
m_prop_d1[fg].emplace_back( ExpDensity(i*log(2.0)/(d1+1.0),1.0) - ExpDensity((i+1)*log(2.0)/(d1+1.0),1.0));
m_prob_d1[fg].emplace_back(ContinuousDecrProba(i, d2));
}
/* 49% of seeds are dispersed between ln(2) and 6.65*ln(2) (less than 5/10000 seeds are lost) according to exponential density law (lambda=1) */
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
/* ln(2)-6.65*ln(2) interval is divide into d2-d1+1 parts */
m_prop_d2[fg].clear();
m_prob_d2[fg].clear();
for (int i=1; i<=d2-d1; i++)
{
m_prop_d2[fg].emplace_back( ExpDensity(log(2.0) + i*6.65*log(2.0)/(d2-d1+1.0),1.0) - ExpDensity(log(2.0) + (i+1)*6.65*log(2.0)/(d2-d1+1.0),1.0));
m_prob_d2[fg].emplace_back(ContinuousDecrProba(i+d1, d2));
}
}
}
/*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/
// TODO (damien#1#): make possible to consider different x and y resolution
void Disp::DoDispersalPacket(unsigned dispOption, int nbCPUs, vector<unsigned> maskCells)
{
omp_set_num_threads( nbCPUs );
#pragma omp parallel for schedule(dynamic) if(nbCPUs>1)
for (unsigned fg=0; fg<m_FGparams->size(); fg++)
{
/* dispOption = 1 (uniform), 2 (expKernel) or 3 (expKernel + proba) */
//unsigned dispOption = m_FGparams->at(fg).getModeDispersed();
if (dispOption==1 || dispOption==2 || dispOption==3)
{
int xt,yt;
vector<int> v1x_select, v2x_select, v1y_select, v2y_select;
vector<float> prop_d1_select, prop_d2_select;
unsigned nbDrawMax = int (max(1,(int)ceil(m_FGdistCircle[fg][0].size()/2.0)));
unsigned seed = chrono::system_clock::now().time_since_epoch().count();
RandomGenerator rng(seed);
/* conversion of dispersal distances in meters into pixels */
int d1 = m_FGparams->at(fg).getDisp50() / m_SeedMapIn->getXres(); // disp50 into pixels
int d2 = m_FGparams->at(fg).getDisp99() / m_SeedMapIn->getXres(); // disp99 into pixels
int dld = m_FGparams->at(fg).getDispLD() / m_SeedMapIn->getXres(); // dispLD into pixels
if (dispOption==2 || dispOption==3)
{
Uni01 random_01(rng);
/* select cell receiving seeds according to a probability decreasing with distance */
for (int id=0; id<(int)m_FGdistCircle[fg][0].size(); id++)
{
int dist_pt = max(abs(m_FGdistCircle[fg][0][id]), abs(m_FGdistCircle[fg][3][id]));
if (dist_pt < m_prob_d1[fg].size())
{
/* get an random number between 0-1 */
/* compare this number to the probability vector value if < then the cell will receive seeds */
if (dispOption == 2 || (dispOption == 3 && random_01() < m_prob_d1[fg][dist_pt]))
{
v1x_select.push_back(m_FGdistCircle[fg][0][id]);
v1y_select.push_back(m_FGdistCircle[fg][3][id]);
prop_d1_select.push_back(m_prop_d1[fg][dist_pt]);
}
}
}
for (int id=0; id<(int)m_FGdistCircle[fg][1].size(); id++)
{
int dist_pt = max(abs(m_FGdistCircle[fg][1][id]), abs(m_FGdistCircle[fg][4][id])) - d1;
if (dist_pt < m_prob_d2[fg].size())
{
/* get an random number between 0-1 */
/* compare this number to the probability vector value if < then the cell will receive seeds */
if (dispOption == 2 || (dispOption == 3 && random_01() < m_prob_d2[fg][dist_pt]))
{
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
v2x_select.push_back(m_FGdistCircle[fg][1][id]);
v2y_select.push_back(m_FGdistCircle[fg][4][id]);
prop_d2_select.push_back(m_prop_d2[fg][dist_pt]);
}
}
}
}
SpatialMap<double, double>* new_SeedMapOut = new SpatialMap<double, double>(m_SeedMapOut->getCoordinates(fg), m_SeedMapOut->getValues(fg));
/* Do dispersal for each cell containing seeds */
for(vector<unsigned>::iterator cell_ID=maskCells.begin(); cell_ID!=maskCells.end(); ++cell_ID)
{ // loop on pixels
unsigned x = *cell_ID % m_SeedMapIn->getXncell();
unsigned y = trunc(*cell_ID / m_SeedMapIn->getXncell());
if ((*m_SeedMapIn)(x,y,fg) > 0)
{ // test of value of each pixel
/* chose p= 100% of seeds into d1 disk and put 50% / area of disk seeds into */
if (d1==0) /* all seeds fall in the same pixel */
{
(*new_SeedMapOut)(x,y) += ( (*m_SeedMapIn)(x,y,fg) * 0.5 );
} else if (d1>0)
{
if (dispOption==1)
{
v1x_select.clear();
v1y_select.clear();
v1x_select = m_FGdistCircle[fg][0];
v1y_select = m_FGdistCircle[fg][3];
}
for (unsigned id=0; id<v1x_select.size(); id++) /* loop on d1 disk */
{
xt = x + v1x_select[id];
yt = y + v1y_select[id];
if (xt>=0 && yt>=0 && xt<(int)m_SeedMapIn->getXncell() && yt<(int)m_SeedMapIn->getYncell())
{
if (dispOption==1)
{
(*new_SeedMapOut)(xt,yt) += (*m_SeedMapIn)(x,y,fg) * 0.5 / (double) v1x_select.size();
} else if (dispOption==2 || dispOption==3)
{
(*new_SeedMapOut)(xt,yt) += (*m_SeedMapIn)(x,y,fg) * prop_d1_select[id];
}
}
} // end of loop on d1 disk
}
/* chose as many as in first disk and into d1 d2 crow and put 49% / area of crown / p seeds into */
if (d2 == 0)
{
(*m_SeedMapOut)(x,y,fg)+= (*m_SeedMapIn)(x,y,fg) * 0.49;
} else if (d2>0)
{
/* seeds will be dispersed in 2 neighbouring cells min disperse in 1 cell and one of its neighbour */
if (dispOption==1)
{
v2x_select.clear();
v2y_select.clear();
UniInt distrib(0,m_FGdistCircle[fg][1].size());
GeneratorUniInt draw_from_distrib(rng,distrib);
for (unsigned nbDraw = 0; nbDraw < nbDrawMax; nbDraw++)
{ /* Draw of cells into crown that will received seeds */
/*!*/
v2x_select.push_back(m_FGdistCircle[fg][1][draw_from_distrib()]);
v2y_select.push_back(m_FGdistCircle[fg][4][draw_from_distrib()]);
/*!*/
}
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
}
for (int id=0; id<(int)v2x_select.size(); id++)
{
xt = x + v2x_select[id];
yt = y + v2y_select[id];
if (xt>=0 && yt>=0 && xt<(int)m_SeedMapIn->getXncell() && yt<(int)m_SeedMapIn->getYncell())
{
/* First cell selected */
if (dispOption==1)
{
(*new_SeedMapOut)(xt,yt) += (*m_SeedMapIn)(x,y,fg) * 0.49 / (nbDrawMax * 2.0);
/* x of its neighbour */
UniInt distrib(0,3);
GeneratorUniInt draw_from_distrib(rng,distrib);
switch(draw_from_distrib())
{
case 0 : xt++;
break;
case 1 : xt--;
break;
case 2 : yt++;
break;
case 3 : yt--;
break;
}
if (xt>=0 && yt>=0 && xt<(int)m_SeedMapIn->getXncell() && yt<(int)m_SeedMapIn->getYncell())
{
(*new_SeedMapOut)(xt,yt) += (*m_SeedMapIn)(x,y,fg) * 0.49 / (nbDrawMax * 2.0);
}
} else if (dispOption==2 || dispOption==3)
{
(*new_SeedMapOut)(xt,yt) +=(*m_SeedMapIn)(x,y,fg) * prop_d2_select[id];
}
}
}
} // end of d50 -> d99 crown dispersal
/* chose 1 cells into d1 d2 crow and put 1% / area of crown / p seeds into */
if (dld == 0)
{
(*m_SeedMapOut)(x,y,fg) += (*m_SeedMapIn)(x,y,fg) * 0.01;
} else if (dld>0)
{
if((int)m_FGdistCircle[fg][2].size()>0)
{
UniInt distrib(0,(int)m_FGdistCircle[fg][2].size() - 1);
GeneratorUniInt draw_from_distrib(rng,distrib);
/*!*/
int LD_draw = draw_from_distrib(); //rand()% vSize;
/*!*/
xt = x + m_FGdistCircle[fg][2][LD_draw];
yt = y + m_FGdistCircle[fg][5][LD_draw];
if (xt>=0 && yt>=0 && xt<(int)m_SeedMapIn->getXncell() && yt<(int)m_SeedMapIn->getYncell())
{
(*new_SeedMapOut)(xt,yt) += (*m_SeedMapIn)(x,y,fg) * 0.01;
}
}
} // end of d99 -> ldd crown dispersal
} // end test if some seed to disperse
} // end loop over pixels
m_SeedMapOut->setValues(fg, new_SeedMapOut->getValues());
} else
{
cout << "No seed dispersal for the PFG "<< fg << "!" << endl;
}
} // end loop over PFGs
}
421422