from numpy import *
from pylab import *

# number of bits in shift register
M = 24
# initialize register with all zeros
bits = [0]*M
# populate the last element of the register
bits[-1] = 1

# load lags data from file
lags = loadtxt("lag_values.txt", dtype = {'names': ('M', 'i'), 'formats': ('i2', 'S9')})

# extract relevant row
trow = lags[M-2]
# extract string of tap indices
data = trow[1]
# initialize taps with all zeros
taps = [0]*M

# set taps
start = 0
while True:
	# find the next comma
	end = data.find(",", start)
	if end == -1:
		# extract the last tap index
		taps[int16(data[start:len(data)]) - 1] = 1
		break
	else:
		# extract tap index
		taps[int16(data[start:end]) - 1] = 1
		# set the new start search value
		start = end + 1

# initialize powers of 2
def f1(x): return 2**x
twos = map(f1, range(M))

# to reduce correlation of LFSR, leap forward M iterations
leap = M
# define random number generator
def rand():
	for i in range(leap):
		# increment the register
		bits.insert(0, mod(dot(bits, taps), 2))
		# push all bits to the right
		del bits[-1]
	# convert from binary to decimal, uniform distribution from 0 to 1
	return float(dot(bits, twos))/(2**M - 1)

# number of walkers
w = 10
# number of steps
N = 1000

# initialize walker positions
x = [[0]*w, [int(rand() > 0.5)*2 - 1 for j1 in range(w)]]

for j2 in range(N-1):
	# compute step, equal probability right or left
	step = [int(rand() > 0.5)*2 - 1 for k1 in range(w)]
	# calculate new position
	new = [x[-1][k2] + step[k2] for k2 in range(w)]
	# update position
	x.append(new)

# initialize error bars
def f2_pls(x): return 1.5*sqrt(x)
root_pls = map(f2_pls, range(N+1))
def f2_mns(x): return -1.5*sqrt(x)
root_mns = map(f2_mns, range(N+1))

# plot walks
fig = figure()
ax = fig.add_subplot(111)
ax.plot(range(N+1), x, 'b', alpha = 0.5)
ax.plot(range(N+1), root_pls, 'k', range(N+1), root_mns, 'k')
title('Random walkers and analytical $3 \sigma$', fontsize = 17)
xlabel('step', fontsize = 20)
ylabel('displacement', fontsize = 20)
savefig('walkers.pdf')
show()