Partitioning a colored grid

Question

I want to gerrymander a grid. What this means is that given an n x m grid with squares colored yellow and red, I want to partition the grid in such a way that yellow will be the majority color in as many partions as possible, as in this image:

All partitions must be continuous, the same number of squares, and all squares will be colored (although it would be awesome if an algorithm could generalize to grids with some squares not colored).

I'm not sure how to even go about 'algorithmitizing' this problem, past brute forcing every possible partition which is hard enough in and of itself and is incredibly inefficient.

What's the best way to accomplish this?

Answer 1

tl;dr: Used simulated annealing, swapping voters between districts. The demo at the bottom lets you do the swap step (Randomly evolve) and optimize for a gerrymandered district (anneal)

General Optimization

We can frame this as optimization problem, where we're trying to maximize the number of districts that red wins and minimize the number of districts that blue wins.

Let's formalize this:

function heuristic(state) {
  return state.districts_that_red_wins - state.districts_that_blue_wins;
}

where state is an assignment of voters to districts.

This will work, but could be improved a little bit. Let's introduce the notion of wasted votes to nudge our optimization in the right direction. We want to maximize wasted blue votes and minimize wasted red votes. I'm arbitrarily weighting them as being 1/10 as important as a district since there are 10 voters per district. That gives us a function to maximize:

function heuristic(state) {
  let {red_wins, blue_wins, tied, wasted_red_votes, wasted_blue_votes} = get_winners(state, voters);
  return red_wins - blue_wins - 0.1 * wasted_red_votes + 0.1 * wasted_blue_votes;
}

You may want to optimize other things e.g. compactness of districts. You can add these to the heursitic function.

Simulated Annealing

Let's pick an optimization algorithm to optimize state. We're subject to some constraints that make it hard to generate random maps that fit these conditions. And I suspect it's impossible to find the best district allocation without brute force, which would be impossible. So let's use an algorithm that lets us iteratively improve our solution.

I like simulated annealing because it's easy to implement and understand, and performs a little bit better than hill-climbing by preventing us from getting stuck in early local optima. My temperature function is simply max(0.8 - iterations/total_iterations, 0). At the beginning, 20% of the time we'll take a new state if and only if it's better; the other 80% we'll take the new state regardless. This slowly becomes more like hill-climbing until we're 80% of the way through our computation budget, then we only change the state if it improves our heuristic score. The choice of 80% is totally arbitrary.

To implement SA, we need an initial state (or a way of generating it). I'm going to use the "Perfect Representation" as the initial state for simplicity, mostly because I don't know how to generate random connected, equally-sized districts. We also need a way of making a small change to the state. I'll discuss this in the next section. Finally, we need a way of scoring the states. Let's use the function from the previous section because it's pretty cheap to compute.

If you're interested, look at the anneal function, or just read the Wikipedia article.

State Evolution

For this problem, given a state, we need to find another similar state that won't change the heuristic score very much to see if we're going in the right direction. I chose to find a pair of points, from two different districts and swap them.

We need to maintain some invariants:

• Keep all districts continuous
• Keep all districts the same size

The second one is easy: always swap points, never (permanently) assign from one district to another. The first is trickier, and we need a brief detour into graph theory. An articulation point (see picture) is a point that cannot be removed without bisecting a graph. For us, this means that we cannot remove articulation points without making a district discontinuous. Once we have a point that can be removed, we need to make sure it's added to a district it's adjacent to. This is pretty simple.

Since we're on a grid and all districts must be continuous, we can just consider the immediate neighbours of a point to determine whether it's an articulation point. If you can't see this, it's not super important, you can use the algorithm that works in general on a graph. I found the grid version easier because it involves no recursion.

See the is_swappable function if you're interested. This is what the "Randomly evolve" button in the demo does.

High-level, our code for evolving our state should look like this:

function evolve_state() {
    randomly pick a source district
    randomly pick a non-articulation point, source_point, from source_district

    for each neighbour of the articulation point
        if the neighbour is in a different district target_district
            temporarily remove source_point from source_district and add it to target_district
            if any articulation point (other than source point), target_point, in target_district is adjacent to source_district
                swap target_point and source_point
                return;
            restore source_point
}

Note: I implemented this in a way that randomly iterates through all source_district, source_point, neighbour, target_district and target_point because I wasn't sure how sparse this would be. If you implement this pseudocode exactly, you will likely need more iterations than I use to converge to a solution.

See evolve_state if you're interested.

Every function I haven't called out is a utility function or for drawing.

Demo

