alkahest coverage - 26.14%

Files
Source

# SMOOTHING
#' @include AllGenerics.R
NULL

# Unweighted smoothing =========================================================
#' @export
#' @rdname smooth_rectangular
#' @aliases smooth_rectangular,numeric,numeric-method
setMethod(
  f = "smooth_rectangular",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, m = 3) {
    assert_length(y, length(x))

    ## Windows
    index <- window_sliding(length(x), m)

    y <- vapply(
      X = index,
      FUN = function(i, data) mean(data[i]),
      FUN.VALUE = numeric(1),
      data = y
    )

    xy <- list(x = x, y = y)
    attr(xy, "method") <- "unweighted smoothing"
    xy
  }
)

#' @export
#' @rdname smooth_rectangular
#' @aliases smooth_rectangular,ANY,missing-method
setMethod(
  f = "smooth_rectangular",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, m) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, m = m)
  }
)

# Weighted smoothing ===========================================================
#' @export
#' @rdname smooth_triangular
#' @aliases smooth_triangular,numeric,numeric-method
setMethod(
  f = "smooth_triangular",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, m = 3) {
    # Validation
    assert_length(y, length(x))
    assert_odd(m)

    # Index
    k <- (m - 1) / 2
    index_k <- seq_len(k)
    index_y <- seq_along(y)
    index_m <- c(index_k, rep_len(k + 1, length(y) - 2 * k), rev(index_k)) - 1

    y <- mapply(
      FUN = function(i, k, data) {
        j <- seq_len(k)
        w <- c(j, k + 1, rev(j))
        index <- seq(from = i - k, to = i + k, by = 1)
        stats::weighted.mean(x = data[index], w = w)
      },
      i = index_y, k = index_m,
      MoreArgs = list(data = y)
    )

    xy <- list(x = x, y = y)
    attr(xy, "method") <- "weighted smoothing"
    xy
  }
)

#' @export
#' @rdname smooth_triangular
#' @aliases smooth_triangular,ANY,missing-method
setMethod(
  f = "smooth_triangular",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, m) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, m = m)
  }
)

# Loess smoothing ==============================================================
#' @export
#' @rdname smooth_loess
#' @aliases smooth_loess,numeric,numeric-method
setMethod(
  f = "smooth_loess",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, span = 0.75, ...) {
    assert_length(y, length(x))

    points <- data.frame(x, y)
    fit <- stats::loess(y ~ x, data = points, span = span, ...)
    bsl <- stats::predict(fit)

    xy <- list(x = x, y = bsl)
    attr(xy, "method") <- "loess smoothing"
    xy
  }
)

#' @export
#' @rdname smooth_loess
#' @aliases smooth_loess,ANY,missing-method
setMethod(
  f = "smooth_loess",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, span = 0.75, ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, span = span, ...)
  }
)

# Savitzky-Golay filter ========================================================
#' @export
#' @rdname smooth_savitzky
#' @aliases smooth_savitzky,numeric,numeric-method
setMethod(
  f = "smooth_savitzky",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, m = 3, p = 2) {
    ## Validation
    assert_length(y, length(x))
    assert_odd(m)

    k <- (m - 1) / 2
    i <- seq(from = -k, to = k, by = 1)
    j <- utils::head(utils::tail(seq_along(y), n = -k), n = -k)
    conv <- coef_savitzky(m, p)

    smoothed <- vapply(
      X = j,
      FUN = function(j, i, conv, data) {
        sum(conv * data[j + i])
      },
      FUN.VALUE = double(1),
      i = i,
      conv = conv,
      data = y
    )
    y[j] <- smoothed

    xy <- list(x = x, y = y)
    attr(xy, "method") <- "Savitzky-Golay filter"
    xy
  }
)

#' @export
#' @rdname smooth_savitzky
#' @aliases smooth_savitzky,ANY,missing-method
setMethod(
  f = "smooth_savitzky",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, m, p) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, m = m, p = p)
  }
)

coef_savitzky <- function(m, p = 2) {
  k <- (m - 1) / 2
  z <- seq(from = -k, to = k, by = 1)
  J <- vapply(X = c(0, p), FUN = function(p, z) z^p, z, FUN.VALUE = double(m))
  (solve(t(J) %*% J) %*% t(J))[1, , drop = TRUE]
}

# Whittaker smoothing ==========================================================
#' @export
#' @rdname smooth_whittaker
#' @aliases smooth_whittaker,numeric,numeric-method
setMethod(
  f = "smooth_whittaker",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, lambda = 1600, d = 2, sparse = FALSE) {
    assert_length(y, length(x))
    m <- length(y)

    if (sparse) {
      assert_Matrix()

      E <- Matrix::Diagonal(m)
      D <- Matrix::diff(E, lag = 1, differences = d)
      B <- Matrix::chol(E + (lambda * Matrix::t(D) %*% D))
      z <- Matrix::solve(B, Matrix::solve(Matrix::t(B), y))
      ## Prior to v1.6, Matrix::solve(a=<Matrix>, b=<vector>) returns a matrix
      z <- as.numeric(z)
    } else {
      E <- diag(m)
      D <- diff(E, lag = 1, differences = d)
      B <- E + (lambda * t(D) %*% D)
      z <- solve(B, y)
    }

    xy <- list(x = x, y = z)
    attr(xy, "method") <- "Whittaker smoothing"
    xy
  }
)

#' @export
#' @rdname smooth_whittaker
#' @aliases smooth_whittaker,ANY,missing-method
setMethod(
  f = "smooth_whittaker",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, lambda = 1600, d = 2, sparse = FALSE) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, lambda = lambda, d = d,
                         sparse = sparse)
  }
)

# Penalized likelihood smoothing ===============================================
#' @export
#' @rdname smooth_likelihood
#' @aliases smooth_likelihood,numeric,numeric-method
setMethod(
  f = "smooth_likelihood",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, lambda, d = 2, SE = FALSE,
                        progress = interactive()) {
    assert_Matrix()
    assert_length(y, length(x))

    m <- length(y)
    n <- length(lambda)

    aic <- rep(0, n)
    mu <- matrix(data = 0, nrow = m, ncol = n)
    se <- matrix(data = 0, nrow = m, ncol = n)
    res <- matrix(data = 0, nrow = m, ncol = n)

    if (progress) pb <- utils::txtProgressBar(min = 0, max = n, style=3)

    n_lambda <- seq_len(n)
    for (i in n_lambda) {
      like <- penalized_likelihood(y, lambda = lambda[i], d = d, SE = SE)

      aic[i] <- like$aic
      mu[, i] <- like$mu
      se[, i] <- like$se
      res[, i] <- like$res

      if (progress) utils::setTxtProgressBar(pb, i)
    }

    if (progress) close(pb)

    k <- which.min(aic)
    mu <- mu[, k]
    se <- se[, k]
    res <- res[, k]
    op_lambda <- lambda[k]
    op_aic <- aic[k]

    xy <- list(x = x, y = mu, aic = aic, se = se, residuals = res,
               optimal_aic = op_aic, optimal_lambda = op_lambda)
    attr(xy, "method") <- "penalized likelihood smoothing"
    xy
  }
)

#' @export
#' @rdname smooth_likelihood
#' @aliases smooth_likelihood,ANY,missing-method
setMethod(
  f = "smooth_likelihood",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, lambda, d = 2, SE = FALSE,
                        progress = interactive()) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, lambda = lambda, d = d,
                         SE = SE, progress = progress)
  }
)

penalized_likelihood <- function(y, lambda, d = 2, SE = FALSE) {

  m <- length(y)
  z <- log(y + 1)
  E <- Matrix::Diagonal(m)
  D <- Matrix::diff(E, differences = d)
  P <- lambda * Matrix::crossprod(D, D)
  mu <- exp(z)

  delta_mu <- 1
  limit <- 1 # Stop condition to prevent infinite loop
  while (delta_mu >= 0.0001 && limit <= 10) {
    W <- Matrix::Diagonal(x = mu)
    z <- Matrix::solve(W + P, y - mu + mu * z)
    new_mu <- exp(as.numeric(z)) # Drop S4 class
    delta_mu <- max(abs(new_mu - mu))
    mu <- new_mu

    limit <- limit + 1
  }

  ## AIC
  H <- Matrix::solve(W + P) %*% W
  ed <- sum(Matrix::diag(H))
  dev <- 2 * sum(y * log((y + 1e-10) / mu))
  aic <- dev + 2 * ed * m / (m - ed)

  ## Standard error
  se <- rep(NA_real_, m)
  if (isTRUE(SE)) {
    HH <- Matrix::tcrossprod(H %*% W, H)
    se <- sqrt(Matrix::diag(HH))
  }

  ## Standardized residuals for counts
  res <- (y - mu) / sqrt(mu)

  list(mu = mu, aic = aic, se = se, res = res)
}

# SLIDING WINDOWS
#' @include AllGenerics.R
NULL

#' @export
#' @rdname window_sliding
#' @aliases window_sliding,integer,integer-method
setMethod(
  f = "window_sliding",
  signature = signature(n = "integer", m = "integer"),
  definition = function(n, m, i = NULL) {
    ## Validation
    assert_odd(m)

    ## Index
    k <- (m - 1) / 2
    x <- if (is.null(i)) seq_len(n) else i

    left <- x - k
    left <- ifelse(left > 0, left, 1)

    right <- x + k
    right <- ifelse(right < n, right, n)

    mapply(
      FUN = seq,
      from = left,
      to = right,
      MoreArgs = list(by = 1),
      SIMPLIFY = FALSE
    )
  }
)

#' @export
#' @rdname window_sliding
#' @aliases window_sliding,numeric,numeric-method
setMethod(
  f = "window_sliding",
  signature = signature(n = "numeric", m = "numeric"),
  definition = function(n, m, i = NULL) {
    n <- as.integer(n)
    m <- as.integer(m)
    methods::callGeneric(n, m, i = i)
  }
)

#' @export
#' @rdname window_tumbling
#' @aliases window_tumbling,integer,integer-method
setMethod(
  f = "window_tumbling",
  signature = signature(n = "integer", m = "integer"),
  definition = function(n, m, drop = FALSE) {
    k <- n %% m
    if (drop && k != 0) {
      n <- n - k
    }

    x <- seq_len(n)
    i <- cut(x, breaks = n / m, labels = FALSE)
    split(x, f = i)
  }
)

#' @export
#' @rdname window_tumbling
#' @aliases window_tumbling,numeric,numeric-method
setMethod(
  f = "window_tumbling",
  signature = signature(n = "numeric", m = "numeric"),
  definition = function(n, m, drop = FALSE) {
    n <- as.integer(n)
    m <- as.integer(m)
    methods::callGeneric(n, m, drop = drop)
  }
)

# BASELINE
#' @include AllGenerics.R
NULL

# Linear =======================================================================
#' @export
#' @rdname baseline_linear
#' @aliases baseline_linear,numeric,numeric-method
setMethod(
  f = "baseline_linear",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, points = range(x)) {
    assert_length(y, length(x))

    ## Find the nearest value
    z <- vapply(X = points, FUN = function(i, x) which_nearest(x, i),
                FUN.VALUE = numeric(1), x = x)
    zi <- which(x >= min(points) & x <= max(points))

    data <- data.frame(x, y)
    fit <- stats::lm(y ~ x, data = data, subset = z)
    pred <- stats::predict(fit, data.frame(x = x[zi]))

    bsl <- rep(NA_real_, length(x))
    bsl[zi] <- pred

    xy <- list(x = x, y = bsl)
    attr(xy, "method") <- "linear baseline"
    xy
  }
)

#' @export
#' @rdname baseline_linear
#' @aliases baseline_linear,ANY,missing-method
setMethod(
  f = "baseline_linear",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, points = range(x)) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, points = points)
  }
)

# Polynomial ===================================================================
#' @export
#' @rdname baseline_polynomial
#' @aliases baseline_polynomial,numeric,numeric-method
setMethod(
  f = "baseline_polynomial",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, d = 3, tolerance = 0.001, stop = 100) {
    assert_length(y, length(x))

    polynom <- cbind(1 / sqrt(length(x)), stats::poly(x, degree = d))

    old_y <- y
    start <- 0
    convergence <- FALSE
    while (!convergence) {
      new_y <- polynom %*% crossprod(polynom, old_y)
      new_y <- pmin(y, new_y)

      criterion <- sum(abs((new_y - old_y) / old_y), na.rm = TRUE)
      convergence <- criterion < tolerance

      old_y <- new_y
      start <- start + 1
      if (start >= stop) {
        warning(tr_("Convergence not reached."), call. = FALSE)
        break
      }
    }

    xy <- list(x = x, y = new_y)
    attr(xy, "method") <- "polynomial baseline"
    xy
  }
)

#' @export
#' @rdname baseline_polynomial
#' @aliases baseline_polynomial,ANY,missing-method
setMethod(
  f = "baseline_polynomial",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, d = 3, tolerance = 0.001, stop = 100) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, d = d, tolerance = tolerance,
                         stop = stop)
  }
)

# Rolling Ball =================================================================
#' @export
#' @rdname baseline_rollingball
#' @aliases baseline_rollingball,numeric,numeric-method
setMethod(
  f = "baseline_rollingball",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, m, s) {
    n <- length(x)
    assert_length(y, n)

    ## Windows
    win_minmax <- window_sliding(n, m)
    win_smooth <- window_sliding(n, s)

    ## Minimize
    T1 <- vapply(
      X = win_minmax,
      FUN = function(i, data) {
        min(data[i])
      },
      FUN.VALUE = numeric(1),
      data = y
    )

    ## Maximize
    T2 <- vapply(
      X = win_minmax,
      FUN = function(i, data) {
        max(data[i])
      },
      FUN.VALUE = numeric(1),
      data = T1
    )

    ## Smooth
    T3 <- vapply(
      X = win_smooth,
      FUN = function(i, data) {
        mean(data[i])
      },
      FUN.VALUE = numeric(1),
      data = T2
    )

    xy <- list(x = x, y = T3)
    attr(xy, "method") <- "rolling ball baseline"
    xy
  }
)

#' @export
#' @rdname baseline_rollingball
#' @aliases baseline_rollingball,ANY,missing-method
setMethod(
  f = "baseline_rollingball",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, m, s) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, m = m, s = s)
  }
)

# Rubberband ===================================================================
#' @export
#' @rdname baseline_rubberband
#' @aliases baseline_rubberband,numeric,numeric-method
setMethod(
  f = "baseline_rubberband",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, noise = 0, spline = TRUE, ...) {
    assert_length(y, length(x))

    ## (chull returns points in clockwise order)
    pts <- grDevices::chull(x, y)

    ## Check that ncol(y) is a position 1,
    ## if not, rotate chull vertex index so that ncol(y) is at position 1
    ## then keep only index from ncol(y) to 1 (i.e. lower part of the hull)
    v_max <- which.max(pts) - 1
    if (v_max > 0) pts <- c(pts[-seq_len(v_max)], pts[seq_len(v_max)])
    pts <- pts[1:which.min(pts)]
    # First and last point must be minima, if not remove them
    if (pts[2] == pts[1] - 1) pts <- pts[-1] # Last point
    pts <- rev(pts) # Sort in ascending order
    if (pts[2] == pts[1] + 1) pts <- pts[-1] # First point

    bsl <- stats::approx(x = x[pts], y = y[pts], xout = x, method = "linear")$y
    if (spline) {
      pts <- which(y <= bsl + noise)
      if (length(pts) > 3) {
        spl <- stats::smooth.spline(x[pts], y[pts], ...)
        tmp <- stats::predict(spl, x, 0)$y
      } else {
        tmp <- stats::spline(x[pts], y[pts], xout = x)$y
      }
    }

    ## Check baseline
    if (anyNA(tmp)) {
      msg <- tr_("Failed to estimate the baseline, please check your parameters.")
      stop(msg, call. = FALSE)
    }

    xy <- list(x = x, y = bsl)
    attr(xy, "method") <- "rubberband baseline"
    xy
  }
)

#' @export
#' @rdname baseline_rubberband
#' @aliases baseline_rubberband,ANY,missing-method
setMethod(
  f = "baseline_rubberband",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, noise = 0, spline = TRUE, ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, noise = noise, spline = spline, ...)
  }
)

# SNIP =========================================================================
#' @export
#' @rdname baseline_snip
#' @aliases baseline_snip,numeric,numeric-method
setMethod(
  f = "baseline_snip",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, LLS = FALSE, decreasing = FALSE, n = 100) {
    assert_length(y, length(x))

    ## LLS operator
    y <- if (LLS) LLS(y) else y

    k <- length(y)
    iter <- if (decreasing) rev(seq_len(n)) else seq_len(n)

    tmp <- y
    for (p in iter) {
      for (i in p:(k - p)) {
        a <- y[i]
        b <- (y[i - p] + y[i + p]) / 2
        tmp[i] <- min(a, b)
      }
      y <- tmp
    }

    ## Inverse LLS operator
    bsl <- if (LLS) inverseLLS(y) else y

    ## Check baseline
    if (anyNA(bsl)) {
      msg <- tr_("Failed to estimate the baseline, please check your parameters.")
      stop(msg, call. = FALSE)
    }

    xy <- list(x = x, y = bsl)
    attr(xy, "method") <- "SNIP"
    xy
  }
)

#' @export
#' @rdname baseline_snip
#' @aliases baseline_snip,ANY,missing-method
setMethod(
  f = "baseline_snip",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, LLS = FALSE, decreasing = FALSE, n = 100) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, LLS = LLS, decreasing = decreasing,
                         n = n)
  }
)

#' LLS operator
#'
#' @param x A [`numeric`] vector.
#' @keywords internal
#' @noRd
LLS <- function(x) {
  log(log(sqrt(x + 1) + 1) + 1)
}
inverseLLS <- function(x) {
  (exp(exp(x) - 1) - 1)^2 - 1
}

# AsLS =========================================================================
#' @export
#' @rdname baseline_asls
#' @aliases baseline_asls,numeric,numeric-method
setMethod(
  f = "baseline_asls",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, p = 0.01, lambda = 10^4, stop = 100) {
    assert_Matrix()
    assert_length(y, length(x))

    m <- length(y)
    E <- Matrix::Diagonal(m)
    D <- Matrix::diff(E, lag = 1, differences = 2)

    w <- rep(1, m) # weights

    start <- 0
    convergence <- FALSE
    while (!convergence) {
      W <- Matrix::Diagonal(x = w)
      C <- Matrix::chol(W + (lambda * Matrix::t(D) %*% D))
      z <- Matrix::solve(C, Matrix::solve(Matrix::t(C), w * y))
      ## Prior to v1.6, Matrix::solve(a=<Matrix>, b=<vector>) returns a matrix
      z <- as.numeric(z)
      w0 <- p * (y > z) + (1 - p) * (y < z)
      if (isTRUE(all.equal(w, w0))) convergence <- TRUE

      w <- w0
      start <- start + 1
      if (start >= stop) {
        warning(tr_("Convergence not reached."), call. = FALSE)
        break
      }
    }

    xy <- list(x = x, y = z)
    attr(xy, "method") <- "AsLS"
    xy
  }
)

