Advent of Code 2020: Day 20 with NetworkX, SciPy cross-correlation in Python
Tags: adventofcode, python
Published: 2020-12-21 11:40:27

I have to admit, this problem kicked my ass

# The Challenge Part 1
[Link to challenge on Advent of Code 2020 website](https://adventofcode.com/2020/day/20)

The challenge is to take same image tiles and arrange them so their borders line up correctly. In image processing, this is a mosaic stitching problem, used in, for example, stitching together telescope images of stars and stitching together a mosaic from a large document (for example in forensic reconstruction of a shredded document). There's also a field of image processing with names like multi-view geometry, photogrammetry, and homography estimation, which relate to doing this kind of image processing but in 3D (which is something I've used in my job for several years).

Unfortunately it turns out that my job knowledge in reconstructing 3D scenes from drone photos didn't come in handy as our domain is mainly based on reconstructing 3D structure, rather than this kind of matching. So we're back to basics.

## Importing the data
The data is similar to Day 17's data, so we import it in a similar way, using the numpy's `view()` method to split up the text. We also use python's destructuring syntax `k, *v, _ = some_array` to break down some_array so that `k` holds the first line, `_` holds the last line, and `v` holds everything inbetween. This lets us break down the input into individual tiles, their titles, and also strip out the newline between them. Sprinkle a little regex to strip out the tile IDs:

Input data:
```
Tile 2311:
..##.#..#.
##..#.....
#...##..#.
####.#...#
##.##.###.
##...#.###
.#.#.#..##
..#....#..
###...#.#.
..###..###

Tile 1951:
#.##...##.
#.####...#
.....#..##
#...######
.##.#....#
.###.#####
###.##.##.
.###....#.
..#.#..#.#
#...##.#..

<snip>
```
Code:
```py
import re
import numpy as np

data = open("input.txt").read().splitlines()
tiles = {int(re.sub(r"Tile (\d+):", r"\1", k)): np.array(v).view('U1').reshape(10, 10) for k, *v, _ in  np.array(data).reshape(-1, 12)}
```
Output (of `tiles`)
```
{2311: array([['.', '.', '#', '#', '.', '#', '.', '.', '#', '.'],
        ['#', '#', '.', '.', '#', '.', '.', '.', '.', '.'],
        ['#', '.', '.', '.', '#', '#', '.', '.', '#', '.'],
        ['#', '#', '#', '#', '.', '#', '.', '.', '.', '#'],
        ['#', '#', '.', '#', '#', '.', '#', '#', '#', '.'],
        ['#', '#', '.', '.', '.', '#', '.', '#', '#', '#'],
        ['.', '#', '.', '#', '.', '#', '.', '.', '#', '#'],
        ['.', '.', '#', '.', '.', '.', '.', '#', '.', '.'],
        ['#', '#', '#', '.', '.', '.', '#', '.', '#', '.'],
        ['.', '.', '#', '#', '#', '.', '.', '#', '#', '#']], dtype='<U1'),
 1951: array([['#', '.', '#', '#', '.', '.', '.', '#', '#', '.'],
        ['#', '.', '#', '#', '#', '#', '.', '.', '.', '#'],
        ['.', '.', '.', '.', '.', '#', '.', '.', '#', '#'],
        ['#', '.', '.', '.', '#', '#', '#', '#', '#', '#'],
        ['.', '#', '#', '.', '#', '.', '.', '.', '.', '#'],
        ['.', '#', '#', '#', '.', '#', '#', '#', '#', '#'],
        ['#', '#', '#', '.', '#', '#', '.', '#', '#', '.'],
        ['.', '#', '#', '#', '.', '.', '.', '.', '#', '.'],
        ['.', '.', '#', '.', '#', '.', '.', '#', '.', '#'],
        ['#', '.', '.', '.', '#', '#', '.', '#', '.', '.']], dtype='<U1'),
```

## Matching the tiles
To match the tiles, we have to pull out each edge and compare it to each other edge of the tile. While the question doesn't say it, but there are no mis-matches! 

To generate the border ID, we simply take a slice out of the 2D array of a tile, and turn it into a string. For example, the right edge of tile 2311 (one of the examples above) can be gotten with:

```py
"".join(tiles[2311][..., -1])
```
Output
```
'...#.##..#'
```

