# When it rains, it pours - August 2016 Community Challenge

### 1. Introduction

This code is my attempt at solving the August 2016 Community Challenge. Coming from a city where it rains cats and dogs on a daily basis this challenge was right up my alley =)

### 2. Algorithm

I used the solution 200_success ♦ outlined in his answer here:

1. Each Cell keeps track of which Basin it belongs to; each Cell is initially assume to be in its own Basin. Each Basin has a sink, or lowest Cell, which acts as a "representative element" of the Basin, as well as a member count. Topography keeps track of all Basins.
2. For each Basin, find lowest of the sink's neighbours. If the lowest is not already a member of this Basin, transfer its cells into the lowest neighbour's Basin, and notify Topography that the higher basin no longer exists.
3. Repeat step 2 until no further action is necessary.
4. Have Topography enumerate the Basins and their counts.

### 3. Input and output

Example 1

Input: rainfall-example-1.txt

Output:

Height Farmland:
[[1 5 2]
[2 4 7]
[3 6 9]]

Basins:
(A)  A  (B)
A   A   B
A   A   A

Letter  Size  Sink
A       7     (0, 0)
B       2     (0, 2)


Example 2

Input: rainfall-example-2.txt

Output:

Height Farmland:
[[1 0 2 5 8]
[2 3 4 7 9]
[3 5 7 8 9]
[1 2 5 4 3]
[3 3 5 2 1]]

Basins:
A   (A)  A  A   A
A    A   A  A   A
B    B   A  C   C
(B)   B   B  C   C
B    B   C  C  (C)

Letter  Size  Sink
A       11    (0, 1)
B        7    (3, 0)
C        7    (4, 4)


Example 3

Input: rainfall-example-3.txt

Output:

Height Farmland:
[[0 2 1 3]
[2 1 0 4]
[3 3 3 3]
[5 5 2 1]]

Basins:
(B)  B   A    A
B   A  (A)   A
B   A   A    C
B   C   C   (C)

Letter  Size  Sink
A       7     (1, 2)
B       5     (0, 0)
C       4     (3, 3)


Example 4

Input: rainfall-example-4.txt

Output:

Height Farmland:
[[ 4 23 25 21 29 16 23 29 12 28]
[ 0 12 26  0 19 23  9 13 11 29]
[26 24 18 21 22  4 29  1  5 28]
[13 15 18  3  6  7 15 15  0  9]
[29 29 23  6 28  1 11  1  3 21]
[ 6  3  0 13 11  0 28  0 25 17]
[20 15  7 24  3  8  5 21 12 23]
[ 0  9 24 12 19 23  9 29 26 21]
[ 1 12 12  2 14  2  0 16  2  6]
[14  5 14  7 26 12 24  6  5 25]
[18 25 20 29 17 23 23  2 24 19]
[ 9  0  6  2 19 19 12 10 18 28]
[ 8 27  7 23 14  9  3 14 18 25]
[ 6 19 13  9  3  0 21  3  2 16]
[ 6  1 14 12 19 22 15  2 19 12]
[17 24 27  8 15 26 16  6  0 27]
[ 0 15  3  4  2 19  0  3 17 19]
[ 3 17 14 19 20 20 25  1  7 19]
[10 13 13 22 27 20 21 28 12  4]
[27 20 19 17 28  0 13  4  1 10]]

Basins:
K     K     M     M     AI   (AI)   AE    A     A     A
(K)    K     M    (M)    M     Y    (AE)   AB    A     A
K     K    (AN)   M     Y    (Y)    AB   (AB)   A     A
(AJ)   AJ    X    (X)    X     C     C     A    (A)    A
B     B     B     X     C     C     C     T     A     A
B     B    (B)    B     C    (C)    C    (T)    T    (AR)
N     B     B     AA   (AA)   C    (AS)   T    (AK)   AK
(N)    N     B     L     AA    F     F     F     P     P
N     N     L    (L)    L     F    (F)    F    (P)    P
N    (AL)   AL    L     L     F     F     I     P     P
H     H     H     U    (AP)   F     I    (I)    I    (AO)
H    (H)    H    (U)    U     G     AF    I     I     I
AH    H     H     U     G     G    (AF)   V     W     W
(AH)   Q     H     G     G    (G)    G     V    (W)    W
Q    (Q)    Q     E     G     G     V    (V)    O    (AQ)
D     Q     AD    E     E     E     J     O    (O)    O
(D)    D    (AD)   E    (E)    J    (J)    J     O     O
D     D     AD    E     E     J     J    (Z)    Z     AG
D     D    (AC)   AC    E     R     R     Z     S    (AG)
D     D     AC   (AM)   R    (R)    R     S    (S)    S

