## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 7,
  fig.height = 4.2,
  dpi  = 96,
  out.width = "100%"
)
set.seed(1)
# the elastic registration methods (tf_register(method = "srvf_mv"),
# tf_register_shape()) require the suggested fdasrvf package; the chunks
# using them are only evaluated if it is installed.
has_fdasrvf <- requireNamespace("fdasrvf", quietly = TRUE)

## ----packages, message = FALSE------------------------------------------------
library(tf)
library(tidyfun)
library(dplyr)

## ----tour-list----------------------------------------------------------------
data(gait)
g <- tfd_mv(list(hip = gait$hip_angle, knee = gait$knee_angle))
g

## ----gait-anatomy-------------------------------------------------------------
tf_ncomp(g)
names(tf_components(g))
dim(g[1:3, c(0.1, 0.5, 0.9)])
class(g[, , component = "hip"])

## ----tour-array---------------------------------------------------------------
arr <- as.matrix(g)            # [curve, arg, component]
dim(arr)
g_arr <- tfd_mv(arr, arg = tf_arg(g))
identical(as.matrix(g_arr), arr)

## ----tour-df------------------------------------------------------------------
df <- as.data.frame(g, unnest = TRUE, long = FALSE)   # id, arg, hip, knee
head(df, 2)
g_df <- tfd_mv(df, id = "id", arg = "arg", value = c("hip", "knee"))

## ----tour-basis---------------------------------------------------------------
g_spline <- tfb_mv(g, basis = "spline", k = list(hip = 5, knee = 12),
                   verbose = FALSE)
g_spline

## ----tour-geom, fig.height = 3.6----------------------------------------------
head(tf_arclength(g))
plot(tf_speed(g), alpha = 0.3, main = expression("pointwise speed " * "||" * f*"'"*(t) * "||"))

## ----gait-facet, fig.height = 3.6---------------------------------------------
plot(g, type = "facet", alpha = 0.4)

## ----gait-traj, fig.height = 4.2----------------------------------------------
plot(g, alpha = 0.4)   # type = "trajectory" by default for d = 2

## ----gait-mean-sd, fig.height = 3.8-------------------------------------------
mu <- mean(g)
s  <- sd(g)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
for (k in seq_len(tf_ncomp(g))) {
  nm <- names(tf_components(g))[k]
  plot(tf_component(g, k), alpha = 0.25,
       ylab = nm, main = nm)
  lines(tf_component(mu, k),              col = "red", lwd = 2)
  lines(tf_component(mu + 2 * s, k),      col = "red", lwd = 1.2, lty = 2)
  lines(tf_component(mu - 2 * s, k),      col = "red", lwd = 1.2, lty = 2)
}
par(op)

## ----gait-arc, fig.height = 4.2-----------------------------------------------
arc <- tf_arclength(g)
extreme <- c(which.min(arc), which.max(arc))
extreme

plot(g, alpha = 0.15)
lines(g[extreme[1]], col = "steelblue", lwd = 2.2)
lines(g[extreme[2]], col = "firebrick", lwd = 2.2)
legend("bottomright", bty = "n",
       legend = c(sprintf("min arc length (%.0f deg)", arc[extreme[1]]),
                  sprintf("max arc length (%.0f deg)", arc[extreme[2]])),
       col = c("steelblue", "firebrick"), lwd = 2)

## ----gait-rung1, fig.width = 8.2, fig.height = 3.6----------------------------
g_unit <- tf_reparam_arclength(g)
plot(g_unit, type = "facet", alpha = 0.3)

## ----gait-rung1-speed, fig.height = 4.0---------------------------------------
sp_raw  <- tf_speed(g)
sp_unit <- tf_speed(g_unit)
yl <- c(0, max(c(unlist(tf_evaluations(sp_raw)),
                 unlist(tf_evaluations(sp_unit))), na.rm = TRUE))

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(sp_raw,  alpha = 0.4, ylim = yl, main = "raw parameterization",
     ylab = "speed (deg / cycle)")
