sit coverage - 84.15%

Files
Source

#' Mean Distance Travelled and Mean Dispersal Distance
#' @family estimates
#'
#' Compute the Mean Distance Travelled (MDT) realised in a SIT experiment.
#' The Mean Dispersal Distance (MDD) is the MDT without accounting for
#' inhomogeneous arrangements of traps. See Details.
#'
#' The mean distance travelled (MDT) is the average distance travelled by all
#' __captured__ individuals, weighted by the relative density of traps.
#'
#' In practice, it is computed as the average distance between the traps
#' and the release point, weighted by the number of individuals from the
#' released population and by the relative density of traps.
#'
#' MDT estimates result in `NaN`, when there are no observations of captures in
#' the corresponding group. This is more likely when estimates are requested by
#' population and age.
#'
#' \deqn{
#'   \text{MDT}_{jk} = \sum_{i = 1}^{n_t} w_i \, n_{ijk} \, d_i \bigg/
#'     \sum_{i = 1}^{n_t} w_i \, n_{ijk}
#' }
#'
#' @inheritParams dispersal
#' @param by Character vector. Either 'population', or 'age', both (default), or
#'   NULL.
#'
#' @return By default (i.e. `by = c('population', 'age')`), a table with MDT
#'   computed for each released population at each _age_ (i.e., number of days
#'   since release). If `by` is either one of the grouping variables, the
#'   results will be presented for each value of the corresponding variable.
#'   Finally, if `by = NULL`, all populations and ages are pooled together and
#'   the function returns a single number.
#' @export
#'
#' @examples
#' sit_mdt(sit_prototype)
#' sit_mdt(sit_prototype, by = 'population')
#' sit_mdt(sit_prototype, by = 'age')
#' sit_mdt(sit_prototype, by = NULL)
sit_mdt <- function(
  x,
  by = c('population', 'age'),
  spatial_adjustment = TRUE,
  # weights = inverse_density(sit_traps(x, stage = 'adult', area = 'sit')),
  following_releases,
  following_days,
  species = NULL
) {

  dispersal_dat <- dispersal(
    x,
    spatial_adjustment = spatial_adjustment,
    following_releases = following_releases,
    following_days = following_days,
    species = species
  )

  ## Weighted mean of distances by counts and trap weights
  mdt <- function(x) {
    stats::weighted.mean(x$dist_m, x$n * x$w)
  }

  ans <- if (is.null(by)) {
    mdt(dispersal_dat)
  } else {
    mdt_values <- by(dispersal_dat, rev(dispersal_dat[by]), mdt, simplify = FALSE)
    cbind(unique(dispersal_dat[by]), mdt = unlist(mdt_values))
  }
  rownames(ans) <- NULL

  ans <- structure(
    ans,
    class = c('sit_mdt', class(ans)),
    dispersal_data = dispersal_dat
  )

  return(ans)
}


#' @importFrom knitr kable
#' @export
print.sit_mdt <- function(x, digits = 2, ...) {
  cli::cli_h1("Mean Distance Travelled")

  if (inherits(x, 'data.frame')) {
    print(
      knitr::kable(
        x,
        format = "simple",
        digits = digits
      )
    )
  } else {
    cat("\n")
    cli::cli_alert_info("Estimated average value: {round(x, digits)}")
  }
  cat("\n")
}


#' @importFrom stats quantile
#' @export
plot.sit_mdt <- function(x, y, ...) {

  dispersal_dat <- attr(x, 'dispersal_data')
  n_surveys <- nrow(dispersal_dat)

  ## Adjusted distances so that MDT = mean(adj_dists)
  adj_dists <- n_surveys *
    with(dispersal_dat, n * w / sum(n * w) * dist_m)

  quant95 <- quantile(adj_dists, probs = c(.025, .975))
  adj_dists_sel <- adj_dists[
    adj_dists > 0 & adj_dists >= quant95[1] & adj_dists <= quant95[2]
  ]

  hist(
    adj_dists_sel,
    main = "Adjusted Distances Travelled",
    xlab = "Distance (in m) adjusted by counts and weights",
    ylab = NULL,
    yaxt = 'n'
  )
  if (inherits(x, 'data.frame')) {
    by <- setdiff(names(x), 'mdt')
    if ('population' %in% by) {
      abline(v = x$mdt, col = x$population)
    } else {
      abline(v = x$mdt, col = x$age)
    }
  } else {
    abline(v = x)
  }
}

#' Age of captured sterile males
#'
#' Compute the _age_ (i.e. number of days since release to capture) of released
#' sterile males.
#'
#'
#' @param surveys A `sit_adult_surveys` object.
#' @param releases A `sit_revents` object with target releases.
#'
#' @return A data.frame, with a similar structure as `surveys` but with captures
#'   of the sterile males from the target releases only, and with additional
#'   variables for the release date and the age of the individual at the moment
#'   of the survey in number of days.
#' @export
#'
#' @examples
#' head(
#'   survey_ages(sit_adult_surveys(sit_prototype), sit_revents(sit_prototype))
#' )
survey_ages <- function(surveys, releases) {
  stopifnot(
    inherits(surveys, 'sit_adult_surveys'),
    inherits(releases, 'sit_revents')
  )

  ## Implicitly filters males, since female surveys necessarily correspond
  ## to a wild population, which is excluded since it has not been released.
  ans <- merge(
    surveys,
    setNames(
      as.data.frame(releases)[, c("colour", "date")],
      c("pop_col", "release_date")
    ),
    by = "pop_col"
  )

  ans$age <- as.numeric(
    difftime(ans$datetime_end, ans$release_date, units = "days")
  )

  return(ans)
}


#' Distances from releases to captures
#'
#' For every survey of sterile adults, return the distance from the release
#' point to the location of the corresponding trap.
#'
#' Populations from areal releases are excluded.
#'
#' @param x A `sit` object.
#'
#' @return A `data.frame` with the survey `id` and the `dist` in `m` from the
#' release point.
#' @keywords internal
#' @importFrom units set_units
#' @export
#'
#' @examples
#' survey_distances_from_release(sit_prototype)
survey_distances_from_release <- function(x) {
  geoms <- x %>%
    dm_select('adult_surveys', 'id', 'trap_id', 'population_id') %>%
    dm_select('release_events', 'id', 'site_id') %>%
    dm_select('release_sites', 'id', released = 'geometry') %>%
    ## This trivial filter selects all release events, but implicitly filters
    ## adult surveys that belong to one of the released populations (i.e.
    ## excluding surveys of wild individuals.)
    ## https://github.com/cynkra/dm/issues/706
    dm_filter(release_events = TRUE) %>%
    dm_filter(release_sites = (release_types(x[['release_sites']]) == 'point')) %>%
    dm_select('traps', 'id', 'type_id', captured = 'geometry') %>%
    dm_flatten_to_tbl(.start = 'adult_surveys', .recursive = TRUE)

  distances <- sf::st_distance(
    geoms$released, geoms$captured, by_element = TRUE
  )
  distances_m <- as.numeric(units::set_units(distances, "m"))
  ans <- data.frame(
    id = geoms$id,
    dist = distances_m
  )

  return(ans)
}


#' Inverse radial density weights
#'
#' Compute a vector of weights for a set of elements (traps), to account for
#' inhomogeneous spatial distribution relative to a reference (release) point.
#'
#' The _radial density_ of traps (i.e. number of traps in a ring ($r$, $r+dr$),
#' divided by the surface area of the ring) is proportional to the radius under
#' a spatially homogeneous distribution of traps, with a proportionality
#' constant of \eqn{2 / R^2},  where \eqn{R} is the maximum radius.
#'
#' This function performs a kernel-density estimation of the radial density
#' of the given set of traps, using a boundary correction to avoid edge effects,
#' and returns weights calculated as the inverse of the relative density, with
#' respect to the expected density under homogeneity.
#'
#' These weights can be used to adjust trap counts by spatial arrangement, in
#' situations where the relationship between counts and the distance from the
#' release point is relevant.
#'
#' @param x Numeric vector of non-negative distances.
#'
#' @return A `inverse_radial_density` object, which is a numeric vector of
#'   weights (one for each trap in `x`, in the same order), with a set of
#'   attributes used in the `plot` method.
#'
#' @importFrom stats density approx
#' @export
#'
#' @examples
#' ## Homogeneous points
#' set.seed(20211129)
#' radial_dists_homog <- sqrt(runif(50)^2 + runif(50)^2)
#' ird <- inverse_radial_density(radial_dists_homog)
#' plot(ird)
inverse_radial_density <- function(
  x
) {

  upper_boundary <- max(x) + mean(diff(sort(x)))

  ## Estimate the distance-density at each trap
  kde <- stats::density(x, from = 0, to = upper_boundary)

  ## Account for border effect
  ## Idea from Venables and Ripley (2002) MASS.
  ## Originally taken from Silverman 1986.
  x_mirror <- c(x, 2 * upper_boundary - x)

  kde2 <- density(
    x_mirror,
    from = 0,
    to = 2 * upper_boundary,
    bw = kde$bw,
    n = 2 * length(kde$x)
  )

  ## Double total density (so it integrates to 2)
  ## and keep the original half
  kde2$y <- 2 * kde2$y
  kde2$x <- kde2$x[seq_along(kde$x)]
  kde2$y <- kde2$y[seq_along(kde$y)]

  ## Interpolate density at the observed distances
  dist_density <- approx(x = kde2$x, y = kde2$y, xout = x)$y

  ## The radial density in a homogeneous arrangement is expected to increase
  ## linearly with the radii, with a slope of: 2/upper_boundary**2
  exp_density_slope <- 2/upper_boundary**2

  # ## Visualise boundary correction
  # plot(kde, col = "darkgrey")
  # lines(kde2)
  # points(x, dist_density)
  # abline(0, 2 / upper_boundary**2, col = 'lightgray')

  relative_density = dist_density / (2 * x / upper_boundary ** 2)

  w = 1 / relative_density

  ans <- structure(
    w,
    distances = x,
    dist_density = dist_density,
    expected_slope = exp_density_slope,
    density = kde2,
    uncorrected_density = kde,
    class = "inverse_radial_density"
  )

  return(ans)
}

#' @export
print.inverse_radial_density <- function(x, ...) {
  print(as.numeric(x))
}

#' @import graphics
#' @export
plot.inverse_radial_density <- function(x, y, ...) {

  ## Visualise boundary correction
  plot(
    attr(x, "uncorrected_density"),
    main = "Estimated radial density",
    col = "darkgrey",
    ...
  )
  lines(attr(x, "density"))
  points(attr(x, "distances"), attr(x, "dist_density"))
  abline(0, attr(x, "expected_slope"), col = 'lightgray')


}

#' Flight Range
#' @family estimates
#'
#' Flight range of sterile individuals at specified levels.
#'
#' @inheritParams dispersal
#' @param pool Logical. Whether to aggregate data from released populations.
#' @param levels Numeric vector with values between 0 and 100. Requested levels
#'   of flight range, in percentage.
#'
#' @return A named numeric vector with one value for each requested `level`.
#' @export
#'
#' @examples
#' sit_flight_range(sit_prototype)
sit_flight_range <- function(
  x,
  pool = FALSE,
  levels = c(50, 90),
  spatial_adjustment = TRUE,
  following_releases,
  following_days,
  species = NULL
) {

  dispersal_dat <- dispersal(
    x,
    spatial_adjustment = spatial_adjustment,
    following_releases = following_releases,
    following_days = following_days,
    species = species
  )

  ## Flight range from a data.frame
  fr_df <- function(x, levels) {
    ## Always adjust, since if spatial_adjustment = FALSE, weights are 1.
    adjusted_counts <- with(x, n * w)
    flight_range(adjusted_counts, x$dist_m, levels = levels)
  }

  ans <- if (pool) {
    ans <- fr_df(dispersal_dat, levels = levels)
    structure(
      data.frame(
        population = "pooled",
        ans
      ),
      fr_data = data.frame(population = "pooled", attr(ans, 'fr_data'))
    )
  } else {
    ans_by_pop <-
      lapply(
        split(dispersal_dat, dispersal_dat$population),
        fr_df,
        levels = levels
      )
    ans <- data.frame(
      population = rep(names(ans_by_pop), each = length(levels)),
      do.call(rbind, ans_by_pop)
    )
    rownames(ans) <- NULL

    fr_dats <- lapply(ans_by_pop, attr, 'fr_data')
    fr_dat <- data.frame(
      population = rep(names(ans_by_pop), times = sapply(fr_dats, nrow)),
      do.call(rbind, fr_dats)
    )
    rownames(fr_dat) <- NULL

    structure(ans, fr_data = fr_dat)
  }

  structure(ans, class = c("sit_flight_range", "data.frame"))

}


#' @importFrom knitr kable
#' @export
print.sit_flight_range <- function(x, digits = 2, ...) {
  cli::cli_h1("Flight range")

  if (inherits(x, 'data.frame')) {
    print(
      knitr::kable(
        x,
        format = "simple",
        digits = digits
      )
    )
  } else {
    cat("\n")
    cli::cli_alert_info("Estimated average value: {round(x, digits)}")
  }
  cat("\n")
}

#' @export
plot.sit_flight_range <- function(x, y, ...) {

  fr_dat <- attr(x, 'fr_data')

  frm <- lm(logd ~ cx*population, data = fr_dat)

  fr_dat$pred <- predict(frm)

  plot(
    logd ~ cx,
    data = fr_dat,
    pch = 19,
    col = fr_dat$population,
    main = "Flight Ranges",
    xlab = "Cumulative captures",
    ylab = "Log-distance (log-m)"
  )

  for (d in split(fr_dat, fr_dat$population)) {
    lines(d$cx, d$pred, col = d$population)
  }

  invisible(fr_dat)
}



#' Diffusion
#' @family estimates
#'
#' Estimate the diffusion coefficient.
#'
#' Assuming that mosquitoes fly randomly (Brownian motion), their dispersion
#' follow [Fick's first
#' law](https://en.wikipedia.org/wiki/Fick%27s_laws_of_diffusion#Fick's_first_law)
#' of diffusion, which postulates that flux goes from regions of high
#' concentration to regions of low concentration, with a magnitude that is
#' proportional to the concentration gradient (spatial derivative). The
#' proportionality constant is the diffusivity $D$ and depend of the species'
#' behaviour for that location, season and weather conditions.
#'
#' Under this model, [the __Mean Squared Displacement__ (MSD) of mosquitoes from
#' their release point at time $t$
#' is](https://en.wikipedia.org/wiki/Fick%27s_laws_of_diffusion#Example_solution_2:_Brownian_particle_and_Mean_squared_displacement)
#' \deqn{\text{MSD}(t) = 4Dt.}
#'
#' From which we can solve for the Diffusion coefficient by using the distances
#' from the release point at which sterile males were captured to estimated the
#' MSD at a given age.
#'
#' This function computes the MSD of the individuals (possibly, by population)
#' at the different observed ages, and derives $D$ as the regression slope.
#' The resulting object can be plotted to see the underlying regression model.
#'
#' @inheritParams dispersal
#' @param pool Logical. Whether to aggregate data from released populations.
#'
#' @return A named numeric vector with one value for each requested `level`.
#' @importFrom stats coef
#' @export
#'
#' @examples
#' (diff_prototype <- sit_diffusion(sit_prototype))
#' plot(diff_prototype)
sit_diffusion <- function(
  x,
  pool = FALSE,
  spatial_adjustment = TRUE,
  following_releases,
  following_days,
  species = NULL
) {

  dispersal_dat <- dispersal(
    x,
    spatial_adjustment = spatial_adjustment,
    following_releases = following_releases,
    following_days = following_days,
    species = species
  )

  if (pool) dispersal_dat$population <- "pooled"

  dispersal_split <- split(
    dispersal_dat,
    list(dispersal_dat$population, dispersal_dat$age)
  )

  ## Diffusion from a single data.frame
  diff_df <- function(x) {
    N <- sum(x$n)
    ## Always adjust, since if spatial_adjustment = FALSE, weights are 1.
    sum_d2 <- with(x, sum(dist_m**2 * n * w))
    MSD <- sum_d2 / N
    data.frame(N, sum_d2, MSD)
  }

  diff_dat <- data.frame(
    population = vapply(
      strsplit(names(dispersal_split), split = "\\."), `[`, character(1), 1
    ),
    age = as.integer(
      vapply(
        strsplit(names(dispersal_split), split = "\\."), `[`, character(1), 2
      )
    ),
    do.call(rbind,lapply(dispersal_split, diff_df)),
    row.names = NULL
  )

  diff_dat <- diff_dat[diff_dat$N > 0, ]
  diff_dat$`4t` <- 4 * diff_dat$age

  diff_split <- split(diff_dat, diff_dat$population)
  estimate_diff <- function(x) {
    unname(coef(lm(MSD ~ 0 + `4t`, data = x)))
  }
  ans <- data.frame(
    population = names(diff_split),
    Dest = vapply(diff_split, estimate_diff, 1),
    row.names = NULL
  )

  structure(
    ans,
    diffusion_data = diff_dat,
    class = c("sit_diffusion", "data.frame")
  )

}


#' @importFrom knitr kable
#' @export
print.sit_diffusion <- function(x, digits = 2, ...) {
  cli::cli_h1("Diffusion coefficient")

  if (inherits(x, 'data.frame')) {
    print(
      knitr::kable(
        x,
        format = "simple",
        digits = digits
      )
    )
  } else {
    cat("\n")
    cli::cli_alert_info("Estimated average value: {round(x, digits)}")
  }
  cat("\n")
}

#' @export
plot.sit_diffusion <- function(x, y, ...) {

  diff_dat <- attr(x, 'diffusion_data')

  diffm <- lm(MSD ~ 0 + `4t`:population, data = diff_dat)

  diff_dat$pred <- predict(diffm)

  plot(
    MSD ~ `4t`,
    data = diff_dat,
    pch = 19,
    col = diff_dat$population,
    main = "Diffusion coefficient",
    xlab = "Scaled time",
    ylab = "MSD"
  )

  for (d in split(diff_dat, diff_dat$population)) {
    lines(d$`4t`, d$pred, col = d$population)
  }

  invisible(diff_dat)
}



#' Flight range
#'
#' Compute the flight range at specified levels.
#'
#'
#' @param x Numeric vector of (possibly adjusted) counts.
#' @param d Numeric vector of distances, of the same length as `x`.
#' @param levels Numeric vector of quantiles, as a percentage, between 0 and
#'   100.
#'
#' @return A numeric vector of the same length as `levels`.
#' @keywords internal
#' @importFrom stats lm predict setNames
#' @export
#'
#' @examples
#' flight_range(rep(50, 100), 100*runif(100), levels = c(50, 90))
flight_range <- function(x, d, levels = c(50, 90)) {

  ## Summarise values of x (catches) by values of d (distance)
  dat <- aggregate(x, by = list(d = d), FUN = sum)
  # ord_d <- order(d)
  # dat <- data.frame(
  #     cx = cumsum(x[ord_d]),
  #     logd = log(d[ord_d] + 1)
  #   )

  dat$cx <- cumsum(dat$x)
  dat$logd <- log(dat$d + 1)

  fr_model <- lm(logd ~ cx, data = dat)

  cum_adj_catch = sum(x) * levels / 100

  fr <- predict(fr_model, newdata = data.frame(cx = cum_adj_catch))

  ans <- structure(
    data.frame(
      level = levels,
      cum_adj_catch = cum_adj_catch,
      FR = exp(fr) - 1
    ),
    fr_data = dat,
    class = c('flight_range', 'data.frame')
  )

  return(ans)
}



#' Dispersal data
#'
#' Compute a `data.frame` with dispersal information. I.e. age and distance
#' from release point for each capture. Possibly adjusting for inhomogeneous
#' trap arrangement.
#'
#' By default, the function uses all adult surveys in the `sit` area of captured
#' individuals released from __point__ releases. You can control this behaviour
#' by passing a custom `sit_revents` object in `following_releases`.
#'
#' Note that the spatial adjustment is performed with the provided traps.
#' Thus, beware of using a subset of traps only.
#'
#' @param x A `sit` object.
#' @param spatial_adjustment Logical. Whether to adjust observations to account
#'   for a irregular spatial arrangement of traps.
# @param weights Numeric vector of length equal to the number of adult traps in
#   the SIT area. Multiplicative adjustment of the catches at each trap. By
#   default, `inverse_density(sit_traps(x, stage = 'adult', area = 'sit'))`.
#' @param following_releases A `sit_release` object with a subset of release
#'   events or missing (default) for all point release events. Use in
#'   combination with `following_days` to filter surveys of the target
#'   populations within a given number of days after release. Note that counts
#'   of wild populations will always be included in the results with a value of
#'   `NA` in `pop_col`.
#' @param following_days Integer or missing (default). Number of days after
#'   releases to return, if `following_releases` is not missing.
#' @param species a character vector of species to be returned. Defaults to
#'   `NULL`, which ignores the species variable.
#'
#' @return A `data.frame` with dispersal data. Specifically, for each survey of
#'   sterile males, their `population` (colour), the `trap`, the `age`, the
#'   capture size `n`, the spatial weights `w` and the distance `dist_m` from
#'   the release point, in metres.
#' @export
#'
#' @examples
#' dispersal(sit_prototype)
#'
#' yellow_release <- sit_revents(sit_prototype)[1,]
#' dispersal(
#'   sit_prototype,
#'   following_releases = yellow_release,
#'   following_days = 7L
#' )
dispersal <- function(
  x,
  spatial_adjustment = TRUE,
  following_releases,
  following_days,
  species = NULL
) {

  stopifnot(inherits(x, 'sit'))
  as <- sit_adult_surveys(
    x,
    area = 'sit',
    following_releases = following_releases,
    following_days = following_days,
    species = species
  )

  ## Extract point releases from x or from following releases if given.
  if (missing(following_releases)) following_releases <- x
  target_releases <- sit_revents(following_releases, type = 'point')

  if (nrow(target_releases) == 0L) {
    stop("No point releases found. Cannot compute dispersal.")
  }

  ## Ages of surveyed __males__ from the populations of the target releases
  ## Excludes observations in as from wild or areal releases.
  sm_ages <- survey_ages(as, target_releases)

  ## Distances from release to capture
  ## Captures from areal releases are excluded.
  dist_release <- survey_distances_from_release(x)

  ## Last age for which there is at least one observation of a sterile male
  max_age <- max_observed_age(sm_ages$n, sm_ages$age)

  ## Include distances and filter ages
  dispersal_dat <- merge(
    sm_ages[sm_ages$age <= max_age, c('id', 'pop_col', 'trap_code', 'age', 'n')],
    dist_release,
    by = 'id'
  )

  ## There should not be any missing value in the distances
  stopifnot(all(!is.na(dispersal_dat$dist)))


  ## Spatial adjustment
  if (spatial_adjustment) {
    trap_pops <- unique(dispersal_dat[, c('pop_col', 'trap_code', 'dist')])

    n_traps_per_release <- table(trap_pops$pop_col)
    if (any(n_traps_per_release < 5)) stop(
      "Too few traps to perform spatial adjustment. ",
      "Please call `sit_mdt` with `sptatial_adjustment = FALSE`."
    )

    ## Compute trap weights, by release
    add_w <- function(x) {
      x$w <- as.numeric(inverse_radial_density(x$dist))
      return(x)
    }
    trap_pops <-
      unsplit(
        lapply(split(trap_pops, trap_pops$pop_col), add_w),
        trap_pops$pop_col
      )
    dispersal_dat <- merge(
      dispersal_dat, trap_pops, by = c('pop_col', 'trap_code', 'dist')
    )
  } else {
    dispersal_dat$w <- 1
  }

  ## Variable selection and renaming
  sel_vars <- c(
    population = 'pop_col',
    trap = 'trap_code',
    age = 'age',
    n = 'n',
    w = 'w',
    dist_m = 'dist'
  )
  ord <- with(dispersal_dat, order(pop_col, age))
  dispersal_dat <- dispersal_dat[ord, sel_vars]  # remove survey 'id'
  names(dispersal_dat) <- names(sel_vars)

  return(dispersal_dat)
}

#' Size of the Wild Population
#'
#' Estimate the size of the wild population, via the Lincoln Index.
#'
#' This provides multiple estimates of the population size, one for each age $t$
#' and for each released population, which can be averaged to compute a final
#' estimate based on all the available data.
#'
#' The calculation needs to be split by released population. First, because
#' survival rates could be different among populations, but also because the
#' number $n_t$ of captured wild males at a given date can be compared against
#' the number of captured sterile males $m_t$ at different ages.
#'
#' @inheritParams sit_flight_range
#'
#' @return A `sit_wild_size` object.
#' @importFrom stats sd
#' @export
#'
#' @examples
#' sit_wild_size(sit_prototype)
sit_wild_size <- function(
  x,
  pool = FALSE,
  species = NULL
) {

  ## Survival
  ## Using spatially-adjusted captures by default here.
  ## But this is irrelevant for the survival parameters.
  surv <- sit_survival(
    x,
    spatial_adjustment = FALSE,  # This should be removed from sit_survival
    pool = pool,
    species = species)[, c("population", "N_released", "SR")]



  ## Adult captures by population and age
  # ## Aggregated captures by population and age
  # age_captures <- attr(surv, "surv_data")

  ## Age surveys from areal releases as well, and without spatial adjustment.
  ages <- survey_ages(sit_adult_surveys(x, area = "sit"), sit_revents(x))

  if (pool) ages$pop_col <- "pooled"

  ## Remove ages for which there are no observations
  ## Last age for which there is at least one observation of a sterile male
  max_age <- max_observed_age(ages$n, ages$age)
  ages <- ages[ages$age <= max_age, ]
  names(ages) <- gsub("^pop_col$", "population", names(ages))

  ages <- aggregate(n ~ population + age + datetime_end, data = ages, sum)


  ## Wild counts aggregated across traps
  asurvs <- sit_adult_surveys(x, area = "sit")
  wild_survs <- asurvs[asurvs$pop_col == "wild" & asurvs$sex == "male", ]

  if (!is.null(species) && species %in% names(wild_survs)) {
    wild_survs <- wild_survs[wild_survs$species %in% species, ]
  }

  wild_survs <- wild_survs[, c("trap_code", "datetime_end", "n")]
  names(wild_survs) <- gsub("^n$", "n_w", names(wild_survs))
  wild_survs <- aggregate(n_w ~ datetime_end, data = wild_survs, sum)

  ages <- merge(ages, surv, by = "population")
  pop_dat <- merge(ages, wild_survs, by = c("datetime_end"))

  ## Lincoln index
  pop_dat$LI <- with(
    pop_dat,
    lincoln_index(
      released = N_released,
      survival_rate = SR,
      age = age,
      marked_recaptured = n,
      wild_catch = n_w
    )
  )

  ## Rename columns
  sel_cols <- c(population = "population", survey_date = "datetime_end",
                age = "age", N_released = "N_released", SR = "SR", n = "n",
                n_wild = "n_w", LI = "LI")
  pop_dat <- setNames(pop_dat[, sel_cols], names(sel_cols))

  ans <- aggregate(LI ~ population, data = pop_dat, function(x) c(p_est = mean(x), p_sd = sd(x)))
  ans <- data.frame(population = ans$population, as.data.frame(ans$LI))

  return(
    structure(
      ans,
      pop_data = pop_dat,
      class = c("sit_wild_size", "data.frame")
    )
  )

}

