## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)

## ----setup--------------------------------------------------------------------
library(rLifting)

if (!requireNamespace("ggplot2", quietly = TRUE)) {
  knitr::opts_chunk$set(eval = FALSE)
  message("'ggplot2' is required to render plots. Vignette code will not run.")
} else {
  library(ggplot2)
}

## ----builtin------------------------------------------------------------------
haar = lifting_scheme("haar")
cdf53 = lifting_scheme("cdf53")
cdf97 = lifting_scheme("cdf97")
dd4 = lifting_scheme("dd4")
print(cdf97)

## ----custom_define------------------------------------------------------------
p_step = lift_step(
  type = "predict",
  coeffs = c(0.5, 0.5),
  start_idx = 0
)

u_step = lift_step(
  type = "update",
  coeffs = c(0.25, 0.25),
  start_idx = -1
)

my_cdf53 = custom_wavelet(
  name = "MyCDF53",
  steps = list(p_step, u_step),
  norm = c(sqrt(2), 1/sqrt(2))
)

print(my_cdf53)

## ----custom_use---------------------------------------------------------------
x = 1:16
# local_linear extension preserves the linear trend at the boundary,
# so all detail coefficients vanish for a degree-1 input (2 VM).
fwd = lwt(
  x, my_cdf53, 
  levels = floor(log2(length(x))),
  extension = "local_linear"
)

# CDF 5/3 has 2 VM: detail coefficients are zero for a linear signal
print(round(fwd$coeffs$d1, 10))

# Perfect reconstruction
rec = ilwt(fwd)
all.equal(as.numeric(x), rec)

## ----diagnose-----------------------------------------------------------------
config = list(
  is_ortho = FALSE, # CDF 5/3 is biorthogonal, not orthogonal
  vm_degrees = c(0, 1), # degree 0 = constant, 1 = ramp
  max_taps = 5 # expected filter support width
)

diag = diagnose_wavelet(
  my_cdf53,
  config = config,
  plot = FALSE # suppress basis plot
)

## ----configs, eval=FALSE------------------------------------------------------
# diagnose_wavelet("haar", list(is_ortho = TRUE,  vm_degrees = 0,    max_taps = 2))
# diagnose_wavelet("db2", list(is_ortho = TRUE,  vm_degrees = 0:1,  max_taps = 4))
# diagnose_wavelet("cdf53", list(is_ortho = FALSE, vm_degrees = 0:1,  max_taps = 5))
# diagnose_wavelet("dd4", list(is_ortho = FALSE, vm_degrees = 0:3,  max_taps = 11))
# diagnose_wavelet("cdf97", list(is_ortho = FALSE, vm_degrees = 0:3,  max_taps = 12))

## ----pipeline, fig.cap="Figure 1: Custom denoising pipeline on a synthetic two-frequency signal. The noisy input (grey) is decomposed with CDF 9/7 at 5 levels; details d1–d3 are denoised with semisoft, d4 with hard, and d5 is left untouched to preserve coarse structure. The reconstructed signal (orange) tracks the underlying truth (black) closely."----
set.seed(42)
n = 512
t_v = seq(0, 1, length.out = n)
pure = sin(2 * pi * 5 * t_v) + 0.5 * sin(2 * pi * 20 * t_v)
noisy = pure + rnorm(n, sd = 0.4)

sch = lifting_scheme("cdf97")

# 1. Forward transform
decomp = lwt(noisy, sch, levels = 5)

# 2. Compute adaptive thresholds
lambdas = compute_adaptive_threshold(decomp, alpha = 0.3, beta = 1.2)

# 3. Custom strategy: semisoft on fine levels, hard on d4; leave d5 untouched
for (lev in 1:3) {
  dname = paste0("d", lev)
  decomp$coeffs[[dname]] = threshold(
    decomp$coeffs[[dname]], lambdas[[dname]], method = "semisoft"
  )
}
decomp$coeffs$d4 = threshold(decomp$coeffs$d4, lambdas$d4, method = "hard")
# d5 untouched to preserve coarse structure

# 4. Inverse transform
reconstructed = ilwt(decomp)

df_pipe = data.frame(
  t = rep(t_v, 3),
  value = c(noisy, pure, reconstructed),
  Signal = rep(c("Noisy", "Truth", "Filtered"), each = n)
)