lines(mean(sp_raw),  col = "firebrick", lwd = 2)
plot(sp_unit, alpha = 0.4, ylim = yl, main = "arc-length parameterization",
     ylab = "speed (deg / cycle)")
lines(mean(sp_unit), col = "firebrick", lwd = 2)
par(op)

## ----gait-rung2, message = FALSE, warning = FALSE, fig.width = 8.2, fig.height = 3.6----
r_knee <- tf_register(g, method = "cc", ref_component = "knee")
plot(tf_aligned(r_knee), type = "facet", alpha = 0.3)

## ----gait-rung2-warps, fig.height = 4.0---------------------------------------
plot(tf_invert(tf_inv_warps(r_knee)), alpha = 0.5,
     main = "reference (knee) warps", ylab = "aligned phase")
abline(0, 1, lty = 3)

## ----gait-rung3, eval = has_fdasrvf, message = FALSE, warning = FALSE, fig.width = 8.2, fig.height = 3.6----
reg_mv <- tf_register(g, method = "srvf_mv", max_iter = 2)
plot(tf_aligned(reg_mv), type = "facet", alpha = 0.3)

## ----gait-rung3-warps, eval = has_fdasrvf, fig.height = 4.0-------------------
plot(tf_invert(tf_inv_warps(reg_mv)), alpha = 0.5,
     main = "srvf_mv warps", ylab = "aligned phase")
abline(0, 1, lty = 3)

## ----gait-rung4, eval = has_fdasrvf, message = FALSE, warning = FALSE, fig.width = 8.2, fig.height = 3.6----
reg_shape <- tf_register_shape(g, max_iter = 2)
plot(tf_aligned(reg_shape), type = "facet", alpha = 0.3)   # note the y-axis: shape space

## ----gait-rung4-traj, eval = has_fdasrvf, message = FALSE, warning = FALSE, fig.width = 8.2, fig.height = 4.2----
op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(g,                     type = "trajectory", alpha = 0.3, main = "original")
plot(tf_aligned(reg_shape), type = "trajectory", alpha = 0.3,
     main = "shape-registered (shape space)")
par(op)

## ----gait-rung4-diag, eval = has_fdasrvf, message = FALSE, warning = FALSE----
round(tf_rotations(reg_shape)[, , 1], 3)   # rotation for subject 1
summary(as.numeric(tf_scales(reg_shape)))  # scales, relative to the template

## ----quant-warp-helper, class.source = "fold-hide"----------------------------
# warp implied by arc-length reparametrization: normalised cumulative arc length
# (raw cycle phase -> arc-length phase, same direction as the registration warps)
arclen_warp <- function(mv) {
  sp  <- tf_speed(mv)
  arg <- as.numeric(tf_arg(sp))
  dom <- tf_domain(sp)
  d   <- diff(arg)
  w <- sapply(tf_evaluations(sp), function(s) {
    cum <- c(0, cumsum((head(s, -1) + tail(s, -1)) / 2 * d))
    dom[1] + diff(dom) * cum / cum[length(cum)]
  })
  tfd(t(w), arg = arg, domain = dom)
}

## ----gait-grid, eval = has_fdasrvf, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 6.6----
cols <- list(
  list(nm = "as-observed", data = g,                     warp = NULL),
  list(nm = "arc-length",  data = g_unit,                warp = arclen_warp(g)),
  list(nm = "ref = knee",  data = tf_aligned(r_knee),    warp = tf_invert(tf_inv_warps(r_knee))),
  list(nm = "srvf_mv",     data = tf_aligned(reg_mv),    warp = tf_invert(tf_inv_warps(reg_mv))),
  list(nm = "shape",       data = tf_aligned(reg_shape), warp = tf_invert(tf_inv_warps(reg_shape)))
)

op <- par(mfrow = c(3, 5), mar = c(2.6, 3.2, 2, 0.8), mgp = c(1.9, 0.6, 0))
for (j in seq_along(cols))
  plot(tf_component(cols[[j]]$data, "hip"),  alpha = 0.3,
       main = cols[[j]]$nm, ylab = if (j == 1) "hip" else "")
