library(ks)
library(colorspace)
library(RColorBrewer)
library(maps)
library(oz)

## colour functions
pt.col <- function(pos=1, alpha=1, ...) {return(paste0(brewer.pal(12, "Paired"), format(as.hexmode(round(alpha*255,0)), width=2))[pos])}
seq.col <- function(n, alpha=1, ...){c("transparent", rev(sequential_hcl(n=n-1, h=290, power=1.1, alpha=alpha, ...)))}
seq.col2 <- function(n, alpha=1, ...) {c("transparent",tail(rev(diverge.col(2*n-1, alpha=alpha, ...)), n-1))}
diverge.col <- function(n, alpha=1, ...){cols <- paste0(brewer.pal(n, "PuOr"), as.hexmode(round(alpha*255,0))); cols[n %/% 2+1] <- grey(1, alpha=0); return(cols)}

## plot parameters
par(oma=c(0,0,0,0), mgp=c(0,0,0), mar=c(0,0,0,0), cex.axis=1.2, cex.lab=1.2)

## WA map
wa.map <- function(xlim, ylim)
{
    if (missing(xlim)) xlim <- c(114.75, 121.5)
    if (missing(ylim)) ylim <- c(-33, -31.5)
    wa.coast <- ozRegion(section=1)
    col <- rainbow_hcl(12, alpha=0.5)[9]
    oz(section=1, xlim=xlim, ylim=ylim)
    polypath(c(wa.coast$lines[[1]]$x, NA, c(wa.coast$rangex+c(-2,2), rev(wa.coast$rangex+c(-2,2)))), c(wa.coast$lines[[1]]$y, NA, rep(wa.coast$rangey+c(-2,2), each=2)), col=col, rule="evenodd")
    box()
}
wa.coast <- ozRegion(section=1)
wa.polygon <- cbind(wa.coast$lines[[1]]$x, wa.coast$lines[[1]]$y)

## load data
data(grevillea)

## Figure 2.8a 
## scattter plot   
wa.map()
points(grevillea, cex=1.2, pch=21, col=grey(0,alpha=0.5), bg=pt.col(4,alpha=0.5))

## Fig 2.8b
## Kernel density estimate with unconstrained bandwidth 
fhat.grevillea <- kde(x=grevillea)
fhat.grevillea <- kde.truncate(fhat.grevillea, wa.polygon)

wa.map()
plot(fhat.grevillea, add=TRUE, drawlabels=FALSE, display="filled.contour", lwd=1, col.fun=seq.col)
plotmixt(mus=c(120.75,-29.5), Sigmas=fhat.grevillea$H, props=1, drawlabels=FALSE, xlim=c(114.75, 121.5), ylim=c(-32.5, -29), add=TRUE)

## Fig 2.8c
## kernel density estimate with diagonal bandwidth 
fhat.grevillea.diag <- kde(x=grevillea, H=Hpi.diag(grevillea))
fhat.grevillea.diag <- kde.truncate(fhat.grevillea.diag, wa.polygon)
					      
wa.map()
plot(fhat.grevillea.diag, add=TRUE, drawlabels=FALSE, xaxs="i", yaxs="i", display="filled.contour", lwd=1, col.fun=seq.col)
plotmixt(mus=c(120.75,-29.5), Sigmas=fhat.grevillea.diag$H, props=1, drawlabels=FALSE, xlim=c(114.75, 121.5), ylim=c(-32.5, -29), add=TRUE)

## Fig 2.10a
## undersmoothed kernel density estimate
fhat.grevillea.us <- kde(grevillea, H=fhat.grevillea$H %*% fhat.grevillea$H)
fhat.grevillea.us <- kde.truncate(fhat.grevillea.us, wa.polygon)

wa.map()
plot(fhat.grevillea.us, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col)

