---
title: "Vector-valued functional data: gait cycles and hurricane tracks"
output:
  rmarkdown::html_vignette:
    code_folding: show
bibliography: references.bib
link-citations: true
vignette: >
  %\VignetteIndexEntry{Vector-valued functional data: gait cycles and hurricane tracks}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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)
```

A `tf` vector represents a sample of functions $f: \mathbb{R} \to \mathbb{R}$
(see the [tf vectors vignette](https://tidyfun.github.io/tidyfun/articles/x01_tf_Vectors.html)
for an introduction to these classes and their operations).
Many real measurement processes, though, produce several coupled signals that
share a common argument: the `(hip, knee)` joint-angle pair sampled across a
gait cycle, the `(longitude, latitude)` position of a hurricane sampled in
time, the `(x, y, z)` body coordinates of a moving animal. The natural object
is a **vector-valued function** $f: \mathbb{R} \to \mathbb{R}^d$, and `tf`
represents a sample of those with the `tf_mv` family: `tfd_mv` for raw
evaluations and `tfb_mv` for basis representations. Internally a `tf_mv` is a
bundle of $d$ ordinary `tf` vectors, so the multivariate methods shown below
reuse the existing univariate numerical kernels by mapping over components.

This article uses two real datasets to put the API through its paces:

* **`tf::gait`** -- regularly sampled hip and knee angles for 39 boys across
  one gait cycle. *Question: how variable is "normal" gait, where in the
  cycle does that variation concentrate, and what are the dominant modes of
  variation?*
* **`dplyr::storms`** -- irregularly sampled position and intensity of every
  Atlantic tropical storm or hurricane from `r min(dplyr::storms$year)` to
  `r max(dplyr::storms$year)`. *Question: do stronger storms travel further
  and move faster than weak ones?*

```{r packages, message = FALSE}
library(tf)
library(tidyfun)
library(dplyr)
```

# 1. A quick tour of vector-valued functional data

Before the case studies, this part introduces the objects: what a vector-valued
function *is*, the two classes that represent samples of them, how to construct
them, and the handful of geometric operations the case studies lean on.

## What is a vector-valued function?

A `tf` vector stores a sample of scalar functions $f: [a, b] \to \mathbb{R}$. A
**vector-valued** function bundles $d$ of them on a *shared* argument,

$$ f(t) = \big(f_1(t), \dots, f_d(t)\big), \qquad t \in [a, b], $$

so that as $t$ runs over the domain the value $f(t)$ traces a *path* through
$\mathbb{R}^d$. `tf` represents a sample of such paths with the `tf_mv` family.
Internally a `tf_mv` *is* a bundle of `d` ordinary `tf` vectors sharing one
argument, so every univariate verb -- evaluation, arithmetic, derivatives,
integration, smoothing, basis fitting -- extends component-wise for free.

Two kinds of variation run through everything below: **amplitude** (how large
the values get) and **phase** (*where in* $t$ the features happen). Much of the
work in functional data analysis [@ramsay2005functional] is telling the two
apart [@marron2015], and Section 2 is largely about exactly that.

## The `tf_mv` classes: `tfd_mv` and `tfb_mv`

Like their univariate cousins `tfd`/`tfb`, vector-valued functions come in two
flavours:

* `tfd_mv` -- *raw evaluations*: each component is a `tfd`, i.e. values on a grid
  plus an interpolation rule.
* `tfb_mv` -- *basis representation*: each component is a `tfb` (spline or FPC),
  i.e. basis coefficients.

Components are named, and a few accessors reach into the bundle: `tf_ncomp()`
(number of components), `tf_components()` (the list of univariate `tf` vectors),
`tf_component(x, k)` (one component, by name or index), and `names()`.

## Constructing `tf_mv` objects

`tfd_mv()` accepts three interchangeable input layouts. We use the built-in
`gait` data [@olshen1989gait; @ramsay2005functional] throughout the gait case
study -- hip and knee angles, in degrees, for 39 boys, each sampled at 20
points across one gait cycle. The most direct
construction is a **named list** of univariate `tf` vectors:

```{r tour-list}
data(gait)
g <- tfd_mv(list(hip = gait$hip_angle, knee = gait$knee_angle))
g
```

The print-out shows `d = 2` components on a shared grid, with a sparkline pair
per subject. It behaves like a vector of 39 bivariate curves; component
accessors and bracket extraction expose the two levels separately -- curves on
rows, argument values on columns, components on the third array dimension:

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

The *same* object can be built from a **3-d array**
`[curve, arg, component]` -- the layout `as.matrix()` returns --

```{r 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)
```

or from a **wide data frame** with one column per component (the long
`(id, arg, component, value)` schema is what `as.data.frame(unnest = TRUE)`
emits; `long = FALSE` gives the wide schema the constructor consumes):

```{r 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"))
```

Each route also takes an optional `domain` and a per-component `evaluator`
(interpolation rule, default `tf_approx_linear`). A `tfd_mv` becomes a basis
object with `tfb_mv()`, choosing a `"spline"` or `"fpc"` basis; basis arguments
apply globally or per component when passed as a named list:

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

## Geometry on the bundle: speed, arc length, reparametrization

Because a `tf_mv` traces a path, the usual curve-geometry quantities are
available, and they return *univariate* summaries of the multivariate object:

* `tf_speed(f)` $= \lVert f'(t)\rVert$ -- pointwise speed, a `tfd`,
* `tf_arclength(f)` -- total path length, one number per curve,
* `tf_reparam_arclength(f)` -- the same path traversed at *unit speed*.

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

## Aligning vector-valued curves: a ladder

When curves differ mainly in *timing*, we usually want to factor that out --
but "factor out" can mean several different things. There is a ladder of
increasingly aggressive operations, distinguished by *what each treats as a
nuisance* [@marron2015]:

| rung | operation | removes | keeps | `tf` entry point |
|---|---|---|---|---|
| 1 | arc-length reparametrization | the clock (parametrization only) | shape **and** size | `tf_reparam_arclength()` |
| 2 | warp from a 1-d reference signal | phase, onto a chosen reference | amplitude, size, orientation | `tf_register(., method = "cc")` |
| 3 | joint multivariate warp | phase, from all components at once | amplitude, size, orientation | `tf_register(., method = "srvf_mv")` |
| 4 | full elastic shape registration | phase **+** rotation **+** scale | shape only | `tf_register_shape()` |

Rung 1 re-labels time within each curve without picking a template; rungs 2 and
3 warp every curve onto a common template (from one reference signal, or from all
components jointly) but leave the values untouched; rung 4 *additionally* rotates
and rescales each curve, landing in a **shape space** where only the geometric
form remains. Section 2 walks all four rungs on the gait data -- and shows where
the top rung is, and is not, the right tool.

# 2. How variable is a "normal" gait cycle?

We built the two-component object `g` above. The two natural views of these 39
curves are the **time-series**
(`type = "facet"`, one panel per component) and the **trajectory in phase
space** (`type = "trajectory"`, the default for two-component objects).

```{r gait-facet, fig.height = 3.6}
plot(g, type = "facet", alpha = 0.4)
```

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

The facet view shows that both joints flex twice per cycle; the trajectory
view collapses time and reveals the characteristic "butterfly" loop traced
out as the leg cycles through stance and swing. Both views suggest the
between-subject spread is far from uniform along the cycle.

## Pointwise mean and standard deviation

`mean()` and `sd()` are vctrs group generics and dispatch component-wise:
they each return a length-1 `tfd_mv`. Plot the raw curves with `plot.tf`,
then overlay the mean +/- 2 sd envelope with `lines.tf` -- all arithmetic
on the components is component-wise on `tf_mv`, so `mu + 2 * s` is itself
a `tfd_mv` and `tf_component(...)` gives the per-axis envelope:

```{r 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)
```

The pointwise sd peaks roughly where the angle itself is changing fastest --
around the stance/swing transitions -- so most of the between-subject
variability lives in *timing* of those transitions, not in the angles
attained at rest.

## Most- and least-varied subjects

`tf_arclength()` measures the path length traced out in `(hip, knee)`-space.
A short orbit = a subject whose `(hip, knee)` excursions are small or who
moves the two joints in lock-step; a long orbit = a subject with
large-amplitude, less synchronised excursions.

```{r 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)
```

## A ladder of registrations

The pointwise sd above conflates two kinds of between-subject variability:
*amplitude* (how far the knee swings) and *phase* (at what fraction of the cycle
heel-strike happens). Removing phase -- "registration" -- can mean progressively
more, depending on what we are willing to treat as a nuisance. Below are four
rungs, applied in turn to the gait sample, each followed by a facet plot of the
result; we quantify them together at the end.

Rungs 3 and 4 (and the rung-comparison at the end) use the elastic
registration routines from the suggested **fdasrvf** package -- those code
chunks are only evaluated if `fdasrvf` is installed.

### Rung 1 -- arc-length reparametrization

The mildest operation re-labels time *within each curve* and uses no template at
all. `tf_reparam_arclength()` traverses every curve at approximately constant
speed in its value space: the `(hip, knee)` path is unchanged -- the same set of
points -- but equal time intervals now cover equal arc length.

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

`tf_speed()` makes the effect explicit. The raw speeds swing wildly (fast
through swing, near-still in stance); after reparametrization they are nearly
flat -- each curve traversed at constant speed:

```{r 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)
```

No curve is moved toward any other: this removes the *parametrization* only,
leaving shape and size untouched and aligning nothing to a common reference.

### Rung 2 -- alignment to a 1-d reference signal

The next rung aligns curves *to each other* by estimating one shared time-warp
per curve and applying it to both components. The warp is driven by a single
univariate **reference signal** -- a component, or any derived quantity. The
choice matters: the knee carries the cycle's sharpest event (heel-strike), so it
locks the alignment best.

```{r 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)
```

The knee component visibly tightens; the hip, phase-coupled to it, tightens only
incidentally. Choosing `ref_component = "hip"` or `= tf_speed` would emphasise
different events -- a 1-d reference forces that modelling choice. The warps
themselves are univariate `tfd`s mapping raw to aligned phase, fanning out around
the swing phase where heel-strike timing varies between subjects:

```{r 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)
```

### Rung 3 -- joint multivariate reparametrization (`srvf_mv`)

`method = "srvf_mv"` removes the choice of reference: it estimates the single
time-warp that best aligns the joint `(hip, knee)` trajectories in the elastic
(square-root-velocity) sense, using all components at once
[@srivastava2011registration; @tucker2013generative]. It is still only a
re-timing -- values are never rotated or rescaled -- and the template is the
multivariate Karcher mean in `tf_template(reg_mv)`.

```{r 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)
```

```{r 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)
```

Because the warp must compromise between both components rather than chase one,
it need not shrink any single channel as hard as a reference-targeted warp -- but
it needs no reference choice and respects both axes symmetrically.

### Rung 4 -- full elastic shape registration

`tf_register_shape()` goes furthest: on top of the time-warp it also fits, per
curve, a **rotation** and a **scale**, and reports the aligned curves in a
centered, normalised **shape space** [@srivastava2011shape; @srivastava2016].
Translation, rotation and size are all quotiented out.

```{r 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
```

The `(hip, knee)` trajectory view makes the collapse vivid -- the original
butterfly loops on the left, the shape-registered curves on the right reduced to
essentially one normalised shape (note the very different axis ranges):

```{r 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)
```

The estimated rotations sit near the identity and the scales near 1 (little
genuine rotation or scaling is present), yet the normalisation alone flattens
the real between-subject differences:

```{r 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
```

**Gait's two axes are not interchangeable**: hip angle and knee angle are
distinct physical quantities in fixed units. Quotienting out rotation mixes them
into meaningless combinations, and quotienting out scale discards amplitude --
which here *is* the signal of interest. Shape registration is the wrong tool for
a bundle of fixed-unit channels; we return to where it *is* right below.

### Quantifying the alignment

Two views across the rungs. First, a visual overview: the two components (rows
*hip*, *knee*) and the estimated time-warp (row *warp*) for each registration,
side by side -- *as-observed*, *arc-length* reparametrized, registered to a 1-d
*reference* (knee), registered jointly (*srvf_mv*), and fully *shape*-registered.
Raw data has no warp; arc-length has no template warp but does reparametrize, so
its warp panel shows the normalised cumulative arc length.

```{r 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)
}
```

```{r 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)
```

Across the warp row, the registration warps fan out around the swing phase where
heel-strike timing varies between subjects; the shape warp is comparable, since
it shares the same time-warp step. Down the `shape` column the components have
been rescaled into normalised shape space -- note the compressed vertical axis --
which is why those curves look so much tighter than the others.

Second, the numbers: the maximum pointwise sd of each component, before and
after. Rungs 1--3 keep the data in its original degree units and are directly
comparable; rung 4 (marked `*`) lives in normalised shape space, so its row is on
a different scale and only illustrates the collapse.

```{r 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)))
}
```

```{r 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
)
```

In degree units the phase rungs chip away at the spread only modestly -- gait's
phase variation is real but subtle. The knee-reference warp shrinks the knee
component most (it targets that channel); `srvf_mv` spreads a gentler correction
across both; and arc-length reparametrization, which aligns nothing to a common
template, barely moves the pointwise sd at all. Rung 4's near-zero entries are
*not* a better alignment -- they are the signal being quotiented away.

## Shape registration and its quotient spaces

Shape registration earns its keep when the components really are *interchangeable
spatial coordinates* and position, orientation and size are all nuisances:
handwriting, gesture or movement paths, object outlines. Here is a synthetic
example -- a single base curve $(t,\, t^2)$, copied three times with each copy
rotated, rescaled and shifted in the plane:

```{r 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)
```

What counts as "the same shape" depends on which transformations we quotient
out, and `tf_register_shape()` lets us choose via the `rotation` and `scale`
flags (translation and the time-warp are always removed). The three quotients
differ in what survives (again, these chunks need the suggested **fdasrvf**
package):

```{r 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)
```

* **rotation-only** (`scale = FALSE`): orientations are aligned but the original
  sizes are kept, so the curves stay nested.
* **scale-only** (`rotation = FALSE`): sizes are equalised but orientations are
  kept, so the curves fan out at one common size.
* **rotation + scale** (the full shape quotient): both are removed, and the three
  congruent curves collapse onto a single shape.

`tf_scales()` reports the per-curve size factors that were removed (relative to
the template; `1` whenever scaling is off):

```{r 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))
```

## Modes of variation via FPC

Back on the gait sample, with the most visible phase variation reduced (we carry
the knee-aligned curves forward), fit a per-component FPC basis
[@ramsay2005functional]. The leading PCs now mostly describe within-component
amplitude modes:

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

```{r 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)
```

The first 2-3 FPCs explain most of the remaining within-component variance,
and the RMSE between the registered and FPC-reconstructed `(hip, knee)`
curves is well below a degree for most subjects.

## Joint modes of variation via MFPCA

The per-component FPC above gives each component its *own* scores: a subject's
hip-mode-1 score and knee-mode-1 score are separate numbers, so the basis
cannot, on its own, express that the two joints co-vary. **Multivariate FPCA**
[@happ2018], via `tfb_mfpc()`, instead produces a *single* score per
subject per mode, paired with a **vector-valued eigenfunction**
$\Psi_m = (\Psi_m^{\text{hip}}, \Psi_m^{\text{knee}})$. One shared score
$s_{im}$ scales both halves at once, so each mode encodes coupled hip-knee
co-variation in a single coordinate system:
$$ f_i(t) \approx \mu(t) + \sum_m s_{im}\,\Psi_m(t). $$

The default `weights = "inverse_variance"` normalises each component to equal
total variance before combining, so the (larger-range) knee does not simply
dominate the hip; `"snr"`, `"equal"`, or a numeric vector are also available.

```{r 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)
```

The multivariate eigenfunctions come back as a `tfd_mv` -- one bivariate "curve"
per mode -- so each MFPC can be inspected component by component. The leading
mode's hip and knee parts move together, which is exactly the coupling the
independent per-component analysis could not represent:

```{r 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)
```

The shared scores live in one coordinate system, so MFPCA summarises each
subject with a single set of `M` numbers rather than one set per component.
That is a genuine trade-off: the independent per-component basis is free to
optimise each component separately, so it reconstructs more accurately for a
given per-component truncation, while MFPCA spends fewer, *coupled* scores and
accepts a little more error in exchange for the shared coordinate system:

```{r 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)
)
```

```{r 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."
))
```

# 3. Atlantic storms as 4-dimensional curves

`dplyr::storms` records four time-varying quantities for every Atlantic
tropical storm or hurricane from `r min(dplyr::storms$year)` to
`r max(dplyr::storms$year)`, sampled every 6 hours along the storm's life:
position (longitude, latitude) and intensity (sustained wind speed, central
pressure). The natural object is a single
**four-component** vector-valued curve per storm, `f_i: [0, T_i] -> R^4`
-- spatial trajectory and intensity life-cycle bundled together.

Two practical points:

* `tracks_km` -- the *spatial* `(x_km, y_km)` view in physical units
  (projecting each storm's `long` / `lat` into a per-storm local-km
  frame at its own mean latitude) -- is what we use for arc length and
  forward speed, because deg/h speeds wrongly upweight east-west motion
  at high latitudes (a longitude degree is ~111 km at the equator but
  only ~78 km at 45 N).
* Most life-cycle comparisons between short and long storms only make
  sense on a **normalised time axis** `phase = t / T_i in [0, 1]`. We
  build a normalised 4-d object on top of the real-time one and use the
  appropriate version per question.

```{r 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+")
  ))
