Wilson's Uniform Spanning Tree Algorithm Animation

Make gif animations of Wilson's uniform spanning tree algorithm and the depth/breadth - first search algorithm.

How to use this program

Run wilspn.py and wait for a few seconds, you will see a .gif file generated in current directory. Enjoy it!

This program can be run with both python2.7+ and python3+. It's written in pure python: no third-party modules nor softwares are needed, just built-in modules struct and random and some built-in functions. I could write it faster by using numpy arrays but I prefer to keep the code being "pure blooded".

The code has been modified many times to make it more readable and improve the efficiency. Since Wilson's algorithm is a random algorithm, its runtime is uncertain. On my old laptop with 1.83GHzx4 Celeron processors the bitrate is about 1000kb/min, or in other words, it takes 30 seconds to generate a file of 500kb.

Note: I have added another script ust.py to show how Wilson's algorithm works. It outputs the final spanning trees as static .png images, not GIF animations. The code is quite short and straight-forward.

What is Wilson's uniform spanning tree algorithm

Consider the following problem:

Let G be a finite, connected graph. How can we choose a random spanning tree among all spanning trees of G with equal probability? (we shall call such a tree an uniform spanning tree, or simply an UST.)

You might say "that's easy, just write a program that lists all spanning trees and then use a random integer to choose one". But for complete graphs K_n, it has n**(n-2) different spanning trees by Cayley's formula, for n=100 this number is 100**98, much more than the number of particles in the universe! (roughly 10**90)

Currently the most efficient algorithm is the one proposed in Wilson's paper

Generating random spanning trees more quickly than the cover time.

It's a random algorithm, which means some times you may be very lucky to get an UST quickly, or you may wait forever. But this algorithm will terminate in finite steps with probability one (note: this does not exclude the possibility of running forever, think about this), and it performs really well in most cases.

The key to understand Wilson's algorithm is the so called loop erased random walk, that is, once the random walk visits a vertex that already exists in its path, the walk immediately erase the loop formed by the two visits and continue the walk from this vertex. Just see the gif animation if you don't understand this, it's obvious to see what "loop erased random walk" means from the animation.

The algorithm runs as follows:

Wilson's algorithm:

1. Choose any vertex v as the root, maintain a tree T, initially T = {v}.

2. For any vertex z that is not in T, start a loop erased random walk from z, until the walk 'hits' T, then add the resulting path of the walk to T.

3. repeat step 2 until all vertices of the graph are in T.

For the proof of the correctness of this algorithm see Wilson's original paper or the book by Russell Lyons and Yuval Peres:

"Probability on Trees and Networks".

How did it come out

This project is motivated by Mike Bostock's wonderful Javascript animation, and also many other nice animations on his website. I learned Wilson's algorithm about 7 years ago and had the idea of writing a python version to produce GIF animations the first sight when I saw Mike's page, but rendering a GIF image which possibly consists of thousands of frames is definitely a formidable task. It's about one year ago when I occasionally touched the GIF89a specification and finally realized the approach of encoding the whole animation into a byte stream.

About the code

For most time you run the script you will get an image contains 2000-3000 frames while the file size is around 200-500KB. How could one do this? There are two key points:

1. Only update the region that has been changed by maintaining a rectangle that defines the position of current frame. If you use some image processing tool like ImageMagick to unpack the example gif image into frames you will soon understand what this means.

2. Write a GIF encoder that has minimal code length 2. The cells in our animation has 4 possible different states hence we need only 4 colors, and 4 colors can be represented by 2 bits. With this small palette and initial code length the file size can be significantly reduced.

Of course you have to understand the GIF89a specification first, espcially its transparent color feature and the LZW encoding algorithm.

The maze is represented by a 2D matrix like (we use a small 7x7 2d array for example)

m m m m m m m
m c w c w c m
m c w c w c m
m c w c w c m
m m m m m m m

Each character represents a square (5x5 pixels by default) in our image. m means "margin", these squares are served as the border of the image and are not used in our animation. The width of the margin can be changed (default to 2 but the example shows only 1). w means this a "wall", and c means this is a "cell". Our grid graph contains only the cells. Walls and margins are not considered to be part of the graph.

As the algorithm runs, the cells become connected and some walls might be set to "in tree" or "in path", but they are still not treated as part of the graph.