#' @importFrom knitr kable
#' @export
print.sit_wild_size <- function(x, digits = 2, ...) {
  cli::cli_h1("Wild male population size")

  if (inherits(x, 'data.frame')) {
    print(
      knitr::kable(
        x,
        format = "simple",
        digits = digits
      )
    )
  } else {
    cat("\n")
    cli::cli_alert_info("Estimated average value: {round(x, digits)}")
  }
  cat("\n")
}

#' @export
plot.sit_wild_size <- function(x, y, ...) {

  pop_dat <- attr(x, 'pop_data')

  boxplot(
    LI ~ population,
    data = pop_dat,
    horizontal = TRUE,
    col = rev(unique(pop_dat$population)),
    main = "Wild male population size",
    xlab = "Lincoln index",
    ylab = "Released population"
  )

  invisible(pop_dat)
}




#' Lincoln index
#'
#' Estimate the size of the wild-male population.
#'
#' A simple estimate is obtained as follows (Thompson 2012, Ch. 18). Let the
#' total captures at a day $t$ be the sum of $m_t$ marked and $n_t$ wild
#' mosquitoes. Assuming that the proportion of marked individuals in the sample
#' equals that in the population of size $P$: \deqn{\frac{m_t}{m_t + n_t} =
#' \frac{M_t}{M_t + P},} where $M_t = R\,S^{a_t}$ is the number of marked
#' individuals captured at time $t$, with $R$ the number of released adults, $S$
#' the daily survival rate and $a_t$ the number of days since release (_age_).
#' I.e., the number of marked individuals at time $t$ is the remaining number
#' from those released that survived for $a_t$ days.
#'
#' The __Lincoln Index__ (_a.k.a._ the __Petersen estimator__) has been used as
#' a simple estimate of the wild __male__ population size, assuming that the
#' survival rate of an individual remains constant. Here we use a _modified_
#' version that corrects for small samples and compensates for daily survival.
#' \deqn{P_t = R\,S^{a_t}\,(n_t + 1) / (m_t + 1).}
#'
#' The values of $R$, $n_t$, $m_t$ and $t$ can be gathered from the adult
#' surveys data. The calculation required the estimation of the survival rate
#' $S$.
#'
#' @param released Numeric. Number of sterile males released.
#' @param survival_rate Numeric. Estimated survival rate.
#' @param age Numeric. Number of days since release.
#' @param marked_recaptured Numeric. Number of re-captured sterile males.
#' @param wild_catch Numeric. Number of wild males captured.
#'
#' @export
#' @references Thompson, Steven K. 2012. Sampling. 3rd ed. Wiley Series in
#'   Probability and Statistics. Hoboken, N.J: Wiley.
#' @examples lincoln_index(
#'   released         = 1e4,
#'   survival_rate    = .80,
#'   age              = 1,
#'   marked_recapture = 15,
#'   wild_catch       = 5
#' )
lincoln_index <- function(
  released,
  survival_rate,
  age,
  marked_recaptured,
  wild_catch
) {
  released * survival_rate^age * (wild_catch + 1) / (marked_recaptured + 1)
}

#' Survival statistics
#'
#' Compute the _Probability of Daily Survival_ (PDS), the _Average Life
#' Expectancy_ `ALE`, _Recapture Rate_ `RR`, _Survival Rate_ `SR`), by
#' population (by default, unless `pool = TRUE`)
#'
#'
#' @inheritParams sit_flight_range
#'
#' @return A `sit_survival` object, which is a data.frame with additional
#' attributes which can be plotted.
#' @export
#'
#' @examples
#' (surv_prototype <- sit_survival(sit_prototype))
#' plot(surv_prototype)
sit_survival <- function(
  x,
  pool = FALSE,
  spatial_adjustment = TRUE,
  # weights = inverse_density(sit_traps(x, stage = 'adult', area = 'sit')),
  following_releases,
  following_days,
  species = NULL
) {

  ## TODO: I think I don't really need the dispersal data, but only counts
  ## aggregated by population and age. Check this.
  dispersal_dat <- dispersal(
    x,
    spatial_adjustment = spatial_adjustment,
    following_releases = following_releases,
    following_days = following_days,
    species = species
  )


  populations <- unique(dispersal_dat$population)
  N <- sit_revents(x)$n[match(populations, sit_revents(x)$colour)]
  if (pool) {
    dispersal_dat$population <- "pooled"
    N <- sum(N)
  }

  surv_dat <- aggregate(n ~ population + age, data = dispersal_dat, FUN = sum)
  surv_dat$log_catch = log(1 + surv_dat$n)

  probability_daily_survival <- function(x) {
    slope <- coef(x)[2]
    exp(slope)
  }

  recapture_rate <- function(x, N) {
    alpha <- coef(x)[1]
    exp(alpha) / (N + exp(alpha))
  }

  survival_rate <- function(x, rr, d = 1) {
    beta <- coef(x)[2]
    exp(beta) / ( (1-rr)^(1/d))
  }

  fm_pop <- lapply(
    split(surv_dat, surv_dat$population),
    function(x) lm(log_catch ~ age, data = x)
  )

  pds <- min(vapply(fm_pop, probability_daily_survival, numeric(1)), 1 - .Machine$double.eps)
  rr <- mapply(recapture_rate, fm_pop, N)

  ans <- data.frame(
    population = names(fm_pop),
    N_released = N,
    PDS = pds,
    ALE = -1 / log(pds),
    RRx1e3  = rr * 1e3,
    SR  = mapply(survival_rate, fm_pop, rr),
    row.names = NULL
  )

  return(
    structure(
      ans,
      surv_data = surv_dat,
      class = c("sit_survival", "data.frame")
    )
  )
}


#' @importFrom knitr kable
#' @export
print.sit_survival <- function(x, digits = 2, ...) {
  cli::cli_h1("Survival")

  if (inherits(x, 'data.frame')) {
    print(
      knitr::kable(
        x,
        format = "simple",
        digits = digits
      )
    )
  } else {
    cat("\n")
    cli::cli_alert_info("Estimated average value: {round(x, digits)}")
  }
  cat("\n")
}

#' @export
plot.sit_survival <- function(x, y, ...) {

  surv_dat <- attr(x, 'surv_data')

  model_formula <- log_catch ~ age*population

  n_populations <- nrow(x)

  if(n_populations == 1L) {
    model_formula <- log_catch ~ age
    surv_dat$population <- "black"
  }

  survm <- lm(model_formula, data = surv_dat)

  surv_dat$pred <- predict(survm)

  plot(
    log_catch ~ age,
    data = surv_dat,
    pch = 19,
    col = surv_dat$population,
    main = "Survival",
    xlab = "Age (days since release)",
    ylab = "Log-catch"
  )

  for (d in split(surv_dat, surv_dat$population)) {
    lines(d$age, d$pred, col = d$population)
  }

  invisible(surv_dat)
}


#' Import or extract release events
#'
#' Provide a `sf` object with a POINT at each record corresponding to a
#' __release point__ or a `data.frame` for a __areal release__, specifying which
#' variables contain the required information. Or, extract the release events
#' used in a `sit` object.
#'
#' @param x Object of class `sf` with POINT elements to be _imported_, or a
#'   `data.frame` with information about areal releases, or
#'   object of class `sit` to _extract_ release events from.
#' @param date Character string representing a date and optionally the time, in
#'   the [RFC 3339 format](https://www.ietf.org/rfc/rfc3339.txt) (a variation of
#'   the [ISO 8601 format](https://en.wikipedia.org/wiki/ISO_8601)), i.e.,
#'   `2021-12-31` or `2019-11-23 15:00`.
#' @param colour Character string. Colour of the release.
#' @param species Optional character. Released species.
#' @param n Numeric. Natural number. Number of individuals released.
#' @param ... Used to pass arguments to specific method.
#'
#' @return A object of class `sit_revents` which can be used in [sit()].
#'   `sit_revents` objects can be concatenated together with `c()`, in order to
#'   import a mixture of point and areal releases into a `sit` object. See
#'   examples.
#'
#' @import sf
#' @export
#'
#' @family importing
#' @seealso [sit()]
#'
#' @examples
#' point_releases <- sit_revents(
#'   sf::st_as_sf(
#'     data.frame(
#'       x = 1:3,
#'       y = 3:1,
#'       date = c("2019-11-25", "2019-12-01", "2019-12-13"),
#'       colour = c("yellow", "red", "blue"),
#'       n = 1e4
#'     ),
#'     coords = c("x", "y")
#'   )
#' )
#'
#' areal_release <- sit_revents(
#'   data.frame(
#'     x = 1, y = 1,
#'     date = "2019-12-21",
#'     colour = "pink",
#'     n = 1e4
#'   )
#' )
#'
#' (c(point_releases, areal_release)) # also a `sit_revent` object
#'
sit_revents <- function(x, ...) UseMethod("sit_revents")


#' @describeIn sit_revents Imports point release data.
#' @export
sit_revents.sf <- function(
  x,
  date = "date",
  colour = "colour",
  n = "n",
  species = NULL,
  ...
) {

  # The features must be all POINTs
  stopifnot(st_geometry_type(x, by_geometry = FALSE) == "POINT")

  ## Build the internal representation
  ans <- build_internal_df(
    x,
    varnames = c(date, colour, n, species),
    internal_names = c("date", "colour", "n", "species"[!is.null(species)])
  )

  ## Interpret dates
  ans$date <- as.POSIXct(ans$date)

  ## Check duplicated releases
  if (any(idx <- duplicated(ans))) {
    stop(
      glue(
        "{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
        "duplicated. Please fix."
      )
    )
  }

  ## Add unique internal numeric identifier of the release event
  ans <- cbind(id = seq.int(nrow(x)), ans)

  ## Table of __unique__ release_sites
  sites <- st_geometry(ans)
  unique_sites <- sites[!duplicated(sites)]
  release_sites <- st_sf(id = seq_along(unique_sites), geometry = unique_sites)

  ## Replace the geometries by the corresponding site_id
  ans <- st_drop_geometry(ans)
  ans$site_id <- as.integer(st_equals(sites, unique_sites))

  ## Reorder columns
  ans <- ans[
    , c("id", "site_id", "date", "colour", "species"[!is.null(species)], "n")
  ]

  attr(ans, "release_sites") <- release_sites
  class(ans) <- c("sit_revents", class(ans))
  return(ans)
}

#' @describeIn sit_revents Imports areal release data.
#' @export
sit_revents.data.frame <- function(
  x,
  date = "date",
  colour = "colour",
  n = "n",
  species = NULL,
  ...
) {

  ## Build the internal representation
  ans <- build_internal_df(
    x,
    varnames = c(date, colour, n, species),
    internal_names = c("date", "colour", "n", "species"[!is.null(species)])
  )

  ## Interpret dates
  ans$date <- as.POSIXct(ans$date)

  ## Check duplicated releases
  if (any(idx <- duplicated(ans))) {
    stop(
      glue(
        "{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
        "duplicated. Please fix."
      )
    )
  }

  ## Add unique internal numeric identifier of the release event
  ans <- cbind(id = seq.int(nrow(x)), ans)

  ## Table of __unique__ release_sites
  ## All areal releases are performed from "the same" release_site
  ## The geometry is a POINT EMPTY.
  release_sites <- st_sf(
    id = 1L,
    geometry = st_sfc(st_point()),
    crs = NA_crs_
  )

  ## Replace the geometries by the corresponding site_id
  ans$site_id <- 1L

  ## Reorder columns
  ans <- ans[
    , c("id", "site_id", "date", "colour", "species"[!is.null(species)], "n")
  ]

  attr(ans, "release_sites") <- release_sites
  class(ans) <- c("sit_revents", class(ans))
  return(ans)
}

#' @describeIn sit_revents Concatenate release events.
#' @export
c.sit_revents <- function(...) {

  lst <- list(...)
  classes_revents <- vapply(lst, inherits, TRUE, 'sit_revents')

  ## All elements need to be of class `sit_revents`
  if (!all(classes_revents)) {
    stop(
      glue(
        "{vars_collapse(
          seq_along(lst)[!classes_revents], 'Argument ', 'Arguments ',
          quote = ''
        )} need to be created with `sit_revents()`."
      )
    )
  }

  ## row-bind release sites
  rs <- lapply(lst, attr, 'release_sites')
  ## Check whether all release sites are areal
  ## I.e. have empty points and NA bboxes
  areals <- vapply(lapply(rs, st_bbox), is.na, TRUE)
  all_areals <- all(areals)
  if (all_areals) {
    ## If all_areals, suppressWarnings related to establishing the bbox
    ## and reset resulting bbox as NA
    release_sites <- suppressWarnings(combine_rows(rs))
    attr(release_sites[[attr(release_sites, "sf_column")]], "bbox") <-
      ## The bbox does not need to have a crs. This causes spurious warnings.
      structure(NA_bbox_, crs = NULL)
  } else {

    ## Check compatibility of CRSs
    crs_list <- lapply(rs, st_crs)
    if (length(crs <- unique(crs_list[!areals])) > 1) {
      stop("Arguments have different coordinate reference systems. ",
           "Please use st_transform() to convert to a common CRS.")
    }

    ## Set areal CRSs to the common CRS
    for (i in which(areals)) {
      st_crs(rs[[i]]) <- crs[[1]]
    }

    release_sites <- combine_rows(rs)
  }


  ## recode site_ids
  element_id <- rep(seq_along(rs), times = vapply(rs, nrow, 1))
  new_codes <- split(seq_along(element_id), element_id)
  for (i in seq_along(lst)) {
    lst[[i]]$site_id <- new_codes[[i]][lst[[i]]$site_id]
  }

  ## row-bind recoded tables of release events
  ## and reset event id
  ans <- combine_rows(lst)

  ## Attach attribute table of release sites
  attr(ans, 'release_sites') <- release_sites

  return(ans)
}


#' @describeIn sit_revents Extracts release events.
#' @param type Character vector. Which _types_ of release events to retrieve.
#'   Either `point` or `areal` or both (default).
#' @export
sit_revents.sit <- function(x, type = c('point', 'areal'), ...) {

  ans <-
    x %>%
    ## Filter release event types
    dm_filter(
      release_sites = (release_types(x[['release_sites']]) %in% .env$type)
    )


  ans <- structure(
    ans[['release_events']],
    release_sites = ans[['release_sites']]
  ) %>%
    set_class('sit_revents')


  return(ans)
}

#' @export
sit_revents.sit_revents <- function(x, type = c('point', 'areal'), ...) {

  sites <- attr(x, "release_sites")
  sites_idx <- release_types(sites) %in% type

  ans <- x[x$site_id %in% sites$id[sites_idx], ]
  attr(ans, "release_sites") <- sites[sites_idx, ]

  return(ans)
}

#' @export
print.sit_revents <- function(
  x,
  ...
) {

  ans <- as.data.frame.sit_revents(x)
  print(ans)
  invisible(ans)
  # ans
}


#' @export
as.data.frame.sit_revents <- function(x, ...) {

  rs <- attr(x, 'release_sites')
  names(rs) <- sub("id", "site_id", names(rs))

  ## Release type
  rs$type <- release_types(rs)

  ## Unclass
  class(x) <- setdiff(class(x), 'sit_revents')

  ans <- merge(
    x, rs, by = "site_id", sort = FALSE
  )

  cols <- c('id', 'type', setdiff(names(ans), c('id', 'type')))
  ans <- ans[order(ans$id), cols]

  return(ans)
}

#' @export
st_as_sf.sit_revents <- function(x, ...) {

  st_as_sf(as.data.frame(x))
}


#' @export
summary.sit_revents <- function(object, ...) {

  rs <- attr(object, 'release_sites')
  n_areal <- sum(sf::st_is_empty(rs))
  n_point <- nrow(rs) - n_areal

  ans <- structure(
    list(
      sit_revents = object,
      release_sites = rs,
      n = c(point = n_point, areal = n_areal)
    ),
    class = "summary.sit_revents"
  )

  return(ans)
}

#' @export
print.summary.sit_revents <- function(x, ...) {

  cat("Number of point releases:", x$n[["point"]], "\n")
  cat("Number of areal releases:", x$n[["areal"]], "\n")
  cat("\n")

  cat("Release events:\n")
  print(x$sit_revents)
  cat("\n")

  cat("Release sites:\n")
  print(x$release_sites)
}

#' Import or extract trap information
#'
#' Provide a `sf` object with a POINT at each record corresponding to a trap
#' location, specifying which variables contain the required information. Or
#' define a set of traps without spatial coordinates. Or, extract the trap data
#' used in a `sit` object.
#'
#' @param x Object of class `sf` with POINT elements to be _imported_, or a
#'   character-coerced vector _trap codes_ to be defined, or an object of class
#'   `sit` to _extract_ trap information from.
#' @param id Character. Variable name with a custom code for the trap. Jointly
#'   with `type`, identifies the trap uniquely.
#' @param type Character. Variable name with __valid__ values of trap types
#'   (i.e., any of `trap_types$label`). For the extraction method, a character
#'   vector of trap types to filter upon.
#' @param area Character. Variable name with values either `control` or `sit`.
#'   For the extraction method, a character vector of areas to filter upon.
#' @param stage Character vector of stages (i.e., `adult` or `egg`) to filter
#'   upon extraction.
#' @param label Optional character. Variable name with custom text associated to
#'   the traps.
#' @param trap_types Table of trap types, created with [sit_trap_types()].
#' @param ... Used to pass arguments to specific method.
#'
#' @return A object of class `sit_traps` which can be used in [sit()].
#' @import sf
#' @export
#'
#' @family importing
#' @seealso [sit()]
#'
#' @examples
#'
#' ## Build spatialised traps from a `sf` object
#' traps_sf <- sf::st_as_sf(
#'   data.frame(
#'     x = 1:3,
#'     y = 3:1,
#'     "Trap.Id" = letters[1:3],
#'     Type = "BGS",
#'     area = "sit"
#'   ),
#'   coords = c("x", "y")
#' )
#' sit_traps(traps_sf, id = "Trap.Id", type = "Type", area = "area")
#'
#' ## Build non-spatialised traps from a vector of codes
#' sit_traps(paste0("C", 1:5), area = 'control', type = 'OVT')
#'
#' ## Retrieve traps in a `sit` object
#' ## sit_traps(sit_prototype)
sit_traps <- function(x, ...) UseMethod("sit_traps")


#' @describeIn sit_traps Imports spatialised trap data
#' @export
sit_traps.sf <- function(
  x,
  id = "id",
  type = "type",
  area = "area",
  label,
  trap_types = sit_trap_types(),
  ...
) {

  ## Fill-in label variable if missing
  if(missing(label)) {
    if (!"label" %in% names(x)) {
      x$label <- NA
    }
    label = "label"
  }

  ## Check trap_types
  if (!inherits(trap_types, 'sit_trap_types')) {
    stop("Please use `sit_trap_types()` to create a valid table of trap types.")
  }


  ## Build the internal representation
  ans <- build_internal_df(
    x,
    varnames = c(id, type, area, label),
    internal_names = c("code", "type_id", "area", "label")
  )

  ## Check duplicated traps
  if (any(idx <- duplicated(ans))) {
    stop(
      glue(
        "{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
        "duplicated. Please fix."
      )
    )
  }

  ## Add unique internal numeric identifier
  ans <- cbind(id = seq.int(nrow(x)), ans)

  ## Check trap types and replace by codes
  obs_types <- unique(ans$type_id)
  xtypes <- setdiff(obs_types, trap_types$label)  # invalid types
  if (length(xtypes) > 0) {
    stop(
      glue(
        "Invalid trap {vars_collapse(xtypes, 'type ', 'types ')} ",
        "in variable {vars_collapse('type', '')}.
        Valid types are the values of 'label' ",
        "in the object passed to `trap_types`, ",
        "i.e.: {vars_collapse(trap_types$label, '', '')}. ",
        "Please verify."
      )
    )
  }
  ans$type_id <- trap_types$id[match(ans$type_id, trap_types$label)]

  attr(ans, "trap_types") <- trap_types
  class(ans) <- c("sit_traps", class(ans))
  return(ans)
}


#' @describeIn sit_traps Imports non-spatialised trap data
#' @export
sit_traps.character <- function(
  x,
  type,
  area,
  trap_types = sit_trap_types(),
  ...
) {

  ## Build a `sf` representation with POINT EMPTY geometry
  dat_sf <- st_sf(
    id   = x,
    type = rep(type, length(x)),
    area = rep(area, length(x)),
    geometry = rep(st_sfc(st_point()), length(x)),
    ...,
    crs = NA_crs_
  )

  ans <- sit_traps.sf(dat_sf, trap_types = trap_types)

  return(ans)
}


#' @describeIn sit_traps Imports non-spatialised trap data
#' @export
sit_traps.numeric <- function(x, ...) {
  sit_traps.character(as.character(x), ...)
}


#' @describeIn sit_traps Combine traps.
#' @export
c.sit_traps <- function(...) {

  lst <- list(...)
  classes_traps <- vapply(lst, inherits, TRUE, 'sit_traps')

  ## All elements need to be of class `sit_traps`
  if (!all(classes_traps)) {
    stop(
      glue(
        "{vars_collapse(
          seq_along(lst)[!classes_traps], 'Argument ', 'Arguments ',
          quote = ''
        )} need to be created with `sit_traps()`."
      )
    )
  }

  ## Combine trap types.
  tt <- lapply(lst, attr, 'trap_types')
  trap_types <- combine_rows(tt, recoder = TRUE)

  ## Recode trap type ids
  for (i in seq_along(lst)) {
    lst[[i]]$type_id <- attr(trap_types, 'recoder')(lst[[i]]$type_id, i)
  }
  attr(trap_types, 'recoder') <- NULL


  ## Check whether all traps are non-spatial
  ## I.e. have empty points and NA bboxes
  non_spatial <- vapply(lapply(lst, st_bbox), is.na, TRUE)
  all_non_spatial <- all(non_spatial)
  if (all_non_spatial) {
    ## If all_non_spatial, suppressWarnings related to establishing the bbox
    ## and reset resulting bbox as NA
    ans <- suppressWarnings(combine_rows(lst))
    attr(ans[[attr(ans, "sf_column")]], "bbox") <-
      ## The bbox does not need to have a crs. This causes spurious warnings.
      structure(NA_bbox_, crs = NULL)
  } else {

    ## Check compatibility of CRSs
    crs_list <- lapply(lst, st_crs)
    if (length(crs <- unique(crs_list[!non_spatial])) > 1) {
      stop("Arguments have different coordinate reference systems. ",
           "Please use st_transform() to convert to a common CRS.")
    }

    ## Set non-spatial CRSs to the common CRS
    for (i in which(non_spatial)) {
      st_crs(lst[[i]]) <- crs[[1]]
    }

    ## row-bind recoded tables of release events
    ## and reset event id
    ans <- combine_rows(lst)
  }

  ## Reset class order
  ## Assignments to 'sf' objects result in class 'sf' becoming first in the list
  ## of classes.
  ans <- set_class(ans, 'sit_traps')

  ## Attach attribute table of release sites, possibly filtering.
  attr(ans, 'trap_types') <- trap_types

  return(ans)
}


#' @describeIn sit_traps Extracts trap data
#' @importFrom rlang .data
#' @export
sit_traps.sit <- function(
  x,
  type  = sit_trap_types(x)$label,
  stage = unique(sit_trap_types(x)$stage),
  area  = unique(x$traps$area),
  ...
) {


  xf <- x %>%
    dm_filter(traps = (.data$area %in% .env$area) ) %>%
    dm_filter(
      trap_types = (.data$stage %in% .env$stage & .data$label %in% .env$type)
    )

  ans <- structure(xf[['traps']], trap_types = xf[['trap_types']])

  ## Reset class order
  ## Assignments to 'sf' objects result in class 'sf' becoming first in the list
  ## of classes.
  ans <- set_class(ans, 'sit_traps')

  return(ans)
}

#' @export
print.sit_traps <- function(
  x,
  ...
) {

  ans <- as.data.frame.sit_traps(x)
  print(ans)
  invisible(ans)
  # ans
}


#' @export
as.data.frame.sit_traps <- function(x, ...) {

  ## trap types renamed into a separate object to prevent issues with
  ## lazy evaluation of the other arguments.
  tt_rn <- structure(
    attr(x, 'trap_types'),
    names = paste("type", names(attr(x, 'trap_types')), sep = "_")
  )

  ## Unclass
  class(x) <- setdiff(class(x), 'sit_traps')

  ans <- merge(
    x, tt_rn, by = "type_id", sort = FALSE
  )[, -1]
  ans <- ans[order(ans$id),]

  return(ans)
}



#' @importFrom stats aggregate
#' @export
summary.sit_traps <- function(object, ...) {

  traps <- as.data.frame.sit_traps(object)

  traps$has_coordinates <- !st_is_empty(st_geometry(traps))

  ans <- structure(
    stats::aggregate(
      # st_drop_geometry(traps),
      data.frame(n = seq.int(nrow(traps))),
      by = list(
        Area = traps$area,
        Type = traps$type_label,
        Coords = traps$has_coordinates
      ),
      length
    ),
    class = c("summary.sit_traps", "data.frame")
  )

  return(ans)
}

#' @export
print.summary.sit_traps <- function(x, ...) {

  cat("Number of traps in each group:\n")
  ans <- as.data.frame(unclass(x))
  print(ans)
  invisible(ans)
}



