R/CJS.R

In hugomflavio/actel: Acoustic Telemetry Data Analysis

#' CJS.R arguments
#' @param replicates A list of arrays containing the respective replicate stations.
#' @param CJS The overall CJS results.
#' @param detections.list A list of the detections split by each target tag.
#' @param spatial The spatial data frame.
#' @param efficiency The efficiency results.
#' @param arrays A list containing information for each array.
#' @param mat,m The presence/absence matrices.
#' @param input A presence/absence matrix.
#' @param movements,moves The movements table.
#' @param status.df The main results table.
#' @param dotmat The matrix of distances between arrays.
#' @param paths A list containing the shortest paths between arrays with distance > 1.
#' @param estimate An estimate of the last array's efficiency, between 0 and 1.
#' @param fixed.efficiency A vector of fixed efficiency estimates \[0, 1\]. `length(fixed.efficiency)` must match `ncol(input)`.
#' @param silent Logical: Should messages be printed? This argument is mainly intended for function calls running within actel's analyses.
#' @name cjs_args
#' @keywords internal
NULL

#' Prepare intra-array matrices for selected arrays
#'
#' @inheritParams cjs_args
#'
#' @return A list of dual matrices for each relevant array.
#'
#' @keywords internal
#'
getDualMatrices <- function(replicates, CJS = NULL, spatial, detections.list) {
  event(type = "debug", "Running getDualMatrices.")
  output <- list()
  for(i in 1:length(replicates)) {
    continue <- TRUE
    if (!is.null(CJS)) {
      if (any(grepl("^max\\.efficiency$", names(CJS))))
        already.calculated <- !is.na(CJS$max.efficiency[names(replicates)[i]])
      else
        already.calculated <- !is.na(CJS$efficiency[names(replicates)[i]])
      if (already.calculated) {
        event(type = c("screen", "warning", "report"),
              "An inter-array efficiency has already been calculated",
              " for array ", names(replicates)[i],".")
        decision <- userInput("Do you want to replace this with an intra-array efficiency estimate?(y/n) ",
                                choices = c("y", "n"), hash = "# replace efficiency?")
        if (!interactive() | decision == "y") {
          event(type = "report",
                "Replacing efficiency estimation.")
          continue <- TRUE
        } else { # nocov start
          event(type = "report",
                "Keeping inter-array efficiency estimation.")
          continue <- FALSE
        } # nocov end
      }
    }
    if (continue) {
      output[[length(output) + 1]] <- dualMatrix(array = names(replicates)[i],
        replicates = replicates[[i]], spatial = spatial, detections.list = detections.list)
      names(output)[[length(output)]] <- names(replicates)[i]
    }
  }
  return(output)
}

#' Incorporate intra-array estimates in the overall CJS object
#'
#' @param m A list of dual matrices
#' @inheritParams cjs_args
#'
#' @return A list containing the updated overall CJS and the intra-array CJS results
#'
#' @keywords internal
#'
includeIntraArrayEstimates <- function(m, efficiency = NULL, CJS = NULL) {
  event(type = "debug", "Running includeIntraArrayEstimates.")
  if (!is.null(efficiency) & !is.null(CJS)) {
    event(type = "stop",
          "Use only one of 'efficiency' or 'CJS' at a time.")
  }
  if (length(m) > 0) {
    intra.CJS <- lapply(m, dualArrayCJS)
    if (!is.null(CJS)) {
      for (i in names(intra.CJS)) {
        CJS$absolutes["estimated", i] <- round(CJS$absolutes["detected", i] / intra.CJS[[i]]$combined.efficiency, 0)
        CJS$efficiency[i] <- intra.CJS[[i]]$combined.efficiency
      }
    }
    if (!is.null(efficiency)) {
      for (i in names(intra.CJS)) {
        efficiency$max.efficiency[i] <- intra.CJS[[i]]$combined.efficiency
        efficiency$min.efficiency[i] <- intra.CJS[[i]]$combined.efficiency
      }
    }
  } else {
    intra.CJS <- NULL
  }
  if (!is.null(CJS))
    return(list(CJS = CJS, intra.CJS = intra.CJS))
  if (!is.null(efficiency))
    return(list(efficiency = efficiency, intra.CJS = intra.CJS))
  return(list(intra.CJS = intra.CJS))
}


#' Assembles CJS tables for all group x release site combinations
#'
#' @param CJS A list of CJS calculated for each group x release site x array combinations
#' @param arrays a list of arrays
#'
#' @return A list containing the CJS results for each group x release site combination.
#'
#' @keywords internal
#'
assembleSplitCJS <- function(mat, CJS, arrays, releases, intra.CJS = NULL) {
  event(type = "debug", "Running assembleSplitCJS.")
  recipient <- lapply(names(CJS), function(i) {
    aux <- releases[releases$Combined == i, ]
    output <- assembleArrayCJS(mat = mat[i], CJS = CJS[[i]], arrays = arrays, releases = aux)[[1]]
    if (!is.null(intra.CJS)) {
      for (i in names(intra.CJS)) {
        output["estimated", i] <- round(output["detected", i] / intra.CJS[[i]]$combined.efficiency, 0)
      }
    }
    output["difference", ] <- output["estimated", ] - output["known", ]
    return(output)
  })
  names(recipient) <- names(CJS)
  return(recipient)
}

