A practical tour of the five most important functions in spinebil for finding, validating, and understanding projection pursuit indices.
Overview
In this vignette we introduce three example datasets; spiral, a sine, and a pipe dataset and then evaluate five key characteristics of projection pursuit indices (PPIs): smoothness, squintability, flexibility, rotation invariance and speed.
We will generate the datasets first, then walk through one section per characteristic with the functions provided in spinebil.
Data Generators
pipe_data(n, p, t); a pipe (circular ring) structure.sin_data(n, p, f); a sine relationship.spiral_data(n, p); an archimedean spiral.
n <- 500
p <- 4
lst <- list(
Pipe = spinebil::pipe_data(n, p),
Sine = spinebil::sin_data(n, p, 1),
Spiral = spinebil::spiral_data(n, p)
)
df_all <- do.call(rbind, lapply(names(lst), function(lbl) {
d <- lst[[lbl]]
data.frame(x = d[[p - 1]], y = d[[p]], structure = lbl)
}))
ggplot2::ggplot(df_all, ggplot2::aes(x, y)) +
ggplot2::geom_point(alpha = 0.6, size = 0.6) +
ggplot2::facet_wrap(~ structure, nrow = 1, scales="free") +
ggplot2::theme(
aspect.ratio = 1,
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
axis.title = ggplot2::element_blank()
)
Now we employ these patterns to evaluate key characteristics of PPIs.
1) compare_smoothing()
Projection pursuit indices evaluated along a tour path can be spiky due to small numerical changes in the projection or to point-level noise. compare_smoothing() provides a principled way to smooth these traces by averaging the index over local perturbations, and to compare different smoothing strategies side-by-side.
It supports two kinds of perturbations:
jitter_angle: randomly jitter the projection plane by a small angle alpha via geodesic moves, then recompute the index.jitter_points: randomly jitter the projected points by a small amount (using base::jitter()), then recompute the index.no_smoothing: the original index with no perturbation.
Averaging across multiple perturbations reduces high-frequency noise and reveals the underlying trend of the index along the tour.
Function Usage
compare_smoothing(
d, # data matrix (n x p)
tPath, # interpolated tour path: list of projection bases (p x 2)
idx, # index function
alphaV = c(0.01, 0.05, 0.1), jitter amounts to compare (for jittering angle or points)
n = 10 # number of evaluations entering mean value calculation
)
d, numeric data matrix with p columns. It will be projected as d %*% basis for each basis in tPath.tPath, list of (p × 2) projection bases. Typically built from a tour history using tourr::save_history() and tourr::interpolate().idx, index function accepting a two-column matrix and returning a single numeric value (e.g., scag_index("stringy")).alphaV, numeric vector of jitter magnitudes to compare.n, number of evaluations entering mean value calculation (larger n → smoother, slower).
Example Usage
d <- as.matrix(spinebil::spiral_data(30,4))
tPath <- tourr::save_history(d, max_bases=2)
tPath <- as.list(tourr::interpolate(tPath, 0.3))
idx <- spinebil::scag_index("stringy")
compS <- spinebil::compare_smoothing(d, tPath, idx, alphaV = c(0.01, 0.05), n=2)
spinebil::plot_smoothing_comparison(compS, lPos = "bottom")
Interpretation
- no_smoothing (red, solid): the raw index trace along the tour. Spiky traces indicate high sensitivity to tiny projection changes.
- jitter_angle (black, dashed): averaging over nearby projection angles. Smooths out high‑frequency variability caused by small plane movements; preserves structure amplitude better when index is robust to minor orientation shifts.
- jitter_points (black, dotted): averaging over point noise. Strong smoothing when the index is sensitive to local point perturbations.
Return value
A tibble with columns:
index_mean, mean index value at each frame (includes the original, unjittered value).t, integer frame index along the tour path.method, one of "jitter_angle", "jitter_points", "no_smoothing".alpha, jitter magnitude used for that row.
2) squint_angle_estimate()
From how far away (in projection space) does the pattern become visible under a chosen index?
squint_angle_estimate(), produces a distribution of squint angles by repeatedly walking from random 2‑D planes toward an assumed optimal plane (the best view of a structure) and recording the first point along each path where a user‑chosen index function exceeds a visibility cutoff.
Interpretation:
- Larger squint angles ⇒ easier to see (the index crosses the cutoff while still far from the optimal plane).
- Smaller squint angles ⇒ harder to see (the index only crosses when we are very close to the optimal plane).
Function usage
squint_angle_estimate(
data, # numeric matrix/data frame (n x p)
indexF, # function: (n x 2) -> numeric scalar
cutoff, # numeric threshold for 'visible'
structure_plane, # 2-D basis (p x 2) representing the optimal view
n = 100, # number of random starts
step_size = 0.01 # interpolation step along the tour path
)
Example Usage
data <- as.matrix(spinebil::spiral_data(50, 4))
indexF <- spinebil::scag_index("stringy")
cutoff <- 0.7
structure_plane <- spinebil::basis_matrix(3,4,4)
spinebil::squint_angle_estimate(data, indexF, cutoff, structure_plane, n=10)
#> [1] 1.678098 1.194991 1.413274 1.552521 1.696683 1.615608 1.790544 1.405470
#> [9] 1.489728 1.701739
indexF (what kind of structure?)
Choose an index that matches the pattern you care about. With scagnostics-style indices, for example:
indexF <- scag_index("stringy") # sine-wave/spiral-like
indexF <- scag_index("skinny") # elongated patterns
cutoff (when do we call it visible?)
Set this data‑driven so results are comparable across datasets and indices.
- Preferred: 95th–99th percentile of the index computed on random projections (a null where no specific structure is targeted).
- You can use
index_noise_threshold() to estimate this cutoff automatically.
structure_plane (where is the best view?)
If you know the two variables that define the structure, construct the basis directly. Otherwise, run a guided tour to maximize the index and use the best basis it returns.
n and step_size
n = 100 typically yields a stable distribution for summaries and plots.step_size = 0.01 is a good accuracy/speed trade‑off. If the threshold crossing looks coarse (the index jumps over the cutoff), try 0.005. If runtime is an issue, use something larger (e.g., 0.02).
Return value
squint_angle_estimate() returns a numeric vector of length n containing the squint-angle estimates under the specified index and cutoff.
3) get_trace()
get_trace() evaluates one or more projection pursuit indices along an interpolated, planned tour path and returns their values at each frame. Plotting these traces reveals whether an index varies smoothly with small changes in the projection (desirable), or exhibits spikes (potentially unstable, overly sensitive, or under‑smoothed).
A smooth trace indicates that small rotations of the view produce small changes in the index,an important property for guided tours and optimisation.
In combination with plot_trace(), you can quickly diagnose index behaviour across a path connecting user‑specified views.
Function Usage
get_trace(
d, # data: matrix/data frame (n x p)
m, # list of projection matrices for the planned tour
index_list, # list of index functions to calculate for each entry -> numeric
index_labels # character vector of labels for the indices
)
d; numeric data with p columns.m; list of (p × 2) projection matrices for the planned tour.index_list; list of functions that take a two‑column matrix and return a numeric scalar (e.g., tourr::holes(), tourr::cmass(), or scag_index("stringy")).index_labels; character vector of the same length and order as index_list; used as column names in the output.
Example Usage
d <- as.matrix(spinebil::spiral_data(100, 4))
m <- list(spinebil::basis_matrix(1,2,4), spinebil::basis_matrix(3,4,4))
index_list <- list(tourr::holes(), tourr::norm_kol(100))
index_labels <- c("holes", "norm kol")
trace <- spinebil::get_trace(d, m, index_list, index_labels)
spinebil::plot_trace(trace)
spinebil::plot_trace(trace, rescY = FALSE)
Return value
A numeric matrix with length(index_labels) + 1 columns and as many rows as interpolation frames. Columns are the index values (named by index_labels) and t (the frame index).
4) profile_rotation()
Does the index value stay the same when the 2‑D data are rotated?
profile_rotation() tests rotation invariance of one or more 2-D projection indices.
Interpretation
- Flat profile → rotation invariant.
- Oscillating profile → orientation dependence.
Function Usage
profile_rotation(
d, # 2-column numeric matrix (the data to rotate)
index_list, # list of functions: (n x 2) -> numeric
index_labels, # character labels for columns
n = 200 # number of rotation steps across [0, 2*pi]
)
d: numeric matrix with 2 columns; the 2-D data to be rotated.index_list: list of functions, each taking an (n × 2) numeric matrix and returning a single numeric index value. Examples: list(tourr::holes(), scag_index("stringy"), mine_indexE("MIC")).index_labels: character vector of labels (one per index in index_list) used as column names in the result.n (default = 200): integer number of rotation steps.
Example Usage
d <- as.matrix(spinebil::sin_data(30, 2))
index_list <- list(tourr::holes(), spinebil::scag_index("stringy"), spinebil::mine_indexE("MIC"))
index_labels <- c("holes", "stringy", "mic")
pRot <- spinebil::profile_rotation(d, index_list, index_labels, n = 50)
spinebil::plot_rotation(pRot)
Interpretation
- Perfect circle (constant radius) → rotation invariant.
- Deformed/flower-like shape → orientation dependence.
Return value
A numeric matrix with n + 1 rows and length(index_labels) + 1 columns:
- One column per index (named by
index_labels) containing the index values at each angle.
- An additional column
alpha giving the corresponding angles in radians.
5) time_sequence()
The cost of evaluating a projection pursuit index can vary with both the data distribution and the projection. time_sequence() times the index on a sequence of projection bases and returns a simple table you can plot or summarise.
Function Usage
time_sequence(
d, # numeric data matrix (n x p)
t, # list of projection matrices (each p x 2); e.g., an interpolated tour path
idx, # index function: (n x 2) -> numeric
pmax # maximum number of projections to evaluate (cut t if longer than pmax)
)
d: numeric matrix of size n × p.t: list of (p × 2) projection bases (e.g., from tourr::basis_random() or an interpolated tour history).idx: index function mapping a two-column matrix to a single numeric value (e.g., scag_index("stringy")).pmax: integer limit; evaluation stops after pmax projections even if t is longer.
Example Usage
d <- as.matrix(spinebil::spiral_data(500, 4))
t <- purrr::map(1:10, ~ tourr::basis_random(4))
idx <- spinebil::scag_index("stringy")
spinebil::time_sequence(d, t, idx, 10)
#> t i
#> 1 0.039 1
#> 2 0.035 2
#> 3 0.115 3
#> 4 0.035 4
#> 5 0.037 5
#> 6 0.037 6
#> 7 0.034 7
#> 8 0.037 8
#> 9 0.034 9
#> 10 0.038 10
Return value
A data frame with two columns:
t: elapsed time in seconds for the index evaluation at that projection.i: the sequence index (1, 2, …) of the projection in t.
