Score a BENDL ensemble

This guide uses one committed Colorado VTD BENDL fixture (data/co_vtd_scoring_10000.bendl). The bundle contains:

  • a 10,000-step, ReCom chain;

  • its 3,158-node dual graph and run metadata;

  • a GeoParquet asset with projected VTD geometry, population, and election columns; and

  • fixture provenance, including the two small connectivity repairs made to the VTD seed.

import tempfile
from pathlib import Path

import numpy as np
import pandas as pd

import gerrytools.scoring as gs
from gerrytools.ben import BendlDecoder, read_geoparquet_asset
/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

Open the BENDL resources

BendlDecoder verifies the checksums, reads the embedded graph, and exposes arbitrary assets. read_geoparquet_asset() reads an embedded GeoParquet directly into a GeoDataFrame. For an ordinary Parquet table, use read_parquet_asset() instead. Assignment positions in the BENDL stream follow the embedded graph’s node order, and the node column records the same order explicitly.

data_dir = Path("data")
bundle_path = data_dir / "co_vtd_scoring_10000.bendl"

bundle = BendlDecoder(bundle_path)
graph = bundle.read_graph()
vtds = read_geoparquet_asset(bundle, "co_vtds_2020.parquet")
pd.Series(
    {
        "VTDs": len(vtds),
        "graph edges": graph.number_of_edges(),
        "districts": vtds["assignment"].nunique(),
        "chain samples": bundle.count_samples(),
        "CRS": vtds.crs.to_string(),
        "assets": ", ".join(bundle.asset_names()),
    },
    name="Colorado scoring fixture",
)
VTDs                                                          3158
graph edges                                                   8825
districts                                                        8
chain samples                                                10000
CRS                                                      EPSG:5070
assets           graph.json, co_vtds_2020.parquet, fixture_meta...
Name: Colorado scoring fixture, dtype: object

Configure an evaluator

The evaluator records borrowed graph and geometry sources when it is configured, then snapshots only the requested resources during the first evaluation. Later evaluations reuse that immutable snapshot and the scoring engine. Do not mutate either source after preparation; adding another metric extends the snapshot only with resources that were not already prepared. The active GeoDataFrame geometry column is reserved for geometry-backed metrics. This example registers representative district, plan, region, election, and compactness metrics:

  • Tally combines all requested numeric columns into one engine pass.

  • PolsbyPopper demonstrates geometry-backed compactness.

  • CutEdges demonstrates unweighted counts and shared-perimeter weights.

  • The three region statistics use counties as fixed regions.

  • TallyByRegion produces a county-by-district table with named values.

Each metric owns its optional result_name, so aliases work naturally with add_metrics(...). Compatible metrics still share engine state, but that implementation detail never changes how results are accessed.

graph.nodes(data=True)[0]
{'path': 'vtd:08001001001',
 'county': '08001',
 'total_pop_20': 3990.0,
 'total_vap_20': 3041.0,
 'bvap_20': 69.0,
 'pres_16_dem': 708.0,
 'pres_16_rep': 240.0,
 'pres_20_dem': 900.0,
 'pres_20_rep': 268.0,
 'pres_24_dem': 717.0,
 'pres_24_rep': 287.0,
 'assignment': 7,
 'area': 2007662.5150314635,
 'perimeter': 6716.882051486405}
evaluator = gs.PlanEvaluator(graph, geometry=vtds, node_id_column="node")
evaluator.add_metrics(
    gs.Tally(
        "total_pop_20",
        "total_vap_20",
        "bvap_20",
        "pres_16_dem",
        "pres_16_rep",
        "pres_20_dem",
        "pres_20_rep",
        "pres_24_dem",
        "pres_24_rep",
        result_name="district_totals",
    ),
    gs.PolsbyPopper(),
    gs.Reock(),
    gs.CutEdges(result_name="cut_edge_count"),
    gs.CutEdges(weight_attr="shared_perim", result_name="cut_edge_perimeter"),
    gs.RegionSplits("county", result_name="county_splits"),
    gs.RegionPieces("county", result_name="county_pieces"),
    gs.RegionParts("county", result_name="county_parts"),
    gs.TallyByRegion(
        "county",
        {"population": "total_pop_20", "BVAP": "bvap_20"},
        include_count=True,
        result_name="county_totals",
    ),
)
<gerrytools.scoring.evaluator.PlanEvaluator at 0x...>

