<jupyter_start><jupyter_text>Short BurstsWe look at the behavior of using short bursts to increase number of opportunity districts for bvap and it projection into the line model<jupyter_code>import matplotlib.pyplot as plt
from gerrychain import (GeographicPartition, Partition, Graph, MarkovChain,
                        proposals, updaters, constraints, accept, Election)
from gerrychain.proposals import recom, propose_random_flip
from gerrychain.tree import recursive_tree_part
from gerrychain.metrics import mean_median, efficiency_gap, polsby_popper, partisan_gini
from functools import (partial, reduce)
import pandas
import geopandas as gp
import numpy as np
import networkx as nx
import pickle
import seaborn as sns
import pprint
import operator
import scipy
from mpl_toolkits.mplot3d import Axes3D
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale, normalize
import random
from nltk.util import bigrams
from nltk.probability import FreqDist
import pylab
## This function takes a name of a shapefile and returns a tuple of the graph
## and its associated dataframe
def build_graph(filename):
    print("Pulling in Graph from Shapefile: " + filename)
    graph = Graph.from_file(filename)
    df = gp.read_file(filename)
    return(graph, df)
def config_markov_chain(initial_part, iters=1000, epsilon=0.05, compactness=True, pop="TOT_POP"):
    # The recom proposal needs to know the ideal population for the districts so
    # that we can improve speed by bailing early on unbalanced partitions.
    ideal_population = sum(initial_part["population"].values()) / len(initial_part)


    # We use functools.partial to bind the extra parameters (pop_col, pop_target,
    # epsilon, node_repeats) of the recom proposal.
    proposal = partial(recom,
                       pop_col=pop,
                       pop_target=ideal_population,
                       epsilon=epsilon,
                       node_repeats=1)

    # To keep districts about as compact as the original plan, we bound the number
    # of cut edges at 2 times the number of cut edges in the initial plan.
    if compactness:
        compactness_bound = constraints.UpperBound(lambda p: len(p["cut_edges"]),
                            2*len(initial_part["cut_edges"]))
        cs = [constraints.within_percent_of_ideal_population(initial_part, epsilon), compactness_bound]
    else:
        cs = [constraints.within_percent_of_ideal_population(initial_part, epsilon)]


    # Configure the MarkovChain.
    return MarkovChain(proposal=proposal, constraints=cs, accept=accept.always_accept,
                       initial_state=initial_part, total_steps=iters)<jupyter_output><empty_output><jupyter_text>Functions to Simulate Short Bursts<jupyter_code>## This function takes a partition and returns the number o
def num_opportunity_dists(part, minority_prec="bvap_prec"):
    dist_precs = part[minority_prec].values()
    return sum(list(map(lambda v: v >= 0.35, dist_precs)))
def short_burst_run(starting_part, score_func, num_bursts=10, num_steps=100, verbose=False,
                    pop="TOT_POP"):
    max_part = (starting_part, num_opportunity_dists(starting_part))
    observed_num_ops = np.zeros((num_bursts, num_steps))

    for i in range(num_bursts):
        if verbose: print("Burst:", i)
        chain = config_markov_chain(max_part[0], iters=num_steps, epsilon=0.1, pop=pop)

        for j, part in enumerate(chain):
            part_score = num_opportunity_dists(part)
            observed_num_ops[i][j] = part_score
            max_part = (part, part_score) if part_score >= max_part[1] else max_part
    
    return (max_part, observed_num_ops)<jupyter_output><empty_output><jupyter_text>Pull in PA data and set up partitions<jupyter_code># graph_PA, df_PA = build_graph("PA_shape/PA_VTD.shp")
graph_PA = pickle.load(open("PA_graph.p", "rb"))
df_PA = pickle.load(open("PA_df.p", "rb"))
PA_updaters = {"population": updaters.Tally("TOT_POP", alias="population"),
               "bvap": updaters.Tally("BLACK_POP", alias="bvap"),
               "vap": updaters.Tally("VAP", alias="vap"),
               "bvap_prec": lambda part: {k: part["bvap"][k] / part["population"][k] for k in part["bvap"]}}
total_population_PA = sum(df_PA.TOT_POP.values)
ideal_population_PA = total_population_PA / 50
seed_part_senate = recursive_tree_part(graph_PA, range(50), pop_col="TOT_POP",
                                pop_target=ideal_population_PA,
                                epsilon=0.01, node_repeats=1)
PA_seed_seante = GeographicPartition(graph_PA, assignment=seed_part_senate,updaters=PA_updaters)
seed_part_house = recursive_tree_part(graph_PA, range(203), pop_col="TOT_POP",
                                pop_target=ideal_population_PA,
                                epsilon=0.01, node_repeats=1)