for (j in seq_along(cols))
  plot(tf_component(cols[[j]]$data, "knee"), alpha = 0.3,
       main = "", ylab = if (j == 1) "knee" else "")
for (j in seq_along(cols)) {
  w <- cols[[j]]$warp
  if (is.null(w)) {
    plot.new(); text(0.5, 0.5, "(identity)", col = "grey55")
  } else {
    plot(w, alpha = 0.4, main = "", ylab = if (j == 1) "warp" else "")
    abline(0, 1, lty = 3)
  }
}
par(op)

## ----gait-register-table-helper, class.source = "fold-hide"-------------------
peak_sd <- function(f, k) {
  arg <- tf_arg(f)
  max(unlist(tf_evaluate(tf_component(sd(f), k), arg)))
}

## ----gait-register-table, eval = has_fdasrvf----------------------------------
rungs <- list(
  "raw"                = g,
  "1 arc-length"       = g_unit,
  "2 reference (knee)" = tf_aligned(r_knee),
  "3 srvf_mv"          = tf_aligned(reg_mv),
  "4 shape (*)"        = tf_aligned(reg_shape)
)
data.frame(
  registration = names(rungs),
  sd_hip       = sapply(rungs, peak_sd, "hip"),
  sd_knee      = sapply(rungs, peak_sd, "knee"),
  row.names    = NULL
)

## ----shape-demo-build---------------------------------------------------------
t      <- seq(0, 1, length.out = 51)
base   <- rbind(t, t^2)
scales <- c(1, 0.7, 1.3)
angles <- c(0, 0.4, -0.25)
offset <- rbind(c(0.2, -0.1), c(0.4, -0.2), c(0.6, -0.3))
beta   <- array(NA_real_, dim = c(3, length(t), 2),
                dimnames = list(c("a", "b", "c"), NULL, c("x", "y")))
for (i in 1:3) {
  rot   <- matrix(c(cos(angles[i]),  sin(angles[i]),
                    -sin(angles[i]), cos(angles[i])), nrow = 2)
  curve <- scales[i] * (rot %*% base) + matrix(offset[i, ], 2, length(t))
  beta[i, , 1] <- curve[1, ]
  beta[i, , 2] <- curve[2, ]
}
shapes <- tfd_mv(beta, arg = t)

## ----shape-demo-quotients, eval = has_fdasrvf, message = FALSE, warning = FALSE, fig.width = 7.4, fig.height = 6.6----
reg_full  <- tf_register_shape(shapes, max_iter = 2, rotation = TRUE,  scale = TRUE)
reg_rot   <- tf_register_shape(shapes, max_iter = 2, rotation = TRUE,  scale = FALSE)
reg_scale <- tf_register_shape(shapes, max_iter = 2, rotation = FALSE, scale = TRUE)

op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(shapes,                asp = 1, col = 1:3, lwd = 2, main = "input: 3 placements")
plot(tf_aligned(reg_rot),   asp = 1, col = 1:3, lwd = 2, main = "rotation-only (keeps size)")
plot(tf_aligned(reg_scale), asp = 1, col = 1:3, lwd = 2, main = "scale-only (keeps orientation)")
plot(tf_aligned(reg_full),  asp = 1, col = 1:3, lwd = 2, main = "rotation + scale (full)")
par(op)

## ----shape-demo-scales, eval = has_fdasrvf------------------------------------
data.frame(curve      = c("a", "b", "c"),
           injected   = scales,
           full       = round(as.numeric(tf_scales(reg_full)),  3),
           rot_only   = round(as.numeric(tf_scales(reg_rot)),   3),
           scale_only = round(as.numeric(tf_scales(reg_scale)), 3))

## ----gait-fpc-helper, class.source = "fold-hide"------------------------------
fpc_var <- function(comp) attr(comp, "score_variance")

## ----gait-fpc, fig.height = 3.8-----------------------------------------------
g_aligned <- tf_aligned(r_knee)
g_b <- tfb_mv(g_aligned, basis = "fpc", verbose = FALSE)

