R/math.R

In GENEAcore: Pre-Processing of 'GENEActiv' Data

#' Weighted sum with logical index
#'
#' @description Computes the weighted sum of `x` using weights `w`, where `idx` is TRUE.
#'
#' @param x Numeric vector of values.
#' @param w Numeric vector of weights, same length as `x`.
#' @param idx Logical vector indicating which elements to include.
#' @keywords internal
#' @export
#' @return A single numeric value. Returns 0 if no TRUE values in `idx`.
#' @examples
#' idx_active <- c(TRUE, TRUE, FALSE, FALSE, FALSE, TRUE)
#' AGSA <- c(0.1, 0.2, 0.03, 0.04, 0.05, 0.6)
#' Duration <- c(100, 90, 80, 70, 60, 50)
#' ActiveVolume <- weighted_sum(AGSA, Duration, idx_active)
weighted_sum <- function(x, w, idx) {
  sum(x[idx] * w[idx], na.rm = TRUE)
}

#' Weighted mean with logical index
#'
#' @description Computes the weighted mean of `x` using weights `w`, where `idx` is TRUE.
#' @param x Numeric vector of values.
#' @param w Numeric vector of weights, same length as `x`.
#' @param idx Logical vector indicating which elements to include.
#' @keywords internal
#' @export
#' @return A single numeric value. Returns NA if the sum of weights is 0.
#' @examples
#' idx_active <- c(TRUE, TRUE, FALSE, FALSE, FALSE, TRUE)
#' AGSA <- c(0.1, 0.2, 0.03, 0.04, 0.05, 0.6)
#' Duration <- c(100, 90, 80, 70, 60, 50)
#' ActiveIntensity <- weighted_mean(AGSA, Duration, idx_active)
weighted_mean <- function(x, w, idx) {
  denominator <- sum(w[idx], na.rm = TRUE)
  if (is.na(denominator) || denominator == 0) {
    return(NA_real_)
  }
  weighted_sum(x, w, idx) / denominator
}


#' Compute Mx and Lx metrics
#' @details Calculates intensity for each x duration window.
#'
#' @param bouts A data frame sorted by time containing at minimum columns for timestamp, duration and an intensity metric.
#' @param time_col Name of the column containing Unix UTC timestamps. Defaults to `"TimeUTC"`.
#' @param duration_col Name of the column containing bout durations. Defaults to `"Duration"`.
#' @param intensity_col Name of the column containing bout intensity values. Defaults to `"AGSAMean"`.
#' @param window Window length in seconds.
#' @param low Set to TRUE for both the maximum (Mx) and minimum (Lx) intensities
#' are returned. If FALSE, only the maximum (Mx) is returned.
#' @returns A list containing the maximum and minimum mean intensity and their
#' corresponding start times for the supplied bouts data at a given window.
#' @keywords internal
#' @export
mx_lx_window <- function(bouts,
                         time_col = "TimeUTC",
                         duration_col = "Duration",
                         intensity_col = "AGSAMean",
                         window,
                         low = TRUE) {
  # Check that columns exist
  required_cols <- c(time_col, duration_col, intensity_col)
  missing_cols <- setdiff(required_cols, names(bouts))

  if (length(missing_cols) > 0) {
    stop(
      sprintf(
        "Column validation failed in `mx_lx_window()`. The following column(s) were not found in `bouts`: %s",
        paste(missing_cols, collapse = ", ")
      )
    )
  }

  # Sort bouts by time
  bouts <- bouts[order(bouts[[time_col]]), ]

  if (nrow(bouts) == 0) {
    return(NULL)
  }

  timeutc <- bouts[[time_col]]
  agsa <- bouts[[intensity_col]]
  duration <- bouts[[duration_col]]

  intensity <- c()

  for (k in 1:nrow(bouts)) {
    tmp <- bouts[timeutc >= timeutc[k] &
      timeutc < (timeutc[k] + window), ]

    intensity[k] <- weighted_mean(tmp[[intensity_col]], tmp[[duration_col]], rep(TRUE, nrow(tmp)))
  }

  i_max <- which.max(intensity)
  out <- list(
    M_intensity = intensity[i_max],
    M_time      = timeutc[i_max]
  )

  if (low) {
    i_min <- which.min(intensity)
    out$L_intensity <- intensity[i_min]
    out$L_time <- timeutc[i_min]
  }

  return(out)
}


#' Compute Cx percentile
#' @details Calculates the xth percentile of cadence for a given percentile.
#'
#' @param bouts A data frame sorted by time containing at minimum a `value_col` and `duration_col`.
#' @param percentile The percentile to compute. Single integer between 0 and 100. Default is 95.
#' @param value_col Name of the column containing values for which percentile is calculated. Defaults to `"StepMean"`.
#' @param duration_col Name of the column containing observation weights used to calculate percentile. Defaults to `"Duration"`.
#' @returns A single numeric value or `NA` if percentile cannot be calculated
#' @keywords internal
#' @export
cx_percentile <- function(bouts,
                          percentile = 95,
                          value_col = "StepMean",
                          duration_col = "Duration") {
  percentile <- as.integer(percentile)
  percentile <- max(0, min(100, percentile))

  # Check that columns exist
  required_cols <- c(value_col, duration_col)
  missing_cols <- setdiff(required_cols, names(bouts))

  if (length(missing_cols) > 0) {
    stop(
      sprintf(
        "Column validation failed in `cx_percentile()`. The following column(s) were not found in `bouts`: %s",
        paste(missing_cols, collapse = ", ")
      )
    )
  }

  bouts <- na.omit(bouts)

  if (nrow(bouts) > 0) {
    bouts <- bouts[order(bouts[[value_col]]), ]
    v <- bouts[[value_col]]
    d <- bouts[[duration_col]]
    cutoff <- (percentile / 100) * sum(d)
    i_cutoff <- which(cumsum(d) >= cutoff)[1]
    return(v[i_cutoff])
  } else {
    return(NA)
  }
}

