Running a Chain With ReCom

In the previous section, we used a proposal function, propose_random_flip, which “flipped” a single node from one district to another to create the next district plan in the Markov chain. This approach evolves slowly from its initial state.

An approach which evolves much faster is to combine two adjacent districts and to then randomly split the combined set of nodes into two new districts.

The GerryChain ReCom class implements proposal functions that do precisely that. You can see how it works in the animation below:

The goal of this section is to become familiar with using the ReCom class.

One of the things we will investigate is how to influence the way ReCom decides to split the combined set of nodes into two new districts. Some states have rules that prioritize keeping some nodes together - for instance, to the extent possible, to not split a county into two districts. GerryChain has implemented so-called “region-aware” settings to accomplish this.

Throughout this guide, we’ll use the toy state of GerryMandria, which needs to be divided into 8 districts.

The legislature of the state of GerryMandria has provided us with the following districting plan:

The corresponding dual graph (showing nodes and edges between nodes) is:

gerrychain works internally with the dual graph of the districting plan, that is, the algorithms in GerryChain are all about nodes and edges.

To make it easier to visualize how ReCom works, the graphics of district plans use the equivalent colored grid layout so that the changing districts and the regions they intersect are easier to compare.

Setup

Let us start by running a simple ReCom chain on this districting plan. Of course, the first thing to do is to import the required packages:

from gerrychain import Partition, Graph, MarkovChain, updaters, constraints, accept
from gerrychain.proposals import ReCom
from gerrychain.constraints import contiguous
from gerrychain.examples import gerrymandria
from functools import partial
import pandas

A Simple ReCom Chain

We begin by instantiating the MarkovChain:

recom_chain = MarkovChain(
    total_steps=40,
    rng=2024,
)

We now need to set the starting state for the chain.

# Use a custom function for this guide to get the initial graph.
graph = gerrymandria()

my_updaters = {
    "population": updaters.Tally("TOTPOP", alias="population"),
}

recom_chain.initial_partition = Partition(graph, assignment="district", updaters=my_updaters)

Now we attach the proposal function

# This should be 8 since each node has population 1 and each district has 8 nodes.
# Note that the key "population" corresponds to the population updater
# that we defined above and not with the population column in the JSON file.
ideal_population = sum(recom_chain.initial_partition["population"].values()) / len(
    recom_chain.initial_partition
)

recom_chain.proposal_fn = ReCom.district_pairs_mst(
    pop_col="TOTPOP",
    pop_target=ideal_population,
    epsilon=0.01,
)

With the proposal function attached, we may now iterate through the chain. For now, we will collect the assignment for each step, so that we can watch the chain work in a fun animation (of course, it would be a bad idea to do this for a chain with a large number of steps).

assignment_list = []

for i, item in enumerate(recom_chain):
    assignment_list.append(item.assignment)

To create the animation, we need to define some functions to do the animation. If the way these functions work is not clear, don’t worry about it - they are just a way to show you how ReCom works.

import numpy as np
import matplotlib.pyplot as plt
import math
import io
from matplotlib.colors import ListedColormap
from PIL import Image
from IPython.display import display, clear_output

DISTRICTR_COLORS = [
    "#0099cd",
    "#ffca5d",
    "#00cd99",
    "#99cd00",
    "#cd0099",
    "#9900cd",
    "#8dd3c7",
    "#bebada",
]


# Define a function to "plot" a GerryChain assignment but as an image
# to be displayed later.
def plot_assignment_to_image(this_assignment, graph):
    # Create a grid (numpy array) with district ids (integers)
    # to use within matplotlib to plot the assignment

    # Assumes a square grid (graph)
    grid_size = math.isqrt(len(this_assignment))

    grid = np.empty((grid_size, grid_size))

    district_labels = set()
    # Get (x,y) position from each node in the graph
    # and set grid value for that position to the district_id
    for node_id, district_id in this_assignment.items():
        x_pos = graph.node_data(node_id)["x"]
        y_pos = graph.node_data(node_id)["y"]
        grid_pos = (y_pos, x_pos)
        grid[grid_pos] = district_id
        district_labels.add(district_id)

    fig, ax = plt.subplots(figsize=(grid_size + 1, grid_size + 1))
    im = ax.imshow(grid, cmap=ListedColormap(DISTRICTR_COLORS[: len(district_labels)]))

    ax.set_xticks(np.arange(-0.5, grid_size, 1), minor=True)
    ax.set_yticks(np.arange(-0.5, grid_size, 1), minor=True)

    ax.grid(which="minor", color="black", linestyle="-", linewidth=1)

    # Hide the numeric labels (axis values)
    ax.set_xticks([])
    ax.set_yticks([])

    buffer = io.BytesIO()
    plt.savefig(buffer, format="png", bbox_inches="tight", pad_inches=0)
    buffer.seek(0)
    image = Image.open(buffer)
    plt.close(fig)
    return image