Evaluate one assignment with lookup

lookup(i) returns one assignment vector without decoding earlier plans. For a one-off score, call a lowercase function with a GeoDataFrame and either an assignment column or assignment vector. The function constructs a temporary evaluator and returns a pandas object or scalar. Geometry scores use the same pattern, for example gs.polsby_popper(vtds, "assignment"); use a persistent evaluator when the statewide geometry will be reused.

single_assignment = bundle.lookup(0)
direct_population = gs.tally(vtds, single_assignment, columns="total_pop_20")
direct_population_deviations = gs.population_deviations(
    graph,
    single_assignment,
    population_attr="total_pop_20",
)
pd.DataFrame(
    {"population": direct_population, "population deviation": direct_population_deviations}
)
population population deviation
district
7 722666.0 0.001319
6 722566.0 0.001180
5 725391.0 0.005094
3 713675.0 -0.011139
2 729030.0 0.010137
0 720710.0 -0.001391
1 717546.0 -0.005775
4 722130.0 0.000576
single = evaluator.evaluate(single_assignment)
single.metrics
('district_totals',
 'polsby_popper',
 'reock',
 'cut_edge_count',
 'cut_edge_perimeter',
 'county_splits',
 'county_pieces',
 'county_parts',
 'county_totals')
single["district_totals"]
metric total_pop_20 total_vap_20 bvap_20 pres_16_dem pres_16_rep pres_20_dem pres_20_rep pres_24_dem pres_24_rep
district
7 722666.0 536948.0 13948.0 125760.0 130363.0 174252.0 158396.0 165019.0 171188.0
6 722566.0 585116.0 11936.0 185903.0 174124.0 256577.0 191231.0 250351.0 183376.0
5 725391.0 556472.0 63191.0 174777.0 130802.0 233782.0 141997.0 211531.0 140139.0
3 713675.0 537315.0 13275.0 116026.0 215761.0 172032.0 254141.0 180383.0 264377.0
2 729030.0 573796.0 8510.0 140719.0 190901.0 186107.0 219678.0 180480.0 219641.0
0 720710.0 585507.0 56523.0 246189.0 63345.0 315388.0 72263.0 280269.0 75377.0
1 717546.0 582115.0 9439.0 242117.0 121363.0 305131.0 128138.0 294416.0 123688.0
4 722130.0 552307.0 42279.0 107379.0 175825.0 161083.0 198763.0 165712.0 199655.0
single["reock"]
district
7    0.438305
6    0.396708
5    0.247015
3    0.404585
2    0.356426
0    0.182698
1    0.641310
4    0.545158
Name: reock, dtype: float64
pd.Series(
    {name: single[name] for name in ("county_splits", "county_pieces", "county_parts")},
    name="county metrics",
)
county_splits    11
county_pieces    83
county_parts     99
Name: county metrics, dtype: int64

With no graph-attribute names supplied, Polsby-Popper derives its measurements from the aligned GeoDataFrame geometry. Supplying any graph-attribute name would instead select the embedded graph and require every compactness measurement to come from its attributes.

single["polsby_popper"]
district
7    0.217284
6    0.238775
5    0.119050
3    0.314329
2    0.313906
0    0.091214
1    0.374271
4    0.584751
Name: polsby_popper, dtype: float64

Evaluate selected assignments with subsample_indices

subsample_indices decodes only the requested zero-based sample indices. Materialize that small selection before passing it to evaluate_many; sample_ids then assigns meaningful, unique labels to the result rows rather than requiring the caller to relabel each table. Every assignment in a batch must use the same district-label set.

sample_indices = [0, 100, 1_000, 9_999]
selected_assignments = list(bundle.subsample_indices(sample_indices))
selected = evaluator.evaluate_many(
    selected_assignments,
    sample_ids=sample_indices,
)
selected["cut_edge_count"]
sample
0       618
100     678
1000    748
9999    710
Name: cut_edge_count, dtype: int64

Stream the complete 10,000-step chain

