Introduction
The creation of soil profile
sketches has been a cornerstone of the {aqp} package since the first
CRAN release in 2010. With the addition of horizon depth annotation in
2017 it became clear that some form of label placement adjustment would
be required to avoid overlapping annotation. Thin horizons and / or
large font sizes almost always lead to overlapping horizon depth labels
(see numbers just to the right of horizon boundaries in the figure
below).
Version 2.0 of the {aqp} package provides two algorithms for avoiding
overlapping annotation, available in plotSPC() for the
adjustment of horizon depth annotation, and for general use via
fixOverlap(). Within the context of plotSPC(),
overlapping horizon depth labels are flagged according to an overlap
threshold computed from label text height on the current graphics
device.
Consider the following example soil profile sketches, created from a
single soil profile template characterized by multiple thin horizons.
Simulation and duplication were used to generate multiple copies of the
three variations on the original template.
library(aqp)
# soil profile template
x <- quickSPC(
"SPC:Oi|Oe|AAAAA|E1|E2|Bhs1Bhs1Bhs1|Bhs2|CCCCCC|RRRRRRRRRR",
interval = 3
)
# set variability in horizon thickness
horizons(x)$thick.sd <- 3
# setup horizon designations
horizons(x)$nm <- factor(x$name, levels = c('Oi', 'Oe', 'A', 'E1', 'E2', 'Bhs1', 'Bhs2', 'C', 'R'))
# simulate 3 realizations of original profile
set.seed(10101)
s <- perturb(x, n = 3, thickness.attr = 'thick.sd', min.thickness = 2)
# create 3 copies of each simulation
s <- duplicate(s, times = 3)
Overlapping horizon depth labels are clear. Note that there are three
copies of each simulated variation on the original template.
Selecting the optimal label-adjustment methods and parameters will
depend on the collection of profiles, horizonation, profile height, font
sizes, and graphics device settings (figure height, resolution).
electrostatic simulation: This is generally the best
method, having deterministic solutions that typically move labels as
little as possible in fewer iterations. It does not do well with many
very thin horizons near the top or bottom of the profile. Very complex
profiles and larger font sizes may require setting q to
values larger than 1. Set q to values between 0.25-0.75 to
minimize adjustment distances. This method cannot (yet) provide a label
adjustment solution in the presence of missing horizon depths or
horizons where top == bottom.
simulated annealing: This method is the most adaptive,
typically capable of providing reasonable solutions for even the most
complex profiles. The non-deterministic solutions can be undesirable
when preparing figures for publication.
Electrostatic Simulation
An electrostatic simulation applies forces of repulsion between
labels that are within a given distance threshold and forces of
attraction to a uniformly spaced sequence. Affected labels are
iteratively perturbed until either no overlap is reported, or a maximum
number of iterations has been reached. Label adjustment solutions are
deterministic and can be controlled by additional arguments to
electroStatics_1D(). The most common adjustment is via
charge density (q). See comparisons below.
Simulated Annealing
This approach makes small adjustments to affected label positions
until overlap is removed, or until a maximum number of iterations is
reached. Rank order and boundary conditions are preserved. The
underlying algorithm is based on simulated annealing. Label placement
solutions are non-deterministic, use set.seed() if
repeatability is important. The “cooling schedule” parameters
T0 and k can be used to tune the algorithm for
specific applications. These additional arguments are passed to
SANN_1D() via fixOverlap().
Note that each solution, within duplicate profiles, is slightly
different.
Quick Comparison
Create a very busy profile with lots of possible overlapping horizon
depth labels.
x <- quickSPC(
"SPC:Oi|Oe|AAA|E1|E2|E3|BhsBhsBhsBhs|Bt1|Bt2|Bt3Bt3|CCCCCC|Ab1|Ab2|2C2C2C2C2C2C|2Cr|2R2R2R2R2R2R2R2R",
interval = 1
)
x$z <- as.numeric(x$hzID)
Set arguments to plotSPC() and define a custom function
to demonstrate various label adjustment settings.
# pretty colors
.bluecolors <- hcl.colors(n = 25, palette = 'Blues')[-25]
# plotSPC arguments
.a <- list(
width = 0.2,
max.depth = 40,
hz.depths = TRUE,
name.style = 'center-center',
cex.names = 1.5,
name = NA,
depth.axis = FALSE,
color = 'z',
show.legend = FALSE,
print.id = FALSE,
col.palette = .bluecolors
)
# set plotSPC default arguments
options(.aqp.plotSPC.args = .a)
# wrapper function to test label collision solutions
testIt <- function(x, ...) {
# make sketches
plotSPC(x, ...)
# a normalized index of label adjustment
.LAI <- get('last_spc_plot', envir = aqp.env)$hz.depth.LAI
.LAI <- ifelse(is.na(.LAI), 0, .LAI)
# annotate with label adjustment index
.txt <- sprintf("LAI: %0.3f", .LAI)
mtext(.txt, side = 1, at = 1, line = -1.5, cex = 0.8)
}
Compare and contrast. The “LAI” is an index of label adjustment,
relative to the original label positions. Larger values imply a higher
degree of adjustment. Note that each simulated annealing (SANN) label
adjustment is different, unless randomness is “controlled” using
set.seed().
par(mar = c(1, 0, 0, 0), mfcol = c(1, 8))
testIt(x, fixLabelCollisions = FALSE)
title('No\nAdjustment', line = -3.5, adj = 0.5)
testIt(x, fixOverlapArgs = list(method = 'S'))
title('SANN\nsim 1', line = -3.5, adj = 0.5)
testIt(x, fixOverlapArgs = list(method = 'S'))
title('SANN\nsim 2', line = -3.5, adj = 0.5)
testIt(x, fixOverlapArgs = list(method = 'S'))
title('SANN\nsim 3', line = -3.5, adj = 0.5)
testIt(x, fixOverlapArgs = list(method = 'E', q = 1.5))
title('Electrostatic\nq = 1.5', line = -3.5, adj = 0.5)
testIt(x, fixOverlapArgs = list(method = 'E', q = 1))
title('Electrostatic\nq = 1', line = -3.5, adj = 0.5)
testIt(x, fixOverlapArgs = list(method = 'E', q = 0.5))
title('Electrostatic\nq = 0.5', line = -3.5, adj = 0.5)
testIt(x, fixOverlapArgs = list(method = 'E', q = 0.25))
title('Electrostatic\nq = 0.25', line = -3.5, adj = 0.5)
General Cases
The label adjustment functionality can be used outside of the context
of plotSPC(). The function overlapMetics() can
be used to determine overlap within a vector of positions based on a
given distance threshold.
x <- c(1, 2, 3, 3.4, 3.5, 5, 6, 10)
overlapMetrics(x, thresh = 0.5)
#> $idx
#> [1] 3 4 5
#>
#> $ov
#> [1] 0.125
The fixOverlap() function will attempt adjustment, given
distance threshold and a method selection. Additional parameters are
passed on to SANN_1D() or electroStatics_1D().
The converged attribute is included in the result to
signify that a successful solution was possible. See the manual page for
fixOverlap() for additional details.
# vector of positions, typically labels but could be profile sketch alignment on the x-axis
s <- c(1, 2, 2.3, 4, 5, 5.5, 7)
# simulated annealing, solution is non-deterministic
fixOverlap(s, thresh = 0.5, method = 'S')
#> 4 iterations
#> [1] 1.000000 2.043921 2.665304 4.000000 5.000000 5.500000 7.000000
#> attr(,"converged")
#> [1] TRUE
# electrostatics-inspired simulation of particles
# solution is deterministic
fixOverlap(s, thresh = 0.5, method = 'E')
#> 2 iterations
#> [1] 1.000000 1.847048 2.419053 4.000000 5.000000 5.500000 7.000000
#> attr(,"converged")
#> [1] TRUE
Compare Algorithms
Define a custom function for comparing the extended output from
fixOverlap(..., trace = TRUE) using both electrostatic
simulation (method = 'E') and simulated annealing
(method = 'S').
evalMethods <- function(x, thresh, q, ...) {
cols <- hcl.colors(n = 9, palette = 'Zissou 1', rev = TRUE)
cols <- colorRampPalette(cols)(length(x))
z <- fixOverlap(x, thresh = thresh, method = 'E', maxIter = 100, trace = TRUE, q = q)
.n <- nrow(z$states)
op <- par(mar = c(0, 2, 2, 0.5), bg = 'black', fg = 'white')
layout(matrix(c(1, 2, 3, 4), ncol = 2, nrow = 2), heights = c(0.33, 0.66))
plot(seq_along(z$cost), z$cost, las = 1, type = 'b', axes = FALSE, cex = 0.66, xlim = c(1, .n))
mtext(text = sprintf("Converged (%s): %s", .n, z$converged), at = 0, side = 3, line = 0, cex = 0.75, font = 3, adj = 0)
matplot(rbind(x, z$states), type = 'l', lty = 1, las = 1, axes = FALSE, col = cols, lwd = 1)
points(x = rep(1, times = length(x)), y = x, cex = 0.66, pch = 16, col = cols)
points(x = rep(.n + 1, times = length(x)), y = z$x, cex = 0.66, pch = 16, col = cols)
text(x = 1, y = x, col = cols, labels = seq_along(x), cex = 0.66, font = 2, pos = 2)
text(x = .n + 1, y = z$x, col = cols, labels = seq_along(x), cex = 0.66, font = 2, pos = 4)
axis(side = 2, at = unique(x), labels = round(unique(x), 1), col.axis = par('fg'), las = 1, cex.axis = 0.6)
title(main = 'Electrostatic Simulation', line = 1, col.main = 'white')
## SANN_1D doesn't always preserve rank ordering
## ->> not designed to use unsorted input
## ->> maybe impossible with ties in x?
z <- fixOverlap(x, thresh = thresh, method = 'S', trace = TRUE, maxIter = 1000)
.n <- nrow(z$states)
plot(seq_along(z$stats), z$stats, las = 1, type = 'b', axes = FALSE, cex = 0.66, xlim = c(1, .n))
mtext(text = sprintf("Converged (%s): %s", .n, z$converged), at = 0, side = 3, line = 0, cex = 0.75, font = 3, adj = 0)
matplot(z$states, type = 'l', lty = 1, las = 1, axes = FALSE, col = cols)
points(x = rep(1, times = length(x)), y = z$states[1, ], cex = 0.66, pch = 16, col = cols)
points(x = rep(.n, times = length(x)), y = z$x, cex = 0.66, pch = 16, col = cols)
text(x = 1, y = z$states[1, ], col = cols, labels = seq_along(x), cex = 0.66, font = 2, pos = 2)
text(x = .n, y = z$x, col = cols, labels = seq_along(x), cex = 0.66, font = 2, pos = 4)
axis(side = 2, at = unique(x), labels = round(unique(x), 1), col.axis = par('fg'), las = 1, cex.axis = 0.6)
title(main = 'Simulated Annealing', line = 1, col.main = 'white')
# reset graphics state
par(op)
layout(1)
}
The top panels represent the objective function, smaller values are
better. Bottom panels illustrate perturbations to elements of the
position vector x. Larger values of charge density
q are required applying the electrostatic simulation method
to increasingly complex problems (e.g. larger thresholds, more overlap).
However, setting q too high will result in chaos and
failure to converge.
# explore effect of charge (q)
# too large -> chaos
x <- c(0, 2, 5, 12, 18, 20, 35, 40, 50, 56, 90, 120, 145, 150)
# just about right, very few perturbations required
evalMethods(x, thresh = 5, q = 1.1)
# ok, but now most label positions are affected
evalMethods(x, thresh = 5, q = 1.8)
# too high, wasting time on more iterations
evalMethods(x, thresh = 5, q = 3)
# far too high, wasting more time with little gain
evalMethods(x, thresh = 5, q = 4)
# chaos and failure to converge
evalMethods(x, thresh = 5, q = 5)
Additional examples to tinker with.
# threshold too large
evalMethods(x, thresh = 10, q = 3)
# large threshold
x <- c(0, 5, 12, 18, 20, 35, 40, 55, 90, 120, 145, 150)
evalMethods(x, thresh = 9, q = 2)
# single iteration enough
x <- c(0, 3, 20, 35, 40, 55, 90, 120, 145, 150)
evalMethods(x, thresh = 6, q = 1)
# clusters
x <- sort(c(0, jitter(rep(10, 3)), jitter(rep(25, 3)), jitter(rep(90, 3)), 150))
evalMethods(x, thresh = 6, q = 3)
evalMethods(x, thresh = 6, q = 2)
## impact of scale / offset
x <- c(0, 5, 12, 18, 20, 35, 40, 50, 120, 145, 150)
# works as expected
evalMethods(x, thresh = 5, q = 1.1)
# works as expected, as long as threshold is scaled
evalMethods(x / 10, thresh = 5 / 10, q = 1.1)
# works as expected, as long as threshold is scaled
evalMethods(x * 10, thresh = 5 * 10, q = 1.1)
# all work as expected, threshold not modified
evalMethods(x + 10, thresh = 5, q = 1.1)
evalMethods(x + 100, thresh = 5, q = 1.1)
evalMethods(x + 1000, thresh = 5, q = 1.1)
# works as expected
x <- c(315, 325, 341, 353, 366, 374, 422)
fixOverlap(x, thresh = 9.7, q = 1, method = 'E')
evalMethods(x, thresh = 9.7, q = 1)
x <- c(1.0075, 1.1200, 1.3450, 1.6450, 1.8700, 1.8825)
fixOverlap(x, thresh = 0.05442329, q = 1)
evalMethods(x, thresh = 0.05442329, q = 1)
Learning More about the SANN Algorithm
Define a custom function for visualizing the extended output from
fixOverlap(..., method = 'S', trace = TRUE).
tracePlot <- function(x, z, cex.axis.labels = 0.85) {
# setup plot device
op <- par(mar = c(4, 4, 1, 1), bg = 'black', fg = 'white')
layout(matrix(c(1,2,3)), widths = 1, heights = c(1,1,2))
# order:
# B: boundary condition violation
# O: rank (order) violation
# +: accepted perturbation
# -: rejected perturbation
cols <- c(grey(0.5), grey(0.85), 'royalblue', 'firebrick')
cols.lines <- hcl.colors(9, 'Zissou 1', rev = TRUE)
cols.lines <- colorRampPalette(cols.lines)(length(x))
# total overlap (objective function) progress
plot(
seq_along(z$stats), z$stats,
type = 'h', las = 1,
xlab = 'Iteration', ylab = 'Total Overlap',
axes = FALSE,
col = cols[as.numeric(z$log)]
)
axis(side = 2, cex.axis = cex.axis.labels, col.axis = 'white', las = 1, line = -2)
mtext('Overlap', side = 2, line = 2, cex = cex.axis.labels, font = 2)
# deviation from original configuration
plot(
seq_along(z$stats), z$ssd,
type = 'h', las = 1,
xlab = 'Iteration', ylab = 'Deviation',
axes = FALSE,
col = cols[as.numeric(z$log)]
)
axis(side = 2, cex.axis = cex.axis.labels, col.axis = 'white', las = 1, line = -2)
mtext('Deviation', side = 2, line = 2, cex = cex.axis.labels, font = 2)
legend('top', legend = c('boundary\nviolation', 'rank\nviolation', 'accepted\nperturbation', 'rejected\nperturbation'), col = cols, bty = 'n', horiz = TRUE, inset = -0.5, lty = 1, lwd = 2, xpd = NA)
# adjustments at each iteration
matplot(
z$states, type = 'l',
lty = 1, las = 1,
xlab = 'Iteration', ylab = 'x-position',
axes = FALSE,
col = cols.lines
)
axis(side = 2, cex.axis = cex.axis.labels, col.axis = 'white', las = 1, at = x, labels = round(x, 1))
axis(side = 4, cex.axis = cex.axis.labels, col.axis = 'white', las = 1, at = z$x, labels = round(z$x, 1), line = -2)
mtext('Position', side = 2, line = 2.5, cex = cex.axis.labels, font = 2)
axis(side = 1, cex.axis = 1, col.axis = 'white', line = 0)
mtext('Iteration', side = 1, line = 2.5, cex = cex.axis.labels, font = 2)
par(op)
layout(1)
}
A relatively challenging example.
x <- c(0, 1, 2, 2.2, 2.8, 3.5, 6, 8, 10, 10.1, 12.8, 13, 14.8, 15, 15.5)
# fix overlap, return debugging information
set.seed(10101)
z <- fixOverlap(x, thresh = 0.73, method = 'S', trace = TRUE)
#> 160 iterations
# check convergence
z$converged
#> [1] TRUE
# inspect algorithm trace
tracePlot(x, z)
# trace log
# B: boundary condition violation
# O: rank (order) violation
# +: accepted perturbation
# -: rejected perturbation
table(z$log)
#>
#> B O + -
#> 22 71 64 2
A very challenging example.
# fix overlap, return debugging information
set.seed(101010)
x <- sort(runif(10, min = 2.5, max = 3.5))
# widen boundary conditions
z <- fixOverlap(x, thresh = 0.2, trace = TRUE, min.x = 0, max.x = 10, maxIter = 2000, adj = 0.05)
#> 1203 iterations
# check convergence
z$converged
#> [1] TRUE
# inspect algorithm trace
tracePlot(x, z)
Cleanup.
# reset plotSPC() options
options(.aqp.plotSPC.args = NULL)
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.