And now let’s check to make sure that it works by plotting an arbitrary assignment - say, assignment_list[35]

saved_district_plan = plot_assignment_to_image(assignment_list[35], graph)
display(saved_district_plan)

Next, we define another utility routine for exploring a sequence of assignments by moving forward and backward through it. This lets you see what GerryChain did at each step.

Note that the plot does not change at every step - this is because sometimes ReCom splits a merged pair of districts in exactly the same way that the districts were before being merged, which looks like nothing happened.

from matplotlib.animation import ArtistAnimation
from IPython.display import HTML


def create_assignment_images(assignment_list, graph, plotting_function=plot_assignment_to_image):
    """Render each district assignment once."""
    return [plotting_function(assignment, graph) for assignment in assignment_list]


def create_assignment_viewer(image_list):
    """Display controls for moving through a sequence of district plans."""
    if not image_list:
        raise ValueError("image_list must contain at least one image")

    fig, ax = plt.subplots(figsize=(5, 5))
    ax.axis("off")
    fig.subplots_adjust(left=0, right=1, bottom=0, top=1)
    frames = [[ax.imshow(image, animated=True)] for image in image_list]
    animation = ArtistAnimation(fig, frames, interval=500)
    plt.close(fig)
    return HTML(animation.to_jshtml())

assignment_images = create_assignment_images(assignment_list, graph)
create_assignment_viewer(assignment_images)

Once Loop Reflect

Use the controls above to move through the chain. If you are running this notebook yourself, the matplotlib player also offers a slider, playback, and speed controls; the rendered documentation hides everything except the slider, Back, and Forward buttons to keep the page simple.

The plan may stay unchanged at some steps. ReCom can reconstruct the same split after merging two districts, and constraints or the acceptance function can also produce a self-loop. These repeated states are part of the Markov chain and can carry statistical importance, so they should generally not be removed from the ensemble.

Region-Aware ReCom

Of course, in the state of GerryMandria, the legislature has decided that it would like to try to keep the municipality of Gerryville together in a single district. In fact, it would really prefer to keep all of the municipalities together if possible. As such, any analysis that you do needs to use an ensemble designed to keep municipalities together. Here is a picture of the municipalities in GerryMandria:

Fortunately, gerrychain has built-in support for region-aware ReCom chains, which create ensembles of districting plans that try to keep particular regions of interest together. And it only takes one extra line of code: we simply update our proposal to include a region_surcharge which increases the importance of the edges within the municipalities.

The way this works under the covers is described later in this guide. For now, all that is important is that ReCom randomizes the creation of each successive district plan by assigning a random weight (between 0.0 and 1.0) to edges, and that adding a surcharge weight to an edge in addition to the random weight will make it more likely that the edge will be selected as a border edge between districts. That is, adding weight to edges that span regions will make it more likely that those edges will be chosen as border edges. The larger the surcharge weight, the more likely it is that edges that span regions will be chosen as border edges.

Let’s create an ensemble of district plans using a large region_surcharge to keep municipalities intact.

# Here we change the proposal function to contain a
# region_surcharge, but leave everything else about the chain the same.
recom_chain.proposal_fn = ReCom.district_pairs_mst(
    pop_col="TOTPOP",
    pop_target=ideal_population,
    epsilon=0.01,
    region_surcharge={"muni": 0.5},
)

assignment_list = []

for i, item in enumerate(recom_chain):
    assignment_list.append(item.assignment)

And then plot them

assignment_images = create_assignment_images(assignment_list, graph)
create_assignment_viewer(assignment_images)

Once Loop Reflect

The interactive viewer above shows the resulting municipality-aware ensemble where no municipality is split between two different districts.

The region_surcharge specifies how much additional weight to assign to edges that span regions. If we reduce the weight from 0.5 (which is a lot) to 0.3 (which is still a good bit), then the code will put less priority on preserving municipalities in districts, and we can see that below.

Surcharges are best read against the random weights they compete with, which are drawn from \([0,1)\). A surcharge of 0.1 only reorders edges whose random weights happen to land within 0.1 of each other, so it is a very weak bias. In this example, values from roughly 0.3 upward are where the effect becomes easier to observe.

# Rerun the region-aware code but with a reduced region_surcharge of 0.3

recom_chain.proposal_fn = ReCom.district_pairs_mst(
    pop_col="TOTPOP",
    pop_target=ideal_population,
    epsilon=0.01,
    region_surcharge={"muni": 0.3},
)

assignment_list = []