```

Before analysing the storms, inspect the irregular object directly. Counts
vary by storm because each lifetime has a different number of 6-hourly
observations, while the long data-frame representation keeps the components
aligned by `(storm, phase)`.

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

## Movement workflow: regularize, then describe

@joo2019 describe movement analysis as a workflow built around
tracking records `(x, y, t)`: clean the fixes, regularize or reconstruct
the path when needed, visualize the track, then extract descriptors such
as speed, heading and turning angles. The storm data follow the same
pattern. The build step above already performed basic preprocessing:
duplicate timestamps were removed, very short tracks were dropped, and
longitude/latitude were projected into local kilometre coordinates.

For diagnostic plots it is useful to regularize a few storms onto a
common lifecycle grid. Differentiating the regularized `(x, y)` curves
gives velocity components; their norm is forward speed, and `atan2(v_y,
v_x)` gives heading. Turning angle is the wrapped difference between
successive headings.

```{r 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)
```

The result is a movement-data diagnostic rather than just a map: storms
with visually similar tracks can differ sharply in when they accelerate,
when they turn, and how abruptly their heading changes during recurvature.

## A single 4-d curve: Hurricane Katrina (2005)

For a single storm the 4-component `tf_mv` is naturally displayed
in facet mode (`type = "facet"`) -- one panel per component, all on the
same normalised cycle phase axis:

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

Wind and pressure are anti-correlated -- pressure dips below 920 mbar as
wind peaks above 150 knots near phase ~0.55 -- and the `(long, lat)`
panels show Katrina sliding northwest across the Gulf and turning north
into Louisiana. Bundling intensity and trajectory as a single `tfd_mv`
lets every operation downstream (subsetting, summaries, arithmetic,
plotting) treat them as one object.

## Track map, faceted by peak intensity

The geographic view uses the `(long, lat)` components only. Faceting by
peak Saffir-Simpson category separates the populations cleanly, and a
coastline backdrop (from `maps::map`, if installed) anchors the geography:

```{r 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)
}
```

```{r 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)
```

Tropical storms / depressions stay in the southern, western basin; cat-3
and cat-4+ tracks reach further north and east -- they live longer, get
caught by the westerlies, and recurve out to sea.

## Scalar features extracted from the 4-d object

Every component of a `tf_mv` is a real `tf` vector and supports the full
univariate API, so per-storm scalar features fall out as one-liners:

```{r 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"
  )
