Nothing
R/demoniche_model.R
In demoniche: demoniche - spatial population dynamics
demoniche_model <-
function(modelname, Niche, Dispersal, repetitions, foldername)
{
BEMDEM <- get(modelname, envir = .GlobalEnv) #assign("BEMDEM", get(modelname))
#load(paste(modelname, ".rda", sep = ""))
# rm(Hmontana)
# require(sp)
require(popbio)
require(lattice)
#### Create vectors ############################################################
Projection <- array(0, # Projection is where all the results go
dim =
c(BEMDEM$no_yrs, length(BEMDEM$stages), nrow(BEMDEM$Niche_ID), length(BEMDEM$years_projections))
,
dimnames =
list(paste( "timesliceyear", 1:BEMDEM$no_yrs, sep="_"), c(paste(BEMDEM$stages)),
BEMDEM$Niche_ID[,"Niche_ID"], paste(BEMDEM$years_projections))
)
eigen_results <- vector(mode = "list" , length(BEMDEM$list_names_matrices))
names(eigen_results) <- unlist(BEMDEM$list_names_matrices)
yrs_total <- BEMDEM$no_yrs * length(BEMDEM$years_projections)
population_sizes <- array(NA, # per matrix!
dim = c(yrs_total, length(BEMDEM$list_names_matrices), repetitions),
dimnames = list(paste("year", 1: yrs_total, sep = ""), BEMDEM$list_names_matrices, paste("rep", 1:repetitions, sep = "_")))
population_results <- array(1:200, # for all matrices!
dim = c(yrs_total, 4, length(BEMDEM$list_names_matrices)),
dimnames = list( paste("year", 1: yrs_total, sep = ""),
c("Meanpop", "SD", "Max", "Min"),
paste(BEMDEM$list_names_matrices)))
metapop_results <- array(NA, # per matrix!
dim = c(yrs_total, length(BEMDEM$list_names_matrices), repetitions),
dimnames = list(paste("year", 1: yrs_total, sep = ""), BEMDEM$list_names_matrices, paste("rep", 1:repetitions, sep = "_")))
simulation_results <- array(NA,
dim = c(length(BEMDEM$list_names_matrices), 7+length(BEMDEM$years_projections)),
dimnames =
list(BEMDEM$list_names_matrices,
c("lambda", "stoch_lambda","mean_perc_ext_final", "initial_population_area", "initial_population", "mean_final_pop", "mean_no_patches_final",
paste("EMA", BEMDEM$years_projections))
) )
EMA <- array(0,
dim = c(repetitions, length(BEMDEM$list_names_matrices), length(BEMDEM$years_projections), 2), dimnames =
list(paste("rep", 1:repetitions, sep = "_"), BEMDEM$list_names_matrices, BEMDEM$years_projections, c("EMA", "No_populations")))
population_Niche <- rep(1, nrow(BEMDEM$Niche_ID))
simulation_results[,"initial_population_area"] <-
sum(BEMDEM$Orig_Populations[,"area_population"])
# original population size
simulation_results[,"initial_population"] <-
round(sum(colSums(BEMDEM$n0_all) * BEMDEM$sumweight), 0)
dir.create(paste(getwd(), "/", foldername, sep = ""), showWarnings = FALSE)
#### Start repetitions #########################################################
for (rx in 1:repetitions) # tx = 1 rx = 1
{
print(paste("Starting projections for repetition:", rx), quote = FALSE)
for(mx in 1:length(BEMDEM$list_names_matrices)) # selects two matrices for the simulations mx =1
{
print(paste("Projecting for scenario/matrix:", (BEMDEM$list_names_matrices)[mx]), quote = FALSE)
yx_tx <- 0 # restart counter
# this selects two matrices, the first basic matrix and the projection one mx =1
Matrix_projection <- cbind(BEMDEM$matrices[,1], (BEMDEM$matrices[, mx]))
# this selects two standard deviation matrices, the first basic matrix and the projection one (column mx)
if(BEMDEM$matrices_var[1] != FALSE) # If matrices are included or not.
{
if(ncol(BEMDEM$matrices_var) > 1){
Matrix_projection_var <- cbind(BEMDEM$matrices_var[,1], (BEMDEM$matrices_var[, mx]))
} else {
Matrix_projection_var <- cbind(BEMDEM$matrices_var[,1], (BEMDEM$matrices_var[, 1]))
}
} else
{
Matrix_projection_var <- FALSE
}
prev_mx <- rep(1, times = yrs_total + 1)
for(tx in 1:length(BEMDEM$years_projections)) # selects which time-slice of the simulation tx = 1
{
if (Niche == TRUE){ # Niche values
population_Niche <- BEMDEM$Niche_values[, tx]
}
for(yx in 1:BEMDEM$no_yrs)
{
yx_tx <- yx_tx + 1 # add one year to counter
####################################### PATCH PROJECTION when tx = 1 ########### tx=1 px=1 tx = 2
################################################################################ yx=2 yx = 1
if(tx == 1 && yx == 1)
{
n0s <- BEMDEM$n0_all[rowSums(BEMDEM$n0_all) > 0,] # take stage vector from intial population, where over zero
n0s_ID <- which(rowSums(BEMDEM$n0_all) > 0) # which populations/rows are above zero
} else { # Second year running the model
if(tx != 1 && yx == 1){ # when going to next time step, yx == 1
n0s <- t(Projection[BEMDEM$no_yrs,, colSums(Projection[BEMDEM$no_yrs,,,tx-1]) > 0, tx-1])
n0s_ID <- which(colSums(Projection[BEMDEM$no_yrs,,,tx-1]) > 0)
BEMDEM$Niche_ID[colSums(Projection[BEMDEM$no_yrs,,,tx-1]) > 0, 2]
} else {
# take population from previous year, same time period
n0s <- t(Projection[yx-1,, colSums(Projection[yx-1,,,tx]) > 0,tx])
n0s_ID <- which(colSums(Projection[yx-1,,,tx]) > 0)
}
} # end if tx == 1 && yx == 1
# n <- c(10,5,5,5,5,5)
population_Niche_short <- population_Niche[n0s_ID]
# run fcn for each population px=1
if (nrow(n0s) > 0) {
for(px in 1:nrow(n0s)) { # px = 1 tx =1
n <- as.vector(n0s[px,])
# source('demoniche_population.r')
# if (sum(n) > 0) {
# selects the right K for this pop and timeperiod
populationmax <-
BEMDEM$populationmax_all[n0s_ID[px],tx]
Projection[yx,,n0s_ID[px],tx] <-
demoniche_population(Matrix_projection = Matrix_projection, Matrix_projection_var = Matrix_projection_var,
n = n, populationmax = populationmax, onepopulation_Niche = population_Niche_short[px],
sumweight = BEMDEM$sumweight, Kweight = BEMDEM$Kweight, prob_scenario = BEMDEM$prob_scenario, noise = BEMDEM$noise, prev_mx = prev_mx,
transition_affected_demogr = BEMDEM$transition_affected_demogr, transition_affected_niche = BEMDEM$transition_affected_niche,
transition_affected_env = BEMDEM$transition_affected_env, env_stochas_type = BEMDEM$env_stochas_type, yx_tx = yx_tx)
# } # end if n > 0
} # end if n0s > 0
} # end px for
# Which many patches persist since last timestep (including seeds)?
metapop_results[yx_tx,mx,rx] <-
length(intersect( which(colSums(Projection[yx,,,tx]) > 1), n0s_ID))
################################################################################
##################### DISPERSAL FOR ALL PATCHES ####################
# Calculate dispersal for one repetition (rx), one matrix (mx), one time period(tx), and one year
# and all populationes.
# load('/noCC_nodispersal/Projection1_meanmatrix.rda')
# print(paste(Projection[yx,1,,tx]))
if(sum(Projection[yx,1,,tx]) > 0) # break
{
if(Dispersal == TRUE)
{
if (Niche == TRUE){
population_Niche <- BEMDEM$Niche_values[, tx]
}
# source("demoniche_dispersal.r")
disp <-
demoniche_dispersal(seeds_per_population = Projection[yx,1,,tx],
fraction_LDD = BEMDEM$fraction_LDD,
dispersal_probabilities = BEMDEM$dispersal_probabilities,
dist_latlong = BEMDEM$dist_latlong, neigh_index = BEMDEM$neigh_index,
fraction_SDD = BEMDEM$fraction_SDD) # T
Projection[yx,1,,tx] <- disp
} # end Dispersal if
} # end Dispersal if
################################################################################
# save the results from one run
population_sizes[yx_tx,mx,rx] <-
sum(rowSums(Projection[yx,,,tx]) * BEMDEM$sumweight)
} # end yx
EMA[rx,mx,tx,1] <- # the minimum population each repetition
min(apply((Projection[,,,tx] * BEMDEM$sumweight), 1, sum))
EMA[rx,mx,tx,2] <- # the number of populations that exist each repetition
sum(colSums(Projection[yx,,,tx]) > 1)
simulation_results[mx, 7+tx] <- mean(EMA[, mx, tx, 1])
} # end tx loop
save(Projection, file = paste(getwd(), "/", foldername,
"/Projection_rep", rx, "_", BEMDEM$list_names_matrices[mx], ".rda", sep = ""))
save(population_sizes, file = paste(getwd(), "/", foldername,
"/population_sizes", ".rda", sep = ""))
# PLOTS of spatial occupancy from the last repetition yx = 1
pop <- data.frame(cbind(BEMDEM$Niche_ID[,2:3],
(colSums(Projection[yx,,,] * BEMDEM$sumweight))))
form <- as.formula(paste(paste(colnames(pop)[-c(1:2)],collapse="+"),"X+Y",sep="~"))
jpeg(filename = paste(getwd(), "/", foldername, "/map_",
BEMDEM$list_names_matrices[mx], ".jpeg", sep = ""))
print(levelplot(form, pop, col.regions=rev(heat.colors(100)), allow.multiple = TRUE,
main = paste("Distribution", foldername, BEMDEM$list_names_matrices[mx], sep = "_")))
dev.off()
} # end mx loop
} # end rx loop
rm(Projection)
#### END OF MODELLING #################################################
####################################################################################
# extra loop to calculate mean, eigenvalues, etc
print("Calculating summary values", quote = FALSE)
for(mx in 1:length(BEMDEM$list_names_matrices)) # mx = 1
{
# for 'eigen_results' object for each matrix.
eigen_results[[mx]] <- c(eigen.analysis(matrix(BEMDEM$matrices[, mx],
ncol = length(BEMDEM$stages), byrow = FALSE)),
LTRE = list(matrix(LTRE(matrix(BEMDEM$matrices[, mx],
ncol = length(BEMDEM$stages), byrow = FALSE),
matrix(BEMDEM$matrices[, 1],
ncol = length(BEMDEM$stages), byrow = FALSE)),
nrow = length(BEMDEM$stages))))
simulation_results[mx,"lambda"] <- eigen_results[[mx]]$lambda1
simulation_results[mx,"stoch_lambda"] <-
stoch.growth.rate(list(matrix(BEMDEM$matrices[, 1],
ncol = length(BEMDEM$stages)), matrix(BEMDEM$matrices[, mx],
ncol = length(BEMDEM$stages))), prob = NULL, maxt = 50000,
verbose = FALSE)$sim
simulation_results[mx, "mean_final_pop"] <-
mean(population_sizes[yrs_total,mx,])
# Final mean percentage of patches extinct
# simulation_results[mx, "mean_perc_ext_final"] <- NA
# (metapop_results[1,mx,] - mean(metapop_results[yrs_total,mx,]))/metapop_results[1,mx,]*100
simulation_results[mx, "mean_no_patches_final"] <-
mean(EMA[,mx,length(BEMDEM$years_projections),2])
for(yx_tx in 1:yrs_total)
{
# mean of all repetitions
population_results[yx_tx, "Meanpop", mx] <- mean(population_sizes[yx_tx,mx,])
# sd of all repetitions
population_results[yx_tx, "SD", mx] <- sd(population_sizes[yx_tx,mx,])
# Min, minimum abundance of all simualtions each year.
population_results[yx_tx, "Min", mx] <- min(population_sizes[yx_tx,mx,])
# Max, maximum abundance of all simualtions each year.
population_results[yx_tx, "Max", mx] <- max(population_sizes[yx_tx,mx,])
}
# jpeg(file = paste(getwd(), "/", foldername, "/mean_",
# BEMDEM$list_names_matrices[mx], ".jpeg", sep = ""))
# plot(population_results[,"Meanpop",mx], ylim = c(0, max(population_results[,"Meanpop",mx])),
# main = paste("meanpopulation", BEMDEM$list_names_matrices[mx], foldername, sep = "_"),type = "l",
# xlab = "Year", ylab = "Population")
# dev.off()
} # end mx loop number 2, for eigenvalues.
######## PLOTS ################
# plot EMAS
jpeg(filename = paste(getwd(), "/", foldername, "/EMAs.jpeg", sep = ""),
width = 580, height = 480)
# matplot(t(simulation_results[, 8:(7+length(BEMDEM$years_projections))]), pch = 15, type = "l")
for(mx in 1:length(BEMDEM$list_names_matrices))
{
par(mar = c(7, 7, 4, 2) + 0.1, cex = 1.5)
plot(simulation_results[mx, 8:(7+length(BEMDEM$years_projections))],
type = 'b', ylim =
range(simulation_results[, 8:c(7+length(BEMDEM$years_projections))])
,
col = mx, xlab = "", ylab = "", axes = FALSE)
par(new = TRUE)
} # end mx loop number 3, for plotting EMAs
axis(1, at = 1:length(BEMDEM$years_projections), labels = FALSE)
text(1:length(BEMDEM$years_projections), y = par("usr")[3] - 5, srt = 45, adj = 1,
labels = (BEMDEM$years_projections), xpd = TRUE)
mtext(1, text = "Time", line = 6, cex = 1.7)
mtext(2, text = "EMA (number of individuals)", line = 6, cex = 1.7)
axis(2, xpd = TRUE, las = 1)
legend("topright", legend = BEMDEM$list_names_matrices, col = 1:mx, fill = 1:mx)
title("EMA for different matrices", cex = 0.9)
dev.off()
#### plot population results
jpeg(filename = paste(getwd(), "/", foldername, "/population_results.jpeg", sep = ""),
width = 580, height = 480)
for(mx in 1:length(BEMDEM$list_names_matrices))
{
par(mar = c(7, 4, 4, 2) + 0.1)
plot(population_results[,"Meanpop",mx], type = 'l', ylim = range(population_results[,1,]),
col = mx, xlab = "", ylab = "EMA", axes = TRUE, lwd = 1.2)
par(new = TRUE)
plot(c(population_results[,"Meanpop",mx] + population_results[,"SD",mx]), type = 'l', lty = 2,
ylim = range(population_results[,1,]),
col = mx, xlab = "", ylab = "EMA", axes = FALSE)
par(new = TRUE)
plot(c(population_results[,"Meanpop",mx] - population_results[,"SD",mx]), type = 'l', lty = 2,
ylim = range(population_results[,1,]),
col = mx, xlab = "", ylab = "EMA", axes = FALSE)
par(new = TRUE)
} # end mx loop number 3, for plotting EMAs
legend("topright", legend = BEMDEM$list_names_matrices, col = 1:mx, fill = 1:mx)
title("Population sizes (+- 1 SD) for different scenarios")
dev.off()
#########################################################
write.table(simulation_results, sep = ";", file =
paste(getwd(), "/", foldername, "/population_results.csv", sep = ""), col.names = NA)
save(simulation_results, file =
paste(getwd(), "/", foldername, "/population_sizes.rda", sep = ""))
save(metapop_results, file =
paste(getwd(), "/", foldername, "/metapop_results.rda", sep = ""))
save(eigen_results, file =
paste(getwd(), "/", foldername, "/eigen_results.rda", sep = ""))
save(EMA, file =
paste(getwd(), "/", foldername, "/EMA.rda", sep = ""))
save(population_results, file = paste(getwd(), "/",
foldername, "/population_results.rda", sep = ""))
print(" All repetitions completed!", quote = FALSE)
return(population_results)
} # end fcn
Try the demoniche package in your browser
Any scripts or data that you put into this service are public.
demoniche documentation built on May 2, 2019, 4:45 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
demoniche_model <-
function(modelname, Niche, Dispersal, repetitions, foldername)
{
BEMDEM <- get(modelname, envir = .GlobalEnv) #assign("BEMDEM", get(modelname))
#load(paste(modelname, ".rda", sep = ""))
# rm(Hmontana)
# require(sp)
require(popbio)
require(lattice)
#### Create vectors ############################################################
Projection <- array(0, # Projection is where all the results go
dim =
c(BEMDEM$no_yrs, length(BEMDEM$stages), nrow(BEMDEM$Niche_ID), length(BEMDEM$years_projections))
,
dimnames =
list(paste( "timesliceyear", 1:BEMDEM$no_yrs, sep="_"), c(paste(BEMDEM$stages)),
BEMDEM$Niche_ID[,"Niche_ID"], paste(BEMDEM$years_projections))
)
eigen_results <- vector(mode = "list" , length(BEMDEM$list_names_matrices))
names(eigen_results) <- unlist(BEMDEM$list_names_matrices)
yrs_total <- BEMDEM$no_yrs * length(BEMDEM$years_projections)
population_sizes <- array(NA, # per matrix!
dim = c(yrs_total, length(BEMDEM$list_names_matrices), repetitions),
dimnames = list(paste("year", 1: yrs_total, sep = ""), BEMDEM$list_names_matrices, paste("rep", 1:repetitions, sep = "_")))
population_results <- array(1:200, # for all matrices!
dim = c(yrs_total, 4, length(BEMDEM$list_names_matrices)),
dimnames = list( paste("year", 1: yrs_total, sep = ""),
c("Meanpop", "SD", "Max", "Min"),
paste(BEMDEM$list_names_matrices)))
metapop_results <- array(NA, # per matrix!
dim = c(yrs_total, length(BEMDEM$list_names_matrices), repetitions),
dimnames = list(paste("year", 1: yrs_total, sep = ""), BEMDEM$list_names_matrices, paste("rep", 1:repetitions, sep = "_")))
simulation_results <- array(NA,
dim = c(length(BEMDEM$list_names_matrices), 7+length(BEMDEM$years_projections)),
dimnames =
list(BEMDEM$list_names_matrices,
c("lambda", "stoch_lambda","mean_perc_ext_final", "initial_population_area", "initial_population", "mean_final_pop", "mean_no_patches_final",
paste("EMA", BEMDEM$years_projections))
) )
EMA <- array(0,
dim = c(repetitions, length(BEMDEM$list_names_matrices), length(BEMDEM$years_projections), 2), dimnames =
list(paste("rep", 1:repetitions, sep = "_"), BEMDEM$list_names_matrices, BEMDEM$years_projections, c("EMA", "No_populations")))
population_Niche <- rep(1, nrow(BEMDEM$Niche_ID))
simulation_results[,"initial_population_area"] <-
sum(BEMDEM$Orig_Populations[,"area_population"])
# original population size
simulation_results[,"initial_population"] <-
round(sum(colSums(BEMDEM$n0_all) * BEMDEM$sumweight), 0)
dir.create(paste(getwd(), "/", foldername, sep = ""), showWarnings = FALSE)
#### Start repetitions #########################################################
for (rx in 1:repetitions) # tx = 1 rx = 1
{
print(paste("Starting projections for repetition:", rx), quote = FALSE)
for(mx in 1:length(BEMDEM$list_names_matrices)) # selects two matrices for the simulations mx =1
{
print(paste("Projecting for scenario/matrix:", (BEMDEM$list_names_matrices)[mx]), quote = FALSE)
yx_tx <- 0 # restart counter
# this selects two matrices, the first basic matrix and the projection one mx =1
Matrix_projection <- cbind(BEMDEM$matrices[,1], (BEMDEM$matrices[, mx]))
# this selects two standard deviation matrices, the first basic matrix and the projection one (column mx)
if(BEMDEM$matrices_var[1] != FALSE) # If matrices are included or not.
{
if(ncol(BEMDEM$matrices_var) > 1){
Matrix_projection_var <- cbind(BEMDEM$matrices_var[,1], (BEMDEM$matrices_var[, mx]))
} else {
Matrix_projection_var <- cbind(BEMDEM$matrices_var[,1], (BEMDEM$matrices_var[, 1]))
}
} else
{
Matrix_projection_var <- FALSE
}
prev_mx <- rep(1, times = yrs_total + 1)
for(tx in 1:length(BEMDEM$years_projections)) # selects which time-slice of the simulation tx = 1
{
if (Niche == TRUE){ # Niche values
population_Niche <- BEMDEM$Niche_values[, tx]
}
for(yx in 1:BEMDEM$no_yrs)
{
yx_tx <- yx_tx + 1 # add one year to counter
####################################### PATCH PROJECTION when tx = 1 ########### tx=1 px=1 tx = 2
################################################################################ yx=2 yx = 1
if(tx == 1 && yx == 1)
{
n0s <- BEMDEM$n0_all[rowSums(BEMDEM$n0_all) > 0,] # take stage vector from intial population, where over zero
n0s_ID <- which(rowSums(BEMDEM$n0_all) > 0) # which populations/rows are above zero
} else { # Second year running the model
if(tx != 1 && yx == 1){ # when going to next time step, yx == 1
n0s <- t(Projection[BEMDEM$no_yrs,, colSums(Projection[BEMDEM$no_yrs,,,tx-1]) > 0, tx-1])
n0s_ID <- which(colSums(Projection[BEMDEM$no_yrs,,,tx-1]) > 0)
BEMDEM$Niche_ID[colSums(Projection[BEMDEM$no_yrs,,,tx-1]) > 0, 2]
} else {
# take population from previous year, same time period
n0s <- t(Projection[yx-1,, colSums(Projection[yx-1,,,tx]) > 0,tx])
n0s_ID <- which(colSums(Projection[yx-1,,,tx]) > 0)
}
} # end if tx == 1 && yx == 1
# n <- c(10,5,5,5,5,5)
population_Niche_short <- population_Niche[n0s_ID]
# run fcn for each population px=1
if (nrow(n0s) > 0) {
for(px in 1:nrow(n0s)) { # px = 1 tx =1
n <- as.vector(n0s[px,])
# source('demoniche_population.r')
# if (sum(n) > 0) {
# selects the right K for this pop and timeperiod
populationmax <-
BEMDEM$populationmax_all[n0s_ID[px],tx]
Projection[yx,,n0s_ID[px],tx] <-
demoniche_population(Matrix_projection = Matrix_projection, Matrix_projection_var = Matrix_projection_var,
n = n, populationmax = populationmax, onepopulation_Niche = population_Niche_short[px],
sumweight = BEMDEM$sumweight, Kweight = BEMDEM$Kweight, prob_scenario = BEMDEM$prob_scenario, noise = BEMDEM$noise, prev_mx = prev_mx,
transition_affected_demogr = BEMDEM$transition_affected_demogr, transition_affected_niche = BEMDEM$transition_affected_niche,
transition_affected_env = BEMDEM$transition_affected_env, env_stochas_type = BEMDEM$env_stochas_type, yx_tx = yx_tx)
# } # end if n > 0
} # end if n0s > 0
} # end px for
# Which many patches persist since last timestep (including seeds)?
metapop_results[yx_tx,mx,rx] <-
length(intersect( which(colSums(Projection[yx,,,tx]) > 1), n0s_ID))
################################################################################
##################### DISPERSAL FOR ALL PATCHES ####################
# Calculate dispersal for one repetition (rx), one matrix (mx), one time period(tx), and one year
# and all populationes.
# load('/noCC_nodispersal/Projection1_meanmatrix.rda')
# print(paste(Projection[yx,1,,tx]))
if(sum(Projection[yx,1,,tx]) > 0) # break
{
if(Dispersal == TRUE)
{
if (Niche == TRUE){
population_Niche <- BEMDEM$Niche_values[, tx]
}
# source("demoniche_dispersal.r")
disp <-
demoniche_dispersal(seeds_per_population = Projection[yx,1,,tx],
fraction_LDD = BEMDEM$fraction_LDD,
dispersal_probabilities = BEMDEM$dispersal_probabilities,
dist_latlong = BEMDEM$dist_latlong, neigh_index = BEMDEM$neigh_index,
fraction_SDD = BEMDEM$fraction_SDD) # T
Projection[yx,1,,tx] <- disp
} # end Dispersal if
} # end Dispersal if
################################################################################
# save the results from one run
population_sizes[yx_tx,mx,rx] <-
sum(rowSums(Projection[yx,,,tx]) * BEMDEM$sumweight)
} # end yx
EMA[rx,mx,tx,1] <- # the minimum population each repetition
min(apply((Projection[,,,tx] * BEMDEM$sumweight), 1, sum))
EMA[rx,mx,tx,2] <- # the number of populations that exist each repetition
sum(colSums(Projection[yx,,,tx]) > 1)
simulation_results[mx, 7+tx] <- mean(EMA[, mx, tx, 1])
} # end tx loop
save(Projection, file = paste(getwd(), "/", foldername,
"/Projection_rep", rx, "_", BEMDEM$list_names_matrices[mx], ".rda", sep = ""))
save(population_sizes, file = paste(getwd(), "/", foldername,
"/population_sizes", ".rda", sep = ""))
# PLOTS of spatial occupancy from the last repetition yx = 1
pop <- data.frame(cbind(BEMDEM$Niche_ID[,2:3],
(colSums(Projection[yx,,,] * BEMDEM$sumweight))))
form <- as.formula(paste(paste(colnames(pop)[-c(1:2)],collapse="+"),"X+Y",sep="~"))
jpeg(filename = paste(getwd(), "/", foldername, "/map_",
BEMDEM$list_names_matrices[mx], ".jpeg", sep = ""))
print(levelplot(form, pop, col.regions=rev(heat.colors(100)), allow.multiple = TRUE,
main = paste("Distribution", foldername, BEMDEM$list_names_matrices[mx], sep = "_")))
dev.off()
} # end mx loop
} # end rx loop
rm(Projection)
#### END OF MODELLING #################################################
####################################################################################
# extra loop to calculate mean, eigenvalues, etc
print("Calculating summary values", quote = FALSE)
for(mx in 1:length(BEMDEM$list_names_matrices)) # mx = 1
{
# for 'eigen_results' object for each matrix.
eigen_results[[mx]] <- c(eigen.analysis(matrix(BEMDEM$matrices[, mx],
ncol = length(BEMDEM$stages), byrow = FALSE)),
LTRE = list(matrix(LTRE(matrix(BEMDEM$matrices[, mx],
ncol = length(BEMDEM$stages), byrow = FALSE),
matrix(BEMDEM$matrices[, 1],
ncol = length(BEMDEM$stages), byrow = FALSE)),
nrow = length(BEMDEM$stages))))
simulation_results[mx,"lambda"] <- eigen_results[[mx]]$lambda1
simulation_results[mx,"stoch_lambda"] <-
stoch.growth.rate(list(matrix(BEMDEM$matrices[, 1],
ncol = length(BEMDEM$stages)), matrix(BEMDEM$matrices[, mx],
ncol = length(BEMDEM$stages))), prob = NULL, maxt = 50000,
verbose = FALSE)$sim
simulation_results[mx, "mean_final_pop"] <-
mean(population_sizes[yrs_total,mx,])
# Final mean percentage of patches extinct
# simulation_results[mx, "mean_perc_ext_final"] <- NA
# (metapop_results[1,mx,] - mean(metapop_results[yrs_total,mx,]))/metapop_results[1,mx,]*100
simulation_results[mx, "mean_no_patches_final"] <-
mean(EMA[,mx,length(BEMDEM$years_projections),2])
for(yx_tx in 1:yrs_total)
{
# mean of all repetitions
population_results[yx_tx, "Meanpop", mx] <- mean(population_sizes[yx_tx,mx,])
# sd of all repetitions
population_results[yx_tx, "SD", mx] <- sd(population_sizes[yx_tx,mx,])
# Min, minimum abundance of all simualtions each year.
population_results[yx_tx, "Min", mx] <- min(population_sizes[yx_tx,mx,])
# Max, maximum abundance of all simualtions each year.
population_results[yx_tx, "Max", mx] <- max(population_sizes[yx_tx,mx,])
}
# jpeg(file = paste(getwd(), "/", foldername, "/mean_",
# BEMDEM$list_names_matrices[mx], ".jpeg", sep = ""))
# plot(population_results[,"Meanpop",mx], ylim = c(0, max(population_results[,"Meanpop",mx])),
# main = paste("meanpopulation", BEMDEM$list_names_matrices[mx], foldername, sep = "_"),type = "l",
# xlab = "Year", ylab = "Population")
# dev.off()
} # end mx loop number 2, for eigenvalues.
######## PLOTS ################
# plot EMAS
jpeg(filename = paste(getwd(), "/", foldername, "/EMAs.jpeg", sep = ""),
width = 580, height = 480)
# matplot(t(simulation_results[, 8:(7+length(BEMDEM$years_projections))]), pch = 15, type = "l")
for(mx in 1:length(BEMDEM$list_names_matrices))
{
par(mar = c(7, 7, 4, 2) + 0.1, cex = 1.5)
plot(simulation_results[mx, 8:(7+length(BEMDEM$years_projections))],
type = 'b', ylim =
range(simulation_results[, 8:c(7+length(BEMDEM$years_projections))])
,
col = mx, xlab = "", ylab = "", axes = FALSE)
par(new = TRUE)
} # end mx loop number 3, for plotting EMAs
axis(1, at = 1:length(BEMDEM$years_projections), labels = FALSE)
text(1:length(BEMDEM$years_projections), y = par("usr")[3] - 5, srt = 45, adj = 1,
labels = (BEMDEM$years_projections), xpd = TRUE)
mtext(1, text = "Time", line = 6, cex = 1.7)
mtext(2, text = "EMA (number of individuals)", line = 6, cex = 1.7)
axis(2, xpd = TRUE, las = 1)
legend("topright", legend = BEMDEM$list_names_matrices, col = 1:mx, fill = 1:mx)
title("EMA for different matrices", cex = 0.9)
dev.off()
#### plot population results
jpeg(filename = paste(getwd(), "/", foldername, "/population_results.jpeg", sep = ""),
width = 580, height = 480)
for(mx in 1:length(BEMDEM$list_names_matrices))
{
par(mar = c(7, 4, 4, 2) + 0.1)
plot(population_results[,"Meanpop",mx], type = 'l', ylim = range(population_results[,1,]),
col = mx, xlab = "", ylab = "EMA", axes = TRUE, lwd = 1.2)
par(new = TRUE)
plot(c(population_results[,"Meanpop",mx] + population_results[,"SD",mx]), type = 'l', lty = 2,
ylim = range(population_results[,1,]),
col = mx, xlab = "", ylab = "EMA", axes = FALSE)
par(new = TRUE)
plot(c(population_results[,"Meanpop",mx] - population_results[,"SD",mx]), type = 'l', lty = 2,
ylim = range(population_results[,1,]),
col = mx, xlab = "", ylab = "EMA", axes = FALSE)
par(new = TRUE)
} # end mx loop number 3, for plotting EMAs
legend("topright", legend = BEMDEM$list_names_matrices, col = 1:mx, fill = 1:mx)
title("Population sizes (+- 1 SD) for different scenarios")
dev.off()
#########################################################
write.table(simulation_results, sep = ";", file =
paste(getwd(), "/", foldername, "/population_results.csv", sep = ""), col.names = NA)
save(simulation_results, file =
paste(getwd(), "/", foldername, "/population_sizes.rda", sep = ""))
save(metapop_results, file =
paste(getwd(), "/", foldername, "/metapop_results.rda", sep = ""))
save(eigen_results, file =
paste(getwd(), "/", foldername, "/eigen_results.rda", sep = ""))
save(EMA, file =
paste(getwd(), "/", foldername, "/EMA.rda", sep = ""))
save(population_results, file = paste(getwd(), "/",
foldername, "/population_results.rda", sep = ""))
print(" All repetitions completed!", quote = FALSE)
return(population_results)
} # end fcn
Try the demoniche package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.