redistricting_frontiers.py

"""
The goal of this software is to run a simulation of a model presented by
Moon Duchin in the math session at the Austin regional site for the Metric
Geometry and Gerrymandering Group


In broad strokes, it takes an input shape of some sort and splits it
in four at some point in the polygon (don't talk to me about points out
of the polygon, you degenerate). It assumes that the polygon has some proportion
of one delegate's worth of population, and that the shapes resulting from the split
have their populations "filled in" with population from the surrounding area.
Also, people from the surrounding areas are universally opposed to people from
within the split polygon.
"""


import argparse
from concurrent.futures import ThreadPoolExecutor, as_completed
from datetime import datetime
from functools import reduce
import random

import matplotlib.pyplot as plt
from shapely.ops import split
from shapely.geometry import (
    Point,
    LineString,
    MultiPoint,
    Polygon
)


def split_polygon_at_point(poly, point):
    """Split a polygon by the vertical and horizontal lines through a point"""
    (x, y) = point.coords[0]
    splitline1 = LineString([(x, 0), (x, 1)])
    splitline2 = LineString([(0, y), (1, y)])
    return reduce(
        lambda x, y: x + y,
        map(lambda x: list(split(x, splitline2)), list(split(poly, splitline1)))
    )


def poly_partisan_balance(population, poly, orig_poly):
    """Calculate the proportion of voters in a polygon that lean toward a party

    It is calculated by the ratio of a polygon's area to the original polygon's
    area, multiplied by the number of representatives worth of population.
    """

    return population * (poly.area / orig_poly.area)


def cracked(n, poly, polys, orig_poly, cracked_range):
    """Calculate whether a polygon is part of cracking

    We have cracking if there are bordering polygons that are both within the
    range set in cracked_range.
    """
    poly_alignment = poly_partisan_balance(n, poly, orig_poly)
    if cracked_range[0] < poly_alignment < cracked_range[1]:
        for other in polys:
            other_alignment = poly_partisan_balance(n, other, orig_poly)
            # poly.relate(other)[4] is the DE-9IM value for the relationship between
            # the two polygons' boundaries. '1' is a line intersection, which is enough to
            # indicate that they border each other.
            if (
                    poly.relate(other)[4] == '1'
                    and cracked_range[0] < other_alignment < cracked_range[1]
            ):
                # We can return once we get a single cracked neighbor, no need to evaluate others
                return True
    return False


def packed(n, poly, orig_poly, threshold):
    """Calculate whether a polygon is packed

    A polygon is packed if its alignment is greater than the threshold
    """

    return poly_partisan_balance(n, poly, orig_poly) > threshold


def run_point(n, startPoly, cracked_range, packed_threshold):
    """Test the results of splitting the polygon at a point

    To sample, check whether the point is inside the polygon. If it is, split
    the polygon at that point, and evaluate the polygons that result for packedness
    and crackedness.
    """

    (xmin, ymin, xmax, ymax) = startPoly.bounds
    x = random.uniform(xmin, xmax)
    y = random.uniform(ymin, ymax)
    p = Point(x, y)
    if not startPoly.contains(p):
        return
    polys = split_polygon_at_point(startPoly, p)
    any_packed = any([packed(n, poly, startPoly, packed_threshold) for poly in polys])
    any_cracked = any(
        [cracked(n, poly, [p for p in polys if p != poly], startPoly, cracked_range) for poly in polys]
    )
    if any_packed and any_cracked:
        color = 'red'
    elif any_packed:
        color = 'orange'
    elif any_cracked:
        color = 'yellow'
    else:
        color = 'green'
    return {'point': (x, y), 'color': color}


