for (i in 1:n)
  {
    ## compute evaluation points 
    eval.x <- gridx[[1]][grid.pts$xmin[i,1]:grid.pts$xmax[i,1]]
    eval.y <- gridx[[2]][grid.pts$xmin[i,2]:grid.pts$xmax[i,2]]
    eval.x.ind <- c(grid.pts$xmin[i,1]:grid.pts$xmax[i,1])
    eval.y.ind <- c(grid.pts$xmin[i,2]:grid.pts$xmax[i,2])
    eval.x.len <- length(eval.x)
    eval.pts <- permute(list(eval.x, eval.y))
    fhat <- dmvnorm(eval.pts, x[i,], H)
    
    ## place vector of density estimate values `fhat' onto grid 'fhat.grid' 
    for (j in 1:length(eval.y))
      fhat.grid[eval.x.ind, eval.y.ind[j]] <- 
        fhat.grid[eval.x.ind, eval.y.ind[j]] + 
          fhat[((j-1) * eval.x.len + 1):(j * eval.x.len)]
  }
  
  fhat.grid <- fhat.grid/n
  gridx1 <- list(gridx[[1]], gridx[[2]]) 
  
  fhat.list <- list(x=x, eval.points=gridx1, estimate=fhat.grid, H=H, gridtype=gridx$gridtype)
  
  return(fhat.list)
}


###############################################################################
# Trivariate kernel density estimate using normal kernels, evaluated over grid
#
# Parameters
# x - data points
# H - bandwidth matrix
# gridsize - number of interval points in grid
# supp - effective support of kernel
#
# Returns
# list with fields
# x - data points
# eval.points - points that KDE is evaluated at
# estimate - KDE evaluated at eval.points 
# H - bandwidth matrix 
###############################################################################


kde.grid.3d <- function(x, H, gridsize, supp, gridx=NULL, grid.pts=NULL, xmin, xmax, gridtype)
{
  ## initialise grid 
  n <- nrow(x)

  if (is.null(gridx))
    gridx <- make.grid.ks(x, matrix.sqrt(H), tol=supp, gridsize=gridsize, xmin=xmin, xmax=xmax, gridtype=gridtype) 
  suppx <- make.supp(x, matrix.sqrt(H), tol=supp)

  if (is.null(grid.pts))
    grid.pts <- find.gridpts(gridx, suppx)    
  fhat.grid <- array(0, dim=c(length(gridx[[1]]), length(gridx[[2]]), 
               length(gridx[[3]])))
  
  for (i in 1:n)
  {
    ## compute evaluation points
    eval.x <- gridx[[1]][grid.pts$xmin[i,1]:grid.pts$xmax[i,1]]
    eval.y <- gridx[[2]][grid.pts$xmin[i,2]:grid.pts$xmax[i,2]]
    eval.z <- gridx[[3]][grid.pts$xmin[i,3]:grid.pts$xmax[i,3]]
    eval.x.ind <- c(grid.pts$xmin[i,1]:grid.pts$xmax[i,1])
    eval.y.ind <- c(grid.pts$xmin[i,2]:grid.pts$xmax[i,2])
    eval.z.ind <- c(grid.pts$xmin[i,3]:grid.pts$xmax[i,3])
    eval.x.len <- length(eval.x)
    eval.pts <- permute(list(eval.x, eval.y))
   
    ## place vector of density estimate values `fhat' onto grid 'fhat.grid' 

    for (k in 1:length(eval.z))
    {
      fhat <- dmvnorm(cbind(eval.pts, eval.z[k]), x[i,], H)
      for (j in 1:length(eval.y))
        fhat.grid[eval.x.ind,eval.y.ind[j], eval.z.ind[k]] <- 
          fhat.grid[eval.x.ind, eval.y.ind[j], eval.z.ind[k]] + 
            fhat[((j-1) * eval.x.len + 1):(j * eval.x.len)]
     }
  }
  
  fhat.grid <- fhat.grid/n

  gridx1 <- list(gridx[[1]], gridx[[2]], gridx[[3]]) 
  fhat.list <- list(x=x, eval.points=gridx1, estimate=fhat.grid, H=H, gridtype=gridx$gridtype)

  return(fhat.list)
}





###############################################################################
# Multivariate kernel density estimate using normal kernels,
# evaluated at each sample point
#
# Parameters
# x - data points
# H - bandwidth matrix
# eval.points - points where to evaluate density estimate
#
# Returns
# list with fields
# x - data points
# eval.points - points that KDE is evaluated at
# estimate - KDE evaluated at eval.points 
# H - bandwidth matrix 
###############################################################################

