# Jaspinder Singh

# Import necessary modules
from numpy import *
from processlibrary import *
from pylab import *
from numpy.fft import fft

# THINGS TO DO EVENTUALLY:
# 1. DEFINE A RANDOM FUNCTION
# 1a. REPLACE ALL RANDOMNESS WITH CALLS TO RANDOM FNC

def dataAppend(state):
	if state == 0:
		return -1
	else:
		return 1

print "1. Unbiased Coin\n"
print "2. Period 1\n"
print "3. Period 2\n"
print "4. Golden Mean\n"
print "5. Even Process\n"
process = int(raw_input("Which process would you like to examine? "))
flips = int(raw_input("How many iterations? "))

data = []
evenCheck = 0

# Initialize the state:
state = random.random()
if state < .5:
	state = 1
else:
	state = 0

for i in range(flips):
	if process == 1:
		state = fairCoin()
		data.append(dataAppend(state))
	elif process == 2:
		state = Period1()
		data.append(dataAppend(state))
	elif process == 3:
		state = Period2(state)
		data.append(dataAppend(state))
	elif process == 4:
		state = GoldenMean(state)
		data.append(dataAppend(state))
	elif process == 5:
		# There is an unresolved problem with the Even process
		if state == 0:
			state = fairCoin()
			if state == 0:
				data.append(-1)
			elif state == 1:
				data.append(state)
		elif state == 1:
			data.append(state)
			state = 0
	else:
		print "You didn't make a valid choice, dummy!\n"

# After the flips are stored to an array, I do a Discrete Fourier Transform
FFTed = (1.0/sqrt(flips))*fft(data)

# To get the power spectrum I simply take the absolute value of each
# term in FFTed and square it
PowerSpectrum = []
for i in range(flips):
	temp = FFTed[i]
	PowerSpectrum.append((temp*conjugate(temp)).real)

# The correlation function is defined as 
# C(n) = 1 / (N - n) sum(0 to N - n) s_m' * s_m'+n

CorrelationFunction = []

# 2-Point Correlation Function
# Useful as a means of cross checking q(n)
# 
#for n in range(flips-1):
#	temp = 0
#	for m in range(flips - n):
#		temp += (data[m] * data[m+n])
#	CorrelationFunction.append((1.0/(flips - n)) * temp)

# Fourier Analysis Correlation Function
freqRange = arange(0, 1, 1/float(flips))
dx = freqRange[1] - freqRange[0]

for n in range(flips):
	integral = 0
	for f in range(flips):
		integral += PowerSpectrum[f]*cos(2*pi*n*freqRange[f])*dx
	CorrelationFunction.append(integral)

# It is more convenient to work with:
# q(n) = 1/2 * [C(n) + 1]

q = []
for n in range(len(CorrelationFunction)):
	q.append(.5 * (CorrelationFunction[n] + 1))

print "q(1): " + str(q[1])
print "q(2): " + str(q[2])
print "q(3): " + str(q[3])
print "q(inf): " + str(average(q))

# Set up the x-axes
n = arange(0, len(CorrelationFunction), 1.0)
f = arange(0, 1.0, 1/float(flips))

# Setup to produce two plots in the same window
#
# Identify the overall plot (aka 'figure')
figure(1)
# Now set up the first graph
subplot(211)
# Give it a name
title('Correlation Function')
ylabel('q(n)')
# Now plot:
#	First, plot (n, C(n)) as a series of blue ('b') lines
plot (n, q, 'b')

# Here's the second graphic, with a new identifier.
subplot(212)
# New title
title('Power Spectrum')
xlabel('f')
ylabel('P(f)')
# And plot it with a black('k') line
plot ( f, PowerSpectrum, 'k')

# To show them
show()