evaluate_stream accepts BEN, XBEN, and finalized BENDL input. It writes bounded, Snappy-compressed Parquet batches and a versioned manifest. A missing output directory is created. New score names are added to an existing run; update=True is only needed to replace a name already stored there.

samples counts expanded chain steps. accepted counts encoded frames written to each table. When consecutive assignments repeat, the Parquet repetitions column preserves their multiplicity.

run_root = Path(tempfile.mkdtemp(prefix="gerrytools-scoring-"))
run_dir = run_root / "scores"
run = evaluator.evaluate_stream(bundle_path, run_dir, progress=True)
run.summary, run.metrics
Evaluating ensemble: 0sample [00:00, ?sample/s]
Evaluating ensemble:   0%|          | 0/10000 [00:00<?, ?sample/s]
Evaluating ensemble:   3%|▎         | 256/10000 [00:00<00:04, 2230.83sample/s]
Evaluating ensemble:   5%|▌         | 512/10000 [00:00<00:04, 2249.56sample/s]
Evaluating ensemble:   8%|▊         | 768/10000 [00:00<00:04, 2227.29sample/s]
Evaluating ensemble:  10%|█         | 1024/10000 [00:00<00:04, 1948.86sample/s]
Evaluating ensemble:  13%|█▎        | 1280/10000 [00:00<00:04, 2043.77sample/s]
Evaluating ensemble:  15%|█▌        | 1536/10000 [00:00<00:04, 2084.71sample/s]
Evaluating ensemble:  18%|█▊        | 1792/10000 [00:00<00:03, 2121.94sample/s]
Evaluating ensemble:  20%|██        | 2048/10000 [00:01<00:04, 1935.74sample/s]
Evaluating ensemble:  23%|██▎       | 2304/10000 [00:01<00:03, 2017.58sample/s]
Evaluating ensemble:  26%|██▌       | 2560/10000 [00:01<00:03, 2071.17sample/s]
Evaluating ensemble:  28%|██▊       | 2816/10000 [00:01<00:03, 2120.58sample/s]
Evaluating ensemble:  31%|███       | 3072/10000 [00:01<00:03, 1951.18sample/s]
Evaluating ensemble:  33%|███▎      | 3328/10000 [00:01<00:03, 2041.25sample/s]
Evaluating ensemble:  36%|███▌      | 3584/10000 [00:01<00:03, 2072.56sample/s]
Evaluating ensemble:  38%|███▊      | 3840/10000 [00:01<00:02, 2114.29sample/s]
Evaluating ensemble:  41%|████      | 4096/10000 [00:02<00:03, 1944.43sample/s]
Evaluating ensemble:  44%|████▎     | 4352/10000 [00:02<00:02, 2016.49sample/s]
Evaluating ensemble:  46%|████▌     | 4608/10000 [00:02<00:02, 2082.46sample/s]
Evaluating ensemble:  49%|████▊     | 4864/10000 [00:02<00:02, 2123.60sample/s]
Evaluating ensemble:  51%|█████     | 5120/10000 [00:02<00:02, 1944.47sample/s]
Evaluating ensemble:  54%|█████▍    | 5376/10000 [00:02<00:02, 2023.00sample/s]
Evaluating ensemble:  56%|█████▋    | 5632/10000 [00:02<00:02, 2078.37sample/s]
Evaluating ensemble:  59%|█████▉    | 5888/10000 [00:02<00:01, 2098.69sample/s]
Evaluating ensemble:  61%|██████▏   | 6144/10000 [00:03<00:01, 1960.03sample/s]
Evaluating ensemble:  64%|██████▍   | 6400/10000 [00:03<00:01, 2027.94sample/s]
Evaluating ensemble:  67%|██████▋   | 6656/10000 [00:03<00:01, 2072.26sample/s]
Evaluating ensemble:  69%|██████▉   | 6912/10000 [00:03<00:01, 2124.25sample/s]
Evaluating ensemble:  72%|███████▏  | 7168/10000 [00:03<00:01, 1954.37sample/s]
Evaluating ensemble:  74%|███████▍  | 7424/10000 [00:03<00:01, 2031.60sample/s]
Evaluating ensemble:  77%|███████▋  | 7680/10000 [00:03<00:01, 2075.91sample/s]
Evaluating ensemble:  79%|███████▉  | 7936/10000 [00:03<00:00, 2107.26sample/s]
Evaluating ensemble:  82%|████████▏ | 8192/10000 [00:04<00:00, 1959.99sample/s]
Evaluating ensemble:  84%|████████▍ | 8448/10000 [00:04<00:00, 2039.89sample/s]
Evaluating ensemble:  87%|████████▋ | 8704/10000 [00:04<00:00, 2090.72sample/s]
Evaluating ensemble:  90%|████████▉ | 8960/10000 [00:04<00:00, 2120.48sample/s]
Evaluating ensemble:  92%|█████████▏| 9216/10000 [00:04<00:00, 1958.83sample/s]
Evaluating ensemble:  95%|█████████▍| 9472/10000 [00:04<00:00, 2023.98sample/s]
Evaluating ensemble:  97%|█████████▋| 9728/10000 [00:04<00:00, 2084.63sample/s]
Evaluating ensemble: 100%|█████████▉| 9984/10000 [00:04<00:00, 2135.83sample/s]
Evaluating ensemble: 100%|██████████| 10000/10000 [00:05<00:00, 1977.28sample/s]
(EvaluationSummary(samples=10000, accepted=10000),
 ('district_totals',
  'polsby_popper',
  'reock',
  'cut_edge_count',
  'cut_edge_perimeter',
  'county_splits',
  'county_pieces',
  'county_parts',
  'county_totals'))