#' @export
#' @rdname baseline_asls
#' @aliases baseline_asls,ANY,missing-method
setMethod(
  f = "baseline_asls",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, p = 0.01, lambda = 10^4, stop = 100) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, lambda = lambda, p = p, stop = stop)
  }
)

# Peak Filling =================================================================
#' @export
#' @rdname baseline_peakfilling
#' @aliases baseline_peakfilling,numeric,numeric-method
setMethod(
  f = "baseline_peakfilling",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, n, m, by = 10, lambda = 1600, d = 2, sparse = FALSE) {
    assert_length(y, length(x))

    ## Number of bucket intervals
    by <- length(x) / by

    ## Exponential decrease in interval width
    win <- m
    if (n != 1) {
      i <- log10(m) * (1 - (0:(n - 2)) / (n - 1))
      win <- ceiling((10)^c(i, 0))
    }

    ## Smoothing
    smoothed <- smooth_whittaker(x, y, lambda = lambda, d = d, sparse = sparse)
    x <- smoothed$x
    y <- smoothed$y

    ## Subsampling
    breaks <- cut(seq_along(x), breaks = by, labels = FALSE)
    mid <- unlist(tapply(X = x, INDEX = breaks, FUN = mean, simplify = FALSE))
    bin <- unlist(tapply(X = y, INDEX = breaks, FUN = min, simplify = FALSE))

    ## Suppression
    iter <- seq_len(n)
    for (i in iter) {
      win_0 <- win[i] # Current half window width

      ## Interval cut-off close to edges
      k <- seq(from = 2, to = by - 1, by = 1)
      j <- by - k + 1
      v <- pmin(k - 1, win_0, by - k)

      ## Point-wise iteration to the right
      a <- mapply(function(k, v, bin) { mean(bin[(k - v):(k + v)]) }, k, v,
                  MoreArgs = list(bin = bin))
      bin[k] <- pmin(a, bin[k]) # Baseline suppression

      ## Point-wise iteration to the left
      a <- mapply(function(j, v, bin) { mean(bin[(j - v):(j + v)]) }, j, v,
                  MoreArgs = list(bin = bin))
      bin[j] <- pmin(a, bin[j]) # Baseline suppression
    }

    ## Stretch
    mid[1] <- x[1]
    mid[by] <- x[length(x)]
    bsl <- stats::approx(mid, bin, x)$y

    xy <- list(x = x, y = bsl)
    attr(xy, "method") <- "4S Peak Filling"
    xy
  }
)

#' @export
#' @rdname baseline_peakfilling
#' @aliases baseline_peakfilling,ANY,missing-method
setMethod(
  f = "baseline_peakfilling",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, n, m, by = 10, lambda = 1600, d = 2, sparse = FALSE) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, n = n, m = m, by = by,
                         lambda = lambda, d = d, sparse = sparse)
  }
)

# HELPERS

## https://michaelchirico.github.io/potools/articles/developers.html
tr_ <- function(...) {
  enc2utf8(gettext(paste0(...), domain = "R-alkahest"))
}

#' Find the Nearest Value in a Vector
#'
#' @param x A [`numeric`] vector.
#' @param value A [`numeric`] value.
#' @return An [`integer`].
#' @keywords internal
#' @noRd
which_nearest <- function(x, value) {
  which.min(abs(x - value))
}

#' Check Odd Numbers
#'
#' @param x A [`numeric`] vector.
#' @return
#'  Throws an error, if any, and returns `x` invisibly otherwise.
#' @keywords internal
#' @noRd
assert_odd <- function(x) {
  arg <- deparse(substitute(x))
  if (any(round(x) %% 2 == 0)) {
    msg <- sprintf(tr_("%s must be an odd integer (%g)."), sQuote(arg), x)
    stop(msg, call. = FALSE)
  }
  invisible(x)
}

#' Check Object Length
#'
#' @param x An object to be checked.
#' @param expected An appropriate expected value.
#' @return
#'  Throws an error, if any, and returns `x` invisibly otherwise.
#' @keywords internal
#' @noRd
assert_length <- function(x, expected) {
  arg <- deparse(substitute(x))
  if (length(x) != expected) {
    str <- tr_("%s must be of length %d; not %d.")
    msg <- sprintf(str, sQuote(arg), expected, length(x))
    stop(msg, call. = FALSE)
  }
  invisible(x)
}

assert_Matrix <- function() {
  if (!requireNamespace("Matrix", quietly = TRUE)) {
    msg <- tr_("The Matrix package is required. Please install it.")
    stop(msg, call. = FALSE)
  }
  invisible()
}

# FIND PEAKS
#' @include AllGenerics.R
NULL

#' @export
#' @rdname peaks_find
#' @aliases peaks_find,numeric,numeric-method
setMethod(
  f = "peaks_find",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, method = "MAD", SNR = 2, m = NULL, ...) {
    ## Validation
    assert_length(y, length(x))
    method <- match.arg(method, several.ok = FALSE)
    if (is.null(m)) {
      m <- as.integer(length(x) * 0.05)
      if (m %% 2 == 0) m <- m + 1
    }

    ## Noise threshold
    threshold <- NULL
    if (SNR != 0) {
      noise <- switch (
        method,
        MAD = stats::mad
      )
      threshold <- noise(y, ...) * SNR
      index_noise <- y < threshold
      y[index_noise] <- 0
    }

    ## Windows
    shape <- diff(sign(diff(y, na.pad = FALSE)))
    win <- window_sliding(length(x), m, i = which(shape < 0) + 1L)

    ## Peaks detection
    pks <- vapply(
      X = win,
      FUN = function(w, k, data) {
        i <- length(w) - k # Middle of the window
        p <- if (all(data[w[-i]] <= data[w[i]])) w[i] else 0
        return(p)
      },
      FUN.VALUE = numeric(1),
      k = (m - 1) / 2,
      data = y
    )

    xy <- list(x = x[pks], y = y[pks])
    attr(xy, "noise") <- threshold
    xy
  }
)

#' @export
#' @rdname peaks_find
#' @aliases peaks_find,ANY,missing-method
setMethod(
  f = "peaks_find",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, method = "MAD", SNR = 2, m = NULL, ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, method = method, SNR = SNR,
                         m = m, ...)
  }
)

#' @export
#' @rdname peaks_fwhm
#' @aliases peaks_fwhm,numeric,numeric-method
setMethod(
  f = "peaks_fwhm",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, center) {
    assert_length(y, length(x))

    i <- which_nearest(x, center)
    peak_height <- y[i]
    half_max <- peak_height / 2

    scale_for_roots <- y - half_max
    root_indices <- which(diff(sign(scale_for_roots)) != 0)

    tmp <- sort(c(root_indices, i))
    k <- which(tmp == i)

    root_left <- root_indices[k - 1]
    root_right <- root_indices[k]

    HWHM_left <- x[i] - x[root_left]
    HWHM_right <- x[root_right] - x[i]

    # xy <- list(x = c(x[root_left], x[root_right]), y = c(half_max, half_max))

    FWHM <- 2 * min(c(HWHM_left, HWHM_right))
    return(FWHM)
  }
)

#' @export
#' @rdname peaks_fwhm
#' @aliases peaks_fwhm,ANY,missing-method
setMethod(
  f = "peaks_fwhm",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, center) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, center = center)
  }
)

# SIGNAL PROCESSING
#' @include AllGenerics.R
NULL

# Bind =========================================================================
#' @export
#' @rdname signal_bind
#' @aliases signal_bind,ANY-method
setMethod(
  f = "signal_bind",
  signature = "ANY",
  definition = function(...) {

    signal <- list(...)
    coords <- lapply(X = signal, FUN = grDevices::xy.coords)

    y <- lapply(X = coords, FUN = `[[`, i = "y")
    ny <- lengths(y)
    if (any(ny != mean(ny))) {
      msg <- tr_("All objects must have the same number of data points.")
      stop(msg, call. = FALSE)
    }
    y <- do.call(rbind, y)

    ## Get names
    subst <- substitute(list(...))[-1]
    arg_names <- vapply(X = subst, FUN = deparse, FUN.VALUE = character(1))
    rownames(y) <- arg_names
    y
  }
)

# Mean =========================================================================
#' @export
#' @rdname signal_mean
#' @aliases signal_mean,ANY-method
setMethod(
  f = "signal_mean",
  signature = "ANY",
  definition = function(...) {

    signal <- list(...)
    coords <- lapply(X = signal, FUN = grDevices::xy.coords)

    x <- lapply(X = coords, FUN = `[[`, i = "x")
    nx <- lengths(x)

    y <- lapply(X = coords, FUN = `[[`, i = "y")
    ny <- lengths(y)

    if (any(nx != mean(nx)) || any(ny != mean(ny))) {
      msg <- tr_("All objects must have the same number of data points.")
      stop(msg, call. = FALSE)
    }

    x <- do.call(cbind, x)
    x <- rowMeans(x)
    y <- do.call(cbind, y)
    y <- rowMeans(y)

    xy <- list(x = x, y = y)
    xy
  }
)

# Subset =======================================================================
#' @export
#' @rdname subset
#' @aliases signal_select,numeric,numeric-method
setMethod(
  f = "signal_select",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, from, to) {
    assert_length(y, length(x))

    ## Find the nearest value
    from <- x[which_nearest(x, from)]
    to <- x[which_nearest(x, to)]
    idx <- which(x >= from & x <= to)

    xy <- list(x = x[idx], y = y[idx])
    xy
  }
)

#' @export
#' @rdname subset
#' @aliases signal_select,ANY,missing-method
setMethod(
  f = "signal_select",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, from, to) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, from = from, to = to)
  }
)

#' @export
#' @rdname subset
#' @aliases signal_slice,numeric,numeric-method
setMethod(
  f = "signal_slice",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, subset) {
    assert_length(y, length(x))

    i <- as.integer(subset)

    if (!all(i > 0) && !all(i < 0)) {
      msg <- tr_("A vector of strictly positive or negative integers is expected.")
      stop(msg, call. = FALSE)
    }

    xy <- list(x = x[i], y = y[i])
    xy
  }
)

#' @export
#' @rdname subset
#' @aliases signal_slice,ANY,missing-method
setMethod(
  f = "signal_slice",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, subset) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, subset = subset)
  }
)

# Shift ========================================================================
#' @export
#' @rdname signal_shift
#' @aliases signal_shift,numeric,numeric-method
setMethod(
  f = "signal_shift",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, lag) {
    assert_length(y, length(x))

    x <- x + lag
    xy <- list(x = x, y = y)
    xy
  }
)

#' @export
#' @rdname signal_shift
#' @aliases signal_shift,ANY,missing-method
setMethod(
  f = "signal_shift",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, lag) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, lag = lag)
  }
)

# Drift ========================================================================
#' @export
#' @rdname signal_drift
#' @aliases signal_drift,numeric,numeric,numeric-method
setMethod(
  f = "signal_drift",
  signature = c(x = "numeric", y = "numeric", lag = "numeric"),
  definition = function(x, y, lag) {
    assert_length(y, length(x))

    y <- y + lag
    xy <- list(x = x, y = y)
    xy
  }
)

#' @export
#' @rdname signal_drift
#' @aliases signal_drift,ANY,missing,ANY-method
setMethod(
  f = "signal_drift",
  signature = c(x = "ANY", y = "missing", lag = "ANY"),
  definition = function(x, lag, subtract = FALSE) {
    xy <- grDevices::xy.coords(x)
    zz <- grDevices::xy.coords(lag)
    lag <- if (subtract) -zz$y else zz$y
    methods::callGeneric(x = xy$x, y = xy$y, lag = lag)
  }
)

# Correct ======================================================================
#' @export
#' @rdname signal_correct
#' @aliases signal_correct,numeric,numeric-method
setMethod(
  f = "signal_correct",
  signature = c(x = "numeric", y = "numeric"),
  definition = function(x, y, method = c("linear", "polynomial", "asls",
                                         "rollingball", "rubberband",
                                         "SNIP", "4S"), ...) {
    ## Validation
    method <- match.arg(method, several.ok = FALSE)

    ## Get method
    f <- switch(
      method,
      asls = baseline_asls,
      linear = baseline_linear,
      polynomial = baseline_polynomial,
      rollingball = baseline_rollingball,
      rubberband = baseline_rubberband,
      SNIP = baseline_snip,
      `4S` = baseline_peakfilling
    )

    ## Baseline estimation
    bsl <- f(x = x, y = y, ...)

    xy <- list(x = x, y = y - bsl$y)
    attr(xy, "method") <- method
    xy
  }
)

#' @export
#' @rdname signal_correct
#' @aliases signal_correct,ANY,missing-method
setMethod(
  f = "signal_correct",
  signature = c(x = "ANY", y = "missing"),
  definition = function(x, method = c("linear", "polynomial", "asls",
                                      "rollingball", "rubberband",
                                      "SNIP", "4S"), ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, method = method, ...)
  }
)

# XRD
#' @include AllGenerics.R
NULL

#' @export
#' @rdname ka2_strip_penalized
#' @aliases ka2_strip_penalized,numeric,numeric-method
setMethod(
  f = "ka2_strip_penalized",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, lambda, wave = c(1.54060, 1.54443), tau = 0.5,
                        nseg = 1, progress = interactive()) {
    assert_Matrix()
    assert_length(y, length(x))

    ## Calculate doublet distance
    wl1 <- wave[1]
    wl2 <- wave[2]
    wr <- wl2 / wl1
    d2r <- pi / 180
    u <- sin(d2r * x / 2) / wl1
    delta <- 2 * asin(wr * sin(d2r * x / 2)) / d2r - x

    m <- length(y)
    n <- length(lambda)
    if (nseg == 1) nseg <- signif(m * 0.50, 1)

    ## Empty matrices to store the results
    aic <- phi <- psi <- rep(0, n)
    Yhat <- Yhat2 <- matrix(data = 0, nrow = m, ncol = n)
    mu <- matrix(data = 0, nrow = m, ncol = n)
    sr <- matrix(data = 0, nrow = m, ncol = n)

    ## Prepare model matrices
    xr <- max(x)
    xl <- min(x) - 1
    xd <- x - delta
    dx <- (xr - xl) / nseg
    knots <- seq(xl - 3 * dx, xr + 3 * dx, by = dx)

    B1 <- sbase(x, lower = xl, upper = xr, n = nseg, deg = 3)
    B2 <- sbase(xd, lower = xl, upper = xr, n = nseg, deg = 3)

    B <- Matrix::rbind2(B1, B2)
    Cb <- Matrix::Diagonal(m)
    C <- Matrix::cbind2(Cb, tau * Cb)

    ## Prepare penalty
    E <- Matrix::Diagonal(ncol(B1))
    D <- Matrix::diff(E, diff = 2)

    ## Initialize
    bst <- Matrix::solve(
      Matrix::crossprod(B1, B1) + 1e-4 * E,
      Matrix::crossprod(B1, log((y + 1) / 2))
    )

    if (progress) pb <- utils::txtProgressBar(min = 0, max = n, style = 3)

    ## Estimate for all lambda
    n_lambda <- seq_len(n)
    for (i in n_lambda) {
      M <- penalized_strip_ka2(y, B, C, D, bst, lambda[i])
      aic[i] <- M$aic
      Yhat[, i] <- M$yhat
      Yhat2[, i] <- M$yhat2
      sr[,i] <- M$res
      mu[, i] <- M$mu
      phi[i] <- M$phi
      psi[i] <- M$psi
      if (progress) utils::setTxtProgressBar(pb, i)
    }

    if (progress) close(pb)

    ## Find optimal model and report
    k <- which.min(aic)
    mu <- mu[, k]
    op_lamb <- lambda[k]
    op_aic <- aic[k]
    yhat <- Yhat[, k]
    yhat2 <- Yhat2[, k]
    sr <- sr[, k]

    xy <- list(x = x, y = yhat, ka2 = yhat2, smoothed = mu)
    xy
  }
)

#' @export
#' @rdname ka2_strip_penalized
#' @aliases ka2_strip_penalized,ANY,missing-method
setMethod(
  f = "ka2_strip_penalized",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, lambda, wave = c(1.54060, 1.54443), tau = 0.5,
                        nseg = 1, progress = interactive()) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, lambda = lambda, wave = wave,
                         tau = tau, nseg = nseg, progress = progress)
  }
)


#' Truncated p-th Power
#'
#' @param x A [`numeric`] vector on which the basis is calculated.
#' @param knot A [`numeric`] scalar giving the truncation point.
#' @param p A [`numeric`] scalar specifying the power for the basis
#'  (e.g. \eqn{p = 3} for cubic TPF).
#' @return
#'  A [`numeric`] vector (piece-wise defined basis functions for \eqn{x > t}).
#' @author P. Eilers
#' @note
#'  Slightly modified from P. Eilers' [JOPS::tpower()].
#' @references
#'  Eilers, P. H. C. and Marx, B. D. (2021). *Practical Smoothing, The Joys of
#'  P-splines.* Cambridge: Cambridge University Press.
#' @keywords internal
#' @noRd
tpower <- function(x, knot, p) {
  (x - knot) ^ p * (x > knot)
}

#' B-spline Design Matrix
#'
#' Computes a B-spline design matrix using evenly spaced knots.
#' @param x A [`numeric`] vector of values at which the B-spline basis functions
#'  are to be evaluated.
#' @param lower A [`numeric`] vector giving the lower limit of the domain of `x`.
#' @param upper A [`numeric`] vector giving the upper limit of the domain of `x`.
#' @param n A [`numeric`] scalar specifying the number of equally sized segments
#'  between `lower` and `upper`.
#' @param deg A [`numeric`] scalar specifyingthe degree of the splines,
#'  usually 1, 2, or 3 (the default).
#' @return
#'  A `matrix` with  `length(x)` rows and `n + deg` columns.
#' @note
#'  Slightly modified from P. Eilers and B. Marx's [JOPS::bbase()].
#' @author P. Eilers, B. Marx
#' @references
#'  Eilers, P. H. C. and Marx, B. D. (2021). *Practical Smoothing, The Joys of
#'  P-splines.* Cambridge: Cambridge University Press.
#' @keywords internal
#' @noRd
bbase <- function(x, lower = min(x), upper = max(x), n = 10, deg = 3) {
  ## Check domain and adjust it if necessary
  if (lower > min(x)) {
    lower <- min(x)
    warning(sprintf(tr_("Left boundary adjusted to %g."), lower), call. = FALSE)
  }
  if (upper <  max(x)) {
    upper <- max(x)
    warning(sprintf(tr_("Right boundary adjusted to %g."), upper), call. = FALSE)
  }

  ## Function for B-spline basis
  dx <- (upper - lower) / n
  knots <- seq(lower - deg * dx, upper + deg * dx, by = dx)
  P <- outer(x, knots, tpower, deg)
  n <- dim(P)[[2L]]
  D <- diff(diag(n), diff = deg + 1) / (gamma(deg + 1) * dx ^ deg)
  B <- (-1) ^ (deg + 1) * P %*% t(D)

  ## Make B-splines exactly zero beyond their end knots
  nb <- ncol(B)
  sk <- knots[(1:nb) + deg + 1]
  mask <- outer(x, sk, "<")
  B <- B * mask

  B
}

