docs: def construct_nonuniform_bins(x_n, nbins):
      """Construct histogram using bins of unequal size to ensure approximately equal population in each bin.
  
      ARGUMENTS
        x_n (1D numpy array of floats) - x_n[n] is data point n
  
      RETURN VALUES
        bin_left_boundary_i (1D numpy array of floats) - data in bin i will satisfy bin_left_boundary_i[i] <= x < bin_left_boundary_i[i+1]
        bin_center_i (1D numpy array of floats) - bin_center_i[i] is the center of bin i
        bin_width_i (1D numpy array of floats) - bin_width_i[i] is the width of bin i
        bin_n (1D numpy array of int32) - bin_n[n] is the bin index (in range(nbins)) of x_n[n]
      """
  
      # Determine number of samples.
      N = x_n.size
  
      # Get indices of elements of x_n sorted in order.
      sorted_indices = x_n.argsort()
  
      # Allocate storage for results.
      bin_left_boundary_i = zeros([nbins+1], float64) 
      bin_right_boundary_i = zeros([nbins+1], float64)
      bin_center_i = zeros([nbins], float64)
      bin_width_i = zeros([nbins], float64)
      bin_n = zeros([N], int32)
      
      # Determine sampled range, adding a little bit to the rightmost range to ensure no samples escape the range.
      x_min = x_n.min()
      x_max = x_n.max()
      x_max += (x_max - x_min) * 1.0e-5
  
      # Determine bin boundaries and bin assignments.
      for bin_index in range(nbins):
          # indices of first and last data points in this span
          first_index = int(float(N) / float(nbins) * float(bin_index))
          last_index = int(float(N) / float(nbins) * float(bin_index+1))
  
          # store left bin boundary
          bin_left_boundary_i[bin_index] = x_n[sorted_indices[first_index]]
  
          # store assignments
          bin_n[sorted_indices[first_index:last_index]] = bin_index                   
  
      # set rightmost boundary
      bin_left_boundary_i[nbins] = x_max
      
      # Determine bin centers and widths
      for bin_index in range(nbins):
          bin_center_i[bin_index] = (bin_left_boundary_i[bin_index] + bin_left_boundary_i[bin_index+1]) / 2.0
          bin_width_i[bin_index] = (bin_left_boundary_i[bin_index+1] - bin_left_boundary_i[bin_index])
  
      # DEBUG
  #    outfile = open('states.out', 'w')
  #    for n in range(N):
  #        outfile.write('%8f %8d\n' % (x_n[n], bin_n[n]))
  #    outfile.close()
          
      return (bin_left_boundary_i, bin_center_i, bin_width_i, bin_n)