## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE, fig.width = 6, fig.height = 4 ) ## ----load-pkg----------------------------------------------------------------- library(adjoin) library(Matrix) ## ----create-grid-------------------------------------------------------------- coords <- as.matrix(expand.grid(x = 1:6, y = 1:6)) # 36 points, 2 columns dim(coords) ## ----first-spatial------------------------------------------------------------ A <- spatial_adjacency(coords, sigma = 1.5, nnk = 8, weight_mode = "heat", include_diagonal = FALSE, normalized = TRUE) cat("size:", nrow(A), "x", ncol(A), "| non-zero entries:", Matrix::nnzero(A), "\n") cat("symmetric:", isSymmetric(A), "\n") ## ----first-spatial-plot, echo = FALSE, fig.cap = "Spatial neighbor graph on a jittered 6 × 6 grid. Thicker, darker lines are higher-weight edges (nearer neighbours); lighter lines are lower-weight (farther neighbours). Point colour encodes degree.", out.width = "70%"---- # Jitter the grid so neighbour distances vary and heat-kernel weights differ set.seed(7) coords_j <- coords + matrix(runif(nrow(coords) * 2, -0.35, 0.35), ncol = 2) # Raw (unnormalized) weights for visual contrast: nearer pairs score higher A_j <- spatial_adjacency(coords_j, sigma = 1.5, nnk = 8, weight_mode = "heat", include_diagonal = FALSE, normalized = FALSE) deg <- as.numeric(Matrix::rowSums(A_j > 0)) pal <- colorRampPalette(c("#deebf7", "#08519c"))(max(deg) + 1) pcol <- pal[deg + 1] plot(coords_j, pch = 16, cex = 1.3, asp = 1, axes = FALSE, xlab = "", ylab = "", main = "Spatial neighbor graph (sigma = 1.5, k \u2264 8, jittered)", col = pcol) box() A_t <- as(A_j, "TsparseMatrix") w_all <- A_t@x[A_t@x > 0] w_min <- min(w_all); w_max <- max(w_all) for (k in seq_along(A_t@x)) { i <- A_t@i[k] + 1L j <- A_t@j[k] + 1L if (i < j) { w <- A_t@x[k] # rescale weight to [0.1, 1] for alpha; heavier = more opaque & thicker alpha_k <- 0.15 + 0.85 * (w - w_min) / (w_max - w_min + 1e-9) lwd_k <- 0.5 + 2.0 * (w - w_min) / (w_max - w_min + 1e-9) segments(coords_j[i, 1], coords_j[i, 2], coords_j[j, 1], coords_j[j, 2], col = adjustcolor("#2171b5", alpha.f = alpha_k), lwd = lwd_k) } } points(coords_j, pch = 16, cex = 1.3, col = pcol) legend("topright", legend = c("fewer neighbors", "more neighbors"), fill = c("#deebf7", "#08519c"), bty = "n", cex = 0.8) ## ----weight-modes------------------------------------------------------------- A_heat <- spatial_adjacency(coords, sigma = 1.5, nnk = 8, weight_mode = "heat", include_diagonal = FALSE, normalized = FALSE, stochastic = FALSE) A_binary <- spatial_adjacency(coords, sigma = 1.5, nnk = 8, weight_mode = "binary", include_diagonal = FALSE, normalized = FALSE, stochastic = FALSE) ## ----weight-mode-plot, echo = FALSE, fig.cap = "Heat weights decay continuously with distance; binary weights are uniform within the neighborhood radius.", fig.width = 7, fig.height = 3---- local({ oldpar <- par(no.readonly = TRUE) tryCatch({ par(mfrow = c(1, 2), mar = c(4, 3, 2.5, 1)) w_heat <- A_heat@x[A_heat@x > 0] hist(w_heat, breaks = 20, col = "#2171b5", border = "white", main = "heat kernel", xlab = "edge weight", ylab = "count") w_bin <- A_binary@x[A_binary@x > 0] hist(w_bin, breaks = 5, col = "#6baed6", border = "white", main = "binary", xlab = "edge weight", ylab = "count") }, finally = { par(oldpar) }) }) ## ----sigma-effect------------------------------------------------------------- A_tight <- spatial_adjacency(coords, sigma = 0.8, nnk = 27, weight_mode = "heat", include_diagonal = FALSE, normalized = TRUE) A_wide <- spatial_adjacency(coords, sigma = 2.5, nnk = 27, weight_mode = "heat", include_diagonal = FALSE, normalized = TRUE) c(tight_edges = Matrix::nnzero(A_tight) / 2L, wide_edges = Matrix::nnzero(A_wide) / 2L) ## ----sigma-effect-plot, echo = FALSE, fig.cap = "Tight sigma (0.8) leaves corner and edge points with few connections; wide sigma (2.5) connects most points richly.", fig.width = 7, fig.height = 3.5---- local({ oldpar <- par(no.readonly = TRUE) tryCatch({ par(mfrow = c(1, 2), mar = c(1, 1, 2.5, 1)) plot_graph <- function(A, coords, title) { deg <- as.numeric(Matrix::rowSums(A > 0)) pal <- colorRampPalette(c("#deebf7", "#08519c"))(max(deg, 1) + 1) cols <- pal[deg + 1] plot(coords, pch = 16, cex = 1.3, asp = 1, axes = FALSE, xlab = "", ylab = "", main = title, col = cols) A_t <- as(A, "TsparseMatrix") for (k in seq_along(A_t@x)) { i <- A_t@i[k] + 1L j <- A_t@j[k] + 1L if (i < j) segments(coords[i, 1], coords[i, 2], coords[j, 1], coords[j, 2], col = "#c6dbef", lwd = 0.8) } points(coords, pch = 16, cex = 1.3, col = cols) } plot_graph(A_tight, coords, "sigma = 0.8 (tight)") plot_graph(A_wide, coords, "sigma = 2.5 (wide)") }, finally = { par(oldpar) }) }) ## ----spatial-smoother--------------------------------------------------------- S <- spatial_smoother(coords, sigma = 1.5, nnk = 8, stochastic = FALSE) set.seed(42) signal_raw <- rnorm(36) signal_smoothed <- as.numeric(S %*% signal_raw) ## ----smooth-plot, echo = FALSE, fig.cap = "One pass of the spatial smoother turns Gaussian noise (left) into a smoothly varying field (right). The same heat kernel defines both neighbor selection and weight decay.", fig.width = 7, fig.height = 3.5---- local({ oldpar <- par(no.readonly = TRUE) tryCatch({ par(mfrow = c(1, 2), mar = c(2, 2, 2.5, 1)) zlim_range <- range(c(signal_raw, signal_smoothed)) heat_pal <- colorRampPalette(c("#d73027", "white", "#4575b4"))(50) image(1:6, 1:6, matrix(signal_raw, 6, 6), zlim = zlim_range, col = heat_pal, main = "Raw signal", xlab = "x", ylab = "y") image(1:6, 1:6, matrix(signal_smoothed, 6, 6), zlim = zlim_range, col = heat_pal, main = "After spatial smooth", xlab = "x", ylab = "y") }, finally = { par(oldpar) }) }) ## ----spatial-laplacian-------------------------------------------------------- L <- spatial_laplacian(coords, dthresh = 4.5, nnk = 8, weight_mode = "binary", normalized = FALSE) # Row sums of a Laplacian are zero by definition range(round(rowSums(L), 10)) ## ----weighted-adj------------------------------------------------------------- set.seed(42) features <- matrix(rnorm(36 * 5), nrow = 36, ncol = 5) # random 5-d features A_sp <- weighted_spatial_adjacency(coords, features, alpha = 0, sigma = 2, nnk = 8) A_mix <- weighted_spatial_adjacency(coords, features, alpha = 0.5, sigma = 2, nnk = 8) A_ft <- weighted_spatial_adjacency(coords, features, alpha = 1, sigma = 2, nnk = 8) ## ----weighted-adj-plot, echo = FALSE, fig.cap = "Adjacency heat maps for three alpha values. Spatial-only (left) shows a clean block structure; feature-only (right) looks noisier because the random features carry no spatial signal.", fig.width = 9, fig.height = 3---- local({ oldpar <- par(no.readonly = TRUE) tryCatch({ par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 0.5)) show_adj <- function(A, title) { image(seq(0, 1, length.out = 36), seq(0, 1, length.out = 36), as.matrix(A), col = colorRampPalette(c("white", "#2c7bb6"))(40), axes = FALSE, xlab = "", ylab = "", main = title) } show_adj(A_sp, "alpha = 0\n(spatial only)") show_adj(A_mix, "alpha = 0.5\n(blend)") show_adj(A_ft, "alpha = 1\n(feature only)") }, finally = { par(oldpar) }) }) ## ----bilateral---------------------------------------------------------------- # Piecewise-constant signal with a hard boundary at x = 3.5 signal_field <- ifelse(coords[, "x"] <= 3, -1, 1) + rnorm(36, sd = 0.2) feature_mat <- matrix(signal_field, ncol = 1) B <- bilateral_smoother(coords, feature_mat, nnk = 8, s_sigma = 1.5, f_sigma = 0.5) signal_bilateral <- as.numeric(B %*% signal_field) signal_spatial <- as.numeric(S %*% signal_field) # plain spatial smoother ## ----bilateral-plot, echo = FALSE, fig.cap = "Bilateral smoother (right) preserves the hard boundary between the two signal regions; a plain spatial smoother (left) blurs across it.", fig.width = 7, fig.height = 3.5---- local({ oldpar <- par(no.readonly = TRUE) tryCatch({ par(mfrow = c(1, 2), mar = c(2, 2, 2.5, 1)) zlim_b <- range(c(signal_spatial, signal_bilateral)) heat_pal <- colorRampPalette(c("#d73027", "white", "#4575b4"))(50) image(1:6, 1:6, matrix(signal_spatial, 6, 6), zlim = zlim_b, col = heat_pal, main = "Spatial smoother\n(blurs boundary)", xlab = "x", ylab = "y") abline(v = 3.5, col = "black", lwd = 2, lty = 2) image(1:6, 1:6, matrix(signal_bilateral, 6, 6), zlim = zlim_b, col = heat_pal, main = "Bilateral smoother\n(preserves boundary)", xlab = "x", ylab = "y") abline(v = 3.5, col = "black", lwd = 2, lty = 2) }, finally = { par(oldpar) }) }) ## ----multi-block-------------------------------------------------------------- coords_small <- as.matrix(expand.grid(x = 1:4, y = 1:4)) # 16-point grid # Pass a list with one entry per block (same grid repeated for two subjects) S2 <- spatial_constraints(list(coords_small, coords_small), nblocks = 2, sigma_within = 1.5, nnk_within = 6, sigma_between = 2.0, nnk_between = 4, shrinkage_factor = 0.15) dim(S2) # 32 x 32: 2 blocks x 16 points ## ----multi-block-plot, echo = FALSE, fig.cap = "32x32 constraint matrix for two 4x4 spatial blocks. Diagonal blocks capture within-block spatial similarity; off-diagonal blocks link corresponding locations across blocks. The leading eigenvalue is 1 by construction.", out.width = "65%"---- image(seq(0, 1, length.out = 32), seq(0, 1, length.out = 32), as.matrix(S2), col = colorRampPalette(c("white", "#2c7bb6"))(50), axes = FALSE, xlab = "", ylab = "", main = "Spatial constraints: 2 blocks x 16 points") abline(v = 0.5, h = 0.5, col = "firebrick", lwd = 1.5) legend("topright", fill = "firebrick", legend = "block boundary", bty = "n", cex = 0.8)