#' Sparse B-spline Basis Matrix
#'
#' Computes a B-spline basis matrix using evenly spaced knots.
#' @param x A [`numeric`] vector of values at which the B-spline basis functions
#'  are to be evaluated.
#' @param lower A [`numeric`] vector giving the lower limit of the domain of `x`.
#' @param upper A [`numeric`] vector giving the upper limit of the domain of `x`.
#' @param n A [`numeric`] scalar specifying the number of equally sized segments
#'  between `lower` and `upper`.
#' @param deg A [`numeric`] scalar specifyingthe degree of the splines,
#'  usually 1, 2, or 3 (the default).
#' @return
#'  A [sparse matrix][Matrix::sparseMatrix()] with `length(x)` rows and
#'  `n + deg` columns.
#' @note
#'  Slightly modified from P. Eilers [JOPS::spbase()].
#' @author P. Eilers, B. Marx
#' @references
#'  Eilers, P. H. C. and Marx, B. D. (2021). *Practical Smoothing, The Joys of
#'  P-splines.* Cambridge: Cambridge University Press.
#' @keywords internal
#' @noRd
sbase <- function(x, lower = min(x), upper = max(x), n = 10, deg = 3) {
  ## Check domain and adjust it if necessary
  if (lower > min(x)) {
    lower <- min(x)
    warning(sprintf(tr_("Left boundary adjusted to %g."), lower), call. = FALSE)
  }
  if (upper <  max(x)) {
    upper <- max(x)
    warning(sprintf(tr_("Right boundary adjusted to %g."), upper), call. = FALSE)
  }

  ## Reduce x to first interval between knots
  dx <- (upper - lower) / n
  ix <- floor((x - lower) / (1.0000001 * dx))
  xx <- (x - lower) - ix * dx

  ## Full basis for reduced x
  Br <- bbase(xx, lower = 0, upper = 0 + dx, n = 1, deg = deg)

  ## Compute proper rows, columns
  m <- length(x)
  p <- ncol(Br)
  i_row <- rep(seq_len(m), each = p)
  i_col <- rep(seq_len(p), times = m) + rep(ix, each = p)

  ## Sparse matrix
  b <- as.vector(t(Br))
  Bs <- Matrix::sparseMatrix(i = i_row, j = i_col, x = b, dims = c(m, n + deg))
  Bs <- Bs[, -ncol(Bs), drop = FALSE] # ???

  Bs
}

#' @param y description
#' @param B description
#' @param C description
#' @param D description
#' @param beta description
#' @param lambda description
#' @author J. J. de Rooi *et al.* (original R code).
#' @references
#'  de Rooi, J. J., van der Pers, N. M., Hendrikx, R. W. A., Delhez, R.,
#'  Böttger A. J. and Eilers, P. H. C. (2014). Smoothing of X-ray diffraction
#'  data and Ka2 elimination using penalized likelihood and the composite link
#'  model. *Journal of Applied Crystallography*, 47: 852-860.
#'  \doi{10.1107/S1600576714005809}
#' @keywords internal
#' @noRd
penalized_strip_ka2 <- function(y, B, C, D, beta, lambda) {
  m <- length(y)

  ## Penalty
  P <- lambda * Matrix::t(D) %*% D

  ## Fit iteratively
  delta_beta <- 1
  limit <- 1 # Stop condition to prevent infinite loop
  while (delta_beta >= 0.0001 && limit <= 30) {
    eta <- B %*% beta
    gam <- exp(eta)
    mu <- as.numeric(C %*% gam)

    M <- Matrix::Diagonal(x = 1 / mu)
    X <- (M %*% C %*% Matrix::Diagonal(x = as.numeric(gam))) %*% B
    W <- Matrix::Diagonal(x = mu)
    G <- Matrix::crossprod(X, W) %*% X

    new_beta <- Matrix::solve(G + P, Matrix::crossprod(X, (y - mu)) + G %*% beta)
    delta_beta <- max(abs(new_beta - beta))
    beta <- new_beta

    limit <- limit + 1
  }

  ## AIC
  H <- Matrix::solve(G + P, G)
  ed <- sum(Matrix::diag(H))
  dev <- 2 * sum(y * log((y + (y == 0)) / mu))
  aic <- dev + 2 * ed

  ## Standardized residuals for counts
  sr <- (y - mu) / sqrt(mu)
  sdr <- sqrt(sum(sr^2) / (m - ed))

  ## L-curve
  phi <- log(sum((D %*% beta)^2))
  psi <- log(sum(sr^2))

  ## estimated signals
  yhat <- gam[1:m]
  yhat2 <- gam[(m + 1):length(gam)] * 0.5

  ## return values
  list(yhat = yhat, yhat2 = yhat2, mu = mu, aic = aic, res = sr,
       beta = beta, phi = phi, psi = psi)
}

# REPLACE
#' @include AllGenerics.R
NULL

# Threshold ====================================================================
#' @export
#' @rdname replace_threshold
#' @aliases replace_threshold,numeric,numeric,function-method
setMethod(
  f = "replace_threshold",
  signature = signature(x = "numeric", y = "numeric", threshold = "function"),
  definition = function(x, y, threshold, value = 0, ...) {
    assert_length(y, length(x))

    threshold <- threshold(y, ...)
    methods::callGeneric(x, y, threshold = threshold, value = value)
  }
)

#' @export
#' @rdname replace_threshold
#' @aliases replace_threshold,ANY,missing,function-method
setMethod(
  f = "replace_threshold",
  signature = signature(x = "ANY", y = "missing", threshold = "function"),
  definition = function(x, threshold, value = 0, ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, threshold = threshold,
                         value = value)
  }
)

#' @export
#' @rdname replace_threshold
#' @aliases replace_threshold,numeric,numeric,numeric-method
setMethod(
  f = "replace_threshold",
  signature = signature(x = "numeric", y = "numeric", threshold = "numeric"),
  definition = function(x, y, threshold, value = 0, ...) {
    assert_length(y, length(x))

    i <- y < threshold
    y[i] <- value

    xy <- list(x = x, y = y)
    xy
  }
)

#' @export
#' @rdname replace_threshold
#' @aliases replace_threshold,ANY,missing,numeric-method
setMethod(
  f = "replace_threshold",
  signature = signature(x = "ANY", y = "missing", threshold = "numeric"),
  definition = function(x, threshold, value = 0, ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, threshold = threshold,
                         value = value)
  }
)

# Negative =====================================================================
#' @export
#' @rdname replace_negative
#' @aliases replace_negative,numeric,numeric-method
setMethod(
  f = "replace_negative",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, value = 0) {
    replace_threshold(x, y, threshold = 0, value = value)
  }
)

#' @export
#' @rdname replace_negative
#' @aliases replace_negative,ANY,missing-method
setMethod(
  f = "replace_negative",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, value = 0) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, value = value)
  }
)

# SCALING
#' @include AllGenerics.R
NULL

# Area =========================================================================
#' @export
#' @rdname rescale_area
#' @aliases rescale_area,numeric,numeric-method
setMethod(
  f = "rescale_area",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, method = c("rectangle", "trapezoid"), ...) {
    ## Validation
    assert_length(y, length(x))
    method <- match.arg(method, several.ok = FALSE)

    ## Get method
    f <- switch(
      method,
      rectangle = integrate_rectangle,
      trapezoid = integrate_trapezoid
    )

    ## Baseline estimation
    auc <- f(x = x, y = y, ...)
    y <- y / auc

    xy <- list(x = x, y = y)
    xy
  }
)

#' @export
#' @rdname rescale_area
#' @aliases rescale_area,ANY,missing-method
setMethod(
  f = "rescale_area",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, method = c("rectangle", "trapezoid"), ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, method = method, ...)
  }
)

# Total ========================================================================
#' @export
#' @rdname rescale_total
#' @aliases rescale_total,numeric,numeric-method
setMethod(
  f = "rescale_total",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, total = 1) {
    assert_length(y, length(x))

    y <- (y * total) / sum(abs(y))
    xy <- list(x = x, y = y)
    xy
  }
)

#' @export
#' @rdname rescale_total
#' @aliases rescale_total,ANY,missing-method
setMethod(
  f = "rescale_total",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, total = 1) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, total = total)
  }
)

# Min-Max ======================================================================
#' @export
#' @rdname rescale_range
#' @aliases rescale_range,numeric,numeric-method
setMethod(
  f = "rescale_range",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, min = 0, max = 1) {
    ## Validation
    assert_length(y, length(x))
    if (min > max) {
      msg <- sprintf(tr_("%s (%g) must be lower than %s (%g)."),
                     sQuote("min"), min, sQuote("max"), max)
      stop(msg, call. = FALSE)
    }

    y <- (y - min(y)) / (max(y) - min(y)) * (max - min) + min
    xy <- list(x = x, y = y)
    xy
  }
)

#' @export
#' @rdname rescale_range
#' @aliases rescale_range,ANY,missing-method
setMethod(
  f = "rescale_range",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, min = 0, max = 1) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, min = min, max = max)
  }
)

#' @export
#' @rdname rescale_range
#' @aliases rescale_min,numeric,numeric-method
setMethod(
  f = "rescale_min",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, min = 0) {
    rescale_range(x, y, min = min, max = max(y))
  }
)

#' @export
#' @rdname rescale_range
#' @aliases rescale_min,ANY,missing-method
setMethod(
  f = "rescale_min",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, min = 0) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, min = min)
  }
)

#' @export
#' @rdname rescale_range
#' @aliases rescale_max,numeric,numeric-method
setMethod(
  f = "rescale_max",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, max = 1) {
    rescale_range(x, y, min = min(y), max = max)
  }
)

#' @export
#' @rdname rescale_range
#' @aliases rescale_max,ANY,missing-method
setMethod(
  f = "rescale_max",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, max = 1) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, max = max)
  }
)

# Transform ====================================================================
#' @export
#' @rdname rescale_transform
#' @aliases rescale_transform,numeric,numeric-method
setMethod(
  f = "rescale_transform",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, f, ...) {
    assert_length(y, length(x))

    y <- f(y, ...)
    xy <- list(x = x, y = y)
    xy
  }
)

#' @export
#' @rdname rescale_transform
#' @aliases rescale_transform,ANY,missing-method
setMethod(
  f = "rescale_transform",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, f, ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, f = f, ...)
  }
)

# SNV ==========================================================================
#' @export
#' @rdname rescale_snv
#' @aliases rescale_snv,numeric,numeric-method
setMethod(
  f = "rescale_snv",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, ...) {
    ## Validation
    assert_length(y, length(x))

    y <- (y - mean(y)) / stats::sd(y)

    xy <- list(x = x, y = y)
    xy
  }
)

#' @export
#' @rdname rescale_snv
#' @aliases rescale_snv,ANY,missing-method
setMethod(
  f = "rescale_snv",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, ...)
  }
)

# GENERIC METHODS

# Baseline =====================================================================
#' Asymmetric Least Squares Smoothing
#'
#' Baseline estimation with asymmetric least squares smoothing.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param p A length-one [`numeric`] vector giving the asymmetry
#'  (\eqn{0.001 < p < 0.1} is a good choice for a signal with positive peaks).
#' @param lambda A length-one [`numeric`] vector giving the smoothing parameter.
#' @param stop An [`integer`] giving the stopping rule (i.e. maximum number of
#'  iterations).
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @references
#'  Eilers, P. H. C. & Boelens, H. F. M. (2005). *Baseline Correction with
#'  Asymmetric Least Squares Smoothing*.
#' @seealso [signal_correct()]
#' @example inst/examples/ex-baseline-asls.R
#' @author P. H. C. Eilers and H. F. M. Boelens (original Matlab code)
#' @docType methods
#' @family baseline estimation methods
#' @aliases baseline_asls-method
setGeneric(
  name = "baseline_asls",
  def = function(x, y, ...) standardGeneric("baseline_asls")
)

#' Linear Baseline Estimation
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param points A [`numeric`] vector specifying the data points to be used in
#'  the fitting process (in `x` unit).
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @seealso [signal_correct()]
#' @example inst/examples/ex-baseline-linear.R
#' @author N. Frerebeau
#' @docType methods
#' @family baseline estimation methods
#' @aliases baseline_linear-method
setGeneric(
  name = "baseline_linear",
  def = function(x, y, ...) standardGeneric("baseline_linear")
)

#' Polynomial Baseline Estimation
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param d An [`integer`] giving the degree of the polynomial. Must be less
#'  than the number of unique points.
#' @param tolerance A [`numeric`] scalar giving the tolerance of difference
#'  between iterations.
#' @param stop An [`integer`] giving the stopping rule (i.e. maximum number of
#'  iterations).
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @seealso [signal_correct()]
#' @references
#'  Lieber, C. A. and Mahadevan-Jansen, A. (2003). Automated Method for
#'  Subtraction of Fluorescence from Biological Raman Spectra. *Applied
#'  Spectroscopy*, 57(11): 1363-67. \doi{10.1366/000370203322554518}.
#' @example inst/examples/ex-baseline-polynomial.R
#' @author N. Frerebeau
#' @docType methods
#' @family baseline estimation methods
#' @aliases baseline_polynomial-method
setGeneric(
  name = "baseline_polynomial",
  def = function(x, y, ...) standardGeneric("baseline_polynomial")
)