kde.points <- function(x, H, eval.points) 
{
  n <- nrow(x)
  Hs <- numeric(0)
  for (i in 1:n)
    Hs <- rbind(Hs, H)
  
  fhat <- dmvnorm.mixt(x=eval.points, mus=x, Sigmas=Hs, props=rep(1, n)/n)

  return(list(x=x, eval.points=eval.points, estimate=fhat, H=H))
}

kde.points.1d <- function(x, h, eval.points, positive=FALSE, adj.positive) 
{
  n <- length(x)

  if (positive)
  {
    if (missing(adj.positive)) adj.positive <- abs(min(x))
    y <- log(x + adj.positive)  ## transform positive data x to real line
    eval.pointsy <- log(eval.points + adj.positive)
  }
  else
  {
    y <- x
    eval.pointsy <- eval.points
  }
  
  fhat <- 0
  for (k in 1:n)
    fhat <- fhat + dnorm(x=eval.pointsy, mean=y[k], sd=h)
  fhat <- fhat/n
  if (positive)
    fhat <- fhat/(eval.points + adj.positive) ##fhat/exp(eval.pointsy)
  
  return(list(x=x, eval.points=eval.points, estimate=fhat, h=h, H=h^2))
}



###############################################################################
# Display kernel density estimate
#
# Parameters
# fhat - output from call to `kde'
###############################################################################

plot.kde <- function(x, drawpoints=FALSE, ...)
{ 
  fhat <- x

  if (is.vector(fhat$x))
    plotkde.1d.v2(fhat, drawpoints=drawpoints, ...)
  else
  {
    d <- ncol(fhat$x)

    if (d==2) 
    {
      plotret <- plotkde.2d.v2(fhat, drawpoints=drawpoints, ...)
      invisible(plotret)
    }
    else if (d==3)
    {
      plotkde.3d(fhat, drawpoints=drawpoints, ...)
       invisible()
    }
    else 
      stop ("Plot function only available for 1, 2 or 3-dimensional data")
  }
}

plotkde.1d.v2 <- function(fhat, xlab, ylab="Density function", add=FALSE,
  drawpoints=TRUE, ptcol="blue", jitter=FALSE, ...) #col="black", ...)
{
  if (missing(xlab)) xlab <- fhat$names
  
  if (add)
    lines(fhat$eval.points, fhat$estimate, xlab=xlab, ylab=ylab, ...)
  else
    plot(fhat$eval.points, fhat$estimate, type="l", xlab=xlab, ylab=ylab, ...) 

  if (drawpoints)
    if (jitter)
      rug(jitter(fhat$x), col=ptcol)
    else
      rug(fhat$x, col=ptcol)

 
  
}


###############################################################################
# Display bivariate kernel density estimate
#
# Parameters 
# fhat - output from 'kde.grid'
# display - "persp" - perspective plot
#         - "slice" - contour plot
#         - "image" image plot
# cont - vector of contours to be plotted
###############################################################################

plotkde.2d.v2 <- function(fhat, display="slice", cont=c(25,50,75), abs.cont,
    xlab, ylab, zlab="Density function", cex=1, pch=1,  
    add=FALSE, drawpoints=TRUE, drawlabels=TRUE, theta=-30, phi=40, d=4,
    ptcol="blue", ...) #shade=0.75, border=NA, persp.col="grey", ...)
{
  disp1 <- substr(display,1,1)
  if (!is.list(fhat$eval.points))
    stop("Need a grid of density estimates")

  if (missing(xlab)) xlab <- fhat$names[1]
  if (missing(ylab)) ylab <- fhat$names[2]

  ##eval1 <- fhat$eval.points[[1]]
  ##eval2 <- fhat$eval.points[[2]]
  
  ## perspective/wireframe plot
  if (disp1=="p")
    plotret <- persp(fhat$eval.points[[1]], fhat$eval.points[[2]], fhat$estimate,
          theta=theta, phi=phi, d=d, xlab=xlab, ylab=ylab, zlab=zlab, ...)
          ##shade=shade, border=border, col=persp.col, ...)
  
  else if (disp1=="s")
  {
    if (!add)
      plot(fhat$x[,1], fhat$x[,2], type="n",xlab=xlab, ylab=ylab, ...)
      
    ## compute contours
    if (missing(abs.cont))
      hts <- contourLevels(fhat, prob=(100-cont)/100)
    else if (is.null(abs.cont))
      hts <- contourLevels(fhat, n.pretty=5)  
    else
      hts <- abs.cont 
   
    ## draw contours         
    for (i in 1:length(hts)) 
    {
      if (missing(abs.cont)) scale <- cont[i]/hts[i]
      else scale <- 1

      contour(fhat$eval.points[[1]], fhat$eval.points[[2]], 
              fhat$estimate*scale, level=hts[i]*scale, add=TRUE, 
              drawlabels=drawlabels, ...)
    }
 
    ## add points 
    if (drawpoints)
      points(fhat$x[,1], fhat$x[,2], col=ptcol, cex=cex, pch=pch)
  }
  ## image plot
  else if (disp1=="i")
  {
    image(fhat$eval.points[[1]], fhat$eval.points[[2]], fhat$estimate, 
            xlab=xlab, ylab=ylab, add=add, ...)
    box()
  }
  else if (disp1=="f")
    filled.contour(fhat$eval.points[[1]], fhat$eval.points[[2]], fhat$estimate, 
                   xlab=xlab, ylab=ylab, ...)

  if (disp1=="p")  invisible(plotret)
  else invisible()
}
  