Letter  Size  Sink
A      12    (3, 8)
B      10    (5, 2)
C       9    (5, 5)
D       9    (16, 0)
E       9    (16, 4)
F       9    (8, 6)
G       9    (13, 5)
H       9    (11, 1)
I       7    (10, 7)
J       6    (16, 6)
K       6    (1, 0)
L       6    (8, 3)
M       6    (1, 3)
N       6    (7, 0)
O       6    (15, 8)
P       6    (8, 8)
Q       5    (14, 1)
R       5    (19, 5)
S       4    (19, 8)
T       4    (5, 7)
U       4    (11, 3)
V       4    (14, 7)
W       4    (13, 8)
X       4    (3, 3)
Y       3    (2, 5)
Z       3    (17, 7)
AA       3    (6, 4)
AB       3    (2, 7)
AC       3    (18, 2)
AE       2    (1, 6)
AF       2    (12, 6)
AG       2    (18, 9)
AH       2    (13, 0)
AI       2    (0, 5)
AJ       2    (3, 0)
AK       2    (6, 8)
AL       2    (9, 1)
AM       1    (19, 3)
AN       1    (2, 2)
AO       1    (10, 9)
AP       1    (10, 4)
AQ       1    (14, 9)
AR       1    (5, 9)
AS       1    (6, 6)


Example 5

Input: rainfall-example-5.txt [20x20 map, height = 1000]

Example 6

Input: rainfall-example-6.txt [map: 55x55, height: 55^2]

• Not happy with using chararray as it seem it is deprecated. I tried using an array with the bool = string, but that threw me an error when I tried to update the array.
• The way I handle strings and str_rep feels wrong.
• The structure of my code feels right, but the classes feels very empty.
• My code struggles with big farmlands for an example for rainfall-example-6.txt. Is this normal, or could the algorithm be improved?
• Useless docstrings?

### 5. Code

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import string
import numpy as np
from numpy import unravel_index

ALPHABETH = string.ascii_uppercase
ALPHABETH_len = len(ALPHABETH)

def num_2_alpha(num):
'''
Converts an arabic number 0, 1, 2.. to it's corresponding letter A, B, C, ....
Example:
0 > A
1 > B
26 > Z
27 > AA
'''
quotient, remainder = divmod(num, ALPHABETH_len)
return quotient*ALPHABETH + ALPHABETH[remainder]

def create_test_file(max_height, shape):
'''
Creates a random height map with dimensions x, y (from shape)
and height from 0 to max_height.
'''
random_integers = np.random.randint(max_height, size=shape*shape)
return random_integers.reshape(shape, shape)

def format_topography(topography_):
'''
Inputs a typography of the farmland formated in a character array
this function formats the typography into a nicer looking string.
Input
[['(A)' 'A' '(B)']
['A' 'A' 'B']
['A' 'A' 'A']]
Output:
(A)  A  (B)
A   A   B
A   A   A
'''
rows, columns = topography_.shape

for i in range(columns):

for i in range(rows):
for j in range(columns):

def basin_2_string(basin_list):
'''
Input is a dictionary of basins where the key is the sink
Example output:
A: 7, Sink: (1, 2)
B: 5, Sink: (0, 0)
C: 4, Sink: (3, 3)
'''

