R/compute.bagplot.r In rainbow: Bagplots, Boxplots and Rainbow Plots for Functional Data

Defines functions compute.bagplot

```compute.bagplot = function (x, y, factor = 3, na.rm = FALSE, approx.limit = 300,
dkmethod = 2, precision = 1, verbose = FALSE, debug.plots = "no")
{
win <- function(dx, dy) {
atan2(y = dy, x = dx)
}
out.of.polygon <- function(xy, pg) {
if (nrow(pg) == 1)
return(xy[, 1] == pg[1] & xy[, 2] == pg[2])
m <- nrow(xy <- matrix(xy, ncol = 2))
n <- nrow(pg)
limit <- -abs(1e-10 * diff(range(pg)))
pgn <- cbind(diff(c(pg[, 2], pg[1, 2])), -diff(c(pg[,
1], pg[1, 1])))
S <- matrix(colMeans(xy), 1, 2)
dxy <- cbind(S[1] - pg[, 1], S[2] - pg[, 2])
S.in.pg <- all(limit < apply(dxy * pgn, 1, sum))
if (!all(limit < apply(dxy * pgn, 1, sum))) {
pg <- pg[n:1, ]
pgn <- -pgn[n:1, ]
}
in.pg <- rep(TRUE, m)
for (j in 1:n) {
dxy <- xy - matrix(pg[j, ], m, 2, byrow = TRUE)
in.pg <- in.pg & limit < (dxy %*% pgn[j, ])
}
return(!in.pg)
}
cut.z.pg <- function(zx, zy, p1x, p1y, p2x, p2y) {
a2 <- (p2y - p1y) / (p2x - p1x)
a1 <- zy/zx
sx <- (p1y - a2 * p1x) / (a1 - a2)
sy <- a1 * sx
sxy <- cbind(sx, sy)
h <- any(is.nan(sxy)) || any(is.na(sxy)) || any(Inf ==
abs(sxy))
if (h) {
if (!exists("verbose"))
verbose <- FALSE
if (verbose)
cat("special")
h <- 0 == (a1 - a2) & sign(zx) == sign(p1x)
sx <- ifelse(h, p1x, sx)
sy <- ifelse(h, p1y, sy)
h <- 0 == (a1 - a2) & sign(zx) != sign(p1x)
sx <- ifelse(h, p2x, sx)
sy <- ifelse(h, p2y, sy)
h <- p1x == p2x & zx != p1x & p1x != 0
sx <- ifelse(h, p1x, sx)
sy <- ifelse(h, zy * p1x / zx, sy)
h <- p1x == p2x & zx != p1x & p1x == 0
sx <- ifelse(h, p1x, sx)
sy <- ifelse(h, 0, sy)
h <- p1x == p2x & zx == p1x & p1x == 0 & sign(zy) ==
sign(p1y)
sx <- ifelse(h, p1x, sx)
sy <- ifelse(h, p1y, sy)
h <- p1x == p2x & zx == p1x & p1x == 0 & sign(zy) !=
sign(p1y)
sx <- ifelse(h, p1x, sx)
sy <- ifelse(h, p2y, sy)
h <- zx == p1x & zy == p1y
sx <- ifelse(h, p1x, sx)
sy <- ifelse(h, p1y, sy)
h <- zx == p2x & zy == p2y
sx <- ifelse(h, p2x, sx)
sy <- ifelse(h, p2y, sy)
h <- zx == 0 & zy == 0
sx <- ifelse(h, 0, sx)
sy <- ifelse(h, 0, sy)
sxy <- cbind(sx, sy)
}
if (!exists("debug.plots"))
debug.plots <- "no"
if (debug.plots == "all") {
segments(sxy[, 1], sxy[, 2], zx, zy, col = "red")
segments(0, 0, sxy[, 1], sxy[, 2], type = "l", col = "green",
lty = 2)
points(sxy, col = "red")
}
return(sxy)
}
find.cut.z.pg <- function(z, pg, center = c(0, 0), debug.plots = "no") {
if (!is.matrix(z))
z <- rbind(z)
if (1 == nrow(pg))
return(matrix(center, nrow(z), 2, TRUE))
n.pg <- nrow(pg)
n.z <- nrow(z)
z <- cbind(z[, 1] - center[1], z[, 2] - center[2])
pgo <- pg
pg <- cbind(pg[, 1] - center[1], pg[, 2] - center[2])
if (!exists("debug.plots"))
debug.plots <- "no"
if (debug.plots == "all") {
plot(rbind(z, pg, 0), bty = "n")
points(z, pch = "p")
lines(c(pg[, 1], pg[1, 1]), c(pg[, 2], pg[1, 2]))
}
apg <- win(pg[, 1], pg[, 2])
apg[is.nan(apg)] <- 0
a <- order(apg)
apg <- apg[a]
pg <- pg[a, ]
az <- win(z[, 1], z[, 2])
segm.no <- apply((outer(apg, az, "<")), 2, sum)
segm.no <- ifelse(segm.no == 0, n.pg, segm.no)
next.no <- 1 + (segm.no%%length(apg))
cuts <- cut.z.pg(z[, 1], z[, 2], pg[segm.no, 1], pg[segm.no,
2], pg[next.no, 1], pg[next.no, 2])
cuts <- cbind(cuts[, 1] + center[1], cuts[, 2] + center[2])
return(cuts)
}
hdepth.of.points <- function(tp, n) {
n.tp <- nrow(tp)
tphdepth <- rep(0, n.tp)
dpi <- 2 * pi - 1e-06
minusplus <- c(rep(-1, n), rep(1, n))
for (j in 1:n.tp) {
dx <- tp[j, 1] - xy[, 1]
dy <- tp[j, 2] - xy[, 2]
a <- win(dx, dy) + pi
h <- a < 10
a <- a[h]
ident <- sum(!h)
init <- sum(a < pi)
a.shift <- (a + pi) %% dpi
minusplus <- c(rep(-1, length(a)), rep(1, length(a)))
h <- cumsum(minusplus[order(c(a, a.shift))])
tphdepth[j] <- init + min(h) + 1
}
tphdepth
}
expand.hull <- function(pg, k) {
resolution <- floor(20 * precision)
pg0 <- xy[hdepth == 1, ]
pg0 <- pg0[chull(pg0[, 1], pg0[, 2]), ]
end.points <- find.cut.z.pg(pg, pg0, center = center,
debug.plots = debug.plots)
lam <- ((0 : resolution)^1) / resolution^1
pg.new <- pg
for (i in 1:nrow(pg)) {
tp <- cbind(pg[i, 1] + lam * (end.points[i, 1] -
pg[i, 1]), pg[i, 2] + lam * (end.points[i, 2] -
pg[i, 2]))
hd.tp <- hdepth.of.points(tp, nrow(xy))
ind <- max(sum(hd.tp >= k), 1)
if (ind < length(hd.tp)) {
tp <- cbind(tp[ind, 1] + lam * (tp[ind + 1, 1] -
tp[ind, 1]), tp[ind, 2] + lam * (tp[ind + 1,
2] - tp[ind, 2]))
hd.tp <- hdepth.of.points(tp, nrow(xy))
ind <- max(sum(hd.tp >= k), 1)
}
pg.new[i, ] <- tp[ind, ]
}
pg.new <- pg.new[chull(pg.new[, 1], pg.new[, 2]), ]
pg.add <- 0.5 * (pg.new + rbind(pg.new[-1, ], pg.new[1,
]))
end.points <- find.cut.z.pg(pg.add, pg0, center = center)
tp <- cbind(pg.add[i, 1] + lam * (end.points[i, 1] -
hd.tp <- hdepth.of.points(tp, nrow(xy))
ind <- max(sum(hd.tp >= k), 1)
if (ind < length(hd.tp)) {
tp <- cbind(tp[ind, 1] + lam * (tp[ind + 1, 1] -
tp[ind, 1]), tp[ind, 2] + lam * (tp[ind + 1,
2] - tp[ind, 2]))
hd.tp <- hdepth.of.points(tp, nrow(xy))
ind <- max(sum(hd.tp >= k), 1)
}
}
pg.new <- pg.new[chull(pg.new[, 1], pg.new[, 2]), ]
}
cut.p.sl.p.sl <- function(xy1, m1, xy2, m2) {
sx <- (xy2[2] - m2 * xy2[1] - xy1[2] + m1 * xy1[1])/(m1 -
m2)
sy <- xy1[2] - m1 * xy1[1] + m1 * sx
if (!is.nan(sy))
return(c(sx, sy))
if (abs(m1) == Inf)
return(c(xy1[1], xy2[2] + m2 * (xy1[1] - xy2[1])))
if (abs(m2) == Inf)
return(c(xy2[1], xy1[2] + m1 * (xy2[1] - xy1[1])))
}
pos.to.pg <- function(z, pg, reverse = FALSE) {
if (reverse) {
int.no <- apply(outer(pg[, 1], z[, 1], ">="), 2,
sum)
zy.on.pg <- pg[int.no, 2] + pg[int.no, 3] * (z[,
1] - pg[int.no, 1])
}
else {
int.no <- apply(outer(pg[, 1], z[, 1], "<="), 2,
sum)
zy.on.pg <- pg[int.no, 2] + pg[int.no, 3] * (z[,
1] - pg[int.no, 1])
}
ifelse(z[, 2] < zy.on.pg, "lower", "higher")
}
find.polygon.center <- function(xy) {
if (length(xy) == 2)
return(xy[1:2])
n <- length(xy[, 1])
mxy <- colMeans(xy)
xy2 <- rbind(xy[-1, ], xy[1, ])
xy3 <- cbind(rep(mxy[1], n), mxy[2])
S <- (xy + xy2 + xy3)/3
F2 <- abs((xy[, 1] - xy3[, 1]) * (xy2[, 2] - xy3[, 2]) -
(xy[, 2] - xy3[, 2]) * (xy2[, 1] - xy3[, 1]))
lambda <- F2/sum(F2)
SP <- colSums(cbind(S[, 1] * lambda, S[, 2] * lambda))
return(SP)
}
xydata <- if (missing(y))
x
else cbind(x, y)
if (is.data.frame(xydata))
xydata <- as.matrix(xydata)
if (any(is.na(xydata))) {
if (na.rm) {
xydata <- xydata[!apply(is.na(xydata), 1, any), ,
drop = FALSE]
print("Warning: NA elements have been removed!!")
}
else {
xy.means <- colMeans(xydata, na.rm = TRUE)
for (j in 1:ncol(xydata)) xydata[is.na(xydata[, j]),
j] <- xy.means[j]
print("Warning: NA elements have been exchanged by mean values!!")
}
}
if (nrow(xydata) < 3) {
print("not enough data points")
return()
}
very.large.data.set <- nrow(xydata) > approx.limit
set.seed(random.seed <- 13)
if (very.large.data.set) {
ind <- sample(seq(nrow(xydata)), size = approx.limit)
xy <- xydata[ind, ]
}
else xy <- xydata
n <- nrow(xy)
points.in.bag <- floor(n/2)
if (verbose)
cat("end of initialization")
prdata <- prcomp(xydata)
is.one.dim <- (min(prdata[[1]])/max(prdata[[1]])) < 1e-04
if (is.one.dim) {
if (verbose)
cat("data set one dimensional")
center <- colMeans(xydata)
res <- list(xy = xy, xydata = xydata, prdata = prdata,
is.one.dim = is.one.dim, center = center)
class(res) <- "bagplot"
return(res)
}
if (verbose)
cat("data not linear")
xym <- apply(xy, 2, mean)
xysd <- apply(xy, 2, sd)
xyxy <- cbind((xy[, 1] - xym[1])/xysd[1], (xy[, 2] - xym[2])/xysd[2])
dx <- (outer(xy[, 1], xy[, 1], "-"))
dy <- (outer(xy[, 2], xy[, 2], "-"))
alpha <- atan2(y = dy, x = dx)
diag(alpha) <- 1000
for (j in 1:n) alpha[, j] <- sort(alpha[, j])
alpha <- alpha[-n, ]
m <- n - 1
if (debug.plots == "all") {
plot(xy, bty = "n")
xdelta <- abs(diff(range(xy[, 1])))
dx <- xdelta * 0.3
for (j in 1:n) {
p <- xy[j, ]
dy <- dx * tan(alpha[, j])
segments(p[1] - dx, p[2] - dy, p[1] + dx, p[2] +
dy, col = j)
text(p[1] - xdelta * 0.02, p[2], j, col = j)
}
}
if (verbose)
print("end of computation of angles")
hdepth <- rep(0, n)
dpi <- 2 * pi - 1e-06
mypi <- pi - 1e-06
minusplus <- c(rep(-1, m), rep(1, m))
for (j in 1:n) {
a <- alpha[, j] + pi
h <- a < 10
a <- a[h]
init <- sum(a < mypi)
a.shift <- (a + pi) %% dpi
minusplus <- c(rep(-1, length(a)), rep(1, length(a)))
h <- cumsum(minusplus[order(c(a, a.shift))])
hdepth[j] <- init + min(h) + 1
}
if (verbose) {
print("end of computation of hdepth:")
print(hdepth)
}
if (debug.plots == "all") {
plot(xy, bty = "n")
xdelta <- abs(diff(range(xy[, 1])))
dx <- xdelta * 0.1
for (j in 1:n) {
a <- alpha[, j] + pi
a <- a[a < 10]
init <- sum(a < pi)
a.shift <- (a + pi) %% dpi
minusplus <- c(rep(-1, length(a)), rep(1, length(a)))
h <- cumsum(minusplus[ao <- (order(c(a, a.shift)))])
no <- which((init + min(h)) == (init + h))[1]
p <- xy[j, ]
dy <- dx * tan(alpha[, j])
segments(p[1] - dx, p[2] - dy, p[1] + dx, p[2] +
dy, col = j, lty = 3)
dy <- dx * tan(c(sort(a), sort(a))[no])
segments(p[1] - 5 * dx, p[2] - 5 * dy, p[1] + 5 *
dx, p[2] + 5 * dy, col = "black")
text(p[1] - xdelta * 0.02, p[2], hdepth[j], col = 1,
cex = 2.5)
}
}
hd.table <- table(sort(hdepth))
d.k <- cbind(dk = rev(cumsum(rev(hd.table))), k = as.numeric(names(hd.table)))
k.1 <- sum(points.in.bag < d.k[, 1])
if (nrow(d.k) > k.1) {
k <- d.k[k.1 + 1, 2]
}
else {
k <- d.k[k.1, 2]
}
if (verbose) {
cat("numbers of members of dk:")
print(hd.table)
}
if (verbose) {
cat("end of computation of k, k=", k)
}
center <- apply(xy[which(hdepth == max(hdepth)), , drop = FALSE],
2, mean)
hull.center <- NULL
if (5 < nrow(xy) && length(hd.table) > 2) {
n.p <- floor(c(32, 16, 8)[1 + (n > 50) + (n > 200)] *
precision)
h <- cands <- xy[rev(order(hdepth))[1:6], ]
cands <- cands[chull(cands[, 1], cands[, 2]), ]
n.c <- nrow(cands)
if (is.null(n.c))
cands <- h
xyextr <- rbind(apply(cands, 2, min), apply(cands, 2,
max))
xydel <- 2 * (xyextr[2, ] - xyextr[1, ])/n.p
h1 <- seq(xyextr[1, 1], xyextr[2, 1], length = n.p)
h2 <- seq(xyextr[1, 2], xyextr[2, 2], length = n.p)
tp <- cbind(matrix(h1, n.p, n.p)[1:n.p^2], matrix(h2,
n.p, n.p, TRUE)[1:n.p^2])
tphdepth <- max(hdepth.of.points(tp, n)) - 1
num <- floor(c(417, 351, 171, 85, 67, 43)[sum(n > c(1,
50, 100, 150, 200, 250))] * precision)
num.h <- floor(num/2)
angles <- seq(0, pi, length = num.h)
ang <- tan(pi/2 - angles)
kkk <- tphdepth
if (verbose) {
cat("max-hdepth found:")
print(kkk)
}
ia <- 1
a <- angles[ia]
xyt <- xyxy %*% c(cos(a), -sin(a))
xyto <- order(xyt)
ind.k <- xyto[kkk]
cutp <- c(xyxy[ind.k, 1], -10)
dxy <- diff(range(xyxy))
pg <- rbind(c(cutp[1], -dxy, Inf), c(cutp[1], dxy, NA))
ind.kk <- xyto[n + 1 - kkk]
cutpl <- c(xyxy[ind.kk, 1], 10)
pgl <- rbind(c(cutpl[1], dxy, Inf), c(cutpl[1], -dxy,
NA))
if (debug.plots == "all") {
plot(xyxy, type = "p", bty = "n")
}
for (ia in seq(angles)[-1]) {
a <- angles[ia]
angtan <- ang[ia]
xyt <- xyxy %*% c(cos(a), -sin(a))
xyto <- order(xyt)
ind.k <- xyto[kkk]
ind.kk <- xyto[n + 1 - kkk]
pnew <- xyxy[ind.k, ]
pnewl <- xyxy[ind.kk, ]
if (debug.plots == "all")
points(pnew[1], pnew[2], col = "red")
if (abs(angtan) > 1e+10) {
pg.no <- sum(pg[, 1] < pnew[1])
cutp <- c(pnew[1], pg[pg.no, 2] + pg[pg.no, 3] *
(pnew[1] - pg[pg.no, 1]))
pg.nol <- sum(pgl[, 1] >= pnewl[1])
cutpl <- c(pnewl[1], pgl[pg.nol, 2] + pgl[pg.nol,
3] * (pnewl[1] - pgl[pg.nol, 1]))
}
else {
pg.inter <- pg[, 2] - angtan * pg[, 1]
pnew.inter <- pnew[2] - angtan * pnew[1]
pg.no <- sum(pg.inter < pnew.inter)
cutp <- cut.p.sl.p.sl(pnew, ang[ia], pg[pg.no,
1:2], pg[pg.no, 3])
pg.interl <- pgl[, 2] - angtan * pgl[, 1]
pnew.interl <- pnewl[2] - angtan * pnewl[1]
pg.nol <- sum(pg.interl > pnew.interl)
cutpl <- cut.p.sl.p.sl(pnewl, angtan, pgl[pg.nol,
1:2], pgl[pg.nol, 3])
}
pg <- rbind(pg[1:pg.no, ], c(cutp, angtan), c(cutp[1] +
dxy, cutp[2] + angtan * dxy, NA))
pgl <- rbind(pgl[1:pg.nol, ], c(cutpl, angtan), c(cutpl[1] -
dxy, cutpl[2] - angtan * dxy, NA))
if (debug.plots == "all") {
points(pnew[1], pnew[2], col = "red")
hx <- xyxy[ind.k, c(1, 1)]
hy <- xyxy[ind.k, c(2, 2)]
segments(hx, hy, c(10, -10), hy + ang[ia] * (c(10,
-10) - hx), lty = 2)
points(cutpl[1], cutpl[2], col = "red")
hx <- xyxy[ind.kk, c(1, 1)]
hy <- xyxy[ind.kk, c(2, 2)]
segments(hx, hy, c(10, -10), hy + ang[ia] * (c(10,
-10) - hx), lty = 2)
}
}
limit <- 1e-10
pg <- pg[c(TRUE, (abs(diff(pg[, 1])) > limit) | (abs(diff(pg[,
2])) > limit)), ]
pgl <- pgl[c(TRUE, (abs(diff(pgl[, 1])) > limit) | (abs(diff(pgl[,
2])) > limit)), ]
pg <- pg[-nrow(pg), ][-1, , drop = FALSE]
pgl <- pgl[-nrow(pgl), ][-1, , drop = FALSE]
indl <- pos.to.pg(pgl, pg)
indu <- pos.to.pg(pg, pgl, TRUE)
sr <- sl <- NULL
if (indu[(npg <- nrow(pg))] == "lower" & indl[1] == "higher") {
rnuml <- which(indl == "lower")[1] - 1
rnumu <- npg + 1 - which(rev(indu == "higher"))[1]
if (is.na(rnuml))
rnuml <- sum(pg[rnumu, 1] < pgl[, 1])
if (is.na(rnumu))
rnumu <- sum(pg[, 1] < pgl[rnuml, 1])
xyl <- pgl[rnuml, ]
xyu <- pg[rnumu, ]
sr <- cut.p.sl.p.sl(xyl[1:2], xyl[3], xyu[1:2], xyu[3])
}
if (indl[(npgl <- nrow(pgl))] == "higher" & indu[1] ==
"lower") {
lnuml <- npgl + 1 - which(rev(indl == "lower"))[1]
lnumu <- which(indu == "higher")[1] - 1
if (is.na(lnuml))
lnuml <- sum(pg[lnumu, 1] < pgl[, 1])
if (is.na(lnumu))
lnumu <- sum(pg[, 1] < pgl[lnuml, 1])
xyl <- pgl[lnuml, ]
xyu <- pg[lnumu, ]
sl <- cut.p.sl.p.sl(xyl[1:2], xyl[3], xyu[1:2], xyu[3])
}
pg <- rbind(pg[indu == "higher", 1:2, drop = FALSE],
sr, pgl[indl == "lower", 1:2, drop = FALSE], sl)
if (debug.plots == "all")
lines(rbind(pg, pg[1, ]), col = "red")
pg <- pg[chull(pg[, 1], pg[, 2]), ]
hull.center <- cbind(pg[, 1] * xysd[1] + xym[1], pg[,
2] * xysd[2] + xym[2])
center <- find.polygon.center(hull.center)
if (verbose) {
cat("hull.center", hull.center)
print(table(tphdepth))
}
}
if (verbose)
cat("center depth:", hdepth.of.points(rbind(center),
n) - 1)
if (verbose) {
print("end of computation of center")
print(center)
}
if (dkmethod == 1) {
xyi <- xy[hdepth >= k, , drop = FALSE]
pdk <- xyi[chull(xyi[, 1], xyi[, 2]), , drop = FALSE]
xyo <- xy[hdepth >= (k - 1), , drop = FALSE]
pdk.1 <- xyo[chull(xyo[, 1], xyo[, 2]), , drop = FALSE]
if (verbose)
cat("hull computed:")
if (debug.plots == "all") {
plot(xy, bty = "n")
h <- rbind(pdk, pdk[1, ])
lines(h, col = "red", lty = 2)
h <- rbind(pdk.1, pdk.1[1, ])
lines(h, col = "blue", lty = 3)
points(center[1], center[2], pch = 8, col = "red")
}
exp.dk <- expand.hull(pdk, k)
exp.dk.1 <- expand.hull(exp.dk, k - 1)
}
else {
num <- floor(c(417, 351, 171, 85, 67, 43)[sum(n > c(1,
50, 100, 150, 200, 250))] * precision)
num.h <- floor(num/2)
angles <- seq(0, pi, length = num.h)
ang <- tan(pi/2 - angles)
kkk <- k
ia <- 1
a <- angles[ia]
xyt <- xyxy %*% c(cos(a), -sin(a))
xyto <- order(xyt)
ind.k <- xyto[kkk]
cutp <- c(xyxy[ind.k, 1], -10)
dxy <- diff(range(xyxy))
pg <- rbind(c(cutp[1], -dxy, Inf), c(cutp[1], dxy, NA))
ind.kk <- xyto[n + 1 - kkk]
cutpl <- c(xyxy[ind.kk, 1], 10)
pgl <- rbind(c(cutpl[1], dxy, Inf), c(cutpl[1], -dxy,
NA))
if (debug.plots == "all") {
plot(xyxy, type = "p", bty = "n")
}
for (ia in seq(angles)[-1]) {
a <- angles[ia]
angtan <- ang[ia]
xyt <- xyxy %*% c(cos(a), -sin(a))
xyto <- order(xyt)
ind.k <- xyto[kkk]
ind.kk <- xyto[n + 1 - kkk]
pnew <- xyxy[ind.k, ]
pnewl <- xyxy[ind.kk, ]
if (debug.plots == "all")
points(pnew[1], pnew[2], col = "red")
if (abs(angtan) > 1e+10) {
pg.no <- sum(pg[, 1] < pnew[1])
cutp <- c(pnew[1], pg[pg.no, 2] + pg[pg.no, 3] *
(pnew[1] - pg[pg.no, 1]))
pg.nol <- sum(pgl[, 1] >= pnewl[1])
cutpl <- c(pnewl[1], pgl[pg.nol, 2] + pgl[pg.nol,
3] * (pnewl[1] - pgl[pg.nol, 1]))
}
else {
pg.inter <- pg[, 2] - angtan * pg[, 1]
pnew.inter <- pnew[2] - angtan * pnew[1]
pg.no <- sum(pg.inter < pnew.inter)
cutp <- cut.p.sl.p.sl(pnew, ang[ia], pg[pg.no,
1:2], pg[pg.no, 3])
pg.interl <- pgl[, 2] - angtan * pgl[, 1]
pnew.interl <- pnewl[2] - angtan * pnewl[1]
pg.nol <- sum(pg.interl > pnew.interl)
cutpl <- cut.p.sl.p.sl(pnewl, angtan, pgl[pg.nol,
1:2], pgl[pg.nol, 3])
}
pg <- rbind(pg[1:pg.no, ], c(cutp, angtan), c(cutp[1] +
dxy, cutp[2] + angtan * dxy, NA))
pgl <- rbind(pgl[1:pg.nol, ], c(cutpl, angtan), c(cutpl[1] -
dxy, cutpl[2] - angtan * dxy, NA))
if (debug.plots == "all") {
points(pnew[1], pnew[2], col = "red")
hx <- xyxy[ind.k, c(1, 1)]
hy <- xyxy[ind.k, c(2, 2)]
segments(hx, hy, c(10, -10), hy + ang[ia] * (c(10,
-10) - hx), lty = 2)
points(cutpl[1], cutpl[2], col = "red")
hx <- xyxy[ind.kk, c(1, 1)]
hy <- xyxy[ind.kk, c(2, 2)]
segments(hx, hy, c(10, -10), hy + ang[ia] * (c(10,
-10) - hx), lty = 2)
}
}
limit <- 1e-10
pg <- pg[c(TRUE, (abs(diff(pg[, 1])) > limit) | (abs(diff(pg[,
2])) > limit)), ]
pgl <- pgl[c(TRUE, (abs(diff(pgl[, 1])) > limit) | (abs(diff(pgl[,
2])) > limit)), ]
pg <- pg[-nrow(pg), ][-1, , drop = FALSE]
pgl <- pgl[-nrow(pgl), ][-1, , drop = FALSE]
indl <- pos.to.pg(pgl, pg)
indu <- pos.to.pg(pg, pgl, TRUE)
sr <- sl <- NULL
if (indu[(npg <- nrow(pg))] == "lower" & indl[1] == "higher") {
rnuml <- which(indl == "lower")[1] - 1
rnumu <- npg + 1 - which(rev(indu == "higher"))[1]
if (is.na(rnuml))
rnuml <- sum(pg[rnumu, 1] < pgl[, 1])
if (is.na(rnumu))
rnumu <- sum(pg[, 1] < pgl[rnuml, 1])
xyl <- pgl[rnuml, ]
xyu <- pg[rnumu, ]
sr <- cut.p.sl.p.sl(xyl[1:2], xyl[3], xyu[1:2], xyu[3])
}
if (indl[(npgl <- nrow(pgl))] == "higher" & indu[1] ==
"lower") {
lnuml <- npgl + 1 - which(rev(indl == "lower"))[1]
lnumu <- which(indu == "higher")[1] - 1
if (is.na(lnuml))
lnuml <- sum(pg[lnumu, 1] < pgl[, 1])
if (is.na(lnumu))
lnumu <- sum(pg[, 1] < pgl[lnuml, 1])
xyl <- pgl[lnuml, ]
xyu <- pg[lnumu, ]
sl <- cut.p.sl.p.sl(xyl[1:2], xyl[3], xyu[1:2], xyu[3])
}
pg <- rbind(pg[indu == "higher", 1:2, drop = FALSE],
sr, pgl[indl == "lower", 1:2, drop = FALSE], sl)
if (debug.plots == "all")
lines(rbind(pg, pg[1, ]), col = "red")
pg <- pg[chull(pg[, 1], pg[, 2]), ]
exp.dk <- cbind(pg[, 1] * xysd[1] + xym[1], pg[, 2] *
xysd[2] + xym[2])
if (kkk > 1)
kkk <- kkk - 1
ia <- 1
a <- angles[ia]
xyt <- xyxy %*% c(cos(a), -sin(a))
xyto <- order(xyt)
ind.k <- xyto[kkk]
cutp <- c(xyxy[ind.k, 1], -10)
dxy <- diff(range(xyxy))
pg <- rbind(c(cutp[1], -dxy, Inf), c(cutp[1], dxy, NA))
ind.kk <- xyto[n + 1 - kkk]
cutpl <- c(xyxy[ind.kk, 1], 10)
pgl <- rbind(c(cutpl[1], dxy, Inf), c(cutpl[1], -dxy,
NA))
if (debug.plots == "all") {
plot(xyxy, type = "p", bty = "n")
}
for (ia in seq(angles)[-1]) {
a <- angles[ia]
angtan <- ang[ia]
xyt <- xyxy %*% c(cos(a), -sin(a))
xyto <- order(xyt)
ind.k <- xyto[kkk]
ind.kk <- xyto[n + 1 - kkk]
pnew <- xyxy[ind.k, ]
pnewl <- xyxy[ind.kk, ]
if (debug.plots == "all")
points(pnew[1], pnew[2], col = "red")
if (abs(angtan) > 1e+10) {
pg.no <- sum(pg[, 1] < pnew[1])
cutp <- c(pnew[1], pg[pg.no, 2] + pg[pg.no, 3] *
(pnew[1] - pg[pg.no, 1]))
pg.nol <- sum(pgl[, 1] >= pnewl[1])
cutpl <- c(pnewl[1], pgl[pg.nol, 2] + pgl[pg.nol,
3] * (pnewl[1] - pgl[pg.nol, 1]))
}
else {
pg.inter <- pg[, 2] - angtan * pg[, 1]
pnew.inter <- pnew[2] - angtan * pnew[1]
pg.no <- sum(pg.inter < pnew.inter)
cutp <- cut.p.sl.p.sl(pnew, ang[ia], pg[pg.no,
1:2], pg[pg.no, 3])
pg.interl <- pgl[, 2] - angtan * pgl[, 1]
pnew.interl <- pnewl[2] - angtan * pnewl[1]
pg.nol <- sum(pg.interl > pnew.interl)
cutpl <- cut.p.sl.p.sl(pnewl, angtan, pgl[pg.nol,
1:2], pgl[pg.nol, 3])
}
pg <- rbind(pg[1:pg.no, ], c(cutp, angtan), c(cutp[1] +
dxy, cutp[2] + angtan * dxy, NA))
pgl <- rbind(pgl[1:pg.nol, ], c(cutpl, angtan), c(cutpl[1] -
dxy, cutpl[2] - angtan * dxy, NA))
if (debug.plots == "all") {
points(pnew[1], pnew[2], col = "red")
hx <- xyxy[ind.k, c(1, 1)]
hy <- xyxy[ind.k, c(2, 2)]
segments(hx, hy, c(10, -10), hy + ang[ia] * (c(10,
-10) - hx), lty = 2)
points(cutpl[1], cutpl[2], col = "red")
hx <- xyxy[ind.kk, c(1, 1)]
hy <- xyxy[ind.kk, c(2, 2)]
segments(hx, hy, c(10, -10), hy + ang[ia] * (c(10,
-10) - hx), lty = 2)
}
}
limit <- 1e-10
pg <- pg[c(TRUE, (abs(diff(pg[, 1])) > limit) | (abs(diff(pg[,
2])) > limit)), ]
pgl <- pgl[c(TRUE, (abs(diff(pgl[, 1])) > limit) | (abs(diff(pgl[,
2])) > limit)), ]
pg <- pg[-nrow(pg), ][-1, , drop = FALSE]
pgl <- pgl[-nrow(pgl), ][-1, , drop = FALSE]
indl <- pos.to.pg(pgl, pg)
indu <- pos.to.pg(pg, pgl, TRUE)
sr <- sl <- NULL
if (indu[(npg <- nrow(pg))] == "lower" & indl[1] == "higher") {
rnuml <- which(indl == "lower")[1] - 1
rnumu <- npg + 1 - which(rev(indu == "higher"))[1]
if (is.na(rnuml))
rnuml <- sum(pg[rnumu, 1] < pgl[, 1])
if (is.na(rnumu))
rnumu <- sum(pg[, 1] < pgl[rnuml, 1])
xyl <- pgl[rnuml, ]
xyu <- pg[rnumu, ]
sr <- cut.p.sl.p.sl(xyl[1:2], xyl[3], xyu[1:2], xyu[3])
}
if (indl[(npgl <- nrow(pgl))] == "higher" & indu[1] ==
"lower") {
lnuml <- npgl + 1 - which(rev(indl == "lower"))[1]
lnumu <- which(indu == "higher")[1] - 1
if (is.na(lnuml))
lnuml <- sum(pg[lnumu, 1] < pgl[, 1])
if (is.na(lnumu))
lnumu <- sum(pg[, 1] < pgl[lnuml, 1])
xyl <- pgl[lnuml, ]
xyu <- pg[lnumu, ]
sl <- cut.p.sl.p.sl(xyl[1:2], xyl[3], xyu[1:2], xyu[3])
}
pg <- rbind(pg[indu == "higher", 1:2, drop = FALSE],
sr, pgl[indl == "lower", 1:2, drop = FALSE], sl)
if (debug.plots == "all")
lines(rbind(pg, pg[1, ]), col = "red")
pg <- pg[chull(pg[, 1], pg[, 2]), ]
exp.dk.1 <- cbind(pg[, 1] * xysd[1] + xym[1], pg[, 2] *
xysd[2] + xym[2])
}
if (nrow(d.k) == k.1 || nrow(d.k) == 1)
lambda <- 0
else {
lambda <- (n/2 - d.k[k.1 + 1, 1])/(d.k[k.1, 1] - d.k[k.1 +
1, 1])
}
if (verbose)
cat("lambda", lambda)
cut.on.pdk.1 <- find.cut.z.pg(exp.dk, exp.dk.1, center = center)
cut.on.pdk <- find.cut.z.pg(exp.dk.1, exp.dk, center = center)
h1 <- (1 - lambda) * exp.dk + lambda * cut.on.pdk.1
h2 <- (1 - lambda) * cut.on.pdk + lambda * exp.dk.1
h <- rbind(h1, h2)
h <- h[!is.nan(h[, 1]) & !is.nan(h[, 2]), ]
hull.bag <- h[chull(h[, 1], h[, 2]), ]
if (verbose)
cat("bag completed:")
if (debug.plots == "all") {
lines(hull.bag, col = "red")
}
hull.loop <- cbind(hull.bag[, 1] - center[1], hull.bag[,
2] - center[2])
hull.loop <- factor * hull.loop
hull.loop <- cbind(hull.loop[, 1] + center[1], hull.loop[,
2] + center[2])
if (verbose)
cat("loop computed")
if (!very.large.data.set) {
pxy.bag <- xydata[hdepth >= k, , drop = FALSE]
pkt.cand <- xydata[hdepth == (k - 1), , drop = FALSE]
pkt.not.bag <- xydata[hdepth < (k - 1), , drop = FALSE]
if (length(pkt.cand) > 0) {
outside <- out.of.polygon(pkt.cand, hull.bag)
if (sum(!outside) > 0)
pxy.bag <- rbind(pxy.bag, pkt.cand[!outside,
])
if (sum(outside) > 0)
pkt.not.bag <- rbind(pkt.not.bag, pkt.cand[outside,
])
}
}
else {
extr <- out.of.polygon(xydata, hull.bag)
pxy.bag <- xydata[!extr, ]
pkt.not.bag <- xydata[extr, , drop = FALSE]
}
if (length(pkt.not.bag) > 0) {
extr <- out.of.polygon(pkt.not.bag, hull.loop)
pxy.outlier <- pkt.not.bag[extr, , drop = FALSE]
if (0 == length(pxy.outlier))
pxy.outlier <- NULL
pxy.outer <- pkt.not.bag[!extr, , drop = FALSE]
}
else {
pxy.outer <- pxy.outlier <- NULL
}
if (verbose)
cat("points of bag, outer points and outlier identified")
hull.loop <- rbind(pxy.outer, hull.bag)
hull.loop <- hull.loop[chull(hull.loop[, 1], hull.loop[,
2]), ]
if (verbose)
cat("end of computation of loop")
res <- list(center = center, hull.center = hull.center, hull.bag = hull.bag,
hull.loop = hull.loop, pxy.bag = pxy.bag, pxy.outer = if (length(pxy.outer) >
0) pxy.outer else NULL, pxy.outlier = if (length(pxy.outlier) >
0) pxy.outlier else NULL, hdepths = hdepth, is.one.dim = is.one.dim,
prdata = prdata, random.seed = random.seed, xy = xy,
xydata = xydata)
if (verbose)
res <- c(res, list(exp.dk = exp.dk, exp.dk.1 = exp.dk.1,
hdepth = hdepth))
class(res) <- "bagplot"
return(res)
}
```

Try the rainbow package in your browser

Any scripts or data that you put into this service are public.

rainbow documentation built on Jan. 29, 2019, 5:09 p.m.