## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 8,
  fig.height = 3.5,
  out.width = "100%"
)
if (!requireNamespace("wavethresh", quietly = TRUE)) {
  knitr::opts_chunk$set(eval = FALSE)
  message(
    "This vignette requires the 'wavethresh' package."
  )
}

## ----pkgs---------------------------------------------------------------------
library(rLifting)
library(wavethresh)

data(BabyECG)
data(BabySS)

## ----raw-plot, fig.height = 4-------------------------------------------------
t_sec  = seq(0, by = 16, length.out = length(BabyECG))
t_hour = t_sec / 3600

sleep_col = c(
  "1" = "#2196F3", "2" = "#90CAF9", 
  "3" = "#FF9800", "4" = "#F44336"
)

plot(
  t_hour, BabyECG, type = "l", col = "grey60", lwd = 0.6,
  xlab = "Time (hours after 21:17)", ylab = "Heart rate (bpm)",
  main = "BabyECG — infant heart rate over one night"
)

for (state in 1:4) {
  idx = which(BabySS == state)
  points(
    t_hour[idx], BabyECG[idx], pch = 20, cex = 0.3,
    col = sleep_col[as.character(state)]
  )
}

legend(
  "topright",
  legend = c("Quiet sleep", "Transitional", "Active sleep", "Awake"),
  col = sleep_col, pch = 20, cex = 0.75, bty = "n"
)

## ----offline-denoise----------------------------------------------------------
ls_cdf53 = lifting_scheme("cdf53")

d_offline = denoise_signal_offline(
  BabyECG, ls_cdf53,
  threshold_method = "universal",
  shrinkage = "semisoft",
  extension = "symmetric"
)

## ----offline-plot-------------------------------------------------------------
plot(
  t_hour, BabyECG, type = "l", col = "grey70", lwd = 0.6,
  xlab = "Time (hours after 21:17)", ylab = "Heart rate (bpm)",
  main = "Offline denoising — cdf53 / universal semisoft"
)

lines(t_hour, d_offline, col = "#1565C0", lwd = 1.5)
legend(
  "topright", legend = c("Raw", "Denoised"), lwd = c(1, 2),
  col = c("grey70", "#1565C0"), bty = "n"
)

## ----offline-timing-----------------------------------------------------------
t_100 = system.time(
  for (i in seq_len(100))
    denoise_signal_offline(
      BabyECG, ls_cdf53,
      threshold_method = "universal",
      shrinkage = "semisoft",
      extension = "symmetric"
    )
)

cat(
  sprintf(
    "Per-call time: %.2f ms  (N = %d)\n",
    t_100["elapsed"] / 100 * 1000, length(BabyECG)
  )
)

## ----offline-diagnostics------------------------------------------------------
lwt_raw = lwt(BabyECG, ls_cdf53, levels = 3L, extension = "symmetric")
sigma_hat = median(abs(lwt_raw$coeffs[[1]])) / 0.6745

resid = BabyECG - d_offline
acf_lag1 = acf(resid, plot = FALSE)$acf[2]

cat(
  sprintf(
    "Estimated noise sigma (MAD, finest level): %.2f bpm\n",
    sigma_hat
  )
)

cat(
  sprintf(
    "Residual sd after denoising: %.2f bpm\n",
    sd(resid)
  )
)

cat(sprintf("Residual ACF at lag 1 : %+.3f\n", acf_lag1))

## ----irr-setup----------------------------------------------------------------
set.seed(42)
n = length(BabyECG)
t_reg = seq(0, by = 16, length.out = n)

keep = sort(sample(n, round(0.70 * n)))  # retain 70 % of samples
t_irr = t_reg[keep]
y_irr = BabyECG[keep]

cat(
  sprintf(
    "Observed: %d of %d samples (%.0f %% retained)\n",
    length(keep), n, length(keep) / n * 100
  )
)

cat(
  sprintf(
    "Gap range: %.0f – %.0f s  (median %.0f s)\n",
    min(diff(t_irr)), max(diff(t_irr)), median(diff(t_irr))
  )
)

## ----irr-rlifting-------------------------------------------------------------
d_irr = denoise_signal_offline(
  y_irr, ls_cdf53,
  threshold_method = "universal",
  shrinkage = "semisoft",
  extension = "symmetric",
  t = t_irr
)

## ----irr-interp---------------------------------------------------------------
y_interp = approx(t_irr, y_irr, xout = t_reg)$y
y_interp[is.na(y_interp)] = mean(y_irr)

d_interp_full = denoise_signal_offline(
  y_interp, ls_cdf53,
  threshold_method = "universal",
  shrinkage = "semisoft",
  extension = "symmetric"
)
d_interp_at_obs = approx(t_reg, d_interp_full, xout = t_irr)$y

## ----irr-heldout--------------------------------------------------------------
dropped = setdiff(seq_len(n), keep)
y_true_dropped = BabyECG[dropped]
t_dropped = t_reg[dropped]

d_irr_at_ho = approx(t_irr, d_irr, xout = t_dropped)$y
d_interp_at_ho = d_interp_full[dropped]
raw_at_ho = approx(t_irr, y_irr, xout = t_dropped)$y

rmse = function(a, b) sqrt(mean((a - b)^2, na.rm = TRUE))

cat(
  sprintf(
    "Held-out RMSE at %d dropped positions (bpm):\n",
    length(dropped)
  )
)

cat(
  sprintf(
    "  rLifting irregular            : %.2f\n",
    rmse(y_true_dropped, d_irr_at_ho)
  )
)

cat(
  sprintf(
    "  Interpolate + offline         : %.2f\n",
    rmse(y_true_dropped, d_interp_at_ho)
  )
)
cat(
  sprintf(
    "  Raw interpolation (no denoise): %.2f\n",
    rmse(y_true_dropped, raw_at_ho)
  )
)