#' Compute Fx metrics
#' @details Calculates frequency of bouts for each x duration window.
#'
#' @param bouts A bouts data frame sorted by time containing at minimum the time and intensity columns.
#' @param time_col Name of the column containing Unix UTC timestamps. Defaults to `"TimeUTC"`.
#' @param duration_col Name of the column containing bout durations. Defaults to `"Duration"`.
#' @param intensity_col Name of the column containing bout intensity values. Defaults to `"AGSAMean"`.
#' @param intensity_threshold Threshold for active classification in g. Defaults to 0.0625 (62.5mg).
#' @param active_duration Duration of active time in seconds. Default is 3600 seconds (60 minutes).
#' @returns A list containing the maximum and minimum bout counts and their
#' corresponding start times for the supplied bouts data at a given window.
#' @keywords internal
#' @export
fx_count <- function(bouts,
                     time_col = "TimeUTC",
                     duration_col = "Duration",
                     intensity_col = "AGSAMean",
                     intensity_threshold = 0.0625,
                     active_duration = 3600) {
  # Check that columns exist
  required_cols <- c(time_col, duration_col, intensity_col)
  missing_cols <- setdiff(required_cols, names(bouts))

  if (length(missing_cols) > 0) {
    stop(
      sprintf(
        "Column validation failed in `fx_count()`. The following column(s) were not found in `bouts`: %s",
        paste(missing_cols, collapse = ", ")
      )
    )
  }

  bouts <- bouts[order(bouts[[time_col]]), ]
  bouts$active <- ifelse(bouts[[intensity_col]] > intensity_threshold, TRUE, FALSE)

  bouts$activecumsum <- cumsum(ifelse(bouts$active, bouts[[duration_col]], 0))
  timeutc <- bouts[[time_col]]

  bout_count <- c()

  if (nrow(bouts) == 0) {
    return(NULL)
  }

  for (k in 1:nrow(bouts)) {
    start_i <- k
    idx <- which(bouts$activecumsum >= (bouts$activecumsum[k] + active_duration))
    if (length(idx) == 0) {
      end_i <- NA
    } else {
      end_i <- min(idx)
    }
    bout_count[k] <- end_i - start_i
  }

  c_max <- which.max(bout_count)
  c_min <- which.min(bout_count)

  out <- list(
    M_count = bout_count[c_max],
    M_time = timeutc[c_max],
    L_count = bout_count[c_min],
    L_time = timeutc[c_min]
  )

  return(out)
}

#' @title Radians MeanDir function
#' @param rotation A numeric vector of radian values
#' @return A single value data.
#' @export
#' @keywords internal
#' @examples
#' tmp1 <- c(1, 3, 2, 6, 4, 5, 3, 9, 10)
#' MeanDir(tmp1)
MeanDir <- function(rotation) {
  sin1 <- sum(sin(rotation), na.rm = T)
  cos1 <- sum(cos(rotation), na.rm = T)
  sin1bar <- (sin1) / length(rotation)
  cos1bar <- (cos1) / length(rotation)

  # Create tan1
  if (sin1bar > 0 & cos1bar > 0) {
    Tan1 <- atan(sin1bar / cos1bar)
  } else if (cos1bar < 0) {
    Tan1 <- atan(sin1bar / cos1bar) + pi
  } else if (sin1bar <= 0 & cos1bar >= 0) {
    Tan1 <- atan(sin1bar / cos1bar) + 2 * pi
  }

  return(Tan1)
}

#' @title Radians CirSD function
#' @param rotation A numeric vector of radian values
#' @return A single value data.
#' @export
#' @keywords internal
#' @examples
#' tmp1 <- c(1, 3, 2, 6, 4, 5, 3, 9, 10)
#' CirSD(tmp1)
CirSD <- function(rotation) {
  sin1 <- sum(sin(rotation), na.rm = T)
  cos1 <- sum(cos(rotation), na.rm = T)
  R1bar <- sqrt((sin1^2 + cos1^2)) / length(rotation)
  CirSD <- sqrt((-2 * log(R1bar)))
  return(CirSD)
}

#' @title Radian functions
#' @param rotation A numeric vector of radian values
#' @param fun Functions to summarise `rotation`
#' @return A single value data.
#' @export
#' @keywords internal
#' @examples
#' tmp1 <- c(1, 3, 2, 6, 4, 5, 3, 9, 10)
#' radian_fun(tmp1, fun = sd)
radian_fun <- function(rotation, fun) {
  if (identical(fun, mean)) {
    return(MeanDir(rotation))
  }
  if (identical(fun, sd)) {
    return(CirSD(rotation))
  }

  out <- fun(rotation)

  if ("Mean" %in% names(out)) out["Mean"] <- MeanDir(rotation)
  if ("SD" %in% names(out)) out["SD"] <- CirSD(rotation)

  out
}