#' Assembles CJS tables for all groups
#'
#' @param CJS A list of CJS calculated for each group x array combinations
#' @inheritParams cjs_args
#'
#' @return A data frame containing the CJS results for each group
#'
#' @keywords internal
#'
assembleGroupCJS <- function(mat, CJS, arrays, releases, intra.CJS = NULL) {
  event(type = "debug", "Running assembleGroupCJS.")
  recipient <- lapply(names(CJS), function(i) {
    link <- grepl(paste0("^", i), names(mat))
    aux <- releases[releases$Group == i, ]
    output <- assembleArrayCJS(mat = mat[link], CJS = CJS[[i]], arrays = arrays, releases = aux)[[1]]
    if (!is.null(intra.CJS)) {
      for (i in names(intra.CJS)) {
        output["estimated", i] <- round(output["detected", i] / intra.CJS[[i]]$combined.efficiency, 0)
      }
    }
    output["difference", ] <- output["estimated", ] - output["known", ]
    return(output)
  })
  names(recipient) <- names(CJS)
  return(recipient)
}

#' Break the detection matrices per array
#'
#' @param type The type of arrays to be matched
#' @param verbose Logical: Should [event()] be used to show and log messages?
#' @inheritParams cjs_args
#'
#' @return A list containing the split matrices for each array
#'
#' @keywords internal
#'
breakMatricesByArray <- function(m, arrays, type = c("peers", "all"), verbose = TRUE) {
  event(type = "debug", "Running breakMatricesByArray.")
  type <- match.arg(type)
  recipient <- list()
  for (i in 1:length(arrays)) {
    if ((type == "peers" & !is.null(arrays[[i]]$after.peers)) | (type == "all" & !is.null(arrays[[i]]$all.after))) {

      # find out relevant arrays
      if (type == "peers")
        a.regex <- paste0("^", c(names(arrays)[i], arrays[[i]]$after.peers), "$", collapse = "|")
      else
        a.regex <- paste0("^", c(names(arrays)[i], arrays[[i]]$all.after), "$", collapse = "|")

      # grab only relevant arrays
      aux  <- lapply(m, function(m_i) m_i[, which(grepl(a.regex, colnames(m_i))), drop = FALSE])

      # Failsafe in case some tags are released at one of the peers
      keep <- unlist(lapply(m, function(m_i) any(grepl(paste0("^", names(arrays)[i], "$"), colnames(m_i)))))
      aux  <- aux[keep]

      # Failsafe in case there is only one column left
      keep <- unlist(lapply(aux, ncol)) > 1
      aux  <- aux[keep]

      # reorder columns if necessary
      aux <- lapply(aux, function(x) {
        peer.cols <- colnames(x)[!grepl(paste0("^", names(arrays)[i], "$"), colnames(x))]
        return(x[, c(names(arrays)[i], peer.cols)])
      })

      # Convert peers to single column and add fake start
      aux <- lapply(aux, function(m) {
        if (ncol(m) > 2) {
          m[, 2:ncol(m)] <- apply(m[, 2:ncol(m)], 2, as.logical)
          m$AnyPeer <- as.numeric(apply(m[-1], 1, any))
        } else {
          colnames(m)[2] <- "AnyPeer"
        }

        # The fake start prevents the CJS functions from breaking if the efficiency of the array is 0
        m$FakeStart <- rep(1, nrow(m))
        return(m[, c(ncol(m), 1, (ncol(m) - 1))])
      })

      # If all peers are 0, the CJS functions will crash. The same happens if the array is all 0's
      own.zero.check <- unlist(lapply(aux, function(x) sum(x[, 2]) == 0))
      peer.zero.check <- unlist(lapply(aux, function(x) sum(x$AnyPeer) == 0))
      zero.check <- all(own.zero.check) | all(peer.zero.check)

      if (zero.check) {
        if (all(own.zero.check) & verbose) {
          event(type = c("screen", "warning", "report"),
                "No tags passed through array ", names(arrays)[i],
                ". Skipping efficiency estimations for this array.")
        } else {
          if (all(peer.zero.check) & verbose)
            event(type = c("screen", "warning", "report"),
                  "No tags passed through any of the efficiency peers",
                  " of array ", names(arrays)[i],
                  ". Skipping efficiency estimations for this array.")
          }
      } else {
        recipient[[length(recipient) + 1]] <- aux
        names(recipient)[length(recipient)] <- names(arrays)[i]
      }
    }
  }
  if (length(recipient) == 0 & verbose) {
    event(type = c("screen", "warning", "report"),
          "None of the arrays has valid efficiency peers.")
    return(NULL)
  } else {
    return(recipient)
  }
}

