multicoin.py

#!/usr/bin/env python3.1
#
#	multicoin.py ... Quantum walk with multiple biased coins
#	
#	Copyright (C) 2010  KennyTM~
#	
#	This program is free software: you can redistribute it and/or modify
#	it under the terms of the GNU General Public License as published by
#	the Free Software Foundation, either version 3 of the License, or
#	(at your option) any later version.
#	
#	This program is distributed in the hope that it will be useful,
#	but WITHOUT ANY WARRANTY; without even the implied warranty of
#	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#	GNU General Public License for more details.
#	
#	You should have received a copy of the GNU General Public License
#	along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

import random
from math import *
from itertools import chain

# Compute theta from s(M-1) .. s1
def findThetaUnbiased(history):
	return pi/4
	
def findThetaBiased(history):
	if history == (-1, -1):
		return 42 * pi/180
	else:
		return pi/4
	
findTheta = findThetaBiased

# Compute the probabilities of flip.
def C(coins, position, amplitude):
	global findTheta
	history = coins[:-1]
	theta = findTheta(history)
	
	s0_up = cos(theta) * amplitude
	s0_down = sin(theta) * 1j * amplitude
	
	if coins[-1] == -1:	# spin-down, flip s0_0 and s0_1
		(s0_up, s0_down) = (s0_down, s0_up)
	
	return [((history + (1,), position), s0_up), ((history + (-1,), position), s0_down)]

# Perform translation
def S(coins, position, amplitude):
	coin = coins[-1]
	position += coin
	return [((coins, position), amplitude)]

# Record history
def R(coins, position, amplitude):
	newCoins = coins[1:] + (coins[0],)
	return [((newCoins, position), amplitude)]

# accumulate
def accumBy(lst):
	d = {}
	for k, v in lst:
		if k in d:
			d[k].append(v)
		else:
			d[k] = [v]
	return {k: fsum(v) for k, v in d.items()}

def accumByComplex(lst):
	(re_part, im_part) = zip( *( ((k, v.real), (k, v.imag)) for k, v in lst) )
	dr = accumBy(re_part)
	di = accumBy(im_part)
	for k, v in di.items():
		dr[k] += v * 1j
	return dr

# Apply an operator to a superposition of states
def evolveWith(states, operator):
	newStateList = chain(*(operator(coins, position, amplitude) for ((coins, position), amplitude) in states.items()))
	return accumByComplex(newStateList)

# Evolution operator
def U(states):
	return evolveWith(evolveWith(evolveWith(states, C), S), R)

# Find position distribution.
def posDistr(states):
	return accumBy((k, v.real*v.real + v.imag*v.imag) for (_, k), v in states.items())

def meanX(xdistr):
	s = fsum(p*q for p, q in xdistr.items())
	n = fsum(xdistr.values())
	return s/n

# actual evolution.
states = {((1,1,1),0): 1.0, ((-1,-1,-1),0): -1.0}
for i in range(250):
	findTheta = [findThetaUnbiased, findThetaUnbiased, findThetaBiased][i % 3]	# AAB

	states = U(states)
	print (i, "\t", meanX(posDistr(states)))