#' Rolling Ball Baseline Estimation
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param m An odd [`integer`] giving the window size (i.e. the number of
#'  adjacent points to be used; see [`window_sliding()`]) for
#'  minimization/maximization.
#' @param s An odd [`integer`] giving the window size (i.e. the number of
#'  adjacent points to be used; see [`window_sliding()`]) for smoothing.
#' @param ... Currently not used.
#' @note
#'  There will be \eqn{(m - 1) / 2} points both at the beginning and at the end
#'  of the data series for which a complete \eqn{m}-width window cannot be
#'  obtained. To prevent data loss, progressively wider/narrower windows are
#'  used at both ends of the data series.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @seealso [signal_correct()]
#' @references
#'  Kneen, M. A. and Annegarn, H. J. (1996). Algorithm for Fitting XRF, SEM and
#'  PIXE X-Ray Spectra Backgrounds. *Nuclear Instruments and Methods in Physics
#'  Research Section B: Beam Interactions with Materials and Atoms*,
#'  109/110: 209-213. \doi{10.1016/0168-583X(95)00908-6}.
#' @example inst/examples/ex-baseline-rollingball.R
#' @author N. Frerebeau
#' @docType methods
#' @family baseline estimation methods
#' @aliases baseline_rollingball-method
setGeneric(
 name = "baseline_rollingball",
 def = function(x, y, ...) standardGeneric("baseline_rollingball")
)

#' Rubberband Baseline Estimation
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param noise A length-one [`numeric`] vector giving the noise level.
#'  Only used if `method` is "`rubberband`".
#' @param spline A [`logical`] scalar: should spline interpolation through the
#'  support points be used instead of linear interpolation?
#'  Only used if `method` is "`rubberband`".
#' @param ... Extra arguments to be passed to [stats::smooth.spline()].
#' @details
#'  A convex envelope of the spectrum is determined and the
#'  baseline is estimated as the part of the convex envelope lying below the
#'  spectrum. Note that the rubber band does not enter the concave regions
#'  (if any) of the spectrum.
#' @note
#'  `baseline_rubberband()` is slightly modified from C. Beleites'
#'  [hyperSpec::spc.rubberband()].
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @seealso [signal_correct()]
#' @example inst/examples/ex-baseline-rubberband.R
#' @author N. Frerebeau
#' @docType methods
#' @family baseline estimation methods
#' @aliases baseline_rubberband-method
setGeneric(
  name = "baseline_rubberband",
  def = function(x, y, ...) standardGeneric("baseline_rubberband")
)

#' SNIP Baseline Estimation
#'
#' Sensitive Nonlinear Iterative Peak clipping algorithm.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param LLS A [`logical`] scalar: should the LLS operator be applied on `x`
#'  before employing SNIP algorithm? Only used if `method` is "`SNIP`".
#' @param decreasing A [`logical`] scalar: should a decreasing clipping window
#'  be used?
#' @param n An [`integer`] value giving the number of iterations.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @seealso [signal_correct()]
#' @references
#'  Morháč, M., Kliman, J., Matoušek, V., Veselský, M. & Turzo, I. (1997).
#'  Background elimination methods for multidimensional gamma-ray spectra.
#'  *Nuclear Instruments and Methods in Physics Research Section A:
#'  Accelerators, Spectrometers, Detectors and Associated Equipment*, 401(1),
#'  p. 113-132. \doi{10.1016/S0168-9002(97)01023-1}
#'
#'  Morháč, M. & Matoušek, V. (2008). Peak Clipping Algorithms for Background
#'  Estimation in Spectroscopic Data. *Applied Spectroscopy*, 62(1), p. 91-106.
#'  \doi{10.1366/000370208783412762}
#'
#'  Ryan, C. G., Clayton, E., Griffin, W. L., Sie, S. H. & Cousens, D. R.
#'  (1988). SNIP, a statistics-sensitive background treatment for the
#'  quantitative analysis of PIXE spectra in geoscience applications.
#'  *Nuclear Instruments and Methods in Physics Research Section B:
#'  Beam Interactions with Materials and Atoms*, 34(3), p. 396-402.
#'  \doi{10.1016/0168-583X(88)90063-8}
#' @example inst/examples/ex-baseline-snip.R
#' @author N. Frerebeau
#' @docType methods
#' @family baseline estimation methods
#' @aliases baseline_snip-method
setGeneric(
  name = "baseline_snip",
  def = function(x, y, ...) standardGeneric("baseline_snip")
)

#' 4S Peak Filling
#'
#' Baseline estimation by iterative mean suppression.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param n An [`integer`] value giving the number of iterations.
#' @param m An odd [`integer`] giving the half window size.
#' @param by A length-one [`numeric`] vector giving the umber of buckets to
#'  divide `x` into.
#' @param lambda An [`integer`] giving the smoothing parameter. The larger
#'  `lambda` is, the smoother the curve (see [smooth_whittaker()]).
#' @param d An [`integer`] specifying the order of the penalty (see
#'  [smooth_whittaker()]).
#' @param sparse A [`logical`] scalar: should sparse matrices be used for
#'  computation (see [smooth_whittaker()])? If `TRUE`, \pkg{Matrix} is required.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @references
#'  Liland, K. H. (2015). 4S Peak Filling - baseline estimation by iterative
#'  mean suppression. *MethodsX*, 2, 135-140. \doi{10.1016/j.mex.2015.02.009}.
#' @seealso [signal_correct()], [smooth_whittaker()]
#' @example inst/examples/ex-baseline-peakfilling.R
#' @author N. Frerebeau
#' @docType methods
#' @family baseline estimation methods
#' @aliases baseline_peakfilling-method
setGeneric(
  name = "baseline_peakfilling",
  def = function(x, y, ...) standardGeneric("baseline_peakfilling")
)

# Integrate ====================================================================
#' Rectangle Rule
#'
#' Approximates the definite integral by using the rectangle rule.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param right A [`logical`] scalar: should the right rule be used instead of
#'  the left rule?
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @example inst/examples/ex-integrate.R
#' @author N. Frerebeau
#' @docType methods
#' @family integration methods
#' @aliases integrate_rectangle-method
setGeneric(
  name = "integrate_rectangle",
  def = function(x, y, ...) standardGeneric("integrate_rectangle")
)

#' Trapezoidal Rule
#'
#' Approximates the definite integral by using the trapezoidal rule.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @example inst/examples/ex-integrate.R
#' @author N. Frerebeau
#' @docType methods
#' @family integration methods
#' @aliases integrate_trapezoid-method
setGeneric(
  name = "integrate_trapezoid",
  def = function(x, y, ...) standardGeneric("integrate_trapezoid")
)

# Peaks ========================================================================
#' Find Peaks
#'
#' Finds local maxima in sequential data.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param method A [`character`] string specifying the method to be used for
#'  background noise estimation (see below).
#' @param SNR An [`integer`] giving the signal-to-noise-ratio for peak detection
#'  (see below).
#' @param m An odd [`integer`] giving the window size (i.e. the number of
#'  adjacent points to be used).
#'  If `NULL`, 5% of the data points is used as the half window size.
#' @param ... Extra parameters to be passed to internal methods.
#' @details
#'  A local maximum has to be the highest one in the given window and has to be
#'  higher than \eqn{SNR \times noise}{SNR * noise} to be recognized as peak.
#'
#'  The following methods are available for noise estimation:
#'  \describe{
#'   \item{`MAD`}{Median Absolute Deviation.}
#'  }
#'
#'  Note that to improve peak detection, it may be helpful to smooth the data
#'  and remove the baseline beforehand.
#' @note
#'  There will be \eqn{(m - 1) / 2} points both at the beginning and at the end
#'  of the data series for which a complete \eqn{m}-width window cannot be
#'  obtained. To prevent data loss, progressively wider/narrower windows are
#'  used at both ends of the data series.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @note
#'  Adapted from Stasia Grinberg's
#'  [`findPeaks`](https://github.com/stas-g/findPeaks) function.
#' @example inst/examples/ex-peaks.R
#' @author N. Frerebeau
#' @docType methods
#' @family peaks detection methods
#' @aliases peaks_find-method
setGeneric(
  name = "peaks_find",
  def = function(x, y, ...) standardGeneric("peaks_find")
)

#' Half-Width at Half-Maximum
#'
#' Estimates the Half-Width at Half-Maximum (FWHM) for a given peak.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param center A [`numeric`] value giving the peak position in `x` units.
#' @param ... Currently not used.
#' @return A [`numeric`] value.
#' @details
#'  It tries to get the smallest possible estimate.
#' @example inst/examples/ex-peaks.R
#' @author N. Frerebeau
#' @docType methods
#' @family peaks detection methods
#' @aliases peaks_fwhm-method
setGeneric(
  name = "peaks_fwhm",
  def = function(x, y, ...) standardGeneric("peaks_fwhm")
)

# Replace ======================================================================
#' Replace Values Below a Given Threshold
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param threshold A [`numeric`] value or a [`function`] that takes a `numeric`
#'  vector as argument and returns a single `numeric` value.
#' @param value A [`numeric`] value to replace values below `threshold`.
#' @param ... Extra parameters to be passed to `threshold`.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @docType methods
#' @family replacement methods
#' @aliases replace_threshold-method
setGeneric(
  name = "replace_threshold",
  def = function(x, y, threshold, ...) standardGeneric("replace_threshold")
)

#' Replace Negative Values
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param value A [`numeric`] value to replace negative values.
#' @param ... Extra parameters to be passed to `threshold`.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @docType methods
#' @family replacement methods
#' @aliases replace_negative-method
setGeneric(
  name = "replace_negative",
  def = function(x, y, ...) standardGeneric("replace_negative")
)

# Resample =====================================================================
#' Bin
#'
#' Averages `x` values and applies a function to the corresponding `y` values.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param by An [`integer`] specifying the binning ratio (i.e. the number of
#'  points to be grouped together; see [`window_tumbling()`]).
#' @param f A [`function`] that takes a `numeric` vector of intensities as
#'  argument and returns a single `numeric` vector. Used to estimate the local
#'  representative value in each bin (defaults to [sum()]; see examples).
#' @param ... Extra parameters to be passed to `f`.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-resample.R
#' @docType methods
#' @family resampling methods
#' @aliases resample_bin-method
setGeneric(
  name = "resample_bin",
  def = function(x, y, ...) standardGeneric("resample_bin")
)

#' Downsample
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param by An [`integer`] specifying the downsampling factor.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-resample.R
#' @docType methods
#' @family resampling methods
#' @aliases resample_down-method
setGeneric(
  name = "resample_down",
  def = function(x, y, ...) standardGeneric("resample_down")
)

# Upsample
#
# @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
# @param by An [`integer`] specifying the upsampling factor.
# @param ... Currently not used.
# @details
#  The signal is upsampled by inserting \eqn{n - 1} zeros between samples.
# @return
#  Returns a [`list`] with two components `x` and `y`.
# @author N. Frerebeau
# @example inst/examples/ex-resample.R
# @docType methods
# @family resampling methods
# @aliases resample_interpolate-method
# setGeneric(
#   name = "resample_up",
#   def = function(x, y, ...) standardGeneric("resample_up")
# )

#' Linearly Interpolate
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param from A length-one [`numeric`] vector giving the starting value of the
#'  sequence where interpolation is to take place.
#' @param to A length-one [`numeric`] vector giving the end value of the
#'  sequence where interpolation is to take place.
#' @param by A length-one [`numeric`] vector specifying the increment of the
#'  sequence.
#' @param ... Extra arguments to be passed to [stats::approx()].
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-resample.R
#' @docType methods
#' @family resampling methods
#' @aliases resample_interpolate-method
setGeneric(
  name = "resample_interpolate",
  def = function(x, y, ...) standardGeneric("resample_interpolate")
)

# Rescale ======================================================================
#' Normalize intensities by AUC
#'
#' Rescales intensities so that the area under the curve (AUC) is equal to 1.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param method A [`character`] string specifying the method for integration.
#'  It must be one of "`rectangle`" or "`trapezoid`".
#'  Any unambiguous substring can be given.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @example inst/examples/ex-rescale.R
#' @author N. Frerebeau
#' @docType methods
#' @family normalization methods
#' @aliases rescale_area-method
setGeneric(
  name = "rescale_area",
  def = function(x, y, ...) standardGeneric("rescale_area")
)

#' Rescale intensities to sum to a specified value
#'
#' Rescales intensities to sum to a specified value.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param total A legnth-one [`numeric`] vector specifying the output total.
#'  Defaults to 1, i.e. normalizes by total intensity.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @example inst/examples/ex-rescale.R
#' @author N. Frerebeau
#' @docType methods
#' @family normalization methods
#' @aliases rescale_total-method
setGeneric(
  name = "rescale_total",
  def = function(x, y, ...) standardGeneric("rescale_total")
)

#' Rescales intensities to have specified minimum and maximum
#'
#' Rescales intensities to have specified minimum and maximum.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param min A legnth-one [`numeric`] vector specifying the output minimum.
#' @param max A legnth-one [`numeric`] vector specifying the output maximum.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @example inst/examples/ex-rescale.R
#' @author N. Frerebeau
#' @docType methods
#' @family normalization methods
#' @aliases rescale_range-method
setGeneric(
  name = "rescale_range",
  def = function(x, y, ...) standardGeneric("rescale_range")
)

#' @rdname rescale_range
#' @aliases rescale_min-method
setGeneric(
  name = "rescale_min",
  def = function(x, y, ...) standardGeneric("rescale_min")
)

#' @rdname rescale_range
#' @aliases rescale_max-method
setGeneric(
  name = "rescale_max",
  def = function(x, y, ...) standardGeneric("rescale_max")
)

#' Transform Intensities
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param f A [`function`] that takes a `numeric` vector of intensities as
#'  argument and returns a `numeric` vector.
#' @param ... Extra arguments to be passed to `f`.
#' @details
#'  Transformation of intensities can be used to improve the identification of
#'  peaks with a low signal-to-noise ratio.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @example inst/examples/ex-transform.R
#' @author N. Frerebeau
#' @docType methods
#' @family normalization methods
#' @aliases rescale_transform-method
setGeneric(
  name = "rescale_transform",
  def = function(x, y, ...) standardGeneric("rescale_transform")
)

#' Standard Normal Variate (SNV) Transformation
#'
#' Subtracts the mean and scales to unit variance.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @references
#'  Barnes, R. J., Dhanoa, M. S. & Lister, S. J. (1989). Standard Normal Variate
#'  Transformation and De-Trending of Near-Infrared Diffuse Reflectance Spectra.
#'  *Applied Spectroscopy*, 43(5): 772-777. \doi{10.1366/0003702894202201}.
#' @example inst/examples/ex-normalize.R
#' @author N. Frerebeau
#' @docType methods
#' @family normalization methods
#' @aliases rescale_snv-method
setGeneric(
  name = "rescale_snv",
  def = function(x, y, ...) standardGeneric("rescale_snv")
)

# Signal =======================================================================
#' Baseline Correction
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param method A [`character`] string specifying the method for baseline
#'  estimation (see details). Any unambiguous substring can be given.
#' @param ... Extra arguments to be passed to `baseline_*()` (see details).
#' @details
#'  Available methods for baseline estimation:
#'  \describe{
#'   \item{`asls`}{Asymmetric Least Squares Smoothing (see [baseline_asls()]).}
#'   \item{`linear`}{Linear baseline estimation (see [baseline_linear()]).}
#'   \item{`polynomial`}{Polynomial baseline estimation (see
#'   [baseline_polynomial()]).}
#'   \item{`rollingball`}{Rolling ball baseline estimation (see
#'   [baseline_rollingball()]).}
#'   \item{`rubberband`}{Rubberband baseline estimation (see
#'   [baseline_rubberband()]).}
#'   \item{`SNIP`}{Sensitive Nonlinear Iterative Peak clipping algorithm
#'   (see [baseline_snip()]).}
#'   \item{`4S`}{4S Peak Filling (see [baseline_peakfilling()]).}
#'  }
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-correct.R
#' @docType methods
#' @family signal processing methods
#' @aliases signal_correct-method
setGeneric(
  name = "signal_correct",
  def = function(x, y, ...) standardGeneric("signal_correct")
)

#' Bind
#'
#' Combines XY objects.
#' @param ... Any object that can be interpreted in a suitable way
#'  (see [grDevices::xy.coords()]).
#' @return
#'  Returns a [`matrix`] of intensities.
#' @author N. Frerebeau
#' @example inst/examples/ex-mean.R
#' @docType methods
#' @family signal processing methods
#' @aliases signal_bind-method
setGeneric(
  name = "signal_bind",
  def = function(...) standardGeneric("signal_bind")
)

#' Mean Intensities
#'
#' @param ... Any object that can be interpreted in a suitable way
#'  (see [grDevices::xy.coords()]).
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-mean.R
#' @docType methods
#' @family signal processing methods
#' @aliases signal_mean-method
setGeneric(
  name = "signal_mean",
  def = function(...) standardGeneric("signal_mean")
)

#' Subset
#'
#' @description
#'  * `signal_select()` allows to subset by values of `x`.
#'  * `signal_slice()` allows to subset by position along `x`.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param subset An [`integer`] vector giving either positive values to keep,
#'  or negative values to drop. The values provided must be either all
#'  positive or all negative (coerced to integer as by [as.integer()]).
#' @param from,to A [`numeric`] value giving the first and last value (in `x`
#'  unit) to be selected.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-subset.R
#' @docType methods
#' @family signal processing methods
#' @name subset
#' @rdname subset
NULL

#' @rdname subset
#' @aliases signal_select-method
setGeneric(
  name = "signal_select",
  def = function(x, y, ...) standardGeneric("signal_select")
)

#' @rdname subset
#' @aliases signal_slice-method
setGeneric(
  name = "signal_slice",
  def = function(x, y, ...) standardGeneric("signal_slice")
)

#' Shift the X Scale
#'
#' Shifts the `x` scale by a given value.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param lag A [`numeric`] vector specifying the offset.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-shift.R
#' @docType methods
#' @family signal processing methods
#' @aliases signal_shift-method
setGeneric(
  name = "signal_shift",
  def = function(x, y, lag, ...) standardGeneric("signal_shift")
)

#' Drift Intensities
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param lag A [`numeric`] vector specifying the offset or any object that can
#'  be interpreted in a suitable way (see [grDevices::xy.coords()])
#' @param subtract A [`logical`] scalar: should `lag` be subtracted to `y`?
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-drift.R
#' @docType methods
#' @family signal processing methods
#' @aliases signal_drift-method
setGeneric(
  name = "signal_drift",
  def = function(x, y, lag, ...) standardGeneric("signal_drift")
)

# Smooth =======================================================================
#' Rectangular Smoothing
#'
#' Unweighted sliding-average or rectangular Smoothing.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param m An odd [`integer`] giving the window size (i.e. the number of
#'  adjacent points to be used; see [`window_sliding()`]).
#' @param ... Currently not used.
#' @details
#'  It replaces each point in the signal with the average of \eqn{m} adjacent
#'  points.
#' @note
#'  There will be \eqn{(m - 1) / 2} points both at the beginning and at the end
#'  of the data series for which a complete \eqn{m}-width window cannot be
#'  obtained. To prevent data loss, progressively wider/narrower windows are
#'  used at both ends of the data series.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-smooth.R
#' @docType methods
#' @family smoothing methods
#' @aliases smooth_rectangular-method
setGeneric(
  name = "smooth_rectangular",
  def = function(x, y, ...) standardGeneric("smooth_rectangular")
)

#' Triangular Smoothing
#'
#' Weighted sliding-average or triangular smoothing.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param m An odd [`integer`] giving the window size (i.e. the number of
#'  adjacent points to be used; see [`window_sliding()`]).
#' @param ... Currently not used.
#' @details
#'  It replaces each point in the signal with the weighted mean of \eqn{m}
#'  adjacent points.
#' @note
#'  There will be \eqn{(m - 1) / 2} points both at the beginning and at the end
#'  of the data series for which a complete \eqn{m}-width window cannot be
#'  obtained. To prevent data loss, progressively wider/narrower windows are
#'  used at both ends of the data series.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-smooth.R
#' @docType methods
#' @family smoothing methods
#' @aliases smooth_triangular-method
setGeneric(
  name = "smooth_triangular",
  def = function(x, y, ...) standardGeneric("smooth_triangular")
)

#' Loess Smoothing
#'
#' Smoothes intensities by loess fitting.
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param span An [`integer`] specifying the degree of smoothing (see
#'  [stats::loess()]).
#' @param ... Extra arguments to be passed to [stats::loess()].
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @author N. Frerebeau
#' @example inst/examples/ex-smooth.R
#' @docType methods
#' @family smoothing methods
#' @aliases smooth_loess-method
setGeneric(
  name = "smooth_loess",
  def = function(x, y, ...) standardGeneric("smooth_loess")
)

#' Savitzky-Golay Filter
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param m An odd [`integer`] giving the window size (i.e. the number of
#'  adjacent points to be used).
#' @param p An [`integer`] giving the degree of the polynomial to be used.
#' @param ... Currently not used.
#' @details
#'  This method is based on the least-squares fitting of polynomials to
#'  segments of \eqn{m} adjacent points.
#' @note
#'  There will be \eqn{(m - 1) / 2} points both at the beginning and at the end
#'  of the data series for which a complete \eqn{m}-width window cannot be
#'  obtained. To prevent data loss, the original \eqn{(m - 1) / 2} points at
#'  both ends of the data series are preserved.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @references
#'  Gorry, P. A. (1990). General Least-Squares Smoothing and Differentiation by
#'  the Convolution (Savitzky-Golay) Method. *Analytical Chemistry*, 62(6),
#'  p. 570-573. \doi{10.1021/ac00205a007}.
#'
#'  Savitzky, A. & Golay, M. J. E. (1964). Smoothing and Differentiation of
#'  Data by Simplified Least Squares Procedures. *Analytical Chemistry*,
#'  36(8), p. 1627-1639. \doi{10.1021/ac60214a047}.
#' @author N. Frerebeau
#' @example inst/examples/ex-smooth.R
#' @docType methods
#' @family smoothing methods
#' @aliases smooth_savitzky-method
setGeneric(
  name = "smooth_savitzky",
  def = function(x, y, ...) standardGeneric("smooth_savitzky")
)

#' Whittaker Smoothing
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param lambda An [`integer`] giving the smoothing parameter. The larger
#'  `lambda` is, the smoother the curve.
#' @param d An [`integer`] specifying the order of the penalty.
#' @param sparse A [`logical`] scalar: should sparse matrices be used for
#'  computation? If `TRUE`, \pkg{Matrix} is required.
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @references
#'  Eilers, P. H. C. (2003). A Perfect Smoother. *Analytical Chemistry*,
#'  75(14): 3631-36. \doi{10.1021/ac034173t}.
#' @author N. Frerebeau
#' @example inst/examples/ex-smooth.R
#' @docType methods
#' @family smoothing methods
#' @aliases smooth_whittaker-method
setGeneric(
  name = "smooth_whittaker",
  def = function(x, y, ...) standardGeneric("smooth_whittaker")
)

#' Penalized Likelihood Smoothing
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param lambda An [`integer`] giving the smoothing parameter. The larger
#'  `lambda` is, the smoother the curve.
#' @param d An [`integer`] specifying the order of the penalty.
#' @param SE A [`logical`] scalar: should standard errors be returned?
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @note
#'  \pkg{Matrix} is required.
#' @references
#'  de Rooi, J. J., van der Pers, N. M., Hendrikx, R. W. A., Delhez, R.,
#'  Böttger A. J. and Eilers, P. H. C. (2014). Smoothing of X-ray diffraction
#'  data and Ka2 elimination using penalized likelihood and the composite link
#'  model. *Journal of Applied Crystallography*, 47: 852-860.
#'  \doi{10.1107/S1600576714005809}
#' @author J. J. de Rooi *et al.* (original R code).
#' @example inst/examples/ex-xrd-ka2.R
#' @docType methods
#' @family smoothing methods
#' @aliases smooth_likelihood-method
setGeneric(
  name = "smooth_likelihood",
  def = function(x, y, ...) standardGeneric("smooth_likelihood")
)

# Strip ========================================================================
#' Strip XRD ka2
#'
#' @param x,y A [`numeric`] vector. If `y` is missing, an attempt is made to
#'  interpret `x` in a suitable way (see [grDevices::xy.coords()]).
#' @param lambda An [`integer`] giving the smoothing parameter. The larger
#'  `lambda` is, the smoother the curve.
#' @param wave A length-two [`numeric`] vector giving the characteristic
#'  wavelengths of the anode material (defaults to copper).
#' @param tau A length-one [`numeric`] vector giving the ratio between
#'  \eqn{\alpha}1 and \eqn{\alpha}2 line intensities (defaults to 1/2).
#' @param nseg A length-one [`numeric`] vector specifying the number of equally
#'  sized segments for B-spline basis matrix computation.
#' @param progress A [`logical`] scalar: should a progress bar be displayed?
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] with two components `x` and `y`.
#' @note
#'  \pkg{Matrix} is required.
#' @author J. J. de Rooi *et al.* (original R code).
#' @references
#'  de Rooi, J. J., van der Pers, N. M., Hendrikx, R. W. A., Delhez, R.,
#'  Böttger A. J. and Eilers, P. H. C. (2014). Smoothing of X-ray diffraction
#'  data and Ka2 elimination using penalized likelihood and the composite link
#'  model. *Journal of Applied Crystallography*, 47: 852-860.
#'  \doi{10.1107/S1600576714005809}
#' @example inst/examples/ex-xrd-ka2.R
#' @docType methods
#' @family specialized tools
#' @aliases ka2_strip_penalized-method
setGeneric(
  name = "ka2_strip_penalized",
  def = function(x, y, ...) standardGeneric("ka2_strip_penalized")
)

# Windows ======================================================================
#' Sliding Windows
#'
#' @param n An [`integer`] giving the length of the data series (will be coerced
#'  with [`as.integer()`] and hence truncated toward zero).
#' @param m An odd [`integer`] giving the window size, i.e. the number of
#'  adjacent points to be used (will be coerced with [`as.integer()`] and hence
#'  truncated toward zero).
#' @param i A vector [`integer`] specifying the indices of the data points for
#'  which windows should be returned. If `NULL` (the default), windows are
#'  evaluated for each data point.
#' @description
#'  There will be \eqn{(m - 1) / 2} points both at the beginning and at the end
#'  of the data series for which a complete \eqn{m}-width window cannot be
#'  obtained. To prevent data loss, progressively wider/narrower windows are
#'  evaluated at both ends of the data series.
#' @param ... Currently not used.
#' @return
#'  Returns a length-\eqn{n} [`list`] of [`integer`] vectors (indices of the
#'  data points in each window).
#' @author N. Frerebeau
#' @example inst/examples/ex-windows.R
#' @docType methods
#' @family moving windows
#' @aliases window_sliding-method
setGeneric(
  name = "window_sliding",
  def = function(n, m, ...) standardGeneric("window_sliding")
)

#' Tumbling Windows
#'
#' @param n An [`integer`] giving the length of the data series (will be coerced
#'  with [`as.integer()`] and hence truncated toward zero).
#' @param m An [`integer`] giving the window size, i.e. the number of
#'  adjacent points to be used (will be coerced with [`as.integer()`] and hence
#'  truncated toward zero).
#' @param drop A [`logical`] scalar: if `m` is not a multiple of `n`, should the
#' last data points be removed so that all windows have the same length?
#' @param ... Currently not used.
#' @return
#'  Returns a [`list`] of [`integer`] vectors (indices of the data points in
#'  each window).
#' @author N. Frerebeau
#' @example inst/examples/ex-windows.R
#' @docType methods
#' @family moving windows
#' @aliases window_tumbling-method
setGeneric(
  name = "window_tumbling",
  def = function(n, m, ...) standardGeneric("window_tumbling")
)

# RESAMPLE
#' @include AllGenerics.R
NULL

# Bin ==========================================================================
#' @export
#' @rdname resample_bin
#' @aliases resample_bin,numeric,numeric-method
setMethod(
  f = "resample_bin",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, by, f = mean, ...) {
    assert_length(y, length(x))

    k <- length(x) %% by
    if (k != 0) {
      x <- utils::tail(x, -k)
      y <- utils::tail(y, -k)
    }
    i <- cut(seq_along(x), breaks = length(x) / by, labels = FALSE)

    mid <- tapply(X = x, INDEX = i, FUN = mean, simplify = FALSE)
    bin <- tapply(X = y, INDEX = i, FUN = f, ..., simplify = FALSE)

    xy <- list(x = unlist(mid), y = unlist(bin))
    xy
  }
)

#' @export
#' @rdname resample_bin
#' @aliases resample_bin,ANY,missing-method
setMethod(
  f = "resample_bin",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, y, by, f = sum) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, by = by, f = f)
  }
)

