The live application below is running completely client-side in your browser using serverless Shinylive (WebAssembly). It provides an interactive Isle Royale Wolf-Moose Platform (IsleRoyaleApp()) featuring 100% offline spatial substrate mapping, continuous micro-step history tracking, side-by-side cowplot age class plots with per-species legends, custom time unit labeling (Steps vs Days), tab-aware sidebar decluttering (show_habitat, show_landmarks, norm, total), and dynamic table discovery (inputApp).
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 880
#| components: [viewer]
library(shiny)
library(bslib)
library(ggplot2)
library(cowplot)
library(dplyr)
library(tidyr)
library(tibble)
library(rlang)
library(stats)
library(graphics)
# --- Auto-Included Data Table: future.moose ---
future.moose <- structure(list(current = c("calf", "yearling", "adult", "senior",
"adult", "death", "predated"), future = c("yearling", "adult",
"senior", "death", "calf", "death", "death"), fid = c(2L, 3L,
4L, 5L, 6L, 7L, 7L), time = c(365L, 365L, 2920L, 2190L, 365L,
0L, 0L), pch = c("0", "1", "2", "3", "G", "D", "D"), color = c("lightgreen",
"green", "darkgreen", "brown", "black", "red", "red"), ageclass = c("calf",
"yearling", "adult", "senior", "adult", NA, NA), event = c("future",
"future", "future", "future", "birth", "death", "death"), init = c(50L,
50L, 400L, 100L, 50L, 0L, 0L)), class = "data.frame", row.names = c(NA,
-7L))
# --- Auto-Included Data Table: future.wolf ---
future.wolf <- structure(list(current = c("pup", "subadult", "adult", "adult",
"death"), future = c("subadult", "adult", "death", "pup", "death"
), fid = 2:6, time = c(365L, 365L, 2555L, 365L, 0L), pch = c("0",
"1", "2", "P", "D"), color = c("orange", "darkorange", "black",
"purple", "red"), ageclass = c("pup", "subadult", "adult", "adult",
NA), event = c("future", "future", "future", "birth", "death"
), init = c(5L, 10L, 15L, 5L, 0L)), class = "data.frame", row.names = c(NA,
-5L))
# --- Auto-Included Data Table: moose.wolf ---
moose.wolf <- structure(list(feed = c(8, 3, 1, 6, 0, 0), offspring = c(1.5,
1, 0.5, 1.2, 0, 0), male = c(NA, NA, NA, NA, NA, NA)), class = "data.frame", row.names = c("calf",
"yearling", "adult", "senior", "death", "predated"))
# --- Auto-Included Data Table: substrate.moose ---
substrate.moose <- structure(list(substrate = c(1.2, 1, 1, 0.8)), class = "data.frame", row.names = c("calf",
"yearling", "adult", "senior"))
# --- Auto-Included Data Table: substrate.substrate ---
substrate.substrate <- structure(list(substrate = 1), class = "data.frame", row.names = "substrate")
# --- Auto-Included Data Table: substrate.wolf ---
substrate.wolf <- structure(list(substrate = c(1, 1, 1)), class = "data.frame", row.names = c("pup",
"subadult", "adult"))
# --- Auto-Included Historical Census Data: wolf_moose ---
wolf_moose <- structure(list(Year = 1980:2019, Wolves = c(50L, 30L, 14L, 23L,
24L, 22L, 20L, 16L, 12L, 11L, 15L, 12L, 12L, 13L, 15L, 16L, 22L,
24L, 14L, 25L, 29L, 19L, 17L, 19L, 29L, 30L, 30L, 21L, 23L, 24L,
19L, 16L, 9L, 8L, 9L, 3L, 2L, 2L, 2L, 14L), Moose = c(664L, 650L,
700L, 900L, 811L, 1062L, 1025L, 1380L, 1653L, 1397L, 1216L, 1313L,
1600L, 1880L, 1800L, 2400L, 1200L, 500L, 700L, 750L, 850L, 900L,
1000L, 900L, 750L, 540L, 385L, 450L, 650L, 530L, 510L, 515L,
750L, 975L, 1050L, 1250L, 1300L, 1600L, 1500L, 2060L)), class = "data.frame", row.names = c(NA,
-40L))
# --- Source: spline.R ---
## $Id: spline.R,v 1.0 2002/12/09 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
###########################################################################################
## Support routines for future.meanvalue
###########################################################################################
spline.rate <- function( meanvalue, x )
{
coeff <- stats::coef( meanvalue )
for( i in 2:ncol( coeff ))
coeff[,i-1] <- i * coeff[,i]
coeff[,ncol( coeff ) ] <- 0
meanvalue$coefficients <- coeff
if( missing( x ))
stats::predict( meanvalue )
else
stats::predict( meanvalue, x )
}
###########################################################################################
spline.deriv <- function( s )
{
s$coefficients <- s$coefficients[,-1]
for( i in seq( 2, ncol( s$coefficients )))
s$coefficients[,i] <- s$coefficients[,i] * i
s
}
###########################################################################################
splinesum <- function( xy, fit = splines::interpSpline( xy$x, xy$y ),
log2 = log( 2 ), tol = 1e-5 )
{
n <- nrow( xy )
## mean time to future event (assume linear off end )
maxx <- xy[n,"x"]
mean.y <- maxx * mean( exp( - stats::predict( fit )$y )) + xy[1,"x"]
rate <- spline.rate( fit, maxx )$y
if( rate > tol )
mean.y <- mean.y + exp( - xy[n,"y"] ) / rate
## median time to future event
adiff <- diff( stats::coef( fit )[,1] )
## undefined if mean value function not monotone
if( !( all( adiff < 0) || all( adiff > 0 )))
median.y <- NA
else {
invmvalue <- splines::backSpline( fit )
median.y <- stats::predict( invmvalue, log2 )$y
if( is.na( median.y ))
median.y <- spline.extrapolate( fit, invmvalue, log2 )
}
c( mean = mean.y, median = median.y )
}
###########################################################################################
summaryshow <- function( xy, fit, col = "black", sums = splinesum( xy, fit ))
{
tmpar <- graphics::par( col = col )
graphics::mtext( paste( "mean =", round( sums[1], 2 )), 3, at = graphics::par("usr")[1], adj = 0 )
medianshow <- if( is.na( sums[2] ))
"curve not monotone"
else
paste( "median =", round( sums[2], 2 ))
graphics::mtext( medianshow, 3, at = graphics::par("usr")[2] / 1.25, adj = 1 )
graphics::par( tmpar )
invisible( sums )
}
###########################################################################################
curve.plot <- function( xy = seq(0,1,by=.25), y = c(0,.1,.5,.9,1),
z = graphics::locator(1,"n"), n=5, action="add",
fit = splines::interpSpline( xy$x, xy$y ),
backfit = TRUE, save.ends = 3, col = c("blue","red"), lwd = 4,
f = function( x ) x, finv = function( x ) x )
{
if( !is.list( xy )) {
xy <- data.frame( x = xy )
xy$y <- y
}
else
xy <- as.data.frame( xy )
if( !match( action, c("refresh","finish"), nomatch = 0 )) {
z <- as.data.frame(z)
tmp <- z$x > max( xy$x )
if( any( tmp )) {
if( all( tmp ))
return( xy )
z$x <- z$x[!tmp]
z$y <- z$y[!tmp]
}
}
remove.points <- function( xy, finv, save.ends = TRUE ) {
## find closest point after standardizing distances
usr <- graphics::par( "usr" )
tmp <- (( xy$x - z$x ) / diff( usr[1:2] )) ^ 2 +
(( xy$y - z$y ) / diff( usr[3:4] )) ^ 2
if( save.ends )
tmp <- tmp[ - c( 1, nrow( xy )) ]
tmp <- save.ends + min( seq( tmp )[ tmp == min( tmp ) ] )
tmpd <- xy[tmp,"y"]
tmpd <- finv( tmpd )
graphics::points( xy[tmp,"x"], tmpd, lwd = lwd, col = "white" )
tmp
}
for( i in 1:n) {
switch( action,
noaction =
return( xy )
,
add = {
graphics::points(z$x,z$y, lwd = lwd )
z$y <- f( z$y )
xy <- rbind(xy,z)
xy <- xy[order(xy$x),]
fit <- splines::interpSpline( xy$x, xy$y )
},
replace = {
graphics::points(z$x,z$y, lwd = lwd )
z$y <- f( z$y )
tmp <- remove.points( xy, finv, save.ends == 3 )
tmpp <- ( tmp > 1 & tmp < length( xy$x ))
if( save.ends != 2 | tmpp )
xy$x[tmp] <- z$x
if( save.ends != 1 | tmpp )
xy$y[tmp] <- z$y
fit <- splines::interpSpline( xy$x, xy$y )
},
delete = {
z$y <- f( z$y )
tmp <- remove.points( xy, finv, save.ends > 0 )
xy <- xy[-tmp,]
fit <- splines::interpSpline( xy$x, xy$y )
},
refresh =, finish = {
xy
}
)
tmpp <- stats::predict( fit )
graphics::lines( tmpp$x, finv( tmpp$y ), col = col[1], lwd = lwd )
adiff <- diff( stats::coef( fit )[,1] )
if( backfit & ( all( adiff < 0) || all( adiff > 0 ))) {
## backspline (not quite a spline) fit to inverse
tmpback <- stats::predict( splines::backSpline( fit ))
graphics::lines( tmpback$y, finv( tmpback$x ), col = col[2], lwd = lwd )
}
}
list( xy = xy[order(xy$x),], fit = fit )
}
###########################################################################################
cdf.lines <- function( data, fig = "mean value", nspline = 8,
conf = c(50,80,90,95), rescale = 1,
col = c("green","blue","red","orange") )
{
rate <- fig == "mean value"
n <- length( data )
data <- sort( data )
prob <- seq( n ) / ( n + 1 )
f <- - rescale * log( 1 - prob )
ylab <- "prob"
if( rate )
ylab <- "cum rate"
else
f <- 1 - exp( -f )
graphics::lines( data, f, lwd = 2 )
conf <- conf / 100
for( i in seq( length( conf ))) {
## lower confidence
tmp <- stats::qbinom( conf[i], n, prob, lower.tail = TRUE ) / ( n + 1 )
tmp <- 1 - exp( rescale * log( 1 - tmp ))
tmpna <- is.na( tmp )
if( any( tmpna ))
tmp[tmpna] <- f[1]
if( rate )
tmp <- - log( 1 - tmp )
graphics::lines( data, tmp, lty = 2, col = col[i] )
## upper confidence
tmp <- stats::qbinom( conf[i], n, prob, lower.tail = FALSE ) / ( n + 1 )
tmp <- 1 - exp( rescale * log( 1 - tmp ))
tmpna <- is.na( tmp )
if( any( tmpna ))
tmp[tmpna] <- f[n]
if( rate )
tmp <- - log( 1 - tmp )
graphics::lines( data, tmp, lty = 2, col = col[i] )
}
}
###########################################################################################
rspline <- function( meantime = 1,
fivepar = c(dispersion = 1, location = 0, intensity = 1, truncation = 0,
rejection = Inf ),
fit = NULL,
meanvalue = fit$meanvalue, invmvalue = fit$invmvalue,
span = Inf )
{
## dispersion = a, location = b, intensity = c
## truncation = -log(1-d), rejection = -log(1-e)
## y = a M^-1( G(d)+cV ) + b if 1 - exp( -cV ) < G(e)
## y = span if 1 - exp( -cV ) >= G(e)
## this is not quite right for truncation, as we know event happened before b
## if 1-exp(-cV) < d, but I am not sure how to pass that information along yet
default <- is.null( meanvalue )
## V ~ exp(1)
## intensity: V/c
## truncation: (G(d)+V)/c
rate <- ( fivepar["truncation"] + rexp( 1 )) / fivepar["intensity"]
## mean value inverse: M^-1( G(d)+V/c )
if( default )
y <- rate
else {
y <- stats::predict( invmvalue, rate )$y
## kludge to linearly extrapolate beyond cubic spline fit
if( is.na( y ))
y <- spline.extrapolate( meanvalue, invmvalue, rate )
}
## dispersion and location: y = a M^-1( cV ) + b
y <- meantime * fivepar["dispersion"] * y + fivepar["location"]
## rejection: y >= e? then set y to span
if( y > fivepar["rejection"] )
y <- span
y
}
###########################################################################################
spline.extrapolate <- function( meanvalue, invmvalue, x )
{
## linear extrapolation of inverse spline beyond upper end
coeff <- stats::coef( invmvalue )
nr <- nrow( coeff ) - 1
xknot <- splines::splineKnots( invmvalue )[nr+(0:1)]
yknot <- splines::splineKnots( meanvalue )[nr+1]
tmpr <- ( xknot[2] - xknot[1] )
slope <- coeff[nr,2] + tmpr * ( 2 * coeff[nr,3] + tmpr * 3 * coeff[nr,4] )
yknot + slope * ( x - xknot[2] )
}
###########################################################################################
### spline.design() is a prototype for designing spline curves
### Ultimately, pieces of spline.design, spline.temp() and spline.meanvalue()
### will be pulled out as subroutines to reduce code overlap
###########################################################################################
spline.design <- function (y = yinit, x = xinit, nspline = 8, xy = data.frame(x = x,
y = y), n = 1, horizontal = FALSE)
{
is.data <- !missing(y)
if (is.data) {
data <- y
if (missing(x))
x <- as.numeric(names(y))
datax <- x
ndata <- length(data)
choose <- round(seq(1, ndata, length = nspline))
xinit <- x <- x[choose]
yinit <- y <- y[choose]
}
else {
tmp <- seq(0, nspline - 1)
if (missing(x))
xinit <- tmp
else xinit <- x
if (missing(y))
yinit <- rep(50, nspline)
else yinit <- y
}
## plot curve and surrounding axes
graphics::par( mfrow = c(1,1), mar = rep(4.1,4))
plotit <- function( xy, fig = "temp", fit = splines::interpSpline( xy$x, xy$y ),
horizontal = FALSE, strip = .25, margin = .1 )
{
switch( fig, {
y <- xy$y
ylim <- range(y)
}
)
xlim <- range(xy$x)
xlim <- xlim + c(-1,1) * margin * diff( xlim )
if( horizontal ) {
if( diff( ylim ) == 0 )
ylim <- ylim * c(.75,1.25)
separator <- ylim[2]
ylim[2] <- ylim[2] + strip * diff( ylim )
}
else {
separator <- xlim[2]
xlim[2] <- xlim[2] + strip * diff( xlim )
}
axt <- c("n","s")
tmpar <- graphics::par( xaxt = axt[1+horizontal], yaxt = axt[2-horizontal] )
plot( xy$x, y, xlim = xlim, ylim = ylim, type="n", xlab = "", ylab = "" )
graphics::par( xaxt = "s", yaxt = "s" )
graphics::points( xy$x, y, lwd = 4 )
graphics::title( fig )
graphics::mtext( "time", 1, 2 )
graphics::mtext( fig, 2, 2 )
if( horizontal ) {
p <- pretty( c(ylim[1],separator) )
graphics::axis( 2, p[ p <= separator ] )
graphics::abline( h = separator, lty = 2 )
}
else {
p <- pretty( c(xlim[1],separator) )
graphics::axis( 1, p[ p <= separator ] )
graphics::abline( v = separator, lty = 2 )
}
curve.plot( xy, n = n, action = "refresh", fit = fit, backfit = FALSE,
save.ends = 0 )
separator
}
## place commands along right strip of plot, highlighting current command
plotcmd <- function( ans, fig, cmds, cmdlocs, usr, col = "green", rest = "black",
horizontal = TRUE )
{
ans <- c( ans, fig )
tmp <- is.na( match( cmds, ans ))
if( any( tmp )) for( i in unique( cmdlocs$adj )) {
tmpi <- tmp & i == cmdlocs$adj
if( any( tmpi ))
graphics::text( cmdlocs$x[tmpi], cmdlocs$y[tmpi], cmds[tmpi], col = rest, adj = i )
}
if( any( !tmp )) for( i in unique( cmdlocs$adj )) {
tmpi <- !tmp & i == cmdlocs$adj
if( any( tmpi ))
graphics::text( cmdlocs$x[tmpi], cmdlocs$y[tmpi], cmds[tmpi], col = col, adj = i )
}
}
cmds <- c("add","delete","replace","","finish","restart","refresh","rescale",
"","data","temp")
newlocs <- if( horizontal )
function( cmds, data = FALSE, usr )
{
if( !data )
cmds <- cmds[ cmds != "data" ]
n <- length( cmds )
blank <- seq( n )[ cmds == "" | cmds == " " ]
tmp <- diff(usr[3:4]) / 20
m <- mean( usr[1:2] )
y <- usr[4] + 0.5 * tmp - c( tmp * seq( blank[1] - 1 ), 0,
tmp * seq( blank[2] - blank[1] - 1 ), 0,
tmp * seq( n - blank[2] ))
x <- c( rep( usr[1], blank[1] - 1 ), mean( m, usr[1] ),
rep( m, blank[2] - blank[1] - 1 ), mean( m, usr[2] ),
rep( usr[2], n - blank[2] ))
adj <- c( rep( 0, blank[1] ),
rep( 0.5, blank[2] - blank[1] ),
rep( 1, n - blank[2] ))
tmp <- data.frame( x = x, y = y, adj = adj )
cmds[blank[2]] <- " "
row.names( tmp ) <- cmds
tmp
}
else
function( cmds, data = FALSE, usr )
{
if( !data )
cmds <- cmds[ cmds != "data" ]
n <- length( cmds )
blank <- seq( n )[ cmds == "" ]
tmp <- diff(usr[3:4]) / 20
m <- mean( usr[3:4] )
tmp <- c( usr[4] - tmp * seq( blank[1] - 1 ),
mean( m, usr[4] ),
m + tmp * ( seq( blank[1] + 1, blank[2] - 1 ) - mean( blank )),
mean( m, usr[3] ),
usr[3] + tmp * seq( n - blank[2] ))
tmp <- data.frame( x = rep( usr[2], n ), y = tmp, adj = rep( 1, n ))
cmds[blank[2]] <- " "
row.names( tmp ) <- cmds
tmp
}
fig <- "temp"
newans <- ans <- "replace"
graphics::par( mar = c(4.1,4.1,3.1,4.1),omi=rep(.25,4))
fit <- splines::interpSpline( xy$x, xy$y )
separator <- plotit( xy, fig, fit, horizontal )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
cmds <- row.names( cmdlocs )
plotcmd( ans, fig, cmds, cmdlocs, usr )
use.data <- FALSE
rescale.data <- 1
repeat {
## get command from plot using cursor
z <- graphics::locator(1,"n")
if(( !horizontal & z$x > separator ) | ( horizontal * z$y > separator )) {
if( horizontal ) { # need to look at both z&y
x <- abs( z$x - cmdlocs$x )
x <- x == min( x )
newans <- cmds[x]
z <- abs( z$y - cmdlocs$y )[x]
newans <- newans[ z == min( z ) ][1]
}
else {
z <- abs(z$y - cmdlocs$y )
newans <- cmds[ z == min( z ) ][1]
}
switch( newans,
finish =, refresh = {
separator <- plotit( xy, fig, fit, horizontal )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
data = {
use.data <- is.data & !use.data
if( is.data & !use.data )
plotit( xy, fig, fit, horizontal )
},
temp = {
fig <- newans
separator <- plotit( xy, fig, fit, horizontal )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
rescale = {
cat( "enter new values followed by RETURN key\n" )
tmpy <- readline( paste( "maximum ", fig, "(",
round( max( xy$y ), 2 ), "):", sep = "" ))
if( tmpy != "" ) {
tmpy <- suppressWarnings(as.numeric( tmpy )) / max( xy$y )
xy$y <- tmpy * xy$y
if( is.data )
rescale.data <- rescale.data * tmpy
}
tmpx <- readline( paste( "maximum time(",
round( max( xy$x ), 2 ), "):", sep = "" ))
if( tmpx != "" ) {
tmpx <- suppressWarnings(as.numeric( tmpx )) / max( xy$x )
xy$x <- tmpx * xy$x
}
fit <- splines::interpSpline( xy$x, xy$y )
separator <- plotit( xy, fig, fit, horizontal )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
restart = {
if( is.data )
rescale.data <- 1
xy <- data.frame( x = xinit, y = yinit )
fit <- splines::interpSpline( xy$x, xy$y )
separator <- plotit( xy, fig, fit, horizontal )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
add =, delete =, replace = {
ans <- newans
}
)
if( use.data ) {
rx <- range( xy$x )
dx <- range( datax )
graphics::lines( rx[1] + ( datax - dx[1] ) * diff( rx ) / diff( dx ), data * rescale.data )
}
plotcmd( ans, fig, cmds, cmdlocs, usr )
}
else {
fit <- curve.plot( xy, n = n, action = ans, z = z, fit = fit, backfit = FALSE,
save.ends = 0 )
xy <- fit$xy
fit <- fit$fit
ans <- "replace"
plotcmd( ans, fig, cmds, cmdlocs, usr )
}
if( newans == "finish" )
break
}
plotcmd( newans, fig, cmds, cmdlocs, "red" )
tmp <- curve.plot( xy, n = n, action = "refresh", backfit = FALSE, save.ends = 0 )
tmp
}
# --- Source: triangle.R ---
## $Id: triangle.R,v 0.9 2002/12/09 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
###########################################################################################
## rtri( n, width )
##
## plot.current( x, species )
## text.current( x, species )
##
###########################################################################################
###########################################################################################
### Tridiagonal Coordinate System S3 Classes & Algebra
###########################################################################################
tricoord <- function(a, b = NULL, c = NULL) {
if (is.data.frame(a) && all(c("a", "b", "c") %in% names(a))) {
res <- a
} else if (is.matrix(a) && ncol(a) == 3) {
res <- as.data.frame(a)
names(res) <- c("a", "b", "c")
} else if (!is.null(b) && !is.null(c)) {
res <- data.frame(a = a, b = b, c = c)
} else if (is.numeric(a) && length(a) == 3) {
res <- data.frame(a = a[1], b = a[2], c = a[3])
} else {
stop("Invalid tricoord input format")
}
class(res) <- c("tricoord", "data.frame")
res
}
`+.tricoord` <- function(e1, e2) {
# Handle vector offsets cleanly
if (is.numeric(e2) && length(e2) == 3) {
tricoord(e1$a + e2[1], e1$b + e2[2], e1$c + e2[3])
} else if (inherits(e2, "tricoord")) {
tricoord(e1$a + e2$a, e1$b + e2$b, e1$c + e2$c)
} else {
stop("Invalid right hand operand for tricoord addition")
}
}
`-.tricoord` <- function(e1, e2) {
if (is.numeric(e2) && length(e2) == 3) {
tricoord(e1$a - e2[1], e1$b - e2[2], e1$c - e2[3])
} else if (inherits(e2, "tricoord")) {
tricoord(e1$a - e2$a, e1$b - e2$b, e1$c - e2$c)
} else {
stop("Invalid right hand operand for tricoord subtraction")
}
}
###########################################################################################
rtri <- function( n, width, tri = matrix(0,3,n), roundoff = TRUE )
{
tri <- as.matrix( tri )
if( n == 1 ) {
xy <- stats::runif( 2, 0, width )
if( roundoff )
xy <- round( xy )
i <- sample( 3, 1 )
i1 <- 1 + i%%3
i2 <- 1 + (i+1)%%3
tri[i1,] <- tri[i1,] + xy[1]
tri[i2,] <- tri[i2,] - xy[2]
tri[i,] <- - ( tri[i1,] + tri[i2,] )
return( tri )
}
else {
xy <- data.frame( x = stats::runif( n, 0, width ),
y = - stats::runif( n, 0, width ))
if( roundoff )
xy <- round( xy )
out <- sample( 3, n, replace = TRUE )
for( i in 1:3 ) {
outi <- out == i
if( any( outi )) {
i1 <- 1 + i%%3
i2 <- 1 + (i+1)%%3
tri[i1,outi] <- tri[i1,outi] + xy$x[outi]
tri[i2,outi] <- tri[i2,outi] + xy$y[outi]
tri[i,outi] <- - ( tri[i1,outi] + tri[i2,outi] )
}
}
}
tri
}
###########################################################################################
car2tri.default <- function(x,y)
car2tri( cbind( x, y ))
car2tri <- function( xy, xmult = ( 2 + sq3 ) / 4, ymult = ( 3 + 2 * sq3 ) / 12,
sq3 = sqrt( 3 ))
{
# if( !is.matrix( xy ))
# xy <- t( as.matrix( xy ))
aa <- xmult * xy[,1] - ymult * xy[,2]
bb <- - xmult * xy[,1] - ymult * xy[,2]
cc <- -( aa + bb )
rbind( a = aa, b = bb, c = cc )
}
###########################################################################################
tri2car.default <- function(aa,bb,cc=-(aa+bb))
tri2car( rbind( aa, bb, cc ))
tri2car <- function(tri, xmult = 2 / ( 2 + sq3 ), ymult = 6 / ( 3 + 2 * sq3 ),
sq3 = sqrt( 3 ))
{
if( inherits(tri, "tricoord") ) {
# If the user passes our S3 tricoord dataframe, map it correctly natively.
x <- ( tri$a - tri$b ) * xmult
y <- - ( tri$a + tri$b ) * ymult
} else {
if( !is.matrix( tri ))
tri <- as.matrix( tri )
x <- ( tri[1,] - tri[2,] ) * xmult
y <- - ( tri[1,] + tri[2,] ) * ymult
}
data.frame( x = x, y = y )
}
###########################################################################################
cardist <- function( xy )
sqrt( xy[,1]^2 + xy[,2]^2 )
###########################################################################################
tridist <- function( tri )
apply( tri, 1, max )
###########################################################################################
gasket <- function( aa, bb )
{
n <- length( aa )
pp <- c(-1,1,0,1)
dda <- diff( aa )
ddb <- diff( bb )
ss <- sign( sign( dda ) - sign( ddb ))
dda <- 2 - abs( dda )
ddb <- 2 - abs( ddb )
aa <- 2 * aa
aa <- c( aa[1], rbind( aa[-n] + pp[dda+1+ss], aa[-1] + pp[dda+1-ss], aa[-1] ))
bb <- 2 * bb
bb <- c( bb[1], rbind( bb[-n] + pp[ddb+1-ss], bb[-1] + pp[ddb+1+ss], bb[-1] ))
data.frame( aa = aa, bb = bb )
}
# --- Source: substrate_triangle.R ---
get_substrate_grid <- function(width, step = 1, orientation = "up") {
pts <- expand.grid(a = seq(0, width - step, by = step),
b = seq(0, width - step, by = step))
if (orientation == "up") {
pts <- subset(pts, a + b <= width - step)
pts$c <- -(pts$a + pts$b)
} else {
pts <- subset(pts, a + b <= width - step)
pts$a <- -pts$a
pts$b <- -pts$b
pts$c <- -(pts$a + pts$b)
}
return(pts)
}
substrate_topology <- function(width = 10, step = 1) {
W <- width - step
# Topology adjacency offsets
list(
fr2 = list(offset = c(0, 0, 0), dir = "down"),
fr1 = list(offset = c(-W, -W, 2*W), dir = "up"),
fr3 = list(offset = c(-W, 0, W), dir = "up"),
fr4 = list(offset = c(0, -W, W), dir = "up"),
tw1 = list(offset = c(-W, 0, W), dir = "down"),
twig = list(offset = c(-W, 0, W), dir = "down"),
tw2 = list(offset = c(-2*W, 0, 2*W), dir = "up"),
lftop = list(offset = c(-2*W, W, W), dir = "down"),
lfbot = list(offset = c(-3*W, W, 2*W), dir = "up")
)
}
create_substrate <- function(topology, width = 10, step = 1) {
W <- width - step
all_points <- data.frame()
labels_df <- data.frame()
poly_df <- data.frame()
for (sub in names(topology)) {
cfg <- topology[[sub]]
grid <- get_substrate_grid(width, step, cfg$dir)
o_a <- cfg$offset[1]
o_b <- cfg$offset[2]
o_c <- cfg$offset[3]
# Needs tricoord and tri2car which are presumably exported/available from R/triangle.R
grid_tri <- tricoord(grid$a, grid$b, grid$c)
grid_tri <- grid_tri + cfg$offset
car_pts <- tri2car(grid_tri)
car_pts$substrate <- sub
all_points <- rbind(all_points, car_pts)
# Determine bounds and midpoints for side labels 1,2,3
if (cfg$dir == "up") {
v_top <- c(o_a, o_b, o_c)
v_br <- c(o_a + W, o_b, o_c - W)
v_bl <- c(o_a, o_b + W, o_c - W)
m1 <- (v_top + v_bl) / 2
m2 <- (v_top + v_br) / 2
m3 <- (v_bl + v_br) / 2
centroid <- (v_top + v_br + v_bl) / 3
p_mat <- cbind(v_top, v_br, v_bl)
} else {
v_bot <- c(o_a, o_b, o_c)
v_tr <- c(o_a, o_b - W, o_c + W)
v_tl <- c(o_a - W, o_b, o_c + W)
m1 <- (v_bot + v_tr) / 2
m2 <- (v_bot + v_tl) / 2
m3 <- (v_tl + v_tr) / 2
centroid <- (v_bot + v_tr + v_tl) / 3
p_mat <- cbind(v_bot, v_tr, v_tl)
}
# Interpolate slightly towards the centroid to put text "just inside" the edges
w_in <- 0.25 # weight towards centroid
l1 <- m1 * (1 - w_in) + centroid * w_in
l2 <- m2 * (1 - w_in) + centroid * w_in
l3 <- m3 * (1 - w_in) + centroid * w_in
mat_l <- cbind(l1, l2, l3)
car_l <- tri2car(mat_l)
car_l$label <- c("1", "2", "3")
car_l$substrate <- sub
labels_df <- rbind(labels_df, car_l)
car_p <- tri2car(p_mat)
car_p$substrate <- sub
poly_df <- rbind(poly_df, car_p)
}
centers <- stats::aggregate(cbind(x,y) ~ substrate, data=all_points, mean)
obj <- list(
points = all_points,
labels = labels_df,
poly = poly_df,
centers = centers,
topology = topology
)
class(obj) <- "substrate"
return(obj)
}
autoplot.substrate <- function(object, ...) {
ggplot2::ggplot() +
# Draw the black boundary lines outlining the substrates exactly over outer dots
ggplot2::geom_polygon(data=object$poly, ggplot2::aes(x=x, y=y, group=substrate), fill=NA, color="black", linewidth=0.7) +
# Plot grid dots
ggplot2::geom_point(data=object$points, ggplot2::aes(x=x, y=y, color=substrate), size=1.5) +
# Plot Substrate Labels (Centers)
ggplot2::geom_text(data=object$centers, ggplot2::aes(x=x, y=y, label=substrate), color="black", fontface="bold", size=5) +
# Axis side numbers
ggplot2::geom_text(data=object$labels, ggplot2::aes(x=x, y=y, label=label), color="darkred", fontface="bold", size=4) +
ggplot2::theme_void() +
ggplot2::coord_fixed() +
ggplot2::ggtitle("Ewing Tridiagonal Substrate Network Mapping")
}
create_hex_overlay <- function(object, step = 1) {
pts <- if (inherits(object, "substrate")) object$points else object
if (is.null(pts) || nrow(pts) == 0) return(data.frame())
xmult <- 2 / (2 + sqrt(3))
ymult <- 6 / (3 + 2 * sqrt(3))
d <- step * sqrt(xmult^2 + ymult^2)
r <- d / sqrt(3)
angles <- (seq(0, 5) * 60 + 30) * pi / 180
dx <- r * cos(angles)
dy <- r * sin(angles)
n_pts <- nrow(pts)
hex_list <- vector("list", n_pts)
for (i in seq_len(n_pts)) {
px <- pts$x[i] + dx
py <- pts$y[i] + dy
sub <- pts$substrate[i]
hex_list[[i]] <- data.frame(
x = px,
y = py,
cell_id = i,
substrate = sub,
stringsAsFactors = FALSE
)
}
do.call(rbind, hex_list)
}
# --- Source: ewing_substrate.R ---
ewing_substrate <- function( community,
species,
headstuff = c( 0, "start", sum( to.plot )),
units = getOrgFeature( community, species[1], "units" ),
right = species[1], adj = c(0,.5,1),
show_sub = NULL,
step = 0,
layout = c("facet", "hex"),
width = 10,
step_density = 1,
rescale = TRUE,
x_var = c("step", "time"),
...)
{
x_var <- match.arg(x_var)
if (inherits(community, "isle_royale_sim")) {
p_map <- autoplot(community$habitat_overlay, show_landmarks = TRUE)
moose_sf <- sf::st_as_sf(community$moose_pop, coords = c("lon", "lat"), crs = sf::st_crs(community$habitat_overlay$layer))
wolf_sf <- sf::st_as_sf(community$wolf_pop, coords = c("lon", "lat"), crs = sf::st_crs(community$habitat_overlay$layer))
hdr_str <- if (x_var == "step") paste0("Step ", community$nstep) else paste0("Time Units ", community$nstep)
p_map <- p_map +
ggplot2::geom_sf(data = moose_sf, color = "#27ae60", shape = 21, fill = NA, stroke = 1.0, size = 0.8, alpha = 0.85) +
ggplot2::geom_sf(data = wolf_sf, color = "#e74c3c", shape = 21, fill = NA, stroke = 1.4, size = 1.5, alpha = 0.95) +
ggplot2::ggtitle(paste0("Isle Royale Substrate Plot (", hdr_str, ")"))
return(p_map)
}
layout <- match.arg(layout)
if (length(species) > 1) {
res_list <- lapply(species, function(sp) {
df <- ewing_substrate(community = community, species = sp, headstuff = headstuff,
units = units, right = right, adj = adj, show_sub = show_sub,
step = step, layout = layout, width = width, step_density = step_density,
rescale = rescale, ...)
if (!is.null(df) && nrow(df) > 0) {
df$species <- sp
}
df
})
res_list <- res_list[!sapply(res_list, is.null)]
if (length(res_list) == 0) return(NULL)
combined <- do.call(rbind, res_list)
attr(combined, "species") <- paste(species, collapse = " & ")
step_val <- if (!is.null(community$step)) community$step else if (!is.null(community$count$step)) community$count$step else if (!is.null(attr(community, "nstep"))) attr(community, "nstep") else step
attr(combined, "step") <- step_val
attr(combined, "layout") <- layout
attr(combined, "width") <- width
attr(combined, "step_density") <- step_density
class(combined) <- c("ewing_substrate", class(combined))
return(combined)
}
## plot current stages for species (except random parasites)
organism <- get.species( community, species )[,-1]
if(is.null(organism)) # species is not in community
return(NULL)
future <- getOrgFuture( community, species, c("color","pch") )
# Substrate names mapping (e.g. fr1, fr2, fr3, fr4, twig, lftop, lfbot)
substrate_feat <- getOrgFeature( community, species, "substrate")
sub_interact <- getOrgInteract(community, substrate_feat, species)
substrates <- rownames(sub_interact)
if (is.null(substrates) || length(substrates) == 0) {
substrates <- names(getOrgInteract(community, substrate_feat, substrate_feat))
}
if (is.null(show_sub)) show_sub <- substrates
position <- paste( "pos", letters[1:3], sep = "." )
if (layout == "hex") {
topo <- substrate_topology(width = width, step = step_density)
n_org <- ncol(organism)
gx <- numeric(n_org)
gy <- numeric(n_org)
sub_indices <- organism["sub.stage", ]
org_sub_names <- substrates[sub_indices]
# Process organisms per substrate patch to rescale local coordinates into substrate surface triangle
unique_subs <- unique(org_sub_names)
for (sub in unique_subs) {
idx <- which(org_sub_names == sub)
target_name <- sub
if (!target_name %in% names(topo)) {
if (target_name == "twig") target_name <- "tw1"
if (target_name == "tw1") target_name <- "twig"
}
cfg <- topo[[target_name]]
pa <- organism["pos.a", idx]
pb <- organism["pos.b", idx]
pc <- organism["pos.c", idx]
if (!is.null(cfg)) {
# Determine substrate surface width (allows substrates of different sizes in future)
W_sub <- if (!is.null(cfg$width)) (cfg$width - step_density) else (width - step_density)
if (rescale) {
amin <- min(pa); amax <- max(pa)
bmin <- min(pb); bmax <- max(pb)
u <- if (amax > amin) (pa - amin) / (amax - amin) else rep(0.5, length(idx))
v <- if (bmax > bmin) (pb - bmin) / (bmax - bmin) else rep(0.5, length(idx))
# 15% inner padding to ensure symbols sit comfortably inside substrate polygon borders
u_m <- 0.15 + 0.70 * u
v_m <- 0.15 + 0.70 * v
a_p <- u_m * W_sub
b_p <- v_m * (W_sub - a_p)
c_p <- -(a_p + b_p)
} else {
a_p <- pa
b_p <- pb
c_p <- pc
}
off <- cfg$offset
if (cfg$dir == "up") {
ga <- a_p + off[1]
gb <- b_p + off[2]
gc <- c_p + off[3]
} else {
ga <- -a_p + off[1]
gb <- -b_p + off[2]
gc <- -c_p + off[3]
}
car <- tri2car(tricoord(ga, gb, gc))
gx[idx] <- car$x
gy[idx] <- car$y
} else {
car <- tri2car(organism[position, idx, drop = FALSE])
gx[idx] <- car$x
gy[idx] <- car$y
}
}
xy <- data.frame(x = gx, y = gy)
} else {
xy <- tri2car( organism[position,] )
}
dat <- dplyr::filter(
dplyr::mutate(
tibble::tibble(xy),
stage = organism["stage",],
substrate = substrates[organism["sub.stage",]],
pchar = factor(as.character( future$pch[.data$stage + 1] ), levels = unique(as.character(future$pch))),
color = as.character( future$color[.data$stage + 1] ),
species = species),
.data$substrate %in% show_sub)
attr(dat, "species") <- species
step_val <- if (!is.null(community$step)) community$step else if (!is.null(community$count$step)) community$count$step else if (!is.null(attr(community, "nstep"))) attr(community, "nstep") else step
attr(dat, "step") <- step_val
attr(dat, "layout") <- layout
attr(dat, "width") <- width
attr(dat, "step_density") <- step_density
class(dat) <- c("ewing_substrate", class(dat))
dat
}
ggplot_ewing_substrate <- function(object,
xlab = "horizontal", ylab = "vertical",
layout = attr(object, "layout"),
width = attr(object, "width"),
step_density = attr(object, "step_density"),
layers = c("poly", "hex", "organisms", "centers", "labels"),
...)
{
if (inherits(object, "ggplot")) return(object)
if (is.null(layout)) layout <- "facet"
if (is.null(width)) width <- 10
if (is.null(step_density)) step_density <- 1
species <- attr(object, "species")
step <- attr(object, "step")
# Allows same color for different pchar, but only one color per pchar.
tmp <- dplyr::arrange(
dplyr::distinct(
dplyr::distinct(object, .data$pchar, .data$color),
.data$pchar, .keep_all = TRUE),
.data$pchar)
col.palate <- tmp$color
names(col.palate) <- as.character(tmp$pchar)
if (layout == "hex") {
topo <- substrate_topology(width = width, step = step_density)
sub_obj <- create_substrate(topo, width = width, step = step_density)
hex_overlay <- create_hex_overlay(sub_obj, step = step_density)
p <- ggplot2::ggplot()
if ("poly" %in% layers && nrow(sub_obj$poly) > 0) {
p <- p + ggplot2::geom_polygon(data = sub_obj$poly, ggplot2::aes(x = x, y = y, group = substrate),
fill = NA, color = "black", linewidth = 0.7)
}
if ("hex" %in% layers && nrow(hex_overlay) > 0) {
p <- p + ggplot2::geom_polygon(data = hex_overlay, ggplot2::aes(x = x, y = y, group = cell_id),
fill = NA, color = "gray75", linewidth = 0.3)
}
if ("organisms" %in% layers && nrow(object) > 0) {
p <- p + ggplot2::geom_text(data = object, ggplot2::aes(x = x, y = y, label = pchar, color = pchar),
fontface = "bold", size = 4) +
ggplot2::scale_color_manual(name = "Stage", values = col.palate) +
ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(label = names(col.palate))))
}
if ("centers" %in% layers && nrow(sub_obj$centers) > 0) {
p <- p + ggplot2::geom_text(data = sub_obj$centers, ggplot2::aes(x = x, y = y, label = substrate),
color = "black", fontface = "bold", size = 4.5)
}
if ("labels" %in% layers && nrow(sub_obj$labels) > 0) {
p <- p + ggplot2::geom_text(data = sub_obj$labels, ggplot2::aes(x = x, y = y, label = label),
color = "darkred", fontface = "bold", size = 3.5)
}
return(p + ggplot2::theme_void() +
ggplot2::theme(plot.margin = ggplot2::margin(2, 2, 2, 2, "pt")) +
ggplot2::coord_fixed() +
ggplot2::ggtitle(paste(species, "on Hex Substrate Grid at", step, "steps")))
}
# Default Facet View by Substrate Component
ggplot2::ggplot(object) +
ggplot2::aes(.data$x, .data$y, label = .data$pchar, col = .data$pchar) +
ggplot2::geom_text() +
ggplot2::facet_wrap(~ substrate) +
ggplot2::xlab(xlab) +
ggplot2::ylab(ylab) +
ggplot2::scale_color_manual(name = "Stage", values = col.palate) +
ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(label = names(col.palate)))) +
ggplot2::ggtitle(paste(species, "on substrate at", step, "steps"))
}
autoplot.ewing_substrate <- function(object, ...)
ggplot_ewing_substrate(object, ...)
# --- Source: ewing_ageclass.R ---
ewing_ageclass <- function(community, substrate = TRUE, total = TRUE,
normalize = TRUE, ...) {
nsim_val <- NULL
nstep_val <- NULL
if (inherits(community, "isle_royale_sim")) {
if (is.null(community$history) || nrow(community$history) == 0) return(NULL)
out <- community$history
if (!substrate) {
out <- dplyr::filter(out, .data$Type != "substrate")
}
if (total) {
tot <- dplyr::mutate(
dplyr::ungroup(
dplyr::summarize(
dplyr::group_by(out, .data$Species, .data$step, .data$time, .data$Type),
Count = sum(.data$Count),
.groups = "drop"
)
),
State = "total"
)
out <- dplyr::bind_rows(out, tot)
}
if (normalize) {
out <- dplyr::ungroup(
dplyr::mutate(
dplyr::group_by(out, .data$Species, .data$State, .data$Type),
Count = {
m <- max(.data$Count, na.rm = TRUE)
if (!is.na(m) && m > 0) .data$Count / m else 0
}
)
)
}
ordered_levels <- unique(c("calf", "yearling", "adult", "senior", "pup", "subadult", "total"))
out$State <- factor(out$State, levels = ordered_levels[ordered_levels %in% unique(out$State)])
attr(out, "nstep") <- community$nstep
attr(out, "units") <- "days"
class(out) <- c("ewing_ageclass", class(out))
return(out)
}
if (inherits(community, "ewing_discrete")) {
nsim_val <- attr(community, "nsim")
nstep_val <- attr(community, "nstep")
community <- community[[1]]
}
count <- readCount(community)
if(!length(count)) return(NULL)
species <- names(count)
if(is.null(species)) return(NULL)
ageclass <- list()
for( i in species ) {
ageclass[[i]] <- levels( getOrgFuture( community, i, "ageclass" ))
}
substrates <- list()
for( i in species ) {
substrates[[i]] <- levels( getOrgInteract( community,, i, "substrate" ))
}
out <- list()
for(i in species) {
out[[i]] <- dplyr::mutate(
tidyr::pivot_longer(
tibble::tibble(
as.data.frame(count[[i]])),
dplyr::any_of(c(ageclass[[i]], substrates[[i]])),
names_to = "State",
values_to = "Count"),
Type = ifelse(.data$State %in% substrates[[i]], "substrate", "ageclass"))
}
out <- dplyr::bind_rows(out, .id = "Species")
if(!substrate) {
out <- dplyr::filter(out, .data$Type != "substrate")
}
if(total) {
tot <- dplyr::mutate(
dplyr::ungroup(
dplyr::summarize(
dplyr::group_by(
out,
.data$Species, .data$step, .data$time, .data$future, .data$Type),
Count = sum(.data$Count),
.groups = "drop")
),
State = "total")
out <- dplyr::bind_rows(out, tot)
}
if(normalize) {
out <- dplyr::ungroup(
dplyr::mutate(
dplyr::group_by(
out,
.data$Species, .data$State, .data$Type),
Count = {
m <- max(.data$Count, na.rm = TRUE)
if (!is.na(m) && m > 0) .data$Count / m else 0
}))
}
subs <- if (substrate) unlist(substrates) else NULL
if (length(species) > 1) {
ordered_levels <- unique(c(ageclass[[species[1]]], "total", unlist(ageclass[species[-1]]), subs))
} else {
ordered_levels <- unique(c(unlist(ageclass), "total", subs))
}
out$State <- factor(out$State, levels = ordered_levels)
attr(out, "nstep") <- if (!is.null(nstep_val)) nstep_val else attr(community, "nstep")
if (is.null(attr(out, "nstep")) && !is.null(out$step)) {
attr(out, "nstep") <- max(out$step, na.rm = TRUE)
}
attr(out, "nsim") <- if (!is.null(nsim_val)) nsim_val else attr(community, "nsim")
attr(out, "units") <- if (inherits(community, "isle_royale_sim")) "days" else tryCatch(getOrgFeature(community, species[1], "units"), error = function(e) "time")
class(out) <- c("ewing_ageclass", class(out))
out
}
ggplot_ewing_ageclass <- function(object, main = NULL, title = NULL, x_var = c("step", "time"), time_unit = NULL, ... )
{
x_var <- match.arg(x_var)
if (is.null(time_unit)) {
time_unit <- attr(object, "units")
}
if (is.null(time_unit) || is.na(time_unit) || time_unit == "NA") {
time_unit <- "time"
}
if (is.null(title)) title <- main
if (is.null(title)) {
nstep <- attr(object, "nstep")
if (is.null(nstep) && !is.null(object$step)) {
nstep <- max(object$step, na.rm = TRUE)
}
nsim <- attr(object, "nsim")
unit_str <- if (x_var == "step") "steps" else time_unit
time_hdr <- if (x_var == "step") "Steps" else paste0(toupper(substring(time_unit, 1, 1)), substring(time_unit, 2))
if (!is.null(nstep)) {
if (!is.null(nsim) && nsim > 1) {
title <- paste0("Age Classes over ", time_hdr, " (", nstep, " ", unit_str, ", nsim = ", nsim, ")")
} else {
title <- paste0("Age Classes over ", time_hdr, " (", nstep, " ", unit_str, ")")
}
} else {
title <- paste0("Age Classes over ", time_hdr)
}
}
species_vec <- unique(as.character(object$Species))
if (length(species_vec) == 0) species_vec <- "Organism"
p_list <- list()
for (sp in species_vec) {
df_sp <- object[object$Species == sp, , drop = FALSE]
if (is.factor(df_sp$State)) {
df_sp$State <- droplevels(df_sp$State)
}
x_col <- if (x_var == "step" && "step" %in% names(df_sp)) "step" else "time"
x_lbl <- if (x_var == "step") "steps" else time_unit
sp_title <- paste(toupper(substring(sp, 1, 1)), substring(sp, 2), " Age Classes", sep = "")
p_sub <- ggplot2::ggplot(df_sp, ggplot2::aes(x = .data[[x_col]], y = .data$Count, col = .data$State, group = .data$State)) +
ggplot2::geom_step(na.rm = TRUE, linewidth = 0.8) +
ggplot2::geom_point(size = 2, na.rm = TRUE) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = sp_title,
x = x_lbl,
y = "Count",
color = "Age Class"
) +
ggplot2::theme(
plot.title = ggplot2::element_text(size = 11, face = "bold", hjust = 0.5),
legend.position = "right"
)
p_list[[sp]] <- p_sub
}
if (length(p_list) == 1) {
return(p_list[[1]] + ggplot2::ggtitle(title))
}
grid_plots <- cowplot::plot_grid(plotlist = p_list, ncol = length(p_list), align = "h")
title_widget <- cowplot::ggdraw() +
cowplot::draw_label(title, fontface = 'bold', x = 0.5, hjust = 0.5, size = 13)
cowplot::plot_grid(title_widget, grid_plots, ncol = 1, rel_heights = c(0.12, 1))
}
autoplot.ewing_ageclass <- function(object, x_var = c("step", "time"), ...) {
x_var <- match.arg(x_var)
ggplot_ewing_ageclass(object, x_var = x_var, ...)
}
# --- Source: ewing_envelope.R ---
ewing_envelope <- function(object, species, item, ordinate = "time", increment = 0.5) {
# Pull out `ordinate` and `item` for each run
pulled <-
tidyr::fill(
dplyr::arrange(
tidyr::pivot_wider(
dplyr::bind_rows(
purrr::map(
object,
function(x) {
dplyr::distinct(
as.data.frame(x[[species]][,c(ordinate, item)]),
.data[[ordinate]],
.keep_all = TRUE)
}),
.id = "run"),
names_from = "run",
values_from = item),
.data[[ordinate]]),
-dplyr::matches(ordinate))
out <- GET::create_curve_set(list(r = as.matrix(pulled)[,1],
obs = as.matrix(pulled[,-1])))
class(out) <- c("ewing_envelope", class(out))
attr(out, "count") <- attr(object, "count")
attr(out, "nstep") <- attr(object, "nstep")
attr(out, "nsim") <- attr(object, "nsim")
attr(out, "species") <- species
attr(out, "item") <- item
attr(out, "ordinate") <- ordinate
out
}
ewing_envelopes <- function(object) {
species <- attr(object, "species")
items <- attr(object, "items")
ordinate <- attr(object, "ordinate")
nstep <- attr(object, "nstep")
count <- attr(object, "count")
nsim <- attr(object, "nsim")
confidence <- (nsim > 2)
envs <- as.list(species)
names(envs) <- species
if(confidence) {
confs <- envs
} else {
confs <- NULL
}
for(specy in species) {
env1 <- as.list(items[[specy]])
if(confidence) {
conf1 <- env1
}
for(item in items[[specy]]) {
env1[[item]] <- ewing_envelope(object, specy, item, ordinate)
if(confidence) {
conf1[[item]] <- GET::fBoxplot(env1[[item]], type = 'area')
}
}
envs[[specy]] <- env1
if(confidence) {
confs[[specy]] <- conf1
}
}
object <- list(env = envs, conf = confs)
class(object) <- c("ewing_envelopes", class(object))
attr(object, "species") <- species
attr(object, "items") <- items
attr(object, "ordinate") <- ordinate
attr(object, "nstep") <- nstep
attr(object, "count") <- count
attr(object, "nsim") <- nsim
attr(object, "confidence") <- confidence
object
}
summary.ewing_envelopes <- function(object, species = NULL, verbose = TRUE, ...) {
# object$conf[[specy]][[item]] is time by 6-num boxplot summary
if(verbose) {
nstep <- attr(object, "nstep")
count <- attr(object, "count")
nsim <- attr(object, "nsim")
cat(nsim, "Runs of ",
nstep, "Steps for",
paste(names(count), count, sep = "=", collapse = ", "), "\n")
}
out <- print(object, species, ...)
if(!is.null(out)) {
out <- dplyr::ungroup(
dplyr::filter(
dplyr::group_by(
out,
.data$species, .data$item),
(.data$r == 0) | (.data$r == max(.data$r))))
}
out
}
print.ewing_envelopes <- function(x, species = NULL, ...) {
# x$conf[[specy]][[item]] is time by 6-num boxplot summary
if(is.null(x$conf)) {
return(NULL)
}
out <- dplyr::bind_rows(
purrr::map(
x$conf,
function(x) {
# somehow get summary across species and items using as.data.frame
x <- x[names(x) != ""]
dplyr::bind_rows(
purrr::map(
x,
as.data.frame),
.id = "item")
}),
.id = "species")
if(!is.null(species)) {
if(species %in% unique(out$species)) {
sp <- species
out <- dplyr::filter(out, species == sp)
}
}
dplyr::mutate(out, dplyr::across(where(is.numeric), function(x) pmax(x,0)))
out
}
ggplot_ewing_envelopes <- function(object, confidence = FALSE, main = "", ...) {
if(inherits(object, "ewing_discrete")) {
object <- ewing_envelopes(object)
}
species <- attr(object, "species")
items <- attr(object, "items")
ordinate <- attr(object, "ordinate")
nstep <- attr(object, "nstep")
count <- attr(object, "count")
nsim <- attr(object, "nsim")
confidence <- confidence & attr(object, "confidence")
patch <- list()
for(specy in species) {
p <- list()
for(item in items[[specy]]) {
if(confidence) {
p[[item]] <- plot(object$conf[[specy]][[item]], main = main) +
ggplot2::labs(x = "time", y = item) +
ggplot2::ggtitle(main) +
ggplot2::ylim(0, NA)
} else {
p[[item]] <- ggplot_ewing_envelope(object$env[[specy]][[item]])
}
}
patch[[specy]] <- cowplot::plot_grid(plotlist = p, nrow = length(p))
}
# NEED TO get attribute count and nstep in here
patch <- cowplot::plot_grid(plotlist = patch, ncol = length(patch))
# Add a title. <https://wilkelab.org/cowplot/articles/plot_grid.html>
title <- cowplot::ggdraw() +
cowplot::draw_label(
paste(nsim, "Runs of ", nstep, "Steps for",
paste(species, count, sep = "=", collapse = ", ")),
x = 0, hjust = 0
) +
ggplot2::theme(
# add margin on the left of the drawing canvas,
# so title is aligned with left edge of first plot
plot.margin = ggplot2::margin(0, 0, 0, 7)
)
cowplot::plot_grid(title, patch, ncol = 1, rel_heights = c(0.1, 1))
}
ggplot_ewing_envelope <- function(object, cols = c("#21908CFF", "#440154FF", "#5DC863FF"),
main = "", ...) {
# Kludge. GET::forder needs at least 3 points; cols can be at most length(object).
lcols <- length(cols)
nsim <- ncol(object$funcs)
if(nsim >= min(3, lcols)) {
A <- GET::forder(object, measure = 'area')
lcols <- min(lcols, length(object))
idx <- order(A)[seq_len(lcols)]
cols <- cols[seq_len(lcols)]
} else {
lcols <- nsim
idx <- seq_len(lcols)
cols <- cols[seq_len(lcols)]
}
item <- attr(object, "item")
species <- attr(object, "species")
ordinate <- attr(object, "ordinate")
if(length(object) >= 50) {
p <- plot(object, idx = idx, col_idx = cols, main = main)
} else {
p <- plot(object)
}
p +
ggplot2::labs(x = ordinate, y = item) +
ggplot2::ggtitle(paste(species, item))
}
autoplot.ewing_envelope <- function(object, ...) {
ggplot_ewing_envelope(object, ...)
}
# --- Source: ewing_snapshot.R ---
ewing_snapshot <- function(object, step = 0, ...)
{
out <- list(step = step,
ageclass = ewing_ageclass(object, ...))
species <- get.species(object)
subs <- list()
for(j in species) {
subs[[j]] <- ewing_substrate(object, j, step = step, ...)
}
out$substrate <- subs
class(out) <- c("ewing_snapshot", "ewing", "list")
out
}
# --- Source: temp.R ---
## $Id: temp.R,v 1.0 2002/12/09 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
###########################################################################################
## initTemp( community, lo.hour, hi.hour ) ### creates Temperature object
## activeTemp( community, lo.hour, hi.hour ) ### updates Temperature object
## checkTime( ) # check if activeTemp needs updating to cover time interval
##
## temp.design( )
## temp.plot( )
###########################################################################################
rescale.temp <- function( low, high, newry = ry, newrx = rx,
ry = range( c( stats::predict( low, low$knots )$y, stats::predict( high, high$knots )$y )))
{
if( max( abs( ry - newry )) > 0 ) {
tmpy <- diff( newry ) / diff( ry )
low$coefficients <- low$coefficients * tmpy
high$coefficients <- high$coefficients * tmpy
tmpy <- newry[1] - ry[1] * tmpy
low$coefficients[,1] <- tmpy + low$coefficients[,1]
high$coefficients[,1] <- tmpy + high$coefficients[,1]
}
rx <- range( low$knots, high$knots )
if( max( abs( rx - newrx )) > 0 ) {
newx <- diff( newrx ) / diff( rx )
low$knots <- newrx[1] + newx * low$knots
high$knots <- newrx[1] + newx * high$knots
tmp <- 1
for( i in seq( 2, ncol( low$coefficients ))) {
tmp <- tmp * newx
low$coefficients[,i] <- low$coefficients[,i] / tmp
high$coefficients[,i] <- high$coefficients[,i] / tmp
}
}
list( low = low, high = high )
}
###########################################################################################
### Temperature data structure
### Min minimum temperature for degree-day computations
### Time list of hours when temperature changes through a day
### Base list of temperature shifts through a day
### list names indicate first day Hour and Base applies
### Low spline fit for daily low temperatures
### High spline fit for daily high temperatures
### DegreeDay spline fit (hour->DD) for currently active days
### Hour ramped backspline fit (DD->hour) for currently active days
### Idea is that DegreeDay and Hour are updated whenever a future event is scheduled
### past the latest future event already scheduled. At that time, the first
### and last time of fit are both adjusted. This can be done without rebuilding
### the whole spline fit by appropriate adjustment of knots and coefficents.
### Steps: (1) combine Hour, Low and High to get knots
### (2) build spline coefficients
###########################################################################################
temp.spline <- function( community, hour, temp, start = 0,
mintemp = getTemp( community, "Min" ), cumulative = TRUE,
mult = getTemp( community, "Unit" ) )
{
if( is.data.frame( hour ) & missing( temp )) {
temp <- hour$temp
hour <- hour$hour
}
## drop low value if ones on either side are low
low <- temp <= mintemp
toolow <- low & c(FALSE,low[-length(low)]) & c(low[-1],FALSE)
hour <- hour[!toolow]
temp <- temp[!toolow]
lh <- length( hour )
## expand single low in middle to two
low <- temp <= mintemp
toolow <- low & c(FALSE,!low[-length(low)]) & c(!low[-1],FALSE)
if( any( toolow )) {
hour <- c(hour,hour[toolow])
temp <- c(temp,temp[toolow])
temp <- temp[ order( hour ) ]
hour <- sort( hour )
}
## truncate on left
low <- temp <= mintemp & c(temp[-1] > mintemp,FALSE)
low1 <- c(FALSE,low[-lh])
if( any( low ))
hour[low] <- hour[low] + ( hour[low1] - hour[low] ) * ( mintemp - temp[low] ) /
( temp[low1] - temp[low] )
## truncate on right
low <- temp <= mintemp & c(FALSE,temp[-lh] > mintemp)
low1 <- c(low[-1],FALSE)
if( any( low ))
hour[low] <- hour[low] + ( hour[low1] - hour[low] ) * ( mintemp - temp[low] ) /
( temp[low1] - temp[low] )
temp[ temp < mintemp ] <- mintemp
temp <- temp[ !duplicated( hour ) ]
hour <- unique( hour )
## linear interpolating spline for temp
lh <- length( hour ) - 1
h0 <- 1:lh
aa <- ( temp[1+h0] - temp[h0] ) / ( hour[1+h0]-hour[h0] )
aa <- c( aa, aa[lh] )
bb <- temp - mintemp
s <- list( knots = hour, coefficients = cbind( bb,aa,0,0) )
dimnames( s$coefficients ) <- list( paste( floor( hour / mult ), round( hour %% mult ),
sep = "." ), c("const","linear","quad","cubic") )
attr(s,"formula") <- temp ~ hour
class( s ) <- c("npolySpline","polySpline","spline")
if( !cumulative )
return( s )
## quadratic spline for cumulative temp
aa <- s$coefficients[,3] <- s$coefficients[,2] / 2 / mult / mult
lh <- length( aa ) - 1
h0 <- 1:lh
h1 <- 1 + h0
hour <- s$knots
bb <- s$coefficients[,2] <- s$coefficients[,1] / mult
cc <- ( hour[h1] - hour[h0] ) * ( bb[h0] + aa[h0] *
( hour[h1] - hour[h0] ))
s$coefficients[,1] <- start + c( 0, cumsum( cc ))
s
}
###########################################################################################
temp.repeat <- function( community, period = range( days ))
{
lodays <- range( getTemp( community, "Low" )$knots )
hidays <- range( getTemp( community, "High" )$knots )
days <- c( max( lodays[1], hidays[1] ), min( lodays[2], hidays[2] ))
lodays <- period[1] >= days[1] & period[1] <= days[2]
if( lodays )
days[1] <- period[1]
hidays <- period[2] <= days[2] & period[2] >= days[1]
if( hidays )
days[2] <- period[2]
if( !( lodays & hidays )) {
stop( paste( "\nSimulation period is outside of Temperature days:\n period =",
paste( round( period ), collapse = "," ), "; days =",
paste( round( days ), collapse = "," ), "\nNeed to run temp.design() and start over!" ))
}
period <- days
days <- seq( days[1], days[2] )
low <- stats::predict( getTemp( community, "Low" ), days )$y
high <- stats::predict( getTemp( community, "High" ), days )$y
periods <- c( as.numeric( names( getTemp( community, "Time" ))), Inf )
temps <- hours <- numeric( )
period[2] <- period[2] + 1
for( i in seq( length( periods ) - 1 )) {
day <- max( period[1], periods[i] )
this.period <- days >= day & days < min( period[2], periods[i+1] )
n <- sum( this.period )
if( n ) {
if( !is.null( getTemp( community, "Time", i ))) {
time <- getTemp( community, "Time", i )
base <- getTemp( community, "Base", i )
}
h <- length( time )
lotemp <- rep( low[this.period], rep(h,n) )
hitemp <- rep( high[this.period], rep(h,n) )
temp <- ( base - min( base )) / diff( range( base ))
temps <- c( temps, rep( temp, n ) * ( hitemp - lotemp ) + lotemp )
hours <- c( hours, rep( time, n ) + rep( getTemp( community, "Unit" ) *
seq( day, day + n - 1 ), rep(h,n) ))
}
}
day <- days == period[2]
if( any( day )) {
temps <- c( temps, low[day] + temp[1] * ( high[day] -
getTemp( community, "Low" )[day] ))
hours <- c( hours, getTemp( community, "Unit" ) * period[2] )
}
data.frame( hour = hours, temp = temps )
}
##########################################################################################
showTemp <- function( community )
{
cat( "Temperature set for days",
paste( range( getTemp( community, "Low" )$knots ), collapse = " to " ), "\n" )
cat( "Daily low temperature range:",
paste( round( range( stats::predict( getTemp( community, "Low" ))$y )),
collapse = " to " ), "\n" )
cat( "Daily high temperature range:",
paste( round( range( stats::predict( getTemp( community, "High" ))$y )),
collapse = " to " ), "\n" )
if( !is.null( getTemp( community, "DegreeDay" ) ))
cat( "Active temperature range:",
paste( round( range( getTemp( community, "DegreeDay", "knots" ) )),
collapse = " to " ), "\n" )
cat( "Run temp.design() to adjust temperature range\n" )
}
##########################################################################################
activeTemp <- function( community,
lo.hour = min( getTemp( community, "DegreeDay", "knots" ) ),
hi.hour = max( getTemp( community, "DegreeDay", "knots" ) ),
degreeday, messages = TRUE )
{
unit <- getTemp( community, "Unit" )
if( missing( hi.hour )) {
hi.hour <- if( length( lo.hour ) > 1 )
lo.hour[2]
else
lo.hour + unit
}
## period is in units of hours
period <- c( floor( lo.hour[1] / unit ), ceiling( hi.hour / unit ))
if( missing( degreeday )) {
degreeday <- if( is.null( getTemp( community, "DegreeDay" ) ))
0
else
getDegreeDay( community, lo.hour[1] )
}
if( is.na( degreeday ))
degreeday <- 0
community <- setTemp( community, "DegreeDay",
temp.spline( community, temp.repeat( community, period ),
start = degreeday ))
community <- setTemp( community, "Hour",
break.backSpline( getTemp( community, "DegreeDay" )))
if(messages) {
tmp <- period * unit
names( tmp ) <- c("lo.hour","hi.hour")
print( tmp )
}
community
}
###########################################################################################
updateTemp <- function( community,
period = range( getTemp( community, "DegreeDay", "knots" ) ))
{
## period is in units of hours
first <- period[1]
last <- period[2]
s <- getTemp( community, "DegreeDay" )
knots <- s$knots
## drop earlier times that are now in the past
drop <- sum( knots < first ) - 1
nk <- length( knots )
change <- drop > 0
if( change ) {
s$knots <- knots[ - seq( drop ) ]
s$coefficients <- s$coefficients[ - seq( drop ), ]
}
## add new days to include last
period <- ceiling( c( knots[nk], last ) / getTemp( community, "Unit" ) )
if( period[1] <= period[2] ) {
change <- TRUE
news <- temp.spline( community, temp.repeat( community, period ),
start = getHour( getTemp( community, "Unit" ) * period[1] ))
nk <- length( s$knots )
if( s$knots[nk] == news$knots[1] ) {
s$knots <- s$knots[-nk]
s$coefficients <- s$coefficients[-nk,]
}
else
s$coefficients[nk,] <- news$coefficients[1,]
s$knots <- c( s$knots, news$knots )
s$coefficients <- rbind( s$coefficients, news$coefficients )
}
if( change ) {
community <- setTemp( community, "DegreeDay", s )
community <- setTemp( community, "Hour", break.backSpline( s ))
}
community
}
###########################################################################################
ramp.backSpline <- function( s )
{
## ramped backspline is a trick to get backSpline when curve is flat in spots
## if tmp <- stats::predict( ramp.backSpline( s ))
## then plot tmp$y versus tmp$x-tmp$y to "recover" original curve.
## problem is that one cannot recover particular x this way!
s$coefficients[,1] <- s$coefficients[,1] + s$knots
s$coefficients[,2] <- s$coefficients[,2] + 1
splines::backSpline( s )
}
###########################################################################################
break.backSpline <- function( tmp )
{
## alternative to ramped backSpline that first removes flat regions (slope 0)
## problem remains that leftover may still have slope 0 at a point
## this can cause anomolous results!
## find flat regions and cut out
tmpc <- diff(tmp$coefficients[,1]) == 0
tmpk <- cumsum( diff(tmp$knots) * ( tmpc ))
tmpk <- c(tmpk,tmpk[length(tmpk)])
tmp$knots <- tmp$knots - tmpk
tmpc <- c(!tmpc,TRUE)
tmp$knots <- tmp$knots[tmpc]
tmp$coefficients <- tmp$coefficients[tmpc,]
tmp$breaks
nb <- dim( tmp$coefficients )
## kludge for backspline: cannot handle slope of zero
tmpn <- tmp$coefficients[,2] == 0
if( any( tmpn )) {
tmpn <- seq( tmpn )[ tmpn ]
tmp$coefficients[ tmpn, 2:nb[2] ] <-
( tmp$coefficients[ pmin( nb[1], tmpn + 1 ), 2:nb[2] ] +
tmp$coefficients[ pmax( 1, tmpn - 1 ), 2:nb[2] ] ) / 2
}
## back spline
tmpb <- splines::backSpline(tmp )
if( any( tmpc )) {
## shift back spline based on breaks
tmpb$coefficients[,1] <- tmpb$coefficients[,1] + tmpk[tmpc]
}
tmpb
}
###########################################################################################
getDegreeDay <- function( community, hour )
{
stats::predict( getTemp( community, "DegreeDay" ), hour )$y
}
###########################################################################################
getHour <- function( community, dd )
{
stats::predict( getTemp( community, "Hour" ), dd )$y
}
###########################################################################################
getTime <- function( community, species, x )
{
if( is.na(x)){
cat("getTime missing value\n")
browser()
}
switch( getOrgFeature( community, species, "units" ),
## organisms on hour basis assumed to be active only 6am-6pm
hr = getDegreeDay( community,
( getTemp( community, "Unit" ) / 2 + x + floor( x / getTemp( community, "Unit" ) )) / 2 ),
DD = x )
}
###########################################################################################
checkTime <- function( community, x, base, units )
{
## NOTE: Sometimes base can be negative!! (reset to 0)
if( units == "hr" ) {
x <- max( x )
# print( c( base = base, x = x, knots = range( getTemp( community, "DegreeDay", "knots" ))))
if( x > max( getTemp( community, "DegreeDay", "knots" )))
community <- activeTemp( community, max( base, 0 ),
x + getTemp( community, "Unit" ), messages = FALSE )
}
community
}
###########################################################################################
transTime <- function( community, org1name, org2name, x,
unit1 = getOrgFeature( community, org1name, "units" ),
unit2 = getOrgFeature( community, org2name, "units" ))
{
if( is.na(x)){
cat("transTime missing value\n")
browser()
}
if( unit1 == unit2 )
return( x )
switch( unit1,
## organisms on hour basis assumed to be active only 6am-6pm
hr = getDegreeDay( community,
( getTemp( community, "Unit" ) / 2 + x + floor( x / getTemp( community, "Unit" ))) / 2 ),
DD = getHour( x ))
}
###########################################################################################
### To do:
### 1. interactive designer for hourly temp fluctation over one day
### 2. interactive designer for daily lows and highs over season: temp.design() DONE
### 3. check future event trees hour vs. DD
### 4. schedule interaction events with hour-DD translation
### 5. make aphytis dormant at night
###########################################################################################
# --- Source: initTemp.R ---
##########################################################################################
### simulation temperature administration
##########################################################################################
initTemp <- function( community, lo.hour = 0, hi.hour = getTemp( community, "Unit" ),
days = TemperaturePar["Days"],
messages = TRUE, datafile = "", ... )
{
if(messages) {
cat( "Initializing Temperature Profile ...\n" )
}
Temperature <- list()
TemperaturePar <- getOrgData(community, "temperature", "par", messages, datafile)
# mydata( "TemperaturePar", getOrgInfo( community, "package" ), messages = messages)
TemperaturePar <- array( TemperaturePar[,"value"],
dimnames = list( row.names( TemperaturePar )))
Temperature$Unit <- TemperaturePar["Unit"]
Temperature$Min <- TemperaturePar["Min"]
## set up daily temperature base
TemperatureBase <- getOrgData(community, "temperature", "base", messages, datafile)
# mydata( "TemperatureBase", getOrgInfo( community, "package" ), messages = messages)
Temperature$Time <- split( TemperatureBase$Time, TemperatureBase$Day )
Temperature$Base <- split( TemperatureBase$Base, TemperatureBase$Day )
community$temp <- Temperature
tmp <- seq( lo.hour / Temperature$Unit,
days + 1 + ( hi.hour / Temperature$Unit ),
length = TemperaturePar["Length"] )
tmp0 <- seq( 0, 1, length = TemperaturePar["Length"] )
tmp1 <- 0.25 * ( TemperaturePar["HighBeg"] - TemperaturePar["LowBeg"] )
Temperature$Low <- splines::interpSpline( tmp, TemperaturePar["LowBeg"] * ( 1 - tmp0 ) +
TemperaturePar["LowEnd"] * tmp0 +
sin( pi * 4 * tmp0 ) * tmp1 )
Temperature$High <- splines::interpSpline( tmp, TemperaturePar["HighBeg"] * ( 1 - tmp0 ) +
TemperaturePar["HighEnd"] * tmp0 +
sin( pi * ( 0.125 + 4 * tmp0 )) * tmp1 )
Temperature$DegreeDay <- NULL
if(messages) {
cat( "Base daily temperature fluctuation:\n" )
}
for( i in names( Temperature$time )) {
cat( "From day", i, ":\n" )
tmp <- Temperature$Base[[i]]
names( tmp ) <- Temperature$Time[[i]]
print( tmp )
}
community$temp <- Temperature
if(messages) {
showTemp( community )
cat( "Initial active temperature:\n" )
}
activeTemp( community, lo.hour, hi.hour, getTemp( community, "Time", 1 )[1],
messages = messages)
}
###########################################################################################
getTemp <- function( community, element, sub )
{
tempelem <- community$temp[[element]]
if( !missing( sub ))
tempelem <- tempelem[[sub]]
tempelem
}
###########################################################################################
setTemp <- function( community, element, value )
{
community$temp[[element]] <- value
community
}
# --- Source: temp.design.R ---
temp.design <- function( community, nspline = 8, n = 1, horizontal = TRUE,
col = c(low="blue",high="red") )
{
low <- getTemp( community, "Low" )
high <- getTemp( community, "High" )
is.data <- !is.null( low )
if( !is.data ) {
tmp <- seq( 0, 60, length = nspline )
low <- interpSpline( tmp, 60 + 0.125 * tmp + sin( 0.25 * tmp ))
high <- interpSpline( tmp, 70 + 0.15 * tmp + sin(( pi / 8 ) + 0.225 * tmp ))
}
## plot curve and surrounding axes
par( mfrow = c(1,1), mar = rep(4.1,4))
cat( "Switch to Graphic Screen to Adjust High and Low Temperatures\n" )
plotit <- function( low, high, fig = fig, horizontal = FALSE, strip = .25, margin = 0 )
{
lowpred <- stats::predict( low, low$knots )
highpred <- stats::predict( high, high$knots )
ylim <- range( c( lowpred$y, highpred$y ))
xlim <- range( c( lowpred$x, highpred$x))
xlim <- xlim + c(-1,1) * margin * diff( xlim )
if( horizontal ) {
if( diff( ylim ) == 0 )
ylim <- ylim * c(.75,1.25)
separator <- ylim[2]
ylim[2] <- ylim[2] + strip * diff( ylim )
}
else {
separator <- xlim[2]
xlim[2] <- xlim[2] + strip * diff( xlim )
}
axt <- c("n","s")
tmpar <- par( xaxt = axt[1+horizontal], yaxt = axt[2-horizontal] )
plot( xlim, ylim, xlim = xlim, ylim = ylim, type="n", xlab = "", ylab = "" )
par( xaxt = "s", yaxt = "s" )
title( fig )
mtext( "day", 1, 2 )
mtext( "temp", 2, 2 )
if( horizontal ) {
p <- pretty( c(ylim[1],separator) )
axis( 2, p[ p <= separator ] )
abline( h = separator, lty = 2 )
}
else {
p <- pretty( c(xlim[1],separator) )
axis( 1, p[ p <= separator ] )
abline( v = separator, lty = 2 )
}
points( highpred$x, highpred$y, lwd = 4 )
curve.plot( highpred, n = n, action = "refresh", fit = high, backfit = FALSE,
save.ends = 3, col = "red", lwd = 2 * ( 1 + ( fig == "high" )))
points( lowpred$x, lowpred$y, lwd = 4 )
curve.plot( lowpred, n = n, action = "refresh", fit = low, backfit = FALSE,
save.ends = 3, col = "blue", lwd = 2 * ( 1 + ( fig == "low" )))
separator
}
## place commands along right strip of plot, highlighting current command
plotcmd <- function( ans, fig, cmds, cmdlocs, usr, col = "green", rest = "black",
horizontal = TRUE, data = FALSE )
{
ans <- c( ans, fig )
if( data )
ans <- c( ans, "data" )
tmp <- is.na( match( cmds, ans ))
if( any( tmp )) for( i in unique( cmdlocs$adj )) {
tmpi <- tmp & i == cmdlocs$adj
if( any( tmpi ))
text( cmdlocs$x[tmpi], cmdlocs$y[tmpi], cmds[tmpi], col = rest, adj = i )
}
if( any( !tmp )) for( i in unique( cmdlocs$adj )) {
tmpi <- !tmp & i == cmdlocs$adj
if( any( tmpi ))
text( cmdlocs$x[tmpi], cmdlocs$y[tmpi], cmds[tmpi], col = col, adj = i )
}
}
cmds <- c("add","delete","replace","rescale","","finish","refresh","restart",
"","data","high","low")
newlocs <- if( horizontal )
function( cmds, data = FALSE, usr )
{
if( !data )
cmds <- cmds[ cmds != "data" ]
n <- length( cmds )
blank <- seq( n )[ cmds == "" | cmds == " " ]
tmp <- diff(usr[3:4]) / 20
m <- mean( usr[1:2] )
y <- usr[4] + 0.5 * tmp - c( tmp * seq( blank[1] - 1 ), 0,
tmp * seq( blank[2] - blank[1] - 1 ), 0,
tmp * seq( n - blank[2] ))
x <- c( rep( usr[1], blank[1] - 1 ), mean( m, usr[1] ),
rep( m, blank[2] - blank[1] - 1 ), mean( m, usr[2] ),
rep( usr[2], n - blank[2] ))
adj <- c( rep( 0, blank[1] ),
rep( 0.5, blank[2] - blank[1] ),
rep( 1, n - blank[2] ))
tmp <- data.frame( x = x, y = y, adj = adj )
cmds[blank[2]] <- " "
row.names( tmp ) <- cmds
tmp
}
else
function( cmds, data = FALSE, usr )
{
if( !data )
cmds <- cmds[ cmds != "data" ]
n <- length( cmds )
blank <- seq( n )[ cmds == "" | cmds == " " ]
tmp <- diff(usr[3:4]) / 20
m <- mean( usr[3:4] )
tmp <- c( usr[4] - tmp * seq( blank[1] - 1 ),
mean( m, usr[4] ),
m + tmp * ( seq( blank[1] + 1, blank[2] - 1 ) - mean( blank )),
mean( m, usr[3] ),
usr[3] + tmp * seq( n - blank[2] ))
tmp <- data.frame( x = rep( usr[2], n ), y = tmp, adj = rep( 1, n ))
cmds[blank[2]] <- " "
row.names( tmp ) <- cmds
tmp
}
newans <- ans <- "replace"
par( mar = c(4.1,4.1,3.1,4.1),omi=rep(.25,4))
fig <- "low"
fit <- low
separator <- plotit( low, high, fig, horizontal )
usr <- par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
cmds <- row.names( cmdlocs )
use.data <- FALSE
plotcmd( ans, fig, cmds, cmdlocs, usr, data = use.data )
rescale.data <- c( range( stats::predict( low, low$knots )$y,
stats::predict( high, high$knots )$y ), range( low$knots, high$knots ))
repeat {
## get command from plot using cursor
z <- locator(1,"n")
if(( !horizontal & z$x > separator ) | ( horizontal & z$y > separator )) {
if( horizontal ) { # need to look at both z&y
x <- abs( z$x - cmdlocs$x )
x <- x == min( x )
newans <- cmds[x]
z <- abs( z$y - cmdlocs$y )[x]
newans <- newans[ z == min( z ) ][1]
}
else {
z <- abs(z$y - cmdlocs$y )
newans <- cmds[ z == min( z ) ][1]
}
switch( newans,
finish =, refresh = {
separator <- plotit( low, high, fig, horizontal )
usr <- par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
data = {
use.data <- is.data & !use.data
if( is.data & !use.data )
plotit( low, high, fig, horizontal )
},
high =, low = {
if( newans != fig )
fig <- newans
separator <- plotit( low, high, fig, horizontal )
usr <- par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
rescale = {
tmpcmds <- cmds
tmpans <- tmpcmds[ cmds == "rescale" ] <-
"rescale: Switch to Character Screen"
plotcmd( tmpans, fig, tmpcmds, cmdlocs, usr, data = use.data )
ry <- range( c( stats::predict( low, low$knots )$y, stats::predict( high, high$knots )$y ))
cat( "\nEnter new minimum/maximum followed by RETURN key\n" )
newry <- ry
show <- c("minimum","maximum")
change <- FALSE
for( i in 1:2 ) {
tmpy <- readline( paste( show[i], " temp (",
round( ry[i], 2 ), "):", sep = "" ))
tmpy <- if( tmpy == "" ) NA
else as.numeric( tmpy )
if( is.na( tmpy ))
tmpy <- ry[i]
else
change <- TRUE
rescale.data[i] <- tmpy
}
tmp <- range( c( low$knots, high$knots ))
for( i in 1:2 ) {
tmpx <- readline( paste( show[i], " time (",
round( tmp[i], 2 ), "):", sep = "" ))
tmpx <- if( tmpx == "" ) NA
else as.numeric( tmpx )
if( is.na( tmpx ))
tmpx <- tmp[i]
else
change <- TRUE
rescale.data[2+i] <- tmpx
}
cat( "Switch to Graphic Screen to Adjust High and Low Temperatures\n" )
if( change ) {
tmp <- rescale.temp( low, high, rescale.data[1:2],
rescale.data[3:4], ry )
low <- tmp$low
high <- tmp$high
}
separator <- plotit( low, high, fig, horizontal )
usr <- par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
restart = {
low <- getTemp( community, "Low" )
high <- getTemp( community, "High" )
if( is.data )
rescale.data <- c( range( stats::predict( low, low$knots )$y,
stats::predict( high, high$knots )$y ), range( low$knots, high$knots ))
separator <- plotit( low, high, fig, horizontal )
usr <- par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
add =, delete =, replace = {
ans <- newans
}
)
if( use.data ) {
tmp <- rescale.temp( getTemp( community, "Low" ), getTemp( community, "High" ),
rescale.data[1:2], rescale.data[3:4] )
for( i in c("low","high") ) {
datax <- tmp[[i]]$knots
lines( datax, stats::predict( tmp[[i]], datax )$y )
}
}
plotcmd( ans, fig, cmds, cmdlocs, usr, data = use.data )
}
else { # modify the curve knots
fit <- get( fig )
if( z$x >= min( fit$knots ) & z$x <= max( fit$knots )) {
fit <- curve.plot( as.data.frame( stats::predict( fit, fit$knots )), n = n, action = ans,
z = z, fit = fit, backfit = FALSE, save.ends = 2, col = col[fig] )
assign( fig, fit$fit )
}
ans <- "replace"
plotcmd( ans, fig, cmds, cmdlocs, usr, data = use.data )
}
if( newans == "finish" )
break
}
plotcmd( newans, fig, cmds, cmdlocs, "red", data = use.data )
community <- setTemp( community, "Low", low )
community <- setTemp( community, "High", high )
activeTemp( community )
}
# --- Source: temp.plot.R ---
temp.plot <- function( community, lo.hour = s$knots[1], hi.hour = max( s$knots ),
length = 201,
col = NULL, derivative = FALSE, ..., printit = FALSE )
{
s <- getTemp( community, "DegreeDay" )
x <- seq( lo.hour, hi.hour, length = length )
## make sure to pick up knots to plot in this region
x <- unique( sort( c( x, s$knots[ s$knots >= lo.hour & s$knots <= hi.hour ] )))
ylab <- "degree-days"
if( derivative ) {
ylab <- "degrees above min"
s$coefficients <- s$coefficients[,-1]
for( i in seq( 2, ncol( s$coefficients )))
s$coefficients[,i] <- s$coefficients[,i] * i
}
y <- stats::predict( s, x )$y
plot( x / getTemp( community, "Unit" ), y, type = "l", xlab = "day",
ylab = ylab, ... )
if( printit )
print( cbind( hi.hour, stats::predict( s, hi.hour )$y ))
if( !is.null( col ))
points( s$knots, coef(s)[,1], col = col )
if( !derivative ) {
s <- getTemp( community, "Hour" )
x <- seq( min( y ), max( y ), length = length )
tmp <- stats::predict( s, x )
lines( tmp$y / getTemp( community, "Unit" ), tmp$x, col = "blue" )
}
}
###########################################################################################
temp.lines <- function( s, mult = 24, col = "red" )
{
x <- seq( s$knots[1], max( s$knots ), length = 51 )
x <- unique( sort( c( x, s$knots )))
p <- stats::predict( s, x )
lines( p$x, p$y, col = col )
if( !is.null( col ))
points( s$knots / mult, coef(s)[,1], col = col )
}
# --- Source: my.R ---
###########################################################################################
## System files
###########################################################################################
my.read <- function(dataname, stringsAsFactors = TRUE)
{
switch(tools::file_ext(dataname),
"txt" =, "tsv" = read.table(dataname, stringsAsFactors = stringsAsFactors),
"csv" = read.csv(dataname, stringsAsFactors = stringsAsFactors),
"xls" = readxl::read_xls(dataname),
"xlsx" = readxl::read_xlsx(dataname))
}
###########################################################################################
my.eval <- function(species, extension, element, checkdata = FALSE )
{
if( !missing( extension ))
species <- paste( species, extension, sep = ".")
if( exists( species ))
organism <- get( species )
else {
if( checkdata ) {
organism <- utils::data( list = species )
if( organism == species )
organism <- NULL
}
else
organism <- NULL
}
if(!( missing(element) | is.null( organism )))
organism <- organism[[element]]
organism
}
###########################################################################################
mydata <- function( dataname, package, restart = FALSE, messages = TRUE )
{
edata <- exists( dataname )
if( restart & edata ) {
remove( list = dataname, pos = 1 )
edata <- !edata
}
if( !edata ) {
utils::data( list = dataname, package = eval( package ))
if(messages) {
cat( "Data", dataname, "loaded\n" )
}
}
else
if(messages) {
cat( "Data", dataname, "already loaded\n" )
}
}
# --- Source: Org.R ---
###########################################################################################
### Organism Features
###########################################################################################
getOrgFeature <- function( community, species, feature = names( OrgFeature ))
{
OrgFeature <- getOrgInfo( community, "Feature" )
if (!is.null(OrgFeature)) {
if( missing( species ))
return( row.names( OrgFeature ))
f <- OrgFeature[ species, feature ]
if( length( feature ) == 1 ) {
if( any( is.na( f )))
return( NA )
f <- as.character( f )
}
else {
f <- apply( f, 2, as.character )
}
f <- c( unlist( f ))
numf <- suppressWarnings(as.numeric( f ))
if( all( !is.na( numf )))
f <- numf
return(f)
}
# Fallback for webR standalone adapter objects
if (missing(species)) return(c("host", "parasite"))
if (missing(feature) || is.null(feature)) return(c(substrate = "substrate", units = "units"))
if (length(feature) == 1) {
if (feature == "substrate") return("substrate")
if (feature == "units") return("units")
return(NA)
}
res <- rep("substrate", length(feature))
names(res) <- feature
res
}
##########################################################################################
getOrgHosts <- function( community, species,
feature = c("offspring","attack","substrate") ###HOST SPECIFIC###
)
{
# This seems overly complicated and adds substrate; maybe function name is wrong
f <- unique( getOrgFeature( community, species, feature ))
# f[match(species, f, nomatch = 0)] # this would only get species
o <- getOrgFeature( community )
o[ match( f, o, nomatch = 0 ) ]
}
###########################################################################################
getOrgFuture <- function( community, species, feature, current,
future = OrgFuture[[species]] )
{
OrgFuture <- getOrgInfo( community, "Future" )
if (!is.null(OrgFuture) && !is.null(OrgFuture[[species]])) {
future <- OrgFuture[[species]]
if( missing( current )) {
if( missing( feature ))
return( future )
future <- future[,feature]
}
else {
if( !is.numeric( current ))
current <- match( current, future$current, nomatch = 0 )
if( missing( feature ))
future <- future[ current, ]
else
future <- future[ current, feature ]
}
if( is.null( future ))
return( NA )
if(is.character( future ))
future <- as.factor(future)
return(future)
}
# Fallback for webR standalone adapter objects
if (!is.null(community$pop[[species]])) {
pch_vec <- community$pop[[species]]$pch
col_vec <- community$pop[[species]]$col
fut <- data.frame(pch = pch_vec, color = col_vec, stringsAsFactors = FALSE)
if (missing(feature)) return(fut)
if (length(feature) == 1) return(fut[[feature]])
return(fut[, feature, drop = FALSE])
}
NULL
}
###########################################################################################
get.interact <- function( community, species, host, avail, event )
{
id <- get.species.element( community, host, "stage", avail )
interact <- getOrgInteract( community, host, species, event )[id]
interact[ is.na( interact ) ] <- 0
interact
}
###########################################################################################
getOrgInteract <- function( community,
org1name = getOrgFeature( community, org2name, "substrate" ),
org2name, event = NULL )
{
OrgInteract <- getOrgInfo( community, "Interact" )
if (!is.null(OrgInteract) && !is.null(OrgInteract[[org1name]][[org2name]])) {
tmp <- OrgInteract[[org1name]][[org2name]]
if( is.null( event ))
return( tmp )
event <- as.character( event )
inter <- tmp[,event]
if( length( event ) == 1 )
names( inter ) <- row.names( tmp )
if(is.character(inter))
inter <- factor(inter)
return(inter)
}
# Fallback for webR standalone adapter objects
sub_names <- if (!is.null(community$sub_names)) community$sub_names else c("fr1", "fr2", "fr3", "fr4", "twig", "lftop", "lfbot")
mat <- matrix(1, nrow = length(sub_names), ncol = 1, dimnames = list(sub_names, "substrate"))
as.data.frame(mat)
}
###########################################################################################
getOrgMeanValue <- function( community, species )
{
OrgMeanValue <- getOrgInfo( community, "MeanValue" )
## The global org$MeanValue[[species]] contains mean value information.
OrgMeanValue[[species]]
}
###########################################################################################
copyOrgInfo <- function( fromname, toname )
{
out <- list()
for( i in c("sim")) {
from <- paste( fromname, i, sep = "." )
if( exists( from )) {
toto <- paste( toname, i, sep = "." )
out[[toto]] <- get( from )
cat( "copied", from, "to", toto, "into list\n" )
}
}
return(out)
}
###########################################################################################
get.alive <- function( community, species, substrate )
{
alive <- getOrgAlive( community, species )
alive <- seq_along( alive )[alive]
alive[ substrate == get.species.element( community, species, "sub.stage", alive ) ]
}
###########################################################################################
getOrgAlive <- function( community, species, element )
{
organism <- get.species( community, species )
## identify dead organisms (free nodes for leftist tree)
tmp <- c( FALSE, apply( organism[c("dist","up","left","right"),-1], 2,
function( x ) any( x > 1 )))
if( !any( tmp )) {
tmpp <- organism["up",1]
if( tmpp > 1 )
tmp[tmpp] <- TRUE
}
if( !missing( element ))
tmp <- organism[ element, tmp ]
tmp
}
###########################################################################################
getOrgAgeClass <- function( community, species, stage = seq_len( nrow( future )),
future = getOrgFuture( community, species ))
{
ageclass <- future$ageclass[stage]
tmp <- !is.na( ageclass )
if( any( tmp ))
ageclass[ !is.na( ageclass ) ]
else
NA
}
###########################################################################################
getOrgSubstrate <- function( community, species, elements = seq_len( nrow( inter )),
substrate = getOrgFeature( community, species, "substrate" ),
inter = getOrgInteract( community, substrate, species ))
{
sites <- inter$substrate[elements]
tmp <- !is.na( sites )
if( any( tmp ))
sites[ !is.na( sites ) ]
else
NA
}
###########################################################################################
sampleOrgSubstrate <- function( community, species, elements = seq_len( nrow( inter )),
substrate.name = getOrgFeature( community, species, "substrate" ),
inter = getOrgInteract( community, substrate.name, species ))
{
if( is.na( substrate.name ))
return( elements )
newsub <- as.matrix( cbind( elements, inter[ elements, levels( factor(inter$substrate) ) ] ))
apply( newsub, 1, function( x, is ) {
ns <- sample( levels( factor(is) ), 1, prob = x[-1] / sum( x[-1] ))
sub <- seq_len( nrow( inter ))[ ns == is ]
if( length( sub ) > 1 ) {
newsub <- getOrgInteract( community, substrate.name, substrate.name )[x[1],sub]
sample( sub, 1, prob = newsub / sum( newsub ))
}
else
sub
}, inter$substrate )
}
##########################################################################################
### simulation organism administration
##########################################################################################
initOrgInfo <- function( package, messages = TRUE, datafile = "", ... )
{
community <- list( pop = list( ))
community$org <- list( )
community$org$package <- package
## Get data
community$org$Feature <- getOrgData(community, "organism", "features",
messages, datafile)
community$pop <- list()
community
}
##########################################################################################
setOrgInfo <- function( community, species, hosts, package, messages = TRUE,
datafile = "", ... )
{
Organism <- community$org
if( is.null( Organism$Future ))
Organism$Future <- list( )
if( is.null( Organism$Interact ))
Organism$Interact <- list( )
for( j in hosts )
if( is.null( Organism$Interact[[j]] ))
Organism$Interact[[j]] <- list( )
## Do not reset MeanValue as it may have important spline fits!
if( is.null( Organism$MeanValue ))
Organism$MeanValue <- list( )
for( i in species ) {
future <- getOrgData(community, "future", i,
messages, datafile)
# Check that future agrees with organism.feature information
subclass <- Organism$Feature[i,"subclass"]
if(!(subclass %in% unique(future$ageclass))) {
stop(paste("Future table", paste("future", i, sep = "."),
"does not include", subclass))
}
level.ageclass <- unique( future$ageclass )
level.ageclass <- as.character( level.ageclass[ !is.na( level.ageclass ) ] )
future$ageclass <- ordered( future$ageclass, level.ageclass )
Organism$Future[[i]] <- future
for( j in hosts )
if( i != j ) {
Organism$Interact[[j]][[i]] <- getOrgData(community, j, i,
messages, datafile)
# Check that interaction agrees with host current stage information
# This is messy!
if(j %in% species) {
if(!all(row.names(Organism$Interact[[j]][[i]]) %in%
c(as.character(Organism$Future[[j]]$current), i))) {
stop(paste("Interaction table", paste(j, i, sep = "."),
"does not match", j, "current stages"))
}
}
}
if( is.null( Organism$MeanValue[[i]] ))
Organism$MeanValue[[i]] <- list( )
else
cat( "Keeping Mean Value information for", i, "if any\n" )
}
for( i in unique( getOrgFeature( community, species, "substrate" ))) {
Organism$Interact[[i]][[i]] <- getOrgData(community, i, i,
messages, datafile)
}
community$org <- Organism
community
}
###########################################################################################
getOrgData <- function(community, left, right,
messages = TRUE, datafile = "")
{
# Get Organism Data from
# package data
# global data supplied by user
# external data file supplied by user
sheet <- paste( left, right, sep = "." )
data_exists <- FALSE
if (is.character(datafile) && length(datafile) == 1 && nzchar(datafile)) {
if (dir.exists(datafile)) {
extensions <- c(".txt", ".tsv", ".csv", ".xls", ".xlsx")
datafile_paths <- file.path(datafile, paste0(sheet, extensions))
exist_idx <- file.exists(datafile_paths)
if (any(exist_idx)) {
datafile <- datafile_paths[exist_idx][1]
data_exists <- TRUE
} else {
data_exists <- FALSE
}
sheet <- ""
} else if (file.exists(datafile)) {
data_exists <- TRUE
}
}
if(!data_exists) {
# Load package data or get user-provided global data.
mydata( sheet, getOrgInfo( community, "package" ), messages = messages)
my.eval( sheet )
} else {
# Read data file from user if provided.
if(sheet == "")
my.read(datafile)
else {
out <- as.data.frame(readxl::read_excel(datafile, sheet = sheet, .name_repair = "none"))
if(names(out)[1] == "") { # first column is actual row names
rownames(out) <- out[[1]]
out[[1]] <- NULL
}
out
}
}
}
###########################################################################################
getOrgNames <- function(datafile = "") {
if(datafile == "") {
c("organism.features", "future.host", "future.parasite",
"substrate.host", "substrate.parasite", "substrate.substrate",
"temperature.base", "temperature.par")
} else {
readxl::excel_sheets(datafile)
}
}
###########################################################################################
getOrgDataSimple <- function(community, dataname, datafile = ""){
out <- getOrgData(
community,
left = stringr::str_remove(dataname, "\\..*"),
right = stringr::str_remove(dataname, ".*\\."),
messages = FALSE, datafile = datafile)
# Kludge to reinstate rownames as a column
if(!identical(rownames(out), as.character(seq_len(nrow(out))))) {
out <- data.frame(rownames = rownames(out), out)
}
out
}
###########################################################################################
getOrgInfo <- function( community, element )
{
community$org[[element]]
}
###########################################################################################
setOrgMeanValue <- function( community, species, stage, mvalue )
{
## The global Organism$MeanValue[[species]] contains mean value information.
community$org$MeanValue[[species]][[stage]] <- mvalue
community
}
# --- Source: organism.features.R ---
NULL
# --- Source: future.R ---
## $Id: future.R,v 1.0 2002/12/11 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
###############################################################################
##
## future.events( community )
##
###############################################################################
###############################################################################
### Get birth and future event
###############################################################################
get.future <- function (community, species,
individuals = get.individual(community, species, id),
id = get.base(community, species))
{
## NOTE: This is the slow routine. For every event, it has to check if there
## is a mean value function and then call rspline.
## the structure future.species is set up to handle competing risks!
future <- getOrgFuture(community, species, c("current", "future", "time"))
future$fid <- match(future$future, future$current)
individuals <- as.matrix(individuals)
rownames <- dimnames(individuals)[[1]]
for (i in seq(ncol(individuals))) {
individual <- individuals[,i]
current <- individual["stage"]
individual["location"] <- individual["time"]
## competing risks based on potential future event times
futures <- future$current == future$current[current]
times <- rep( individual["time"], sum( futures ))
cur <- seq(nrow(future))[futures]
for( j in seq( sum( futures ))) {
meantime <- future[cur[j], "time"]
if( meantime > 0 ) {
for.stage <- as.character(future$current[future$fid[ cur[j] ]])
times[j] <- rspline( meantime, individual,
getOrgMeanValue(community, species)[[for.stage]])
}
}
individual["time"] <- min( times )
current <- cur[ times == individual["time"] ][1]
individual["future"] <- future$fid[current]
if (individual["time"] == individual["rejection"])
individual["future"] <- future$fid[future$current == "death"]
individuals[,i] <- individual
}
individuals
}
###############################################################################
event.death <- function( community, species,
id = get.base( community, species ))
{
## remove dead individual from leftist tree
community <- put.species( community, species,
leftist.remove( get.species( community, species ), id ))
## free up individual for reuse
community <- put.base( community, species, id )
community
}
###############################################################################
update_mintime <- function( object, species, ... )
{
base <- get.base( object, species )
mintime <- max( getCount( object, species, "mintime" ),
getTime( object, species,
get.species.element( object, species, "time", base )))
setCount( object, species, list( base = base, mintime = mintime ))
}
###############################################################################
set.birth <- function( community, species, neworg )
{
## merge immediate new births (if any)
newbirths <- ncol( neworg )
if( newbirths > 0 ) {
neworg[c("dist","left","right","up"),] <- 1
community <- checkTime( community, neworg["time",],
getCount( community, species, "mintime" ),
getOrgFeature( community, species, "units" ))
oldbase <- getCount( community, species, "base" )
tmp <- leftist.birth( get.species( community, species ), neworg,
getCount( community, species, "free" ))
community <- put.species( community, species, tmp$tree )
community <- put.base( community, species, free = tmp$free )
community <- updateCounts( community, species, newbirths )
}
community
}
# --- Source: future.host.R ---
NULL
# --- Source: future.meanvalue.R ---
###########################################################################################
## init.meanvalue( organism, stage )
## spline.meanvalue( x, y )
## spline.meanvalue( data = data )
##
## five.show( )
## five.plot( )
##
###########################################################################################
# Curve Designing routines -- under development
###########################################################################################
future.meanvalue <- function( community, species, event = future$current[1],
data )
{
future <- getOrgFuture( community, species )
mvalue <- getOrgMeanValue( community, species )[[event]]
if( missing( data )) {
if( !is.null( mvalue )) {
mvalue <- stats::predict( mvalue$meanvalue, mvalue$meanvalue$knots )
mvalue <- spline.meanvalue( mvalue$x, mvalue$y )$fit
}
else
mvalue <- spline.meanvalue( )$fit
}
else
mvalue <- spline.meanvalue( data = data )$fit
setOrgMeanValue( community, species, event, mvalue )
}
###########################################################################################
spline.meanvalue <- function( x = xinit, y = yinit, data, nspline = 8,
xy = data.frame( x = x, y = y ),
tol = 1e-5, n = 1 )
{
is.data <- !missing( data )
if( !is.data ) {
tmp <- - log( 1 - seq( 0, 1 - exp( -5 ), length = nspline ))
if( missing( x ))
xinit <- tmp
else
xinit <- x
if( missing( y ))
yinit <- tmp
else
yinit <- y
data <- NULL
}
else {
xinit <- sort( data )
ndata <- length( data )
yinit <- seq( ndata ) / ( 1 + ndata )
choose <- round( seq( 1, ndata, length = nspline ))
xinit <- xinit[choose]
yinit <- - log( 1 - yinit[choose] )
}
fs <- list( probability = function( x ) { - log( 1 - x ) } )
finvs <- list( probability = function( x ) { 1 - exp( - x ) } )
for( i in c("mean value","rate","density") )
fs[[i]] <- finvs[[i]] <- function( x ) x
## plot curve and surrounding axes
graphics::par( mfrow = c(1,1), mar = rep(4.1,4))
plotit <- function( xy, fig = "mean value", fit = splines::interpSpline( xy$x, xy$y ))
{
switch( fig, {
y <- xy$y
ylim <- range(c(0,y))
},
probability = {
y <- 1 - exp( - xy$y )
ylim <- range(c(0,y))
},
rate = {
y <- spline.rate( fit, xy$x )$y
tmp <- spline.rate( fit )
ylim <- range(c(0,tmp$y))
},
density = {
y <- spline.rate( fit, xy$x )$y * exp( - stats::predict( fit, xy$x )$y )
tmp <- spline.rate( fit )
tmp$y <- tmp$y * exp( - stats::predict( fit )$y )
ylim <- range(c(0,tmp$y))
}
)
plot(xy$x,y,xlim=1.25*range(c(0,xy$x)), ylim = ylim,
type="n", xlab = "", ylab = "" )
graphics::points( xy$x, y, lwd = 4 )
graphics::title( fig )
graphics::mtext( "time", 1, 2 )
graphics::mtext( fig, 2, 2 )
graphics::abline( v = max( xy$x ), lty = 2 )
switch( fig, {
if( fig == "probability" ) {
tmp <- c(.1,.2,.5,1:10)
ltmp <- 1-exp(-tmp)
graphics::mtext( "mean value", 4, 2 )
}
else {
tmp <- c(seq(0,.9,,by=.1),.95,.98,.99,.999)
ltmp <- -log(1-tmp)
graphics::mtext( "probability", 4, 2 )
}
usr <- graphics::par("usr")
s <- ltmp <= usr[4] & ltmp >= usr[3]
tmpar <- graphics::par( cex = .75 )
graphics::axis(4,ltmp[s],tmp[s])
graphics::par( tmpar )
summaryshow( xy, fit, "white" )
tmp <- curve.plot( xy, n = n, action = "refresh", fit = fit,
f = fs[[fig]], finv = finvs[[fig]] )
summaryshow( xy, tmp$fit )
},
rate =, density = {
graphics::lines( tmp$x, tmp$y )
summaryshow( xy, fit )
}
)
}
## place commands along right strip of plot, highlighting current command
plotcmd <- function( ans, fig, cmds, cmdlocs, usr, col = "green", rest = "black",
data = FALSE )
{
ans <- c( ans, fig )
if( data )
ans <- c( ans, "data" )
tmp <- is.na( match( cmds, ans ))
if( any( tmp ))
graphics::text( rep(usr[2],sum(tmp)), cmdlocs[tmp], cmds[tmp], col = rest, adj = 1 )
if( any( !tmp ))
graphics::text( rep(usr[2],sum(!tmp)), cmdlocs[!tmp], cmds[!tmp], col = col, adj = 1 )
}
fig <- "mean value"
newans <- ans <- "replace"
graphics::par( mar = c(4.1,4.1,3.1,4.1),omi=rep(.25,4))
fit <- splines::interpSpline( xy$x, xy$y )
sums <- plotit( xy, fig, fit )
cmds <- c("refresh","add","delete","replace","rescale","shrink to 1","finish","restart",
"","data","mean value","probability","rate","density")
newlocs <- function( cmds, data = FALSE, usr )
{
if( !data )
cmds <- cmds[ cmds != "data" ]
n <- length( cmds )
blank <- seq( n )[ cmds == "" ]
tmp <- diff(usr[3:4]) / 20
tmp <- c( usr[4] - tmp * seq( blank - 1 ), mean( usr[3:4] ),
usr[3] + tmp * seq( n - blank ))
names( tmp ) <- cmds
tmp
}
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
cmds <- names( cmdlocs )
use.data <- FALSE
plotcmd( ans, fig, cmds, cmdlocs, usr )
rescale.data <- 1
repeat {
## get command from plot using cursor
z <- graphics::locator(1,"n")
if( z$x > max( xy$x )) {
z <- abs(z$y - cmdlocs )
newans <- cmds[z==min(z)][1]
switch( newans,
finish =, refresh = {
sums <- plotit( xy, fig, fit )
},
data = {
use.data <- is.data & !use.data
if( is.data & !use.data & match( fig, c("mean value","probability"),
nomatch = 0 ))
sums <- plotit( xy, fig, fit )
},
"mean value" =, probability =, rate =, density = {
fig <- newans
sums <- plotit( xy, fig, fit )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
"shrink to 1" = {
while( abs( sums[1] - 1 ) > tol ) {
xy$y <- xy$y * sums[1]
if( is.data )
rescale.data <- rescale.data * sums[1]
fit <- splines::interpSpline( xy$x, xy$y )
sums <- splinesum( xy, fit, tol )
}
sums <- plotit( xy, fig, fit )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
rescale = {
cat( "enter new values followed by RETURN key\n" )
tmpy <- readline( paste( "maximum mean value(",
round( max( xy$y ), 2 ), "):", sep = "" ))
if( tmpy != "" ) {
tmpy <- as.numeric( tmpy ) / max( xy$y )
xy$y <- tmpy * xy$y
if( is.data )
rescale.data <- rescale.data * tmpy
}
tmpx <- readline( paste( "maximum time(",
round( max( xy$x ), 2 ), "):", sep = "" ))
if( tmpx != "" ) {
tmpx <- as.numeric( tmpx ) / max( xy$x )
xy$x <- tmpx * xy$x
if( is.data )
data <- data * tmpx
}
fit <- splines::interpSpline( xy$x, xy$y )
sums <- plotit( xy, fig, fit )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
restart = {
if( is.data )
rescale.data <- 1
xy <- data.frame( x = xinit, y = yinit )
fit <- splines::interpSpline( xy$x, xy$y )
sums <- plotit( xy, fig, fit )
usr <- graphics::par("usr")
cmdlocs <- newlocs( cmds, data = is.data, usr = usr)
},
add =, delete =, replace = {
ans <- newans
}
)
if( use.data & match( fig, c("mean value","probability"), nomatch = 0 ))
cdf.lines( data, fig, rescale = rescale.data )
plotcmd( ans, fig, cmds, cmdlocs, usr, data = use.data )
}
else {
if( is.na( match( fig, c("rate","density") ))) {
summaryshow( xy, fit, "white", sums )
fit <- curve.plot( xy, n = n, action = ans, z = z, fit = fit,
f = fs[[fig]], finv = finvs[[fig]] )
xy <- fit$xy
fit <- fit$fit
sums <- summaryshow( xy, fit )
ans <- "replace"
plotcmd( ans, fig, cmds, cmdlocs, usr, data = use.data )
}
}
if( newans == "finish" )
break
}
plotcmd( newans, fig, cmds, cmdlocs, "red", data = use.data )
summaryshow( xy, fit, "white", sums )
tmp <- curve.plot( xy, n = n, action = "refresh",
f = fs[[fig]], finv = finvs[[fig]] )
tmp$meanvalue <- tmp$fit
tmp$fit <- NULL
tmp$invmvalue <- splines::backSpline( tmp$meanvalue )
sums <- summaryshow( xy, tmp$meanvalue )
tmp$mean <- sums[1]
tmp$median <- sums[2]
tmp
}
# --- Source: substrate.host.R ---
NULL
# --- Source: host.parasite.R ---
NULL
# --- Source: simdata.R ---
NULL
# --- Source: redscale.R ---
NULL
# --- Source: init.simulation.R ---
init.simulation <- function( package = "ewing",
count = 200,
interact = FALSE,
messages = TRUE,
...)
{
community <- initOrgInfo( package, messages = messages, ... )
community <- initTemp( community, messages = messages, ... )
species <- getOrgFeature( community )[1:2]
hosts <- getOrgHosts( community, species )
if(messages) {
cat( "Creating simulation organism set using species:\n",
paste( species, collapse = ", " ), "\n\n" )
}
community <- setOrgInfo( community, species, hosts, package,
messages = messages, ... )
if(messages) {
cat( "\n" )
}
species <- unique( species )
num <- numeric( length( species ))
names( num ) <- species
count <- rep_len(count, length(species))
names(count) <- species
for( i in species ) {
num[i] <- reuse <- count[i]
if(messages | (interact & interactive())) {
cat( paste( "Initialize ", i, " at size ", reuse, sep = "" ))
}
if(interact & interactive()) {
cat(" :")
r <- readline( )
if( r != "" & is.na( pmatch( substring( r, 1, 1 ), c("y","Y") )))
reuse <- suppressWarnings(as.numeric( r ))
if( is.na( reuse ))
reuse <- count[i]
}
if( reuse ) {
if(messages) {
cat( "...\n" )
}
community <- init.population( community, i, n = reuse, messages = messages, ... )
num[i] <- reuse
}
}
class( community ) <- c("ewing", "list")
attr(community, "count") <- count
community
}
# --- Source: init.population.R ---
init.population <- function( community, species, n = 200, width = 100,
units = getOrgFeature( community, species, "units" ),
timeit = FALSE,
reject = Inf,
position = rtri( n, width ),
colnames = c(leftistnames,paramnames,posnames,eventnames),
init.stage = istage,
init.weight = getOrgFuture( community, species, "init" ),
messages = TRUE,
...)
{
leftistnames <- c("dist","left","right","up")
paramnames <- c("dispersion","location","intensity","truncation","rejection")
posnames <- paste("pos",letters[1:3], sep = ".")
eventnames <- c("time","stage","future","offspring","sex","sub.stage","sub.future")
organism <- matrix( 0, length( colnames ), n+1,
dimnames = list( colnames, NULL ))
organism["time",1] <- Inf
## 5-parameter initialization
organism[c("dispersion","intensity"),-1] <- 1
organism["rejection",-1] <- if( reject == Inf )
rep( Inf, n )
else
reject * stats::rexp( n )
## triangular coordinates
organism[posnames,-1] <- position
## substrate
substrate.name <- getOrgFeature( community, species, "substrate" )
if( !is.na( substrate.name )) {
substrate <- getOrgInteract( community, substrate.name, species, "init" )
organism["sub.stage",-1] <- organism["sub.future",-1] <- sample( length( substrate ),
n, replace = TRUE, prob = substrate / sum( substrate ))
}
## randomly generate events proportional to future time units
nstage <- length( init.weight )
istage <- sample( nstage, n, replace = TRUE, prob = init.weight / sum( init.weight ))
init.stage <- array( init.stage, n )
organism["stage",-1] <- init.stage
## schedule future events
if( timeit )
proctime <- proc.time()
organism[,-1] <- get.future( community, species, organism[,-1] )
if( timeit ) {
tmp <- proc.time() - proctime
cat( "future time: user=", tmp[1], "system=", tmp[2], "total=", tmp[3], "\n" )
}
## create leftist tree
if( timeit )
proctime <- proc.time()
community <- put.species( community, species, leftist.create( organism ))
if( timeit ) {
tmp <- proc.time() - proctime
cat( "leftist time: user=", tmp[1], "system=", tmp[2], "total=", tmp[3], "\n" )
}
## mean number of offspring
if( timeit )
proctime <- proc.time()
#***This is where parasite is crashing--no offspring?**
community <- initOffspring( community, species )
if( timeit ) {
tmp <- proc.time() - proctime
cat( "offspring time: user=", tmp[1], "system=", tmp[2], "total=", tmp[3], "\n" )
}
if(messages) {
cat( "Initializing events for", species, "with", ncol( organism ) - 1, "individuals\n" )
}
community
}
# --- Source: init.timing.R ---
###########################################################################################
### Timing of simulation run
###########################################################################################
init.timing <- function( community )
{
## initialize timing
events <- NULL
for( species in get.species( community ))
events <- c( events, levels( getOrgFuture( community, species, "event" )))
events <- sort( unique( events ))
tmp <- c("total",events,"refresh","other")
cpu <- matrix( 0, 3, length( tmp ),
dimnames = list( c("user","system","total"), tmp ))
community$cpu <- cpu
community <- set.timing( community, "total" )
community
}
###########################################################################################
set.timing <- function( community, string, flag = -1 ) {
if( !is.null( community$cpu ))
community$cpu[,string] <- community$cpu[,string] + flag * proc.time()[1:3]
community
}
###########################################################################################
fini.timing <- function( community )
{
if( !is.null( community$cpu )) {
community <- set.timing( community, "total", 1 )
community$cpu[,"other"] <- community$cpu[,"total"] - apply( community$cpu[,-1], 1, sum )
}
community
}
# --- Source: initCount.R ---
###########################################################################################
### Simulation count object administration
###########################################################################################
initCount <- function( community, species, debugit = FALSE, file = NULL, append = FALSE,
messages = TRUE, ... )
{
if(messages) {
cat( "initial" )
for( i in species)
cat( ":", i, sum( apply( get.species( community, i ), 2,
function(x) !all(x[c("dist","left","right","up")]==1))) - 1 )
cat( "\n" )
}
old_counts <- if (append) getCount( community, , "counts" ) else NULL
old_step <- if (append) getCount( community, , "step" ) else 0
count <- list()
## leftist tree counters
count$mintime <- numeric( length( species ))
count$base <- numeric( length( species ))
names( count$base ) <- names( count$mintime ) <- species
count$free <- list()
## initialize lists to keep track of events
count$events <- count$countage <- count$countsub <- count$nameage <- count$namesub <- list()
## Set up hour to degreeday spline based on range of hours if any
simmin <- c(hr=Inf,DD=Inf)
for( i in species ) {
if( is.null( get.species( community, i )))
stop( paste( "Missing species", i ))
count$free[[i]] <- 1
count$base[i] <- get.base( community, i )
units <- getOrgFeature( community, i, "units" )
tmp <- get.individual( community, i )["time"]
if( tmp < simmin[units] )
simmin[units]<- tmp
}
if( max( simmin ) < Inf ) {
community <- activeTemp( community, simmin["hr"], , simmin["DD"], messages = messages )
}
esums <- c("initial","during","final")
tmpfn <- function( counter )
{
rownames <- levels( counter )
array( 0, length( rownames ), dimnames = list( rownames ))
}
subclass <- getOrgFeature( community, species, "subclass" )
names( subclass ) <- species
for( i in species ) {
species.time <- get.individual( community, i )["time"]
count$mintime[i] <- getTime( community, i, species.time )
future <- getOrgFuture( community, i )
## possible future events
count$events[[i]] <- matrix( 0, nrow( future ), length( esums ),
dimnames = list( as.character( future$current ), esums ))
## current record of future events
count$countage[[i]] <- tmpfn( getOrgFuture( community, i, "ageclass" ))
count$countsub[[i]] <- tmpfn( getOrgInteract( community,, i, "substrate" ))
if( species.time < Inf ) {
## count by age groups
stage <- getOrgAlive( community, i, "stage" )
if( length( stage )) {
classes <- getOrgAgeClass( community, i, stage )
tmp <- table( classes )
count$countage[[i]][ names( tmp ) ] <- tmp
}
## count by substrate
substage <- getOrgAlive( community, i, "sub.stage" )
## only for individuals of class = subclass[i]
substage <- substage[ subclass[i] == getOrgAgeClass( community, i, stage ) ]
if( length( substage )) {
classes <- getOrgInteract( community,, i, "substrate" )
tmp <- table( classes[substage] )
count$countsub[[i]][ names( tmp ) ] <- tmp
}
}
}
count$debug <- debugit
# If file is NULL, then don't write to file; keep counts internal
count$file <- file
if (append && !is.null(old_counts)) {
count$counts <- old_counts
count$step <- old_step
}
community$count <- count
## Put counts in file
community <- putCount( community, append )
## tally events at start of simulation
setEvents( community, "initial" )
}
###########################################################################################
getCount <- function( community, species, element )
{
count <- community$count[[element]]
if( !missing( species ))
count <- count[[species]]
count
}
###########################################################################################
set.step <- function( community, step )
setCount( community,, list( step = step ))
###########################################################################################
setCount <- function( community, species, elements )
{
count <- community$count
for( i in names( elements )) {
if( missing( species ))
count[[i]] <- elements[[i]]
else
count[[i]][[species]] <- elements[[i]]
}
community$count <- count
community
}
# --- Source: future.events.R ---
future.events <- function( community,
nstep = 4000,
species = get.species( community ),
refresh = nstep / 20, cex = 0.5,
substrate.plot = TRUE, extinct = TRUE,
timeit = TRUE, debugit = FALSE,
messages = TRUE, append = NULL, ... )
{
## Integrity check of dataset, and initialization of tallies.
if( missing( community ))
stop( "Must specify a community." )
if (is.null(append)) {
append <- !is.null(community$count$counts) && nrow(community$count$counts) > 0
}
if( debugit ) cat( "initialization\n" )
community <- initCount( community, species, debugit = debugit, file = NULL,
append = append, messages = messages )
if( timeit )
community <- init.timing( community )
mintime <- getCount( community, , "mintime" )
species.now <- species[ mintime == min( mintime ) ][1]
future <- getOrgFuture( community, species.now )
# Set up list for plot information.
p <- list()
pstep <- 0
if (!is.null(community$plot)) {
p <- community$plot
pstep <- length(p)
}
start_step <- if (append && !is.null(community$count$step)) community$count$step else 0
## for nstep steps schedule future events and process immediate events
for( istep in seq( nstep )) {
## stop if any extinct and extinct flag on, or all extinct
omintime <- mintime
mintime <- getCount( community, , "mintime" )
if( debugit ) print( mintime )
if( min( mintime ) < min( omintime )) {
cat( "time reversal!\n" ) # should not happen
browser()
}
tmp <- mintime == Inf
if( any( tmp )) {
if( extinct | all( tmp )) {
for( i in names( mintime )[tmp] )
cat( "***", i, "is extinct ***\n" )
if( plotit )
plot.ewing( community, substrate = substrate.plot, cex = cex, ...)
break
}
}
## each species is always sorted so 1st element is next future event
species.prev <- species.now
species.now <- species[ mintime == min( mintime ) ][1]
individual <- get.individual( community, species.now )
if( is.na( individual["time"] ) | individual["time"] == Inf ) {
cat( "No more finite future events. End of simulation.\n" )
break
}
if( all( individual[c("dist","left","right","up")] == 1 )) {
cat( individual["time"], ": last", species.now, "alive",
getCount( community, species.now, "base" ), "\n" )
}
if( species.now != species.prev )
future <- getOrgFuture( community, species.now )
## make future event the current stage
current <- individual["stage"]
stage <- as.character( future$current[current] )
if(!length(stage))
stop(paste("no stage", current))
community <- updateCount( community, species.now, individual, stage == "death",
start_step + istep )
individual["stage"] <- current
individual["sub.stage"] <- individual["sub.future"]
community <- put.individual( community, species.now, individual )
if( debugit ) {
cat( species.now, istep, "base",
getCount( community, species.now, "base" ), "\n" )
print( c( step = istep,
countage = sum( getCount( community, species.now, "countage" )),
countsub = sum( getCount( community, species.now, "countsub" )),
round( individual["time"], 2 ), current, stage ))
}
## processing of immediate, pending and future events
## this is the main show--all the rest is setup
event.type <- as.character( future[ current, "event" ] )
if( debugit ) cat( "do", event.type, istep, stage,
as.character( future$current[ individual["future"] ] ),
individual["time"], "\n" )
community <- set.timing( community, event.type )
if( event.type == "death" )
community <- event.death( community, species.now )
else {
if( event.type != "future" ) {
## this routine could be user supplied
## generic routines are event.birth, event.attack
event.parsed <- get( paste( "event", event.type, sep = "." ))
community <- event.parsed( community, species.now )
}
community <- event.future( community, species.now )
}
community <- set.timing( community, event.type, 1 )
## refresh plot
if( refresh & ! ( istep %% refresh )) {
community <- set.timing( community, "refresh" )
# Save the ewing_ageclass and ewing_substrate objects.
pstep <- pstep + 1
p[[pstep]] <- ewing_snapshot(community, start_step + istep, ...)
if(messages) {
cat( "refresh", istep )
for( j in get.species( community ))
cat( ":", j, sum( getCount( community, j, "countage" ), na.rm = TRUE ))
cat( "\n" )
}
community <- set.timing( community, "refresh", 1 )
}
## periodic browser if in debug mode
if( debugit ) {
cat( "done", istep, "\n" )
if( refresh & !( istep %% refresh )) {
cat( "Type \"c\" to continue or \"Q\" to quit.\n" )
browser()
}
}
}
community <- set.timing( community, "refresh" )
## end of main loop on future events
if( debugit ) cat("done\n")
if( sum( getOrgAlive( community, species.now )) > 1 ) {
## tally events at end of simulation
community <- setEvents( community, "final" )
}
community <- set.timing( community, "refresh", 1 )
community <- fini.timing( community )
community$step <- start_step + nstep
attr(community, "nstep") <- start_step + nstep
if( !refresh | (nstep%%refresh)) {
pstep <- pstep + 1
p[[pstep]] <- ewing_snapshot(community, start_step + istep, ...)
}
community$plot <- p
community
}
# --- Source: event.R ---
## $Id: event.R,v 1.0 2002/12/11 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
###############################################################################
##
## event.birth( community, species )
## event.attack( community, species )
##
## Issues to resolve:
## 2. search strategy for predator/parasite/parasitoid
## 3. generic calls to do.x
###############################################################################
event.birth <- function( community, species )
{
## now only one offspring at a time,
## but could depend on individual
offspring <- get.offspring( community, species)
## update parent based on anticipated offspring
community <- parent.birth( community, species, offspring )
if( offspring ) {
## get new births
newbirths <- get.birth( community, species, offspring )
## merge births into community
community <- set.birth( community, species, newbirths )
}
community
}
###############################################################################
parent.birth <- function( community, species, offspring )
{
## get individual record
individual <- get.individual( community, species )
individual["offspring"] <- individual["offspring"] - offspring
## female starves when egg load depleted (less than or equal to 0)
if( individual["offspring"] <= 0 )
individual["stage"] <- set.future( community, species, "starved" )
community <- put.individual( community, species, individual )
}
###############################################################################
get.birth <- function( community, species, offspring )
{
## get individual record
individual <- get.individual( community, species )
## matrix of new offspring for community
newbirths <- matrix( individual, length( individual ), offspring )
dimnames( newbirths ) <- list( names( individual ), NULL )
if( offspring ) {
## assumes newborn is stage 1, and next stage is 2
newbirths["stage",] <- 1
newbirths["future",] <- 2
## set up as unlinked node for leftist tree
newbirths[c("dist","left","right","up"),] <- 1
## disperse offspring across substrate types
newbirths <- event.move( community, species, newbirths )
}
if( getCount( community,, "debug" ))
cat( round( individual["time"] ),
getOrgFeature( community, species, "units" ),
":", species, "offspring", individual["offspring"], "\n" )
if( offspring ) {
## get future events for new organisms
get.future( community, species, newbirths)
}
else
community
}
###############################################################################
get.deplete <- function( community, species )
{
## get individual record of attacker
individual <- get.individual( community, species )
## Deplete reserves based on time spent searching for host.
individual["offspring"] <- individual["offspring"] -
( individual["time"] - individual["location"] ) /
getOrgFeature( community, species, "deplete" )
community <- put.individual( community, species, individual )
}
###############################################################################
set.future <- function( community, species, stage )
{
current <- getOrgFuture( community, species, "current" )
seq( length( current ))[ current == stage ]
}
###############################################################################
### Interaction Events (only attack for now)
###############################################################################
event.attack <- function( community, species )
{
## dyadic event: attack of host by adult parasitoid
## deplete individual based on time spent searching for host
community <- get.deplete( community, species )
## get individual record of attacker
individual <- get.individual( community, species )
## find a host if parasite has offspring reserve left
if( individual["offspring"] > 0 ) {
## get name of host for attacker
host <- getOrgFeature( community, species, "attack" )
## get attack parasite and event types
attack <- get.attack( community, species, individual )
## find host located on the same substrate
found <- event.find( community, species, host, attack["event"] )
if( length( found )) {
## host-parasite interaction
event.parsed <- get( paste( "event", attack["event"], sep="." ))
community <- event.parsed( community, species, host, found )
event.parsed <- get( paste( "host", attack["parasite"], sep="." ))
community <- event.parsed( community, species, host, found )
individual <- get.individual( community, species )
}
}
if( individual["offspring"] > 0 ) {
## parasite moves along substrate
individual <- event.move( community, species, individual )
}
else {
## parasite dies if it does not feed enough
individual["stage"] <- set.future( community, species, "starved" )
}
## put updated individual back in community
community <- put.individual( community, species, individual )
}
###############################################################################
get.attack <- function( community, species, individual )
{
## get parasite type ("ecto" or "endo") and current event ("feed" or "ovip")
parasite <- getOrgFeature( community, species, "parasite" )
event <- getOrgFuture( community, species, "current", individual["future"] )
event <- as.character( event )
if( parasite=="endo" )
event <- "ovip"
else if( individual["offspring"] < 1 ) {
## must feed if depleted
event <- "feed"
}
c( event = event, parasite = parasite )
}
###############################################################################
host.ecto <- function( community, species, host, dead )
{
## ectoparasites effectively kill their host
## get individual doing the attack
individual <- get.individual( community, species )
## get host individual that is attacked
hostindiv <- get.individual( community, host, dead )
## set host time to now, which may involve hr-DD translation
hostindiv["time"] <- transTime( community, species, host, individual["time"] )
## schedule immediate death of host
hostindiv["future"] <- set.future( community, host, "death" )
## update host record in community
community <- put.individual( community, host, hostindiv, dead )
## update leftist tree and mintime
community <- put.species( community, host,
leftist.update( get.species( community, host ), dead ))
update_mintime( community, host )
}
###############################################################################
host.endo <- function( community, species, host, dead, harm )
{
## endoparasites reduces capacity of host (assumed by half here)
## get host individual that is attacked
hostindiv <- get.individual( community, host, harm )
## schedule harm for hosts (reduce egg capacity by half)
hostindiv["offspring"] <- floor( hostindiv["offspring"] / 2 )
if( hostindiv["offspring"] == 0 ) {
## get individual doing the attack
individual <- get.individual( community, species )
## set host time to now, which may involve hr-DD translation
hostindiv["time"] <- transTime( community, species, host, individual["time"] )
## schedule immediate death of host
hostindiv["future"] <- set.future( community, host, "death" )
}
## update host record in community
community <- put.individual( community, host, hostindiv, harm )
if( hostindiv["offspring"] == 0 ) {
## update leftist tree and mintime
community <- put.species( community, host,
leftist.update( get.species( community, host ), dead ))
community <- update_mintime( community, host )
}
community
}
###############################################################################
event.feed <- function( community, species, host, dead )
{
## feed: adult parasite feeds on host
interact <- get.interact( community, species, host, dead, "feed" )
## host-parasite interaction: feeding
if( interact ) {
individual <- get.individual( community, species )
individual["offspring"] <- individual["offspring"] + interact
community <- put.individual( community, species, individual )
}
community
}
###############################################################################
event.ovip <- function( community, species, host, dead, gender=TRUE )
{
## ovip: adult lays egg in the host to emerge later as adult
offspring <- get.offspring( community, species )
interact <- get.interact( community, species, host, dead, "ovip" )
## host-parasite interaction: feeding
if( interact ) {
## update parent individual, depleting energy after egg laying
community <- parent.birth( community, species, offspring )
## get new births
newbirths <- get.birth( community, species, offspring )
## gender preference for offspring
if( get.interact( community, species, host, dead, "male" ) < stats::runif( 1 ) |
!gender ){
## set offspring for female eggs based on dead host
newbirths["offspring", ] <- set.offspring( community, species, host, dead )
}
else {
## produce male and put on queue for immediate death
newbirths <- set.male( community, species, host, newbirths )
}
community <- set.birth( community, species, newbirths )
}
community
}
###############################################################################
set.male <- function( community, species, host, newbirths )
{
## produce male offspring, which is queued for immediate death
newbirths["future",1] <- set.future( community, species, "male" )
newbirths["time",1] <- get.individual( community, species )["time"]
newbirths
}
# --- Source: event.future.R ---
event.future <- function( community, species )
{
## schedule future event based on current stage
individual <- get.future( community, species )
## move if appropriate
individual <- event.move( community, species, individual )
## update time translation if needed
community <- checkTime( community, individual["time"],
getCount( community, species, "mintime"),
getOrgFeature( community, species, "units" ))
## put updated individual back in community
community <- put.individual( community, species, individual )
## reprioritize the leftist tree if time has changed
if( individual["time"] > individual["location"] ) {
## and time is longer than next scheduled time
if( individual["time"] >
min( get.species.element( community, species, "time",
individual[c("left","right")] ))) {
## remove individual from leftist tree
community <- put.species( community, species,
leftist.update( get.species( community, species )))
}
community <- update_mintime( community, species )
}
community
}
# --- Source: Events.R ---
###########################################################################################
updateEvents <- function( community, species, event, increment = 1 )
{
community$count$events[[species]][event,"during"] <-
community$count$events[[species]][event,"during"] + increment
community
}
###########################################################################################
setEvents <- function( community, period )
{
count <- community$count
for( species in get.species( community )) {
current <- getOrgFuture( community, species, "current" )
events <- rep( 0, length( current ))
names( events ) <- as.character( current )
stage <- getOrgAlive( community, species, "stage" )
if( length( stage )) {
tmp <- tapply( stage, current[stage], length )
tmp[ is.na( tmp ) ] <- 0
events[ names( tmp ) ] <- tmp
}
count$events[[species]][,period] <- events
}
community$count <- count
community
}
# --- Source: leftist.R ---
## $Id: leftist.R,v 1.0 2002/12/11 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
###########################################################################################
## This version has list structure for individuals.
## In addition it will keep species separate.
###########################################################################################
## Triply linked leftist trees
###########################################################################################
leftist.tree <- function( x,
tree = data.frame( root = empty ),
key = "time",
empty = c( time = NA, dist = 1, left = 1, right = 1, up = 1 ))
{
node <- empty
node[key] <- Inf
node[c("right","left")] <- 2
node["dist"] <- 0
if( missing( tree ))
names( tree ) <- NULL
node <- empty
node[key] <- x[1]
base <- 2
tree[[base]] <- node
for( i in 1 + seq( 2, length( x ))) {
node <- empty
node[key] <- x[i-1]
if( node[key] <= tree[[base]][key] ) {
## insert at root
node["left"] <- base
tree[[base]]["up"] <- i
base <- i
tree[[base]] <- node
}
else {
tree[[i]] <- node
## insert in tree (e.g. merge two trees)
tree <- as.data.frame( leftist.merge( as.matrix( tree ), i, base ))
base <- tree["up",1]
}
}
tree[[1]][c("right","left","up")] <- base
tree
}
###########################################################################################
leftist.create <- function( tree, key = "time" )
{
## assume 1st element of tree is for base
## and keys are in place already
base <- 2
tree[c("right","left","up"),] <- 1
tree["dist",-1] <- 1
ntree <- ncol( tree )
for( i in seq( 3, ntree )) {
if( tree[key,i] <= tree[key,base] ) {
## insert at root
tree["left",i] <- base
tree["up",base] <- i
base <- i
}
else {
## insert in tree (e.g. merge two trees)
tree <- leftist.merge( tree, i, base )
base <- tree["up",1]
}
}
## set 1st element to point to base
tree[c("right","left","up"),1] <- base
tree
}
###########################################################################################
leftist.merge <- function( tree, P = 1, Q = 1, R = 1, key = "time" )
{
while( TRUE ) {
if( P == 1 ) {
P <- Q
Q <- 1
}
if( Q == 1 ) {
D <- tree["dist",P]
while( R > 1 ) {
Q <- tree["right",R]
temp <- tree["left",R]
if( tree["dist",temp] < D ) {
D <- tree["dist",temp] + 1
tree["right",R] <- tree["left",R]
tree["left",R] <- P
}
else {
D <- D + 1
tree["right",R] <- P
}
tree["up",P] <- R
tree["dist",R] <- D
P <- R
R <- Q
}
tree["up",P] <- 1
tree[c("left","right","up"),1] <- P
return( tree )
}
## merge two right lists
if( tree[key,P] <= tree[key,Q] ) {
temp <- tree["right",P]
tree["right",P] <- R
tree["up",R] <- P
R <- P
P <- temp
}
else {
temp <- tree["right",Q]
tree["right",Q] <- R
tree["up",R] <- Q
R <- Q
Q <- temp
}
}
tree["up",1] <- P
tree[c("left","right","up"),P] <- 1
tree
}
###########################################################################################
leftist.remove <- function( tree, P )
{
oldbase <- tree["up",1]
up <- tree["up",P]
## merge the subtrees below P
tree <- leftist.merge( tree, tree["left",P], tree["right",P] )
## reset node to empty
tree[c("up","left","right","dist"),P] <- 1
## return if base node removed
if( oldbase == P )
return( tree )
## make up node leftist
if( tree["left",up] == P )
tree["left",up] <- tree["right",up]
tree["right",up] <- 1
tree["dist",up] <- 1
## traverse back up the tree to make sure it is leftist to base
upup <- tree["up",up]
left <- tree["left",upup]
right <- tree["right",upup]
while( tree["dist",left] < tree["dist",right] ) {
tree["dist",upup] <- tree["dist",left] + 1
tree["left",upup] <- right
tree["right",upup] <- left
upup <- tree["up",upup]
left <- tree["left",upup]
right <- tree["right",upup]
}
## merge down and up trees
leftist.merge( tree, oldbase, tree["up",1] )
}
###########################################################################################
leftist.birth <- function( organism, neworg, free )
{
newbase <- organism["up",1]
norganism <- ncol( neworg )
j <- norganism
nfree <- free[1]
while( j > 0 & nfree > 1 ) {
newbabe <- free[nfree]
organism[,newbabe] <- neworg[,j]
organism <- leftist.merge( organism, newbase, newbabe )
newbase <- organism["up",1]
nfree <- nfree - 1
free[1] <- nfree
j <- j - 1
}
if( j > 0 ) {
norganism <- ncol( organism )
js <- 1:j
organism <- cbind( organism, neworg[,js] )
for( i in norganism + js ) {
organism <- leftist.merge( organism, newbase, i )
newbase <- organism["up",1]
}
}
list( tree = organism, base = newbase, free = free )
}
###########################################################################################
leftist.free <- function( free, id )
{
free[1] <- tmp <- free[1] + 1
free[tmp] <- id
free
}
###########################################################################################
leftist.update <- function( tree, P = tree["up",1] )
{
tree <- leftist.remove( tree, P )
leftist.merge( tree, tree["up",1], P )
}
# --- Source: community.R ---
get.organisms <- function(datafile = "") {
org <- list(species = c("host", "parasite"), substrates = "substrate")
if(datafile != "") {
if(tools::file_ext(datafile) %in% c("xls","xlsx")){
sheets <- readxl::excel_sheets(datafile)
species <- stringr::str_remove(
sheets[stringr::str_detect(sheets, "future\\.")],
"future\\.")
substrates <- unique(stringr::str_remove(
sheets[stringr::str_detect(sheets, paste("\\.", species, sep = "", collapse = "|")) &
!stringr::str_detect(sheets, paste(c("future", species), "\\.", sep = "", collapse = "|"))],
"\\..*"))
org <- list(species = species, substrates = substrates)
}
}
org
}
###########################################################################################
get.species <- function( community, species ) {
if( missing( species ))
return( names( community$pop ))
if( is.numeric( species ))
species <- names( community$pop )[species]
if( is.null( species ) || !species %in% names( community$pop ))
return( NULL )
ans <- community$pop[[species]]
if (is.list(ans) && !is.matrix(ans) && !is.data.frame(ans) && !is.null(ans$org)) {
ans <- ans$org
}
if (is.matrix(ans) || is.data.frame(ans)) {
if (!"up" %in% rownames(ans) && (is.null(colnames(ans)) || colnames(ans)[1] != "dummy")) {
dummy <- ans[, 1, drop = FALSE]
colnames(dummy) <- "dummy"
return(cbind(dummy, ans))
}
}
ans
}
###########################################################################################
get.species.element <- function( community, species, rows, cols )
community$pop[[species]][rows,cols]
###########################################################################################
put.species <- function( community, species, value )
{
community$pop[[species]] <- value
community
}
###############################################################################
put.individual <- function( community, species, individual,
id = get.base( community, species ))
{
community$pop[[species]][,id] <- individual
community
}
###############################################################################
get.individual <- function( community, species,
id = get.base( community, species ))
community$pop[[species]][,id]
###############################################################################
get.base <- function( community, species )
community$pop[[species]]["up",1]
# --- Source: c.ewing.R ---
c.ewing <- function(...) {
communities <- list(...)
if(length(communities) < 2) {
if(length(communities) == 1)
return(communities[[1]])
return(NULL)
}
community <- communities[[1]]
# *** need to verify that all communities have save structure
# Elements org, temp should be identical
for(comi in seq(2, length(communities))) {
# Element pop
for(species in names(community$pop)) {
community$pop[[species]] <-
cbind(community$pop[[species]], communities[[comi]]$pop[[species]])
}
# Element cpu
community$cpu <- community$cpu + communities[[comi]]$cpu
}
# Element plot is more complicated as it contains items for plots
# Probably want some form of appendX functions
}
# --- Source: move.R ---
## $Id: move.R,v 1.0 2002/12/09 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
###############################################################################
event.move <- function( community, species, individual )
{
individual <- as.matrix( individual )
## move only if individual is in stage that moves
if( all( is.move( community, species, individual ))) {
## move among substrates?
individual["sub.future",] <-
sampleOrgSubstrate( community, species, individual["sub.stage",] )
## move to new position in substrate
position <- paste( "pos", letters[1:3], sep = "." )
individual[position,] <- rtri( ncol( individual ), 10, individual[position,] )
}
if( ncol( individual ) == 1 )
individual <- individual[,1]
individual
}
###############################################################################
is.move <- function( community, species, individual )
{
!is.na( match( getOrgFeature( community, species, "move" ),
getOrgFuture( community, species, "current" )[ individual["stage",] ] ))
}
###############################################################################
event.find <- function( community, species, host, event )
{
individual <- get.individual( community, species )
substrate <- individual["sub.stage"]
## pending event: need to find available hosts on substrate
avail <- get.alive( community, host, substrate )
navail <- length( avail )
if( !navail )
return( avail )
## preferences based on schedule
interact <- get.interact( community, species, host, avail, event )
if( length( event ) > 1 )
interact <- apply( interact, 1, sum )
sinteract <- sum( interact )
if( sinteract ) {
if( navail > 1 ) {
found <- sample( avail, 1, prob = interact / sinteract )
}
else
found <- avail
found
}
else
numeric(0)
}
# --- Source: offspring.R ---
## $Id: init.R,v 1.0 2002/12/09 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## Offspring Information
##########################################################################################
getOffspring <- function( community, species,
offspring = getOrgFeature( community, species, "offspring" ))
{
if( is.na( offspring ))
return( 0 )
if( is.numeric( offspring ))
return( offspring )
getOrgInteract( community, offspring, species, "offspring" )
}
##########################################################################################
initOffspring <- function( community, species )
{
hostname <- getOrgFeature( community, species, "offspring" )
## find if there is offspring load based on host
orgoffspring <- getOffspring( community, species, hostname )
norganism <- sum( getOrgAlive( community, species ))
if( length( orgoffspring ) == 1 ) {
## mean offspring does not depend on any host
offspring <- stats::rpois( norganism, orgoffspring )
}
else {
orgoffspring <- orgoffspring[ orgoffspring > 0 ]
if(!length(orgoffspring))
return(community)
## figure out initial offspring load based on host distribution
## mean offspring depends on host stages and events
host <- get.species( community, hostname )
if( is.null( host ))
stop( paste( "Host", hostname, "not initiated yet" ))
## get weights of host stages in terms of future event times
host <- host[ , getOrgAlive( community, hostname ) ]
## find host stages that are preferred by parasite
## need to take subset of current that are actually in host
hoststages <- match( names( orgoffspring ), getOrgFuture( community, hostname )$current,
nomatch = 0 )
host <- as.matrix( host[ , !is.na( match( host["stage",], hoststages )) ] )
if( ncol( host ) == 0 )
return( community )
tmp <- !is.na( match( hoststages, host["stage",] ))
hoststages <- hoststages[tmp]
orgoffspring <- orgoffspring[tmp]
if(!length(orgoffspring))
return(community)
dd <- tapply( host["time",], host["stage",], sum )
dd[ as.character( hoststages[
is.na( match( hoststages, names( dd ))) ] ) ] <- 0
dd[ is.na( dd ) ] <- 0
sdd <- sum( dd )
if( length( dd ) > 1 & sdd > 0)
offspring <- as.vector( sample( orgoffspring, norganism, replace = TRUE,
prob = dd / sdd ))
else
offspring <- rep( ( sdd > 0 ) * orgoffspring[1], norganism )
offspring[ is.na( offspring ) ] <- 0
}
organism <- get.species( community, species )
organism["offspring",-1] <- offspring
put.species( community, species, organism )
}
###############################################################################
get.offspring <- function( community, species )
{
individual <- get.individual( community, species )
if( individual["offspring"] > 0 )
1
else
0
}
###########################################################################################
set.offspring <- function( community, species, host, dead )
{
stage <- get.species.element( community, host, "stage", dead )
current <- getOrgFuture( community, host, "current", stage )
offspring <- getOrgInteract( community, host, species, "offspring")
offspring <- as.vector( offspring[ as.character( current ) ] )
offspring[ is.na( offspring ) ] <- 0
offspring
}
# --- Source: sim.R ---
## $Id: sim.R,v 0.9 2002/12/09 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
### migrate count$count* to writeCount and then retrieve with readCount
###########################################################################################
updateCount <- function(community, species, individual, is.death = FALSE, step) {
community <- updateEvents(community, species, individual["future"])
if (!missing(step)) {
community <- setCount(community, elements = list(step = step))
}
## move individual through age classes or drop if it dies
countage <- getCount( community, species, "countage" )
ageclass <- getOrgAgeClass( community, species, individual["stage"] )
if( !is.na( ageclass )) {
ageclass <- as.character( ageclass )
countage[ageclass] <- max(countage[ageclass] - 1, 0)
}
ageclass <- getOrgAgeClass( community, species, individual["future"] )
if( !is.na( ageclass )) {
ageclass <- as.character( ageclass )
countage[ageclass] <- countage[ageclass] + 1
}
## move individual across substrate elements or drop if it dies
countsub <- getCount( community, species, "countsub" )
substrate <- getOrgFeature( community, species, "substrate" )
elements <- getOrgInteract( community, substrate, species, "substrate" )
element <- elements[ individual["sub.stage"] ]
subclass <- getOrgFeature( community, species, "subclass" )
include <- subclass == getOrgAgeClass( community, species, individual[c("stage","future")] )
# need 2 values; if only one, replicate it.
include <- rep_len(include, 2)
## leave old substrate
if (!is.na(element) && include[1]) {
element <- as.character( element )
countsub[element] <- max(countsub[element] - 1, 0)
}
if (!is.death && include[2]) {
## move to new substrate
newsub <- individual["sub.future"]
if( !is.na( newsub )) {
element <- elements[newsub]
if( !is.na( element )) {
element <- as.character( element )
countsub[element] <- countsub[element] + 1
}
}
}
## record counts
community <- setCount( community, species,
list( countage = countage, countsub = countsub ))
community <- writeCount( community, species, individual["time"], individual["future"],
countage, countsub )
}
###########################################################################################
updateCounts <- function( community, species, newbirths )
{
community <- updateEvents( community, species, 1, newbirths )
countage <- getCount( community, species, "countage" )
countsub <- getCount( community, species, "countsub" )
subclass <- getOrgFeature( community, species, "subclass" )
ageclass <- getOrgAgeClass( community, species, 1 )
individual <- get.individual( community, species )
if( !is.na( ageclass[1] )) {
ageclass <- as.character( ageclass )
countage[ageclass] <- newbirths + countage[ageclass]
if( subclass == ageclass ) {
substrate <- getOrgSubstrate( community, species, individual["sub.stage"] )
if( !is.na( substrate[1] )) {
substrate <- as.character( substrate )
countsub[substrate] <- newbirths + countsub[substrate]
}
}
}
community <- setCount( community, species,
list( countage = countage, countsub = countsub ))
community <- writeCount( community, species, individual["time"], 1, countage, countsub )
}
###########################################################################################
put.base <- function( community, species, id,
free = getCount( community, species, "free" ))
{
base <- get.base( community, species )
if( !missing( id ))
free <- leftist.free( free, id )
mintime <- max( getCount( community, species, "mintime" ),
getTime( community, species,
get.species.element( community, species, "time", base )))
setCount( community, species, list( base = base, free = free, mintime = mintime ))
}
# --- Source: simple.R ---
simpleServer <- function(id) {
shiny::moduleServer(id, function(input, output, session) {
ns <- session$ns
})
}
simpleInput <- function(id) {
ns <- shiny::NS(id)
shiny::tagList(
shiny::h4("Simulation Settings"),
shiny::sliderInput(ns("steps"),
label = "Simulation steps:",
min = 1000,
max = 10000,
value = 1000,
step = 500),
shiny::radioButtons(ns("nsim"),
"Number of Simulations",
c(1,10,20,50,100,200),
1, inline = TRUE),
shiny::actionButton(ns("go"), "Start Simulation"),
shiny::HTML("<hr style='height:1px;border:none;color:#333;background-color:#333;' />"),
shiny::h4("Save Files"),
shiny::uiOutput(ns("version"))
)
}
simpleOutput <- function(id) {
ns <- shiny::NS(id)
shiny::tagList(
shiny::radioButtons(ns("button"), "", c("Plots", "Input Data"),
"Plots", inline = TRUE),
)
}
simpleApp <- function(title = "Population Ethology") {
ui <- shiny::fluidPage(
shiny::titlePanel(title),
shiny::sidebarLayout(
shiny::sidebarPanel(
simpleInput("simple")
),
shiny::mainPanel(
simpleOutput("simple")
)))
server <- function(input, output, server) {
simpleServer("simple")
}
shiny::shinyApp(ui = ui, server = server)
}
# --- Source: ring.R ---
## $Id: ring.R,v 1.0 2002/12/09 yandell@stat.wisc.edu Exp $
##
## Functions for Bland Ewing's modeling.
##
## Copyright (C) 2000,2001,2002 Brian S. Yandell.
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
###########################################################################################
## Doubly linked rings
###########################################################################################
ring.add <- function( ring = data.frame( root = c( key = NA, left = 1 + nx, right = 2 )),
x )
{
nx <- length( x )
n1 <- ring["left","root"]
ring[[ as.character( n1 + 1 ) ]] <- c(
key = x[1],
left = suppressWarnings(as.numeric( ring["left","root"] )),
right = 1 )
if( nx > 1 ) for( i in seq( 2, nx ))
{
ring[[ as.character( n1 + i ) ]] <- c(
key = x[i],
left = n1+i-1,
right = 1 )
ring[ "right", as.character( n1 + i - 1 ) ] <- n1+i
}
ring["left","root"] <- n1 + nx
ring
}
###########################################################################################
ring.remove <- function( ring, P )
{
aP <- as.character( P )
if( is.na( match( aP, names( ring ))))
{
cat( paste( "Warning:", aP, "not found in ring\n" ))
return( ring )
}
left <- ring["left",aP]
right <- ring["right",aP]
ring[ "right", as.character( left ) ] <- right
ring[ "left", as.character( right ) ] <- left
ring[[aP]] <- NULL
ring
}
# --- Source: sierpinski.R ---
sierpinski <- function( stage = 5, reset = TRUE )
{
if( reset )
tmpar <- graphics::par( pty = "s", bty = "n", xaxt = "n", yaxt = "n", omi = rep(0,4),
mar = rep(0,4) )
aa <- 0:1
bb <- - aa
for( i in seq( stage )) {
tmp <- gasket( aa, bb )
aa <- tmp$aa
bb <- tmp$bb
tri <- tri2car.default( aa, bb )
r <- range( unlist( tri ))
plot( r, r, type = "n", xlab = "", ylab = "" )
graphics::lines(tri )
graphics::mtext( paste( "(", letters[1+i], ") Gasket of Order ", i, sep = "" ), 3, -2 )
# graphics::mtext( paste( "Gasket of order", i ), 3, -2 )
}
if( reset )
graphics::par( tmpar )
invisible( tri )
}
# --- Source: fileCount.R ---
putCount <- function( community, append = FALSE )
{
species <- get.species( community )
countage <- getCount( community,, "countage")
countsub <- getCount( community,, "countsub")
countbase <- getCount( community,, "base")
step_now <- if (append) (getCount( community,, "step" ) %||% 0) else 0
cnames <- c("step", "time", "future",
paste( "count",
seq( max( unlist( lapply( countage, length )) +
unlist( lapply( countsub, length )))),
sep = "" ))
cnums <- list()
for( i in species ) {
cnums[[i]] <- c(
step_now,
get.species.element( community, i, c("time","stage"), countbase[i] ),
countage[[i]],
countsub[[i]])
}
file <- getCount( community,, "file" )
if(!is.null(file) && is.character(file) && length(file) > 0) {
if( !( file.exists( file ) & append ))
cat( "species", cnames, "\n", file = file )
for( i in species ) {
cat( i, cnums[[i]], "\n", file = file, append = TRUE )
}
community
} else {
if(append) {
counts <- getCount( community,, "counts")
} else {
counts <- NULL
}
newcounts <- matrix(NA, length(species), length(cnames),
dimnames = list(species, cnames))
for(i in species) {
newcounts[i, seq_along(cnums[[i]])] <- cnums[[i]]
}
newcounts <- tibble::tibble(
data.frame(
species = species,
newcounts))
counts <- dplyr::bind_rows(
counts,
newcounts)
setCount( community,, list(counts = counts))
}
}
###########################################################################################
writeCount <- function( community, species, time, future, countage, countsub)
{
nstep <- getCount( community,, "step" )
cnums <- c(nstep, time, future, countage, countsub)
file <- getCount( community,, "file" )
if(!is.null(file) && is.character(file) && length(file) > 0) {
cat( species, cnums, "\n", file = file, append = TRUE )
community
} else {
counts <- getCount( community,, "counts")
cnames <- names(counts)[-1]
newcounts <- matrix(NA, 1, length(cnames),
dimnames = list(species, cnames))
newcounts[1, seq_along(cnums)] <- cnums
newcounts <- tibble::tibble(
data.frame(
species = species,
newcounts))
counts <- dplyr::bind_rows(
counts,
newcounts)
setCount( community,, list(counts = counts))
}
}
###########################################################################################
readCount <- function( community, species = unique(counts$species) )
{
file <- getCount( community,, "file" )
if(!is.null(file)) {
counts <- utils::read.table( file, header = TRUE, fill = TRUE )
} else {
counts <- getCount( community,, "counts")
}
count <- list()
for( i in species ) {
colnames <- c( levels( getOrgFuture( community, i, "ageclass" )),
levels( getOrgInteract( community,, i, "substrate" )))
count[[i]] <- as.matrix( counts[ counts$species == i, seq( 2, 4 + length( colnames )) ] )
dimnames( count[[i]] ) <- list( count[[i]][,"step"],
c( "step", "time", "future", colnames ))
}
count
}
# --- Source: count.join.R ---
count.join <- function( ... )
{
x <- list( ... )
numnum <- list()
for( i in seq( length( x ))) {
for( j in names( x[[i]] )) {
if( is.null( numnum[[j]] ))
numnum[[j]] <- x[[i]][[j]]
else
numnum[[j]] <- cbind( numnum[[j]], x[[i]][[j]] )
}
}
numnum
}
# --- Source: summary.ewing.R ---
summary.ewing <- function(object, ...) {
out <- list()
out$package <- object$org$package
out$species <- names(object$org$Future)
out$interact <- names(object$org$Interact)
out$meanvalue <- list()
for(i in out$species) {
out$meanvalue <- names(object$org$MeanValue[[i]])
}
if(length(object$pop)) {
out$stage <- list()
for(i in out$species)
out$stage[[i]] <-
table(object$org$Future[[i]]$current[getOrgAlive(object, i, "stage")])
}
if(!is.null(object$count)) {
out$events <- object$count$events
for( i in seq( length(out$events)))
out$events[[i]] <- apply(out$events[[i]], 2, function(x) {
tmp <- sum(x, na.rm = TRUE)
if(tmp > 0)
c(round(100 * x / tmp, 1), total = tmp)
else
c(x, total = 0)
})
}
out$cpu <- signif(object$cpu, 4)
class(out) <- c("summary.ewing", class(out))
out
}
print.summary.ewing <- function(x, ...) {
cat("Data initialization package:", x$package, "\n")
cat("Community species:", paste(x$species, collapse = ", " ), "\n")
cat("Community hosts:", paste(x$interact, collapse = ", "), "\n")
cat("Mean Value curves by species:")
mv <- FALSE
for(i in x$species) {
meanvalue <- x$meanvalue[[i]]
mv <- mv | !is.null(meanvalue)
if( !is.null(meanvalue)) {
cat("\n ", i, ":", paste( meanvalue, collapse = ", " ), "\n")
}
}
if(!mv)
cat(" none\n")
if(!is.null(x$stage)) {
cat("\nSimulation community has following counts:\n",
paste(x$species, lapply(x$stage, sum), sep = "=", collapse = ", "),
"\n")
print(x$stage)
}
if( !is.null( x$temp )) {
# ** later
}
if(!is.null(x$events)) {
print(x$events)
}
if(!is.null(x$cpu)) {
cat( "CPU timing by event in simulation\n" )
print(x$cpu)
}
}
# --- Source: summary_simobj.R ---
summary_simobj <- function(object) {
out <- sapply(object, function(x) {
if(is.list(x)) {
x <- names(x)
}
paste(x, collapse=",")
})
paste(paste(names(out), out, sep = ": "), collapse = "<br>")
}
# --- Source: ewing_discrete.R ---
ewing_discrete <- function(nsim, verbose = FALSE, ...) {
sims <- seq_len(nsim)
object <- as.list(sims)
names(object) <- sims
for(i in sims) {
if(verbose) cat(".")
object[[i]] <- ewing_discrete1(...)
}
make_ewing_discrete(object)
}
ewing_discrete1 <- function(siminit = init.simulation(interact = FALSE,
messages = FALSE, ...),
increment = 0.5, ...)
{
# Make sure increment is 1,2,5 x power of 10
incr <- pretty(increment)
increment <- incr[which.min(abs(incr - increment))[1]]
out <- future.events(siminit, refresh = 1000,
plotit = FALSE, messages = FALSE, ...)
attrs <- attributes(out)
# Get age classes used later for summaries and plots
items <- purrr::map(out$org$Future, function(x) levels(factor(x$ageclass)))
out <- readCount(out)
out <- purrr::map(
out,
function(x) {
purrr::map_df(
dplyr::distinct(
purrr::map_df(
dplyr::mutate(
as.data.frame(x),
time = ifelse(.data$step == 0, 0,
increment * ceiling(.data$time / increment))),
rev),
.data$time, .keep_all = TRUE),
rev)
})
attr(out, "count") <- attrs$count
attr(out, "nstep") <- attrs$nstep
attr(out, "items") <- items
out
}
summary.ewing_discrete <- function(object, ...) {
summary(ewing_envelopes(object), ...)
}
# --- Source: make_ewing_discrete.R ---
make_ewing_discrete <- function(object) {
nsim <- length(object)
class(object) <- c("ewing_discrete", class(object))
attr(object, "species") <- species <- names(object[[1]])
attr(object, "ordinate") <- "time"
attr(object, "count") <- attr(object[[1]], "count")
attr(object, "nstep") <- attr(object[[1]], "nstep")
attr(object, "items") <- attr(object[[1]], "items")
attr(object, "nsim") <- nsim
object
}
# --- Source: ggplot_ewing.R ---
ggplot_ewing <- function(object, step = 0, ageclass = TRUE,
substrate = !ageclass, ...)
{
if(!inherits(object, "ewing_snapshot")) {
object <- ewing_snapshot(object, step, ...)
}
step <- object$step
p <- list()
i <- 0
if(ageclass) {
i <- i + 1
p[[i]] <- ggplot2::autoplot(object$ageclass, ...)
}
if(substrate) {
species <- names(object$substrate)
for(j in species) {
i <- i + 1
p[[i]] <- ggplot2::autoplot(object$substrate[[j]], ...)
}
}
if(length(p) == 1) p <- p[[1]]
p
}
autoplot.ewing <- function(object, ...) {
ggplot_ewing(object, ...)
}
plot.ewing <- function(x, ...) {
ggplot_ewing(x, ...)
}
# --- Source: ggplot_current.R ---
ggplot_current <- function( x,
species,
col = as.character( future$color[stage] ),
headstuff = c( 0, "start"),
units = getOrgFeature( x, species, "units" ),
right = species, adj = c(0,.5,1),
position = paste( "pos", letters[1:3], sep = "." ),
pch = as.character( future$pch[stage] ), cex = 0.5,
stage = organism["stage",],
xlab = "horizontal", ylab = "vertical",
future = getOrgFuture( x, species, c("color","pch") ),
facet = TRUE, ...)
{
## plot current stages for species (except random parasites)
organism <- get.species( x, species )[,-1]
if(is.null(organism))
return(NULL)
tri_coord <- tri2car( organism[position,] )
tri_coord$col <- NA
tri_coord$label <- pch
values <- future$color
names(values) <- future$pch
values <- unique(values)
# Facet by Substrate.
# This needs to use generic function to get substrate names.
tmp <- names(x$org$Interact$substrate$substrate)
tri_coord$substrate <- tmp[organism["sub.stage",]]
p <- ggplot2::ggplot(tri_coord) +
ggplot2::aes(x, y, col = label, label = label) +
ggplot2::geom_text(size=3) +
ggplot2::xlab(xlab) +
ggplot2::ylab(ylab) +
ggplot2::ggtitle(paste(units, "future event", right),
subtitle = paste(headstuff, collapse = " ")) +
ggplot2::scale_color_manual(values = values)
if(facet) {
p <- p + ggplot2::facet_wrap(~substrate)
}
p
}
# --- Source: plot_null.R ---
plot_null <- function(msg = "no data") {
ggplot2::ggplot(data.frame(x = 1, y = 1),
ggplot2::aes(.data$x, .data$y, label = msg)) +
ggplot2::geom_text(size = 10) +
ggplot2::theme_void()
}
# --- Source: step_controls.R ---
step_size_choices <- c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000)
step_size_slider <- function(inputId, label = "Steps per click:", selected = 50) {
idx <- match(selected, step_size_choices)
if (is.na(idx)) idx <- 6
sl <- shiny::sliderInput(inputId, label, min = 1, max = length(step_size_choices), value = idx, step = 1, ticks = TRUE)
sl$children[[2]]$attribs[['data-values']] <- paste(step_size_choices, collapse = ",")
sl
}
parse_step_size <- function(val) {
if (is.null(val)) return(50)
num <- round(as.numeric(val))
if (is.na(num)) return(50)
# Direct large values (e.g., 20, 50, 100, 200, 500, 1000, 2000)
if (num %in% step_size_choices && num > 10) {
return(num)
}
# 0-based JavaScript index from ion.rangeSlider (0 to 10)
if (num >= 0 && num < length(step_size_choices)) {
return(step_size_choices[num + 1])
}
# Fallback for direct value
if (num %in% step_size_choices) {
return(num)
}
50
}
axisUnitInput <- function(id, time_label = "Time", selected = "step") {
ns <- shiny::NS(id)
choices <- c("Steps" = "step")
choices[time_label] <- "time"
shiny::radioButtons(ns("x_var"), "Display Units:",
choices = choices,
selected = selected, inline = TRUE)
}
axisUnitServer <- function(id) {
shiny::moduleServer(id, function(input, output, session) {
shiny::reactive({ if (!is.null(input$x_var)) input$x_var else "step" })
})
}
ageClassControlInput <- function(id, time_label = "Time") {
ns <- shiny::NS(id)
shiny::tagList(
axisUnitInput(ns("axis_unit"), time_label = time_label),
shiny::checkboxInput(ns("norm"), "Normalize Plot", TRUE),
shiny::checkboxInput(ns("total"), "Include Total in Plot", TRUE)
)
}
ageClassControlServer <- function(id) {
shiny::moduleServer(id, function(input, output, session) {
x_var <- axisUnitServer("axis_unit")
list(
x_var = x_var,
norm = shiny::reactive({ if (!is.null(input$norm)) input$norm else TRUE }),
total = shiny::reactive({ if (!is.null(input$total)) input$total else TRUE })
)
})
}
# --- Source: initParApp.R ---
initParApp <- function(title = "Population Ethology") {
ui <- bslib::page_sidebar(
title = title,
sidebar = bslib::sidebar(
initParInput("init_par")),
initParUI("init_par"),
initParOutput("init_par")
)
server <- function(input, output, server) {
init_par <- initParServer("init_par")
}
shiny::shinyApp(ui = ui, server = server)
}
initParServer <- function(id, simres = shiny::reactiveVal(NULL), datafile = shiny::reactiveVal("")) {
shiny::moduleServer(id, function(input, output, session) {
ns <- session$ns
species <- shiny::reactive({
get.organisms()$species
})
output$sppsize <- shiny::renderUI({
shiny::req(species()) # "host", "parasite"
lapply(species(), function(x) {
shiny::sliderInput(ns(x),
label = paste0("Number of ", x, "s:"),
min = 0,
max = 500,
value = 100,
step = 20)
})
})
datanames <- shiny::reactive({
getOrgNames()
})
output$inputfiles <- shiny::renderUI({
shiny::tagList(
shiny::selectInput(ns("dataname"), "", datanames(), "organism.features"),
DT::dataTableOutput(ns("org_table")))
})
output$org_table <- DT::renderDataTable({
getOrgDataSimple(simres(), shiny::req(input$dataname), datafile())
}, escape = FALSE, options = list(scrollX = TRUE, pageLength = 10))
# Show parameters
output$show_par <- shiny::renderUI({
nlist <- names(input)
# Remove any internal inputs, which have numbers.
glist <- grep("[0-9]", names(input))
if(length(glist))
nlist <- nlist[-glist]
# Construct output string.
out <- paste0("inputs: ", paste(nlist, collapse = ", "))
for(i in nlist) {
out <- paste(out, "<br>",
paste(i, input[[i]], sep = " = "))
}
shiny::HTML(out)
})
# Return.
input
})
}
initParInput <- function(id) {
ns <- shiny::NS(id)
shiny::uiOutput(ns("sppsize"))
}
initParUI <- function(id) {
ns <- shiny::NS(id)
shiny::uiOutput(ns("show_par"))
}
initParOutput <- function(id) {
ns <- shiny::NS(id)
shiny::uiOutput(ns("inputfiles"))
}
# --- Source: initApp.R ---
initApp <- function(title = "Population Ethology") {
ui <- bslib::page_sidebar(
title = title,
sidebar = bslib::sidebar(
initParInput("init_par")),
initOutput("init"),
substrateOutput("substrate")
)
server <- function(input, output, server) {
init_par <- initParServer("init_par")
siminit <- initServer("init", init_par)
substrateServer("substrate", siminit)
}
shiny::shinyApp(ui = ui, server = server)
}
initServer <- function(id, init_par) {
shiny::moduleServer(id, function(input, output, session) {
ns <- session$ns
siminit <- shiny::reactive({
init.simulation(count = as.numeric(c(shiny::req(init_par$host),
shiny::req(init_par$parasite))))
})
output$init <- shiny::renderUI({
out <- summary_simobj(shiny::req(siminit()))
shiny::HTML(out)
})
# Return.
siminit
})
}
initOutput <- function(id) {
ns <- shiny::NS(id)
shiny::uiOutput(ns("init"))
}
# --- Source: substrateApp.R ---
substrateApp <- function(title = "Substrate Organism Movement Explorer") {
ui <- bslib::page_sidebar(
title = title,
sidebar = bslib::sidebar(
initParInput("init_par"),
shiny::hr(),
substrateInput("substrate")),
substrateOutput("substrate")
)
server <- function(input, output, server) {
init_par <- initParServer("init_par")
siminit <- initServer("init", init_par)
substrateServer("substrate", siminit)
}
shiny::shinyApp(ui = ui, server = server)
}
substrateInput <- function(id) {
ns <- shiny::NS(id)
shiny::tagList(
shiny::div(
style = "font-size: 0.85rem;",
shiny::h4("Substrate Display & Stepping", style = "font-size: 1rem; font-weight: 600; margin-bottom: 8px;"),
shiny::checkboxGroupInput(ns("show_species"), "Species to Display:",
choices = c("Host" = "host", "Parasite" = "parasite"),
selected = c("host", "parasite"),
inline = TRUE),
shiny::radioButtons(ns("species_mode"), "Species Mode:",
choices = c("Overlay (1 Map)" = "overlay", "Separate (Adjacent Maps)" = "separate"),
selected = "overlay", inline = TRUE),
shiny::div(style = "border-top: 1px solid rgba(0,0,0,0.1); margin: 6px 0;"),
shiny::radioButtons(ns("layout"), "Layout View:",
choices = c("Hex Substrate Overlay" = "hex", "Faceted Substrates" = "facet"),
selected = "hex", inline = TRUE),
shiny::div(style = "border-top: 1px solid rgba(0,0,0,0.1); margin: 6px 0;"),
shiny::span("Simulation Stepping:", style = "font-weight: 600; color: #1a73e8; display: block; margin-bottom: 4px;"),
shiny::div(
style = "display: flex; gap: 4px; margin-bottom: 8px;",
shiny::actionButton(ns("step1"), "+1 Step", class = "btn-sm btn-outline-primary flex-fill"),
shiny::actionButton(ns("step10"), "+10 Steps", class = "btn-sm btn-outline-primary flex-fill"),
shiny::actionButton(ns("step100"), "+100 Steps", class = "btn-sm btn-outline-primary flex-fill"),
shiny::actionButton(ns("reset_sim"), "Reset", class = "btn-sm btn-outline-secondary")
),
shiny::div(style = "border-top: 1px solid rgba(0,0,0,0.1); margin: 6px 0;"),
shiny::checkboxGroupInput(ns("layers"), "Display Layers:",
choices = c("Substrate Boundaries" = "poly",
"Hex Grid Overlay" = "hex",
"Organisms" = "organisms",
"Substrate Names" = "centers",
"Side Numbers" = "labels"),
selected = c("poly", "hex", "organisms", "centers", "labels")),
shiny::div(style = "border-top: 1px solid rgba(0,0,0,0.1); margin: 6px 0;"),
shiny::div(
style = "display: flex; gap: 8px;",
shiny::numericInput(ns("width"), "Radius:", value = 10, min = 2, max = 30, step = 1),
shiny::numericInput(ns("step_density"), "Step Density:", value = 1, min = 0.5, max = 5, step = 0.5)
)
)
)
}
substrateServer <- function(id, simres, width = 10, step_density = 1) {
shiny::moduleServer(id, function(input, output, session) {
ns <- session$ns
current_sim <- shiny::reactiveVal(NULL)
shiny::observeEvent(simres(), {
current_sim(simres())
})
shiny::observeEvent(input$step1, {
sim <- current_sim()
if (!is.null(sim)) {
res <- future.events(sim, nstep = 1, plotit = FALSE)
current_sim(res)
}
})
shiny::observeEvent(input$step10, {
sim <- current_sim()
if (!is.null(sim)) {
res <- future.events(sim, nstep = 10, plotit = FALSE)
current_sim(res)
}
})
shiny::observeEvent(input$step100, {
sim <- current_sim()
if (!is.null(sim)) {
res <- future.events(sim, nstep = 100, plotit = FALSE)
current_sim(res)
}
})
shiny::observeEvent(input$reset_sim, {
current_sim(simres())
})
available_species <- shiny::reactive({
sim <- current_sim()
if (!is.null(sim) && !is.null(sim$pop)) names(sim$pop) else NULL
})
selected_species <- shiny::reactive({
avail <- available_species()
if (is.null(avail)) return(NULL)
sel <- input$show_species
if (is.null(sel) || length(sel) == 0) avail else intersect(sel, avail)
})
sppplot <- shiny::reactive({
spp <- selected_species()
shiny::req(spp)
sim <- current_sim()
shiny::req(sim)
layout_val <- if (!is.null(input$layout)) input$layout else "hex"
mode_val <- if (!is.null(input$species_mode)) input$species_mode else "overlay"
w_val <- if (!is.null(input$width)) input$width else width
sd_val <- if (!is.null(input$step_density)) input$step_density else step_density
layers_val <- if (!is.null(input$layers)) input$layers else c("poly", "hex", "organisms", "centers", "labels")
sim_single <- if (inherits(sim, "ewing_discrete") && is.list(sim) && length(sim) > 0) sim[[1]] else sim
if (inherits(sim_single, "ewing")) {
if (mode_val == "overlay" && layout_val == "hex") {
sub_data <- ewing_substrate(sim_single, spp, layout = layout_val, width = w_val, step_density = sd_val)
if (!is.null(sub_data)) {
p_obj <- ggplot_ewing_substrate(sub_data, layout = layout_val, width = w_val, step_density = sd_val, layers = layers_val)
list(p_obj)
} else {
list()
}
} else {
p <- lapply(spp, function(x) {
sub_data <- ewing_substrate(sim_single, x, layout = layout_val, width = w_val, step_density = sd_val)
if (!is.null(sub_data)) {
p_obj <- ggplot_ewing_substrate(sub_data, layout = layout_val, width = w_val, step_density = sd_val, layers = layers_val)
p_obj
} else {
NULL
}
})
p[!sapply(p, is.null)]
}
} else {
list()
}
})
output$sppPlot <- shiny::renderPlot({
plots <- sppplot()
if (!is.null(plots) && length(plots) > 0) {
cowplot::plot_grid(plotlist = plots, ncol = length(plots), align = "h")
} else {
ggplot2::ggplot() + ggplot2::theme_void() + ggplot2::ggtitle("No active species selected to plot")
}
})
output$substrate_plot <- shiny::renderUI({
plots <- sppplot()
h_px <- 500
shiny::plotOutput(ns("sppPlot"), height = paste0(h_px, "px"))
})
# Return current simulation state for downstream composition
current_sim
})
}
substrateOutput <- function(id) {
ns <- shiny::NS(id)
shiny::uiOutput(ns("substrate_plot"))
}
# --- Source: distPlotApp.R ---
distPlotApp <- function(title = "Population Ethology") {
ui <- bslib::page_sidebar(
title = title,
sidebar = bslib::sidebar(
initParInput("init_par"),
distPlotInput("dist_plot"),
futureInput("future")),
distPlotOutput("dist_plot")
)
server <- function(input, output, server) {
init_par <- initParServer("init_par")
siminit <- initServer("init", init_par)
simres <- futureServer("future", siminit)
distPlotServer("dist_plot", simres)
}
shiny::shinyApp(ui = ui, server = server)
}
distPlotServer <- function(id, simres, x_var = NULL, total = NULL, norm = NULL) {
shiny::moduleServer(id, function(input, output, session) {
ns <- session$ns
dist_plot <- shiny::reactive({
sim <- if (is.reactive(simres)) simres() else simres
shiny::req(sim)
tot_val <- if (is.reactive(total)) total() else if (!is.null(total)) total else if (!is.null(input$total)) input$total else TRUE
norm_val <- if (is.reactive(norm)) norm() else if (!is.null(norm)) norm else if (!is.null(input$norm)) input$norm else TRUE
xv_val <- if (is.reactive(x_var)) x_var() else if (!is.null(x_var)) x_var else if (!is.null(input$x_var)) input$x_var else "step"
if (is.null(tot_val)) tot_val <- TRUE
if (is.null(norm_val)) norm_val <- TRUE
if (is.null(xv_val) || !xv_val %in% c("step", "time")) xv_val <- "step"
object <- tryCatch({
ewing_ageclass(sim, total = tot_val, normalize = norm_val)
}, error = function(e) NULL)
if (is.null(object)) return(plot_null("no simulation"))
ggplot2::autoplot(object, x_var = xv_val)
})
output$dist_plot <- shiny::renderPlot({
dist_plot()
})
# Return.
dist_plot
})
}
distPlotInput <- function(id) {
ns <- shiny::NS(id)
shiny::tagList(
shiny::checkboxInput(ns("norm"), "Normalize Plot", TRUE),
shiny::checkboxInput(ns("total"), "Include Total", TRUE))
}
distPlotOutput <- function(id) {
ns <- shiny::NS(id)
shiny::plotOutput(ns("dist_plot"), height = "400px")
}
# --- Source: multApp.R ---
multApp <- function(title = "Population Ethology") {
ui <- shiny::fluidPage(
shiny::titlePanel(title),
shiny::sidebarLayout(
shiny::sidebarPanel(
multInput("ewing")
),
shiny::mainPanel(
multOutput("ewing")
)))
server <- function(input, output, server) {
multServer("ewing")
}
shiny::shinyApp(ui = ui, server = server)
}
multServer <- function(id) {
shiny::moduleServer(id, function(input, output, session) {
ns <- session$ns
simres <- shiny::reactive({
nsim <- as.integer(shiny::req(input$nsim))
shiny::withProgress(
message = paste('Ewing Discrete', nsim, 'Simulations ...'),
value = 0,
{
out <- as.list(seq_len(nsim))
inc <- 1 / nsim
for(i in seq_len(nsim)) {
shiny::incProgress(inc)
out[[i]] <- ewing_discrete1(
count = as.numeric(c(input$host, input$parasite)),
nstep = input$steps)
}
}
)
})
distplot <- shiny::reactive({
if(inherits(simres(), "ewing")) {
ggplot2::autoplot(ewing_ageclass(simres(), total = input$total,
normalize = input$norm))
} else {
NULL
}
})
output$distPlot <- shiny::renderPlot({
distplot()
})
# *** This is not right. Need to get each species name here and in `ewing_substrate`
# species <- ewing:::getOrgFeature(simres)
# gives list but includes substrates.
# can figure out what substrate goes to species with
# ewing:::getOrgFeature(simres, species[i], "substrate")
# if it is NA (or "NA"), then that is a substrate.
# So cycle through species generating plots.
# put as much in `ewing_substrate` as possible.
species <- shiny::reactive({
get.organisms(datafile())$species
})
substrates <- shiny::reactive({
get.organisms(datafile())$substrates
})
output$sppsize <- shiny::renderUI({
shiny::req(species())
lapply(species(), function(x) {
shiny::sliderInput(ns(x),
label = paste0("Number of ", x, "s:"),
min = 0,
max = 500,
value = 100,
step = 20)
})
})
sppplot <- shiny::reactive({
shiny::req(species(), simres())
if(inherits(simres(), "ewing")) {
if(!is.null(simres())) {
p <- lapply(species(), function(x) {
p <- ggplot2::autoplot(ewing_substrate(simres(), x))
if(inherits(p, "ggplot"))
p <- p + ggplot2::ggtitle(paste(x, "on", substrates()[1]))
p
})
if(any(unlist(purrr::map(p, is.null))))
p <- NULL
p
}
} else {
ggplot2::ggplot()
}
})
output$sppPlot <- shiny::renderPlot({
if(!is.null(sppplot())) {
spp <- length(species())
cowplot::plot_grid(plotlist = sppplot(), nrow = spp)
} else {
ggplot2::ggplot()
}
})
envdata <- shiny::reactive({
shiny::req(simres())
if(inherits(simres(), "ewing_discrete")) {
ewing_envelopes(simres())
} else {
NULL
}
})
envelopePlot <- shiny::reactive({
shiny::req(envdata())
nsim <- as.integer(shiny::req(input$nsim))
conf <- (nsim >= 10) & input$conf
if(inherits(simres(), "ewing_discrete")) {
ggplot_ewing_envelopes(envdata(), conf)
} else {
NULL
}
})
output$envPlot <- shiny::renderPlot({
envelopePlot()
})
output$plots <- shiny::renderUI({
nsim <- as.integer(shiny::req(input$nsim), simres())
if(nsim == 1) {
shiny::req(species())
shiny::tagList(
shiny::plotOutput(ns("distPlot"), height = "400px"),
shiny::plotOutput(ns("sppPlot"), height = paste0(200 * length(species()), "px")))
} else {
shiny::plotOutput(ns("envPlot"))
}
})
data <- reactive({
nsim <- as.integer(shiny::req(input$nsim))
species <- shiny::req(input$species)
if(nsim == 1) {
readCount(simres())[[species]]
} else {
shiny::req(envdata())
print(envdata(), species = species)
}
})
params <- shiny::reactive({
nsim <- shiny::req(input$nsim)
paste(shiny::req(input$host), shiny::req(input$parasite),
shiny::req(input$steps), nsim, sep = "_")
})
output$downloadRun <- shiny::downloadHandler(
filename = function() {
paste0(paste(shiny::req(input$outfile), shiny::req(input$species), params(), sep = "_"), ".csv") },
content = function(file) {
utils::write.csv(data(), file, row.names = FALSE)
}
)
output$downloadPlot <- shiny::downloadHandler(
filename = function() {
paste0(paste(shiny::req(input$plotfile), params(), sep = "_"), ".pdf") },
content = function(file) {
grDevices::pdf(file, width = 9)
nsim <- as.integer(shiny::req(input$nsim))
if(nsim == 1) {
print(distplot())
for(i in species()) {
print(sppplot()[[i]])
}
} else {
print(envelopePlot())
}
grDevices::dev.off()
}
)
datanames <- shiny::reactive({
getOrgNames(datafile())
})
output$inputfiles <- shiny::renderUI({
shiny::tagList(
shiny::selectInput(ns("dataname"), "", datanames(), "organism.features"),
DT::renderDataTable({
getOrgDataSimple(simres(),shiny::req(input$dataname), datafile())
}, escape = FALSE,
options = list(scrollX = TRUE, pageLength = 10)))
})
datafile <- shiny::reactive({
if(shiny::isTruthy(input$datafile)) {
input$datafile$datapath
} else {
""
}
})
output$outs <- shiny::renderUI({
switch(input$button,
Plots = shiny::uiOutput(ns("plots")),
"Input Data" = shiny::uiOutput(ns("inputfiles")))
})
output$plottype <- shiny::renderUI({
if(input$nsim == 1) {
shiny::tagList(
shiny::checkboxInput(ns("norm"), "Normalize Plot", TRUE),
shiny::checkboxInput(ns("total"), "Include Total", TRUE))
} else {
shiny::checkboxInput(ns("conf"), "Confidence band", FALSE)
}
})
output$version <- shiny::renderText({
paste("Ewing package version ", utils::packageVersion("ewing"))
})
})
}
multInput <- function(id) {
ns <- shiny::NS(id)
shiny::tagList(
shiny::h4("Simulation Settings"),
shiny::uiOutput(ns("sppsize")),
shiny::sliderInput(ns("steps"),
label = "Simulation steps:",
min = 1000,
max = 10000,
value = 1000,
step = 500),
shiny::radioButtons(ns("nsim"),
"Number of Simulations",
c(1,10,20,50,100,200),
1, inline = TRUE),
shiny::fileInput(ns("datafile"), "Optional XLSX Input Data File",
multiple = FALSE,
accept = c(".xls", ".xlsx")),
shiny::actionButton(ns("go"), "Start Simulation"),
shiny::HTML("<hr style='height:1px;border:none;color:#333;background-color:#333;' />"),
shiny::uiOutput(ns("plottype")),
shiny::h4("Save Files"),
shiny::fluidRow(
shiny::column(6, shiny::textInput(ns("outfile"), "Species Table", "mysim")),
shiny::column(3, shiny::selectInput(ns("species"), "", c("host", "parasite"), "host")),
shiny::column(3, shiny::downloadButton(ns("downloadRun"), "Table"))),
shiny::fluidRow(
shiny::column(9, shiny::textInput(ns("plotfile"), "Plot File", "myplot")),
shiny::column(3, shiny::downloadButton(ns("downloadPlot"), "Plots"))),
shiny::HTML("<hr style='height:1px;border:none;color:#333;background-color:#333;' />"),
shiny::HTML("See <a href='https://github.com/byandell/ewing'>ewing package on github</a>"),
shiny::uiOutput(ns("version"))
)
}
multOutput <- function(id) {
ns <- shiny::NS(id)
shiny::tagList(
shiny::radioButtons(ns("button"), "", c("Plots", "Input Data"),
"Plots", inline = TRUE),
shiny::uiOutput(ns("outs"))
)
}
# --- Source: inputApp.R ---
inputApp <- function(title = "Input Data Explorer") {
ui <- bslib::page_sidebar(
title = title,
sidebar = bslib::sidebar(
inputAppInput("input_app")
),
bslib::card(
inputAppOutput("input_app")
)
)
server <- function(input, output, session) {
inputAppServer("input_app")
}
shiny::shinyApp(ui = ui, server = server)
}
discover_dataset_tables <- function(datafile = "", sim = NULL) {
found <- character(0)
# 1. Inspect datafile directory
d_path <- if (is.character(datafile) && datafile != "") datafile else if (!is.null(sim) && !is.null(sim$datafile)) sim$datafile else ""
if (d_path != "" && file.exists(d_path)) {
if (dir.exists(d_path)) {
files <- list.files(d_path, pattern = "\\.(txt|csv|rds)$", full.names = FALSE)
if (length(files) > 0) {
found <- tools::file_path_sans_ext(files)
}
} else if (grepl("\\.xlsx$", d_path, ignore.case = TRUE)) {
sheets <- tryCatch(readxl::excel_sheets(d_path), error = function(e) character(0))
if (length(sheets) > 0) found <- sheets
}
}
# 2. Inspect sim$datasets if present
if (!is.null(sim) && !is.null(sim$datasets)) {
found <- unique(c(found, names(sim$datasets)))
}
# 3. Default fallback choices if nothing found
if (length(found) == 0) {
found <- c(
"organism.features", "future.host", "future.parasite",
"substrate.host", "substrate.parasite", "substrate.substrate",
"host.parasite", "temperature.base", "temperature.par"
)
}
unique(found)
}
inputAppInput <- function(id, choices = NULL) {
ns <- shiny::NS(id)
default_choices <- if (!is.null(choices)) choices else c(
"organism.features", "future.host", "future.parasite",
"substrate.host", "substrate.parasite", "substrate.substrate",
"host.parasite", "temperature.base", "temperature.par"
)
shiny::tagList(
shiny::selectInput(ns("dataname"), "Select Dataset Table:",
choices = default_choices,
selected = default_choices[1]
)
)
}
inputAppOutput <- function(id) {
ns <- shiny::NS(id)
shiny::tableOutput(ns("org_table"))
}
inputAppServer <- function(id, simres = shiny::reactiveVal(NULL), datafile = shiny::reactiveVal("")) {
shiny::moduleServer(id, function(input, output, session) {
ns <- session$ns
# Dynamically update select choices based on folder / simulation contents
shiny::observe({
sim <- if (is.reactive(simres)) simres() else simres
dfile <- if (is.reactive(datafile)) datafile() else datafile
discovered <- discover_dataset_tables(dfile, sim)
if (length(discovered) > 0) {
current_sel <- input$dataname
sel <- if (!is.null(current_sel) && current_sel %in% discovered) current_sel else discovered[1]
shiny::updateSelectInput(session, "dataname", choices = discovered, selected = sel)
}
})
output$org_table <- shiny::renderTable(
{
name <- input$dataname %||% "organism.features"
sim <- if (is.reactive(simres)) simres() else simres
dfile <- if (is.reactive(datafile)) datafile() else datafile
res <- NULL
# 0. Check direct file in datafile directory if points to folder
if (is.character(dfile) && dfile != "" && dir.exists(dfile)) {
txt_path <- file.path(dfile, paste0(name, ".txt"))
if (file.exists(txt_path)) res <- tryCatch(utils::read.table(txt_path, header = TRUE, sep = "\t", stringsAsFactors = FALSE), error = function(e) NULL)
if (is.null(res)) {
csv_path <- file.path(dfile, paste0(name, ".csv"))
if (file.exists(csv_path)) res <- tryCatch(utils::read.csv(csv_path, stringsAsFactors = FALSE), error = function(e) NULL)
}
}
# 1. Check if dataset is stored in sim$datasets (e.g. injected in webR demo)
if (is.null(res) && !is.null(sim) && !is.null(sim$datasets) && !is.null(sim$datasets[[name]])) {
res <- sim$datasets[[name]]
}
# 2. Extract dynamically via getOrgDataSimple or getOrg* package routines
if (is.null(res) && !is.null(sim) && inherits(sim, "ewing")) {
sim_single <- if (inherits(sim, "ewing_discrete")) sim[[1]] else sim
res <- tryCatch({
getOrgDataSimple(sim_single, name, datafile = dfile)
}, error = function(e) NULL)
}
# 3. Dynamic fallback to extracting from sim_single$org state structures
if ((is.null(res) || !is.data.frame(res) || nrow(res) == 0) && !is.null(sim) && inherits(sim, "ewing")) {
sim_single <- if (inherits(sim, "ewing_discrete")) sim[[1]] else sim
res <- tryCatch({
left <- stringr::str_remove(name, "\\..*")
right <- stringr::str_remove(name, ".*\\.")
if (left == "organism" && right == "features") {
if (!is.null(sim_single$org$Feature)) as.data.frame(sim_single$org$Feature) else NULL
} else if (left == "future") {
if (!is.null(sim_single$org$Future[[right]])) sim_single$org$Future[[right]] else getOrgFuture(sim_single, right)
} else if (!is.null(sim_single$org$Interact[[left]][[right]])) {
sim_single$org$Interact[[left]][[right]]
} else if (!is.null(sim_single$org[[left]][[right]])) {
sim_single$org[[left]][[right]]
} else {
NULL
}
}, error = function(e) NULL)
}
if (is.null(res) || !is.data.frame(res) || nrow(res) == 0) {
res <- data.frame(Info = paste("Dataset", name, "is not available in current simulation instance."))
}
res
},
striped = TRUE,
hover = TRUE,
bordered = TRUE
)
})
}
# --- Source: habitat.R ---
get_habitat_features <- function(watershed_obj,
categories = c("lakes", "waterways", "forests", "bogs"),
use_cache = TRUE) {
huc_layer <- watershed_obj$layer
if (is.null(huc_layer)) {
stop("Invalid watershed_obj: layer component is required.")
}
# Check for pre-created cached data for Isle Royale
cache_dir <- system.file("extdata/isle_royale", package = "ewing")
if (cache_dir == "") cache_dir <- "inst/extdata/isle_royale"
cache_file <- file.path(cache_dir, "isle_royale_features.rds")
if (use_cache && file.exists(cache_file)) {
cached_sf <- tryCatch(readRDS(cache_file), error = function(e) NULL)
if (!is.null(cached_sf) && inherits(cached_sf, "sf")) {
cached_sf <- sf::st_transform(cached_sf, sf::st_crs(huc_layer))
clipped <- suppressWarnings(sf::st_intersection(cached_sf, huc_layer))
if (nrow(clipped) > 0) return(clipped)
}
}
# Query OpenStreetMap if osmdata package is installed
if (!requireNamespace("osmdata", quietly = TRUE)) {
warning("The 'osmdata' package is not installed. Returning fallback habitat geometries.")
return(get_fallback_habitat_features(huc_layer))
}
bbox <- sf::st_bbox(sf::st_transform(huc_layer, 4326))
bbox_str <- paste(bbox["ymin"], bbox["xmin"], bbox["ymax"], bbox["xmax"], sep = ",")
old_url <- osmdata::get_overpass_url()
osmdata::set_overpass_url("https://lz4.overpass-api.de/api/interpreter")
features_list <- list()
old_s2 <- sf::sf_use_s2()
sf::sf_use_s2(FALSE)
on.exit({
sf::sf_use_s2(old_s2)
osmdata::set_overpass_url(old_url)
}, add = TRUE)
# Extract Inland Lakes & Waterbodies
if ("lakes" %in% categories) {
ql_lakes <- paste0(
"[out:xml][timeout:30];\n(\n",
" natural[\"water\"](", bbox_str, ");\n",
" waterway[\"riverbank\"](", bbox_str, ");\n",
");\nout body;\n>;\nout skel qt;\n"
)
res_lakes <- tryCatch(osmdata::osmdata_sf(ql_lakes), error = function(e) NULL)
if (!is.null(res_lakes) && !is.null(res_lakes$osm_polygons) && nrow(res_lakes$osm_polygons) > 0) {
poly <- res_lakes$osm_polygons
poly$habitat_type <- "Lake/Pond"
features_list$lakes <- poly[, c("habitat_type", "geometry")]
}
}
# Extract Waterways & Beaver Ponds
if ("waterways" %in% categories) {
ql_water <- paste0(
"[out:xml][timeout:30];\n(\n",
" waterway[\"stream\"](", bbox_str, ");\n",
" waterway[\"river\"](", bbox_str, ");\n",
" waterway[\"drain\"](", bbox_str, ");\n",
");\nout body;\n>;\nout skel qt;\n"
)
res_water <- tryCatch(osmdata::osmdata_sf(ql_water), error = function(e) NULL)
if (!is.null(res_water) && !is.null(res_water$osm_lines) && nrow(res_water$osm_lines) > 0) {
lines <- res_water$osm_lines
lines$habitat_type <- "Waterway"
features_list$waterways <- lines[, c("habitat_type", "geometry")]
}
}
# Extract Cool Shaded Forests
if ("forests" %in% categories) {
ql_forest <- paste0(
"[out:xml][timeout:30];\n(\n",
" landuse[\"forest\"](", bbox_str, ");\n",
" natural[\"wood\"](", bbox_str, ");\n",
");\nout body;\n>;\nout skel qt;\n"
)
res_forest <- tryCatch(osmdata::osmdata_sf(ql_forest), error = function(e) NULL)
if (!is.null(res_forest) && !is.null(res_forest$osm_polygons) && nrow(res_forest$osm_polygons) > 0) {
poly <- res_forest$osm_polygons
poly$habitat_type <- "Forest"
features_list$forests <- poly[, c("habitat_type", "geometry")]
}
}
# Extract Bogs & Wetlands
if ("bogs" %in% categories) {
ql_bogs <- paste0(
"[out:xml][timeout:30];\n(\n",
" natural[\"wetland\"](", bbox_str, ");\n",
" wetland[\"bog\"](", bbox_str, ");\n",
" wetland[\"marsh\"](", bbox_str, ");\n",
");\nout body;\n>;\nout skel qt;\n"
)
res_bogs <- tryCatch(osmdata::osmdata_sf(ql_bogs), error = function(e) NULL)
if (!is.null(res_bogs) && !is.null(res_bogs$osm_polygons) && nrow(res_bogs$osm_polygons) > 0) {
poly <- res_bogs$osm_polygons
poly$habitat_type <- "Bog/Wetland"
features_list$bogs <- poly[, c("habitat_type", "geometry")]
}
}
if (length(features_list) == 0) {
return(get_fallback_habitat_features(huc_layer))
}
combined <- do.call(rbind, features_list)
combined <- sf::st_transform(combined, sf::st_crs(huc_layer))
combined <- suppressWarnings(sf::st_make_valid(combined))
clipped <- suppressWarnings(sf::st_intersection(combined, huc_layer))
if (nrow(clipped) == 0) {
return(get_fallback_habitat_features(huc_layer))
}
return(clipped)
}
get_fallback_habitat_features <- function(huc_layer) {
crs_target <- sf::st_crs(huc_layer)
# Landmark reference coordinates (WGS84)
# Washington Creek Windigo (-89.146, 47.923)
# Ojibway Lake (-88.618, 48.113)
# Feldtmann Lake (-88.961, 47.876)
# Hidden Lake (-88.490, 48.151)
mk_poly <- function(lon, lat, dx = 0.015, dy = 0.01) {
pts <- matrix(c(
lon - dx, lat - dy,
lon + dx, lat - dy,
lon + dx, lat + dy,
lon - dx, lat + dy,
lon - dx, lat - dy
), ncol = 2, byrow = TRUE)
sf::st_polygon(list(pts))
}
p_ojibway <- mk_poly(-88.618, 48.113, 0.012, 0.008)
p_feldtmann <- mk_poly(-88.961, 47.876, 0.018, 0.010)
p_hidden <- mk_poly(-88.490, 48.151, 0.010, 0.006)
p_wash_forest <- mk_poly(-89.146, 47.923, 0.025, 0.015)
p_bog <- mk_poly(-88.750, 48.020, 0.020, 0.012)
geom_sfc <- sf::st_sfc(p_ojibway, p_feldtmann, p_hidden, p_wash_forest, p_bog, crs = 4326)
df <- data.frame(
habitat_type = c("Lake/Pond", "Lake/Pond", "Lake/Pond", "Forest", "Bog/Wetland"),
stringsAsFactors = FALSE
)
fallback_sf <- sf::st_sf(df, geometry = geom_sfc)
fallback_sf <- sf::st_transform(fallback_sf, crs_target)
clipped <- suppressWarnings(sf::st_intersection(fallback_sf, huc_layer))
return(clipped)
}
get_moose_landmarks <- function(watershed_obj, use_cache = TRUE) {
huc_layer <- watershed_obj$layer
cache_dir <- system.file("extdata/isle_royale", package = "ewing")
if (cache_dir == "") cache_dir <- "inst/extdata/isle_royale"
cache_file <- file.path(cache_dir, "isle_royale_landmarks.rds")
if (use_cache && file.exists(cache_file)) {
cached_pts <- tryCatch(readRDS(cache_file), error = function(e) NULL)
if (!is.null(cached_pts) && inherits(cached_pts, "sf")) {
return(sf::st_transform(cached_pts, sf::st_crs(huc_layer)))
}
}
df <- data.frame(
name = c(
"Washington Creek (Windigo)",
"Ojibway Lake",
"Feldtmann Lake",
"Hidden Lake (Tobin Harbor)"
),
location = c("Windigo", "Ojibway", "Feldtmann", "Tobin Harbor"),
description = c(
"Feeding area along stream & forest cover",
"Aquatic vegetation feeding lake",
"Major southwest inland lake habitat",
"Aquatic plant feeding area near Tobin Harbor"
),
lon = c(-89.146, -88.618, -88.961, -88.490),
lat = c(47.923, 48.113, 47.876, 48.151),
stringsAsFactors = FALSE
)
pts_sf <- sf::st_as_sf(df, coords = c("lon", "lat"), crs = 4326)
pts_sf <- sf::st_transform(pts_sf, sf::st_crs(huc_layer))
return(pts_sf)
}
add_habitat_hex_overlay <- function(hex_obj, habitat_sf = NULL, landmarks_sf = NULL, features = NULL, landmarks = NULL) {
if (is.null(habitat_sf)) habitat_sf <- features
if (is.null(landmarks_sf)) landmarks_sf <- landmarks
if (is.character(habitat_sf) && file.exists(habitat_sf)) {
habitat_sf <- tryCatch(readRDS(habitat_sf), error = function(e) NULL)
}
if (is.character(landmarks_sf) && file.exists(landmarks_sf)) {
landmarks_sf <- tryCatch(readRDS(landmarks_sf), error = function(e) NULL)
}
if (is.null(habitat_sf)) {
habitat_sf <- get_habitat_features(hex_obj)
}
if (is.null(landmarks_sf)) {
landmarks_sf <- get_moose_landmarks(hex_obj)
}
hex_mesh <- hex_obj$hex_overlay
# Calculate habitat preference score per hexagon
# Base weight = 1.0; bonus for lake (+2.0), waterway (+1.5), forest (+1.0), bog (+1.8)
scores <- numeric(length(hex_mesh))
types_list <- character(length(hex_mesh))
if (!is.null(habitat_sf) && nrow(habitat_sf) > 0) {
inter <- safe_st_intersects(hex_mesh, habitat_sf)
for (i in seq_along(inter)) {
indices <- inter[[i]]
if (length(indices) > 0) {
sub_types <- habitat_sf$habitat_type[indices]
score <- 1.0
if ("Lake/Pond" %in% sub_types) score <- score + 2.0
if ("Waterway" %in% sub_types) score <- score + 1.5
if ("Bog/Wetland" %in% sub_types) score <- score + 1.8
if ("Forest" %in% sub_types) score <- score + 1.0
scores[i] <- score
types_list[i] <- paste(unique(sub_types), collapse = ", ")
} else {
scores[i] <- 1.0
types_list[i] <- "Upland/Open"
}
}
} else {
scores[] <- 1.0
types_list[] <- "General"
}
hex_sf <- sf::st_sf(
hex_id = seq_along(hex_mesh),
habitat_score = scores,
habitat_type = types_list,
geometry = hex_mesh
)
res <- hex_obj
res$habitat_sf <- habitat_sf
res$landmarks_sf <- landmarks_sf
res$hex_habitat_sf <- hex_sf
class(res) <- c("habitat_hex_overlay", class(hex_obj))
return(res)
}
autoplot.habitat_hex_overlay <- function(object, show_landmarks = TRUE, ...) {
p <- ggplot2::ggplot() +
# Underlying island boundary
ggplot2::geom_sf(data = object$layer, fill = "#eef4f8", color = "#2c3e50", linewidth = 0.8)
# Render habitat layers if available
if (!is.null(object$habitat_sf) && nrow(object$habitat_sf) > 0) {
p <- p + ggplot2::geom_sf(
data = object$habitat_sf,
ggplot2::aes(fill = .data$habitat_type),
alpha = 0.5, color = NA
)
}
# Hexagonal substrate overlay colored by habitat suitability score
if (!is.null(object$hex_habitat_sf)) {
p <- p + ggplot2::geom_sf(
data = object$hex_habitat_sf,
ggplot2::aes(color = .data$habitat_score),
fill = NA, linewidth = 0.6
) +
ggplot2::scale_color_viridis_c(option = "viridis", name = "Habitat Weight")
}
# Render Moose Sighting Landmarks
if (show_landmarks && !is.null(object$landmarks_sf) && nrow(object$landmarks_sf) > 0) {
p <- p +
ggplot2::geom_sf(data = object$landmarks_sf, color = "#d35400", size = 3, shape = 18) +
ggplot2::geom_sf_text(
data = object$landmarks_sf,
ggplot2::aes(label = .data$name),
color = "#900c3f", size = 3, fontface = "bold", vjust = -0.7
)
}
title_txt <- "Isle Royale Moose Habitat & Substrate Overlay Model"
if (!is.null(object$feature_name) && object$feature_name != "") {
title_txt <- paste0(title_txt, " (", object$feature_name, ")")
}
p +
ggplot2::theme_minimal() +
ggplot2::ggtitle(title_txt) +
ggplot2::labs(
x = "Longitude", y = "Latitude",
caption = "Habitats: Inland Lakes, Beaver Ponds/Waterways, Shaded Forests & Bogs"
)
}
add_leaflet_habitat_overlay <- function(map, object) {
if (is.null(map) || is.null(object)) return(map)
# Render habitat suitability polygons/hexagons
if (!is.null(object$hex_habitat_sf)) {
hex_wgs <- sf::st_transform(object$hex_habitat_sf, 4326)
map <- map |>
leaflet::addPolygons(
data = hex_wgs,
color = "#e74c3c",
weight = 1,
fillOpacity = 0.15,
popup = paste0("<b>Hex ID:</b> ", hex_wgs$hex_id,
"<br/><b>Habitat Score:</b> ", hex_wgs$habitat_score,
"<br/><b>Habitat Types:</b> ", hex_wgs$habitat_type),
group = "Habitat Substrate Mesh"
)
}
# Render Moose Sighting Area Markers
if (!is.null(object$landmarks_sf)) {
lm_wgs <- sf::st_transform(object$landmarks_sf, 4326)
map <- map |>
leaflet::addCircleMarkers(
data = lm_wgs,
color = "#d35400",
radius = 7,
fillOpacity = 0.9,
popup = paste0("<b>Moose Sighting Area:</b> ", lm_wgs$name,
"<br/><b>Description:</b> ", lm_wgs$description),
group = "Moose Sighting Areas"
)
}
return(map)
}
# --- Source: watershed.R ---
get_watershed <- function(huc_id, feature_name = NULL, huc_layer = NULL) {
# Standardize huc_id input (support comma-separated string or character vector)
if (is.character(huc_id) && length(huc_id) == 1 && grepl(",", huc_id)) {
huc_id <- trimws(unlist(strsplit(huc_id, ",")))
}
# Get HUC12 sf object (queries USGS WBD if layer not provided)
if (is.null(huc_layer)) {
huc_layer <- nhdplusTools::get_huc(id = huc_id, type = "huc12")
}
if (is.null(huc_layer) || nrow(huc_layer) == 0) {
stop("Invalid HUC12 ID or could not retrieve watershed data from USGS.")
}
# Track individual component HUCs if multi-HUC
individual_hucs <- huc_layer
actual_huc_ids <- if ("huc12" %in% names(huc_layer)) unique(huc_layer$huc12) else huc_id
if (!is.null(feature_name)) {
if (!requireNamespace("osmdata", quietly = TRUE)) {
stop("The 'osmdata' package is required to filter geographic features by name. Install it using install.packages('osmdata')")
}
# Try querying OpenStreetMap nominatim for the geographic feature boundary
feature_geom <- tryCatch({
osmdata::getbb(feature_name, format_out = "sf_polygon", limit = 1)
}, error = function(e) {
warning(paste("osmdata could not find a valid polygon for feature:", feature_name, "- Generating whole HUC region instead."))
return(NULL)
})
if (!is.null(feature_geom)) {
# Formats may vary (list containing polygon/multipolygon, or an sf object directly)
feature_sf <- NULL
if (inherits(feature_geom, "sf") || inherits(feature_geom, "sfc")) {
feature_sf <- feature_geom
} else if (is.list(feature_geom)) {
if (!is.null(feature_geom$multipolygon)) {
feature_sf <- feature_geom$multipolygon
} else if (!is.null(feature_geom$polygon)) {
feature_sf <- feature_geom$polygon
} else if (length(feature_geom) > 0) {
feature_sf <- feature_geom[[1]]
}
}
if (!is.null(feature_sf)) {
# Project to HUC's CRS and spatially intersect to restrict bounds
feature_sf <- sf::st_transform(feature_sf, sf::st_crs(huc_layer))
clipped_layer <- suppressWarnings(sf::st_intersection(huc_layer, feature_sf))
# If intersection yields valid geometries, use it; otherwise fallback to full HUC
if (nrow(clipped_layer) > 0) {
huc_layer <- clipped_layer
} else {
warning(paste("Feature", feature_name, "does not overlap with specified HUC region - Generating whole HUC region instead."))
}
} else {
warning(paste("osmdata could not extract a valid polygon for feature:", feature_name, "- Generating whole HUC region instead."))
}
}
}
# Aggregate multiple HUC12 polygons into a single combined region geometry via sf::st_union
unified_layer <- if (nrow(huc_layer) > 1) {
suppressWarnings(sf::st_union(huc_layer))
} else {
huc_layer
}
# Calculate geographic centroid of the final geometry
centroid <- suppressWarnings(sf::st_centroid(sf::st_geometry(unified_layer)))
coords <- sf::st_coordinates(centroid)
list(
huc_id = actual_huc_ids,
feature_name = feature_name,
lon = as.numeric(coords[1, "X"]),
lat = as.numeric(coords[1, "Y"]),
individual_hucs = individual_hucs,
layer = unified_layer
)
}
safe_st_intersects <- function(x, y) {
tryCatch(
sf::st_intersects(x, y),
error = function(e) {
old_s2 <- sf::sf_use_s2(FALSE)
on.exit(sf::sf_use_s2(old_s2), add = TRUE)
x_val <- tryCatch(sf::st_make_valid(x), error = function(e2) x)
y_val <- tryCatch(sf::st_make_valid(y), error = function(e2) y)
sf::st_intersects(x_val, y_val)
}
)
}
add_watershed_hex_overlay <- function(huc_info, hex_diameter = 0.01) {
huc_layer <- huc_info$layer
# Use sf::st_make_grid with square = FALSE to mathematically build a spatial hex mesh atop the bounding box
hex_mesh <- sf::st_make_grid(huc_layer, square = FALSE, cellsize = c(hex_diameter, hex_diameter))
# Filter the generated mesh to only retain hexagons crossing the actual geographical feature
# (lengths > 0 signifies the hexagon touches the island geometry)
hex_overlay <- hex_mesh[lengths(safe_st_intersects(hex_mesh, huc_layer)) > 0]
huc_info$hex_overlay <- hex_overlay
huc_info$hex_diameter <- hex_diameter
class(huc_info) <- "watershed_hex_overlay"
return(huc_info)
}
autoplot.watershed_hex_overlay <- function(object, ...) {
huc_str <- if (length(object$huc_id) > 1) {
paste0(length(object$huc_id), " Combined HUC12s")
} else {
paste("Watershed:", object$huc_id)
}
title_txt <- paste("Geographic Hexagonal Grid (", huc_str, ")",
"\nHexagon Extent Diameter:", object$hex_diameter)
if (!is.null(object$feature_name) && object$feature_name != "") {
title_txt <- paste0(title_txt, " - Restricted to: ", object$feature_name)
}
p <- ggplot2::ggplot()
# Display individual component HUC boundaries if multi-HUC
if (!is.null(object$individual_hucs) && nrow(object$individual_hucs) > 1) {
p <- p + ggplot2::geom_sf(data = object$individual_hucs, fill = NA, color = "purple", linetype = "dashed", linewidth = 0.4)
}
p +
# Represents underlying restricted outline
ggplot2::geom_sf(data = object$layer, fill = "lightblue", alpha = 0.3, color = "blue", linewidth = 0.7) +
# Overlay our spatial hexagons
ggplot2::geom_sf(data = object$hex_overlay, fill = NA, color = "darkred", linewidth = 0.7) +
ggplot2::theme_minimal() +
ggplot2::ggtitle(title_txt) +
ggplot2::labs(x = "Longitude", y = "Latitude")
}
discover_watershed_features <- function(huc_id, feature_types = c("natural", "waterway", "leisure")) {
# Get HUC12 sf object to establish the tight bounding limit
huc_layer <- nhdplusTools::get_huc(id = huc_id, type = "huc12")
if (is.null(huc_layer) || nrow(huc_layer) == 0) {
stop("Invalid HUC12 ID or could not retrieve watershed data from USGS.")
}
if (!requireNamespace("osmdata", quietly = TRUE)) {
stop("The 'osmdata' package is required to execute dynamic feature discovery. Install it using install.packages('osmdata')")
}
# Map coordinates to WGS84 for the public Overpass API interface
bbox <- sf::st_bbox(sf::st_transform(huc_layer, 4326))
bbox_str <- paste(bbox["ymin"], bbox["xmin"], bbox["ymax"], bbox["xmax"], sep = ",")
# Construct raw Overpass QL to execute a UNION (OR) query natively
ql_union_body <- paste(
sapply(feature_types, function(key) {
paste0(
" node[\"", key, "\"](", bbox_str, ");\n",
" way[\"", key, "\"](", bbox_str, ");\n",
" relation[\"", key, "\"](", bbox_str, ");\n"
)
}), collapse = ""
)
ql_query <- paste0(
"[out:xml][timeout:60];\n(\n",
ql_union_body,
");\n",
"out body;\n>;\nout skel qt;\n"
)
discovered_names <- c()
# Temporarily reroute to the lz4 alternative Overpass Mirror to dodge IP bans on main branch
old_url <- osmdata::get_overpass_url()
osmdata::set_overpass_url("https://lz4.overpass-api.de/api/interpreter")
tryCatch({
# Temporarily disable standard S2 spherical tracking to bypass OSM topological boundary crashes
old_s2 <- sf::sf_use_s2()
sf::sf_use_s2(FALSE)
# Send the raw, unified union query into the osmdata parser
osm_res <- osmdata::osmdata_sf(ql_query)
# Extract polygons (Lakes, Parks)
if (!is.null(osm_res$osm_polygons) && "name" %in% colnames(osm_res$osm_polygons)) {
poly_sf <- osm_res$osm_polygons[!is.na(osm_res$osm_polygons$name), ]
if (nrow(poly_sf) > 0) {
poly_sf <- sf::st_transform(poly_sf, sf::st_crs(huc_layer))
poly_sf <- suppressWarnings(sf::st_make_valid(poly_sf))
inter <- suppressWarnings(lengths(sf::st_intersects(poly_sf, huc_layer))) > 0
discovered_names <- c(discovered_names, poly_sf$name[inter])
}
}
# Extract multipolygons (Great Lakes, large reserves)
if (!is.null(osm_res$osm_multipolygons) && "name" %in% colnames(osm_res$osm_multipolygons)) {
mpoly_sf <- osm_res$osm_multipolygons[!is.na(osm_res$osm_multipolygons$name), ]
if (nrow(mpoly_sf) > 0) {
mpoly_sf <- sf::st_transform(mpoly_sf, sf::st_crs(huc_layer))
mpoly_sf <- suppressWarnings(sf::st_make_valid(mpoly_sf))
inter <- suppressWarnings(lengths(sf::st_intersects(mpoly_sf, huc_layer))) > 0
discovered_names <- c(discovered_names, mpoly_sf$name[inter])
}
}
# Extract linestrings (Dams, Rivers)
if (!is.null(osm_res$osm_lines) && "name" %in% colnames(osm_res$osm_lines)) {
line_sf <- osm_res$osm_lines[!is.na(osm_res$osm_lines$name), ]
if (nrow(line_sf) > 0) {
line_sf <- sf::st_transform(line_sf, sf::st_crs(huc_layer))
line_sf <- suppressWarnings(sf::st_make_valid(line_sf))
inter <- suppressWarnings(lengths(sf::st_intersects(line_sf, huc_layer))) > 0
discovered_names <- c(discovered_names, line_sf$name[inter])
}
}
# Reset standard mapping state
sf::sf_use_s2(old_s2)
}, error = function(e) {
if (exists("old_s2")) sf::sf_use_s2(old_s2)
warning(paste("Raw OSM Extraction timeout or failure:", e$message))
})
# Return sanitized unique list
clean_names <- unique(discovered_names)
clean_names <- sort(clean_names[!is.na(clean_names) & trimws(clean_names) != ""])
# Reset the API routing mirror
osmdata::set_overpass_url(old_url)
return(clean_names)
}
# --- Source: isle_royale_sim.R ---
init_isle_royale_sim <- function(year = 1980,
n_moose = NULL,
n_wolves = NULL,
hex_diameter = 0.01,
datafile = "",
features_rds = NULL,
landmarks_rds = NULL) {
# Load historical benchmark time series data
csv_path <- system.file("doc/isle_royale/wolf_moose.csv", package = "ewing")
if (csv_path == "" || !file.exists(csv_path)) csv_path <- "inst/doc/isle_royale/wolf_moose.csv"
hist_data <- NULL
if (file.exists(csv_path)) {
hist_data <- utils::read.csv(csv_path, stringsAsFactors = FALSE)
}
# Lookup baseline population counts for start year if not explicitly provided
if (!is.null(hist_data) && year %in% hist_data$Year) {
row_match <- hist_data[hist_data$Year == year, ]
if (is.null(n_moose)) n_moose <- as.numeric(row_match$Moose[1])
if (is.null(n_wolves)) n_wolves <- as.numeric(row_match$Wolves[1])
}
if (is.null(n_moose)) n_moose <- 664
if (is.null(n_wolves)) n_wolves <- 50
# 1. Initialize Base Spatial Geography & Habitat Overlay
huc_info <- get_watershed("041800000101", feature_name = "Isle Royale")
hex_obj <- add_watershed_hex_overlay(huc_info, hex_diameter = hex_diameter)
habitat_overlay <- add_habitat_hex_overlay(hex_obj, features = features_rds, landmarks = landmarks_rds)
# 2. Setup Configuration Data Directory
if (datafile == "") {
pkg_dir <- system.file("extdata/isle_royale", package = "ewing")
if (pkg_dir != "" && dir.exists(pkg_dir)) {
datafile <- pkg_dir
} else if (dir.exists("inst/extdata/isle_royale")) {
datafile <- "inst/extdata/isle_royale"
}
}
# 3. Initialize ewing Community Core
community <- NULL
tryCatch({
community <- init.simulation(package = "ewing", count = c(n_moose, n_wolves), datafile = datafile, messages = FALSE)
}, error = function(e) {
# Soft fallback if ewing initialization runs in isolated test environment
community <<- list(pop = list())
})
# 4. Spatially Sample Initial Positions Weighted by Substrate Habitat Score
hex_sf <- habitat_overlay$hex_habitat_sf
probs <- hex_sf$habitat_score / sum(hex_sf$habitat_score)
sampled_moose_hex <- sample(seq_len(nrow(hex_sf)), size = n_moose, replace = TRUE, prob = probs)
sampled_wolf_hex <- sample(seq_len(nrow(hex_sf)), size = n_wolves, replace = TRUE, prob = probs)
# Extract hex centroids as point locations
centroids <- sf::st_centroid(sf::st_geometry(hex_sf))
moose_pts <- centroids[sampled_moose_hex]
wolf_pts <- centroids[sampled_wolf_hex]
moose_coords <- sf::st_coordinates(moose_pts)
wolf_coords <- sf::st_coordinates(wolf_pts)
moose_df <- data.frame(
id = paste0("M", seq_len(n_moose)),
species = "Moose",
ageclass = sample(c("calf", "yearling", "adult", "senior"), size = n_moose, replace = TRUE, prob = c(0.15, 0.15, 0.55, 0.15)),
hex_id = sampled_moose_hex,
lon = moose_coords[, 1],
lat = moose_coords[, 2],
stringsAsFactors = FALSE
)
wolf_df <- data.frame(
id = paste0("W", seq_len(n_wolves)),
species = "Wolf",
ageclass = sample(c("pup", "subadult", "adult"), size = n_wolves, replace = TRUE, prob = c(0.20, 0.25, 0.55)),
hex_id = sampled_wolf_hex,
lon = wolf_coords[, 1],
lat = wolf_coords[, 2],
stringsAsFactors = FALSE
)
# Initialize Population History Tallies (for Dist Plots)
m_counts <- table(factor(moose_df$ageclass, levels = c("calf", "yearling", "adult", "senior")))
w_counts <- table(factor(wolf_df$ageclass, levels = c("pup", "subadult", "adult")))
hist_df <- data.frame(
step = 0,
time = 0,
Species = c(rep("moose", 4), rep("wolf", 3)),
State = c(names(m_counts), names(w_counts)),
Type = "ageclass",
Count = c(as.numeric(m_counts), as.numeric(w_counts)),
stringsAsFactors = FALSE
)
res <- list(
community = community,
habitat_overlay = habitat_overlay,
start_year = year,
moose_pop = moose_df,
wolf_pop = wolf_df,
historical_data = hist_data,
history = hist_df,
nstep = 0
)
class(res) <- "isle_royale_sim"
return(res)
}
run_isle_royale_sim <- function(sim_obj, nstep = 1000, refresh = 10, ...) {
if (!inherits(sim_obj, "isle_royale_sim")) {
stop("Input must be an object of class 'isle_royale_sim'.")
}
hex_sf <- sim_obj$habitat_overlay$hex_habitat_sf
n_hex <- nrow(hex_sf)
scores <- hex_sf$habitat_score
# Precompute neighbor adjacency list and centroid coordinates
adj_list <- suppressWarnings(sf::st_touches(hex_sf))
centroids <- sf::st_centroid(sf::st_geometry(hex_sf))
cent_coords <- sf::st_coordinates(centroids)
# Helper to step individuals to adjacent hexes weighted by habitat_score
move_pop <- function(pop_df, move_prob = 0.5) {
if (is.null(pop_df) || nrow(pop_df) == 0) return(pop_df)
n <- nrow(pop_df)
new_hex <- pop_df$hex_id
for (i in seq_len(n)) {
if (stats::runif(1) < move_prob) {
curr_h <- pop_df$hex_id[i]
nbrs <- adj_list[[curr_h]]
cands <- c(curr_h, nbrs)
w <- scores[cands]
p <- w / sum(w)
new_h <- sample(cands, 1, prob = p)
new_hex[i] <- new_h
}
}
pop_df$hex_id <- new_hex
# Calculate coords with small spatial jitter
hex_diam <- sim_obj$habitat_overlay$hex_diameter
if (is.null(hex_diam)) hex_diam <- 0.01
jitter_x <- stats::rnorm(n, mean = 0, sd = hex_diam * 0.15)
jitter_y <- stats::rnorm(n, mean = 0, sd = hex_diam * 0.15)
pop_df$lon <- cent_coords[new_hex, 1] + jitter_x
pop_df$lat <- cent_coords[new_hex, 2] + jitter_y
return(pop_df)
}
# Helper for dynamic life stage progression (aging, births, natural deaths)
update_demographics <- function(moose_df, wolf_df) {
# 1. Moose Life Stage Transitions & Mortality
if (!is.null(moose_df) && nrow(moose_df) > 0) {
n_m <- nrow(moose_df)
# Aging
calf_idx <- which(moose_df$ageclass == "calf")
if (length(calf_idx) > 0) {
trans <- calf_idx[stats::runif(length(calf_idx)) < 0.03]
if (length(trans) > 0) moose_df$ageclass[trans] <- "yearling"
}
yearling_idx <- which(moose_df$ageclass == "yearling")
if (length(yearling_idx) > 0) {
trans <- yearling_idx[stats::runif(length(yearling_idx)) < 0.03]
if (length(trans) > 0) moose_df$ageclass[trans] <- "adult"
}
adult_idx <- which(moose_df$ageclass == "adult")
if (length(adult_idx) > 0) {
trans <- adult_idx[stats::runif(length(adult_idx)) < 0.015]
if (length(trans) > 0) moose_df$ageclass[trans] <- "senior"
# Reproduction (Calf births from adults)
n_births <- sum(stats::runif(length(adult_idx)) < 0.012)
if (n_births > 0) {
parent_hexes <- sample(moose_df$hex_id[adult_idx], n_births, replace = TRUE)
max_id <- suppressWarnings(max(as.numeric(gsub("[^0-9]", "", moose_df$id)), na.rm = TRUE))
if (!is.finite(max_id)) max_id <- nrow(moose_df)
new_calves <- data.frame(
id = paste0("M", max_id + seq_len(n_births)),
species = "Moose",
ageclass = "calf",
hex_id = parent_hexes,
lon = cent_coords[parent_hexes, 1] + stats::rnorm(n_births, 0, 0.001),
lat = cent_coords[parent_hexes, 2] + stats::rnorm(n_births, 0, 0.001),
stringsAsFactors = FALSE
)
moose_df <- rbind(moose_df, new_calves)
}
}
# Natural mortality for seniors
senior_idx <- which(moose_df$ageclass == "senior")
if (length(senior_idx) > 0) {
deaths <- senior_idx[stats::runif(length(senior_idx)) < 0.025]
if (length(deaths) > 0) moose_df <- moose_df[-deaths, ]
}
}
# 2. Wolf Life Stage Transitions & Mortality
if (!is.null(wolf_df) && nrow(wolf_df) > 0) {
pup_idx <- which(wolf_df$ageclass == "pup")
if (length(pup_idx) > 0) {
trans <- pup_idx[stats::runif(length(pup_idx)) < 0.04]
if (length(trans) > 0) wolf_df$ageclass[trans] <- "subadult"
}
sub_idx <- which(wolf_df$ageclass == "subadult")
if (length(sub_idx) > 0) {
trans <- sub_idx[stats::runif(length(sub_idx)) < 0.04]
if (length(trans) > 0) wolf_df$ageclass[trans] <- "adult"
}
w_adult_idx <- which(wolf_df$ageclass == "adult")
if (length(w_adult_idx) > 0) {
# Reproduction (Pup births)
n_w_births <- sum(stats::runif(length(w_adult_idx)) < 0.015)
if (n_w_births > 0) {
parent_hexes <- sample(wolf_df$hex_id[w_adult_idx], n_w_births, replace = TRUE)
max_w_id <- suppressWarnings(max(as.numeric(gsub("[^0-9]", "", wolf_df$id)), na.rm = TRUE))
if (!is.finite(max_w_id)) max_w_id <- nrow(wolf_df)
new_pups <- data.frame(
id = paste0("W", max_w_id + seq_len(n_w_births)),
species = "Wolf",
ageclass = "pup",
hex_id = parent_hexes,
lon = cent_coords[parent_hexes, 1] + stats::rnorm(n_w_births, 0, 0.001),
lat = cent_coords[parent_hexes, 2] + stats::rnorm(n_w_births, 0, 0.001),
stringsAsFactors = FALSE
)
wolf_df <- rbind(wolf_df, new_pups)
}
# Natural adult wolf mortality
w_deaths <- w_adult_idx[stats::runif(length(w_adult_idx)) < 0.01]
if (length(w_deaths) > 0) wolf_df <- wolf_df[-w_deaths, ]
}
}
list(moose = moose_df, wolf = wolf_df)
}
# Micro-step execution loop (recording history continuously as changes happen)
step_chunk <- max(1, min(10, refresh))
steps_left <- nstep
while (steps_left > 0) {
chunk <- min(step_chunk, steps_left)
steps_left <- steps_left - chunk
# 1. Underlying ewing simulation engine step if present
if (!is.null(sim_obj$community) && inherits(sim_obj$community, "ewing")) {
sim_obj$community <- tryCatch({
future.events(sim_obj$community, nstep = chunk, refresh = chunk, plotit = FALSE, ...)
}, error = function(e) sim_obj$community)
}
# 2. Update Spatial Movement
sim_obj$moose_pop <- move_pop(sim_obj$moose_pop, move_prob = 0.6)
sim_obj$wolf_pop <- move_pop(sim_obj$wolf_pop, move_prob = 0.75)
# 3. Process Wolf Predation on Vulnerable Moose
if (!is.null(sim_obj$wolf_pop) && !is.null(sim_obj$moose_pop) && nrow(sim_obj$wolf_pop) > 0 && nrow(sim_obj$moose_pop) > 0) {
wolf_hexes <- unique(sim_obj$wolf_pop$hex_id)
vulnerable_idx <- which(sim_obj$moose_pop$hex_id %in% wolf_hexes & sim_obj$moose_pop$ageclass %in% c("calf", "senior"))
if (length(vulnerable_idx) > 0) {
n_predated <- min(length(vulnerable_idx), max(1, round(nrow(sim_obj$wolf_pop) * 0.02)))
pred_remove <- sample(vulnerable_idx, n_predated)
sim_obj$moose_pop <- sim_obj$moose_pop[-pred_remove, ]
}
}
# 4. Process Life Stage Transitions & Demographics
demog_res <- update_demographics(sim_obj$moose_pop, sim_obj$wolf_pop)
sim_obj$moose_pop <- demog_res$moose
sim_obj$wolf_pop <- demog_res$wolf
sim_obj$nstep <- sim_obj$nstep + chunk
# 5. Append Micro-step History Tally (for continuous Dist Plot step curves)
m_curr <- table(factor(sim_obj$moose_pop$ageclass, levels = c("calf", "yearling", "adult", "senior")))
w_curr <- table(factor(sim_obj$wolf_pop$ageclass, levels = c("pup", "subadult", "adult")))
step_hist <- data.frame(
step = sim_obj$nstep,
time = sim_obj$nstep,
Species = c(rep("moose", 4), rep("wolf", 3)),
State = c(names(m_curr), names(w_curr)),
Type = "ageclass",
Count = c(as.numeric(m_curr), as.numeric(w_curr)),
stringsAsFactors = FALSE
)
sim_obj$history <- rbind(sim_obj$history, step_hist)
}
return(sim_obj)
}
ggplot_isle_royale_sim <- function(x, ...) {
# 1. Spatial Landscape Plot
p_map <- autoplot(x$habitat_overlay, show_landmarks = TRUE)
# Prepare spatial organism data frames
moose_sf <- sf::st_as_sf(x$moose_pop, coords = c("lon", "lat"), crs = sf::st_crs(x$habitat_overlay$layer))
wolf_sf <- sf::st_as_sf(x$wolf_pop, coords = c("lon", "lat"), crs = sf::st_crs(x$habitat_overlay$layer))
p_map <- p_map +
ggplot2::geom_sf(data = moose_sf, color = "#27ae60", shape = 21, fill = NA, stroke = 1.0, size = 0.8, alpha = 0.85) +
ggplot2::geom_sf(data = wolf_sf, color = "#e74c3c", shape = 21, fill = NA, stroke = 1.4, size = 1.5, alpha = 0.95) +
ggplot2::ggtitle(paste0("Isle Royale Wolf-Moose Spatial Landscape (Year ", x$start_year, ")"))
# 2. Historical Benchmark Trajectory Plot
if (!is.null(x$historical_data) && nrow(x$historical_data) > 0) {
df_hist <- x$historical_data
p_hist <- ggplot2::ggplot(df_hist, ggplot2::aes(x = .data$Year)) +
ggplot2::geom_line(ggplot2::aes(y = .data$Moose, color = "Historical Moose"), linewidth = 1.0) +
ggplot2::geom_line(ggplot2::aes(y = .data$Wolves * 40, color = "Historical Wolves (x40)"), linewidth = 1.0, linetype = "dashed") +
ggplot2::scale_color_manual(
name = "Empirical Time Series",
values = c("Historical Moose" = "#27ae60", "Historical Wolves (x40)" = "#e74c3c")
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Isle Royale Empirical Census Benchmarks (1980 - 2019)",
x = "Year", y = "Moose Abundance (Wolves x40)"
)
if (requireNamespace("cowplot", quietly = TRUE)) {
return(cowplot::plot_grid(p_map, p_hist, ncol = 1, rel_heights = c(1.2, 1)))
}
}
return(p_map)
}
plot.isle_royale_sim <- function(x, ...) {
print(ggplot_isle_royale_sim(x, ...))
}
# --- Source: IsleRoyaleApp.R ---
IsleRoyaleApp <- function(title = "Isle Royale Wolf-Moose Simulation Platform") {
ui <- bslib::page_sidebar(
title = title,
sidebar = bslib::sidebar(
width = 340,
IsleRoyaleInput("isleroyale")
),
IsleRoyaleOutput("isleroyale")
)
server <- function(input, output, session) {
IsleRoyaleServer("isleroyale")
}
shiny::shinyApp(ui = ui, server = server)
}
IsleRoyaleInput <- function(id) {
ns <- shiny::NS(id)
# Fetch available years from wolf_moose.csv if available
csv_path <- system.file("doc/isle_royale/wolf_moose.csv", package = "ewing")
if (csv_path == "" || !file.exists(csv_path)) csv_path <- "inst/doc/isle_royale/wolf_moose.csv"
years <- 1980:2019
if (file.exists(csv_path)) {
df <- tryCatch(utils::read.csv(csv_path), error = function(e) NULL)
if (!is.null(df) && "Year" %in% names(df)) years <- df$Year
}
shiny::tagList(
shiny::h4("Isle Royale Controls", style = "font-size: 1.0rem; font-weight: 600;"),
shiny::selectInput(ns("start_year"), "Historical Baseline Year:", choices = years, selected = 1980),
shiny::sliderInput(ns("n_moose"), "Initial Moose Count:", min = 50, max = 3000, value = 664, step = 50),
shiny::sliderInput(ns("n_wolves"), "Initial Wolf Count:", min = 0, max = 60, value = 50, step = 2),
shiny::sliderInput(ns("hex_diameter"), "Hexagon Extent Diameter (Degrees):", min = 0.005, max = 0.03, value = 0.01, step = 0.001),
# Conditional Spatial Overlay Controls (shown ONLY on Substrate Plot & Census Benchmarks)
shiny::conditionalPanel(
condition = sprintf("input['%s'] == 'Substrate Plot' || input['%s'] == 'Census Benchmarks'", ns("tabset"), ns("tabset")),
shiny::div(style = "border-top: 1px solid rgba(0,0,0,0.1); margin: 6px 0;"),
shiny::h4("Map Overlay Options", style = "font-size: 0.9rem; font-weight: 600; margin-bottom: 4px;"),
shiny::checkboxInput(ns("show_habitat"), "Overlay Moose Habitat Features", value = TRUE),
shiny::checkboxInput(ns("show_landmarks"), "Show Moose Sighting Landmarks", value = TRUE)
),
# Conditional Substrate Plot Axis Units
shiny::conditionalPanel(
condition = sprintf("input['%s'] == 'Substrate Plot'", ns("tabset")),
axisUnitInput(ns("substrate_axis"), time_label = "Days")
),
# Conditional Age Classes Display Controls (shown ONLY on Age Classes tab)
shiny::conditionalPanel(
condition = sprintf("input['%s'] == 'Age Classes'", ns("tabset")),
shiny::div(style = "border-top: 1px solid rgba(0,0,0,0.1); margin: 6px 0;"),
shiny::h4("Age Classes Options", style = "font-size: 0.9rem; font-weight: 600; margin-bottom: 4px;"),
ageClassControlInput(ns("age_ctrls"), time_label = "Days")
),
shiny::HTML("<hr style='margin: 10px 0;'/>"),
shiny::actionButton(ns("reset_sim"), "Reset Simulation", class = "btn-warning", style = "margin-bottom: 6px; width: 100%;"),
step_size_slider(ns("step_size"), "Steps per click:", selected = 200),
shiny::actionButton(ns("step_sim"), "Run Simulation Steps", class = "btn-primary", style = "width: 100%;"),
shiny::br(), shiny::br(),
shiny::uiOutput(ns("status"))
)
}
IsleRoyaleOutput <- function(id) {
ns <- shiny::NS(id)
shiny::tagList(
shiny::tabsetPanel(
id = ns("tabset"),
type = "tabs",
shiny::tabPanel(
"Substrate Plot",
shiny::plotOutput(ns("substrate_plot"), height = "650px")
),
shiny::tabPanel(
"Age Classes",
distPlotOutput(ns("dist_plot"))
),
shiny::tabPanel(
"Census Benchmarks",
shiny::plotOutput(ns("autoplot"), height = "650px")
),
shiny::tabPanel(
"Live Demographics",
shiny::br(),
shiny::h4("Live Organism Population Summary (Updates Live on Simulation Step)"),
shiny::tableOutput(ns("summary_table"))
),
shiny::tabPanel(
"Input Data",
shiny::br(),
inputAppInput(ns("input_data")),
inputAppOutput(ns("input_data"))
)
)
)
}
IsleRoyaleServer <- function(id) {
shiny::moduleServer(id, function(input, output, session) {
ns <- session$ns
status_msg <- shiny::reactiveVal("")
sim_state <- shiny::reactiveVal(NULL)
sub_x_var <- axisUnitServer("substrate_axis")
age_ctrls <- ageClassControlServer("age_ctrls")
# Auto-update initial population slider defaults when baseline year changes
shiny::observeEvent(input$start_year, {
yr <- as.numeric(input$start_year)
csv_path <- system.file("doc/isle_royale/wolf_moose.csv", package = "ewing")
if (csv_path == "" || !file.exists(csv_path)) csv_path <- "inst/doc/isle_royale/wolf_moose.csv"
if (file.exists(csv_path)) {
df <- tryCatch(utils::read.csv(csv_path), error = function(e) NULL)
if (!is.null(df) && "Year" %in% names(df) && yr %in% df$Year) {
row_match <- df[df$Year == yr, ]
shiny::updateSliderInput(session, "n_moose", value = as.numeric(row_match$Moose[1]))
shiny::updateSliderInput(session, "n_wolves", value = as.numeric(row_match$Wolves[1]))
}
}
})
# Initialize / Reset Simulation with Default Initial Steps (Matching sysetholApp behavior)
shiny::observeEvent(list(input$reset_sim, input$start_year), {
yr <- as.numeric(input$start_year)
nm <- input$n_moose
nw <- input$n_wolves
hd <- input$hex_diameter
shiny::req(yr, nm, nw, hd)
sim <- init_isle_royale_sim(
year = yr,
n_moose = nm,
n_wolves = nw,
hex_diameter = hd
)
# Run initial default simulation steps (matching sysetholApp.R)
steps <- parse_step_size(input$step_size)
if (is.null(steps) || steps <= 0) steps <- 200
sim <- run_isle_royale_sim(sim, nstep = steps)
sim_state(sim)
status_msg(paste0("<div style='color:green;'><b>Simulation Initialized:</b> Year ", yr, " with ", nm, " Moose and ", nw, " Wolves (Executed ", steps, " initial steps).</div>"))
}, ignoreNULL = FALSE)
# Step Simulation Execution
shiny::observeEvent(input$step_sim, {
sim <- sim_state()
shiny::req(sim)
steps <- parse_step_size(input$step_size)
if (is.null(steps) || steps <= 0) steps <- 200
updated_sim <- run_isle_royale_sim(sim, nstep = steps)
sim_state(updated_sim)
status_msg(paste0("<div style='color:blue;'><b>Executed ", steps, " steps:</b> Total Steps = ", updated_sim$nstep, "</div>"))
})
# Render Offline Isle Royale Substrate Plot (0 API calls!)
output$substrate_plot <- shiny::renderPlot({
sim <- sim_state()
shiny::req(sim)
ewing_substrate(sim, x_var = sub_x_var())
})
# Compose Dist Plot Module for Age Classes
distPlotServer("dist_plot", simres = sim_state, x_var = age_ctrls$x_var, total = age_ctrls$total, norm = age_ctrls$norm)
# Render ggplot autoplot (Dual-panel benchmark plot)
output$autoplot <- shiny::renderPlot({
sim <- sim_state()
shiny::req(sim)
ggplot_isle_royale_sim(sim)
})
# Render Summary Table
output$summary_table <- shiny::renderTable({
sim <- sim_state()
shiny::req(sim)
moose_counts <- table(sim$moose_pop$ageclass)
wolf_counts <- table(sim$wolf_pop$ageclass)
df_moose <- data.frame(
Species = "Moose",
AgeClass = names(moose_counts),
Count = as.numeric(moose_counts),
stringsAsFactors = FALSE
)
df_wolf <- data.frame(
Species = "Wolf",
AgeClass = names(wolf_counts),
Count = as.numeric(wolf_counts),
stringsAsFactors = FALSE
)
rbind(df_moose, df_wolf)
})
# Compose Input Data App Module for Input Table Inspection
pkg_dir <- system.file("extdata/isle_royale", package = "ewing")
if (pkg_dir == "" || !dir.exists(pkg_dir)) pkg_dir <- "inst/extdata/isle_royale"
inputAppServer("input_data", simres = sim_state, datafile = shiny::reactiveVal(pkg_dir))
# Output status message
output$status <- shiny::renderUI({
shiny::HTML(status_msg())
})
})
}
# --- Launch Application ---
IsleRoyaleApp()Launch the interactive application natively in R using exported package functions:
library(ewing)
# Launch interactive offline Shiny app
IsleRoyaleApp()Below is an executable simulation workflow displaying spatial individual organism distributions on the Isle Royale substrate map, step-by-step demographic tallies, and historical 40-year empirical census benchmarks (1980–2019).
library(ewing)
library(ggplot2)
# 1. Initialize Isle Royale spatial simulation with 1980 baseline counts (664 Moose, 50 Wolves)
sim <- init_isle_royale_sim(year = 1980)
# 2. Run simulation steps (e.g. 200 days)
sim <- run_isle_royale_sim(sim, nstep = 200)
# 3. Visualize spatial landscape & benchmark trajectories
ggplot_isle_royale_sim(sim)
# Plot side-by-side per-species age class dynamics over time (days)
autoplot(ewing_ageclass(sim), x_var = "time")
sf layers (inst/extdata/isle_royale/) with 0 API calls.show_habitat, show_landmarks, norm, total, Steps vs Days) strictly when their target tab is active.step_controls, distPlotApp, inputApp, and substrateApp.