ggplot(df_pipe, aes(x = t, y = value, colour = Signal)) +
  geom_line(alpha = 0.7, linewidth = 0.4) +
  labs(
    title = "Custom denoising pipeline (CDF 9/7, 5 levels)",
    x = "Time", y = "Amplitude"
  ) +
  scale_colour_manual(
    values = c(
      "Noisy" = "grey70", "Truth" = "black",
      "Filtered" = "#D55E00"
    )
  ) +
  theme_minimal()

## ----boundary_ll_k, eval=FALSE------------------------------------------------
# # Compare 2-point vs 4-point extrapolation on a noisy boundary
# res_k2 = denoise_signal_offline(
#   noisy, sch,
#   extension = "local_linear", ll_k = 2L
# )
# 
# res_k4 = denoise_signal_offline(
#   noisy, sch,
#   extension = "local_linear", ll_k = 4L
# )

## ----boundary_compare, fig.cap="Figure 2: Left-boundary zoom of an offline denoising on a trending signal (3·t + noise + sinusoid) using CDF 9/7 at 4 levels. symmetric (blue) and local_linear (green) reflect or extrapolate the boundary; one_sided (pink) uses only existing samples. On a smooth trending signal in offline mode all three are visually similar; differences become important in causal mode."----
set.seed(42)
n_bnd = 256
t_bnd = seq(0, 1, length.out = n_bnd)
trend = 3 * t_bnd
sig_b = trend + 0.3 * rnorm(n_bnd) + 0.5 * sin(2 * pi * 8 * t_bnd)
sch_b = lifting_scheme("cdf97")

res_sym = denoise_signal_offline(
  sig_b, sch_b, levels = 4,
  extension = "symmetric"
)

res_ll = denoise_signal_offline(
  sig_b, sch_b, levels = 4,
  extension = "local_linear"
)
res_os = denoise_signal_offline(
  sig_b, sch_b, levels = 4,
  extension = "one_sided"
)

df_bnd = data.frame(
  t = rep(t_bnd, 4),
  value = c(sig_b, res_sym, res_ll, res_os),
  Signal = rep(
    c("Noisy", "symmetric", "local_linear", "one_sided"),
    each = n_bnd
  )
)

ggplot(df_bnd, aes(x = t, y = value, colour = Signal, linewidth = Signal)) +
  geom_line(alpha = 0.85) +
  scale_linewidth_manual(
    values = c(
      "Noisy" = 0.3, "symmetric" = 0.7, "local_linear" = 0.9, "one_sided" = 0.9
    )
  ) +
  scale_colour_manual(values = c(
    "Noisy" = "grey75",
    "symmetric" = "#0072B2",
    "local_linear" = "#009E73",
    "one_sided" = "#CC79A7"
  )) +
  coord_cartesian(xlim = c(0, 0.15)) +
  labs(
    title = "Left-boundary zoom: symmetric vs local_linear vs one_sided",
    subtitle = "Linearly trending signal, CDF 9/7, 4 levels",
    x = "Time", y = "Amplitude"
  ) +
  theme_minimal()

## ----boundary_causal, fig.cap="Figure 3: Causal-stream denoising of a trending signal with sinusoidal component (Haar, window = 63). symmetric (blue) reflects the most recent sample inward at each step; one_sided (pink) renormalises the filter against the available window. Both follow the trend; the differences emerge in regions where the filter window contains a transition."----
set.seed(1)
n_c = 512
t_c = seq(0, 2, length.out = n_c)
x_c = 2 * t_c + 0.4 * rnorm(n_c) + sin(2 * pi * 4 * t_c)

sch_c = lifting_scheme("haar")
proc_sym = new_wavelet_stream(
  sch_c, window_size = 63, levels = 3,
  extension = "symmetric"
)

proc_os = new_wavelet_stream(
  sch_c, window_size = 63, levels = 3,
  extension = "one_sided"
)

out_sym = out_os = numeric(n_c)
for (i in seq_len(n_c)) {
  out_sym[i] = proc_sym(x_c[i])
  out_os[i]  = proc_os(x_c[i])
}

df_c = data.frame(
  t = rep(t_c, 3),
  value = c(x_c, out_sym, out_os),
  Signal = rep(c("Noisy", "symmetric", "one_sided"), each = n_c)
)

