## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  cache = TRUE,
  collapse = TRUE,
  comment = "#>"
)
library(tf)
oldpar <- par(no.readonly = TRUE)
par(las = 1)
alpha_palette <- c(
  "#00000080", "#DF536B80", "#61D04F80", "#2297E680", "#28E2E580",
  "#CD0BBC80", "#F5C71080", "#9E9E9E80", "#00000080", "#DF536B80"
)
have_srvf <- rlang::is_installed("fdasrvf")
have_cc <- TRUE

## ----eval=FALSE, echo = TRUE, class.source = "fold-show"----------------------
# # One-shot registration (returns tf_registration object):
# reg <- tf_register(x, method = "...")
# tf_aligned(reg)   # registered/aligned curves
# tf_inv_warps(reg)     # estimated inverse warping functions (observed → aligned time)
# tf_template(reg)  # template used
# summary(reg)      # alignment diagnostics
# plot(reg)          # 3-panel diagnostic plot
# 
# # Or step by step:
# warps <- tf_estimate_warps(x, method = "...")
# x_registered <- tf_align(x, warps)

## ----warp-illustration, echo = FALSE, fig.width = 6, fig.height = 6-----------
# Synthetic illustration of different warp types
s <- seq(0, 1, length.out = 201)

# Template: 
template_synth <- tfd(exp(-((s - 0.45)^2) / (2 * 0.08^2)) + cos(2*pi*s) + s, arg = s)

# Warp 1: affine shift + scale (not domain-preserving — extends beyond [0,1])
h_shift <- tfd(0.8 * s + 0.1, arg = s)

# Warp 2: smooth diffeomorphism (domain-preserving)
h_smooth_raw <- s - 0.25 * sin(pi * s)
h_smooth <- tfd(
  (h_smooth_raw - h_smooth_raw[1]) / (h_smooth_raw[length(s)] - h_smooth_raw[1]),
  arg = s)

# Warp 3: wiggly diffeomorphism (domain-preserving)
h_wiggly_raw <- .1 * pbeta(s, .1, .1) + .5 * pbeta(s, 25, 25) + 
  .2 * pbeta(s, 4, 20) + .2 * pbeta(s, 100, 10)
h_wiggly <- tfd(h_wiggly_raw, arg = s)  

# Compose template with warps: m(h(s))
warped_shift <- tf_warp(template_synth, h_shift) |> suppressWarnings()
warped_smooth <- tf_warp(template_synth, h_smooth)
warped_wiggly <- tf_warp(template_synth, h_wiggly)

warp_cols <- c("#DF536B", "#2297E6", "#61D04F")

layout(matrix(1:9, 3, 3))
par(mar = c(4, 4, 3, 1))

warp_names <- c("Affine (shift & scale)", "Smooth elastic", "Wiggly elastic")
warps <- list(h_shift, h_smooth, h_wiggly)
warped <- list(warped_shift, warped_smooth, warped_wiggly)

# Column 1: blank, template, blank
plot.new()
plot(template_synth, lwd = 2.5,
     main = "Template m(s)", xlab = "s", ylab = "m(s)")
plot.new()

# Column 2: one warp per row
for (i in 1:3) {
  plot(s, s, type = "l", lty = 3, lwd = 1,
       main = paste0("h(s): ", warp_names[[i]]),
       xlab = "s", ylab = "h(s)", ylim = c(-0.05, 1.15))
  do.call(graphics::lines, list(warps[[i]], col = warp_cols[[i]], lwd = 2.5))
}

# Column 3: template (dashed) + warped curve per row
for (i in 1:3) {
  plot(template_synth, lwd = 2, lty = 3,
       main = paste0("x(t) = m(h(s)): ", warp_names[[i]]),
       xlab = "t", ylab = "x(t)")
  do.call(graphics::lines, list(warped[[i]], col = warp_cols[[i]], lwd = 2.5))
}