Now for a demo. :) (Uses Lodash for utility functions and Mithril for DOM manipulation)

If you want to play around with this, it might be easier to use my Plunker: http://plnkr.co/edit/Bho4qhQBKRShXWX8fHmt.

const RED = 'R';
const BLUE = 'B';
const VOTERS_PER_DISTRICT = 10;
const neighbours = [{x: 1, y: 0}, {x: 0, y: 1}, {x: -1, y: 0}, {x: 0, y: -1}];

/* UTILITY FUNCTIONS */

/**
 Create a generator that starts at a random point p, 0 <= p < max
 The generator will iterate over values p, p+1, ... max, 0, ... p-1
*/
function* cyclic_generator(max) {
    let start = _.random(max);

    for (let i=0; i<max; i++) {
        yield (i + start) % max;
    }
}

/**
  Return grid[x][y] if x and y are within grid. Otherwise return undefined
*/
function grid_get(grid, x, y) {
  if(_.isUndefined(grid[x])) {
    return undefined;
  }
  else {
    return grid[x][y];
  }
}

/** Generates a 2d array red and blue voters */
function generate_voters() {
  return _.times(5, x => _.times(10, () => {return {vote: x > 2 ? RED : BLUE, district_vote: 0xffffff}}))
}

/** Generate an initial state */ 
function generate_initial_state() {
  return _.range(5).map(x => _.range(10).map(y => {return {x, y}}));
}

/**
  Randomly swap two squares in the grid between two districts. 
  The new square to be added must be connected to the district, and the 
  old square must not break another district in two
 */
function evolve_state(state) {
  state = _.cloneDeep(state);
  
  // Create a grid with the district number
  let point_to_district = _.range(5).map(x => _.range(10).map(y => -1));
  state.forEach((district, i) => district.forEach(({x, y}) => point_to_district[x][y] = i));

  // swap a point from source_district to target_district.
  // then swap a point from target_district to source_district. 
  for(let source_district_idx of cyclic_generator(state.length)) {
    let source_articulation_points = state[source_district_idx].filter(point => is_swappable(point_to_district, point, source_district_idx));
    for(let source_point_idx of cyclic_generator(source_articulation_points.length)) {
      let source_point = source_articulation_points[source_point_idx];
      
      for(let neighbour_idx of cyclic_generator(4)) {
        let neighbour = neighbours[neighbour_idx];
        let target_district_idx = grid_get(point_to_district, source_point.x + neighbour.x, source_point.y + neighbour.y);
        if (_.isUndefined(target_district_idx) || target_district_idx == source_district_idx) {
          continue;
        }
        
        // swap the source point
        point_to_district[source_point.x][source_point.y] = target_district_idx;
        _.remove(state[source_district_idx], ({x, y}) => x == source_point.x && y == source_point.y); 
        // we don't add the point the the target array yet because we don't want to swap that point back
        
        // try to find a point in target_district that we can move to source_district
        let target_articulation_points = state[target_district_idx].filter(point => is_swappable(point_to_district, point, target_district_idx));
        for(let target_point_idx of cyclic_generator(target_articulation_points.length)) {
          let target_point = target_articulation_points[target_point_idx];
          for(let n of neighbours) {
            if(grid_get(point_to_district, target_point.x + n.x, target_point.y + n.y) === source_district_idx) {
              
              // found a point that we can swap!
              // console.log('swapping points!', source_point, target_point);
              
              _.remove(state[target_district_idx], ({x, y}) => x == target_point.x && y == target_point.y);
              state[target_district_idx].push(source_point);
              state[source_district_idx].push(target_point);
              return state;
            }
          }
        }
        
        // unswap source point since we were unsuccessful
        point_to_district[source_point.x][source_point.y] = source_district_idx;
        state[source_district_idx].push(source_point);
      }
    }
  }
  throw 'Could not find any states to swap' // this should never happen, since there will always be the option of reversing the previous step
}

