library(ks)
library(colorspace)
library(RColorBrewer)

## 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))}
seq.col3d <- function(n, alpha=1, ...){rev(sequential_hcl(n=n, h=290, power=1.1, c.=c(80,70), l=c(30,50), alpha=alpha, ...))}
heat.col <- function(n, alpha=1, ...){heat_hcl(n, h=c(0, -100), l=c(75, 40), c=c(40, 80), power=1, alpha=alpha, ...)}
terrain.col <- function(n, alpha=1, ...){terrain_hcl(n, c=c(65, 0), l=c(45, 95), power=c(1/3, 1.5), alpha=alpha, ...)}
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)}

## prepare data
data(tempb)
tempb <- tempb[, c("tmin", "tmax")]

## plot parameters
par(oma=c(0,0,0,0)+0.1, mgp=c(1.8,0.5,0), mar=c(2.9,2.9,0,0)+0.1, cex.axis=1.2, cex.lab=1.2)
xlab <- "Min temperature (deg C)" 
ylab <- "Max temperature (deg C)" 
xlim <- c(-5.5, 27.5)
ylim <- c(5.5, 42.5)

## Figure 1.1
## scattter plot   
plot(tempb, cex=0.7, pch=16, col=pt.col(4, alpha=0.1), asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab) 

## Fig 2.1a
## histogram 
hist.tempb <- histde(tempb)

