Map of Islands (and a river)

Question

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:

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

Answer 1

C + PicoSAT, 2345 995 952 bytes

#include<picosat.h>
#define f(i,a)for(i=a;i;i--)
#define g(a)picosat_add(E,a)
#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:

3. 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.
4. 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.
5. 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).
6. 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.
7. 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

Answer 2

GLPK/GNU MathProg, 1286 1227 1086 998 947 bytes

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 (golfed)

param M;param N;param P;set X:=0..N-1;set Y:=0..M-1;set Z:=0..P-1;set V:=X cross Y;set E:=setof{(p,q)in V,(r,s)in V:((abs(p-r)=1)&&(q=s))||((p=r)&&(abs(q-s)=1))}(p,q,r,s);param I{X,Y,Z},default 0;param t{z in Z}:=sum{x in X,y in Y}I[x,y,z];var u{X,Y},binary;var v{X,Y,Z},integer;var w{E,Z}>=0;minimize o:sum{(p,q,r,s)in E,z in Z}w[p,q,r,s,z];a{(x,y)in V,z in Z}:I[x,y,z]+sum{(p,q,r,s)in E:x=r&&y=s}w[p,q,x,y,z]-sum{(p,q,r,s)in E:x=p&&y=q}w[x,y,r,s,z]=v[x,y,z];b{(x,y)in V,z in Z}:v[x,y,z]*99>=I[x,y,z];c{(x,y)in V}:u[x,y]+sum{z in Z}v[x,y,z]=1;d{z in Z}:sum{(x,y)in V}v[x,y,z]=t[z];e{(x,y)in V:x>0&&y>0}:3>=u[x,y]+u[x-1,y-1]+u[x-1,y]+u[x,y-1];f{(x,y)in V}:sum{(p,q,r,s)in E:x=r&&y=s}u[p,q]>=u[x,y];g{(x,y)in V,z in Z:t[z]>1}:sum{(p,q,r,s)in E:x=r&&y=s}v[p,q,z]>=v[x,y,z];h{(x,y)in V,z in Z}:sum{(p,q,r,s)in E,oz in Z:x=p&&y=q&&z!=oz}v[r,s,oz]<=4*P*(1-v[x,y,z]);solve;for{y in M-1..0 by -1}{for{x in X}printf if u[x,y]=1 then"."else"#";printf"\n";}

Code (ungolfed)

# SETS
param M > 0, integer; # height
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 I

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 I{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} I[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. I_are_connected{(x, y) in V, z in Z}:

    I[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 I[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. I_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. I_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 I:=
# 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 I:=
# 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 I:=
# 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).

Answer 3

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
  c.add(t[2]);n=D(Q[_],m);U(k,n)
  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.