---
title: "Descriptive covariance and landmark applications"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Descriptive covariance and landmark applications}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", error = FALSE,
                     fig.width = 7, fig.height = 4)
library(Riemann)
```

The package contains covariance matrices and landmark configurations that can
illustrate complete fitted workflows. These examples are descriptive. Source
provenance and the meaning of observations limit the scientific claims they can
support. The separate `write-JSS-Riemann/` writing workspace contains
`replication/code.R`, which separately regenerates full data
validation, model objects, figures, tables, and an output manifest.

## Event-related covariance matrices

A reconstruction of the [historical pyRiemann 0.2.7 recipe](https://github.com/pyRiemann/pyRiemann/blob/176e766f540bd4c7846f38573165fc3d27fc69ca/examples/ERP/plot_embedding_EEG.py)
matched all 216 `ERP` matrices, with maximum relative Frobenius error
`5.56e-12` and exact label order, without rescaling, reordering, or sign changes.
These are test epochs from the single participant in the MNE sample audiovisual
experiment. The recipe selects 305 MEG sensors and fits supervised xDAWN filters
on 72 separate training epochs. Its 32 covariance coordinates combine 16
class-prototype and 16 filtered trial features; they are not original sensor
channels. See `?ERP` for filtering, epoch, split, and covariance settings.
The [public MNE sample record](https://doi.org/10.18112/openneuro.ds000248.v1.2.4)
declares CC0. The reconstruction used recorded contemporary dependencies, so
the numerical match is not a claim of byte identity or recovery of the original
runtime. These descriptive examples do not treat 216 epochs as independent
participants, and a later split of the stored matrices does not recreate the
label-dependent feature-learning step or establish new-participant performance.

```{r erp-validation}
data(ERP)
dim(ERP$covariance)
table(ERP$label)
x <- wrap.spd(ERP$covariance)
spectral <- t(vapply(x$data, function(a) {
  values <- eigen(a, symmetric = TRUE, only.values = TRUE)$values
  c(minimum = min(values), maximum = max(values), ratio = min(values) / max(values))
}, numeric(3)))
stopifnot(all(spectral[, "minimum"] > 0))
apply(spectral, 2, range)
```

No ridge, scaling, or guessed physical-unit conversion is applied. Very small
absolute eigenvalues do not themselves imply poor conditioning; the dimensionless
minimum/maximum ratio is recorded separately.

Fit two geometries to the same observations and retain their complete references.
The mean diagnostics report numerical convergence.

```{r erp-fits}
airm <- riem.pga(x, ndim = 2, geometry = "affine_invariant", maxiter = 300, eps = 1e-7)
loge <- riem.pga(x, ndim = 2, geometry = "log_euclidean", maxiter = 300, eps = 1e-7)
stopifnot(airm$diagnostics$converged, loge$diagnostics$converged)
data.frame(
  geometry = c("affine_invariant", "log_euclidean"),
  mean_objective = c(airm$diagnostics$mean$objective, loge$diagnostics$mean$objective),
  tangent_variance_retained = c(sum(airm$explained.variance), sum(loge$explained.variance))
)
```

The two-dimensional displays retain only
`r sprintf('%.1f%%', 100 * sum(airm$explained.variance))` and
`r sprintf('%.1f%%', 100 * sum(loge$explained.variance))` of their respective tangent
variances. They are partial views of a higher-dimensional sample; their appearance
cannot establish classification performance or a neurophysiological difference.
Objectives under different geometries are different estimands, so their numerical
magnitudes alone do not rank the geometries.

```{r erp-display}
colors <- c("#3565a0", "#bd5845", "#43886c", "#84669b")
old <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(airm$embed, col = colors[as.integer(ERP$label)], pch = 19, cex = 0.6,
     main = "Affine-invariant", xlab = "Tangent PC1", ylab = "Tangent PC2")
plot(loge$embed, col = colors[as.integer(ERP$label)], pch = 19, cex = 0.6,
     main = "Log-Euclidean", xlab = "Tangent PC1", ylab = "Tangent PC2")
legend("bottomleft", levels(ERP$label), col = colors, pch = 19, bty = "n", cex = 0.7)
par(old)
```

The same saved object can transform observations one at a time. This is a software
reuse check, not an independent test set.

```{r erp-reuse}
one_matrix <- wrap.spd(list(ERP$covariance[, , 1]))
stopifnot(max(abs(predict(loge, one_matrix) - loge$embed[1, , drop = FALSE])) < 1e-8)
reconstructed <- riem.reconstruct(loge, predict(loge, one_matrix))
stopifnot(min(eigen(reconstructed$data[[1]], symmetric = TRUE)$values) > 0)
```

For a modest downstream task, fit descriptive clusters in the declared metric.
The label/cluster table is a cross-tabulation on the same observations, not held-out
prediction accuracy. No population-level uncertainty estimate is justified here.

```{r erp-clustering}
set.seed(91401)
clusters <- riem.kmeans(x, k = 4, geometry = "log_euclidean", nstart = 2,
                        maxiter = 100, mean.eps = 1e-8)
