## -----------------------------------------------------------------------------
#| label: setup
#| include: false
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library("Coreset")


## -----------------------------------------------------------------------------
#| label: install
#| eval: false
# install.packages("Coreset")


## -----------------------------------------------------------------------------
#| label: quickstart
# Load the `eurodist` dist object
data("eurodist")

# Set a seed for a reproducible selection
set.seed(1)

# Select 4 maximally dispersed cities
ffPick <- FarFirst(4L, eurodist)

# View distances between selected cities
as.matrix(eurodist)[ffPick, ffPick]

# Quickly extract the minimum distance for a given selection
MinDist(eurodist, ffPick)


## -----------------------------------------------------------------------------
#| label: gonzalez-random
set.seed(1)
FarFirst(6L, eurodist)   # default: best of three random starts


## -----------------------------------------------------------------------------
#| label: gonzalez-ensemble
set.seed(1)
# Fewer starts provides a faster run, but the solution found may be inferior
FarFirst(6L, eurodist, nSeeds = 2L)


## -----------------------------------------------------------------------------
#| label: gonzalez-seeds
ffPick <- FarFirst(6L, eurodist, strategy = c("diameter", "anti_medoid"))
MinDist(eurodist, ffPick)

attr(ffPick, "winning_strategy")


## -----------------------------------------------------------------------------
#| label: gonzalez-oracle
data("USArrests")
arrestTypes <- USArrests[, c("Murder", "Assault", "Rape")]
StateDist <- function(i) {
  diffs <- sweep(arrestTypes, 2, unlist(arrestTypes[i, ]), "-")
  sqrt(rowSums(diffs ^ 2))
}
idx <- FarFirst(4L, StateDist, N = nrow(arrestTypes), strategy = 1L)
arrestTypes[idx, ]


## -----------------------------------------------------------------------------
#| label: dropadd
daPick <- DropAdd(6L, eurodist, plateau = 500L)

daPick
labels(eurodist)[daPick]


## -----------------------------------------------------------------------------
#| label: grasp
set.seed(0)
grPick <- Grasp(6L, eurodist, plateau = 50L)
grPick
labels(eurodist)[grPick]
attr(grPick, "pr_calls")   # path-relinking calls performed


## -----------------------------------------------------------------------------
#| label: sim-data
set.seed(1) # Seed selected such that FarFirst < DropAdd < Grasp
pts <- matrix(rnorm(100), ncol = 2)   # 50 points, 2 dimensions
d50 <- dist(pts)
k   <- 8L


## -----------------------------------------------------------------------------
#| label: compare-run
set.seed(1) # Seed selected such that FarFirst < DropAdd < Grasp
ffPick <- FarFirst(k, d50)
da50Pick <- DropAdd(k, d50, plateau = 500L)
gr50Pick <- Grasp(k, d50, plateau = 50L)


## -----------------------------------------------------------------------------
#| label: compare-scores
scores <- c(
  FarFirst  = attr(ffPick, "score"),
  DropAdd = attr(da50Pick, "score"),
  Grasp   = attr(gr50Pick, "score")
)
round(scores, 3)


## -----------------------------------------------------------------------------
#| label: compare-plot
#| fig-width: 6.5
#| fig-height: 6
#| fig-cap: "Selections returned by each method on 50 random 2-D points
#|   (k = 8). Coloured symbols mark selected points; grey circles are the
#|   full candidate set."

methods <- list(
  FarFirst  = ffPick,
  DropAdd = da50Pick,
  Grasp   = gr50Pick
)
cols <- c(FarFirst = "#E41A1C", DropAdd = "#377EB8", Grasp = "#4DAF4A")
pchs <- c(FarFirst = 1L, DropAdd = 3L, Grasp = 4L)

plot(pts, pch = 1L, col = "grey75", asp = 1L,
     xlab = "x", ylab = "y", frame.plot = FALSE,
     main = "Coreset method comparison")

for (nm in names(methods)) {
  sel <- methods[[nm]]
  points(pts[sel, 1L], pts[sel, 2L], pch = pchs[nm], col = cols[nm],
         cex = 1.6)
}