# Downsample ===================================================================
#' @export
#' @rdname resample_down
#' @aliases resample_down,numeric,numeric-method
setMethod(
  f = "resample_down",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, by) {
    assert_length(y, length(x))

    i <- seq(from = 1, to = length(x), by = by)
    xy <- list(x = x[i], y = y[i])
    xy
  }
)

#' @export
#' @rdname resample_down
#' @aliases resample_down,ANY,missing-method
setMethod(
  f = "resample_down",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, y, by) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, by = by)
  }
)

# Interpolate ==================================================================
#' @export
#' @rdname resample_interpolate
#' @aliases resample_interpolate,numeric,numeric-method
setMethod(
  f = "resample_interpolate",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, from, to, by, ...) {
    assert_length(y, length(x))

    ## New x scale
    x_scale <- seq(from = from, to = to, by = by)

    ## Interpolate
    new_data <- stats::approx(x = x, y = y, xout = x_scale, ...)

    xy <- list(x = new_data$x, y = new_data$y)
    xy
  }
)

#' @export
#' @rdname resample_interpolate
#' @aliases resample_interpolate,ANY,missing-method
setMethod(
  f = "resample_interpolate",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, y, from, to, by, ...) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, from = from, to = to, by = by, ...)
  }
)

# SIGNAL INTEGRATION
#' @include AllGenerics.R
NULL

#' @export
#' @rdname integrate_rectangle
#' @aliases integrate_rectangle,numeric,numeric-method
setMethod(
  f = "integrate_rectangle",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y, right = FALSE) {
    assert_length(y, length(x))

    h <- if (right) utils::tail(y, -1) else utils::head(y, -1)
    sum(h * diff(x))
  }
)

#' @export
#' @rdname integrate_rectangle
#' @aliases integrate_rectangle,ANY,missing-method
setMethod(
  f = "integrate_rectangle",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x, right = FALSE) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y, right = right)
  }
)

#' @export
#' @rdname integrate_trapezoid
#' @aliases integrate_trapezoid,numeric,numeric-method
setMethod(
  f = "integrate_trapezoid",
  signature = signature(x = "numeric", y = "numeric"),
  definition = function(x, y) {
    assert_length(y, length(x))

    sum((utils::head(y, -1) + utils::tail(y, -1)) * diff(x) / 2)
  }
)

#' @export
#' @rdname integrate_trapezoid
#' @aliases integrate_trapezoid,ANY,missing-method
setMethod(
  f = "integrate_trapezoid",
  signature = signature(x = "ANY", y = "missing"),
  definition = function(x) {
    xy <- grDevices::xy.coords(x)
    methods::callGeneric(x = xy$x, y = xy$y)
  }
)

1		# SMOOTHING
2		#' @include AllGenerics.R
3		NULL
4
5		# Unweighted smoothing =========================================================
6		#' @export
7		#' @rdname smooth_rectangular
8		#' @aliases smooth_rectangular,numeric,numeric-method
9		setMethod(
10		f = "smooth_rectangular",
11		signature = signature(x = "numeric", y = "numeric"),
12		definition = function(x, y, m = 3) {
13	!	assert_length(y, length(x))
14
15		## Windows
16	!	index <- window_sliding(length(x), m)
17
18	!	y <- vapply(
19	!	X = index,
20	!	FUN = function(i, data) mean(data[i]),
21	!	FUN.VALUE = numeric(1),
22	!	data = y
23		)
24
25	!	xy <- list(x = x, y = y)
26	!	attr(xy, "method") <- "unweighted smoothing"
27	!	xy
28		}
29		)
30
31		#' @export
32		#' @rdname smooth_rectangular
33		#' @aliases smooth_rectangular,ANY,missing-method
34		setMethod(
35		f = "smooth_rectangular",
36		signature = signature(x = "ANY", y = "missing"),
37		definition = function(x, m) {
38	!	xy <- grDevices::xy.coords(x)
39	!	methods::callGeneric(x = xy$x, y = xy$y, m = m)
40		}
41		)
42
43		# Weighted smoothing ===========================================================
44		#' @export
45		#' @rdname smooth_triangular
46		#' @aliases smooth_triangular,numeric,numeric-method
47		setMethod(
48		f = "smooth_triangular",
49		signature = signature(x = "numeric", y = "numeric"),
50		definition = function(x, y, m = 3) {
51		# Validation
52	!	assert_length(y, length(x))
53	!	assert_odd(m)
54
55		# Index
56	!	k <- (m - 1) / 2
57	!	index_k <- seq_len(k)
58	!	index_y <- seq_along(y)
59	!	index_m <- c(index_k, rep_len(k + 1, length(y) - 2 * k), rev(index_k)) - 1
60
61	!	y <- mapply(
62	!	FUN = function(i, k, data) {
63	!	j <- seq_len(k)
64	!	w <- c(j, k + 1, rev(j))
65	!	index <- seq(from = i - k, to = i + k, by = 1)
66	!	stats::weighted.mean(x = data[index], w = w)
67		},
68	!	i = index_y, k = index_m,
69	!	MoreArgs = list(data = y)
70		)
71
72	!	xy <- list(x = x, y = y)
73	!	attr(xy, "method") <- "weighted smoothing"
74	!	xy
75		}
76		)
77
78		#' @export
79		#' @rdname smooth_triangular
80		#' @aliases smooth_triangular,ANY,missing-method
81		setMethod(
82		f = "smooth_triangular",
83		signature = signature(x = "ANY", y = "missing"),
84		definition = function(x, m) {
85	!	xy <- grDevices::xy.coords(x)
86	!	methods::callGeneric(x = xy$x, y = xy$y, m = m)
87		}
88		)
89
90		# Loess smoothing ==============================================================
91		#' @export
92		#' @rdname smooth_loess
93		#' @aliases smooth_loess,numeric,numeric-method
94		setMethod(
95		f = "smooth_loess",
96		signature = signature(x = "numeric", y = "numeric"),
97		definition = function(x, y, span = 0.75, ...) {
98	!	assert_length(y, length(x))
99
100	!	points <- data.frame(x, y)
101	!	fit <- stats::loess(y ~ x, data = points, span = span, ...)
102	!	bsl <- stats::predict(fit)
103
104	!	xy <- list(x = x, y = bsl)
105	!	attr(xy, "method") <- "loess smoothing"
106	!	xy
107		}
108		)
109
110		#' @export
111		#' @rdname smooth_loess
112		#' @aliases smooth_loess,ANY,missing-method
113		setMethod(
114		f = "smooth_loess",
115		signature = signature(x = "ANY", y = "missing"),
116		definition = function(x, span = 0.75, ...) {
117	!	xy <- grDevices::xy.coords(x)
118	!	methods::callGeneric(x = xy$x, y = xy$y, span = span, ...)
119		}
120		)
121
122		# Savitzky-Golay filter ========================================================
123		#' @export
124		#' @rdname smooth_savitzky
125		#' @aliases smooth_savitzky,numeric,numeric-method
126		setMethod(
127		f = "smooth_savitzky",
128		signature = signature(x = "numeric", y = "numeric"),
129		definition = function(x, y, m = 3, p = 2) {
130		## Validation
131	!	assert_length(y, length(x))
132	!	assert_odd(m)
133
134	!	k <- (m - 1) / 2
135	!	i <- seq(from = -k, to = k, by = 1)
136	!	j <- utils::head(utils::tail(seq_along(y), n = -k), n = -k)
137	!	conv <- coef_savitzky(m, p)
138
139	!	smoothed <- vapply(
140	!	X = j,
141	!	FUN = function(j, i, conv, data) {
142	!	sum(conv * data[j + i])
143		},
144	!	FUN.VALUE = double(1),
145	!	i = i,
146	!	conv = conv,
147	!	data = y
148		)
149	!	y[j] <- smoothed
150
151	!	xy <- list(x = x, y = y)
152	!	attr(xy, "method") <- "Savitzky-Golay filter"
153	!	xy
154		}
155		)
156
157		#' @export
158		#' @rdname smooth_savitzky
159		#' @aliases smooth_savitzky,ANY,missing-method
160		setMethod(
161		f = "smooth_savitzky",
162		signature = signature(x = "ANY", y = "missing"),
163		definition = function(x, m, p) {
164	!	xy <- grDevices::xy.coords(x)
165	!	methods::callGeneric(x = xy$x, y = xy$y, m = m, p = p)
166		}
167		)
168
169		coef_savitzky <- function(m, p = 2) {
170	!	k <- (m - 1) / 2
171	!	z <- seq(from = -k, to = k, by = 1)
172	!	J <- vapply(X = c(0, p), FUN = function(p, z) z^p, z, FUN.VALUE = double(m))
173	!	(solve(t(J) %% J) %% t(J))[1, , drop = TRUE]
174		}
175
176		# Whittaker smoothing ==========================================================
177		#' @export
178		#' @rdname smooth_whittaker
179		#' @aliases smooth_whittaker,numeric,numeric-method
180		setMethod(
181		f = "smooth_whittaker",
182		signature = signature(x = "numeric", y = "numeric"),
183		definition = function(x, y, lambda = 1600, d = 2, sparse = FALSE) {
184	!	assert_length(y, length(x))
185	!	m <- length(y)
186
187	!	if (sparse) {
188	!	assert_Matrix()
189
190	!	E <- Matrix::Diagonal(m)
191	!	D <- Matrix::diff(E, lag = 1, differences = d)
192	!	B <- Matrix::chol(E + (lambda * Matrix::t(D) %*% D))
193	!	z <- Matrix::solve(B, Matrix::solve(Matrix::t(B), y))
194		## Prior to v1.6, Matrix::solve(a=<Matrix>, b=<vector>) returns a matrix
195	!	z <- as.numeric(z)
196		} else {
197	!	E <- diag(m)
198	!	D <- diff(E, lag = 1, differences = d)
199	!	B <- E + (lambda * t(D) %*% D)
200	!	z <- solve(B, y)
201		}
202
203	!	xy <- list(x = x, y = z)
204	!	attr(xy, "method") <- "Whittaker smoothing"
205	!	xy
206		}
207		)
208
209		#' @export
210		#' @rdname smooth_whittaker
211		#' @aliases smooth_whittaker,ANY,missing-method
212		setMethod(
213		f = "smooth_whittaker",
214		signature = signature(x = "ANY", y = "missing"),
215		definition = function(x, lambda = 1600, d = 2, sparse = FALSE) {
216	!	xy <- grDevices::xy.coords(x)
217	!	methods::callGeneric(x = xy$x, y = xy$y, lambda = lambda, d = d,
218	!	sparse = sparse)
219		}
220		)
221
222		# Penalized likelihood smoothing ===============================================
223		#' @export
224		#' @rdname smooth_likelihood
225		#' @aliases smooth_likelihood,numeric,numeric-method
226		setMethod(
227		f = "smooth_likelihood",
228		signature = signature(x = "numeric", y = "numeric"),
229		definition = function(x, y, lambda, d = 2, SE = FALSE,
230		progress = interactive()) {
231	!	assert_Matrix()
232	!	assert_length(y, length(x))
233
234	!	m <- length(y)
235	!	n <- length(lambda)
236
237	!	aic <- rep(0, n)
238	!	mu <- matrix(data = 0, nrow = m, ncol = n)
239	!	se <- matrix(data = 0, nrow = m, ncol = n)
240	!	res <- matrix(data = 0, nrow = m, ncol = n)
241
242	!	if (progress) pb <- utils::txtProgressBar(min = 0, max = n, style=3)
243
244	!	n_lambda <- seq_len(n)
245	!	for (i in n_lambda) {
246	!	like <- penalized_likelihood(y, lambda = lambda[i], d = d, SE = SE)
247
248	!	aic[i] <- like$aic
249	!	mu[, i] <- like$mu
250	!	se[, i] <- like$se
251	!	res[, i] <- like$res
252
253	!	if (progress) utils::setTxtProgressBar(pb, i)
254		}
255
256	!	if (progress) close(pb)
257
258	!	k <- which.min(aic)
259	!	mu <- mu[, k]
260	!	se <- se[, k]
261	!	res <- res[, k]
262	!	op_lambda <- lambda[k]
263	!	op_aic <- aic[k]
264
265	!	xy <- list(x = x, y = mu, aic = aic, se = se, residuals = res,
266	!	optimal_aic = op_aic, optimal_lambda = op_lambda)
267	!	attr(xy, "method") <- "penalized likelihood smoothing"
268	!	xy
269		}
270		)
271
272		#' @export
273		#' @rdname smooth_likelihood
274		#' @aliases smooth_likelihood,ANY,missing-method
275		setMethod(
276		f = "smooth_likelihood",
277		signature = signature(x = "ANY", y = "missing"),
278		definition = function(x, lambda, d = 2, SE = FALSE,
279		progress = interactive()) {
280	!	xy <- grDevices::xy.coords(x)
281	!	methods::callGeneric(x = xy$x, y = xy$y, lambda = lambda, d = d,
282	!	SE = SE, progress = progress)
283		}
284		)
285
286		penalized_likelihood <- function(y, lambda, d = 2, SE = FALSE) {
287
288	!	m <- length(y)
289	!	z <- log(y + 1)
290	!	E <- Matrix::Diagonal(m)
291	!	D <- Matrix::diff(E, differences = d)
292	!	P <- lambda * Matrix::crossprod(D, D)
293	!	mu <- exp(z)
294
295	!	delta_mu <- 1
296	!	limit <- 1 # Stop condition to prevent infinite loop
297	!	while (delta_mu >= 0.0001 && limit <= 10) {
298	!	W <- Matrix::Diagonal(x = mu)
299	!	z <- Matrix::solve(W + P, y - mu + mu * z)
300	!	new_mu <- exp(as.numeric(z)) # Drop S4 class
301	!	delta_mu <- max(abs(new_mu - mu))
302	!	mu <- new_mu
303
304	!	limit <- limit + 1
305		}
306
307		## AIC
308	!	H <- Matrix::solve(W + P) %*% W
309	!	ed <- sum(Matrix::diag(H))
310	!	dev <- 2 * sum(y * log((y + 1e-10) / mu))
311	!	aic <- dev + 2 * ed * m / (m - ed)
312
313		## Standard error
314	!	se <- rep(NA_real_, m)
315	!	if (isTRUE(SE)) {
316	!	HH <- Matrix::tcrossprod(H %*% W, H)
317	!	se <- sqrt(Matrix::diag(HH))
318		}
319
320		## Standardized residuals for counts
321	!	res <- (y - mu) / sqrt(mu)
322
323	!	list(mu = mu, aic = aic, se = se, res = res)
324		}

1		# SLIDING WINDOWS
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname window_sliding
7		#' @aliases window_sliding,integer,integer-method
8		setMethod(
9		f = "window_sliding",
10		signature = signature(n = "integer", m = "integer"),
11		definition = function(n, m, i = NULL) {
12		## Validation
13	3x	assert_odd(m)
14
15		## Index
16	3x	k <- (m - 1) / 2
17	3x	x <- if (is.null(i)) seq_len(n) else i
18
19	3x	left <- x - k
20	3x	left <- ifelse(left > 0, left, 1)
21
22	3x	right <- x + k
23	3x	right <- ifelse(right < n, right, n)
24
25	3x	mapply(
26	3x	FUN = seq,
27	3x	from = left,
28	3x	to = right,
29	3x	MoreArgs = list(by = 1),
30	3x	SIMPLIFY = FALSE
31		)
32		}
33		)
34
35		#' @export
36		#' @rdname window_sliding
37		#' @aliases window_sliding,numeric,numeric-method
38		setMethod(
39		f = "window_sliding",
40		signature = signature(n = "numeric", m = "numeric"),
41		definition = function(n, m, i = NULL) {
42	3x	n <- as.integer(n)
43	3x	m <- as.integer(m)
44	3x	methods::callGeneric(n, m, i = i)
45		}
46		)
47
48		#' @export
49		#' @rdname window_tumbling
50		#' @aliases window_tumbling,integer,integer-method
51		setMethod(
52		f = "window_tumbling",
53		signature = signature(n = "integer", m = "integer"),
54		definition = function(n, m, drop = FALSE) {
55	2x	k <- n %% m
56	2x	if (drop && k != 0) {
57	1x	n <- n - k
58		}
59
60	2x	x <- seq_len(n)
61	2x	i <- cut(x, breaks = n / m, labels = FALSE)
62	2x	split(x, f = i)
63		}
64		)
65
66		#' @export
67		#' @rdname window_tumbling
68		#' @aliases window_tumbling,numeric,numeric-method
69		setMethod(
70		f = "window_tumbling",
71		signature = signature(n = "numeric", m = "numeric"),
72		definition = function(n, m, drop = FALSE) {
73	2x	n <- as.integer(n)
74	2x	m <- as.integer(m)
75	2x	methods::callGeneric(n, m, drop = drop)
76		}
77		)