stopifnot(clusters$converged, identical(predict(clusters, x), clusters$cluster))
clusters$starts
table(event_label = ERP$label, descriptive_cluster = clusters$cluster)
```

## Gorilla skull shape

The arrays contain 29 male and 30 female configurations with eight corresponding
landmarks in two dimensions. Their exact upstream mapping was checked against
`shapes` 1.2.8: reorder `gorm.dat` and `gorf.dat` by
`c(1, 5, 4, 3, 2, 8, 7, 6)` on the landmark dimension, retaining coordinate and
individual order. The [shapes manual](https://cran.r-project.org/web/packages/shapes/refman/shapes.html)
attributes these data to Paul O'Higgins and cites O'Higgins and Dryden (1993),
*Journal of Human Evolution* 24, 183–205. The separate writing workspace includes the executable
`maintenance/verify-gorilla-provenance.R` array-correspondence verification script.

The workflow analyzes shape after removing location and scale. Its quotient uses
all orthogonal transformations, including reflections. This choice is explicit:
orientation-preserving shape or original skull size would ask a different question.

```{r gorilla-fit}
data(gorilla)
as_list <- function(a) lapply(seq_len(dim(a)[3]), function(i) a[, , i])
raw <- c(as_list(gorilla$male), as_list(gorilla$female))
group <- factor(c(rep("male", 29), rep("female", 30)), levels = c("male", "female"))
shapes <- wrap.landmark(raw)
shape_fit <- riem.pga(shapes, ndim = 2, geometry = "shape_orthogonal",
                      maxiter = 300, eps = 1e-8)
stopifnot(shape_fit$diagnostics$converged)
summary(shape_fit)
```

The retained tangent variance is
`r sprintf('%.1f%%', 100 * sum(shape_fit$explained.variance))`. Sex labels describe
this sample; the example does not perform a population test or estimate a causal
or biological effect.

```{r gorilla-display}
plot(shape_fit$embed, col = c("#3565a0", "#bd5845")[as.integer(group)], pch = 19,
     xlab = "Tangent PC1", ylab = "Tangent PC2")
legend("topleft", levels(group), col = c("#3565a0", "#bd5845"), pch = 19, bty = "n")
```

Verify that the model uses its trained frame when new batches are wrapped.
This check covers positive scaling, translation, rotation, reflection, and
single-observation prediction.

```{r gorilla-invariance}
angle <- 0.6
rotation <- matrix(c(cos(angle), sin(angle), -sin(angle), cos(angle)), 2)
transformed <- lapply(raw, function(a) sweep(2.3 * a %*% rotation, 2, c(30, -5), "+"))
reflected <- lapply(raw, function(a) a %*% diag(c(-1, 1)))
checks <- c(
  location_scale_rotation = max(abs(predict(shape_fit, wrap.landmark(transformed)) - shape_fit$embed)),
  reflection = max(abs(predict(shape_fit, wrap.landmark(reflected)) - shape_fit$embed)),
  singleton = max(abs(do.call(rbind, lapply(raw, function(a) {
    predict(shape_fit, wrap.landmark(list(a)))
  })) - shape_fit$embed))
)
stopifnot(all(checks < 1e-7))
checks
```

Reconstruct a small variation along the first component. Landmark numbers indicate
stored correspondences; no anatomical names or outlines are inferred.

```{r gorilla-mode}
scores <- matrix(0, 3, 2)
scores[, 1] <- c(-2, 0, 2) * sqrt(shape_fit$variance[1])
mode <- riem.reconstruct(shape_fit, scores)$data
plot(mode[[2]], asp = 1, pch = 19,
     xlim = range(vapply(mode, function(a) a[, 1], numeric(8))),
     ylim = range(vapply(mode, function(a) a[, 2], numeric(8))),
     xlab = "Trained preshape coordinate 1",
     ylab = "Trained preshape coordinate 2")
points(mode[[1]], col = "#3565a0", pch = 1)
points(mode[[3]], col = "#bd5845", pch = 2)
text(mode[[2]], labels = seq_len(8), pos = 3, cex = 0.7)
legend("bottomleft", c("-2 tangent SD", "reference", "+2 tangent SD"),
       col = c("#3565a0", "black", "#bd5845"), pch = c(1, 19, 2), bty = "n")
```

These reconstructions remain in the supported regular local shape neighborhood.
Singular alignments and reconstructions outside that chart are explicit errors,
not additional validated cases. Full numerical records and both component-mode
plots are generated by the separate replication entry point.

```{r session}
sessionInfo()
```