Partition a map of water flows

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.

Assume the elevation maps are square. Input will begin with a line with one integer, S, the height (and width) of the map. The next S lines will each contain a row of the map, each with S integers – the elevations of the S cells in the row. Some farmers have small land plots such as the examples below, while some have larger plots. However, in no case will a farmer have a plot of land larger than S = 5000.

Your code should output a space-separated list of the basin sizes, in descending order. (Trailing spaces are ignored.)

A few examples are below.

Input:

3
1 5 2
2 4 7
3 6 9


Output: 7 2

The basins, labeled with A’s and B’s, are:

A A B
A A B
A A A


Input:

1
10


Output: 1

There is only one basin in this case.

Input:

5
1 0 2 5 8
2 3 4 7 9
3 5 7 8 9
1 2 5 4 2
3 3 5 2 1


Output: 11 7 7

The basins, labeled with A’s, B’s, and C’s, are:

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


Input:

4
0 2 1 3
2 1 0 4
3 3 3 3
5 5 2 1


Output: 7 5 4

The basins, labeled with A’s, B’s, and C’s, are:

A A B B
A B B B
A B B C
A C C C

• I edited your question to make it more appropriate for this site. Earlier, it was a programming question / code review. Now it is in the form of the challenge. This site is for releasing code challenges / problems to the community for them to attempt. Note: you still require a winning criteria: shortest code (code-golf) is recommended. – Justin Jan 23 '14 at 7:15
• @OP If you want an answer for your original question instead of a number of alternative golfed solutions, I suggest asking it again on Stack Overflow (or maybe Code Review?) – Gareth Jan 23 '14 at 8:24
• @JanDvorak I think the original question before editing might be okay on Code Review (there wasn't any golfing involved to start with)? You're probably right about SO though. – Gareth Jan 23 '14 at 8:29
• @JanDvorak I think just edit it out and make it a valid code-golf – Justin Jan 23 '14 at 9:53
• I have posted the problem on code review - codereview.stackexchange.com/questions/39895/… – AnkitSablok Jan 23 '14 at 18:17

Mathematica

The basin size list can be gotten by

WatershedComponents[
CornerNeighbors -> False,
Method -> "Basins"
] // Reverse@Sort@Part[Tally[Flatten@#], All, 2] &


where m is the given input data. To display a matrix like the ones in the question one can replace // Reverse@Sort@Part[Tally[Flatten@#], All, 2] & with /. {1 -> "A", 2 -> "B", 3 -> "C"} // MatrixForm or one can display it as an image instead using //ImageAdjust//Image.

• Don't leave us hanging! The sorted basin size list would use BinCounts[] & Sort[], right? – Scott Leadley May 3 '14 at 15:27
• @ScottLeadley I didn't realize it was the basin size list that was requested, thank you for pointing that out. I've fixed the answer (although the last part can probably be made a lot shorter). – user11030 May 4 '14 at 14:42

JavaScript - 673 707 730 751

e=[],g=[],h=[],m=[],q=[];function r(){a=s,b=t;function d(d,A){n=a+d,p=b+A;c>e[n][p]&&(u=!1,v>e[n][p]&&(v=e[n][p],w=n,k=p))}c=e[a][b],u=!0,v=c,w=a,k=b;0!=a&&d(-1,0);a!=l&&d(1,0);0!=b&&d(0,-1);b!=l&&d(0,1);g[a][b]=w;h[a][b]=k;return u}function x(a,b,d){function c(a,b,c,k){g[a+b][c+k]==a&&h[a+b][c+k]==c&&(d=x(a+b,c+k,d))}d++;0!=a&&c(a,-1,b,0);a!=l&&c(a,1,b,0);0!=b&&c(a,0,b,-1);b!=l&&c(a,0,b,1);return d}y=$EXEC('cat "'+$ARG[0]+'"').split("\n");l=y[0]-1;for(z=-1;z++<l;)e[z]=y[z+1].split(" "),g[z]=[],h[z]=[];for(s=-1;s++<l;)for(t=-1;t++<l;)r()&&m.push([s,t]);for(z=m.length-1;0<=z;--z)s=m[z][0],t=m[z][1],q.push(x(s,t,0));print(q.sort(function(a,b){return b-a}).join(" "));

Test results (using Nashorn):

$for i in A B C D; do jjs -scripting minlm.js -- "test$i"; done
7 2
1
11 7 7
7 5 4
$ There would probably be stack problems for maps of size 5000 (but that's an implementation detail :). The unminified source in all it's fugliness: // lm.js - find the local minima // Globalization of variables. /* The map is a 2 dimensional array. Indices for the elements map as: [0,0] ... [0,n] ... [n,0] ... [n,n] Each element of the array is a structure. The structure for each element is: Item Purpose Range Comment ---- ------- ----- ------- h Height of cell integers s Is it a sink? boolean x X of downhill cell (0..maxIndex) if s is true, x&y point to self y Y of downhill cell (0..maxIndex) Debugging only: b Basin name ('A'..'A'+# of basins) Use a separate array-of-arrays for each structure item. The index range is 0..maxIndex. */ var height = []; var sink = []; var downhillX = []; var downhillY = []; //var basin = []; var maxIndex; // A list of sinks in the map. Each element is an array of [ x, y ], where // both x & y are in the range 0..maxIndex. var basinList = []; // An unordered list of basin sizes. var basinSize = []; // Functions. function isSink(x,y) { var myHeight = height[x][y]; var imaSink = true; var bestDownhillHeight = myHeight; var bestDownhillX = x; var bestDownhillY = y; /* Visit the neighbors. If this cell is the lowest, then it's the sink. If not, find the steepest downhill direction. This would be the place to test the assumption that "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." But right now, we'll take that on faith. */ function visit(deltaX,deltaY) { var neighborX = x+deltaX; var neighborY = y+deltaY; if (myHeight > height[neighborX][neighborY]) { imaSink = false; if (bestDownhillHeight > height[neighborX][neighborY]) { bestDownhillHeight = height[neighborX][neighborY]; bestDownhillX = neighborX; bestDownhillY = neighborY; } } } if (x !== 0) { // upwards neighbor exists visit(-1,0); } if (x !== maxIndex) { // downwards neighbor exists visit(1,0); } if (y !== 0) { // left-hand neighbor exists visit(0,-1); } if (y !== maxIndex) { // right-hand neighbor exists visit(0,1); } downhillX[x][y] = bestDownhillX; downhillY[x][y] = bestDownhillY; return imaSink; } function exploreBasin(x,y,currentSize) {//,basinName) { // This cell is in the basin. //basin[x][y] = basinName; currentSize++; /* Visit all neighbors that have this cell as the best downhill path and add them to the basin. */ function visit(x,deltaX,y,deltaY) { if ((downhillX[x+deltaX][y+deltaY] === x) && (downhillY[x+deltaX][y+deltaY] === y)) { currentSize = exploreBasin(x+deltaX,y+deltaY,currentSize); //,basinName); } return 0; } if (x !== 0) { // upwards neighbor exists visit(x,-1,y,0); } if (x !== maxIndex) { // downwards neighbor exists visit(x,1,y,0); } if (y !== 0) { // left-hand neighbor exists visit(x,0,y,-1); } if (y !== maxIndex) { // right-hand neighbor exists visit(x,0,y,1); } return currentSize; } // Read map from file (1st argument). var lines =$EXEC('cat "' + $ARG[0] + '"').split('\n'); maxIndex = lines.shift() - 1; for (var i = 0; i<=maxIndex; i++) { height[i] = lines.shift().split(' '); // Create all other 2D arrays. sink[i] = []; downhillX[i] = []; downhillY[i] = []; //basin[i] = []; } // Everyone decides if they are a sink. Create list of sinks (i.e. roots). for (var x=0; x<=maxIndex; x++) { for (var y=0; y<=maxIndex; y++) { if (sink[x][y] = isSink(x,y)) { // This node is a root (AKA sink). basinList.push([x,y]); } } } //for (var i = 0; i<=maxIndex; i++) { print(sink[i]); } // Each root explores it's basin. //var basinName = 'A'; for (var i=basinList.length-1; i>=0; --i) { // i-- makes Closure Compiler sad var x = basinList[i][0]; var y = basinList[i][1]; basinSize.push(exploreBasin(x,y,0)); //,basinName)); //basinName = String.fromCharCode(basinName.charCodeAt() + 1); } //for (var i = 0; i<=maxIndex; i++) { print(basin[i]); } // Done. print(basinSize.sort(function(a, b){return b-a}).join(' '));  I got better minimization results by breaking up the element objects into separate arrays, globalizing everywhere possible and embracing side-effects. NSFW. The effects of code minimization: • 4537 bytes, unminified • 1180 bytes, packer • 855 bytes, packer + hand optimizations (1 character global names) • 751 bytes, Google Closure Compiler with ADVANCED_OPTIMIZATIONS (NB, it elided a vestigial "return 0" as dead code) • 730 bytes, reckless hand optimization (I'm not changing the unminified source, so NSFW) • 707 bytes, more reckless hand optimization (remove all references to sink[]); • 673 bytes, remove all "var"s, drop Nashorn -strict flag I could have achieved close to 700 bytes without editing the minimized code if I'd been willing to modify the original source. But I didn't because I think leaving it as-is gives an interesting view from the starting point. • You can shorten var e=[],g=[],h=[],l,m=[],q=[] to e=g=h=l=m=q=[]. You can probably get rid of other uses of the var keyword too if you aren't shadowing any global variables. – nyuszika7h May 4 '14 at 17:22 • @nyuszika7h No can do. e=g=h=l=m=q=[] would have them all using a pointer to the same array. And Nashorn requires the var. – Scott Leadley May 4 '14 at 19:55 • @nyuszika7h You kicked my out of my rut. I dropped Nashorn -strict and deleted all the "var"s. – Scott Leadley May 4 '14 at 22:02 Python: 276 306 365 bytes This is my first golf attempt. Suggestions are appreciated! edit: imports and closing files take too many characters! So does storing files in variables and nested list comprehension. t=map(int,open('a').read().split());n=t.pop(0);q=n*n;r,b,u=range(q),[1]*q,1 while u!=0: u=0 for j in r: d=min((t[x],x)for x in [j,j-1,j+1,j-n,j+n]if int(abs(j/n-x/n))+abs(j%n-x%n)<=1 and x in r)[1] if j-d:u|=b[j];b[d]+=b[j];b[j]=0 for x in sorted(b)[::-1]:print x or '',  fully commented (2130 bytes...) from math import floor with open('a') as f: l = f.read() terrain = map(int,l.split()) # read in all the numbers into an array (treating the 2D array as flattened 1D) n = terrain.pop(0) # pop the first value: the size of the input valid_indices = range(n*n) # 0..(n*n)-1 are the valid indices of this grid water=[1]*(n*n) # start with 1 unit of water at each grid space. it will trickle down and sum in the basins. updates=1 # keep track of whether each iteration included an update # helper functions def dist(i,j): # returns the manhattan (L1) distance between two indices row_dist = abs(floor(j/n) - floor(i/n)) col_dist = abs(j % n - i % n) return row_dist + col_dist def neighbors(j): # returns j plus up to 4 valid neighbor indices possible = [j,j-1,j+1,j-n,j+n] # validity criteria: neighbor must be in valid_indices, and it must be one space away from j return [x for x in possible if dist(x,j)<=1 and x in valid_indices] def down(j): # returns j iff j is a sink, otherwise the minimum neighbor of j # (works by constructing tuples of (value, index) which are min'd # by their value, then the [1] at the end returns its index) return min((terrain[i],i) for i in neighbors(j))[1] while updates!=0: # break when there are no further updates updates=0 # reset the update count for this iteration for j in valid_indices: # for each grid space, shift its water d =down(j) if j!=d: # only do flow if j is not a sink updates += water[j] # count update (water[j] is zero for all non-sinks when the sinks are full!) water[d] += water[j] # move all of j's water into the next lowest spot water[j] = 0 # indicate that all water has flown out of j # at this point, water is zeros everywhere but the sinks. # the sinks have a value equal to the size of their watershed. # so, sorting water and printing nonzero answers gives us the result we want! water = sorted(water)[::-1] # [::-1] reverses the array (high to low) nonzero_water = [w for w in water if w] # 0 evaulates to false. print " ".join([str(w) for w in nonzero_water]) # format as a space-separated list  • Please don't golf a year. 365 chars is too nice. :P – tomsmeding May 5 '14 at 8:01 • I got it down to 306! I need those extra 59 days of vacation time. – wrongu May 5 '14 at 13:33 • You should be able to just do open('a').read(), I think. – MrLemon May 5 '14 at 14:45 JavaScript (ECMAScript 6) - 226 Characters s=S.split(/\s/);n=s.shift(k=[]);u=k.a;t=s.map((v,i)=>[v,i,1]);t.slice().sort(X=(a,b)=>a[0]-b[0]).reverse().map(v=>{i=v[1];p=[v,i%n?t[i-1]:u,t[i-n],(i+1)%n?t[i+1]:u,t[+n+i]].sort(X)[0];p==v?k.push(v[2]):p[2]+=v[2]});k.join(' ')  Explanation s=S.split(/\s/); // split S into an array using whitespace as the boundary. n=s.shift(); // remove the grid size from s and put it into n. k=[]; // an empty array to hold the position of the sinks. u=k.a; // An undefined variable t=s.map((v,i)=>[v,i,1]); // map s to an array of: // - the elevation // - the position of this grid square // - the number of grid squares which have flowed into // this grid square (initially 1). X=(a,b)=>a[0]-b[0]; // A comparator function for sorting. t.slice() // Take a copy of t .sort(X) // Then sort it by ascending elevation .reverse() // Reverse it to be sorted in descending order .map(v=>{ // For each grid square (starting with highest elevation) i=v[1]; // Get the position within the grid p=[v,i%n?t[i-1]:u,t[i-n],(i+1)%n?t[i+1]:u,t[+n+i]] // Create an array of the grid square and 4 adjacent // squares (or undefined if off the edge of the grid) .sort(X) // Then sort by ascending elevation [0]; // Then get the square with the lowest elevation. p==v // If the current grid square has the lowest elevation ?k.push(v[2]) // Then add the number of grid square which have // flowed into it to k :p[2]+=v[2]}); // Else flow the current grid square into its lowest // neighbour. k.join(' ') // Output the sizes of the block with space separation.  Previous Version - 286 Characters s=S.split(/\s/);n=s.shift()*1;k=[];u=k[1];t=s.map((v,i)=>({v:v,p:i,o:[]}));for(i in t){t[p=[t[i],i%n?t[i-1]:u,t[i-n],(+i+1)%n?t[+i+1]:u,t[+i+n]].sort((a,b)=>(a.v-b.v))[0].p].o.push([i]);p==i&&k.push([i])}k.map(x=>{while(x[L="length"]<(x=[].concat(...x.map(y=>t[y].o)))[L]);return x[L]})  Assumes that the input is in a variable S; Explanation s=S.split(/\s/); // split S into an array using whitespace as the boundary. n=s.shift()*1; // remove the grid size from s and put it into n. k=[]; // an empty array to hold the position of the sinks. u=k[1]; // Undefined t=s.map((v,i)=>({v:v,p:i,o:[]})); // map s to an Object with attributes: // - v: the elevation // - p: the position of this grid square // - o: an array of positions of neighbours which // flow into this grid square. for(i in t){ // for each grid square p=[t[i],i%n?t[i-1]:u,t[i-n],(+i+1)%n?t[+i+1]:u,t[+i+n]] // start with an array containing the objects // representing that grid square and its 4 neighbours // (or undefined for those neighbours which are // outside the grid) .sort((a,b)=>(a.v-b.v)) // then sort that array in ascending order of elevation [0].p // then get the first array element (with lowest // elevation) and get the position of that grid square. t[p].o.push([i]); // Add the position of the current grid square to the // array of neighbours which flow into the grid square // we've just found. p==i&&k.push([i]) // Finally, if the two positions are identical then // we've found a sink so add it to the array of sinks (k) } k.map(x=>{ // For each sink start with an array, x, containing the // position of the sink. while(x.length<(x=[].concat(...x.map(y=>t[y].o))).length); // Compare x to the concatenation of x with all the // positions of grid squares which flow into squares // in x and loop until it stops growing. return x.length // Then return the number of grid squares. })  Test S="3\n1 5 2\n2 4 7\n3 6 9"; s=S.split(/\s/);n=s.shift()*1;k=[];u=k[1];t=s.map((v,i)=>({v:v,p:i,o:[]}));for(i in t){t[p=[t[i],i%n?t[i-1]:u,t[i-n],(+i+1)%n?t[+i+1]:u,t[+i+n]].sort((a,b)=>(a.v-b.v))[0].p].o.push([i]);p==i&&k.push([i])}k.map(x=>{while(x[L="length"]<(x=[].concat(...x.map(y=>t[y].o)))[L]);return x[L]})  Outputs: [7, 2] S="5\n1 0 2 5 8\n2 3 4 7 9\n3 5 7 8 9\n1 2 5 4 2\n3 3 5 2 1" s=S.split(/\s/);n=s.shift()*1;k=[];u=k[1];t=s.map((v,i)=>({v:v,p:i,o:[]}));for(i in t){t[p=[t[i],i%n?t[i-1]:u,t[i-n],(+i+1)%n?t[+i+1]:u,t[+i+n]].sort((a,b)=>(a.v-b.v))[0].p].o.push([i]);p==i&&k.push([i])}k.map(x=>{while(x[L="length"]<(x=[].concat(...x.map(y=>t[y].o)))[L]);return x[L]})  Outputs: [11, 7, 7] S="4\n0 2 1 3\n2 1 0 4\n3 3 3 3\n5 5 2 1" s=S.split(/\s/);n=s.shift()*1;k=[];u=k[1];t=s.map((v,i)=>({v:v,p:i,o:[]}));for(i in t){t[p=[t[i],i%n?t[i-1]:u,t[i-n],(+i+1)%n?t[+i+1]:u,t[+i+n]].sort((a,b)=>(a.v-b.v))[0].p].o.push([i]);p==i&&k.push([i])}k.map(x=>{while(x[L="length"]<(x=[].concat(...x.map(y=>t[y].o)))[L]);return x[L]})  Outputs: [5, 7, 4] • To my eye, the arrow (=>) function definitions are much clearer. – Scott Leadley May 9 '14 at 18:49 Julia, 315 function f(a,i,j) z=size(a,1) n=filter((x)->0<x[1]<=z&&0<x[2]<=z,[(i+1,j),(i-1,j),(i,j-1),(i,j+1)]) v=[a[b...] for b in n] all(v.>a[i,j]) && (return i,j) f(a,n[indmin(v)]...) end p(a)=prod(["$n " for n=(b=[f(a,i,j) for i=1:size(a,1),j=1:size(a,2)];sort([sum(b.==s) for s=unique(b)],rev=true))])


Just a recursive function that either determines the current cell is a sink or finds the drain, then call that on every set of indices. Didn't bother to do the input part since I wasn't going to win anyway, and that part isn't fun.

import Data.List
m=map
q[i,j]=[-1..1]>>= \d->[[i+d,j],[i,j+d]]
x%z=m(\i->snd.fst.minimum.filter((elemq i).snd)$zip(zip z[0..])x)x g(n:z)=iterate(\v->m(v!!)v)(sequence[[1..n],[1..n]]%z)!!(n*n) main=interact$unwords.m show.reverse.sort.m length.group.sort.g.m read.words


Might be still some code to be golf'd here.

& runhaskell 19188-Partition.hs <<INPUT
> 5
> 1 0 2 5 8
> 2 3 4 7 9
> 3 5 7 8 9
> 1 2 5 4 2
> 3 3 5 2 1
INPUT
11 7 7


Explanation

Basic idea: For each cell (i,j) find the lowest cell in the "neighborhood". This gives a graph [(i,j)(mi,mj)]. If a cell is the lowest cell itself, then (i,j) == (mi,mj).

This graph can be iterated: For each a → b in the graph, replace it with a → c where b → c is in the graph. When this iteration yields no more changes, then each cell in the graph points at the lowest cell it will flow to.

To golf this, several changes have been made: First, coordinates are represented as a list of length 2, rather than a pair. Second, once the neighbors have been found, cells are represented by their index into a linear array of the cells, not 2D coordinates. Third, as there are n*n cells, after n*n iterations, the graph must be stable.

Ungolf'd

type Altitude = Int     -- altitude of a cell

type Coord = Int        -- single axis coordinate: 1..n
type Coords = [Coord]   -- 2D location, a pair of Coord
-- (Int,Int) would be much more natural, but Coords are syntehsized
-- later using sequence, which produces lists

type Index = Int        -- cell index
type Graph = [Index]    -- for each cell, the index of a lower cell it flows to

neighborhood :: Coords -> [Coords]                              -- golf'd as q
neighborhood [i,j] = concatMap (\d -> [[i+d,j], [i,j+d]]) [-1..1]
-- computes [i-1,j] [i,j-1] [i,j] [i+1,j] [i,j+1]
-- [i,j] is returned twice, but that won't matter for our purposes

flowsTo :: [Coords] -> [Altitude] -> Graph                      -- golf'd as (%)
flowsTo cs vs = map lowIndex cs
where
lowIndex is = snd . fst                          -- take just the Index of
. minimum                          -- the lowest of
. filter (inNeighborhood is . snd) -- those with coords nearby
$gv -- from the data inNeighborhood :: Coords -> Coords -> Bool inNeighborhood is ds = ds elem neighborhood is gv :: [((Altitude, Index), Coords)] -- the altitudes paired with their index and coordinates gv = zip (zip vs [0..]) cs flowInput :: [Int] -> Graph -- golf'd as g flowInput (size:vs) = iterate step (flowsTo coords vs) !! (size * size) where coords = sequence [[1..size],[1..size]] -- generates [1,1], [1,2] ... [size,size] step :: Graph -> Graph step v = map (v!!) v -- follow each arc one step main' :: IO () main' = interact$
unwords . map show      -- counts a single line of text
. reverse . sort        -- counts from hi to lo
. map length            -- for each common group, get the count
. group . sort          -- order cells by common final cell index
. flowInput             -- compute the final cell index graph
. map read . words      -- all input as a list of Int

• It'd be great if you could explain what's going on here. – Not that Charles May 8 '14 at 15:58
• @Charles - done! – MtnViewMark May 12 '14 at 14:42

Ruby, 216

r=[]
M=gets('').split.map &:to_i
N=M.shift
g=M.map{1}
M.sort.reverse.map{|w|t=[c=M.index(w),c%N<0?c:c-1,c%N<N-1?c+1:c,c+N,c-N].min_by{|y|M[y]&&y>=0?M[y]:M.max}
M[c]+=1
t!=c ?g[t]+=g[c]:r<<g[c]}
$><<r.sort.reverse*' '  It's a slightly different approach, only invoking "flow" on each square once (performance depends what the performance of Array::index is). It goes from the highest elevation to the lowest, emptying out one cell at a time into its lowest neighbor and marking the cell done (by adding 1 to the elevation) when it's done. Commented and spaced: results=[] ELEVATIONS = gets('').split.map &:to_i # ELEVATIONS is the input map MAP_SIZE = ELEVATIONS.shift # MAP_SIZE is the first line of input watershed_size = ELEVATIONS.map{1} # watershed_size is the size of the watershed of each cell ELEVATIONS.sort.reverse.map { |water_level| # target_index is where the water flows to. It's the minimum elevation of the (up to) 5 cells: target_index = [ current_index = ELEVATIONS.index(water_level), # this cell (current_index % MAP_SIZE) < 0 ? current_index : current_index-1, # left if possible (current_index % MAP_SIZE) >= MAP_SIZE-1 ? current_index : current_index+1, # right if possible current_index + MAP_SIZE, # below current_index - MAP_SIZE # above ].min_by{ |y| # if y is out of range, use max. Else, use ELEVATIONS[y] (ELEVATIONS[y] && y>=0) ? ELEVATIONS[y] : ELEVATIONS.max } # done with this cell. # increment the elevation to mark done since it no longer matters ELEVATIONS[current_index] += 1 # if this is not a sink (target_index != current_index) ? # add my watershed size to the target's watershed_size[target_index] += watershed_size[current_index] # else, push my watershed size onto results : results << watershed_size[current_index]}  Changelog: 216 - better way to deselect out-of-bounds indices 221 - turns out, "11" comes before "2"... revert to to_i, but save some space on our getses. 224 - Why declare s, anyway? And each => map 229 - massive golfing - sort the elevations first into s (and thereby drop the while clause), use min_by instead of sort_by{...}[0], don't bother to_i for elevations, use flat_map, and shrink select{} block 271 - moved watershed size into new array and used sort_by 315 - moved results to array which gave all sorts of benefits, and shortened neighbor index listing. also gained one char in index lambda. 355 - first commit Python - 470 447 445 393 392 378 376 375 374 369 bytes I can't stop myself! Not a winning solution, but I had a lot of fun creating it. This version does not assume the input to be stored anywhere and instead reads it from stdin. Maximum recursion depth = longest distance from a point to it's sink. def f(x,m=[],d=[],s=[]): n=[e[a]if b else 99for a,b in(x-1,x%z),(x+1,x%z<z-1),(x-z,x/z),(x+z,x/z<z-1)];t=min(n) if t<e[x]:r=f(x+(-1,1,-z,z)[n.index(t)])[0];s[r]+=x not in m;m+=[x] else:c=x not in d;d+=[x]*c;r=d.index(x);s+=[1]*c return r,s z,e=input(),[] exec'e+=map(int,raw_input().split());'*z for x in range(z*z):s=f(x)[1] print' '.join(map(str,sorted(s)[::-1]))  I don't have time to explain it today, but here's the ungolfed code: It actually is quite different than the original code. I read S lines from the stdin, split, map to ints and flatten the lists to get the flattened field. Then I loop through all tiles (let me call them tiles) once. The flow-function checks the neighboring tiles and picks the one with the smallest value. If it's smaller than value of the current tile, move to it and recurse. If not, the current tile is a sink and new basin is created. The return value of the recursion is the id of the basin. # --- ORIGINAL SOURCE --- # lowest neighboring cell = unique and next # neihboring cells all higher = sink and end basinm = [] # list of the used tiles basins = {} # list of basin sizes basinf = [] # tuples of basin sinks field = [] # 2d-list representing the elevation map size = 0 def flow(x, y): global basinf, basinm print "Coordinate: ", x, y nearby = [] nearby += [field[y][x-1] if x > 0 else 99] nearby += [field[y][x+1] if x < size-1 else 99] nearby += [field[y-1][x] if y > 0 else 99] nearby += [field[y+1][x] if y < size-1 else 99] print nearby next = min(nearby) if next < field[y][x]: i = nearby.index(next) r = flow(x+(-1,1,0,0)[i], y+(0,0,-1,1)[i]) if (x,y) not in basinm: basins[r] += 1 basinm += [(x,y)] else: c = (x,y) not in basinf if c: basinf += [(x,y)] r = basinf.index((x,y)) if c: basins[r] = 1 return r size = input() field = [map(int,raw_input().split()) for _ in range(size)] print field for y in range(size): for x in range(size): flow(x, y) print print ' '.join(map(str,sorted(basins.values(),reverse=1)))  JavaScript (ES6) 190 203 Edit A little more ES6ish (1 year later...) Define a function with input rows as a string, including newlines, return output as string with trailing blanks F=l=>{[s,...m]=l.split(/\s+/);for(j=t=[];k=j<s*s;t[i]=-~t[i])for(i=j++;k;i+=k)k=r=0,[for(z of[-s,+s,i%s?-1:+s,(i+1)%s?1:+s])(q=m[z+i]-m[i])<r&&(k=z,r=q)];return t.sort((a,b)=>b-a).join(' ')} // Less golfed U=l=>{ [s,...m] = l.split(/\s+/); for (j=t=[]; k=j<s*s; t[i]=-~t[i]) for(i=j++; k; i+=k) k=r=0, [for(z of [-s,+s,i%s?-1:+s,(i+1)%s?1:+s]) (q=m[z+i]-m[i]) < r && (k=z,r=q)]; return t.sort((a,b)=>b-a).join(' ') } // TEST out=x=>O.innerHTML += x + '\n'; out(F('5\n1 0 2 5 8\n 2 3 4 7 9\n 3 5 7 8 9\n 1 2 5 4 2\n 3 3 5 2 1'))// "11 7 7" out(F('4\n0 2 1 3\n2 1 0 4\n3 3 3 3\n5 5 2 1')) //"7 5 4" <pre id=O></pre> Perl 6, 419 404 Newlines added for clarity. You can safely remove them. my \d=$*IN.lines[0];my @a=$*IN.lines.map(*.trim.split(" "));my @b;my$i=0;my $j=0; for @a {for @$_ {my $c=$_;my $p=$i;my $q=$j;my &y={@a[$p+$_[0]][$q+$_[1]]//Inf};
loop {my @n=(0,1),(1,0);push @n,(-1,0) if $p;push @n,(0,-1) if$q;my \o=@n.sort(
&y)[0];my \h=y(o);last if h>$c;$c=h;$p+=o[0];$q+=o[1]};@b[$i][$j]=($p,$q);++$j};$j=0;++$i};say join " ",bag(@b.map(*.flat).flat.map(~*)).values.sort: {$^b <=>$^a}  Old solution: my \d=$*IN.lines[0];my @a=$*IN.lines.map(*.trim.split(" "));my @b;my$i=0;my $j=0; for @a {for @$_ {
my $c=$_;my $p=$i;my $q=$j;
loop {my @n=(0,1),(1,0);@n.push: (-1,0) if $p;@n.push: (0,-1) if$q;
my \o=@n.sort({@a[$p+$_[0]][$q+$_[1]]//Inf})[0];
my \h=@a[$p+o[0]][$q+o[1]];last if h>$c;$c=h;$p+=o[0];$q+=o[1]};@b[$i][$j]=($p,$q);++$j};$j=0;++$i}; say join " ",bag(@b.map(*.flat.flat).flat.map(~*)).values.sort: {$^b <=>\$^a}


And yet I get beaten by Python and JavaScript solutions.