diff --git a/Analytical-pckg.R b/Analytical-pckg.R new file mode 100644 index 0000000..cb06597 --- /dev/null +++ b/Analytical-pckg.R @@ -0,0 +1,373 @@ +multiRun=function(PATH,nrep=10,b0,d0,m0,inc.b,inc.d,step,radius,dens.t,config,N0,lands=land,tm=100) +{ + oldpath <- getwd() + dir.create(PATH, recursive=TRUE) + setwd(PATH) + for(i in 1:nrep) + { + TWoLife(taxa.basal=b0, + taxa.morte=d0, + move=m0, + incl.birth=inc.b, + incl.death=inc.d, + passo=step, + raio=radius, + density.type=dens.t, + ini.config=config, + ##################### + N=N0, + AngVis=360, + death.mat=1, + landscape=lands, + tempo=tm, + out.code=i) + } + #Nt.reps=Nt.data(wdir=PATH,intervals=1,tmax=tm) + #plot.Nt(data=Nt.reps,growth=b0-d0,sum.incl=inc.b+inc.d) + #return(Nt.reps) + setwd(oldpath) +} + +# land <- Landscape(cover=1,type="b",cell.size=100) +# Dcoef=10000*0.8/4 +# sim1=multiRun(PATH="~/Desktop/Comb-0023", +# nrep=20, +# b0=0.4, +# d0=0.1, +# m0=0.8, +# inc.b=100450, +# inc.d=0, +# step=100, +# radius=1460, +# dens.t=1, +# config=0, +# N0=50, +# lands=land, +# tm=2000 +# ) + + +######################### +# Ancillary functions 1 +######################### +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) +} + +#################### +# Get N(t) data +#################### +Nt.data=function(wdir,intervals=1,tmax){ + oldpath <- getwd() + setwd(wdir) + files=dir() + Nt.arr=array(0,c(tmax+1,2,length(files))) + for(i in 1:length(files)) + { + t1=read.table(files[i],sep=" ") + Nt=as.data.frame(table(t1[,1])) + if(is.na(t1[dim(t1)[1],2])) # condicao que indica se houve extincao. + { + Nt[dim(Nt)[1],2]=0 # Corrije o data.frame Nt para que N = 0 quando t (coluna 1) = t de extincao. Isso se deve ao fato + # do table levar em conta o tempo de extincao como um fator representado uma vez na planilha da dados. + } + Nt[,1]=as.double(levels(Nt[,1])) + if(max(Nt[,1])<=tmax){miss=data.frame(Var1=seq(0,tmax,by=intervals))} else{miss=data.frame(Var1=seq(0,floor(max(Nt[,1])),by=intervals))} + Nt=merge(miss,Nt,all.x=T,all.y=T) + + index=sort(which(is.na(Nt[,2])),decreasing=T) + if(length(index>0)){ + if(index[1]==dim(Nt)[1]) + { + Nt[index[1],2]=0 + for(j in index[-1]) + { + Nt[j,2]=Nt[j+1,2] + } + } else{ + for(j in index) + { + Nt[j,2]=Nt[j+1,2] + } + } + } + Nt.arr[,1,i]=0:tmax + Nt.arr[,2,i]=Nt[1:(tmax+1),2] + } +return(Nt.arr) +setwd(oldpath) +} +# Nt.test=Nt.data("~/Desktop/Comb-0004",tmax=50) + +#################### +# Plot N(t) data +#################### +plot.Nt=function(dataSim=sim1,growth,sum.incl=0,land.area=10^8,name="N(t)",radius=0) +{ + plot(dataSim[,1,1],dataSim[,2,1],type="n",xlab="t",ylab="N (population size)",ylim=c(0,max(dataSim[,2,])+10), + cex=1.5, main = name) + for(i in 1:dim(dataSim)[3]) + { + lines(dataSim[,1,i],dataSim[,2,i],col="gray50") + } + if(sum.incl==0) + { + curve(dataSim[1,2,1]*exp(growth*x),add=T,col="red",lwd=4) + legend(dim(dataSim)[1]/4,max(dataSim[,2,]),legend=c("Observed","Predicted","Replicates"),lty=1,lwd=2,col=c(1,2,"gray50"),bty="n") + lines(0:(dim(dataSim)[1]-1),apply(dataSim[,2,which(dataSim[dim(dataSim)[1],2,]!=0)],1,mean),col=1,lwd=4) + } else + { + K=growth*land.area/sum.incl + lines(0:(dim(dataSim)[1]-1),apply(dataSim[,2,which(dataSim[dim(dataSim)[1],2,]!=0)],1,mean),col=1,lwd=4) + curve(K/(1+((K/dataSim[1,2,1])-1)*exp(-growth*x)),add=T,col="red",lwd=4) + legend(dim(dataSim)[1]/4,max(dataSim[,2,])/2,legend=c("Observed","Predicted","Replicates"),lty=1,lwd=2,col=c(1,2,"gray50"),bty="n") + } + +} +# plot.Nt(Nt.test,growth=0,0,10^8) +# plot.Nt(growth=0.3,sum.incl=100450) +# sim1=Nt.data(wdir="~/Desktop/Comb-0022",tmax=2000) + +######################### +# Ancillary functions 2 +######################### +displacement=function(dado) +{ + disp=sqrt(dado[,3]^2+dado[,4]^2) + return(disp) +} +dist.front=function(data,ind=10) +{ + vec2=sort(data,decreasing=T) + if(length(data)>=ind){front=vec2[ind]} else {front=vec2[length(data)]} + return(front) +} + +######################### +# Get velocity(t) data +######################### +velocity=function(work.dir="~/Desktop/Comb-0003",tmax=50,nrep=20,intervals=1) +{ + oldpath <- getwd() + setwd(work.dir) + files=dir() + vels.arr=matrix(0,tmax,(nrep+1)) + vels.arr[,1]=1:tmax + for(i in 1:nrep) + { + t2=read.table(files[i],sep=" ") + if(is.na(t2[dim(t2)[1],2])) # condicao que indica se houve extincao. + { + t2=t2[-(dim(t2)[1]),] # Corrije o data.frame Nt para que N = 0 quando t (coluna 1) = t de extincao. Isso se deve ao fato + # do table levar em conta o tempo de extincao como um fator representado uma vez na planilha da dados. + } + vec1=displacement(t2) + new.t2=data.frame(t2[,1],vec1) + front=aggregate(new.t2[,2],by = list(new.t2[,1]),FUN=dist.front,ind=10) + if(dim(front)[1]1){lines(object[,1],apply(object[,-1],1,mean),lwd=4)} # mean replicates (observed) + if(b0==0 && d0==0){curve(2*sqrt(Dcoef/x),col=2,add=T,lwd=3)} else {abline(h=2*sqrt(growth*Dcoef),col=2,lwd=3)} # expected +} + +# teste=velocity("~/Desktop/Comb-0014",nrep=100,tmax=50) +# Dcoef=10000*0.8/4 +# plot.vel(teste,Dcoef=Dcoef,nrep=20,growth=0.03) + +#################################### +# Plot XY density distribution data +#################################### +plot.XYdistrib=function(FILE="output-00001.txt",Dcoef,n.tests=5,tmax=50) +{ + times=sort(sample(1:tmax,n.tests,replace=F)) + t3=read.table(FILE,sep=" ") + + par(mfrow=c(n.tests,2)) + for(i in times) + { + #x + hist(t3[which(t3[,1]==i),3],freq=F,breaks=13) + curve(exp(-(x^2)/(4*Dcoef*i))/sqrt(4*pi*Dcoef*i),col=2,add=T) + #y + hist(t3[which(t3[,1]==i),4],freq=F,breaks=13) + curve(exp(-(x^2)/(4*Dcoef*i))/sqrt(4*pi*Dcoef*i),col=2,add=T) + } +} + +############################ +# Plot sd(x) & sd(y) +############################ +plot.sdXY=function(FILE="output-00001.txt",Dcoef,tmax=50,name="Standard Deviation") +{ + ## Teste desvio padrao + t3=read.table(FILE,sep=" ") + #par(mfrow=c(2,1)) + sd.x=aggregate(t3[,3],list(t3[,1]),sd) + sd.y=aggregate(t3[,4],list(t3[,1]),sd) + # sd x + plot(sd.x[,1],sd.x[,2],xlim=c(0,tmax),type="n",xlab="Time",ylab= "Std. Deviation of X positions",main=name) + lines(sd.x[,1],sd.x[,2],lwd=3) + curve(sqrt(2*Dcoef*x),col=2,add=T,lwd=3) + #sd y + plot(sd.y[,1],sd.y[,2],xlim=c(0,tmax),type="n",xlab="Time",ylab= "Std. Deviation of Y positions",main=name) + lines(sd.y[,1],sd.y[,2],lwd=3) + curve(sqrt(2*Dcoef*x),col=2,add=T,lwd=3) + #par(mfrow=c(1,1)) +} + +# Dcoef=10000*0.8/4 +# plot.XYdistrib(Dcoef=Dcoef) + +################ +# Plot XY data +################ + +# t3=read.table("output-00001.txt",sep=" ") +# for(i in 1:20){ +# plot(t3[which(t3[,1]==i),3],t3[which(t3[,1]==i),4],pch=16)} +# +# head(t3) +# t3 +# 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[,],pop[,], pch=4, col=2) +# } + + +########################### +# Critical Patch Size Test +########################### + +#params=c(0.007,0.004,0.8,0,0,500,0,0,0,20,0,1000) +multi.run.cP=function(PATH="~/Desktop/TWoTests/Critical_Patch/Comb-01",nrep=20,b0=0.007,d0=0.004,m0=0.8,inc.b=0,inc.d=0, + step=500,radius=0,dens.t=0,config=0,N0=20,tm=1000,lands.range=c(-5000,5000,1000)) +{ + oldpath <- getwd() + dir.create(PATH) + setwd(PATH) + Dcoef=m0*(step^2)/4 + A.crit=2*(pi^2)*Dcoef/(b0-d0) + land.sides=ceiling(sqrt(A.crit))+seq(lands.range[1],lands.range[2],by=lands.range[3]) + + for(j in 1:length(land.sides)) + { + lands=Landscape(cover=1,type="b",cell.size=land.sides[j]/10,numb.cells = 10) + dir.create(paste(j,land.sides[j],sep="_")) + setwd(paste(PATH,land.sides[j],sep="/")) + + for(i in 1:nrep) + { + TWoLife(taxa.basal=b0, + taxa.morte=d0, + move=m0, + incl.birth=inc.b, + incl.death=inc.d, + passo=step, + raio=radius, + density.type=dens.t, + ini.config=config, + ##################### + N=N0, + AngVis=360, + death.mat=1, + landscape=lands, + tempo=tm, + out.code=i) + } + setwd(PATH) + } + return(land.sides) + setwd(oldpath) +} + +##multi.run.cP(nrep=50) +##multi.run.cP(PATH="TWoTests/Critical_Patch/Comb-03",lands.range=c(2000,17000,1000),nrep=50) + +data.Acrit=function(PATH="~/Desktop/TWoTests/Critical_Patch/Comb-01",nreps=50,tmaxi=1000) +{ + directories=paste(PATH,dir(PATH),sep="/") + Nt.end=matrix(0,nreps,length(directories)) + for(i in 1:length(directories)) + { + teste=Nt.data(wdir=directories[i],tmax=tmaxi) + Nt.end[,i]=teste[tmaxi+1,2,] + } + return(list(Nt.end,teste)) +} +#testes=data.Acrit("~/Desktop/TWoTests/Critical_Patch/Comb-02") + +plot.Acrit=function(object,move=0.8,step=500,growth=0.003,N0=20,range=c(-17000,17000,1000),tmaxi=1000,Ncrit=0) +{ + Dcoef=move*(step^2)/4 + A.crit=2*(pi^2)*Dcoef/growth + sides=ceiling(sqrt(A.crit)+seq(range[1],range[2],by=range[3])) + areas=sides^2 + + par(mfrow=c(1,2)) + plot(rep(areas,each=dim(object)[1]),as.numeric(object),pch=1,xlab="Patch Size",ylab="Population size (Nt =1000)", + main=paste("Reaction-diffusion: Exponential Growth \n growth = ",growth,", D = ", Dcoef, ", N0 = ", N0,sep=""), + cex.main=0.8) + lines(areas,apply(Nt.end,2,mean),col=2) + abline(v=areas[ceiling(length(areas)/2)],h=N0*exp(growth*tmaxi),lty=c(1,3),col=c(4,1)) + + + plot(rep(areas,each=dim(object)[1]),as.numeric(object),type="n",ylim=c(0,1.1),xlab="Patch Size", + ylab="Extinction Probability for t = 1000", + main=paste("Reaction-diffusion: Exponential Growth \n growth = ",growth,", D = ", Dcoef, ", N0 = ", N0,sep=""), + cex.main=0.8) + for(i in Ncrit) + { + indexs=which(object<=i) + occurence=object + occurence[indexs]=0 + occurence[-indexs]=1 + points(areas,1-apply(occurence,2,sum)/dim(occurence)[1],col=i+1) + } + abline(v=areas[ceiling(length(areas)/2)],h=1,lty=c(1,3),col=c(2,1)) + legend(x =areas[ceiling(length(areas)/2)+2] ,y =0.6 ,"Critical Patch Size",bty="n",lty=3,col=1,cex=0.8) +} + +#plot.Acrit(teste[[1]]) + +################################# +###### Fragmented Landscape Tests +################################# diff --git a/TWoLife.R b/TWoLife.R index b0feef5..03b6d08 100644 --- a/TWoLife.R +++ b/TWoLife.R @@ -73,55 +73,80 @@ TWoLife <- function ( density.type=0, death.mat=7, landscape, - tempo=20) + 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") } - saida.C <- .C("TWoLife", as.double(raio), as.integer(N),as.double(AngVis), as.double(passo), - as.double(move), as.double(taxa.basal), as.double(taxa.morte), - as.double(incl.birth), as.double(incl.death), as.integer(landscape$numb.cells), - as.double(landscape$cell.size), as.integer(landscape$land.shape), as.integer(density.type), - as.double(death.mat), as.integer(landscape$bound.condition), as.integer(landscape$scape), - as.double(tempo), as.integer(0), - as.double(rep(0, 5000)), as.double(rep(0,5000)) ## verificar se precisa definir o tamanho e se isto nao dará problemas + 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[[18]] - x <- saida.C[[19]] - y <- saida.C[[20]] + 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)) } -## Um teste rapido -land <- Landscape(cover=0.95, type="b") -## Uma rodada: coordenadas dos sobreviventes apos t=20 -teste <- TWoLife(raio=0.1, - N=80, - AngVis=360, - passo=5, - move=3, - taxa.basal=2, - taxa.morte=0.1, - incl.birth=0.5/0.01, - incl.death=0, - density.type=0, - death.mat=700, - landscape=land, - tempo=30) +# ## 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) -TWoPlot <- function(pop, land) { - 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 = "gray20" - } else { - color = c("gray70", "gray20") - } - image(s, s, matrix(land$scape,ncol=n), col=color) - points(pop, pch=4, col=2) -} -TWoPlot(teste, 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 @@ -137,11 +162,11 @@ TWoPlot(teste, land) # } ## 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) -} +# 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.cpp b/TWoLife.cpp index ffcf8e2..1064099 100644 --- a/TWoLife.cpp +++ b/TWoLife.cpp @@ -1,10 +1,13 @@ #include +#include #include #include #include #include #include "paisagem.cpp" #include "individuo.cpp" +#include +#include using namespace std; @@ -14,28 +17,61 @@ using namespace std; /** Funcao que recebe os parâmetros da simulação e a executa, -// baseada no texto em 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! -// -// Veja a descrição da classe \ref paisagem para o significado dos parâmetros -// */ + // baseada no texto em 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! + // + // Veja a descrição da classe \ref paisagem para o significado dos parâmetros + // */ 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 * bound_condition, int * scape, double * tempo, int * nPop, double * x, double * y) + 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) { - GetRNGstate(); /* (mudar para doxygen): este comando chama o engine de numeros aleatorios do R - Por causa dela nossa biblioteca nao eh standalone */ + // This sequence creates an attribute containing the output file name. The template is output-00000.txt. + string fileNAME = "output-00000.txt"; + stringstream tmps; + tmps<tempo_do_mundo < tempo[0] && floresta->conta_individuos() > 0) { - floresta->update(); + move[0], taxa_basal[0], taxa_morte[0], incl_b[0], + incl_d[0], numb_cells[0], cell_size[0], land_shape[0], + density_type[0], death_mat[0], inipos[0], bound_condition[0], + scape); + + ofstream outputSIM; // ofstream for the output file + outputSIM.open(fileNAME.c_str()); + for(unsigned int i=0; iconta_individuos();i++) + { + outputSIM << floresta->tempo_do_mundo << " " << floresta->get_individuos(i)->get_id() << " " << floresta->get_individuos(i)->get_x() << " " << floresta->get_individuos(i)->get_y() << endl; } + + while (floresta->tempo_do_mundo < tempo[0] && floresta->conta_individuos() > 0) + { + double t_ant = floresta->tempo_do_mundo; + int lowerInd=floresta->update(); + if(t_ant < (int)floresta->tempo_do_mundo) + { + for(unsigned int i=0; iconta_individuos();i++) + { + outputSIM << (int)floresta->tempo_do_mundo << " " << floresta->get_individuos(i)->get_id() << " " << floresta->get_individuos(i)->get_x() << " " << floresta->get_individuos(i)->get_y() << endl; + } + } + floresta->realiza_acao(lowerInd); + } + if(floresta->conta_individuos()==0){outputSIM << floresta->tempo_do_mundo << " " << "NA" << " " << "NA" << " " << "NA" << endl;} + 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(); @@ -44,3 +80,5 @@ extern "C" void TWoLife (double * raio, int * N, double * angulo_visada, double delete floresta; PutRNGstate(); } + + diff --git a/Tests.R b/Tests.R new file mode 100644 index 0000000..218bd74 --- /dev/null +++ b/Tests.R @@ -0,0 +1,325 @@ +## source("TWoLife.R") ## calls the C++ sources with parallel. +source("Analytical-pckg.R") +########################################## +########## Set of Parameteres Combinations +########################################## +b0s=c(rep(0.04,6),# only growth tests: exponential + rep(rep(c(0.4,0.7,0.4),each=2),2),#only growth: logistic - global density + rep(0.4,6),rep(0.7,2),#only growth: logistic - local density + rep(0,3),# only diffusion + rep(c(0.04,0.07,0.04),2), # reaction-diffusion: Skellam model (exponential growth) + rep(c(rep(0.4,3),0.7),2), # reaction-diffusion: FIsher-Kolmogorov (logistic growth, global density) + rep(c(rep(0.4,3),0.7),2)) # reaction-diffusion: FIsher-Kolmogorov (logistic growth, local density) +d0s=c(rep(0.01,6),# only growth tests: exponential + rep(rep(c(0.1,0.4,0.1),each=2),2),#only growth: logistic - global density + rep(0.1,6),rep(0.4,2),#only growth: logistic - local density + rep(0,3),# only diffusion + rep(c(0.01,0.04,0.01),2), # reaction-diffusion: Skellam model (exponential growth) + rep(c(rep(0.1,3),0.4),2), # reaction-diffusion: FIsher-Kolmogorov (logistic growth, global density) + rep(c(rep(0.1,3),0.4),2)) # reaction-diffusion: FIsher-Kolmogorov (logistic growth, local density) +m0s=c(rep(0,26), # only growth + 0.8,0.8,0.4, # only diffusion + rep(c(rep(0.8,2),1.6),2), # reaction-diffusion: Skellam model (exponential growth) + rep(0.8,16)) +incl.bs=c(rep(0,6), # only growth tests: exponential + rep(c(rep(10^5,4),rep(0.5*10^5,2)),2), # only growth: logistic - global density + rep(100500,2),rep(1750,2),rep(100500,4), # only growth: logistic - local density + 0,0,0, # only diffusion + rep(0,6), # reaction-diffusion: Skellam model (exponential growth) + rep(c(rep(100500,2),50250,100500),2), # reaction-diffusion: FIsher-Kolmogorov (logistic growth, global density) + rep(c(rep(100500,2),50250,100500),2)) # reaction-diffusion: FIsher-Kolmogorov (logistic growth, local density) +incl.ds=c(rep(0,6), # only growth tests: exponential + rep(c(rep(0,4),rep(0.5*10^5,2)),2), # only growth: logistic - global density + rep(0,8), # only growth: logistic - local density + 0,0,0, # only diffusion + rep(0,6), # reaction-diffusion: Skellam model (exponential growth) + rep(c(rep(0,2),50250,0),2), # reaction-diffusion: FIsher-Kolmogorov (logistic growth, global density) + rep(c(rep(0,2),50250,0),2)) # reaction-diffusion: FIsher-Kolmogorov (logistic growth, local density) +steps=c(rep(0,26), # only growth + rep(100,3), # only diffusion + rep(100,6), # reaction-diffusion: Skellam model (exponential growth) + rep(100,16)) # reaction-diffusion: Fisher-Kolmogorov model (logistic growth; global and local density) +radii=c(rep(0,18), # only growth: exponential + logistic global density + 750,1500,rep(750,6), # only growth: logistic local density + rep(0,3), # only diffusion + rep(0,6), # reaction-diffusion: Skellam model (exponential growth) + rep(0,8), # reaction-diffusion: Fisher-Kolmogorov model (logistic growth; global density) + rep(c(750,1500,750,750),2)) # reaction-diffusion: Fisher-Kolmogorov model (logistic growth; local density) +d.types=c(rep(0,18), # only growth: exponential + logistic global density + rep(1,8), # only growth: logistic local density + rep(0,3), # only diffusion + rep(0,6), # reaction-diffusion: Skellam model (exponential growth) + rep(0,8), # reaction-diffusion: Fisher-Kolmogorov model (logistic growth; global density) + rep(1,8)) +i.configs=c(rep(0:2,2), + rep(c(0,1),each=6), + rep(0,4),rep(1,2),0,0, + 0,2,0, + rep(c(0,1),each=3), + rep(c(0,1),each=4), + rep(c(0,1),each=4)) +N0s=c(rep(c(50,10),each=3), + rep(c(50,500),6), + rep(50,3),500,rep(50,4), + rep(500,3), + rep(50,6), + rep(50,16)) +d.matrix=rep(0,51) +times=c(rep(c(50,100),each=3), + rep(2000,20), + rep(100,3), + rep(50,6), + rep(2000,16)) +#n.replicates=c() + +cb1=paste("0",1:9,sep="") +cb2=paste(10:51) +CODE=c(cb1,cb2) + +# Parameter combinations table for testing TwoLife +params=data.frame(b0=b0s,d0=d0s,movem.rate=m0s,incl.birth=incl.bs,incl.death=incl.ds,step.length=steps,radius=radii, + dens.type=d.types,initial.config=i.configs,N0=N0s,death.matrix=d.matrix,tmax=times) +# params$b0[51] +# dim(params) +## Create a directory for the output files, here named TWotests +land <- Landscape(cover=1,type="b",cell.size=100) +combs=1:51 +for (i in combs) +{ + multiRun(PATH=paste("./TWoTests/Comb-",CODE[i],sep=""), + nrep=20, + b0=params$b0[i], + d0=params$d0[i], + m0=params$movem.rate[i], + inc.b=params$incl.birth[i], + inc.d=params$incl.death[i], + step=params$step.length[i], + radius=params$radius[i], + dens.t=params$dens.type[i], + config=params$initial.config[i], + N0=params$N0[i], + lands=land, + tm=params$tmax[i]) +} + +################################### +################################### +################################### + + + +###################### +##### Plot main titles +###################### +titles.exp=paste("Simple Birth-Death (COMB ", CODE[1:6],") \n", expression(r), + paste(" = ", params[1:6,1]-params[1:6,2],", ini.config = ",params[1:6,9],sep=""), sep="") +titles.log1=paste("General Birth-Death - Logistic Growth (COMB ", CODE[7:18],") \n", expression(b0), + paste(" =", params[7:18,1],", "), expression(slope(b)),paste(" =", params[7:18,4],", "), + expression(d0), paste(" =", params[7:18,2],", "),expression(slope(d)),paste(" =", params[7:18,5],"\n"), + paste("ini.config = ",params[7:18,9],", "),expression(N0),paste(" =", params[7:18,10]),sep="") +titles.log2= paste("General Birth-Death - Logistic Growth (COMB ", CODE[19:26],") \n", + expression(b0), paste(" =", params[19:26,1],", "), expression(slope(b)),paste(" =", params[19:26,4],", "), + expression(d0), paste(" =", params[19:26,2],", "), expression(slope(d)),paste(" =", params[19:26,5],"\n"), + "movement rate = ", params[19:26,3], ", step length = ", params[19:26,6], ", R = ", params[19:26,7], + ", ini.config = ",params[19:26,9],sep="") +titles.mov=paste("Simple Random Walk (COMB ", CODE[27:29],") \n", "movement rate = ", params[27:29,3], + ", step length = ", params[27:29,6], ", ini.config = ",params[27:29,9], sep="") +titles.skellam= paste("Reaction-Diffusion with Exponential Growth (COMB ", CODE[30:35],") \n", + expression(b0), paste(" =", params[30:35,1],", "), expression(d0), paste(" =", params[30:35,2],"\n"), + "movement rate = ", params[30:35,3], ", step length = ", params[30:35,6], + ", ini.config = ",params[30:35,9],sep="") +titles.fisher1= paste("Reaction-Diffusion with Logistic Growth (COMB ", CODE[36:43],") \n", + expression(b0), paste(" =", params[36:43,1],", "), expression(slope(b)),paste(" =", params[36:43,4],", "), + expression(d0), paste(" =", params[36:43,2],", "), expression(slope(d)),paste(" =", params[36:43,5],"\n"), + "movement rate = ", params[36:43,3], ", step length = ", params[36:43,6], + ", ini.config = ",params[36:43,9],sep="") +titles.fisher2= paste("Reaction-Diffusion with Logistic Growth (COMB ", CODE[44:51],") \n", + expression(b0), paste(" =", params[44:51,1],", "), expression(slope(b)),paste(" =", params[44:51,4],", "), + expression(d0), paste(" =", params[44:51,2],", "), expression(slope(d)),paste(" =", params[44:51,5],"\n"), + "movement rate = ", params[44:51,3], ", step length = ", params[44:51,6], ", R = ", params[44:51,7], + ", ini.config = ",params[44:51,9],sep="") +titles=c(titles.exp,titles.log1,titles.log2,titles.mov,titles.skellam,titles.fisher1,titles.fisher2) +titles +########### + + +############### +####### Plots +## First create a directory for the plots (./TWoResults here) +############### +paths=paste("./TWoTests/Comb-",CODE,sep="") +paths +ind1=c(1,4,2,5,3,6) + +pdf("./TWoResults/Nt/Nt_Simple-Birth-Death.pdf",width = 8,height = 12) +par(mfrow=c(3,2)) +for(i in ind1) +{ + teste=Nt.data(wdir = paths[i],tmax=params$tmax[i]) + setwd("./TWoResults/Nt") + plot.Nt(teste,growth=params$b0[i]-params$d0[i],sum.incl = params$incl.birth[i]+params$incl.death[i],name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() + +ind=c(19:23,25) + +pdf("./TWoResults/Nt/Nt_General-Birth-Death-gd2.pdf",width = 8,height = 12) +par(mfrow=c(3,2)) +for(i in ind) +{ + teste=Nt.data(wdir = paths[i],tmax=params$tmax[i]) + setwd("./TWoResults/Nt") + plot.Nt(teste,growth=params$b0[i]-params$d0[i],sum.incl = params$incl.birth[i]+params$incl.death[i],name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() + +ind=c(19:23,25) + +pdf("./TWoResults/Nt/Nt_General-Birth-Death-ld.pdf",width = 8,height = 12) +par(mfrow=c(3,2)) +for(i in ind) +{ + teste=Nt.data(wdir = paths[i],tmax=params$tmax[i]) + setwd("./TWoResults/Nt") + plot.Nt(teste,growth=params$b0[i]-params$d0[i],sum.incl = params$incl.birth[i]+params$incl.death[i],name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() + + + +################# +################# +ind=27:29 +pdf("./TWoResults/Velocity/Vel_Simple-Random-Walk.pdf",width = 8,height = 8) +par(mfrow=c(2,2)) +for(i in ind) +{ + teste=velocity(work.dir = paths[i],tmax = params$tmax[i],nrep = 50) + setwd("./TWoResults/Velocity") + plot.vel(teste,Dcoef = (params$step.length[i]^2)*params$movem.rate[i]/4, growth=params$b0[i]-params$d0[i], + name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() +####################### +pdf("./TWoResults/Distribution/Dist_Simple-Random-Walk.pdf",width = 8,height = 12) +par(mfrow=c(length(ind),2)) +for(i in ind) +{ + Dcoefs=(params$step.length[i]^2)*params$movem.rate[i]/4 + plot.sdXY(FILE = paste(paths[i],"/output-00001.txt",sep=""),Dcoef = Dcoefs,tmax = params$tmax[i],name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() + +############### +############### +ind=c(30,33,31,34,32,35) +pdf("./TWoResults/Velocity/Vel_Skellam-b.pdf",width = 8,height = 12) +par(mfrow=c(3,2)) +for(i in ind) +{ + Dcoefs=(params$step.length[i]^2)*params$movem.rate[i]/4 + growths=params$b0[i]-params$d0[i] + teste=velocity(work.dir = paths[i],tmax = params$tmax[i],nrep = 50) + setwd("./TWoResults/Velocity") + plot.vel(teste,Dcoef = Dcoefs, b0=params$b0[i],d0=params$d0[i], yrange = c(-1,30),name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() +### +pdf("./TWoResults/Nt/Nt_Skellam.pdf",width = 8,height = 12) +par(mfrow=c(3,2)) +for(i in ind) +{ + teste=Nt.data(wdir = paths[i],tmax=params$tmax[i]) + setwd("./TWoResults/Nt") + plot.Nt(teste,growth=params$b0[i]-params$d0[i],sum.incl = params$incl.birth[i]+params$incl.death[i],name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() + + +ind=c(36:43) +pdf("./TWoResults/Velocity/Vel_Fisher-gd-b.pdf",width = 8,height = 16) +par(mfrow=c(4,2)) +for(i in ind) +{ + Dcoefs=(params$step.length[i]^2)*params$movem.rate[i]/4 + growths=params$b0[i]-params$d0[i] + teste=velocity(work.dir = paths[i],tmax = params$tmax[i],nrep = 20) + setwd("./TWoResults/Velocity") + plot.vel(teste,Dcoef = Dcoefs, b0=params$b0[i],d0=params$d0[i], yrange = c(-100,100),name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() +pdf("./TWoResults/Nt/Nt_Fisher-gd.pdf",width = 8,height = 12) +par(mfrow=c(4,2)) +for(i in ind) +{ + teste=Nt.data(wdir = paths[i],tmax=params$tmax[i]) + setwd("./TWoResults/Nt") + plot.Nt(teste,growth=params$b0[i]-params$d0[i],sum.incl = params$incl.birth[i]+params$incl.death[i],name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() + +ind=c(44:47) +pdf("./TWoResults/Velocity/Vel_Fisher-ld.pdf",width = 8,height = 8) +par(mfrow=c(2,2)) +for(i in ind) +{ + teste=velocity(work.dir = paths[i],tmax = params$tmax[i],nrep = 20) + setwd("./TWoResults/Velocity") + plot.vel(teste,Dcoef = (params$step.length[i]^2)*params$movem.rate[i]/4, growth=params$b0[i]-params$d0[i],name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() +pdf("./TWoResults/Nt/Nt_Fisher-ld.pdf",width = 8,height = 8) +par(mfrow=c(2,2)) +for(i in ind) +{ + teste=Nt.data(wdir = paths[i],tmax=params$tmax[i]) + setwd("./TWoResults/Nt") + plot.Nt(teste,growth=params$b0[i]-params$d0[i],sum.incl = params$incl.birth[i]+params$incl.death[i],name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() + +ind=c(48:51) +pdf("./TWoResults/Velocity/Vel_Fisher-ld-c2b.pdf",width = 8,height = 8) +par(mfrow=c(2,2)) +for(i in ind) +{ + Dcoefs=(params$step.length[i]^2)*params$movem.rate[i]/4 + growths=params$b0[i]-params$d0[i] + teste=velocity(work.dir = paths[i],tmax = params$tmax[i],nrep = 20) + setwd("./TWoResults/Velocity") + plot.vel(teste,Dcoef = Dcoefs, b0=params$b0[i],d0=params$d0[i], yrange = c(-150+2*sqrt(growths*Dcoefs),50+2*sqrt(growths*Dcoefs)), + name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() +pdf("./TWoResults/Nt/Nt_Fisher-ld-c2.pdf",width = 8,height = 8) +par(mfrow=c(2,2)) +for(i in ind) +{ + teste=Nt.data(wdir = paths[i],tmax=params$tmax[i]) + setwd("./TWoResults/Nt") + plot.Nt(teste,growth=params$b0[i]-params$d0[i],sum.incl = params$incl.birth[i]+params$incl.death[i], + name = titles[i]) +} +par(mfrow=c(1,1)) +dev.off() + +##################### +##################### +testes=data.Acrit("./TWoTests/Critical_Patch/Comb-02") + +pdf("./TWoResults/Acrit/Skellam_Comb-02.pdf",width = 8,height = 5) +plot.Acrit(teste[[1]]) +dev.off() + diff --git a/individuo.cpp b/individuo.cpp index a41b1bd..ea7dc49 100644 --- a/individuo.cpp +++ b/individuo.cpp @@ -101,13 +101,13 @@ individuo::individuo(const individuo& rhs) * - 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). */ -void individuo::update() +void individuo::update(double dens) { - double densi = (this->lisViz.size()+1)/(M_PI*(this->raio*this->raio)); // densidade inclui n de vizinhos + o individuo + double densi = dens; // densidade inclui n de vizinhos + o individuo if(this->tipo_habitat==0) { this->birth = 0; - // 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. + // 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 diff --git a/individuo.hpp b/individuo.hpp index 56d944b..a06ef47 100644 --- a/individuo.hpp +++ b/individuo.hpp @@ -100,8 +100,10 @@ class individuo inline const double get_y() const {return this->y;} /** Retorna o raio de percep√ßa√£o do indiv√≠duo */ const double get_raio() const {return this->raio;} - /** Retorna o tiop de densidadeque afeta o indiv√≠duo */ + /** Retorna o tipo de densidade que afeta o indiv√≠duo */ 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) */ + const int NBHood_size() const {return this->lisViz.size()+1;} // outros metodos publicos /** Retorna o tempo sorteado para o próximo evento acontecer com este indivíduo. @@ -115,7 +117,7 @@ class individuo /** Atualiza a taxa de nascimento e/ou de morte baseado na densidade de indiv√≠duos dentro do raio de percep√ß√£o e sorteia o tempo * do pr√≥ximo evento baseado nas taxas de nascimento, morte e movimenta√ß√£o * \sa \ref individuo::get_tempo */ - void update(); + void update(double dens); /** Faz com que o indiv√≠duo ande um passo, do 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. */ diff --git a/paisagem.cpp b/paisagem.cpp index e63e650..119163f 100644 --- a/paisagem.cpp +++ b/paisagem.cpp @@ -5,14 +5,15 @@ 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 bound_condition, int scape[]): + 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) + boundary_condition(bound_condition), + initialPos(inipos) { for(int i=0;inumb_cells;i++) { @@ -23,7 +24,7 @@ paisagem::paisagem(double raio, int N, double angulo_visada, double passo, doubl } } - // Calculo do raio dependendo do tipo de densidade + // 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); @@ -39,30 +40,77 @@ void paisagem::populating(double raio, int N, double angulo_visada, double passo { individuo::reset_id(); // reinicia o contador de id dos individuos // Considerar diferentes possibilidades de posições iniciais. TBI. - 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==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)); + } } + 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() +int paisagem::update() { if(this->popIndividuos.size()>0) { @@ -79,14 +127,52 @@ void paisagem::update() // aleatorias sao chamadas sempre na mesma ordem (garante reprodutibilidade) for(unsigned int i=0; ipopIndividuos.size(); i++) { - this->popIndividuos[i]->update(); //e atualiza o individuo i da populacao + double dsty=this->calcDensity(popIndividuos[i]); + this->popIndividuos[i]->update(dsty); //e atualiza o individuo i da populacao } + + // 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(); + } + } + this->tempo_do_mundo = this->tempo_do_mundo+menor_tempo; + return menor; + } +} + +void paisagem::realiza_acao(int lower) //TODO : criar matriz de distancias como atributo do mundo e atualiza-la apenas quanto ao individuos afetado nesta funcao) +{ + int acao = this->popIndividuos[lower]->sorteia_acao(); - this->realiza_acao();//escolhe tempo, indica e faz + switch(acao) //0 eh morte, 1 eh nascer, 2 eh andar + { + 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(); + this->apply_boundary(popIndividuos[lower]); + break; } } - //para quando tiver especies, a definir... //int paisagem::conta_especies() //{ @@ -109,43 +195,6 @@ void paisagem::update() // } //} -void paisagem::realiza_acao() //TODO : criar matriz de distancias como atributo do mundo e atualiza-la apenas quanto ao individuos afetado nesta funcao) -{ - 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(); - } - } - - this->tempo_do_mundo = this->tempo_do_mundo+menor_tempo; - int acao = this->popIndividuos[menor]->sorteia_acao(); - - switch(acao) //0 eh morte, 1 eh nascer, 2 eh andar - { - case 0: - delete this->popIndividuos[menor]; - this->popIndividuos.erase(this->popIndividuos.begin()+menor); - break; - - case 1: - individuo* chosen; - //Novo metodo para fazer copia do individuo: - chosen = new individuo(*this->popIndividuos[menor]); - this->popIndividuos.push_back(chosen); - break; - - case 2: - this->popIndividuos[menor]->anda(); - this->apply_boundary(popIndividuos[menor]); - break; - } -} // 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 @@ -182,7 +231,7 @@ void paisagem::apply_boundary(individuo * const ind) //const { for(unsigned int i=0; ipopIndividuos[i]->get_id()==(int)ind->get_id()) + if(this->popIndividuos[i]->get_id()==(int)ind->get_id()) //DUVIDA: porque tem int? { delete this->popIndividuos[i]; this->popIndividuos.erase(this->popIndividuos.begin()+i); @@ -232,6 +281,120 @@ double paisagem::calcDist(const individuo* a1, const individuo* a2) const //Viro } } +// 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; +} + + /* 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 @@ -270,9 +433,13 @@ void paisagem::atualiza_habitat(individuo * const ind) const // 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+this->numb_cells/2; + hy= ((double)ind->get_y()/this->cell_size)*(-1)+this->numb_cells/2; ind->set_habitat(this->landscape[hx][hy]); } diff --git a/paisagem.hpp b/paisagem.hpp index 0a1d196..73d2b34 100644 --- a/paisagem.hpp +++ b/paisagem.hpp @@ -2,6 +2,7 @@ #define PAISAGEM_H #include "individuo.hpp" #include +#include #include #include ///** Dimensão maxima do mapa, usada para rotinas de fragmentação TBI */ @@ -26,6 +27,7 @@ class paisagem const int landscape_shape; // paramâmetro do modelo relacionado à forma da paisagem (por enquanto, 1 = "quadrada" ou 0 = "circular") const int boundary_condition; // tipo de condição de contorno (0 = absortiva, 1 = periódica, 2 = reflexiva) int landscape[dim][dim];//[linha][coluna] temporariamente substituido or valor fixo + const int initialPos; //metodos privados void populating( @@ -48,17 +50,18 @@ class paisagem const double incl_d, /** Constante de acrescimo na taxa basal de morte na matriz em relação à essa taxa basal no habitat */ const double death_m, - /** Tipo de densidade ("g" = GLOBAL, "l" = LOCAL) */ + /** Tipo de densidade (0 = GLOBAL, 1 = LOCAL) */ const int dens_type ); void atualiza_vizinhos(individuo * const ind) const;//contabilizador de vizinhos void atualiza_habitat(individuo * const i) const;//vai informar o individuo em que tipo de habitat ele esta - void realiza_acao();//vai pegar os tempos de cada individuo e informa qual foi o escolhido e manda ele fazer + //int define_tempo(); void apply_boundary(individuo * const ind); //const; // metodo para aplicação da condicao de contorno public: + //vector popIndividuos; /** Contador de quanto tempo já transcorreu no mundo simulado */ double tempo_do_mundo; @@ -94,14 +97,17 @@ class paisagem const int density_type, /** Constante de acrescimo na taxa basal de morte na matriz em relação à essa taxa basal no habitat */ const double death_mat, + /** How individuals are initially set into the landscpae **/ + const int inipos, /** Condição de contorno */ const int bound_condition, /** Vetor de cobertura de habitat na paisagem */ - int scape[] + int scape[] ); //construtor /** Atualiza as variáveis de todos os indivíduos (ver individuo::set_vizinhos, individuo::set_habitat e individuo::update) e escolhe uma ação para ser executada. Executa a ação e atualiza o tempo do mundo de acordo \sa \ref Introdução */ - void update();//atualizador + int update();//atualizador + void realiza_acao(int lower);//vai pegar os tempos de cada individuo e informa qual foi o escolhido e manda ele fazer /** Retorna o número total de indivíduos na paisagem */ const int conta_individuos() const{return popIndividuos.size();} /** Retorna um vetor contendo todos os indivíduos na paisagem */ @@ -112,6 +118,8 @@ class paisagem const double get_tamanho() const {return this->tamanho;} double calcDist(const individuo* a1, const individuo* a2) const; + + double calcDensity(const individuo* ind1) const; /** Retorna false se o indivíduo estava no ambiente no passod 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 */