ggplot(df_c, aes(x = t, y = value, colour = Signal, linewidth = Signal)) +
  geom_line(alpha = 0.85) +
  scale_linewidth_manual(
    values = c("Noisy" = 0.3, "symmetric" = 0.8, "one_sided" = 0.8)
  ) +
  scale_colour_manual(
    values = c(
      "Noisy" = "grey75",
      "symmetric" = "#0072B2",
      "one_sided" = "#CC79A7"
    )
  ) +
  labs(
    title = "Causal stream: symmetric vs one_sided (Haar, window = 63)",
    subtitle = "Trending signal, differences emerge at window boundaries",
    x = "Time", y = "Amplitude"
  ) +
  theme_minimal()

## ----irr_custom---------------------------------------------------------------
# 4-point cubic interpolating predict (sum = 1 -> degree = 3)
p_cubic = lift_step(
  "predict",
  coeffs = c(-1/16, 9/16, 9/16, -1/16),
  start_idx = -1
)

u_cubic = lift_step(
  "update",
  coeffs = c(-1/32, 9/32, 9/32, -1/32),
  start_idx = -1
)

my_dd4 = custom_wavelet(
  "MyDD4", list(p_cubic, u_cubic), c(sqrt(2), 1/sqrt(2))
)

# Confirm the degree was inferred correctly (3 means cubic Lagrange)
my_dd4$steps[[1]]$degree

## ----irr_denoise, fig.cap="Figure 4: Custom 4-point cubic wavelet applied to a smooth signal on an irregular grid. The noisy input (grey dots) is sampled at non-uniform time positions drawn from a log-normal-ish process; the underlying signal (black dashed) is recovered by `denoise_signal_offline` with the custom wavelet and `t = t_irr` (orange). The predict step uses Lagrange interpolation evaluated at each sample's physical position rather than fixed coefficients."----
set.seed(20260523)
n_irr = 256

# Irregular grid: cumulative spacing with high variance
t_irr = cumsum(c(0, abs(rnorm(n_irr - 1, mean = 1, sd = 0.6))))

# Smooth underlying signal
pure_irr = sin(t_irr / 8) + 0.3 * cos(t_irr / 3)
noisy_irr = pure_irr + rnorm(n_irr, sd = 0.3)

# Denoise using the custom wavelet on the irregular grid
clean_irr = denoise_signal_offline(
  noisy_irr, my_dd4,
  levels = 4, t = t_irr,
  shrinkage = "semisoft", alpha = 0.3, beta = 1.2
)

df_irr = data.frame(
  t = rep(t_irr, 3),
  value = c(noisy_irr, pure_irr, clean_irr),
  Signal = factor(
    rep(c("Noisy", "Truth", "Denoised (MyDD4)"), each = n_irr),
    levels = c("Noisy", "Truth", "Denoised (MyDD4)")
  )
)

ggplot(df_irr, aes(x = t, y = value, colour = Signal)) +
  geom_point(
    data = subset(df_irr, Signal == "Noisy"),
    size = 0.6, alpha = 0.5
  ) +
  geom_line(
    data = subset(df_irr, Signal == "Truth"),
    linewidth = 0.5, linetype = "dashed"
  ) +
  geom_line(
    data = subset(df_irr, Signal == "Denoised (MyDD4)"),
    linewidth = 0.8
  ) +
  scale_colour_manual(values = c(
    "Noisy" = "grey60",
    "Truth" = "black",
    "Denoised (MyDD4)" = "#D55E00"
  )) +
  theme_minimal() +
  labs(
    title = "Custom wavelet on an irregular grid",
    subtitle = sprintf(
      "n = %d samples; spacing CV = %.0f%%",
      n_irr, 100 * sd(diff(t_irr)) / mean(diff(t_irr))
    ),
    x = "Physical time t", y = "Amplitude"
  )

## ----basis, fig.height=3, fig.cap="Figure 5: Scaling and wavelet basis functions of CDF 9/7 obtained by cascade iteration of the lifting steps. Useful for verifying that a custom wavelet produces reasonable basis functions and for diagnosing coefficient specification errors."----
plot(lifting_scheme("cdf97"))