## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----builtin------------------------------------------------------------------
library(rLifting)

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

haar   = lifting_scheme("haar")
cdf97  = lifting_scheme("cdf97")
print(haar)
print(cdf97)

## ----custom_define------------------------------------------------------------
# Step 1: Predict (P) — coeffs = c(0.5, 0.5), starting at index 0
p_step = lift_step(
  type = "predict",
  coeffs = c(0.5, 0.5),
  start_idx = 0
)

# Step 2: Update (U) — coeffs = c(0.25, 0.25), starting at index -1
u_step = lift_step(
  type = "update",
  coeffs = c(0.25, 0.25),
  start_idx = -1
)

# Combine into a scheme
my_linear = custom_wavelet(
  name = "LinearManual",
  steps = list(p_step, u_step),
  norm = c(sqrt(2), 1/sqrt(2))  # biorthogonal normalization
)

print(my_linear)

## ----custom_use---------------------------------------------------------------
# A linear ramp: detail coefficients should be zero
x = 1:16
fwd = lwt(x, my_linear, levels = floor(log2(length(x))))

# Detail at level 1: all zeros for a linear signal
print(round(fwd$coeffs$d1, 10))

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

## ----diagnose-----------------------------------------------------------------
diag = diagnose_wavelet(
  my_linear,
  config = list(
    is_ortho  = FALSE,   # CDF 5/3 is biorthogonal, not orthogonal
    vm_degrees = c(0, 1), # test vanishing moments for degree 0 and 1
    max_taps  = 5         # expected filter support width
  ),
  plot = FALSE  # suppress basis plot in vignette
)

## ----individual_tests---------------------------------------------------------
# Test perfect reconstruction separately
pr = validate_perfect_reconstruction(my_linear)
cat("Perfect reconstruction:", pr$passed, "\n")
cat("Max error:", format(pr$max_error, scientific = TRUE), "\n")

## ----threshold_demo, fig.width=6, fig.height=4--------------------------------
x = seq(-3, 3, length.out = 201)
lambda = 1.0

df_th = data.frame(
  x = rep(x, 3),
  y = c(threshold(x, lambda, "hard"),
        threshold(x, lambda, "soft"),
        threshold(x, lambda, "semisoft")),
  Method = rep(c("Hard", "Soft", "Semisoft"), each = length(x))
)

library(ggplot2)
ggplot(df_th, aes(x = x, y = y, colour = Method)) +
  geom_line(linewidth = 0.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dotted", colour = "grey60") +
  geom_vline(xintercept = c(-lambda, lambda), linetype = "dashed", colour = "grey40") +
  labs(title = "Thresholding functions",
       subtitle = expression(paste(lambda, " = 1")),
       x = "Input coefficient", y = "Output coefficient") +
  theme_minimal() +
  scale_colour_manual(values = c("Hard" = "#E69F00", "Soft" = "#56B4E9",
                                  "Semisoft" = "#009E73"))

## ----adaptive-----------------------------------------------------------------
set.seed(42)
n = 256
signal = sin(seq(0, 4*pi, length.out = n)) + rnorm(n, sd = 0.3)

sch = lifting_scheme("cdf97")
decomp = lwt(signal, sch, levels = 4)

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

## ----pipeline, fig.width=7, fig.height=4--------------------------------------
set.seed(42)
n = 512
t_vec = seq(0, 1, length.out = n)
pure = sin(2 * pi * 5 * t_vec) + 0.5 * sin(2 * pi * 20 * t_vec)
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 thresholding: apply semisoft to detail levels,
#    but leave the coarsest level (d5) untouched
for (lev in 1:4) {
  dname = paste0("d", lev)
  decomp$coeffs[[dname]] = threshold(
    decomp$coeffs[[dname]], lambdas[[dname]], method = "semisoft"
  )
}
# d5 is not thresholded — preserving low-frequency structure

# 4. Inverse transform
reconstructed = ilwt(decomp)

# 5. Compare
df_pipe = data.frame(
  t = rep(t_vec, 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()

## ----basis, fig.width=6, fig.height=4-----------------------------------------
plot(lifting_scheme("cdf97"))