basin_string = 'Letter  Size  Sink\n'
for i, basin in enumerate(basin_list):
basin_string += '''{:>{}} {} {:{}d} {} {}
return basin_string

def create_height_map(filename):
'''
Input: a filename with a textfile formated as
1 5 2
2 4 7
3 6 9
Ouputs a matrix of the height map.
'''
file = open(filename, 'r')
matrix = np.matrix([map(int, line.split()) for line in file])
file.close()
return matrix

def create_matrix_map(class_type, shape):
'''
Creates a dictionary where the keys are (x, y) coordinates.
This is used to store / index the basins and the cells
'''
length, width = shape
matrix = dict()
for i in xrange(length):
for j in xrange(width):
matrix[(i, j)] = class_type([i, j], [length, width])
return matrix

def neighboors_list(coordinates, shape):
'''
Makes a list of all the neighboors to a point in the height map
'''
length, width = shape
x, y = coordinates
neighboors = []

if x < 0 or x == length or y < 0 or y == width:
ValueError("The coordinates lies outside the matrix")
if x > 0:
neighboors.append((x-1, y))
if x + 1 < length:
neighboors.append((x+1, y))
if y > 0:
neighboors.append((x, y-1))
if y + 1 < width:
neighboors.append((x, y+1))
return neighboors

def create_basins_and_cells(height_map, shape):
'''
This creates the basins and the cells using the following algorithm:
1. Each Cell keeps track of which Basin it belongs to; each Cell is initially assume to be in its own Basin.
Each Basin has a sink, or lowest Cell, which acts as a "representative element" of the Basin,
as well as a member count. Topography keeps track of all Basins.
2. For each Basin, find lowest of the sink's neighbours. If the lowest is not already a member of this Basin,
transfer its cells into the lowest neighbour's Basin, and notify Topography that the higher basin no longer exists.
3. Repeat step 2 until no further action is necessary.
'''

basins = create_matrix_map(Basin, shape)
cells = create_matrix_map(Cell, shape)

topography_changed = True
while topography_changed:
topography_changed = False

for old_basin_coords in basins:
sink_coords = basins[old_basin_coords].sink
lowest_neighboor_height = height_map[sink_coords]
lowest_neighboor_coords = sink_coords

for neighboor in cells[sink_coords].neighboors:
if height_map[neighboor] < lowest_neighboor_height:
lowest_neighboor_height = height_map[neighboor]
lowest_neighboor_coords = neighboor

if lowest_neighboor_coords not in basins[old_basin_coords].cells:
topography_changed = True
new_basin_coords = cells[lowest_neighboor_coords].basin

for cell_coords in basins[old_basin_coords].cells:
basins[new_basin_coords].cells.append(tuple(cell_coords))
cells[cell_coords].basin = new_basin_coords

basins[new_basin_coords].size = len(
basins[new_basin_coords].cells)
del basins[old_basin_coords]
break

basin_list = sorted(basins.values(), key=lambda x: x.size, reverse=True)
return basin_list, cells

def create_topography(basin_list, shape):
'''
Enumerates the basins, in practice this creates the topography of the farmland
'''
character_length = len(num_2_alpha(len(basin_list)))
topography = np.chararray(shape, itemsize=character_length+2)
for i, basin in enumerate(basin_list):
letter = num_2_alpha(i)
for coords in basin.cells:
if coords == basin.sink:
topography[coords] = '('+letter+')'
else:
topography[coords] = letter

class Cell():

'''
Simple class describing a single cell in the height map.
'''

def __init__(self, coordinates, shape):
self.coordinates = tuple(coordinates)
self.neighboors = neighboors_list(coordinates, shape)
self.basin = tuple(coordinates)

class Basin():

'''
Simple class describing the basins in the height map.
'''

def __init__(self, coordinates, shape):
self.coordinates = tuple(coordinates)
self.sink = tuple(coordinates)
self.cells = [tuple(coordinates)]
self.size = 1

def __repr__(self):
return '''<Size: {}, Sink: {}, Cells: {}>\n
'''.format(self.size, self.sink, self.cells)

class Topography():

'''
A group of farmers has some elevation data, and we're going to help them understand
how rainfall flows over their farmland.

We'll represent the land as a two-dimensional array of altitudes and use the following
model, based on the idea that water flows downhill:

If a cell’s four neighboring cells all have higher altitudes, we call this cell a sink; water collects in sinks.

Otherwise, water will flow to the neighboring cell with the lowest altitude. If a cell is not a sink,
you may assume it has a unique lowest neighbor and that this neighbor will be lower than the cell.

