diff --git a/Manuscript-scripts/Appendix S9/Appendix S9_main text_plots.r b/Manuscript-scripts/Appendix S8/Appendix S8_main text_plots.r similarity index 92% rename from Manuscript-scripts/Appendix S9/Appendix S9_main text_plots.r rename to Manuscript-scripts/Appendix S8/Appendix S8_main text_plots.r index 3b1c1e2..c9b708d 100644 --- a/Manuscript-scripts/Appendix S9/Appendix S9_main text_plots.r +++ b/Manuscript-scripts/Appendix S8/Appendix S8_main text_plots.r @@ -9,7 +9,7 @@ # Either place this rdata file in the default working directory, or use setwd() # # to set a working directory. # # # -# Figs. 3 & 4 take about 10 minutes to run in 8 cores # +# Figs. 2 & S2 take about 10 minutes to run in 8 cores # ################################################################################# library(ks) # 2D kernel density @@ -24,8 +24,8 @@ step.length <- 10 # number of cores to use in parallel simulations; adjust as necessary n.cores <- 8 -if(!file.exists("Appendix S9_resistance-rasters.rdata")) { - stop("Can't find data file.\n************************************************\nPlease copy 'Appendix S9_resistance-rasters.rdata' to the folder\n", getwd()) +if(!file.exists("Appendix S8_resistance-rasters.rdata")) { + stop("Can't find data file.\n************************************************\nPlease copy 'Appendix S8_resistance-rasters.rdata' to the folder\n", getwd()) } ## Load rasters corresponding to how different animals might see the landscape: @@ -35,7 +35,7 @@ if(!file.exists("Appendix S9_resistance-rasters.rdata")) { # more often in forested areas. e.g. otter # - aquatic animal: Moves exclusively in water. e.g. fish -load("Appendix S9_resistance-rasters.rdata") +load("Appendix S8_resistance-rasters.rdata") # this is a list with 3 elements: terr, amph, fish, corresponding to the # descriptions above @@ -119,14 +119,14 @@ plot.relocs <- function(resist, relocs, n.rep = NA , xlim = NA, ylim = NA ## START MENU choice <- menu(c( - "Fig.3 - Influence of landscape" - , "Fig.4 - Influence of perceptual range" - , "Fig.2 - Basic plots" + "Fig.2 - Influence of landscape" + , "Fig.S2 - Influence of perceptual range" + , "Fig.1 - Basic plots" , "Just plot last run" )) switch(choice , { - output.filename <- "Figure-3-influence-landscape" + output.filename <- "Figure-2-influence-landscape" simulated.rasters <- rasters # create species levy.walker <- species((state.RW() * 100) + (state.CRW(0.95) * 500) @@ -173,7 +173,7 @@ switch(choice , { ) }, { - output.filename <- "Figure-4-influence-percep-range" + output.filename <- "Figure-S2-influence-percep-range" species <- lapply(c(5000, 2000, 500), function(pr) { species(list( state(0, perceptualRange("cir", pr), step.length) @@ -245,7 +245,7 @@ switch(choice , { relocs <- mapply(my.simulate, simulated.rasters, n.rep, species , list(c(691391, 4629949)), SIMPLIFY = FALSE) - tiff("Figure-2-basic-illustration.tif", width = 10 * 2, height = 10 * 2 + tiff("Figure-1-basic-illustration.tif", width = 10 * 2, height = 10 * 2 , unit = "cm", comp = "lzw", res = 1000) par(mar = c(0.1, 0.1, 2.5, 0.1)) layout(matrix(c(seq_along(relocs)), nrow = 2, byrow = TRUE)) @@ -266,7 +266,7 @@ switch(choice , { }) -## Make Fig.3 & 4, and the full res plots in appendix +## Make Fig.2 & S2, and the full res plots (not shown) if(choice %in% c(1, 2)) { cat("Computing kernel densities...\n"); flush.console() if(n.rep > 1) { diff --git a/Manuscript-scripts/Appendix S9/Appendix S9_resistance-rasters.rdata b/Manuscript-scripts/Appendix S8/Appendix S8_resistance-rasters.rdata similarity index 100% rename from Manuscript-scripts/Appendix S9/Appendix S9_resistance-rasters.rdata rename to Manuscript-scripts/Appendix S8/Appendix S8_resistance-rasters.rdata diff --git a/Manuscript-scripts/Appendix S10/Appendix S10_model_parameterization.r b/Manuscript-scripts/Appendix S9/Appendix S9_model_parameterization.r similarity index 91% rename from Manuscript-scripts/Appendix S10/Appendix S10_model_parameterization.r rename to Manuscript-scripts/Appendix S9/Appendix S9_model_parameterization.r index 38d6632..99275e9 100644 --- a/Manuscript-scripts/Appendix S10/Appendix S10_model_parameterization.r +++ b/Manuscript-scripts/Appendix S9/Appendix S9_model_parameterization.r @@ -1,6 +1,6 @@ ################################################################################ # EXAMPLES OF INPUT PARAMETER APPROXIMATION FROM REAL DATASETS # -# This script produces the plots presented in Appendix S2 of the paper # +# This script produces the plots presented in Appendix S3 & S4 of the paper # #------------------------------------------------------------------------------# # IMPORTANT NOTE: Run it with the command source('filename'), # # NOT with copy/paste! # @@ -12,19 +12,8 @@ library(adehabitatLT) # for L library(moveHMM) # for moveHMM simulation library(TeachingDemos) # for subplots inside plot -# number of generations to run optimization algorithm -# NOTE: we don't need 1000 generations here, as the algorithm generally converges -# faster. You can reduce this to e.g. 500 -generations <- seq(5, 1000, by = 5) # number of times each candidate solution is simulated to compute average fitness nrepetitions <- 6 -# simulate at 50 times higher frequency than real data -downsample <- 50 -# use 7-bin histograms for turning angle and step length ditributions during -# optimization, for solution fitness evaluation -nbins.hist <- c(7, 7, 0) -# use the logarithm of the step lengths for computing the histograms above -log.step.length <- TRUE # Just a custom function to plot trajectories and their histograms of turning # angles and step lengths inset at the corners @@ -38,7 +27,6 @@ movementplot <- function(relocs, downsample = 1, nbins = c(4, 5) plot(relocs, type = "l", asp = 1, axes = FALSE) } else { plot(relocs, type="l", asp=1, col="#aaaaaa", axes = FALSE) -# points(relocs, col=c("#aaaaaa", "red", "blue")[relocs[, 3] + 1], pch=19, cex=0.1) lines(relocs.st$relocs, lwd=0.5) } } else { @@ -82,6 +70,18 @@ switch(menu(c( "Optimization convergence plots with real and simulated data" , "Validation plots (ability to recover true parameter values after downsampling)" )), { + # number of generations to run optimization algorithm + # NOTE: we don't need 1000 generations here, as the algorithm generally converges + # faster. You can reduce this to e.g. 500 + generations <- seq(5, 1000, by = 5) + # simulate at 50 times higher frequency than real data + downsample <- 50 + # use 7-bin histograms for turning angle and step length ditributions during + # optimization, for solution fitness evaluation + nbins.hist <- c(7, 7, 0) + # use the logarithm of the step lengths for computing the histograms above + log.step.length <- TRUE + switch(menu(c( "Elk data (Morales et al., 2004)" , "moveHMM simulation" @@ -180,8 +180,8 @@ switch(menu(c( # NOTE: because the optimization algorithm is not yet fully optimized for speed, # running this procedure takes about 22 h! The code is provided just to illustrate # how to do it. - if(!file.exists("Appendix S10_otter-realdata.rdata")) { - stop("Can't find data file.\n************************************************\nPlease copy 'Appendix S10_otter-realdata.rdata' to the folder\n", getwd()) + if(!file.exists("Appendix S9_otter-realdata.rdata")) { + stop("Can't find data file.\n************************************************\nPlease copy 'Appendix S9_otter-realdata.rdata' to the folder\n", getwd()) } cat("NOTE: this options takes >20 hours to run, see script comments.") @@ -191,7 +191,7 @@ switch(menu(c( # equally well with 50 downsample <- 25 - load("Appendix S10_otter-realdata.rdata") + load("Appendix S9_otter-realdata.rdata") filename <- "otter" tmp <- sampleMovement(real.data) max.step.length <- (max(tmp$stat[, "steplengths"]) / downsample) * 2 @@ -222,7 +222,7 @@ switch(menu(c( , resistance = resistance, coords = startCoord , resol = downsample, generations = generations , nrepetitions = nrepetitions, nbins.hist = nbins.hist - , step.hist.log = log.step.length, trace = TRUE, TA.variation = FALSE + , step.hist.log = log.step.length, trace = TRUE , aggregate = TRUE) }) # ... And that's all it takes to approximate parameters from real data. @@ -319,10 +319,13 @@ switch(menu(c( # Is the optimization able to recover true parameters after a 1:50 downsampling? downsample <- 50 - # Use 5-bin histograms for turning angle distributions and 12-bin histograms + # Use 5-bin histograms for turning angle distributions and 8-bin histograms # for step length ditributions during optimization, for solution fitness evaluation nbins.hist <- c(5, 8, 0) + # use the logarithm of the step lengths for computing the histograms above + log.step.length <- TRUE + # define simulation parameters that we'll use for validation # one state CRW, two state RW-CRW and three state RW-CRW-CRW simulation.parameters <- list( @@ -340,7 +343,7 @@ switch(menu(c( # note that the performance of the optimization in recovering the true parameters # depends on the particular trajectory that was used as "real data". Hence, multiple # runs may produce different validation plots. - for(sp in seq_along(species.models)[3]) { + for(sp in seq_along(species.models)) { smodel <- species.models[sp] truepars <- simulation.parameters[[sp]] diff --git a/Manuscript-scripts/Appendix S10/Appendix S10_otter-realdata.rdata b/Manuscript-scripts/Appendix S9/Appendix S9_otter-realdata.rdata similarity index 100% rename from Manuscript-scripts/Appendix S10/Appendix S10_otter-realdata.rdata rename to Manuscript-scripts/Appendix S9/Appendix S9_otter-realdata.rdata