v_hip  <- fpc_var(tf_component(g_b, "hip"));  v_hip  <- v_hip  / sum(v_hip)
v_knee <- fpc_var(tf_component(g_b, "knee")); v_knee <- v_knee / sum(v_knee)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
barplot(head(v_hip,  6), names.arg = 1:6, main = "hip: variance share",
        ylab = "proportion", col = "grey70")
barplot(head(v_knee, 6), names.arg = 1:6, main = "knee: variance share",
        ylab = "proportion", col = "grey70")
par(op)

# residual: how well does the low-rank FPC basis approximate the curves?
g_round <- vctrs::vec_cast(g_b, g_aligned)
resid   <- g_aligned - g_round
rmse_per_subject <- sqrt(unlist(lapply(tf_evaluations(resid), function(m) {
  mean(as.matrix(m[, -1L, drop = FALSE])^2)
})))
summary(rmse_per_subject)

## ----gait-mfpc----------------------------------------------------------------
g_m <- tfb_mfpc(g_aligned, pve = 0.95)
g_m

scores <- tf_mfpc_scores(g_m)          # 39 x M, shared across both components
nu     <- attr(g_m, "mfpc")$evalues    # joint (multivariate) eigenvalues
ve     <- nu / sum(nu)                  # variance share per shared mode
dim(scores)
round(head(ve, 4), 3)

## ----gait-mfpc-efun, fig.width = 8.2, fig.height = 3.6------------------------
psi    <- tf_mfpc_efunctions(g_m)
k_show <- min(3L, length(psi))
op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(tf_component(psi, "hip")[seq_len(k_show)], col = seq_len(k_show), lwd = 2,
     main = "MFPC eigenfunctions: hip", ylab = "loading")
plot(tf_component(psi, "knee")[seq_len(k_show)], col = seq_len(k_show), lwd = 2,
     main = "MFPC eigenfunctions: knee", ylab = "loading")
legend("topright", bty = "n", lwd = 2, col = seq_len(k_show),
       legend = paste("MFPC", seq_len(k_show)), cex = 0.85)
par(op)

## ----gait-mfpc-rmse, class.source = "fold-hide"-------------------------------
rmse_mv <- function(approx) {
  resid <- g_aligned - vctrs::vec_cast(approx, g_aligned)
  sqrt(mean(unlist(lapply(tf_evaluations(resid), function(m) {
    as.matrix(m[, -1L, drop = FALSE])^2
  }))))
}
n_indep <- sum(vapply(tf_components(g_b),
                      function(co) length(attr(co, "score_variance")), integer(1)))
data.frame(
  representation = c("independent FPC", "joint MFPCA"),
  stored_scores  = c(n_indep, attr(g_m, "mfpc")$npc),
  rmse           = round(c(rmse_mv(g_b), rmse_mv(g_m)), 3)
)

## ----gait-mfpc-recap, echo = FALSE, results = "asis"--------------------------
cat(sprintf(
  "The first %d shared modes already capture %.0f%% of the joint variance. ",
  min(2L, length(ve)), 100 * sum(head(ve, 2))
))
cat(sprintf(
  "New subjects can be projected onto this fitted basis with `tf_rebase()`, which re-scores them jointly rather than component by component."
))

## ----storms-build, class.source = "fold-hide"---------------------------------
KM_PER_DEG <- 111.32

storms_clean <- storms |>
  mutate(
    storm_id = paste(name, year),
    ts       = as.POSIXct(ISOdate(year, month, day, hour), tz = "UTC")
  ) |>
  group_by(storm_id) |>
  distinct(ts, .keep_all = TRUE) |>      # dedupe duplicate-timestamp rows
  mutate(
    t_hours = as.numeric(ts - min(ts), units = "hours"),
    phase   = if (max(t_hours) > 0) t_hours / max(t_hours) else 0,
    ref_lat = mean(lat),
    x_km    = (long - mean(long)) * KM_PER_DEG * cos(ref_lat * pi / 180),
    y_km    = (lat  - ref_lat)    * KM_PER_DEG
  ) |>
  filter(n() >= 16, !is.na(wind), !is.na(pressure)) |>  # >= 4 days
  ungroup()