## ----irr-roughness------------------------------------------------------------
rough = function(x) sd(diff(x, differences = 2), na.rm = TRUE)

cat(sprintf("Roughness (sd of 2nd diff) at observed positions:\n"))
cat(sprintf("rLifting irregular: %.3f\n", rough(d_irr)))
cat(sprintf("Interp + offline: %.3f\n", rough(d_interp_at_obs)))

## ----irr-plot-----------------------------------------------------------------
par(mfrow = c(1, 2))

t_irr_h = t_irr / 3600

plot(
  t_irr_h, y_irr, type = "l", col = "grey70", lwd = 0.5,
  xlab = "Time (h)", ylab = "Heart rate (bpm)",
  main = "rLifting — irregular grid"
)

lines(t_irr_h, d_irr, col = "#1565C0", lwd = 1.5)
legend(
  "topright", legend = c("Observed", "Denoised"),
  lwd = c(1, 2), col = c("grey70", "#1565C0"), bty = "n", cex = 0.8
)

plot(
  t_irr_h, y_irr, type = "l", col = "grey70", lwd = 0.5,
  xlab = "Time (h)", ylab = "Heart rate (bpm)",
  main = "Interpolate → denoise"
)

lines(t_irr_h, d_interp_at_obs, col = "#C62828", lwd = 1.5)
legend(
  "topright", legend = c("Observed", "Interp + denoise"),
  lwd = c(1, 2), col = c("grey70", "#C62828"), bty = "n", cex = 0.8
)

par(mfrow = c(1, 1))

## ----irr-timing---------------------------------------------------------------
t_irr_100 = system.time(
  for (i in seq_len(100))
    denoise_signal_offline(
      y_irr, ls_cdf53,
      threshold_method = "universal",
      shrinkage = "semisoft",
      extension = "symmetric", t = t_irr
    )
)
cat(
  sprintf(
    "Irregular-grid per-call: %.2f ms\n",
    t_irr_100["elapsed"] / 100 * 1000
  )
)

## ----stream-setup-------------------------------------------------------------
ls_haar = lifting_scheme("haar")

stream = new_wavelet_stream(
  ls_haar,
  extension = "symmetric",
  threshold_method = "universal",
  shrinkage = "semisoft"
)

## ----stream-run---------------------------------------------------------------
out_stream = numeric(length(BabyECG))
t_stream = system.time(
  for (i in seq_along(BabyECG))
    out_stream[i] = stream(BabyECG[i])
)

cat(
  sprintf(
    "Total for %d samples : %.1f ms\n",
    length(BabyECG), t_stream["elapsed"] * 1000
  )
)

cat(
  sprintf(
    "Per-sample latency   : %.1f µs\n",
    t_stream["elapsed"] / length(BabyECG) * 1e6
  )
)

## ----stream-convergence-------------------------------------------------------
plot(
  t_hour, BabyECG, type = "l", col = "grey70", lwd = 0.6,
  xlab = "Time (hours after 21:17)", ylab = "Heart rate (bpm)",
  main = "Streaming denoising — haar / universal semisoft"
)

lines(t_hour, out_stream, col = "#2E7D32", lwd = 1.5)
abline(v = t_hour[256], lty = 2, col = "grey40")
text(t_hour[256] + 0.05, 175, "window full\n(~68 min)", cex = 0.75, adj = 0)

legend(
  "topright",
  legend = c("Raw", "Streaming denoised", "Warm-up boundary"),
  lwd = c(1, 2, 1), lty = c(1, 1, 2),
  col = c("grey70", "#2E7D32", "grey40"),
  bty = "n", cex = 0.8
)

## ----summary-table, echo = FALSE----------------------------------------------
t_off_100 = system.time(
  for (i in seq_len(100))
    denoise_signal_offline(
      BabyECG, ls_cdf53,
      threshold_method = "universal",
      shrinkage = "semisoft",
      extension = "symmetric"
    )
)

t_cas_100 = system.time(
  for (i in seq_len(100))
    denoise_signal_causal(
      BabyECG, ls_haar,
      threshold_method = "universal", shrinkage = "semisoft",
      extension = "symmetric", window_size = 256
    )
)

latency_us = t_stream["elapsed"] / length(BabyECG) * 1e6

ms_off = round(t_off_100["elapsed"] / 100 * 1000, 2)
ms_irr = round(t_irr_100["elapsed"] / 100 * 1000, 2)
ms_cas = round(t_cas_100["elapsed"] / 100 * 1000, 2)
ms_str = round(t_stream["elapsed"] * 1000, 1)
us_off = round(ms_off / length(BabyECG) * 1000, 2)
us_irr = round(ms_irr / length(y_irr)   * 1000, 2)
us_cas = round(ms_cas / length(BabyECG) * 1000, 2)
us_str = round(latency_us, 1)

summary_tbl = data.frame(
  Mode = c(
    "Offline (regular)", "Offline (irregular)",
    "Causal", "Stream"
  ),
  N = c(
    length(BabyECG), length(y_irr),
    length(BabyECG), length(BabyECG)
  ),
  `Per-call (ms)` = c(ms_off, ms_irr, ms_cas, ms_str),
  `Per-sample (µs)` = c(us_off, us_irr, us_cas, us_str),
  Wavelet = c("cdf53", "cdf53", "haar", "haar"),
  Irregular = c("no", "yes (30 % gap)", "no", "no"),
  check.names = FALSE
)

knitr::kable(
  summary_tbl, align = "lrrrrr",
  caption = paste0(
    "Computational summary — BabyECG ",
    "(N = 2 048, 16 s sampling)"
  )
)