PlanEvaluator

Browse tutorial data

PlanEvaluator is the reusable scoring interface. It aligns a graph with an optional GeoDataFrame, prepares the resources required by a collection of metrics, and then applies those metrics to one plan, several plans, GerryChain partitions, or an encoded ensemble.

This guide uses one 6-by-6 grid throughout. The example is small enough to inspect directly, but it includes the same pieces as a precinct-level analysis: a dual graph, projected polygons, population and election columns, region labels, and two districting plans.

import tempfile
from math import pi
from pathlib import Path
from time import perf_counter

import geopandas as gpd
import matplotlib.pyplot as plt
import networkx as nx
import pandas as pd
from binary_ensemble import BenEncoder
from gerrychain import Partition
from shapely.geometry import box

import gerrytools.scoring as gs
from gerrytools.ben import BendlDecoder, read_geoparquet_asset

# NetworkX preserves this node order when the coordinate labels are replaced by integers.
graph = nx.convert_node_labels_to_integers(
    nx.grid_2d_graph(6, 6),
    label_attribute="grid_position",
)

records = []
for node, data in graph.nodes(data=True):
    row, column = data.pop("grid_position")
    records.append(
        {
            "node": node,
            "row": row,
            "column": column,
            "population": 90 + 2 * row + column,
            "democratic": 35 + 7 * column + 2 * row,
            "republican": 65 - 2 * column,
            "county": "west" if column < 3 else "east",
            "district": (row // 3) * 2 + column // 3,
            "geometry": box(column, row, column + 1, row + 1),
        }
    )

units = gpd.GeoDataFrame(records, crs="EPSG:3857")

for unit in units.itertuples():
    graph.nodes[unit.node].update(
        population=unit.population,
        democratic=unit.democratic,
        republican=unit.republican,
        county=unit.county,
        district=unit.district,
        area=1.0,
        boundary_perim=float((unit.row in (0, 5)) + (unit.column in (0, 5))),
    )

nx.set_edge_attributes(graph, 1.0, "shared_perim")
units.head()

/home/docs/checkouts/readthedocs.org/user_builds/gerrytools/envs/latest/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

	node	row	column	population	democratic	republican	county	district	geometry
0	0	0	0	90	35	65	west	0	POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))
1	1	0	1	91	42	63	west	0	POLYGON ((2 0, 2 1, 1 1, 1 0, 2 0))
2	2	0	2	92	49	61	west	0	POLYGON ((3 0, 3 1, 2 1, 2 0, 3 0))
3	3	0	3	93	56	59	east	1	POLYGON ((4 0, 4 1, 3 1, 3 0, 4 0))
4	4	0	4	94	63	57	east	1	POLYGON ((5 0, 5 1, 4 1, 4 0, 5 0))

Align graph and geometry

The graph defines the unit order used by assignment vectors and supplies the topology for graph metrics. The GeoDataFrame supplies projected geometry and, when present, is authoritative for ordinary node and region columns. node_id_column="node" makes the alignment explicit; without it, the evaluator would match graph nodes against the GeoDataFrame index.

The graph also carries area, boundary_perim, and shared_perim attributes. They are not needed for the first evaluator, which derives compactness measurements from units, but they let us compare the two preparation paths later.

Assignment order

A mapping is reordered by graph node identifier. A sequence is interpreted in the graph’s construction-time node order. Prefer a mapping when assignments come from a table whose row order may have changed.

from gerrytools.plotting import GeoPlot

gp = GeoPlot(units)
gp.add_districting_plan_layer("district", edgecolor="black")
gp.show(pad_inches=-0.15, bbox_inches="tight")

Prepare a reusable metric collection

add_geometry() must be called before the first metric because geometry changes which resources a metric can request. add_metrics() validates the whole registration batch before adding it. Column existence, numeric values, topology, and geometry contracts are checked later during lazy preparation, when the evaluator first scores a plan.

The metrics below exercise several result shapes. Tally produces multiple values for every district, Seats and CutEdges produce one value for the plan, and the two compactness metrics produce one value per district. result_name= controls the result key without changing the metric’s calculation.

