#which_vars indexes all of the variables you want to plot implausibility for, IN ORDER
impList_tinydancer <- function(which_vars, var_names, ImpData, resolution=c(15,15), cutoff=3, SpaceLower=0){
  combGrid <- expand.grid(which_vars[-length(which_vars)],which_vars[-1])
  badRows <- c()
  for(i in 1:length(combGrid[,1])){
    if(combGrid[i,1] >= combGrid[i,2])
      badRows <- c(badRows,i)
  }
  if(!is.null(badRows))
    combGrid <- combGrid[-badRows,]
  combGrid <- combGrid[do.call(order,combGrid),]
  gridList <- lapply(which_vars[-length(which_vars)], function(k) combGrid[which(combGrid[,1]==k),])
  parallel::mclapply(gridList, function(e) lapply(1:length(e[,1]), function(k) 
    imp_panel_tinydancer(x1var=var_names[e[k,1]], x2var=var_names[e[k,2]], 
                         ImpData = ImpData, resolution=resolution, 
                         cutoff=cutoff, SpaceLower = SpaceLower)),mc.cores = 8)
}
imp_panel_tinydancer <- function(x1var, x2var, ImpData, resolution=c(10,10), cutoff=3, SpaceLower=0){
  #find pixel boundaries
  x1RHboundaries <- seq(from=SpaceLower,to=1,length=resolution[1]+1)[-1]
  x2RHboundaries <- seq(from=SpaceLower,to=1,length=resolution[2]+1)[-1]
  x1LHboundaries <- c(SpaceLower,x1RHboundaries[-resolution[1]])
  x2LHboundaries <- c(SpaceLower,x2RHboundaries[-resolution[2]])
  whereX1 <- which(colnames(ImpData)==x1var)
  whereX2 <- which(colnames(ImpData)==x2var)
  #create pointers to each pixel
  n <- prod(resolution)
  xs <- 1:resolution[1]
  ys <- 1:resolution[2]
  xyGrid <- expand.grid(xs,ys)
  #What proportion of sample points are in each pixel?
  num_points <- nrow(ImpData)
  ImpOnePixel <- function(index){
    tx1Indices <- which(ImpData[,whereX1] > x1LHboundaries[xyGrid[index,1]] & ImpData[,whereX1] < x1RHboundaries[xyGrid[index,1]])
    tx2Indices <- which(ImpData[,whereX2] > x2LHboundaries[xyGrid[index,2]] & ImpData[,whereX2] < x2RHboundaries[xyGrid[index,2]])
    tIndices <- tx1Indices[which(tx1Indices %in% tx2Indices)]
    N <- length(tIndices)
    if(N<1){
      #No NROY points behind this pixel, return max imp above the cutoff and 0 as the relative density of the sample
      return(c(cutoff+0.25,0))
    }
    else{
      Timps <- ImpData[tIndices,ncol(ImpData),drop=F]
      c(min(Timps),nrow(Timps)/num_points)
    }
  }
  VectOnePixel <- Vectorize(ImpOnePixel)
  VectOnePixel(1:n)
}

#This code contains functions for drawing NROY implausibility matrices. The tricky part in using it is having the list of NROY densities and Implausibility matrices in the right format. #Separate files will include code for generating one of these given emulators and data etc. There will also be an example list of lists and separate code for drawing general maps in layout form.
rain <- function(n){
  rainbow(n,st=1-0.01/12,en=2/6,alpha=1)[n:1]
}
#library(rgl)
library(shape)
rain <- function(n){
  intpalette(c("darkgreen","green","yellow","orange","red"),numcol=n)
}

breakVec <- sort(c(200,3,2.5,2,1.5,1,0.5,0))
cold1 <- function(n,fra1=0.1,fra2=0.25,fra3=0.3,n1=1){
  n4 <- round(n*fra3,0)
  n2 <- round(n*fra2,0)
  n3 <- n-n1-n2-n4
  co1 <- rep(hsv(0.68,0.25,0.55),n1)
  co2 <- rainbow(n2,st=3/6,en=3.001/6,alpha=seq(0,0.9,len=n2))
  co3 <- rainbow(n3,st=3/6,en=1,alpha=seq(0.7,0.8,len=n3))
  co4 <- rainbow(n4,st=0,en=1/6,alpha=seq(0.8,1,len=n4))
  c(co1,co2,co3,co4)
}

