# read shapefile and convert into a Polygon
bolsachica <- sf::st_read(system.file("shapes/FTB_lines.shp", package="animalEKF")[1])
bolsachica <- sf::st_polygonize(bolsachica)
island <- sf::st_read(system.file("shapes/FTB_island.shp", package="animalEKF")[1])
# the actual available room for movement is the area in the water, subtracting the island inside
bolsachica <- sf::st_difference(sf::st_geometry(bolsachica), sf::st_geometry(island))
# sample 5 points approximately equally spaced within the shapefile, as region centroids
regions <- sf::st_sample(x=sf::st_geometry(bolsachica), size=5, type="regular", exact=TRUE)
# extract the coordinates
regions <- as.data.frame(sf::st_coordinates(regions)[, c("X","Y")])
#define Voronoi tessellation tile in which to start shark paths
vortess <- deldir::deldir(x=regions[,"X"], y=regions[,"Y"], wlines="tess",
plotit=FALSE, suppressMsge=TRUE)
# convert these to a set of Polygons, and choose one of them as the starting polygon
vtiles <- sf::st_as_sf(tess2spat(vortess))
sf::st_crs(vtiles) <- sf::st_crs(bolsachica)
# extract only the 3rd tile
# note, want to have the simulation paths spread out, so a given draw may result in
# cramped and thus hard to estimate paths
starting_polygon <- sf::st_sfc(vtiles[[1]][[3]], crs=sf::st_crs(vtiles))
starting_polygon <- sf::st_intersection(sf::st_geometry(bolsachica),
sf::st_geometry(starting_polygon))
#define list of transition matrices between behaviors
tmat_list <- list(matrix(c(8, 2, 2, 4), ncol=2, byrow=TRUE),
matrix(c(1.5*5, 1.5*1, 3, 3), ncol=2, byrow=TRUE),
matrix(c(7, 1, 1, 7), ncol=2, byrow=TRUE))
#generate 4-shark simulated trajectory with 100 regular steps of length 120 seconds.
#Sharks 3 and 4 will be interacting with the others, but 1 and 2 will not.
nsharks <- 4
#simulate trajectory
#setting gen_irreg=TRUE generates an irregular trajectory from the regular-step one
#with the log-normal specified in dt_lnorm_mu and dt_lnorm_sd
#sim_4sharks$di would contain the irregular dataset
#otherwise, say you wanted to try different interpolations, you can use the same regular
#step from sim_trajectory_joint and then interpolate separately with interp_trajectory_joint.
#make simulated trajectories all start in the same area so they will be close enough to be
#interacting, for the purposes of this exercise
#note that the simulation may time out trying to draw points in this starting polygon that end
#up in the shapefile boundary
nsteps_sim <- 100
reg_dt <- 120
sim_4sharks <- sim_trajectory_joint(area_map=sf::st_geometry(bolsachica), centroids=regions,
transition_matrices=tmat_list, nsharks=nsharks,
mu0_pars=list(alpha=c(-4 ,-1.6), beta=c(0,0)),
var0_pars=list(alpha=c(1,0.25), beta=c(1,.25)),
N=nsteps_sim, nstates=2, reg_dt=reg_dt,
gen_irreg=FALSE, one_d=FALSE,
starting_polygon=starting_polygon, interact=TRUE,
interact_pars=list(interacting_sharks=c(3:4),
time_radius=60*30, spat_radius=150,
min_num_neibs=10,
eta_mu=c(2,1), rho_sd=c(0.75, 0.75)),
time_dep_trans=FALSE,
dt_lnorm_mu=log(120), dt_lnorm_sd=0.4)
#plot trajectories
shark_names <- dimnames(sim_4sharks$d)[[ 3 ]]
shark_colors <- 2:5
names(shark_colors) <- shark_names
sp::plot(bolsachica, main="Full trajectories")
deldir::plot.deldir(vortess, wlines="tess", add=TRUE)
for (ss in shark_names) {
lines(sim_4sharks$d[,c("X","Y"), ss], col=shark_colors[ss])
}
#now interpolate to uneven steps with lognormal mean log(120) (so they are on
#average the same as the regular steps and sd=0.4
#d is the regular step, di is irregular
#if want to interpolate separately. Otherwise just set gen_irreg=TRUE above
#this is so you can interpolate a dataset not generated by sim_trajectory_joint
#if gen_irreg=TRUE in sim_trajectory_joint,
#interp_ds will be returned as the 'di' object
interp_ds <- interp_trajectory_joint(d=sim_4sharks$d, nstates=2,
one_d=FALSE,
dt_lnorm_mu=log(reg_dt),
dt_lnorm_sd=0.4,
centroids=regions)
#now plot observed ones, may differ
sp::plot(bolsachica, main="Observed trajectories")
deldir::plot.deldir(vortess, wlines="tess", add=TRUE)
for (ss in shark_names) {
lines(interp_ds[ interp_ds$tag == ss ,c("X","Y")], col=shark_colors[ss])
}
# \donttest{
#try to recover EKF with steps at the original 120 seconds
#use the original simulated transition and foraging probabilities for comparison
#intial values for some parameters
tau_pars_init <- c(8, 14, 10,1) #2
sigma_pars_init <- c(5, 8, 8, 3)
#measurement error
bmat <- matrix(c(1, -0.3, -0.3, 1), ncol=2)
Errvar_init1 <-5*20*bmat
Errvar_init2 <- 15*20*bmat
#particle error
Particle_err_init <- 0.5*20*bmat
# only estimate movement on first 5 steps
# for better results, npart should be set higher, like 150 or more
nsteps_estimate <- 5
npart <- 15
#again, if you use gen_irreg=TRUE in sim_trajectory_joint,
#the input 'd' argument should be sim_4sharks$di or interp_ds
#NOTE: user should set output_plot=TRUE to see PDF,
#for purposes of package testing we set it to FALSE
# if FALSE, plots will still appear in the console
ekf_interp_mod <- EKF_interp_joint(d=interp_ds, npart=npart,
area_map=bolsachica,
state_favor=c(1,2),
centroids=regions,
sigma_pars=sigma_pars_init,
tau_pars=tau_pars_init,
Errvar0=list(Errvar_init1, Errvar_init2),
Particle_errvar0=Particle_err_init,
mu0_pars=list(alpha=c(-4 ,-1.3), beta=c(0,0)),
truncate=TRUE,
neff_sample=0.75, dirichlet_init=c(8,2,2,4),
smoothing=TRUE, fix_smoothed_behaviors=FALSE,
time_dep_trans=FALSE, resamp_full_hist=FALSE,
nstates=2, reg_dt=reg_dt, interact=TRUE,
maxStep=nsteps_estimate, update_eachstep=TRUE,
compare_with_known=TRUE,
known_trans_prob=sim_4sharks$true_transition_prob,
known_foraging_prob=sim_4sharks$true_foraging_prob,
known_regular_step_ds=sim_4sharks$d_ds,
output_plot=FALSE)
#simulate one-dimensional movement for 1 robot (shark)
#here we use gen_irreg=TRUE instead of generating a separate interpolation object
one_d <- sim_trajectory_joint(centroids=NULL, N=nsteps_sim,
mu0_pars=list(alpha=c(4, 9)),
var0_pars=list(alpha=c(1, 1)),
transition_matrices=tmat_list[[ 1 ]], nstates=2,
reg_dt=reg_dt, gen_irreg=TRUE, one_d=TRUE,
dt_lnorm_mu=log(120), dt_lnorm_sd=0.55)
#measurement error
bmat <- matrix(1)
Errvar_init1 <-1*bmat
Errvar_init2 <-3*bmat
#particle error
Particle_err_init <- 0.1*bmat
ekf_1d <- EKF_1d_interp_joint(d=one_d$di, npart=npart, maxStep=nsteps_estimate,
state_favor=c(1,1), nstates=2, lowvarsample=TRUE,
neff_sample=1, time_dep_trans=FALSE, reg_dt=reg_dt,
max_int_wo_obs=15, resamp_full_hist=FALSE,
alpha0_pars=list(mu0=c(4, 9), V0=c(0.25, 0.25)),
sigma_pars=sigma_pars_init,
Errvar0=list(Errvar_init1, Errvar_init2),
Particle_errvar0=Particle_err_init,
compare_with_known=TRUE,
known_trans_prob=one_d$true_transition_prob,
known_foraging_prob=one_d$true_foraging_prob,
known_regular_step_ds=one_d$d_ds, update_eachstep=TRUE,
smoothing=TRUE, output_plot=FALSE)
# }
Run the code above in your browser using DataLab