#-- required library hexbin --
 source("http://bioconductor.org/biocLite.R")
 biocLite("hexbin")
 library("hexbin")
#--------------------

hexplot=
function (x, y, style = "lattice",     # either of the two: "lattice", "colorscale"
          xlim = range(x), ylim = range(y),
          xbins = 30, ybins = 30, shape = 1, 
          legend = TRUE,
          minarea = 0.04, maxarea = 1, 
          mincnt = 1, maxcnt = 1,      # default is max(hex@count),
          leg.cnt = 0,
          colramp = function(n) LinGray(n, beg = 90, end = 15),            
          xlab = "", ylab = "", 
          xaxt = c("s", "n"), yaxt = c("s", "n"), main = "")
{
    if (!is.vector(x)) 
        stop("x must be a vector of x-coordinates")
    if (!is.vector(y)) 
        stop("y must be a vector of y-coordinates")
    if (length(x) != length(y)) 
        stop("x and y coordinates must be vectors of the same length")
    if (!is.logical(legend)) 
        stop("Legend is a logical parameter")
    if ((xbins <= 0) | !is.numeric(xbins))
        stop("xbins must be a positive number")
    if ((ybins <= 0) | !is.numeric(ybins))
        stop("ybins must be a positive number")
    if (minarea < 0) 
        stop("Minimum area must be non-negative")
    if (maxarea > 1) 
        warning("Maximum area should be <= 1, this leads to overlapping hexagons")
    if (minarea > maxarea) 
        stop("Minarea must be <= maxarea")
    style.types = c("lattice", "colorscale")
    if (!is.element(style, style.types))
        stop("Style not defined") 
    if (!is.vector(leg.cnt))
        stop("Leg.cnt must be a vector of the count cuts")
        else leg.cnt = leg.cnt[order(leg.cnt)]
    if (!is.vector(leg.cnt))
        stop("Leg.cnt must be a vector of the count cuts")
        else leg.cnt = leg.cnt[order(leg.cnt)]

    xaxt <- match.arg(xaxt)
    yaxt <- match.arg(yaxt)

good = (x>=xlim[1]) & (x<=xlim[2]) & (y>=ylim[1]) & (y<=ylim[2])
x = x[good]
y = y[good]
h = ybins*shape / xbins
sh = ybins / xbins
hex = hexbin (x, y, xbins = xbins, shape = sh, xbnds = xlim, ybnds = ylim) 
 if (maxcnt == 1) 
     maxcnt = max(hex@count)
 if (length(leg.cnt) == 1)
     leg.cnt = pretty(mincnt:maxcnt, n = min(14,(maxcnt - mincnt + 1)))

B = L = 4 
T = R = 2
 if (h > 1) {
  B = B*h
  T = T*h
 } else {
  R = R/h
  L = L/h
 }
    min = max(mincnt, min(hex@count))
    low = max(which(leg.cnt <= min), 1)        # lower & upper cnts used
    up = max(which(leg.cnt < min(maxcnt, max(hex@count))), low)
    bnds = leg.cnt[low : (up + 1)]
    bnds[1] = min
    plus = ifelse(is.na(leg.cnt[up + 1]), 1, 0)


if (legend == TRUE) {        ### legend ###

layout(mat = matrix(c(2,1),nrow=1,byrow=TRUE), widths = c(10,3), heights = 10*h, respect = TRUE)
  ysize = diff(ylim)                         
  xsize = diff(xlim) * (3/(8.2 - 1*h - 2*(1-h)*(h<1) ))    
par(mar = c(B,0,T,0))
plot(c(0,xsize), c(0,ysize), bty = "n", type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "")  
             
    n = length(bnds)
    spacing <- ysize/(n + 3)     
    midx <- xsize/3
    textx <- (2 * xsize)/3

 ### >>> colorscale <<<

    switch(style, colorscale = {              
     dy = (min(xsize/2, (sqrt(3) * spacing))) / (4 * sqrt(3)) 
     dx = dy * (diff(hex@xbnds)/diff(hex@ybnds)) * hex@shape * sqrt(3)
     hexC = hexcoords(dx, dy, n = 1, sep = NULL)

        pen <- colramp(n)
        tx <- hexC$x + midx
        for (i in seq(length = n - 2)) {
            polygon(tx, hexC$y + i * spacing, col = pen[i])                             
            text(labels = as.character(bnds[i]), x = textx, y = (i - 0.5) * spacing)     
        }
        polygon(tx, hexC$y + (n-1) * spacing, col = pen[n-1])                              
        if (plus == 1)
          text(labels = paste(as.character(bnds[n-1]),"+"), x = textx, y = (n - 1.5) * spacing)        
        if (plus == 0) { 
          text(labels = as.character(bnds[n-1]), x = textx, y = (n - 1.5) * spacing)  
          text(labels = as.character(bnds[n]), x = textx, y = (n - 0.5) * spacing) 
        }       
        text(labels = "Counts", x = xsize/2, y = (n+1) * spacing)     
      
 ### >>> lattice <<<
        
    }, lattice = {                            
     sx <- hex@xbins/diff(hex@xbnds)
     sy <- (hex@xbins * hex@shape)/diff(hex@ybnds)
     inner <- 0.5
     outer <- (2 * inner)/sqrt(3)
     dx <- inner/sx
     dy <- outer/(2 * sy)
     hexC <- hexcoords(dx, dy, n = 1, sep = NULL)

       pen = 1
       nc = length(bnds)          
       bnds[nc] = ifelse(is.na(bnds[nc]), min(maxcnt, max(hex@count)), bnds[nc])
       bnds.range = bnds[nc] - bnds[1]
       if (bnds.range > 0) {
        rbnds = (bnds - bnds[1]) / bnds.range
        } else 
        rbnds = c(0, 1)

        radius <- sqrt(minarea + (maxarea - minarea) * rbnds)
        nc6 = rep(6, nc)
        tx <- rep.int(hexC$x, nc) * rep.int(radius, nc6) + midx
        ty <- rep.int(hexC$y, nc) * rep.int(radius, nc6) 

       for (i in seq(length = nc - 2)) {
            polygon(tx[(6*(i-1) + 1) : (6*i)], ty[(6*(i-1) + 1) : (6*i)] + i * spacing, col = pen)      
            text(labels = as.character(bnds[i]), x = textx, y = (i - 0.5) * spacing)   
            }
       polygon(tx[(6*(nc-2) + 1) : (6*(nc-1))], ty[(6*(nc-2) + 1) : (6*(nc-1))] + (nc-1) * spacing, col = pen)                           
        if (plus == 1)
           text(labels = paste(as.character(bnds[n-1]),"+"), x = textx, y = (n - 1.5) * spacing)        
        if (plus == 0) { 
           text(labels = as.character(bnds[n-1]), x = textx, y = (n - 1.5) * spacing)  
           text(labels = as.character(bnds[n]), x = textx, y = (n - 0.5) * spacing) 
       }       
       text(labels = "Counts", x = xsize/2, y = (n + 1) * spacing)                    
     })
} 
                         ### hexagons ###

par(mar = c(B,L,T,R))  
plot(0, 0, type = "n", xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, 
    xaxt = xaxt, yaxt = yaxt, main = main)
     
    cnt <- hex@count
    xbins <- hex@xbins  # *
    shape <- hex@shape
    tmp <- hcell2xy(hex)   # Computes x and y coordinates from hexagon cell id's 
    good <- mincnt <= cnt & cnt <= maxcnt
    xnew <- tmp$x[good]
    ynew <- tmp$y[good]
    cnt <- cnt[good]     # data cleanin'
    rcnt <- {            # relative counts
        if (maxcnt == mincnt) 
            rep.int(1, length(cnt))
        else (cnt - mincnt)/(maxcnt - mincnt)                 
        }
    
 ### >>> colorscsale <<<

    switch(style, colorscale = {
      nc = length(bnds)                     # number of categories + 1
      bnds[nc] = ifelse(is.na(bnds[nc]), min(maxcnt, max(hex@count)), bnds[nc])
      bnds.range = bnds[nc] - bnds[1]
      if (bnds.range > 0) {
       rbnds = (bnds - bnds[1]) / bnds.range
       rbnds[1] <- rbnds[1] - 1e-06
       rbnds[nc] <- rbnds[nc] + 1e-06
       colgrp <- cut(rcnt, rbnds, labels = FALSE)
       if (any(is.na(colgrp))) 
          colgrp <- ifelse(is.na(colgrp), 0, colgrp)
       clrs <- colramp(nc)
       pen <- clrs[colgrp]
      } else
       pen = rep(colramp(2)[1], length(rcnt))
    })

    sx <- xbins/diff(hex@xbnds)
    sy <- (xbins * shape)/diff(hex@ybnds)
    inner <- 0.5
    outer <- (2 * inner)/sqrt(3)
    dx <- inner/sx
    dy <- outer/(2 * sy)
    rad <- sqrt(dx^2 + dy^2)
    hexC <- hexcoords(dx, dy, sep = NULL)   # function for computing hexagons' coord.s; == 6 values
                                            # dx and dy = horizontal and vertical width of the hexagon(s).
    switch(style, colorscale = {
       hp = hexpolygon(xnew, ynew, hexC)     
       newvx = as.vector(rbind(matrix(hp$x, nrow=6), rep(NA, length(hp$x)/6)))
       newvy = as.vector(rbind(matrix(hp$y, nrow=6), rep(NA, length(hp$y)/6)))
       hp$x = newvx 
       hp$y = newvy
       polygon(hp, col = pen, border = 0)    
    })

 ### >>>  lattice  <<<

    area <- pmin(minarea + rcnt * (maxarea - minarea), maxarea)    # doesn't need to use pmin
    radius <- sqrt(area)
    n <- length(radius)
    n6 <- rep.int(6:6, n)
    pltx <- rep.int(hexC$x, n) * rep.int(radius, n6) + rep.int(xnew, n6)
    plty <- rep.int(hexC$y, n) * rep.int(radius, n6) + rep.int(ynew, n6)

    switch(style, lattice = {
           newpltx = as.vector(rbind(matrix(pltx, nrow = 6), rep(NA, length(pltx)/6)))
           newplty = as.vector(rbind(matrix(plty, nrow = 6), rep(NA, length(plty)/6)))
           polygon(newpltx, newplty, col = 1)  
    })
}

### ----------------- end -----------------------