###############################################################################
## Display trivariate kernel density estimate
###############################################################################


plotkde.3d <- function(fhat, cont=c(25,50,75), abs.cont, colors,
                       alphavec, size=3, ptcol="blue", add=FALSE, 
                       xlab, ylab, zlab, drawpoints=FALSE, alpha=1, ...)

{
  require(rgl)
  require(misc3d)
 
  if (missing(abs.cont))
    hts <- contourLevels(fhat, prob=(100-cont)/100)
  else
    hts <- abs.cont 
  nc <- length(hts)
  
  if (missing(colors))
    colors <- rev(heat.colors(nc))

  if (missing(xlab)) xlab <- fhat$names[1]
  if (missing(ylab)) ylab <- fhat$names[2]
  if (missing(zlab)) zlab <- fhat$names[3]
    #if (is.null(x.names)) zlab <- "z" else zlab <- x.names[3]
  
  if (missing(alphavec))
    alphavec <- seq(0.1,0.5,length=nc)

  
  ##if (!add) clear3d()
  
  bg3d(col="white")
  
  if (drawpoints)
    plot3d(fhat$x[,1],fhat$x[,2],fhat$x[,3], size=size, col=ptcol, alpha=alpha, xlab=xlab, ylab=ylab, zlab=zlab, add=add, ...)
  else
    plot3d(fhat$x[,1],fhat$x[,2],fhat$x[,3], type="n", xlab=xlab, ylab=ylab, zlab=zlab, add=add, ...)
  
  for (i in 1:nc)
    contour3d(fhat$estimate, level=hts[nc-i+1], x=fhat$eval.points[[1]],
              y=fhat$eval.points[[2]], z=fhat$eval.points[[3]], add=TRUE,
              color=colors[i], alpha=alphavec[i], ...)
}




###############################################################################
### Contour levels 
###############################################################################

## create S3 generic 
contourLevels <- function(x, ...){  
  UseMethod("contourLevels")  
}   

contourLevels.kde <- function(x, prob, cont, nlevels=5, ...)
{ 
  fhat <- x
  if (is.vector(fhat$x))
  {
    d <- 1; n <- length(fhat$x); H <- as.matrix(fhat$H)
    bgridsize <- length(fhat$estimate)
  }
  else
  {
    d <- ncol(fhat$x); n <-nrow(fhat$x); H <- fhat$H
    bgridsize <- dim(fhat$estimate)
  }

  ## for large sample sizes, use binned approx. 
  if (n >= 5e3 & d <= 4 & fhat$binned)
  {
    bin.par <- binning(fhat$x, bgridsize=bgridsize, H=H, supp=3.7)
    dobs <- rep(fhat$estimate, round(bin.par$counts,0))
    dobs <- dobs[dobs>0]
  }
  else
    dobs <- kde(x=fhat$x, H=fhat$H, eval.points=fhat$x)$estimate 
  
  if (missing(prob) & missing(cont))
    hts <- pretty(dobs, n=nlevels) 
  
  if (!missing(prob) & missing(cont))
    hts <- quantile(dobs, prob=prob)
  
  if (missing(prob) & !missing(cont))
    hts <- quantile(dobs, prob=(100-cont)/100)
  
  return(hts)
}