for i, item in enumerate(recom_chain):
    assignment_list.append(item.assignment)

assignment_images = create_assignment_images(assignment_list, graph)
create_assignment_viewer(assignment_images)

Once Loop Reflect

In this case, the district plans do a better job of keeping municipalities intact than plans generated without a surcharge, but they do not always keep municipalities intact.

Multiple Regions at Once

Now, the legislature of GerryMandria has decided that it would also like to try to keep the water districts intact in addition to keeping the municipalities intact.

Here is a picture of the water districts in GerryMandria:

The picture reflects the fact that a river cuts through the middle of the state, which means it is not going to be possible to keep all of the water districts together and all of the municipalities together in one plan. However, we can try to keep the water districts together as much as possible, and then, within those water districts, try to be sensitive to the boundaries of the municipalities. Again, this only requires us to edit the region_surcharge parameter of the proposal.

We reduce the surcharge for municipalities, “muni”, from 0.3 to 0.2, and add a surcharge for water districts, “water_dist”, of 0.8, making preserving water districts more important than preserving municipalities.

Since we are trying to be sensitive to multiple regions, we will increase the length of the example chain.

recom_chain.proposal_fn = ReCom.district_pairs_mst(
    pop_col="TOTPOP",
    pop_target=ideal_population,
    epsilon=0.01,
    # reduce "muni" surcharge, add "water_dist" surcharge
    region_surcharge={"muni": 0.2, "water_dist": 0.8},
)

# We update the total number of steps and the RNG seed
recom_chain.total_steps = 200
recom_chain.rng = 2024

assignment_list = []

for i, item in enumerate(recom_chain):
    if (i + 1) % 100 == 0:
        print(f"\rRunning step {i + 1}...", end="\r", flush=True)
    assignment_list.append(item.assignment)

Running step 100...

Running step 200...

And plot the last 40 assignments. We will also change the image display function so that it is easier to see the influence of the region surcharges on each of the partitions.

You also see that the chain tends to repeat assignments more frequently than previously. This is to be expected since surcharging the edges of the spanning tree fundamentally alters the kind of trees that are likely to be drawn.

def plot_plan_with_regions(this_assignment, graph):
    """Plot a district plan with both kinds of region outlined on top of it.

    Municipalities are boxed in black and water districts outlined in white. A municipality whose
    box holds a single color was kept intact; one holding two colors was split between two
    districts.
    """
    grid_size = math.isqrt(len(this_assignment))
    districts = np.empty((grid_size, grid_size))
    water = np.empty((grid_size, grid_size), dtype=object)
    muni = np.empty((grid_size, grid_size), dtype=object)
    for node_id, district_id in this_assignment.items():
        data = graph.node_data(node_id)
        districts[data["y"], data["x"]] = district_id
        water[data["y"], data["x"]] = data["water_dist"]
        muni[data["y"], data["x"]] = data["muni"]

    fig, ax = plt.subplots(figsize=(grid_size + 1, grid_size + 1))
    ax.imshow(districts, cmap=ListedColormap(DISTRICTR_COLORS), zorder=0)
    ax.set_xticks([])
    ax.set_yticks([])

    for muni_id in sorted(set(muni.flat), key=int):
        nodes = [
            (x_pos, y_pos)
            for y_pos in range(grid_size)
            for x_pos in range(grid_size)
            if muni[y_pos, x_pos] == muni_id
        ]
        left = min(x_pos for x_pos, _ in nodes) - 0.5
        right = max(x_pos for x_pos, _ in nodes) + 0.5
        top = min(y_pos for _, y_pos in nodes) - 0.5
        bottom = max(y_pos for _, y_pos in nodes) + 0.5
        ax.plot(
            [left, right, right, left, left],
            [top, top, bottom, bottom, top],
            color="black",
            linewidth=0.8,
            zorder=3,
        )

    # A white line wherever two neighboring nodes fall in different water districts.
    for y_pos in range(grid_size):
        for x_pos in range(grid_size):
            if x_pos + 1 < grid_size and water[y_pos, x_pos] != water[y_pos, x_pos + 1]:
                ax.plot(
                    [x_pos + 0.5, x_pos + 0.5],
                    [y_pos - 0.5, y_pos + 0.5],
                    color="white",
                    linewidth=5,
                    zorder=2,
                )
            if y_pos + 1 < grid_size and water[y_pos, x_pos] != water[y_pos + 1, x_pos]:
                ax.plot(
                    [x_pos - 0.5, x_pos + 0.5],
                    [y_pos + 0.5, y_pos + 0.5],
                    color="white",
                    linewidth=5,
                    zorder=2,
                )
    buffer = io.BytesIO()
    plt.savefig(buffer, format="png", bbox_inches="tight", pad_inches=0)
    buffer.seek(0)
    image = Image.open(buffer)
    plt.close(fig)
    return image

