Metropolized Forest RECOM
Overview
Metropolized Forest Recom Autry et al. (2021) is an extension of the Recom (or Recombination) algorithm developed by DeFord, Duchin, and Solomon (DeFord et al. 2021). It modifies the original, whose stationary measure is unknown, to an algorithm that is computationally feasible to use as a proposal for a Metropolis-Hastings scheme. This gives the user better knowledge of what measure the ensemble of maps has been drawn from and how they prioritize various legal and cultural norms in redistricting.
Additionally, the multiscale version has been used with great success to create redistrictings that preserve particular geographic structures. In our experience it works better than weights or score functions. This is especially true when one wants approach the minimum number of county splits possible.
There have been a number of different code bases that implemented this algorithm, including ones in C++, Python and Julia. We currently recommend the Julia codebase that is official Julia package and hence can be installed using the Julia Package manager.
More information can be found at the Git Repository.
Installation
Metropolized Forest Recom is written in the Julia programming language and registered as an official Julia package. It can be installed using the Julia package manager. Generally it is recommended that you run the code in an local environment. This can be done with the activate command. After activating the environment, we need to install the needed packages.
using Pkg
Pkg.activate(".")
Pkg.add("RandomNumbers")
Pkg.add("MetropolizedForestRecom")Once this installation is complete, one can activate the local environment and load the needed packages as follows.
using Pkg
Pkg.activate(".")
using RandomNumbers
using MetropolizedForestRecomA one-level example on NC
Download the Adjacency File
The first step is to download the needed JSON adjacency file. For this example it is call NC_pct21.json and was discussed more fully in this Section on JSON Adjacency Files. You can download the file to current directory using the following commandline at the terminal
curl https://raw.githubusercontent.com/jonmjonm/QGDocs/refs/heads/main/Geo/Adjacency/NC_pct21.json -o NC_pct21.jsonInitialize the Graph
# manually read in base_graph so that we can set a unique identifier that
# combines county and precinct id
pctGraphPath = joinpath(".", "NC_pct21.json")
nodeData = Set(["county", "prec_id", "pop2020cen", "area", "border_length"]);
base_graph = BaseGraph(pctGraphPath, "pop2020cen", inc_node_data=nodeData,
area_col="area", node_border_col="border_length",
edge_perimeter_col="length", edge_weights="connections");
# combines county and precinct id
for ii = 1:length(base_graph.node_attributes)
county = base_graph.node_attributes[ii]["county"]
prec_id = base_graph.node_attributes[ii]["prec_id"]
name = county*"_"*prec_id
base_graph.node_attributes[ii]["county_and_prec_id"] = name
end
# now that the field "county_and_prec_id" is set, we can use it to create the
# graph object that we will sample on
graph = Graph(base_graph, "county_and_prec_id");Add Constraints and Initial Random Partition
# add constraints
num_dists=14
pop_dev = 0.02
constraints = initialize_constraints()
add_constraint!(constraints, PopulationConstraint(graph, num_dists, pop_dev))
# create initial random partition
rng_seed = 110934571
rng = PCG.PCGStateOneseq(UInt64, rng_seed)
partition = Partition(graph, constraints, num_dists; rng=rng);Build Metropolis-Hastings Proposal Chain
proposal = build_forest_recom2(constraints)Specify Target Measure
gamma=0.0;
alpha=1.1;
iso=0.45;
measure = Measure(gamma, alpha); # spanning forest measure;
# first number is exponent on trees, second on linking edges
push_measure!(measure, get_isoperimetric_score, iso); # add iso parametric scoreSet Output
output_file_path = joinpath("output", "nc",
"metropolizedForestRECOM_1level_gamma"*string(gamma)*"_iso"*string(iso)*".jsonl.gz")
writer = Writer(measure, constraints, partition, output_file_path)
push_writer!(writer, get_log_spanning_trees)
push_writer!(writer, get_log_spanning_forests)
push_writer!(writer, get_isoperimetric_scores)Run Metropolized MCMC
steps = 100
output_freq=1
println("starting mcmc")
run_metropolis_hastings!(partition, proposal, measure, steps, rng,
writer=writer, output_freq=output_freq);Putting Everything Together
using Pkg
Pkg.activate(".")
using RandomNumbers
using MetropolizedForestRecom
# manually read in base_graph so that we can set a unique identifier that
# combines county and precinct id
pctGraphPath = joinpath(".", "NC_pct21.json")
nodeData = Set(["county", "prec_id", "pop2020cen", "area", "border_length"]);
base_graph = BaseGraph(pctGraphPath, "pop2020cen", inc_node_data=nodeData,
area_col="area", node_border_col="border_length",
edge_perimeter_col="length", edge_weights="connections");
# combines county and precinct id
for ii = 1:length(base_graph.node_attributes)
county = base_graph.node_attributes[ii]["county"]
prec_id = base_graph.node_attributes[ii]["prec_id"]
name = county*"_"*prec_id
base_graph.node_attributes[ii]["county_and_prec_id"] = name
end
# now that the field "county_and_prec_id" is set, we can use it to create the
# graph object that we will sample on
graph = Graph(base_graph, "county_and_prec_id");
# add constraints
num_dists=14
pop_dev = 0.02
constraints = initialize_constraints()
add_constraint!(constraints, PopulationConstraint(graph, num_dists, pop_dev))
# create initial random partition
rng_seed = 110934571
rng = PCG.PCGStateOneseq(UInt64, rng_seed)
partition = Partition(graph, constraints, num_dists; rng=rng);
# build proposal
proposal = build_forest_recom2(constraints)
# set output file and data
output_file_path = joinpath("output", "nc",
"metropolizedForestRECOM_1level_gamma"*string(gamma)*"_iso"*string(iso)*".jsonl.gz")
writer = Writer(measure, constraints, partition, output_file_path)
push_writer!(writer, get_log_spanning_trees)
push_writer!(writer, get_log_spanning_forests)
push_writer!(writer, get_isoperimetric_scores)
# run the MCMC Chain
steps = 100
output_freq=1
println("starting mcmc")
run_metropolis_hastings!(partition, proposal, measure, steps, rng,
writer=writer, output_freq=output_freq);A Two-level Example on NC
Now we will run MCMC on a phase space that consists of hierarchical graphs with the top level being counties and the finer level being precincts.
Loading Packages
Assuming we have already created a local environment as in the previous example. Then we can load the packages as follows.
import Pkg
Pkg.activate("./runMetropolizedRecomEnv")
Pkg.instantiate()
using RandomNumbers
using MetropolizedForestRecomInitialize the Graph
We now initialize a two-level graph. We specify the two layers by ["county", "prec_id"]. The graph object that is instantiation of a 2-level graph.
pctGraphPath = joinpath(".", "NC_pct21.json")
nodeData = Set(["county", "prec_id", "pop2020cen", "area", "border_length"]);
graph = Graph(pctGraphPath, "pop2020cen", ["county", "prec_id"];
inc_node_data=nodeData, area_col="area",
node_border_col="border_length", edge_perimeter_col="length",
edge_weights="connections")Constraints
We now create the set of constraints. We include the same constraint as before on the population deviation. But we also add a constraint that ensures pieces of districts in a county are contiguous and a constraint on the maximum number of split counties.
num_dists = 14
pop_dev = 0.02
constraints = initialize_constraints()
add_constraint!(constraints, PopulationConstraint(graph, num_dists, pop_dev))
# ensures pieces of districts in a county are contiguous; currently required
add_constraint!(constraints, ConstrainDiscontinuousTraversals(graph))
# fixes the maximum number of split counties
add_constraint!(constraints, MaxCoarseNodeSplits(num_dists+1)) Generate Initial Partition
Now we generate a random initial 2-level partition.
rng_seed = 454190
rng = PCG.PCGStateOneseq(UInt64, rng_seed)
partition = Partition(graph, constraints, num_dists; rng=rng);Proposal
Next, we instantiate the proposal chain which will be a multilevel Forest Recom. This proposal chain first mearges two districts. Then it loosely draws a uniform spanning tree on the top-level (which is counties in this case). Then it picks which edge to remove. This implies which county must be split. One then draws a uniform spanning tree on the finer level (which is precincts in this case). More details can be found in (Autry et al. 2021).
proposal = build_forest_recom2(constraints)Target Measure
iso=0.45
gamma=0.0
pop_dev = 0.02
measure = Measure(gamma, 1.0)
# to add elements to the measure, e.g.
push_measure!(measure, get_isoperimetric_score, iso )Setup output
output_file_path = joinpath("output", "nc",
"metropolizedForestRECOM_2level_gamma"*string(gamma)*"_iso"*string(iso)*".jsonl.gz")
writer = Writer(measure, constraints, partition, output_file_path)
push_writer!(writer, get_log_spanning_forests)
push_writer!(writer, get_isoperimetric_scores)Run Markov Chain
steps = 1000
output_freq=10
println("starting mcmc")
run_metropolis_hastings!(partition, proposal, measure, steps, rng,
writer=writer, output_freq=output_freq);