#' Combine the individual CJS's of each array into a single table
#'
#' @param CJS A list of CJS calculated for each array
#' @inheritParams cjs_args
#'
#' @return A list containing the CJS absolute numbers and efficiency estimates
#'
#' @keywords internal
#'
assembleArrayCJS <- function(mat, CJS, arrays, releases, silent = TRUE) {
  event(type = "debug", "Running assembleArrayCJS.")
  # Compile final objects
  absolutes <- matrix(nrow = 5, ncol = length(arrays))
  colnames(absolutes) <- names(arrays)
  rownames(absolutes) <- c("detected", "here plus on peers", "not here but on peers", "known", "estimated")
  absolutes <- as.data.frame(absolutes)
  efficiency <- rep(NA, length(arrays))
  names(efficiency) <- names(arrays)
  for (i in 1:length(CJS)) {
    absolutes[, names(CJS)[i]] <- CJS[[i]]$absolutes[, 2]
    efficiency[names(CJS)[i]] <- CJS[[i]]$efficiency[2]
  }
  # Fix max estimated based on arrays before, if needed
  for (i in names(arrays)) {
    # if an estimation is made and previous arrays exist
    if (!is.na(absolutes["estimated", i]) && !is.null(arrays[[i]]$before)) {
      # if estimations were made for all the previous arrays
      if (all(!is.na(absolutes["estimated", arrays[[i]]$before]))) {
        # if tags were included in the system in the target array
        if (any(releases$Array == i))
          # the max is the sum of the estimated for the previous arrays plus the released tags.
          the.max <- sum(absolutes["estimated", arrays[[i]]$before]) + sum(releases$n[releases$Array == i])
        else
          # the max is the sum of the estimated for the previous arrays.
          the.max <- sum(absolutes["estimated", arrays[[i]]$before])
        #if the estimated is too high, correct it.
        if (absolutes["estimated", i] > the.max)
          absolutes["estimated", i] <- the.max
      }
    }
  }
  # Final fixes
  fix.zero <- na.as.false(efficiency == 0)
  fix.peers <- is.na(efficiency)
  for (i in 1:length(arrays)) {
    if (fix.zero[i]) {
      # fix estimated for arrays with 0 efficiency. # if any of the before is NA, this will be NA too.
      if (!silent) {
        event(type = c("screen", "warning", "report"),
              "Array '", names(arrays)[i], "' has 0 efficiency. Attempting",
              " to round estimated based on peers, if possible.")
      }
      if (!is.null(arrays[[i]]$before))
        absolutes["estimated", i] <- sum(absolutes["estimated", arrays[[i]]$before])
    }
    if (fix.peers[i]) {
      # fix absolutes for arrays with no peers
      absolutes["detected", i] <- sum(unlist(lapply(mat, function(m) {
        if (any(colnames(m) == names(arrays)[i]))
          return(sum(m[, names(arrays)[i]]))
      })))
      # Fix knowns for arrays with no peers
      absolutes["known", i] <- absolutes["detected", i]
    }
  }
  return(list(absolutes = absolutes, efficiency = efficiency))
}

#' Assemble detection matrices
#'
#' @inheritParams cjs_args
#'
#' @return A list containing the detection matrices split by groups and release sites
#'
#' @keywords internal
#'
assembleMatrices <- function(spatial, movements, status.df, arrays, paths, dotmat) {
  event(type = "debug", "Running assembleMatrices.")
  temp <- efficiencyMatrix(movements = movements, arrays = arrays, paths = paths, dotmat = dotmat)
  output <- lapply(temp, function(x) {
    # include transmitters that were never detected
    x <- includeMissing(x = x, status.df = status.df)

    # sort the rows by the same order as status.df (I think these two lines are not needed, but leaving them in just in case)
    link <- sapply(status.df$Transmitter, function(i) grep(paste0("^", i, "$"), rownames(x)))
    x <- x[link, ]

    # split by group*release site combinations
    aux <- split(x, paste0(status.df$Group, ".", status.df$Release.site))

    # Re-order
    the.order <- c()
    for (i in unique(status.df$Group)) {
      the.order <- c(the.order, paste0(i, ".", unique(spatial$release.sites$Standard.name)))
    }
    aux <- aux[order(match(names(aux), the.order))]

    unique.release.arrays <- unique( # only keep each name once
                              unlist( # turn output into a string
                                sapply(spatial$release.sites$Array, function(x) {
                                    unlist(strsplit(x, "|", fixed = TRUE)) # break arrays by '|'
                                  })
                                )
                              )

    # If the release sites start in different arrays, trim the matrices as needed
    if (length(unique.release.arrays) > 1) {
      for(i in 1:length(aux)){ # for each matrix, find the corresponding release site.
        the_release_site <- sapply(spatial$release.sites$Standard.name, function(x) grepl(paste0("\\.", x, "$"), names(aux)[i]))
        if(sum(the_release_site) > 1) { # if there is more than one matching release site, stop.
          event(type = "stop",
                "Multiple release sites match the matrix name.",
                " Make sure that the release sites' names are not contained",
                " within the animal groups or within themselves.\n")
        }
        # else, find which is the first column to keep. This is tricky for multi-branch sites...
        the.col <- min(which(grepl(spatial$release.sites$Array[the_release_site], colnames(aux[[i]]))))
        # then keep only the relevant columns
        aux[[i]] <- aux[[i]][, c(1, the.col:ncol(aux[[i]]))]
      }
    }
    return(aux)
  })
  return(output)
}