## ----quickstart, class.source = "fold-show"-----------------------------------
pinch_small <- pinch[1:10]
template_affine <- pinch[7] |> tf_smooth(f= .2)

reg_shift <- tf_register(
  pinch_small,
  method = "affine",
  template = template_affine,
  type = "shift"
)
reg_shift
summary(reg_shift)
pinch_inv_warp_shift <- tf_inv_warps(reg_shift)
pinch_small_reg <- tf_aligned(reg_shift)
pinch_template_shift <- tf_template(reg_shift)

## ----quickstart-plot, echo = FALSE, class.source = "fold-show", fig.width = 6, fig.height = 3----
# plot(reg_shift) gives a compact diagnostic; below reuses the extracted
# components for a customized version of the same view:
layout(t(1:3))
plot(
  pinch_small,
  main = "Original Curves",
  xlab = "Observed Time t [sec]",
  ylab = "Force [N]",
  lwd = 1.5,
  col = alpha_palette
)
lines(pinch_template_shift, lwd = 2, lty = 3, col = rgb(0,0,.5))
plot(
  pinch_inv_warp_shift,
  points = FALSE,
  main = "Inverse Warping Functions",
  xlab = "Observed Time t [sec]",
  ylab = "Registered Time s [sec]",
  lwd = 1.5,
  col = alpha_palette
)
abline(0, 1, lty = 3)
plot(
  pinch_small_reg,
  main = "Registered Curves",
  xlab = "Registered Time s [sec]",
  ylab = "Force [N]",
  lwd = 1.5,
  col = alpha_palette, 
  points = FALSE
)
lines(pinch_template_shift, lwd = 2, lty = 3, col = rgb(0,0,.5))

## ----template-choice, fig.width = 6, fig.height = 4---------------------------
# Gaussian bumps with large shifts:
s <- seq(-4, 6, length.out = 201)
mus <- c(-2, -1, 0, 1, 2)
bumps <- tfd(t(sapply(mus, \(mu) dnorm(s, mu, sd = 0.5))), arg = s)

# Pointwise mean is smeared — not a good template:
bumps_mean <- mean(bumps)
# MBD median picks the most central observed curve:
bumps_median <- median(bumps, depth = "MBD")

# Register with default (mean) vs MBD median template:
reg_mean <- tf_register(bumps, method = "affine", type = "shift",
                        template = bumps_mean)
reg_median <- tf_register(bumps, method = "affine", type = "shift",
                          template = bumps_median)

layout(matrix(1:6, 2, 3, byrow = TRUE))
plot(bumps, col = alpha_palette, lwd = 1.5,
     main = "Original + mean template")
lines(bumps_mean, lwd = 3, lty = 2)
plot(tf_inv_warps(reg_mean), col = alpha_palette, lwd = 1.5, points = FALSE,
     main = "Inverse warps (mean template)",
     xlab = "Observed time t", ylab = "Aligned time s")
abline(0, 1, lty = 3)
plot(tf_aligned(reg_mean), col = alpha_palette, lwd = 1.5, points = FALSE,
     main = "Aligned to mean")
lines(bumps_mean, lwd = 3, lty = 2)


plot(bumps, col = alpha_palette, lwd = 1.5,
     main = "Original + MBD median template")
lines(bumps_median, lwd = 3, lty = 2, col = "red3")
plot(tf_inv_warps(reg_median), col = alpha_palette, lwd = 1.5, points = FALSE,
     main = "Inverse warps (median template)",
     xlab = "Observed time t", ylab = "Aligned time s")
abline(0, 1, lty = 3)
plot(tf_aligned(reg_median), col = alpha_palette, lwd = 1.5, points = FALSE,
     main = "Aligned to MBD median")
lines(bumps_median, lwd = 3, lty = 2, col = "red3")

## ----diagnostics, class.source = "fold-show"----------------------------------
x <- pinch[1:10]

# Register with affine shift warps:
reg_aff <- tf_register(x, method = "affine", type = "shift_scale")