Cells that drain into the same sink – directly or indirectly – are said to be part of the same basin.

Your challenge is to partition the map into basins. In particular, given a map of elevations,
your code should partition the map into basins and output the sizes of the basins, in descending order.
'''

def __init__(self, filename):
self.height_map = create_height_map(filename)
self.shape = self.height_map.shape
self.basins, self.cells = create_basins_and_cells(
self.height_map, self.shape)
self.topography = create_topography(self.basins, self.shape)
self.basin_sizes = [basin.size for basin in self.basins]

def __str__(self):
string = '''
Height Farmland:
{}

Basins:
{}
{}'''.format(
str(self.height_map),
format_topography(self.topography),
basin_2_string(self.basins))
return string

if __name__ == '__main__':

n = 50
seperator = '-'*n
print seperator
for i in range(1, 5):
print Topography('rainfall-example-{}.txt'.format(i))
print seperator
print
# np.savetxt('rainfall-example-4.txt', create_test_file(30, (20, 10)), delimiter = ' ', fmt='%d')
# np.savetxt('rainfall-example-5.txt', create_test_file(15, (30, 25)), delimiter = ' ', fmt='%d')


The classes feel empty because you are basically just using fancy dictionaries and a bunch of functions. Don't get me wrong - that's basically all a class is, but if you're going to commit to using classes then you really do want to wrap up functionality inside of them. Pretty much all of your functions should become instance methods of one of your classes.

I'm going to overview your code as written, from the entry point to the end.

The biggest problem that jumps out at me is that your create_test_file is

1. Poorly named - a better name is create_test_matrix
2. Broken. It doesn't guarantee that every cell will either be a sink or have a unique least neighbor.

There are three things I'd do to solve it.

1. Write a method called create_test_matrix that makes the matrix
2. Change create_test_file to actually write a valid test file
3. Make the matrix creation an iterative process - generate new values left-to-right, top-to-bottom (or whatever orientation you prefer). Make sure that they don't violate the constraint of the already created values, and move on.

Now that we're actually generating valid test files, let's look at the rest of your code.

It smells funny to me that your Topography class takes a filename. I'd expect a Topography to take data, and provide a classmethod that loads it from a file. This also means that with testing you don't need to go through the intermediate file. As an aside, your implementation is buggy (and your test files are invalid). They should have the size of the matrix as the first line, and your test files don't have that (and your code doesn't handle it).

This classmethod will end up basically being create_height_map so I'll pull the contents of that function out and put it right in there. While doing that, I'll also use a context manager to safely open and close the file. I'll also rename the file object variable so it doesn't mask the builtin file. I'll also use a list comprehension instead of map, as it's more Pythonic and easier to read. I forget if np.matrix can take a generator object, but if so then I'd use a generator comprehension instead of a list comprehension.

@classmethod
def from_file(cls, filename):
with open(filename, 'r') as rainfall_data:
data_iter = iter(rainfall_data)
# skip the shape
next(data_iter)
rainfall_matrix = np.matrix(
[[int(data) for data in line.split()]
for line in data_iter]
)


Next I'm going to pull create_basins_and_cells into Topography, and create_matrix_map as well.

I'll start with create_matrix_map. This can just be a dictionary comprehension. I'm also going to exchange the lists for tuples, as you should generally prefer immutability.

@staticmethod
def create_matrix_map(cls, shape):
length, width = shape
return {(i, j): cls((i, j), shape) for i in xrange(length) for j in xrange(width)}


I made it a static method because it doesn't require any intrinsic attribute of the class or instance. I still give it the cls first argument instead of class_type because that's more understandable for me - YMMV.

Next I'll dive into create_basins_and_cells. I'm not a huge fan of your algorithm - it seems that you should be able to iterate over the cells in a more intelligent manner. When I solved this myself I sorted the cells by descending height, and that lent itself to a very clean algorithm. I've also run into performance issues with my solution, however, so I'm going to propose something new here (in pseudocode because I haven't quite worked out the code yet).

The goal is to be able to iterate over all of the cells in an arbitrary order without having to do it way too many times. To do this I think we want to combine the "create basins and cells" stage with the "sort into the correct basins" stage. This should probably look something like this (pseudocode, I haven't completely figured out the actual code yet).

for row in topography
for cell in row
for each neighbor of cell
if the neighbor feeds into this cell and has already been loaded
add it to this cell's basin
if the cell is a sink
add it to an array of sinks
else

sizes = []
for sink in sinks
count the number of connected components to that sink, append it to sizes

print sizes sorted in descending order, space separated


This is still pretty rough, but I think we can speed things up by only making one pass through the data. We can also load the data into memory in stages, and if we store the results of these calculations on the other cells then I think we'll avoid having to do too much extra work per-cell (that's another thing you could do with your computation - a lot of your calculations are repeated and could be stored on the cells). I hope to come back to this answer this weekend and fix up this part and add the actual code.

I'll just touch on some of your string stuff to finish up. Overall, these should all become instance methods of your different classes. You should make sure each class has its own __str__, __repr__, and __unicode__ methods. Then when you finally print the topography you don't have to worry about having reasonable string representations of the basins and the cells, all you ahve to do is worry about padding.

Also, your best friend will become textwrap.dedent - you don't have to do that ugly unindent when you use multiline strings.

Minor automatic bookkeeping note: You do

basins[new_basin_coords].size = len(
basins[new_basin_coords].cells)


but since that relationship (A Basin's size is always the length of its .cells) is fixed, you should make it automatic in code on the Basin object:

@property
def size(self):
return len(self.cells)


Also, you're right that your classes are a bit bare. Try making them self-contained. The only thing that calls create_topography is Topography.__init__ so perhaps move it to be Topography.__init_topography. Similar with create_height_map and create_basins_and_cells and neighbors_list and basin_2_string and format_topography...

## Use a better algorithm

Since it hasn't been pointed out yet, the algorithm is faulty. The problem can easily be seen with simple input in which the entire region is one square basin:

Height Farmland:
[[2 2 2 2 2]
[2 1 1 1 2]
[2 1 0 1 2]
[2 1 1 1 2]
[2 2 2 2 2]]

Basins:
(H)   D    A    B   (I)
D   (D)   A   (B)   B
A    A   (A)   A    A
E   (E)   A   (C)   C
(F)   E    A    C   (G)

Letter  Size  Sink
A       9     (2, 2)
B       3     (1, 3)
C       3     (3, 3)
D       3     (1, 1)
E       3     (3, 1)
F       1     (4, 0)
G       1     (4, 4)
H       1     (0, 0)
I       1     (0, 4)


One way to fix that would be to consider, for each square, "which way would water run?" If we assume, as the problem states, that water always runs down the steepest gradient, it's simply a matter of creating a path from each point on the grid to wherever the water no longer flows.

The docstring for num_2_alpha says that 26 would be converted to Z and 27 would yield AA but that's clearly not the case. In fact, since numbering starts at 0, 25 would yield Z and 26 and 27 would yield AA and AB respectively.

## In Python 2.7, prefer xrange to range

It's a minor point, but the current code uses range several places in which xrange would likely offer better performance. See this question for a discussion of range and xrange.

This neat little challenge has two tricky aspects ('Knackpunkte', as they say in German):

• identification of the sink for a given cell based on a local neighbourhood of 4 cells
• identification of the connected components

Correctness of a solution hinges on getting these two things right - the rest is sleepwalking. Hence it is important make these two things separate focus points of the code, so that they can be coded (and verified/tested) separately.

It would be harmful if parts of the related logic where scattered high and wide over dozens of artifacts (functions, structures, classes) and hundreds of source code lines. For example, I still haven't understood some important details of the original Python solution, and although I suspect that it isn't correct I cannot say so for certain.

So here's a proposal for an alternative approach that should allow better separation of concerns and much greater simplicity.

## Disjoint-Set Union

The simplest and most efficient means of identifying the sets of connected components ('basins') is of course the Disjoint Set Forest, which is just a humble array with a tiny function. See Disjoint-set data structure in the Wikipedia, or the related HackerEarth tutorial.

The classic form of the function that finds the set representative ('root') for a given set member looks thus:

static int root (int[] a, int i)
{
for (int j = a[i]; j != i; )
a[i] = j = a[i = j];

return i;
}


While looking for the root, this function also performs a bit of path compression, to speed up further searches. This is not the full path compression that would make the algorithm provably optimal (as explained in the Wikipedia article) since this would require an auxiliary stack or a really humungous call stack. But it does improve things at negligible cost.

Here I'm proposing a variant where the set root slots are used to keep track of membership counts (stored negative), instead of pointing at themselves. This requires a slight modification of the root finding code: a lag of one step needs to be introduced since roots are no longer pointing at themselves and hence their slot contents can no longer be blindly stored to other slots. As a side effect, this modification eliminates the final redundant store of the classic version:

static int root (int[] a, int i)
{
int j = a[i];

if (j >= 0)
for (int k; (k = a[i = j]) >= 0; j = k)
a[i] = k;

return i;
}


In the classic form, two sets identified by arbitrary members x and y can be joined via a simple assignment:

a[root(x)] = root(y);


Having the membership counts available makes it possible to do Union by Rank, i.e. always appending the smaller set to the larger in order to keep things compact and to minimise the workload for the path compression. This ensures that the structure performs well under heavy loads (despite the not-quite optimal path compression) and so it is worth it to create a separate function for joining sets:

static void join (int[] a, int i, int j)
{
if ((i = root(a, i)) != (j = root(a, j)))
{
if (a[i] <= a[j])  // the roots hold negative membership counts
{
a[i] += a[j];
a[j] = i;
}
else
{
a[j] += a[i];
a[i] = j;
}
}
}


That's the connected component analysis done and dusted, with simplicity and (almost) perfect efficiency.

## Analysing the Local Neighbourhood of a Cell

Now on to the neighbourhood/sink analysis. The sink for each cell is determined by its immediate neighbours at the four points of the compass; the sink thus determined is final - there is no way for it to change, ever. Basin membership can vary over the course of processing but not what's the sink for any given cell. Hence it is possible to perform this analysis in a single sweep through the input, keeping only three rows in memory at a time (previous, current, and next).

An easy hack for simplifying matters is to pad the outer edge of the map with harmless dummy values, to reduce the number of cases that need to be handled separately (edge cases, literally). Since we are looking for minima, the padding cells are transparent if they have the maximum possible value for their data type. Encoding coordinates as row * column_count + column reduces clutter and maps the coordinates to a compact range of indices that can be used directly with the disjoint set forest.

The meat of the processing code is this loop which goes through the current row, finds the sink for each cell and updates the disjoint set forest dsf[] which represents basin membership:

for (int col = 1; col <= col_count; ++col)
{
var n_w = Math.Min(prev_row[col], curr_row[col - 1]);
var s_e = Math.Min(next_row[col], curr_row[col + 1]);
var min = Math.Min(n_w, s_e);

if (curr_row[col] > min)
{
int curr = (row - 1) * col_count + (col - 1), sink = curr;

if (min == prev_row[col])
sink -= col_count;
else if (min == next_row[col])
sink += col_count;
else if (min == curr_row[col - 1])
sink -= 1;
else
sink += 1;

join(dsf, curr, sink);
}
}


Column and row indices are base 1 because of the padding, hence the necessity for subtracting 1 when encoding the coordinates. Minimum-finding and determination of related coordinates is done separately to get clearer and much more compact code, which is necessary because C# lacks Python's elegance regarding things like tuple assignment etc.

The rest is elementary and not worth talking about. See CodeReview_Rainfall.cs on PasteBin, with conditional using directives so that it can be used unchanged with LINQPad (which can be downloaded for free). For completeness' sake, here's how to print the counts in descending order:

var basins = new HashSet<int>(Enumerable.Range(0, dsf.Length).Select(i => root(dsf, i)));
var counts = basins.Select(sink => -dsf[sink]).OrderBy(x => -x);

Console.WriteLine(string.Join(" ", counts));


Just to show that languages like Python and Perl no longer hold the monopoly on getting stuff done with a minimum of fuss. ;-)

• It seems to me that this might be better as a new C# question rather than a Python answer. – Edward Aug 21 '16 at 13:21