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-1 2024-05-20 pi-weer > > 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: 95 bsms_n> (a <- bsims_animate(p, vocal_rate=phi, duration=dur)) bSims events 1 km x 1 km stratification: HER total abundance: 95 duration: 10 min bsms_n> (o <- bsims_detect(a, tau=tau)) bSims detections 1 km x 1 km stratification: HER total abundance: 95 duration: 10 min detected: 28 heard bsms_n> (x <- bsims_transcribe(o, tint=tbr, rint=rbr)) bSims transcript 1 km x 1 km stratification: HER total abundance: 95 duration: 10 min detected: 28 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 1 0 0 50-100m 1 0 1 100-150m 1 0 1 150+m 8 3 1 bsms_n> get_table(x, "visits") 0-3min 3-5min 5-10min 0-50m 1 1 1 50-100m 1 1 2 100-150m 1 1 2 150+m 8 6 10 bsms_n> head(get_events(a)) x y t v a i 1 -4.216416 -1.665967 0.004256148 1 128 17 2 -1.819787 -1.397067 0.035912504 1 193 30 3 1.657029 2.333562 0.046456345 1 170 73 4 -2.014351 -4.179592 0.061841253 1 197 15 5 2.405320 2.377107 0.106170326 1 185 64 6 4.120898 2.317694 0.134605668 1 108 77 bsms_n> plot(get_events(a)) bsms_n> head(get_detections(o)) x y t v a d f i j 7 -0.1930718 -4.4713025 0.1919153 1 168 4.4754690 NA 47 47 8 1.0642734 0.3348885 0.2772081 1 0 1.1157187 NA 58 58 9 1.1458289 -1.3176072 0.3460845 1 301 1.7461422 NA 59 59 10 0.4426139 -0.3861752 0.3521154 1 43 0.5873997 NA 51 51 16 0.4426139 -0.3861752 0.4417134 1 20 0.5873997 NA 51 51 17 0.9592812 -3.2834556 0.4871696 1 78 3.4207164 NA 60 60 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: 59 duration: 10 min detected: 8 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 0 0 0 50-100m 1 0 0 100-150m 0 0 0 150+m 0 0 0 [[2]] 0-3min 3-5min 5-10min 0-50m 0 0 0 50-100m 1 0 0 100-150m 0 0 0 150+m 0 0 0 [[3]] 0-3min 3-5min 5-10min 0-50m 1 0 0 50-100m 1 0 0 100-150m 1 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.58 0.00 0.58 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 1.06 bsms_n> system.time(bb <- b$replicate(B, cl=cl)) user system elapsed 0.00 0.02 0.47 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 -4.879637 -3.18272168 NA 2 2 H -3.405989 4.22705952 NA 3 3 H -4.998252 -4.83519284 NA 4 4 H -3.487108 -0.09482567 NA 5 5 H -3.993120 3.29131962 NA 6 6 H -3.078560 0.74322227 NA gt_nst> head(get_nests(e)) i s x y g 1 1 H -4.879637 -3.18272168 G1 2 2 H -3.405989 4.22705952 G1 3 3 H -4.998252 -4.83519284 G1 4 4 H -3.487108 -0.09482567 G1 5 5 H -3.993120 3.29131962 G1 6 6 H -3.078560 0.74322227 G1 gt_nst> head(get_nests(d)) i s x y g 1 1 H -4.879637 -3.18272168 G1 2 2 H -3.405989 4.22705952 G1 3 3 H -4.998252 -4.83519284 G1 4 4 H -3.487108 -0.09482567 G1 5 5 H -3.993120 3.29131962 G1 6 6 H -3.078560 0.74322227 G1 gt_nst> head(get_nests(x)) i s x y g 1 1 H -4.879637 -3.18272168 G1 2 2 H -3.405989 4.22705952 G1 3 3 H -4.998252 -4.83519284 G1 4 4 H -3.487108 -0.09482567 G1 5 5 H -3.993120 3.29131962 G1 6 6 H -3.078560 0.74322227 G1 gt_nst> ## abundance gt_nst> get_abundance(p) [1] 102 gt_nst> get_abundance(e) [1] 102 gt_nst> get_abundance(d) [1] 102 gt_nst> get_abundance(x) [1] 102 gt_nst> ## density gt_nst> get_density(p) [1] 1.02 gt_nst> get_density(e) [1] 1.02 gt_nst> get_density(d) [1] 1.02 gt_nst> get_density(x) [1] 1.02 gt_nst> ## events gt_nst> head(get_events(e)) x y t v a i 1 -3.495883 -1.4240042 0.008108801 1 26 7 2 1.119322 2.3768926 0.023272952 1 151 68 3 -2.726107 -0.4036950 0.027281035 1 50 23 4 -3.564113 0.3150039 0.046934498 1 180 19 5 2.791033 2.0182472 0.050672594 1 217 92 6 3.035678 3.7355000 0.072153976 1 220 91 gt_nst> head(get_events(d)) x y t v a f i 1 -3.495883 -1.4240042 0.008108801 1 26 NA 7 2 1.119322 2.3768926 0.023272952 1 151 NA 68 3 -2.726107 -0.4036950 0.027281035 1 50 NA 23 4 -3.564113 0.3150039 0.046934498 1 180 NA 19 5 2.791033 2.0182472 0.050672594 1 217 NA 92 6 3.035678 3.7355000 0.072153976 1 220 NA 91 gt_nst> head(get_events(x)) x y t v a f i 1 -3.495883 -1.4240042 0.008108801 1 26 NA 7 2 1.119322 2.3768926 0.023272952 1 151 NA 68 3 -2.726107 -0.4036950 0.027281035 1 50 NA 23 4 -3.564113 0.3150039 0.046934498 1 180 NA 19 5 2.791033 2.0182472 0.050672594 1 217 NA 92 6 3.035678 3.7355000 0.072153976 1 220 NA 91 gt_nst> ## detections gt_nst> head(get_detections(d)) x y t v a d f i j 2 1.1193218 2.3768926 0.02327295 1 151 2.6272609 NA 68 68 11 -1.4891123 3.4400339 0.13386602 1 82 3.7485048 NA 45 45 15 2.2483441 -0.3216313 0.32955540 1 334 2.2712326 NA 78 78 16 0.7687933 -2.4843894 0.33329143 1 21 2.6006218 NA 70 70 32 -0.8324822 -0.5091484 0.58833475 1 160 0.9758374 NA 55 55 44 0.1923586 -3.1665576 0.76250671 1 103 3.1723948 NA 60 60 gt_nst> head(get_detections(x)) x y t v a d f i j 2 1.1193218 2.3768926 0.02327295 1 151 2.6272609 NA 68 68 11 -1.4891123 3.4400339 0.13386602 1 82 3.7485048 NA 45 45 15 2.2483441 -0.3216313 0.32955540 1 334 2.2712326 NA 78 78 16 0.7687933 -2.4843894 0.33329143 1 21 2.6006218 NA 70 70 32 -0.8324822 -0.5091484 0.58833475 1 160 0.9758374 NA 55 55 44 0.1923586 -3.1665576 0.76250671 1 103 3.1723948 NA 60 60 gt_nst> get_table(x, "removal") 0-3min 3-5min 5-10min 0-50m 0 0 0 50-100m 4 0 0 100-150m 2 0 0 150+m 10 2 0 gt_nst> get_table(x, "visits") 0-3min 3-5min 5-10min 0-50m 0 0 0 50-100m 4 3 4 100-150m 2 0 2 150+m 10 7 10 ---------- bSims example: rlnorm2 ---------- rlnrm2> summary(rlnorm2(10^6, 1.3, 0.5)) # mean ~ 1.3 Min. 1st Qu. Median Mean 3rd Qu. Max. 0.09709 0.81885 1.14811 1.30026 1.60863 12.23329 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.260816 0.7468329 rmvn> rmvn(2, c(a=0, b=0), diag(1, 2, 2)) a b [1,] 1.735134 0.44255782 [2,] -1.891226 -0.09748579 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 4.60 0.34 6.56