docs: 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