assignments_to_plot = assignment_list[-40:]
assignment_images = create_assignment_images(
    assignments_to_plot, graph, plotting_function=plot_plan_with_regions
)
create_assignment_viewer(assignment_images)

Once Loop Reflect

The combined picture above shows that the chain has done a pretty good job of keeping the water districts together (the regions bounded by thick, white lines) while also being sensitive to the municipalities (the regions bounded by black boxes). Most municipalities that lie inside a water district are kept whole, and where the two goals conflict, the municipalities tend to be the ones that give way, which is what the larger "water_dist" surcharge asked for.

How the Region-Aware Implementation Works

ReCom works by combining all of the nodes in two adjacent districts and then randomly splitting that set of nodes back into two districts. It does this by creating what is known as a spanning tree.

In graph theory, a “tree” is a special kind of graph. Like a graph, it is composed of nodes and edges, but a tree contains no cycles; that is, you cannot follow edges and end up back where you started. A spanning tree for a given graph is a tree that contains all of the nodes in the graph (but typically not all of the edges). The characteristic of a spanning tree that is important for GerryChain is that you can create two connected subtrees out of a spanning tree by just cutting one edge. All of the nodes on one side of the cut edge will be connected and all of the nodes on the other side of the cut edge will be connected. This is exactly what we want ReCom to do - split a collection of nodes into two connected subsets.

The ReCom variant that we have been using in this guide, district_pairs_mst, works by randomly assigning weights from the interval \([0,1]\) to the edges in the combined set of nodes from the two adjacent districts. It then builds a spanning tree by greedily selecting the lowest weighted edge next. This has the effect of randomizing which spanning tree is created from the universe of all possible spanning trees.

When working with region-aware ReCom chains, an additional surcharge weight is applied to the edges in the graph that span different regions specified by the region_surcharge dictionary. This has the effect of making surcharged edges less likely to be selected next when building a spanning tree. As a result, nodes without a surcharge are more likely to be selected first and therefore close to each other in the spanning tree. When the spanning tree is cut, these un-surcharged nodes will tend to be in the same district.

So in our example above, we have region_surcharge={"muni": 0.2, "water_dist": 0.8}, which means that the edges that span different municipalities will be upweighted by 0.2 and the edges that span different water districts will be upweighted by 0.8. We then create a minimum spanning tree by greedily selecting the lowest-weight edges via Kruskal’s algorithm. The surcharges on the edges help ensure that the algorithm picks the edges interior to the region before it picks the edges that bridge different regions.

This makes it more likely that each region is largely contained in a connected subtree attached to a bridge node. Thus, when we make a cut, the regions attached to the bridge node are more likely to be (mostly) preserved in the subtree on either side of the cut.

In the implementation of gerrychain.tree.bipartition_tree we further bias this choice by deterministically cutting bridge edges first (when possible). In the event that multiple types of regions are specified, the surcharges are added together, and edges are selected first by the number of types of regions that they span, and then by the surcharge added to those weights. So, if we have a region surcharge dictionary of {"a": 1, "b": 4, "c": 2} then edges which bridge all three regions, “a”, “b”, and “c”, would be cut first (total weight of 7), and then edges which bridge regions, “b” and “c” (total weight of 6), etc. So, the order in which edges would be selected would be:

• (“a”, “b”, “c”)

• (“b”, “c”)

• (“a”, “b”)

• (“a”, “c”)

• (“b”)

• (“c”)

• (“a”)

• the maximum-weight fallback

Choosing the Cut Yourself

If the user wants different behaviour, such as choosing an edge at random, they can alter the cut_choice_fn used by bipartition_tree_fn to apply a criterion other than the default described above.

Because the build_recom_proposal_fn expects a function reference for the value of its bipartition_tree_fn, we need to create a new function using Python’s functools.partial() that will bind the cut_choice_fn to be used by the bipartition_tree_fn. GerryChain supplies the chain-owned RNG automatically:

from gerrychain.tree import bipartition_tree
from functools import partial


def choose_random_cut(cuts, *, rng):
    """Choose one balanced cut uniformly at random."""
    return rng.choice(cuts)


my_new_bipartition_tree_fn = partial(
    bipartition_tree,
    cut_choice_fn=choose_random_cut,
)

To those of you for whom functools.partial() is a new concept:

Using functools.partial

The functools.partial function allows us to create a new function from an existing function by binding the values of some of the arguments. For example, we might have a function to make a colored square:

from PIL import Image