run.frames.head()
sample_offset repetitions
accepted
0 0 1
1 1 1
2 2 1
3 3 1
4 4 1

Apply array formulas to streamed tallies

If you usually load a large DataFrame with pd.read_parquet(...), think of run.read("district_totals") as the equivalent operation for one scored metric. It eagerly loads that metric and returns an ordinary pandas Series or DataFrame with meaningful index and column labels. The underlying Parquet column names and table layout remain storage details.

A streamed chain can store a self-loop once with a repetition count instead of writing several identical rows. By default, run.read(...) returns one row per accepted frame. Setting expand_repetitions=True repeats those rows so the sample index covers all 10,000 original chain steps in this example. This is convenient when downstream pandas code expects one unweighted row per step, but the expanded result can use substantially more memory.

As with any large DataFrame load, an eager read must fit in memory. GerryTools warns when the predicted peak reaches 2 GiB and raises EvaluationMemoryError at 8 GiB. If that happens, iterate over run.iter_batches(...); each batch has the same semantic pandas layout without loading the whole metric at once. Use allow_large=True only when the machine deliberately has enough memory. To keep the peak predictable, result reads decode Parquet columns serially, so very wide metrics may load more slowly than a maximally parallel Parquet read.

After loading a metric, .to_numpy() provides the arrays used by scoring.formulas. For a typical sample-by-district DataFrame, the array shape is (samples, districts): formulas operate on every leading sample or batch axis and treat the last axis as districts. Formulas that combine elections expect an array shaped (..., elections, districts).

tallies = run.read("district_totals", expand_repetitions=True)
districts = tallies["total_pop_20"].columns

dem_2020 = tallies["pres_20_dem"].to_numpy()
rep_2020 = tallies["pres_20_rep"].to_numpy()
vote_shares_2020 = gs.formulas.district_vote_shares(dem_2020, rep_2020)
wins_2020 = gs.formulas.district_wins(dem_2020, rep_2020)