1		# BASELINE
2		#' @include AllGenerics.R
3		NULL
4
5		# Linear =======================================================================
6		#' @export
7		#' @rdname baseline_linear
8		#' @aliases baseline_linear,numeric,numeric-method
9		setMethod(
10		f = "baseline_linear",
11		signature = c(x = "numeric", y = "numeric"),
12		definition = function(x, y, points = range(x)) {
13	!	assert_length(y, length(x))
14
15		## Find the nearest value
16	!	z <- vapply(X = points, FUN = function(i, x) which_nearest(x, i),
17	!	FUN.VALUE = numeric(1), x = x)
18	!	zi <- which(x >= min(points) & x <= max(points))
19
20	!	data <- data.frame(x, y)
21	!	fit <- stats::lm(y ~ x, data = data, subset = z)
22	!	pred <- stats::predict(fit, data.frame(x = x[zi]))
23
24	!	bsl <- rep(NA_real_, length(x))
25	!	bsl[zi] <- pred
26
27	!	xy <- list(x = x, y = bsl)
28	!	attr(xy, "method") <- "linear baseline"
29	!	xy
30		}
31		)
32
33		#' @export
34		#' @rdname baseline_linear
35		#' @aliases baseline_linear,ANY,missing-method
36		setMethod(
37		f = "baseline_linear",
38		signature = c(x = "ANY", y = "missing"),
39		definition = function(x, points = range(x)) {
40	!	xy <- grDevices::xy.coords(x)
41	!	methods::callGeneric(x = xy$x, y = xy$y, points = points)
42		}
43		)
44
45		# Polynomial ===================================================================
46		#' @export
47		#' @rdname baseline_polynomial
48		#' @aliases baseline_polynomial,numeric,numeric-method
49		setMethod(
50		f = "baseline_polynomial",
51		signature = c(x = "numeric", y = "numeric"),
52		definition = function(x, y, d = 3, tolerance = 0.001, stop = 100) {
53	!	assert_length(y, length(x))
54
55	!	polynom <- cbind(1 / sqrt(length(x)), stats::poly(x, degree = d))
56
57	!	old_y <- y
58	!	start <- 0
59	!	convergence <- FALSE
60	!	while (!convergence) {
61	!	new_y <- polynom %*% crossprod(polynom, old_y)
62	!	new_y <- pmin(y, new_y)
63
64	!	criterion <- sum(abs((new_y - old_y) / old_y), na.rm = TRUE)
65	!	convergence <- criterion < tolerance
66
67	!	old_y <- new_y
68	!	start <- start + 1
69	!	if (start >= stop) {
70	!	warning(tr_("Convergence not reached."), call. = FALSE)
71	!	break
72		}
73		}
74
75	!	xy <- list(x = x, y = new_y)
76	!	attr(xy, "method") <- "polynomial baseline"
77	!	xy
78		}
79		)
80
81		#' @export
82		#' @rdname baseline_polynomial
83		#' @aliases baseline_polynomial,ANY,missing-method
84		setMethod(
85		f = "baseline_polynomial",
86		signature = c(x = "ANY", y = "missing"),
87		definition = function(x, d = 3, tolerance = 0.001, stop = 100) {
88	!	xy <- grDevices::xy.coords(x)
89	!	methods::callGeneric(x = xy$x, y = xy$y, d = d, tolerance = tolerance,
90	!	stop = stop)
91		}
92		)
93
94		# Rolling Ball =================================================================
95		#' @export
96		#' @rdname baseline_rollingball
97		#' @aliases baseline_rollingball,numeric,numeric-method
98		setMethod(
99		f = "baseline_rollingball",
100		signature = c(x = "numeric", y = "numeric"),
101		definition = function(x, y, m, s) {
102	!	n <- length(x)
103	!	assert_length(y, n)
104
105		## Windows
106	!	win_minmax <- window_sliding(n, m)
107	!	win_smooth <- window_sliding(n, s)
108
109		## Minimize
110	!	T1 <- vapply(
111	!	X = win_minmax,
112	!	FUN = function(i, data) {
113	!	min(data[i])
114		},
115	!	FUN.VALUE = numeric(1),
116	!	data = y
117		)
118
119		## Maximize
120	!	T2 <- vapply(
121	!	X = win_minmax,
122	!	FUN = function(i, data) {
123	!	max(data[i])
124		},
125	!	FUN.VALUE = numeric(1),
126	!	data = T1
127		)
128
129		## Smooth
130	!	T3 <- vapply(
131	!	X = win_smooth,
132	!	FUN = function(i, data) {
133	!	mean(data[i])
134		},
135	!	FUN.VALUE = numeric(1),
136	!	data = T2
137		)
138
139	!	xy <- list(x = x, y = T3)
140	!	attr(xy, "method") <- "rolling ball baseline"
141	!	xy
142		}
143		)
144
145		#' @export
146		#' @rdname baseline_rollingball
147		#' @aliases baseline_rollingball,ANY,missing-method
148		setMethod(
149		f = "baseline_rollingball",
150		signature = c(x = "ANY", y = "missing"),
151		definition = function(x, m, s) {
152	!	xy <- grDevices::xy.coords(x)
153	!	methods::callGeneric(x = xy$x, y = xy$y, m = m, s = s)
154		}
155		)
156
157		# Rubberband ===================================================================
158		#' @export
159		#' @rdname baseline_rubberband
160		#' @aliases baseline_rubberband,numeric,numeric-method
161		setMethod(
162		f = "baseline_rubberband",
163		signature = c(x = "numeric", y = "numeric"),
164		definition = function(x, y, noise = 0, spline = TRUE, ...) {
165	!	assert_length(y, length(x))
166
167		## (chull returns points in clockwise order)
168	!	pts <- grDevices::chull(x, y)
169
170		## Check that ncol(y) is a position 1,
171		## if not, rotate chull vertex index so that ncol(y) is at position 1
172		## then keep only index from ncol(y) to 1 (i.e. lower part of the hull)
173	!	v_max <- which.max(pts) - 1
174	!	if (v_max > 0) pts <- c(pts[-seq_len(v_max)], pts[seq_len(v_max)])
175	!	pts <- pts[1:which.min(pts)]
176		# First and last point must be minima, if not remove them
177	!	if (pts[2] == pts[1] - 1) pts <- pts[-1] # Last point
178	!	pts <- rev(pts) # Sort in ascending order
179	!	if (pts[2] == pts[1] + 1) pts <- pts[-1] # First point
180
181	!	bsl <- stats::approx(x = x[pts], y = y[pts], xout = x, method = "linear")$y
182	!	if (spline) {
183	!	pts <- which(y <= bsl + noise)
184	!	if (length(pts) > 3) {
185	!	spl <- stats::smooth.spline(x[pts], y[pts], ...)
186	!	tmp <- stats::predict(spl, x, 0)$y
187		} else {
188	!	tmp <- stats::spline(x[pts], y[pts], xout = x)$y
189		}
190		}
191
192		## Check baseline
193	!	if (anyNA(tmp)) {
194	!	msg <- tr_("Failed to estimate the baseline, please check your parameters.")
195	!	stop(msg, call. = FALSE)
196		}
197
198	!	xy <- list(x = x, y = bsl)
199	!	attr(xy, "method") <- "rubberband baseline"
200	!	xy
201		}
202		)
203
204		#' @export
205		#' @rdname baseline_rubberband
206		#' @aliases baseline_rubberband,ANY,missing-method
207		setMethod(
208		f = "baseline_rubberband",
209		signature = c(x = "ANY", y = "missing"),
210		definition = function(x, noise = 0, spline = TRUE, ...) {
211	!	xy <- grDevices::xy.coords(x)
212	!	methods::callGeneric(x = xy$x, y = xy$y, noise = noise, spline = spline, ...)
213		}
214		)
215
216		# SNIP =========================================================================
217		#' @export
218		#' @rdname baseline_snip
219		#' @aliases baseline_snip,numeric,numeric-method
220		setMethod(
221		f = "baseline_snip",
222		signature = c(x = "numeric", y = "numeric"),
223		definition = function(x, y, LLS = FALSE, decreasing = FALSE, n = 100) {
224	!	assert_length(y, length(x))
225
226		## LLS operator
227	!	y <- if (LLS) LLS(y) else y
228
229	!	k <- length(y)
230	!	iter <- if (decreasing) rev(seq_len(n)) else seq_len(n)
231
232	!	tmp <- y
233	!	for (p in iter) {
234	!	for (i in p:(k - p)) {
235	!	a <- y[i]
236	!	b <- (y[i - p] + y[i + p]) / 2
237	!	tmp[i] <- min(a, b)
238		}
239	!	y <- tmp
240		}
241
242		## Inverse LLS operator
243	!	bsl <- if (LLS) inverseLLS(y) else y
244
245		## Check baseline
246	!	if (anyNA(bsl)) {
247	!	msg <- tr_("Failed to estimate the baseline, please check your parameters.")
248	!	stop(msg, call. = FALSE)
249		}
250
251	!	xy <- list(x = x, y = bsl)
252	!	attr(xy, "method") <- "SNIP"
253	!	xy
254		}
255		)
256
257		#' @export
258		#' @rdname baseline_snip
259		#' @aliases baseline_snip,ANY,missing-method
260		setMethod(
261		f = "baseline_snip",
262		signature = c(x = "ANY", y = "missing"),
263		definition = function(x, LLS = FALSE, decreasing = FALSE, n = 100) {
264	!	xy <- grDevices::xy.coords(x)
265	!	methods::callGeneric(x = xy$x, y = xy$y, LLS = LLS, decreasing = decreasing,
266	!	n = n)
267		}
268		)
269
270		#' LLS operator
271		#'
272		#' @param x A [`numeric`] vector.
273		#' @keywords internal
274		#' @noRd
275		LLS <- function(x) {
276	!	log(log(sqrt(x + 1) + 1) + 1)
277		}
278		inverseLLS <- function(x) {
279	!	(exp(exp(x) - 1) - 1)^2 - 1
280		}
281
282		# AsLS =========================================================================
283		#' @export
284		#' @rdname baseline_asls
285		#' @aliases baseline_asls,numeric,numeric-method
286		setMethod(
287		f = "baseline_asls",
288		signature = c(x = "numeric", y = "numeric"),
289		definition = function(x, y, p = 0.01, lambda = 10^4, stop = 100) {
290	!	assert_Matrix()
291	!	assert_length(y, length(x))
292
293	!	m <- length(y)
294	!	E <- Matrix::Diagonal(m)
295	!	D <- Matrix::diff(E, lag = 1, differences = 2)
296
297	!	w <- rep(1, m) # weights
298
299	!	start <- 0
300	!	convergence <- FALSE
301	!	while (!convergence) {
302	!	W <- Matrix::Diagonal(x = w)
303	!	C <- Matrix::chol(W + (lambda * Matrix::t(D) %*% D))
304	!	z <- Matrix::solve(C, Matrix::solve(Matrix::t(C), w * y))
305		## Prior to v1.6, Matrix::solve(a=<Matrix>, b=<vector>) returns a matrix
306	!	z <- as.numeric(z)
307	!	w0 <- p * (y > z) + (1 - p) * (y < z)
308	!	if (isTRUE(all.equal(w, w0))) convergence <- TRUE
309
310	!	w <- w0
311	!	start <- start + 1
312	!	if (start >= stop) {
313	!	warning(tr_("Convergence not reached."), call. = FALSE)
314	!	break
315		}
316		}
317
318	!	xy <- list(x = x, y = z)
319	!	attr(xy, "method") <- "AsLS"
320	!	xy
321		}
322		)
323
324		#' @export
325		#' @rdname baseline_asls
326		#' @aliases baseline_asls,ANY,missing-method
327		setMethod(
328		f = "baseline_asls",
329		signature = c(x = "ANY", y = "missing"),
330		definition = function(x, p = 0.01, lambda = 10^4, stop = 100) {
331	!	xy <- grDevices::xy.coords(x)
332	!	methods::callGeneric(x = xy$x, y = xy$y, lambda = lambda, p = p, stop = stop)
333		}
334		)
335
336		# Peak Filling =================================================================
337		#' @export
338		#' @rdname baseline_peakfilling
339		#' @aliases baseline_peakfilling,numeric,numeric-method
340		setMethod(
341		f = "baseline_peakfilling",
342		signature = c(x = "numeric", y = "numeric"),
343		definition = function(x, y, n, m, by = 10, lambda = 1600, d = 2, sparse = FALSE) {
344	!	assert_length(y, length(x))
345
346		## Number of bucket intervals
347	!	by <- length(x) / by
348
349		## Exponential decrease in interval width
350	!	win <- m
351	!	if (n != 1) {
352	!	i <- log10(m) * (1 - (0:(n - 2)) / (n - 1))
353	!	win <- ceiling((10)^c(i, 0))
354		}
355
356		## Smoothing
357	!	smoothed <- smooth_whittaker(x, y, lambda = lambda, d = d, sparse = sparse)
358	!	x <- smoothed$x
359	!	y <- smoothed$y
360
361		## Subsampling
362	!	breaks <- cut(seq_along(x), breaks = by, labels = FALSE)
363	!	mid <- unlist(tapply(X = x, INDEX = breaks, FUN = mean, simplify = FALSE))
364	!	bin <- unlist(tapply(X = y, INDEX = breaks, FUN = min, simplify = FALSE))
365
366		## Suppression
367	!	iter <- seq_len(n)
368	!	for (i in iter) {
369	!	win_0 <- win[i] # Current half window width
370
371		## Interval cut-off close to edges
372	!	k <- seq(from = 2, to = by - 1, by = 1)
373	!	j <- by - k + 1
374	!	v <- pmin(k - 1, win_0, by - k)
375
376		## Point-wise iteration to the right
377	!	a <- mapply(function(k, v, bin) { mean(bin[(k - v):(k + v)]) }, k, v,
378	!	MoreArgs = list(bin = bin))
379	!	bin[k] <- pmin(a, bin[k]) # Baseline suppression
380
381		## Point-wise iteration to the left
382	!	a <- mapply(function(j, v, bin) { mean(bin[(j - v):(j + v)]) }, j, v,
383	!	MoreArgs = list(bin = bin))
384	!	bin[j] <- pmin(a, bin[j]) # Baseline suppression
385		}
386
387		## Stretch
388	!	mid[1] <- x[1]
389	!	mid[by] <- x[length(x)]
390	!	bsl <- stats::approx(mid, bin, x)$y
391
392	!	xy <- list(x = x, y = bsl)
393	!	attr(xy, "method") <- "4S Peak Filling"
394	!	xy
395		}
396		)
397
398		#' @export
399		#' @rdname baseline_peakfilling
400		#' @aliases baseline_peakfilling,ANY,missing-method
401		setMethod(
402		f = "baseline_peakfilling",
403		signature = c(x = "ANY", y = "missing"),
404		definition = function(x, n, m, by = 10, lambda = 1600, d = 2, sparse = FALSE) {
405	!	xy <- grDevices::xy.coords(x)
406	!	methods::callGeneric(x = xy$x, y = xy$y, n = n, m = m, by = by,
407	!	lambda = lambda, d = d, sparse = sparse)
408		}
409		)

1		# HELPERS
2
3		## https://michaelchirico.github.io/potools/articles/developers.html
4		tr_ <- function(...) {
5	3x	enc2utf8(gettext(paste0(...), domain = "R-alkahest"))
6		}
7
8		#' Find the Nearest Value in a Vector
9		#'
10		#' @param x A [`numeric`] vector.
11		#' @param value A [`numeric`] value.
12		#' @return An [`integer`].
13		#' @keywords internal
14		#' @noRd
15		which_nearest <- function(x, value) {
16	3x	which.min(abs(x - value))
17		}
18
19		#' Check Odd Numbers
20		#'
21		#' @param x A [`numeric`] vector.
22		#' @return
23		#' Throws an error, if any, and returns `x` invisibly otherwise.
24		#' @keywords internal
25		#' @noRd
26		assert_odd <- function(x) {
27	5x	arg <- deparse(substitute(x))
28	5x	if (any(round(x) %% 2 == 0)) {
29	1x	msg <- sprintf(tr_("%s must be an odd integer (%g)."), sQuote(arg), x)
30	1x	stop(msg, call. = FALSE)
31		}
32	4x	invisible(x)
33		}
34
35		#' Check Object Length
36		#'
37		#' @param x An object to be checked.
38		#' @param expected An appropriate expected value.
39		#' @return
40		#' Throws an error, if any, and returns `x` invisibly otherwise.
41		#' @keywords internal
42		#' @noRd
43		assert_length <- function(x, expected) {
44	28x	arg <- deparse(substitute(x))
45	28x	if (length(x) != expected) {
46	1x	str <- tr_("%s must be of length %d; not %d.")
47	1x	msg <- sprintf(str, sQuote(arg), expected, length(x))
48	1x	stop(msg, call. = FALSE)
49		}
50	27x	invisible(x)
51		}
52
53		assert_Matrix <- function() {
54	!	if (!requireNamespace("Matrix", quietly = TRUE)) {
55	!	msg <- tr_("The Matrix package is required. Please install it.")
56	!	stop(msg, call. = FALSE)
57		}
58	!	invisible()
59		}