# summary() gives a quick quantitative overview:
summary(reg_aff)

## ----diagnostics-plot, fig.width = 6, fig.height = 3, class.source = "fold-show"----
# plot() provides a 3-panel diagnostic:
plot(reg_aff)

## ----diagnostics-features, class.source = "fold-show"-------------------------
# Are global peak locations more concentrated after alignment?
peak_before <- tf_where(x, value == max(value)) |> as.numeric()
peak_after_aff <- tf_where(tf_aligned(reg_aff), value == max(value)) |> as.numeric()

data.frame(
  metric = c("sd_peak_before", "sd_peak_after_affine"),
  value = c(sd(peak_before), sd(peak_after_aff))
)

## ----diagnostics-optional-srvf, eval = have_srvf, include = have_srvf, class.source = "fold-show"----
# Compare with SRVF (non-linear warps):
reg_srvf <- tf_register(x, method = "srvf")

# Side-by-side summaries:
summary(reg_srvf)

## ----diagnostics-optional-srvf-compare, eval = have_srvf, include = have_srvf, class.source = "fold-show"----
peak_after_srvf <- tf_where(tf_aligned(reg_srvf), value == max(value)) |> as.numeric()

data.frame(
  metric = c("sd_peak_before", "sd_peak_after_affine", "sd_peak_after_srvf"),
  value = c(sd(peak_before), sd(peak_after_aff), sd(peak_after_srvf))
)

## ----diagnostics-optional-srvf-plot, eval = have_srvf, include = have_srvf, fig.width = 6, fig.height = 4----
layout(matrix(1:6, 2, 3, byrow = TRUE))
# Affine:
plot(x, main = "Original", col = alpha_palette, lwd = 1.5)
lines(tf_template(reg_aff), lwd = 2, lty = 2)
plot(tf_inv_warps(reg_aff), main = "Affine Inverse Warps",
     col = alpha_palette, lwd = 1.5,  points = FALSE,
     xlab = "Observed time t", ylab = "Aligned time s")
abline(0, 1, lty = 3)
plot(tf_aligned(reg_aff), main = "Affine Aligned",
     col = alpha_palette, lwd = 1.5, points = FALSE)
lines(tf_template(reg_aff), lwd = 2, lty = 2)
# SRVF:
plot(x, main = "Original", col = alpha_palette, lwd = 1.5)
lines(tf_template(reg_srvf), lwd = 2, lty = 2)
plot(tf_inv_warps(reg_srvf), main = "SRVF Inverse Warps",
     col = alpha_palette, lwd = 1.5,   points = FALSE,
     xlab = "Observed time t", ylab = "Aligned time s")
abline(0, 1, lty = 3)
plot(tf_aligned(reg_srvf), main = "SRVF Aligned",
     col = alpha_palette, lwd = 1.5, points = FALSE)
lines(tf_template(reg_srvf), lwd = 2, lty = 2)

## ----pinch-compare, class.source = "fold-show"--------------------------------
pinch_small <- pinch[1:10]

reg_aff <- tf_register(pinch_small, method = "affine", type = "shift_scale")
inv_warp_aff <- tf_inv_warps(reg_aff)

# tf_landmarks_extrema needs smoothed inputs,
#     otherwise it tends to detect lots of spurious features:
pinch_small |> tf_landmarks_extrema(which = "max") |> head()
#     ... so in this case, we simply use the global maximum for each curve:
(peak_locs <- pinch_small |> tf_where(value == max(value)) |> unlist() |> as.matrix())
reg_lm <- tf_register(pinch_small, method = "landmark", landmarks = peak_locs)
inv_warp_lm <- tf_inv_warps(reg_lm)

# ... but using more than one peak location as the only landmark will get us
# better alignment here - use times where curve first & last exceeds 3 as well:
pinch_landmarks  <- cbind(
  start = pinch_small |> tf_where(value > 3, return = "first") |> unlist(),
  peak  = peak_locs,
  end   = pinch_small |> tf_where(value > 3, return = "last") |> unlist()
)
pinch_landmarks |> head()
reg_lm2 <- tf_register(pinch_small, method = "landmark", landmarks = pinch_landmarks)
inv_warp_lm2 <- tf_inv_warps(reg_lm2)

