Domino shuffling algorithm

Implement the domino shuffling algorithm that we discussed in class, to generate a large random perfect matching. For reference, you can look at the paper "Generalized Domino Shuffling" by Propp, or the paper "Alternating Sign Matrices and Domino Tilings II" by Elkies, Kuperberg, Larsen, Propp.

I recommend storing perfect matchings on an order-$n$ aztec diamond tilted by 45 degrees, in a $2n$ by $2n$ array of booleans (or you can use zeros and ones). So, for instance, the two perfect matchings on the order $1$ aztec diamond would then be

10   01
01   10  

and the eight perfect matchings on the order 2 aztec diamond would be

1010  1001  1010  1001  1001  0101  1001  0101  
0001  0010  0001  0010  0100  1000  0100  1000 
1000  1000  0100  0100  0010  0010  0001  0001  
0101  0101  1001  1001  1001  1001  1010  1010  

With these conventions: "northbound" edges are located in cell $(2i, 2j)$, westbound in $(2i+1,2j)$, etc (I am using matrix coordinates).

Aztec Diamond Picture

Write code that draws a picture of a perfect matching on an Aztec Diamond.

I suggest actually drawing the "domino tiling" description, rather than the perfect matching description, so that you can show a bigger example; use four different colors to represent north, south, east, west dominos.

Note: history's solution is to use red, green, blue, yellow for the four different types of dominos. You can do this if you want. However it might be better to pick almost any other colors; people with color blindness (including some of the authors on the first Aztec Diamond paper) can't tell some of these colors apart.

Here's an order-2 example which I did by hand. Note: this uses cartesian coordinates (x,y), and "North" is the direction (-1, -1) - I don't particularly care which conventions you use for your picture, and I'll leave it to you to decide which of the above tilings it is.

import matplotlib.pyplot as plt
north_color = '#053363'
south_color = '#175b91'
east_color = '#40b0bf'
west_color = '#fcb735'
north_dominos = [[(0,2),(1,3),(3,1),(2,0)], [(2,4),(3,5),(5,3),(4,2)]]
east_dominos = [[(0,4),(1,3),(3,5),(2,6)]]
west_dominos = [[(4,0),(3,1),(5,3),(6,2)]]
south_dominos = [[(1,3),(2,4),(4,2),(3,1)], [(3,5),(4,6),(6,4),(5,3)]]

for d in north_dominos:
    plt.fill(*zip(*d), color=north_color)
for d in east_dominos:
    plt.fill(*zip(*d), color=east_color)
for d in west_dominos:
    plt.fill(*zip(*d), color=west_color)
for d in south_dominos:
    plt.fill(*zip(*d), color=south_color)
plt.axis('equal')
plt.show()

Kasteleyn Matrix

(a) Compute a Kasteleyn matrix for the order 3 aztec diamond graph, by choosing some explicit ordering of the shaded and unshaded vertices. It should be some 24 by 24 matrix. Verify that its determinant is $\pm 64$.

(b) There are 36 edges in the order 3 aztec diamond graph, as we learned in Problem 1. For each, compute the probability that the edge is in a randomly chosen perfect matching. Note: 36 is much less than $24^2$ - so definitely not all the entries in $K^{-1}$ correspond to such probabilities! Write the answers down in a 6 by 6 matrix.