/*
  Return whether a point can be removed from a district without creating disjoint districts. 
  In graph theory, points that cannot be removed are articulation points. 
  For a general algorithm, see: https://stackoverflow.com/questions/15873153/explanation-of-algorithm-for-finding-articulation-points-or-cut-vertices-of-a-gr
  
  My version takes advantage of the fact that we're on a grid and that all the districts must be continuous, 
  so we can consider only the immediate neighbours of a point.
*/
function is_swappable(grid, p, district) {
  
  // if the the point is not even in this district, it makes no sense for this to consider this point at all
  if(grid[p.x][p.y] != district) {
    return false;
  }
  
  // if two opposite edges are part of this district, this is an articulation point
  // .x.  x is an articulation point
  // Exception:
  // .x.  x is not an articulation point
  // ...
  if (grid_get(grid, p.x+1, p.y) === district && grid_get(grid, p.x-1, p.y) === district && grid_get(grid, p.x, p.y+1) !== district && grid_get(grid, p.x, p.y-1) !== district) {
    return false;
  }
  if (grid_get(grid, p.x, p.y+1) === district && grid_get(grid, p.x, p.y-1) === district && grid_get(grid, p.x+1, p.y) !== district && grid_get(grid, p.x-1, p.y) !== district) {
    return false;
  }
  
  // check if any corners are missing:
  // .x  x is not an articulation point      .x   x is an articulation point
  // ..                                       .
  for(let i = 0; i < 4; i++) {
    let nx = neighbours[i].x;
    let ny = neighbours[i].y;
    let nx2 = neighbours[(i+1)%4].x;
    let ny2 = neighbours[(i+1)%4].y;
    
    if (grid_get(grid, p.x+nx, p.y+ny) === district && grid_get(grid, p.x+nx2, p.y+ny2) === district && grid_get(grid, p.x+nx+nx2, p.y+ny+ny2) !== district) {
      return false;
    }
  }
  return true;
}

/** Count how many districts each party wins */
function get_winners(state, voters) {
  let red_wins = 0;
  let blue_wins = 0;
  let tied = 0;
  
  let wasted_red_votes= 0; // see https://en.wikipedia.org/wiki/Wasted_vote
  let wasted_blue_votes = 0;
  
  state.forEach(district => {
    let counts = _.countBy(district.map(({x, y}) => voters[x][y].vote))
    
    if ((counts[BLUE] || 0) > (counts[RED] || 0)) {
      blue_wins++;
      wasted_blue_votes += (counts[BLUE] || 0) - VOTERS_PER_DISTRICT / 2 - 1;
      wasted_red_votes += (counts[RED] || 0);
    }
    else if ((counts[RED] || 0) > (counts[BLUE] || 0)) {
      red_wins++;
      wasted_red_votes += (counts[RED] || 0) - VOTERS_PER_DISTRICT / 2 - 1;
      wasted_blue_votes += (counts[BLUE] || 0);
    }
    else {
      tied++;
    }
  });
  return {red_wins, blue_wins, tied, wasted_red_votes, wasted_blue_votes};
}

/* GUI */

/* Display a grid showing which districts each party won */
function render_districts(state, voters) {
  let red_districts = 0;
  let blue_districts = 0;
  let grey_districts = 0;
  
  // Color each district
  state.forEach(district => {
    let counts = _.countBy(district.map(({x, y}) => voters[x][y].vote))

    let district_color;
    if ((counts[BLUE] || 0) > (counts[RED] || 0)) {
      district_color = 'blue' + blue_districts++;
    }
    else if ((counts[RED] || 0) > (counts[BLUE] || 0)) {
      district_color = 'red' + red_districts++;
    }
    else {
      district_color = 'grey' + grey_districts++;
    }
    
    district.map(({x, y}) => voters[x][y].district_color = district_color);
  });
  
  return m('table', [
    m('tbody', voters.map(row => 
      m('tr', row.map(cell => m('td', {'class': cell.district_color}, cell.vote)))
    ))
  ]);
}

/** Score a state with four criteria: 
  - maximize number of red districts
  - minimize number of blue districts
  - minimize number of red voters in districts that red wins
  - maximize number of blue voters in districts that blue wins
  
  The first two criteria are arbitrarily worth 10x more than the latter two
  The latter two are to nudge the final result toward the correct solution
 */
function heuristic(state) {
  let {red_wins, blue_wins, tied, wasted_red_votes, wasted_blue_votes} = get_winners(state, voters);
  return red_wins - blue_wins - 0.1 * wasted_red_votes + 0.1 * wasted_blue_votes;
}