#' Analytical CJS model
#'
#' Computes an analytical CJS model for a presence/absence matrix.
#'
#' @references Perry et al (2012), 'Using mark-recapture models to estimate survival from telemetry data'. url: <https://www.researchgate.net/publication/256443823_Using_mark-recapture_models_to_estimate_survival_from_telemetry_data>
#'
#' @inheritParams cjs_args
#'
#' @examples
#' # prepare a dummy presence/absence matrix
#' x <- matrix(c(TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE), ncol = 3)
#' colnames(x) <- c("Release", "Array1", "Array2")
#'
#' # run CJS
#' simpleCJS(x)
#'
#' @return A list containing:
#' \itemize{
#'   \item \code{absolutes} A data frame with the absolute number of tags detected and missed,
#'   \item \code{efficiency} A vector of calculated array detection efficiencies,
#'   \item \code{survival} A matrix of calculated survivals,
#'   \item \code{lambda} A combined detection efficiency * survival estimate for the last array.
#' }
#'
#' @export
#'
simpleCJS <- function(input, estimate = NULL, fixed.efficiency = NULL, silent = TRUE){
  event(type = "debug", "Running simpleCJS")
  if (!is.matrix(input) & !is.data.frame(input)) {
    event(type = "stop",
          "input must be a matrix or data frame containing only 0's and 1's.")
  }

  # stop if there is weird data in the input
  if (any(input != 0 & input != 1)) {
    event(type = "stop",
          "input must be a matrix or data frame containing only 0's and 1's.")
  }

  # stop if both estimate and fixed efficiency are present
  if (!is.null(estimate) & !is.null(fixed.efficiency)) {
    event(type = "stop",
          "Please choose only one of 'estimate' or 'fixed.efficiency'.")
  }

  if (any(input[, 1] == 0)) {
    event(type = "stop",
          "The first column of the input should only contain 1's",
          " (i.e. release point).")
  }

  # Only check below if estimate is set
  if (!is.null(estimate)) {
    # stop if multiple estimate values are provided.
    if (!is.null(estimate) && length(estimate) != 1)
      event(type = "stop",
            "Please use only one value for estimate.")
    # stop if estimate exceeds 1
    if (!is.null(estimate) && (estimate < 0 | estimate > 1))
      event(type = "stop",
            "'estimate' must be between 0 and 1.")
    # all good
  }

  # Only check below if fixed.efficiency was set
  if (!is.null(fixed.efficiency)) {
    # stop if there are not enough efficiency values
    if (length(fixed.efficiency) != ncol(input)) {
     event(type = "stop",
           "Fixed efficiency was set but its length is not the same as",
           " the number of columns in the input.")
    }
    # stop if any efficiency exceeds 1
    if (any(fixed.efficiency > 1, na.rm = TRUE)) {
      event(type = "stop", "Fixed efficiency estimates must be",
            " between 0 and 1.")
    }
    # Fake estimate to avoid further changes below.
    if (!is.na(tail(fixed.efficiency, 1)))
     estimate <- tail(fixed.efficiency, 1)
    # all good
    if(!silent) {
      event(type = "screen",
            "M: Running CJS with fixed efficiency estimates.")
    }
  }

  # Start the calculations

  # S: Survival probability between arrays
  # m: tags detected at i
  # r: tags detected at i and on peers
  # z: tags NOT detected at i but detected on peers
  # p: probability of detection (efficiency)
  # M: animals estimated to be alive at i

  S <- rep(NA, ncol(input) - 1)
  r <- z <- p <- m <- M <- rep(NA, ncol(input))
  for(i in 1:(ncol(input) - 1)){
    # number of tags detected at i and on peers (r)
    tr <- input[input[, i] == 1, (i + 1):ncol(input)]
    if (is.data.frame(tr) | is.matrix(tr))
      r[i] = sum(apply(tr, 1, function(f) any(f == 1)))
    if (is.vector(tr))
      r[i] = sum(tr)
    # number of tags NOT detected at i but detected on peers (z)
    tz <- input[input[, i] == 0, (i + 1):ncol(input)]
    if (is.data.frame(tz) | is.matrix(tz))
      z[i] = sum(apply(tz, 1, function(f) any(f == 1)))
    if (is.vector(tz))
      z[i] = sum(tz)
    # probability of detection (p)
    if (is.null(fixed.efficiency)) {
      if (r[i] == 0 & z[i] == 0)
        p[i] <- 0
      else
        p[i] <- r[i] / (r[i] + z[i])
    } else {
      p[i] <- fixed.efficiency[i]
    }
    # number of detected tags at i (m)
    m[i] = sum(input[, i])
    # number of animals estimated alive at i (M)
    # Failsafe for array with 0 efficiency. Issues warning.
    if (p[i] == 0 | m[i] == 0) {
      if(p[i] == 0 & !silent) {
        event(type = c("screen", "warning"),
              "Array '", colnames(input)[i],
              "' had 0% efficiency. Skipping survival estimation.")
      }
      if(m[i] == 0 & !silent) {
        event(type = c("screen", "warning"),
              "No tags were detected at array '", colnames(input)[i],
              "'. Skipping survival estimation.")
      }
      M[i] = M[i - 1]
      S[i - 1] = -999
      if (i == (ncol(input)-1))
        l <- 0
    } else {
      M[i] = round(m[i] / p[i], 0)
      # Fix M if the fixed efficiency value draws it too low.
      if (!is.null(fixed.efficiency) && M[i] < sum(m[i], z[i])){
        # Prevent fairly common warning from showing up while inside Actel.
        if(!silent) {
          event(type = c("screen", "warning"),
                "The fixed efficiency caused a too low estimate at iteration ",
                i,". Forcing higher estimate.")
        }
        M[i] = sum(m[i], z[i])
      }
      # Correct M if the estimate at i is bigger than the estimate at i - 1
      if (i > 1 && M[i] > M[i - 1])
        M[i] = M[i - 1]
      # survival probability from i - 1 to i (S)
      if (i > 1)
        S[i - 1] <- M[i] / M[i - 1]
      # lambda (last S * last P)
      if (i == (ncol(input)-1))
        l <- r[i] / m[i]
    }
  }
  m[ncol(input)] <- sum(input[, ncol(input)])

  # Complete final values if estimate is present
  if (!is.null(estimate)) {
    if (estimate == 0) {
      S[ncol(input) - 1] <- -999
    } else {
      S[ncol(input) - 1] <- l / estimate
      if (S[ncol(input) - 1] > 1)
        S[ncol(input) - 1] = 1
    }
    M[ncol(input)] <- round(m[ncol(input)] / estimate, 0)
    if (M[ncol(input)] > M[ncol(input) - 1])
      M[ncol(input)] = M[ncol(input) - 1]
    p[ncol(input)] = estimate
  }

  # prepare printings
  k <- sapply(1:length(m), function(i) sum(c(m[i], z[i]), na.rm = TRUE))
  absolutes <- matrix(c(m, r, z, k, M), ncol = ncol(input), byrow = TRUE)
  rownames(absolutes) <- c("detected", "here plus on peers", "not here but on peers", "known", "estimated")
  colnames(absolutes) <- colnames(input)
  names(p) <- colnames(input)
  link <- is.na(S) | S >= 0
  survival <- matrix(S[link], nrow = sum(link))
  if (any(!link))
    the.names <- colnames(input)[-(which(!link) + 1)]
  else
    the.names <- colnames(input)
  maxcharA <- max(nchar(the.names[-length(the.names)]))
  maxcharB <- max(nchar(the.names[-1]))
  the.rows <- c()
  for (i in 1:(length(the.names) - 1)) {
    if(nchar(the.names[i] != maxcharA))
      centeringA <- paste(rep(" ", maxcharA - nchar(the.names[i])), collapse = "")
    else
      centeringA <- ""
    if(nchar(the.names[i] != maxcharB))
      centeringB <- paste(rep(" ", maxcharB - nchar(the.names[i + 1])), collapse = "")
    else
      centeringB <- ""
    the.rows[length(the.rows) + 1] <- paste0(centeringA, the.names[i], " -> ", the.names[i + 1], centeringB, " =")
  }
  rownames(survival) <- the.rows
  colnames(survival) <- ""

  return(list(absolutes = absolutes, efficiency = p, survival = survival, lambda = l))
}

