Stack Puzzle: Finding "patterns" in 2D arrays using Convolution Operations
Stack Puzzle: Finding "patterns" in 2D arrays using Convolution Operations
#stackoverflow #python #numpy #convolution #datascience
Around the month of January from 2021, I solved an interesting problem on Stack Overflow1 so in this post, allow me to share the problem with you as a coding puzzle, and let’s explore an interesting solution to it.
Puzzle: Given an imperfect 2D lattice/array with 1s and 0s forming a repeating pattern, find the “gaps” in the structure where some of the 1s are missing.
Here are a few examples of what the 2D lattices and the “gaps” in them look like -
The yellow points represent 1s and the rest represent 0s. Notice that in a few places where the regular structure should persist, a few 1s are missing. The goal of this puzzle is to write code that can identify these “gaps”. Here is the starter code to create these 2D arrays in python -
import numpy as np
arr = np.array([[1, 0, 0, 0, 1, 0, 0, 0, 0], #one gap here
[0, 0, 1, 0, 0, 0, 1, 0, 0],
[1, 0, 0, 0, 1, 0, 0, 0, 1],
[0, 0, 1, 0, 0, 0, 0, 0, 0], #one gap here
[1, 0, 0, 0, 1, 0, 0, 0, 1]])NOTE: There is one assumption here. There should exist at least one complete pattern in the lattice/array. Not every instance of the pattern can have a gap in it!
And, let’s say you are given a function to generate such “impure” lattices on demand, so you can test your algorithm properly!
import random
import numpy as np
import matplotlib.pyplot as plt
def generate_random_lattice(gaps=0.1):
#Generate random tile pattern
n,m = random.randrange(3,9,2), random.randrange(3,9,2)
tile = np.zeros((n,m)).astype(int)
tile[::tile.shape[0]-1, ::tile.shape[1]-1] = 1 #update corners
tile[tile.shape[0]//2, tile.shape[1]//2] = 1 #update center
tile = tile[:-1, :-1] #remove bottom edge
#Create pure lattice
x,y = np.random.randint(2,10), np.random.randint(2,10)
lattice = np.tile(tile, (x,y))
#Add impurities / gaps
ones_shape = lattice[lattice.nonzero()].shape
lattice[lattice.nonzero()] = np.random.binomial(n=1, p=1-gaps, size=ones_shape)
#Plot lattice
plt.imshow(lattice)
return lattice
arr = generate_random_lattice()
Do think about how you would solve this puzzle algorithmically, or better yet, attempt it before scrolling down to the (potential) solution!
There are a million ways to solve this problem, but I would like to explore one that involves an operation which is is very frequently used in the world of Data Science and Machine Learning: Convolutions.
Solution: Using Convolutions
The folks with a background in ML/Data Science would already know what the term “Convolution” is, as it is foundational to Computer Vision and Image Processing. Convolution Neural Networks (CNNs) are a class of neural network architectures and at their core, they use the operation named Convolution.
I will not go a lot deeper into what the Convolution operation or CNNs are but here is a quick overview.
What is a 2D Convolution?
A 2D convolution operation involves a matrix/image and a filter. The filter is moved/convolved over the matrix and a dot product is performed, which is returned in the output matrix. In an implementation, since these are independent operations, these can be done in parallel in a vectorized manner.
We can use other operations instead of a dot product as well! The idea here is that if the pattern in the filter is found in the image exactly, it should return a very high value. The analysis of the output of the operation reveals what the filter and the operation is designed to discover.
There is additional complexity here due to how padding can be leveraged etc, but I will skip that for this post. There is a lot of great material on this online, for example, this one. Here is a visual representation of what this operation looks like (image source: wikipedia)
An interesting thing here is that Convolutions are used for identifying “patterns” in an image, but they can also be used to find a lack of a given pattern, which is the case in this coding puzzle.
With this knowledge, let’s design an algorithm of how we can detect and mark these gaps given an “imperfect lattice” or an array with a repeating pattern with some imperfections.
Here is the pseudo-code -
Identify the smallest convolution window size by analyzing the proximity of 1s, row-wise and column-wise.
Estimate the number of such repeating patterns in the array over both axes.
Fetch all the convolution windows by iterating over the array.
Take the sum of each convolution window to search for “gaps”. The largest sum would indicate a complete pattern, while anything less would indicate incomplete patterns.
Fetch a copy of the complete pattern and subtract it (with broadcasting) from each convolution window to get the “gaps” as -1’s.
Check the location of the -1’s and mark them in the original array.
Let’s begin implementing it. First, let’s start by generating an example to work with.
arr = generate_random_lattice() #using the function above
print(arr.shape)
(19, 17)
Let’s start by analyzing the proximity of 1s over the 2 axes. When working on anything related to “distance between elements” or “span of a specific element” in arrays, I find using numpy.cumsum() is quite common. Because we have gaps in the lattice, we can get the pattern over the 2 axes by just taking a sum as below -
print(arr.sum(1))
print(arr.sum(0))
[4 0 0 3 0 0 5 0 0 3 0 0 5 0 0 4 0 0 4] #pattern of x-2-x-2-x
[4 0 2 0 4 0 3 0 4 0 3 0 4 0 2 0 2] #pattern of x-1-x-1-xWe can see that this shows us the structure of the pattern over both axes while “ignoring” the gaps in the lattice. Let’s use this to get the shape of the convolution window that we want to fetch from the array.
xdims = np.unique((arr.sum(1)!=0).cumsum(), return_counts=True)[1].max()*2+1
ydims = np.unique((arr.sum(0)!=0).cumsum(), return_counts=True)[1].max()*2+1
pattern_shape = (xdims, ydims)
print(pattern_shape)
(7, 5)With this, now we know that the smallest pattern shape that we can create convolution windows is (7,5). Now, we could create a perfect pattern of this shape and just use that to find imperfections in our lattice, but where is the fun in that? Let’s first find out how many sliding window operations are needed over both axes, for a pattern of this shape to repeat and fit into the overall array.
num_xpatterns = int(arr.shape[0]/(pattern_shape[0]-1))
num_ypatterns = int(arr.shape[1]/(pattern_shape[1]-1))
print((num_xpatterns, num_ypatterns))
(3, 4)This is straightforward. Since the pattern is not completely repeating, but shares an edge, we need to remove the last row/column of pixels from the pattern shape and then calculate how many times we need to move. Here, we need to repeat the (7,5) window 3 times to the right and 4 times downwards.
Let’s now break the original array into these 3 X 4 = 12 matrices of (7,5) shape sliding windows. I will use numpy.lib.stride_tricks2 for this.
#Calculating the stride and shape that needs to be taken with stride_tricks
shp = (num_xpatterns, num_ypatterns, xdims, ydims)
strd = (arr.strides[0]*(xdims-1), arr.strides[1]*(ydims-1), arr.strides[0], arr.strides[1])
#Generating rolling windows/convolutions over the image to separate the patterns.
convolve_pattern = np.lib.stride_tricks.as_strided(arr, shape=shp, strides=strd)
convolve_pattern.shape
(3, 4, 7, 5)As expected, we get a tensor of the shape (3,4,7,5) which means 3X4 matrices of shape (7,5). Let’s fetch one of them to check what it looks like.
convolve_pattern[0,-1,:,:]
array([[1, 0, 0, 0, 0], #<- one gap here
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[1, 0, 0, 0, 1]])This corresponds to this window of the array -
Great next, we just calculate the sum of each of these windows to find the gaps!
pattern_sums = convolve_pattern.sum(axis=(-1,-2))
pattern_sums
array([[4, 5, 5, 4],
[5, 5, 5, 4],
[5, 5, 5, 4]])So, in this (3,4) matrix, the values where the sum is 5 are examples of perfect patterns, while the ones < 5 are patterns with gaps. Let’s fetch a perfect pattern from this view.
pattern_sums = convolve_pattern.sum(axis=(-1,-2))
idx = np.unravel_index(np.argmax(pattern_sums), pattern_sums.shape)
truth_pattern = convolve_pattern[idx]
print(truth_pattern)
array([[1, 0, 0, 0, 1],
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[1, 0, 0, 0, 1]])Nice, we are almost there! Now, since we have this (7,5) pattern, we can subtract it from the convolved view of the array (3,4,7,5) with broadcasting to get the gap locations!
gaps = convolve_pattern - truth_pattern[None, None, :, :]
gaps[0,-1,:,:]
array([[ 0, 0, 0, 0, -1], #<- gap marked with -1
[ 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0]])Then we can use the beautiful concept of numpy views vs copies3, to overwrite the 0s of these gap locations on the memory of the original matrix with -1s. Overwriting the locations on the (3,4,7,5) view also overwrites the original location in the (19,17) array since they share the same memory! This is the power of understanding how to work with views vs copies!
for i in np.argwhere(gaps==-1):
convolve_pattern[tuple(i)]=-1Thats it! Here is the complete code and the function in action!
#IDENTIFYING IMPERFECT PATTERNS IN A LATTICE
import numpy as np
import matplotlib.pyplot as plt
def find_gaps(x, verbose=False):
#Identifying the size and shape of the repeating pattern
xdims = np.unique((x.sum(1)!=0).cumsum(), return_counts=True)[1].max()*2+1
ydims = np.unique((x.sum(0)!=0).cumsum(), return_counts=True)[1].max()*2+1
pattern_shape = (xdims, ydims)
#Calculating number of rolling windows that exist with that pattern
num_xpatterns = int(x.shape[0]/(pattern_shape[0]-1))
num_ypatterns = int(x.shape[1]/(pattern_shape[1]-1))
#Calculating the stride and shape that needs to be taken with stride_tricks
shp = (num_xpatterns, num_ypatterns, xdims, ydims)
strd = (x.strides[0]*(xdims-1), x.strides[1]*(ydims-1), x.strides[0], x.strides[1])
#Generating rolling windows/convolutions over the image to separate the patterns.
convolve_pattern = np.lib.stride_tricks.as_strided(x, shape=shp, strides=strd)
#Assuming at least 1 untouched pattern exists, finding that pure pattern
pattern_sums = convolve_pattern.sum(axis=(-1,-2))
idx = np.unravel_index(np.argmax(pattern_sums), pattern_sums.shape)
truth_pattern = convolve_pattern[idx]
#Printing Debugging info
if verbose==True:
print('x shape:',x.shape)
print('pattern shape:',pattern_shape)
print('convolved shape:',convolve_pattern.shape)
print('')
print('pattern sums')
print(pattern_sums)
print('')
print('true pattern')
print(truth_pattern)
#Identifying the gaps by subtracting the convolved image with the truth pattern
gaps = convolve_pattern - truth_pattern[None, None, :, :]
#Setting the gaps as -1 directly into the location of memory of the original image
for i in np.argwhere(gaps==-1):
convolve_pattern[tuple(i)]=-1
plt.imshow(x)
return x
arr = generate_random_lattice() #Generate lattice
o = find_gaps(x, verbose=False) #Find gaps in lattice
Thanks for reading all the way to the end! My goal here was to demonstrate the power of Convolution operations and Numpy in solving unique scenarios. Hope this was entertaining and educational for you, the reader!
https://stackoverflow.com/questions/65296608
https://numpy.org/doc/stable/reference/generated/numpy.lib.stride_tricks.as_strided.html
https://numpy.org/doc/stable/user/basics.copies.html