Skip to content

Commit

Permalink
Updated manuscript scripts.
Browse files Browse the repository at this point in the history
  • Loading branch information
miguel-porto committed Jan 31, 2018
1 parent 0cde278 commit fa2d35d
Show file tree
Hide file tree
Showing 4 changed files with 33 additions and 30 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
Expand All @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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))
Expand All @@ -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) {
Expand Down
Original file line number Diff line number Diff line change
@@ -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! #
Expand All @@ -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
Expand All @@ -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 {
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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.")
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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(
Expand All @@ -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]]

Expand Down

0 comments on commit fa2d35d

Please sign in to comment.