#' Estimate competitiveness of sterile males
#' @family estimates
#'
#' Computes the Fried Index, which estimates the competitiveness based on
#' observed fertility rates and the sterile-wild male ratio.
#'
#' It is strongly recommended to estimate competitiveness based on surveys from
#' areal releases. The calculation requires the estimation of the sterile-wild
#' male ratio from surveys pooled across traps and days, which requires a
#' roughly homogeneous distribution of sterile males.
#'
#' @param x A `sit` object.
#' @param following_releases A `sit_release` object with a subset of release
#'   events or missing (default) for all release events. Use in combination with
#'   `following_days` to filter surveys of the target populations within a given
#'   number of days after release. Note that counts of wild populations will
#'   always be included in the results with a value of `NA` in `pop_col`.
#' @param following_days Integer or missing (default). Number of days after
#'   releases to return, if `following_releases` is not missing.
#' @param residual_fertility Numeric. Default: 0.05.
#' @param species Character. Filter results for a given species or ignore the
#'   species if missing.
#'
#' @return A numeric value.
#' @export
#'
#' @examples
#'  sit_competitiveness(sit_prototype)
sit_competitiveness <- function(
  x,
  following_releases = sit_revents(x, type = "areal"),
  following_days = 7,
  residual_fertility = 0.01,
  species = NULL
) {

  ## Use only surveys within 7 days since the areal releases and consider
  ## only sterile males from those releases.
  adult_areal <- sit_adult_surveys(
    x,
    area = "sit",
    following_releases = following_releases,
    following_days = following_days,
    species = species
  )

  if (!nrow(adult_areal[adult_areal$pop_col == "wild", ]) > 0 |
      !nrow(adult_areal[adult_areal$pop_col != "wild", ]) > 0 ) {
    message(
      "Cannot compute sterile male competitiveness. ",
      "Surveys from adult sterile and wild males following an areal release missing."
    )
    return(NULL)
  }

  swmr <- sterile_wild_male_ratio(adult_areal, pool = TRUE)

  egg_control <- sit_egg_surveys(x, area = "control", species = species)

  if (!nrow(egg_control) > 0) {
    message(
      "Cannot compute sterile male competitiveness. ",
      "Surveys from ovitraps in the control area are missing."
    )
    return(NULL)
  }

  natural_fertility <- fertility_rate(
    egg_control,
    pool = TRUE
  )

  egg_areal <- sit_egg_surveys(
    x,
    area = "sit",
    following_releases = following_releases,
    following_days = following_days,
    species = species
  )

  if (!nrow(egg_areal) > 0) {
    message(
      "Cannot compute sterile male competitiveness. ",
      "Surveys from ovitraps following an areal release missing."
    )
    return(NULL)
  }

  sit_fertility <- fertility_rate(egg_areal, pool = TRUE)

  res <- fried_index(
    sterile_wild_mr    = swmr,
    natural_fertility  = natural_fertility,
    sit_fertility      = sit_fertility,
    residual_fertility = residual_fertility
  )

  comp_table <- data.frame(
    Component = c(
      "Sterile-wild male ratio", "Natural fertility",
      "Observed fertility in SIT area", "Sterile fertility (assumed)"
    ),
    Value = c(swmr, natural_fertility, sit_fertility, residual_fertility)
  )

  ans <- structure(res, components = comp_table, class = 'sit_competitiveness')
  return(ans)
}


#' @importFrom knitr kable
#' @export
print.sit_competitiveness <- function(x, digits = 2, ...) {
  cli::cli_h1("Sterile male competitiveness")
  cat("\n")
  cli::cli_alert_info("Estimated value: {round(x, digits)}")
  print(
    knitr::kable(
      attr(x, "components"),
      format = "simple",
      digits = digits
    )
  )
  cat("\n")
}

#' Fried index
#'
#' Index of the mating competitiveness of sterile males relative to wild males.
#'
#' The competitiveness \eqn{\gamma} of the sterile male individuals is defined
#' as their relative capacity to mate with a wild female, compared to a wild
#' male.
#'
#' Thus, in a homogeneously mixed population with \eqn{M_s} sterile males and
#' \eqn{M_w} wild males, the probability that a mating occurs with a sterile
#' individual is
#' \deqn{
#'   p_s = \frac{\gamma M_s}{M_w + \gamma M_s} =
#'       \frac{\gamma R_{sw}}{1 + \gamma R_{sw}},
#' }
#' where \eqn{R_{sw} = M_s/M_w} is the sterile-wild male ratio.
#' At a given sterile-wild male ratio \eqn{R_{sw} = M_s/M_w} we observe a
#' fertility rate \eqn{H_s} in the field.
#'
#' Assuming a residual fertility rate \eqn{H_{rs}} for sterile males and a
#' natural fertility rate \eqn{H_w} for wild males, the observed fertility rate
#' \eqn{H_s} in the field is:
#' \deqn{
#'    H_s = p_s H_{rs} + (1-p_s)H_w =
#'        \frac{\gamma R_{sw}}{1 + \gamma R_{sw}} H_{rs} +
#'            \frac{1}{1 + \gamma R_{sw}} H_w.
#'  }

#' @param sterile_wild_mr Non-negative number. Ratio of sterile to wild males
#'   in the target population.
#' @param natural_fertility Number between 0 and 1. Proportion of fertile eggs
#'   in a natural (wild) population.
#' @param sit_fertility Number between 0 and 1. Proportion of fertile eggs
#'   in the target population.
#' @param residual_fertility Number between 0 and 1. Proportion of fertile eggs
#'   in a completely sterile population. Also called _residual_ fertility.
#'
#' @return Non-negative number.
#' @export
#'
#' @examples
#' fried_index(
#'   sterile_wild_mr = 0.3,
#'   natural_fertility = .9,
#'   sit_fertility = .71,
#'   residual_fertility = .05
#' )
fried_index <- function(
  sterile_wild_mr,
  natural_fertility,
  sit_fertility,
  residual_fertility
) {

  if (!any(sterile_wild_mr > 0)) {
    stop(
      "Cannot estimate competitiveness with a estimated wild population of 0."
    )
  }

  if (any(sit_fertility == residual_fertility, na.rm = TRUE)) {
    stop(
      "Cannot estimate competitiveness when the estimated fertility equals",
      " the sterile fertility."
    )
  }

  (natural_fertility - sit_fertility) /
    (sit_fertility - residual_fertility) /
    sterile_wild_mr

}

#' Fertility rate
#'
#' Observed fertility rate from egg surveys.
#'
#' The `sit` method uses all egg surveys in the study. For a `sit` object `x`,
#' `fertility_rate(x)` is thus equivalent to
#' `fertility_rate(sit_egg_surveys(x))`.
#'
#' @param x Either a `sit` or a `sit_egg_surveys` object.
#' @param species Character. Filter results for a given species or ignore the
#'   species if missing.
#' @param pool Logical. If `TRUE` pools fertile and sterile counts from all traps
#'   and survey dates. Otherwise (default) yields results by trap and survey
#'   date.
#'
#' @return If `pool = TRUE`, a number. Otherwise, a data.frame with counts of
#'   fertile and sterile eggs by trap and survey period, with the corresponding
#'   proportion of fertile eggs.
#' @export
#' @examples
#' fertility_rate(sit_prototype, pool = TRUE)
fertility_rate <- function(x, species = NULL, pool = FALSE) {
  UseMethod('fertility_rate')
}

#' @export
fertility_rate.sit <- function(x, species, pool) {

  as_sit <- sit_egg_surveys(x)
  fertility_rate.sit_egg_surveys(as_sit, species, pool)

}

#' @importFrom stats reshape
#' @export
fertility_rate.sit_egg_surveys <- function(x, species, pool) {

  ## Aggregate counts by trap, survey date and fertility status
  x_long <- as.data.frame(x[, -1])

  x_wide <- stats::reshape(
    x_long,
    direction = "wide",
    idvar = c("trap_code", "datetime_start", "datetime_end"),
    timevar = "fertile",
    v.names = "n",
    varying = list(c("fertile", "sterile"))
  )
  attr(x_wide, "reshapeWide") <- NULL

  if (pool) {
    tot_fertile <- sum(x_wide$fertile)
    total <- tot_fertile + sum(x_wide$sterile)
    return(tot_fertile / total)
  }

  x_wide$fertility_rate = x_wide$fertile / (x_wide$fertile + x_wide$sterile)

  return(x_wide)

}


#' Sterile-Wild male ratio
#'
#' Observed ratio between sterile and wild male individuals from adult surveys.
#'
#' The `sit` method uses all adult surveys from the study (`sit`) area. For a
#' `sit` object `x`, `sterile_wild_male_ratio(x)` is thus equivalent to
#' `sterile_wild_male_ratio(sit_adult_surveys(x, area = 'sit'))`.
#'
#' @param x Either a `sit` or a `sit_adult_surveys` object.
#' @param pool Logical. If `TRUE` pools sterile and wild counts from all traps
#'   and survey dates. Otherwise (default) yields results by trap and survey
#'   date.
#'
#' @return If `pool = TRUE`, a number. Otherwise, a data.frame with counts of
#'   sterile and wild males by trap and survey period, with the corresponding
#'   ratio.
#' @export
#' @examples
#' sterile_wild_male_ratio(sit_prototype, pool = TRUE)
sterile_wild_male_ratio <- function(x, pool) {
  UseMethod('sterile_wild_male_ratio')
}

#' @export
sterile_wild_male_ratio.sit <- function(x, pool = FALSE) {

  as_sit <- sit_adult_surveys(x, area = "sit")
  sterile_wild_male_ratio.sit_adult_surveys(as_sit, pool)

}

#' @importFrom stats reshape
#' @export
sterile_wild_male_ratio.sit_adult_surveys <- function(x, pool = FALSE) {

  x_males <- x[x$sex == "male", ]

  ## Aggregate counts from all sterile populations together
  x_males$group <- ifelse(x_males$pop_col == 'wild', 'wild', 'sterile')
  x_group <- aggregate(
    x_males['n'],
    by = list(
      trap_code = x_males$trap_code,
      datetime_start = x_males$datetime_start,
      datetime_end = x_males$datetime_end,
      group = x_males$group
    ),
    FUN = "sum",
    na.rm = TRUE
  )

  x_wide <- stats::reshape(
    x_group,
    direction = "wide",
    idvar = c("trap_code", "datetime_start", "datetime_end"),
    timevar = "group"
  )
  x_wide <- setNames(x_wide, gsub("^n\\.", "", names(x_wide)))
  attr(x_wide, "reshapeWide") <- NULL

  ## Replace NAs due to missing records of some group at some trap and date by 0
  x_wide$sterile[is.na(x_wide$sterile)] <- 0
  x_wide$wild[is.na(x_wide$wild)] <- 0

  if (pool) {
    return(sum(x_wide$sterile) / sum(x_wide$wild))
  }

  x_wide$swmr = x_wide$sterile / x_wide$wild

  return(x_wide)

}

#' Build a `sit` object
#'
#' Gathers data about traps, release events and field survey data into
#' a relational data model.
#'
#' @param traps A `sit_traps` object built with `sit_traps()`.
#' @param release_events A `sit_revents` object built with `sit_revents()`.
#' @param adult_surveys A `sit_adult_surveys` object built with `sit_adult_surveys()`.
#' @param egg_surveys A `sit_egg_surveys` object built with `sit_egg_surveys()`.
#'
#' @return A `sit` object that can be queried with functions from the package.
# @seealso [sit_sterile_male_competitivenes()], [sit_mdt()],
#   [sit_flight_range()], [sit_diffusion()], [sit_survival()].
#' @import dm
#' @export
#'
#' @examples
#' sit(
#'   sit_traps(fake_traps),
#'   c(sit_revents(fake_rpoints), sit_revents(fake_rareal)),
#'   sit_adult_surveys(fake_adults),
#'   sit_egg_surveys(fake_eggs)
#' )
sit <- function(traps, release_events, adult_surveys = NULL, egg_surveys = NULL) {

  # # Test example from Readme.
  # traps = sit_traps(fake_traps)
  # release_events = c(sit_revents(fake_rpoints), sit_revents(fake_rareal))
  # adult_surveys = sit_adult_surveys(fake_adults)
  # egg_surveys = sit_egg_surveys(fake_eggs)

  ## Extract embedded tables from attributes
  trap_types <- attr(traps, "trap_types")
  release_sites <- attr(release_events, "release_sites")
  attr(traps, "trap_types") <- NULL
  attr(release_events, "release_sites") <- NULL

  ## Remove the corresponding classes
  class(traps) <- setdiff(class(traps), 'sit_traps')
  class(release_events) <- setdiff(class(release_events), 'sit_revents')


  if (is.null(adult_surveys)) {
    ## If not specified, use a empty table with no observations
    adult_surveys <- sit_adult_surveys(x = sit::fake_adults[c(),])
  }

  if (is.null(egg_surveys)) {
    ## If not specified, use a empty table with no observations
    egg_surveys <- sit_egg_surveys(x = sit::fake_eggs[c(),])
  }

  ## Check whether variable 'trap_type' is required in adult data to uniquely
  ## identify traps
  non_unique_trap_codes <- any(
    duplicated(unique(st_drop_geometry(traps[c("code", "type_id")]))$code)
  )
  if (non_unique_trap_codes) {
    multiple_adult_traps <- sum(trap_types$stage == 'adult') > 1
    if (multiple_adult_traps) {
      adults_has_type <- 'trap_type' %in% names(adult_surveys)
      if (!adults_has_type) {
        stop(
          glue(
            "Variable `trap_type` is required in adult data, ",
            "in order to uniquely identify the correct trap."
          )
        )
      }
    } else {
      # Only one adult trap. Use it as trap type with message.
      adult_surveys$type <- trap_types$label[trap_types$stage == "adult"]
    }
    # If traps are identified by code and type, need to match both
    # variables using a merge.
    trap_labels <- merge(
      st_drop_geometry(traps)[, c("id", "code", "type_id")],
      trap_types[, c("id", "label")],
      by.x = "type_id",
      by.y = "id"
    )

    a_trap_code.idx <- merge(
      setNames(
        adult_surveys[, c("id", "trap_code", "trap_type")],
        c("survey_id", "code", "label")
      ),
      trap_labels,
      by = c("code", "label"),
      sort = FALSE
    )
    a_trap_code.idx <- a_trap_code.idx$id[order(a_trap_code.idx$survey_id)]
  } else {
    # Unique trap codes: simply match codes
    a_trap_code.idx <- match(adult_surveys$trap_code, traps$code)
  }

  ## Check trap codes in adult data
  if (any(idx <- is.na(a_trap_code.idx))) {
    stop(
      glue(
        "Trap {vars_collapse(unique(adult_surveys$trap_code[idx]), 'code ', 'codes ')} ",
        "in adult field data not found in the imported traps."
      )
    )
  }

  ## Check trap types in adult data
  stages <- trap_types$stage[traps$type_id[a_trap_code.idx]]
  if (any(idx <- stages != "adult")) {
    stop(
      glue(
        "Trap {vars_collapse(unique(adult_surveys$trap_code[idx]), 'code ', 'codes ')} ",
        "not of adult type. Please verify."
      )
    )
  }

  ## Check population colours in adult data
  ## a_cols.idx matches colours and NAs for wild populations
  ## cols.idx is the reduced set of non-wild colours
  a_cols.idx <- match(adult_surveys$pop_col, release_events$colour)
  if (any(idx <- is.na(a_cols.idx) & adult_surveys$pop_col != "wild")) {
    stop(
      glue(
        "Release with ",
        "{vars_collapse(unique(adult_surveys$pop_col[idx]), 'colour ', 'colours ')} ",
        "not found. Please verify."
      )
    )
  }

  ## Build foreign keys
  adult_surveys$trap_id <- traps$id[a_trap_code.idx]
  adult_surveys$population_id <- release_events$id[a_cols.idx]


  ## Select and reorder columns
  avars.idx <- c(
    grep("id", names(adult_surveys)),
    grep("_code|_col|id", names(adult_surveys), invert = TRUE)
  )
  adult_surveys <- adult_surveys[, avars.idx]


  ## Check trap codes in egg data
  e_trap_code.idx <- match(egg_surveys$trap_code, traps$code)
  if (any(idx <- is.na(e_trap_code.idx))) {
    stop(
      glue(
        "Trap {vars_collapse(egg_surveys$trap_code[idx], 'code ', 'codes ')} ",
        "in egg field data not found in the imported traps."
      )
    )
  }

  ## Check trap types in egg data
  stages <- trap_types$stage[traps$type_id[e_trap_code.idx]]
  if (any(idx <- stages != "egg")) {
    stop(
      glue(
        "Trap {vars_collapse(egg_surveys$trap_code[idx], 'code ', 'codes ')} ",
        "not of egg type. Please verify."
      )
    )
  }

  ## Build foreign keys
  egg_surveys$trap_id <- traps$id[e_trap_code.idx]

  ## Select and reorder columns
  evars.idx <- c(
    grep("id", names(egg_surveys)),
    grep("_code|id", names(egg_surveys), invert = TRUE)
  )
  egg_surveys <- egg_surveys[, evars.idx]



  ## Build a relational data model with dm
  ## https://cynkra.github.io/dm/articles/howto-dm-df.html
  ans <- dm(
    adult_surveys,
    egg_surveys,
    release_events,
    release_sites,
    traps,
    trap_types
  )

  ## Specify primary keys and check uniqueness
  ans <- dm_add_pk(ans, adult_surveys, "id", check = TRUE)
  ans <- dm_add_pk(ans, egg_surveys, "id", check = TRUE)
  ans <- dm_add_pk(ans, release_events, "id", check = TRUE)
  ans <- dm_add_pk(ans, release_sites, "id", check = TRUE)
  ans <- dm_add_pk(ans, traps, "id", check = TRUE)
  ans <- dm_add_pk(ans, trap_types, "id", check = TRUE)

  ## Specify foreign keys and check correspondence
  ans <- dm_add_fk(ans, egg_surveys, "trap_id", traps, check = FALSE)
  ans <- dm_add_fk(ans, adult_surveys, "trap_id", traps, check = FALSE)
  ans <- dm_add_fk(
    ans,
    adult_surveys,
    "population_id",
    release_events,
    ## Cannot check since there can be NA values corresponding to wild
    ## populations, which are not present in table release_events.
    check = FALSE
  )
  ans <- dm_add_fk(ans, release_events, "site_id", release_sites, check = FALSE)

  ## Avoid checking for consistency. Problematic since traps is `sf`
  ## https://github.com/cynkra/dm/issues/705
  ans <- dm_add_fk(ans, traps, "type_id", trap_types, check = FALSE)


  ans <- dm_set_colors(
    ans,
    peru           = traps,
    darkmagenta    = release_events,
    darkolivegreen = c("egg_surveys", "adult_surveys")
  )
  class(ans) <- c("sit", class(ans))

  # dm_draw(ans, view_type = "all")

  ans <- validate_sit(ans)

  return(ans)
}


#' @keywords internal
#' @export
as_dm.sit <- function(...) {
  ## Remove sit class
  x <- list(...)[[1]]
  class(x) <- setdiff(class(x), "sit")
  return(x)
}

#' @import cli
#' @export
print.sit <- function(x, ...) {

  cli::cat_rule("sit", col = "violet")

  revents <- as.data.frame(sit_revents(x))
  revents_by_site <- table(revents$site_id)
  revents_by_type <- table(revents$type)

  revents_site_chr <- paste0(
    sum(revents_by_site),
    " (", paste(revents_by_site, collapse = " + "), ")"
    )

  # cli::cat_line("Release events: ", revents_site_chr)
  cli::cat_line("Release events: ", describe_numbers(revents_by_type))
  cli::cat_line("Release sites : ", length(revents_by_site))

  traps <- st_drop_geometry(as.data.frame(sit_traps(x)))
  traps_type <- table(traps$type_label)
  traps_stage_area <- table(traps$type_stage, traps$area)

  cli::cat_line("Adult traps   : ", describe_numbers(traps_stage_area["adult", ]))
  cli::cat_line("Egg traps     : ", describe_numbers(traps_stage_area["egg", ]))
  cli::cat_line("Trap types    : ", describe_numbers(traps_type))

  surveys <- c(adult = nrow(sit_adult_surveys(x)), egg = nrow(sit_egg_surveys(x)))

  cli::cat_line("Surveys       : ", describe_numbers(surveys))
}


#' @import cli
#' @importFrom stats setNames
#' @export
summary.sit <- function(object, ...) {

  cli::cat_rule("Releases", col = "violet")
  print(summary(sit_revents(object)))
  cli::cat_line()

  cli::cat_rule("Traps", col = "orange")
  print(summary(sit_traps(object)))
  cli::cat_line()

  cli::cat_rule("Adult surveys (total captures)", col = "green")

  as <- sit_adult_surveys(object)
  as$pop_col[is.na(as$pop_col)] <- 'wild'
  as$pop_col <- factor(as$pop_col, c(sit_revents(object)$colour, 'wild'))
  ## Add trap area
  as <- merge(
    as,
    setNames(
      st_drop_geometry(
        as.data.frame(sit_traps(object))
      )[, c('code', 'area')],
      c('trap_code', 'area')
    )
      ,
    by = 'trap_code'
  )
  print(
    tapply(
      seq_len(nrow(as)),
      list(as$pop_col, as$sex, as$area),
      function(.) sum(as$n[.])
    )
  )
  cli::cat_line()

  cli::cat_rule("Egg surveys (total captures)", col = "green")
  es <- sit_egg_surveys(object)
  es$fertile <- factor(es$fertile, labels = c('Sterile', 'Fertile'))
  ## Add trap area
  es <- merge(
    es,
    setNames(
      st_drop_geometry(
        as.data.frame(sit_traps(object))
      )[, c('code', 'area')],
      c('trap_code', 'area')
    )
    ,
    by = 'trap_code'
  )
  print(
    tapply(
      seq_len(nrow(es)),
      list(es$area, es$fertile),
      function(.) sum(es[., "n"])
    )
  )
  cli::cat_line()

  invisible(object)
}


#' @export
plot.sit <- function(x, y, ...) {

  ## Only plot areas and trap types for which there are coordinates
  summary_traps <- summary(sit_traps(x))
  target_traps <- summary_traps[summary_traps$Coords, ]

  ## Target only the sit area
  target_traps <- target_traps[target_traps$Area == "sit"]

  ## Plot release points and sit-area traps by type
  traps <- sit_traps(x, area = "sit", type = target_traps$Type)
  ## Merge with trap_types attribute to gain access to type_label
  ## Despite the name, returns a sf object
  traps <- as.data.frame(traps)

  rpoints <- st_as_sf(sit_revents(x, type = "point"))

  plot(traps["type_label"], pch = 20, col = "darkgray", reset = FALSE,
       main = "Traps and release points in the sit area")
  plot(rpoints["colour"], add = TRUE, pch = 19, col = rpoints$colour)

}

#' Collapse variables
#'
#' Collapse a vector of variables names into a properly formatted string that
#' can be used in messages.
#'
#' @param x Vector coerced to character. Variable names to be concatenated.
#' @param intro1 Character. String included at the beginning if `x` is of
#'   length 1.
#' @param intros Character. String included at the beginning if `x` is
#'   of length greater than 1.
#' @param quote Character. Symbol used to quote variable names.
#' @param sep Character. Symbol used to separate variable names.
#' @param last Character. String used to separate the last two items if `x` has
#'   at least 2 items.
#'
#' If `x` is an empty character vector, returns another empty character vector.
#'
#' @return
#' Character string.
#' @import glue
#' @keywords internal
#'
#' @examples
#' sit:::vars_collapse(paste0("var", 1:4))
#' sit:::vars_collapse("single_var")
#' sit:::vars_collapse("single_var", intro1 = "", quote = "")
vars_collapse <- function(
  x,
  intro1 = "Variable ",
  intros = "Variables ",
  quote = "'",
  sep = ", ",
  last = " and "
) {

  if (length(x) == 0L) return(x)

  vars_str <- glue_collapse(
    paste0(quote, x, quote), sep = sep, last = last
  )
  ans <- glue(
    "{variables}{vars_str}",
    variables = ifelse(length(x) > 1, intros, intro1)
  )
  return(unclass(ans))
}

#' Are given variable names present in a data.frame
#'
#' Verify that a vector of variable names are all actual variables in
#' a `data.frame`.
#'
#' @param x `data.frame`-like object
#' @param names_given Character vector of variable names.
#'
#' @return NULL if all variables exist in the `data.frame`. Otherwise, stops
#'   with a informative error message.
#' @import glue
#' @keywords internal
#'
#' @examples
#' sit:::check_variable_names(data.frame(a = 1), "a")
check_variable_names <- function(x, names_given) {
  if (! all(matching <- names_given %in% names(x))) {
    vars <- vars_collapse(names_given[!matching])
    stop(glue("{vars} not found in import data. Please verify."))
  }
}


#' Build the internal representation of a data set
#'
#' This function selects the relevant variables from the data set
#' and renames them according to a given set of internal names.
#'
#' @param x `data.frame`.
#' @param varnames Character vector. Variable names from `x`.
#' @param internal_names Character vector. Internal names.
#'
#' @return A `data.frame` with the internal variable names in the expected
#'   order.
#' @keywords internal
#'
#' @examples
#' xx <- data.frame(a = 1, b = 2, c = 3)
#' sit:::build_internal_df(xx, c("a", "c"), c("A", "C"))
build_internal_df <- function(x, varnames, internal_names) {

  ## Check that `varnames` exist in `x`
  check_variable_names(x, varnames)

  ## Remove unused columns
  ans <- x[, varnames]

  ## Rename variables using internal names
  ## The subsetting in match is necessary in the case of a `sf` object
  ## which has an additional column `geometry`.
  names(ans)[match(varnames, names(ans))] <- internal_names

  return(ans)
}


#' Test if a number is whole
#'
#' @param x numeric.
#' @param tol numeric. Tolerance.
#'
#' @return Logical.
#' @keywords internal
#'
#' @examples
#' sit:::is.wholenumber(5)
#' sit:::is.wholenumber(5.00002)
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) {
  abs(x - round(x)) < tol
}


#' Combine rows
#'
#' Remove duplicate rows (not taking into account id) and reset the id.
#'
#' @param x a list of `data.frame`s. All must share the same variables.
#' @param id_col Numeric. Column number of the identification variable.
#'
#' @return An object of the same class of the first element in `x`.
#' @keywords internal
#'
#' @examples
#' dat1 <- data.frame(id = 1:3, x = letters[1:3])
#' dat2 <- data.frame(id = 1:2, x = letters[2:3])
#' sit:::combine_rows(list(dat1, dat2))
combine_rows <- function(x, id_col = 1, recoder = FALSE) {

  ## Bind rows
  xbind <- do.call('rbind', x)

  ## Remove duplicates, disregarding the id
  ans <- unique(xbind[-id_col])

  ## Reset id
  ans$id <- seq.int(nrow(ans))

  ## Reorder columns
  ans <- ans[, c(ncol(ans), seq.int(ncol(ans) - 1))]

  if (recoder) {
    ## Recoding map
    rec_l <- lapply(
      x, function(.)
      merge(., ans, by = intersect(names(.)[-id_col], names(ans)[-1]))
    )

    ## Recoding function: recode a vector x from item i of the initial list.
    rec_f <- function(x, i) {
      rec_l[[i]]$id.y[match(x, rec_l[[i]]$id.x)]
    }

    attr(ans, 'recoder') <- rec_f
  }

  return(ans)
}


#' Compute survey periods
#'
#' Given the survey date and either the activation or duration, compute the
#' starting and ending survey dates and times.
#'
#' The function checks the precedence of the resulting starting and ending
#' times.
#'
#' @param x A `data.frame` object.
#' @param input A Character string. Either 'activation'  or 'duration'.
#'
#' @return The same data.frame with the input variables replaced by
#'   datetime_start and datetime_end.
#' @keywords internal
#'
#' @examples
#' td <- data.frame(
#'   id = 1:2, survey = c('2021-11-15', '2021-11-01'),
#'   activation = c('2021-11-10', '2021-10-30'), duration = c(7, 5)
#' )
#'
#' ## Note the difference between the two inputs due to daylight saving time
#' ## (which depends on the time zone).
#' sit:::compute_periods(td, input = "activation")
#' sit:::compute_periods(td, input = "duration")
compute_periods <- function(x, input) {

  if (input == 'duration' && !all(idx <- x$duration > 0)) {
    stop("Argument `duration` must be non-negative.")
  }

  x$datetime_end <- as.POSIXct(x$survey)
  if (input == 'activation') {
    x$datetime_start <- as.POSIXct(x$activation)
  } else {
    stopifnot(input == 'duration')
    x$datetime_start <- x$datetime_end -
      as.difftime(x$duration, units = "days")
  }

  if (!all(x$datetime_end - x$datetime_start > 0)) {
    stop("Dates of `activation` must precede survey dates.")
  }

  x$survey <- x$duration <- x$activation <- NULL

  return(x)
}


