R Under development (unstable) (2024-05-20 r86569 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #remotes::install_github("psolymos/bSims") > > library(bSims) Loading required package: intrval Loading required package: mefa4 Loading required package: Matrix mefa4 0.3-11 2024-05-20 Attaching package: 'mefa4' The following object is masked from 'package:intrval': %notin% Loading required package: MASS Loading required package: deldir deldir 2.0-4 Nickname: "Idol Comparison" The syntax of deldir() has changed since version 0.0-10. Read the help!!!. bSims 0.3-2 2024-05-21 chip-chewy > > help_pages <- c( + "acceptreject", + "bsims_init", + "dist_fun2", + "events", + "expand_list", + "get_nests", + "rlnorm2", + "rmvn") > > for (i in help_pages) { + cat("\n\n---------- bSims example:", i, "----------\n\n") + eval(parse(text=paste0("example('", i, + "', package='bSims', run.dontrun=TRUE, run.donttest=TRUE)"))) + } ---------- bSims example: acceptreject ---------- accptr> ## complete spatial randomness accptr> plot(acceptreject(100), asp=1) accptr> ## more systematic accptr> distance <- seq(0,1,0.01) accptr> f <- function(d) accptr+ (1-exp(-d^2/0.1^2) + dlnorm(d, 0.2)/dlnorm(exp(0.2-1),0.2)) / 2 accptr> op <- par(mfrow = c(1, 2)) accptr> plot(distance, f(distance), type="l") accptr> plot(acceptreject(100, f, m=1), asp=1) accptr> par(op) accptr> ## more clustered accptr> f <- function(d) accptr+ exp(-d^2/0.1^2) + 0.5*(1-exp(-d^2/0.4^2)) accptr> op <- par(mfrow = c(1, 2)) accptr> plot(distance, f(distance), type="l") accptr> plot(acceptreject(100, f, m=1), asp=1) accptr> par(op) ---------- bSims example: bsims_init ---------- bsms_n> phi <- 0.5 bsms_n> tau <- 1:3 bsms_n> dur <- 10 bsms_n> rbr <- c(0.5, 1, 1.5, Inf) bsms_n> tbr <- c(3, 5, 10) bsms_n> (l <- bsims_init(10, 0.5, 1)) bSims landscape 1 km x 1 km stratification: HER bsms_n> (p <- bsims_populate(l, 1)) bSims population 1 km x 1 km stratification: HER total abundance: 103 bsms_n> (a <- bsims_animate(p, vocal_rate=phi, duration=dur)) bSims events 1 km x 1 km stratification: HER total abundance: 103 duration: 10 min bsms_n> (o <- bsims_detect(a, tau=tau)) bSims detections 1 km x 1 km stratification: HER total abundance: 103 duration: 10 min detected: 36 heard bsms_n> (x <- bsims_transcribe(o, tint=tbr, rint=rbr)) bSims transcript 1 km x 1 km stratification: HER total abundance: 103 duration: 10 min detected: 36 heard 1st event detected by breaks: [0, 3, 5, 10 min] [0, 50, 100, 150, Inf m] bsms_n> plot(x) bsms_n> get_table(x, "removal") 0-3min 3-5min 5-10min 0-50m 2 0 0 50-100m 2 2 0 100-150m 0 1 0 150+m 7 1 3 bsms_n> get_table(x, "visits") 0-3min 3-5min 5-10min 0-50m 2 1 2 50-100m 2 3 2 100-150m 0 2 2 150+m 7 6 14 bsms_n> head(get_events(a)) x y t v a i 1 4.4005662 -3.8552792 0.009417043 1 20 79 2 0.3751943 -3.2730626 0.093780628 1 149 56 3 -3.7321867 3.9509588 0.107290643 1 281 34 4 4.2688294 0.4080248 0.176989985 1 315 82 5 2.3275329 -4.2081446 0.177479506 1 104 76 6 0.8343926 2.8678146 0.185008359 1 108 64 bsms_n> plot(get_events(a)) bsms_n> head(get_detections(o)) x y t v a d f i j 27 -1.59181186 1.8860926 0.5395908 1 56 2.4680378 NA 2 2 33 0.37519431 -3.2730626 0.6932617 1 140 3.2944968 NA 56 56 35 0.38709520 0.5869883 0.7327568 1 255 0.7031344 NA 48 48 36 -0.28239990 0.7625467 0.7385949 1 139 0.8131588 NA 58 58 63 0.07796123 -2.2904200 1.3754761 1 221 2.2917465 NA 46 46 68 0.76679501 1.4415011 1.4420373 1 123 1.6327584 NA 59 59 bsms_n> plot(get_detections(o), "time") bsms_n> plot(get_detections(o), "distance") bsms_n> ## wrapper function for all the bsims_* layers bsms_n> b <- bsims_all(road=1, density=0.5, tint=tbr, rint=rbr) bsms_n> ## alternatively: supply a list bsms_n> #settings <- list(road=1, density=0.5, tint=tbr, rint=rbr) bsms_n> #b <- bsims_all(settings) bsms_n> b$settings() $road [1] 1 $density [1] 0.5 $tint [1] 3 5 10 $rint [1] 0.5 1.0 1.5 Inf bsms_n> b$new() bSims transcript 1 km x 1 km stratification: HR total abundance: 46 duration: 10 min detected: 4 heard 1st event detected by breaks: [0, 3, 5, 10 min] [0, 50, 100, 150, Inf m] bsms_n> bb <- b$replicate(3) bsms_n> lapply(bb, get_table) [[1]] 0-3min 3-5min 5-10min 0-50m 1 0 0 50-100m 2 0 0 100-150m 0 0 0 150+m 1 0 0 [[2]] 0-3min 3-5min 5-10min 0-50m 2 0 0 50-100m 2 0 0 100-150m 0 0 0 150+m 0 0 0 [[3]] 0-3min 3-5min 5-10min 0-50m 0 0 0 50-100m 0 0 0 100-150m 0 0 0 150+m 0 0 0 bsms_n> ## No test: bsms_n> ## parallel simulations bsms_n> library(parallel) bsms_n> b <- bsims_all(density=0.5) bsms_n> B <- 4 # number of runs bsms_n> nc <- 2 # number of cores bsms_n> ## sequential bsms_n> system.time(bb <- b$replicate(B, cl=NULL)) user system elapsed 0.66 0.01 0.67 bsms_n> ## parallel clusters bsms_n> cl <- makeCluster(nc) bsms_n> ## note: loading the package is optional bsms_n> system.time(clusterEvalQ(cl, library(bSims))) user system elapsed 0.00 0.00 0.88 bsms_n> system.time(bb <- b$replicate(B, cl=cl)) user system elapsed 0.01 0.00 0.50 bsms_n> stopCluster(cl) bsms_n> ## parallel forking bsms_n> if (.Platform$OS.type != "windows") { bsms_n+ system.time(bb <- b$replicate(B, cl=nc)) bsms_n+ } bsms_n> ## End(No test) bsms_n> bsms_n> bsms_n> ---------- bSims example: dist_fun2 ---------- dst_f2> tau <- c(1, 2, 3, 2, 1) dst_f2> d <- seq(0, 4, 0.01) dst_f2> dist_fun <- function(d, tau) exp(-(d/tau)^2) # half normal dst_f2> #dist_fun <- function(d, tau) exp(-d/tau) # exponential dst_f2> #dist_fun <- function(d, tau) 1-exp(-(d/tau)^-2) # hazard rate dst_f2> dst_f2> b <- c(0.5, 1, 1.5, 2) # boundaries dst_f2> op <- par(mfrow=c(2, 1)) dst_f2> plot(d, dist_fun2(d, tau[1], dist_fun), type="n", dst_f2+ ylab="P(detection)", xlab="Distance", axes=FALSE, dst_f2+ main="Sound travels from left to right") dst_f2> axis(1) dst_f2> axis(2) dst_f2> for (i in seq_len(length(b)+1)) { dst_f2+ x1 <- c(0, b, 4)[i] dst_f2+ x2 <- c(0, b, 4)[i+1] dst_f2+ polygon(c(0, b, 4)[c(i, i, i+1, i+1)], c(0, 1, 1, 0), dst_f2+ border=NA, dst_f2+ col=c("darkolivegreen1", "burlywood1", "lightgrey", dst_f2+ "burlywood1", "darkolivegreen1")[i]) dst_f2+ } dst_f2> lines(d, dist_fun2(d, tau[1], dist_fun)) dst_f2> lines(d, dist_fun2(d, tau[2], dist_fun)) dst_f2> lines(d, dist_fun2(d, tau[3], dist_fun)) dst_f2> lines(d, dist_fun2(d, tau, dist_fun, b), col=2, lwd=3) dst_f2> plot(rev(d), dist_fun2(d, tau[1], dist_fun), type="n", dst_f2+ ylab="P(detection)", xlab="Distance", axes=FALSE, dst_f2+ main="Sound travels from right to left") dst_f2> axis(1) dst_f2> axis(2) dst_f2> for (i in seq_len(length(b)+1)) { dst_f2+ x1 <- c(0, b, 4)[i] dst_f2+ x2 <- c(0, b, 4)[i+1] dst_f2+ polygon(c(0, b, 4)[c(i, i, i+1, i+1)], c(0, 1, 1, 0), dst_f2+ border=NA, dst_f2+ col=c("darkolivegreen1", "burlywood1", "lightgrey", dst_f2+ "burlywood1", "darkolivegreen1")[i]) dst_f2+ } dst_f2> lines(rev(d), dist_fun2(d, tau[1], dist_fun)) dst_f2> lines(rev(d), dist_fun2(d, tau[2], dist_fun)) dst_f2> lines(rev(d), dist_fun2(d, tau[3], dist_fun)) dst_f2> lines(rev(d), dist_fun2(d, tau, dist_fun, rev(4-b)), col=2, lwd=3) dst_f2> par(op) ---------- bSims example: events ---------- events> timetoevent(0, 10) numeric(0) events> timetoevent(Inf, 10) [1] 0 events> rr <- 1 events> tt <- timetoevent(rr, 10) events> op <- par(mfrow=c(1,2)) events> plot(ecdf(tt)) events> curve(1-exp(-rr*x), add=TRUE, col=2) # cdf events> plot(stepfun(sort(tt), 0:length(tt)/length(tt)), ylab="F(x)") events> curve(1-exp(-rr*x), add=TRUE, col=2) # cdf events> par(op) events> e <- events(movement=1, duration=60) events> mx <- max(abs(e[,1:2])) events> plot(e[,1:2], col="grey", type="l", asp=1, events+ xlim=2*c(-mx, mx), ylim=2*c(-mx, mx)) events> points(e[,1:2], col=e$v+1) events> abline(h=0, v=0, lty=2) events> legend("topright", pch=21, col=1:2, horiz=TRUE, events+ legend=c("movement", "vocalization")) ---------- bSims example: expand_list ---------- expnd_> b <- expand_list( expnd_+ movement = c(0, 1, 2), expnd_+ rint = list(c(0.5, 1, 1.5, Inf)), # in a list to keep as one expnd_+ xy_fun = list(NULL, function(z) z)) expnd_> b[[1]] $movement [1] 0 $rint [1] 0.5 1.0 1.5 Inf $xy_fun NULL expnd_> str(b) List of 6 $ :List of 3 ..$ movement: num 0 ..$ rint : num [1:4] 0.5 1 1.5 Inf ..$ xy_fun : NULL $ :List of 3 ..$ movement: num 1 ..$ rint : num [1:4] 0.5 1 1.5 Inf ..$ xy_fun : NULL $ :List of 3 ..$ movement: num 2 ..$ rint : num [1:4] 0.5 1 1.5 Inf ..$ xy_fun : NULL $ :List of 3 ..$ movement: num 0 ..$ rint : num [1:4] 0.5 1 1.5 Inf ..$ xy_fun :function (z) .. ..- attr(*, "srcref")= 'srcref' int [1:8] 11 23 11 35 23 35 11 11 .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' $ :List of 3 ..$ movement: num 1 ..$ rint : num [1:4] 0.5 1 1.5 Inf ..$ xy_fun :function (z) .. ..- attr(*, "srcref")= 'srcref' int [1:8] 11 23 11 35 23 35 11 11 .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' $ :List of 3 ..$ movement: num 2 ..$ rint : num [1:4] 0.5 1 1.5 Inf ..$ xy_fun :function (z) .. ..- attr(*, "srcref")= 'srcref' int [1:8] 11 23 11 35 23 35 11 11 .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' ---------- bSims example: get_nests ---------- gt_nst> phi <- 0.5 # singing rate gt_nst> tau <- 1:3 # EDR by strata gt_nst> dur <- 10 # simulation duration gt_nst> tbr <- c(3, 5, 10) # time intervals gt_nst> rbr <- c(0.5, 1, 1.5, Inf) # counting radii gt_nst> l <- bsims_init(10, 0.5, 1)# landscape gt_nst> p <- bsims_populate(l, 1) # population gt_nst> e <- bsims_animate(p, # events gt_nst+ vocal_rate=phi, duration=dur) gt_nst> d <- bsims_detect(e, # detections gt_nst+ tau=tau) gt_nst> x <- bsims_transcribe(d, # transcription gt_nst+ tint=tbr, rint=rbr) gt_nst> ## next locations gt_nst> head(get_nests(p)) i s x y g 1 1 H -2.477539 -2.1279875 NA 2 2 H -4.823014 -4.4605843 NA 3 3 H -3.696612 0.5544238 NA 4 4 H -2.161887 -2.9572063 NA 5 5 H -2.025115 0.1865945 NA 6 6 H -4.325523 2.3188061 NA gt_nst> head(get_nests(e)) i s x y g 1 1 H -2.477539 -2.1279875 G1 2 2 H -4.823014 -4.4605843 G1 3 3 H -3.696612 0.5544238 G1 4 4 H -2.161887 -2.9572063 G1 5 5 H -2.025115 0.1865945 G1 6 6 H -4.325523 2.3188061 G1 gt_nst> head(get_nests(d)) i s x y g 1 1 H -2.477539 -2.1279875 G1 2 2 H -4.823014 -4.4605843 G1 3 3 H -3.696612 0.5544238 G1 4 4 H -2.161887 -2.9572063 G1 5 5 H -2.025115 0.1865945 G1 6 6 H -4.325523 2.3188061 G1 gt_nst> head(get_nests(x)) i s x y g 1 1 H -2.477539 -2.1279875 G1 2 2 H -4.823014 -4.4605843 G1 3 3 H -3.696612 0.5544238 G1 4 4 H -2.161887 -2.9572063 G1 5 5 H -2.025115 0.1865945 G1 6 6 H -4.325523 2.3188061 G1 gt_nst> ## abundance gt_nst> get_abundance(p) [1] 97 gt_nst> get_abundance(e) [1] 97 gt_nst> get_abundance(d) [1] 97 gt_nst> get_abundance(x) [1] 97 gt_nst> ## density gt_nst> get_density(p) [1] 0.97 gt_nst> get_density(e) [1] 0.97 gt_nst> get_density(d) [1] 0.97 gt_nst> get_density(x) [1] 0.97 gt_nst> ## events gt_nst> head(get_events(e)) x y t v a i 1 4.1480370 1.9555012 0.02277008 1 196 75 2 -2.6825213 -0.9192479 0.04395041 1 201 24 3 -1.6709891 1.8654063 0.07025847 1 161 32 4 -3.8615339 -1.4014938 0.07230658 1 256 18 5 -2.0747972 -1.3641240 0.11111091 1 57 14 6 -0.4757805 1.7088403 0.15397749 1 171 50 gt_nst> head(get_events(d)) x y t v a f i 1 4.1480370 1.9555012 0.02277008 1 196 NA 75 2 -2.6825213 -0.9192479 0.04395041 1 201 NA 24 3 -1.6709891 1.8654063 0.07025847 1 161 NA 32 4 -3.8615339 -1.4014938 0.07230658 1 256 NA 18 5 -2.0747972 -1.3641240 0.11111091 1 57 NA 14 6 -0.4757805 1.7088403 0.15397749 1 171 NA 50 gt_nst> head(get_events(x)) x y t v a f i 1 4.1480370 1.9555012 0.02277008 1 196 NA 75 2 -2.6825213 -0.9192479 0.04395041 1 201 NA 24 3 -1.6709891 1.8654063 0.07025847 1 161 NA 32 4 -3.8615339 -1.4014938 0.07230658 1 256 NA 18 5 -2.0747972 -1.3641240 0.11111091 1 57 NA 14 6 -0.4757805 1.7088403 0.15397749 1 171 NA 50 gt_nst> ## detections gt_nst> head(get_detections(d)) x y t v a d f i j 5 -2.0747972 -1.3641240 0.1111109 1 57 2.4830662 NA 14 14 6 -0.4757805 1.7088403 0.1539775 1 171 1.7738383 NA 50 50 7 1.3013181 -1.7660829 0.1568836 1 14 2.1937360 NA 70 70 21 1.3013181 -1.7660829 0.6148041 1 95 2.1937360 NA 70 70 22 -0.2353686 0.7201898 0.6426921 1 358 0.7576752 NA 57 57 26 -0.7589418 -0.7406259 0.6900792 1 92 1.0604336 NA 40 40 gt_nst> head(get_detections(x)) x y t v a d f i j 5 -2.0747972 -1.3641240 0.1111109 1 57 2.4830662 NA 14 14 6 -0.4757805 1.7088403 0.1539775 1 171 1.7738383 NA 50 50 7 1.3013181 -1.7660829 0.1568836 1 14 2.1937360 NA 70 70 22 -0.2353686 0.7201898 0.6426921 1 358 0.7576752 NA 57 57 38 0.8103048 -2.7748346 0.9619932 1 269 2.8907267 NA 69 69 42 0.3815950 1.2397293 1.0016899 1 20 1.2971289 NA 51 51 gt_nst> get_table(x, "removal") 0-3min 3-5min 5-10min 0-50m 0 0 0 50-100m 1 0 0 100-150m 1 0 0 150+m 10 2 2 gt_nst> get_table(x, "visits") 0-3min 3-5min 5-10min 0-50m 0 1 0 50-100m 1 1 1 100-150m 1 3 4 150+m 10 7 11 ---------- bSims example: rlnorm2 ---------- rlnrm2> summary(rlnorm2(10^6, 1.3, 0.5)) # mean ~ 1.3 Min. 1st Qu. Median Mean 3rd Qu. Max. 0.09236 0.81925 1.14759 1.29986 1.60840 13.34982 rlnrm2> exp(log(1.3) - 0.5^2/2) # ~ median [1] 1.147246 ---------- bSims example: rmvn ---------- rmvn> rmvn(0, c(a=0, b=0), diag(1, 2, 2)) a b rmvn> rmvn(1, c(a=0, b=0), diag(1, 2, 2)) a b [1,] -1.216723 1.127792 rmvn> rmvn(2, c(a=0, b=0), diag(1, 2, 2)) a b [1,] -1.0016664 1.033442 [2,] 0.8402779 2.541317 rmvn> sapply(0:10, function(n) dim(rmvn(n, c(a=0, b=0), diag(1, 2, 2)))) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 0 1 2 3 4 5 6 7 8 9 10 [2,] 2 2 2 2 2 2 2 2 2 2 2 > > cat("\n\n---------- Parsing Shiny apps:", "----------\n\n") ---------- Parsing Shiny apps: ---------- > files <- list.files(system.file("shiny", package="bSims"), + pattern="\\.R$") > for (file in files) { + cat(" * file:", file) + tmp <- parse(file.path( + system.file("shiny", package="bSims"), file)) + cat(" -- OK\n") + } * file: bsimsH.R -- OK * file: bsimsHER.R -- OK * file: distfunH.R -- OK * file: distfunHER.R -- OK > > proc.time() user system elapsed 5.09 0.28 6.92