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))}
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)
xlim <- c(-5.6,5.6)
ylim <- c(-5.6,5.6)

## dumbbell density random sample
mus <- rbind(c(-2,2), c(0,0), c(2,-2))
Sigmas <- rbind(diag(2), 0.8*invvech(c(1, -0.9, 1)), diag(2))
cwt <- 3/11
props <- c((1-cwt)/2, cwt, (1-cwt)/2)
samp <- 1000
set.seed(81928192)
x <- rmvnorm.mixt(n=samp, mus=mus, Sigmas=Sigmas, props=props)

## Fig 2.10a
## target dumbbell density 
target.db <- plotmixt(mus=mus, Sigmas=Sigmas, props=props, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE)

plot(target.db, drawlabels=FALSE, axes=FALSE, lwd=1, display="filled.contour", col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="")
box()

## Fig 2.10b
## kernel density estimate with unconstrained bandwidth
fhat.db <- kde(x=x, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2]))
db.single <- plotmixt(mus=c(3,3), Sigmas=fhat.db$H, props=1, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE)

plot(fhat.db, drawlabels=FALSE, axes=FALSE, lwd=1, display="filled.contour", col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="")
plot(db.single, add=TRUE, drawlabels=FALSE, lwd=1)
box()

## Fig 2.10c
## kernel density estimate with diagonal bandwidth
fhat.db.diag <- kde(x=x, H=Hpi.diag(x), xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2]))
db.single.diag <- plotmixt(mus=c(3,3), Sigmas=fhat.db.diag$H, props=1, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE)

plot(fhat.db.diag, drawlabels=FALSE, axes=FALSE, lwd=1, display="filled.contour", col.fun=seq.col, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="")
plot(db.single.diag, add=TRUE, drawlabels=FALSE, lwd=1)
box()

## Fig 5.4a
## target dumbbell density gradient 
target1.db <- plotmixt(mus=mus, Sigmas=Sigmas, props=props, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE, deriv.order=1)

plot(target1.db, display="quiver", thin=6, lwd=2, transf=1/12, scale=0.8, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="")
box()

## Fig 5.4b
## kernel density gradient estimate with unconstrained bandwidth
fhat1.db <- kdde(x=x, deriv.order=1, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2]))

plot(fhat1.db, display="quiver", thin=6, lwd=2, transf=1/12, scale=0.8, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="")
box()

## Fig 5.4c
## kernel density gradient estimate with diagonal bandwidth
fhat1.db.diag <- kdde(x=x, H=Hpi.diag(x, deriv.order=1), binned=TRUE, deriv.order=1, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2]))

png.create(file="fig5-4c-dumbbell-kdde.png")
plot(fhat1.db.diag, display="quiver", thin=6, lwd=2, transf=1/12, scale=0.8, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="")
box()

## Fig 5.4d
## target dumbbell density Hessian
target2.db <- plotmixt(mus=mus, Sigmas=Sigmas, props=props, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE, deriv.order=2)
target2.db.curv <- kcurv(target2.db)

png.create(file="fig5-4d-dumbbell-kdde.png")
plot(target2.db.curv, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col2, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="")
box()

## Fig 5.4e
## kernel density Hessian estimate with unconstrained bandwidth
fhat2.db <- kdde(x=x, deriv.order=2, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2]))
fhat2.db.curv <- kcurv(fhat2.db)

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

## Fig 5.4f
## kernel density Hessian estimate with diagonal bandwidth
fhat2.db.diag <- kdde(x=x, H=Hpi.diag(x, deriv.order=2), binned=TRUE, deriv.order=2, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2]))
fhat2.db.curv.diag <- kcurv(fhat2.db.diag)

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