## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE,
  fig.width = 7.2,
  fig.height = 3.2,
  fig.align = "center",
  dpi = 110
)
library(fibermargin)

field_palette <- function(labels) {
  classes <- sort(unique(as.character(labels[!is.na(labels)])))
  stats::setNames(grDevices::hcl.colors(length(classes), "Dark 3"), classes)
}

plot_field <- function(xy, labels, main = NULL, cex = 0.55) {
  palette <- field_palette(labels)
  graphics::plot(
    xy[, 1L], xy[, 2L],
    col = palette[as.character(labels)],
    pch = 16, cex = cex, asp = 1,
    xlab = colnames(xy)[1L], ylab = colnames(xy)[2L],
    main = main
  )
}

plot_mask <- function(mask, main = NULL, palette = field_palette(mask)) {
  values <- matrix(
    match(as.character(mask), names(palette)),
    nrow = nrow(mask), ncol = ncol(mask)
  )
  graphics::image(
    seq_len(nrow(mask)), seq_len(ncol(mask)), values,
    col = unname(palette), useRaster = TRUE, axes = FALSE,
    xlab = "", ylab = "", main = main, asp = 1
  )
}

## ----api-index----------------------------------------------------------------
api <- data.frame(
  Function = c(
    "refine_spatial_labels", "clean_categorical_mask",
    "corrupt_categorical_mask", "evaluate_mask_cleaning",
    "spatial_benchmark", "evaluate_spatial_refinement",
    "benchmark_spatial_refiners", "simulate_spatial_clusters",
    "simulate_spatial_domains", "simulate_complex_spatial_domains",
    "simulate_gradient_regions", "simulate_volumetric_domains",
    "available_spatial_benchmarks", "load_spatial_benchmark"
  ),
  Purpose = c(
    "Refine labels at irregular 2D/3D coordinates",
    "Clean a categorical matrix or 3D array",
    "Generate reproducible mask errors",
    "Score a cleaned mask",
    "Validate and package a benchmark",
    "Score coordinate-indexed refinement",
    "Run identical inputs through multiple methods",
    "Simulate separated Gaussian-like clusters",
    "Simulate structured 2D or curved-layer 3D domains",
    "Simulate held-out complex geometries",
    "Simulate the A-B-C gradient design",
    "Simulate variable-density 3D domains",
    "List real-data availability and licensing",
    "Load a bundled frozen real-data scenario"
  ),
  stringsAsFactors = FALSE
)
knitr::kable(api)

## ----quick-start--------------------------------------------------------------
gradient <- simulate_gradient_regions(
  n = 1200L,
  minority = 0.10,
  dimensions = 2L,
  samples = 2L,
  density_profile = "moderate",
  seed = 12L
)

refined_gradient <- refine_spatial_labels(
  xy = gradient$xy,
  labels = gradient$labels,
  samples = gradient$samples,
  workers = 1L
)

gradient_metrics <- evaluate_spatial_refinement(
  truth = gradient$truth,
  initial = gradient$labels,
  refined = refined_gradient,
  boundary = gradient$boundary,
  regions = gradient$area,
  sparse = gradient$sparse,
  method = "FiberMargin"
)

gradient_metrics[, c(
  "method", "initial_accuracy", "accuracy", "ari",
  "correction_recall", "damage_rate"
)]

## ----quick-start-plot, fig.height=3.0-----------------------------------------
old_par <- graphics::par(no.readonly = TRUE)
graphics::par(mfrow = c(1, 3), mar = c(3, 3, 2, 1))
plot_field(gradient$xy, gradient$truth, "Reference", cex = 0.35)
plot_field(gradient$xy, gradient$labels, "Observed", cex = 0.35)
plot_field(gradient$xy, refined_gradient, "FiberMargin", cex = 0.35)
graphics::par(old_par)

## ----diagnostics--------------------------------------------------------------
diagnostic_names <- c(
  "candidate", "margin_score", "required", "repair_margin",
  "atlas_dispersion", "isolation", "changed"
)

diagnostics <- data.frame(
  observed = head(as.character(gradient$labels), 8L),
  refined = head(as.character(refined_gradient), 8L),
  candidate = head(as.character(attr(refined_gradient, "candidate")), 8L),
  margin_score = head(attr(refined_gradient, "margin_score"), 8L),
  required = head(attr(refined_gradient, "required"), 8L),
  repair_margin = head(attr(refined_gradient, "repair_margin"), 8L),
  changed = head(attr(refined_gradient, "changed"), 8L)
)
diagnostics

setdiff(diagnostic_names, names(attributes(refined_gradient)))