def make_color_square(red_val, green_val, blue_val):
    img = Image.new('RGB', (100, 100), color = (red_val, green_val, blue_val))
    return img

And we can then use this to make a new function that always makes a blue square:

make_blue_square = partial(make_color_square, red_val=0, green_val=0)

make_color_square(red_val=255, green_val=0, blue_val=0).show() # Makes a red square
make_blue_square(blue_val=255).show() # Makes a blue square

To build the updated proposal function, we will need to make use of the more verbose build_recom_proposal_fn function. The functions attached to the ReCom namespace that we used before are intentionally slim to make them easy to use, but there are times when the limited set of parameters that they expose are insufficient.

from gerrychain.tree import bipartition_tree
from gerrychain.proposals import build_recom_proposal_fn


recom_chain.proposal_fn = build_recom_proposal_fn(
    pop_col="TOTPOP",
    pop_target=ideal_population,
    epsilon=0.01,
    region_surcharge={
        "muni": 0.2,
        "water_dist": 0.8,
    },
    bipartition_tree_fn=my_new_bipartition_tree_fn,
)

recom_chain.rng = 2027

assignment_list = []

for i, item in enumerate(recom_chain):
    if i % 100 == 0:
        print(f"\rRunning step {i}...", end="", flush=True)
    assignment_list.append(item.assignment)

Running step 0...
Running step 100...

And plot the last 40 assignments.

assignments_to_plot = assignment_list[-40:]
assignment_images = create_assignment_images(assignments_to_plot, graph)
create_assignment_viewer(assignment_images)

Once Loop Reflect

Surcharges should be interpreted relative to the base random edge weights in \([0,1)\). Since the standard random() method of Python’s built-in random module outputs values in \([0,1)\), surcharges closer to 1 will have a stronger effect on region preservation compared to surcharges closer to 0. Surcharges above 1 are allowed (as shown below), but they do not have any stronger effect than a surcharge of 1.

What to Do if the Chain Gets Stuck

Sometimes, either because of the constraints that you have imposed or because of the shape of the graph that you are working with, a ReCom chain can get stuck and will throw an error. For example, if we try to be a bit too demanding of the region-aware chain given above and ask for a plan that effectively never splits a municipality nor a water district, then the chain will get stuck and throw an error. Here is the setup:

from gerrychain import Partition, Graph, MarkovChain, updaters, constraints, accept
from gerrychain.proposals import build_recom_proposal_fn
from gerrychain.tree import bipartition_tree
from gerrychain.constraints import contiguous
from gerrychain.examples import gerrymandria
from functools import partial


recom_chain = MarkovChain(
    total_steps=20,
    rng=0,
)

my_updaters = {
    "population": updaters.Tally("TOTPOP"),
    "cut_edges": updaters.cut_edges
}

graph = gerrymandria()
recom_chain.initial_partition = Partition(
    graph,
    assignment="district",
    updaters=my_updaters
)

ideal_population = sum(initial_partition["population"].values()) / len(initial_partition)

recom_chain.proposal_fn = build_recom_proposal_fn(
    pop_col="TOTPOP",
    pop_target=ideal_population,
    epsilon=0.01,
    region_surcharge={
        "muni": 2.0,
        "water_dist": 2.0,
    },
    bipartition_tree_fn=partial(
        bipartition_tree,
        max_attempts=100,
        warn_attempts=50,
    ),
)

assignment_list = []
for item in recom_chain:
    assignment_list.append(item.assignment)

This deterministic example emits a warning and then raises an error:

BipartitionWarning:
Failed to find a balanced cut after 50 attempts.

RuntimeError: Could not find a possible cut after 100 attempts.

Here, max_attempts is the total number of balanced-cut searches for the selected pair of districts. node_repeats is the number of additional roots tried on each spanning tree, so a tree is searched node_repeats + 1 times. With the default memoized cut finder, each search already examines every edge in the tree. Re-rooting cannot make an uncuttable tree cuttable, so the default node_repeats=0 redraws immediately, and setting it above zero with the memoized finder emits a warning. Positive values remain useful with the contraction finder or a custom cut-edge finder whose result can depend on the root choice.

Pair Reselection

The hard regional surcharges can still make a particular district pair difficult to split. Pair reselection lets ReCom try another adjacent pair after exhausting the budget. It helps when some other adjacent pair is easier to split; if epsilon is tight across the whole plan, every pair is hard and reselection mostly delays the error:

recom_chain.proposal_fn = build_recom_proposal_fn(
    pop_col="TOTPOP",
    pop_target=ideal_population,
    epsilon=0.01,
    region_surcharge={
        "muni": 2.0,
        "water_dist": 2.0,
    },
    bipartition_tree_fn=partial(
        bipartition_tree,
        max_attempts=100,
        warn_attempts=50,
        allow_pair_reselection=True,
    ),
)
recom_chain.total_steps = 20
recom_chain.rng = 0