## ----pinch-compare-plot, fig.width = 6, fig.height = 6------------------------
layout(t(matrix(1:12, 4, 3)))
par(cex.main = 0.8)
plot.new()
plot(inv_warp_aff, main = "Affine Inverse Warps", ylab = "",
     points = FALSE, col = alpha_palette, lwd = 1.5)
abline(0, 1, lty = 3)
plot(tf_aligned(reg_aff), main = "Affine Registered",
     col = alpha_palette, lwd = 1.5, points = FALSE)
plot(tf_aligned(reg_aff), main = "Affine Registered",
     type = "lasagna", col = hcl.colors(12, rev = TRUE))

plot(pinch_small, main = "Original",
     col = alpha_palette, lwd = 1.5)
plot(inv_warp_lm, main = "Landmark Inverse Warps \n (Peak only)", ylab = "",
     points = FALSE, col = alpha_palette, lwd = 1.5)
abline(0, 1, lty = 3)
plot(tf_aligned(reg_lm), main = "Landmark Registered\n (Peak only)", col = alpha_palette, lwd = 1.5)
plot(tf_aligned(reg_lm), main = "Landmark Registered\n (Peak only)", type = "lasagna", col = hcl.colors(12, rev = TRUE))

plot.new()
plot(inv_warp_lm2, main = "Landmark Inverse Warps \n (Start + Peak + End)", ylab = "",
     points = FALSE, col = alpha_palette, lwd = 1.5)
abline(0, 1, lty = 3)
plot(tf_aligned(reg_lm2), main = "Landmark Registered\n (Start + Peak + End)", col = alpha_palette, lwd = 1.5)
plot(tf_aligned(reg_lm2), main = "Landmark Registered\n (Start + Peak + End)", type = "lasagna", col = hcl.colors(12, rev = TRUE))

## ----pinch-compare-srvf-cc, eval = have_srvf && have_cc, class.source = "fold-show"----
reg_srvf <- tf_register(pinch_small, method = "srvf")
inv_warp_srvf <- tf_inv_warps(reg_srvf)
reg_cc_unpen <- tf_register(pinch_small, method = "cc", max_iter = 10, nbasis = 10, crit = 1)
inv_warp_cc_unpen <- tf_inv_warps(reg_cc_unpen)
reg_cc_pen <- tf_register(pinch_small, method = "cc", lambda = 1e-4, max_iter = 20)
inv_warp_cc_pen <- tf_inv_warps(reg_cc_pen)

## ----pinch-compare-srvf-cc-summary, eval = have_srvf && have_cc, class.source = "fold-show"----
summary(reg_cc_unpen)
# ... ouch! max slope almost 100 and only 40% amplitude variance reduction ....

summary(reg_srvf)
summary(reg_cc_pen)

## ----pinch-compare-srvf-cc-plot, eval = have_srvf && have_cc, fig.width = 6, fig.height = 6----
layout(t(matrix(1:12, 4, 3)))
par(cex.main = 0.8)
plot.new()
plot(inv_warp_srvf, main = "SRVF Inverse Warps", ylab = "",
     points = FALSE, col = alpha_palette, lwd = 1.5)
abline(0, 1, lty = 3)
plot(tf_aligned(reg_srvf), main = "SRVF Registered",
     col = alpha_palette, lwd = 1.5, points = FALSE)
plot(tf_aligned(reg_srvf), main = "SRVF Registered",
     type = "lasagna", col = hcl.colors(12, rev = TRUE))

plot(pinch_small, main = "Original",
     col = alpha_palette, lwd = 1.5)
plot(inv_warp_cc_unpen, main = "CC Inverse Warps (unpen.)", ylab = "",
     points = FALSE, col = alpha_palette, lwd = 1.5)