PA_seed_house = GeographicPartition(graph_PA, assignment=seed_part_house,updaters=PA_updaters)
enacted_senate = GeographicPartition(graph_PA, assignment="SSD", updaters=PA_updaters)
num_opportunity_dists(enacted_senate)<jupyter_output><empty_output><jupyter_text>Pull in AR data and set up partitions<jupyter_code>graph_AR, df_AR = build_graph("AR_shape/AR.shp")
AR_updaters = {"population": updaters.Tally("TOTPOP", alias="population"),
               "bvap": updaters.Tally("BVAP", alias="bvap"),
               "vap": updaters.Tally("VAP", alias="vap"),
               "bvap_prec": lambda part: {k: part["bvap"][k] / part["vap"][k]
                                          for k in part["bvap"]}}
AR_enacted_senate = GeographicPartition(graph_AR, assignment="SSD", updaters=AR_updaters)
num_opportunity_dists(AR_enacted_senate)
total_population_AR = sum(df_AR.TOTPOP.values)
ideal_population_AR = total_population_AR / 35
senate_seed = recursive_tree_part(graph_AR, range(35), pop_col="TOTPOP",
                                pop_target=ideal_population_AR,
                                epsilon=0.01, node_repeats=1)
AR_seed_senate = GeographicPartition(graph_AR, assignment=senate_seed,updaters=AR_updaters)
num_opportunity_dists(AR_seed_senate)
house_seed = recursive_tree_part(graph_AR, range(100), pop_col="TOTPOP",
                                pop_target=ideal_population_AR,
                                epsilon=0.01, node_repeats=1)
AR_seed_house = GeographicPartition(graph_AR, assignment=house_seed,updaters=AR_updaters)<jupyter_output><empty_output><jupyter_text>How does burst length affect findings<jupyter_code>plt.figure(figsize=(10,8))
plt.xlim(-.5, 8)
plt.xlabel("Number of opportunity districts")
plt.ylabel("Steps")
plt.title("PA short bursts of different lengths")

total_steps = 1000

for color, len_burst in [("k", 1), ("b", 5), ("r", 25), ("g", 50), ("y", 100),
                         ("cyan", 1000)]:
    print(len_burst)
    num_bursts =  int(total_steps / len_burst)
    _, observations = short_burst_run(enacted_senate, num_opportunity_dists,
                                      num_bursts=num_bursts, num_steps=len_burst)

    for i in range(num_bursts):
        plt.plot(observations[i], range(len_burst*i, len_burst*(i+1)),
                 color=color, alpha=0.5, marker=".", markevery=[0,-1])
    plt.plot([], color=color, label=("Burst_len " + str(len_burst)))
    
plt.legend()
plt.show()
plt.figure(figsize=(10,8))
plt.xlim(-.5, 8)
plt.xlabel("Number of opportunity districts")
plt.ylabel("Steps")
plt.title("AR short bursts of different lengths")

total_steps = 1000

for color, len_burst in [("k", 1), ("b", 5), ("r", 25), ("g", 50), ("y", 100),
                         ("cyan", 1000)]:
    print(len_burst)
    num_bursts =  int(total_steps / len_burst)
    _, observations = short_burst_run(AR_seed_senate, num_opportunity_dists,
                                      num_bursts=num_bursts, num_steps=len_burst, pop="TOTPOP")

    for i in range(num_bursts):
        plt.plot(observations[i], range(len_burst*i, len_burst*(i+1)),
                 color=color, alpha=0.5, marker=".", markevery=[0,-1])
    plt.plot([], color=color, label=("Burst_len " + str(len_burst)))
    
plt.legend()
plt.show()<jupyter_output>1
5
25
50
100
1000<jupyter_text>Run on PA state senate and visualizing results<jupyter_code>num_bursts = 50
num_steps = 10
max_part, observed_num_ops = short_burst_run(enacted_senate, num_opportunity_dists,
                                             num_bursts=num_bursts, num_steps=num_steps)
_, observed_num_ops_flat = short_burst_run(enacted_senate, num_opportunity_dists,
                                             num_bursts=1, num_steps=500)
plt.hist(observed_num_ops.flatten())
plt.hist(observed_num_ops_flat.flatten())
for i, row in enumerate(observed_num_ops):
    plt.figure()
    plt.title("Observations at burst " + str(i + 1))
    plt.xlabel("Number of opportunity districts")
    plt.ylabel("Frequency")
    plt.xlim(0,10)
    bins = np.linspace(0, 10, 11)
    sns.distplot(row, kde=False, bins=bins)
    #plt.savefig("plots/PA_senate_10_100_seed_3/short_burst_" + str(i+1) + "_hist.png")
names = list(map(lambda i: "Burst " + str(i), range(1,num_bursts + 1)))
plt.figure(figsize=(10,8))
plt.xlim(-.5, 8)
plt.xlabel("Number of opportunity districts")
plt.ylabel("Steps")
for i in range(num_bursts):
    plt.plot(observed_num_ops[i], range(num_steps*i, num_steps*(i+1)))#, label=names[i])
# plt.legend()
# plt.savefig("plots/burst_walk_100_100.png")
plt.figure(figsize=(10,8))
colors = ["red", "orange", "yellow", "lime", "green",
          "cyan", "blue", "indigo", "purple", "fuchsia"]
