# Introduction

For many centuries, there has been a certain river that has never been mapped. The Guild of Cartographers want to produce a map of the river, however, they have never managed to succeed -- for some reason, all the cartographers they have sent to map the river have been eaten by wild animals in the area. A different approach is required.

# Input Description

The area is a rectangular grid of cells of length m and width n. The cell in the bottom left would be 0,0, and the cell in the top right would be m-1,n-1. m and n are provided in the input as a tuple m,n.

By using long distance geographical sounding techniques the location of islands around the river have been identified. The size of the islands (i.e. how many cells the island occupies) have also been identified but the shape has not. We supply this information in a tuple s,x,y where s is the size of the island, and x and y are the x and y positions of one particular cell of that island. Each tuple in the input is space separated, so here is an example input:

7,7 2,0,0 2,3,1 2,6,1 2,4,3 2,2,4 8,0,6 1,2,6 3,4,6


To illustrate more clearly, here are is the input on a graph:

 y 6|8 1 3
5|
4|  2
3|    2
2|
1|   2  2
0|2
=======
0123456
x


# Output Description

Output a map using ASCII characters to represent parts of the area. Each cell will either be # (land) or . (water). The map should follow these rules:

1. Definition. An island is a orthogonally contiguous group of land cells that is bounded entirely by river cells and/or the border of the area.
2. Definition. A river is an orthogonally contiguous group of water cells that is bounded entirely by land cells and/or the border of the area, and does not contain "lakes" (2x2 areas of water cells).
3. Rule. The map shall contain exactly one river.
4. Rule. Each numbered cell in the input shall be part of an island containing exactly s cells.
5. Rule. Every island in the map shall contain exactly one of the numbered cells in the input.
6. Rule. There exists a single unique map for every input.

Here is the output of the example input:

#.#.##.
#....#.
#.##...
##..##.
###....
...##.#
##....#


Here is another input and output.

Input:

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


Output:

#.#.#
#.#.#
.....
###.#
.....