#' Calculate CJS for each group.release combination
#'
#' @inheritParams cjs_args
#'
#' @return A list containing CJS results split by group and release site
#'
#' @keywords internal
#'
mbSplitCJS <- function(mat, fixed.efficiency = NULL) {
  event(type = "debug", "Running mbSplitCJS.")
  recipient <- lapply(mat, function(m) {
    if (is.na(fixed.efficiency[grepl(paste0("^", colnames(m[[1]])[2], "$"), names(fixed.efficiency))]))
      array.efficiency <- NULL
    else
      array.efficiency <- c(1, fixed.efficiency[grepl(paste0("^", colnames(m[[1]])[2], "$"), names(fixed.efficiency))], 1)
    output <- lapply(m, function(mm) {
      simpleCJS(mm, fixed.efficiency = array.efficiency, silent = TRUE)
    })
    return(output)
  })
  groups <- unique(unlist(lapply(recipient, names)))
  output <- list()
  for (group in groups) {
    output[[length(output) + 1]] <- lapply(recipient, function(r) {
      if (any(grepl(group, names(r))))
        return(r[[which(grepl(group, names(r)))]])
      else
        return(NULL)
    })
  }
  output <- lapply(output, function(x) {
    keep <- !unlist(lapply(x, is.null))
    return(x[keep])
  })
  names(output) <- groups
  return(output)
}

