I'm trying to write a program (with javascript or python) that samples a random lozenge tiling of a hexagon with Propp - Wilson's coupling from the past algorithm. I'm quite clear of the framework of the algorithm, but I don't know which is the most efficient way to encode a tiling, or whether different encoding methods would result in non-uniform samplers. 


I have at least 3 possible way to encode a lozenge tiling of a hexagon:


1. View it as a plane partition, i.e. piling cubes in a 3D room. This is quite intuitive, but not easy to implement the add - remove step.
2. View it as a non-intersecting path system. Each path is further represented by a 0-1 array.
3. encode it as a interlacing array (which I have not understand it yet). But I have seen at least 3 people that mentioned this approach. For example [here](http://faculty.virginia.edu/petrov//blog/2015/02/18/Glauber-Simulation/)


So my question is: which is the best way to encode a tiling, or can anyone explain the 3rd approach?



implementing Propp-Wilson's Coupling From the Past on Lozenge Tilings of a Hexagon
You should represent it as a plane partition. But you shouldn't represent a plane partition as a 3D array of 0/1 for absent/present. Most people think of a partition as a set of natural numbers, not a Ferrers diagram. The Ferrers diagram has more symmetry, but usually the set (or decreasing sequence) of numbers is the better choice. Similarly, Wikipedia *defines* a [plane partition](https://en.wikipedia.org/wiki/Plane_partition) as a 2D array of natural numbers, subject to constraints that correspond to the Ferrers diagram not falling over. This is well suited to a computer: it is a compact representation; the compatibilities are inequalities between adjacent entries; adding a box corresponds to incrementing a number. 


CFTP requires that you have a Gibbs sampling Markov chain. The way Gibbs sampling works is that you erase part of the diagram, compute all possible ways of filling it in and uniformly choose from this small set. You could probably create a Gibbs sampling Markov chain from any of your representations whose stationary distribution would be uniform.


But CFTP uses more structure, specifically a lattice structure, a partial order such that sup and inf exist. Plane partitions have this structure by product partial order — one partition is bigger than another if all of its components are bigger. If you want to use a different representation, you must check that it has this structure.


Actually, there is a slightly different representation that might be better and has been popular in the past. In the case of regular partitions, it turns the partition 45 degrees, so that the gravity pulling boxes into the corner is now the $y$ axis. The $x$-axis parameterizes diagonal stacks of boxes, touching corner to corner. We record an array index by $x$ of the heights of these stacks. This restores the symmetry of the axes. Similarly, we can represent a (boxed) plane partition as a hexagonal grid of heights, where height is $x+y+z$.