1		# FIND PEAKS
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname peaks_find
7		#' @aliases peaks_find,numeric,numeric-method
8		setMethod(
9		f = "peaks_find",
10		signature = signature(x = "numeric", y = "numeric"),
11		definition = function(x, y, method = "MAD", SNR = 2, m = NULL, ...) {
12		## Validation
13	!	assert_length(y, length(x))
14	!	method <- match.arg(method, several.ok = FALSE)
15	!	if (is.null(m)) {
16	!	m <- as.integer(length(x) * 0.05)
17	!	if (m %% 2 == 0) m <- m + 1
18		}
19
20		## Noise threshold
21	!	threshold <- NULL
22	!	if (SNR != 0) {
23	!	noise <- switch (
24	!	method,
25	!	MAD = stats::mad
26		)
27	!	threshold <- noise(y, ...) * SNR
28	!	index_noise <- y < threshold
29	!	y[index_noise] <- 0
30		}
31
32		## Windows
33	!	shape <- diff(sign(diff(y, na.pad = FALSE)))
34	!	win <- window_sliding(length(x), m, i = which(shape < 0) + 1L)
35
36		## Peaks detection
37	!	pks <- vapply(
38	!	X = win,
39	!	FUN = function(w, k, data) {
40	!	i <- length(w) - k # Middle of the window
41	!	p <- if (all(data[w[-i]] <= data[w[i]])) w[i] else 0
42	!	return(p)
43		},
44	!	FUN.VALUE = numeric(1),
45	!	k = (m - 1) / 2,
46	!	data = y
47		)
48
49	!	xy <- list(x = x[pks], y = y[pks])
50	!	attr(xy, "noise") <- threshold
51	!	xy
52		}
53		)
54
55		#' @export
56		#' @rdname peaks_find
57		#' @aliases peaks_find,ANY,missing-method
58		setMethod(
59		f = "peaks_find",
60		signature = signature(x = "ANY", y = "missing"),
61		definition = function(x, method = "MAD", SNR = 2, m = NULL, ...) {
62	!	xy <- grDevices::xy.coords(x)
63	!	methods::callGeneric(x = xy$x, y = xy$y, method = method, SNR = SNR,
64	!	m = m, ...)
65		}
66		)
67
68		#' @export
69		#' @rdname peaks_fwhm
70		#' @aliases peaks_fwhm,numeric,numeric-method
71		setMethod(
72		f = "peaks_fwhm",
73		signature = signature(x = "numeric", y = "numeric"),
74		definition = function(x, y, center) {
75	1x	assert_length(y, length(x))
76
77	1x	i <- which_nearest(x, center)
78	1x	peak_height <- y[i]
79	1x	half_max <- peak_height / 2
80
81	1x	scale_for_roots <- y - half_max
82	1x	root_indices <- which(diff(sign(scale_for_roots)) != 0)
83
84	1x	tmp <- sort(c(root_indices, i))
85	1x	k <- which(tmp == i)
86
87	1x	root_left <- root_indices[k - 1]
88	1x	root_right <- root_indices[k]
89
90	1x	HWHM_left <- x[i] - x[root_left]
91	1x	HWHM_right <- x[root_right] - x[i]
92
93		# xy <- list(x = c(x[root_left], x[root_right]), y = c(half_max, half_max))
94
95	1x	FWHM <- 2 * min(c(HWHM_left, HWHM_right))
96	1x	return(FWHM)
97		}
98		)
99
100		#' @export
101		#' @rdname peaks_fwhm
102		#' @aliases peaks_fwhm,ANY,missing-method
103		setMethod(
104		f = "peaks_fwhm",
105		signature = signature(x = "ANY", y = "missing"),
106		definition = function(x, center) {
107	1x	xy <- grDevices::xy.coords(x)
108	1x	methods::callGeneric(x = xy$x, y = xy$y, center = center)
109		}
110		)

1		# SIGNAL PROCESSING
2		#' @include AllGenerics.R
3		NULL
4
5		# Bind =========================================================================
6		#' @export
7		#' @rdname signal_bind
8		#' @aliases signal_bind,ANY-method
9		setMethod(
10		f = "signal_bind",
11		signature = "ANY",
12		definition = function(...) {
13
14	!	signal <- list(...)
15	!	coords <- lapply(X = signal, FUN = grDevices::xy.coords)
16
17	!	y <- lapply(X = coords, FUN = `[[`, i = "y")
18	!	ny <- lengths(y)
19	!	if (any(ny != mean(ny))) {
20	!	msg <- tr_("All objects must have the same number of data points.")
21	!	stop(msg, call. = FALSE)
22		}
23	!	y <- do.call(rbind, y)
24
25		## Get names
26	!	subst <- substitute(list(...))[-1]
27	!	arg_names <- vapply(X = subst, FUN = deparse, FUN.VALUE = character(1))
28	!	rownames(y) <- arg_names
29	!	y
30		}
31		)
32
33		# Mean =========================================================================
34		#' @export
35		#' @rdname signal_mean
36		#' @aliases signal_mean,ANY-method
37		setMethod(
38		f = "signal_mean",
39		signature = "ANY",
40		definition = function(...) {
41
42	!	signal <- list(...)
43	!	coords <- lapply(X = signal, FUN = grDevices::xy.coords)
44
45	!	x <- lapply(X = coords, FUN = `[[`, i = "x")
46	!	nx <- lengths(x)
47
48	!	y <- lapply(X = coords, FUN = `[[`, i = "y")
49	!	ny <- lengths(y)
50
51	!	if (any(nx != mean(nx)) \|\| any(ny != mean(ny))) {
52	!	msg <- tr_("All objects must have the same number of data points.")
53	!	stop(msg, call. = FALSE)
54		}
55
56	!	x <- do.call(cbind, x)
57	!	x <- rowMeans(x)
58	!	y <- do.call(cbind, y)
59	!	y <- rowMeans(y)
60
61	!	xy <- list(x = x, y = y)
62	!	xy
63		}
64		)
65
66		# Subset =======================================================================
67		#' @export
68		#' @rdname subset
69		#' @aliases signal_select,numeric,numeric-method
70		setMethod(
71		f = "signal_select",
72		signature = c(x = "numeric", y = "numeric"),
73		definition = function(x, y, from, to) {
74	1x	assert_length(y, length(x))
75
76		## Find the nearest value
77	1x	from <- x[which_nearest(x, from)]
78	1x	to <- x[which_nearest(x, to)]
79	1x	idx <- which(x >= from & x <= to)
80
81	1x	xy <- list(x = x[idx], y = y[idx])
82	1x	xy
83		}
84		)
85
86		#' @export
87		#' @rdname subset
88		#' @aliases signal_select,ANY,missing-method
89		setMethod(
90		f = "signal_select",
91		signature = c(x = "ANY", y = "missing"),
92		definition = function(x, from, to) {
93	1x	xy <- grDevices::xy.coords(x)
94	1x	methods::callGeneric(x = xy$x, y = xy$y, from = from, to = to)
95		}
96		)
97
98		#' @export
99		#' @rdname subset
100		#' @aliases signal_slice,numeric,numeric-method
101		setMethod(
102		f = "signal_slice",
103		signature = c(x = "numeric", y = "numeric"),
104		definition = function(x, y, subset) {
105	1x	assert_length(y, length(x))
106
107	1x	i <- as.integer(subset)
108
109	1x	if (!all(i > 0) && !all(i < 0)) {
110	!	msg <- tr_("A vector of strictly positive or negative integers is expected.")
111	!	stop(msg, call. = FALSE)
112		}
113
114	1x	xy <- list(x = x[i], y = y[i])
115	1x	xy
116		}
117		)
118
119		#' @export
120		#' @rdname subset
121		#' @aliases signal_slice,ANY,missing-method
122		setMethod(
123		f = "signal_slice",
124		signature = c(x = "ANY", y = "missing"),
125		definition = function(x, subset) {
126	1x	xy <- grDevices::xy.coords(x)
127	1x	methods::callGeneric(x = xy$x, y = xy$y, subset = subset)
128		}
129		)
130
131		# Shift ========================================================================
132		#' @export
133		#' @rdname signal_shift
134		#' @aliases signal_shift,numeric,numeric-method
135		setMethod(
136		f = "signal_shift",
137		signature = c(x = "numeric", y = "numeric"),
138		definition = function(x, y, lag) {
139	1x	assert_length(y, length(x))
140
141	1x	x <- x + lag
142	1x	xy <- list(x = x, y = y)
143	1x	xy
144		}
145		)
146
147		#' @export
148		#' @rdname signal_shift
149		#' @aliases signal_shift,ANY,missing-method
150		setMethod(
151		f = "signal_shift",
152		signature = c(x = "ANY", y = "missing"),
153		definition = function(x, lag) {
154	1x	xy <- grDevices::xy.coords(x)
155	1x	methods::callGeneric(x = xy$x, y = xy$y, lag = lag)
156		}
157		)
158
159		# Drift ========================================================================
160		#' @export
161		#' @rdname signal_drift
162		#' @aliases signal_drift,numeric,numeric,numeric-method
163		setMethod(
164		f = "signal_drift",
165		signature = c(x = "numeric", y = "numeric", lag = "numeric"),
166		definition = function(x, y, lag) {
167	2x	assert_length(y, length(x))
168
169	2x	y <- y + lag
170	2x	xy <- list(x = x, y = y)
171	2x	xy
172		}
173		)
174
175		#' @export
176		#' @rdname signal_drift
177		#' @aliases signal_drift,ANY,missing,ANY-method
178		setMethod(
179		f = "signal_drift",
180		signature = c(x = "ANY", y = "missing", lag = "ANY"),
181		definition = function(x, lag, subtract = FALSE) {
182	2x	xy <- grDevices::xy.coords(x)
183	2x	zz <- grDevices::xy.coords(lag)
184	2x	lag <- if (subtract) -zz$y else zz$y
185	2x	methods::callGeneric(x = xy$x, y = xy$y, lag = lag)
186		}
187		)
188
189		# Correct ======================================================================
190		#' @export
191		#' @rdname signal_correct
192		#' @aliases signal_correct,numeric,numeric-method
193		setMethod(
194		f = "signal_correct",
195		signature = c(x = "numeric", y = "numeric"),
196		definition = function(x, y, method = c("linear", "polynomial", "asls",
197		"rollingball", "rubberband",
198		"SNIP", "4S"), ...) {
199		## Validation
200	!	method <- match.arg(method, several.ok = FALSE)
201
202		## Get method
203	!	f <- switch(
204	!	method,
205	!	asls = baseline_asls,
206	!	linear = baseline_linear,
207	!	polynomial = baseline_polynomial,
208	!	rollingball = baseline_rollingball,
209	!	rubberband = baseline_rubberband,
210	!	SNIP = baseline_snip,
211	!	`4S` = baseline_peakfilling
212		)
213
214		## Baseline estimation
215	!	bsl <- f(x = x, y = y, ...)
216
217	!	xy <- list(x = x, y = y - bsl$y)
218	!	attr(xy, "method") <- method
219	!	xy
220		}
221		)
222
223		#' @export
224		#' @rdname signal_correct
225		#' @aliases signal_correct,ANY,missing-method
226		setMethod(
227		f = "signal_correct",
228		signature = c(x = "ANY", y = "missing"),
229		definition = function(x, method = c("linear", "polynomial", "asls",
230		"rollingball", "rubberband",
231		"SNIP", "4S"), ...) {
232	!	xy <- grDevices::xy.coords(x)
233	!	methods::callGeneric(x = xy$x, y = xy$y, method = method, ...)
234		}
235		)

1		# XRD
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname ka2_strip_penalized
7		#' @aliases ka2_strip_penalized,numeric,numeric-method
8		setMethod(
9		f = "ka2_strip_penalized",
10		signature = signature(x = "numeric", y = "numeric"),
11		definition = function(x, y, lambda, wave = c(1.54060, 1.54443), tau = 0.5,
12		nseg = 1, progress = interactive()) {
13	!	assert_Matrix()
14	!	assert_length(y, length(x))
15
16		## Calculate doublet distance
17	!	wl1 <- wave[1]
18	!	wl2 <- wave[2]
19	!	wr <- wl2 / wl1
20	!	d2r <- pi / 180
21	!	u <- sin(d2r * x / 2) / wl1
22	!	delta <- 2 * asin(wr * sin(d2r * x / 2)) / d2r - x
23
24	!	m <- length(y)
25	!	n <- length(lambda)
26	!	if (nseg == 1) nseg <- signif(m * 0.50, 1)
27
28		## Empty matrices to store the results
29	!	aic <- phi <- psi <- rep(0, n)
30	!	Yhat <- Yhat2 <- matrix(data = 0, nrow = m, ncol = n)
31	!	mu <- matrix(data = 0, nrow = m, ncol = n)
32	!	sr <- matrix(data = 0, nrow = m, ncol = n)
33
34		## Prepare model matrices
35	!	xr <- max(x)
36	!	xl <- min(x) - 1
37	!	xd <- x - delta
38	!	dx <- (xr - xl) / nseg
39	!	knots <- seq(xl - 3 * dx, xr + 3 * dx, by = dx)
40
41	!	B1 <- sbase(x, lower = xl, upper = xr, n = nseg, deg = 3)
42	!	B2 <- sbase(xd, lower = xl, upper = xr, n = nseg, deg = 3)
43
44	!	B <- Matrix::rbind2(B1, B2)
45	!	Cb <- Matrix::Diagonal(m)
46	!	C <- Matrix::cbind2(Cb, tau * Cb)
47
48		## Prepare penalty
49	!	E <- Matrix::Diagonal(ncol(B1))
50	!	D <- Matrix::diff(E, diff = 2)
51
52		## Initialize
53	!	bst <- Matrix::solve(
54	!	Matrix::crossprod(B1, B1) + 1e-4 * E,
55	!	Matrix::crossprod(B1, log((y + 1) / 2))
56		)
57
58	!	if (progress) pb <- utils::txtProgressBar(min = 0, max = n, style = 3)
59
60		## Estimate for all lambda
61	!	n_lambda <- seq_len(n)
62	!	for (i in n_lambda) {
63	!	M <- penalized_strip_ka2(y, B, C, D, bst, lambda[i])
64	!	aic[i] <- M$aic
65	!	Yhat[, i] <- M$yhat
66	!	Yhat2[, i] <- M$yhat2
67	!	sr[,i] <- M$res
68	!	mu[, i] <- M$mu
69	!	phi[i] <- M$phi
70	!	psi[i] <- M$psi
71	!	if (progress) utils::setTxtProgressBar(pb, i)
72		}
73
74	!	if (progress) close(pb)
75
76		## Find optimal model and report
77	!	k <- which.min(aic)
78	!	mu <- mu[, k]
79	!	op_lamb <- lambda[k]
80	!	op_aic <- aic[k]
81	!	yhat <- Yhat[, k]
82	!	yhat2 <- Yhat2[, k]
83	!	sr <- sr[, k]
84
85	!	xy <- list(x = x, y = yhat, ka2 = yhat2, smoothed = mu)
86	!	xy
87		}
88		)
89
90		#' @export
91		#' @rdname ka2_strip_penalized
92		#' @aliases ka2_strip_penalized,ANY,missing-method
93		setMethod(
94		f = "ka2_strip_penalized",
95		signature = signature(x = "ANY", y = "missing"),
96		definition = function(x, lambda, wave = c(1.54060, 1.54443), tau = 0.5,
97		nseg = 1, progress = interactive()) {
98	!	xy <- grDevices::xy.coords(x)
99	!	methods::callGeneric(x = xy$x, y = xy$y, lambda = lambda, wave = wave,
100	!	tau = tau, nseg = nseg, progress = progress)
101		}
102		)
103
104
105		#' Truncated p-th Power
106		#'
107		#' @param x A [`numeric`] vector on which the basis is calculated.
108		#' @param knot A [`numeric`] scalar giving the truncation point.
109		#' @param p A [`numeric`] scalar specifying the power for the basis
110		#' (e.g. \eqn{p = 3} for cubic TPF).
111		#' @return
112		#' A [`numeric`] vector (piece-wise defined basis functions for \eqn{x > t}).
113		#' @author P. Eilers
114		#' @note
115		#' Slightly modified from P. Eilers' [JOPS::tpower()].
116		#' @references
117		#' Eilers, P. H. C. and Marx, B. D. (2021). *Practical Smoothing, The Joys of
118		#' P-splines.* Cambridge: Cambridge University Press.
119		#' @keywords internal
120		#' @noRd
121		tpower <- function(x, knot, p) {
122	!	(x - knot) ^ p * (x > knot)
123		}
124
125		#' B-spline Design Matrix
126		#'
127		#' Computes a B-spline design matrix using evenly spaced knots.
128		#' @param x A [`numeric`] vector of values at which the B-spline basis functions
129		#' are to be evaluated.
130		#' @param lower A [`numeric`] vector giving the lower limit of the domain of `x`.
131		#' @param upper A [`numeric`] vector giving the upper limit of the domain of `x`.
132		#' @param n A [`numeric`] scalar specifying the number of equally sized segments
133		#' between `lower` and `upper`.
134		#' @param deg A [`numeric`] scalar specifyingthe degree of the splines,
135		#' usually 1, 2, or 3 (the default).
136		#' @return
137		#' A `matrix` with `length(x)` rows and `n + deg` columns.
138		#' @note
139		#' Slightly modified from P. Eilers and B. Marx's [JOPS::bbase()].
140		#' @author P. Eilers, B. Marx
141		#' @references
142		#' Eilers, P. H. C. and Marx, B. D. (2021). *Practical Smoothing, The Joys of
143		#' P-splines.* Cambridge: Cambridge University Press.
144		#' @keywords internal
145		#' @noRd
146		bbase <- function(x, lower = min(x), upper = max(x), n = 10, deg = 3) {
147		## Check domain and adjust it if necessary
148	!	if (lower > min(x)) {
149	!	lower <- min(x)
150	!	warning(sprintf(tr_("Left boundary adjusted to %g."), lower), call. = FALSE)
151		}
152	!	if (upper < max(x)) {
153	!	upper <- max(x)
154	!	warning(sprintf(tr_("Right boundary adjusted to %g."), upper), call. = FALSE)
155		}
156
157		## Function for B-spline basis
158	!	dx <- (upper - lower) / n
159	!	knots <- seq(lower - deg * dx, upper + deg * dx, by = dx)
160	!	P <- outer(x, knots, tpower, deg)
161	!	n <- dim(P)[[2L]]
162	!	D <- diff(diag(n), diff = deg + 1) / (gamma(deg + 1) * dx ^ deg)
163	!	B <- (-1) ^ (deg + 1) * P %*% t(D)
164
165		## Make B-splines exactly zero beyond their end knots
166	!	nb <- ncol(B)
167	!	sk <- knots[(1:nb) + deg + 1]
168	!	mask <- outer(x, sk, "<")
169	!	B <- B * mask
170
171	!	B
172		}
173
174		#' Sparse B-spline Basis Matrix
175		#'
176		#' Computes a B-spline basis matrix using evenly spaced knots.
177		#' @param x A [`numeric`] vector of values at which the B-spline basis functions
178		#' are to be evaluated.
179		#' @param lower A [`numeric`] vector giving the lower limit of the domain of `x`.
180		#' @param upper A [`numeric`] vector giving the upper limit of the domain of `x`.
181		#' @param n A [`numeric`] scalar specifying the number of equally sized segments
182		#' between `lower` and `upper`.
183		#' @param deg A [`numeric`] scalar specifyingthe degree of the splines,
184		#' usually 1, 2, or 3 (the default).
185		#' @return
186		#' A [sparse matrix][Matrix::sparseMatrix()] with `length(x)` rows and
187		#' `n + deg` columns.
188		#' @note
189		#' Slightly modified from P. Eilers [JOPS::spbase()].
190		#' @author P. Eilers, B. Marx
191		#' @references
192		#' Eilers, P. H. C. and Marx, B. D. (2021). *Practical Smoothing, The Joys of
193		#' P-splines.* Cambridge: Cambridge University Press.
194		#' @keywords internal
195		#' @noRd
196		sbase <- function(x, lower = min(x), upper = max(x), n = 10, deg = 3) {
197		## Check domain and adjust it if necessary
198	!	if (lower > min(x)) {
199	!	lower <- min(x)
200	!	warning(sprintf(tr_("Left boundary adjusted to %g."), lower), call. = FALSE)
201		}
202	!	if (upper < max(x)) {
203	!	upper <- max(x)
204	!	warning(sprintf(tr_("Right boundary adjusted to %g."), upper), call. = FALSE)
205		}
206
207		## Reduce x to first interval between knots
208	!	dx <- (upper - lower) / n
209	!	ix <- floor((x - lower) / (1.0000001 * dx))
210	!	xx <- (x - lower) - ix * dx
211
212		## Full basis for reduced x
213	!	Br <- bbase(xx, lower = 0, upper = 0 + dx, n = 1, deg = deg)
214
215		## Compute proper rows, columns
216	!	m <- length(x)
217	!	p <- ncol(Br)
218	!	i_row <- rep(seq_len(m), each = p)
219	!	i_col <- rep(seq_len(p), times = m) + rep(ix, each = p)
220
221		## Sparse matrix
222	!	b <- as.vector(t(Br))
223	!	Bs <- Matrix::sparseMatrix(i = i_row, j = i_col, x = b, dims = c(m, n + deg))
224	!	Bs <- Bs[, -ncol(Bs), drop = FALSE] # ???
225
226	!	Bs
227		}
228
229		#' @param y description
230		#' @param B description
231		#' @param C description
232		#' @param D description
233		#' @param beta description
234		#' @param lambda description
235		#' @author J. J. de Rooi et al. (original R code).
236		#' @references
237		#' de Rooi, J. J., van der Pers, N. M., Hendrikx, R. W. A., Delhez, R.,
238		#' Böttger A. J. and Eilers, P. H. C. (2014). Smoothing of X-ray diffraction
239		#' data and Ka2 elimination using penalized likelihood and the composite link
240		#' model. Journal of Applied Crystallography, 47: 852-860.
241		#' \doi{10.1107/S1600576714005809}
242		#' @keywords internal
243		#' @noRd
244		penalized_strip_ka2 <- function(y, B, C, D, beta, lambda) {
245	!	m <- length(y)
246
247		## Penalty
248	!	P <- lambda * Matrix::t(D) %*% D
249
250		## Fit iteratively
251	!	delta_beta <- 1
252	!	limit <- 1 # Stop condition to prevent infinite loop
253	!	while (delta_beta >= 0.0001 && limit <= 30) {
254	!	eta <- B %*% beta
255	!	gam <- exp(eta)
256	!	mu <- as.numeric(C %*% gam)
257
258	!	M <- Matrix::Diagonal(x = 1 / mu)
259	!	X <- (M %% C %% Matrix::Diagonal(x = as.numeric(gam))) %*% B
260	!	W <- Matrix::Diagonal(x = mu)
261	!	G <- Matrix::crossprod(X, W) %*% X
262
263	!	new_beta <- Matrix::solve(G + P, Matrix::crossprod(X, (y - mu)) + G %*% beta)
264	!	delta_beta <- max(abs(new_beta - beta))
265	!	beta <- new_beta
266
267	!	limit <- limit + 1
268		}
269
270		## AIC
271	!	H <- Matrix::solve(G + P, G)
272	!	ed <- sum(Matrix::diag(H))
273	!	dev <- 2 * sum(y * log((y + (y == 0)) / mu))
274	!	aic <- dev + 2 * ed
275
276		## Standardized residuals for counts
277	!	sr <- (y - mu) / sqrt(mu)
278	!	sdr <- sqrt(sum(sr^2) / (m - ed))
279
280		## L-curve
281	!	phi <- log(sum((D %*% beta)^2))
282	!	psi <- log(sum(sr^2))
283
284		## estimated signals
285	!	yhat <- gam[1:m]
286	!	yhat2 <- gam[(m + 1):length(gam)] * 0.5
287
288		## return values
289	!	list(yhat = yhat, yhat2 = yhat2, mu = mu, aic = aic, res = sr,
290	!	beta = beta, phi = phi, psi = psi)
291		}

