% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/perturb.R
\name{perturb}
\alias{perturb}
\alias{permute_profile}
\title{Perturb soil horizon depths using boundary distinctness}
\usage{
perturb(
  p,
  n = 100,
  id = NULL,
  thickness.attr = NULL,
  boundary.attr = NULL,
  min.thickness = 1,
  max.depth = NULL,
  new.idname = "pID"
)
}
\arguments{
\item{p}{A SoilProfileCollection}

\item{n}{Number of new profiles to generate (default: \code{100}) per profile in \code{p}}

\item{id}{a vector of profile IDs with length equal to (\code{n}). Overrides use of \code{seq_len(n)} as default profile ID values.}

\item{thickness.attr}{Horizon variance attribute containing numeric "standard deviations" reflecting horizon thickness}

\item{boundary.attr}{Horizon variance attribute containing numeric "standard deviations" reflecting boundary transition distinctness}

\item{min.thickness}{Minimum thickness of permuted horizons (default: \code{1})}

\item{max.depth}{Depth below which horizon depths are not perturbed (default: \code{NULL})}

\item{new.idname}{New column name to contain unique profile ID (default: \code{pID})}
}
\value{
a SoilProfileCollection with \code{n} realizations of each profile in \code{p}
}
\description{
"Perturbs" the \strong{boundary between horizons} or the \strong{thickness of horizons} using a standard deviation specified as a horizon-level attribute. This is selected using either \code{boundary.attr} or \code{thickness.attr} to specify the column name.

The boundary standard deviation corresponds roughly to the concept of "horizon boundary distinctness." In contrast, the \emph{horizon thickness} standard deviation corresponds roughly to the "variation in horizon thickness" so it may be determined from several similar profiles that have a particular layer "in common."
}
\details{
Imagine a Normal curve with mean centered on the vertical (depth axis) at a representative value (RV) horizon bottom depth or thickness. By the Empirical Rule for Normal distribution, two "standard deviations" above or below that "central" mean value represent 95\% of the "typical volume" of that horizon or boundary.

\code{perturb()} can leverage semi-quantitative (ordered factor) levels of boundary distinctness/topography for the upper and lower boundary of individual horizons. A handy function for this is \code{\link[=hzDistinctnessCodeToOffset]{hzDistinctnessCodeToOffset()}}. The \code{boundary.attr} is arguably easier to parameterize from a single profile description or "Form 232" where \emph{horizon boundary distinctness} classes (based on vertical distance of transition) are conventionally recorded for each layer.

Alternately, \code{perturb()} can be parameterized using standard deviation in thickness of layers derived from a group. Say, the variance parameters are defined from a set of pedons correlated to a particular series or component, and the template "seed" profile is, for example, the Official Series Description or the Representative Component Pedon.
}
\examples{

### THICKNESS

# load sample data and convert into SoilProfileCollection
data(sp3)
depths(sp3) <- id ~ top + bottom

# select a profile to use as the basis for simulation
s <- sp3[3,]

# reset horizon names
s$name <- paste('H', seq_along(s$name), sep = '')

# simulate 25 new profiles
horizons(s)$hz.sd <- 2 # constant standard deviation
sim.1 <- perturb(s, n = 25, thickness.attr = "hz.sd")

# simulate 25 new profiles using different SD for each horizon
horizons(s)$hz.sd <- c(1, 2, 5, 5, 5, 10, 3)
sim.2 <- perturb(s, n = 25, thickness.attr = "hz.sd")

# plot
par(mfrow = c(2, 1), mar = c(0, 0, 0, 0))
plot(sim.1)
mtext(
  'SD = 2',
  side = 2,
  line = -1.5,
  font = 2,
  cex = 0.75
)
plot(sim.2)
mtext(
  'SD = c(1, 2, 5, 5, 5, 10, 3)',
  side = 2,
  line = -1.5,
  font = 2,
  cex = 0.75
)

# aggregate horizonation of simulated data
# note: set class_prob_mode=2 as profiles were not defined to a constant depth
sim.2$name <- factor(sim.2$name)
a <- slab(sim.2, ~ name, cpm=2)

# convert to long format for plotting simplicity
library(data.table)
a.long <- data.table::melt(data.table::as.data.table(a),
               id.vars = c('top', 'bottom'),
                measure.vars = levels(sim.2$name))

# plot horizon probabilities derived from simulated data
# dashed lines are the original horizon boundaries
library(lattice)

xyplot(
  top ~ value,
  groups = variable,
  data = a.long,
  subset = value > 0,
  ylim = c(100,-5),
  type = c('l', 'g'),
  asp = 1.5,
  ylab = 'Depth (cm)',
  xlab = 'Probability',
  auto.key = list(
    columns = 4,
    lines = TRUE,
    points = FALSE
  ),
  panel = function(...) {
    panel.xyplot(...)
    panel.abline(h = s$top, lty = 2, lwd = 2)
  }
)

### BOUNDARIES

# example with sp1 (using boundary distinctness)
data("sp1")
depths(sp1) <- id ~ top + bottom

# specify "standard deviation" for boundary thickness
#   consider a normal curve centered at boundary RV depth
# lookup table: ~maximum thickness of boundary distinctness classes, divided by 3
bound.lut <- c('V'=0.5,'A'=2,'C'=5,'G'=15,'D'=45) / 3

## V          A          C          G          D
## 0.1666667  0.6666667  1.6666667  5.0000000 15.0000000

sp1$bound_sd <- bound.lut[sp1$bound_distinct]

# hold any NA boundary distinctness constant
sp1$bound_sd[is.na(sp1$bound_sd)] <- 0

quantile(sp1$bound_sd, na.rm = TRUE)
p <- sp1[3]

# assume boundary sd is 1/12 midpoint of horizon depth
# (i.e. general relationship: SD increases (less well known) with depth)
sp1 <- transform(sp1, midpt = (bottom - top) / 2 + top, bound_sd = midpt / 12)
quantile(sp1$bound_sd)

perturb(p, boundary.attr = "bound_sd", n = 10)


### Custom IDs

ids <- sprintf("\%s-\%03d", profile_id(p), 1:10) 
perturb(p, boundary.attr = "bound_sd", id = ids)


}
\seealso{
\code{\link[=random_profile]{random_profile()}} \code{\link[=hzDistinctnessCodeToOffset]{hzDistinctnessCodeToOffset()}}
}
\author{
D.E. Beaudette, A.G. Brown
}