names = list(map(lambda i: "Burst " + str(i), range(1,num_bursts + 1)))
             
plt.hist(list(observed_num_ops), color=colors, label=names, align='mid', rwidth=10)
plt.xlim(1,7)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Reprojecting bursts onto the line  AR<jupyter_code>_, observed_num_ops_AR = short_burst_run(AR_seed_senate, num_opportunity_dists,
                                             num_bursts=1000, num_steps=10, pop="TOTPOP")
transitions_AR = []
for row in observed_num_ops_AR:
    transitions_AR.extend(bigrams(row))
fdist_AR = FreqDist(transitions_AR)
probs = {}
N = np.matrix(observed_num_ops_AR).astype(int)
dim = N.max()
A = np.zeros((dim, dim))
for k, v in fdist_AR.items():
    p = v / fdist_AR.N()
    probs[k] = p
    A[int(k[0]) - 1][int(k[1]) - 1] = p
A_normed = normalize(A, norm="l1")
G = nx.from_numpy_array(A_normed, create_using=nx.DiGraph)
mapping = {n: n+1 for n in G.nodes}
G1 = nx.relabel_nodes(G, mapping)
AR_line_proj=dict([((u,v,), round(d['weight'],4))
                 for u,v,d in G1.edges(data=True)])
AR_line_proj<jupyter_output><empty_output><jupyter_text>PA<jupyter_code>transitions = []
for row in observed_num_ops:
    transitions.extend(bigrams(row))
fdist = FreqDist(transitions)
fdist
probs = {}
N = np.matrix(observed_num_ops).astype(int)
dim = N.max()
A = np.zeros((dim, dim))
for k, v in fdist.items():
    p = v / fdist.N()
    probs[k] = p
    A[int(k[0]) - 1][int(k[1]) - 1] = p
A_normed = normalize(A, norm="l1")
A_normed[2]
G = nx.from_numpy_array(A_normed, create_using=nx.DiGraph)
G.edges(data=True)
mapping = {n: n+1 for n in G.nodes}
G1 = nx.relabel_nodes(G, mapping)
PA_line_proj=dict([((u,v,), round(d['weight'],4))
                 for u,v,d in G1.edges(data=True)])
PA_line_proj
plt.figure()
pos=nx.shell_layout(G1)
nx.draw(G1, pos=pos, edge_cmap=plt.cm.Reds)
nx.draw_networkx_edge_labels(G1,pos=pos,edge_labels=edge_labels)
plt.show()
H = nx.DiGraph()
for (n1,n2) in probs:
    H.add_edges_from([(n1, n2)], weight=probs[n1,n2])

edge_labels=dict([((u,v,), round(d['weight'],4))
                 for u,v,d in H.edges(data=True)])
edge_labels
plt.figure()
pos=nx.shell_layout(H)
nx.draw(H, pos=pos, edge_cmap=plt.cm.Reds)
nx.draw_networkx_edge_labels(H,pos=pos,edge_labels=edge_labels)
plt.show()

from networkx.drawing.nx_agraph import write_dot
plt.figure(figsize=(10,10))
nx.draw_networkx(
    G1, pos=nx.circular_layout(G1),
    node_shape='o', node_size=1000, font_size=10,
    edge_color='#555555', width=3.0
)
nx.draw_networkx_edge_labels(
    G1, pos=nx.circular_layout(G1), 
    edge_labels=edge_labels,
    font_size=10, label_pos=0.25, rotate=False
)
# plt.axis('off')
write_dot(G,'graph.dot')
from __future__ import division
import matplotlib as mpl

G = H
pos = nx.layout.spring_layout(G)

# node_sizes = [3 + 10 * i for i in range(len(G))]
M = G.number_of_edges()
edge_colors = range(2, M + 2)
edge_alphas = [round(d['weight'],4)*2 for u,v,d in H.edges(data=True)]

nodes = nx.draw_networkx_nodes(G, pos, node_color='blue')
edges = nx.draw_networkx_edges(G, pos, arrowstyle='->',
                               arrowsize=10, edge_color=edge_colors,
                               edge_cmap=plt.cm.Blues, width=2)
# set alpha value for each edge
for i in range(M):
    edges[i].set_alpha(edge_alphas[i])

pc = mpl.collections.PatchCollection(edges, cmap=plt.cm.Blues)
pc.set_array(edge_colors)
plt.colorbar(pc)

ax = plt.gca()
ax.set_axis_off()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12763664357197846, 'fraction_numerical': 0.018475750577367205, 'mean_word_length': 3.95272588638963, 'pattern_counts': {'":': 8, '<': 32, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '497248', 'license': 'MIT', 'n_tokens_mistral': 4924, 'n_tokens_neox': 4572, 'n_words': 913, 'path': 'misc/simple_random_walks_and_short_bursts.ipynb', 'repo_name': 'vrdi/shortbursts-gingles', 'size': '12990'}