#' Colours for concepts in the data model.
#'
#' @param x Character vector. Concepts among `traps`, `release_events` or
#'   `surveys`. If missing, all three are returned in that order.
#'
#' @return A named vector of colours.
#' @keywords internal
#'
#' @examples
#' sit:::sit_colours()
#' sit:::sit_colours('release_events')
#' sit:::sit_colours(c('surveys', 'traps', 'surveys'))
sit_colours <- function(x) {
  cols <- c(
    traps = "peru", release_events = "darkmagenta", surveys = "darkolivegreen"
  )

  if (missing(x)) x <- names(cols)

  return(cols[x])
}

#' Assign a S3 class
#'
#' Append a class to the beginning of an object's list of classes.
#'
#' If the object already inherited from the class, it is moved to the beginning
#' of the class list. Otherwise, the new class is added.
#'
#' @param x Any object.
#' @param classname Character string. The class name to be assigned.
#'
#' @keywords internal
#' @return The same object `x` with the updated `class` attribute.
#'
#' @examples
#' sit:::set_class(7, 'prime')
#' sit:::set_class( sit:::set_class(7, 'prime'), 'numeric')
set_class <- function(x, classname) {

  class(x) <- c(classname, setdiff(class(x), classname))
  return (x)
}


#' Print a comma-separated vector of terms
#'
#' Borrowed from [dm].
#'
#' @param x Character vector.
#' @param max_commas Integer. Maximum number of terms to print.
#' @param capped Logical. If TRUE do not show the total number of elements for
#'   long lists
#' @param fun One-sided formula. Function to apply to the vector. Must return
#'   another vector.
#'
#' @keywords internal
#' @import cli
#' @import rlang
#' @examples
#' sit:::commas(letters[1:6])
#' sit:::commas(letters[1:9])
#' sit:::commas(letters[1:9], capped = TRUE)
#' sit:::commas(letters[1:9], fun = toupper)
commas <- function(x, max_commas = 6L, capped = FALSE, fun = identity) {

  fun <- rlang::as_function(fun)

  if (rlang::is_empty(x)) {
    x <- ""
  } else if (length(x) > max_commas) {
    cap <- if (!capped) paste0(" (", length(x), " total)")

    x <- c(
      fun(x[seq_len(max_commas - 1)]),
      paste0(cli::symbol$ellipsis, cap)
    )
  } else {
    x <- fun(x)
  }

  glue_collapse(x, sep = ", ")
}


#' Description of numeric vectors
#'
#' Textual description of numeric vectors as a sum of components for descriptive
#' purposes.
#'
#' If the vector is named, use the names in parentheses besides the
#' corresponding number.
#'
#' @param x Numeric vector. Optionally named.
#' @param collapse Character string. Symbol separating numbers.
#' @param zero.print Logical. If TRUE, omit zero-valued entries.
#'
#' @keywords internal
#' @examples
#' sit:::describe_numbers(1:4)
#' sit:::describe_numbers(1:4, collapse = ", ")
#' sit:::describe_numbers(setNames(1:3, letters[1:3]))
describe_numbers <- function(x, collapse = " + ", zero.print = FALSE) {

  if (missing(x)) x <- numeric()

  if (!is.null(names(x))) {
    xn <- paste0("(", names(x), ")")
    ans <- paste(x, xn)
  } else {
    ans <- as.character(x)
  }

  if (!zero.print) {
    ans <- ans[x > 0]
  }

  ans <- paste(ans, collapse = collapse)

  return(ans)
}


#' Are times within a period
#'
#' Test whether some date-times fall within a set of time periods.
#'
#' For each date-time in `x` returns `TRUE` if it falls within __any__ of the
#' temporal intervals defined by the corresponding elements in `starts` and
#' `ends`.
#'
#' The interval is open by the left and closed by the right. I.e., it includes
#' the right boundary but not the left.
#'
#' If `starts` and `end` are empty vectors (e.g. `as.POSIXct(character())`),
#' return `TRUE` for each element in `x`.
#'
#' @param x Vector of POSIXct date-times.
#' @param starts Vector of POSIXct date-time starting times.
#' @param periods Vector of POSIXct date-time ending times.
#'
#' @return Logical vector of the same size as `x`.
#' @keywords internal
#'
#' @examples
#' test_times <- as.POSIXct(
#'   c('2021-05-04', '2021-06-28', '2012-10-05', '2021-10-05')
#' )
#' spring_2021 <- as.POSIXct(c('2021-03-21', '2021-06-20'))
#' fall_2021 <- as.POSIXct(c('2021-09-21', '2021-12-20'))
#' starts <- c(spring_2021[1], fall_2021[1])
#' ends <- c(spring_2021[2], fall_2021[2])
#' sit:::is_date_within(test_times, starts, ends)
#'
#' # Left-open and right-closed
#' day <- as.POSIXct('2021-05-04')
#' next_day <- as.POSIXct('2021-05-05')
#' sit:::is_date_within(day, day, next_day)      # FALSE
#' sit:::is_date_within(next_day, day, next_day) # TRUE
is_date_within <- function(x, starts, ends) {

  stopifnot(
    inherits(x, 'POSIXct'),
    inherits(starts, 'POSIXct'),
    inherits(ends, 'POSIXct'),
    (n_periods <- length(starts)) == length(ends)
  )

  if (length(x) == 0L) return(logical())
  if (length(starts) == 0L) return(rep(TRUE, length(x)))

  stopifnot(all(starts <= ends))

  ## Matrix of size n_periods x n_times
  conditions <-
    matrix(
      vapply(x, `>`, rep(TRUE, length(starts)), starts) &
        vapply(x, `<=`, rep(TRUE, length(ends)), ends),
      n_periods, length(x)
    )

  return(apply(conditions, 2, any))
}


#' Release types
#'
#' @param x `sf` of `sfc` object, representing the table of release sites.
#' @return Logical vector with as many elements as rows in `x`.
#' @keywords internal
#' @examples
#' sit:::release_types(sit_prototype[['release_sites']])
release_types <- function(x) {
  ifelse(st_is_empty(x), 'areal', 'point')
}


#' Get time periods following releases
#'
#' Retrieve the time periods from the release dates to a given number of
#' days after the corresponding release dates.
#'
#' If both arguments are missing, returns empty vectors.
#' Otherwise, specifying `following_releases` with missing `following_days`
#' results in an error, whereas the converse is ignored with a warning.
#'
#' @param following_releases A `sit_revents` object.
#' @param following_days Numeric. Number of days targeted. Can be fractional.
#'
#' @return A list with vectors `starts` and `ends` of `POSIXct` dates.
#' @keywords internal
#'
#' @examples
#' sit:::get_periods(sit_revents(sit_prototype), 7.5)
get_periods <- function(following_releases, following_days) {
  if (!missing(following_releases)) {
    if (missing(following_days)) {
      stop("Use of argument following_releases requires ",
           "following_days, which is currently missing.")
    }
    stopifnot(inherits(following_releases, 'sit_revents'))
    stopifnot(is.numeric(following_days), all(following_days > 0))

    starts <- as.POSIXct(following_releases$date)
    ends <- starts + as.difftime(following_days, units = "days")
  } else {
    ## missing following_releases
    if (!missing(following_days)) {
      warning("Argument `following_days` ignored.")
    }
    ## Using empty periods in is_date_within returns TRUE
    starts <- ends <- as.POSIXct(character())
  }

  return(list(starts = starts, ends = ends))
}


#' Maximum observed age
#'
#' Maximum age observed in a set of captures.
#'
#' @param captures Number of observed individuals.
#' @param ages Ages of the corresponding captures.
#'
#' @return Number.
#' @keywords internal
#'
#' @examples
#' sit:::max_observed_age(captures = rep(1, 5), ages = c(1, 2, 8, 3, 0))
#' sit:::max_observed_age(captures = c(1, 1, 0, 1, 1), ages = c(1, 2, 8, 3, 0))
max_observed_age <- function(captures, ages) {

  stopifnot(
    length(captures) == length(ages),
    is.numeric(captures),
    is.numeric(ages)
  )

  if (!length(captures) > 0) return(NA)

  counts_by_age <-
    tapply(captures, ages, FUN = sum, na.rm = TRUE)

  max_age.idx <- max(which(counts_by_age > 0))
  ans <- as.numeric(names(counts_by_age)[max_age.idx])

  return(ans)
}



#' Print a table of variables and their description
#'
#' @param ... A even number of character strings. Variable names followed by
#'   their description. Variables are printed in mono-spaced font.
#' @param header Logical. Include the header "Variable" "Description".
#' @param caption Character string. Table caption.
#'
#' @return A kable object.
#' @keywords internal
#' @export
#'
#' @examples
#' print_variable_table("Var1", "Desc1", "Var2", "Desc2")
#' print_variable_table(
#'   "Var1", "Desc1",
#'   "Var2", "Desc2",
#'   header = FALSE
#' )
print_variable_table <- function(
  ..., header = TRUE, caption = NULL) {

  ans <- as.data.frame(matrix(c(...), byrow = TRUE, ncol = 2))

  ## Use monospace font for variables
  ans[, 1] <- paste0("`", ans[, 1], "`")

  ## Set header
  ans <- setNames(
    ans,
    switch(header, c("Variable", "Description"), NULL)
  )

  kableExtra::kbl(
    ans,
    booktabs = TRUE,
    toprule = '',
    midrule = '',
    bottomrule = '',
    vline = '',
    escape = FALSE,
    centering = FALSE
  ) %>%
    kableExtra::kable_styling(full_width = FALSE, position = "left") %>%
    kableExtra::column_spec(2, width = "4in")

}

#' Import or retrieve adult surveys
#'
#' Import survey field data from adult traps, or retrieve the data from a `sit`
#' object.
#'
#' The argument `type` is mandatory only if needed to uniquely identify the trap
#' (i.e., more than one adult trap type is declared with [sit_trap_types()] and
#' the same code has been used for traps of different types). If the user used
#' unique codes for each trap, which will be checked internally, this is
#' redundant and will be ignored in any case.
#'
#' __Dates__ must be provided in the [RFC 3339
#' format](https://www.ietf.org/rfc/rfc3339.txt) (a variation of the [ISO 8601
#' format](https://en.wikipedia.org/wiki/ISO_8601)), i.e., `2021-12-31` or
#' `2019-11-23 15:00`.
#'
#' @param x A `data.frame` with field survey data in
#'   [tidy](https://www.openscapes.org/blog/2020/10/12/tidy-data/) format.
#' @param trap Identification code of the trap. This is the user’s code
#'   corresponding to the `id` argument in [sit_traps()].
#' @param type Character string. Trap type.
#' @param survey Character. Date and optionally the time of the survey. See
#'   details on the date format.
#' @param activation Character. Date and optionally the time of the survey
#'   activation. See details on the date format. Either `activation` or
#'   `duration` are required.
#' @param duration Numeric. Number of days (integer or fractional) that the trap
#'   has been functioning.
#' @param population Character. Either `wild` or any of the colour codes in the
#'   release events imported with [sit_revents()].
#' @param species Character. Only required if the counts of wild individuals are
#'   further categorised by species. Make sure to use the same nomenclature as
#'   in [sit_revents()]. In the extraction method, a character vector of species
#'   to be returned. Defaults to `NULL`, which ignores the species variable.
#' @param sex Character. Either `male` or `female`.
#' @param n Mandatory. An integer number of the surveyed individuals in the
#'   corresponding group.
#' @param ... Used to pass arguments to specific method.
#'
#' @return A object of class `sit_adult_surveys` which can be used in [sit()].
#' @importFrom utils head tail
#' @export
#'
#' @family importing
#' @seealso [sit()]
#'
#' @examples
#' ad_surv <- data.frame(
#'   trap = 1:3,
#'   survey = c("2021-04-01", "2021-08-21", "2021-08-21"),
#'   duration = rep(7L, 3),
#'   population = c("blue", "yellow", "yellow"),
#'   species = "aeg",
#'   sex = c("male", "male", "female"),
#'   n = 1:3
#' )
#' sit_adult_surveys(ad_surv)
sit_adult_surveys <- function(x, ...) {
  UseMethod('sit_adult_surveys')
}


#' @describeIn sit_adult_surveys Imports field survey data from adult traps.
#' @export
sit_adult_surveys.data.frame <- function(
  x,
  trap = "trap",
  type = "type",
  survey = "survey",
  activation = "activation",
  duration = "duration",
  population = "population",
  species = "species",
  sex = "sex",
  n = "n",
  ...
) {
  # # Debug
  # args <- names(formals(sit_adult_surveys.data.frame)[c(-1, -11)])
  # for (arg in args) assign(arg, arg)

  ## The requirement of the `type` argument can only be verified when building
  ## the `sit` object, by looking at the codes and types of the available traps.

  ## Check here whether there are variables in the data set corresponding
  ## to optional arguments and set them to NULL if not.
  opt_args <- c(type, activation, duration, species)
  for (arg in opt_args[!opt_args %in% names(x)]) {
    assign(arg, NULL)
  }

  ## Exactly one of `activation` or `duration` are required
  if ((l_ad <- length(c(activation, duration))) != 1L) {
    if (l_ad == 0L) {
      stop("Either activation or duration are required arguments. ",
           "Both are currently missing.")
    }
    if (l_ad == 2L) {
      stop("Only one of `activation` or `duration` needed. ",
           "Please remove one.")
    }
  }

  ## Build the internal representation
  ## The `population_id` can only be completed once merging with
  ## a `sit_revents` into a `sit` object.
  inames <- c(
    "trap_code", "trap_type"[!is.null(type)], "survey",
    "activation"[!is.null(activation)],
    "duration"[!is.null(duration)], "pop_col",
    "species"[!is.null(species)], "sex", "n"
  )
  # Variable names reported through non-NULL arguments
  # vnames <- unlist(lapply(tail(head(names(formals()), -1), -1), get,
  #                         pos = environment()))
  # Alternatively:
  vnames <- unname(
    unlist(
      as.list(
        environment())[tail(head(names(formals()), -1), -1)]
    )
  )

  ans <- build_internal_df(
    x[!is.na(x$n), ],  # remove surveys with missing counts
    varnames = vnames,
    internal_names = inames
  )

  ## Check duplicated surveys
  if (any(idx <- duplicated(ans))) {
    stop(
      glue(
        "{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
        "duplicated. Please fix."
      )
    )
  }

  ## Compute starting and ending dates and times
  ans <- compute_periods(
    ans,
    input = ifelse(is.null(activation), "duration", "activation")
  )

  ## Survey id
  ans$id <- if (nrow(ans) > 0) {
    seq.int(nrow(ans))
  } else {
    integer()
  }

  ## Reorder columns
  ans <- ans[
    , c("id",
        "trap_code",          # Changed into `trap_id` within `sit()`
        "trap_type"[!is.null(type)],  # Possibly needed for `trap_id`
        "pop_col",            # Changed into `population_id` within `sit()`
        "datetime_start",
        "datetime_end",
        "species"[!is.null(species)],
        "sex",
        "n"
        )
  ]

  class(ans) <- c("sit_adult_surveys", class(ans))
  return(ans)
}


#' @param area Character vector. Either `control`, `sit` or both (default).
#' @param trap_type Character vector. Any subset of
#'   `unique(sit_trap_types(x)$label)`.
#' @param following_releases A `sit_release` object with a subset of release
#'   events or missing (default) for all release events. Use in combination with
#'   `following_days` to filter surveys of the target populations within a given
#'   number of days after release. Note that counts of wild populations will
#'   always be included in the results with a value of `NA` in `pop_col`.
#' @param following_days Integer or missing (default). Number of days after
#'   releases to return, if `following_releases` is not missing.
#' @describeIn sit_adult_surveys Retrieve adult survey data, possibly filtering
#'   results.
#' @importFrom rlang .data .env
#' @importFrom tidyr everything
#' @export
sit_adult_surveys.sit <- function(
  x,
  area = c('control', 'sit'),
  trap_type = sit_trap_types(x)$label,
  following_releases,
  following_days,
  species = NULL,
  ...) {

  target_periods <- get_periods(following_releases, following_days)

  ## Filter trap areas, trap types, and dates
  x_sub <-
    x %>%
    dm_filter(traps = (.data$area %in% .env$area)) %>%
    dm_filter(trap_types = (.data$label %in% .env$trap_type)) %>%
    dm_filter(
      adult_surveys =
        is_date_within(
          .data$datetime_end,
          .env$target_periods$starts,
          .env$target_periods$ends
        )
    )

  # ## Filters for population if required
  # ## Can't use dm_filter, as it removes surveys of wild population
  # ## in adult_surveys. Filter the final object manually below.
  # ## https://github.com/cynkra/dm/issues/706
  # if (!missing(following_releases)) {
  #   x_sub <- dm_filter(
  #     x_sub,
  #     'release_events',
  #     .data$colour %in% following_releases$colour
  #   )
  #
  # }

  ## Filters for species if required
  if (!is.null(species)) {
    if (!'species' %in% names(x[['adult_surveys']])) {
      stop(
        "There is no variable `species` in the adult surveys."
      )
    }
    x_sub <- dm_filter(
      x_sub,
      adult_surveys = (.data$species %in% .env$species)
    )
  }

  ans <- x_sub %>%
    # dm_get_tables() %>% `[`(c("traps", "trap_types"))
    ## Select and rename columns from relevant tables
    dm_select('release_events', .data$id, pop_col = .data$colour) %>%
    dm_select('traps', .data$id, trap_code = .data$code) %>%
    dm_flatten_to_tbl('adult_surveys') %>%
    select(-.data$trap_id, -.data$population_id, -.data$geometry) %>%
    select(.data$id, .data$trap_code, .data$pop_col, everything()) %>%
    set_class('sit_adult_surveys')

  ## Return wild population as 'wild' rather than NA
  ans$pop_col[is.na(ans$pop_col)] <- 'wild'

  ## Filters for population if required. Always keep wild populations.
  if (!missing(following_releases)) {
    idx <- ans$pop_col %in% c(following_releases$colour, 'wild')
    ans <- ans[idx, ]
  }


  return(ans)
}


#' Import or extract trap types
#'
#' Extracts the table of trap types used in a `sit` experiment, or defines
#' a table of trap types.
#'
#' __Returns__ the table of trap types used in the `sit` object `x`.
#'
#' __Defines__ a table of trap types by passing a `data.frame` with variables
#' `name`, `label`, `stage` (either `egg` or `adult`) and optionally
#' `description`. Any other variable is ignored with a warning.
#'
#' __Returns__ the default table of trap types by calling `sit_trap_types()`
#' without any argument.
#'
#' Check the default trap types with `sit_trap_types()`.
#' You can store the table into an object and make edits as for a regular
#' `data.frame`. E.g. add another trap type, or edit labels or descriptions.
#' You can use the edited object as an input to the `trap_types` arguments in
#' [sit_traps()].
#'
#' @param x A `sit` object, a `data.frame`-like object or empty.
#'
#' @return Table of trap types used in the `sit` experiment `x`. If `x` is a
#'   `data.frame`, returns it as a table of trap types after some verifications.
#'   If `x` is missing, returns the current table of trap types.
#' @export
#'
#' @examples
#' sit_trap_types()  # default trap types.
#'
#' # Define new trap types
#' my_traps <- data.frame(name = "My New Trap", label = "MNT", stage = "adult")
#' sit_trap_types(my_traps)
sit_trap_types <- function(x) {

  if (missing(x)) {
    ## Return default table of trap types
    return(default_trap_types())
  }

  UseMethod("sit_trap_types")

}


#' @describeIn sit_trap_types Extract trap data
#' @export
sit_trap_types.sit <- function(x) {

  return(x[['trap_types']])
}

#' @describeIn sit_trap_types Import trap data
#' @export
sit_trap_types.data.frame <- function(x) {

  mandatory_vars <- c("name", "label", "stage")

  ## id is optional, but ignored
  ## description is optional and used
  optional_vars <- c("id", "description")

  ## Check mandatory variables
  if (any(miss <- (!mandatory_vars %in% names(x)))) {
    stop(
      glue(
        "Required variables ",
        "{vars_collapse(mandatory_vars[miss], intros = '')}",
        " not found in {substitute(x)}. Please verify."
      )
    )
  }

  ## Check valid values in stages
  xtages <- setdiff(unique(x$stage), c("egg", "adult"))
  if (length(xtages) > 0) {
    stop(
      glue(
        "Trap-type ",
        "{vars_collapse(xtages, intro1 = 'stage ', intros = 'stages ')} ",
        "not valid. Only 'egg' and 'adult' are valid values. Case sensitively."
      )
    )
  }

  ## Warning about ignoring id
  if ('id' %in% names(x)) {
    warning(
      glue(
        "Variable 'id' from {substitute(x)} ignored. ",
        "Trap-type ids are overriden internally."
      )
    )
  }

  ## Warning about ignoring additional variables
  extra_vars <- setdiff(names(x), c(mandatory_vars, optional_vars))
  if (length(extra_vars) > 0) {
    warning(
      glue(
        "{vars_collapse(extra_vars)} from {substitute(x)} ignored."
      )
    )
  }

  ## Make empty description if not provided
  if (!'description' %in% names(x)) {
    x$description <- NA
  }

  ans <- cbind(
    id = seq.int(nrow(x)),
    as.data.frame(x[, c(mandatory_vars, "description")])
  )
  class(ans) <- c("sit_trap_types", class(ans))

  return(ans)
}

#' @export
sit_trap_types.default <- function(x) {

  stop(glue(
    "{deparse(substitute(x))} must either be a 'sit' object to extract",
    " trap types from, a data.frame to define trap types",
    " or missing to get the currently defined trap types."
  ))
}

#' Import or retrieve egg surveys
#'
#' Import survey field data from egg traps, or retrieve the data from a `sit`
#' object.
#'
#' The argument `type` is mandatory only if needed to uniquely identify the trap
#' (i.e., more than one egg trap type is declared with [sit_trap_types()] and
#' the same code has been used for traps of different types). If the user used
#' unique codes for each trap, which will be checked internally, this is
#' redundant and will be ignored in any case.
#'
#' __Dates__ must be provided in the [RFC 3339
#' format](https://www.ietf.org/rfc/rfc3339.txt) (a variation of the [ISO 8601
#' format](https://en.wikipedia.org/wiki/ISO_8601)), i.e., `2021-12-31` or
#' `2019-11-23 15:00`.
#'
#' @param x A `data.frame` with field survey data in
#'   [tidy](https://www.openscapes.org/blog/2020/10/12/tidy-data/) format.
#' @param trap Identification code of the trap. This is the user’s code
#'   corresponding to the `id` argument in [sit_traps()].
#' @param type Character string. Trap type.
#' @param survey Character. Date and optionally the time of the survey. See
#'   details on the date format.
#' @param activation Character. Date and optionally the time of the survey
#'   activation. See details on the date format. Either `activation` or
#'   `duration` are required.
#' @param duration Numeric. Number of days (integer or fractional) that the trap
#'   has been functioning.
#' @param species Character.
#' @param fertile Logical. Whether `n` represents counts of __fertile__ eggs
#'   (`TRUE`) of __sterile__ eggs (`FALSE`).
#' @param n Mandatory. An integer number of the surveyed individuals in the
#'   corresponding group.
#' @param ... Used to pass arguments to specific method.
#'
#' @return A object of class `sit_egg_surveys` which can be used in [sit()].
#' @importFrom utils head tail
#' @export
#'
#' @family importing
#' @seealso [sit()]
#'
#' @examples
#' egg_surv <- data.frame(
#'   trap = 1:3,
#'   survey = c("2021-04-01", "2021-08-21", "2021-08-21"),
#'   duration = rep(7L, 3),
#'   fertile = c(TRUE, FALSE, TRUE),
#'   n = 1:3
#' )
#' sit_egg_surveys(egg_surv)
sit_egg_surveys <- function(x, ...) {
  UseMethod('sit_egg_surveys')
}


#' @describeIn sit_egg_surveys Imports field survey data from egg traps.
#' @export
sit_egg_surveys.data.frame <- function(
  x,
  trap = "trap",
  type = "type",
  survey = "survey",
  activation = "activation",
  duration = "duration",
  species = "species",
  fertile = "fertile",
  n = "n",
  ...
) {
  # # Debug
  # args <- names(formals(sit_egg_surveys.data.frame)[c(-1, -11)])
  # for (arg in args) assign(arg, arg)

  ## The requirement of the `type` argument can only be verified when building
  ## the `sit` object, by looking at the codes and types of the available traps.

  ## Check here whether there are variables in the data set corresponding
  ## to optional arguments and set them to NULL if not.
  opt_args <- c(type, activation, duration, species)
  for (arg in opt_args[!opt_args %in% names(x)]) {
    assign(arg, NULL)
  }

  ## Exactly one of `activation` or `duration` are required
  if ((l_ad <- length(c(activation, duration))) != 1L) {
    if (l_ad == 0L) {
      stop("Either activation or duration are required arguments. ",
           "Both are currently missing.")
    }
    if (l_ad == 2L) {
      stop("Only one of `activation` or `duration` needed. ",
           "Please remove one.")
    }
  }

  ## Build the internal representation
  inames <- c(
    "trap_code", "trap_type"[!is.null(type)], "survey",
    "activation"[!is.null(activation)],
    "duration"[!is.null(duration)],
    "species"[!is.null(species)], "fertile", "n"
  )
  # Variable names reported through non-NULL arguments
  # vnames <- unlist(lapply(tail(head(names(formals()), -1), -1), get,
  #                         pos = environment()))
  # Alternatively:
  vnames <- unname(
    unlist(
      as.list(
        environment())[tail(head(names(formals()), -1), -1)]
    )
  )

  ans <- build_internal_df(x, varnames = vnames, internal_names = inames)

  ## Check duplicated surveys
  if (any(idx <- duplicated(ans))) {
    stop(
      glue(
        "{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
        "duplicated. Please fix."
      )
    )
  }

  ## Compute starting and ending dates and times
  ans <- compute_periods(
    ans,
    input = ifelse(is.null(activation), "duration", "activation")
  )

  ## Survey id
  ans$id <- if (nrow(ans) > 0) {
    seq.int(nrow(ans))
  } else {
    integer()
  }

  ## Reorder columns
  ans <- ans[
    , c("id",
        "trap_code",          # Changed into `trap_id` within `sit()`
        "trap_type"[!is.null(type)],  # Possibly needed for `trap_id`
        "datetime_start",
        "datetime_end",
        "species"[!is.null(species)],
        "fertile",
        "n"
        )
  ]

  class(ans) <- c("sit_egg_surveys", class(ans))
  return(ans)
}


