#' accumulate
#'
#' bootstrap wrapper for recorder function
#'
#' @param n is the number of accumulated arrivals
#' @param ri is the rate in, a float
#' @param ro is the rate out, a float
#' @param ia is a vector of inherited ages, integers
#' @param it is the number of iterations of the bootstrap, an integer
#' @return the median deposit ages of `n` accumulations after `it` bootstraps
#' @export
accumulate <- function(n, ri, ro, ia, it = 100) {
ar <- array(0, c(n, it))
for (i in 1:it) {
b <- recorder(n, ri, ro, ia)
ar[1:length(b), i] <- b
}
vec <- apply(ar, 1, median)
return(vec)
}
#' accumulate1
#'
#' bootstrap wrapper for recorder function
#' adds minimum age
#'
#' @param n is the number of accumulated arrivals
#' @param ri is the rate in, a float
#' @param ro is the rate out, a float
#' @param ia is a vector of inherited ages, integers
#' @param mn is the minimum age
#' @param it is the number of iterations of the bootstrap, an integer
#' @return the median deposit ages of `n` accumulations after `it` bootstraps
#' @seealso recorder1
#' @export
accumulate1 <- function(n, ri, ro, ia, mn, it = 100) {
ar <- array(0, c(n, it))
for (i in 1:it) {
b <- recorder1(n, ri, ro, ia, mn)
ar[1:length(b), i] <- b
}
vec <- apply(ar, 1, median)
return(vec)
}
#' accumulater
#'
#' bootstrap wrapper for accumulate function
#'
#' @param n is the number of accumulated arrivals
#' @param ti_rng is the range of rates in, a vector of floats
#' @param to_rng is the range of rates out, a vector of floats
#' @param ia is a vector of inherited ages, integers
#' @param obs is the vector of observed estimated mean charcoal ages
#' @param batch is the batch size of random draws
#' @param it is the number of iterations of the bootstrap, an integer
#' @return the median deposit ages of `n` accumulations after `it` bootstraps
#' @seealso accumulate
#' @seealso recorder
#' @export
accumulater <- function(n, ti_rng, to_rng, ia, obs, batch = 10, it = 100) {
go <- runif(batch, min(to_rng), max(to_rng))
gi <- runif(batch, min(ti_rng), max(ti_rng))
res <- mcmapply(
function(x, y, a, b, c) accumulate(a, x, y, b, c),
x = gi, y = go,
MoreArgs = list(a = n, b = ia, c = it)
)
gof <- apply(res, 2, function(x) cdf_gof(x, obs))
gof <- matrix(unlist(gof), ncol = 2, byrow = TRUE)
df <- data.frame(
ti = gi,
to = go,
ks = gof[, 1],
kp = gof[, 2]
)
return(df)
}
#' accumulater1
#'
#' bootstrap wrapper for accumulate function
#' adds minimum age
#'
#' @param n is the number of accumulated arrivals
#' @param ti_rng is the range of rates in, a vector of floats
#' @param to_rng is the range of rates out, a vector of floats
#' @param ia is a vector of inherited ages, integers
#' @param obs is the vector of observed estimated mean charcoal ages
#' @param batch is the batch size of random draws
#' @param it is the number of iterations of the bootstrap, an integer
#' @return the median deposit ages of `n` accumulations after `it` bootstraps
#' @seealso accumulate
#' @seealso recorder
#' @export
accumulater1 <- function(n, ti_rng, to_rng, ia, obs, batch = 10, it = 100) {
go <- runif(batch, min(to_rng), max(to_rng))
gi <- runif(batch, min(ti_rng), max(ti_rng))
mn <- min(obs)
res <- mcmapply(
function(x, y, a, b, c, d) accumulate1(a, x, y, b, c, d),
x = gi, y = go,
MoreArgs = list(a = n, b = ia, c = mn, d = it)
)
gof <- apply(res, 2, function(x) cdf_gof(x, obs))
gof <- matrix(unlist(gof), ncol = 2, byrow = TRUE)
df <- data.frame(
ti = gi,
to = go,
ks = gof[, 1],
kp = gof[, 2]
)
return(df)
}
#' boot_accum
#'
#' bootstrap utility for estimating confidence intervals
#' wraps accumulater()
#'
#' @param n is the number of accumulated arrivals
#' @param ti_rng is the range of rates in, a vector of floats
#' @param to_rng is the range of rates out, a vector of floats
#' @param ia is a vector of inherited ages, integers
#' @param obs is the vector of observed estimated mean charcoal ages
#' @param batch is the batch size of random draws
#' @param it is the number of iterations of the bootstrap, an integer
#' @return the result with the lowest ks test value
#' @seealso accumulater
#' @seealso accumulate
#' @seealso recorder
#' @export
boot_accum <- function(n, ti_rng, to_rng, ia, obs, batch = 1000, it = 100) {
robs <- sample(obs, length(obs), replace = TRUE)
rec <- accumulater(n, ti_rng, to_rng, ia, robs, batch, it)
res <- rec[rec$ks == min(rec$ks), ]
res[1, ]
}
#' recorder
#'
#' models arrivals and removals from a reservoir
#'
#' @param n is the number of accumulated arrivals
#' @param ri is the rate in, a float
#' @param ro is the rate out, a float
#' @param ia is a vector of inherited ages, integers
#' @return deposit ages of accumulated arrivals, a numeric vector
#' @export
recorder <- function(n, ri, ro, ia) {
im <- 0
k <- 1
rec <- 0
om <- 0
too_old <- 50000
while (length(rec) < n &
im < too_old
) {
om <- om + rexp(1, ro)
while (
im < om &
length(rec) <= n
) {
im <- im + rexp(1, ri)
rec[k] <- im
k <- k + 1
}
m <- length(rec[rec < om])
if (m >= 1 ) {
s <- round(runif(1, 1, m))
rec <- rec[-s]
k <- k - 1
}
}
# rec <- rec[1:k]
rec <- sort(max(rec) - rec)
rc <- 0
for (i in seq_along(rec)) {
rc[i] <- rec[i] + sample(ia, 1)
}
return(sort(rc - min(rc)))
}
#' recorder1
#'
#' models arrivals and removals from a reservoir
#' adds a minimum age
#'
#' @param n is the number of accumulated arrivals
#' @param ri is the rate in, a float
#' @param ro is the rate out, a float
#' @param ia is a vector of inherited ages, integers
#' @param mn is a minimum age
#' @return deposit ages of accumulated arrivals, a numeric vector
#' @export
recorder1 <- function(n, ri, ro, ia, mn) {
im <- 0
k <- 1
rec <- 0
om <- 0
too_old <- 50000
while (length(rec) < n &
im < too_old
) {
om <- om + rexp(1, ro)
while (
im < om &
length(rec) <= n
) {
im <- im + rexp(1, ri)
rec[k] <- im
k <- k + 1
}
m <- length(rec[rec < om])
if (m >= 1 ) {
s <- round(runif(1, 1, m))
rec <- rec[-s]
k <- k - 1
}
}
# rec <- rec[1:k]
rec <- sort(max(rec) - rec)
rc <- 0
for (i in seq_along(rec)) {
rc[i] <- rec[i] + sample(ia, 1)
}
return(sort(rc - min(rc) + mn))
}
#' rec_vol
#'
#' function to record volumes
#'
#' @param node is a spatial object (creek node)
#' @param n is an integer representing the number of years to accumulate
#' @param vi is the volume input rate, a float
#' @param vo is the volume output rate, a float
#' @param return matrix of volumes and levels
#' @export
rec_vol <- function(node, n, vi, vo) {
t <- 0
age <- 0
vol <- 0
k <- 0
while (t < n) {
k <- k + 1
t <- t + rexp(1, node$ri)
age[k] <- t
vol[k] <- rexp(1, vi)
}
df <- data.frame(
t = age,
r = 'input',
vol = vol,
lvl = vol / as.numeric(node$valley_width)
)
t <- 0
age <- 0
vol <- 0
k <- 0
while (t < n) {
k <- k + 1
t <- t + rexp(1, node$ro)
age[k] <- t
vol[k] <- rexp(1, vo)
}
df <- rbind(df, data.frame(
t = age,
r = 'output',
vol = vol,
lvl = vol / as.numeric(node$valley_width)))
df <- df[order(df$t),]
vol <- 0
lvl <- 0
for (i in 1:nrow(df)) {
if (df$r[i] == 'input') {
vol <- vol + df$vol[i]
lvl <- lvl + df$lvl[i]
}
if (df$r[i] == 'output') {
if (df$vol[i] >= vol) {
vol <- 0
}
if (df$vol[i] < vol) {
vol <- vol - df$vol[i]
}
if (df$lvl[i] >= lvl) {
lvl <- 0
}
if (df$lvl[i] < lvl) {
lvl <- lvl - df$lvl[i]
}
}
}
data.frame(vol = vol, lvl = lvl)
}
#' fit_volumes
#'
#' function to fit volumes to creek nodes
#'
#' @param nodes is a spatial object (creek nodes)
#' @param n is an integer representing years to simulate accumulation
#' @param vi is the volumetric input rate, a float
#' @param vo is the volumetric output rate, a float
#' @param it is an integer representing the number of iterations per node
#' @return a data.frame with cols c('vol', 'lvl')
#' @export
fit_volumes <- function(nodes, n, vi, vo, it = 10) {
vol <- array(0, c(nrow(nodes), it))
lvl <- array(0, c(nrow(nodes), it))
for (i in 1:it) {
for (j in 1:nrow(nodes)) {
rec <- rec_vol(nodes[j, ], n, vi, vo)
vol[j, i] <- rec[1,1]
lvl[j, i] <- rec[1,2]
}
}
data.frame(
vol = apply(vol, 1, median),
lvl = apply(lvl, 1, median))
}
#' select_rate_pairs
#'
#' function to return rate pairs given a preference matrix
#'
#' @param ar is the preference matrix (a pdf superimposed on a grid)
#' @param xs is the vector of values along the x-axis of `ar`
#' @param ys is the vector of values along the y-axis of `ar`
#' @return a pair of floats c(x, y)
#' @export
select_rate_pairs <- function(ar, xs, ys) {
vec <- as.vector(ar)
ardim <- dim(ar)
rvec <- rep(1:ardim[1], ardim[2])
cvec <- 0
for (i in 1:ardim[2]) {
cvec <- c(cvec, rep(i, ardim[1]))
}
cvec <- cvec[-1]
cdf <- cumsum(vec)
xid <- rvec[cdf == min(cdf[cdf > runif(1)])]
yid <- cvec[cdf == min(cdf[cdf > runif(1)])]
x <- runif(1, xs[xid], xs[xid+1])
y <- runif(1, ys[yid], ys[yid+1])
c(x, y)
}
#' weight_array
#'
#' Returns preference matrix given topology matrix.
#' The scale parameter deflates preference weighting from zero to one,
#' inflates preference weighting above one,
#' and crashes below zero.
#'
#' @param ar is a topology matrix from weight_grid()
#' @param scl is the scale parameter
#' @return a preference matrix
#' @seealso weight_grid
#' @export
weight_array <- function(ar, scl = 1) {
vec <- as.vector(ar)
wmax <- max(vec)
p <- (1 - vec) * scl * (1 / length(vec))
p <- p / sum(p)
array(p, dim(ar))
}
#' weight_grid
#'
#' Returns topology of coordinate on grid given record of test values.
#' Types for `test`:
#' Kolmogorov-Smirnov statistic = 'ks'
#' Kuiper statistic = 'kp'
#' `xy_grid` determines the size of the outgoing topology array.
#' Types for `by`:
#' Unit cross-sectional area in m2 = 'vol'
#' Average cross-sectional depth per node = 'lvl'
#' `scl` ranges [0,1] and deflates the effects of neighbors on each other.
#'
#' @param pos is the coordinate pair of data in the record
#' @param xy_grid is the topology array to replace
#' @param rec is the record of test values
#' @param test is a string naming the test type c('ks', 'kp)
#' @param by is a string naming the test varaible c('vol', 'lvl')
#' @param scl is a scaling parameter, a float
#' @return the estimated test value at `pos`
#' @export
weight_grid <- function(pos, xy_grid, rec, test, by, scl = 1) {
val <- 0
if (test == 'ks') val <- 3
if (test == 'kp') val <- 5
if (by == 'lvl') val <- val + 1
x_dist <- xy_grid[, 1] - pos[1]
y_dist <- xy_grid[, 2] - pos[2]
dist <- sqrt(x_dist^2 + y_dist^2)
dmax <- max(dist) + 1
wts <- (dmax - dist) / dmax * scl
sum((wts * rec[, val])) / sum(wts)
}
#' update_search_array
#'
#' Returns preference topology given a record of test values.
#' Types for `test`:
#' Kolmogorov-Smirnov statistic = 'ks'
#' Kuiper statistic = 'kp'
#' Types for `by`:
#' Unit cross-sectional area in m2 = 'vol'
#' Average cross-sectional depth per node = 'lvl'
#'
#' @param rec is a data.frame of test values with cols c(x, y, ks, kp)
#' @param xs is a vector of values for the x-axis of topology array
#' @param ys is a vector of values for the y-axis of topology array
#' @param ar is a preference topology to update
#' @param test is a string naming the test type c('ks', 'kp)
#' @param by is a string naming the test varaible c('vol', 'lvl')
#' @param d_scl is the scale value for weight_grid()
#' @param w_scl is the scale value for weight_array()
#' @return a preference topology array
update_search_array <- function(rec, xs, ys, ar, test, by,
d_scl = 1, w_scl = 0.9) {
xy_grid <- matrix(0, nrow = nrow(rec), ncol = 2)
for (i in 1:nrow(rec)) {
x_bin <- 0
k <- 1
while (!x_bin) {
k <- k + 1
if (rec[i, 1] <= xs[k] & rec[i, 1] > xs[k-1]) x_bin <- k
}
y_bin <- 0
k <- 1
while (!y_bin) {
k <- k + 1
if (rec[i, 2] <= ys[k] & rec[i, 2] > ys[k-1]) y_bin <- k
}
xy_grid[i, ] <- c(x_bin, y_bin)
}
for (i in 1:nrow(ar)) {
for (j in 1:ncol(ar)) {
ar[i, j] <- weight_grid(c(i,j), xy_grid, rec, test, by, d_scl)
}
}
weight_array(ar, w_scl)
}
#' search_topology
#'
#'
#'
#' @param so is a spatial object (creek nodes)
#' @param x_range is a vector of values for the x-axis of the preference array
#' @param y_range is a vector of values for the y-axis of the preference array
#' @param rec is a data.frame of test values with cols c(x, y, ks, kp)
#' @param bt is the number of times to simulate accumulation over `so`
#' @param it is the number of times to simulate accumulation per node of `so`
#' @param n is the number of years to accumulate deposits
#' @param test is a string naming the test type c('ks', 'kp)
#' @param by is a string naming the test varaible c('vol', 'lvl')
#' @param grid is the number of intervals to divide x and y ranges into
#' @param d_scl is the scale value for weight_grid()
#' @param w_scl is the scale value for weight_array()
#' @return a record of test value results at randomly generated coords in `x_range` and `y_range`
#' @export
search_topology <- function(so,
x_range,
y_range,
rec = 0,
bt = 10,
it = 20,
n = 10000,
test = 'ks',
by = 'vol',
grid = 100,
d_scl = 1,
w_scl = 3,
linked = FALSE,
backfill = TRUE
) {
x_dist <- max(x_range) - min(x_range)
y_dist <- max(y_range) - min(y_range)
xs <- seq(
min(x_range),
max(x_range),
x_dist / grid
)
ys <- seq(
min(y_range),
max(y_range),
y_dist / grid
)
ar <- array(1/(grid^2), c(grid, grid))
if (length(rec) == 1) vec <- 0
if (length(rec) > 1) {
vec <- rec[1, ]
if (nrow(rec) > 1) {
for (i in 2:nrow(rec)) {
vec <- c(vec, rec[i, ])
}
}
}
begin <- Sys.time()
for (i in 1:bt) {
out <- 0
if (length(vec) > 1) {
ar <- update_search_array(
matrix(unlist(vec), ncol = 7, byrow = TRUE),
xs, ys, ar, test, by, d_scl, w_scl)
present_topology(ar, xs, ys)
}
rates <- select_rate_pairs(ar, xs, ys)
if (linked) {
res <- fit_volumes1(so, n, rates[1], rates[2], it)
boot <- res[[1]]
out <- res[[2]]
}
if (backfill) {
res <- backfill(so, n, rates[1], rates[2], it)
boot <- res[[1]]
out <- res[[2]]
}
if (!linked & !backfill) boot <- fit_volumes(so, n, rates[1], rates[2], it)
vol_gof <- cdf_gof(boot[, 1], so$xsec_area)
lvl_gof <- cdf_gof(boot[, 2], so$lvl)
if (length(vec) > 1) vec <- c(vec, rates, vol_gof, lvl_gof, out)
if (length(vec) == 1) vec <- c(rates, vol_gof, lvl_gof, out)
now <- Sys.time()
dif <- now - begin
pace <- dif / i
rem <- (pace * bt) - dif
print(paste0('Percent Done: ', round(i/bt*100, 2), '%'))
print(paste0('Time Elapsed: ', round(as.numeric(dif, units = 'hours'), 2), ' hours'))
print(paste0('Time Remaining: ', round(as.numeric(rem, units = 'hours'), 2), ' hours'))
}
end <- Sys.time()
print(end - begin)
matrix(unlist(vec), ncol = 7, byrow = TRUE)
}
#' present_topology
#'
#' prints a visualization of the topology array
#'
#' @param ar is a topology array
#' @param xs is a vector of values for the x-axis of the array
#' @param ys is a vector of values for the y-axis of the array
#' @return prints a visualization of the topology array
#' @export
present_topology <- function(ar, xs, ys) {
xvec <- rep(1:nrow(ar), ncol(ar))
yvec <- c(
matrix(rep(1:ncol(ar), nrow(ar)), ncol = ncol(ar), byrow = TRUE)
)
vec <- 1 - (c(ar) / max(ar))
print(scatter3D(xs[xvec], ys[yvec], vec, pch = 20, ticktype = 'detailed',
phi = 40, theta = 320))
return()
}
#' weight_by_delprob
#'
#' converts MB delivery probabilities to Rose-Lancaster optimized delivery probabilites
#'
#' @param delprob is a vector of delivery probabilities, f64
#' @param weight is a matrix with cols c(cdf, delprob)
#' @return vector of optimized delivery probabilities
#' @export
weight_by_delprob <- function(delprob, weight = df_wt) {
flag <- 0
k <- 1
ln <- nrow(weight)
while (k <= ln & flag == 0) {
if (weight[k, 2] >= delprob) {
flag <- k
}
if (weight[k, 2] < delprob & k < ln) {
k <- k + 1
}
if (weight[k, 2] < delprob & k == ln) {
flag <- k
}
}
res <- 0
if (flag > 0) {
res <- weight[flag, 1]
}
if (flag == 0) {
res <- NA
}
res
}
#' record_volume
#'
#' Record acculumation at a node. Takes arrivals and emits departures.
#' Nodes are linked, with deposits departing downstream.
#'
#' @param node is a row from a spatial dataframe (creek nodes)
#' @param n is the number of years to simulate accumulation
#' @param vi is the volumetric input rate
#' @param vo is the volumetric output rate
#' @param arrivals is a dataframe describing upstream arrivals
#' @return a list with elements c(accumulation, departures)
#' @seealso fit_volumes1
#' @export
record_volume <- function(node, n, vi, vo,
arrivals = 0) {
departures <- 0
if (length(arrivals) != 1) {
rolls <- runif(nrow(arrivals))
if (sum(rolls > node$rdp) >= 1) {
departures <- arrivals[rolls > node$rdp, ]
arrivals <- arrivals[rolls <= node$rdp, ]
}
arrivals$lvl <- arrivals$vol / as.numeric(node$valley_width)
}
t <- 0
age <- 0
vol <- 0
k <- 0
while (t < n) {
k <- k + 1
t <- t + rexp(1, node$ri)
age[k] <- t
vol[k] <- rexp(1, vi)
}
df <- data.frame(
t = age,
r = 'input',
vol = vol,
lvl = vol / as.numeric(node$valley_width)
)
if (length(arrivals) != 1) df <- rbind(df, arrivals)
t <- 0
age <- 0
vol <- 0
k <- 0
while (t < n) {
k <- k + 1
t <- t + rexp(1, node$ro)
age[k] <- t
vol[k] <- rexp(1, vo)
}
dep <- data.frame(
t = age,
r = 'output',
vol = vol,
lvl = vol / as.numeric(node$valley_width))
if (length(departures) == 1) departures <- dep
if (length(departures) != 1) departures <- rbind(departures, dep)
departures$r <- 'input'
df <- rbind(df, dep)
df <- df[order(df$t),]
vol <- 0
lvl <- 0
for (i in 1:nrow(df)) {
if (df$r[i] == 'input') {
vol <- vol + df$vol[i]
lvl <- lvl + df$lvl[i]
}
if (df$r[i] == 'output') {
if (df$vol[i] >= vol) {
vol <- 0
}
if (df$vol[i] < vol) {
vol <- vol - df$vol[i]
}
if (df$lvl[i] >= lvl) {
lvl <- 0
}
if (df$lvl[i] < lvl) {
lvl <- lvl - df$lvl[i]
}
}
}
res <- data.frame(vol = vol, lvl = lvl)
list(res, departures)
}
#' fit_volumes1
#'
#' Function to fit volume rate to observed volumes.
#' Creek nodes are linked, passing deposits downstream.
#'
#' @param nodes is a spatial object (creek nodes)
#' @param n is the number of years to simulate accumulation, an integer
#' @param vi is the volumetric input rate
#' @param vo is the volumetric output rate
#' @param it is the number of iterations to simulate accumulation, an integer
#' @return the median accumulation record after `it` iterations in volume and level
#' @seealso record_volume
#' @export
fit_volumes1 <- function(nodes, n, vi, vo, it = 10) {
vol <- array(0, c(nrow(nodes), it))
lvl <- array(0, c(nrow(nodes), it))
out <- 0
nodes <- nodes[order(nodes$ToMouth_km, decreasing = TRUE), ]
arrivals <- data.frame(t = 0, r = 'input', vol = 0, lvl = 0)
res <- mclapply(1:it, function(x) voluminous(nodes, n, vi, vo, arrivals))
for (i in 1:it) {
rec <- res[[i]]
vol[ , i] <- rec[[1]]
lvl[ , i] <- rec[[2]]
dep <- rec[[3]]
dep <- dep[dep$t < n, ]
out[i] <- sum(dep$vol) / n
}
df <- data.frame(
vol = apply(vol, 1, median),
lvl = apply(lvl, 1, median)
)
return(list(df, mean(out)))
}
#' voluminous
#'
#' Wrapper for record_volume.
#'
#' @param nodes is a spatial object (creek nodes)
#' @param n is the number of years to simulate accumulation, an integer
#' @param vi is the volumetric input rate
#' @param vo is the volumetric output rate
#' @param arrivals is the record of arrivals, a data.frame with cols (t, r, vol, lvl)
#' @return a list with elements c(volumes, levels, arrivals)
#' @export
voluminous <- function(nodes, n, vi, vo, arrivals) {
vol <- 0
lvl <- 0
for (i in 1:nrow(nodes)) {
res <- record_volume(nodes[i, ], n, vi, vo, arrivals)
rec <- res[[1]]
arrivals <- res[[2]]
vol[i] <- rec[1,1]
lvl[i] <- rec[1,2]
}
return(list(vol = vol, lvl = lvl, arrivals = arrivals))
}
#' voluminous1
#'
#' Simulates accumulation record, using a linked-bucket model
#' with backfilling based on average elevation of bank deposits.
#'
#' @param nodes is a spatial object (creek nodes)
#' @param n is the number of years to simulate accumulation, an integer
#' @param vi is the volumetric input rate
#' @param vo is the volumetric output rate
#' @param b_scl is the proportion of volume that settles out of deposits at obstructions
#' @return a list with elements c(volumes, levels, arrivals)
#' @seealso eventer
#' @export
voluminous1 <- function(nodes, n, vi, vo, b_scl) {
nodes <- nodes[order(nodes$ToMouth_km), ]
vol <- array(0, nrow(nodes))
lvl <- array(0, nrow(nodes))
elev <- nodes$elev
arr <- mapply(function(x, y, z, n, m) as.matrix(eventer(x, y, z, n, m)),
x = nodes$ri, y = nodes$NODE_ID,
MoreArgs = list(z = vi, n = n, m = 0))
rec <- arr[[1]]
for (i in 2:length(arr)) {
rec <- rbind(rec, arr[[i]])
}
dep <- mapply(function(x, y, z, n, m) as.matrix(eventer(x, y, z, n, m)),
x = nodes$ro, y = nodes$NODE_ID,
MoreArgs = list(z = vo, n = n, m = 1))
for (i in seq_along(dep)) {
rec <- rbind(rec, dep[[i]])
}
rec <- rec[order(rec[, 1]), ]
# print('record:')
# print(head(rec[rec[,1] > 0, ], 50))
ids <- nodes$NODE_ID
departures <- data.frame(age = 0, vol = 0, type = 0, id = 0)
for (i in 1:nrow(rec)) {
# print(paste0('on record ', i))
k <- sum((ids == rec[i, 4]) * 1:length(ids))
# print(paste0('on index ', k))
# print(paste0('is input: ', rec[i, 3] == 0))
if (rec[i, 3] == 0) {
vol[k] <- vol[k] + rec[i, 2]
lev <- lvl[k] + rec[i, 2] /
as.numeric(nodes$valley_width[k])
lvl[k] <- lev
elev[k] <- elev[k] + lev
}
# print(paste0('is output: ', rec[i, 3] == 1))
if (rec[i, 3] == 1) {
scoop <- 0
settle <- 0
if (rec[i, 2] > vol[k]) {
scoop <- vol[k]
}
if (rec[i, 2] <= vol[k]) {
scoop <- rec[i, 2]
}
# print(paste0('scoop is ', scoop))
og_k <- k
vol[k] <- vol[k] - scoop
lvl[k] <- lvl[k] - scoop / as.numeric(nodes$valley_width[k])
elev[k] <- elev[k] - scoop / as.numeric(nodes$valley_width[k])
going <- TRUE
idl <- length(ids)
while (k <= idl & going == TRUE) {
if (k == idl) {
rec[i, 2] <- scoop
departures <- rbind(departures, rec[i, ])
going <- FALSE
# print('off and away')
}
if (k < idl) {
k <- k + 1
# print(paste0('going to ', k))
if (lvl[k-1] < lvl[k]) {
settle <- scoop * b_scl
scoop <- scoop - settle
vol[k] <- vol[k] + settle
lvl[k] <- lvl[k] + settle / as.numeric(nodes$valley_width[k])
elev[k] <- elev[k] + settle / as.numeric(nodes$valley_width[k])
# print('never got started')
}
if (lvl[k-1] >= lvl[k]) {
r <- runif(1)
if (r <= nodes$rdp[k]) {
settle <- scoop * b_scl
scoop <- scoop - settle
vol[k] <- vol[k] + settle
lvl[k] <- lvl[k] + settle / as.numeric(nodes$valley_width[k])
elev[k] <- elev[k] + settle / as.numeric(nodes$valley_width[k])
# print(paste0('new vol is: ', vol[k]))
}
}
}
}
}
}
return(list(vol = vol, lvl = lvl, elev = elev, departures = departures)
)
}
#' eventer
#'
#' Utility for voluminous().
#' Generates an event data.frame with cols (age, vol, type, id),
#' given a temporal and volumetric event rate, over a given period of time.
#'
#' @param rt is the temporal event rate
#' @param id is the creek node number, an integer.
#' @param vl is the volumetric event rate
#' @param n is the period of time to simulate events
#' @param lab is the type label, which sounds like a string but is numeric!
#' @return a data.frame of events with cols (age, vol, type, id)
#' @export
eventer <- function(rt, id, vl, n, lab) {
t <- 0
age <- 0
vol <- 0
k <- 1
while(t < n) {
t <- t + rexp(1, rt)
if (t <= n) {
age[k] <- t
vol[k] <- rexp(1, vl)
k <- k + 1
}
}
return(
data.frame(
age = age,
vol = vol,
type = rep(lab, length(vol)),
id = rep(id, length(vol))
)
)
}
#' backfill
#'
#' Function to fit volume rate to observed volumes.
#' Creek nodes are linked, passing deposits downstream.
#' Deposits backfill based upon average elevation
#'
#' @param nodes is a spatial object (creek nodes)
#' @param n is the number of years to simulate accumulation, an integer
#' @param vi is the volumetric input rate
#' @param vo is the volumetric output rate
#' @param it is the number of iterations to simulate accumulation, an integer
#' @param b_scl is the proportion of volume that settles out of deposits at obstructions
#' @return the median accumulation record after `it` iterations in volume and level
#' @seealso record_volume
#' @export
backfill <- function(nodes, n, vi, vo, b_scl, it = 10) {
vol <- array(0, c(nrow(nodes), it))
lvl <- array(0, c(nrow(nodes), it))
elev <- array(0, c(nrow(nodes), it))
out <- 0
nodes <- nodes[order(nodes$ToMouth_km, decreasing = TRUE), ]
arrivals <- data.frame(t = 0, r = 'input', vol = 0, lvl = 0)
res <- lapply(1:it, function(x) voluminous1(nodes, n, vi, vo, b_scl))
for (i in 1:it) {
rec <- res[[i]]
vol[ , i] <- rec[[1]]
lvl[ , i] <- rec[[2]]
elev[ , i] <- rec[[3]]
dep <- rec[[4]]
dep <- dep[dep[,1] < n, ]
out[i] <- sum(dep$vol) / n
}
df <- data.frame(
vol = apply(vol, 1, median),
lvl = apply(lvl, 1, median),
elev = apply(elev, 1, median)
)
return(list(record = df, output = mean(out)))
}
#' slope_to_elev
#'
#' Converts slopes to elevations based on distance from a base.
#'
#' @param slope is a numeric vector of slopes at creek nodes
#' @param dist is a numeric vector of distances to mouth at creek nodes
#' @param base is the starting elevation
#' @return a vector of elevations
#' @export
slope_to_elev <- function(slope, dist, base = 0) {
elev <- base
for (i in seq_along(slope)) {
if (i == 1) {
elev[i] <- base
}
if (i > 1) {
elev[i] <- elev[i-1] + ((dist[i] - dist[i-1]) * slope[i-1] * 1000)
}
}
elev
}
#' n_dimensional_search
#'
#' For finding a needle in a haystack.
#' Uses backfill() to simulate accumulation.
#'
#' @param so is a spatial object (creek nodes)
#' @param ti_rng is the range of temporal input rates to search
#' @param to_rng is the range of temporal output rates to search
#' @param vi_rng is the range of volumetric input rates to search
#' @param vo_rng is the range of volumetric input rates to search
#' @param ti_scl is the range of temporal input scales to search
#' @param to_scl is the range of temporal output scales to search
#' @param ks_scl is the range of streampower scales to search
#' @param dp_scl is the range of deposit probability scales to search
#' @param b_scl is the proportion of volume that settles out of deposits at obstructions
#' @param it is the number of iterations to perform at each point
#' @param n is the number of years to simulate accumulation, an integer
#' @seealso backfill
#' @return a list with elements (record, mean output rate)
#' @export
n_dimensional_search <- function(so,
ti_rng,
to_rng,
vi_rng,
vo_rng,
ti_scl = 1,
to_scl = 1,
ks_scl = 1,
dp_scl = 1,
m_scl = 1,
b_scl = 1,
type = 0,
it = 100,
n = 10000,
batch = 10) {
go <- runif(batch, min(to_rng), max(to_rng))
gi <- runif(batch, go, max(ti_rng))
vo <- runif(batch, min(vo_rng), max(vo_rng))
vi <- runif(batch, min(vi_rng), max(vi_rng))
to_scale <- runif(batch, min(to_scl), max(to_scl))
ks_scale <- runif(batch, min(ks_scl), max(ks_scl))
ti_scale <- runif(batch, min(ti_scl), max(ti_scl))
dp_scale <- runif(batch, min(dp_scl), max(dp_scl))
m_scale <- runif(batch, min(m_scl), max(m_scl))
b_scale <- runif(batch, min(b_scl), max(b_scl))
if (type == 0) type <- sample(c('ks', 'dp', 'mixed'), 1)
res <- mcmapply(
function (x, a, b, c, d, e, f, g, h, j, k, l, m, p) backfill(
rater(x, a, b, c, d, e, f, g, h), l, j, k, m, p),
a = gi, b = go, c = ti_scale, d = to_scale, e = ks_scale, f = dp_scale, g = m_scale, j = vi, k = vo, m = b_scale,
MoreArgs = list(x = so, h = type, l = n, p = it)
)
vol_gof <- array(0, c(batch, 2))
lvl_gof <- array(0, c(batch, 2))
elev_gof <- array(0, c(batch, 2))
out <- array(0, batch)
vor <- array(0, batch)
for (i in 1:(length(res)/2)) {
bt <- res[[i + (i-1)]]
vol_gof[i, ] <- as.numeric(cdf_gof(as.numeric(unlist(bt[1])), so$xsec_area))
lvl_gof[i, ] <- as.numeric(cdf_gof(as.numeric(unlist(bt[2])), so$xsec_area / as.numeric(so$valley_width)))
elev_gof[i, ] <- as.numeric(cdf_gof(as.numeric(unlist(bt[3])), slope_to_elev(so$slope, so$ToMouth_km)))
out[i] <- res[[i + (i-1) + 1]]
}
data.frame(
ti = gi,
to = go,
vi = vi,
vo = vo,
ti_scl = ti_scale,
to_scl = to_scale,
ks_scl = ks_scale,
dp_scl = dp_scale,
m_scl = m_scale,
b_scl = b_scale,
ks_vol = vol_gof[ ,1],
kp_vol = vol_gof[ ,2],
ks_lvl = lvl_gof[ ,1],
kp_lvl = lvl_gof[ ,2],
ks_elev = elev_gof[, 1],
kp_elev = elev_gof[ ,2],
out = out,
type = type,
stringsAsFactors = FALSE
)
}
#' rater
#'
#' Derives variables related to temporal input and output rates.
#' The average rate per node is the global rate divided by the number of sampled nodes in sa_nodes.
#' Types of `type` include c('ks', 'dp', 'mixed'):
#' Streampower coefficient = 'ks'.
#' Inverse delivery probability = 'dp'.
#' Both = 'mixed'.
#'
#' @param so is a spatial object (creek nodes)
#' @param gi is a global temporal input rate
#' @param go is a global temporal output rate
#' @param ti_scl is the range of temporal input scales to search
#' @param to_scl is the range of temporal output delivery scales to search
#' @param ks_scl is the range of temporal output streampower scales to search
#' @param dp_scl is the range of deposit probability scales to search
#' @param m_scl is the range of distance probability scales to search
#' @param type is a string indicated the type of output weight to apply
#' @param sa_nodes is the spatial object of sampled nodes
#' @return spatial object with rates calculated
#' @seealso n_dimensional_search
#' @export
rater <- function(so, gi, go, ti_scl, to_scl, ks_scl, dp_scl, m_scl, type, sa = sa_nodes) {
so <- so[order(so$ToMouth_km, decreasing = FALSE), ]
wt <- unlist(lapply(so$DebrisFlow, function(x) weight_by_delprob(x)))
dp <- so$DebrisFlow * wt
rdp <- dp * dp_scl
ldp_k <- exp(log(dp) - mean(log(dp)))
# wdp <- log(dp) - min(log(dp))
lvl <- so$xsec_area / as.numeric(so$valley_width)
elev <- slope_to_elev(so$slope, so$ToMouth_km)
pwr_df <- data.frame(CA = log(so$contr_area),
slope = log(so$slope))
pwr_mod <- lm(CA ~ slope, data = pwr_df)
pwr_prd <- predict(pwr_mod, newdata = data.frame(slope = log(so$slope)))
pwr_k <- log(so$contr_area) - pwr_prd
wdp <- -pwr_k
wdp <- wdp - min(wdp)
mth <- so$ToMouth_km - min(so$ToMouth_km)
mth <- mth / max(mth)
mth <- 1 - (mth - m_scl * mth^2)
to_ave <- go / nrow(sa)
ti_ave <- gi / nrow(sa)
to <- 0
if (type == 'ks') to <- to_ave * exp(pwr_k) * ks_scl
if (type == 'dp') to <- to_ave * (1 / ldp_k) * to_scl
if (type == 'mixed') to <- to_ave * exp(pwr_k) * (1 / ldp_k) * ks_scl * to_scl
to[to <= 0] <- min(to[to > 0])
to <- to * sum(to_ave * nrow(so)) / sum(to)
ti <- ti_ave * wdp * mth * ti_scl
ti[ti <= 0] <- min(ti[ti > 0])
ti <- ti * sum(ti_ave * nrow(so)) / sum(ti)
so$dp <- dp
so$rdp <- rdp
so$ri <- ti
so$ro <- to
so$lvl <- lvl
so$elev <- elev
so
}
#' voluminous2
#'
#' Simulates accumulation record, using a linked-bucket model
#' with backfilling based on average elevation of bank deposits.
#'
#' @param nodes is a spatial object (creek nodes)
#' @param n is the number of years to simulate accumulation, an integer
#' @param vi is the volumetric input rate
#' @param vo is the volumetric output rate
#' @param tp is the turbulent deposition probability
#' @param lp is the laminar deposition probability
#' @return a list with elements c(volumes, levels, arrivals)
#' @seealso eventer
#' @export
voluminous2 <- function(nodes, n, vi, vo, tp = 0.1, lp = 0.5) {
nodes <- nodes[order(nodes$ToMouth_km, decreasing = TRUE), ]
vol <- array(0, nrow(nodes)) # unit cross-sectional area m2
lvl <- array(0, nrow(nodes)) # average valley depth m
lam <- array(0, nrow(nodes)) # laminar or turbulent boolean
elev <- nodes$elev
recs <- 1:nrow(nodes)
recs <- lapply(recs, function(x) data.frame(id = x, vol = 0, age = 0, type = 0))
# assemble record of arrivals and departures
arr <- mapply(function(x, y, z, n, m) as.matrix(eventer(x, y, z, n, m)),
x = nodes$ri, y = nodes$NODE_ID,
MoreArgs = list(z = vi, n = n, m = 0))
rec <- arr[[1]]
for (i in 2:length(arr)) {
rec <- rbind(rec, arr[[i]])
}
dep <- mapply(function(x, y, z, n, m) as.matrix(eventer(x, y, z, n, m)),
x = nodes$ro, y = nodes$NODE_ID,
MoreArgs = list(z = vo, n = n, m = 1))
for (i in seq_along(dep)) {
rec <- rbind(rec, dep[[i]])
}
rec <- rec[order(rec[, 1]), ]
ids <- nodes$NODE_ID
departures <- data.frame(id =0, vol = 0, age = 0, type = 0)
# process arrivals and departures
for (i in 1:nrow(rec)) {
# k is id in ids of node on which rec occurs
k <- sum((ids == rec[i, 4]) * 1:length(ids))
# deposit input
if (rec[i, 3] == 0) {
vol[k] <- vol[k] + rec[i, 2]
lev <- lvl[k] + rec[i, 2] /
as.numeric(nodes$valley_width[k])
lvl[k] <- lev
recs[[k]] <- rbind(recs[[k]], c(k, rec[i, 2], rec[i, 1], 0))
# mark upstream portions laminar (inundated)
# if their elevation is below the height of the deposit
thresh <- elev[k] + lvl[k]
l <- k
if (l > 1) l <- l - 1
while ((elev[l] + lvl[l]) <= thresh & l > 1) {
lam[l] <- 1
l <- l - 1
}
}
# withdraw output
if (rec[i, 3] == 1 & vol[k] > 0) {
scoop <- 0
if (rec[i, 2] > vol[k]) {
scoop <- vol[k]
}
if (rec[i, 2] <= vol[k]) {
scoop <- rec[i, 2]
}
vol[k] <- vol[k] - scoop
lvl[k] <- lvl[k] - scoop / as.numeric(nodes$valley_width[k])
going <- TRUE
idl <- length(ids)
# update node volume record
slug <- data.frame(id = k, vol = 0, age = 0, type = 1)
cup <- 0
nrec <- recs[[k]]
while (sum(slug[, 2]) < scoop) {
strat <- sample(1:nrow(nrec), 1)
cup <- scoop - sum(slug[, 2])
if (cup <= nrec[strat, 2]) { # scoop takes portion of strat
slug <- rbind(slug, c(k, cup, nrec[strat, 3], 1))
nrec[strat, 2] <- nrec[strat, 2] - cup
}
if (cup > nrec[strat, 2]) { # scoop consumes strat
slug <- rbind(slug, c(k, nrec[strat, 2], nrec[strat, 3], 1))
}
}
recs[[k]] <- nrec
slug <- slug[slug[, 2] > 0, ] # drop zero volumes from slug
# send output downstream
while (k <= idl & going == TRUE) {
# case: end of study area
if (k == idl) {
departures <- rbind(departures, slug)
going <- FALSE
}
# case: going downstream
if (k < idl) {
k <- k + 1
r <- runif(1)
settle <- r <= tp # turbulent probability
if (lam[k]) settle <- r <= lp # laminar probability
if (settle) {
lump <- sum(slug[, 2])
vol[k] <- vol[k] + lump
lvl[k] <- lvl[k] + lump / as.numeric(nodes$valley_width[k])
recs[[k]] <- rbind(recs[[k]], slug)
going <- FALSE
}
}
}
}
}
rex <- 0
for (i in seq_along(recs)) {
if (i == 1) rex <- recs[[i]]
if (i > 1 ) rex <- rbind(rex, recs[[i]])
}
return(list(vol = vol, lvl = lvl, recs = rex, dep = departures)
)
}
#' backfill1
#'
#' Function to fit volume rate to observed volumes.
#' Creek nodes are linked, passing deposits downstream.
#' Deposits backfill based upon average elevation
#'
#' @param nodes is a spatial object (creek nodes)
#' @param n is the number of years to simulate accumulation, an integer
#' @param vi is the volumetric input rate
#' @param vo is the volumetric output rate
#' @param tp is the turbulent deposition probability
#' @param lp is the laminar deposition probability
#' @param it is the number of iterations to simulate accumulation, an integer
#' @return the median accumulation record after `it` iterations in volume and level
#' @seealso record_volume
#' @export
backfill1 <- function(nodes, n, vi, vo, tp, lp, it = 10) {
vol <- array(0, c(nrow(nodes), it))
lvl <- array(0, c(nrow(nodes), it))
out <- 0
nodes <- nodes[order(nodes$ToMouth_km, decreasing = TRUE), ]
arrivals <- data.frame(t = 0, r = 'input', vol = 0, lvl = 0)
res <- lapply(1:it, function(x) voluminous2(nodes, n, vi, vo, tp, lp))
for (i in 1:it) {
rec <- res[[i]]
vol[ , i] <- rec[[1]]
lvl[ , i] <- rec[[2]]
recs <- rec[[3]]
df_vol <- sum(recs$vol[recs$type == 0])
df_pct <- df_vol / sum(recs$vol)
dep <- rec[[4]]
dep <- dep[dep[,3] < n, ]
out[i] <- sum(dep$vol) / n
}
df <- data.frame(
vol = apply(vol, 1, median),
lvl = apply(lvl, 1, median)
)
return(list(record = df, output = mean(out), df_pct = df_pct))
}
#' ndim_search
#'
#' For finding a needle in a haystack.
#' Uses backfill() to simulate accumulation.
#'
#' @param so is a spatial object (creek nodes)
#' @param ti_rng is the range of temporal input rates to search
#' @param to_rng is the range of temporal output rates to search
#' @param vi_rng is the range of volumetric input rates to search
#' @param vo_rng is the range of volumetric input rates to search
#' @param tp is the turbulent deposition probability
#' @param lp is the laminar deposition probability
#' @param it is the number of iterations to perform at each point
#' @param n is the number of years to simulate accumulation, an integer
#' @seealso backfill
#' @return a list with elements (record, mean output rate)
#' @export
ndim_search <- function(so,
ti_rng,
to_rng,
vi_rng,
vo_rng,
tp = 0.1,
lp = 0.5,
it = 100,
n = 10000,
batch = 10) {
go <- runif(batch, min(to_rng), max(to_rng))
gi <- runif(batch, go, max(ti_rng))
vo <- runif(batch, min(vo_rng), max(vo_rng))
vi <- runif(batch, min(vi_rng), max(vi_rng))
tps <- runif(batch, min(tp), max(tp))
lps <- runif(batch, min(lp), max(lp))
res <- mcmapply(
function (x, a, b, c, d, e, f, g, h) backfill1(
rater1(x, a, b), c, d, e, f, g, h),
a = gi, b = go, d = vi, e = vo, f = tps, g = lps,
MoreArgs = list(x = so, c = n, h = it)
)
vol_gof <- array(0, c(batch, 2))
lvl_gof <- array(0, c(batch, 2))
out <- array(0, batch)
df_pct <- array(0, batch)
for (i in 1:(length(res)/3)) {
bt <- res[[3 * (i-1) + 1]]
vol_gof[i, ] <- as.numeric(cdf_gof(as.numeric(unlist(bt[1])), so$xsec_area))
lvl_gof[i, ] <- as.numeric(cdf_gof(as.numeric(unlist(bt[2])), so$xsec_area / as.numeric(so$valley_width)))
out[i] <- res[[3 * (i-1) + 2]]
df_pct[i] <- res[[3 * (i-1) + 3]]
}
data.frame(
ti = gi,
to = go,
vi = vi,
vo = vo,
tp = tps,
lp = lps,
ks_vol = vol_gof[ ,1],
kp_vol = vol_gof[ ,2],
ks_lvl = lvl_gof[ ,1],
kp_lvl = lvl_gof[ ,2],
out = out,
df_pct = df_pct,
stringsAsFactors = FALSE
)
}
#' rater1
#'
#' Derives variables related to temporal input and output rates.
#' The average rate per node is the global rate divided by the number of sampled nodes in sa_nodes.
#'
#' @param so is a spatial object (creek nodes)
#' @param gi is a global temporal input rate
#' @param go is a global temporal output rate
#' @param sa_nodes is the spatial object of sampled nodes
#' @return spatial object with rates calculated
#' @seealso n_dimensional_search
#' @export
rater1 <- function(so, gi, go, sa = sa_nodes) {
so <- so[order(so$ToMouth_km, decreasing = FALSE), ]
wt <- unlist(lapply(so$DebrisFlow, function(x) weight_by_delprob(x)))
dp <- so$DebrisFlow * wt
dp_wt <- log(dp) - min(log(dp))
lvl <- so$xsec_area / as.numeric(so$valley_width)
elev <- slope_to_elev(so$slope, so$ToMouth_km)
pwr_df <- data.frame(CA = log(so$contr_area),
slope = log(so$slope))
pwr_mod <- lm(CA ~ slope, data = pwr_df)
pwr_prd <- predict(pwr_mod, newdata = data.frame(slope = log(so$slope)))
pwr_k <- log(so$contr_area) - pwr_prd
ca_wt <- so$contr_area / max(so$contr_area)
to_ave <- go / nrow(sa)
ti_ave <- gi / nrow(sa)
to <- to_ave * exp(pwr_k)
to[to <= 0] <- min(to[to > 0])
to <- to * sum(to_ave * nrow(so)) / sum(to)
ti <- ti_ave * dp_wt * ca_wt / exp(pwr_k)
ti[ti <= 0] <- min(ti[ti > 0])
ti <- ti * sum(ti_ave * nrow(so)) / sum(ti)
so$k <- pwr_k
so$dp <- dp
so$ri <- ti
so$ro <- to
so$lvl <- lvl
so$elev <- elev
so
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.