def plot_points(points):
    """Put the points in a plot"""
    non_null_points = [x for x in points if x is not None]
    (xmin, ymin, xmax, ymax) = MultiPoint([x['point'] for x in non_null_points]).bounds
    fig, ax = plt.subplots()
    colors = set([x['color'] for x in non_null_points])
    for c in colors:
        filtered = [p for p in non_null_points if p['color'] == c]
        xs = [p['point'][0] for p in filtered]
        ys = [p['point'][1] for p in filtered]
        ax.scatter(xs, ys, color=c, s=1)
    ax.set_xlim(xmin, xmax)
    ax.set_ylim(ymin, ymax)
    plt.show()


def make_random_line():
    """Create a line between two random points"""
    return LineString([
        (random.random(), random.random()),
        (random.random(), random.random())
    ])


def weird_shape(n_lines, buffer_distance):
    """Weird shape generates a weird blobby polygon from buffered random lines"""
    polys = [make_random_line().buffer(buffer_distance) for _ in range(n_lines)]
    return reduce(lambda x, y: x.union(y), polys).simplify(tolerance=0.05)


# TODO
# collect results instead of plotting them, to run many trials

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument(
        'N',
        help='How much population the region should have as a ratio over one delegate',
        type=float
    )
    parser.add_argument(
        'shape',
        help='What the starting shape should be',
        choices=['circle', 'rectangle', 'strange', 'creepyfingers'],
        default='circle'
    )
    blob_generation_lines_help = (
        'How many lines to use to generate a strange shape. Ignored if the shape is not "strange"'
    )
    blob_generation_buffer_help = (
        'How far to buffer from each blob generation line. If not provided and generation lines is'
        ', will default to 1 / blob-generation-lines'
    )
    parser.add_argument('--num-runs', '-r', type=int, default=10000,
                        help='How many points to sample')
    parser.add_argument('--cracked-min', type=float, default=0.3,
                        help='Floor of the cracked districts threshold')
    parser.add_argument('--cracked-max', type=float, default=0.42,
                        help='Ceiling of the cracked districts threshold')
    parser.add_argument('--packed-threshold', type=float, default=0.7,
                        help='Ceiling of the cracked districts threshold')
    parser.add_argument('--blob-generation-lines', '-b', default=10, type=int,
                        help=blob_generation_lines_help)
    parser.add_argument('--blob-generation-buffer', type=float, default=0.15,
                        help=blob_generation_buffer_help)
    parser.add_argument('--seed', type=int, default=1,
                        help='Random seed to use for reproducibility')
    args = parser.parse_args()
    if args.seed != parser.get_default(args.seed):
        random.seed(args.seed)
    blob_generation_lines = args.blob_generation_lines
    calculate_generation_buffer = (
        blob_generation_lines != parser.get_default('blob_generation_lines')
        and args.blob_generation_buffer == parser.get_default('blob_generation_buffer')
        and args.shape == 'strange'
    )
    if calculate_generation_buffer:
        blob_generation_buffer = 1 / float(blob_generation_lines)
    else:
        blob_generation_buffer = args.blob_generation_buffer
    cracked_range = (args.cracked_min, args.cracked_max)
    packed_threshold = args.packed_threshold
    if args.shape == 'circle':
        orig_poly = Point(0.5, 0.5).buffer(0.5)
    elif args.shape == 'rectangle':
        orig_poly = Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
    elif args.shape == 'strange':
        orig_poly = weird_shape(blob_generation_lines, blob_generation_buffer)
    # elif args.shape == 'creepyfingers':
    else:
        # Start with two connected points
        # draw a buffered line between them
        # draw a random number of lines outward from each point
        # maybe continue based on some depth parameter
        raise NotImplementedError('chill bro')
    print(datetime.now().isoformat())
    points = []
    with ThreadPoolExecutor(max_workers=4) as executor:
        point_futures = [
            executor.submit(run_point, args.N, orig_poly, cracked_range, packed_threshold)
            for _ in range(args.num_runs)
        ]
        for future in as_completed(point_futures):
            data = future.result()
            points.append(data)
    print(datetime.now().isoformat())
    plot_points([x.result() for x in point_futures])


if __name__ == '__main__':
    main()

requirements.txt

matplotlib==2.1.1
shapely==1.6.4.post1