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)
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]