#' @param area Character vector. Either `control`, `sit` or both (default).
#' @param trap_type Character vector. Any subset of
#'   `unique(sit_trap_types(x)$label)`.
#' @param following_releases A `sit_release` object with a subset of release
#'   events or missing (default) for all release events. Use in combination with
#'   `following_days` to filter surveys of the target populations within a given
#'   number of days after release. Note that counts of wild populations will
#'   always be included in the results with a value of `NA` in `pop_col`.
#' @param following_days Integer or missing (default). Number of days after
#'   releases to return, if `following_releases` is not missing.
#' @describeIn sit_egg_surveys Retrieve egg survey data, possibly filtering
#'   results.
#' @importFrom rlang .data .env
#' @importFrom tidyr everything
#' @export
sit_egg_surveys.sit <- function(
  x,
  area = c('control', 'sit'),
  trap_type = sit_trap_types(x)$label,
  following_releases,
  following_days,
  species = NULL,
  ...) {

  target_periods <- get_periods(following_releases, following_days)

  ## Filter trap areas, trap types, and dates
  x_sub <-
    x %>%
    dm_filter(traps = .data$area %in% .env$area) %>%
    dm_filter(trap_types = .data$label %in% .env$trap_type) %>%
    dm_filter(
      egg_surveys =
        is_date_within(
          .data$datetime_end,
          .env$target_periods$starts,
          .env$target_periods$ends
        )
    )

  ## Filters for species if required
  if (!is.null(species)) {
    if (!'species' %in% names(x[['egg_surveys']])) {
      stop(
        "There is no variable `species` in the egg surveys."
      )
    }
    x_sub <- dm_filter(
      x_sub,
      egg_surveys = (.data$species %in% .env$species)
    )
  }

  ans <- x_sub %>%
    # dm_get_tables() %>% `[`(c("traps", "trap_types"))
    ## Select and rename columns from relevant tables
    dm_select('release_events', .data$id, pop_col = .data$colour) %>%
    dm_select('traps', .data$id, trap_code = .data$code) %>%
    dm_flatten_to_tbl('egg_surveys') %>%
    select(-.data$trap_id, -.data$geometry) %>%
    select(.data$id, .data$trap_code, everything()) %>%
    set_class('sit_egg_surveys')

  return(ans)
}


#' Default table of trap types
#' @return A `data.frame` with the table of default trap types.
#' @keywords internal
default_trap_types <- function() {
  ans <- data.frame(
    id = 1:3,
    name = c("Ovitrap", "BG-Sentinel", "Human Landing Catch"),
    label = c("OVT", "BGS", "HLC"),
    stage = c("egg", "adult", "adult"),
    description = NA
  )

  class(ans) <- c("sit_trap_types", class(ans))

  return(ans)
}

#' Validator
#'
#' Check the internal consistency of a `sit` object.
#'
#' Performs the following checks:
#'
#' - Consistency of survey dates. Surveys of sterile males must precede their
#' release date. _Dummy_ surveys with 0 counts will be silently removed.
#'
#' @param x An object.
#'
#' @return Returns the `sit`, invisibly, after finishing all checks.
#'
#' @keywords internal
#' @importFrom rlang .data .env
#' @export
validate_sit <- function(x) {

  stopifnot(inherits(x, 'sit'))

  ## Surveys of sterile males must precede their release date.
  ages <- survey_ages(sit_adult_surveys(x), sit_revents(x))

  inconsistent_dates <- ages$age < 0 & ages$n > 0
  if ( (sid <- sum(inconsistent_dates)) > 0) {
    stop(
      "Inconsistent dates. ",
      "There are ", sid,
      " counts of sterile males surveyed at a date earlier than their release.")
  }

  ## Remove dummy surveys
  dummy_surveys <- ages$age <= 0  # must be n = 0
  if (any(dummy_surveys)) {
    dummy_surv_ids <- ages$id[dummy_surveys]
    # x <- dm_filter(x, adult_surveys = (!.data$id %in% .env$dummy_surv_ids))
    x <- dm_mutate_tbl(
      x,
      adult_surveys = x[["adult_surveys"]][
        -match(dummy_surv_ids, x[["adult_surveys"]][["id"]]),
      ]
    )
    class(x) <- c("sit", class(x))
  }

  invisible(x)
}

1		#' Survival statistics
2		#'
3		#' Compute the _Probability of Daily Survival_ (PDS), the _Average Life
4		#' Expectancy_ `ALE`, _Recapture Rate_ `RR`, _Survival Rate_ `SR`), by
5		#' population (by default, unless `pool = TRUE`)
6		#'
7		#'
8		#' @inheritParams sit_flight_range
9		#'
10		#' @return A `sit_survival` object, which is a data.frame with additional
11		#' attributes which can be plotted.
12		#' @export
13		#'
14		#' @examples
15		#' (surv_prototype <- sit_survival(sit_prototype))
16		#' plot(surv_prototype)
17		sit_survival <- function(
18		x,
19		pool = FALSE,
20		spatial_adjustment = TRUE,
21		# weights = inverse_density(sit_traps(x, stage = 'adult', area = 'sit')),
22		following_releases,
23		following_days,
24		species = NULL
25		) {
26
27		## TODO: I think I don't really need the dispersal data, but only counts
28		## aggregated by population and age. Check this.
29	4x	dispersal_dat <- dispersal(
30	4x	x,
31	4x	spatial_adjustment = spatial_adjustment,
32	4x	following_releases = following_releases,
33	4x	following_days = following_days,
34	4x	species = species
35		)
36
37
38	3x	populations <- unique(dispersal_dat$population)
39	3x	N <- sit_revents(x)$n[match(populations, sit_revents(x)$colour)]
40	3x	if (pool) {
41	1x	dispersal_dat$population <- "pooled"
42	1x	N <- sum(N)
43		}
44
45	3x	surv_dat <- aggregate(n ~ population + age, data = dispersal_dat, FUN = sum)
46	3x	surv_dat$log_catch = log(1 + surv_dat$n)
47
48	3x	probability_daily_survival <- function(x) {
49	5x	slope <- coef(x)[2]
50	5x	exp(slope)
51		}
52
53	3x	recapture_rate <- function(x, N) {
54	5x	alpha <- coef(x)[1]
55	5x	exp(alpha) / (N + exp(alpha))
56		}
57
58	3x	survival_rate <- function(x, rr, d = 1) {
59	5x	beta <- coef(x)[2]
60	5x	exp(beta) / ( (1-rr)^(1/d))
61		}
62
63	3x	fm_pop <- lapply(
64	3x	split(surv_dat, surv_dat$population),
65	3x	function(x) lm(log_catch ~ age, data = x)
66		)
67
68	3x	pds <- min(vapply(fm_pop, probability_daily_survival, numeric(1)), 1 - .Machine$double.eps)
69	3x	rr <- mapply(recapture_rate, fm_pop, N)
70
71	3x	ans <- data.frame(
72	3x	population = names(fm_pop),
73	3x	N_released = N,
74	3x	PDS = pds,
75	3x	ALE = -1 / log(pds),
76	3x	RRx1e3 = rr * 1e3,
77	3x	SR = mapply(survival_rate, fm_pop, rr),
78	3x	row.names = NULL
79		)
80
81	3x	return(
82	3x	structure(
83	3x	ans,
84	3x	surv_data = surv_dat,
85	3x	class = c("sit_survival", "data.frame")
86		)
87		)
88		}
89
90
91		#' @importFrom knitr kable
92		#' @export
93		print.sit_survival <- function(x, digits = 2, ...) {
94	3x	cli::cli_h1("Survival")
95
96	3x	if (inherits(x, 'data.frame')) {
97	3x	print(
98	3x	knitr::kable(
99	3x	x,
100	3x	format = "simple",
101	3x	digits = digits
102		)
103		)
104		} else {
105	!	cat("\n")
106	!	cli::cli_alert_info("Estimated average value: {round(x, digits)}")
107		}
108	3x	cat("\n")
109		}
110
111		#' @export
112		plot.sit_survival <- function(x, y, ...) {
113
114	3x	surv_dat <- attr(x, 'surv_data')
115
116	3x	model_formula <- log_catch ~ age*population
117
118	3x	n_populations <- nrow(x)
119
120	3x	if(n_populations == 1L) {
121	2x	model_formula <- log_catch ~ age
122	2x	surv_dat$population <- "black"
123		}
124
125	3x	survm <- lm(model_formula, data = surv_dat)
126
127	3x	surv_dat$pred <- predict(survm)
128
129	3x	plot(
130	3x	log_catch ~ age,
131	3x	data = surv_dat,
132	3x	pch = 19,
133	3x	col = surv_dat$population,
134	3x	main = "Survival",
135	3x	xlab = "Age (days since release)",
136	3x	ylab = "Log-catch"
137		)
138
139	3x	for (d in split(surv_dat, surv_dat$population)) {
140	5x	lines(d$age, d$pred, col = d$population)
141		}
142
143	3x	invisible(surv_dat)
144		}
145

1		#' Import or extract trap types
2		#'
3		#' Extracts the table of trap types used in a `sit` experiment, or defines
4		#' a table of trap types.
5		#'
6		#' __Returns__ the table of trap types used in the `sit` object `x`.
7		#'
8		#' __Defines__ a table of trap types by passing a `data.frame` with variables
9		#' `name`, `label`, `stage` (either `egg` or `adult`) and optionally
10		#' `description`. Any other variable is ignored with a warning.
11		#'
12		#' __Returns__ the default table of trap types by calling `sit_trap_types()`
13		#' without any argument.
14		#'
15		#' Check the default trap types with `sit_trap_types()`.
16		#' You can store the table into an object and make edits as for a regular
17		#' `data.frame`. E.g. add another trap type, or edit labels or descriptions.
18		#' You can use the edited object as an input to the `trap_types` arguments in
19		#' [sit_traps()].
20		#'
21		#' @param x A `sit` object, a `data.frame`-like object or empty.
22		#'
23		#' @return Table of trap types used in the `sit` experiment `x`. If `x` is a
24		#' `data.frame`, returns it as a table of trap types after some verifications.
25		#' If `x` is missing, returns the current table of trap types.
26		#' @export
27		#'
28		#' @examples
29		#' sit_trap_types() # default trap types.
30		#'
31		#' # Define new trap types
32		#' my_traps <- data.frame(name = "My New Trap", label = "MNT", stage = "adult")
33		#' sit_trap_types(my_traps)
34		sit_trap_types <- function(x) {
35
36	162x	if (missing(x)) {
37		## Return default table of trap types
38	37x	return(default_trap_types())
39		}
40
41	125x	UseMethod("sit_trap_types")
42
43		}
44
45
46		#' @describeIn sit_trap_types Extract trap data
47		#' @export
48		sit_trap_types.sit <- function(x) {
49
50	114x	return(x[['trap_types']])
51		}
52
53		#' @describeIn sit_trap_types Import trap data
54		#' @export
55		sit_trap_types.data.frame <- function(x) {
56
57	10x	mandatory_vars <- c("name", "label", "stage")
58
59		## id is optional, but ignored
60		## description is optional and used
61	10x	optional_vars <- c("id", "description")
62
63		## Check mandatory variables
64	10x	if (any(miss <- (!mandatory_vars %in% names(x)))) {
65	1x	stop(
66	1x	glue(
67	1x	"Required variables ",
68	1x	"{vars_collapse(mandatory_vars[miss], intros = '')}",
69	1x	" not found in {substitute(x)}. Please verify."
70		)
71		)
72		}
73
74		## Check valid values in stages
75	9x	xtages <- setdiff(unique(x$stage), c("egg", "adult"))
76	9x	if (length(xtages) > 0) {
77	1x	stop(
78	1x	glue(
79	1x	"Trap-type ",
80	1x	"{vars_collapse(xtages, intro1 = 'stage ', intros = 'stages ')} ",
81	1x	"not valid. Only 'egg' and 'adult' are valid values. Case sensitively."
82		)
83		)
84		}
85
86		## Warning about ignoring id
87	8x	if ('id' %in% names(x)) {
88	1x	warning(
89	1x	glue(
90	1x	"Variable 'id' from {substitute(x)} ignored. ",
91	1x	"Trap-type ids are overriden internally."
92		)
93		)
94		}
95
96		## Warning about ignoring additional variables
97	8x	extra_vars <- setdiff(names(x), c(mandatory_vars, optional_vars))
98	8x	if (length(extra_vars) > 0) {
99	1x	warning(
100	1x	glue(
101	1x	"{vars_collapse(extra_vars)} from {substitute(x)} ignored."
102		)
103		)
104		}
105
106		## Make empty description if not provided
107	8x	if (!'description' %in% names(x)) {
108	6x	x$description <- NA
109		}
110
111	8x	ans <- cbind(
112	8x	id = seq.int(nrow(x)),
113	8x	as.data.frame(x[, c(mandatory_vars, "description")])
114		)
115	8x	class(ans) <- c("sit_trap_types", class(ans))
116
117	8x	return(ans)
118		}
119
120		#' @export
121		sit_trap_types.default <- function(x) {
122
123	1x	stop(glue(
124	1x	"{deparse(substitute(x))} must either be a 'sit' object to extract",
125	1x	" trap types from, a data.frame to define trap types",
126	1x	" or missing to get the currently defined trap types."
127		))
128		}

1		#' Size of the Wild Population
2		#'
3		#' Estimate the size of the wild population, via the Lincoln Index.
4		#'
5		#' This provides multiple estimates of the population size, one for each age $t$
6		#' and for each released population, which can be averaged to compute a final
7		#' estimate based on all the available data.
8		#'
9		#' The calculation needs to be split by released population. First, because
10		#' survival rates could be different among populations, but also because the
11		#' number $n_t$ of captured wild males at a given date can be compared against
12		#' the number of captured sterile males $m_t$ at different ages.
13		#'
14		#' @inheritParams sit_flight_range
15		#'
16		#' @return A `sit_wild_size` object.
17		#' @importFrom stats sd
18		#' @export
19		#'
20		#' @examples
21		#' sit_wild_size(sit_prototype)
22		sit_wild_size <- function(
23		x,
24		pool = FALSE,
25		species = NULL
26		) {
27
28		## Survival
29		## Using spatially-adjusted captures by default here.
30		## But this is irrelevant for the survival parameters.
31	!	surv <- sit_survival(
32	!	x,
33	!	spatial_adjustment = FALSE, # This should be removed from sit_survival
34	!	pool = pool,
35	!	species = species)[, c("population", "N_released", "SR")]
36
37
38
39		## Adult captures by population and age
40		# ## Aggregated captures by population and age
41		# age_captures <- attr(surv, "surv_data")
42
43		## Age surveys from areal releases as well, and without spatial adjustment.
44	!	ages <- survey_ages(sit_adult_surveys(x, area = "sit"), sit_revents(x))
45
46	!	if (pool) ages$pop_col <- "pooled"
47
48		## Remove ages for which there are no observations
49		## Last age for which there is at least one observation of a sterile male
50	!	max_age <- max_observed_age(ages$n, ages$age)
51	!	ages <- ages[ages$age <= max_age, ]
52	!	names(ages) <- gsub("^pop_col$", "population", names(ages))
53
54	!	ages <- aggregate(n ~ population + age + datetime_end, data = ages, sum)
55
56
57		## Wild counts aggregated across traps
58	!	asurvs <- sit_adult_surveys(x, area = "sit")
59	!	wild_survs <- asurvs[asurvs$pop_col == "wild" & asurvs$sex == "male", ]
60
61	!	if (!is.null(species) && species %in% names(wild_survs)) {
62	!	wild_survs <- wild_survs[wild_survs$species %in% species, ]
63		}
64
65	!	wild_survs <- wild_survs[, c("trap_code", "datetime_end", "n")]
66	!	names(wild_survs) <- gsub("^n$", "n_w", names(wild_survs))
67	!	wild_survs <- aggregate(n_w ~ datetime_end, data = wild_survs, sum)
68
69	!	ages <- merge(ages, surv, by = "population")
70	!	pop_dat <- merge(ages, wild_survs, by = c("datetime_end"))
71
72		## Lincoln index
73	!	pop_dat$LI <- with(
74	!	pop_dat,
75	!	lincoln_index(
76	!	released = N_released,
77	!	survival_rate = SR,
78	!	age = age,
79	!	marked_recaptured = n,
80	!	wild_catch = n_w
81		)
82		)
83
84		## Rename columns
85	!	sel_cols <- c(population = "population", survey_date = "datetime_end",
86	!	age = "age", N_released = "N_released", SR = "SR", n = "n",
87	!	n_wild = "n_w", LI = "LI")
88	!	pop_dat <- setNames(pop_dat[, sel_cols], names(sel_cols))
89
90	!	ans <- aggregate(LI ~ population, data = pop_dat, function(x) c(p_est = mean(x), p_sd = sd(x)))
91	!	ans <- data.frame(population = ans$population, as.data.frame(ans$LI))
92
93	!	return(
94	!	structure(
95	!	ans,
96	!	pop_data = pop_dat,
97	!	class = c("sit_wild_size", "data.frame")
98		)
99		)
100
101		}
102
103		#' @importFrom knitr kable
104		#' @export
105		print.sit_wild_size <- function(x, digits = 2, ...) {
106	!	cli::cli_h1("Wild male population size")
107
108	!	if (inherits(x, 'data.frame')) {
109	!	print(
110	!	knitr::kable(
111	!	x,
112	!	format = "simple",
113	!	digits = digits
114		)
115		)
116		} else {
117	!	cat("\n")
118	!	cli::cli_alert_info("Estimated average value: {round(x, digits)}")
119		}
120	!	cat("\n")
121		}
122
123		#' @export
124		plot.sit_wild_size <- function(x, y, ...) {
125
126	!	pop_dat <- attr(x, 'pop_data')
127
128	!	boxplot(
129	!	LI ~ population,
130	!	data = pop_dat,
131	!	horizontal = TRUE,
132	!	col = rev(unique(pop_dat$population)),
133	!	main = "Wild male population size",
134	!	xlab = "Lincoln index",
135	!	ylab = "Released population"
136		)
137
138	!	invisible(pop_dat)
139		}
140
141
142
143
144		#' Lincoln index
145		#'
146		#' Estimate the size of the wild-male population.
147		#'
148		#' A simple estimate is obtained as follows (Thompson 2012, Ch. 18). Let the
149		#' total captures at a day $t$ be the sum of $m_t$ marked and $n_t$ wild
150		#' mosquitoes. Assuming that the proportion of marked individuals in the sample
151		#' equals that in the population of size $P$: \deqn{\frac{m_t}{m_t + n_t} =
152		#' \frac{M_t}{M_t + P},} where $M_t = R\,S^{a_t}$ is the number of marked
153		#' individuals captured at time $t$, with $R$ the number of released adults, $S$
154		#' the daily survival rate and $a_t$ the number of days since release (_age_).
155		#' I.e., the number of marked individuals at time $t$ is the remaining number
156		#' from those released that survived for $a_t$ days.
157		#'
158		#' The __Lincoln Index__ (_a.k.a._ the __Petersen estimator__) has been used as
159		#' a simple estimate of the wild __male__ population size, assuming that the
160		#' survival rate of an individual remains constant. Here we use a _modified_
161		#' version that corrects for small samples and compensates for daily survival.
162		#' \deqn{P_t = R\,S^{a_t}\,(n_t + 1) / (m_t + 1).}
163		#'
164		#' The values of $R$, $n_t$, $m_t$ and $t$ can be gathered from the adult
165		#' surveys data. The calculation required the estimation of the survival rate
166		#' $S$.
167		#'
168		#' @param released Numeric. Number of sterile males released.
169		#' @param survival_rate Numeric. Estimated survival rate.
170		#' @param age Numeric. Number of days since release.
171		#' @param marked_recaptured Numeric. Number of re-captured sterile males.
172		#' @param wild_catch Numeric. Number of wild males captured.
173		#'
174		#' @export
175		#' @references Thompson, Steven K. 2012. Sampling. 3rd ed. Wiley Series in
176		#' Probability and Statistics. Hoboken, N.J: Wiley.
177		#' @examples lincoln_index(
178		#' released = 1e4,
179		#' survival_rate = .80,
180		#' age = 1,
181		#' marked_recapture = 15,
182		#' wild_catch = 5
183		#' )
184		lincoln_index <- function(
185		released,
186		survival_rate,
187		age,
188		marked_recaptured,
189		wild_catch
190		) {
191	!	released * survival_rate^age * (wild_catch + 1) / (marked_recaptured + 1)
192		}

1		#' Import or extract release events
2		#'
3		#' Provide a `sf` object with a POINT at each record corresponding to a
4		#' __release point__ or a `data.frame` for a __areal release__, specifying which
5		#' variables contain the required information. Or, extract the release events
6		#' used in a `sit` object.
7		#'
8		#' @param x Object of class `sf` with POINT elements to be _imported_, or a
9		#' `data.frame` with information about areal releases, or
10		#' object of class `sit` to _extract_ release events from.
11		#' @param date Character string representing a date and optionally the time, in
12		#' the [RFC 3339 format](https://www.ietf.org/rfc/rfc3339.txt) (a variation of
13		#' the [ISO 8601 format](https://en.wikipedia.org/wiki/ISO_8601)), i.e.,
14		#' `2021-12-31` or `2019-11-23 15:00`.
15		#' @param colour Character string. Colour of the release.
16		#' @param species Optional character. Released species.
17		#' @param n Numeric. Natural number. Number of individuals released.
18		#' @param ... Used to pass arguments to specific method.
19		#'
20		#' @return A object of class `sit_revents` which can be used in [sit()].
21		#' `sit_revents` objects can be concatenated together with `c()`, in order to
22		#' import a mixture of point and areal releases into a `sit` object. See
23		#' examples.
24		#'
25		#' @import sf
26		#' @export
27		#'
28		#' @family importing
29		#' @seealso [sit()]
30		#'
31		#' @examples
32		#' point_releases <- sit_revents(
33		#' sf::st_as_sf(
34		#' data.frame(
35		#' x = 1:3,
36		#' y = 3:1,
37		#' date = c("2019-11-25", "2019-12-01", "2019-12-13"),
38		#' colour = c("yellow", "red", "blue"),
39		#' n = 1e4
40		#' ),
41		#' coords = c("x", "y")
42		#' )
43		#' )
44		#'
45		#' areal_release <- sit_revents(
46		#' data.frame(
47		#' x = 1, y = 1,
48		#' date = "2019-12-21",
49		#' colour = "pink",
50		#' n = 1e4
51		#' )
52		#' )
53		#'
54		#' (c(point_releases, areal_release)) # also a `sit_revent` object
55		#'
56	117x	sit_revents <- function(x, ...) UseMethod("sit_revents")
57
58
59		#' @describeIn sit_revents Imports point release data.
60		#' @export
61		sit_revents.sf <- function(
62		x,
63		date = "date",
64		colour = "colour",
65		n = "n",
66		species = NULL,
67		...
68		) {
69
70		# The features must be all POINTs
71	37x	stopifnot(st_geometry_type(x, by_geometry = FALSE) == "POINT")
72
73		## Build the internal representation
74	37x	ans <- build_internal_df(
75	37x	x,
76	37x	varnames = c(date, colour, n, species),
77	37x	internal_names = c("date", "colour", "n", "species"[!is.null(species)])
78		)
79
80		## Interpret dates
81	37x	ans$date <- as.POSIXct(ans$date)
82
83		## Check duplicated releases
84	36x	if (any(idx <- duplicated(ans))) {
85	1x	stop(
86	1x	glue(
87	1x	"{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
88	1x	"duplicated. Please fix."
89		)
90		)
91		}
92
93		## Add unique internal numeric identifier of the release event
94	35x	ans <- cbind(id = seq.int(nrow(x)), ans)
95
96		## Table of __unique__ release_sites
97	35x	sites <- st_geometry(ans)
98	35x	unique_sites <- sites[!duplicated(sites)]
99	35x	release_sites <- st_sf(id = seq_along(unique_sites), geometry = unique_sites)
100
101		## Replace the geometries by the corresponding site_id
102	35x	ans <- st_drop_geometry(ans)
103	35x	ans$site_id <- as.integer(st_equals(sites, unique_sites))
104
105		## Reorder columns
106	35x	ans <- ans[
107	35x	, c("id", "site_id", "date", "colour", "species"[!is.null(species)], "n")
108		]
109
110	35x	attr(ans, "release_sites") <- release_sites
111	35x	class(ans) <- c("sit_revents", class(ans))
112	35x	return(ans)
113		}
114
115		#' @describeIn sit_revents Imports areal release data.
116		#' @export
117		sit_revents.data.frame <- function(
118		x,
119		date = "date",
120		colour = "colour",
121		n = "n",
122		species = NULL,
123		...
124		) {
125
126		## Build the internal representation
127	17x	ans <- build_internal_df(
128	17x	x,
129	17x	varnames = c(date, colour, n, species),
130	17x	internal_names = c("date", "colour", "n", "species"[!is.null(species)])
131		)
132
133		## Interpret dates
134	17x	ans$date <- as.POSIXct(ans$date)
135
136		## Check duplicated releases
137	17x	if (any(idx <- duplicated(ans))) {
138	1x	stop(
139	1x	glue(
140	1x	"{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
141	1x	"duplicated. Please fix."
142		)
143		)
144		}
145
146		## Add unique internal numeric identifier of the release event
147	16x	ans <- cbind(id = seq.int(nrow(x)), ans)
148
149		## Table of __unique__ release_sites
150		## All areal releases are performed from "the same" release_site
151		## The geometry is a POINT EMPTY.
152	16x	release_sites <- st_sf(
153	16x	id = 1L,
154	16x	geometry = st_sfc(st_point()),
155	16x	crs = NA_crs_
156		)
157
158		## Replace the geometries by the corresponding site_id
159	16x	ans$site_id <- 1L
160
161		## Reorder columns
162	16x	ans <- ans[
163	16x	, c("id", "site_id", "date", "colour", "species"[!is.null(species)], "n")
164		]
165
166	16x	attr(ans, "release_sites") <- release_sites
167	16x	class(ans) <- c("sit_revents", class(ans))
168	16x	return(ans)
169		}
170
171		#' @describeIn sit_revents Concatenate release events.
172		#' @export
173		c.sit_revents <- function(...) {
174
175	18x	lst <- list(...)
176	18x	classes_revents <- vapply(lst, inherits, TRUE, 'sit_revents')
177
178		## All elements need to be of class `sit_revents`
179	18x	if (!all(classes_revents)) {
180	1x	stop(
181	1x	glue(
182	1x	"{vars_collapse(
183	1x	seq_along(lst)[!classes_revents], 'Argument ', 'Arguments ',
184	1x	quote = ''
185	1x	)} need to be created with `sit_revents()`."
186		)
187		)
188		}
189
190		## row-bind release sites
191	17x	rs <- lapply(lst, attr, 'release_sites')
192		## Check whether all release sites are areal
193		## I.e. have empty points and NA bboxes
194	17x	areals <- vapply(lapply(rs, st_bbox), is.na, TRUE)
195	17x	all_areals <- all(areals)
196	17x	if (all_areals) {
197		## If all_areals, suppressWarnings related to establishing the bbox
198		## and reset resulting bbox as NA
199	1x	release_sites <- suppressWarnings(combine_rows(rs))
200	1x	attr(release_sites[[attr(release_sites, "sf_column")]], "bbox") <-
201		## The bbox does not need to have a crs. This causes spurious warnings.
202	1x	structure(NA_bbox_, crs = NULL)
203		} else {
204
205		## Check compatibility of CRSs
206	16x	crs_list <- lapply(rs, st_crs)
207	16x	if (length(crs <- unique(crs_list[!areals])) > 1) {
208	1x	stop("Arguments have different coordinate reference systems. ",
209	1x	"Please use st_transform() to convert to a common CRS.")
210		}
211
212		## Set areal CRSs to the common CRS
213	15x	for (i in which(areals)) {
214	13x	st_crs(rs[[i]]) <- crs[[1]]
215		}
216
217	15x	release_sites <- combine_rows(rs)
218		}
219
220
221		## recode site_ids
222	16x	element_id <- rep(seq_along(rs), times = vapply(rs, nrow, 1))
223	16x	new_codes <- split(seq_along(element_id), element_id)
224	16x	for (i in seq_along(lst)) {
225	30x	lst[[i]]$site_id <- new_codes[[i]][lst[[i]]$site_id]
226		}
227
228		## row-bind recoded tables of release events
229		## and reset event id
230	16x	ans <- combine_rows(lst)
231
232		## Attach attribute table of release sites
233	16x	attr(ans, 'release_sites') <- release_sites
234
235	16x	return(ans)
236		}
237
238
239		#' @describeIn sit_revents Extracts release events.
240		#' @param type Character vector. Which _types_ of release events to retrieve.
241		#' Either `point` or `areal` or both (default).
242		#' @export
243		sit_revents.sit <- function(x, type = c('point', 'areal'), ...) {
244
245	59x	ans <-
246	59x	x %>%
247		## Filter release event types
248	59x	dm_filter(
249	59x	release_sites = (release_types(x[['release_sites']]) %in% .env$type)
250		)
251
252
253	59x	ans <- structure(
254	59x	ans[['release_events']],
255	59x	release_sites = ans[['release_sites']]
256		) %>%
257	59x	set_class('sit_revents')
258
259
260	59x	return(ans)
261		}
262
263		#' @export
264		sit_revents.sit_revents <- function(x, type = c('point', 'areal'), ...) {
265
266	4x	sites <- attr(x, "release_sites")
267	4x	sites_idx <- release_types(sites) %in% type
268
269	4x	ans <- x[x$site_id %in% sites$id[sites_idx], ]
270	4x	attr(ans, "release_sites") <- sites[sites_idx, ]
271
272	4x	return(ans)
273		}
274
275		#' @export
276		print.sit_revents <- function(
277		x,
278		...
279		) {
280
281	1x	ans <- as.data.frame.sit_revents(x)
282	1x	print(ans)
283	1x	invisible(ans)
284		# ans
285		}
286
287
288		#' @export
289		as.data.frame.sit_revents <- function(x, ...) {
290
291	32x	rs <- attr(x, 'release_sites')
292	32x	names(rs) <- sub("id", "site_id", names(rs))
293
294		## Release type
295	32x	rs$type <- release_types(rs)
296
297		## Unclass
298	32x	class(x) <- setdiff(class(x), 'sit_revents')
299
300	32x	ans <- merge(
301	32x	x, rs, by = "site_id", sort = FALSE
302		)
303
304	32x	cols <- c('id', 'type', setdiff(names(ans), c('id', 'type')))
305	32x	ans <- ans[order(ans$id), cols]
306
307	32x	return(ans)
308		}
309
310		#' @export
311		st_as_sf.sit_revents <- function(x, ...) {
312
313	!	st_as_sf(as.data.frame(x))
314		}
315
316
317		#' @export
318		summary.sit_revents <- function(object, ...) {
319
320	1x	rs <- attr(object, 'release_sites')
321	1x	n_areal <- sum(sf::st_is_empty(rs))
322	1x	n_point <- nrow(rs) - n_areal
323
324	1x	ans <- structure(
325	1x	list(
326	1x	sit_revents = object,
327	1x	release_sites = rs,
328	1x	n = c(point = n_point, areal = n_areal)
329		),
330	1x	class = "summary.sit_revents"
331		)
332
333	1x	return(ans)
334		}
335
336		#' @export
337		print.summary.sit_revents <- function(x, ...) {
338
339	1x	cat("Number of point releases:", x$n[["point"]], "\n")
340	1x	cat("Number of areal releases:", x$n[["areal"]], "\n")
341	1x	cat("\n")
342
343	1x	cat("Release events:\n")
344	1x	print(x$sit_revents)
345	1x	cat("\n")
346
347	1x	cat("Release sites:\n")
348	1x	print(x$release_sites)
349		}