# spatial (km, real time): the compact data-frame constructor.
tracks_km <- storms_clean |>
  transmute(storm_id, t_hours, x = x_km, y = y_km) |>
  tfd_mv(id = "storm_id", arg = "t_hours", value = c("x", "y"))

# Equivalent spelling if the component tfd vectors are already available.
tracks_km_from_components <- tfd_mv(list(
  x = tfd(storms_clean, id = "storm_id", arg = "t_hours", value = "x_km"),
  y = tfd(storms_clean, id = "storm_id", arg = "t_hours", value = "y_km")
))

# full 4-d, on normalised life-cycle phase
tracks4 <- tfd_mv(list(
  long = tfd(storms_clean, id = "storm_id", arg = "phase", value = "long"),
  lat  = tfd(storms_clean, id = "storm_id", arg = "phase", value = "lat"),
  wind = tfd(storms_clean, id = "storm_id", arg = "phase", value = "wind"),
  pres = tfd(storms_clean, id = "storm_id", arg = "phase", value = "pressure")
))
tracks4

# storm-level metadata
peak <- storms |>
  group_by(name, year) |>
  summarise(peak_cat = suppressWarnings(max(category, na.rm = TRUE)),
            .groups = "drop") |>
  mutate(peak_cat = if_else(is.finite(peak_cat), peak_cat, 0),
         peak_cat = as.integer(peak_cat),
         storm_id = paste(name, year))

df <- tibble::tibble(
  storm_id  = names(tracks4),
  track     = tracks4,
  track_km  = tracks_km
) |>
  left_join(peak, by = "storm_id") |>
  mutate(strength = factor(
    pmin(peak_cat, 4),
    levels = 0:4,
    labels = c("TS/TD", "Cat 1", "Cat 2", "Cat 3", "Cat 4+")
  ))

## ----storms-irregular---------------------------------------------------------
head(tf_count(tracks4))
as.data.frame(tracks4[1:2], unnest = TRUE) |> head()

## ----storms-movement-workflow, fig.width = 8.2, fig.height = 6.6--------------
movement_ids <- c("Andrew 1992", "Katrina 2005", "Sandy 2012")
movement_ids <- movement_ids[movement_ids %in% df$storm_id]
movement_rows <- match(movement_ids, df$storm_id)
movement_grid <- seq(0, 1, length.out = 81)
movement_duration <- vapply(
  tf_arg(df$track_km[movement_rows]),
  \(t) max(t) - min(t),
  numeric(1)
)

movement_eval <- lapply(seq_along(movement_rows), function(k) {
  tf_evaluate(
    df$track_km[movement_rows[k]],
    arg = movement_grid * movement_duration[k]
  )[[1]]
})
movement_x <- do.call(rbind, lapply(movement_eval, `[[`, "x"))
movement_y <- do.call(rbind, lapply(movement_eval, `[[`, "y"))

movement_track <- tfd_mv(list(
  x = tfd(movement_x, arg = movement_grid, domain = c(0, 1)),
  y = tfd(movement_y, arg = movement_grid, domain = c(0, 1))
), domain = c(0, 1))
names(movement_track) <- movement_ids

movement_velocity <- tf_derive(movement_track)
movement_speed <- tf_norm(movement_velocity) / movement_duration

vx <- as.matrix(tf_component(movement_velocity, "x"),
                arg = movement_grid, interpolate = TRUE)
vy <- as.matrix(tf_component(movement_velocity, "y"),
                arg = movement_grid, interpolate = TRUE)
heading_mat <- atan2(vy, vx) * 180 / pi
heading <- tfd(heading_mat, arg = movement_grid, domain = c(0, 1))
names(heading) <- movement_ids