1		# REPLACE
2		#' @include AllGenerics.R
3		NULL
4
5		# Threshold ====================================================================
6		#' @export
7		#' @rdname replace_threshold
8		#' @aliases replace_threshold,numeric,numeric,function-method
9		setMethod(
10		f = "replace_threshold",
11		signature = signature(x = "numeric", y = "numeric", threshold = "function"),
12		definition = function(x, y, threshold, value = 0, ...) {
13	1x	assert_length(y, length(x))
14
15	1x	threshold <- threshold(y, ...)
16	1x	methods::callGeneric(x, y, threshold = threshold, value = value)
17		}
18		)
19
20		#' @export
21		#' @rdname replace_threshold
22		#' @aliases replace_threshold,ANY,missing,function-method
23		setMethod(
24		f = "replace_threshold",
25		signature = signature(x = "ANY", y = "missing", threshold = "function"),
26		definition = function(x, threshold, value = 0, ...) {
27	1x	xy <- grDevices::xy.coords(x)
28	1x	methods::callGeneric(x = xy$x, y = xy$y, threshold = threshold,
29	1x	value = value)
30		}
31		)
32
33		#' @export
34		#' @rdname replace_threshold
35		#' @aliases replace_threshold,numeric,numeric,numeric-method
36		setMethod(
37		f = "replace_threshold",
38		signature = signature(x = "numeric", y = "numeric", threshold = "numeric"),
39		definition = function(x, y, threshold, value = 0, ...) {
40	3x	assert_length(y, length(x))
41
42	3x	i <- y < threshold
43	3x	y[i] <- value
44
45	3x	xy <- list(x = x, y = y)
46	3x	xy
47		}
48		)
49
50		#' @export
51		#' @rdname replace_threshold
52		#' @aliases replace_threshold,ANY,missing,numeric-method
53		setMethod(
54		f = "replace_threshold",
55		signature = signature(x = "ANY", y = "missing", threshold = "numeric"),
56		definition = function(x, threshold, value = 0, ...) {
57	1x	xy <- grDevices::xy.coords(x)
58	1x	methods::callGeneric(x = xy$x, y = xy$y, threshold = threshold,
59	1x	value = value)
60		}
61		)
62
63		# Negative =====================================================================
64		#' @export
65		#' @rdname replace_negative
66		#' @aliases replace_negative,numeric,numeric-method
67		setMethod(
68		f = "replace_negative",
69		signature = signature(x = "numeric", y = "numeric"),
70		definition = function(x, y, value = 0) {
71	1x	replace_threshold(x, y, threshold = 0, value = value)
72		}
73		)
74
75		#' @export
76		#' @rdname replace_negative
77		#' @aliases replace_negative,ANY,missing-method
78		setMethod(
79		f = "replace_negative",
80		signature = signature(x = "ANY", y = "missing"),
81		definition = function(x, value = 0) {
82	1x	xy <- grDevices::xy.coords(x)
83	1x	methods::callGeneric(x = xy$x, y = xy$y, value = value)
84		}
85		)

1		# SCALING
2		#' @include AllGenerics.R
3		NULL
4
5		# Area =========================================================================
6		#' @export
7		#' @rdname rescale_area
8		#' @aliases rescale_area,numeric,numeric-method
9		setMethod(
10		f = "rescale_area",
11		signature = signature(x = "numeric", y = "numeric"),
12		definition = function(x, y, method = c("rectangle", "trapezoid"), ...) {
13		## Validation
14	1x	assert_length(y, length(x))
15	1x	method <- match.arg(method, several.ok = FALSE)
16
17		## Get method
18	1x	f <- switch(
19	1x	method,
20	1x	rectangle = integrate_rectangle,
21	1x	trapezoid = integrate_trapezoid
22		)
23
24		## Baseline estimation
25	1x	auc <- f(x = x, y = y, ...)
26	1x	y <- y / auc
27
28	1x	xy <- list(x = x, y = y)
29	1x	xy
30		}
31		)
32
33		#' @export
34		#' @rdname rescale_area
35		#' @aliases rescale_area,ANY,missing-method
36		setMethod(
37		f = "rescale_area",
38		signature = signature(x = "ANY", y = "missing"),
39		definition = function(x, method = c("rectangle", "trapezoid"), ...) {
40	1x	xy <- grDevices::xy.coords(x)
41	1x	methods::callGeneric(x = xy$x, y = xy$y, method = method, ...)
42		}
43		)
44
45		# Total ========================================================================
46		#' @export
47		#' @rdname rescale_total
48		#' @aliases rescale_total,numeric,numeric-method
49		setMethod(
50		f = "rescale_total",
51		signature = signature(x = "numeric", y = "numeric"),
52		definition = function(x, y, total = 1) {
53	1x	assert_length(y, length(x))
54
55	1x	y <- (y * total) / sum(abs(y))
56	1x	xy <- list(x = x, y = y)
57	1x	xy
58		}
59		)
60
61		#' @export
62		#' @rdname rescale_total
63		#' @aliases rescale_total,ANY,missing-method
64		setMethod(
65		f = "rescale_total",
66		signature = signature(x = "ANY", y = "missing"),
67		definition = function(x, total = 1) {
68	1x	xy <- grDevices::xy.coords(x)
69	1x	methods::callGeneric(x = xy$x, y = xy$y, total = total)
70		}
71		)
72
73		# Min-Max ======================================================================
74		#' @export
75		#' @rdname rescale_range
76		#' @aliases rescale_range,numeric,numeric-method
77		setMethod(
78		f = "rescale_range",
79		signature = signature(x = "numeric", y = "numeric"),
80		definition = function(x, y, min = 0, max = 1) {
81		## Validation
82	5x	assert_length(y, length(x))
83	5x	if (min > max) {
84	1x	msg <- sprintf(tr_("%s (%g) must be lower than %s (%g)."),
85	1x	sQuote("min"), min, sQuote("max"), max)
86	1x	stop(msg, call. = FALSE)
87		}
88
89	4x	y <- (y - min(y)) / (max(y) - min(y)) * (max - min) + min
90	4x	xy <- list(x = x, y = y)
91	4x	xy
92		}
93		)
94
95		#' @export
96		#' @rdname rescale_range
97		#' @aliases rescale_range,ANY,missing-method
98		setMethod(
99		f = "rescale_range",
100		signature = signature(x = "ANY", y = "missing"),
101		definition = function(x, min = 0, max = 1) {
102	1x	xy <- grDevices::xy.coords(x)
103	1x	methods::callGeneric(x = xy$x, y = xy$y, min = min, max = max)
104		}
105		)
106
107		#' @export
108		#' @rdname rescale_range
109		#' @aliases rescale_min,numeric,numeric-method
110		setMethod(
111		f = "rescale_min",
112		signature = signature(x = "numeric", y = "numeric"),
113		definition = function(x, y, min = 0) {
114	2x	rescale_range(x, y, min = min, max = max(y))
115		}
116		)
117
118		#' @export
119		#' @rdname rescale_range
120		#' @aliases rescale_min,ANY,missing-method
121		setMethod(
122		f = "rescale_min",
123		signature = signature(x = "ANY", y = "missing"),
124		definition = function(x, min = 0) {
125	2x	xy <- grDevices::xy.coords(x)
126	2x	methods::callGeneric(x = xy$x, y = xy$y, min = min)
127		}
128		)
129
130		#' @export
131		#' @rdname rescale_range
132		#' @aliases rescale_max,numeric,numeric-method
133		setMethod(
134		f = "rescale_max",
135		signature = signature(x = "numeric", y = "numeric"),
136		definition = function(x, y, max = 1) {
137	2x	rescale_range(x, y, min = min(y), max = max)
138		}
139		)
140
141		#' @export
142		#' @rdname rescale_range
143		#' @aliases rescale_max,ANY,missing-method
144		setMethod(
145		f = "rescale_max",
146		signature = signature(x = "ANY", y = "missing"),
147		definition = function(x, max = 1) {
148	2x	xy <- grDevices::xy.coords(x)
149	2x	methods::callGeneric(x = xy$x, y = xy$y, max = max)
150		}
151		)
152
153		# Transform ====================================================================
154		#' @export
155		#' @rdname rescale_transform
156		#' @aliases rescale_transform,numeric,numeric-method
157		setMethod(
158		f = "rescale_transform",
159		signature = signature(x = "numeric", y = "numeric"),
160		definition = function(x, y, f, ...) {
161	1x	assert_length(y, length(x))
162
163	1x	y <- f(y, ...)
164	1x	xy <- list(x = x, y = y)
165	1x	xy
166		}
167		)
168
169		#' @export
170		#' @rdname rescale_transform
171		#' @aliases rescale_transform,ANY,missing-method
172		setMethod(
173		f = "rescale_transform",
174		signature = signature(x = "ANY", y = "missing"),
175		definition = function(x, f, ...) {
176	1x	xy <- grDevices::xy.coords(x)
177	1x	methods::callGeneric(x = xy$x, y = xy$y, f = f, ...)
178		}
179		)
180
181		# SNV ==========================================================================
182		#' @export
183		#' @rdname rescale_snv
184		#' @aliases rescale_snv,numeric,numeric-method
185		setMethod(
186		f = "rescale_snv",
187		signature = signature(x = "numeric", y = "numeric"),
188		definition = function(x, y, ...) {
189		## Validation
190	1x	assert_length(y, length(x))
191
192	1x	y <- (y - mean(y)) / stats::sd(y)
193
194	1x	xy <- list(x = x, y = y)
195	1x	xy
196		}
197		)
198
199		#' @export
200		#' @rdname rescale_snv
201		#' @aliases rescale_snv,ANY,missing-method
202		setMethod(
203		f = "rescale_snv",
204		signature = signature(x = "ANY", y = "missing"),
205		definition = function(x, ...) {
206	1x	xy <- grDevices::xy.coords(x)
207	1x	methods::callGeneric(x = xy$x, y = xy$y, ...)
208		}
209		)

1		# RESAMPLE
2		#' @include AllGenerics.R
3		NULL
4
5		# Bin ==========================================================================
6		#' @export
7		#' @rdname resample_bin
8		#' @aliases resample_bin,numeric,numeric-method
9		setMethod(
10		f = "resample_bin",
11		signature = signature(x = "numeric", y = "numeric"),
12		definition = function(x, y, by, f = mean, ...) {
13	1x	assert_length(y, length(x))
14
15	1x	k <- length(x) %% by
16	1x	if (k != 0) {
17	1x	x <- utils::tail(x, -k)
18	1x	y <- utils::tail(y, -k)
19		}
20	1x	i <- cut(seq_along(x), breaks = length(x) / by, labels = FALSE)
21
22	1x	mid <- tapply(X = x, INDEX = i, FUN = mean, simplify = FALSE)
23	1x	bin <- tapply(X = y, INDEX = i, FUN = f, ..., simplify = FALSE)
24
25	1x	xy <- list(x = unlist(mid), y = unlist(bin))
26	1x	xy
27		}
28		)
29
30		#' @export
31		#' @rdname resample_bin
32		#' @aliases resample_bin,ANY,missing-method
33		setMethod(
34		f = "resample_bin",
35		signature = signature(x = "ANY", y = "missing"),
36		definition = function(x, y, by, f = sum) {
37	1x	xy <- grDevices::xy.coords(x)
38	1x	methods::callGeneric(x = xy$x, y = xy$y, by = by, f = f)
39		}
40		)
41
42		# Downsample ===================================================================
43		#' @export
44		#' @rdname resample_down
45		#' @aliases resample_down,numeric,numeric-method
46		setMethod(
47		f = "resample_down",
48		signature = signature(x = "numeric", y = "numeric"),
49		definition = function(x, y, by) {
50	1x	assert_length(y, length(x))
51
52	1x	i <- seq(from = 1, to = length(x), by = by)
53	1x	xy <- list(x = x[i], y = y[i])
54	1x	xy
55		}
56		)
57
58		#' @export
59		#' @rdname resample_down
60		#' @aliases resample_down,ANY,missing-method
61		setMethod(
62		f = "resample_down",
63		signature = signature(x = "ANY", y = "missing"),
64		definition = function(x, y, by) {
65	1x	xy <- grDevices::xy.coords(x)
66	1x	methods::callGeneric(x = xy$x, y = xy$y, by = by)
67		}
68		)
69
70		# Interpolate ==================================================================
71		#' @export
72		#' @rdname resample_interpolate
73		#' @aliases resample_interpolate,numeric,numeric-method
74		setMethod(
75		f = "resample_interpolate",
76		signature = signature(x = "numeric", y = "numeric"),
77		definition = function(x, y, from, to, by, ...) {
78	!	assert_length(y, length(x))
79
80		## New x scale
81	!	x_scale <- seq(from = from, to = to, by = by)
82
83		## Interpolate
84	!	new_data <- stats::approx(x = x, y = y, xout = x_scale, ...)
85
86	!	xy <- list(x = new_data$x, y = new_data$y)
87	!	xy
88		}
89		)
90
91		#' @export
92		#' @rdname resample_interpolate
93		#' @aliases resample_interpolate,ANY,missing-method
94		setMethod(
95		f = "resample_interpolate",
96		signature = signature(x = "ANY", y = "missing"),
97		definition = function(x, y, from, to, by, ...) {
98	!	xy <- grDevices::xy.coords(x)
99	!	methods::callGeneric(x = xy$x, y = xy$y, from = from, to = to, by = by, ...)
100		}
101		)

1		# SIGNAL INTEGRATION
2		#' @include AllGenerics.R
3		NULL
4
5		#' @export
6		#' @rdname integrate_rectangle
7		#' @aliases integrate_rectangle,numeric,numeric-method
8		setMethod(
9		f = "integrate_rectangle",
10		signature = signature(x = "numeric", y = "numeric"),
11		definition = function(x, y, right = FALSE) {
12	4x	assert_length(y, length(x))
13
14	4x	h <- if (right) utils::tail(y, -1) else utils::head(y, -1)
15	4x	sum(h * diff(x))
16		}
17		)
18
19		#' @export
20		#' @rdname integrate_rectangle
21		#' @aliases integrate_rectangle,ANY,missing-method
22		setMethod(
23		f = "integrate_rectangle",
24		signature = signature(x = "ANY", y = "missing"),
25		definition = function(x, right = FALSE) {
26	3x	xy <- grDevices::xy.coords(x)
27	3x	methods::callGeneric(x = xy$x, y = xy$y, right = right)
28		}
29		)
30
31		#' @export
32		#' @rdname integrate_trapezoid
33		#' @aliases integrate_trapezoid,numeric,numeric-method
34		setMethod(
35		f = "integrate_trapezoid",
36		signature = signature(x = "numeric", y = "numeric"),
37		definition = function(x, y) {
38	1x	assert_length(y, length(x))
39
40	1x	sum((utils::head(y, -1) + utils::tail(y, -1)) * diff(x) / 2)
41		}
42		)
43
44		#' @export
45		#' @rdname integrate_trapezoid
46		#' @aliases integrate_trapezoid,ANY,missing-method
47		setMethod(
48		f = "integrate_trapezoid",
49		signature = signature(x = "ANY", y = "missing"),
50		definition = function(x) {
51	1x	xy <- grDevices::xy.coords(x)
52	1x	methods::callGeneric(x = xy$x, y = xy$y)
53		}
54		)