@@ -0,0 +1,2336 @@
+%0A# encoding: utf-8%0A%22%22%22Routines for sampling from the Haar measures of the classical compact%0Agroups. Algorithms taken from http://arxiv.org/abs/math-ph/0609050.%0A%0ATODO Symplectic groups are missing%0A%22%22%22%0A%0Afrom __future__ import division, print_function%0Aimport numpy as np%0Afrom scipy.linalg import qr, eigvals%0A%0A%0Adef goe(dim, randn=np.random.randn):%0A    %22%22%22Returns a sample from the Gaussian orthogonal ensemble of given%0A    dimension.  (i.e. the haar measure on U(dim)).%0A%0A    :param int dim: Dimension%0A    :param randn: Function to create real N(0,1) distributed random variables.%0A        It should take the shape of the output as numpy.random.randn does%0A        (default: numpy.random.randn)%0A    %22%22%22%0A    z = randn(dim, dim)%0A    q, r = qr(z)%0A    d = np.diagonal(r)%0A    ph = d / np.abs(d)%0A    return q * ph%0A%0A%0Adef gue(dim, randn=np.random.randn):%0A    %22%22%22Returns a sample from the Gaussian unitary ensemble of given dimension.%0A    (i.e. the haar measure on U(dim)).%0A%0A    :param int dim: Dimension%0A    :param randn: Function to create real N(0,1) distributed random variables.%0A        It should take the shape of the output as numpy.random.randn does%0A        (default: numpy.random.randn)%0A    %22%22%22%0A    z = (randn(dim, dim) + 1j * randn(dim, dim)) / np.sqrt(2.0)%0A    q, r = qr(z)%0A    d = np.diagonal(r)%0A    ph = d / np.abs(d)%0A    return q * ph%0A%0A%0A#############%0A#  Tesing   #%0A#############%0ALEVEL_SPACING_DF = %7B'goe': lambda s: np.pi / 2 * s * np.exp(-np.pi / 4 * s**2),%0A                    'gue': lambda s: 32 / np.pi**2 * s**2 * np.exp(-4 / np.pi * s**2)%0A                    %7D%0A%0A%0Adef _test_ensemble(dim, ensemble, samples=1000):%0A    import matplotlib%0A    matplotlib.use('TkAgg')%0A    from matplotlib import pyplot as pl%0A    from tools.helpers import Progress%0A%0A    eigdists = %5B%5D%0A    get_sample = globals()%5Bensemble%5D%0A    for _ in Progress(xrange(samples)):%0A        u = get_sample(dim)%0A        eigs = np.sort(np.angle(eigvals(u)))%0A        eigdists += list(eigs%5B1:%5D - eigs%5B:-1%5D)%0A%0A    eigdists = np.asarray(eigdists) / (np.sum(eigdists) / len(eigdists))%0A    pl.hist(eigdists, bins=50, normed=True)%0A    dist = lambda s: 32 / np.pi**2 * s**2 * np.exp(-4 / np.pi * s**2)%0A    s = np.linspace(0, 4, 100)%0A    pl.plot(s, dist(s))%0A    pl.show()%0A%0Aif __name__ == '__main__':%0A    _test_ensemble(50, 'gue', samples=10000)%0A    _test_ensemble(50, 'goe', samples=10000)%0A

Add the random matrix module