abline(0, 1, lty = 3)
plot(tf_aligned(reg_cc_unpen), main = "CC Registered (unpen.)",
     col = alpha_palette, lwd = 1.5, points = FALSE)
plot(tf_aligned(reg_cc_unpen), main = "CC Registered (unpen.)",
     type = "lasagna", col = hcl.colors(12, rev = TRUE))

plot.new()
plot(inv_warp_cc_pen, main = expression("CC Inverse Warps (" * lambda * " = 1e-4)"),
     ylab = "", points = FALSE, col = alpha_palette, lwd = 1.5)
abline(0, 1, lty = 3)
plot(tf_aligned(reg_cc_pen), main = expression("CC Registered (" * lambda * " = 1e-4)"),
     col = alpha_palette, lwd = 1.5, points = FALSE)
plot(tf_aligned(reg_cc_pen), main = expression("CC Registered (" * lambda * " = 1e-4)"),
     type = "lasagna", col = hcl.colors(12, rev = TRUE))

## ----growth-prep, class.source = "fold-show"----------------------------------
growth <- tf::growth |> dplyr::filter(gender == "female")

# Raw velocity via finite differences — noisy, only 30 midpoints from 31 measurements:
growth$raw_vel <- tf_derive(growth$height)

# Smooth velocity via spline representation on a much denser grid, then derive analytically:
growth$smooth_vel <- tfb(
  growth$height,
  k = 15,
  bs = "tp",
  arg = seq(1, 18, l = 80),
  global = TRUE
  # family = gaussian(link = "log") # ensures positivity of velocity
) |>
  tf_derive()

## ----growth-raw-vs-smooth, eval = have_srvf, include = have_srvf, class.source = "fold-show"----
# SRVF on raw (noisy) velocity:
reg_raw_obj <- tf_register(growth$raw_vel, method = "srvf")
inv_warp_raw <- tf_inv_warps(reg_raw_obj)
reg_raw <- tf_aligned(reg_raw_obj)

# SRVF on smooth velocity:
reg_smooth_obj <- tf_register(growth$smooth_vel, method = "srvf")
inv_warp_smooth <- tf_inv_warps(reg_smooth_obj)
reg_smooth <- tf_aligned(reg_smooth_obj)

## ----growth-raw-vs-smooth-plot, eval = have_srvf, include = have_srvf, fig.width = 6, fig.height = 4----
reg_brks <- range(c(tf_evaluations(reg_raw), tf_evaluations(reg_smooth))) |>
  (\(x) seq(x[1], x[2], l = 13))()

layout(t(matrix(1:8, 4, 2)))
par(cex.main = 0.8)

plot(growth$raw_vel, main = "Raw Velocity (finite diff.)",
     xlab = "Age [years]", ylab = "cm/year", lwd = 0.8, ylim = c(0, 30))
plot(inv_warp_raw, points = FALSE, main = "SRVF Inverse Warps (raw)",
     xlab = "Chronological Age", ylab = "Registered Age", lwd = 0.8)
abline(0, 1, lty = 3)
plot(reg_raw, main = "SRVF Registered (raw)",
     xlab = "Registered Age [years]", ylab = "cm/year", lwd = 0.8, ylim = c(0, 30))
plot(reg_raw, main = "SRVF Registered (raw)",
     type = "lasagna", col = hcl.colors(12, rev = TRUE), breaks = reg_brks)

plot(growth$smooth_vel, main = "Smooth Velocity (tfb deriv.)",
     xlab = "Age [years]", ylab = "cm/year", lwd = 0.8, ylim = c(0, 30))
plot(inv_warp_smooth, points = FALSE, main = "SRVF Inverse Warps (smooth)",
     xlab = "Chronological Age", ylab = "Registered Age", lwd = 0.8)
abline(0, 1, lty = 3)
plot(reg_smooth, main = "SRVF Registered (smooth)",
     xlab = "Registered Age [years]", ylab = "cm/year", lwd = 0.8, ylim = c(0, 30))