## ----summary-diagnostics------------------------------------------------------
summary_names <- c(
  "workers", "dimensions_used", "labels_changed", "changed_fraction",
  "classes_before", "classes_after", "removed_classes", "sample_sizes"
)
lapply(summary_names, function(name) attr(refined_gradient, name))

## ----diagnostic-contract------------------------------------------------------
identical(
  as.logical(attr(refined_gradient, "changed")),
  as.character(refined_gradient) != as.character(gradient$labels)
)
attr(refined_gradient, "workers")

## ----binary-specialization----------------------------------------------------
binary_data <- simulate_spatial_clusters(
  n = 700L,
  dimensions = 2L,
  k = 2L,
  samples = 1L,
  noise = 0.12,
  seed = 21L
)
binary_refined <- refine_spatial_labels(
  binary_data$xy, binary_data$labels,
  samples = binary_data$samples, workers = 1L
)

c(
  input_accuracy = mean(binary_data$labels == binary_data$truth),
  refined_accuracy = mean(binary_refined == binary_data$truth),
  changed_fraction = mean(binary_refined != binary_data$labels)
)

## ----constant-z---------------------------------------------------------------
flat_2d <- refine_spatial_labels(
  gradient$xy, gradient$labels, gradient$samples, workers = 1L
)
flat_3d <- refine_spatial_labels(
  cbind(gradient$xy, z = 0), gradient$labels,
  gradient$samples, workers = 1L
)
identical(flat_2d, flat_3d)
attr(flat_3d, "dimensions_used")

## ----genuine-3d---------------------------------------------------------------
volume_example <- simulate_volumetric_domains(
  n = 1800L,
  shape = "folded_layers",
  samples = 2L,
  noise = 0.18,
  seed = 23L
)
volume_refined <- refine_spatial_labels(
  volume_example$xy, volume_example$labels,
  volume_example$samples, workers = 2L
)
attr(volume_refined, "dimensions_used")

## ----specimen-isolation-------------------------------------------------------
source_rows <- which(gradient$samples == levels(gradient$samples)[1L])
overlapping_xy <- rbind(
  gradient$xy[source_rows, , drop = FALSE],
  gradient$xy[source_rows, , drop = FALSE]
)
overlapping_labels <- factor(
  rep(as.character(gradient$labels[source_rows]), 2L),
  levels = levels(gradient$labels)
)
overlapping_samples <- rep(
  c("section_1", "section_2"), each = length(source_rows)
)
joint <- refine_spatial_labels(
  overlapping_xy, overlapping_labels,
  samples = overlapping_samples, workers = 2L
)

separate <- unsplit(
  lapply(split(seq_len(nrow(overlapping_xy)), overlapping_samples), function(rows) {
    refine_spatial_labels(
      overlapping_xy[rows, , drop = FALSE],
      overlapping_labels[rows],
      workers = 1L
    )
  }),
  overlapping_samples
)

identical(as.character(joint), as.character(separate))

## ----worker-determinism, eval=FALSE-------------------------------------------
# one_worker <- refine_spatial_labels(
#   gradient$xy, gradient$labels, gradient$samples, workers = 1L
# )
# four_workers <- refine_spatial_labels(
#   gradient$xy, gradient$labels, gradient$samples, workers = 4L
# )
# attr(one_worker, "workers") <- NULL
# attr(four_workers, "workers") <- NULL
# identical(one_worker, four_workers)

## ----construct-mask-----------------------------------------------------------
mask_size <- 48L
mask_grid <- expand.grid(row = seq_len(mask_size), column = seq_len(mask_size))
radius <- sqrt(
  (mask_grid$row - (mask_size + 1) / 2)^2 +
    (mask_grid$column - (mask_size + 1) / 2)^2
)

reference_mask <- matrix("background", mask_size, mask_size)
reference_mask[radius < 17] <- "ring"
reference_mask[radius < 10] <- "core"
reference_mask[1:4, 1:4] <- NA_character_

mask_corruptions <- lapply(
  c("impulse", "boundary", "patch"),
  function(mechanism) {
    corrupt_categorical_mask(
      reference_mask, mechanism = mechanism,
      rate = 0.15, seed = 100L
    )
  }
)
names(mask_corruptions) <- c("impulse", "boundary", "patch")

vapply(
  mask_corruptions,
  function(x) mean(x != reference_mask, na.rm = TRUE),
  numeric(1L)
)

## ----clean-mask---------------------------------------------------------------
initial_mask <- mask_corruptions$boundary
cleaned_mask <- clean_categorical_mask(initial_mask, workers = 1L)

