import numpy
import matplotlib.pyplot as plt

r1, r2, r3, r4, r5 = numpy.load('random.npy')
N = r1.size

# PDF Histo
print "PDF histogram"
plt.figure()
plt.subplot(221)
plt.hist(r1, 100, histtype='step', normed=True)
plt.hist(r2, 100, histtype='step', normed=True)
plt.hist(r3, 100, histtype='step', normed=True)
plt.hist(r4, 100, histtype='step', normed=True)
plt.hist(r5, 100, histtype='step', normed=True)

# Chi square
print "chi square"
def chisquare(rs, sample_size=10):
    N = rs.size
    n = N/sample_size
    means = numpy.empty(n)
    for i in range(n):
        means[i] = rs[i*sample_size:(i+1)*sample_size].mean()
    return means

plt.subplot(222)
def gaussian(x, mu, sigma):
    return 1./(sigma * numpy.sqrt(2.0*numpy.pi)) * numpy.exp(-(x-mu)**2 / (2.0*sigma**2))

x = numpy.linspace(0.0, 1.0, 100)
plt.plot(x, gaussian(x, 0.5, 1.0/numpy.sqrt(120.0)), linewidth=2, color='r')

plt.hist(chisquare(r1), 100, histtype='step', normed=True)
plt.hist(chisquare(r2), 100, histtype='step', normed=True)
plt.hist(chisquare(r3), 100, histtype='step', normed=True)
plt.hist(chisquare(r4), 100, histtype='step', normed=True)
plt.hist(chisquare(r5), 100, histtype='step', normed=True)

plt.legend(("expected", "r1", "r2", "r3", "r4", "r5"))

# Fourier analysis
print "Fourier analysis"
def fourier(rs):
    return numpy.abs(numpy.fft.rfft(rs))/rs.size

plt.subplot(223)
plt.semilogy(fourier(r1))
plt.semilogy(fourier(r2))
plt.semilogy(fourier(r3))
plt.semilogy(fourier(r4))
plt.semilogy(fourier(r5))

# Triplet correlation
print "triplet correlation"
def triplet_corr_error(rs, p=250):
    corr = numpy.empty(p-1)

    N = rs.size
    for k in range(1,p):
        c = rs[p:]*rs[:-p]*rs[p-k:-k]
        corr[k-1] = c.sum() / (N-p)
        
    error = abs(corr - 0.5**3)

    return error

plt.subplot(224)
plt.semilogy(triplet_corr_error(r1), '-')
plt.semilogy(triplet_corr_error(r2), '-')
plt.semilogy(triplet_corr_error(r3), '-')
plt.semilogy(triplet_corr_error(r4), '-')
plt.semilogy(triplet_corr_error(r5), '-')
plt.show()