evaluator = (
    gs.PlanEvaluator(graph)
    .add_geometry(units, node_id_column="node")
    .add_metrics(
        gs.Tally(
            "population",
            "democratic",
            "republican",
            result_name="district_totals",
        ),
        gs.Seats(
            "democratic",
            "republican",
            result_name="democratic_seats",
        ),
        gs.PolsbyPopper(),
        gs.Reock(),
        gs.CutEdges(result_name="cut_edges"),
    )
)
evaluator.metrics

('district_totals', 'democratic_seats', 'polsby_popper', 'reock', 'cut_edges')

The metrics property reports logical result names in evaluation order. Registration does not perform any scoring. On the first call to evaluate(), evaluate_many(), or evaluate_stream(), the evaluator aligns the requested columns, snapshots them, prepares shared Rust state, and checks metric-specific contracts. Later calls reuse that prepared state.

Evaluate one plan

evaluate() accepts an assignment mapping, an assignment sequence, or a GerryChain Partition. Here the mapping is built from the explicit node column, so its meaning does not depend on the current GeoDataFrame row order.

PlanEvalResult is a read-only mapping. A district metric returns a Series, a multi-column district metric returns a DataFrame, and a plan metric returns a scalar. This keeps each result in its natural pandas shape instead of forcing unrelated metrics into one rectangular table.

assignment = units.set_index("node")["district"].to_dict()
result = evaluator.evaluate(assignment)

result["district_totals"].join(result["polsby_popper"].rename("Polsby-Popper")).join(
    result["reock"].rename("Reock")
)

	population	democratic	republican	Polsby-Popper	Reock
district
0	837.0	396.0	567.0	0.785398	0.63662
1	864.0	585.0	513.0	0.785398	0.63662
2	891.0	450.0	567.0	0.785398	0.63662
3	918.0	639.0	513.0	0.785398	0.63662

pd.Series(
    {
        "Democratic seats": result["democratic_seats"],
        "cut edges": result["cut_edges"],
        "Polsby-Popper array shape": result.array("polsby_popper").shape,
    }
)

Democratic seats                  2
cut edges                        12
Polsby-Popper array shape    (1, 4)
dtype: object

array(name) exposes the same prepared result as an immutable NumPy view. It is useful when a later calculation expects arrays rather than pandas labels; ordinary inspection is usually clearer through result[name].

Compare several plans in memory

evaluate_many() is intended for a manageable collection of selected plans. The second plan below uses four connected, stepped districts with the same nine-unit district sizes as the quadrant plan. Plotting both assignments makes the compactness results easier to interpret.