Where to learn the GIF89a specification

The best and possibly the only reference you will need is What's in a GIF. You should read these three articles word by word before you can encode your gif image bits by bits.

Things you should know when implementing your GIF encoder with python

In this section I'll give a not-so-detailed introduction to the GIF89a specification, you should always refer to What's in a GIF when you have difficulties understanding my words.

Structure of a GIF File

1. Always begins with 6 bytes GIF89a.
2. Then follows immediately with 7 bytes - the logical screen descriptor, which specify the width & height of the image and the size of the global color table.
3. Then follows immediately with the global color table.
4. Then comes with the actual data of the frames. The data of each frame is further divided into 3 parts:
   1. firstly the graphics control block, which specifies the delay and transparent color of this frame.
   2. then comes with the image descriptor, which specifies the relative position of this frame in the image window as well as the size of the local color table. We will not use any local color table in our program.
   3. Finally the LZW-encoded pixel data of the frame.
5. The image file ends with a single byte 0x3B.

Note: the description above does not apply to all GIF images, there can be some variations, for example:

1. There may be no global color table (so you have to specify a local color table for each frame).

2. You may insert a background image immediately after the global color table. You do not need to specify the graphcs control block for this background image. When combined with the transparent feature one can make effects like "painting on this background image". (we used this feature in the program)

   

   The above image shows what happens if a background image is missing: the region that has not been occupied by any frame will be set to transparent color.

3. The file may end without the byte 0x3B, most decoders can still decode the image correctly.

For beginners I suggest you do not care about these exceptional cases.

Now we explain the details of each part.

The head GIF89a

In python you may write it like

struct.pack('6s', b'GIF89a')

Why b'GIF89a' not simply 'GIF89a'? This is for compatibility with python3+ since python3's default encoding is unicode whereas python2's is ascii.

The logical screen descriptor

Example:

struct.pack('<2H3B', width, height, 0b10010001, 0, 0)

Here you shoud note the format string < (little endian).

The byte 0b10010001 is called a packed field. Let's read it from the left:

1. The first bit 1 means we have a global color table.
2. The next 3 bits specify the color depth. You need not care about what they mean, modern decoders like firefox and chrome do not use them.
3. The 5-th bit is the sort flag and is not used today, always set to 0.
4. The ending 3 bits represent an integer x in range 0-7, x specifies the size of the global color table is 2**(x+1).

The final two bytes are of little importance and won't be explained here.

Global color table

Example: if we want to use 4 colors red, green, blue, yellow, then it can be written as

bytearray([255, 0, 0, 0, 255, 0, 0, 0, 255, 255, 255, 0])

It's important that the number of colors in this array must match the size of the global color table specified in the packed byte in logical screen descriptor.

Graphics control block

Example:

struct.pack("<4BH2B", 0x21, 0xF9, 4, 0b00000101, delay, trans_index, 0)

The 3 beginning bytes 0x21, 0xF9, 4 are fixed and it's used to inform the decoder "Hey, I'm a graphics control block, see the next 4 bytes!"

The next byte 0b00000101 is also a packed field, read from left:

1. The first 3 bits are useless and are always 0.
2. The next 3 bits are called desposal method, they represent an integer x in range 0-7 and x specifies how we should dispose this frame after it's displayed.

   1. x=0 means it's undefined, decoders will use default 1 in this case.

   2. x=1 is the default, it means leave this frame here. So if the next frame is not overlapped with this frame, then these two frames will both be displayed. Otherwise the overlapped area in this frame will be covered by the next frame. See the example below:

      

   3. x=2 means remove this frame and clear the window to background image/color, see the following example:

      

      You can see each frame is remove immediately after it's displayed, so one can only see one frame at a time.

   4. x=3 also means remove this frame, but restore the image to the previous frame.

   5. 4-7 are unused.

The next 2 bytes delay specifies the delay of the frame in centisecond, so delay=3 means "last for 0.03 second".

The next byte trans_index specifies the transparent color in this frame, the pixels in this frame using this color are transparent: you can see through them to see the previous frame (if it is reserved) or the background image.

Image descriptor block

Example:

struct.pack('<B4HB', 0x2C, left, top, width, height, 0)

Quite straight-forward to understand. The last byte is 0 since we do not need any local color table.