wrap_degrees <- function(x) ((x + 180) %% 360) - 180
turning_mat <- t(apply(heading_mat, 1, function(x) {
  c(NA_real_, wrap_degrees(diff(x)))
}))
turning <- tfd(turning_mat, arg = movement_grid, domain = c(0, 1))
names(turning) <- movement_ids

movement_cols <- c("#1b9e77", "#d95f02", "#7570b3")[seq_along(movement_ids)]
xy <- as.matrix(movement_track, interpolate = TRUE)

op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(range(xy[, , "x"], na.rm = TRUE), range(xy[, , "y"], na.rm = TRUE),
     type = "n", asp = 1, xlab = "x (km)", ylab = "y (km)",
     main = "regularized tracks")
lines(movement_track, col = movement_cols, lwd = 2)
legend("topleft", bty = "n", lwd = 2, col = movement_cols,
       legend = movement_ids, cex = 0.85)
plot(movement_speed, col = movement_cols, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "km / h", main = "forward speed")
plot(heading, col = movement_cols, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "degrees", main = "heading")
plot(turning, col = movement_cols, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "degrees", main = "turning angle")
par(op)

## ----storms-katrina, fig.width = 8, fig.height = 5.4--------------------------
katrina <- df |> filter(storm_id == "Katrina 2005") |> pull(track)
plot(katrina, type = "facet", lwd = 2)

## ----storms-plot-helpers, class.source = "fold-hide"--------------------------
have_maps <- requireNamespace("maps", quietly = TRUE)
pal <- c("grey50", "#fed976", "#feb24c", "#fd8d3c", "#e31a1c")

draw_coast <- function(xlim, ylim) {
  if (!have_maps) return(invisible())
  m <- maps::map("world", plot = FALSE,
                 xlim = xlim, ylim = ylim, fill = FALSE)
  lines(m$x, m$y, col = "grey55", lwd = 0.6)
}

## ----storms-map, fig.width = 8.2, fig.height = 6.8----------------------------
xlim <- range(unlist(tf_evaluations(tf_component(df$track, "long"))), na.rm = TRUE)
ylim <- range(unlist(tf_evaluations(tf_component(df$track, "lat"))),  na.rm = TRUE)

op <- par(mfrow = c(2, 3), mar = c(3.4, 3.4, 2, 1), mgp = c(2, 0.7, 0))
for (k in seq_along(levels(df$strength))) {
  lev   <- levels(df$strength)[k]
  trks  <- df$track[df$strength == lev]
  spatial <- tfd_mv(list(
    long = tf_component(trks, "long"),
    lat  = tf_component(trks, "lat")
  ))
  plot(range(xlim), range(ylim), type = "n",
       xlab = "long", ylab = "lat",
       main = sprintf("%s (n = %d)", lev, length(trks)))
  draw_coast(xlim, ylim)
  lines(spatial, col = pal[k], alpha = 0.4, lwd = 1)
}
plot.new()
legend("center", bty = "n", lwd = 2, col = pal, legend = levels(df$strength),
       title = "peak intensity", cex = 1.05)
par(op)

## ----storms-scalars-----------------------------------------------------------
df <- df |> mutate(
  path_km    = tf_arclength(track_km),
  duration   = vapply(tf_arg(track_km), \(t) max(t) - min(t), numeric(1)),
  mean_speed = path_km / duration,            # km/h, lifetime average
  peak_wind  = vapply(tf_evaluations(tf_component(track, "wind")), max, numeric(1)),
  min_pres   = vapply(tf_evaluations(tf_component(track, "pres")), min, numeric(1))
)

df |>
  group_by(strength) |>
  summarise(
    n = dplyr::n(),
    median_path_km     = round(median(path_km)),
    median_speed_kmh   = round(median(mean_speed), 1),
    median_peak_wind   = median(peak_wind),
    median_min_pres    = median(min_pres),
    .groups = "drop"
  )

## ----storms-lifecycle-prep, class.source = "fold-hide"------------------------
phase_grid <- seq(0, 1, length.out = 41)