plot(reg_smooth, main = "SRVF Registered (smooth)",
     type = "lasagna", col = hcl.colors(12, rev = TRUE), breaks = reg_brks)

## ----growth-penalization, eval = have_srvf, include = have_srvf, class.source = "fold-show"----
# Penalized SRVF on the raw velocity:
reg_raw_pen_obj <- tf_register(growth$raw_vel, method = "srvf", lambda = 0.1)
inv_warp_raw_pen <- tf_inv_warps(reg_raw_pen_obj)
reg_raw_pen <- tf_aligned(reg_raw_pen_obj)

## ----growth-penalization-plot, eval = have_srvf, include = have_srvf, fig.width = 6, fig.height = 3----
layout(t(matrix(1:4, 4, 1)))
par(cex.main = 0.8)
plot(growth$raw_vel, main = "Raw Velocity",
     xlab = "Age [years]", ylab = "cm/year", lwd = 0.8)
plot(inv_warp_raw_pen, points = FALSE,
     main = expression("Penalized Warps (raw, " * lambda * " = 0.1)"),
     xlab = "Chronological Age", ylab = "Registered Age", lwd = 0.8)
abline(0, 1, lty = 3)
plot(reg_raw_pen,
     main = expression("Penalized Registered (raw, " * lambda * " = 0.1)"),
     xlab = "Registered Age [years]", ylab = "cm/year", lwd = 0.8)
plot(reg_raw_pen, main = "Penalized Registered (raw)",
     type = "lasagna", col = hcl.colors(12, rev = TRUE), breaks = reg_brks)

## ----growth-landmarks, class.source = "fold-show"-----------------------------
# End of rapid infant growth: 
#    less than 2/3 of max early childhood growth (1-3) velocity before age 5
growth_slows <- growth$smooth_vel |> tf_zoom(begin = 1, end = 5) |> 
  tf_where(value < 0.66 * max(value[1:10]), return = "first") |> unlist()

# Pubertal peak: maximum velocity in age 8--17
growth_peaks <- growth$smooth_vel |> tf_zoom(begin = 8, end = 17) |> 
  tf_where(value == max(value)) |> unlist()

# Pre-pubertal trough: minimum velocity between age 5 and 1 year before the peak
growth_troughs <- growth$smooth_vel |> tf_zoom(begin = 5, end = growth_peaks - 1) |> 
  tf_where(value == min(value)) |> unlist() 

# Build landmark matrix and register:
growth_lm <- cbind(slowdown = growth_slows, trough = growth_troughs, peak = growth_peaks)
reg_lm_obj <- tf_register(growth$smooth_vel, method = "landmark", landmarks = growth_lm)
inv_warp_lm <- tf_inv_warps(reg_lm_obj)
reg_lm <- tf_aligned(reg_lm_obj)

## ----growth-landmarks-plot, fig.width = 6, fig.height = 3---------------------
if (!exists("reg_brks")) {
  reg_brks <- range(c(tf_evaluations(growth$smooth_vel), tf_evaluations(reg_lm))) |>
    (\(x) seq(x[1], x[2], l = 13))()
}
layout(t(1:4))
plot(growth$smooth_vel, main = "Smooth Velocity",
     xlab = "Age [years]", ylab = "cm/year", lwd = 0.8)
plot(inv_warp_lm, points = FALSE, main = "Landmark Inverse Warps",
     xlab = "Chronological Age", ylab = "Registered Age", lwd = 0.8)
abline(0, 1, lty = 3)
plot(reg_lm, main = "Landmark Registered",
     xlab = "Registered Age [years]", ylab = "cm/year", lwd = 0.8)
plot(reg_lm, main = "Landmark Registered",
     type = "lasagna", col = hcl.colors(12, rev = TRUE), breaks = reg_brks)

## ----restore-par, include = FALSE---------------------------------------------
par(oldpar)