{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Imports\n", "\n", "from gerrychain import Graph, GeographicPartition, Partition, Election\n", "from gerrychain.updaters import Tally, cut_edges\n", "import geopandas as gpd\n", "import numpy as np\n", "from gerrychain.random import random\n", "\n", "from gerrychain import MarkovChain\n", "from gerrychain.constraints import single_flip_contiguous, contiguous\n", "from gerrychain.proposals import propose_random_flip, propose_chunk_flip\n", "from gerrychain.accept import always_accept\n", "from gerrychain.metrics import polsby_popper\n", "from gerrychain import constraints\n", "\n", "\n", "from gerrychain.metrics import partisan_gini\n", "\n", "\n", "\n", "\n", "from gerrychain.tree import recursive_tree_part\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# read in data\n", "\n", "# set to \"True\" for North Carolina, \"False for Pennsylvania\"\n", "nc = True\n", "\n", "if nc:\n", " df = gpd.read_file(\"https://github.com/mggg-states/NC-shapefiles/raw/master/NC_VTD.zip\")\n", " num_dist = 13\n", " sen_elec = Election(\"SENElec\", {\"Democratic\": \"EL14G_US_1\", \"Republican\": \"EL14G_USS_\"})\n", " \n", " def_assn = \"newplan\"\n", "\n", "else:\n", " df = gpd.read_file(\"https://github.com/mggg-states/PA-shapefiles/raw/master/PA_VTDs.zip\")\n", " num_dist = 18\n", " sen_elec = Election(\"SENElec\", {\"Republican\": \"T16SENR\", \"Democratic\": \"T16SEND\"})\n", " \n", " def_assn = \"TS\"\n", " \n", " \n", " \n", "graph = Graph.from_geodataframe(df, reproject=False,ignore_errors=True)\n", "graph.add_data(df,list(df)) \n", "\n", "\n", "pop = sum(df[\"TOTPOP\"])\n", "\n", "\n", "updaters = {\n", " \"population\": Tally(\"TOTPOP\", alias=\"population\"),\n", " \"cut_edges\": cut_edges,\n", " \"polsby_popper\": polsby_popper,\n", " sen_elec.name : sen_elec\n", "}\n", "initial_partition = GeographicPartition(graph, assignment=def_assn, updaters = updaters)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# function definitions\n", "\n", "\n", "def no_worse_pg(partition):\n", " if random.random() < .001: return True\n", " return abs(partisan_gini(partition[\"SENElec\"])) <= abs(partisan_gini(partition.parent[\"SENElec\"]))\n", "\n", "\n", "\n", "def no_worse_pp(partition):\n", " if random.random() < .001: return True\n", " return np.mean(list(partition[\"polsby_popper\"].values())) >= np.mean(list(partition.parent[\"polsby_popper\"].values())) \n", "\n", "\n", "\n", "def step2temp(step,maxsteps):\n", " \n", " if step/maxsteps < .33:\n", " return 0\n", " else:\n", " return 4* (((step/maxsteps) - .33)/(1-.33))\n", "\n", "\n", "\n", "class MH_Annealer:\n", " def __init__(self, temp, maxsteps):\n", " self.counter = 0 # counts number of steps in the chain, increment on accept\n", " self.temp = temp # function which eats a count and gives the current temp\n", " self.maxsteps = maxsteps\n", " \n", " \n", " \n", " def __call__(self, partition):\n", " if self.temp(self.counter, self.maxsteps) == 0 : \n", " self.counter +=1 \n", " return True\n", " else:\n", " prob = get_score(partition) / (get_score(partition.parent) * self.temp(self.counter, self.maxsteps))\n", " \n", " \n", " if prob > random.random():\n", " self.counter += 1\n", " return True\n", " else:\n", " return False\n", "\n", " \n", "def sometimes_chunk(partition):\n", " r = random.random()\n", " if r < .65: return propose_random_flip(partition)\n", " else: return propose_chunk_flip(partition)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# run the chain\n", "\n", "pps = []\n", "pgs = []\n", "\n", "samples = []\n", "\n", "\n", "# increase these for larger sample, will take much longer to run!\n", "steps = 3\n", "#chainsteps == 150000 in the paper, can be adjusted if desired\n", "chainsteps = 20000\n", "\n", "\n", "for i in range(steps):\n", " \n", " samp = []\n", " \n", " accept = MH_Annealer(temp = step2temp, maxsteps = chainsteps)\n", " new_plan = recursive_tree_part(graph, range(num_dist), pop/num_dist, \"TOTPOP\", .02, 1)\n", " new_plan = GeographicPartition(graph, assignment=new_plan, updaters = updaters)\n", " \n", " \n", " def pareto_improves(partition):\n", " \n", " if count > 9.5*chainsteps/10: return (no_worse_pg(partition) and no_worse_pp(partition))\n", " return no_worse_pp(partition)\n", " \n", " t = 1- (i+1)/(steps)\n", " \n", " if random.random()**16 <= t:\n", " return no_worse_pp(partition)\n", " else:\n", " return no_worse_pg(partition)\n", " \n", " \n", " \n", " \n", " \n", " chain = MarkovChain(\n", " proposal=sometimes_chunk,\n", " constraints=[\n", " constraints.within_percent_of_ideal_population(initial_partition, 0.05),\n", " contiguous\n", " ],\n", " accept=pareto_improves,\n", " initial_state=new_plan,\n", " total_steps=chainsteps,\n", " )\n", "\n", " samp.append(new_plan)\n", " count = 0\n", " for p in chain:\n", " count+=1\n", " if count %10000 ==0:\n", " samp.append(p)\n", " pgs.append(partisan_gini(p[\"SENElec\"]))\n", " pps.append(np.mean(list(p[\"polsby_popper\"].values())))\n", " samples.append(samp)" ] }, { "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.0" } }, "nbformat": 4, "nbformat_minor": 2 }