assignment_list = [item.assignment for item in recom_chain]

This chain completes all 20 steps because it can move on from an unsuitable district pair.

When Epsilon Is the Real Problem

Surcharges are not the only way to get stuck, and, more often than not, a chain will stall if the bound for \(\varepsilon\) is too tight.

It is worth being precise about what “stuck” means here. A merged pair always has at least one valid cut. ReCom only ever merges adjacent districts, so a spanning tree of the first district joined to a spanning tree of the second by a single edge between them is itself a spanning tree of the merged pair, and cutting that edge restores the two districts you started from. Those were population-balanced already, so that cut is valid by construction.

The difficulty is that such trees can be rare. Cutting a spanning tree splits it at whole nodes, so the reachable populations are fixed by the node populations, and a tight epsilon may admit only a few of them. Each attempt draws a fresh tree, so if only a small fraction of trees happen to contain an acceptable edge, all max_attempts draws can miss.

A small example makes this concrete. Suppose a merged pair holds 16 nodes of population 1000 each, so pop_target is 8000. Every candidate cut yields a multiple of 1000, so the reachable populations are 7000, 8000, 9000, and so on:

epsilon	acceptable window	cuts that fit
0.01	7920 – 8080	8000
0.1	7200 – 8800	8000
0.15	6800 – 9200	7000, 8000, 9000

At epsilon=0.01 only an exactly even split qualifies, so a tree is usable only when one of its edges happens to separate 8 nodes from the other 8. Widening to epsilon=0.1 changes nothing, because 7000 and 9000 are still outside the window. At epsilon=0.15 any 7/8/9 split works, so far more trees qualify and the search succeeds quickly.

So when a chain stalls, compare your tolerance against the granularity of your data:

• Compare epsilon * pop_target with typical node populations. If the window is narrower than a single node, only near-perfect splits qualify and cuttable trees are scarce. Real precinct data usually has enough variety to hide this, but small or highly uniform graphs do not.

• Relax epsilon as a diagnostic. If a chain that fails at epsilon=0.01 runs at epsilon=0.05, the tolerance was the binding constraint, not the graph or the surcharges.

• More attempts help, but with diminishing returns. Raising max_attempts improves the odds of drawing a cuttable tree, since one always exists. When acceptable cuts are very rare, though, you pay a lot of time for that chance, and widening epsilon is usually the more effective lever.

• Watch for interactions. Tight surcharges and a tight epsilon compound: the surcharges restrict where a tree is likely to be cut, and epsilon restricts which of those cuts count.

Seed generation is a different matter. recursive_tree_part splits a region that has no pre-existing valid division, so a tight epsilon can genuinely be infeasible rather than merely unlikely.

A Real-World Example

In this example, we’ll use GerryChain to analyze Pennsylvania’s 2011 congressional districting plan. We’ll compare the partisan vote shares in the 2011 plan to those in an ensemble of districting plans generated by our ReCom chain.

Imports

As always, the first step is to import everything we need

import matplotlib.pyplot as plt
from gerrychain import (
    GeographicPartition,
    Partition,
    Graph,
    MarkovChain,
    updaters,
    constraints,
    accept,
    Election,
)
import pandas

Setting Up the Markov Chain

Download PA File

pa_chain = MarkovChain(
    total_steps=1000,
    rng=2024,
)

We’ll create our graph using the example Pennsylvania JSON file.

graph = Graph.from_json("./PA_VTDs.json")

Then we instantiate the initial partition with the standard population updater.

# Population updater, for computing how close to equality the district
# populations are. "TOTPOP" is the population column from our shapefile.
my_updaters = {"population": updaters.Tally("TOT_POP", alias="population")}

pa_chain.initial_partition = GeographicPartition(
    graph, assignment="2011_PLA_1", updaters=my_updaters
)

The GeographicPartition class comes with built-in area and perimeter updaters. We do not use them here since (i) the JSON file that we are working with does not have geometric information and (ii) geometric updaters tend to slow the chain quite considerably (and this is just an example), but they would allow us to compute compactness scores like Polsby-Popper that depend on these measurements.

For this example, we are interested in tracking election outcomes in a few elections. To this end, we will add some Election updaters to the chain.

Remember that the data used by our updaters exists as node_data in our graph, and the node_data in our graph came from the file that was used to create our graph.