However, we also need to be able to handle this border when the tile is flipped over. We can achieve this using python's extended slice notation to slice backwards:
```py
"".join(tiles[2311][..., -1][::-1])
```
Output
```
'#..##.#...'
```

The right edge is the index `(..., -1)` with `...` (python's Ellipses object) meaning "every value along this axis", and `-1` meaning "the last index on this axis". We can store all four edges in a tuple for convenience:

```py
edges = (
    (..., -1 ),
    (0  , ...),
    (...,  0 ),
    (-1 , ...),
)
```

And we can iterate over them in a For loop. Since we also need their flipped forms, we can either use a double-for loop, or the itertools.product() function:

```py
from itertools import product
for (i, j), f in product(edges, (1, -1)):
    border_id = "".join(tile[i, j][::f])
```

This loop generates the 8 combinations of each of the 4 edges, with each of the 2 flip directions.

For each border_id, we'll stick it in a dictionary along with the tile name, so that next time the border_id is seen again, we can match up the two tiles and create a graph of all the matches. I'm again using the `networkx` library, a library that I've used extensively through this series:

```py
import networkx as nx
from itertools import product

edges = (
    (..., -1 ),
    (0  , ...),
    (...,  0 ),
    (-1 , ...),
)

borders = {}
graph = nx.Graph()
for tile_name, tile in tiles.items():
    for (i, j), f in product(edges, (1, -1)):
        border_id = "".join(tile[i, j][::f])
            
        if border_id in borders:
            graph.add_edge(borders[border_id], tile_name)
        else:
            borders[border_id] = tile_name
```

We can visualize the resulting graph, and actually this accidentally solves the puzzle by directly plotting what the arrangement should look like:

```py
pos = nx.nx_agraph.graphviz_layout(graph)
nx.draw(graph, pos=pos, node_color="lightgrey", edge_color="grey", node_size=50, with_labels=True)
```
(from my puzzle input)

![SPOILER_unknown](https://dev-to-uploads.s3.amazonaws.com/i/e8jv0rztcpdfja858yb7.png)

The challenge asks for the IDs of the four corners multiplied together. We can actually just read them off this graph!

But instead, we can also find the corners by looking at the nodes of the graph with only two edges:

```py
from math import prod
print("prod", prod(tile for tile in graph if len(graph[tile]) == 2))
```

The full code for this part of the challenge:
```py
import numpy as np
import networkx as nx
from itertools import product
from math import prod
import re

data = open("input.txt").read().splitlines()
tiles = {int(re.sub(r"Tile (\d+):", r"\1", k)): np.array(v).view('U1').reshape(10, 10) for k, *v, _ in  np.array(data).reshape(-1, 12)}

edges = (
    (..., -1 ),
    (0  , ...),
    (...,  0 ),
    (-1 , ...),
)

borders = {}
graph = nx.Graph()
for tile_name, tile in tiles.items():
    for (i, j), f in product(edges, (1, -1)):
        border_id = "".join(tile[i, j][::f])
            
        if border_id in borders:
            graph.add_edge(borders[border_id], tile_name)
        else:
            borders[border_id] = tile_name
            
print("prod", prod(tile for tile in graph if len(graph[tile]) == 2))
```

# The Challenge Part 2
Oh boy, this is a big one. This challenge contains multiple steps. Now that we have our graph of tiles, we must first re-assemble those tiles into an image. Once assembled, borders must be trimmed off. Then we must find occurrences of a sea monster inside the image. Then, we find how many "#" characters make up the remainder of the image...

## Reassembling the tiles into an image
This is by far the most tedious step. Since so far we know which tile connects to which other tile, and on which borders, in theory we can work out exactly which way up each tile needs to be in order to join up seamlessly. Unfortunately this requires too much finesse in geometry, and I can't quite wrap my head around the transformation needed to work out which direction, rotation, and flip to draw a tile given knowledge of which pair of borders match, and also the global rotation and flip of the tile it's connected to.

So instead, I'm just going to plop down a tile, and then brute-force match each of the 4 edges with the 8 possible orientations of the next tile, knowing that only one of these will match. It's still linear-time solve, just on average 16 times slower.

To do this, I need to be able to iterate over the 8 orientations of the tiles, so I create this generator

```py
def orient(tile):
    for i in range(4):
        yield np.rot90(tile, i)
        yield np.fliplr(np.rot90(tile, i))
```

Next, because I only want to compare an edge in one tile against its corresponding edge in another tile, I create a flip lookup table, that gives me the flipped edge.

Recall `edge` looks like this:
```py
edges = (
    (..., -1 ),
    (0  , ...),
    (...,  0 ),
    (-1 , ...),
)
```
It turns out, that to transform each of these edges to its opposite one, I simply need to swap any `-1` for `0` and `0` for `-1`. I can do this using the following lookup table (as odd a syntax as it looks):

```py
flip = {
    ...: ...,
    -1 :  0 ,
     0 : -1 ,
}
```
This way, I can compare tiles like so:
```
for i, j in edges:
  if str(old_tile[i, j]) == str(new_tile[flip[i], flip[j]]):
```

Next, once a match is found, I need to store the data. Because I'm not orienting the first tile, and I'm not selecting any particular corner, the grid can grow out in any direction. It would be more memory-efficient to find a corner, and orient it correctly so that the grid grows right and down to fill an exact chunk of memory.

But I'm lazy. Since the data takes up 12x12, if I create a large enough space, and start in the middle, I can guarantee that I don't exceed the bounds whichever direction the grid grows in. So I create a 25x25x10x10 4-dimensional numpy grid to store the 10x10 tiles in a 2D arrangement.

To store the data, I also need to find the direction in the grid to store it in. If the matched edge is the right-side, I need to draw the new tile into the space to the right of the old one. If the matched edge is below, then the new tile needs to go below. It turns out I can convert our edge array into a direction vector with another lookup table:

```py
dire = {
    ...:  0 ,
    -1 :  1 ,
     0 : -1 ,
}
```

For example, the right-edge given by (..., -1) maps to vector (0, 1) using this lookup table.

The full code looks like:

```py
from math import sqrt

dim = int(sqrt(len(tiles)))

first_tile = next(iter(tiles))
result = np.full([2*dim+1, 2*dim+1, 10, 10], "")
result[dim][dim] = tiles[first_tile]
location = {first_tile: np.array([dim, dim])}

for left, right in nx.dfs_edges(graph, first_tile):
    if left in location:
        old, new = left, right
    else:
        old, new = right, left
        
    old_loc = location[old]
    old_tile = result[old_loc[0], old_loc[1]]
    for i, j in edges:
        for oriented_tile in orient(tiles[new]):
            if str(old_tile[i, j]) == str(oriented_tile[flip[i], flip[j]]):
                new_loc = old_loc + np.array([dire[i], dire[j]])
                result[new_loc[0], new_loc[1]] = oriented_tile
                location[new] = new_loc
                break
```
The output here is a giant 24x24x10x10 numpy array with the assembled grid of tiles somewhere in the middle.

## Clipping the borders and forming the full picture
The challenge requires us to remove the borders, so each tile goes from a 10x10 array to an 8x8 one. Since we constructed the results as a 25x25 grid of 10x10 arrays, we can easily do this in one go:

```py
clipped = result[:,:,1:-1,1:-1]
```

To form the full picture, we need to collapse our 4D array into a 2D one. Unfortunately some care is needed as we have to collapse the axes in the right order to form the image.
```py
flat = np.swapaxes(clipped, 1, 2).reshape((2*dim+1)*8, -1)
```

Next, the 25x25 grid contains a lot of extra space, I can remove these by just selecting everything that's empty. This results in a 1D array however, but that can be re-shaped easily:

```py
flat_filtered = flat [~(flat == '')].reshape(dim*8, -1)
```
Output (displayed here with the 3x3 demo data):
```
array([['.', '.', '.', '#', '#', '#', '.', '.', '.', '#', '#', '.', '.', '.', '#', '.', '.', '.', '#', '.', '.', '#', '#', '#'],
       ['.', '#', '.', '#', '#', '#', '.', '.', '#', '#', '.', '.', '#', '#', '.', '.', '#', '#', '#', '#', '.', '#', '#', '.'],
       ['#', '.', '#', '#', '.', '.', '#', '.', '.', '#', '.', '.', '.', '#', '.', '.', '#', '#', '#', '#', '.', '.', '.', '#'],
       ['#', '#', '#', '#', '#', '.', '.', '#', '#', '#', '#', '#', '.', '.', '.', '#', '#', '#', '.', '.', '.', '.', '#', '#'],
       ['#', '.', '.', '#', '#', '#', '#', '.', '.', '.', '#', '.', '#', '.', '#', '.', '#', '#', '#', '.', '#', '#', '#', '.'],
       ['.', '.', '#', '.', '#', '.', '.', '#', '.', '.', '#', '.', '#', '.', '#', '.', '#', '#', '#', '#', '.', '#', '#', '#'],
       ['.', '#', '#', '#', '#', '.', '#', '#', '#', '.', '#', '.', '.', '.', '#', '#', '#', '.', '#', '.', '.', '#', '.', '#'],
       ['.', '#', '.', '#', '.', '#', '#', '#', '.', '#', '#', '.', '#', '#', '.', '#', '.', '.', '#', '.', '#', '#', '.', '.'],
       ['#', '#', '#', '.', '#', '.', '.', '.', '#', '.', '.', '#', '.', '#', '#', '.', '#', '#', '#', '#', '#', '#', '.', '.'],
       ['.', '#', '.', '#', '.', '.', '.', '.', '#', '.', '#', '#', '.', '#', '.', '.', '.', '#', '#', '#', '.', '#', '#', '.'],
       ['.', '.', '.', '#', '.', '.', '#', '.', '.', '#', '.', '#', '.', '#', '#', '.', '.', '#', '#', '#', '.', '#', '#', '#'],
       ['#', '#', '.', '.', '#', '#', '.', '#', '.', '.', '.', '#', '.', '.', '.', '#', '.', '#', '.', '#', '.', '#', '.', '.'],
       ['#', '.', '#', '#', '#', '#', '.', '.', '.', '.', '#', '#', '.', '.', '#', '#', '#', '#', '#', '#', '#', '#', '.', '#'],
       ['#', '#', '#', '.', '#', '.', '#', '.', '.', '.', '#', '.', '#', '#', '#', '#', '#', '#', '.', '#', '.', '.', '#', '#'],
       ['#', '.', '#', '#', '#', '#', '.', '.', '#', '.', '#', '#', '#', '#', '.', '#', '.', '#', '.', '#', '#', '#', '.', '.'],
       ['#', '.', '.', '#', '.', '#', '#', '.', '.', '#', '.', '.', '#', '#', '#', '.', '#', '.', '#', '#', '.', '.', '.', '.'],
       ['.', '#', '#', '#', '#', '.', '.', '.', '#', '.', '.', '#', '.', '.', '.', '.', '.', '#', '.', '.', '.', '.', '.', '.'],
       ['.', '.', '.', '.', '#', '.', '.', '#', '.', '.', '.', '#', '#', '.', '.', '#', '.', '#', '.', '#', '#', '#', '.', '.'],
       ['.', '.', '.', '#', '#', '#', '#', '#', '#', '#', '#', '.', '#', '.', '.', '.', '.', '#', '#', '#', '#', '#', '.', '#'],
       ['#', '#', '.', '#', '.', '.', '.', '.', '#', '.', '#', '#', '.', '#', '#', '#', '#', '.', '.', '.', '#', '.', '#', '#'],
       ['#', '#', '#', '.', '#', '#', '#', '#', '#', '.', '.', '.', '#', '.', '#', '#', '#', '#', '#', '.', '#', '.', '.', '#'],
       ['#', '#', '.', '#', '#', '.', '#', '#', '#', '.', '#', '.', '#', '.', '.', '#', '#', '#', '#', '#', '#', '.', '.', '.'],
       ['#', '#', '#', '.', '.', '.', '.', '#', '.', '#', '.', '.', '.', '.', '#', '.', '.', '#', '.', '.', '.', '.', '.', '.'],
       ['.', '#', '.', '#', '.', '.', '#', '.', '#', '#', '.', '.', '.', '#', '.', '#', '#', '.', '.', '#', '#', '#', '#', '#']], dtype='<U1')
```
But actually we want 1s and 0s to for the next step:
```py
image = (flat_filtered == "#").astype(int)
```
Output
```
array([[0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1],
       [0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0],
       [1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1],
       [1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1],
       [1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0],
       [0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1],
       [0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1],
       [0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0],
       [1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0],
       [0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1],
       [1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0],
       [1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1],
       [1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1],
       [1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0],
       [1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0],
       [0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0],
       [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1],
       [1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1],
       [1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1],
       [1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
       [1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1]])
```

## Importing the monster
Importing the monster is done similarly to importing the tiles. It's saved to a text file, and imported with numpy. We also need to know the total number of "#" symbols in it for later.

```py
monster = (np.array(open("monster.txt").read().splitlines()).view('U1') == "#").reshape((3, -1)).astype(int)
monster_value = monster.sum()
```

## Finding the monster
Finally, in order to find the monster, we just need to run a 2D Cross-correlation between the monster image (in all its different possible orientations) and the main image. Scipy has a nice cross-correlation function. It essentially scans the monster image across the main image, and sums up where they match. We expect a match to equal the sum of "#" in the monster, which we calculated earlier.

```py
from scipy.ndimage import correlate
correlate(image, monster, mode="constant") == monster_value
```
This returns a True or False for every position where a monster is detected. Unfortunately, we get all False because the image is in the wrong orientation. We have to loop over all orientations of the image (or monster), which we can do with the `orient()` function used earlier.

```py
total_monsters = sum((correlate(image, oriented_monster, mode="constant") == monster_value).sum() for oriented_monster in orient(monster))
```

Finally, once we have the number of monsters, we can count the `#` in the image and remove the number of monsters from it

```py
print("remaining", image.sum() - total_monsters*monster_value)
```

The full code is (note: this reproduces part of the code from Part 1):
```py
import numpy as np
import networkx as nx
from itertools import product
from math import prod, sqrt
import re
from scipy.ndimage import correlate

data = open("input.txt").read().splitlines()
tiles = {int(re.sub(r"Tile (\d+):", r"\1", k)): np.array(v).view('U1').reshape(10, 10) for k, *v, _ in  np.array(data).reshape(-1, 12)}

edges = (
    (..., -1 ),
    (0  , ...),
    (...,  0 ),
    (-1 , ...),
)

borders = {}
graph = nx.Graph()
for tile_name, tile in tiles.items():
    for (i, j), f in product(edges, (1, -1)):
        border_id = "".join(tile[i, j][::f])
            
        if border_id in borders:
            graph.add_edge(borders[border_id], tile_name)
        else:
            borders[border_id] = tile_name
            
def orient(tile):
    for i in range(4):
        yield np.rot90(tile, i)
        yield np.fliplr(np.rot90(tile, i))

flip = {
    ...: ...,
    -1 :  0 ,
     0 : -1 ,
}

dire = {
    ...:  0 ,
    -1 :  1 ,
     0 : -1 ,
}

dim = int(sqrt(len(tiles)))

first_tile = next(iter(tiles))
result = np.full([2*dim+1, 2*dim+1, 10, 10], "")
result[dim][dim] = tiles[first_tile]
location = {first_tile: np.array([dim, dim])}

for left, right in nx.dfs_edges(graph, first_tile):
    if left in location:
        old, new = left, right
    else:
        old, new = right, left
        
    old_loc = location[old]
    old_tile = result[old_loc[0], old_loc[1]]
    for i, j in edges:
        for oriented_tile in orient(tiles[new]):
            if str(old_tile[i, j]) == str(oriented_tile[flip[i], flip[j]]):
                new_loc = old_loc + np.array([dire[i], dire[j]])
                result[new_loc[0], new_loc[1]] = oriented_tile
                location[new] = new_loc
                break

flat = np.swapaxes(result[:,:,1:-1,1:-1], 1, 2).reshape((2*dim+1)*8, -1)
flat_filtered  = flat[~(flat == '')].reshape(dim*8, -1)
image = (flat_filtered == "#").astype(int)

monster = (np.array(open("monster.txt").read().splitlines()).view('U1') == "#").reshape((3, -1)).astype(int)
monster_value = monster.sum()

total_monsters = sum((correlate(image, oriented_monster, mode="constant") == monster_value).sum() for oriented_monster in orient(monster))

print("remaining", image.sum() - total_monsters*monster_value)
```

Phew, this was so far the longest and most tedious challenge yet.
Onward!