# forward speed on normalised time: re-arg each storm's km-speed by
# t / T_i, so all storms share phase domain [0, 1].
speed_km  <- tf_speed(df$track_km)
durations <- df$duration
speed_long <- do.call(rbind, lapply(seq_along(speed_km), function(i) {
  data.frame(
    id    = names(speed_km)[i],
    phase = tf_arg(speed_km[i])[[1]] / durations[i],
    value = tf_evaluations(speed_km[i])[[1]]
  )
}))
speed_phase <- tfd(speed_long, id = "id", arg = "phase", value = "value",
                   domain = c(0, 1))

# per-stratum mean curve on the regular phase grid, packaged as length-G tfd
stratum_mean_tfd <- function(comp, grp, grid = phase_grid) {
  mat <- as.matrix(comp, arg = grid, interpolate = TRUE)
  means <- vapply(levels(grp), function(g) {
    rows <- which(grp == g)
    if (length(rows) < 2L) return(rep(NA_real_, length(grid)))
    colMeans(mat[rows, , drop = FALSE], na.rm = TRUE)
  }, numeric(length(grid)))
  out <- tfd(t(means), arg = grid, domain = c(0, 1))
  names(out) <- levels(grp)
  out
}

wind_avg  <- stratum_mean_tfd(tf_component(df$track, "wind"), df$strength)
pres_avg  <- stratum_mean_tfd(tf_component(df$track, "pres"), df$strength)
speed_avg <- stratum_mean_tfd(speed_phase,                     df$strength)

## ----storms-lifecycle, fig.width = 8.2, fig.height = 6.8----------------------
op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(wind_avg, col = pal, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "wind (knots)",
     main = "mean sustained wind")
legend("topright", bty = "n", lwd = 2, col = pal, legend = levels(df$strength),
       cex = 0.9, title = "peak intensity")
plot(pres_avg,  col = pal, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "pressure (mbar)",
     main = "mean central pressure")
plot(speed_avg, col = pal, lwd = 2, alpha = 1,
     xlab = "lifecycle phase", ylab = "forward speed (km/h)",
     main = "mean forward speed")
plot.new()
text(0.5, 0.7, "wind and pressure peak\nnear lifecycle phase 0.5;",
     adj = 0.5, cex = 1.05)
text(0.5, 0.35, "forward speed peaks\nlate, during recurvature",
     adj = 0.5, cex = 1.05)
par(op)

## ----storms-tfb, fig.width = 8.2, fig.height = 4.8, warning = FALSE-----------
top6 <- df |> arrange(desc(path_km)) |> slice(1:6) |> pull(storm_id)

# fit a per-component spline basis on (long, lat, wind, pres) over [0, 1]
tb <- tfb_mv(df$track[df$storm_id %in% top6], k = 12, verbose = FALSE)

# pull just the (long, lat) components out of the 4-d objects so plot.tf_mv
# defaults to the trajectory (long, lat) view
raw_xy <- df$track[df$storm_id %in% top6, , c("long", "lat")]
sm_xy  <- tb[, , c("long", "lat")]

op <- par(mfrow = c(1, 2), mar = c(3.6, 3.6, 2, 1), mgp = c(2, 0.7, 0))
plot(range(xlim), range(ylim), type = "n",
     xlab = "long", ylab = "lat", main = "raw observations")
draw_coast(xlim, ylim)
lines(raw_xy, col = 1:6, lwd = 1.6)
plot(range(xlim), range(ylim), type = "n",
     xlab = "long", ylab = "lat", main = "spline-smoothed")
draw_coast(xlim, ylim)
lines(sm_xy, col = 1:6, lwd = 1.6)
par(op)

## ----storms-tfb-intensity, fig.width = 8.2, fig.height = 3.4------------------
op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(tf_component(tb, "wind"), col = 1:6, lwd = 2,
     xlab = "lifecycle phase", ylab = "wind (knots)",
     main = "smoothed wind")
plot(tf_component(tb, "pres"), col = 1:6, lwd = 2,
     xlab = "lifecycle phase", ylab = "pressure (mbar)",
     main = "smoothed pressure")
par(op)