/**
 Optimization routine to find the maximum of prob_fcn.
 prob_fcn: function to maximize. should take state as its argument
 transition: how to generate another state from the previous state
 initialize_state: a function that returns an initial state
 iters: number of iterations to run
 
 Stolen from my repo here: https://github.com/c2huc2hu/automated-cryptanalysis/blob/master/part3.js
*/
function anneal(prob_fcn, transition, initialize_state, seeds=1, iters=1000) {
    let best_result = initialize_state();

    for(let i=0; i<seeds; i++) {
        let curr_state = initialize_state();
        let curr_cost = prob_fcn(curr_state);

        // perform annealing. do a few extra steps with temp=0 to refine the final solution
        for(let j=0; j<iters; j++) {
            let candidate_state = transition(curr_state);
            let candidate_cost = prob_fcn(candidate_state);
            temp = 0.8 - j / iters;

            if(candidate_cost >= curr_cost || Math.random() < temp) {
                curr_state = candidate_state;
                curr_cost = candidate_cost;
            }
        }

        if(prob_fcn(curr_state) > prob_fcn(best_result)) {
            best_result = curr_state;
        }
    }

    return best_result;
}

let voters = generate_voters();
let state = generate_initial_state();

// main rendering code: this code renders the UI
m.mount(document.getElementById('actions'), {view: function() {
  return m('div', [
    m('button', {onclick: () => state = generate_initial_state()}, 'Reset'),
    m('button', {onclick: () => state = evolve_state(state)}, 'Randomly evolve'), // randomly evolves

    m('br'),
    m('label', {'for': 'radio-blue'}, 'Gerrymander for blue'),
    m('input', {type: 'radio', name: 'heuristic', value: 'blue', id: 'radio-blue'}),
    m('label', {'for': 'radio-red'}, 'Gerrymander for red'),
    m('input', {type: 'radio', name: 'heuristic', value: 'red', id: 'radio-red'}),
    m('br'),
    m('label', {'for': 'anneal-steps'}, 'Anneal steps: '),
    m('input', {id: 'anneal-steps', type: 'number', value: '500'}),
    m('button', {onclick: function() {
      let minimize = document.getElementById('radio-red').checked;
      let _heuristic = minimize ? heuristic : state => -heuristic(state)
      
      let new_state = anneal(_heuristic, evolve_state, generate_initial_state, 1, parseInt(document.getElementById('anneal-steps').value));
      if(_heuristic(new_state) > _heuristic(state)) {
        state = new_state;
      }
      else {
        console.log('found no better solutions')
      }
    }}, 'Anneal!'),
  ]);
}});

// This renders the grid
m.mount(document.getElementById('grid'), {
  view: function() {
    return render_districts(state, voters)
  }
});

// state = anneal(heuristic, evolve_state, generate_initial_state, 5, 1000);
document.getElementById('radio-red').checked = true;
m.redraw();

/* Layout */

table {
  border: solid 1px black;
}

td {
  padding: 5px;
  border: solid 1px black;
}

button {
  margin: 10px;
}

p {
  max-width: 500px;
}

/* Colour classes. In hindsight, this wasn't a good idea */
.red0 {
  background-color: red;
}
.red1 {
  background-color: darkred;
}
.red2 {
  background-color: pink;
}
.red3 {
  background-color: deeppink;
}
.red4 {
  background-color: lightsalmon;
}

.blue0 {
  background-color: aqua;
}
.blue1 {
  background-color: cadetblue;
}
.blue2 {
  background-color: steelblue;
}
.blue3 {
  background-color: royalblue;
}
.blue4 {
  background-color: midnightblue;
}

.grey0 {
  background-color: lightgrey;
}
.grey1 {
  background-color: silver;
}
.grey2 {
  background-color: darkgray;
}
.grey3 {
  background-color: gray;
}
.grey4 {
  background-color: dimgray;
}

<!DOCTYPE html>
<html>

  <head>
    <script data-require="lodash.js@4.17.4" data-semver="4.17.4" src="https://cdn.jsdelivr.net/npm/lodash@4.17.4/lodash.min.js"></script>
    <script data-require="mithril@1.0.1" data-semver="1.0.1" src="https://cdnjs.cloudflare.com/ajax/libs/mithril/1.0.1/mithril.js"></script>
    <link rel="stylesheet" href="style.css" />
  </head>

  <body>
    <h1>Gerrymandering simulation</h1>
    
    <p>
      There are two parties, red and blue (chosen because they contrast well).
      Each person will always vote a certain way, and is marked with R or B in the table. 
      People are divided into districts, shown here as groups of people marked in a single colour.
    </p>
    
    <p>
      Use the buttons below to divide up districts.
      The reset button will restore the initial state.
      The randomly-evolve button will swap two people between districts
      The anneal button will optimize for your chosen party. 
      You should limit the number of steps to ~1000 or your browser will appear to hang.
      In general, it is sufficient to run a few seeds for 500 iterations. 
    </p>
    
    <div id="grid"></div>
    
    <div id="actions"></div>
    
    <script src="script.js"></script>
  </body>

</html>