elections_updaters = {
    "SEN10": Election("SEN10", {"Democratic": "SEN10D", "Republican": "SEN10R"}),
    "SEN12": Election("SEN12", {"Democratic": "USS12D", "Republican": "USS12R"}),
    "SEN16": Election("SEN16", {"Democratic": "T16SEND", "Republican": "T16SENR"}),
    "PRES12": Election("PRES12", {"Democratic": "PRES12D", "Republican": "PRES12R"}),
    "PRES16": Election("PRES16", {"Democratic": "T16PRESD", "Republican": "T16PRESR"}),
}

pa_chain.add_updaters(elections_updaters)

Now we add the proposal function

ideal_population = sum(pa_chain.initial_partition["population"].values()) / len(
    pa_chain.initial_partition
)

pa_chain.proposal_fn = build_recom_proposal_fn(
    pop_col="TOT_POP",
    pop_target=ideal_population,
    epsilon=0.02,
)

To keep districts about as compact as the original plan, we would like to constrain the number of cut edges between all of the districts. We can do this using the gerrychain.constraints.UpperBound constraint, and, as a general heuristic, we’ll bound the number of cut edges by twice the number of cut edges in the initial plan.

def cut_edges_length(partition):
    return len(partition["cut_edges"])


pa_chain.add_constraint(
    constraints.UpperBound(cut_edges_length, 2 * len(pa_chain.initial_partition["cut_edges"]))
)

Coding note

We can simplify the calling of this compactness bound using lambda functions.

compactness_bound = constraints.UpperBound(
  lambda p: len(p["cut_edges"]),
  2*len(initial_partition["cut_edges"])
)

The use of lambda functions tends to be a more advanced coding technique, but the benefit is that we do not need to define a new function for each constraint that we want to use, and they can make the code more readable.

Running the Chain

Now we’ll run the chain, putting the sorted Democratic vote percentages directly into a pandas DataFrame for analysis and plotting. The DataFrame will have a row for each state of the chain. The first column of the DataFrame will hold the lowest Democratic vote share among the districts in each partition in the chain, the second column will hold the second-lowest Democratic vote share, and so on.

# This might take a few minutes.

data = pandas.DataFrame(sorted(partition["SEN12"].percents("Democratic")) for partition in pa_chain)

If you are wondering what the for loop inside of the parentheses is doing, please see this note. If you install the tqdm package, you can see a progress bar as the chain runs by running this code instead

data = pandas.DataFrame(
    sorted(partition["SEN12"].percents("Democratic"))
    for partition in chain.with_progress_bar()
)

Create a Plot

Now we’ll create a box plot to help visualize the data report.

fig, ax = plt.subplots(figsize=(8, 6))

# Draw 50% line
ax.axhline(0.5, color="#cccccc")

# Draw boxplot
data.boxplot(ax=ax, positions=range(len(data.columns)), grid=False, whis=(1, 99))

# Draw initial plan's Democratic vote %s (.iloc[0] gives the first row)
plt.plot(data.iloc[0], "ro")

# Annotate
ax.set_title("Comparing the 2011 plan to an ensemble")
ax.set_ylabel("Democratic vote % (Senate 2012)")
ax.set_xlabel("Sorted districts")
ax.set_ylim(0, 1)
ax.set_yticks([0, 0.25, 0.5, 0.75, 1])

plt.show()

The exact plot is reproducible when the input data, GerryChain version, configuration, and seed are all the same. Different seeds produce different valid random samples. This short chain is an illustration of the workflow, not evidence by itself for a legal or statistical conclusion.

Multi-Member ReCom

Every chain so far has drawn single-member districts, where each district elects one representative. GerryChain also supports multi-member districts, where a district elects a fixed number of representatives, and where different districts may elect different numbers.

Why Would Districts Elect Different Numbers of Members?

The main use for this is election modeling, for bodies like city councils and school boards that fill several seats from a single district.

Uneven member counts arise more naturally than they might first appear. Such a body is often required to have an odd number of members so that votes cannot tie, and when that total does not divide evenly among the districts, uneven districts are forced. A seven-member council drawn into three districts has to split the seats 3-2-2, because there is no even option.

How Member Counts Work in the Code

A key distinction between the standard ReCom and Multi-Member ReCom is that in the Multi-Member variant, member counts attach to district labels. So, if you declare districts 1 and 2 to both have 1 member, district 3 to have 2 members, and district 4 to have 4 members, this will remain invariant for all steps of the chain.

This is why population is balanced per member rather than per district. Each district is balanced against pop_target * members_per_district[label], where pop_target is the target for a single member, so a four-member district should hold roughly four times the population of a one-member district rather than the same amount.

That also means the equal-population constraint used earlier in this guide is not appropriate here: it would try to force every district to the same size regardless of its seat count.

Setting Up the Initial Partition

There are two things worth keeping separate:

1. building a starting plan that already satisfies the member counts, and