c(
  void_cells_preserved = all(is.na(cleaned_mask) == is.na(initial_mask)),
  changed_cells = sum(attr(cleaned_mask, "changed"), na.rm = TRUE)
)
dim(attr(cleaned_mask, "repair_margin"))

## ----mask-plot, fig.height=2.7------------------------------------------------
mask_colors <- field_palette(reference_mask)
old_par <- graphics::par(no.readonly = TRUE)
graphics::par(mfrow = c(1, 3), mar = c(1, 1, 2, 1))
plot_mask(reference_mask, "Reference", mask_colors)
plot_mask(initial_mask, "Boundary corruption", mask_colors)
plot_mask(cleaned_mask, "FiberMargin", mask_colors)
graphics::par(old_par)

## ----sampled-mask, eval=FALSE-------------------------------------------------
# mask_samples <- matrix(
#   ifelse(col(reference_mask) <= mask_size / 2, "left", "right"),
#   nrow = mask_size, ncol = mask_size
# )
# cleaned_by_section <- clean_categorical_mask(
#   initial_mask, samples = mask_samples, workers = 2L
# )

## ----voxel-mask, eval=FALSE---------------------------------------------------
# voxel_reference <- array("A", dim = c(20, 20, 8))
# voxel_reference[, 8:14, ] <- "B"
# voxel_reference[, 15:20, ] <- "C"
# voxel_initial <- corrupt_categorical_mask(
#   voxel_reference, mechanism = "boundary", rate = 0.10, seed = 9L
# )
# voxel_cleaned <- clean_categorical_mask(voxel_initial, workers = 2L)

## ----evaluate-mask------------------------------------------------------------
mask_metrics <- evaluate_mask_cleaning(
  reference = reference_mask,
  initial = initial_mask,
  cleaned = cleaned_mask,
  method = "FiberMargin"
)

mask_metrics[, c(
  "method", "initial_accuracy", "accuracy", "ari",
  "mean_iou", "mean_boundary_iou", "rare_class_iou",
  "correction_recall", "damage_rate"
)]

## ----repair-damage-identity---------------------------------------------------
with(
  mask_metrics,
  c(
    reported = accuracy,
    decomposed = initial_accuracy * (1 - damage_rate) +
      (1 - initial_accuracy) * correction_recall,
    numerical_error = repair_identity_error
  )
)

## ----spatial-benchmark--------------------------------------------------------
custom_benchmark <- spatial_benchmark(
  xy = gradient$xy,
  labels = gradient$labels,
  truth = gradient$truth,
  samples = gradient$samples,
  boundary = gradient$boundary,
  regions = gradient$area,
  sparse = gradient$sparse,
  name = "two_gradient_tissues",
  metadata = list(seed = 12L, purpose = "vignette")
)

class(custom_benchmark)
names(custom_benchmark)

## ----evaluate-spatial---------------------------------------------------------
manual_metrics <- evaluate_spatial_refinement(
  truth = custom_benchmark$truth,
  initial = custom_benchmark$labels,
  refined = refined_gradient,
  boundary = custom_benchmark$boundary,
  regions = custom_benchmark$regions,
  sparse = custom_benchmark$sparse,
  elapsed = 0,
  method = "FiberMargin"
)

manual_metrics[, c(
  "accuracy", "accuracy_gain", "ari", "macro_recall", "worst_recall",
  "boundary_accuracy", "sparse_region_accuracy",
  "correction_recall", "damage_rate", "changed_precision"
)]

## ----benchmark-refiners-------------------------------------------------------
jagged_small <- simulate_spatial_domains(
  n = 600L,
  pattern = "jagged_stripes",
  k = 4L,
  noise = 0.18,
  noise_type = "boundary",
  seed = 31L
)

comparison <- benchmark_spatial_refiners(
  data = list(
    gradient = custom_benchmark,
    jagged = jagged_small
  ),
  methods = list(
    FiberMargin = function(xy, labels, samples = NULL) {
      refine_spatial_labels(xy, labels, samples, workers = 1L)
    },
    IdentityExample = function(xy, labels, ...) labels
  ),
  include_initial = TRUE,
  seed = 5L,
  on_error = "stop"
)

comparison[, c(
  "dataset", "method", "accuracy", "ari",
  "correction_recall", "damage_rate", "seconds"
)]

## ----simulate-clusters--------------------------------------------------------
cluster_sim <- simulate_spatial_clusters(
  n = 600L,
  dimensions = 3L,
  k = 4L,
  samples = 2L,
  noise = 0.10,
  seed = 2L
)
c(dimensions = ncol(cluster_sim$xy), classes = nlevels(cluster_sim$truth))
table(cluster_sim$samples)

