Bottleneck_Basic_Recom.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bottleneck distance visualisation of PA congressional plan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminary Steps\n",
    "\n",
    "Import some packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "\n",
    "import random\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial\n",
    "import networkx as nx\n",
    "\n",
    "from gerrychain import MarkovChain\n",
    "from gerrychain.constraints import (\n",
    "    Validator,\n",
    "    single_flip_contiguous,\n",
    "    within_percent_of_ideal_population,\n",
    ")\n",
    "from gerrychain.proposals import propose_random_flip\n",
    "from gerrychain.accept import always_accept\n",
    "from gerrychain.updaters import Election, Tally, cut_edges\n",
    "from gerrychain.partition import Partition\n",
    "from gerrychain.proposals import recom\n",
    "from gerrychain.metrics import mean_median, efficiency_gap\n",
    "\n",
    "import csv\n",
    "import os\n",
    "from functools import partial\n",
    "import json\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from gerrychain import (\n",
    "    Election,\n",
    "    Graph,\n",
    "    MarkovChain,\n",
    "    Partition,\n",
    "    accept,\n",
    "    constraints,\n",
    "    updaters,\n",
    ")\n",
    "\n",
    "from gerrychain.tree import recursive_tree_part"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define functions to create district adjacency graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def district_pairs(part):\n",
    "    pairs = []\n",
    "    for x, y in part['cut_edges']:\n",
    "        pairs.append((part.assignment[x],part.assignment[y]))\n",
    "    return set(pairs)\n",
    "\n",
    "def adjacency_graph_cut_edges(part):\n",
    "    edges = district_pairs(part)\n",
    "    adjacency_graph = nx.Graph()\n",
    "    adjacency_graph.add_nodes_from(list({x for (x,y) in edges}))\n",
    "    adjacency_graph.add_edges_from(list(edges))\n",
    "    return adjacency_graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import some data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop_col = \"TOT_POP\"\n",
    "election_names = [\n",
    "    \"PRES12\",\n",
    "    \"PRES16\",\n",
    "    \"SENW101216\",\n",
    "]\n",
    "election_columns = [\n",
    "    [\"PRES12D\", \"PRES12R\"],\n",
    "    [\"T16PRESD\", \"T16PRESR\"],\n",
    "    [\"W101216D\", \"W101216R\"],\n",
    "]\n",
    "graph = Graph.from_file(\"./PA_VTD/PA_VTD.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define updaters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "updaters1 = {\n",
    "    \"population\": updaters.Tally(\"TOT_POP\", alias=\"population\"),\n",
    "    \"cut_edges\": cut_edges,\n",
    "}\n",
    "\n",
    "elections = [\n",
    "    Election(\n",
    "        election_names[i],\n",
    "        {\"Democratic\": election_columns[i][0], \"Republican\": election_columns[i][1]},\n",
    "    )\n",
    "    for i in range(len(election_names))\n",
    "]\n",
    "\n",
    "election_updaters = {election.name: election for election in elections}\n",
    "\n",
    "updaters1.update(election_updaters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get some \"real\" plans"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will look at some 'premade' partitions. In particular, `partition_2011` is the real districting plan from 2011."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(\"Sample node: \", graph.nodes[0])\n",
    "#fix some strings which should be numbers\n",
    "for n in graph.nodes():\n",
    "    graph.nodes[n]['538CPCT__1'] = int(graph.nodes[n]['538CPCT__1'])\n",
    "    graph.nodes[n]['538DEM_PL'] = int(graph.nodes[n]['538DEM_PL'])\n",
    "    graph.nodes[n]['538GOP_PL'] = int(graph.nodes[n]['538GOP_PL'])\n",
    "    graph.nodes[n]['8THGRADE_1'] = int(graph.nodes[n]['8THGRADE_1'])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "partition_2011 = Partition(graph, \"2011_PLA_1\", updaters1)\n",
    "partition_GOV = Partition(graph, \"GOV\", updaters1)\n",
    "partition_TS = Partition(graph, \"TS\", updaters1)\n",
    "partition_REMEDIAL = Partition(graph, \"REMEDIAL_P\", updaters1)\n",
    "partition_CPCT = Partition(graph, \"538CPCT__1\", updaters1)\n",
    "partition_DEM = Partition(graph, \"538DEM_PL\", updaters1)\n",
    "partition_GOP = Partition(graph, \"538GOP_PL\", updaters1)\n",
    "partition_8th = Partition(graph, \"8THGRADE_1\", updaters1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ideal popuation:  704718.2777777778\n"
     ]
    }
   ],
   "source": [
    "ideal_population = sum(partition_2011[\"population\"].values())/len(partition_2011)\n",
    "print(\"Ideal popuation: \", ideal_population)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Recom Plans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_partition = Partition(graph, \"GOV\", updaters1)\n",
    "\n",
    "proposal = partial(\n",
    "    recom, pop_col=\"TOT_POP\", pop_target=ideal_population, epsilon=0.02, node_repeats=2\n",
    ")\n",
    "\n",
    "compactness_bound = constraints.UpperBound(\n",
    "    lambda p: len(p[\"cut_edges\"]), 2 * len(initial_partition[\"cut_edges\"])\n",
    ")\n",
    "\n",
    "recom_steps = 10000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "chain = MarkovChain(\n",
    "    proposal=proposal,\n",
    "    constraints=[\n",
    "        constraints.within_percent_of_ideal_population(initial_partition, 0.02),\n",
    "        compactness_bound,  # single_flip_contiguous#no_more_discontiguous\n",
    "    ],\n",
    "    accept=accept.always_accept,\n",
    "    initial_state=initial_partition,\n",
    "    total_steps=recom_steps,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Made Recom 20\n",
      "Made Recom 40\n",
      "Made Recom 60\n",
      "Made Recom 80\n",
      "Made Recom 100\n",
      "Made Recom 120\n",
      "Made Recom 140\n",
      "Made Recom 160\n",
      "Made Recom 180\n",
      "Made Recom 200\n",
      "Made Recom 220\n",
      "Made Recom 240\n",
      "Made Recom 260\n",
      "Made Recom 280\n",
      "Made Recom 300\n",
      "Made Recom 320\n",
      "Made Recom 340\n",
      "Made Recom 360\n",
      "Made Recom 380\n",
      "Made Recom 400\n",
      "Made Recom 420\n",
      "Made Recom 440\n",
      "Made Recom 460\n",
      "Made Recom 480\n",
      "Made Recom 500\n",
      "Made Recom 520\n",
      "Made Recom 540\n",
      "Made Recom 560\n",
      "Made Recom 580\n",
      "Made Recom 600\n",
      "Made Recom 620\n",
      "Made Recom 640\n",
      "Made Recom 660\n",
      "Made Recom 680\n",
      "Made Recom 700\n",
      "Made Recom 720\n",
      "Made Recom 740\n",
      "Made Recom 760\n",
      "Made Recom 780\n",
      "Made Recom 800\n",
      "Made Recom 820\n",
      "Made Recom 840\n",
      "Made Recom 860\n",
      "Made Recom 880\n",
      "Made Recom 900\n",
      "Made Recom 920\n",
      "Made Recom 940\n",
      "Made Recom 960\n",
      "Made Recom 980\n",
      "Made Recom 1000\n",
      "Made Recom 1020\n",
      "Made Recom 1040\n",
      "Made Recom 1060\n",
      "Made Recom 1080\n",
      "Made Recom 1100\n",
      "Made Recom 1120\n",
      "Made Recom 1140\n",
      "Made Recom 1160\n",
      "Made Recom 1180\n",
      "Made Recom 1200\n",
      "Made Recom 1220\n",
      "Made Recom 1240\n",
      "Made Recom 1260\n",
      "Made Recom 1280\n",
      "Made Recom 1300\n",
      "Made Recom 1320\n",
      "Made Recom 1340\n",
      "Made Recom 1360\n",
      "Made Recom 1380\n",
      "Made Recom 1400\n",
      "Made Recom 1420\n",
      "Made Recom 1440\n",
      "Made Recom 1460\n",
      "Made Recom 1480\n",
      "Made Recom 1500\n",
      "Made Recom 1520\n",
      "Made Recom 1540\n",
      "Made Recom 1560\n",
      "Made Recom 1580\n",
      "Made Recom 1600\n",
      "Made Recom 1620\n",
      "Made Recom 1640\n",
      "Made Recom 1660\n",
      "Made Recom 1680\n",
      "Made Recom 1700\n",
      "Made Recom 1720\n",
      "Made Recom 1740\n",
      "Made Recom 1760\n",
      "Made Recom 1780\n",
      "Made Recom 1800\n",
      "Made Recom 1820\n",
      "Made Recom 1840\n",
      "Made Recom 1860\n",
      "Made Recom 1880\n",
      "Made Recom 1900\n",
      "Made Recom 1920\n",
      "Made Recom 1940\n",
      "Made Recom 1960\n",
      "Made Recom 1980\n",
      "Made Recom 2000\n",
      "Made Recom 2020\n",
      "Made Recom 2040\n",
      "Made Recom 2060\n",
      "Made Recom 2080\n",
      "Made Recom 2100\n",
      "Made Recom 2120\n",
      "Made Recom 2140\n",
      "Made Recom 2160\n",
      "Made Recom 2180\n",
      "Made Recom 2200\n",
      "Made Recom 2220\n",
      "Made Recom 2240\n",
      "Made Recom 2260\n",
      "Made Recom 2280\n",
      "Made Recom 2300\n",
      "Made Recom 2320\n",
      "Made Recom 2340\n",
      "Made Recom 2360\n",
      "Made Recom 2380\n",
      "Made Recom 2400\n",
      "Made Recom 2420\n",
      "Made Recom 2440\n",
      "Made Recom 2460\n",
      "Made Recom 2480\n",
      "Made Recom 2500\n",
      "Made Recom 2520\n",
      "Made Recom 2540\n",
      "Made Recom 2560\n",
      "Made Recom 2580\n",
      "Made Recom 2600\n",
      "Made Recom 2620\n",
      "Made Recom 2640\n",
      "Made Recom 2660\n",
      "Made Recom 2680\n",
      "Made Recom 2700\n",
      "Made Recom 2720\n",
      "Made Recom 2740\n",
      "Made Recom 2760\n",
      "Made Recom 2780\n",
      "Made Recom 2800\n",
      "Made Recom 2820\n",
      "Made Recom 2840\n",
      "Made Recom 2860\n",
      "Made Recom 2880\n",
      "Made Recom 2900\n",
      "Made Recom 2920\n",
      "Made Recom 2940\n",
      "Made Recom 2960\n",
      "Made Recom 2980\n",
      "Made Recom 3000\n",
      "Made Recom 3020\n",
      "Made Recom 3040\n",
      "Made Recom 3060\n",
      "Made Recom 3080\n",
      "Made Recom 3100\n",
      "Made Recom 3120\n",
      "Made Recom 3140\n",
      "Made Recom 3160\n",
      "Made Recom 3180\n",
      "Made Recom 3200\n",
      "Made Recom 3220\n",
      "Made Recom 3240\n",
      "Made Recom 3260\n",
      "Made Recom 3280\n",
      "Made Recom 3300\n",
      "Made Recom 3320\n",
      "Made Recom 3340\n",
      "Made Recom 3360\n",
      "Made Recom 3380\n",
      "Made Recom 3400\n",
      "Made Recom 3420\n",
      "Made Recom 3440\n",
      "Made Recom 3460\n",
      "Made Recom 3480\n",
      "Made Recom 3500\n",
      "Made Recom 3520\n",
      "Made Recom 3540\n",
      "Made Recom 3560\n",
      "Made Recom 3580\n",
      "Made Recom 3600\n",
      "Made Recom 3620\n",
      "Made Recom 3640\n",
      "Made Recom 3660\n",
      "Made Recom 3680\n",
      "Made Recom 3700\n",
      "Made Recom 3720\n",
      "Made Recom 3740\n",
      "Made Recom 3760\n",
      "Made Recom 3780\n",
      "Made Recom 3800\n",
      "Made Recom 3820\n",
      "Made Recom 3840\n",
      "Made Recom 3860\n",
      "Made Recom 3880\n",
      "Made Recom 3900\n",
      "Made Recom 3920\n",
      "Made Recom 3940\n",
      "Made Recom 3960\n",
      "Made Recom 3980\n",
      "Made Recom 4000\n",
      "Made Recom 4020\n",
      "Made Recom 4040\n",
      "Made Recom 4060\n",
      "Made Recom 4080\n",
      "Made Recom 4100\n",
      "Made Recom 4120\n",
      "Made Recom 4140\n",
      "Made Recom 4160\n",
      "Made Recom 4180\n",
      "Made Recom 4200\n",
      "Made Recom 4220\n",
      "Made Recom 4240\n",
      "Made Recom 4260\n",
      "Made Recom 4280\n",
      "Made Recom 4300\n",
      "Made Recom 4320\n",
      "Made Recom 4340\n",
      "Made Recom 4360\n",
      "Made Recom 4380\n",
      "Made Recom 4400\n",
      "Made Recom 4420\n",
      "Made Recom 4440\n",
      "Made Recom 4460\n",
      "Made Recom 4480\n",
      "Made Recom 4500\n",
      "Made Recom 4520\n",
      "Made Recom 4540\n",
      "Made Recom 4560\n",
      "Made Recom 4580\n",
      "Made Recom 4600\n",
      "Made Recom 4620\n",
      "Made Recom 4640\n",
      "Made Recom 4660\n",
      "Made Recom 4680\n",
      "Made Recom 4700\n",
      "Made Recom 4720\n",
      "Made Recom 4740\n",
      "Made Recom 4760\n",
      "Made Recom 4780\n",
      "Made Recom 4800\n",
      "Made Recom 4820\n",
      "Made Recom 4840\n",
      "Made Recom 4860\n",
      "Made Recom 4880\n",
      "Made Recom 4900\n",
      "Made Recom 4920\n",
      "Made Recom 4940\n",
      "Made Recom 4960\n",
      "Made Recom 4980\n",
      "Made Recom 5000\n",
      "Made Recom 5020\n",
      "Made Recom 5040\n",
      "Made Recom 5060\n",
      "Made Recom 5080\n",
      "Made Recom 5100\n",
      "Made Recom 5120\n",
      "Made Recom 5140\n",
      "Made Recom 5160\n",
      "Made Recom 5180\n",
      "Made Recom 5200\n",
      "Made Recom 5220\n",
      "Made Recom 5240\n",
      "Made Recom 5260\n",
      "Made Recom 5280\n",
      "Made Recom 5300\n",
      "Made Recom 5320\n",
      "Made Recom 5340\n",
      "Made Recom 5360\n",
      "Made Recom 5380\n",
      "Made Recom 5400\n",
      "Made Recom 5420\n",
      "Made Recom 5440\n",
      "Made Recom 5460\n",
      "Made Recom 5480\n",
      "Made Recom 5500\n",
      "Made Recom 5520\n",
      "Made Recom 5540\n",
      "Made Recom 5560\n",
      "Made Recom 5580\n",
      "Made Recom 5600\n",
      "Made Recom 5620\n",
      "Made Recom 5640\n",
      "Made Recom 5660\n",
      "Made Recom 5680\n",
      "Made Recom 5700\n",
      "Made Recom 5720\n",
      "Made Recom 5740\n",
      "Made Recom 5760\n",
      "Made Recom 5780\n",
      "Made Recom 5800\n",
      "Made Recom 5820\n",
      "Made Recom 5840\n",
      "Made Recom 5860\n",
      "Made Recom 5880\n",
      "Made Recom 5900\n",
      "Made Recom 5920\n",
      "Made Recom 5940\n",
      "Made Recom 5960\n",
      "Made Recom 5980\n",
      "Made Recom 6000\n",
      "Made Recom 6020\n",
      "Made Recom 6040\n",
      "Made Recom 6060\n",
      "Made Recom 6080\n",
      "Made Recom 6100\n",
      "Made Recom 6120\n",
      "Made Recom 6140\n",
      "Made Recom 6160\n",
      "Made Recom 6180\n",
      "Made Recom 6200\n",
      "Made Recom 6220\n",
      "Made Recom 6240\n",
      "Made Recom 6260\n",
      "Made Recom 6280\n",
      "Made Recom 6300\n",
      "Made Recom 6320\n",
      "Made Recom 6340\n",
      "Made Recom 6360\n",
      "Made Recom 6380\n",
      "Made Recom 6400\n",
      "Made Recom 6420\n",
      "Made Recom 6440\n",
      "Made Recom 6460\n",
      "Made Recom 6480\n",
      "Made Recom 6500\n",
      "Made Recom 6520\n",
      "Made Recom 6540\n",
      "Made Recom 6560\n",
      "Made Recom 6580\n",
      "Made Recom 6600\n",
      "Made Recom 6620\n",
      "Made Recom 6640\n",
      "Made Recom 6660\n",
      "Made Recom 6680\n",
      "Made Recom 6700\n",
      "Made Recom 6720\n",
      "Made Recom 6740\n",
      "Made Recom 6760\n",
      "Made Recom 6780\n",
      "Made Recom 6800\n",
      "Made Recom 6820\n",
      "Made Recom 6840\n",
      "Made Recom 6860\n",
      "Made Recom 6880\n",
      "Made Recom 6900\n",
      "Made Recom 6920\n",
      "Made Recom 6940\n",
      "Made Recom 6960\n",
      "Made Recom 6980\n",
      "Made Recom 7000\n",
      "Made Recom 7020\n",
      "Made Recom 7040\n",
      "Made Recom 7060\n",
      "Made Recom 7080\n",
      "Made Recom 7100\n",
      "Made Recom 7120\n",
      "Made Recom 7140\n",
      "Made Recom 7160\n",
      "Made Recom 7180\n",
      "Made Recom 7200\n",
      "Made Recom 7220\n",
      "Made Recom 7240\n",
      "Made Recom 7260\n",
      "Made Recom 7280\n",
      "Made Recom 7300\n",
      "Made Recom 7320\n",
      "Made Recom 7340\n",
      "Made Recom 7360\n",
      "Made Recom 7380\n",
      "Made Recom 7400\n",
      "Made Recom 7420\n",
      "Made Recom 7440\n",
      "Made Recom 7460\n",
      "Made Recom 7480\n",
      "Made Recom 7500\n",
      "Made Recom 7520\n",
      "Made Recom 7540\n",
      "Made Recom 7560\n",
      "Made Recom 7580\n",
      "Made Recom 7600\n",
      "Made Recom 7620\n",
      "Made Recom 7640\n",
      "Made Recom 7660\n",
      "Made Recom 7680\n",
      "Made Recom 7700\n",
      "Made Recom 7720\n",
      "Made Recom 7740\n",
      "Made Recom 7760\n",
      "Made Recom 7780\n",
      "Made Recom 7800\n",
      "Made Recom 7820\n",
      "Made Recom 7840\n",
      "Made Recom 7860\n",
      "Made Recom 7880\n",
      "Made Recom 7900\n",
      "Made Recom 7920\n",
      "Made Recom 7940\n",
      "Made Recom 7960\n",
      "Made Recom 7980\n",
      "Made Recom 8000\n",
      "Made Recom 8020\n",
      "Made Recom 8040\n",
      "Made Recom 8060\n",
      "Made Recom 8080\n",
      "Made Recom 8100\n",
      "Made Recom 8120\n",
      "Made Recom 8140\n",
      "Made Recom 8160\n",
      "Made Recom 8180\n",
      "Made Recom 8200\n",
      "Made Recom 8220\n",
      "Made Recom 8240\n",
      "Made Recom 8260\n",
      "Made Recom 8280\n",
      "Made Recom 8300\n",
      "Made Recom 8320\n",
      "Made Recom 8340\n",
      "Made Recom 8360\n",
      "Made Recom 8380\n",
      "Made Recom 8400\n",
      "Made Recom 8420\n",
      "Made Recom 8440\n",
      "Made Recom 8460\n",
      "Made Recom 8480\n",
      "Made Recom 8500\n",
      "Made Recom 8520\n",
      "Made Recom 8540\n",
      "Made Recom 8560\n",
      "Made Recom 8580\n",
      "Made Recom 8600\n",
      "Made Recom 8620\n",
      "Made Recom 8640\n",
      "Made Recom 8660\n",
      "Made Recom 8680\n",
      "Made Recom 8700\n",
      "Made Recom 8720\n",
      "Made Recom 8740\n",
      "Made Recom 8760\n",
      "Made Recom 8780\n",
      "Made Recom 8800\n",
      "Made Recom 8820\n",
      "Made Recom 8840\n",
      "Made Recom 8860\n",
      "Made Recom 8880\n",
      "Made Recom 8900\n",
      "Made Recom 8920\n",
      "Made Recom 8940\n",
      "Made Recom 8960\n",
      "Made Recom 8980\n",
      "Made Recom 9000\n",
      "Made Recom 9020\n",
      "Made Recom 9040\n",
      "Made Recom 9060\n",
      "Made Recom 9080\n",
      "Made Recom 9100\n",
      "Made Recom 9120\n",
      "Made Recom 9140\n",
      "Made Recom 9160\n",
      "Made Recom 9180\n",
      "Made Recom 9200\n",
      "Made Recom 9220\n",
      "Made Recom 9240\n",
      "Made Recom 9260\n",
      "Made Recom 9280\n",
      "Made Recom 9300\n",
      "Made Recom 9320\n",
      "Made Recom 9340\n",
      "Made Recom 9360\n",
      "Made Recom 9380\n",
      "Made Recom 9400\n",
      "Made Recom 9420\n",
      "Made Recom 9440\n",
      "Made Recom 9460\n",
      "Made Recom 9480\n",
      "Made Recom 9500\n",
      "Made Recom 9520\n",
      "Made Recom 9540\n",
      "Made Recom 9560\n",
      "Made Recom 9580\n",
      "Made Recom 9600\n",
      "Made Recom 9620\n",
      "Made Recom 9640\n",
      "Made Recom 9660\n",
      "Made Recom 9680\n",
      "Made Recom 9700\n",
      "Made Recom 9720\n",
      "Made Recom 9740\n",
      "Made Recom 9760\n",
      "Made Recom 9780\n",
      "Made Recom 9800\n",
      "Made Recom 9820\n",
      "Made Recom 9840\n",
      "Made Recom 9860\n",
      "Made Recom 9880\n",
      "Made Recom 9900\n",
      "Made Recom 9920\n",
      "Made Recom 9940\n",
      "Made Recom 9960\n",
      "Made Recom 9980\n",
      "Made Recom 10000\n"
     ]
    }
   ],
   "source": [
    "recom_partitions = []\n",
    "recom_labels = []\n",
    "t=0\n",
    "for part in chain:\n",
    "    t+=1\n",
    "    #recom_partitions.append(part)\n",
    "    rand = random.random()\n",
    "    if rand < 0.01:\n",
    "        recom_partitions.append(part)\n",
    "        recom_labels.append(\"GOVRecom\" + str(t))\n",
    "    if t % 20 == 0:\n",
    "        print(\"Made Recom \" + str(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "second_partition = Partition(graph, \"2011_PLA_1\", updaters1)\n",
    "\n",
    "second_chain = MarkovChain(\n",
    "    proposal=proposal,\n",
    "    constraints=[\n",
    "        constraints.within_percent_of_ideal_population(second_partition, 0.02),\n",
    "        compactness_bound,  # single_flip_contiguous#no_more_discontiguous\n",
    "    ],\n",
    "    accept=accept.always_accept,\n",
    "    initial_state=second_partition,\n",
    "    total_steps=recom_steps,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Made Recom 20\n",
      "Made Recom 40\n",
      "Made Recom 60\n",
      "Made Recom 80\n",
      "Made Recom 100\n",
      "Made Recom 120\n",
      "Made Recom 140\n",
      "Made Recom 160\n",
      "Made Recom 180\n",
      "Made Recom 200\n",
      "Made Recom 220\n",
      "Made Recom 240\n",
      "Made Recom 260\n",
      "Made Recom 280\n",
      "Made Recom 300\n",
      "Made Recom 320\n",
      "Made Recom 340\n",
      "Made Recom 360\n",
      "Made Recom 380\n",
      "Made Recom 400\n",
      "Made Recom 420\n",
      "Made Recom 440\n",
      "Made Recom 460\n",
      "Made Recom 480\n",
      "Made Recom 500\n",
      "Made Recom 520\n",
      "Made Recom 540\n",
      "Made Recom 560\n",
      "Made Recom 580\n",
      "Made Recom 600\n",
      "Made Recom 620\n",
      "Made Recom 640\n",
      "Made Recom 660\n",
      "Made Recom 680\n",
      "Made Recom 700\n",
      "Made Recom 720\n",
      "Made Recom 740\n",
      "Made Recom 760\n",
      "Made Recom 780\n",
      "Made Recom 800\n",
      "Made Recom 820\n",
      "Made Recom 840\n",
      "Made Recom 860\n",
      "Made Recom 880\n",
      "Made Recom 900\n",
      "Made Recom 920\n",
      "Made Recom 940\n",
      "Made Recom 960\n",
      "Made Recom 980\n",
      "Made Recom 1000\n",
      "Made Recom 1020\n",
      "Made Recom 1040\n",
      "Made Recom 1060\n",
      "Made Recom 1080\n",
      "Made Recom 1100\n",
      "Made Recom 1120\n",
      "Made Recom 1140\n",
      "Made Recom 1160\n",
      "Made Recom 1180\n",
      "Made Recom 1200\n",
      "Made Recom 1220\n",
      "Made Recom 1240\n",
      "Made Recom 1260\n",
      "Made Recom 1280\n",
      "Made Recom 1300\n",
      "Made Recom 1320\n",
      "Made Recom 1340\n",
      "Made Recom 1360\n",
      "Made Recom 1380\n",
      "Made Recom 1400\n",
      "Made Recom 1420\n",
      "Made Recom 1440\n",
      "Made Recom 1460\n",
      "Made Recom 1480\n",
      "Made Recom 1500\n",
      "Made Recom 1520\n",
      "Made Recom 1540\n",
      "Made Recom 1560\n",
      "Made Recom 1580\n",
      "Made Recom 1600\n",
      "Made Recom 1620\n",
      "Made Recom 1640\n",
      "Made Recom 1660\n",
      "Made Recom 1680\n",
      "Made Recom 1700\n",
      "Made Recom 1720\n",
      "Made Recom 1740\n",
      "Made Recom 1760\n",
      "Made Recom 1780\n",
      "Made Recom 1800\n",
      "Made Recom 1820\n",
      "Made Recom 1840\n",
      "Made Recom 1860\n",
      "Made Recom 1880\n",
      "Made Recom 1900\n",
      "Made Recom 1920\n",
      "Made Recom 1940\n",
      "Made Recom 1960\n",
      "Made Recom 1980\n",
      "Made Recom 2000\n",
      "Made Recom 2020\n",
      "Made Recom 2040\n",
      "Made Recom 2060\n",
      "Made Recom 2080\n",
      "Made Recom 2100\n",
      "Made Recom 2120\n",
      "Made Recom 2140\n",
      "Made Recom 2160\n",
      "Made Recom 2180\n",
      "Made Recom 2200\n",
      "Made Recom 2220\n",
      "Made Recom 2240\n",
      "Made Recom 2260\n",
      "Made Recom 2280\n",
      "Made Recom 2300\n",
      "Made Recom 2320\n",
      "Made Recom 2340\n",
      "Made Recom 2360\n",
      "Made Recom 2380\n",
      "Made Recom 2400\n",
      "Made Recom 2420\n",
      "Made Recom 2440\n",
      "Made Recom 2460\n",
      "Made Recom 2480\n",
      "Made Recom 2500\n",
      "Made Recom 2520\n",
      "Made Recom 2540\n",
      "Made Recom 2560\n",
      "Made Recom 2580\n",
      "Made Recom 2600\n",
      "Made Recom 2620\n",
      "Made Recom 2640\n",
      "Made Recom 2660\n",
      "Made Recom 2680\n",
      "Made Recom 2700\n",
      "Made Recom 2720\n"
     ]
    }
   ],
   "source": [
    "second_recom_partitions = []\n",
    "second_recom_labels = []\n",
    "t=0\n",
    "for part in chain:\n",
    "    t+=1\n",
    "    #recom_partitions.append(part)\n",
    "    rand = random.random()\n",
    "    if rand < 0.01:\n",
    "        second_recom_partitions.append(part)\n",
    "        second_recom_labels.append(\"2011Recom\" + str(t))\n",
    "    if t % 20 == 0:\n",
    "        print(\"Made Recom \" + str(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "third_partition = Partition(graph, \"538DEM_PL\", updaters1)\n",
    "\n",
    "third_chain = MarkovChain(\n",
    "    proposal=proposal,\n",
    "    constraints=[\n",
    "        constraints.within_percent_of_ideal_population(third_partition, 0.02),\n",
    "        compactness_bound,  # single_flip_contiguous#no_more_discontiguous\n",
    "    ],\n",
    "    accept=accept.always_accept,\n",
    "    initial_state=third_partition,\n",
    "    total_steps=recom_steps,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "third_recom_partitions = []\n",
    "third_recom_labels = []\n",
    "t=0\n",
    "for part in chain:\n",
    "    t+=1\n",
    "    #recom_partitions.append(part)\n",
    "    rand = random.random()\n",
    "    if rand < 0.01:\n",
    "        third_recom_partitions.append(part)\n",
    "        third_recom_labels.append(\"538DemRecom\" + str(t))\n",
    "    if t % 20 == 0:\n",
    "        print(\"Made Recom \" + str(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create tree plans and add them to the list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_partitions = []\n",
    "\n",
    "tree_plans = 20 #set this higher later...\n",
    "\n",
    "for i in range(tree_plans):\n",
    "    print('Finished tree plan', i)\n",
    "    cddict = recursive_tree_part(graph, range(18), ideal_population, \"TOT_POP\", .01, 2)\n",
    "    tree_partitions.append(Partition(graph, cddict, updaters1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "partition_list = [partition_2011, partition_GOV, partition_TS,\n",
    "                  partition_REMEDIAL, partition_CPCT, partition_DEM,\n",
    "                  partition_GOP, partition_8th]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = ['2011', 'GOV', 'TS', 'REMEDIAL', 'CPCT', 'DEM', 'GOP', '8th']+[\"Tree\" + str(k + 1) for k in range(tree_plans)];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels += recom_labels\n",
    "labels += second_recom_labels\n",
    "labels += third_recom_labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(labels)\n",
    "print(len(labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_partition_list = partition_list + tree_partitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_partition_list += recom_partitions\n",
    "full_partition_list += second_recom_partitions\n",
    "full_partition_list += third_recom_partitions\n",
    "#full_partition_list = [partition_GOV] + recom_partitions\n",
    "#labels = ['GOV'] + [\"Recom\" + str(k+1) for k in range(recom_steps)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_partitions = len(full_partition_list)\n",
    "print(num_partitions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing Adjacency Graphs and Barcodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gudhi as gd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for k in range(num_partitions):\n",
    "    part = full_partition_list[k]\n",
    "\n",
    "    adjacency_graph = adjacency_graph_cut_edges(part) #2 districts have edge between them if share boundary\n",
    "\n",
    "    spCpx = gd.SimplexTree()\n",
    "    for edge in adjacency_graph.edges:\n",
    "        spCpx.insert(list(edge))\n",
    "\n",
    "    Democratic_voter_share = part['PRES12'].percents('Democratic')\n",
    "\n",
    "    zero_skeleton = spCpx.get_skeleton(0) #all the vertices\n",
    "\n",
    "    for j in range(len(zero_skeleton)):\n",
    "        spCpx.assign_filtration(zero_skeleton[j][0], filtration=Democratic_voter_share[j])\n",
    "\n",
    "    spCpx.make_filtration_non_decreasing() # when have boundary of complex, add it?\n",
    "\n",
    "    BarCodes = spCpx.persistence()\n",
    "\n",
    "    fig = plt.figure(figsize = (10,5))\n",
    "\n",
    "    ax = fig.add_subplot(1,2,1)\n",
    "    nx.draw_networkx(adjacency_graph, node_color=Democratic_voter_share, cmap=plt.cm.Blues)\n",
    "    plt.title('Partition '+labels[k])\n",
    "    ax.axis('off')\n",
    "\n",
    "    fig.add_subplot(1,2,2)\n",
    "    gd.plot_persistence_diagram(BarCodes);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bottleneck Distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "partitions_spCpx = []\n",
    "\n",
    "for k in range(num_partitions):\n",
    "    part = full_partition_list[k]\n",
    "\n",
    "    adjacency_graph = adjacency_graph_cut_edges(part)\n",
    "\n",
    "    spCpx = gd.SimplexTree()\n",
    "    for edge in adjacency_graph.edges:\n",
    "        spCpx.insert(list(edge))\n",
    "    \n",
    "    Democratic_voter_share = part['PRES12'].percents('Democratic')\n",
    "    \n",
    "    zero_skeleton = spCpx.get_skeleton(0)\n",
    "\n",
    "    for j in range(len(zero_skeleton)):\n",
    "        spCpx.assign_filtration(zero_skeleton[j][0], filtration=Democratic_voter_share[j])\n",
    "\n",
    "    spCpx.make_filtration_non_decreasing()\n",
    "\n",
    "    partitions_spCpx.append(spCpx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create a distance matrix containing the bottleneck distance between any pair of barcodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distMat = np.zeros((num_partitions,num_partitions))\n",
    "\n",
    "for j in range(num_partitions):\n",
    "    for k in range(num_partitions):\n",
    "        spCpx0 = partitions_spCpx[j]\n",
    "        spCpx1 = partitions_spCpx[k]\n",
    "        \n",
    "        spCpx0.persistence()\n",
    "        spCpx1.persistence()\n",
    "\n",
    "        I0 = spCpx0.persistence_intervals_in_dimension(0)\n",
    "        I1 = spCpx1.persistence_intervals_in_dimension(0)\n",
    "\n",
    "        distMat[j,k] = gd.bottleneck_distance(I0,I1)\n",
    "        \n",
    "# Symmetrize to remove any numerical errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distMat = np.maximum(distMat, distMat.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We view the distance matrix as an image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(distMat,cmap = 'hot');\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand the shape of the 'space of districting plans', we can use Multi-Dimensional Scaling. This algorithm looks for the set of points in $\\mathbb{R}^2$ (or $\\mathbb{R}^3$) whose distance matrix is as close as possible to the distance matrix we just computed. For a more precise description, go here: https://en.wikipedia.org/wiki/Multidimensional_scaling\n",
    "\n",
    "The result gives us a visualization of how similar the districting plans are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = len(partition_list)*['r'] + len(tree_partitions)*['b']+len(recom_partitions)*['g'] + len(second_recom_partitions)*['m'] + len(third_recom_partitions)*['c']\n",
    "#colors = ['r']+len(recom_partitions)*['g']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Import a package containing the MDS algorithm and set options for the algorithm\n",
    "from sklearn import manifold\n",
    "mds = manifold.MDS(n_components=2, dissimilarity=\"precomputed\")\n",
    "\n",
    "# Compute MDS and extract the coordinates of the points\n",
    "results = mds.fit(distMat)\n",
    "coords = results.embedding_\n",
    "\n",
    "z = coords[:,0]\n",
    "y = coords[:,1]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(z, y, c = colors)\n",
    "\n",
    "for i, txt in enumerate(labels):\n",
    "    ax.annotate(txt, (z[i], y[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import a package containing the MDS algorithm and set options for the algorithm\n",
    "from sklearn import manifold\n",
    "mds = manifold.MDS(n_components=2, dissimilarity=\"precomputed\")\n",
    "\n",
    "# Compute MDS and extract the coordinates of the points\n",
    "results = mds.fit(distMat)\n",
    "coords = results.embedding_\n",
    "\n",
    "z = coords[:,0]\n",
    "y = coords[:,1]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(z, y, c = colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use MDS to project the data to 3D. It takes a little more work to make the plot. You can try running this a few times; MDS doesn't have a unique answer and the algorithm involves some randomness, so you will get something different every time you run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "mds = manifold.MDS(n_components=3, dissimilarity=\"precomputed\")\n",
    "results = mds.fit(distMat)\n",
    "coords = results.embedding_\n",
    "\n",
    "fig = plt.figure(figsize = (7,7))\n",
    "ax = Axes3D(fig)\n",
    "\n",
    "for i, txt in enumerate(labels):\n",
    "    ax.scatter(coords[i,0],coords[i,1],coords[i,2], c = colors[i]) \n",
    "    ax.text(coords[i,0],coords[i,1],coords[i,2],  '%s' % (txt), size=10, zorder=1) \n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "mds = manifold.MDS(n_components=3, dissimilarity=\"precomputed\")\n",
    "results = mds.fit(distMat)\n",
    "coords = results.embedding_\n",
    "\n",
    "fig = plt.figure(figsize = (7,7))\n",
    "ax = Axes3D(fig)\n",
    "\n",
    "for i, txt in enumerate(labels):\n",
    "    ax.scatter(coords[i,0],coords[i,1],coords[i,2], c = colors[i]) \n",
    "\n",
    "\n",
    "plt.show()\n",
    "plt.savefig('PennsylvaniaRecomExperiment.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.cluster import hierarchy\n",
    "\n",
    "linkage = hierarchy.linkage(distMat, 'single')\n",
    "\n",
    "\n",
    "plt.figure(figsize = (10,10))\n",
    "hierarchy.dendrogram(linkage, labels = colors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}