• Note: this is the same question as a Nurikabe solver. – absinthe May 30 '15 at 4:02
• Can we take input in any convenient format, or should we stick to the one in the question? – Erik the Outgolfer Jun 28 '18 at 18:27
• this is also problem 4 from the 2012 Dyalog competition – ngn Jun 29 '18 at 11:46
• @ngn Since when is "post a cryptographic hash"... usual? (but I suppose it's allowed when a challenge explicitly allow it) – user202729 Jul 6 '18 at 10:21
• here's a bookmarklet for puzzle-nurikabe.com - it converts the current puzzle to a valid input for this challenge and shows it in red just below the board: javascript:(_=>{var t=Game.nurikabe().task,m=t.length,n=t[0].length,s=[m,n];for(var i=0;i<m;i++)for(var j=0;j<n;j++)if(t[i][j]>=0)s+=' '+[t[i][j],i,j];puzzleContainerDiv.insertAdjacentHTML('beforeend','<hr><tt style=color:red>'+s+'</tt><hr>')})();void(0) – ngn Jul 12 '18 at 8:25

# C + PicoSAT, 2345995 952 bytes

#include<picosat.h>
#define f(i,a)for(i=a;i;i--)
#define b calloc(z+1,sizeof z)
#define e(a,q)if(a)A[q]^A[p]?l[q]++||(j[++k]=q):s[q]||(i[q]=p,u(q));
z,F,v,k,n,h,p,q,r,C,*x,*A,*i,*l,*s,*j,*m;u(p){s[m[++n]=p]=1;e(p%F-1,p-1)e(p%F,p+1)e(p>F,p-F)e(p<=F*v-F,p+F)}t(){f(q,k)l[j[q]]=0;f(q,n)s[m[q]]=0;k=n=0;i[p]=-1;u(p);}main(){void*E=picosat_init();if(scanf("%d,%d",&F,&v)-2)abort();z=F*v;for(x=b;scanf("%d,%d,%d",&r,&p,&q)==3;g(p),g(0))x[p=F-p+q*F]=r;f(p,F*v-F)if(p%F)g(p),g(p+1),g(p+F),g(p+F+1),g(0);for(A=b,i=b,l=b,s=b,j=b,m=b;!C;){picosat_sat(E,C=h=-1);f(p,F*v)A[p]=picosat_deref(E,p)>0,i[p]=0;f(p,F*v)if(x[p])if(i[q=p]){for(g(-q);i[q]+1;)q=i[q],g(-q);g(C=0);}else if(t(),r=n-x[p]){f(q,r<0?k:n)g(r<0?j[q]:-m[q]);g(C=0);}f(p,F*v)if(!i[p])if(t(),A[p]){g(-++z);f(q,k)g(j[q]);g(C=0);f(q,n)g(-m[q]),g(z),g(0);}else{C&=h++;f(q,k)g(-j[q]);g(++z);g(++z);g(0);f(q,F*v)g(s[q]-z),g(q),g(0);}}f(p,F*v)putchar(A[p]?35:46),p%F-1||puts("");}


Try it online!

(Warning: this TIO link is a 30 kilobyte URL that contains a minified copy of PicoSAT 965, so you might not be able to load it in some browsers, but it loads in at least Firefox and Chrome.)

### How it works

We initialize the SAT solver with a variable for each cell (land or water), and only the following constraints:

1. Every numbered cell is land.
2. Every 2×2 rectangle has at least one land.

The rest of the constraints are difficult to encode directly into SAT, so instead we run the solver to get a model, run a sequence of depth-first searches to find the connected regions of this model, and add additional constraints as follows:

1. For every numbered cell in a land region that’s too big, add a constraint that there should be at least one water cell among the current land cells in that region.
2. For every numbered cell in a land region that’s too small, add a constraint that there should be at least one land cell among the current water cells bordering that region.
3. For every numbered cell in the same land region as another numbered cell, add a constraint that there should be at least one water cell along the path of current land cells between them (found by walking the parent pointers left over from the depth-first search).
4. For every land region including no numbered cells, add constraints that either
• all of those current land cells should be water, or
• at least one of the current water cells bordering that region should be land.
5. For every water region, add constraints that either
• all of those current water cells should be land, or
• every cell other than those current water cells should be land, or
• at least one of the current land cells bordering that region should be water.

Taking advantage of the incremental interface to the PicoSAT library, we can immediately rerun the solver including the added constraints, preserving all the previous inferences made by the solver. PicoSAT gives us a new model, and we continue iterating the above steps until the solution is valid.

This is remarkably effective; it solves 15×15 and 20×20 instances in a tiny fraction of a second.

(Thanks to Lopsy for suggesting this idea of interactively constraining connected regions in an incremental SAT solver, a while back.)

### Some larger test cases from puzzle-nurikabe.com

A random page of 15×15 hard puzzles (5057541, 5122197, 5383030, 6275294, 6646970, 6944232):

15,15 1,5,1 3,9,1 5,4,2 1,6,2 2,11,2 2,2,3 3,9,3 2,4,4 1,10,4 5,12,4 3,1,5 1,3,5 3,8,5 1,13,5 5,5,6 1,12,6 1,2,8 2,9,8 1,1,9 2,6,9 6,11,9 3,13,9 5,2,10 2,4,10 4,10,10 1,5,11 2,12,11 2,3,12 2,8,12 5,10,12 1,5,13 1,9,13 1,6,14 1,8,14
15,15 4,2,0 2,5,0 1,3,1 2,14,2 1,3,3 2,11,3 1,13,3 1,5,4 11,7,4 1,9,4 1,4,5 1,8,5 2,10,5 12,14,5 3,5,6 1,4,7 2,10,7 3,9,8 4,0,9 1,4,9 1,6,9 3,10,9 1,5,10 1,7,10 8,9,10 1,1,11 10,3,11 2,11,11 6,0,12 1,11,13 2,9,14 1,12,14
15,15 2,2,0 8,10,0 2,3,1 2,14,2 2,3,3 3,5,3 3,9,3 2,11,3 5,13,3 6,0,4 3,7,4 3,3,5 2,11,5 2,6,6 1,8,6 1,4,7 2,10,7 1,6,8 2,8,8 5,3,9 2,11,9 2,7,10 7,14,10 2,1,11 4,3,11 2,5,11 1,9,11 2,11,11 2,0,12 4,6,13 1,11,13 3,4,14 1,12,14
15,15 2,0,0 2,4,0 3,6,1 2,10,1 1,13,1 2,5,2 2,12,2 3,0,3 2,2,3 4,7,3 2,9,3 1,14,3 1,4,4 1,8,4 2,12,5 4,2,6 3,4,6 1,14,6 7,7,7 1,10,8 2,12,8 3,2,9 2,14,9 2,0,10 2,6,10 1,10,10 2,5,11 4,7,11 2,12,11 1,14,11 3,2,12 3,9,12 1,1,13 2,4,13 3,8,13 2,10,14 5,14,14
15,15 1,3,0 1,14,0 3,7,1 3,10,1 2,13,1 3,1,2 4,5,2 2,12,3 3,3,4 1,8,4 1,1,5 3,5,5 1,9,5 5,13,5 3,3,6 1,8,6 2,2,7 2,12,7 1,6,8 1,8,8 2,11,8 2,1,9 4,5,9 2,9,9 2,13,9 2,6,10 4,11,10 1,2,11 3,9,12 2,13,12 3,1,13 2,4,13 3,7,13 1,0,14
15,15 2,8,0 2,4,1 2,7,1 1,10,1 6,4,3 1,1,4 12,5,4 3,11,4 5,13,4 3,10,5 3,0,6 1,6,6 2,8,6 4,13,7 2,3,8 1,6,8 3,8,8 2,14,8 2,4,9 5,1,10 4,3,10 1,9,10 6,13,10 3,8,11 1,10,11 3,4,13 2,7,13 3,10,13 1,6,14 1,14,14


A random page of 20×20 normal puzzles (536628, 3757659):

20,20 1,0,0 3,2,0 2,6,0 1,13,0 3,9,1 3,15,1 2,7,2 3,13,2 3,0,3 2,3,3 3,18,3 3,5,4 2,9,4 2,11,4 2,16,4 1,0,5 2,7,5 1,10,5 1,19,5 3,2,6 1,11,6 2,17,6 2,0,7 3,4,7 5,6,7 2,9,7 4,13,7 3,15,7 1,3,8 1,10,8 1,14,9 2,18,9 3,1,10 2,4,10 1,8,10 1,10,10 3,12,10 3,16,10 1,9,11 1,17,11 2,19,11 2,0,12 2,2,12 1,4,12 4,6,12 2,13,12 2,15,12 1,14,13 2,17,13 1,3,14 2,5,14 4,7,14 2,15,14 3,0,15 1,2,15 2,13,15 3,18,15 3,7,16 7,10,16 1,17,16 2,0,17 2,3,17 2,5,17 3,11,17 3,15,17 1,0,19 1,2,19 1,4,19 2,6,19 5,8,19 1,11,19 1,13,19 3,15,19 2,18,19
20,20 1,0,0 1,4,0 5,8,0 1,17,0 1,19,0 2,17,2 3,6,3 2,10,3 2,12,3 4,14,3 6,0,4 3,4,4 4,7,4 1,11,4 1,18,4 1,6,5 3,12,5 4,15,5 4,4,6 2,16,6 2,19,6 6,0,7 3,10,7 2,12,8 2,17,8 3,3,9 2,5,9 4,8,9 2,10,9 3,0,10 1,2,10 5,14,10 2,16,10 2,19,10 7,7,11 3,12,12 2,17,12 2,2,13 4,4,13 3,6,13 4,14,13 3,0,14 1,3,14 1,5,14 3,16,14 1,2,15 1,9,15 2,11,15 5,13,15 3,19,15 1,4,16 3,6,16 1,3,17 1,12,17 1,14,17 1,16,17 6,0,19 2,2,19 3,5,19 2,7,19 5,9,19 1,11,19 2,13,19 1,15,19 4,17,19


# GLPK, (non-competing) 2197 bytes

This is a non-competing entry, as

• I don't follow the input data format (as GLPK can only read input data from files) and
• GLPK seems not to be available on RIO.

I'll save a still ungolfed, yet functional version here. Depending on feedback, I might try to shorten it + add an explanation if there's interest. So far, the constraint names serve as in-place docs.

The main idea is to encode the "contiguous islands" constraint by introducing a preserved flow variable that has a pre-specified source at the hint location. The decision variable is_island then plays the role of placeable sinks. By minimizing the total sum of this flow, the islands are forced to remain connected in the optimum. The other constraints rather directly implement the various rules. What puzzles me that I still seem to need island_fields_have_at_least_one_neighbor.

Performance of this formulation is not great, though. By directly eating all the complexity with a large amount of constraints, the solver takes close to 15 seconds for the 15 x 15 example.

## Code (still ungolfed)

# SETS
param M > 0, integer; # length
param N > 0, integer; # width
param P > 0, integer; # no. of islands

set X := 0..N-1;  # set of x coords
set Y := 0..M-1;  # set of y coords
set Z := 0..P-1;  # set of islands

set V := X cross Y;
set E within V cross V := setof{
(sx, sy) in V, (tx, ty) in V :

((abs(sx - tx) = 1) and (sy = ty)) or
((sx = tx) and (abs(sy - ty) = 1))
}
(sx, sy, tx, ty);

# PARAMETERS
param islands{x in X, y in Y, z in Z}, integer, default 0;
param island_area{z in Z} := sum{x in X, y in Y} islands[x, y, z];

# VARIABLES
var is_river{x in X, y in Y}, binary;
var is_island{x in X, y in Y, z in Z}, binary;
var flow{(sx, sy, tx, ty) in E, z in Z} >= 0;

# OBJECTIVE
minimize obj: sum{(sx, sy, tx, ty) in E, z in Z} flow[sx, sy, tx, ty, z];

# CONSTRAINTS
s.t. islands_are_connected{(x, y) in V, z in Z}:

islands[x, y, z]
- is_island[x, y, z]
+ sum{(sx, sy, tx, ty) in E: x = tx and y = ty} flow[sx, sy, x, y, z]
- sum{(sx, sy, tx, ty) in E: x = sx and y = sy} flow[x, y, tx, ty, z]
= 0;

s.t. island_contains_hint_location{(x, y) in V, z in Z}:

is_island[x, y, z] >= if islands[x, y, z] > 0 then 1 else 0;

s.t. each_square_is_river_or_island{(x, y) in V}:

is_river[x, y] + sum{z in Z} is_island[x, y, z] = 1;

s.t. islands_match_hint_area{z in Z}:

sum{(x, y) in V} is_island[x, y, z] = island_area[z];

s.t. river_has_no_lakes{(x,y) in V: x > 0 and y > 0}:

3 >= is_river[x, y] + is_river[x - 1, y - 1]
+ is_river[x - 1, y] + is_river[x, y - 1];

s.t. river_squares_have_at_least_one_neighbor{(x, y) in V}:

sum{(sx, sy, tx, ty) in E: x = tx and y = ty} is_river[sx, sy]
>= is_river[x, y];

s.t. island_fields_have_at_least_one_neighbor{(x, y) in V, z in Z: island_area[z] > 1}:

sum{(sx, sy, tx, ty) in E: x = tx and y = ty} is_island[sx, sy, z]
>= is_island[x, y, z];

s.t. islands_are_separated_by_water{(x, y) in V, z in Z}:

sum{(sx, sy, tx, ty) in E, oz in Z: x = sx and y = sy and z != oz} is_island[tx, ty, oz]
<= 4 * P * (1 - is_island[x, y, z]);

solve;

for{y in M-1..0 by -1}
{
for {x in X}
{
printf "%s", if is_river[x, y] = 1 then "." else "#";
}
printf "\n";
}


## Problem data

### 5 x 5

data;
param M := 5;
param N := 5;
param P := 5;
param islands :=
# x,y,z,area
0,1,0,3
4,1,1,1
0,4,2,2
2,4,3,2
4,4,4,2;
end;


### 7 x 7

data;
param M := 7;
param N := 7;
param P := 8;
param islands :=
# x,y,z,area
0,0,0,2
3,1,1,2
6,1,2,2
4,3,3,2
2,4,4,2
0,6,5,8
2,6,6,1
4,6,7,3;
end;


### 15 x 15

data;
param M := 15;
param N := 15;
param P := 34;
param islands :=
# x,y,   z,area
5,  1,   0, 1
9,  1,   1, 3
4,  2,   2, 5
6,  2,   3, 1
11, 2,   4, 2
2,  3,   5, 2
9,  3,   6, 3
4,  4,   7, 2
10, 4,   8, 1
12, 4,   9, 5
1,  5,  10, 3
3,  5,  11, 1
8,  5,  12, 3
13, 5,  13, 1
5,  6,  14, 5
12, 6,  15, 1
2,  8,  16, 1
9,  8,  17, 2
1,  9,  18, 1
6,  9,  19, 2
11, 9,  20, 6
13, 9,  21, 3
2,  10, 22, 5
4,  10, 23, 2
10, 10, 24, 4
5,  11, 25, 1
12, 11, 26, 2
3,  12, 27, 2
8,  12, 28, 2
10, 12, 29, 5
5,  13, 30, 1
9,  13, 31, 1
6,  14, 32, 1
8,  14  33, 1;
end;


## Usage

Have glpsol installed, model as one file (e.g. river.mod), input data in separate file(s) (e.g. 7x7.mod). Then:

glpsol -m river.mod -d 7x7.mod


This plus some patience will print the solution in the specified output format (together with "some" diagnostic output).

• I think this answer should be considered to be competing as it is possible for other people to verify that it works. The differences in IO format should be covered by the assumption that any reasonable IO format should be accepted. – fəˈnɛtɪk Jul 13 '18 at 21:04
• @fəˈnɛtɪk Agreed. It wasn’t eligible for ngn’s bounty that just ended, which specifically required an answer runnable on TIO, but that isn’t a requirement of the question itself. – Anders Kaseorg Jul 13 '18 at 22:38
• Given that I haven't started golfing it, I will not consider it competing (yet). Once I'm sure that I have pruned all redundant constraints, I'll one-char all declarations. – ojdo Jul 16 '18 at 9:11

# Python 3, 1295 bytes

This is a python only solution. It uses no libraries and loads the board through standard input. Further explanation to come. This is my first attempt at such a large golf. There is a link to the commented and none-golfed code at the bottom.

L,S,O,R,F=len,set,None,range,frozenset
U,N,J,D,I=S.update,F.union,F.isdisjoint,F.difference,F.intersection
def r(n,a,c):
U(c,P)
if L(I(N(Q[n],C[n]),a))<2:return 1
w=D(P,N(a,[n]));e=S();u=S([next(iter(w))])
while u:n=I(Q[u.pop()],w);U(u,D(n,e));U(e,n)
return L(e)==L(w)
def T(a,o,i,c,k):
s,p,m=a
for _ in o:
t=s,p,N(m,[_]);e=D(o,[_])
if t[2] in c:o=e;continue
if not J(i,n)or not r(_,N(m,i),k):o=e
elif s==L(t[2]):yield t
else:yield from T(t,N(e,n),i,c,k)
s,*p=input().split()
X,Y=eval(s)
A=[]
l=1,-1,0,0
P=F((x,y)for y in R(Y)for x in R(X))
exec("Q%sl,l[::-1]%s;C%s(1,1,-1,-1),l[:2]*2%s"%(('={(x,y):F((x+i,y+j)for i,j in zip(',')if X>x+i>-1<y+j<Y)for x,y in P}')*2))
for a in p:a,x,y=eval(a);k=x,y;A+=[(a,k,F([k]))]
A.sort(reverse=1)
k=F(a[1]for a in A)
p=[O]*L([a for a in A if a[0]!=1])
g,h=p[:],p[:]
i=0
while 1:
if g[i]is O:h[i]=S();f=O;g[i]=T(A[i],Q[A[i][1]],D(N(k,*p[:i]),[A[i][1]]),S(),h[i])
try:p[i]=g[i].send(f)[2]
except:
f=I(N(k,*p[:i]),h[i]);g[i]=p[i]=O;i-=1
while J(p[i],f):g[i]=p[i]=O;i-=1
else:
i+=1
if i==L(p):
z=N(k,*p)
if not any(J(z,F(zip([x,x+1]*2,[y,y,y+1,y+1])))for x in R(X-1)for y in R(Y-1)):break
for c in h:U(c,z)
b=[X*['.']for i in R(Y)]
for x,y in z:b[y][x]='#'
for l in b[::-1]:print(''.join(l))


Try it online!

Look at the un-golfed code.