## ----simulate-domains---------------------------------------------------------
domain_sim <- simulate_spatial_domains(
  n = 800L,
  pattern = "thin_layers",
  k = 5L,
  noise = 0.20,
  dimensions = 2L,
  samples = 2L,
  noise_type = "boundary",
  feature_scale = 0.8,
  density_profile = "strong",
  seed = 3L
)
table(domain_sim$truth)

## ----simulate-complex---------------------------------------------------------
complex_sim <- simulate_complex_spatial_domains(
  n = 700L,
  shape = "tubular_network",
  density_profile = "extreme",
  noise_type = "patch",
  noise = 0.20,
  samples = 1L,
  k = 5L,
  seed = 4L
)
c(dimensions = ncol(complex_sim$xy), observations = nrow(complex_sim$xy))

## ----simulate-gradient-custom-------------------------------------------------
gradient_3d <- simulate_gradient_regions(
  n = 700L,
  minority = 0.05,
  dimensions = 3L,
  samples = 2L,
  seed = 5L,
  curvature = 0.06,
  density_profile = c(1, 2, 0.5, 3)
)
table(gradient_3d$area, gradient_3d$truth)

## ----simulate-volume----------------------------------------------------------
volume_sim <- simulate_volumetric_domains(
  n = 700L,
  shape = "folded_layers",
  acquisition = "class_imbalanced",
  noise_type = "boundary",
  noise = 0.20,
  k = 5L,
  samples = 2L,
  seed = 6L
)

c(
  dimensions = ncol(volume_sim$xy),
  observations = nrow(volume_sim$xy),
  samples = nlevels(volume_sim$samples)
)
volume_sim$region_counts

## ----simulator-gallery, fig.height=5.8----------------------------------------
old_par <- graphics::par(no.readonly = TRUE)
graphics::par(mfrow = c(2, 2), mar = c(3, 3, 2, 1))
plot_field(cluster_sim$xy, cluster_sim$truth, "Spatial clusters", cex = 0.45)
plot_field(domain_sim$xy, domain_sim$truth, "Thin layers", cex = 0.45)
plot_field(complex_sim$xy, complex_sim$truth, "Tubular network", cex = 0.45)
plot_field(volume_sim$xy, volume_sim$truth, "3D folded layers: x-y", cex = 0.45)
graphics::par(old_par)

## ----available-benchmarks-----------------------------------------------------
benchmark_status <- available_spatial_benchmarks()
benchmark_status[c(
  "dataset", "included", "observations", "classes", "scenarios"
)]

## ----load-dlpfc---------------------------------------------------------------
dlpfc_case <- load_spatial_benchmark("dlpfc", scenario = 1L)
c(
  observations = nrow(dlpfc_case$xy),
  classes = nlevels(dlpfc_case$truth),
  specimens = nlevels(dlpfc_case$samples),
  input_accuracy = mean(dlpfc_case$labels == dlpfc_case$truth)
)
dlpfc_case$scenario

## ----load-crc, eval=FALSE-----------------------------------------------------
# crc_case <- load_spatial_benchmark(
#   "crc", scenario = "CRC_random_25_r1", seed = 1040001L
# )
# crc_refined <- refine_spatial_labels(
#   crc_case$xy, crc_case$labels,
#   samples = crc_case$samples,
#   workers = 4L
# )
# evaluate_spatial_refinement(
#   crc_case$truth, crc_case$labels, crc_refined,
#   boundary = crc_case$boundary,
#   sparse = crc_case$sparse,
#   method = "FiberMargin"
# )

## ----end-to-end-template, eval=FALSE------------------------------------------
# # One row per spatial location.
# xy <- as.matrix(my_data[, c("x", "y")])
# labels <- factor(my_data$annotation)
# samples <- factor(my_data$section_id)
# 
# # Repair does not receive the reference labels.
# refined <- refine_spatial_labels(
#   xy = xy,
#   labels = labels,
#   samples = samples,
#   workers = 4L
# )
# 
# # Inspect edits before replacing an annotation column.
# audit <- data.frame(
#   observed = labels,
#   candidate = attr(refined, "candidate"),
#   refined = refined,
#   changed = attr(refined, "changed"),
#   margin_score = attr(refined, "margin_score"),
#   required = attr(refined, "required"),
#   repair_margin = attr(refined, "repair_margin"),
#   isolation = attr(refined, "isolation")
# )
# 
# my_data$fibermargin_annotation <- refined

## ----end-to-end-evaluation, eval=FALSE----------------------------------------
# evaluate_spatial_refinement(
#   truth = reference_labels,
#   initial = labels,
#   refined = refined,
#   boundary = boundary_locations,
#   sparse = rare_or_sparse_locations,
#   method = "FiberMargin"
# )