pd.concat(
    {
        "Democratic two-party share": pd.DataFrame(
            vote_shares_2020[:5],
            columns=districts,
        ),
        "Democratic win": pd.DataFrame(
            wins_2020[:5],
            columns=districts,
        ),
    },
    axis="columns",
).rename_axis(index="sample", columns=["quantity", "district"])
quantity Democratic two-party share Democratic win
district 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7
sample
0 0.813587 0.704253 0.458634 0.403667 0.447644 0.622126 0.572962 0.523833 True True False False False True True True
1 0.813587 0.704253 0.458634 0.506262 0.447644 0.505746 0.572962 0.523833 True True False True False True True True
2 0.813587 0.704253 0.458634 0.447325 0.507371 0.505746 0.572962 0.523833 True True False False True True True True
3 0.813587 0.704253 0.458634 0.447325 0.507371 0.507775 0.572962 0.522139 True True False False True True True True
4 0.813587 0.704253 0.458634 0.447325 0.574186 0.507775 0.508156 0.522139 True True False False True True True True
partisan_scores_2020 = pd.DataFrame(
    {
        "seats": gs.formulas.seats(dem_2020, rep_2020),
        "overall_vote_share": gs.formulas.overall_vote_share(dem_2020, rep_2020),
        "efficiency_gap": gs.formulas.efficiency_gap(dem_2020, rep_2020),
        "simplified_efficiency_gap": gs.formulas.simplified_efficiency_gap(
            dem_2020,
            rep_2020,
        ),
        "mean_median": gs.formulas.mean_median(dem_2020, rep_2020),
        "partisan_bias_equal": gs.formulas.partisan_bias(
            dem_2020,
            rep_2020,
            turnout_model="equal",
        ),
        "partisan_bias_observed": gs.formulas.partisan_bias(
            dem_2020,
            rep_2020,
            turnout_model="observed",
        ),
        "partisan_gini_equal": gs.formulas.partisan_gini(
            dem_2020,
            rep_2020,
            turnout_model="equal",
        ),
        "partisan_gini_observed": gs.formulas.partisan_gini(
            dem_2020,
            rep_2020,
            turnout_model="observed",
        ),
    }
)
partisan_scores_2020.agg(["mean", "std", "min", "max"]).T
mean std min max
seats 5.141000 6.311567e-01 4.000000 7.000000
overall_vote_share 0.569383 1.110279e-16 0.569383 0.569383
efficiency_gap 0.008446 8.127099e-02 -0.164526 0.250681
simplified_efficiency_gap 0.003859 7.889459e-02 -0.138766 0.236234
mean_median -0.007717 2.029084e-02 -0.065593 0.057738
partisan_bias_equal -0.018775 7.366252e-02 -0.250000 0.125000
partisan_bias_observed -0.018038 7.324102e-02 -0.250000 0.125000
partisan_gini_equal 0.031613 1.400150e-02 0.001302 0.084842
partisan_gini_observed 0.031715 1.393516e-02 0.001708 0.084842

The fixture contains three presidential elections. Stacking them gives arrays with shape (plans, elections, districts) for cross-election summaries.

years = (2016, 2020, 2024)
dem_elections = np.stack(
    [tallies[f"pres_{year % 100:02d}_dem"].to_numpy() for year in years],
    axis=1,
)
rep_elections = np.stack(
    [tallies[f"pres_{year % 100:02d}_rep"].to_numpy() for year in years],
    axis=1,
)

wins_by_district = gs.formulas.party_wins_by_district(
    dem_elections,
    rep_elections,
)
pd.DataFrame(wins_by_district[:5], columns=districts)
district 0 1 2 3 4 5 6 7
0 3 3 0 0 0 3 3 1
1 3 3 0 2 0 1 3 1
2 3 3 0 0 2 1 3 1
3 3 3 0 0 2 1 3 1
4 3 3 0 0 3 1 2 1
cross_election_scores = pd.DataFrame(
    {
        "competitive_contests": gs.formulas.competitive_contests(
            dem_elections,
            rep_elections,
            vote_share_margin=0.05,
        ),
        "swing_districts": gs.formulas.swing_districts(
            dem_elections,
            rep_elections,
        ),
        "democratic_districts": gs.formulas.party_districts(
            dem_elections,
            rep_elections,
        ),
        "republican_districts": gs.formulas.opposition_party_districts(
            dem_elections,
            rep_elections,
        ),
        "aggregate_democratic_seats": gs.formulas.aggregate_seats(
            dem_elections,
            rep_elections,
        ),
        "mean_signed_seat_vote_gap": gs.formulas.mean_signed_seat_vote_gap(
            dem_elections,
            rep_elections,
        ),
        "mean_absolute_seat_vote_gap": gs.formulas.mean_absolute_seat_vote_gap(
            dem_elections,
            rep_elections,
        ),
    }
)
cross_election_scores.agg(["mean", "std", "min", "max"]).T
mean std min max
competitive_contests 7.309400 2.081131 2.000000 15.000000
swing_districts 0.901200 0.740333 0.000000 4.000000
democratic_districts 4.262600 0.604878 3.000000 6.000000
republican_districts 2.836200 0.631514 1.000000 4.000000
aggregate_democratic_seats 14.399400 1.527846 11.000000 19.000000
mean_signed_seat_vote_gap 0.049081 0.063660 -0.092561 0.240773
mean_absolute_seat_vote_gap 0.079987 0.042679 0.046305 0.240773