#' Calculate CJS for each group for each array
#'
#' @inheritParams cjs_args
#'
#' @return A list containing CJS results split by group.
#'
#' @keywords internal
#'
mbGroupCJS <- function(mat, status.df, fixed.efficiency = NULL) {
  event(type = "debug", "Running mbGroupCJS.")
  output <- list()
  for (i in 1:length(unique(status.df$Group))) {
    output[[i]] <- lapply(mat, function(m_i) {
      if (is.na(fixed.efficiency[grepl(paste0("^", colnames(m_i[[1]])[2], "$"), names(fixed.efficiency))]))
        array.efficiency <- NULL
      else
        array.efficiency <- c(1, fixed.efficiency[grepl(paste0("^", colnames(m_i[[1]])[2], "$"), names(fixed.efficiency))], 1)
      link <- grepl(paste0("^", unique(status.df$Group)[i]), names(m_i))
      if (any(link)) {
        if(sum(link) == 1)
          return(simpleCJS(m_i[[which(link)]], fixed.efficiency = array.efficiency, silent = TRUE))
        else
          return(combineCJS(m_i[link], fixed.efficiency = array.efficiency, silent = TRUE))
      } else {
        return(NULL)
      }
    })
    names(output[[i]]) <- names(mat)
  }
  output <- lapply(output, function(x) x[!unlist(lapply(x, is.null))])
  names(output) <- unique(status.df$Group)
  return(output)
}

#' Compile inter-array detection matrix
#'
#' @inheritParams cjs_args
#'
#' @return A matrix of detection histories per tag.
#'
#' @keywords internal
#'
efficiencyMatrix <- function(movements, arrays, paths, dotmat) {
  event(type = "debug", "Running efficiencyMatrix.")
  max.ef <- as.data.frame(matrix(ncol = length(arrays) + 1, nrow = length(movements)))
  colnames(max.ef) <- c("Release", names(arrays))
  rownames(max.ef) <- names(movements)
  max.ef[is.na(max.ef)] = 0
  max.ef$Release = 1
  min.ef <- max.ef

  capture <- lapply(names(movements), function(tag) {
    max.aux <- c(1, rep(0, length(arrays)))
    names(max.aux) <- c("Release", names(arrays))
    one.way <- oneWayMoves(movements = movements[[tag]], arrays = arrays)
    if (!is.null(one.way)) {
      max.aux[match(one.way$Array, names(max.aux))] <- 1
      min.aux <- max.aux
      if (nrow(one.way) > 1) {
        aux <- countArrayFailures(moves = one.way, paths = paths, dotmat = dotmat)
        aux <- unique(aux[grepl("unsure", names(aux))])
        if (!is.null(aux))
          max.aux[match(aux, names(max.aux))] <- 1
      }
      max.ef[tag, ] <<- max.aux
      min.ef[tag, ] <<- min.aux
    }
  })
  return(list(maxmat = max.ef, minmat = min.ef))
}

#' Trim movements table to contain only uni-directional movements
#'
#' @inheritParams cjs_args
#'
#' @return A data frame with the uni-directional movements for the target tags.
#'
#' @keywords internal
#'
oneWayMoves <- function(movements, arrays) {
  event(type = "debug", "Running oneWayMoves.")
  if (nrow(movements) > 1) {
    while (TRUE) {
      aux <- data.frame(from = movements$Array[-nrow(movements)], to = movements$Array[-1])
      link <- apply(aux, 1, function(x) x[2] %in% arrays[[x[1]]]$all.after)
      if (any(link)) {
        if (all(link))
          return(movements)
        else
          movements <- movements[c(TRUE, link), ]
      } else {
        return(movements[1, , drop = FALSE])
      }
    }
  } else {
    return(movements)
  }
}

