Add BVAP-biased ReCom

BVAP-biased ReCom.ipynb

@@ -0,0 +1,178 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "import matplotlib.pyplot as plt\n",
+    "from gerrychain import (GeographicPartition, Partition, Graph, MarkovChain,\n",
+    "                        proposals, updaters, constraints, accept, Election)\n",
+    "from gerrychain.proposals import recom\n",
+    "from gerrychain.metrics import mean_median\n",
+    "from bvap_recom import bvap_biased_recom\n",
+    "from functools import partial\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt; plt.style.use('ggplot')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "graph = Graph.from_file(\"data/VA_precincts/VA_precincts.shp\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for node in graph.nodes:\n",
+    "    total_pop = graph.nodes[node]['TOTPOP']\n",
+    "    bvap = graph.nodes[node]['BVAP']\n",
+    "    graph.nodes[node]['nBVAP'] = total_pop - bvap\n",
+    "    graph.nodes[node]['population'] = total_pop"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# from https://gerrychain.readthedocs.io/en/latest/user/recom.html\n",
+    "my_updaters = {\n",
+    "    \"population\": updaters.Tally(\"TOTPOP\", alias=\"population\"),\n",
+    "    \"cut_edges\": updaters.cut_edges,\n",
+    "}\n",
+    "elections = [\n",
+    "    Election(\"BVAP\", {\"BVAP\": \"BVAP\", \"nBVAP\": \"nBVAP\"})\n",
+    "]\n",
+    "election_updaters = {election.name: election for election in elections}\n",
+    "my_updaters.update(election_updaters)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import bvap_recom\n",
+    "from importlib import reload\n",
+    "bvap_recom = reload(bvap_recom)\n",
+    "\n",
+    "initial_partition = GeographicPartition(graph, assignment=\"HDIST_11\", updaters=my_updaters)\n",
+    "ideal_population = sum(initial_partition[\"population\"].values()) / len(initial_partition)\n",
+    "pop_tolerance = 0.15\n",
+    "\n",
+    "proposal = partial(bvap_recom.bvap_biased_recom,\n",
+    "                   pop_col=\"TOTPOP\",\n",
+    "                   bvap_col=\"BVAP\",\n",
+    "                   pop_target=ideal_population,\n",
+    "                   epsilon=pop_tolerance,\n",
+    "                   node_repeats=5,\n",
+    "                   target_districts=[86, 87],\n",
+    "                   p_target=0.25\n",
+    "                  )\n",
+    "\n",
+    "compactness_bound = constraints.UpperBound(\n",
+    "    lambda p: len(p[\"cut_edges\"]),\n",
+    "    2*len(initial_partition[\"cut_edges\"])\n",
+    ")\n",
+    "\n",
+    "pop_constraint = constraints.within_percent_of_ideal_population(initial_partition, pop_tolerance)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chain = MarkovChain(\n",
+    "    proposal=proposal,\n",
+    "    constraints=[\n",
+    "        pop_constraint,\n",
+    "        compactness_bound\n",
+    "    ],\n",
+    "    accept=accept.always_accept,\n",
+    "    initial_state=initial_partition,\n",
+    "    total_steps=1000\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bvap_86 = []\n",
+    "bvap_87 = []\n",
+    "for run in range(1):\n",
+    "    print('--- Run', run, '---')\n",
+    "    for idx, part in enumerate(chain):\n",
+    "        if idx % 25 == 0: print(idx)\n",
+    "        bvaps = sorted(part['BVAP'].percents('BVAP'))\n",
+    "        bvap_86.append(bvaps[86])\n",
+    "        bvap_87.append(bvaps[87])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.scatter(bvap_86, bvap_87)\n",
+    "plt.xlabel('BVAP 86')\n",
+    "plt.ylabel('BVAP 87')\n",
+    "plt.title(f'BVAP-biased ReCom ({len(bvap_86)} runs)')\n",
+    "plt.savefig('results/bvap_biased_recom/bvap_86_vs_87.png', dpi=300)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

bvap_recom.py

@@ -0,0 +1,59 @@
+from random import random, choice
+from collections import defaultdict
+from gerrychain.tree import recursive_tree_part
+
+def bvap_biased_recom(partition, pop_col, bvap_col, pop_target, epsilon,
+                      node_repeats, target_districts, p_target):
+    # TODO: Do we want to consider other denominators? (Voting-age population?)
+    bvaps = defaultdict(int)
+    pops = defaultdict(int)
+    for node in partition.graph.nodes:
+        bvaps[partition.assignment[node]] = partition.graph.nodes[node][bvap_col]
+        pops[partition.assignment[node]] = partition.graph.nodes[node][pop_col]
+    bvap_ratios = {
+        district_idx: bvaps[district_idx] / district_pop
+        for district_idx, district_pop in pops.items()
+    }
+    by_bvap = sorted(bvap_ratios.items(), key=lambda kv: kv[1])
+    bvap_mapping = {
+        bvap_idx: orig[0]
+        for bvap_idx, orig in enumerate(by_bvap)
+    }
+    edges = sorted(partition['cut_edges'], key=lambda k: k[0])
+    if random() < p_target:
+        # With probability p_target, we uniformly choose a target
+        # district to mix. Target indices are given by ordinal BVAP
+        # (index 0 refers to the district with the loewst BVAP).
+        target = choice(target_districts)
+        # Get the district index of the district with the ordinal
+        # BVAP value given by ``target``.
+        district_idx = bvap_mapping[target]
+        # print('Merging district', district_idx, '[', target, ']')
+        district_edges = [edge for edge in edges
+                          if partition.assignment[edge[0]] == district_idx]
+        edge = choice(district_edges)
+    else:
+        # Otherwise, we just run a normal ReCom step.
+        edge = choice(edges)
+    parts_to_merge = (partition.assignment[edge[0]],
+                      partition.assignment[edge[1]])
+
+    # (copied from original ReCom proposal)
+    subgraph = partition.graph.subgraph(
+        partition.parts[parts_to_merge[0]] | partition.parts[parts_to_merge[1]]
+    )
+
+    try:
+        flips = recursive_tree_part(
+            subgraph,
+            parts_to_merge,
+            pop_col=pop_col,
+            pop_target=pop_target,
+            epsilon=epsilon,
+            node_repeats=node_repeats,
+        )
+        return partition.flip(flips)
+    except KeyError as ex:
+        # TODO: Do something more robust here
+        print('ReCom fails!')
+        return partition

data/VA_precincts/VA_precincts.cpg

@@ -0,0 +1 @@
+ISO-8859-1

data/VA_precincts/VA_precincts.prj

@@ -0,0 +1 @@
+PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",37],PARAMETER["standard_parallel_2",39.5],PARAMETER["latitude_of_origin",36],PARAMETER["central_meridian",-79.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]