1		#' Import or extract trap information
2		#'
3		#' Provide a `sf` object with a POINT at each record corresponding to a trap
4		#' location, specifying which variables contain the required information. Or
5		#' define a set of traps without spatial coordinates. Or, extract the trap data
6		#' used in a `sit` object.
7		#'
8		#' @param x Object of class `sf` with POINT elements to be _imported_, or a
9		#' character-coerced vector _trap codes_ to be defined, or an object of class
10		#' `sit` to _extract_ trap information from.
11		#' @param id Character. Variable name with a custom code for the trap. Jointly
12		#' with `type`, identifies the trap uniquely.
13		#' @param type Character. Variable name with __valid__ values of trap types
14		#' (i.e., any of `trap_types$label`). For the extraction method, a character
15		#' vector of trap types to filter upon.
16		#' @param area Character. Variable name with values either `control` or `sit`.
17		#' For the extraction method, a character vector of areas to filter upon.
18		#' @param stage Character vector of stages (i.e., `adult` or `egg`) to filter
19		#' upon extraction.
20		#' @param label Optional character. Variable name with custom text associated to
21		#' the traps.
22		#' @param trap_types Table of trap types, created with [sit_trap_types()].
23		#' @param ... Used to pass arguments to specific method.
24		#'
25		#' @return A object of class `sit_traps` which can be used in [sit()].
26		#' @import sf
27		#' @export
28		#'
29		#' @family importing
30		#' @seealso [sit()]
31		#'
32		#' @examples
33		#'
34		#' ## Build spatialised traps from a `sf` object
35		#' traps_sf <- sf::st_as_sf(
36		#' data.frame(
37		#' x = 1:3,
38		#' y = 3:1,
39		#' "Trap.Id" = letters[1:3],
40		#' Type = "BGS",
41		#' area = "sit"
42		#' ),
43		#' coords = c("x", "y")
44		#' )
45		#' sit_traps(traps_sf, id = "Trap.Id", type = "Type", area = "area")
46		#'
47		#' ## Build non-spatialised traps from a vector of codes
48		#' sit_traps(paste0("C", 1:5), area = 'control', type = 'OVT')
49		#'
50		#' ## Retrieve traps in a `sit` object
51		#' ## sit_traps(sit_prototype)
52	58x	sit_traps <- function(x, ...) UseMethod("sit_traps")
53
54
55		#' @describeIn sit_traps Imports spatialised trap data
56		#' @export
57		sit_traps.sf <- function(
58		x,
59		id = "id",
60		type = "type",
61		area = "area",
62		label,
63		trap_types = sit_trap_types(),
64		...
65		) {
66
67		## Fill-in label variable if missing
68	40x	if(missing(label)) {
69	39x	if (!"label" %in% names(x)) {
70	36x	x$label <- NA
71		}
72	39x	label = "label"
73		}
74
75		## Check trap_types
76	40x	if (!inherits(trap_types, 'sit_trap_types')) {
77	1x	stop("Please use `sit_trap_types()` to create a valid table of trap types.")
78		}
79
80
81		## Build the internal representation
82	39x	ans <- build_internal_df(
83	39x	x,
84	39x	varnames = c(id, type, area, label),
85	39x	internal_names = c("code", "type_id", "area", "label")
86		)
87
88		## Check duplicated traps
89	38x	if (any(idx <- duplicated(ans))) {
90	1x	stop(
91	1x	glue(
92	1x	"{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
93	1x	"duplicated. Please fix."
94		)
95		)
96		}
97
98		## Add unique internal numeric identifier
99	37x	ans <- cbind(id = seq.int(nrow(x)), ans)
100
101		## Check trap types and replace by codes
102	37x	obs_types <- unique(ans$type_id)
103	37x	xtypes <- setdiff(obs_types, trap_types$label) # invalid types
104	37x	if (length(xtypes) > 0) {
105	1x	stop(
106	1x	glue(
107	1x	"Invalid trap {vars_collapse(xtypes, 'type ', 'types ')} ",
108	1x	"in variable {vars_collapse('type', '')}.
109	1x	Valid types are the values of 'label' ",
110	1x	"in the object passed to `trap_types`, ",
111	1x	"i.e.: {vars_collapse(trap_types$label, '', '')}. ",
112	1x	"Please verify."
113		)
114		)
115		}
116	36x	ans$type_id <- trap_types$id[match(ans$type_id, trap_types$label)]
117
118	36x	attr(ans, "trap_types") <- trap_types
119	36x	class(ans) <- c("sit_traps", class(ans))
120	36x	return(ans)
121		}
122
123
124		#' @describeIn sit_traps Imports non-spatialised trap data
125		#' @export
126		sit_traps.character <- function(
127		x,
128		type,
129		area,
130		trap_types = sit_trap_types(),
131		...
132		) {
133
134		## Build a `sf` representation with POINT EMPTY geometry
135	3x	dat_sf <- st_sf(
136	3x	id = x,
137	3x	type = rep(type, length(x)),
138	3x	area = rep(area, length(x)),
139	3x	geometry = rep(st_sfc(st_point()), length(x)),
140		...,
141	3x	crs = NA_crs_
142		)
143
144	3x	ans <- sit_traps.sf(dat_sf, trap_types = trap_types)
145
146	3x	return(ans)
147		}
148
149
150		#' @describeIn sit_traps Imports non-spatialised trap data
151		#' @export
152		sit_traps.numeric <- function(x, ...) {
153	1x	sit_traps.character(as.character(x), ...)
154		}
155
156
157		#' @describeIn sit_traps Combine traps.
158		#' @export
159		c.sit_traps <- function(...) {
160
161	4x	lst <- list(...)
162	4x	classes_traps <- vapply(lst, inherits, TRUE, 'sit_traps')
163
164		## All elements need to be of class `sit_traps`
165	4x	if (!all(classes_traps)) {
166	1x	stop(
167	1x	glue(
168	1x	"{vars_collapse(
169	1x	seq_along(lst)[!classes_traps], 'Argument ', 'Arguments ',
170	1x	quote = ''
171	1x	)} need to be created with `sit_traps()`."
172		)
173		)
174		}
175
176		## Combine trap types.
177	3x	tt <- lapply(lst, attr, 'trap_types')
178	3x	trap_types <- combine_rows(tt, recoder = TRUE)
179
180		## Recode trap type ids
181	3x	for (i in seq_along(lst)) {
182	6x	lst[[i]]$type_id <- attr(trap_types, 'recoder')(lst[[i]]$type_id, i)
183		}
184	3x	attr(trap_types, 'recoder') <- NULL
185
186
187		## Check whether all traps are non-spatial
188		## I.e. have empty points and NA bboxes
189	3x	non_spatial <- vapply(lapply(lst, st_bbox), is.na, TRUE)
190	3x	all_non_spatial <- all(non_spatial)
191	3x	if (all_non_spatial) {
192		## If all_non_spatial, suppressWarnings related to establishing the bbox
193		## and reset resulting bbox as NA
194	1x	ans <- suppressWarnings(combine_rows(lst))
195	1x	attr(ans[[attr(ans, "sf_column")]], "bbox") <-
196		## The bbox does not need to have a crs. This causes spurious warnings.
197	1x	structure(NA_bbox_, crs = NULL)
198		} else {
199
200		## Check compatibility of CRSs
201	2x	crs_list <- lapply(lst, st_crs)
202	2x	if (length(crs <- unique(crs_list[!non_spatial])) > 1) {
203	1x	stop("Arguments have different coordinate reference systems. ",
204	1x	"Please use st_transform() to convert to a common CRS.")
205		}
206
207		## Set non-spatial CRSs to the common CRS
208	1x	for (i in which(non_spatial)) {
209	1x	st_crs(lst[[i]]) <- crs[[1]]
210		}
211
212		## row-bind recoded tables of release events
213		## and reset event id
214	1x	ans <- combine_rows(lst)
215		}
216
217		## Reset class order
218		## Assignments to 'sf' objects result in class 'sf' becoming first in the list
219		## of classes.
220	2x	ans <- set_class(ans, 'sit_traps')
221
222		## Attach attribute table of release sites, possibly filtering.
223	2x	attr(ans, 'trap_types') <- trap_types
224
225	2x	return(ans)
226		}
227
228
229		#' @describeIn sit_traps Extracts trap data
230		#' @importFrom rlang .data
231		#' @export
232		sit_traps.sit <- function(
233		x,
234		type = sit_trap_types(x)$label,
235		stage = unique(sit_trap_types(x)$stage),
236		area = unique(x$traps$area),
237		...
238		) {
239
240
241	18x	xf <- x %>%
242	18x	dm_filter(traps = (.data$area %in% .env$area) ) %>%
243	18x	dm_filter(
244	18x	trap_types = (.data$stage %in% .env$stage & .data$label %in% .env$type)
245		)
246
247	18x	ans <- structure(xf[['traps']], trap_types = xf[['trap_types']])
248
249		## Reset class order
250		## Assignments to 'sf' objects result in class 'sf' becoming first in the list
251		## of classes.
252	18x	ans <- set_class(ans, 'sit_traps')
253
254	18x	return(ans)
255		}
256
257		#' @export
258		print.sit_traps <- function(
259		x,
260		...
261		) {
262
263	1x	ans <- as.data.frame.sit_traps(x)
264	1x	print(ans)
265	1x	invisible(ans)
266		# ans
267		}
268
269
270		#' @export
271		as.data.frame.sit_traps <- function(x, ...) {
272
273		## trap types renamed into a separate object to prevent issues with
274		## lazy evaluation of the other arguments.
275	16x	tt_rn <- structure(
276	16x	attr(x, 'trap_types'),
277	16x	names = paste("type", names(attr(x, 'trap_types')), sep = "_")
278		)
279
280		## Unclass
281	16x	class(x) <- setdiff(class(x), 'sit_traps')
282
283	16x	ans <- merge(
284	16x	x, tt_rn, by = "type_id", sort = FALSE
285	16x	)[, -1]
286	16x	ans <- ans[order(ans$id),]
287
288	16x	return(ans)
289		}
290
291
292
293		#' @importFrom stats aggregate
294		#' @export
295		summary.sit_traps <- function(object, ...) {
296
297	2x	traps <- as.data.frame.sit_traps(object)
298
299	2x	traps$has_coordinates <- !st_is_empty(st_geometry(traps))
300
301	2x	ans <- structure(
302	2x	stats::aggregate(
303		# st_drop_geometry(traps),
304	2x	data.frame(n = seq.int(nrow(traps))),
305	2x	by = list(
306	2x	Area = traps$area,
307	2x	Type = traps$type_label,
308	2x	Coords = traps$has_coordinates
309		),
310	2x	length
311		),
312	2x	class = c("summary.sit_traps", "data.frame")
313		)
314
315	2x	return(ans)
316		}
317
318		#' @export
319		print.summary.sit_traps <- function(x, ...) {
320
321	1x	cat("Number of traps in each group:\n")
322	1x	ans <- as.data.frame(unclass(x))
323	1x	print(ans)
324	1x	invisible(ans)
325		}
326

1
2		#' Estimate competitiveness of sterile males
3		#' @family estimates
4		#'
5		#' Computes the Fried Index, which estimates the competitiveness based on
6		#' observed fertility rates and the sterile-wild male ratio.
7		#'
8		#' It is strongly recommended to estimate competitiveness based on surveys from
9		#' areal releases. The calculation requires the estimation of the sterile-wild
10		#' male ratio from surveys pooled across traps and days, which requires a
11		#' roughly homogeneous distribution of sterile males.
12		#'
13		#' @param x A `sit` object.
14		#' @param following_releases A `sit_release` object with a subset of release
15		#' events or missing (default) for all release events. Use in combination with
16		#' `following_days` to filter surveys of the target populations within a given
17		#' number of days after release. Note that counts of wild populations will
18		#' always be included in the results with a value of `NA` in `pop_col`.
19		#' @param following_days Integer or missing (default). Number of days after
20		#' releases to return, if `following_releases` is not missing.
21		#' @param residual_fertility Numeric. Default: 0.05.
22		#' @param species Character. Filter results for a given species or ignore the
23		#' species if missing.
24		#'
25		#' @return A numeric value.
26		#' @export
27		#'
28		#' @examples
29		#' sit_competitiveness(sit_prototype)
30		sit_competitiveness <- function(
31		x,
32		following_releases = sit_revents(x, type = "areal"),
33		following_days = 7,
34		residual_fertility = 0.01,
35		species = NULL
36		) {
37
38		## Use only surveys within 7 days since the areal releases and consider
39		## only sterile males from those releases.
40	3x	adult_areal <- sit_adult_surveys(
41	3x	x,
42	3x	area = "sit",
43	3x	following_releases = following_releases,
44	3x	following_days = following_days,
45	3x	species = species
46		)
47
48	3x	if (!nrow(adult_areal[adult_areal$pop_col == "wild", ]) > 0 \|
49	3x	!nrow(adult_areal[adult_areal$pop_col != "wild", ]) > 0 ) {
50	!	message(
51	!	"Cannot compute sterile male competitiveness. ",
52	!	"Surveys from adult sterile and wild males following an areal release missing."
53		)
54	!	return(NULL)
55		}
56
57	3x	swmr <- sterile_wild_male_ratio(adult_areal, pool = TRUE)
58
59	3x	egg_control <- sit_egg_surveys(x, area = "control", species = species)
60
61	3x	if (!nrow(egg_control) > 0) {
62	!	message(
63	!	"Cannot compute sterile male competitiveness. ",
64	!	"Surveys from ovitraps in the control area are missing."
65		)
66	!	return(NULL)
67		}
68
69	3x	natural_fertility <- fertility_rate(
70	3x	egg_control,
71	3x	pool = TRUE
72		)
73
74	3x	egg_areal <- sit_egg_surveys(
75	3x	x,
76	3x	area = "sit",
77	3x	following_releases = following_releases,
78	3x	following_days = following_days,
79	3x	species = species
80		)
81
82	3x	if (!nrow(egg_areal) > 0) {
83	!	message(
84	!	"Cannot compute sterile male competitiveness. ",
85	!	"Surveys from ovitraps following an areal release missing."
86		)
87	!	return(NULL)
88		}
89
90	3x	sit_fertility <- fertility_rate(egg_areal, pool = TRUE)
91
92	3x	res <- fried_index(
93	3x	sterile_wild_mr = swmr,
94	3x	natural_fertility = natural_fertility,
95	3x	sit_fertility = sit_fertility,
96	3x	residual_fertility = residual_fertility
97		)
98
99	3x	comp_table <- data.frame(
100	3x	Component = c(
101	3x	"Sterile-wild male ratio", "Natural fertility",
102	3x	"Observed fertility in SIT area", "Sterile fertility (assumed)"
103		),
104	3x	Value = c(swmr, natural_fertility, sit_fertility, residual_fertility)
105		)
106
107	3x	ans <- structure(res, components = comp_table, class = 'sit_competitiveness')
108	3x	return(ans)
109		}
110
111
112		#' @importFrom knitr kable
113		#' @export
114		print.sit_competitiveness <- function(x, digits = 2, ...) {
115	!	cli::cli_h1("Sterile male competitiveness")
116	!	cat("\n")
117	!	cli::cli_alert_info("Estimated value: {round(x, digits)}")
118	!	print(
119	!	knitr::kable(
120	!	attr(x, "components"),
121	!	format = "simple",
122	!	digits = digits
123		)
124		)
125	!	cat("\n")
126		}
127
128		#' Fried index
129		#'
130		#' Index of the mating competitiveness of sterile males relative to wild males.
131		#'
132		#' The competitiveness \eqn{\gamma} of the sterile male individuals is defined
133		#' as their relative capacity to mate with a wild female, compared to a wild
134		#' male.
135		#'
136		#' Thus, in a homogeneously mixed population with \eqn{M_s} sterile males and
137		#' \eqn{M_w} wild males, the probability that a mating occurs with a sterile
138		#' individual is
139		#' \deqn{
140		#' p_s = \frac{\gamma M_s}{M_w + \gamma M_s} =
141		#' \frac{\gamma R_{sw}}{1 + \gamma R_{sw}},
142		#' }
143		#' where \eqn{R_{sw} = M_s/M_w} is the sterile-wild male ratio.
144		#' At a given sterile-wild male ratio \eqn{R_{sw} = M_s/M_w} we observe a
145		#' fertility rate \eqn{H_s} in the field.
146		#'
147		#' Assuming a residual fertility rate \eqn{H_{rs}} for sterile males and a
148		#' natural fertility rate \eqn{H_w} for wild males, the observed fertility rate
149		#' \eqn{H_s} in the field is:
150		#' \deqn{
151		#' H_s = p_s H_{rs} + (1-p_s)H_w =
152		#' \frac{\gamma R_{sw}}{1 + \gamma R_{sw}} H_{rs} +
153		#' \frac{1}{1 + \gamma R_{sw}} H_w.
154		#' }
155
156		#' @param sterile_wild_mr Non-negative number. Ratio of sterile to wild males
157		#' in the target population.
158		#' @param natural_fertility Number between 0 and 1. Proportion of fertile eggs
159		#' in a natural (wild) population.
160		#' @param sit_fertility Number between 0 and 1. Proportion of fertile eggs
161		#' in the target population.
162		#' @param residual_fertility Number between 0 and 1. Proportion of fertile eggs
163		#' in a completely sterile population. Also called _residual_ fertility.
164		#'
165		#' @return Non-negative number.
166		#' @export
167		#'
168		#' @examples
169		#' fried_index(
170		#' sterile_wild_mr = 0.3,
171		#' natural_fertility = .9,
172		#' sit_fertility = .71,
173		#' residual_fertility = .05
174		#' )
175		fried_index <- function(
176		sterile_wild_mr,
177		natural_fertility,
178		sit_fertility,
179		residual_fertility
180		) {
181
182	8x	if (!any(sterile_wild_mr > 0)) {
183	1x	stop(
184	1x	"Cannot estimate competitiveness with a estimated wild population of 0."
185		)
186		}
187
188	7x	if (any(sit_fertility == residual_fertility, na.rm = TRUE)) {
189	1x	stop(
190	1x	"Cannot estimate competitiveness when the estimated fertility equals",
191	1x	" the sterile fertility."
192		)
193		}
194
195	6x	(natural_fertility - sit_fertility) /
196	6x	(sit_fertility - residual_fertility) /
197	6x	sterile_wild_mr
198
199		}
200
201		#' Fertility rate
202		#'
203		#' Observed fertility rate from egg surveys.
204		#'
205		#' The `sit` method uses all egg surveys in the study. For a `sit` object `x`,
206		#' `fertility_rate(x)` is thus equivalent to
207		#' `fertility_rate(sit_egg_surveys(x))`.
208		#'
209		#' @param x Either a `sit` or a `sit_egg_surveys` object.
210		#' @param species Character. Filter results for a given species or ignore the
211		#' species if missing.
212		#' @param pool Logical. If `TRUE` pools fertile and sterile counts from all traps
213		#' and survey dates. Otherwise (default) yields results by trap and survey
214		#' date.
215		#'
216		#' @return If `pool = TRUE`, a number. Otherwise, a data.frame with counts of
217		#' fertile and sterile eggs by trap and survey period, with the corresponding
218		#' proportion of fertile eggs.
219		#' @export
220		#' @examples
221		#' fertility_rate(sit_prototype, pool = TRUE)
222		fertility_rate <- function(x, species = NULL, pool = FALSE) {
223	9x	UseMethod('fertility_rate')
224		}
225
226		#' @export
227		fertility_rate.sit <- function(x, species, pool) {
228
229	2x	as_sit <- sit_egg_surveys(x)
230	2x	fertility_rate.sit_egg_surveys(as_sit, species, pool)
231
232		}
233
234		#' @importFrom stats reshape
235		#' @export
236		fertility_rate.sit_egg_surveys <- function(x, species, pool) {
237
238		## Aggregate counts by trap, survey date and fertility status
239	9x	x_long <- as.data.frame(x[, -1])
240
241	9x	x_wide <- stats::reshape(
242	9x	x_long,
243	9x	direction = "wide",
244	9x	idvar = c("trap_code", "datetime_start", "datetime_end"),
245	9x	timevar = "fertile",
246	9x	v.names = "n",
247	9x	varying = list(c("fertile", "sterile"))
248		)
249	9x	attr(x_wide, "reshapeWide") <- NULL
250
251	9x	if (pool) {
252	8x	tot_fertile <- sum(x_wide$fertile)
253	8x	total <- tot_fertile + sum(x_wide$sterile)
254	8x	return(tot_fertile / total)
255		}
256
257	1x	x_wide$fertility_rate = x_wide$fertile / (x_wide$fertile + x_wide$sterile)
258
259	1x	return(x_wide)
260
261		}
262
263
264		#' Sterile-Wild male ratio
265		#'
266		#' Observed ratio between sterile and wild male individuals from adult surveys.
267		#'
268		#' The `sit` method uses all adult surveys from the study (`sit`) area. For a
269		#' `sit` object `x`, `sterile_wild_male_ratio(x)` is thus equivalent to
270		#' `sterile_wild_male_ratio(sit_adult_surveys(x, area = 'sit'))`.
271		#'
272		#' @param x Either a `sit` or a `sit_adult_surveys` object.
273		#' @param pool Logical. If `TRUE` pools sterile and wild counts from all traps
274		#' and survey dates. Otherwise (default) yields results by trap and survey
275		#' date.
276		#'
277		#' @return If `pool = TRUE`, a number. Otherwise, a data.frame with counts of
278		#' sterile and wild males by trap and survey period, with the corresponding
279		#' ratio.
280		#' @export
281		#' @examples
282		#' sterile_wild_male_ratio(sit_prototype, pool = TRUE)
283		sterile_wild_male_ratio <- function(x, pool) {
284	8x	UseMethod('sterile_wild_male_ratio')
285		}
286
287		#' @export
288		sterile_wild_male_ratio.sit <- function(x, pool = FALSE) {
289
290	2x	as_sit <- sit_adult_surveys(x, area = "sit")
291	2x	sterile_wild_male_ratio.sit_adult_surveys(as_sit, pool)
292
293		}
294
295		#' @importFrom stats reshape
296		#' @export
297		sterile_wild_male_ratio.sit_adult_surveys <- function(x, pool = FALSE) {
298
299	8x	x_males <- x[x$sex == "male", ]
300
301		## Aggregate counts from all sterile populations together
302	8x	x_males$group <- ifelse(x_males$pop_col == 'wild', 'wild', 'sterile')
303	8x	x_group <- aggregate(
304	8x	x_males['n'],
305	8x	by = list(
306	8x	trap_code = x_males$trap_code,
307	8x	datetime_start = x_males$datetime_start,
308	8x	datetime_end = x_males$datetime_end,
309	8x	group = x_males$group
310		),
311	8x	FUN = "sum",
312	8x	na.rm = TRUE
313		)
314
315	8x	x_wide <- stats::reshape(
316	8x	x_group,
317	8x	direction = "wide",
318	8x	idvar = c("trap_code", "datetime_start", "datetime_end"),
319	8x	timevar = "group"
320		)
321	8x	x_wide <- setNames(x_wide, gsub("^n\\.", "", names(x_wide)))
322	8x	attr(x_wide, "reshapeWide") <- NULL
323
324		## Replace NAs due to missing records of some group at some trap and date by 0
325	8x	x_wide$sterile[is.na(x_wide$sterile)] <- 0
326	8x	x_wide$wild[is.na(x_wide$wild)] <- 0
327
328	8x	if (pool) {
329	6x	return(sum(x_wide$sterile) / sum(x_wide$wild))
330		}
331
332	2x	x_wide$swmr = x_wide$sterile / x_wide$wild
333
334	2x	return(x_wide)
335
336		}

