generate sorted random numbers

sorted_random.py

"""
Generate sorted random numbers.
method from: http://repository.cmu.edu/cgi/viewcontent.cgi?article=3483&context=compsci

Has same interface as random module except it takes N as the first argument
for the number of numbers to generate.
"""
from random import random as rand, seed

def random(N):
    """
    Generate N sorted uniform random numbers between 0 and 1
    >>> seed(42)
    >>> r = list(random(10))
    >>> all(r0 < r1 for r0, r1 in zip(r, r[1:]))
    True
    >>> len(r)
    10
    """
    curmax = 1.0
    for I in range(N, 0, -1):
        curmax *= rand() ** (1.0 / float(I))
        yield (1.0 - curmax)


def randint(N, imin, imax, include_end=False):
    """
    Generate N sorted uniform random numbers between imin and imax.
    >>> seed(42)
    >>> list(randint(6, 100, 200))
    [107, 155, 167, 180, 183, 188]
    """
    # include end-point like random.randint()
    span = imax - imin + int(bool(include_end))
    for r in random(N):
        yield int(imin + r * span)

if "__main__" == __name__:
    import doctest
    doctest.testmod()

    import numpy as np
    import seaborn as sns
    from matplotlib import pyplot as plt

    seed()
    N = 10000
    chrom_size = 500000
    starts = np.array(list(randint(N, 0, chrom_size)))

    assert np.all(starts[1:] >= starts[:-1])
    print starts.min(), starts.max()

    plt.hist(starts, 15)
    plt.show()

Once we get a breather, we should sit down and implement this in bedtools.