Derive population, demographic, and compactness scores

Population and demographic functions use district tallies from the same streamed table. Schwartzberg compactness is derived from the scoring-engine Polsby-Popper output.

population = tallies["total_pop_20"].to_numpy()
voting_age_population = tallies["total_vap_20"].to_numpy()
black_voting_age_population = tallies["bvap_20"].to_numpy()

population_deviations = gs.formulas.population_deviations(population)
bvap_share = gs.formulas.demographic_shares(
    black_voting_age_population,
    voting_age_population,
)
population_scores = pd.DataFrame(
    {
        "max_absolute_deviation": gs.formulas.max_absolute_population_deviation(
            population,
            relative_to_ideal=True,
        ),
        "maximum_deviation": gs.formulas.max_population_deviation(
            population,
            relative_to_ideal=True,
        ),
        "districts_above_40_percent_BVAP": gs.formulas.districts_above_threshold(
            black_voting_age_population,
            voting_age_population,
            threshold=0.4,
        ),
    }
)
population_scores.agg(["mean", "std", "min", "max"]).T
mean std min max
max_absolute_deviation 0.044438 0.005187 0.011139 0.049998
maximum_deviation 0.082233 0.011228 0.019261 0.099901
districts_above_40_percent_BVAP 0.000000 0.000000 0.000000 0.000000
polsby = run.read("polsby_popper", expand_repetitions=True).to_numpy()
schwartzberg = gs.formulas.schwartzberg(polsby)

pd.concat(
    {
        "population deviation": pd.DataFrame(
            population_deviations[:5],
            columns=districts,
        ),
        "BVAP share": pd.DataFrame(
            bvap_share[:5],
            columns=districts,
        ),
        "Schwartzberg": pd.DataFrame(
            schwartzberg[:5],
            columns=districts,
        ),
    },
    axis="columns",
)
population deviation BVAP share Schwartzberg
district 0 1 2 3 4 5 6 7 0 1 ... 6 7 0 1 2 3 4 5 6 7
0 -0.001391 -0.005775 0.010137 -0.011139 0.000576 0.005094 0.00118 0.001319 0.096537 0.016215 ... 0.020399 0.025976 3.311071 1.634582 1.784843 1.783643 1.307720 2.898250 2.046469 2.145292
1 -0.001391 -0.005775 0.010137 -0.021843 0.000576 0.015798 0.00118 0.001319 0.096537 0.016215 ... 0.020399 0.025976 3.311071 1.634582 1.784843 1.893948 1.307720 1.777326 2.046469 2.145292
2 -0.001391 -0.005775 0.010137 0.013213 -0.034479 0.015798 0.00118 0.001319 0.096537 0.016215 ... 0.020399 0.025976 3.311071 1.634582 1.784843 2.503794 2.622338 1.777326 2.046469 2.145292
3 -0.001391 -0.005775 0.010137 0.013213 -0.034479 0.017385 0.00118 -0.000268 0.096537 0.016215 ... 0.020399 0.115688 3.311071 1.634582 1.784843 2.503794 2.622338 2.733609 2.046469 1.812951
4 -0.001391 -0.005775 0.010137 0.013213 -0.025059 0.017385 -0.00824 -0.000268 0.096537 0.016215 ... 0.027823 0.115688 3.311071 1.634582 1.784843 2.503794 2.143320 2.733609 2.695370 1.812951

5 rows × 24 columns

Inspect region-by-district tallies

TallyByRegion exposes values for each fixed region and proposed district. For one plan, regions form the row index; the first column level contains meaningful value names, and the second contains ordered district labels.