1		#' Build a `sit` object
2		#'
3		#' Gathers data about traps, release events and field survey data into
4		#' a relational data model.
5		#'
6		#' @param traps A `sit_traps` object built with `sit_traps()`.
7		#' @param release_events A `sit_revents` object built with `sit_revents()`.
8		#' @param adult_surveys A `sit_adult_surveys` object built with `sit_adult_surveys()`.
9		#' @param egg_surveys A `sit_egg_surveys` object built with `sit_egg_surveys()`.
10		#'
11		#' @return A `sit` object that can be queried with functions from the package.
12		# @seealso [sit_sterile_male_competitivenes()], [sit_mdt()],
13		# [sit_flight_range()], [sit_diffusion()], [sit_survival()].
14		#' @import dm
15		#' @export
16		#'
17		#' @examples
18		#' sit(
19		#' sit_traps(fake_traps),
20		#' c(sit_revents(fake_rpoints), sit_revents(fake_rareal)),
21		#' sit_adult_surveys(fake_adults),
22		#' sit_egg_surveys(fake_eggs)
23		#' )
24		sit <- function(traps, release_events, adult_surveys = NULL, egg_surveys = NULL) {
25
26		# # Test example from Readme.
27		# traps = sit_traps(fake_traps)
28		# release_events = c(sit_revents(fake_rpoints), sit_revents(fake_rareal))
29		# adult_surveys = sit_adult_surveys(fake_adults)
30		# egg_surveys = sit_egg_surveys(fake_eggs)
31
32		## Extract embedded tables from attributes
33	20x	trap_types <- attr(traps, "trap_types")
34	20x	release_sites <- attr(release_events, "release_sites")
35	20x	attr(traps, "trap_types") <- NULL
36	20x	attr(release_events, "release_sites") <- NULL
37
38		## Remove the corresponding classes
39	20x	class(traps) <- setdiff(class(traps), 'sit_traps')
40	20x	class(release_events) <- setdiff(class(release_events), 'sit_revents')
41
42
43	20x	if (is.null(adult_surveys)) {
44		## If not specified, use a empty table with no observations
45	6x	adult_surveys <- sit_adult_surveys(x = sit::fake_adults[c(),])
46		}
47
48	20x	if (is.null(egg_surveys)) {
49		## If not specified, use a empty table with no observations
50	8x	egg_surveys <- sit_egg_surveys(x = sit::fake_eggs[c(),])
51		}
52
53		## Check whether variable 'trap_type' is required in adult data to uniquely
54		## identify traps
55	20x	non_unique_trap_codes <- any(
56	20x	duplicated(unique(st_drop_geometry(traps[c("code", "type_id")]))$code)
57		)
58	20x	if (non_unique_trap_codes) {
59	2x	multiple_adult_traps <- sum(trap_types$stage == 'adult') > 1
60	2x	if (multiple_adult_traps) {
61	2x	adults_has_type <- 'trap_type' %in% names(adult_surveys)
62	2x	if (!adults_has_type) {
63	1x	stop(
64	1x	glue(
65	1x	"Variable `trap_type` is required in adult data, ",
66	1x	"in order to uniquely identify the correct trap."
67		)
68		)
69		}
70		} else {
71		# Only one adult trap. Use it as trap type with message.
72	!	adult_surveys$type <- trap_types$label[trap_types$stage == "adult"]
73		}
74		# If traps are identified by code and type, need to match both
75		# variables using a merge.
76	1x	trap_labels <- merge(
77	1x	st_drop_geometry(traps)[, c("id", "code", "type_id")],
78	1x	trap_types[, c("id", "label")],
79	1x	by.x = "type_id",
80	1x	by.y = "id"
81		)
82
83	1x	a_trap_code.idx <- merge(
84	1x	setNames(
85	1x	adult_surveys[, c("id", "trap_code", "trap_type")],
86	1x	c("survey_id", "code", "label")
87		),
88	1x	trap_labels,
89	1x	by = c("code", "label"),
90	1x	sort = FALSE
91		)
92	1x	a_trap_code.idx <- a_trap_code.idx$id[order(a_trap_code.idx$survey_id)]
93		} else {
94		# Unique trap codes: simply match codes
95	18x	a_trap_code.idx <- match(adult_surveys$trap_code, traps$code)
96		}
97
98		## Check trap codes in adult data
99	19x	if (any(idx <- is.na(a_trap_code.idx))) {
100	1x	stop(
101	1x	glue(
102	1x	"Trap {vars_collapse(unique(adult_surveys$trap_code[idx]), 'code ', 'codes ')} ",
103	1x	"in adult field data not found in the imported traps."
104		)
105		)
106		}
107
108		## Check trap types in adult data
109	18x	stages <- trap_types$stage[traps$type_id[a_trap_code.idx]]
110	18x	if (any(idx <- stages != "adult")) {
111	1x	stop(
112	1x	glue(
113	1x	"Trap {vars_collapse(unique(adult_surveys$trap_code[idx]), 'code ', 'codes ')} ",
114	1x	"not of adult type. Please verify."
115		)
116		)
117		}
118
119		## Check population colours in adult data
120		## a_cols.idx matches colours and NAs for wild populations
121		## cols.idx is the reduced set of non-wild colours
122	17x	a_cols.idx <- match(adult_surveys$pop_col, release_events$colour)
123	17x	if (any(idx <- is.na(a_cols.idx) & adult_surveys$pop_col != "wild")) {
124	1x	stop(
125	1x	glue(
126	1x	"Release with ",
127	1x	"{vars_collapse(unique(adult_surveys$pop_col[idx]), 'colour ', 'colours ')} ",
128	1x	"not found. Please verify."
129		)
130		)
131		}
132
133		## Build foreign keys
134	16x	adult_surveys$trap_id <- traps$id[a_trap_code.idx]
135	16x	adult_surveys$population_id <- release_events$id[a_cols.idx]
136
137
138		## Select and reorder columns
139	16x	avars.idx <- c(
140	16x	grep("id", names(adult_surveys)),
141	16x	grep("_code\|_col\|id", names(adult_surveys), invert = TRUE)
142		)
143	16x	adult_surveys <- adult_surveys[, avars.idx]
144
145
146		## Check trap codes in egg data
147	16x	e_trap_code.idx <- match(egg_surveys$trap_code, traps$code)
148	16x	if (any(idx <- is.na(e_trap_code.idx))) {
149	1x	stop(
150	1x	glue(
151	1x	"Trap {vars_collapse(egg_surveys$trap_code[idx], 'code ', 'codes ')} ",
152	1x	"in egg field data not found in the imported traps."
153		)
154		)
155		}
156
157		## Check trap types in egg data
158	15x	stages <- trap_types$stage[traps$type_id[e_trap_code.idx]]
159	15x	if (any(idx <- stages != "egg")) {
160	1x	stop(
161	1x	glue(
162	1x	"Trap {vars_collapse(egg_surveys$trap_code[idx], 'code ', 'codes ')} ",
163	1x	"not of egg type. Please verify."
164		)
165		)
166		}
167
168		## Build foreign keys
169	14x	egg_surveys$trap_id <- traps$id[e_trap_code.idx]
170
171		## Select and reorder columns
172	14x	evars.idx <- c(
173	14x	grep("id", names(egg_surveys)),
174	14x	grep("_code\|id", names(egg_surveys), invert = TRUE)
175		)
176	14x	egg_surveys <- egg_surveys[, evars.idx]
177
178
179
180		## Build a relational data model with dm
181		## https://cynkra.github.io/dm/articles/howto-dm-df.html
182	14x	ans <- dm(
183	14x	adult_surveys,
184	14x	egg_surveys,
185	14x	release_events,
186	14x	release_sites,
187	14x	traps,
188	14x	trap_types
189		)
190
191		## Specify primary keys and check uniqueness
192	14x	ans <- dm_add_pk(ans, adult_surveys, "id", check = TRUE)
193	14x	ans <- dm_add_pk(ans, egg_surveys, "id", check = TRUE)
194	14x	ans <- dm_add_pk(ans, release_events, "id", check = TRUE)
195	14x	ans <- dm_add_pk(ans, release_sites, "id", check = TRUE)
196	14x	ans <- dm_add_pk(ans, traps, "id", check = TRUE)
197	14x	ans <- dm_add_pk(ans, trap_types, "id", check = TRUE)
198
199		## Specify foreign keys and check correspondence
200	14x	ans <- dm_add_fk(ans, egg_surveys, "trap_id", traps, check = FALSE)
201	14x	ans <- dm_add_fk(ans, adult_surveys, "trap_id", traps, check = FALSE)
202	14x	ans <- dm_add_fk(
203	14x	ans,
204	14x	adult_surveys,
205	14x	"population_id",
206	14x	release_events,
207		## Cannot check since there can be NA values corresponding to wild
208		## populations, which are not present in table release_events.
209	14x	check = FALSE
210		)
211	14x	ans <- dm_add_fk(ans, release_events, "site_id", release_sites, check = FALSE)
212
213		## Avoid checking for consistency. Problematic since traps is `sf`
214		## https://github.com/cynkra/dm/issues/705
215	14x	ans <- dm_add_fk(ans, traps, "type_id", trap_types, check = FALSE)
216
217
218	14x	ans <- dm_set_colors(
219	14x	ans,
220	14x	peru = traps,
221	14x	darkmagenta = release_events,
222	14x	darkolivegreen = c("egg_surveys", "adult_surveys")
223		)
224	14x	class(ans) <- c("sit", class(ans))
225
226		# dm_draw(ans, view_type = "all")
227
228	14x	ans <- validate_sit(ans)
229
230	14x	return(ans)
231		}
232
233
234		#' @keywords internal
235		#' @export
236		as_dm.sit <- function(...) {
237		## Remove sit class
238	1x	x <- list(...)[[1]]
239	1x	class(x) <- setdiff(class(x), "sit")
240	1x	return(x)
241		}
242
243		#' @import cli
244		#' @export
245		print.sit <- function(x, ...) {
246
247	1x	cli::cat_rule("sit", col = "violet")
248
249	1x	revents <- as.data.frame(sit_revents(x))
250	1x	revents_by_site <- table(revents$site_id)
251	1x	revents_by_type <- table(revents$type)
252
253	1x	revents_site_chr <- paste0(
254	1x	sum(revents_by_site),
255	1x	" (", paste(revents_by_site, collapse = " + "), ")"
256		)
257
258		# cli::cat_line("Release events: ", revents_site_chr)
259	1x	cli::cat_line("Release events: ", describe_numbers(revents_by_type))
260	1x	cli::cat_line("Release sites : ", length(revents_by_site))
261
262	1x	traps <- st_drop_geometry(as.data.frame(sit_traps(x)))
263	1x	traps_type <- table(traps$type_label)
264	1x	traps_stage_area <- table(traps$type_stage, traps$area)
265
266	1x	cli::cat_line("Adult traps : ", describe_numbers(traps_stage_area["adult", ]))
267	1x	cli::cat_line("Egg traps : ", describe_numbers(traps_stage_area["egg", ]))
268	1x	cli::cat_line("Trap types : ", describe_numbers(traps_type))
269
270	1x	surveys <- c(adult = nrow(sit_adult_surveys(x)), egg = nrow(sit_egg_surveys(x)))
271
272	1x	cli::cat_line("Surveys : ", describe_numbers(surveys))
273		}
274
275
276		#' @import cli
277		#' @importFrom stats setNames
278		#' @export
279		summary.sit <- function(object, ...) {
280
281	1x	cli::cat_rule("Releases", col = "violet")
282	1x	print(summary(sit_revents(object)))
283	1x	cli::cat_line()
284
285	1x	cli::cat_rule("Traps", col = "orange")
286	1x	print(summary(sit_traps(object)))
287	1x	cli::cat_line()
288
289	1x	cli::cat_rule("Adult surveys (total captures)", col = "green")
290
291	1x	as <- sit_adult_surveys(object)
292	1x	as$pop_col[is.na(as$pop_col)] <- 'wild'
293	1x	as$pop_col <- factor(as$pop_col, c(sit_revents(object)$colour, 'wild'))
294		## Add trap area
295	1x	as <- merge(
296	1x	as,
297	1x	setNames(
298	1x	st_drop_geometry(
299	1x	as.data.frame(sit_traps(object))
300	1x	)[, c('code', 'area')],
301	1x	c('trap_code', 'area')
302		)
303		,
304	1x	by = 'trap_code'
305		)
306	1x	print(
307	1x	tapply(
308	1x	seq_len(nrow(as)),
309	1x	list(as$pop_col, as$sex, as$area),
310	1x	function(.) sum(as$n[.])
311		)
312		)
313	1x	cli::cat_line()
314
315	1x	cli::cat_rule("Egg surveys (total captures)", col = "green")
316	1x	es <- sit_egg_surveys(object)
317	1x	es$fertile <- factor(es$fertile, labels = c('Sterile', 'Fertile'))
318		## Add trap area
319	1x	es <- merge(
320	1x	es,
321	1x	setNames(
322	1x	st_drop_geometry(
323	1x	as.data.frame(sit_traps(object))
324	1x	)[, c('code', 'area')],
325	1x	c('trap_code', 'area')
326		)
327		,
328	1x	by = 'trap_code'
329		)
330	1x	print(
331	1x	tapply(
332	1x	seq_len(nrow(es)),
333	1x	list(es$area, es$fertile),
334	1x	function(.) sum(es[., "n"])
335		)
336		)
337	1x	cli::cat_line()
338
339	1x	invisible(object)
340		}
341
342
343		#' @export
344		plot.sit <- function(x, y, ...) {
345
346		## Only plot areas and trap types for which there are coordinates
347	!	summary_traps <- summary(sit_traps(x))
348	!	target_traps <- summary_traps[summary_traps$Coords, ]
349
350		## Target only the sit area
351	!	target_traps <- target_traps[target_traps$Area == "sit"]
352
353		## Plot release points and sit-area traps by type
354	!	traps <- sit_traps(x, area = "sit", type = target_traps$Type)
355		## Merge with trap_types attribute to gain access to type_label
356		## Despite the name, returns a sf object
357	!	traps <- as.data.frame(traps)
358
359	!	rpoints <- st_as_sf(sit_revents(x, type = "point"))
360
361	!	plot(traps["type_label"], pch = 20, col = "darkgray", reset = FALSE,
362	!	main = "Traps and release points in the sit area")
363	!	plot(rpoints["colour"], add = TRUE, pch = 19, col = rpoints$colour)
364
365		}

1		#' Collapse variables
2		#'
3		#' Collapse a vector of variables names into a properly formatted string that
4		#' can be used in messages.
5		#'
6		#' @param x Vector coerced to character. Variable names to be concatenated.
7		#' @param intro1 Character. String included at the beginning if `x` is of
8		#' length 1.
9		#' @param intros Character. String included at the beginning if `x` is
10		#' of length greater than 1.
11		#' @param quote Character. Symbol used to quote variable names.
12		#' @param sep Character. Symbol used to separate variable names.
13		#' @param last Character. String used to separate the last two items if `x` has
14		#' at least 2 items.
15		#'
16		#' If `x` is an empty character vector, returns another empty character vector.
17		#'
18		#' @return
19		#' Character string.
20		#' @import glue
21		#' @keywords internal
22		#'
23		#' @examples
24		#' sit:::vars_collapse(paste0("var", 1:4))
25		#' sit:::vars_collapse("single_var")
26		#' sit:::vars_collapse("single_var", intro1 = "", quote = "")
27		vars_collapse <- function(
28		x,
29		intro1 = "Variable ",
30		intros = "Variables ",
31		quote = "'",
32		sep = ", ",
33		last = " and "
34		) {
35
36	1x	if (length(x) == 0L) return(x)
37
38	24x	vars_str <- glue_collapse(
39	24x	paste0(quote, x, quote), sep = sep, last = last
40		)
41	24x	ans <- glue(
42	24x	"{variables}{vars_str}",
43	24x	variables = ifelse(length(x) > 1, intros, intro1)
44		)
45	24x	return(unclass(ans))
46		}
47
48		#' Are given variable names present in a data.frame
49		#'
50		#' Verify that a vector of variable names are all actual variables in
51		#' a `data.frame`.
52		#'
53		#' @param x `data.frame`-like object
54		#' @param names_given Character vector of variable names.
55		#'
56		#' @return NULL if all variables exist in the `data.frame`. Otherwise, stops
57		#' with a informative error message.
58		#' @import glue
59		#' @keywords internal
60		#'
61		#' @examples
62		#' sit:::check_variable_names(data.frame(a = 1), "a")
63		check_variable_names <- function(x, names_given) {
64	156x	if (! all(matching <- names_given %in% names(x))) {
65	4x	vars <- vars_collapse(names_given[!matching])
66	4x	stop(glue("{vars} not found in import data. Please verify."))
67		}
68		}
69
70
71		#' Build the internal representation of a data set
72		#'
73		#' This function selects the relevant variables from the data set
74		#' and renames them according to a given set of internal names.
75		#'
76		#' @param x `data.frame`.
77		#' @param varnames Character vector. Variable names from `x`.
78		#' @param internal_names Character vector. Internal names.
79		#'
80		#' @return A `data.frame` with the internal variable names in the expected
81		#' order.
82		#' @keywords internal
83		#'
84		#' @examples
85		#' xx <- data.frame(a = 1, b = 2, c = 3)
86		#' sit:::build_internal_df(xx, c("a", "c"), c("A", "C"))
87		build_internal_df <- function(x, varnames, internal_names) {
88
89		## Check that `varnames` exist in `x`
90	149x	check_variable_names(x, varnames)
91
92		## Remove unused columns
93	148x	ans <- x[, varnames]
94
95		## Rename variables using internal names
96		## The subsetting in match is necessary in the case of a `sf` object
97		## which has an additional column `geometry`.
98	148x	names(ans)[match(varnames, names(ans))] <- internal_names
99
100	148x	return(ans)
101		}
102
103
104		#' Test if a number is whole
105		#'
106		#' @param x numeric.
107		#' @param tol numeric. Tolerance.
108		#'
109		#' @return Logical.
110		#' @keywords internal
111		#'
112		#' @examples
113		#' sit:::is.wholenumber(5)
114		#' sit:::is.wholenumber(5.00002)
115		is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) {
116	4x	abs(x - round(x)) < tol
117		}
118
119
120		#' Combine rows
121		#'
122		#' Remove duplicate rows (not taking into account id) and reset the id.
123		#'
124		#' @param x a list of `data.frame`s. All must share the same variables.
125		#' @param id_col Numeric. Column number of the identification variable.
126		#'
127		#' @return An object of the same class of the first element in `x`.
128		#' @keywords internal
129		#'
130		#' @examples
131		#' dat1 <- data.frame(id = 1:3, x = letters[1:3])
132		#' dat2 <- data.frame(id = 1:2, x = letters[2:3])
133		#' sit:::combine_rows(list(dat1, dat2))
134		combine_rows <- function(x, id_col = 1, recoder = FALSE) {
135
136		## Bind rows
137	37x	xbind <- do.call('rbind', x)
138
139		## Remove duplicates, disregarding the id
140	37x	ans <- unique(xbind[-id_col])
141
142		## Reset id
143	37x	ans$id <- seq.int(nrow(ans))
144
145		## Reorder columns
146	37x	ans <- ans[, c(ncol(ans), seq.int(ncol(ans) - 1))]
147
148	37x	if (recoder) {
149		## Recoding map
150	3x	rec_l <- lapply(
151	3x	x, function(.)
152	3x	merge(., ans, by = intersect(names(.)[-id_col], names(ans)[-1]))
153		)
154
155		## Recoding function: recode a vector x from item i of the initial list.
156	3x	rec_f <- function(x, i) {
157	6x	rec_l[[i]]$id.y[match(x, rec_l[[i]]$id.x)]
158		}
159
160	3x	attr(ans, 'recoder') <- rec_f
161		}
162
163	37x	return(ans)
164		}
165
166
167		#' Compute survey periods
168		#'
169		#' Given the survey date and either the activation or duration, compute the
170		#' starting and ending survey dates and times.
171		#'
172		#' The function checks the precedence of the resulting starting and ending
173		#' times.
174		#'
175		#' @param x A `data.frame` object.
176		#' @param input A Character string. Either 'activation' or 'duration'.
177		#'
178		#' @return The same data.frame with the input variables replaced by
179		#' datetime_start and datetime_end.
180		#' @keywords internal
181		#'
182		#' @examples
183		#' td <- data.frame(
184		#' id = 1:2, survey = c('2021-11-15', '2021-11-01'),
185		#' activation = c('2021-11-10', '2021-10-30'), duration = c(7, 5)
186		#' )
187		#'
188		#' ## Note the difference between the two inputs due to daylight saving time
189		#' ## (which depends on the time zone).
190		#' sit:::compute_periods(td, input = "activation")
191		#' sit:::compute_periods(td, input = "duration")
192		compute_periods <- function(x, input) {
193
194	55x	if (input == 'duration' && !all(idx <- x$duration > 0)) {
195	1x	stop("Argument `duration` must be non-negative.")
196		}
197
198	53x	x$datetime_end <- as.POSIXct(x$survey)
199	53x	if (input == 'activation') {
200	3x	x$datetime_start <- as.POSIXct(x$activation)
201		} else {
202	50x	stopifnot(input == 'duration')
203	50x	x$datetime_start <- x$datetime_end -
204	50x	as.difftime(x$duration, units = "days")
205		}
206
207	53x	if (!all(x$datetime_end - x$datetime_start > 0)) {
208	1x	stop("Dates of `activation` must precede survey dates.")
209		}
210
211	52x	x$survey <- x$duration <- x$activation <- NULL
212
213	52x	return(x)
214		}
215
216
217		#' Colours for concepts in the data model.
218		#'
219		#' @param x Character vector. Concepts among `traps`, `release_events` or
220		#' `surveys`. If missing, all three are returned in that order.
221		#'
222		#' @return A named vector of colours.
223		#' @keywords internal
224		#'
225		#' @examples
226		#' sit:::sit_colours()
227		#' sit:::sit_colours('release_events')
228		#' sit:::sit_colours(c('surveys', 'traps', 'surveys'))
229		sit_colours <- function(x) {
230	3x	cols <- c(
231	3x	traps = "peru", release_events = "darkmagenta", surveys = "darkolivegreen"
232		)
233
234	1x	if (missing(x)) x <- names(cols)
235
236	3x	return(cols[x])
237		}
238
239		#' Assign a S3 class
240		#'
241		#' Append a class to the beginning of an object's list of classes.
242		#'
243		#' If the object already inherited from the class, it is moved to the beginning
244		#' of the class list. Otherwise, the new class is added.
245		#'
246		#' @param x Any object.
247		#' @param classname Character string. The class name to be assigned.
248		#'
249		#' @keywords internal
250		#' @return The same object `x` with the updated `class` attribute.
251		#'
252		#' @examples
253		#' sit:::set_class(7, 'prime')
254		#' sit:::set_class( sit:::set_class(7, 'prime'), 'numeric')
255		set_class <- function(x, classname) {
256
257	168x	class(x) <- c(classname, setdiff(class(x), classname))
258	168x	return (x)
259		}
260
261
262		#' Print a comma-separated vector of terms
263		#'
264		#' Borrowed from [dm].
265		#'
266		#' @param x Character vector.
267		#' @param max_commas Integer. Maximum number of terms to print.
268		#' @param capped Logical. If TRUE do not show the total number of elements for
269		#' long lists
270		#' @param fun One-sided formula. Function to apply to the vector. Must return
271		#' another vector.
272		#'
273		#' @keywords internal
274		#' @import cli
275		#' @import rlang
276		#' @examples
277		#' sit:::commas(letters[1:6])
278		#' sit:::commas(letters[1:9])
279		#' sit:::commas(letters[1:9], capped = TRUE)
280		#' sit:::commas(letters[1:9], fun = toupper)
281		commas <- function(x, max_commas = 6L, capped = FALSE, fun = identity) {
282
283	11x	fun <- rlang::as_function(fun)
284
285	11x	if (rlang::is_empty(x)) {
286	1x	x <- ""
287	10x	} else if (length(x) > max_commas) {
288	6x	cap <- if (!capped) paste0(" (", length(x), " total)")
289
290	6x	x <- c(
291	6x	fun(x[seq_len(max_commas - 1)]),
292	6x	paste0(cli::symbol$ellipsis, cap)
293		)
294		} else {
295	4x	x <- fun(x)
296		}
297
298	11x	glue_collapse(x, sep = ", ")
299		}
300
301
302		#' Description of numeric vectors
303		#'
304		#' Textual description of numeric vectors as a sum of components for descriptive
305		#' purposes.
306		#'
307		#' If the vector is named, use the names in parentheses besides the
308		#' corresponding number.
309		#'
310		#' @param x Numeric vector. Optionally named.
311		#' @param collapse Character string. Symbol separating numbers.
312		#' @param zero.print Logical. If TRUE, omit zero-valued entries.
313		#'
314		#' @keywords internal
315		#' @examples
316		#' sit:::describe_numbers(1:4)
317		#' sit:::describe_numbers(1:4, collapse = ", ")
318		#' sit:::describe_numbers(setNames(1:3, letters[1:3]))
319		describe_numbers <- function(x, collapse = " + ", zero.print = FALSE) {
320
321	2x	if (missing(x)) x <- numeric()
322
323	12x	if (!is.null(names(x))) {
324	6x	xn <- paste0("(", names(x), ")")
325	6x	ans <- paste(x, xn)
326		} else {
327	6x	ans <- as.character(x)
328		}
329
330	12x	if (!zero.print) {
331	12x	ans <- ans[x > 0]
332		}
333
334	12x	ans <- paste(ans, collapse = collapse)
335
336	12x	return(ans)
337		}
338
339
340		#' Are times within a period
341		#'
342		#' Test whether some date-times fall within a set of time periods.
343		#'
344		#' For each date-time in `x` returns `TRUE` if it falls within __any__ of the
345		#' temporal intervals defined by the corresponding elements in `starts` and
346		#' `ends`.
347		#'
348		#' The interval is open by the left and closed by the right. I.e., it includes
349		#' the right boundary but not the left.
350		#'
351		#' If `starts` and `end` are empty vectors (e.g. `as.POSIXct(character())`),
352		#' return `TRUE` for each element in `x`.
353		#'
354		#' @param x Vector of POSIXct date-times.
355		#' @param starts Vector of POSIXct date-time starting times.
356		#' @param periods Vector of POSIXct date-time ending times.
357		#'
358		#' @return Logical vector of the same size as `x`.
359		#' @keywords internal
360		#'
361		#' @examples
362		#' test_times <- as.POSIXct(
363		#' c('2021-05-04', '2021-06-28', '2012-10-05', '2021-10-05')
364		#' )
365		#' spring_2021 <- as.POSIXct(c('2021-03-21', '2021-06-20'))
366		#' fall_2021 <- as.POSIXct(c('2021-09-21', '2021-12-20'))
367		#' starts <- c(spring_2021[1], fall_2021[1])
368		#' ends <- c(spring_2021[2], fall_2021[2])
369		#' sit:::is_date_within(test_times, starts, ends)
370		#'
371		#' # Left-open and right-closed
372		#' day <- as.POSIXct('2021-05-04')
373		#' next_day <- as.POSIXct('2021-05-05')
374		#' sit:::is_date_within(day, day, next_day) # FALSE
375		#' sit:::is_date_within(next_day, day, next_day) # TRUE
376		is_date_within <- function(x, starts, ends) {
377
378	101x	stopifnot(
379	101x	inherits(x, 'POSIXct'),
380	101x	inherits(starts, 'POSIXct'),
381	101x	inherits(ends, 'POSIXct'),
382	101x	(n_periods <- length(starts)) == length(ends)
383		)
384
385	10x	if (length(x) == 0L) return(logical())
386	67x	if (length(starts) == 0L) return(rep(TRUE, length(x)))
387
388	24x	stopifnot(all(starts <= ends))
389
390		## Matrix of size n_periods x n_times
391	24x	conditions <-
392	24x	matrix(
393	24x	vapply(x, `>`, rep(TRUE, length(starts)), starts) &
394	24x	vapply(x, `<=`, rep(TRUE, length(ends)), ends),
395	24x	n_periods, length(x)
396		)
397
398	24x	return(apply(conditions, 2, any))
399		}
400
401
402		#' Release types
403		#'
404		#' @param x `sf` of `sfc` object, representing the table of release sites.
405		#' @return Logical vector with as many elements as rows in `x`.
406		#' @keywords internal
407		#' @examples
408		#' sit:::release_types(sit_prototype[['release_sites']])
409		release_types <- function(x) {
410	108x	ifelse(st_is_empty(x), 'areal', 'point')
411		}
412
413
414		#' Get time periods following releases
415		#'
416		#' Retrieve the time periods from the release dates to a given number of
417		#' days after the corresponding release dates.
418		#'
419		#' If both arguments are missing, returns empty vectors.
420		#' Otherwise, specifying `following_releases` with missing `following_days`
421		#' results in an error, whereas the converse is ignored with a warning.
422		#'
423		#' @param following_releases A `sit_revents` object.
424		#' @param following_days Numeric. Number of days targeted. Can be fractional.
425		#'
426		#' @return A list with vectors `starts` and `ends` of `POSIXct` dates.
427		#' @keywords internal
428		#'
429		#' @examples
430		#' sit:::get_periods(sit_revents(sit_prototype), 7.5)
431		get_periods <- function(following_releases, following_days) {
432	92x	if (!missing(following_releases)) {
433	16x	if (missing(following_days)) {
434	1x	stop("Use of argument following_releases requires ",
435	1x	"following_days, which is currently missing.")
436		}
437	15x	stopifnot(inherits(following_releases, 'sit_revents'))
438	15x	stopifnot(is.numeric(following_days), all(following_days > 0))
439
440	15x	starts <- as.POSIXct(following_releases$date)
441	15x	ends <- starts + as.difftime(following_days, units = "days")
442		} else {
443		## missing following_releases
444	76x	if (!missing(following_days)) {
445	1x	warning("Argument `following_days` ignored.")
446		}
447		## Using empty periods in is_date_within returns TRUE
448	76x	starts <- ends <- as.POSIXct(character())
449		}
450
451	91x	return(list(starts = starts, ends = ends))
452		}
453
454
455		#' Maximum observed age
456		#'
457		#' Maximum age observed in a set of captures.
458		#'
459		#' @param captures Number of observed individuals.
460		#' @param ages Ages of the corresponding captures.
461		#'
462		#' @return Number.
463		#' @keywords internal
464		#'
465		#' @examples
466		#' sit:::max_observed_age(captures = rep(1, 5), ages = c(1, 2, 8, 3, 0))
467		#' sit:::max_observed_age(captures = c(1, 1, 0, 1, 1), ages = c(1, 2, 8, 3, 0))
468		max_observed_age <- function(captures, ages) {
469
470	15x	stopifnot(
471	15x	length(captures) == length(ages),
472	15x	is.numeric(captures),
473	15x	is.numeric(ages)
474		)
475
476	1x	if (!length(captures) > 0) return(NA)
477
478	13x	counts_by_age <-
479	13x	tapply(captures, ages, FUN = sum, na.rm = TRUE)
480
481	13x	max_age.idx <- max(which(counts_by_age > 0))
482	13x	ans <- as.numeric(names(counts_by_age)[max_age.idx])
483
484	13x	return(ans)
485		}
486
487
488
489		#' Print a table of variables and their description
490		#'
491		#' @param ... A even number of character strings. Variable names followed by
492		#' their description. Variables are printed in mono-spaced font.
493		#' @param header Logical. Include the header "Variable" "Description".
494		#' @param caption Character string. Table caption.
495		#'
496		#' @return A kable object.
497		#' @keywords internal
498		#' @export
499		#'
500		#' @examples
501		#' print_variable_table("Var1", "Desc1", "Var2", "Desc2")
502		#' print_variable_table(
503		#' "Var1", "Desc1",
504		#' "Var2", "Desc2",
505		#' header = FALSE
506		#' )
507		print_variable_table <- function(
508		..., header = TRUE, caption = NULL) {
509
510	!	ans <- as.data.frame(matrix(c(...), byrow = TRUE, ncol = 2))
511
512		## Use monospace font for variables
513	!	ans[, 1] <- paste0("`", ans[, 1], "`")
514
515		## Set header
516	!	ans <- setNames(
517	!	ans,
518	!	switch(header, c("Variable", "Description"), NULL)
519		)
520
521	!	kableExtra::kbl(
522	!	ans,
523	!	booktabs = TRUE,
524	!	toprule = '',
525	!	midrule = '',
526	!	bottomrule = '',
527	!	vline = '',
528	!	escape = FALSE,
529	!	centering = FALSE
530		) %>%
531	!	kableExtra::kable_styling(full_width = FALSE, position = "left") %>%
532	!	kableExtra::column_spec(2, width = "4in")
533
534		}

