## ----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)