county_totals = single["county_totals"]
display(county_totals.head())
metric count population BVAP
district 7 6 5 3 2 0 1 4 7 6 ... 1 4 7 6 5 3 2 0 1 4
county
08001 230 2 14 9 0 0 0 0 456673.0 2892.0 ... 0.0 0.0 9896.0 45.0 5724.0 160.0 0.0 0.0 0.0 0.0
08003 0 0 0 0 8 0 0 0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 277.0 0.0 0.0 0.0
08005 0 0 372 22 0 9 0 0 0.0 0.0 ... 0.0 0.0 0.0 0.0 56525.0 2721.0 0.0 399.0 0.0 0.0
08007 0 0 0 0 8 0 0 0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 105.0 0.0 0.0 0.0
08009 0 0 0 9 0 0 0 0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 39.0 0.0 0.0 0.0 0.0

5 rows × 24 columns

Save and compare several runs

Give each assignment stream its own score directory. evaluate_stream creates the directory and returns its open EnsembleEvalResult. This tutorial has one BENDL fixture, so the example reuses it for both paths; an actual comparison would pass a different BEN or BENDL file to each call.

stats_dir = run_root / "stats"
run_1_path = stats_dir / "ben_file_1_scores"
run_2_path = stats_dir / "ben_file_2_scores"

totals = gs.PlanEvaluator(graph).add_metric(
    gs.Tally(
        "total_pop_20",
        "total_vap_20",
        result_name="district_totals",
    )
)
run_1_scores = totals.evaluate_stream(bundle_path, run_1_path)
run_2_scores = totals.evaluate_stream(bundle_path, run_2_path)
run_1_scores.summary, run_2_scores.summary
(EvaluationSummary(samples=10000, accepted=10000),
 EvaluationSummary(samples=10000, accepted=10000))

The returned object is already open. Use read() to extract a saved logical score; callers do not need to find or combine the physical Parquet files themselves.

run_1_totals = run_1_scores.read("district_totals")
run_2_totals = run_2_scores.read("district_totals")
run_1_totals.head()
metric total_pop_20 total_vap_20
district 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7
accepted
0 720710.0 717546.0 729030.0 713675.0 722130.0 725391.0 722566.0 722666.0 585507.0 582115.0 573796.0 537315.0 552307.0 556472.0 585116.0 536948.0
1 720710.0 717546.0 729030.0 705950.0 722130.0 733116.0 722566.0 722666.0 585507.0 582115.0 573796.0 537128.0 552307.0 556659.0 585116.0 536948.0
2 720710.0 717546.0 729030.0 731250.0 696830.0 733116.0 722566.0 722666.0 585507.0 582115.0 573796.0 560099.0 529336.0 556659.0 585116.0 536948.0
3 720710.0 717546.0 729030.0 731250.0 696830.0 734261.0 722566.0 721521.0 585507.0 582115.0 573796.0 560099.0 529336.0 549047.0 585116.0 544560.0
4 720710.0 717546.0 729030.0 731250.0 703629.0 734261.0 715767.0 721521.0 585507.0 582115.0 573796.0 560099.0 570189.0 549047.0 544263.0 544560.0

A new result name is added without a flag. Reopen the run afterward so the object sees the new manifest, then read both the old totals and the new cut-edge score.

edges = gs.PlanEvaluator(graph).add_metric(gs.CutEdges())
edges.evaluate_stream(bundle_path, run_1_path)
edges.evaluate_stream(bundle_path, run_2_path)

run_1_scores = gs.EnsembleEvalResult.open(run_1_path)
run_2_scores = gs.EnsembleEvalResult.open(run_2_path)
old_totals = run_1_scores.read("district_totals")
new_cut_edges = run_1_scores.read("cut_edges")
display(old_totals.head())
new_cut_edges.head()
metric total_pop_20 total_vap_20
district 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7
accepted
0 720710.0 717546.0 729030.0 713675.0 722130.0 725391.0 722566.0 722666.0 585507.0 582115.0 573796.0 537315.0 552307.0 556472.0 585116.0 536948.0
1 720710.0 717546.0 729030.0 705950.0 722130.0 733116.0 722566.0 722666.0 585507.0 582115.0 573796.0 537128.0 552307.0 556659.0 585116.0 536948.0
2 720710.0 717546.0 729030.0 731250.0 696830.0 733116.0 722566.0 722666.0 585507.0 582115.0 573796.0 560099.0 529336.0 556659.0 585116.0 536948.0
3 720710.0 717546.0 729030.0 731250.0 696830.0 734261.0 722566.0 721521.0 585507.0 582115.0 573796.0 560099.0 529336.0 549047.0 585116.0 544560.0
4 720710.0 717546.0 729030.0 731250.0 703629.0 734261.0 715767.0 721521.0 585507.0 582115.0 573796.0 560099.0 570189.0 549047.0 544263.0 544560.0
accepted
0    618
1    603
2    635
3    658
4    662
Name: cut_edges, dtype: int64