legend_labels <- lapply(seq_along(methods), function(i) {
  bquote(.(names(methods)[i]) ~ (T[k] == .(sprintf("%.3f", scores[i]))))
})

legend("topleft", legend = as.expression(legend_labels), pch = pchs, col = cols,
       pt.bg = cols, pt.cex = 1.4, bty = "n")


## -----------------------------------------------------------------------------
#| label: install-highs
#| eval: false
# install.packages("highs")


## -----------------------------------------------------------------------------
#| label: exact-data
set.seed(1L)
pts30 <- matrix(rnorm(60L), ncol = 2L)
d30   <- dist(pts30)


## -----------------------------------------------------------------------------
#| label: exact
#| eval: !expr requireNamespace("highs", quietly = TRUE)
# exPick <- ExactMaxMin(6L, d30, maxSeconds = 30L)
# 
# attr(exPick, "proven")      # TRUE  ⟹  objective is the global optimum
# attr(exPick, "score")
# 
# # Compare to the greedy heuristic on the same instance
# ff30Pick <- FarFirst(6L, d30)
# c(exact    = attr(exPick, "score"),
#   farFirst = MinDist(d30, ff30Pick))


## -----------------------------------------------------------------------------
#| label: MinDist
MinDist(d50, ffPick)                               # from dist
MinDist(as.matrix(d50), ffPick)                    # from square matrix
MinDist(points = pts, idx = ffPick)                # from coordinates


## -----------------------------------------------------------------------------
#| label: exact-maxsum
#| eval: !expr requireNamespace("highs", quietly = TRUE)
# smPick <- ExactMaxSum(6L, d30, maxSeconds = 30L)
# 
# attr(smPick, "proven")      # TRUE  ⟹  objective is the global optimum
# attr(smPick, "score")       # total pairwise distance within the selection


## -----------------------------------------------------------------------------
#| label: maxmean-data
set.seed(1)
affinity <- matrix(runif(30L * 30L, min = -10, max = 10), nrow = 30L)
affinity <- (affinity + t(affinity)) / 2   # symmetric
diag(affinity) <- 0


## -----------------------------------------------------------------------------
#| label: maxmean
set.seed(1)
mmPick <- MaxMean(affinity, maxSeconds = 2)
mmPick
attr(mmPick, "size")              # the algorithm chose this subset size
attr(mmPick, "score")             # achieved mean-dispersion objective


## -----------------------------------------------------------------------------
#| label: meandist
MeanDist(affinity, mmPick)        # matches attr(mmPick, "score")
MeanDist(affinity, 1:30)          # the full set scores lower


## -----------------------------------------------------------------------------
#| label: kcentre
centres <- KCentre(4L, eurodist)
labels(eurodist)[centres]
centres


## -----------------------------------------------------------------------------
#| label: kcentre-radius
ff <- FarFirst(4L, eurodist, strategy = "peripheral")
c(KCentre  = KCentreRadius(eurodist, centres),
  FarFirst = KCentreRadius(eurodist, ff))


## -----------------------------------------------------------------------------
#| label: exact-kcentre
#| eval: !expr requireNamespace("highs", quietly = TRUE)
# kc <- ExactKCentre(4L, eurodist)
# kc
# attr(kc, "proven")      # TRUE  ⟹  radius is the global covering optimum


## -----------------------------------------------------------------------------
#| label: dispersion-vs-covering
#| eval: !expr requireNamespace("highs", quietly = TRUE)
# disp <- ExactMaxMin(4L, eurodist)
# labels(eurodist)[disp]  # dispersion: pushed to the extremes
# labels(eurodist)[kc]    # covering: pulled toward the interior


## -----------------------------------------------------------------------------
#| label: maxentropy
mePick <- MaxEntropy(4L, eurodist)
labels(eurodist)[mePick]
attr(mePick, "score")     # the achieved log-determinant
attr(mePick, "exact")      # TRUE if certified by exact enumeration


## -----------------------------------------------------------------------------
#| label: maxentropy-negmass
attr(mePick, "negMass")