#' Find and list arrays which failed during the movements of the tags
#'
#' @param moves the valid array movements
#' @inheritParams res_efficiency
#' @inheritParams dotPaths
#'
#' @return NULL if no arrays failed, or a list containing information on the arrays that failed
#'
#' @keywords internal
#'
countArrayFailures <- function(moves, paths, dotmat) {
  event(type = "debug", "Running countArrayFailures.")
  x <- lapply(1:(nrow(moves) - 1), function(i) {
    A <- moves$Array[i]
    B <- moves$Array[i + 1]
    if (A != B & dotmat[A, B] != 1)
      blameArrays(from = A, to = B, paths = paths)
    else
      NULL
  })
  return(unlist(x))
}

#' Find which arrays to blame for a jump in movement events
#'
#' @param from The array where the tag started
#' @param to The array where the tag was next detected
#' @inheritParams res_efficiency
#'
#' @return A list containing information on the arrays that failed
#'
#' @keywords internal
#'
blameArrays <- function(from, to, paths) {
  event(type = "debug", "Running blameArrays.")
  the.paths <- paths[[paste0(from, "_to_", to)]]
  if (is.null(the.paths)) {
    event(type = "stop",
          "Either 'from' is not connected to 'to', or both are neighbours.")
  }
  if (length(the.paths) == 1) {
    return(list(known = unique(unlist(strsplit(the.paths, " -> ")))))
  } else {
    aux <- strsplit(the.paths, " -> ")
    aux <- lapply(aux, unique)
    combined.unsure <- unique(unlist(aux))
    ocurrences <- table(unlist(aux))
    if (any(ocurrences == length(aux))) {
      combined.knowns <- names(ocurrences)[ocurrences == length(aux)]
      combined.unsure <- combined.unsure[!combined.unsure %in% combined.knowns]
    } else {
      combined.knowns <- NULL
    }
    return(list(known = combined.knowns, unsure = combined.unsure))
  }
}

#' Include tags that were never detected
#'
#' @param x an efficiency matrix
#' @inheritParams cjs_args
#'
#' @return A matrix of detection histories per tag, including tags that were never detected.
#'
#' @keywords internal
#'
includeMissing <- function(x, status.df){
  event(type = "debug", "Running includeMissing.")
  missing_transmitters <- status.df$Transmitter[!status.df$Transmitter %in% rownames(x)]
  x[missing_transmitters, ] <- 0
  x[missing_transmitters, 1] <- 1
  return(x)
}

#' Compile detection matrix for last array
#'
#' @inheritParams cjs_args
#'
#' @return A matrix of detection histories per tag for the last array.
#'
#' @keywords internal
#'
dualMatrix <- function(array, replicates, spatial, detections.list){
  event(type = "debug", "Running dualMatrix.")
  all.stations <- spatial$stations$Standard.name[spatial$stations$Array == array]
  original <- all.stations[!all.stations %in% replicates]
  efficiency <- as.data.frame(matrix(ncol = 2, nrow = length(detections.list)))
  colnames(efficiency) <- c("original","replicates")
  rownames(efficiency) <- names(detections.list)
  for (i in 1:length(detections.list)) {
    efficiency[i, "original"] <- any(!is.na(match(original, detections.list[[i]]$Standard.name)))
    efficiency[i, "replicates"] <- any(!is.na(match(replicates, detections.list[[i]]$Standard.name)))
  }
  colnames(efficiency) <- c("R1", "R2")
  return(efficiency)
}

#' Calculate estimated last-array efficiency
#'
#' @references Perry et al (2012), 'Using mark-recapture models to estimate survival from telemetry data'. url: <https://www.researchgate.net/publication/256443823_Using_mark-recapture_models_to_estimate_survival_from_telemetry_data>
#'
#' @inheritParams cjs_args
#'
#' @examples
#' # prepare a dummy presence/absence matrix
#' x <- matrix(c(TRUE, TRUE, TRUE, TRUE, FALSE, TRUE), ncol = 2)
#' colnames(x) <- c("R1", "R2")
#'
#' # run CJS
#' dualArrayCJS(x)
#'
#' @return A list containing:
#' \itemize{
#'   \item \code{absolutes} A matrix with the absolute number of tags detected at each replicate and at both,
#'   \item \code{single.efficiency} A vector of calculated array detection efficiencies for each of the replicates,
#'   \item \code{combined.efficiency} The value of the combined detection efficiency for the array.
#' }
#'
#' @export
#'
dualArrayCJS <- function(input){
  event(type = "debug", "Running dualArrayCJS.")
  is_a_table <- inherits(input, "matrix") | inherits(input, "data.frame")
  has_2_cols <- is_a_table && ncol(input) == 2
  is_logical <- has_2_cols && all(apply(input, 2, is.logical))
  if (!is_logical) {
    event(type = "stop",
          "Please provide a data.frame or matrix of TRUE/FALSE's",
          " with two columns")
  }
  if (is.null(colnames(input)))
    colnames(input) <- c("R1", "R2")
  # r: tags detected at i and on other pair
  # z: tags NOT detected at i but detected on other pair
  # p: probability of detection (efficiency)
  r <- z <- p <- rep(NA, 2)
  for(i in 1:2){
    # number of tags detected at i and elsewhere (r)
    tr <- input[input[, i] == 1, -i]
    if (is.vector(tr))
      r[i] = sum(tr)
    # number of tags NOT detected at i but detected elsewhere (z)
    tz <- input[input[, i] == 0, -i]
    if (is.vector(tz))
      z[i] = sum(tz)
    # probability of detection (p)
    p[i] <- r[i] / (r[i] + z[i])
  }
  combined.p <- 1 - prod(1 - p)
  absolutes <- matrix(c(apply(input, 2, sum), r[1]), nrow = 3)
  colnames(absolutes) <- ""
  rownames(absolutes) <- c(paste0("detected at ", colnames(input)[1], ": "), paste0("detected at ", colnames(input)[2], ": "), "detected at both: ")
  names(p) <- colnames(input)
  return(list(absolutes = absolutes, single.efficiency = p, combined.efficiency = combined.p))
}