1		#' Import or retrieve adult surveys
2		#'
3		#' Import survey field data from adult traps, or retrieve the data from a `sit`
4		#' object.
5		#'
6		#' The argument `type` is mandatory only if needed to uniquely identify the trap
7		#' (i.e., more than one adult trap type is declared with [sit_trap_types()] and
8		#' the same code has been used for traps of different types). If the user used
9		#' unique codes for each trap, which will be checked internally, this is
10		#' redundant and will be ignored in any case.
11		#'
12		#' __Dates__ must be provided in the [RFC 3339
13		#' format](https://www.ietf.org/rfc/rfc3339.txt) (a variation of the [ISO 8601
14		#' format](https://en.wikipedia.org/wiki/ISO_8601)), i.e., `2021-12-31` or
15		#' `2019-11-23 15:00`.
16		#'
17		#' @param x A `data.frame` with field survey data in
18		#' [tidy](https://www.openscapes.org/blog/2020/10/12/tidy-data/) format.
19		#' @param trap Identification code of the trap. This is the user’s code
20		#' corresponding to the `id` argument in [sit_traps()].
21		#' @param type Character string. Trap type.
22		#' @param survey Character. Date and optionally the time of the survey. See
23		#' details on the date format.
24		#' @param activation Character. Date and optionally the time of the survey
25		#' activation. See details on the date format. Either `activation` or
26		#' `duration` are required.
27		#' @param duration Numeric. Number of days (integer or fractional) that the trap
28		#' has been functioning.
29		#' @param population Character. Either `wild` or any of the colour codes in the
30		#' release events imported with [sit_revents()].
31		#' @param species Character. Only required if the counts of wild individuals are
32		#' further categorised by species. Make sure to use the same nomenclature as
33		#' in [sit_revents()]. In the extraction method, a character vector of species
34		#' to be returned. Defaults to `NULL`, which ignores the species variable.
35		#' @param sex Character. Either `male` or `female`.
36		#' @param n Mandatory. An integer number of the surveyed individuals in the
37		#' corresponding group.
38		#' @param ... Used to pass arguments to specific method.
39		#'
40		#' @return A object of class `sit_adult_surveys` which can be used in [sit()].
41		#' @importFrom utils head tail
42		#' @export
43		#'
44		#' @family importing
45		#' @seealso [sit()]
46		#'
47		#' @examples
48		#' ad_surv <- data.frame(
49		#' trap = 1:3,
50		#' survey = c("2021-04-01", "2021-08-21", "2021-08-21"),
51		#' duration = rep(7L, 3),
52		#' population = c("blue", "yellow", "yellow"),
53		#' species = "aeg",
54		#' sex = c("male", "male", "female"),
55		#' n = 1:3
56		#' )
57		#' sit_adult_surveys(ad_surv)
58		sit_adult_surveys <- function(x, ...) {
59	88x	UseMethod('sit_adult_surveys')
60		}
61
62
63		#' @describeIn sit_adult_surveys Imports field survey data from adult traps.
64		#' @export
65		sit_adult_surveys.data.frame <- function(
66		x,
67		trap = "trap",
68		type = "type",
69		survey = "survey",
70		activation = "activation",
71		duration = "duration",
72		population = "population",
73		species = "species",
74		sex = "sex",
75		n = "n",
76		...
77		) {
78		# # Debug
79		# args <- names(formals(sit_adult_surveys.data.frame)[c(-1, -11)])
80		# for (arg in args) assign(arg, arg)
81
82		## The requirement of the `type` argument can only be verified when building
83		## the `sit` object, by looking at the codes and types of the available traps.
84
85		## Check here whether there are variables in the data set corresponding
86		## to optional arguments and set them to NULL if not.
87	29x	opt_args <- c(type, activation, duration, species)
88	29x	for (arg in opt_args[!opt_args %in% names(x)]) {
89	57x	assign(arg, NULL)
90		}
91
92		## Exactly one of `activation` or `duration` are required
93	29x	if ((l_ad <- length(c(activation, duration))) != 1L) {
94	2x	if (l_ad == 0L) {
95	1x	stop("Either activation or duration are required arguments. ",
96	1x	"Both are currently missing.")
97		}
98	1x	if (l_ad == 2L) {
99	1x	stop("Only one of `activation` or `duration` needed. ",
100	1x	"Please remove one.")
101		}
102		}
103
104		## Build the internal representation
105		## The `population_id` can only be completed once merging with
106		## a `sit_revents` into a `sit` object.
107	27x	inames <- c(
108	27x	"trap_code", "trap_type"[!is.null(type)], "survey",
109	27x	"activation"[!is.null(activation)],
110	27x	"duration"[!is.null(duration)], "pop_col",
111	27x	"species"[!is.null(species)], "sex", "n"
112		)
113		# Variable names reported through non-NULL arguments
114		# vnames <- unlist(lapply(tail(head(names(formals()), -1), -1), get,
115		# pos = environment()))
116		# Alternatively:
117	27x	vnames <- unname(
118	27x	unlist(
119	27x	as.list(
120	27x	environment())[tail(head(names(formals()), -1), -1)]
121		)
122		)
123
124	27x	ans <- build_internal_df(
125	27x	x[!is.na(x$n), ], # remove surveys with missing counts
126	27x	varnames = vnames,
127	27x	internal_names = inames
128		)
129
130		## Check duplicated surveys
131	27x	if (any(idx <- duplicated(ans))) {
132	1x	stop(
133	1x	glue(
134	1x	"{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
135	1x	"duplicated. Please fix."
136		)
137		)
138		}
139
140		## Compute starting and ending dates and times
141	26x	ans <- compute_periods(
142	26x	ans,
143	26x	input = ifelse(is.null(activation), "duration", "activation")
144		)
145
146		## Survey id
147	26x	ans$id <- if (nrow(ans) > 0) {
148	20x	seq.int(nrow(ans))
149		} else {
150	6x	integer()
151		}
152
153		## Reorder columns
154	26x	ans <- ans[
155	26x	, c("id",
156	26x	"trap_code", # Changed into `trap_id` within `sit()`
157	26x	"trap_type"[!is.null(type)], # Possibly needed for `trap_id`
158	26x	"pop_col", # Changed into `population_id` within `sit()`
159	26x	"datetime_start",
160	26x	"datetime_end",
161	26x	"species"[!is.null(species)],
162	26x	"sex",
163	26x	"n"
164		)
165		]
166
167	26x	class(ans) <- c("sit_adult_surveys", class(ans))
168	26x	return(ans)
169		}
170
171
172		#' @param area Character vector. Either `control`, `sit` or both (default).
173		#' @param trap_type Character vector. Any subset of
174		#' `unique(sit_trap_types(x)$label)`.
175		#' @param following_releases A `sit_release` object with a subset of release
176		#' events or missing (default) for all release events. Use in combination with
177		#' `following_days` to filter surveys of the target populations within a given
178		#' number of days after release. Note that counts of wild populations will
179		#' always be included in the results with a value of `NA` in `pop_col`.
180		#' @param following_days Integer or missing (default). Number of days after
181		#' releases to return, if `following_releases` is not missing.
182		#' @describeIn sit_adult_surveys Retrieve adult survey data, possibly filtering
183		#' results.
184		#' @importFrom rlang .data .env
185		#' @importFrom tidyr everything
186		#' @export
187		sit_adult_surveys.sit <- function(
188		x,
189		area = c('control', 'sit'),
190		trap_type = sit_trap_types(x)$label,
191		following_releases,
192		following_days,
193		species = NULL,
194		...) {
195
196	59x	target_periods <- get_periods(following_releases, following_days)
197
198		## Filter trap areas, trap types, and dates
199	59x	x_sub <-
200	59x	x %>%
201	59x	dm_filter(traps = (.data$area %in% .env$area)) %>%
202	59x	dm_filter(trap_types = (.data$label %in% .env$trap_type)) %>%
203	59x	dm_filter(
204	59x	adult_surveys =
205	59x	is_date_within(
206	59x	.data$datetime_end,
207	59x	.env$target_periods$starts,
208	59x	.env$target_periods$ends
209		)
210		)
211
212		# ## Filters for population if required
213		# ## Can't use dm_filter, as it removes surveys of wild population
214		# ## in adult_surveys. Filter the final object manually below.
215		# ## https://github.com/cynkra/dm/issues/706
216		# if (!missing(following_releases)) {
217		# x_sub <- dm_filter(
218		# x_sub,
219		# 'release_events',
220		# .data$colour %in% following_releases$colour
221		# )
222		#
223		# }
224
225		## Filters for species if required
226	59x	if (!is.null(species)) {
227	5x	if (!'species' %in% names(x[['adult_surveys']])) {
228	1x	stop(
229	1x	"There is no variable `species` in the adult surveys."
230		)
231		}
232	4x	x_sub <- dm_filter(
233	4x	x_sub,
234	4x	adult_surveys = (.data$species %in% .env$species)
235		)
236		}
237
238	58x	ans <- x_sub %>%
239		# dm_get_tables() %>% `[`(c("traps", "trap_types"))
240		## Select and rename columns from relevant tables
241	58x	dm_select('release_events', .data$id, pop_col = .data$colour) %>%
242	58x	dm_select('traps', .data$id, trap_code = .data$code) %>%
243	58x	dm_flatten_to_tbl('adult_surveys') %>%
244	58x	select(-.data$trap_id, -.data$population_id, -.data$geometry) %>%
245	58x	select(.data$id, .data$trap_code, .data$pop_col, everything()) %>%
246	58x	set_class('sit_adult_surveys')
247
248		## Return wild population as 'wild' rather than NA
249	58x	ans$pop_col[is.na(ans$pop_col)] <- 'wild'
250
251		## Filters for population if required. Always keep wild populations.
252	58x	if (!missing(following_releases)) {
253	9x	idx <- ans$pop_col %in% c(following_releases$colour, 'wild')
254	9x	ans <- ans[idx, ]
255		}
256
257
258	58x	return(ans)
259		}
260

1		#' Import or retrieve egg surveys
2		#'
3		#' Import survey field data from egg traps, or retrieve the data from a `sit`
4		#' object.
5		#'
6		#' The argument `type` is mandatory only if needed to uniquely identify the trap
7		#' (i.e., more than one egg trap type is declared with [sit_trap_types()] and
8		#' the same code has been used for traps of different types). If the user used
9		#' unique codes for each trap, which will be checked internally, this is
10		#' redundant and will be ignored in any case.
11		#'
12		#' __Dates__ must be provided in the [RFC 3339
13		#' format](https://www.ietf.org/rfc/rfc3339.txt) (a variation of the [ISO 8601
14		#' format](https://en.wikipedia.org/wiki/ISO_8601)), i.e., `2021-12-31` or
15		#' `2019-11-23 15:00`.
16		#'
17		#' @param x A `data.frame` with field survey data in
18		#' [tidy](https://www.openscapes.org/blog/2020/10/12/tidy-data/) format.
19		#' @param trap Identification code of the trap. This is the user’s code
20		#' corresponding to the `id` argument in [sit_traps()].
21		#' @param type Character string. Trap type.
22		#' @param survey Character. Date and optionally the time of the survey. See
23		#' details on the date format.
24		#' @param activation Character. Date and optionally the time of the survey
25		#' activation. See details on the date format. Either `activation` or
26		#' `duration` are required.
27		#' @param duration Numeric. Number of days (integer or fractional) that the trap
28		#' has been functioning.
29		#' @param species Character.
30		#' @param fertile Logical. Whether `n` represents counts of __fertile__ eggs
31		#' (`TRUE`) of __sterile__ eggs (`FALSE`).
32		#' @param n Mandatory. An integer number of the surveyed individuals in the
33		#' corresponding group.
34		#' @param ... Used to pass arguments to specific method.
35		#'
36		#' @return A object of class `sit_egg_surveys` which can be used in [sit()].
37		#' @importFrom utils head tail
38		#' @export
39		#'
40		#' @family importing
41		#' @seealso [sit()]
42		#'
43		#' @examples
44		#' egg_surv <- data.frame(
45		#' trap = 1:3,
46		#' survey = c("2021-04-01", "2021-08-21", "2021-08-21"),
47		#' duration = rep(7L, 3),
48		#' fertile = c(TRUE, FALSE, TRUE),
49		#' n = 1:3
50		#' )
51		#' sit_egg_surveys(egg_surv)
52		sit_egg_surveys <- function(x, ...) {
53	59x	UseMethod('sit_egg_surveys')
54		}
55
56
57		#' @describeIn sit_egg_surveys Imports field survey data from egg traps.
58		#' @export
59		sit_egg_surveys.data.frame <- function(
60		x,
61		trap = "trap",
62		type = "type",
63		survey = "survey",
64		activation = "activation",
65		duration = "duration",
66		species = "species",
67		fertile = "fertile",
68		n = "n",
69		...
70		) {
71		# # Debug
72		# args <- names(formals(sit_egg_surveys.data.frame)[c(-1, -11)])
73		# for (arg in args) assign(arg, arg)
74
75		## The requirement of the `type` argument can only be verified when building
76		## the `sit` object, by looking at the codes and types of the available traps.
77
78		## Check here whether there are variables in the data set corresponding
79		## to optional arguments and set them to NULL if not.
80	29x	opt_args <- c(type, activation, duration, species)
81	29x	for (arg in opt_args[!opt_args %in% names(x)]) {
82	86x	assign(arg, NULL)
83		}
84
85		## Exactly one of `activation` or `duration` are required
86	29x	if ((l_ad <- length(c(activation, duration))) != 1L) {
87	2x	if (l_ad == 0L) {
88	1x	stop("Either activation or duration are required arguments. ",
89	1x	"Both are currently missing.")
90		}
91	1x	if (l_ad == 2L) {
92	1x	stop("Only one of `activation` or `duration` needed. ",
93	1x	"Please remove one.")
94		}
95		}
96
97		## Build the internal representation
98	27x	inames <- c(
99	27x	"trap_code", "trap_type"[!is.null(type)], "survey",
100	27x	"activation"[!is.null(activation)],
101	27x	"duration"[!is.null(duration)],
102	27x	"species"[!is.null(species)], "fertile", "n"
103		)
104		# Variable names reported through non-NULL arguments
105		# vnames <- unlist(lapply(tail(head(names(formals()), -1), -1), get,
106		# pos = environment()))
107		# Alternatively:
108	27x	vnames <- unname(
109	27x	unlist(
110	27x	as.list(
111	27x	environment())[tail(head(names(formals()), -1), -1)]
112		)
113		)
114
115	27x	ans <- build_internal_df(x, varnames = vnames, internal_names = inames)
116
117		## Check duplicated surveys
118	27x	if (any(idx <- duplicated(ans))) {
119	1x	stop(
120	1x	glue(
121	1x	"{vars_collapse(which(idx), 'Row ', 'Rows ', quote = '')} ",
122	1x	"duplicated. Please fix."
123		)
124		)
125		}
126
127		## Compute starting and ending dates and times
128	26x	ans <- compute_periods(
129	26x	ans,
130	26x	input = ifelse(is.null(activation), "duration", "activation")
131		)
132
133		## Survey id
134	24x	ans$id <- if (nrow(ans) > 0) {
135	16x	seq.int(nrow(ans))
136		} else {
137	8x	integer()
138		}
139
140		## Reorder columns
141	24x	ans <- ans[
142	24x	, c("id",
143	24x	"trap_code", # Changed into `trap_id` within `sit()`
144	24x	"trap_type"[!is.null(type)], # Possibly needed for `trap_id`
145	24x	"datetime_start",
146	24x	"datetime_end",
147	24x	"species"[!is.null(species)],
148	24x	"fertile",
149	24x	"n"
150		)
151		]
152
153	24x	class(ans) <- c("sit_egg_surveys", class(ans))
154	24x	return(ans)
155		}
156
157
158		#' @param area Character vector. Either `control`, `sit` or both (default).
159		#' @param trap_type Character vector. Any subset of
160		#' `unique(sit_trap_types(x)$label)`.
161		#' @param following_releases A `sit_release` object with a subset of release
162		#' events or missing (default) for all release events. Use in combination with
163		#' `following_days` to filter surveys of the target populations within a given
164		#' number of days after release. Note that counts of wild populations will
165		#' always be included in the results with a value of `NA` in `pop_col`.
166		#' @param following_days Integer or missing (default). Number of days after
167		#' releases to return, if `following_releases` is not missing.
168		#' @describeIn sit_egg_surveys Retrieve egg survey data, possibly filtering
169		#' results.
170		#' @importFrom rlang .data .env
171		#' @importFrom tidyr everything
172		#' @export
173		sit_egg_surveys.sit <- function(
174		x,
175		area = c('control', 'sit'),
176		trap_type = sit_trap_types(x)$label,
177		following_releases,
178		following_days,
179		species = NULL,
180		...) {
181
182	30x	target_periods <- get_periods(following_releases, following_days)
183
184		## Filter trap areas, trap types, and dates
185	30x	x_sub <-
186	30x	x %>%
187	30x	dm_filter(traps = .data$area %in% .env$area) %>%
188	30x	dm_filter(trap_types = .data$label %in% .env$trap_type) %>%
189	30x	dm_filter(
190	30x	egg_surveys =
191	30x	is_date_within(
192	30x	.data$datetime_end,
193	30x	.env$target_periods$starts,
194	30x	.env$target_periods$ends
195		)
196		)
197
198		## Filters for species if required
199	30x	if (!is.null(species)) {
200	5x	if (!'species' %in% names(x[['egg_surveys']])) {
201	1x	stop(
202	1x	"There is no variable `species` in the egg surveys."
203		)
204		}
205	4x	x_sub <- dm_filter(
206	4x	x_sub,
207	4x	egg_surveys = (.data$species %in% .env$species)
208		)
209		}
210
211	29x	ans <- x_sub %>%
212		# dm_get_tables() %>% `[`(c("traps", "trap_types"))
213		## Select and rename columns from relevant tables
214	29x	dm_select('release_events', .data$id, pop_col = .data$colour) %>%
215	29x	dm_select('traps', .data$id, trap_code = .data$code) %>%
216	29x	dm_flatten_to_tbl('egg_surveys') %>%
217	29x	select(-.data$trap_id, -.data$geometry) %>%
218	29x	select(.data$id, .data$trap_code, everything()) %>%
219	29x	set_class('sit_egg_surveys')
220
221	29x	return(ans)
222		}
223

1		#' Default table of trap types
2		#' @return A `data.frame` with the table of default trap types.
3		#' @keywords internal
4		default_trap_types <- function() {
5	40x	ans <- data.frame(
6	40x	id = 1:3,
7	40x	name = c("Ovitrap", "BG-Sentinel", "Human Landing Catch"),
8	40x	label = c("OVT", "BGS", "HLC"),
9	40x	stage = c("egg", "adult", "adult"),
10	40x	description = NA
11		)
12
13	40x	class(ans) <- c("sit_trap_types", class(ans))
14
15	40x	return(ans)
16		}

1		#' Validator
2		#'
3		#' Check the internal consistency of a `sit` object.
4		#'
5		#' Performs the following checks:
6		#'
7		#' - Consistency of survey dates. Surveys of sterile males must precede their
8		#' release date. _Dummy_ surveys with 0 counts will be silently removed.
9		#'
10		#' @param x An object.
11		#'
12		#' @return Returns the `sit`, invisibly, after finishing all checks.
13		#'
14		#' @keywords internal
15		#' @importFrom rlang .data .env
16		#' @export
17		validate_sit <- function(x) {
18
19	15x	stopifnot(inherits(x, 'sit'))
20
21		## Surveys of sterile males must precede their release date.
22	15x	ages <- survey_ages(sit_adult_surveys(x), sit_revents(x))
23
24	15x	inconsistent_dates <- ages$age < 0 & ages$n > 0
25	15x	if ( (sid <- sum(inconsistent_dates)) > 0) {
26	!	stop(
27	!	"Inconsistent dates. ",
28	!	"There are ", sid,
29	!	" counts of sterile males surveyed at a date earlier than their release.")
30		}
31
32		## Remove dummy surveys
33	15x	dummy_surveys <- ages$age <= 0 # must be n = 0
34	15x	if (any(dummy_surveys)) {
35	1x	dummy_surv_ids <- ages$id[dummy_surveys]
36		# x <- dm_filter(x, adult_surveys = (!.data$id %in% .env$dummy_surv_ids))
37	1x	x <- dm_mutate_tbl(
38	1x	x,
39	1x	adult_surveys = x[["adult_surveys"]][
40	1x	-match(dummy_surv_ids, x[["adult_surveys"]][["id"]]),
41		]
42		)
43	1x	class(x) <- c("sit", class(x))
44		}
45
46	15x	invisible(x)
47		}