def hexbin(data, x_bin_size, y_bin_size):
    """Bin 2D data into hexagons.

       Tessalate the data plane with hexagons the centers of which are
       x_bin_size apart along the x axis and y_bin_size along the y axis,
       i.e., the distance between neighbouring centers is
       sqrt(xbin_size**2 + ybin_size**2).

       Returns a dictionary with the keys designating 2-dimensional indices
       of the hexagons and the values being lists of x- and y-coordinate of
       the centroid of the data falling into the given hexagon, and the
       number of data points in the hexagon, i.e., [cx, cy, n].
       """
    hexbinner = {}
    for x, y in data:
        xbin = x // x_bin_size
        ybin = y // y_bin_size
        if (xbin + ybin) % 2:
            if (x - x_bin_size * xbin) ** 2 + (y - y_bin_size * (ybin + 1)
                ) ** 2 < (x - x_bin_size * (xbin + 1)) ** 2 + (y - 
                y_bin_size * ybin) ** 2:
                hdata = hexbinner.setdefault((xbin, ybin + 1), [0, 0, 0])
            else:
                hdata = hexbinner.setdefault((xbin + 1, ybin), [0, 0, 0])
        elif (x - x_bin_size * xbin) ** 2 + (y - y_bin_size * ybin) ** 2 < (x -
            x_bin_size * (xbin + 1)) ** 2 + (y - y_bin_size * (ybin + 1)) ** 2:
            hdata = hexbinner.setdefault((xbin, ybin), [0, 0, 0])
        else:
            hdata = hexbinner.setdefault((xbin + 1, ybin + 1), [0, 0, 0])
        hdata[0] = hdata[0] + x
        hdata[1] = hdata[1] + y
        hdata[2] += 1
    for hdata in hexbinner.values():
        hdata[0] /= hdata[2]
        hdata[1] /= hdata[2]
    return hexbinner
