from numpy import *
from pylab import *

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

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

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

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

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

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

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

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

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

# 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()