#' Combine multiple CJS models
#'
#' @param ... The detection matrices to be joined, or a list containing detection matrices
#' @inheritParams cjs_args
#'
#' @return A list containing:
#' \itemize{
#'   \item \code{absolutes} A data frame with the absolute number of tags detected and missed,
#'   \item \code{efficiency} A vector of calculated array detection efficiencies,
#'   \item \code{survival} A matrix of calculated survivals,
#'   \item \code{lambda} A combined detection efficiency * survival estimate for the last array.
#' }
#'
#' @keywords internal
#'
combineCJS <- function(..., estimate = NULL, fixed.efficiency = NULL, silent = FALSE){
  event(type = "debug", "Running combineCJS.")
  # stop if both estimate and fixed efficiency are present
  if (!is.null(estimate) & !is.null(fixed.efficiency)) {
    event(type = "stop",
          "Please choose only one of 'estimate' or 'fixed.efficiency'.")
  }

  # stop if multiple estimate values are provided.
  if (!is.null(estimate) && length(estimate) != 1) {
    event(type = "stop",
          "Please use only one value for estimate.")
  }

  # stop if any efficiency is outside [0, 1]
  not_null <- !is.null(fixed.efficiency)
  above_1 <- not_null && any(fixed.efficiency > 1, na.rm = TRUE)
  below_0 <- not_null && any(fixed.efficiency < 0, na.rm = TRUE)
  if (above_1 | below_0) {
    event(type = "stop",
          "Fixed efficiency estimates must be between 0 and 1.")
  }

  not_null <- !is.null(estimate)
  above_1 <- not_null && any(estimate > 1, na.rm = TRUE)
  below_0 <- not_null && any(estimate < 0, na.rm = TRUE)
  if (above_1 | below_0) {
    event(type = "stop",
          "'estimate' must be between 0 and 1.")
  }

  # hack to figure the number of arguments in the dots, because I could not find a better way around it.
  arg.names <- names(match.call())
  if (is.null(arg.names))
    discount <- 1
  else
    discount <- 1 + sum(!is.na(match(c("estimate", "fixed.efficiency", "silent"), arg.names)))

  rm(arg.names)

  # Prepare input
  if ((length(match.call()) - discount) == 1) {
    if(is.list(...)) {
      input <- (...)
      if (length(input) == 1) {
        event(type = "stop",
              "Input appears to contain a list with only one element.")
      }
    } else {
      event(type = "stop",
            "Only one object provided but it is not a list.")
    }
  } else {
    input <- list(...)
  }

  check <- lapply(input, function(x) !is.matrix(x) & !is.data.frame(x))
  if(any(unlist(check))) {
    event(type = "stop",
          "Not all objects provided are matrices or data frames.",
          " Please use either one list of matrices/data frames or",
          " multiple matrices/data frames.")
  }

  # stop if not all matrices finish in the same array
  the.last <- unlist(lapply(input, function(x) tail(colnames(x),1)))
  if (length(unique(the.last)) > 1) {
    event(type = "stop",
          "The last array is not the same in all input matrices.")
  }

  ncols <- sapply(input, ncol)
  groups <- rev(sort(unique(ncols)))

  if (!is.null(fixed.efficiency)){
    # stop if there are not enough efficiency values
    if (length(fixed.efficiency) != max(sapply(input, ncol))) {
      event(type = "stop",
            "Fixed efficiency was set but its length is not the same as the",
            " maximum number of columns in the input.")
    } else {
      if (!silent) {
        event(type = "screen",
              "M: Running CJS with fixed efficiency values.")
      }
    }
  }

  combinedmatrices <- as.data.frame(data.table::rbindlist(input))
  return(simpleCJS(combinedmatrices, estimate = estimate, fixed.efficiency = fixed.efficiency))
}