2. interpreting what comes out of the chain.

The starting plan must be contiguous and already consistent with the member counts, and the keys of members_per_district have to match the partition’s district labels exactly. Here we combine GerryMandria’s eight columns into four districts electing one, one, two, and four members.

from gerrychain.proposals import MultiMemberReCom

multi_member_chain = MarkovChain(
    total_steps=40,
    rng=42,
)

multi_member_graph = gerrymandria()

# Each district label carries a fixed member count for the whole run.
members_per_district = {"1": 1, "2": 1, "3": 2, "4": 4}

# Columns 0-7 of the grid are grouped so that each district starts out with
# 8 people per member: 1 column for districts "1" and "2", 2 for "3", 4 for "4".
column_to_district_label = {
    0: "1",
    1: "2",
    2: "3",
    3: "3",
    4: "4",
    5: "4",
    6: "4",
    7: "4",
}
multi_member_assignment = {
    node: column_to_district_label[multi_member_graph.node_data(node)["x"]]
    for node in multi_member_graph.node_indices
}

multi_member_chain.initial_partition = Partition(
    multi_member_graph,
    assignment=multi_member_assignment,
    updaters={
        "population": updaters.Tally("TOTPOP", alias="population"),
        "cut_edges": updaters.cut_edges,
    },
)

initial_multi_member_partition = multi_member_chain.initial_partition
print("Raw population by district:")
print(f"\t{dict(initial_multi_member_partition['population'])}")
print("Population per member:")

pop_per_member = {
    part: initial_multi_member_partition["population"][part] / members_per_district[part]
    for part in sorted(initial_multi_member_partition.parts)
}
print(f"\t{pop_per_member}")

Raw population by district:
	{'1': 8, '2': 8, '3': 16, '4': 32}
Population per member:
	{'1': 8.0, '2': 8.0, '3': 8.0, '4': 8.0}

The two printouts show the distinction that matters. The raw populations are 8, 8, 16, and 32, which look wildly unbalanced. Divided by each district’s member count, they are all 8, which is exactly the balance a multi-member plan is supposed to have.

Attaching the Proposal and the Constraint

Like in the single member case, pop_target is the target population per member of the legislative body. So, to properly instantiate the Multi-Member chain, we also need to indicate the number of members each district is intended to have by passing the same members_per_district dictionary constructed above.

pop_target = sum(initial_multi_member_partition["population"].values()) / sum(
    members_per_district.values()
)
print(f"Target population per member: {pop_target}")

multi_member_chain.proposal_fn = MultiMemberReCom.district_pairs_mst(
    pop_col="TOTPOP",
    pop_target=pop_target,
    epsilon=0.01,
    members_per_district=members_per_district,
)

Target population per member: 8.0

Interpreting the Output

Now we can run the chain.

populations_per_member = []
assignment_list = []
for partition in multi_member_chain:
    assignment_list.append(partition.assignment)
    populations_per_member.append(
        {
            part: population / members_per_district[part]
            for part, population in partition["population"].items()
        }
    )

final_partition = partition
print("Final raw population by district:")
print(
    f"\t{ {part: final_partition['population'][part] for part in sorted(final_partition.parts)} }"
)
print("Final population per member:")
print(f"\t{ {part: populations_per_member[-1][part] for part in sorted(final_partition.parts)} }")

Final raw population by district:
	{'1': 8, '2': 8, '3': 16, '4': 32}
Final population per member:
	{'1': 8.0, '2': 8.0, '3': 8.0, '4': 8.0}

The raw populations still come out 8, 8, 16, and 32, while the per-member figures are all 8, even though ReCom has been redrawing the boundaries the whole time. District "4" is a different set of nodes than it was at the start, but it is still the four-member district and it still holds four members’ worth of people. We can see this in action in the animation below:

assignment_images = create_assignment_images(assignment_list, multi_member_graph)
create_assignment_viewer(assignment_images)

Once Loop Reflect

Recap

A ReCom proposal can be used directly, made region-aware with surcharges, or customized with a cut-selection function. If a difficult district pair exhausts the cut-search budget, pair reselection lets the chain try another adjacent pair. For bodies that fill several seats from one district, MultiMemberReCom attaches a fixed member count to each district label and balances population per member.

To build on this example, here are some possible next steps:

• Add, remove, or tweak the constraints

• Perform a similar analysis on a different districting plan for Pennsylvania

• Perform a similar analysis on a different state

• Compute partisan symmetry scores like Efficiency Gap or Mean-Median, and create a histogram of the scores of the ensemble.

• Perform the same analysis using a different election than the 2012 Senate election

• Collect Democratic vote percentages for all the elections we set up, instead of just the 2012 Senate election.