def take2D(histogram, x, y, bins_x, bins_y):
    """
    Take the value from a two-dimensional histogram from the bin corresponding to (x, y).

    Parameters:
    -----------
    histogram : The values in the histogram (n,m) (ADW: is this ordering right?)
    x : the x-value to take from the hist
    y : the y-value to take from the hist
    bins_x : the xbin edges, including upper edge (n-dim)
    bins_y : the ybin edges, including upper edge (m-dim)
    """
    histogram = np.array(histogram)
    
    if np.isscalar(x):
        x = [x]
    if np.isscalar(y):
        y = [y]

    bins_x[-1] += 1.e-10 * (bins_x[-1] - bins_x[-2]) # Numerical stability
    bins_y[-1] += 1.e-10 * (bins_y[-1] - bins_y[-2])

    #return np.take(histogram, (histogram.shape[1] * (np.digitize(y, bins_y) - 1)) + (np.digitize(x, bins_x) - 1))

    # Return np.nan for entries which are outside the binning range on either axis
    index = (histogram.shape[1] * (np.digitize(y, bins_y) - 1)) + (np.digitize(x, bins_x) - 1)
    index_clipped = np.clip(index, 0, (histogram.shape[0] * histogram.shape[1]) - 1)
    val = np.take(histogram, index_clipped)

    outlier_x = np.logical_or(x < bins_x[0], x > bins_x[-1])
    outlier_y = np.logical_or(y < bins_y[0], y > bins_y[-1])
    outlier = np.logical_or(outlier_x, outlier_y)
    val[outlier] = np.nan

    return val