diff --git a/.gitignore b/.gitignore index 9aba099..aa558ed 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ *.o html/ TWoTests/ -.Rhistory \ No newline at end of file +.Rhistory +.Rproj.user diff --git a/New Functions.R b/New Functions.R new file mode 100644 index 0000000..5e7759f --- /dev/null +++ b/New Functions.R @@ -0,0 +1,86 @@ +indspec<-function(N, TNW , BIC , Method = 2 ){ + + + + WIC<- TNW - BIC + + + ind<- matrix(nrow = N, ncol = 2) + colnames(ind) <- c("genotype_Mean", "width_sd") + + + + x<- runif(1,0,1-TNW) + + + #Random + + if(Method == 0 || Method == "random" || Method == "Random"){ + + for (i in 1:N) { + + for (i in 1:N) { + ind[i,1]<- runif(1,(x+WIC/2),(x+TNW-WIC/2)) + } + } + + } + + + #Sistematical + + if(Method == 1 || Method =="sistematical" || Method == "Sistematical"){ + + ind[,1] <-seq(from = x+WIC/2, to= x+TNW-WIC/2, length.out = N) + + } + + #Normal + + if(Method == 2 || Method == "Normal" || Method == "Normal"){ + + for (i in 1:N) { + + ind[i]<- rnorm(1,(x+TNW/2),BIC/2) + while (ind[i,1]x+TNW-WIC/2) { + ind[i]<- rnorm(1,(x+TNW/2),BIC/2) + } + } + } + + ind[,2]<-WIC/2 + + return(ind) +} + +show_ind_fitness_curves<-function(nInd,TNW,BIC, Method=2){ + + + inds<- indspec(nInd,TNW,BIC, Method) + + x<- seq(from= 0, to= 1, length.out = 1000) # Generates the z_scores + y = matrix(nrow=1000, ncol=nInd) # creates a matrix for storing the Distributions + + for (i in 1:nrow(inds)) { + + y[,i] = dnorm(x,inds[i,1], inds[i,2]) + } + + Y=apply(y, 1, sum) + + plot(Y, col="darkblue", type = "line") + lines(nInd*dnorm(x,mean = mean(inds[,1]),TNW/2),col="grey") # Plots the Distribution of individual optimum trait values + lines(Y/nInd, type = "line", col="purple") + for (i in 1:nrow(inds)) { + lines(y[,i]) + } + + info<-list() + + info[[1]]<-inds + info[[2]]<-Y/nInd + + names(info)<-c("inds","pop") + + return(info) +} diff --git a/TWoLife.R b/TWoLife.R index 9af9b66..358ce87 100644 --- a/TWoLife.R +++ b/TWoLife.R @@ -1,3 +1,6 @@ + +setwd("~/Desktop/TWoLife") + # Função "wrapper" para a chamada em C: # # Os passos abaixo foram adaptados de http://users.stat.umn.edu/~geyer/rc/ @@ -12,162 +15,256 @@ dyn.load("TWoLife.so") ## carrega os source resultantes como biblioteca dinamica # numb.cells represents both the lenght AND width of the landscape, so numb.cells=100 creates a 100x100 landscape # Land.shape can be 0 = XXX or 1 = XXX. # Bound.condition can be 0 = XXX or 1 = XXX. -Landscape <- function (numb.cells = 100, cell.size = 1, land.shape = 1, type=c("random","blob"), bound.condition=0, cover=1) { - type=match.arg(type) - if(cover < 0 || cover > 1) { - stop("Error creating landscape. Cover must be between 0 and 1") - } - # scape represents the actual landscape - scape <- rep(1, numb.cells*numb.cells) - if(cover < 1) { - NtoRemove=round((1-cover)*numb.cells*numb.cells); - if(type=="random") { - while(NtoRemove>0) - { - i=round(runif(1,0,numb.cells-1)); - j=round(runif(1,0,numb.cells-1)); - # tests to see if this point has already been removed - if(scape[1+numb.cells*j+i] == 1) { - NtoRemove = NtoRemove - 1 - scape[1+numb.cells*j+i] = 0 - } - } - } - if(type=="blob") { - i=round(runif(1,0,numb.cells-1)); - j=round(runif(1,0,numb.cells-1)); - while(NtoRemove>0) - { - # tests to see if this point has already been removed - if(scape[1+numb.cells*j+i] == 1) { - NtoRemove = NtoRemove - 1 - scape[1+numb.cells*j+i] = 0 - } - # Draft a new point to be removed (random walk!) - if(sample(1:2,1) == 1) { - i = i + (-1)**sample(1:2,1) - } else { - j = j + (-1)**sample(1:2,1) - } - if(i == -1) { i=numb.cells-1} - if(i == numb.cells) { i=1} - if(j == -1) { j=numb.cells-1} - if(j == numb.cells) { j=1} - } - } - } - land <- list(numb.cells = numb.cells, cell.size=cell.size, land.shape=land.shape, type=type, bound.condition=bound.condition, cover=cover, scape=scape) - class(land) <- "landscape" - return(land) + +indspec<-function(N, TNW , BIC , Method = 2 ){ + + + + WIC<- TNW - BIC + + + ind<- matrix(nrow = N, ncol = 2) + colnames(ind) <- c("genotype_Mean", "width_sd") + + + #Random + + if(Method == 0 || Method == "random" || Method == "Random"){ + + for (i in 1:N) { + + for (i in 1:N) { + ind[i,1]<- runif(1,(1-(TNW-(WIC/2))),(1-(WIC/2))) + } + } + + } + + + #Sistematical + + if(Method == 1 || Method =="sistematical" || Method == "Sistematical"){ + + ind[,1] <-seq(from = (1-(TNW-(WIC/2))), to= (1-(WIC/2)), length.out = N) + + } + + #Normal + + if(Method == 2 || Method == "Normal" || Method == "Normal"){ + + for (i in 1:N) { + + ind[i]<- rnorm(1,(1-(TNW/2)),BIC/2) + while (ind[i,1]<(1-(TNW-(WIC/2))) || ind[i,1]>(1-(WIC/2))) { + ind[i]<- rnorm(1,(1-(TNW/2)),BIC/2) + } + } + } + + ind[,2]<-WIC/2 + + return(ind) } -TWoLife <- function ( - raio=0.1, - N=80, - AngVis=360, - passo=5, - move=0.5, - taxa.basal=0.6, - taxa.morte=0.1, - incl.birth=0.5/0.01, - incl.death=0, - density.type=0, - death.mat=7, - landscape, - tempo=20, - ini.config=0, - out.code=1) -{ - if(class(landscape) != "landscape") { - stop("Error in function TWoLife: you must provide a valid landscape. See ?Landscape") - } - if(raio>landscape$numb.cells*landscape$cell.size/2) - {stop("Error in function TWoLife: the radius must be lower than or equal to the half of landscape side (radius <= numb.cells*cell.size/2)")} - - saida.C <- .C("TWoLife", - as.double(raio),# 1 - as.integer(N),# 2 - as.double(AngVis),# 3 - as.double(passo),# 4 - as.double(move),# 5 - as.double(taxa.basal),# 6 - as.double(taxa.morte),# 7 - as.double(incl.birth),# 8 - as.double(incl.death),# 9 - as.integer(landscape$numb.cells),# 10 - as.double(landscape$cell.size),# 11 - as.integer(landscape$land.shape),# 12 - as.integer(density.type),# 13 - as.double(death.mat), # 14 - as.integer(ini.config), #15 - as.integer(landscape$bound.condition), #16 - as.integer(landscape$scape), #17 - as.double(tempo), #18 - as.integer(0), # 19 - as.double(rep(0, 5000)), # 20 - as.double(rep(0,5000)), # 21 - as.integer(out.code) - ## verificar se precisa definir o tamanho e se isto nao dará problemas (dois ultimos argumentos) - ) - n <- saida.C[[19]] - x <- saida.C[[20]] - y <- saida.C[[21]] - x <- x[1:n]; y <- y[1:n] - return(data.frame(x=x,y=y)) +frag<- function(directory,TNW =1 , numb.cells = 100, cell.size = 1, land.shape = 1, type=c("random","blob"), bound.condition=0, cover=1){ + + library(vegan) + + land<- Landscape(numb.cells, cell.size, land.shape, type, bound.condition, cover) + temp<- numeric(length(land$scape)) + + setwd(directory) + files<-list.files() + + for (i in 1:length(files)) { + + frag=read.table(files[i]) + + count=0 + for (j in 1:sqrt(length(land$scape))) { + for (k in 1:sqrt(length(land$scape))) { + count = count +1 + temp[count]= temp[count]+frag[j,k] + } + } + } + + #land$scape<-decostand(temp, method = "range") + land$scape<-1-(decostand(temp, method = "range")*(TNW)) + + return (land) } -# ## Um teste rapido -# land <- Landscape(cover=1,type="b",cell.size=100) -# # ## Uma rodada: coordenadas dos sobreviventes apos t=20 -# teste <- TWoLife(raio=1560, -# N=10, -# AngVis=360, -# passo=10, -# move=0, -# taxa.basal=0.2, -# taxa.morte=0, -# incl.birth=1529.076, -# incl.death=0, -# density.type=1, -# death.mat=1, -# landscape=land, -# tempo=30, -# ini.config=1, -# out.code=234) +eventos<-function(arquivo, npop){ + nlines <- as.integer(unlist(strsplit(system(paste("wc -l", arquivo), intern=TRUE), split=" "))[1]) + dados = file(arquivo, "r") + dpaisagem = readLines(dados, n=9) + + + IDmax <- 0 + indID <- as.integer(rep(0, nlines)) + indM <- as.double(rep(-1, nlines)) + + for(i in 1:npop) + { + lin = readLines(dados, n=1) + lin<-strsplit(lin, " ") + line <- unlist(lin) + + IDmax <- IDmax + 1 + indM[IDmax]<- line[6] + + + } + + + lin = readLines(dados, n=1) + + while(lin != "EOF") + { + lin<-strsplit(lin, " ") + line <- unlist(lin) + acao <-strtoi(line[2]) + ID<-strtoi(line[3]) + + if(acao == 1) + { + + indID[ID] <- indID[ID]+1 + + IDmax <- IDmax + 1 + indM[IDmax]<- line[7] + + } + lin = readLines(dados, n=1) + } + + close(dados) + + return(resp<-cbind(indM[1:IDmax],indID[1:IDmax])) + +} + +Landscape <- function (numb.cells = 100, cell.size = 1, land.shape = 1, type=c("random","blob"), bound.condition=0, cover=1) { + type=match.arg(type) + if(cover < 0 || cover > 1) { + stop("Error creating landscape. Cover must be between 0 and 1") + } + # scape represents the actual landscape + scape <- rep(1, numb.cells*numb.cells) + if(cover < 1) { + NtoRemove=round((1-cover)*numb.cells*numb.cells); + if(type=="random") { + while(NtoRemove>0) + { + i=round(runif(1,0,numb.cells-1)); + j=round(runif(1,0,numb.cells-1)); + # tests to see if this point has already been removed + if(scape[1+numb.cells*j+i] == 1) { + NtoRemove = NtoRemove - 1 + scape[1+numb.cells*j+i] = 0 + } + } + } + if(type=="blob") { + i=round(runif(1,0,numb.cells-1)); + j=round(runif(1,0,numb.cells-1)); + while(NtoRemove>0) + { + # tests to see if this point has already been removed + if(scape[1+numb.cells*j+i] == 1) { + NtoRemove = NtoRemove - 1 + scape[1+numb.cells*j+i] = 0 + } + # Draft a new point to be removed (random walk!) + if(sample(1:2,1) == 1) { + i = i + (-1)**sample(1:2,1) + } else { + j = j + (-1)**sample(1:2,1) + } + if(i == -1) { i=numb.cells-1} + if(i == numb.cells) { i=1} + if(j == -1) { j=numb.cells-1} + if(j == numb.cells) { j=1} + } + } + } + land <- list(numb.cells = numb.cells, cell.size=cell.size, land.shape=land.shape, type=type, bound.condition=bound.condition, cover=cover, scape=scape) + class(land) <- "landscape" + return(land) +} -# TWoPlot <- function(pop, land, col1="gray20", col2="gray70") { -# n = land$numb.cells -# s <- seq(-n*land$cell.size/2, n*land$cell.size/2, length=n) # creates the x- and y- sequences for image -# if (sum(land$scape) == n*n) { -# color = col1 -# } else { -# color = c(col2, col1) -# } -# image(s, s, matrix(land$scape,ncol=n), col=color) -# points(pop, pch=4, col=2) -# } -# TWoPlot(teste, land) -#plot(teste1, xlim=c(-100,100), ylim=c(-100,100)) -#dim(teste1) -## Tamanho de populacao apos t=6 de 100 repeticoes -#pop.size<- numeric() -#for (i in 1:20) -# { -# pop.size[i] = -# nrow( -# TWoLife(raio=0.1, N=80, AngVis=360, passo=5, move=0.1, taxa.basal=0.6, -# taxa.morte=0.1, -# incl.birth=0.5/0.01, incl.death=0, numb.cells=200, cell.size=1, land.shape=1, -# density.type=0, death.mat=7,bound.condition=0, cover=1, tempo=6)) -# } +TWoLife <- function ( + raio=0.1, + N=80, + AngVis=360, + passo=5, + taxa.move=0.5, + taxa.basal=0.6, + taxa.morte=0.1, + incl.birth=0.5/0.01, + incl.death=0, + density.type=0, + death.mat=7, + move.mat=7, + landscape, + tempo=20, + ini.config=0, + out.code="1", + genotype_means=(rep(1, N)), + width_sds=(rep(0, N)), + points=(rep(0, N)), + Null= FALSE, + initialPosX= NULL, + initialPosY= NULL +) +{ + if(class(landscape) != "landscape") { + stop("Error in function TWoLife: you must provide a valid landscape. See ?Landscape") + } + #if(raio>landscape$numb.cells*landscape$cell.size/2) + #{stop("Error in function TWoLife: the radius must be lower than or equal to the half of landscape side (radius <= numb.cells*cell.size/2)")} + + saida.C <- .C("TWoLife", + as.double(raio),# 1 + as.integer(N),# 2 + as.double(AngVis),# 3 + as.double(passo),# 4 + as.double(taxa.move),# 5 + as.double(taxa.basal),# 6 + as.double(taxa.morte),# 7 + as.double(incl.birth),# 8 + as.double(incl.death),# 9 + as.integer(landscape$numb.cells),# 10 + as.double(landscape$cell.size),# 11 + as.integer(landscape$land.shape),# 12 + as.integer(density.type),# 13 + as.double(death.mat), # 14 + as.double(move.mat), # 15 + as.integer(ini.config), #16 + as.integer(landscape$bound.condition), #17 + as.double(landscape$scape), #18 + as.double(tempo), #19 + as.double(initialPosX), #20 + as.double(initialPosY), # 21 + as.double(genotype_means), # 22 + as.double(width_sds), # 23 + as.integer(points), # 24 + as.logical (Null), # 25 + as.integer(0), #26 + as.double(rep(0, 500000)), # 27 + as.double(rep(0,500000)), # 28 + as.double(rep(0, 500000)), # 29 + as.character(out.code)# 30 + ## verificar se precisa definir o tamanho e se isto nao dará problemas (dois ultimos argumentos) + ) + n <- saida.C[[26]] + x <- saida.C[[27]] + y <- saida.C[[28]] + m <- saida.C[[29]] + x <- x[1:n]; y <- y[1:n] ; m <- m[1:n]; + return(data.frame(x=x,y=y, m=m)) +} -## esperado: capacidade de suporte -# Support <- function(taxa.basal=0.6, taxa.morte=0.1, incl.birth=0.5/0.01, -# incl.death=0, numb.cells=200, cell.size=2) { -# densi.max = (taxa.basal-taxa.morte)/(incl.birth+incl.death) -# return ((numb.cells*cell.size)^2 * densi.max) -# } -## Media das simulacoes -#print(pop.size - Support()) -#print(mean(pop.size - Support())) diff --git a/TWoLife.Rproj b/TWoLife.Rproj new file mode 100644 index 0000000..8e3c2eb --- /dev/null +++ b/TWoLife.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX diff --git a/TWoLife.cpp b/TWoLife.cpp index 6002cf8..642b665 100644 --- a/TWoLife.cpp +++ b/TWoLife.cpp @@ -4,109 +4,208 @@ #include #include #include -#include "paisagem.cpp" -#include "individuo.cpp" +#include "paisagem.cpp" +#include "individuo.cpp" #include #include using namespace std; -/** \file TWoLife.cpp \brief Arquivo usado na integração R/C - * - * Este arquivo define a função TWoLife em C, que será chamada a partir do R */ +/** \file TWoLife.cpp \brief File used for R/C++ interaction + * + * This file defines the TWoLife function in C, to be be called by R */ - -/** Funcao que recebe os parâmetros da simulação e a executa, - // baseada no texto em http://users.stat.umn.edu/~geyer/rc/ +/** The function recieves the inputed parameters and executes the simulation + // Based on the text avaiable on http://users.stat.umn.edu/~geyer/rc/ // - // NOTE que funções em C++ precisam ser indicadas como extern "C" - // para poderem ser acessadas facilmente pela interface do R! + // Note that functions in C++ need to be indicated with the " extern "C" " + // as to be easely accessd by the R interface // - // Veja a descrição da classe \ref paisagem para o significado dos parâmetros + // See the description off the \ref paisagem class for parameter meenings // */ -extern "C" void TWoLife (double * raio, int * N, double * angulo_visada, double * passo, double * move, - double * taxa_basal, double * taxa_morte, double * incl_b, double * incl_d, - int * numb_cells, double * cell_size, int * land_shape, int * density_type, - double * death_mat, int * inipos, int * bound_condition, int * scape, double * tempo, int * nPop, double * x, double * y, int * outCode) -{ - // This sequence creates an attribute containing the output file name. The template is output-00000.txt. - string fileNAME = "output-00000.txt"; - stringstream tmps; - tmps<get_tamanho() << endl; - outputSIM << "Tamanho pixel: " << cell_size[0] << endl; - outputSIM << "Pixels lado: " << numb_cells[0] << endl; - outputSIM << "N fragmentos: " << floresta->get_numb_patches() << endl; - outputSIM << "Area fragmentos: "; - double sum = 0; - //outputSIM << floresta->get_patch_area(0) << " "; //Imprimir area matriz? - for (int i = 1; i <= floresta->get_numb_patches(); i++) - { - sum += floresta->get_patch_area(i); - outputSIM << floresta->get_patch_area(i) << " "; - } - outputSIM << "\nProporção de habitat: " << sum/(floresta->get_tamanho()*floresta->get_tamanho())<< "\n\n\n\n"; - - for(unsigned int i=0; iconta_individuos();i++) - { - outputSIM << floresta->tempo_do_mundo << " " << floresta->get_individuos(i)->get_id() << " " << floresta->get_individuos(i)->get_patch() << " " << floresta->get_individuos(i)->get_x() << " " << floresta->get_individuos(i)->get_y() << endl; - } - while (floresta->tempo_do_mundo < tempo[0] && floresta->conta_individuos() > 0) - { - int ind_neo = floresta->sorteia_individuo(); - int acao = floresta->sorteia_acao(ind_neo); - floresta->atualiza_tempo(ind_neo); +/* + Function to be called by R + All the parameters are Pointers! + +Input "inputs": + Paran double* raio - Radius of density dependance influence + Paran int* N - Number of individuals at the start of the simulation + Paran double* angulo_visada - Angle used for orientation when dispersing + Paran double* passo - The Lenght distance of a dispersal event (random walk) or maximum dispersal radius (Selection step) + Paran double* move - The rate at which the individuals disperse + Paran double* taxa_basal - The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours) + Paran double * taxa_morte - The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) + Paran double * incl_b - The slope of the birth density dependence function + Paran double * incl_d - The slope of the death density dependence function + Paran int * numb_cells - Number of pixels in a square landscape side (if land_shape == 1) + Paran double * cell_size - Resolution Lenght of a square pixel side + Paran int * land_shape - Shape of the landscape (0= circle, 1= square) + Paran int * density_type - Density type (0 = global, 1 = local/within a individual radius) + Paran double * death_mat - Constant that indicates how many times higher the death rate should be on non-habitat pixels + Paran double * move_mat - Constant that indicates how many times lower the movement rate should be on non-habitat pixels + Paran int * inipos - The initial postion of individuals (0 = origin, 1 = random, 2 = normaly distributed with mean on origin, 3 = inputed initial coordinates) + Paran int * bound_condition - The boundary condition type affects how individuals interact with the edges of the landscape (0 = absortive, 1 = periodical (pacman),2 = reflexive) + Paran int * scape - Vector containing the evironmental values of the landscape pixels (in the bynary case: 0= matrix, 1= habitat) + Paran double * tempo - Defines the maximum number of events to happen until the end of the simulation + Paran double * initialPosX - Vector containing the x coordinates initial individuals + Paran double * initialPosY - Vector containing the y coordinates initial individuals + Paran double * genotype_means - Vector containing the genotypical trait means of the initial individuals + Paran double * width_sds - Vector containing the standard deviations of the initial individuals + Paran int * points - Vector containing the numer of points initial individuals sample when selecting habitats + Paran bool * Null - Booling value switching simuluations to neutral state (all individuals acting as an average individual) + + Output "inputs": + Paran int * nPop - Variable for storing the number of individuals at the end of the simulation + Paran double * x - Vector to store the final x coodinate of the living individuals at the end of the simulation + double * y - Vector to store the final y coodinate of the living individuals at the end of the simulation + Paran char ** outCode - A value/name used for identifying the simulatio run + */ - int ID_neo = floresta->get_individuos(ind_neo)->get_id(); - double x_neo = floresta->get_individuos(ind_neo)->get_x(); - double y_neo = floresta->get_individuos(ind_neo)->get_y(); - int patch_neo = floresta->get_individuos(ind_neo)->get_patch(); - - bool emigrou = floresta->realiza_acao(acao, ind_neo); - floresta->update(); - if(acao==2 && !emigrou) - { - x_neo = floresta->get_individuos(ind_neo)->get_x(); - y_neo = floresta->get_individuos(ind_neo)->get_y(); - patch_neo = floresta->get_individuos(ind_neo)->get_patch(); - } - else if(acao==2 && emigrou) - acao = 3; - - outputSIM << floresta->tempo_do_mundo << " " << acao << " " << ID_neo << " " << patch_neo << " " << x_neo << " " << y_neo << endl; - } - outputSIM<< "EOF\n"; - outputSIM.close(); //end of output file - - - *nPop = floresta->conta_individuos(); - for (int i =0; i < *nPop; i ++) { - x[i] = floresta->get_individuos(i)->get_x(); - y[i] = floresta->get_individuos(i)->get_y(); - } //DUVIDA: porque x[i] e y[i] nao tem asterisco antes? - delete floresta; - PutRNGstate(); +extern "C" void TWoLife (double * raio, + int * N, + double * angulo_visada, + double * passo, + double * taxa_move, + double * taxa_basal, + double * taxa_morte, + double * incl_b, + double * incl_d, + int * numb_cells, + double * cell_size, + int * land_shape, + int * density_type, + double * death_mat, + double * move_mat, + int * inipos, + int * bound_condition, + double * scape, + double * tempo, + double * initialPosX, + double * initialPosY, + double * genotype_means, + double * width_sds, + int * points, + bool * Null,// Output "inputs" + int * nPop, + double * x, + double * y, + double * m, + char ** outCode) +{ + // Names the output file name + stringstream tmps; // Creates an object of the stringstream to recieve the output name + tmps<get_tamanho() << endl; // Outputs the lenght of the landscape side + outputSIM << "Tamanho pixel: " << cell_size[0] << endl; // Outputs the lenght of a pixel side + outputSIM << "Pixels lado: " << numb_cells[0] << endl; // Outputs the number of pixels on the landscapeside + outputSIM << "N fragmentos: " << floresta->get_numb_patches() << endl; // Outputs the number of fragments on the landscape + outputSIM << "Area fragmentos: "; // Start Outputing the area of each of the fragments + + + double sum = 0;// Creates a accumulator variable + //outputSIM << floresta->get_patch_area(0) << " "; //Imprimir area matriz? + for (int i = 1; i <= floresta->get_numb_patches(); i++) // Passes through each patch + { + sum += floresta->get_patch_area(i); // Sums a patch area to the acumulator + outputSIM << floresta->get_patch_area(i) << " "; // Outputs the area of each fragment + } + + outputSIM << "\nProporção de habitat: " << sum/(floresta->get_tamanho()*floresta->get_tamanho())<< "\n\n\n\n"; // Outputs the proportion of the landscape that is habitable + + + + + for(unsigned int i=0; iconta_individuos();i++) // Passes through each starting individual + { + // Outputs all the starting individuals and their informations + outputSIM << floresta->tempo_do_mundo << " " << floresta->get_individuos(i)->get_id() << " " << floresta->get_individuos(i)->get_patch() << " " << floresta->get_individuos(i)->get_x() << " " << floresta->get_individuos(i)->get_y() << " " << floresta->get_individuos(i)->get_genotype_mean() << endl; + } + + /** Performs the simulation */ + while (floresta->tempo_do_mundo < tempo[0] && floresta->conta_individuos() > 0 && floresta->conta_individuos() < 2000) // Performs the simulation until the maximum desired time is reached + { + int ind_neo = floresta->sorteia_individuo(); // Drafts an individual to perform an action + int acao = floresta->sorteia_acao(ind_neo); // Drafts an action for the individual to perform + floresta->atualiza_tempo(ind_neo); // Updstes the world Clock + + int ID_neo = floresta->get_individuos(ind_neo)->get_id(); // Obtains the ID of the drafted individual + double x_neo = floresta->get_individuos(ind_neo)->get_x(); // Obtains the X coordinate of the drafted individual + double y_neo = floresta->get_individuos(ind_neo)->get_y(); // Obtains the Y coordinate of the drafted individual + int patch_neo = floresta->get_individuos(ind_neo)->get_patch(); // Obtains the curent patch of the drafted individual + double m_neo = floresta->get_individuos(ind_neo)->get_genotype_mean(); + + bool emigrou = floresta->realiza_acao(acao, ind_neo); // Makes the drafted individual perform the chosen action + + if(acao==2 && emigrou) //Checks if the action performed was a migration event AND the drafted individual left the landscape boundary + floresta->update(0, ind_neo); // Updates the drafted individual migration event as a death event + else // Comtemplates all other cases + floresta->update(acao, ind_neo); // Updates the individual according to the performed action + + if(acao==2 && !emigrou) //Checks if the action performed was a migration event AND the drafted individuals new location is within the landscape boundary + { + x_neo = floresta->get_individuos(ind_neo)->get_x(); // Obtains the individuals current x coordinate + y_neo = floresta->get_individuos(ind_neo)->get_y(); // Obtains the individuals current y coordinate + patch_neo = floresta->get_individuos(ind_neo)->get_patch(); // Obtains the individuals current patch + } + else if(acao==2 && emigrou) //Checks if the action performed was a migration event AND the drafted individual left the landscape boundary + acao = 3; // Sets the output entry of the individual as an emigration event + + // Outputs the actions performed by the drafted individual + outputSIM << floresta->tempo_do_mundo << " " << acao << " " << ID_neo << " " << patch_neo << " " << x_neo << " " << y_neo << " " << m_neo << endl; + } + + + outputSIM<< "EOF\n";// Outbuts a warning stating the simulation was run to completion + outputSIM.close(); //Closes the output file + + //Outputs the remaining individuals + *nPop = floresta->conta_individuos(); // Checks how many individuals are left + for (int i =0; i < *nPop; i ++) { // Passes through each remaining individual + x[i] = floresta->get_individuos(i)->get_x(); // Outputs its x coordinate + y[i] = floresta->get_individuos(i)->get_y(); // Outputs its y coordinate + m[i] = floresta->get_individuos(i)->get_genotype_mean(); + } + + delete floresta; // Deletes the object of the landscape/paisagem class + + PutRNGstate(); //Ends the random seed based on the R (For stand-alone one could use "srand(time(0));") } diff --git a/individuo.cpp b/individuo.cpp index 0f6bc11..deba091 100644 --- a/individuo.cpp +++ b/individuo.cpp @@ -8,157 +8,302 @@ using namespace std; +// Initializes max ID unsigned long individuo::MAXID = 0; -/** \brief Exceção no indivíduo. - * - * Esta classe implementa uma exceção ocorrida na classe indivíduo. Atualmente, é usada para detectar valores impossíveis - * passados no construtor, mas vai ser expandida para compreender mais casos (\ref TBI). */ +/** \brief Individual exception + Displays warning messages for impossible values + (\ref TBI). */ class individuo_exception: public exception { - virtual const char* what() const throw() - { - return "Aconteceu uma exceção na classe indivíduo!"; - } + virtual const char* what() const throw() + { + return "Aconteceu uma exceção na classe indivíduo!"; + } } myex; -individuo::individuo(double x, double y, int especie, double taxa_morte, - double orientacao, double angulo_visada, - double passo, double move, double raio, - double taxa_basal, int semente, //retirar int semente - double incl_b, double incl_d, - double death_mat, int dens_type): - id(++MAXID), // pega o próximo ID livre - x(x), - y(y), - especie(especie), - taxa_morte(taxa_morte), - move(move), - passo(passo), - orientacao(orientacao), - ang_visada(angulo_visada), - raio(raio), - taxa_basal(taxa_basal), - semente(semente), - incl_birth(incl_b), - incl_death(incl_d), - const_d_matrix(death_mat), - dens_type(dens_type) +// Class Constructor +individuo::individuo(double x, + double y, + int especie, + double taxa_morte, + double orientacao, + double angulo_visada, + double passo, + double taxa_move, + double raio, + double taxa_basal, + int semente, + double incl_b, + double incl_d, + double death_mat, + double move_mat, + int dens_type, + vector genotype, + vector width, + int points): +// This section is similar to regular atribution and is run before brackets +id(++MAXID), // pega o próximo ID livre +x(x), +y(y), +especie(especie), +taxa_morte(taxa_morte), +orientacao(orientacao), +ang_visada(angulo_visada), +passo(passo), +taxa_move(taxa_move), +move(taxa_move), +raio(raio), +taxa_basal(taxa_basal), +semente(semente), +incl_birth(incl_b), +incl_death(incl_d), +const_d_matrix(death_mat), +const_m_matrix(move_mat), +dens_type(dens_type), +points(points) { - if (taxa_morte < 0) throw myex; - if (passo < 0) throw myex; - if (move < 0) throw myex; - if (raio < 0) throw myex; - if (taxa_basal < 0) throw myex; - - if(incl_birth!=0 && incl_death!=0) - { - this->densi_max = (taxa_basal-taxa_morte)/(incl_b+incl_d); - this->birth_death_eq = taxa_morte+incl_d*((taxa_basal-taxa_morte)/(incl_b+incl_d)); - } + // Checks for impossible values + if (taxa_morte < 0) throw myex; + if (passo < 0) throw myex; + if (taxa_move < 0) throw myex; + if (raio < 0) throw myex; + if (taxa_basal < 0) throw myex; + + if(incl_birth!=0 && incl_death!=0) // Checks if the birth and death rates are density dependant + { + this->densi_max = (taxa_basal-taxa_morte)/(incl_b+incl_d); // Sets the maximum density + this->birth_death_eq = taxa_morte+incl_d*((taxa_basal-taxa_morte)/(incl_b+incl_d)); // Sets the point where rates are equal + } + + this->genotype_mean = genotype; //Sets the individuals genotype value (or values) + + this->width_sd = width; //Sets the individuals width value (or values) + + this->rdn_noise = 0; //Sets random noises (mutation) to zero (for future implementation) + + for (int i=0; igenotype_mean.size(); i++) { // passes by the genotype means of the individual + + this->env_optimum.push_back(this->rdn_noise+this->genotype_mean[i]); // Generates phenotype by adding a random value to the genotype value + } } /** \brief função para reprodução assexuada de um indivíduo - * + * * Faz uma cópia do indivíduo pai. Recebe um ponteiro dereferenciado (nomeado rhs) e usa isto para * Criar novo individuo com os mesmos valores de atributo do pai, exceto * - id (veja \ref individuo::get_id) * - vizinhos (veja \ref individuo::set_vizinhos) - * - tempo para evento (veja \ref individuo::get_tempo) - * Usa notacao :atributo(valor) ao inves de atribuicão. + * - tempo para evento (veja \ref individuo::get_tempo) + * Usa notacao :atributo(valor) ao inves de atribuicão. * Funcao chamada por paisagem::realiza_acao() quando a ação é um nascimentto * @param rhs ponteiro dereferenciado para o pai */ -individuo::individuo(const individuo& rhs) - :id(++MAXID), - x(rhs.x), - y(rhs.y), - especie(rhs.especie), - taxa_morte(rhs.taxa_morte), - move(rhs.move), - passo(rhs.passo), - ang_visada(rhs.ang_visada), - raio(rhs.raio), - taxa_basal(rhs.taxa_basal), - tipo_habitat(rhs.tipo_habitat), - semente(rhs.semente), - incl_birth(rhs.incl_birth), - incl_death(rhs.incl_death), - const_d_matrix(rhs.const_d_matrix), - dens_type(rhs.dens_type), - birth_death_eq(rhs.birth_death_eq) - -{ //precisamos dessa chave e da que fecha ela? - + +/** Copy constructor, used for generating new individuals by assexual reproduction + All the characteristics of the parent individual will be copyed, exept: + - id (veja \ref individuo::get_id/ will be set to new value) + -list of Neighbours (veja \ref individuo::set_vizinhos/ will be updated) + - time until nex event (veja \ref individuo::get_tempo/ will be drafted) + Paran const individuo& rhs - Parent individual */ +individuo::individuo(const individuo& rhs): +id(++MAXID), +x(rhs.x), +y(rhs.y), +especie(rhs.especie), +taxa_morte(rhs.taxa_morte), +ang_visada(rhs.ang_visada), +passo(rhs.passo), +taxa_move(rhs.taxa_move), +move(rhs.taxa_move), +raio(rhs.raio), +taxa_basal(rhs.taxa_basal), +semente(rhs.semente), +incl_birth(rhs.incl_birth), +incl_death(rhs.incl_death), +const_d_matrix(rhs.const_d_matrix), +const_m_matrix(rhs.const_m_matrix), +dens_type(rhs.dens_type), +tipo_habitat(rhs.tipo_habitat), +birth_death_eq(rhs.birth_death_eq) + +{ + this->orientacao = runif(0,360); // Sets an ramdom orientation + + if (rhs.genotype_mean.size()==1) // Checks if this is not a null simulation + { + this->rdn_noise= 0; //Sets random noises (recombination/ mutation) to zero (for future implementation) + //this->genotype_mean.push_back(rhs.genotype_mean[0] + runif(-1.0,1.0)); + this->genotype_mean.push_back(rhs.genotype_mean[0]+rdn_noise); //Sets the individuals genotype value + this->width_sd= rhs.width_sd; //Sets the individuals width value (or values) + //this->rdn_noise= 0; //Sets random noises (recombination/ mutation) to zero (for future implementation) + this->env_optimum.push_back(this->rdn_noise+this->genotype_mean[0]); // Generates the genotype of the offspring by adding a random value to its parents genotype value + + } + else // Checks if this is a null simulation + { + + this->rdn_noise= 0; + + for (int i=0; igenotype_mean.push_back(this->rdn_noise+rhs.genotype_mean[i]); // Generates the genotype of the offspring by adding a random value to its parents genotype value + } + + //this->genotype_mean= rhs.genotype_mean; + this->width_sd = rhs.width_sd; //Sets the individuals width values + this->rdn_noise = 0; //Sets random noises (recombination/ mutation) to zero (for future implementation) + + for (int i=0; igenotype_mean.size(); i++) { // passes by the genotype means of the individual + + this->env_optimum.push_back(this->rdn_noise+this->genotype_mean[i]); // Generates phenotype by ading a random value to the genotype values + } + + } + this->points=rhs.points; //Sets the points equal to the parents } -/** \brief Método de atualização dos indivíduos +/** \brief Método de atualização dos indivíduos * Esta função é a camada de atualização dos indivíduos * A cada execução deste método: * - É atualizada a taxa de nascimento de acordo com a densidade de vizinhos no raio de vizinhanca obtido com paisagem::atualiza_vizinhos * - Sorteia o tempo de acordo com as novas taxas * As taxas de morte e movimentação no momento fixas. Mas tambem serão funções da densidade de vizinhos (\ref TBI). */ + +/** Function that updates both death and birth rates of individuals based on the number of individuals within the area of their density dependance radius +It them calls a function to draft the time needed for execution of that individual based on its birth, death and dispersal rates*/ void individuo::update(double dens) { - double densi = dens; // densidade inclui n de vizinhos + o individuo - if(this->tipo_habitat==0) - { - this->birth = 0; - // ToDo: Implementar aqui modelo mais geral para mortalidade na matriz. Aqui a denso dependencia é igual à do habitat, só muda a mortalidade basal que é maior que no habitat. - this->death = this->const_d_matrix*this->taxa_morte+this->incl_death*densi; - } - else - { - this->birth = this->taxa_basal-this->incl_birth*densi; - this->death = this->taxa_morte+this->incl_death*densi; - if(this->birth<0){this->birth=0;} - } - - this->sorteiaTempo(); + double densi = dens; // Density of neighbour individuals (includes focal one) + + + if (this->width_sd[0] == 0) { // Checks if individuals are completely specialists + if(this->tipo_habitat==0) // Checks if the individual is currently on the matrix + { + this->birth = 0; // Sets birth rate to 0 + // ToDo: Implementar aqui modelo mais geral para mortalidade na matriz. Aqui a denso dependencia é igual à do habitat, só muda a mortalidade basal que é maior que no habitat. + this->death = this->const_d_matrix*this->taxa_morte+this->incl_death*densi; // Sets the higer mortality rate atributed to nonhabita patches + this->move = this->const_m_matrix*this->taxa_move; // Sets the higher/lower movement rates atributed to nonhabita patches + } + else + { + this->birth = this->taxa_basal-this->incl_birth*densi; // Computes the actual birth rate on habitat patch (that is influenced by the density of neighbours) + this->death = this->taxa_morte+this->incl_death*densi; // Computes the actual death rate on habitat patch (that is influenced by the density of neighbours) + this->move = this->taxa_move; // Sets the standard movement rates + } + } + else{ + + //this->move = this->taxa_move; + + this->birth = this->taxa_basal-this->incl_birth*densi; // Computes the actual birth rate on habitat patch (that is influenced by the density of neighbours) + this->death = ((this->const_d_matrix*this->taxa_morte)-((dnorm_sum(this->tipo_habitat, this->env_optimum, this->width_sd)/dnorm(this->env_optimum[0],this->env_optimum[0],this-> width_sd[0]))*((this->const_d_matrix*this->taxa_morte)-this->taxa_morte))); // Computes the actual death rate on habitat patch (that is influenced by the suitability of its current habitat) + + } + + if(this->birth<0) //Checks if the birth rate is lower than possible + { + this->birth=0; // Sets to the lowest possible value to Birth + + } + + this->sorteiaTempo(); // Calls the function to draft the time needed to execute the next action } - +// Function that drafts one of the three possible actions (acounting for their respective taxes) and returns the drafted action to the landscape void individuo::sorteiaTempo() { - this->tempo_evento = rexp(1.0/(this->move+this->birth+this->death)); + this->tempo_evento = rexp(1.0/(this->move+this->birth+this->death)); //Sets the time to the individual } - +// Function that drafts one of the three possible actions (acounting for their respective taxes) and returns the drafted action to the landscape int individuo::sorteia_acao() { - vector probscum; - double total = this->death+this->birth+this->move; - probscum.push_back(this->death/total); - probscum.push_back(probscum[0]+(this->birth/total)); - probscum.push_back(probscum[1]+(this->move/total)); - double evento; - evento = runif(0.0,1.0); - - int decisao; - for(unsigned int i=0; i probscum; //creates a vector for storing the probabilities + double total = this->death+this->birth+this->move; // Sums all rates + probscum.push_back(this->death/total); // stores the death event probability at the first position + probscum.push_back(probscum[0]+(this->birth/total)); // stores the birth event probability at the secon position + probscum.push_back(probscum[1]+(this->move/total)); // stores the dispersal event probability at the third position + double evento; // temporary variable for storing a random value + evento = runif(0.0,1.0); // Samples between 0 and 1 + + int decisao; // temporary variable for storing dracted event + for(unsigned int i=0; ievento)//encontrar o primeiro maior valor + if(probscum[i]>evento)// Finds the first action probability higher than the drafted number { - decisao = i; - return decisao;//retorna a decisao tomada para o mundo (0 = morte, 1 = nascer, 2 = andar) - break; - } - } - return probscum.size()-1; + decisao = i; //Stores the drated action + return decisao; // Returns the drafted action to the landscape (0 = death, 1 = birth, 2 = dispersal) + break; + } + } + return probscum.size()-1; // returns sorted action } -void individuo::anda(bool aleatorio) +//Function that dislocates the XY corrdinates of an individual by a "passo" lenght distance and oriented by the angle the individuals is currently facing added of a random component (angulo_vizada) (if angulo_vizada== 360 -> Random Walk) +void individuo::anda() { - //a cada movimentacao ele muda de direcao; random walk dentro de um angulo de visao - if (aleatorio) { - this->orientacao = runif(-180.0,180.0);//random way point para sortear uma direcao qualquer - } else { - this->orientacao+= runif(-ang_visada/2.0, ang_visada/2.0);//random way point - } - double oriRad=this->orientacao*M_PI/180.0;//tranforma em radianos para poder calcular as distancias das posicoes x e y - double dx= cos(oriRad)*this->passo; - double dy= sin(oriRad)*this->passo; - this->x+=dx; - this->y+=dy; + this->orientacao+= runif(-this->ang_visada/2.0, this->ang_visada/2.0);//random way point to draft any direction + if(this->orientacao < 0) // Checks for impossible (negative)values + this->orientacao +=360; // Transforms in absolute dgrees + else if(this->orientacao >= 360) // Checks for impossible (big)values + this->orientacao -= 360; // Fixes spilage + double oriRad=this->orientacao*M_PI/180.0;// Tranforms to radians to calculate the new XY coordinates + + double dx= cos(oriRad)*this->passo; // Calculates the new x coordinate + double dy= sin(oriRad)*this->passo; // Calculates the new y coordinate + this->x+=dx; // Sets new x coordinate + this->y+=dy; // Sets new y coordinate } + +//Function that returns the value of the probability density function for the normal distribution +double individuo::dnorm(double x ,double mean, double sd){ + + return (1/(sd*sqrt(2*M_PI))*(exp(-1*(pow(x-mean, 2)/(2*pow(sd, 2)))))); + +} + +//Function that returns the value the sum of several probability density function for the normal distribution +double individuo::dnorm_sum( double x ,vector mean, vector sd){ + + double probcumsum=0; // Temporary variable for storing the sum + + for (int i=0; ipoints]; // Temporary vector for storing the scores + double cumsum=0, choice=0, score=0; // Temporary variable for storing the drafting components + int final_pos= 0; // Temporary variable for storing the chosen index + + for (int i=0; ipoints; i++) { // Passes by the points + + scores[i] = exp(dnorm_sum(possibilitities[i][2], this->genotype_mean, this->width_sd)); // Atributes a score to the coordinate based on the probability density function of the migrating individual + cumsum += scores[i]; // Sums that score to the total + } + + choice= runif(0.0,1.0); // Drafts between 0 and 1 + + for (int i=0; ipoints; i++) { // Passes by the points + + score += (scores[i]/cumsum); // Computes the comparative value of the ranks + + if (score > choice) { // Checks the comparative value of the ranks has exceeded the drafted value + + final_pos= i; // Choses the index of the selected coordinate + break; + } + } + + this->x = possibilitities[final_pos][0]; // Sets new x coordinate + this->y = possibilitities[final_pos][1]; // Sets new y coordinate +} + diff --git a/individuo.hpp b/individuo.hpp index 9f60dd5..7f0cb82 100644 --- a/individuo.hpp +++ b/individuo.hpp @@ -4,154 +4,233 @@ using namespace std; -/** \brief A classe individuo representa um agente da simulação. +/** \brief The individual class represents an agent of the simulation * - * Esta classe contém informações referentes a cada indivíduo, incluindo sua localização e estado, uma lista contendo ponteiros para os vizinhos próximos, tamanho do passo de movimentação, etc. Esta classe NÃO contém métodos que sejam de responsabilidade do ambiente, como o método de atualizar os vizinhos */ + * This class conatins informations pertinent to an individual, including it's location, state, a vector containing pointers to close neighbours, lenght distance of a dispersal event, rate of dispersal, etc. This class DOES NOT contain methods of landscape/paisagem responsability, mutch like the update of the neighbour vector of each individual. */ class individuo { private: - //propriedades privadas - /** Identificador único de cada indivíduo */ + + // Private properties + + /** Unique individual identifier */ const unsigned long id; - /** Identificador máximo já usado */ + /** Current maximum individual identifier */ static unsigned long MAXID; - /** Posição X da localização do indivíduo */ + /** The X coordinate of the individual */ double x; - /** Posição Y da localização do indivíduo */ - double y; - /** Espécie a qual o indivíduo pertence. Cada espécie será representada por um número único, TBI*/ + /** The Y coordinate of the individual */ + double y; + /** The species identifier number of the individual (TBI)*/ const int especie; - /** Taxa de natalidade final, após a consideração de todos os efeitos */ + /** Final/current birth rate of the individual (all effects considered) */ double birth; - /** Taxa intrínseca de mortalidade quando no habitat e na ausência de outros indivíduos na vizinhança local */ + /** The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) */ const double taxa_morte; - /** Taxa de movimentação */ - const double move; - /** Tamanho do passo dado a cada movimentação */ + /** The basal Dispersal rate of the individual */ + const double taxa_move; + /** Final/current Dispersal rate of the individual (all effects considered) */ + double move; + /** The Lenght distance of a dispersal event (random walk) or maximum dispersal radius (Selection step) */ const double passo; - /** Ângulo de orientação sorteado para a movimentação */ + /** Angle that an individual is curently facing */ double orientacao; - /** Ângulo de visada */ + /** Angle used for orientation when dispersing */ const double ang_visada; - /** Tempo sorteado para o próximo evento acontecer com este indivíduo */ + /** Drafted time require for the individual to execute an action */ double tempo_evento; - /** Densidade máxima suportada pelo indivíduo*/ + /** Maximum density the individual is capable of bearing*/ double densi_max; - /** Raio de percepção da densidade */ + /** Radius of density dependant influence */ double raio; - /** Vetor de vizinhanca */ - vector lisViz;//vetor de vizinhanca ----//CAMADA PERCEPTIVA - /** Taxa intrínseca de natalidade quando no habitat e na ausência de outros indivíduos na vizinhança local */ + /** The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours)*/ const double taxa_basal; - /** Identificdor do tipo de habitat do pixel correspondente à atual posição do indivíduo (0 = matriz; 1 = habitat) */ - int tipo_habitat;//CAMADA PERCEPTIVA - /** Identificdor do fragmento correspondente à atual posição do indivíduo (0 = matriz; 1, 2, ... = fragmentos de habitat) */ + /** Identfier of the Type of habitat the individual is currently on (0 = matriz; 1 = habitat) */ + double tipo_habitat;//CAMADA PERCEPTIVA + /** Identfier of the patch of habitat the individual is currently on */ int patch_label; - /** Semente para gerar os numeros aleatorios */ + /** Seed for generating random numers for the individual */ const int semente; - /** Taxa de mortalidade final, após a consideração de todos os efeitos */ - double death; - /** Taxa de natalidade e mortalidade quando N = K ou dens = dens_max */ - double birth_death_eq; - /** Inclinação da curva de denso-depedência para natalidade */ - const double incl_birth; - /** Inclinação da curva de denso-depedência para mortalidade */ - const double incl_death; - /** Constante que indica quantas vezes a mortalidade basal na matriz eh maior que no habitat */ - const double const_d_matrix;// Pensar em como considerar diferentes mortalidades (constantes) em diferentes tipos de matriz - /** Forma como a densidade é calculada (0 = GLOBAL, 1 = LOCAL)*/ - const int dens_type; - //metodos privados - /** Gera um número aleatório, segundo distribuição exponencial, correspondente ao tempo para a ocorrência do próximo evento, levando em consideração as taxas de mortalidade, natalidade e movimentação */ + /** Final/current death rate of the individual (all effects considered) */ + double death; + /** Birth and death rates when they are at equilibrium (populational K) */ + double birth_death_eq; + /** The slope of the birth density dependance function*/ + const double incl_birth; + /** The slope of the death density dependance function */ + const double incl_death; + /** Constant that indicates how many times higher the death rate should be on non-habitat pixels */ + const double const_d_matrix; + /** Constant that indicates how many times lower the movement rate should be on non-habitat pixels*/ + const double const_m_matrix; + /** Density type (0 = global, 1 = local/within a individual radius)*/ + const int dens_type; + /** The genetical optimum environmental value for a individual (where its mortality rate is the lowest)*/ + vector env_optimum; + /** The phenotipical optimum environmental value for a individual (where its mortality rate is the lowest) */ + vector genotype_mean; + /** The standard deviation of environmental usage by a individual, how generalist it is*/ + vector width_sd; + /**The randomly selected displacement value (This is used to dislocate the fitness distribution of the individual) */ + double rdn_noise; + /** The number of samplining corrdinates drafted at each dispersal event */ + int points; + + + // Private Methods + /** Fucntion that generates random numbers, following a exponential distribution, corresponding to the time needed to execute the next action, taking in account the birth, death and dispersal rates. + void sorteiaTempo(); */ void sorteiaTempo(); - + public: - /** Construtor da classe individuo. Deve ser chamado pela paisagem para posicionar os - * indivíduos que já estão na paisagem no início da simulação. */ + /** Vector for storing the neightbours of the individual */ + vector lisViz;//vetor de vizinhanca ----//CAMADA PERCEPTIVA + /** Constructor of the individual/individuo class + Must be called by the landscape/paisagem class to position individuals at the star of the simualtion. */ + individuo( - /** Posição X da localização do indivíduo */ - double x, - /** Posição Y da localização do indivíduo */ - double y, - /** Espécie a qual o indivíduo pertence. Cada espécie é representada por um número único. \ref TBI */ - const int especie, - /** Taxa intrínseca de mortalidade quando no habitat e na ausência de outros indivíduos na vizinhança local */ - const double taxa_morte, - /** Ângulo de orientação sorteado para a movimentação, veja \ref individuo::anda */ - double orientacao, - /** Ângulo de visada, veja \ref individuo::anda */ - const double angulo_visada, - /** Tamanho do passo dado a cada movimentação, veja \ref individuo::anda */ - const double passo, - /** Taxa de movimentação */ - const double move, - /** Raio de percepção da densidade */ - const double raio, - /** Taxa intrínseca de natalidade quando no habitat e na ausência de outros indivíduos na vizinhança local */ - const double taxa_basal, - /** Semente para gerar os numeros aleatorios */ - const int semente, - /** Inclinação da curva de denso-depedência para natalidade */ - const double incl_birth, - /** Inclinação da curva de denso-depedência para mortalidade */ - const double incl_death, - /** Constante que indica quantas vezes a mortalidade basal na matriz é maior que no habitat */ - const double death_mat, - /** Forma como a densidade é calculada (0 = GLOBAL, 1 = LOCAL)*/ - const int dens_type); - /** Construtor de cópia, usado para gerar novos indivíduos por reprodução assexuada. - * Todos os dados do individuo pai serão copiados, com a exceção de: - * - id (veja \ref individuo::get_id) - * - vizinhos (veja \ref individuo::set_vizinhos) - * - tempo para evento (veja \ref individuo::get_tempo) - * */ - individuo(/** Indivíduo pai */ const individuo& rhs); - - /** Reinicia o contador de indivíduos */ + /** The X coordinate of the individual */ + double x, + /** The Y coordinate of the individual */ + double y, + /** The species identifier number of the individual \ref TBI */ + const int especie, + /** The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) */ + const double taxa_morte, + /** Angle that an individual is curently facing \ref individuo::anda */ + double orientacao, + /** Angle used for orientation when dispersing \ref individuo::anda */ + const double angulo_visada, + /** The Lenght distance of a dispersal event (random walk) or maximum dispersal radius (Selection step) \ref individuo::anda */ + const double passo, + /** The basal Dispersal rate of the individual */ + const double taxa_move, + /** Radius of density dependant influence */ + const double raio, + /** The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours) */ + const double taxa_basal, + /** Seed for generating random numers for the individual */ + const int semente, + /** The slope of the birth density dependance function */ + const double incl_birth, + /** The slope of the death density dependance function */ + const double incl_death, + /** Constant that indicates how many times higher the death rate should be on non-habitat pixels */ + const double death_mat, + /** Constant that indicates how many times lower the movement rate should be on non-habitat pixels */ + const double move_mat, + /** Density type (0 = global, 1 = local/within a individual radius)*/ + const int dens_type, + /** The phenotipical optimum environmental value for a individual (where its mortality rate is the lowest) */ + vector genotype, + /** The standard deviation of environmental usage by a individual, how generalist it is*/ + vector width, + /** The number of samplining corrdinates drafted at each dispersal event */ + int points + ); + + /** Copy constructor, used for generating new individuals by assexual reproduction + All the characteristics of the parent individual will be copyed, exept: + - id (veja \ref individuo::get_id/ will be set to new value) + -list of Neighbours (veja \ref individuo::set_vizinhos/ will be updated) + - time until nex event (veja \ref individuo::get_tempo/ will be drafted) + Paran const individuo& rhs - Parent individual */ + individuo(/** Indivíduo pai */ const individuo& rhs); + + /** Function that sets the maximum id to 0 */ static void reset_id() {MAXID = 0;} - /** Retorna o identificador único deste indivíduo */ - const unsigned long get_id() const {return this->id;} - /** Atualiza a lista de vizinhos deste indivíduo. Deve ser chamada a cada passo de tempo pela \ref paisagem */ + + /** Function that returns the ID of an individual */ + const unsigned long get_id() const {return this->id;} + + /** Function that passes the imputed vector of neighbours of an individual (individuals within a radius distance) to the individuals. + Must be called at each time step by the landscape + Paran const vector lis - The list of individuals within aradius dinstace of the focal individual*/ void set_vizinhos (/** Lista dos vizinhos */ const vector lis){this->lisViz=lis;} - /** Atualiza o tipo de hábitat no qual o indivíduo está. Deve ser chamada a cada passo de tempo pela \ref paisagem. */ - void set_habitat (const int tipo){this->tipo_habitat=tipo;} - /** Atualiza o fragmento do indivíduo */ + + /** Function that updates the habitat type of the pixeld the individual is currently on + Mut be called at each time step by the landscape (0= matrix 1= habitat) + Paran const int tipo - Pixel adress on the landscape matrix */ + void set_habitat (const double tipo){this->tipo_habitat=tipo;} + + /** + Function that updates the fragment identfier of the pixel the individual is currently on + Paran const int label - Pixel adress on the patch id matrix */ void set_patch (const int label){this->patch_label=label;} - /** Atualiza a posição X do indivíduo */ + + /** Function that updates the X coordinate of the individual + Paran double i - The new x Coordinate*/ void set_x(/** Nova posição */double i){this->x =i;} - /** Atualiza a posição Y do indivíduo */ + + /** Function that updates the y coordinate of the individual + Paran double i - The new y Coordinate */ void set_y(/** Nova posição */double i){this->y =i;} - /** Retorna a posição X do indivíduo */ + + /** Function that returns the x coordinate of the individual */ inline const double get_x() const {return this->x;} - /** Retorna a posição Y do indivíduo */ + + /** Function that returns the y coordinate of the individual */ inline const double get_y() const {return this->y;} - /** Retorna o raio de percepção do indivíduo */ + + /** Function that returns the density dependance radius of the individual */ const double get_raio() const {return this->raio;} - /** Retorna o tipo de densidade que afeta o indivíduo */ + + /** Function that returns the density type afecting the individual (0= global, 1=local) */ const int get_densType() const {return this->dens_type;} - /** Returns the number of individuals inside the neighbourhood of the individual (it includes the focal individual) */ + + /** Function that returns the number of individuals within the readius of the focal individuals for density calculations (it also includes the focal individual) */ const int NBHood_size() const {return this->lisViz.size()+1;} - /** Retorna o fragmento no qual o indivíduo está */ + + /** Function that retuns the environmental value of the pixel the individual is currently on */ const int get_patch() const {return this->patch_label;} - - // outros metodos publicos - /** Retorna o tempo sorteado para o próximo evento acontecer com este indivíduo. - * \sa - * \ref individuo::update - * \ref paisagem::update */ + + /** Function that returns the step lengh of the individual*/ + const double get_passo() const {return this->passo;} + + /** Function that returns the amounts of sampling coordinates an individual has */ + const int get_points() const {return this->points;} + + /** Function that returns the genotype mean the individual */ + const double get_genotype_mean() const {return this->genotype_mean[0];} + + // Other Public Methods + /** Function that Returns the drafted time for executing an action by the individual based on its birth, death and dispersal rates. + * \sa + * \ref individuo::update + * \ref paisagem::update */ const double get_tempo(){return this->tempo_evento;} - /** Sorteia uma das três ações possíveis de acordo com suas respectivas taxas e retorna a ação sorteada para a paisagem - * (0 = morte, 1 = nascimento, 2 = movimentação) */ + + /** Function that drafts one of the three possible actions (0 = death, 1 = birth, 2 = dispersal) (acounting for their respective taxes) and returns the drafted action to the landscape int sorteia_acao(); */ int sorteia_acao(); - /** Atualiza a taxa de natalidade e/ou de mortalidade baseado na densidade de indivíduos dentro do raio de percepção e sorteia o tempo - * do próximo evento baseado nas taxas de natalidade, de mortalidade e de movimentação - * \sa \ref individuo::get_tempo */ + + /** Function that updates both death and birth rates of individuals based on the number of individuals within the area of their density dependance radius + It them calls a function to draft the time needed for execution of that individual based on its birth, death and dispersal rates + Paran double dens - The number of individuals pondered by the area of the density dependance radius + * \sa \ref individuo::get_tempo */ void update(double dens); - /** Faz com que o indivíduo ande um passo, de tamanho passo. A orientação na qual o indivíduo vai andar é a orientação atual - * (definida no construtor como orientacao) mais um ângulo aleatório dentro do ângulo de visada (angulo_visada). A definição de um - * ângulo de visada de 360 graus equivale a uma caminhada aleatória */ - void anda(/** Passe aleatorio = true para forçar uma caminhada aleatória */ bool aleatorio = true); + + /** Function that dislocates the XY corrdinates of an individual by a "passo" lenght distance and oriented by the angle the individuals is currently facing added of a random component (angulo_vizada) (if angulo_vizada== 360 -> Random Walk) + Paran bool aleatorio = true - */ + void anda(); + + /** Function that returns the value of the probability density function for the normal distribution + Paran double x - Quantile of interest + Paran double mean - Mean value of the distribution + Paran double sd - Standad deviation of the distribution*/ + double dnorm(double x ,double mean, double sd); + + /** Function that returns the value the sum of several probability density function for the normal distribution + Paran double x - Quantile of interest + Paran vector mean - vector containing theMean values of the distribution + Paran vector sd - vector containing the Standad deviations of the distribution*/ + double dnorm_sum( double x ,vector mean, vector sd); + + /** Function that selects from within a given set of coordinates based on the dispersing individuals preference ranks via a softmax function + Paran double possibilitities[][3] - Vector conatining X and Y coordinates (first two cols), and the envinronmental values of that coordinats pixel. (third col) */ + void habitat_selection(double possibilitities[][3]); + }; //Falta mexer no doxygen dos construtores e da update #endif // INDIVIDUO_H + diff --git a/paisagem.cpp b/paisagem.cpp index c8b7f2c..27b0ab4 100644 --- a/paisagem.cpp +++ b/paisagem.cpp @@ -3,210 +3,346 @@ #include #include -paisagem::paisagem(double raio, int N, double angulo_visada, double passo, double move, double taxa_basal, - double taxa_morte, double incl_b, double incl_d, int numb_cells, double cell_size, int land_shape, - int density_type, double death_mat, int inipos, int bound_condition, int scape[]): - tamanho(numb_cells*cell_size), - N(N), - tempo_do_mundo(0), - numb_cells(numb_cells), - cell_size(cell_size), - landscape_shape(land_shape), - boundary_condition(bound_condition), - initialPos(inipos) +// Class Constructor +paisagem::paisagem(double raio, + int N, + double angulo_visada, + double passo, + double taxa_move, + double taxa_basal, + double taxa_morte, + double incl_b, + double incl_d, + int numb_cells, + double cell_size, + int land_shape, + int density_type, + double death_mat, + double move_mat, + int inipos, + int bound_condition, + double scape[], + double initialPosX[], + double initialPosY[], + double genotype_means[], + double width_sds[], + int points [], + bool Null): +// This section is similar to regular atribution and is run before brackets +tamanho(numb_cells*cell_size), +N(N), +tempo_do_mundo(0), +numb_cells(numb_cells), +cell_size(cell_size), +landscape_shape(land_shape), +boundary_condition(bound_condition), +initialPos(inipos) { - for(int i=0;inumb_cells;i++) - { - for(int j=0;jnumb_cells;j++) - { - // transforma o vetor scape recebido do construtor em uma matriz - this->landscape[i][j]=scape[j*numb_cells+i]; - } - } - - //Atribui valores à matriz patches, que determina o fragmento a que cada pixel pertence - int component = 0; - for (int k = 0; k < this->numb_cells; ++k) - for (int l = 0; l < this->numb_cells; ++l) - if (!this->patches[k][l] && this->landscape[k][l]) find_patches(k, l, ++component); - this->numb_patches = component; - - this->patch_area = new double[this->numb_patches+1]; - - for (unsigned int j = 0; jpatch_area[j] = 0; - - for(int i = 0; i < this->numb_cells; i++) - for(int j = 0; j < this->numb_cells; j++) - this->patch_area[patches[i][j]] += 1; - for (unsigned int j = 0; jpatch_area[j] = this->patch_area[j]*this->cell_size*this->cell_size; - - // Calculo do raio dependendo do tipo de densidade. 0 = global, 1 = local (raio), 2 = kernel. - if(density_type==0) - { - raio = this->tamanho/sqrt(M_PI); - } - /* Coloca os indivíduos na paisagem por meio da função populating() */ - this->populating(raio,N,angulo_visada,passo,move,taxa_basal,taxa_morte,incl_b,incl_d,death_mat,density_type); - - for(unsigned int i=0; ipopIndividuos.size(); i++) - { - this->atualiza_vizinhos(this->popIndividuos[i]);//atualiza os vizinhos - this->atualiza_habitat(this->popIndividuos[i]);//retorna o tipo de habitat - this->atualiza_patch(this->popIndividuos[i]); + for(int i=0;inumb_cells;i++) // Passes through the rows + { + for(int j=0;jnumb_cells;j++) // Passes through the columns + { + + this->landscape[i][j]=scape[j*numb_cells+i]; // Re-shapes the scape vector into a matrix + this->patches[i][j]=0; // Initializas the patches vector giving it 0 as a null value } - - - for(unsigned int i=0; ipopIndividuos.size(); i++) - { - double dsty=this->calcDensity(popIndividuos[i]); - this->popIndividuos[i]->update(dsty); //e atualiza o individuo i da populacao - } - + } + + // Atributes values to the landscape pixels, determining the patch they belong to + int component = 0; // Counter variable used for identyfing the patches + for (int k = 0; k < this->numb_cells; ++k) // Passes through the rows + for (int l = 0; l < this->numb_cells; ++l) // Passes through the columns + if (!this->patches[k][l] && this->landscape[k][l]) find_patches(k, l, ++component); // Checks if each landscape pixel has already been assigned a patch identfication value AND if it is a non-matrix patch + this->numb_patches = component; // Stores the number of habitat patches in the landscape + + this->patch_area = new double[this->numb_patches+1]; //Points to a new vector object used for storing the area of all habitats and the matrix + + + for (unsigned int j = 0; jpatch_area[j] = 0; // Initializes the vector giving it 0 as a null value + + for(int i = 0; i < this->numb_cells; i++) //Passes through the rows + for(int j = 0; j < this->numb_cells; j++) //Passes through the columns + this->patch_area[patches[i][j]] += 1; // Indexes the patch_area vector by the patch identfication numbers of each pixel in the landscape so that it adds the number o pixels of each patch + for (unsigned int j = 0; jpatch_area[j] = this->patch_area[j]*this->cell_size*this->cell_size; // multiplies the number of pixels of each patch by the area of a pixel + + + // Modification of the radius for the global density case + if(density_type==0) // Checks if the density type is set to zero (global) + { + raio = this->tamanho/sqrt(M_PI); // Changes the perception radius to the equivalent of a radius of a circle with the same area of a square with the lenght of the landscape + } + + // Inserts N individuals in the landscape through the populating() function + this->populating(raio,N,angulo_visada,passo,taxa_move,taxa_basal,taxa_morte,incl_b,incl_d,death_mat,move_mat,density_type,initialPosX, initialPosY,genotype_means,width_sds,points, Null); + + for(unsigned int i=0; ipopIndividuos.size(); i++) // Passes through each individuals + { + this->set_vizinhos(this->popIndividuos[i]); // Calls a function of the individual/individuo class to update the neighbours of a individual (the individuals within a radius distance of the focal individual) + this->atualiza_habitat(this->popIndividuos[i]); // Calls a function of the individual/individuo class to update the habitat type the individual is currently on + this->atualiza_patch(this->popIndividuos[i]); // Calls a function of the individual/individuo class to update the patch identfier of the patch he is currently on + } + + + for(unsigned int i=0; ipopIndividuos.size(); i++) // Passes through each individual + { + double dsty=this->calcDensity(popIndividuos[i]); // Calls a function of the individual/individuo class to determine the density of individuals within a radius distance area from a focal individual + this->popIndividuos[i]->update(dsty); // Calls a function of the individual/individuo class to update the characteristics of each individual + } + } -void paisagem::populating(double raio, int N, double angulo_visada, double passo, double move, double taxa_basal, - double taxa_morte, double incl_b, double incl_d, double death_m, - int dens_type) +// Function responsible for calling the constructor of the individual/individuo class to create N individuals at the start of the simulation +void paisagem::populating(double raio, + int N, + double angulo_visada, + double passo, + double taxa_move, + double taxa_basal, + double taxa_morte, + double incl_b, + double incl_d, + double death_m, + double move_m, + int dens_type, + double initialPosX[], + double initialPosY[], + double genotype_means[], + double width_sds[], + int points[], + bool Null) { - individuo::reset_id(); // reinicia o contador de id dos individuos - // Considerar diferentes possibilidades de posições iniciais. TBI. - if(this->initialPos==0) - { - for(int i=0; iN; i++) - { - this->popIndividuos.push_back(new individuo( - 0,//posicao x - 0,//posicao y - 0,//especie - taxa_morte,//taxa de morte - runif(0,360),// orientacao - angulo_visada,//angulo de visada - passo,//tamanho do passo - move, //taxa de movimentacao - raio,//tamanho do raio - taxa_basal,// taxa máxima de nascimento - 99, // semente de numero aleatorio - incl_b, - incl_d, - death_m, - dens_type)); - // como o popAgentes eh um ponteiro de vetores, ao adicionar enderecos das variaveis, usamos os new. Dessa forma fica mais rapido - //pois podemos acessar apenas o endereco e nao ficar guardando todos os valores - } - } - if(this->initialPos==1) // Random initial positions (initialPos==1) - { - for(int i=0; iN; i++) - { - this->popIndividuos.push_back(new individuo( - runif(this->tamanho/(-2),this->tamanho/2), //posicao x - runif(this->tamanho/(-2),this->tamanho/2), //posicao y - 0,//especie - taxa_morte,//taxa de morte - runif(0,360),// orientacao - angulo_visada,//angulo de visada - passo,//tamanho do passo - move, //taxa de movimentacao - raio,//tamanho do raio - taxa_basal,// taxa máxima de nascimento - 99, // semente de numero aleatorio - incl_b, - incl_d, - death_m, - dens_type)); - } + individuo::reset_id(); // Restars the id counter of the individuals + + vector genotype, width; // Creates vectors for storing either a genotype/wisght value (normal case), or all of them (null case) + + if (Null==FALSE) {// Checks if this is not a Null simulation run + + genotype.resize(1), // Resizes de genotype vectors for containing only a value + width.resize(1); // Resizes de width vectors for containing only a value + } + else{ // Checks if this a Null simulation run + + for (int i=0; iinitialPos==0) // Checks if the initialPos is set to 0 (origin) + { + for(int i=0; iN; i++) // Goess through the amount of initial individuals selected + { + if (Null==FALSE) { // Checks if this is not a Null simulation run + genotype[0] = genotype_means[i]; // Stores the genotype value of the current individual + width[0] = width_sds[i]; // Stores the width value of the current individual + } + + //OBS: As popAgents is a pointer of vectors, when adding variable addresses we use "new". This way should be faster, as we can acsess only the adress instead of keep storing the values + + // Creates an individual in the origin + this->popIndividuos.push_back(new individuo( + 0,//X Coordinate + 0,//Y Coordinate + 0,//Species ID + taxa_morte,//The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) + runif(0,360),// Angle that an individual is curently facing + angulo_visada,//Angle used for orientation when dispersing + passo,//The Lenght distance of a dispersal event + taxa_move, // The rate at which the individuals disperse + raio,// Density dependance radius + taxa_basal,// The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours) + 99, // Randomn seed + incl_b, //The slope of the birth density dependece function + incl_d, //The slope of the death density dependece function + death_m, //Constant that indicates how many times higher the death rate should be on non-habitat pixels + move_m, // Constant that indicates how many times lower the movement rate should be on non-habitat pixels + dens_type, //Density type (0 = global, 1 = local/within a individual radius) + genotype, // Individuals genotype trait mean, + width, // Individuals genotype trait width, + points[i])); // Number of sampling points + + } + } + if(this->initialPos==1) // Checks if the initialPos is set to 1 (Random initial positions) + { + for(int i=0; iN; i++) // Goess through the amount of initial individuals selected + { + if (Null==FALSE) { // Checks if this is not a Null simulation run + genotype[0] = genotype_means[i]; // Stores the genotype value of the current individual + width[0] = width_sds[i]; // Stores the width value of the current individual + } + + // Creates an individual at random locations + this->popIndividuos.push_back(new individuo( + runif(this->tamanho/(-2),this->tamanho/2), //X Coordinate + runif(this->tamanho/(-2),this->tamanho/2), //Y Coordinate + 0,//Species ID + taxa_morte,//The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) + runif(0,360),// Angle that an individual is curently facing + angulo_visada, //Angle used for orientation when dispersing, + passo,//The Lenght distance of a dispersal event + taxa_move, // The rate at which the individuals disperse + raio, // Density dependance radius + taxa_basal, // The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours) + 99, // Randomn seed + incl_b, //The slope of the birth density dependece function + incl_d, //The slope of the death density dependece function + death_m, //Constant that indicates how many times higher the death rate should be on non-habitat pixels + move_m, // Constant that indicates how many times lower the movement rate should be on non-habitat pixels + dens_type, //Density type (0 = global, 1 = local/within a individual radius) + genotype, // Individuals genotype trait mean, + width, // Individuals genotype trait width, + points[i])); // Number of sampling points + } + } + if(this->initialPos==2) // Checks if the initialPos is set to 2, Random initial positions with normal distribution. + { + for(int i=0; iN; i++) // Goess through the amount of initial individuals selected + { + if (Null==FALSE) { // Checks if this is not a Null simulation run + genotype[0] = genotype_means[i]; // Stores the genotype value of the current individual + width[0] = width_sds[i]; // Stores the width value of the current individual + } + + // Creates an individual at Random initial positions with normal distribution. + this->popIndividuos.push_back(new individuo( + rnorm(0,sqrt(taxa_move)*passo),//X Coordinate + rnorm(0,sqrt(taxa_move)*passo),//Y Coordinate + 0,//Species ID + taxa_morte,//The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) + runif(0,360),// Angle that an individual is curently facing + angulo_visada,//Angle used for orientation when dispersing + passo,//The Lenght distance of a dispersal event + taxa_move, // The rate at which the individuals disperse + raio,// Density dependance radius + taxa_basal,// The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours) + 99, // Randomn seed + incl_b, //The slope of the birth density dependece function + incl_d, //The slope of the death density dependece function + death_m, //Constant that indicates how many times higher the death rate should be on non-habitat pixels + move_m, // Constant that indicates how many times lower the movement rate should be on non-habitat pixels + dens_type, //Density type (0 = global, 1 = local/within a individual radius) + genotype, // Individuals genotype trait mean, + width, // Individuals genotype trait width, + points[i])); // Number of sampling points + } + } + + if(this->initialPos==3) // Checks if the initialPos is set to 3, Sets given initial positions to individuals + { + + for(int i=0; iN; i++) // Goess through the amount of initial individuals selected + { + if (Null==FALSE) { // Checks if this is not a Null simulation run + genotype[0] = genotype_means[i]; // Stores the genotype value of the current individual + width[0] = width_sds[i]; // Stores the width value of the current individual + } + + // Creates an individual at given coordinate. + this->popIndividuos.push_back(new individuo( + initialPosX[i],//X Coordinate + initialPosY[i],//Y Coordinatey + 0,//Species ID + taxa_morte,//The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) + runif(0,360), // Angle that an individual is curently facing + angulo_visada,//Angle used for orientation when dispersing + passo,//The Lenght distance of a dispersal event + taxa_move, // The rate at which the individuals disperse + raio,// Density dependance radius + taxa_basal,// The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours) + 99, // Randomn seed + incl_b, //The slope of the birth density dependece function + incl_d, //The slope of the death density dependece function + death_m, //Constant that indicates how many times higher the death rate should be on non-habitat pixels + move_m, // Constant that indicates how many times lower the movement rate should be on non-habitat pixels + dens_type, //Density type (0 = global, 1 = local/within a individual radius) + genotype, // Individuals genotype trait mean, + width, // Individuals genotype trait width, + points[i])); // Number of sampling points + } } - if(this->initialPos==2) // Random initial positions with normal distribution. TBI: tornar os parametros da rnorm livres - { - for(int i=0; iN; i++) - { - this->popIndividuos.push_back(new individuo( - rnorm(0,sqrt(move)*passo),//posicao x - rnorm(0,sqrt(move)*passo),//posicao y - 0,//especie - taxa_morte,//taxa de morte - runif(0,360),// orientacao - angulo_visada,//angulo de visada - passo,//tamanho do passo - move, //taxa de movimentacao - raio,//tamanho do raio - taxa_basal,// taxa máxima de nascimento - 99, // semente de numero aleatorio - incl_b, - incl_d, - death_m, - dens_type)); - } - } } -void paisagem::update() +// Updater function +void paisagem::update(int acao, int ind) { - if(this->popIndividuos.size()>0) + if(this->popIndividuos.size()>0) // checks if there are any currently alive individuals { - // Este for loop pode ser paralelizado, pois o que acontece com cada individuo eh independente - #ifdef PARALLEL - #pragma omp parallel for - #endif - for(unsigned int i=0; ipopIndividuos.size(); i++) + switch (acao) // Switches between possible actions (0= death, 1= birth, 2= dispersal) { - this->atualiza_vizinhos(this->popIndividuos[i]);//atualiza os vizinhos - this->atualiza_habitat(this->popIndividuos[i]);//retorna o tipo de habitat - this->atualiza_patch(this->popIndividuos[i]); + case 0: //0= death + break; + case 1: //1= birth + this->atualiza_habitat(this->popIndividuos[this->popIndividuos.size()-1]); // updates in the individual the environmental value of the pixel corresponent to its current coordinate + this->atualiza_patch(this->popIndividuos[this->popIndividuos.size()-1]); // updates the individual the patch identification value of the pixel corresponent to its current coordinate + break; + case 2: // 2= dispersal + this->atualiza_habitat(this->popIndividuos[ind]); // updates the individual the environmental value of the pixel corresponent to its current coordinate + this->atualiza_patch(this->popIndividuos[ind]); // updates the individual the patch identification value of the pixel corresponent to its current coordinate + break; } - // Este loop não é parelelizado, APESAR de ser independente, para garantir que as funcoes - // aleatorias sao chamadas sempre na mesma ordem (garante reprodutibilidade) - for(unsigned int i=0; ipopIndividuos.size(); i++) + + for(unsigned int i=0; ipopIndividuos.size(); i++) //Goes through all currently alive individuals { - double dsty=this->calcDensity(popIndividuos[i]); - this->popIndividuos[i]->update(dsty); //e atualiza o individuo i da populacao + double dsty=this->calcDensity(popIndividuos[i]); // returns the dentity of individuals within a radius distance area from each individual + this->popIndividuos[i]->update(dsty); // Updates the rates of each individual } - - } + } + } +// Function that drafts an individual int paisagem::sorteia_individuo() { - // time for next event and simulation time update - int menor=0; - double menor_tempo = this->popIndividuos[0]->get_tempo(); - - for(unsigned int i=1; ipopIndividuos.size(); i++) - { - if(this->popIndividuos[i]->get_tempo()popIndividuos[i]->get_tempo(); - } - } - return menor; + // time for next event and simulation time update + int menor=0; // Stores the first individual of the vector as the one with curently the least required time. + double menor_tempo = this->popIndividuos[0]->get_tempo(); // Drafts a time from a exponential equation used as the needed time to execute a action + + for(unsigned int i=1; ipopIndividuos.size(); i++) //Goess through all individuals + { + if(this->popIndividuos[i]->get_tempo()popIndividuos[i]->get_tempo(); // Stores the current lower time + } + } + return menor; // Returns the individual with the lowest time } +//Function that executes the selected action. It also returns a positive boolean value if the individual disperses out of the landscape boundary. bool paisagem::realiza_acao(int acao, int lower) //TODO : criar matriz de distancias como atributo do mundo e atualiza-la apenas quanto ao individuos afetado nesta funcao) { - bool emigrou = false; - switch(acao) //0 eh morte, 1 eh nascer, 2 eh andar + bool emigrou = false; // Initializes a boolean variable used to record if an individual has left the boundaries of the landscape as false + switch(acao) //(0= death, 1= birth, 2= dispersal) { - case 0: - delete this->popIndividuos[lower]; - this->popIndividuos.erase(this->popIndividuos.begin()+lower); - break; - - case 1: - individuo* chosen; - //Novo metodo para fazer copia do individuo: - chosen = new individuo(*this->popIndividuos[lower]); - this->popIndividuos.push_back(chosen); - break; - - case 2: - this->popIndividuos[lower]->anda(); - emigrou = this->apply_boundary(popIndividuos[lower]); - break; + case 0: //0= death + this->atualiza_vizinhos(acao, lower); + delete this->popIndividuos[lower]; // Deletes the object in the chosen position + this->popIndividuos.erase(this->popIndividuos.begin()+lower); // Erases the position that object occupied in the vector of individuals + break; + + case 1: // 1= birth + individuo* chosen; // Declares a new object of the individual class + //Novo metodo para fazer copia do individuo: + chosen = new individuo(*this->popIndividuos[lower]); // Calls a function to create a new individuals with similar characteristics to the parent one + this->popIndividuos.push_back(chosen); // Stores the new individual at the postition of the landscape individuals vector + this->atualiza_vizinhos(acao, lower); + break; + + case 2: // 2= birth + + emigrou = this->walk(lower); // Performs movement action and returs if the invidual left the landscape boundaries + + if(emigrou) //checks if an invidual left the landscape boundaries + this->realiza_acao(0, lower); //Treats an emigrating individual as if it died + else if(!emigrou) //checks if an invidual left the landscape boundaries + this->atualiza_vizinhos(acao, lower); //Apply the boundary condition if the individual disperses out of the landscape boundary Changes the boolean value to pisitive if a individual went ou of the landscape boundaries + break; } - return emigrou; + return emigrou; // Returs true if an invidual left the landscape boundaries } //para quando tiver especies, a definir... @@ -232,276 +368,448 @@ bool paisagem::realiza_acao(int acao, int lower) //TODO : criar matriz de distan //} -// metodo para condicao de contorno, argumento é um ponteiro para um individuo -//TODO: conferir se a combinacao x , y da condicao esta gerando o efeito desejado -//TBI: condicao periodica do codigo antigo feito com Garcia. Verificar se estah correta FE:Verifiquei e tinha um problema -// (veja p. ex. um unico individuo apenas se movimentando) +// Funtion to apply the desired boundary conditions when a dispersal event leaves an individual out of the landscape boundarys bool paisagem::apply_boundary(individuo * const ind) //const { - bool emigrou = false; - double rad = (double)ind->get_raio(); - switch(this->boundary_condition) - { - - case 0: - if(this->landscape_shape==0) - { - if(rad*rad < (double)ind->get_x()*(double)ind->get_x()+(double)ind->get_y()*(double)ind->get_y()) - { - for(unsigned int i=0; ipopIndividuos[i]->get_id()==(int)ind->get_id()) - { - delete this->popIndividuos[i]; - this->popIndividuos.erase(this->popIndividuos.begin()+i); - } - } - emigrou = true; - } - } - - if(this->landscape_shape==1) - { - if((double)ind->get_x()>=this->numb_cells*this->cell_size/2 || //>= porque na paisagem quadrado as bordas mais distantes de 0 iniciariam um proximo pixel que estaria fora da paisagem. Ou teriamos que assumir que esses pixels mais extremos tenha uma área maior, o que daria um trabalho adicional para implementar uma situação irreal. - (double)ind->get_x()<(this->numb_cells*this->cell_size/2)*(-1)|| - (double)ind->get_y()>this->numb_cells*this->cell_size/2 || - (double)ind->get_y()<=(this->numb_cells*this->cell_size/-2)) - { - for(unsigned int i=0; ipopIndividuos[i]->get_id()==(int)ind->get_id()) //DUVIDA: porque tem int? - { - delete this->popIndividuos[i]; - this->popIndividuos.erase(this->popIndividuos.begin()+i); - } - } - emigrou = true; - } - } - break; - - case 1: - if(ind->get_x()<(this->numb_cells*this->cell_size/2)*(-1)) - ind->set_x(this->tamanho+ind->get_x()); - if(ind->get_x()>=this->numb_cells*this->cell_size/2) - ind->set_x(ind->get_x()-this->tamanho); - if(ind->get_y()<(this->numb_cells*this->cell_size/2)*(-1)) - ind->set_y(this->tamanho+ind->get_y()); - if(ind->get_y()>=this->numb_cells*this->cell_size/2 ) - ind->set_y(ind->get_y()-this->tamanho); - break; - } - return emigrou; - /* TBI - case 2: reflexiva - */ + bool emigrou = false; // Initializes a boolean variable used to record if an individual has left the boundaries of the landscape as false + double rad = (double)ind->get_raio(); // Obtains the radius of the individual that dispersed + + switch(this->boundary_condition) // Switches between conditions (0= absortive, 1= periodical (pacman),2= reflexive) + { + + case 0: //0= absortive + if(this->landscape_shape==0) // Checks if the landscape shape is set to 0 (circle) + { + if(rad*rad < (double)ind->get_x()*(double)ind->get_x()+(double)ind->get_y()*(double)ind->get_y()) + // Checks if the individuals is out of bounds + { + emigrou = true; // Sets emigration to true + } + } + + if(this->landscape_shape==1) // Checks if the landscape shape is set to 1 (square) + { + if((double)ind->get_x()>=this->numb_cells*this->cell_size/2 || //Checks if the x coordinate is higher than the maximum value OR + // obs: >= porque na paisagem quadrado as bordas mais distantes de 0 iniciariam um proximo pixel que estaria fora da paisagem. Ou teriamos que assumir que esses pixels mais extremos tenha uma área maior, o que daria um trabalho adicional para implementar uma situação irreal. + (double)ind->get_x()<(this->numb_cells*this->cell_size/2)*(-1)|| // if the x coordinate is lower than the minimum value OR + (double)ind->get_y()>this->numb_cells*this->cell_size/2 || // if the y coordinate is higher than the maximum value Or + (double)ind->get_y()<=(this->numb_cells*this->cell_size/-2)) // if the x coordinate is lower than the minimum value + { + emigrou = true; // Sets emigration to true + } + } + break; + + case 1: // 1= periodical + if(ind->get_x()<(this->numb_cells*this->cell_size/2)*(-1)) //Checks if the x coordinate is higher than the maximum value + ind->set_x(this->tamanho+ind->get_x()); // Changes the x coordinate to represent the periodical condition (left corner to right corner) + if(ind->get_x()>=this->numb_cells*this->cell_size/2) //Checks if the x coordinate is lower than the minimum value + ind->set_x(ind->get_x()-this->tamanho); // Changes the x coordinate to represent the periodical condition (lright corner to left corner) + if(ind->get_y()<(this->numb_cells*this->cell_size/2)*(-1)) //Checks if the y coordinate is higher than the maximum value + ind->set_y(this->tamanho+ind->get_y()); // Changes the y coordinate to represent the periodical condition (bottom corner to top corner) + if(ind->get_y()>=this->numb_cells*this->cell_size/2 ) //Checks if the y coordinate is lower than the maximum value + ind->set_y(ind->get_y()-this->tamanho); // Changes the y coordinate to represent the periodical condition (up corner to bottom corner) + break; + + case 2: // 2= reflexive + if(ind->get_x() < -this->tamanho/2) //Checks if the x coordinate is lower than the minimum value + ind->set_x( -this->tamanho/2 + abs(this->tamanho/2 - ind->get_x()) ); // Sets the x coordinate to the boundary minus its excess lengh + if(ind->get_x() >= this->tamanho/2) //Checks if the x coordinate is higher than the maximum value + ind->set_x( this->tamanho/2 - abs(this->tamanho/2 - ind->get_x()) );// Sets the x coordinate to the boundary minus its excess lengh + if(ind->get_y() < -this->tamanho/2) //Checks if the y coordinate is lower than the minimum value + ind->set_y( -this->tamanho/2 + abs(this->tamanho/2 - ind->get_y()) );// Sets the y coordinate to the boundary minus its excess lengh + if(ind->get_y() >= this->tamanho/2) //Checks if the y coordinate is higher than the maximum value + ind->set_x( this->tamanho/2 - abs(this->tamanho/2 - ind->get_y()) );// Sets the y coordinate to the boundary minus its excess lengh + break; + } + return emigrou; // Returns a positive boolean value if the individual emigrated + } - +// Function that calculate the distance between two individuals double paisagem::calcDist(const individuo* a1, const individuo* a2) const //Virou método da paisagem pois não estava conseguindo usar as propriedades privadas da paisagem nessa função (p. ex. this->boundary_condition). Tive que tirar o primeiro const. { - switch(this->boundary_condition) - { - case 0: - return sqrt(((double)a1->get_x()-(double)a2->get_x())*((double)a1->get_x()-(double)a2->get_x())+((double)a1->get_y()-(double)a2->get_y())*((double)a1->get_y()-(double)a2->get_y())); - break; - - case 1: - double x1=a1->get_x(); - double x2=a2->get_x(); - double y1=a1->get_y(); - double y2=a2->get_y(); - double dx= x1 > x2 ? x1 - x2 : x2 - x1; // escolhe o menor entre x1-x2 e x2-x1 - double dy= y1 > y2 ? y1 - y2 : y2 - y1; - dx = dx > this->tamanho - dx ? this->tamanho - dx : dx; // escolhe o menor entre tam-dx e dx - dy = dy > this->tamanho - dy ? this->tamanho - dy : dy; - return sqrt(dx*dx + dy*dy); - break; - } + switch(this->boundary_condition) //switches between boundary conditoion cases (0= absortive, 1= periodical (pacman),2= reflexive) + { + case 0: // Euclidian Distance + return sqrt(((double)a1->get_x()-(double)a2->get_x())*((double)a1->get_x()-(double)a2->get_x())+((double)a1->get_y()-(double)a2->get_y())*((double)a1->get_y()-(double)a2->get_y())); + break; + + case 1: // The Euclidian distance acounting for the periodical effect (individuals on diferent edges may be close to one another) + double x1=a1->get_x(); // Gets the x coordinate of the first individuals + double x2=a2->get_x(); // Gets the x coordinate of the second individuals + double y1=a1->get_y(); // Gets the y coordinate of the first individuals + double y2=a2->get_y(); // Gets the y coordinate of the second individuals + double dx= x1 > x2 ? x1 - x2 : x2 - x1; // choses the lower between x1-x2 e x2-x1 + double dy= y1 > y2 ? y1 - y2 : y2 - y1; //choses the lower between y1-y2 e y2-y1 + dx = dx > this->tamanho - dx ? this->tamanho - dx : dx; // choses the lower between tam-dx e dx + dy = dy > this->tamanho - dy ? this->tamanho - dy : dy; //choses the lower between tam-dy e dy + return sqrt(dx*dx + dy*dy); // Computes and returns the distance + break; + } } // A function to calculate de density of individuals according to density type (global or local) and considering landscape boundary effects in the calculation. double paisagem::calcDensity(const individuo* ind1) const { - double density; - density=ind1->NBHood_size()/(M_PI*(ind1->get_raio()*ind1->get_raio())); - - // Functions for local density calculation - - /* 1. Circular area defining a region in which denso-dependence occurs: landscape boundary effects. - In this case, density is the number of individuals inside the circle divided by circle area. - This is the same calculation as for global density, except by the cases in which landscape boundary affects - the area of the circle. - */ - - // Condition giving the boundary effects cases - if(ind1->get_densType()==1) - { - if(ind1->get_x()*ind1->get_x()>((this->tamanho/2)-ind1->get_raio())*((this->tamanho/2)-ind1->get_raio()) || - ind1->get_y()*ind1->get_y()>((this->tamanho/2)-ind1->get_raio())*((this->tamanho/2)-ind1->get_raio())) - { - // temporary objects - double modIx = fabs(ind1->get_x()); - double modIy = fabs(ind1->get_y()); - double XYmax = this->tamanho/2; - vectorsecX; - vectorsecY; - - // Functions for adjusted local density calculation, according to the specific case - // 1) - if(modIx>XYmax-ind1->get_raio() && modIy<=XYmax-ind1->get_raio()) - { - secX.push_back(XYmax); - secX.push_back(XYmax); - secY.push_back(modIy+sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); - secY.push_back(modIy-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); - - double distSecs = secY[0]-secY[1]; - double height = XYmax - modIx; - double theta = acos(1-(distSecs*distSecs/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians - double adjArea = M_PI*ind1->get_raio()*ind1->get_raio() - (theta*ind1->get_raio()*ind1->get_raio()/2 - (distSecs*height/2)); - - density = ind1->NBHood_size()/adjArea; - - } - - // 2) - if(modIx<=XYmax-ind1->get_raio() && modIy>XYmax-ind1->get_raio()) - { - secY.push_back(XYmax); - secY.push_back(XYmax); - secX.push_back(modIx+sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); - secX.push_back(modIx-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); - - double distSecs = secX[0]-secX[1]; - double height = XYmax - modIy; - double theta = acos(1-(distSecs*distSecs/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians - double adjArea = M_PI*ind1->get_raio()*ind1->get_raio() - (theta*ind1->get_raio()*ind1->get_raio()/2 - (distSecs*height/2)); - - density = ind1->NBHood_size()/adjArea; - } - - - if(modIx>XYmax-ind1->get_raio() && modIy>XYmax-ind1->get_raio()) - { - - // 3) - if((modIx-XYmax)*(modIx-XYmax)+(modIy-XYmax)*(modIy-XYmax)>ind1->get_raio()*ind1->get_raio()) - { - secX.push_back(modIx+sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); - secY.push_back(XYmax); - secX.push_back(modIx-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); - secY.push_back(XYmax); - secX.push_back(XYmax); - secY.push_back(modIy+sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); - secX.push_back(XYmax); - secY.push_back(modIy-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); - - double distSecs = sqrt((secX[3]-secX[1])*(secX[3]-secX[1])+(secY[1]-secY[3])*(secY[1]-secY[3])); - double distSecs2 = sqrt((secX[2]-secX[0])*(secX[2]-secX[0])+(secY[0]-secY[2])*(secY[0]-secY[2])); - double theta = acos(1-(distSecs*distSecs/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians - double phi = acos(1-(distSecs2*distSecs2/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians - double adjArea = (2*M_PI-theta)*ind1->get_raio()*ind1->get_raio()/2 + phi*ind1->get_raio()*ind1->get_raio()/2 + (secX[0]-secX[1])*(XYmax-modIy)/2 + (secY[2]-secY[3])*(XYmax-modIx)/2; - - density = ind1->NBHood_size()/adjArea; - - } - // 4) - else - { - secX.push_back(modIx-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); - secY.push_back(XYmax); - secX.push_back(XYmax); - secY.push_back(modIy-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); - - double distSecs = sqrt((secX[1]-secX[0])*(secX[1]-secX[0])+(secY[0]-secY[1])*(secY[0]-secY[1])); - double theta = acos(1-(distSecs*distSecs/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians - double adjArea = theta*ind1->get_raio()*ind1->get_raio()/2 + (XYmax-secX[0])*(XYmax-modIy) + (XYmax-secY[1])*(XYmax-modIx); - - density = ind1->NBHood_size()/adjArea; - } - - } - } - } - - /* 2. Density kernel (TBI). - - if(ind1->get_densType()==2) {} - - */ - return density; + double density; // Creates a temporary variable for storing density + + if((M_PI*(ind1->get_raio()*ind1->get_raio()))>(this->tamanho*this->tamanho)){ // check if the individual radius would be higher than the landscape + + density=ind1->NBHood_size()/(this->tamanho*this->tamanho); // Computes the standard case of density using the total area od the world + } + + else{ + density=ind1->NBHood_size()/(M_PI*(ind1->get_raio()*ind1->get_raio())); // Computes the standard case of density (when an individual radius area is completely within the landscape) + } + + // Functions for local density calculation + + /* 1. Circular area defining a region in which denso-dependence occurs: landscape boundary effects. + In this case, density is the number of individuals inside the circle divided by circle area. + This is the same calculation as for global density, except by the cases in which landscape boundary affects + the area of the circle. + */ + + // Condition giving the boundary effects cases + if(ind1->get_densType()==1) // checks if the density type is set to 1 (local/ within a radius distance) + { + if(ind1->get_x()*ind1->get_x()>((this->tamanho/2)-ind1->get_raio())*((this->tamanho/2)-ind1->get_raio()) || // Checks if the x coordinate is within a dadius distance from the edge OR + ind1->get_y()*ind1->get_y()>((this->tamanho/2)-ind1->get_raio())*((this->tamanho/2)-ind1->get_raio())) // if the y coordinate is within a dadius distance from the edge + { // OBS: The absolute values retain their distance from the edge, so the compuatations from all quadrants can be done with a single set of equations + + // temporary objects + double modIx = fabs(ind1->get_x()); // Stores the the absolute value of x + double modIy = fabs(ind1->get_y()); // Stores the the absolute value of y + double XYmax = this->tamanho/2; // Stores the lenght of quadrant side + vectorsecX; + vectorsecY; + + // Functions for adjusted local density calculation, according to the specific case + // 1) + if(modIx>XYmax-ind1->get_raio() && modIy<=XYmax-ind1->get_raio()) // Checks if the x coordinate is within a dadius distance from the edge AND if the y coordinate is not + { + secX.push_back(XYmax); + secX.push_back(XYmax); + secY.push_back(modIy+sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); + secY.push_back(modIy-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); + + double distSecs = secY[0]-secY[1]; + double height = XYmax - modIx; + double theta = acos(1-(distSecs*distSecs/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians + double adjArea = M_PI*ind1->get_raio()*ind1->get_raio() - (theta*ind1->get_raio()*ind1->get_raio()/2 - (distSecs*height/2)); + + density = ind1->NBHood_size()/adjArea; + + } + + // 2) + if(modIx<=XYmax-ind1->get_raio() && modIy>XYmax-ind1->get_raio()) // Checks if the y coordinate is within a dadius distance from the edge AND if the x coordinate is not + { + secY.push_back(XYmax); + secY.push_back(XYmax); + secX.push_back(modIx+sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); + secX.push_back(modIx-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); + + double distSecs = secX[0]-secX[1]; + double height = XYmax - modIy; + double theta = acos(1-(distSecs*distSecs/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians + double adjArea = M_PI*ind1->get_raio()*ind1->get_raio() - (theta*ind1->get_raio()*ind1->get_raio()/2 - (distSecs*height/2)); + + density = ind1->NBHood_size()/adjArea; + } + + + if(modIx>XYmax-ind1->get_raio() && modIy>XYmax-ind1->get_raio()) + { + + // 3) + if((modIx-XYmax)*(modIx-XYmax)+(modIy-XYmax)*(modIy-XYmax)>ind1->get_raio()*ind1->get_raio()) + { + secX.push_back(modIx+sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); + secY.push_back(XYmax); + secX.push_back(modIx-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); + secY.push_back(XYmax); + secX.push_back(XYmax); + secY.push_back(modIy+sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); + secX.push_back(XYmax); + secY.push_back(modIy-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); + + double distSecs = sqrt((secX[3]-secX[1])*(secX[3]-secX[1])+(secY[1]-secY[3])*(secY[1]-secY[3])); + double distSecs2 = sqrt((secX[2]-secX[0])*(secX[2]-secX[0])+(secY[0]-secY[2])*(secY[0]-secY[2])); + double theta = acos(1-(distSecs*distSecs/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians + double phi = acos(1-(distSecs2*distSecs2/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians + double adjArea = (2*M_PI-theta)*ind1->get_raio()*ind1->get_raio()/2 + phi*ind1->get_raio()*ind1->get_raio()/2 + (secX[0]-secX[1])*(XYmax-modIy)/2 + (secY[2]-secY[3])*(XYmax-modIx)/2; + + density = ind1->NBHood_size()/adjArea; + + } + // 4) + else + { + secX.push_back(modIx-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIy)*(XYmax-modIy)))); + secY.push_back(XYmax); + secX.push_back(XYmax); + secY.push_back(modIy-sqrt(ind1->get_raio()*ind1->get_raio()-((XYmax-modIx)*(XYmax-modIx)))); + + double distSecs = sqrt((secX[1]-secX[0])*(secX[1]-secX[0])+(secY[0]-secY[1])*(secY[0]-secY[1])); + double theta = acos(1-(distSecs*distSecs/(2*ind1->get_raio()*ind1->get_raio()))); // angle in radians + double adjArea = theta*ind1->get_raio()*ind1->get_raio()/2 + (XYmax-secX[0])*(XYmax-modIy) + (XYmax-secY[1])*(XYmax-modIx); + + density = ind1->NBHood_size()/adjArea; + } + + } + } + } + + /* 2. Density kernel (TBI). + + if(ind1->get_densType()==2) {} + + */ + return density; // Returns the individuals curent density } /* - Sempre adicione const aos argumentos de métodos quando o método não - deve alterá-los. Previne vários erros e pode otimizar compilação -*/ + Sempre adicione const aos argumentos de métodos quando o método não + deve alterá-los. Previne vários erros e pode otimizar compilação + */ -void paisagem::atualiza_vizinhos(individuo * const ag1) const //acessando os vizinhos dos agentes +// Function that updates an individuals neighbourhood list +void paisagem::set_vizinhos(individuo * const ag1) const //acessando os vizinhos dos agentes { - vector listViz; - if(ag1->get_densType()==0) //dens_type poderia voltar como propriedade da paisagem. Facilitariam as coisas. Como muitas propriedades e métodos deste código, elas podem ser interpretadas das duas formas (como do individuo ou como da paisagem). O que está dando confusão é que estamos fazendo um IBM, mas para algumas situações estamos querendo simular dinâmicas cujas variáveis de interesse são propriedades populacionais e não do indivíduo. Se aceito, limar o método get_densType() do individuo.h. - { - for(unsigned int j=0; jpopIndividuos[j]; - if(ag1==ag2) continue; - listViz.push_back(ag2); - } - } - else - { - double rad = (double)ag1->get_raio(); - for(unsigned int j=0; jpopIndividuos[j]; - if(ag1==ag2) continue; - double d=this->calcDist(ag1,ag2); - if(d<=rad) {listViz.push_back(ag2);} - } - } - ag1->set_vizinhos(listViz); + vector listViz; // Creates a vector of individuals to store the neighbours + if(ag1->get_densType()==0) // Checks if the density type is set to 0 (global) + //dens_type poderia voltar como propriedade da paisagem. Facilitariam as coisas. Como muitas propriedades e métodos deste código, elas podem ser interpretadas das duas formas (como do individuo ou como da paisagem). O que está dando confusão é que estamos fazendo um IBM, mas para algumas situações estamos querendo simular dinâmicas cujas variáveis de interesse são propriedades populacionais e não do indivíduo. Se aceito, limar o método get_densType() do individuo.h. + { + for(unsigned int j=0; jpopIndividuos[j]; //Stores a copy of each individual + if(ag1==ag2) continue; //Except if it is the focal individual itself + listViz.push_back(ag2); // Stores the copy of each other individual in a neighbours vector + } + } + else //Checks if the density type is not set to 0 (1= local) + { + double rad = (double)ag1->get_raio(); // Gets the radius of the focal individual + for(unsigned int j=0; jpopIndividuos[j]; //Stores a copy of each individual + if(ag1==ag2) continue; //Except if it is the focal individual itself + double d=this->calcDist(ag1,ag2); // Checks if the copyed individual is within a radius distance + if(d<=rad) {listViz.push_back(ag2);} // Stores the copy of each other individual within a radius distance into a neighbours vector + } + } + ag1->set_vizinhos(listViz); // Passes to an individual its neighbourhood list + +} +// Function for updating individuals neighbours depending on their last action +void paisagem::atualiza_vizinhos(int acao, int ind) const //acessando os vizinhos dos agentes +{ + vector listViz; // Temporary object for string neighbours + switch(acao) // Switches between actions (0= death, 1= birth, 2= dispersal) + { + case 0: // Performs the neighbourhood update in case of a death event + for(unsigned int i=0; ipopIndividuos[ind]->NBHood_size()-1; i++) // Passes through each of the drafted individuals neighbours + { + for(unsigned int j=0; jpopIndividuos[ind]->lisViz[i]->NBHood_size()-1; j++)// Passes through each of the drafted individuals neighbours neighbours + { + if(this->popIndividuos[ind]->lisViz[i]->lisViz[j]->get_id() == this->popIndividuos[ind]->get_id()) // Checks if the current neighbours neighbour has the id of the drafted individual + { + this->popIndividuos[ind]->lisViz[i]->lisViz.erase(this->popIndividuos[ind]->lisViz[i]->lisViz.begin()+j); // Erases the drafted individual from the neighbourhood of its neighbours + } + } + } + break; // Breaks switch + + case 1:// Performs the neighbourhood update in case of a birth event + for(unsigned int i=0; ipopIndividuos[ind]->NBHood_size()-1; i++)// Passes through each of the drafted individuals neighbours + { + this->popIndividuos[ind]->lisViz[i]->lisViz.push_back(this->popIndividuos[this->popIndividuos.size()-1]);// Adds the last born individual to its neighbours neighbourhood + listViz.push_back(this->popIndividuos[ind]->lisViz[i]);// Stores the drafted individuals neighbours + } + + this->popIndividuos[ind]->lisViz.push_back(this->popIndividuos[this->popIndividuos.size()-1]);// Adds the last born individual to its parent neighbourhood + listViz.push_back(this->popIndividuos[ind]);// Adds the last born individual to its own neighbourhood + this->popIndividuos[this->popIndividuos.size()-1]->set_vizinhos(listViz);// Calls a funtion that stores the last born individuals neighbours in its neighbourhood + break; // Breaks switch + + case 2:// Performs the neighbourhood update in case of a migration event + + if(this->popIndividuos[ind]->get_densType()!=0){// Checks if the density type is not set to global + + for(unsigned int i=0; ipopIndividuos[ind]->NBHood_size()-1; i++) // Passes through each of the drafted individuals neighbours + { + for(unsigned int j=0; jpopIndividuos[ind]->lisViz[i]->NBHood_size()-1; j++)// Passes through each of the drafted individuals neighbours neighbours + { + if(this->popIndividuos[ind]->lisViz[i]->lisViz[j]->get_id() == this->popIndividuos[ind]->get_id()) // Checks if the current neighbours neighbour has the id of the drafted individual + { + this->popIndividuos[ind]->lisViz[i]->lisViz.erase(this->popIndividuos[ind]->lisViz[i]->lisViz.begin()+j); // Erases the drafted individual from the neighbourhood of its neighbours + } + } + } + + double rad = (double) this->popIndividuos[ind]->get_raio(); // Obtains the drafted individuals radius + for(unsigned int k=0; kpopIndividuos.size(); k++) // Passes through each living individual + { + if(k==ind) // Checks if the compared individual isnt itself + continue; + if(this->calcDist(this->popIndividuos[ind], this->popIndividuos[k]) <= rad) // Checks if the compared individual is within a radius distance + { + listViz.push_back(this->popIndividuos[k]); // Stores the drafted individuals neighbours + this->popIndividuos[k]->lisViz.push_back(this->popIndividuos[ind]); // Stores the drafted individual on its neighbours neighbourhood + } + } + this->popIndividuos[ind]->set_vizinhos(listViz); // Calls a funtion that stores the neighbours in the individuals neighbourhood + break; // Breaks switch + + } + } } +// Function that updates if an individual is currently on habitat or matrix void paisagem::atualiza_habitat(individuo * const ind) const { - // Tinha um IF com landscape_shape que eu removi. Não entendi como a paisagem ser circular - // interfere em ser habitat ou não: isso deve interferir na apply_boundary apenas, certo? - // Also: Tinha uma inversão do y que eu também não entendi e removi - // A.C. 10.07.13 - - // Um termo (-1) foi removido erroneamente por A.C.. Para o hy, o sentido em que o número de células aumenta é o - //contrário do sentido em que as coordenadas aumentam. Portanto a multiplicação por - 1 é necessária. - // M.A. 12.09.14 - int hx,hy; - hx= (double)ind->get_x()/this->cell_size+this->numb_cells/2; - hy= ((double)ind->get_y()/this->cell_size)*(-1)+this->numb_cells/2; - ind->set_habitat(this->landscape[hx][hy]); + // Tinha um IF com landscape_shape que eu removi. Não entendi como a paisagem ser circular + // interfere em ser habitat ou não: isso deve interferir na apply_boundary apenas, certo? + // Also: Tinha uma inversão do y que eu também não entendi e removi + // A.C. 10.07.13 + + // Um termo (-1) foi removido erroneamente por A.C.. Para o hy, o sentido em que o número de células aumenta é o + //contrário do sentido em que as coordenadas aumentam. Portanto a multiplicação por - 1 é necessária. + // M.A. 12.09.14 + + int hx,hy; //Initializes variables for storing pixels coordinats + hx= (double)ind->get_x()/this->cell_size+this->numb_cells/2; // Finds the pixel X coodenate + hy= ((double)ind->get_y()/this->cell_size)*(-1)+this->numb_cells/2; //Finds the pixel y coodenate + ind->set_habitat(this->landscape[hx][hy]); // Sets the patch label of the individual to the label of the patch it is currently on } +// Function to update the patch the individual is currently in void paisagem::atualiza_patch(individuo * const ind) const { - int hx,hy; - hx= (double)ind->get_x()/this->cell_size+this->numb_cells/2; - hy= ((double)ind->get_y()/this->cell_size)*(-1)+this->numb_cells/2; - ind->set_patch(this->patches[hx][hy]); + int hx,hy; //Initializes variables for storing pixels coordinats + hx= (double)ind->get_x()/this->cell_size+this->numb_cells/2; // Finds the pixel X coodenate + hy= ((double)ind->get_y()/this->cell_size)*(-1)+this->numb_cells/2; //Finds the pixel y coodenate + ind->set_patch(this->patches[hx][hy]); // Sets the patch label of the individual to the label of the patch it is currently on } +// Retursive function that finds and identfy each landscape path void paisagem::find_patches(int x, int y, int current_label) { - if (x < 0 || x == this->numb_cells) return; // out of bounds - if (y < 0 || y == this->numb_cells) return; // out of bounds - if (this->patches[x][y] || !this->landscape[x][y]) return; // already labeled or not marked with 1 in m - - // mark the current cell - this->patches[x][y] = current_label; + if (x < 0 || x == this->numb_cells) return; // Checks if the x coordinate is out of bounds + if (y < 0 || y == this->numb_cells) return; // Checks if the y coordinate is out of bounds + if (this->patches[x][y] || !this->landscape[x][y]) return; // Checks if the curret pixel is already labeled AND if is a habitat pixel + + // marks the current cell with its respective patch ID value + this->patches[x][y] = current_label; + + // recursively mark the neighbors + find_patches(x + 1, y, current_label); + find_patches(x, y + 1, current_label); + find_patches(x - 1, y, current_label); + find_patches(x, y - 1, current_label); + + // 8-rule + find_patches(x + 1, y + 1, current_label); + find_patches(x - 1, y + 1, current_label); + find_patches(x - 1, y - 1, current_label); + find_patches(x + 1, y - 1, current_label); +} - // recursively mark the neighbors - find_patches(x + 1, y, current_label); - find_patches(x, y + 1, current_label); - find_patches(x - 1, y, current_label); - find_patches(x, y - 1, current_label); +// Function that decides upon, and calls other functions to perform, the desired method of dispersion (Random walk or habitat selection). In the habitat selection case this function also samples x Points within the individuals "passo" radius and the indivuduals relative fitness on that location. +bool paisagem::walk(int lower){ + + bool emigrou =FALSE; // Initilizes boolean variable that mars if the indiviual emigrated + + if(this->popIndividuos[lower]->get_points()>0) // Checks if the the individuals sample their step radius + { + double possibilitities[this->popIndividuos[lower]->get_points()][3]; // Creates a vector for storing the possible migration locations + + double choice, dist; // Temporary variables for storing angle and lenght + int hx,hy; // Temporary variables for storing the possible coordinates + + for (int i=0; ipopIndividuos[lower]->get_points(); i++) // Passes by each possible location + { + + choice=runif(0.0,360.0); // Samples an angle + dist=runif(0.0,this->popIndividuos[lower]->get_passo()); // Samples a distance + + possibilitities[i][0]=this->popIndividuos[lower]->get_x()+cos(choice)*dist; // Stores the possible x + possibilitities[i][1]=this->popIndividuos[lower]->get_y()+sin(choice)*dist; // Stores the possible y + + + + switch(this->boundary_condition){ // switches between baoundary conditions (0= absortive, 1= periodical (pacman),2= reflexive) + + case 0: // 0= absortive + + if(possibilitities[i][0]<(this->numb_cells*this->cell_size/2)*(-1) || //Checks if the x coordinate is higher than the maximum value OR + possibilitities[i][0]>=this->numb_cells*this->cell_size/2 || // if the x coordinate is lower than the minimum value OR + possibilitities[i][1]<(this->numb_cells*this->cell_size/2)*(-1) || // if the y coordinate is higher than the maximum value Or + possibilitities[i][1]>=this->numb_cells*this->cell_size/2)// if the x coordinate is lower than the minimum value + { + + possibilitities[i][2] = -10000; // Sets the out of boundaries locations rank to a very small number + } + + case 1: // 1= periodical (pacman) + if(possibilitities[i][0]<(this->numb_cells*this->cell_size/2)*(-1)) //Checks if the x coordinate is higher than the maximum value OR + possibilitities[i][0]=(this->tamanho+possibilitities[i][0]); // Transtports coordinate to the other edge + if(possibilitities[i][0]>=this->numb_cells*this->cell_size/2) // if the x coordinate is lower than the minimum value OR + possibilitities[i][0]=(possibilitities[i][0]-this->tamanho); // Transtports coordinate to the other edge + if(possibilitities[i][1]<(this->numb_cells*this->cell_size/2)*(-1))// if the y coordinate is higher than the maximum value Or + possibilitities[i][1]=(this->tamanho+possibilitities[i][1]); // Transtports coordinate to the other edge + if(possibilitities[i][1]>=this->numb_cells*this->cell_size/2 ) // if the x coordinate is lower than the minimum value + possibilitities[i][1]=(possibilitities[i][1]-this->tamanho); // Transtports coordinate to the other edge + + hx= (double)possibilitities[i][0]/this->cell_size+this->numb_cells/2; // Discovers the x pixel index of the coordinate + hy= ((double)possibilitities[i][1]/this->cell_size)*(-1)+this->numb_cells/2; // Discovers the y pixel index of the coordinate + + possibilitities[i][2] = this->landscape[hx][hy]; // Stores the habitat value of the current pixel + + break; + + case 2:// 2= reflexive + // Prevents out of boundary corrdinates from being sampled + while (possibilitities[i][0]<(this->numb_cells*this->cell_size/2)*(-1) || //Checks if the x coordinate is higher than the maximum value OR + possibilitities[i][0]>=this->numb_cells*this->cell_size/2 || // if the x coordinate is lower than the minimum value OR + possibilitities[i][1]<(this->numb_cells*this->cell_size/2)*(-1) || // if the y coordinate is higher than the maximum value Or + possibilitities[i][1]>=this->numb_cells*this->cell_size/2) { // if the x coordinate is lower than the minimum value + + choice=runif(0.0,360.0); // re-Samples an angle + dist=runif(0.0,this->popIndividuos[lower]->get_passo()); // re-Samples a distance + + possibilitities[i][0]=this->popIndividuos[lower]->get_x()+cos(choice)*dist;// Stores the possible x + possibilitities[i][1]=this->popIndividuos[lower]->get_y()+sin(choice)*dist;// Stores the possible y + } + hx= (double)possibilitities[i][0]/this->cell_size+this->numb_cells/2; // Discovers the x pixel index of the coordinate + hy= ((double)possibilitities[i][1]/this->cell_size)*(-1)+this->numb_cells/2; // Discovers the y pixel index of the coordinate + + possibilitities[i][2] = this->landscape[hx][hy]; // Stores the habitat value of the current pixel + + break; + } + + } + + this->popIndividuos[lower]->habitat_selection(possibilitities);// Passes the possible locations for the individual to chose from + + if (this->boundary_condition==0) { + emigrou = this->apply_boundary(popIndividuos[lower]); + } + + } + + else // Checks if the the individuals dont sample their step radius (ramdom) + { + this->popIndividuos[lower]->anda(); // Calls a function that changes the X and Y coordinates of the individuals via randomwalk + emigrou = this->apply_boundary(popIndividuos[lower]); // Applyes boundary condition + } + return emigrou; } + + diff --git a/paisagem.hpp b/paisagem.hpp index 13b5eb0..0020e9f 100644 --- a/paisagem.hpp +++ b/paisagem.hpp @@ -5,160 +5,240 @@ #include #include #include -///** Dimensão maxima do mapa, usada para rotinas de fragmentação TBI */ +///** Maximum landscape side lenght, used in fragmentation routines TBI */ #define dim 10000 /*Em algum momento pode ser alterado para um argumento do - construtor. No momento não é prioritário. */ +construtor. No momento não é prioritário. */ //Isto aqui não está aparecendo no doxygen using namespace std; /** \brief A classe paisagem implementa o mundo onde os agentes estão. * - * Esta classe é responsável pela criação dos indivíduos, por manter o relógio do mundo, e por estabelecer a comunicação entre os indivíduos \sa \ref paisagem::update */ + * The paisagem/landscape class is responsible for generating the domain were the agents(individuals) reside + * This class is responsable for creating individuals, updating the world clock, and stablishes comunication between individuals + * \sa \ref paisagem::update */ + class paisagem { private: - - //propriedades privadas - /** Tamanho do lado da paisagem */ + + // Private properties + /** Lenght of a square landscape side */ const double tamanho; - /** Número de indivíduos no começo da simulação */ - const unsigned long N; - /** Vetor dos indivíduos da populacao */ - vector popIndividuos; - /** Número de pixels do lado da paisagem */ + /** Number of individuals at the start of the simulation */ + const unsigned long N; + /** Vector for storing the individuals of the population */ + vector popIndividuos; + /** Number of pixels in a square landscape side */ const int numb_cells; - /** Tamanho do lada de cada pixel da imagem usada para a classificação da paisagem */ + /** Resolution Lenght of a square pixel side */ const double cell_size; - /** Forma da paisagem (0 = circular, 1 = quadrada)*/ - const int landscape_shape; - /** Tipo de condição de contorno (0 = absortiva, 1 = periódica, 2 = reflexiva) */ - const int boundary_condition; - /** Matriz de pixels da paisagem, assume os valores 0 quando o pixel está na matriz e 1 quando está no habitat */ - int landscape[dim][dim];//[linha][coluna] temporariamente substituido or valor fixo - /** Matriz com a determinação do fragmento a que cada pixel pertence (0 para matriz; 1, 2, 3, ... para fragmentos */ - int patches[dim][dim]; - /** Número de fragmentos da paisagem, desconsiderando-se a matriz. */ - int numb_patches; - /** Ponteiro para o vetor contendo a área de cada fragmento */ - double* patch_area; - /** Posição dos indivíduos no início da simulação (0 = origem; 1 = aleatória com distribuição uniforme na paisagem; 2 = aleatória com distribuição normal na paisagem)*/ - const int initialPos; - - //metodos privados - void populating( - /** Raio no qual os indivíduos percebem vizinhos. Deve ser menor que o tamanho do mundo */ - const double raio, - /** Número de indivíduos no começo da simulação */ - const int N, - /** Ângulo de visada dos indivíduos da população */ - const double angulo_visada, - /** Passo de caminhada dos indivíduos da população */ - const double passo, - /** Taxa de movimentação dos individuos da população */ - const double move, - /** Taxa de nascimento de um indivíduo no habitat favorável e sem vizinhos */ - const double taxa_basal, - /** Taxa de mortalidade dos indivíduos */ - const double taxa_morte, - /** Inclinação da curva de denso-depedência para natalidade */ - const double incl_b, - /** Inclinação da curva de denso-depedência para mortalidade */ - const double incl_d, - /** Constante que indica quantas vezes a mortalidade basal na matriz é maior que no habitat */ - const double death_m, - /** Tipo de densidade (0 = GLOBAL, 1 = LOCAL) */ - const int dens_type - ); - - /** Atualiza a lista de vizinhos de um indivíduo */ - void atualiza_vizinhos(individuo * const ind) const; - /** Informa o indivíduo o tipo de habitat (matriz ou habitat) correspondente à sua atual posição, atualizando \ref individuo::tipo_habitat */ + /** Shape of the landscape (0= circle, 1= square)*/ + const int landscape_shape; + /** The boundary condition type affects how individuals interact with the edges of the landscape (0= absortive, 1= periodical (pacman),2= reflexive) */ + const int boundary_condition; + /** Matrix containing the evironmental values of the landscape pixels (0= matrix, 1= habitat)*/ + double landscape[dim][dim]; + /** Matriz coatining the fragment ID that each pixel resides (0 for matriz; 1, 2, 3, ... for fragments) */ + int patches[dim][dim]; + /** number of non-matrix patches */ + int numb_patches; + /** Pointer for a vector storing the area of each patch */ + double* patch_area; + /** The initial postion of individuals (0 = origin, 1 = random, 2 = normaly distributed with mean on origin, 3 = inputed initial coordinates)*/ + const int initialPos; + + //Private methods + + // Function responsible for calling the constructor of the individual/individuo class to create N individuals at the start of the simulation + + void populating( + /** Radius of density dependant influence */ + const double raio, + /** Number of individuals at the start of the simulation */ + const int N, + /** Angle used for orientation when dispersing */ + const double angulo_visada, + /** The Lenght distance of a dispersal event (random walk) or maximum dispersal radius (Selection step) */ + const double passo, + /** The rate at which the individuals disperse */ + const double taxa_move, + /** The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours) */ + const double taxa_basal, + /** The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) */ + const double taxa_morte, + /** The slope of the birth density dependance function */ + const double incl_b, + /** The slope of the death density dependance function*/ + const double incl_d, + /** Constant that indicates how many times higher the death rate should be on non-habitat pixels */ + const double death_m, + /** Constant that indicates how many times lower the movement rate should be on non-habitat pixels */ + const double move_m, + /** Density type (0 = global, 1 = local/within a individual radius) */ + const int dens_type, + /** Vector containing the x coordinates initial individuals */ + double initialPosX[], + /** Vector containing the y coordinates initial individuals */ + double initialPosY[], + /** Vector containing the genotypical trait means of the initial individuals */ + double genotype_means[], + /** Vector containing the standard deviations of the initial individuals */ + double width_sds[], + /** Vector containing the numer of points initial individuals sample when selecting habitats */ + int points[], + /** Booling value switching simuluations to neutral state (all individuals acting as an average individual) */ + bool Null); + + /** Creates an individuals neighbourhood list (each of the individuals within a radius distance of the focal individual), or all of them (global) + Paran individuo * const ind - an object of the individual/individuo class */ + void set_vizinhos(individuo * const ind) const; + + /** Updates an individuals neighbourhood list (each of the individuals within a radius distance of the focal individual), or all of them (global) + Paran individuo * const ind - an object of the individual/individuo class + Paran int acao - an variable containing the action performed by the individual */ + void atualiza_vizinhos(int acao, int ind) const; + + /** Updates in the individual the environmental value of the pixel corresponent to its current coordinate (if binary: 0= matrix, 1= habitat) + Paran individuo * const ind - an object of the individual/individuo class + \ref individuo::tipo_habitat */ void atualiza_habitat(individuo * const ind) const; - /** Informa o indivíduo o fragmento correspondente à sua atual posição, atualizando \ref individuo::patch_label */ + + /** Updates in the individual the patch identification value of the pixel corresponent to its current coordinate + Paran individuo * const ind - an object of the individual/individuo class + \ref individuo::patch_label */ void atualiza_patch(individuo * const ind) const; - /** Aplica a condição de contorno após a movimentação */ - bool apply_boundary(individuo * const ind); //const; // metodo para aplicação da condicao de contorno - - + + /** Aplies the boundary condition after a dispersal event + Paran individuo * const ind - an object of the individual/individuo class */ + bool apply_boundary(individuo * const ind); + + public: - //vector popIndividuos; - - /** Contador de quanto tempo já transcorreu no mundo simulado */ + + /** The world counter used for storing how much time has already passed */ double tempo_do_mundo; - - //metodos publicos - /** Construtor da classe paisagem */ + + //Public methods + /** Constructor of the landscape/paisagem class */ paisagem( - /** Raio no qual os indivíduos percebem vizinhos. Deve ser menor que o tamanho do mundo */ - const double raio, - /** Número de indivíduos no começo da simulação */ - const int N, - /** Ângulo de visada dos indivíduos */ - const double angulo_visada, - /** Passo de caminhada dos indivíduos */ - const double passo, - /** Taxa de movimentação dos indivíduos */ - const double move, - /** Taxa de nascimento de um indivíduo no habitat favorável e sem vizinhos */ - const double taxa_basal, - /** Taxa de morte dos indivíduos */ - const double taxa_morte, - /** Inclinação da curva de denso-depedência para natalidade */ - const double incl_b, - /** Inclinação da curva de denso-depedência para mortalidade */ - const double incl_d, - /** Número de pixels do lado da paisagem */ - const int numb_cells, - /** Tamanho do lada de cada pixel da imagem usada para a classificação da paisagem */ - const double cell_size, - /** Forma da paisagem (0 = circular, 1 = quadrada)*/ - const int land_shape, - /** Tipo de densidade ("g" = GLOBAL, "l" = LOCAL) */ - const int density_type, - /** Constante que indica quantas vezes a mortalidade basal na matriz é maior que no habitat */ - const double death_mat, - /** Posição dos indivíduos no início da simulação (0 = origem; 1 = aleatória com distribuição uniforme na paisagem; 2 = aleatória com distribuição normal na paisagem)*/ - const int inipos, - /** Condição de contorno (0 = absortiva, 1 = periódica, 2 = reflexiva)*/ - const int bound_condition, - /** Vetor de cobertura de habitat na paisagem */ - int scape[] - ); //construtor - - /** Atualiza as */ - void update();//atualizador - /** Seleciona o individuo com menor tempo até o proximo evento para realizar uma ação */ - int sorteia_individuo(); - /** Após a seleção do indivíduo que realizará a ação, sorteia uma das três ações possíveis de acordo com suas respectivas taxas e retorna a ação sorteada para a paisagem - * (0 = morte, 1 = nascimento, 2 = movimentação) */ - int sorteia_acao(const int lower){return this->popIndividuos[lower]->sorteia_acao();} - /** Realiza a ação sorteado do indivíduo selecionado. Além disso informa se o indivíduo sai da paisagem por emigração */ - bool realiza_acao(int acao, int lower);//vai pegar os tempos de cada individuo e informa qual foi o escolhido e manda ele fazer */ - /** Atualiza o contador de tempo, somando o tempo para o evento do indivíduo selecionado */ - void atualiza_tempo(const int lower){this->tempo_do_mundo = this->tempo_do_mundo + this->popIndividuos[lower]->get_tempo();} - /** Retorna o número total de indivíduos na paisagem */ + /** Radius of density dependant influence */ + const double raio, + /** Number of individuals at the start of the simulation */ + const int N, + /** Angle used for orientation when dispersing */ + const double angulo_visada, + /** The Lenght distance of a dispersal event (random walk) or maximum dispersal radius (Selection step) */ + const double passo, + /** The rate at which the individuals disperse */ + const double taxa_move, + /** The basal birth rate (The rate at which the individuals give birth on a habitat patch without neigbours) */ + const double taxa_basal, + /** The basal death rate (The rate at which the individuals die on a habitat patch without neigbours) */ + const double taxa_morte, + /** The slope of the birth density dependance function */ + const double incl_b, + /** The slope of the death density dependance functione */ + const double incl_d, + /** Number of pixels in a square landscape side */ + const int numb_cells, + /** Resolution Lenght of a square pixel side*/ + const double cell_size, + /** Shape of the landscape (0= circle, 1= square)*/ + const int land_shape, + /** Density type (0 = global, 1 = local/within a individual radius) */ + const int density_type, + /** Constant that indicates how many times higher the death rate should be on non-habitat pixels*/ + const double death_mat, + /** Constant that indicates how many times lower the movement rate should be on non-habitat pixels */ + const double move_mat, + /** The initial postion of individuals (0 = origin, 1 = random, 2 = normaly distributed with mean on origin)*/ + const int inipos, + /** The boundary condition type affects how individuals interact with the edges of the landscape (0= absortive, 1= periodical (pacman),2= reflexive)*/ + const int bound_condition, + /** Vector containing the evironmental values of the landscape pixels (0= matrix, 1= habitat) */ + double scape[], + /** Vector containing the x coordinates initial individuals */ + double initialPosX[], + /** Vector containing the y coordinates initial individuals */ + double initialPosY[], + /** Vector containing the genotypical trait means of the initial individuals */ + double genotype_means[], + /** Vector containing the standard deviations of the initial individuals */ + double width_sds[], + /** Vector containing the numer of points initial individuals sample when selecting habitats */ + int points[], + /** Vector containing the numer of points initial individuals sample when selecting habitats */ + bool Null + ); + + /** Function that calls other functions of the individual/individuo class to update the vector of individuals of the landscape/paisagem object (atualiza_vizinhos, atualiza_habitat, atualiza_patch, update()) + Paran int acao - A value representing the selected action (0= death, 1= birth, 2= dispersal) + Paranint ind - Index of the Individual performing the action + */ + void update(int acao, int ind); + + /** Function that uses the get_tempo function of the individual/individuo class to draft times for each individual within the landscape and returns the individual with the least amount of time required for the execution of the next action */ + int sorteia_individuo(); + + /** Function that recieves the selected individual and calls a function member of the individual/individuao class to randomly select one of the three possible actions for that individual to perform (0= death, 1= birth, 2= dispersal) + Paran const int lower - The postion of the individual with the lowest drafted time */ + int sorteia_acao(const int lower){return this->popIndividuos[lower]->sorteia_acao();} + + /** Function that executes the selected action. It also returns a positive boolean value if the individual disperses out of the landscape boundary + Paran int acao - A value representing the selected action (0= death, 1= birth, 2= dispersal) + Paran int lower - The postion of the individual with the lowest drafted time */ + bool realiza_acao(int acao, int lower); + + /** Function that updates the world clock, adding time required for the action execution by the selection individual + const int lower - The postion of the individual with the lowest drafted time */ + void atualiza_tempo(const int lower){this->tempo_do_mundo = this->tempo_do_mundo + this->popIndividuos[lower]->get_tempo();} + + /** Function that calls a function of the individual/individuo class to return the total number of individuals currently in the landscape */ const int conta_individuos() const{return popIndividuos.size();} - /** Retorna um vetor contendo todos os indivíduos na paisagem */ + + /** Function that returs an individual from within the landscape vector of individuals + Paran int i - The position(i vector) of the individual to be returned */ individuo* get_individuos(int i) const {return popIndividuos[i];} - /** Retorna o número de espécies existentes na paisagem \ref TBI */ + + /** Function that returns the curent number of species in the landscape + \ref TBI */ const int conta_especies() const; - /** Retorna o tamanho da paisagem (definido no construtor) */ + + /** Function that returns the lenght of square landscape side */ const double get_tamanho() const {return this->tamanho;} - /** Calcula a distância entre dois indivíduos */ - double calcDist(const individuo* a1, const individuo* a2) const; - /** Calcula a densidade de vizinhos de um indivíduo */ - double calcDensity(const individuo* ind1) const; - - /** Retorna false se o indivíduo estava no ambiente no passo de tempo 0, e true se ele nasceu durante a simulação. - * Usado para pintar os indivíduos nascidos de um cor diferente dos individuos originais */ + + /** Function that computes and returns the distance between a pair of individuals + Paran const individuo* a1 - First individual + Paran const individuo* a2 - Second individual */ + double calcDist(const individuo* a1, const individuo* a2) const; + + /** Function that computes and returns the density of individuals within a radius distance area from a focal individual + Paran const individuo* a1 - Focal individual */ + double calcDensity(const individuo* ind1) const; + + /** Function that returs a negative boolean value if the inserted individual was present at the start of the simulation and a positive value is the individuals was born during the simulation run time (used for colouring original individuals diferently) + Paran individuo * const ind - Individual to be tested */ const bool nascido(individuo * const ind) const {return ind->get_id() > this->N;} - /** Função recursiva utilizada para encontrar os freagmentos da paisagem */ + + /** Recursive function used to find and identify the witch pixels belong to each of the landscape patches + Paran int x - The initial value of the x coordinate + Paran int y - The initial value of the y coordinate + Paran int current_label - The initial path identification number */ void find_patches(int x, int y, int current_label); - /** Retorna o número de fragmentos da paisagem */ + + /** Function that returns the number of fragments on the landscape */ int get_numb_patches(){return numb_patches;} - /** Retorna a área de um dado fragmento da paisagem */ + + /** Function that computes and returns the area of a fragment in the landscape + Paran int i - The patch number */ double get_patch_area(int i) const {return this->patch_area[i];} - + + /** Function that decides upon, and calls other functions to perform, the desired method of dispersion (Random walk or habitat selection). + In the habitat selection case this function also samples x Points within the individuals "passo" radius and the indivuduals relative fitness on that location. + const int lower - The index of the individual with the lowest drafted time */ + bool walk(int lower); + }; #endif // PAISAGEM_H +