```

Median path length grows strongly from tropical-storm to cat-4+ storms,
whereas lifetime-average forward speed increases more modestly. The
canonical wind / pressure gap between categories is visible.

## Intensity and forward-speed life-cycles per category

Plotting all forward-speed and intensity time courses at once is hopeless
(hundreds of irregular noodles). What we actually want is the
*mean curve per intensity stratum* on the normalised lifecycle phase
axis. Evaluate each component on a common phase grid, average within each
stratum, and re-wrap the result as a length-`G` univariate `tfd` so we
can use the standard `plot.tf` / `lines.tf` machinery:

```{r 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)
```

```{r 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)
```

Three clean stories on a single `tf_mv` object: intensity peaks
mid-life-cycle (around `phase ~ 0.5`) and the gap between categories is
roughly uniform across the cycle; central pressure mirrors wind; forward
speed peaks much later (around `phase ~ 0.75`) and the strongest storms
accelerate the most -- the canonical recurvature signature.

## Smoothing on normalised time

Fitting a basis representation on raw real-time data has a tractable
domain-extrapolation problem -- the longest storm runs 600+ hours, but
most tracks end after 100-200 hours, and a shared basis over [0, 600]
extrapolates wildly. On the **normalised phase** version this just
disappears, because every storm spans `[0, 1]`:

```{r 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)
```

The smoothed trajectories pick up the gross recurvature shape without the
6-hourly sampling jitter and -- thanks to the normalised time axis --
without any extrapolation pathology. The same `tb` object also carries
clean smoothed wind / pressure curves for each storm; for example:

```{r 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)
```

# 4. Recap

The two case studies above exercised exactly the same surface:

* construction from per-component `tf` vectors (`tfd_mv(list(...))`), from a 3-d
  `[curve, arg, component]` array, or from a long/wide data frame;
* `tf_components()`, `$component`, `tf_ncomp()` for component access;
* facet vs. trajectory plotting (`plot(..., type = ...)`);
* vctrs-native arithmetic / `mean()` / `sd()` returning a length-1 `tf_mv`;
* `tf_arclength()` and `tf_speed()` as geometric primitives on the bundle;
* `tfb_mv()` for smoothing (`basis = "spline"`) or per-component PC
  decomposition (`basis = "fpc"`), and `tfb_mfpc()` for joint multivariate
  FPCA with a single set of shared scores per curve [@happ2018];
* the alignment ladder: `tf_reparam_arclength()` for constant-speed (shape-only)
  reparametrization, `tf_register()` for shared-warp alignment -- from a single
  reference component (`method = "cc"`) or jointly from all components
  (`method = "srvf_mv"`) -- and `tf_register_shape()` for full elastic shape
  registration, whose `tf_rotations()` / `tf_scales()` accessors expose the
  estimated rotations and (template-relative) scale factors;
* full dplyr / tibble integration via the underlying vctrs proxy.

The common pattern is simple: keep coupled signals bundled as one object, use
component accessors when a scalar or univariate summary is needed, and let the
multivariate methods reuse the existing univariate kernels component-wise.

# References