The directory name is included in every filename, so a copied Parquet file still identifies its run:

stats/
├── ben_file_1_scores/
│   ├── manifest__ben_file_1_scores.json
│   ├── cut_edges__ben_file_1_scores.parquet
│   └── district_totals/
│       ├── total_pop_20_tallies__ben_file_1_scores.parquet
│       └── total_vap_20_tallies__ben_file_1_scores.parquet
└── ben_file_2_scores/
    ├── manifest__ben_file_2_scores.json
    ├── cut_edges__ben_file_2_scores.parquet
    └── district_totals/
        ├── total_pop_20_tallies__ben_file_2_scores.parquet
        └── total_vap_20_tallies__ben_file_2_scores.parquet

Use update=True only when replacing a result with the same name. Scores not registered on the replacement evaluator stay as they are.

replacement = gs.PlanEvaluator(graph).add_metric(
    gs.Tally("total_pop_20", result_name="district_totals")
)
replacement.evaluate_stream(bundle_path, run_1_path, update=True)

run_1_scores = gs.EnsembleEvalResult.open(run_1_path)
updated_totals = run_1_scores.read("district_totals")
preserved_cut_edges = run_1_scores.read("cut_edges")
updated_totals.head(), preserved_cut_edges.head()
(district         0         1         2         3         4         5  \
 accepted                                                               
 0         720710.0  717546.0  729030.0  713675.0  722130.0  725391.0   
 1         720710.0  717546.0  729030.0  705950.0  722130.0  733116.0   
 2         720710.0  717546.0  729030.0  731250.0  696830.0  733116.0   
 3         720710.0  717546.0  729030.0  731250.0  696830.0  734261.0   
 4         720710.0  717546.0  729030.0  731250.0  703629.0  734261.0   
 
 district         6         7  
 accepted                      
 0         722566.0  722666.0  
 1         722566.0  722666.0  
 2         722566.0  722666.0  
 3         722566.0  721521.0  
 4         715767.0  721521.0  ,
 accepted
 0    618
 1    603
 2    635
 3    658
 4    662
 Name: cut_edges, dtype: int64)

evaluate_stream checks that the sample count, accepted-frame count, and district IDs fit the existing directory. You are responsible for passing the same assignment stream when adding or replacing scores. Reopen an EnsembleEvalResult after changing its directory; an object opened earlier keeps the old manifest in memory.

Result contracts

evaluate returns PlanEvalResult; indexing it by metric name returns a scalar, Series, or DataFrame with semantic labels. evaluate_many returns ManyPlanEvalResult; its first axis uses sample_ids when supplied. Region results use (sample, region) rows for many plans and region rows for one plan, with (metric, district) columns in both cases. array(name) is available when a canonical immutable NumPy view is preferable.

evaluate_stream returns EnsembleEvalResult after writing Parquet score tables plus a JSON manifest. Tallies and region tallies use one table per requested attribute; other metrics use one table. Every filename ends with the run-directory name. frames exposes accepted-frame offsets and repetition counts. read() restores the same logical pandas shapes at either accepted-frame or expanded-sample resolution, while raw() remains available for physical Parquet access. Large results can instead be processed with iter_batches(), iter_raw_batches(), or iter_frame_batches() without materializing a whole table. Array formulas deliberately remain separate, so they work with in-memory results, streamed tables, or arrays produced elsewhere.

Sequence assignments must follow graph-node order. Mapping assignments are aligned by node identifier; both forms are captured by the exported gerrytools.scoring.Assignment alias. Geometry must be projected for area and distance metrics. Consult each metric or formula docstring for its formula, tie convention, turnout model, and literature references.