## Fig 2.10b
## oversmoothed kernel density estimate
fhat.grevillea.os <- kde(grevillea, H=matrix.sqrt(fhat.grevillea$H))
fhat.grevillea.os <- kde.truncate(fhat.grevillea.os, wa.polygon)

wa.map()
plot(fhat.grevillea.os, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col)

## Fig 2.11a
## kernel density estimate with normal scale bandwidth
fhat.grevillea.ns <- kde(grevillea, H=Hns(grevillea))
fhat.grevillea.ns <- kde.truncate(fhat.grevillea.ns, wa.polygon)

wa.map()
plot(fhat.grevillea.ns, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col)

## Fig 2.11b
## kernel density estimate with normal mixture bandwidth
fhat.grevillea.nm <- kde(x=grevillea, H=Hnm(grevillea))
fhat.grevillea.nm <- kde.truncate(fhat.grevillea.nm, wa.polygon)

wa.map()
plot(fhat.grevillea.nm, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col)

## Fig 2.11c
## kernel density estimate with UCV bandwidth
fhat.grevillea.ucv <- kde(grevillea,  H=Hucv(jitter(grevillea)))
fhat.grevillea.ucv <- kde.truncate(fhat.grevillea.ucv, wa.polygon)

wa.map()
plot(fhat.grevillea.ucv, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col)

## Fig 2.11d
## kernel density estimate with BCV bandwidth
fhat.grevillea.bcv <- kde(grevillea, H=Hbcv(grevillea))
fhat.grevillea.bcv <- kde.truncate(fhat.grevillea.bcv, wa.polygon)

wa.map()
plot(fhat.grevillea.bcv, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col)

## Fig 2.11e
## kernel density estimate with plug-in bandwidth
wa.map()
plot(fhat.grevillea, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col)

## Fig 2.11f
## kernel density estimate with SCV bandwidth
fhat.grevillea.scv <- kde(grevillea, H=Hscv(grevillea))
fhat.grevillea.scv <- kde.truncate(fhat.grevillea.scv, wa.polygon)

wa.map()
plot(fhat.grevillea.scv, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col)

## Fig 5.3a
## kernel density gradient estimate with unconstrained bandwidth 
fhat1.grevillea <- kdde(x=grevillea, deriv.order=1)
fhat1.grevillea <- kdde.truncate(fhat1.grevillea, wa.polygon)

wa.map()
plot(fhat1.grevillea, display="quiver", add=TRUE, thin=6, lwd=2, transf=1/8, asp=1)

## Fig 5.3b
## Kernel density gradient estimate with diagonal bandwidth 
fhat1.grevillea.diag <- kdde(x=grevillea, deriv.order=1, H=Hpi.diag(grevillea, deriv.order=1))
fhat1.grevillea.diag <- kdde.truncate(fhat1.grevillea.diag, wa.polygon)

wa.map()
plot(fhat1.grevillea.diag, display="quiver", add=TRUE, thin=6, lwd=2, transf=1/8, asp=1)

## Fig 5.3c
## kernel density curvature estimate with unconstrained bandwidth 
fhat2.grevillea <- kdde(x=grevillea, deriv.order=2)
fhat2.grevillea.curv <- kcurv(fhat2.grevillea)

wa.map()
plot(fhat2.grevillea.curv, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col2)

## Fig 5.3d
## kernel density curvature estimate with diagonal bandwidth 
fhat2.grevillea.diag <- kdde(x=grevillea, deriv.order=2, H=Hpi.diag(grevillea, deriv.order=2))
fhat2.grevillea.curv.diag <- kcurv(fhat2.grevillea.diag)

wa.map()
plot(fhat2.grevillea.curv.diag, display="filled.contour2", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col2)

## Fig 5.6a
## kernel density gradient estimate with normal scale bandwidth
fhat1.grevillea.ns <- kdde(x=grevillea, deriv.order=1, H=Hns(grevillea, deriv.order=1))
fhat1.grevillea.ns <- kdde.truncate(fhat1.grevillea.ns, wa.polygon)