imp.layout01 <- function(ImpList,VarNames,VariableDensity=TRUE,newPDF=FALSE,the.title=NULL,newPNG=FALSE,newJPEG=FALSE,newEPS=FALSE,Points=NULL, ...){
  if(newPDF){
    if(is.null(the.title))
      stop("Require a Title for the PDF")
    else
      pdf(the.title,compress=FALSE, ...)
  }
  else if(newPNG){
    if(is.null(the.title))
      stop("Require a Title for the PDF")
    else
      png(the.title)
  }
  else if(newJPEG){
    if(is.null(the.title))
      stop("Require a JPEG title")
    else
      jpeg(the.title,quality=100)
  }
  else if(newEPS){
    if(is.null(the.title))
      stop("Require an eps title")
    else
      postscript(the.title)
  }
  else{
    #	dev.new()
  }
  tBig <- 0
  n <- length(VarNames)
  n1 <- n-1
  for(j in 1:n1){
    for(i in 1:(max(n1-j+1))){
      tBig <- max(tBig,max(ImpList[[j]][[i]][2,]))
    }
  }
  Mainlevs <- pretty(c(0,tBig),100)
  Maincols <- cold1(length(Mainlevs)-1,n1=1)
  mar <- par("mar")
  tmar <- mar
  tmar[4] <- tmar[2]
  tmar[2] <- 1
  tmar <- tmar/4
  tmar[4] <- 3
  tmar[1] <- 0.3
  tmar[3] <- tmar[1]
  key <- rep(1,n)
  key2 <- rep(2,n)
  numMatrix <- diag(3:(n+2))
  currentNum <- n+3
  for(i in 1:(n-1)){
    for(j in (i+1):n){
      numMatrix[i,j] <- currentNum
      currentNum <- currentNum + 1
      numMatrix[j,i] <- currentNum
      currentNum <- currentNum + 1
    }
  }
  layout(cbind(key2,numMatrix,key),widths=c(lcm(1.6),rep(0.3,n),lcm(1.6)))
  par(mar=tmar,las=1)
  plot.new()
  plot.window(c(0,1),range(Mainlevs),xaxs="i",yaxs="i")
  rect(0,Mainlevs[-length(Mainlevs)],1,Mainlevs[-1],col=Maincols)
  axis(4)
  par(mar=c(0.3,2.2,0.3,0.8))
  plot.new()
  plot.window(c(0,1),c(0,3.45),xaxs="i",yaxs="i")
  bV <- breakVec
  bV[8] <- 3.5
  rect(0,bV[-length(bV)],1,bV[-1],col=rain(7))
  axis(2)
  par(mar=rep(0.1,4),usr=c(0,1,0,1))
  for(k in 1:n){
    plot.new()
    plot.window(c(0,1),c(0,1))
    rect(0,0,1,1)
    text(0.5,0.5,labels=VarNames[k],cex=10/n)
    if(k==1){
      arrows(0,0,0,1.01,col=1,code=2,lwd=0.8,length=0.1)
      arrows(0,0,1.01,0,col=1,code=2,lwd=0.8,length=0.1)
    }
  }
  for(i in 1:(n-1)){
    for(j in (i+1):n){
      zmatod <- ImpList[[i]][[j-(i-1)-1]][2,]
      tmpn <- sqrt(length(zmatod))
      dim(zmatod) <- c(tmpn,tmpn)
      x1 <- seq(from=1,by=1,to=tmpn)/tmpn - (1/(2*tmpn))
      plot.new()
      plot.window(c(0,1),c(0,1))
      #.Internal(filledcontour(x1,x1,t(zmatod),Mainlevs,Maincols))
      .filled.contour(x1,x1,t(zmatod),Mainlevs,Maincols)
      rect(0,0,1,1)
      if(!is.null(Points)){
        points(Points[,j],Points[,i],pch=c(16:(16-length(Points[,1]))),col=1,cex=1.2)
        #points(Points[,j],Points[,i],pch=c(16,17),col=1,cex=1.2)
      }
      if(VariableDensity){
        plot.new()
        plot.window(c(0,1),c(0,1))
        levs <- pretty(range(zmatod),50)
        cols <- cold1(length(levs)-1)
        #.Internal(filledcontour(x1,x1,zmatod,levs,cols))
        .filled.contour(x1,x1,t(zmatod),levs,cols)
        rect(0,0,1,1)
      }
      else{
        zmatimp <- ImpList[[i]][[j-(i-1)-1]][1,]
        tmpn <- sqrt(length(zmatimp))
        dim(zmatimp) <- c(tmpn,tmpn)
        x1 <- seq(from=1,by=1,to=tmpn)/tmpn - (1/(2*tmpn))
        plot.new()
        #.Internal(filledcontour(x1,x1,zmatimp,breakVec,rain(100)))
        .filled.contour(x1,x1,t(zmatimp),breakVec,rain(7))
        rect(0,0,1,1)
        if(!is.null(Points)){
          points(Points[,j],Points[,i],pch=c(16:(16-length(Points[,1]))),col=1,cex=1.2)
          #points(Points[,j],Points[,i],pch=c(16,17),col=1,cex=1.2)
        }
      }
    }       
  }
}