plot(hist.tempb, col.fun=seq.col, col.rect=grey(0,alpha=0.5), lty.rect=1, nbreaks=5, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(hist.tempb$eval.points[[1]][4], hist.tempb$eval.points[[2]][5], col=1, pch=16, cex=1.6)

## Fig 2.1b
## effect of anchor point placement
hist.tempb2 <- histde(tempb, adj=0.5)

plot(hist.tempb2, col.fun=seq.col, col.rect=grey(0,alpha=0.5), lty.rect=1, nbreaks=5, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(hist.tempb2$eval.points[[1]][4], hist.tempb2$eval.points[[2]][5], col=1, pch=16, cex=1.6)

## Fig 2.2a
## effect of binwidth - undersmoothing
hist.tempb.us <- histde(tempb, binw=hist.tempb$binw/3)

plot(hist.tempb.us, col.fun=seq.col, col.rect=grey(0,alpha=0.5), lty.rect=1, nbreaks=5, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.2b
## effect of binwidth - oversmoothing
hist.tempb.os <- histde(tempb, binw=hist.tempb$binw*3)

plot(hist.tempb.os, col.fun=seq.col, col.rect=grey(0,alpha=0.5), lty.rect=1, nbreaks=5, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.3a
## kernel functions on tempb subset
tempb.sub <- tempb[c(1608,2344,3820,5083,8184,10434,13470,13632,14044,21234),]
nsub <- nrow(tempb.sub)
H <- 10*Hpi(tempb)

for (i in 1:nsub)
{
    fhat.kern1 <- plotmixt(mus=as.numeric(tempb.sub[i,]), Sigmas=H, props=1, draw=FALSE, asp=1, xlim=xlim, ylim=ylim)
    plot(fhat.kern1, add=i>1, drawlabels=FALSE, lwd=1, col=grey(0.1, alpha=0.9), cont=seq(1,nsub,by=3)/nsub*45, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
}
points(tempb.sub, pch=21, bg=pt.col(4))

## Fig 2.3b
## kernel density estimate on tempb subset
fhat.tempb.sub <- kde(x=tempb.sub, H=H, xmin=c(-10.231763, -2.377353), xmax=c(29.03176, 48.55315))

plot(fhat.tempb.sub, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(tempb.sub, pch=21, bg=pt.col(4))

## Fig 2.4a
## kernel density estimate - filled contour plot
fhat.tempb <- kde(x=tempb)
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.4b
## kernel density estimate - perspective plot
plot(fhat.tempb, display="persp", phi=20, theta=20, col.fun=function(n){seq.col(n, alpha=1)[-1]}, border=grey(0, alpha=0.2), thin=3, xlab="Min temp.", ylab="Max temp.", zlab="")

## Fig 2.5a
## filled contour plot - 9 probability contour levels
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col, cont=seq(10,90,by=10), asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.5b
## filled contour plot - 9 linear contour levels
ct.lev <- seq(0.0009, 0.0059, length=9)
	  
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, abs.cont=ct.lev, col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.6a
## filled contour plot - 50% probability contour level
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, cont=50, col=seq.col(9)[c(1,6)], asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.6b
## filled contour plot - 95% probability contour level
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, cont=95, col=seq.col(9)[c(1,2)], asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.7a
## filled contour plot - heat colour scheme 
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=heat.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.7b
## filled contour plot - terrain colour scheme 
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1 col.fun=terrain.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.8a
## unconstrained bandwidth matrix
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, asp=1, col.fun=seq.col, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
plotmixt(mus=c(26,38), Sigmas=fhat.tempb$H, props=1, drawlabels=FALSE, add=TRUE)

## Fig 2.8b
## diagonal bandwidth matrix
plot(fhat.tempb.diag, display="filled.contour", drawlabels=FALSE, lwd=1 col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, ylab=ylab)
plotmixt(mus=c(26,38), Sigmas=fhat.tempb.diag$H, props=1, drawlabels=FALSE, add=TRUE)

## Fig 5.1a
## kernel density gradient estimate 
fhat1.tempb <- kdde(x=tempb, deriv.order=1)

## contour plot - partial derivative f^(1,0)
plot(fhat1.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=1, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.1b
## contour plot - partial derivative f^(0,1)
plot(fhat1.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=2, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.1c
## quiver plot
plot(fhat1.tempb, display="quiver", transf=1/12, thin=6, lwd=2, col.fun=diverge.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.2a
## kernel density Hessian estimate
fhat2.tempb <- kdde(x=tempb, deriv.order=2)
                                                  
## partial derivative f^(2,0)
plot(fhat2.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=1, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.2b
## partial derivative f^(1,1)
plot(fhat2.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=2, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.2c
## partial derivative f^(0,2)
plot(fhat2.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=4, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.2d
## kernel summary curvature estimate
fhat2.tempb.curv <- kcurv(fhat2.tempb)

plot(fhat2.tempb.curv, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col2, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 6.1a
## modal region - 25% contour of density estimate
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, asp=1, cont=25, col.fun=seq.col, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 6.1b
## modal region - 25% contour of summary curvature 
plot(fhat2.tempb.curv, display="filled.contour", drawlabels=FALSE, cont=25, col.fun=seq.col2, lwd=1, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()                                               

## Fig 6.3a
## Devroye-Wise support estimate
supp.tempb.dw <- dwsupp(tempb, H=fhat.tempb$H, xmin=c(-20,-10), xmax=c(35,55))
chull.tempb.dw <- ksupp(supp.tempb.dw, cont=99, convex.hull=TRUE)

plot(fhat.tempb, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(tempb, cex=0.4, pch=16, col=pt.col(4,alpha=0.1))
plot(supp.tempb.dw, add=TRUE, drawlabels=FALSE, lwd=2, col=1, abs.cont=0.5, asp=1)
polygon(chull.tempb.dw, lwd=2, lty=3, border=grey(0,alpha=0.5))

## Fig 6.3b
## kernel support estimate
chull.tempb9995 <- ksupp(fhat.tempb, cont=99.95, convex.hull=TRUE)
chull.tempb <- fhat.tempb$x[chull(fhat.tempb$x),]

plot(fhat.tempb, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(tempb, cex=0.4, pch=16, col=pt.col(4,alpha=0.1))
plot(fhat.tempb, add=TRUE, drawlabels=FALSE, lwd=2, col=1, abs.cont=contourLevels(fhat.tempb, cont=c(95, 99.5, 99.95)))
polygon(chull.tempb9995, lwd=2, lty=2, border=seq.col(3)[3])

## Fig 6.7
## kernel mean shift clustering
## warning: kms may take hours to execute on tempb
ms.tempb <- kms(tempb, verbose=TRUE, min.clust.size=0.05*nrow(tempb))

plot(fhat.tempb, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab)
plot(ms.tempb, pch=16, col=pt.col(c(4,10,2), alpha=0.1), add=TRUE)
points(ms.tempb$mode, pch=21, col=1, bg=pt.col(7))

## Fig 6.11
## kernel feature significance
fs.tempb <- kfs(tempb)

plot(fs.tempb, display="filled.contour", col=pt.col(7), lwd=2, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
plot(fhat2.tempb.curv, cont=25, add=TRUE, drawlabels=FALSE, col=pt.col(10), lwd=5, lty="12")

## Fig 8.2a
## coarse estimation grid size
fhat.tempb.gs1 <- kde(x=tempb, bgridsize=rep(11,2))

plot(fhat.tempb, display="slice", drawlabels=FALSE, col=pt.col(10,0.5), lwd=3, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
plot(fhat.tempb.gs1, display="slice", drawlabels=FALSE, add=TRUE, col=1,lwd=3, lty=3)

## Fig 8.2b
## fine estimation grid 
fhat.tempb.gs2 <- kde(x=tempb, bgridsize=rep(301,2))

plot(fhat.tempb, display="slice", drawlabels=FALSE, col=pt.col(10,0.5), lwd=3, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
plot(fhat.tempb.gs2, display="slice", drawlabels=FALSE, add=TRUE, col=1,lwd=3, lty=3)