stepped = {unit.node: (unit.row * 6 + unit.column) // 9 for unit in units.itertuples()}

figure, axes = plt.subplots(1, 2, figsize=(8, 4))
for ax, plan in zip(axes, (assignment, stepped), strict=True):
    units["district"] = pd.Series(plan)
    gp = GeoPlot(units)
    gp.add_districting_plan_layer("district", edgecolor="white")
    gp.bind_to_ax(ax)

comparison = evaluator.evaluate_many(
    [assignment, stepped],
    sample_ids=["quadrants", "stepped"],
    track_uniqueness=True,
)

pd.DataFrame(
    {
        "Democratic seats": comparison["democratic_seats"],
        "cut edges": comparison["cut_edges"],
        "mean Polsby-Popper": comparison["polsby_popper"].mean(axis=1),
        "mean Reock": comparison["reock"].mean(axis=1),
    }
)

	Democratic seats	cut edges	mean Polsby-Popper	mean Reock
sample
quadrants	2	12	0.785398	0.636620
stepped	2	20	0.441786	0.286479

The sample IDs become the result index. District-valued metrics become sample-by-district DataFrames, while plan-valued metrics become Series. track_uniqueness=True also records counts that ignore district labels and district ordering. It does not filter repeated plans; it adds summary information.

evaluate_many() materializes its assignments and results in memory. That is convenient for a small comparison like this one. Use evaluate_stream() for a long BEN, XBEN, or BENDL recording.

comparison.summary

EvaluationSummary(samples=2, accepted=2, unique_plans=2, unique_districts=8)

Choose geometry or graph compactness

With geometry attached and no *_attr arguments, Polsby-Popper and Schwartzberg derive unit areas, total perimeters, shared boundaries, and rook adjacency from the aligned GeoDataFrame. Reock and the hull-based metrics likewise prepare their geometric resources from the polygons. This path has a larger one-time preparation cost, but it keeps the measurements tied to one projected geometry source.

Supplying any compactness *_attr argument selects graph-backed scoring. Every compactness measurement then comes from graph attributes, including attributes whose names retain their defaults. This avoids deriving measurements from polygons, but all node and edge attributes must have been computed from compatible geometry in the same projected CRS.

graph_compactness = (
    gs.PlanEvaluator(graph, geometry=units, node_id_column="node")
    .add_metric(
        gs.PolsbyPopper(
            area_attr="area",
            boundary_perimeter_attr="boundary_perim",
            shared_perimeter_attr="shared_perim",
        )
    )
    .evaluate(assignment)
)

pd.concat(
    {
        "derived from geometry": result["polsby_popper"],
        "read from graph attributes": graph_compactness["polsby_popper"],
    },
    axis=1,
)

	derived from geometry	read from graph attributes
district
0	0.785398	0.785398
1	0.785398	0.785398
2	0.785398	0.785398
3	0.785398	0.785398

The values agree because this example’s graph attributes were computed from the same unit squares. For total-perimeter mode, perimeter_attr names each unit’s complete perimeter. For boundary mode, boundary_perimeter_attr names only the exterior portion and shared_perimeter_attr supplies the edge lengths needed to reconstruct the total.

Compare Colorado VTD geometry and graph attributes

The grid makes the source-selection rules visible, but it does not show the preparation cost of a real precinct layer. The Colorado fixture used in the BENDL tutorial contains 3,158 VTDs, an embedded dual graph, and an enacted congressional assignment. After downloading the tutorial data, place the data directory beside this notebook.

The two evaluators below score the same assignment. One receives the VTD GeoDataFrame and derives all Polsby-Popper measurements during preparation. The other reads area, exterior perimeter, and shared perimeter directly from graph attributes.

bundle = BendlDecoder(Path("data/co_vtd_scoring_10000.bendl"))
co_graph = bundle.read_graph()
co_vtds = read_geoparquet_asset(bundle, "co_vtds_2020.parquet")
co_assignment = co_vtds.set_index("node")["assignment"].to_dict()

pd.Series(
    {
        "VTDs": len(co_vtds),
        "graph edges": co_graph.number_of_edges(),
        "districts": co_vtds["assignment"].nunique(),
        "CRS": co_vtds.crs.to_string(),
    },
    name="Colorado fixture",
)

VTDs                3158
graph edges         8825
districts              8
CRS            EPSG:5070
Name: Colorado fixture, dtype: object

Prepare a consistent graph representation

Graph-backed compactness assumes that every measurement came from compatible geometry in one projected CRS. The fixture’s original perimeter attributes are retained for provenance, so this example writes a separate set with _5070 names from the embedded EPSG:5070 VTD polygons.

For each graph edge, the shared perimeter is the length of the common polygon boundary. A unit’s exterior perimeter is its complete polygon perimeter minus the shared lengths on all incident edges. Values within floating-point tolerance of zero are clamped to zero. In production these attributes would normally be prepared once with the graph and reused by every analysis; this cell keeps their origin inspectable.

co_graph_attributes = co_graph.copy()
co_geometry = co_vtds.set_index("node").geometry

nx.set_node_attributes(
    co_graph_attributes,
    co_geometry.area.to_dict(),
    "area_5070",
)

shared_perimeter = {
    (left, right): co_geometry[left].boundary.intersection(co_geometry[right].boundary).length
    for left, right in co_graph_attributes.edges
}
nx.set_edge_attributes(
    co_graph_attributes,
    shared_perimeter,
    "shared_perim_5070",
)

incident_shared = dict.fromkeys(co_graph_attributes, 0.0)
for (left, right), length in shared_perimeter.items():
    incident_shared[left] += length
    incident_shared[right] += length

boundary_perimeter = {
    node: co_geometry[node].length - incident_shared[node] for node in co_graph_attributes
}
assert min(boundary_perimeter.values()) > -1e-6
boundary_perimeter = {node: max(0.0, length) for node, length in boundary_perimeter.items()}
nx.set_node_attributes(
    co_graph_attributes,
    boundary_perimeter,
    "boundary_perim_5070",
)

Measure first and repeated evaluations

Preparation is lazy, so timing construction alone would miss the expensive work. The first geometry-backed evaluation decodes polygons, computes unit measurements and shared boundaries, and prepares the scoring state. The second evaluation reuses that state and measures only plan normalization and scoring.

The graph-backed evaluator still validates and snapshots its attributes on the first call, but it does not reconstruct measurements from polygon geometry. Exact timings depend on the machine; the important comparison is between the first and repeated calls for the same evaluator. The graph attribute construction above is intentionally outside these timings because persisted graph attributes are an input to that evaluator, not part of its preparation.

co_geometry_evaluator = gs.PlanEvaluator(
    co_graph, geometry=co_vtds, node_id_column="node"
).add_metric(gs.PolsbyPopper())
co_graph_evaluator = gs.PlanEvaluator(co_graph_attributes).add_metric(
    gs.PolsbyPopper(
        area_attr="area_5070",
        boundary_perimeter_attr="boundary_perim_5070",
        shared_perimeter_attr="shared_perim_5070",
    )
)


def timed_evaluation(evaluator, plan):
    start = perf_counter()
    result = evaluator.evaluate(plan)
    return result, perf_counter() - start


co_geometry_result, geometry_first = timed_evaluation(
    co_geometry_evaluator,
    co_assignment,
)
_, geometry_repeated = timed_evaluation(
    co_geometry_evaluator,
    co_assignment,
)
co_graph_result, graph_first = timed_evaluation(
    co_graph_evaluator,
    co_assignment,
)
_, graph_repeated = timed_evaluation(
    co_graph_evaluator,
    co_assignment,
)

pd.DataFrame(
    {
        "first evaluation (seconds)": [geometry_first, graph_first],
        "repeated evaluation (seconds)": [geometry_repeated, graph_repeated],
    },
    index=["VTD geometry", "graph attributes"],
).rename_axis("measurement source")

	first evaluation (seconds)	repeated evaluation (seconds)
measurement source
VTD geometry	2.579989	0.001599
graph attributes	0.031458	0.001649

Inspect the score components

Polsby-Popper is \(4\pi A/P^2\), where \(A\) is district area and \(P\) is district perimeter. The graph representation makes those aggregates easy to inspect. District area is the sum of its unit areas. District perimeter combines unit exterior boundaries with the shared boundary of every graph edge cut by the plan; a cut edge contributes its length to the perimeter of both adjacent districts.

The table recomputes the formula from those components and places both evaluator results beside it. Small floating-point differences are expected because the geometry path prepares its own boundary representation, but the two versions should agree at the displayed precision.

district_area = dict.fromkeys(co_assignment.values(), 0.0)
district_perimeter = dict.fromkeys(co_assignment.values(), 0.0)

for node, district in co_assignment.items():
    district_area[district] += co_graph_attributes.nodes[node]["area_5070"]
    district_perimeter[district] += co_graph_attributes.nodes[node]["boundary_perim_5070"]

for (left, right), length in shared_perimeter.items():
    left_district = co_assignment[left]
    right_district = co_assignment[right]
    if left_district != right_district:
        district_perimeter[left_district] += length
        district_perimeter[right_district] += length

co_components = pd.DataFrame(
    {
        "area": pd.Series(district_area),
        "perimeter": pd.Series(district_perimeter),
        "geometry score": co_geometry_result["polsby_popper"],
        "graph score": co_graph_result["polsby_popper"],
    }
).sort_index()
co_components["formula check"] = 4 * pi * co_components["area"] / co_components["perimeter"].pow(2)
co_components

	area	perimeter	geometry score	graph score	formula check
0	4.029343e+08	2.356082e+05	0.091214	0.091214	0.091214
1	3.002974e+10	1.004124e+06	0.374271	0.374271	0.374271
2	1.296021e+11	2.277776e+06	0.313906	0.313906	0.313906
3	8.368600e+10	1.829108e+06	0.314329	0.314329	0.314329
4	3.813207e+09	2.862627e+05	0.584751	0.584751	0.584751
5	8.858860e+08	3.057946e+05	0.119050	0.119050	0.119050
6	1.866913e+10	9.912245e+05	0.238775	0.238775	0.238775
7	2.515528e+09	3.814221e+05	0.217284	0.217284	0.217284

The geometry cost is paid once per evaluator, not once per plan. Reusing co_geometry_evaluator for additional enacted maps or an ensemble keeps the prepared VTD measurements and scoring state. Adding a compatible metric can rebuild the engine while retaining already prepared geometry resources; constructing a new evaluator starts preparation from scratch.

Reuse the evaluator as GerryChain updaters

to_updaters() returns a mapping accepted directly by Partition. The first GerryTools metric requested from a partition evaluates the complete registered collection. GerryChain caches that combined result on the partition, so subsequent GerryTools metrics do not rescore it.

The updater mapping is a snapshot of the metric names registered when to_updaters() is called. Call it again if metrics are added later.

partition = Partition(
    graph,
    assignment,
    updaters=evaluator.to_updaters(),
)

pd.Series(
    {
        "Democratic seats": partition["democratic_seats"],
        "cut edges": partition["cut_edges"],
        "districts": len(partition["district_totals"]),
    }
)

Democratic seats     2
cut edges           12
districts            4
dtype: int64

Extend a prepared evaluator

add_metric() registers one additional metric and returns the same evaluator. Adding a metric after evaluation invalidates the scoring engine; the next evaluation rebuilds the expanded metric collection. Compatible resource snapshots are retained, so already prepared geometry and columns do not need to be aligned again.

Polsby-Popper and Schwartzberg share one area-and-perimeter state when their four *_attr arguments match. Schwartzberg is therefore a useful example of extending the evaluator without requesting a new geometric resource.

evaluator.add_metric(gs.Schwartzberg())
extended = evaluator.evaluate(assignment)

pd.concat(
    [
        extended["polsby_popper"],
        extended["schwartzberg"],
    ],
    axis=1,
)

	polsby_popper	schwartzberg
district
0	0.785398	1.128379
1	0.785398	1.128379
2	0.785398	1.128379
3	0.785398	1.128379

Stream an encoded ensemble

evaluate_stream(source, output_dir) reads BEN, XBEN, or finalized BENDL assignments in native batches and writes a Parquet result directory. A missing directory is created. New score names are added to an existing run, while matching names require update=True. batch_size controls assignment frames per engine batch, while max_samples limits logical samples after encoded repetitions are expanded.

A BENDL graph can verify exact node order. Raw BEN and XBEN files do not carry graph metadata, so their vectors must already follow the evaluator’s graph-node order. The tiny BEN file below uses the same two plans evaluated above.

temporary = tempfile.TemporaryDirectory()
source = Path(temporary.name) / "plans.ben"
output_dir = Path(temporary.name) / "scores"

with BenEncoder(source, variant="standard") as stream:
    stream.write([assignment[node] for node in graph.nodes])
    stream.write([stepped[node] for node in graph.nodes])

run = evaluator.evaluate_stream(
    source,
    output_dir=output_dir,
    batch_size=2,
    track_uniqueness=True,
)

pd.Series(
    {
        "samples": run.summary.samples,
        "accepted frames": run.summary.accepted,
        "unique plans": run.summary.unique_plans,
        "metrics": ", ".join(run.metrics),
    }
)

samples                                                            2
accepted frames                                                    2
unique plans                                                       2
metrics            district_totals, democratic_seats, polsby_popp...
dtype: object

run.read("polsby_popper")

district	0	1	2	3
accepted
0	0.785398	0.785398	0.785398	0.785398
1	0.441786	0.441786	0.441786	0.441786

EnsembleEvalResult.read() reconstructs the same logical pandas shape used by in-memory evaluation. Pass return_type="series" or return_type="dataframe" to request that shape and narrow the static return type. Series and one-column DataFrame results convert between shapes. For a large result, iter_batches() yields bounded semantic batches instead. Repeated encoded frames remain compressed unless expand_repetitions=True is requested.

The temporary directory remains available while this notebook kernel is running, so run, its manifest, and its Parquet tables can be inspected in additional cells.

Related

• Overview

• Convenience Functions

• Scoring BENDL files

• Working with GerryChain

• Scoring API