wa.map()
plot(fhat1.grevillea.ns, display="quiver", add=TRUE, thin=6, lwd=2, transf=1/8, asp=1)

## Fig 5.6b
## kernel density gradient estimate with UCV bandwidth
fhat1.grevillea.ucv <- kdde(x=grevillea, deriv.order=1, H=Hucv(jitter(grevillea), deriv.order=1))
fhat1.grevillea.ucv <- kdde.truncate(fhat1.grevillea.ucv, wa.polygon)

wa.map()
plot(fhat1.grevillea.ucv, display="quiver", add=TRUE, thin=6, lwd=2, transf=1/8, asp=1)

## Fig 5.6c
## kernel density gradient estimate with plug-in bandwidth
wa.map()
plot(fhat1.grevillea, display="quiver", add=TRUE, thin=6, lwd=2, transf=1/8, asp=1)

## Fig 5.6d
## kernel density gradient estimate with SCV bandwidth
fhat1.grevillea.scv <- kdde(x=grevillea, deriv.order=1, H=Hscv(grevillea, deriv.order=1))
fhat1.grevillea.scv <- kdde.truncate(fhat1.grevillea.scv, wa.polygon)

wa.map()
plot(fhat1.grevillea.scv, display="quiver", add=TRUE, thin=6, lwd=2, transf=1/8, asp=1)

## Fig 5.7a
## kernel density curvature estimate with normal scale bandwidth
fhat2.grevillea.ns <- kdde(x=grevillea, deriv.order=2, H=Hns(grevillea, deriv.order=2))
fhat2.grevillea.ns <- kdde.truncate(fhat2.grevillea.ns, wa.polygon)
fhat2.grevillea.curv.ns <- kcurv(fhat2.grevillea.ns)

wa.map()
plot(fhat2.grevillea.curv.ns, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col2)

## Fig 5.7b
## kernel density curvature estimate with UCV bandwidth
fhat2.grevillea.ucv <- kdde(x=grevillea, deriv.order=2, H=Hucv(jitter(grevillea), deriv.order=2))
fhat2.grevillea.ucv <- kdde.truncate(fhat2.grevillea.ucv, wa.polygon)
fhat2.grevillea.curv.ucv <- kcurv(fhat2.grevillea.ucv)

wa.map()
plot(fhat2.grevillea.curv.ucv, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col2)

## Fig 5.7c
## kernel density curvature estimate with plug-in bandwidth
wa.map()
plot(fhat2.grevillea.curv, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1, col.fun=seq.col2)

## Fig 5.7d
## kernel density curvature estimate with SCV bandwidth
fhat2.grevillea.scv <- kdde(x=grevillea, deriv.order=2, H=Hscv(grevillea, deriv.order=2))
fhat2.grevillea.scv <- kdde.truncate(fhat2.grevillea.scv, wa.polygon)
fhat2.grevillea.curv.scv <- kcurv(fhat2.grevillea.scv)

wa.map()
plot(fhat2.grevillea.curv.scv, display="filled.contour", add=TRUE, drawlabels=FALSE, lwd=1,  col.fun=seq.col2)

## Fig 6.4
## kernel support/range estimate
chull.grevillea <- fhat.grevillea$x[chull(fhat.grevillea$x),]
chull.grevillea95 <- ksupp(fhat.grevillea, cont=95, convex.hull=TRUE)
chull.grevillea99 <- ksupp(fhat.grevillea, cont=99, convex.hull=TRUE)

wa.map()
points(grevillea, cex=1, pch=21, col=grey(0,alpha=0.5), bg=pt.col(4,alpha=0.5))
plot(fhat.grevillea, display="slice", cont=99, add=TRUE, drawlabels=FALSE, lwd=2)
polygon(chull.grevillea99, lwd=2, lty=2, border=seq.col(3)[3])
polygon(chull.grevillea, lwd=2, lty=3, border=grey(0,alpha=0.5))