Reputation: 67
The Hitori puzzle consists of an nxn grid filled with numbers 1..n, some of them repeated in rows and/or columns. The goal is to cover some of the repeated numbers so that
i) there are no repeated numbers in any row or column,
ii) no covered square is vertically or horizontally adjacent to another covered square,
iii) and the uncovered squares are vertically or horizontally connected (in a single graph), ie there are no isolated islands of uncovered cells.
Below is a MiniZinc model and a small dataset.
When this model is run using the given dataset, it produces the following output.
Running hitori.mzn, hitori.dzn
79msec
hitorimask hitorisoln
0 1 0 1 1 . 5 . 4 3
1 1 1 1 0 4 2 5 1 .
1 0 1 0 1 5 . 1 . 4
0 1 0 1 1 . 4 . 5 2
1 1 1 1 0 2 1 4 3 .
% time elapsed: 86msec
----------
hitorimask hitorisoln
0 1 0 1 0 . 5 . 4 .
1 1 1 1 1 4 2 5 1 3
1 0 1 0 1 5 . 1 . 4
0 1 0 1 1 . 4 . 5 2
1 1 1 1 0 2 1 4 3 .
% time elapsed: 86msec
----------
==========
Finished in 79msec.
The first “solution” shows two regions of connected uncovered cells. The second solution is the correct one, where all the open cells are connected.
My last constraint in the MiniZinc model is a first stab at the connectedness condition. It simply avoids single cell islands. If one comments it out, 33 solutions are generated, many with single isolated cells.
Question: How can I model the connected condition in MiniZinc?
I have thought of modelling the connections of a candidate solution as an n^2 by n^2 Markov transition matrix and then riasing it to the n^2 power. I have not actually tried it, but it seems inefficient and awefully inelegant.
Suggestions will be appreciated.
Thank you.
% hitori.mzn
% Given square a grid of sidelength n
% The grid is filled with numbers 1..n, some of them repeated in rows and columns
% Goal: Cover some of the repeated numbers so that
% i) there are no duplicates in any row of column
% ii) no covered square is vertically of horizontally adjacent to another coverde square
% iii) the uncovered blocks are connected
% data instance defined in *.dzn
int: nsize ;
array[SDX,SDX] of SDX: hitoritabl ;
set of int: SDX = 1..nsize ;
set of int: SDX0 = 0..nsize ;
set of int: SDX2 = 2..nsize ;
% primary unknowns, 0 or 1 whether whether a square is covered or open
array[SDX,SDX] of var 0..1 : hitorimask ;
% modelling variable to display the solution, 0 is used to indicate a covered cell, ie a blacked cell
array[SDX,SDX] of var SDX0: hitorisoln ;
constraint forall( ir in SDX, jc in SDX ) ( hitorisoln[ir,jc] = hitorimask[ir,jc] * hitoritabl[ir,jc] ) ;
% model constraints
% use primary variables to construct quantity on which nonrepeat constraint can be imposed
include "alldifferent_except_0.mzn" ;
constraint forall ( ir in SDX) ( alldifferent_except_0 ( [ hitorisoln[ir,jc]| jc in SDX ] ) ) ; % no repeats in a row
constraint forall ( jc in SDX) ( alldifferent_except_0 ( [ hitorisoln[ir,jc]| ir in SDX ] ) ) ; % no repeats in a column
% impose condition on non-adjacent masked cells
% no black square directly below anothe black square
constraint forall (ir in SDX2, jc in SDX) ( 0 < hitorimask[ir-1,jc] + hitorimask[ir,jc] ) ;
% no black square directly to the right of anothe black square
constraint forall (ir in SDX, jc in SDX2) ( 0 < hitorimask[ir,jc-1] + hitorimask[ir,jc] ) ;
% this constraint simply excludes islands of singleton open cells
constraint forall (ir in SDX, jc in SDX)
( if hitorimask[ir,jc]=1 then
( if ir > 1 then hitorimask[ir-1,jc ] else 0 endif
+ if jc < nsize then hitorimask[ir ,jc+1] else 0 endif
+ if ir < nsize then hitorimask[ir+1,jc ] else 0 endif
+ if jc > 1 then hitorimask[ir ,jc-1] else 0 endif
> 0
)
endif ) ;
%
% still need a condition of connectedness of open squares
%
solve satisfy ;
output ["\nhitorimask hitorisoln \n" ] ;
output [
if jc <= nsize then "" ++ show_int(3, hitorimask[ir,jc] ) else "" endif
++ if jc = nsize then " " else "" endif
++ if nsize < jc /\ jc <= 2*nsize
then if fix(hitorisoln[ir,jc-nsize])>0
then show_int(3, hitorisoln[ir,jc-nsize] )
else " ."
endif
else ""
endif
++ if jc = 2*nsize then "\n" else "" endif
| ir in SDX, jc in 1..2*nsize ] ;
data file
% hitori.dzn
nsize = 5 ; % size of square grid
hitoritabl = [|
4, 5, 4, 4, 3 |
4, 2, 5, 1, 3 |
5, 1, 1, 3, 4 |
5, 4, 2, 5, 2 |
2, 1, 4, 3, 1 |] ;
Upvotes: 1
Views: 259
Reputation: 111
This second approach to ensure connectivity among all unshaded cells is based on flow conservation:
2. flow conservation method
Edit:
adding PoC (Python CP-SAT):
performance is slower than method (3)
### FLOW CONSERVATION
from ortools.sat.python import cp_model
import numpy as np
def check_connectivity_flow(matrix):
model = cp_model.CpModel()
rows, cols = matrix.shape
# Identify all True cells and count them.
true_cells = [(r, c) for r in range(rows) for c in range(cols) if matrix[r, c]]
num_true = len(true_cells)
# Pick first True cell as a root
root = true_cells[0]
# Allowed edges: adjacent in same row (col diff 1) or same column (row diff 1)
allowed_edges = []
for (r, c) in true_cells:
for (dr, dc) in [(0, 1), (1, 0), (0, -1), (-1, 0)]:
nr, nc = r + dr, c + dc
if (nr, nc) in true_cells:
allowed_edges.append(((r, c), (nr, nc)))
# Flow variables for each allowed directed edge
flow_vars = {}
for (u, v) in allowed_edges:
flow_vars[(u, v)] = model.new_int_var(0, num_true - 1, f"flow_{u}_{v}")
flow_vars[(v, u)] = model.new_int_var(0, num_true - 1, f"flow_{v}_{u}")
flow_direct = model. new_bool_var(f"flow_direct_{u}_{v}")
# Prevent bidirectional flow
model.add(flow_vars[(u, v)] > 0).only_enforce_if(flow_direct)
model.add(flow_vars[(v, u)] == 0).only_enforce_if(flow_direct)
model.add(flow_vars[(v, u)] > 0).only_enforce_if(flow_direct.Not())
model.add(flow_vars[(u, v)] == 0).only_enforce_if(flow_direct.Not())
# Flow conservation
for cell in true_cells:
out_flow = sum(flow_vars[(cell, v)] for (u, v) in flow_vars if u == cell)
in_flow = sum(flow_vars[(u, cell)] for (u, v) in flow_vars if v == cell)
if cell == root:
# Root must send out exactly num_true - 1 units of flow
model.add(out_flow == num_true - 1)
# Optionally, force in_flow == 0
model.add(in_flow == 0)
else:
# Every other cell must have net inflow 1
model.add(in_flow == out_flow + 1)
# Solve the model
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.ResponseStats())
print(f"Input mask:\n{matrix}")
# If the model is infeasible, connectivity cannot be established.
if status != cp_model.OPTIMAL and status != cp_model.FEASIBLE:
return False
# Flow difference matrix for each cell
flow_diff_matrix = np.zeros((rows, cols), dtype=int)
in_flow_matrix = np.zeros((rows, cols), dtype=int)
out_flow_matrix = np.zeros((rows, cols), dtype=int)
for cell in true_cells:
in_flow = sum(solver.value(flow_vars[(u, cell)]) for (u, v) in flow_vars if v == cell)
out_flow = sum(solver.value(flow_vars[(cell, v)]) for (u, v) in flow_vars if u == cell)
flow_diff_matrix[cell[0], cell[1]] = in_flow - out_flow
in_flow_matrix[cell[0], cell[1]] = in_flow
out_flow_matrix[cell[0], cell[1]] = out_flow
print(f"\nflow conservation:\n{flow_diff_matrix}")
print(f"\nin_flow:\n{in_flow_matrix}")
print(f"\nout_flow:\n{out_flow_matrix}")
return True
# Examples
small_mask = np.array([
[True, False, False],
[True, False, True ],
[True, False, True ]
])
true_mask = np.array([
[0, 1, 0, 1, 0],
[1, 1, 1, 1, 1],
[1, 0, 1, 0, 1],
[0, 1, 0, 1, 1],
[1, 1, 1, 1, 0]
])
false_mask = np.array([
[0, 1, 0, 1, 1],
[1, 1, 1, 1, 0],
[1, 0, 1, 0, 1],
[0, 1, 0, 1, 1],
[1, 1, 1, 1, 0]
])
bigger_mask = np.array([
[1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1],
[0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0],
[1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1],
[1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0],
[0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0],
[1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0],
[1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1],
[1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0],
[1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1],
[0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1],
[1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1],
[1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1],
[0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1],
[1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1],
[1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1]
])
print(f"\n> Is Connected: {check_connectivity_flow(false_mask)}\n***************\n")
print(f"\n> Is Connected: {check_connectivity_flow(bigger_mask)}\n***************\n")
Output:
CpSolverResponse summary:
status: INFEASIBLE
objective: NA
best_bound: NA
integers: 0
booleans: 0
conflicts: 0
branches: 0
propagations: 0
integer_propagations: 0
restarts: 0
lp_iterations: 0
walltime: 0.000847509
usertime: 0.000847564
deterministic_time: 0
gap_integral: 0
Input mask:
[[0 1 0 1 1]
[1 1 1 1 0]
[1 0 1 0 1]
[0 1 0 1 1]
[1 1 1 1 0]]
> Is Connected: False
***************
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
integers: 263
booleans: 131
conflicts: 0
branches: 262
propagations: 512
integer_propagations: 2161
restarts: 262
lp_iterations: 126
walltime: 0.0671798
usertime: 0.0671798
deterministic_time: 0.00821366
gap_integral: 0
solution_fingerprint: 0x611a630ef34b9365
Input mask:
[[1 1 1 1 0 1 1 0 1 1 1 0 1 1 0 1 1]
[0 1 1 0 1 0 1 1 0 1 1 1 1 0 1 1 0]
[1 0 1 1 1 1 1 0 1 1 0 1 0 1 0 1 1]
[1 1 1 0 1 0 1 1 0 1 1 1 1 1 1 1 0]
[0 1 0 1 0 1 0 1 1 1 0 1 1 1 0 1 1]
[1 1 1 1 1 1 1 1 0 1 0 1 1 0 1 1 0]
[0 1 1 0 1 1 0 1 1 0 1 0 1 1 0 1 0]
[1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 1 0]
[1 0 1 1 1 0 1 0 1 0 1 0 1 1 0 1 1]
[1 1 0 1 0 1 1 1 1 1 0 1 0 1 1 1 0]
[1 0 1 1 1 1 0 1 0 1 1 1 1 0 1 0 1]
[0 1 0 1 0 1 1 1 1 0 1 0 1 1 0 1 1]
[1 1 1 0 1 0 1 0 1 1 0 1 1 1 1 0 1]
[1 1 0 1 1 1 0 1 0 1 1 1 0 1 0 1 1]
[0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 1]
[1 1 0 1 1 0 1 1 1 0 1 1 0 1 1 1 1]
[1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 0 1]]
flow conservation:
[[-201 1 1 1 0 1 1 0 1 1 1 0 1 1 0 1 1]
[ 0 1 1 0 1 0 1 1 0 1 1 1 1 0 1 1 0]
[ 1 0 1 1 1 1 1 0 1 1 0 1 0 1 0 1 1]
[ 1 1 1 0 1 0 1 1 0 1 1 1 1 1 1 1 0]
[ 0 1 0 1 0 1 0 1 1 1 0 1 1 1 0 1 1]
[ 1 1 1 1 1 1 1 1 0 1 0 1 1 0 1 1 0]
[ 0 1 1 0 1 1 0 1 1 0 1 0 1 1 0 1 0]
[ 1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 1 0]
[ 1 0 1 1 1 0 1 0 1 0 1 0 1 1 0 1 1]
[ 1 1 0 1 0 1 1 1 1 1 0 1 0 1 1 1 0]
[ 1 0 1 1 1 1 0 1 0 1 1 1 1 0 1 0 1]
[ 0 1 0 1 0 1 1 1 1 0 1 0 1 1 0 1 1]
[ 1 1 1 0 1 0 1 0 1 1 0 1 1 1 1 0 1]
[ 1 1 0 1 1 1 0 1 0 1 1 1 0 1 0 1 1]
[ 0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 1]
[ 1 1 0 1 1 0 1 1 1 0 1 1 0 1 1 1 1]
[ 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 0 1]]
in_flow:
[[ 0 201 131 1 0 1 2 0 1 6 4 0 2 1 0 2 1]
[ 0 69 197 0 1 0 4 1 0 7 3 4 3 0 1 4 0]
[ 1 0 196 91 90 87 86 0 1 8 0 4 0 1 0 6 1]
[ 2 103 104 0 1 0 81 80 0 54 45 44 16 18 16 15 0]
[ 0 100 0 1 0 1 0 79 57 56 0 23 17 10 0 8 1]
[ 1 99 89 88 86 24 9 21 0 1 0 11 16 0 1 6 0]
[ 0 16 9 0 69 90 0 11 1 0 1 0 15 1 0 4 0]
[ 5 6 0 1 0 89 96 9 0 1 4 5 13 0 1 3 0]
[ 4 0 1 4 1 0 95 0 1 0 1 0 7 6 0 2 1]
[ 3 1 0 5 0 45 94 48 18 16 0 1 0 5 4 2 0]
[ 1 0 1 8 9 44 0 29 0 15 14 12 10 0 1 0 1]
[ 0 1 0 1 0 34 33 59 58 0 1 0 10 9 0 1 3]
[ 3 6 1 0 1 0 1 0 57 56 0 2 2 8 1 0 4]
[ 2 7 0 17 19 20 0 1 0 55 21 20 0 5 0 1 6]
[ 0 10 11 16 0 21 26 28 32 33 0 17 12 22 8 0 7]
[ 2 3 0 4 1 0 4 2 4 0 1 4 0 21 20 10 9]
[ 1 0 1 2 0 1 2 0 2 1 0 2 1 0 1 0 1]]
out_flow:
[[201 200 130 0 0 0 1 0 0 5 3 0 1 0 0 1 0]
[ 0 68 196 0 0 0 3 0 0 6 2 3 2 0 0 3 0]
[ 0 0 195 90 89 86 85 0 0 7 0 3 0 0 0 5 0]
[ 1 102 103 0 0 0 80 79 0 53 44 43 15 17 15 14 0]
[ 0 99 0 0 0 0 0 78 56 55 0 22 16 9 0 7 0]
[ 0 98 88 87 85 23 8 20 0 0 0 10 15 0 0 5 0]
[ 0 15 8 0 68 89 0 10 0 0 0 0 14 0 0 3 0]
[ 4 5 0 0 0 88 95 8 0 0 3 4 12 0 0 2 0]
[ 3 0 0 3 0 0 94 0 0 0 0 0 6 5 0 1 0]
[ 2 0 0 4 0 44 93 47 17 15 0 0 0 4 3 1 0]
[ 0 0 0 7 8 43 0 28 0 14 13 11 9 0 0 0 0]
[ 0 0 0 0 0 33 32 58 57 0 0 0 9 8 0 0 2]
[ 2 5 0 0 0 0 0 0 56 55 0 1 1 7 0 0 3]
[ 1 6 0 16 18 19 0 0 0 54 20 19 0 4 0 0 5]
[ 0 9 10 15 0 20 25 27 31 32 0 16 11 21 7 0 6]
[ 1 2 0 3 0 0 3 1 3 0 0 3 0 20 19 9 8]
[ 0 0 0 1 0 0 1 0 1 0 0 1 0 0 0 0 0]]
> Is Connected: True
***************
Addendum
This is the result for a 15x15 Hitori on Colab:
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
integers: 356
booleans: 46
conflicts: 0
branches: 92
propagations: 1058
integer_propagations: 2398
restarts: 92
lp_iterations: 350
walltime: 0.10211
usertime: 0.102111
deterministic_time: 0.00997731
gap_integral: 0
solution_fingerprint: 0xa5547dd0b45d7440
Input:
[ 15 3 1 4 15 13 1 6 15 5 7 2 1 7 11 ]
[ 8 13 6 13 4 14 13 7 11 12 10 1 15 5 1 ]
[ 6 1 15 8 6 15 4 6 13 14 15 3 5 6 9 ]
[ 11 12 3 12 12 1 12 15 14 12 5 12 8 2 4 ]
[ 6 8 3 1 13 6 5 9 3 15 3 4 10 11 9 ]
[ 10 6 11 10 7 9 5 5 3 6 2 15 10 4 8 ]
[ 12 11 7 2 12 6 15 4 7 9 13 7 3 12 10 ]
[ 3 11 4 11 5 11 9 11 6 2 11 14 2 10 11 ]
[ 6 5 6 9 6 12 6 14 6 10 15 6 11 6 3 ]
[ 4 14 2 13 15 5 11 3 7 8 1 10 5 9 3 ]
[ 11 2 14 3 8 10 11 1 11 7 3 11 13 3 6 ]
[ 2 4 12 15 1 12 6 12 8 12 11 5 12 14 12 ]
[ 5 13 12 14 13 8 7 11 13 4 9 13 2 13 15 ]
[ 14 15 5 12 9 4 12 10 2 11 1 13 1 12 7 ]
[ 13 7 13 6 4 3 2 4 10 13 14 3 4 15 4 ]
Mask (unshaded == 1):
[ 1 1 1 1 0 1 0 1 0 1 0 1 0 1 1 ]
[ 1 0 1 0 1 1 1 1 1 1 1 1 1 1 0 ]
[ 0 1 1 1 1 0 1 0 1 1 0 1 1 0 1 ]
[ 1 0 1 0 1 1 0 1 1 0 1 0 1 1 1 ]
[ 1 1 0 1 1 0 1 1 0 1 1 1 1 1 0 ]
[ 0 1 1 1 1 1 0 1 1 0 1 1 0 1 1 ]
[ 1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 ]
[ 1 0 1 1 1 0 1 0 1 1 0 1 0 1 0 ]
[ 0 1 0 1 0 1 0 1 0 1 1 0 1 1 1 ]
[ 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 ]
[ 0 1 1 0 1 1 0 1 0 1 0 1 1 1 1 ]
[ 1 1 0 1 1 0 1 1 1 0 1 1 0 1 0 ]
[ 1 0 1 1 0 1 1 1 0 1 1 0 1 1 1 ]
[ 1 1 1 1 1 1 0 1 1 1 0 1 1 0 1 ]
[ 0 1 0 1 0 1 1 0 1 1 1 0 1 1 0 ]
Solution:
[ 15 3 1 4 _ 13 _ 6 _ 5 _ 2 _ 7 11 ]
[ 8 _ 6 _ 4 14 13 7 11 12 10 1 15 5 _ ]
[ _ 1 15 8 6 _ 4 _ 13 14 _ 3 5 _ 9 ]
[ 11 _ 3 _ 12 1 _ 15 14 _ 5 _ 8 2 4 ]
[ 6 8 _ 1 13 _ 5 9 _ 15 3 4 10 11 _ ]
[ _ 6 11 10 7 9 _ 5 3 _ 2 15 _ 4 8 ]
[ 12 11 _ 2 _ 6 15 4 _ 9 13 7 3 _ 10 ]
[ 3 _ 4 11 5 _ 9 _ 6 2 _ 14 _ 10 _ ]
[ _ 5 _ 9 _ 12 _ 14 _ 10 15 _ 11 6 3 ]
[ 4 14 2 13 15 5 11 3 7 8 1 10 _ 9 _ ]
[ _ 2 14 _ 8 10 _ 1 _ 7 _ 11 13 3 6 ]
[ 2 4 _ 15 1 _ 6 12 8 _ 11 5 _ 14 _ ]
[ 5 _ 12 14 _ 8 7 11 _ 4 9 _ 2 13 15 ]
[ 14 15 5 12 9 4 _ 10 2 11 _ 13 1 _ 7 ]
[ _ 7 _ 6 _ 3 2 _ 10 13 14 _ 4 15 _ ]
If we don’t need to see the inflows and outflows of the unshaded cells separately, we can simplify the model.
By directly aggregating all flows of a cell at once, we can reduce the numbers of decision variables and constraints.This simplified model matches more with method (3), both in terms of complexity and performance.
Stats for the same 15x15 Hitori (simplified model):
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
integers: 36
booleans: 0
conflicts: 0
branches: 0
propagations: 0
integer_propagations: 0
restarts: 0
lp_iterations: 39
walltime: 0.0686117
usertime: 0.0686119
deterministic_time: 0.001577
gap_integral: 0
solution_fingerprint: 0x3c10622297f28acb
Stats for the 20x20 Hitori (simplified model):
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
integers: 33
booleans: 0
conflicts: 0
branches: 0
propagations: 0
integer_propagations: 0
restarts: 0
lp_iterations: 28
walltime: 0.11172
usertime: 0.111721
deterministic_time: 0.00226062
gap_integral: 0
solution_fingerprint: 0x2274178dfe825f1e
Upvotes: 1
Reputation: 111
I figured out 2 approaches to formulate the connectedness constraint for all unshaded (not covered) cells.
The first approach is based on transitive relation among the unshaded cells:
1. transitivity method
Edit:
The (1) transitivity method struggles to scale effectively … therefore, replacing it with
the (3) min. distance method, which outperforms the (2) flow conservation method as well:
3. minimum distance method
PoC (Python CP-SAT) to check a given matrix whether the True/1-cells form a connected area:
note: it was easier for me to setup colab with ortools than with minizinc and additional solvers, and I am also more familiar with CP-SAT
### MIN DISTANCE
from ortools.sat.python import cp_model
import numpy as np
import math
def check_connectivity_min_dist(matrix):
model = cp_model.CpModel()
rows, cols = matrix.shape
# Max distance is n * n
max_distance = rows * cols
# Identify all True cells and count them.
true_cells = [(r, c) for r in range(rows) for c in range(cols) if matrix[r, c]]
# Create a distance variable for each cell (including non-colored cells)
distance_vars = {}
for (r, c) in true_cells:
distance_vars[(r, c)] = model.new_int_var(0, max_distance, f"distance_{r}_{c}")
# Set distance of root cell (first true cell) to 0
root = true_cells[0]
# Function to get neighbors for any cell
def get_neighbors(r, c):
neighbors = []
for dr, dc in [(0, 1), (1, 0), (0, -1), (-1, 0)]:
nr, nc = r + dr, c + dc
if 0 <= nr < rows and 0 <= nc < cols:
neighbors.append((nr, nc))
return neighbors
# Distance Constraints:
for (r, c) in true_cells:
if (r, c) == root:
model.Add(distance_vars[(r, c)] == 0)
else:
min_neighbor_distance = model.new_int_var(0, max_distance, f"min_neighbor_distance_{r}_{c}")
neighbors = [distance_vars[nr, nc] for nr, nc in get_neighbors(r, c) if matrix[nr, nc]]
model.add_min_equality(target=min_neighbor_distance, exprs=neighbors)
model.add(distance_vars[(r, c)] == 1 + min_neighbor_distance)
model.add(distance_vars[(r, c)] < max_distance)
# Solve the model
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.ResponseStats())
print(f"Input mask:\n{matrix}")
# If the model is infeasible, connectivity cannot be established.
if status != cp_model.OPTIMAL and status != cp_model.FEASIBLE:
return False
# Distance matrix
distance_matrix = np.full((rows, cols), -1, dtype=int)
for (r, c) in true_cells:
distance_matrix[r, c] = solver.Value(distance_vars[(r, c)])
print(f"\nDistances:\n{distance_matrix}")
return True
# Examples
small_mask = np.array([
[True, False, False],
[True, False, True ],
[True, False, True ]
])
true_mask = np.array([
[0, 1, 0, 1, 0],
[1, 1, 1, 1, 1],
[1, 0, 1, 0, 1],
[0, 1, 0, 1, 1],
[1, 1, 1, 1, 0]
])
false_mask = np.array([
[0, 1, 0, 1, 1],
[1, 1, 1, 1, 0],
[1, 0, 1, 0, 1],
[0, 1, 0, 1, 1],
[1, 1, 1, 1, 0]
])
bigger_mask = np.array([
[1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1],
[0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0],
[1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1],
[1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0],
[0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0],
[1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0],
[1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1],
[1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0],
[1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1],
[0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1],
[1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1],
[1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1],
[0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1],
[1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1],
[1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1]
])
print(f"\n> Is Connected: {check_connectivity_min_dist(false_mask)}\n***************\n")
print(f"\n> Is Connected: {check_connectivity_min_dist(bigger_mask)}\n***************\n")
Output:
CpSolverResponse summary:
status: INFEASIBLE
objective: NA
best_bound: NA
integers: 0
booleans: 0
conflicts: 0
branches: 0
propagations: 0
integer_propagations: 0
restarts: 0
lp_iterations: 0
walltime: 0.000979077
usertime: 0.000979146
deterministic_time: 0
gap_integral: 0
Input mask:
[[0 1 0 1 1]
[1 1 1 1 0]
[1 0 1 0 1]
[0 1 0 1 1]
[1 1 1 1 0]]
> Is Connected: False
***************
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
integers: 0
booleans: 0
conflicts: 0
branches: 0
propagations: 0
integer_propagations: 0
restarts: 0
lp_iterations: 0
walltime: 0.0105172
usertime: 0.0105172
deterministic_time: 6.5e-07
gap_integral: 0
solution_fingerprint: 0xd22cc995cb38921f
Input mask:
[[1 1 1 1 0 1 1 0 1 1 1 0 1 1 0 1 1]
[0 1 1 0 1 0 1 1 0 1 1 1 1 0 1 1 0]
[1 0 1 1 1 1 1 0 1 1 0 1 0 1 0 1 1]
[1 1 1 0 1 0 1 1 0 1 1 1 1 1 1 1 0]
[0 1 0 1 0 1 0 1 1 1 0 1 1 1 0 1 1]
[1 1 1 1 1 1 1 1 0 1 0 1 1 0 1 1 0]
[0 1 1 0 1 1 0 1 1 0 1 0 1 1 0 1 0]
[1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 1 0]
[1 0 1 1 1 0 1 0 1 0 1 0 1 1 0 1 1]
[1 1 0 1 0 1 1 1 1 1 0 1 0 1 1 1 0]
[1 0 1 1 1 1 0 1 0 1 1 1 1 0 1 0 1]
[0 1 0 1 0 1 1 1 1 0 1 0 1 1 0 1 1]
[1 1 1 0 1 0 1 0 1 1 0 1 1 1 1 0 1]
[1 1 0 1 1 1 0 1 0 1 1 1 0 1 0 1 1]
[0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 1]
[1 1 0 1 1 0 1 1 1 0 1 1 0 1 1 1 1]
[1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 0 1]]
Distances:
[[ 0 1 2 3 -1 11 10 -1 18 17 18 -1 20 21 -1 23 24]
[-1 2 3 -1 7 -1 9 10 -1 16 17 18 19 -1 23 22 -1]
[ 8 -1 4 5 6 7 8 -1 16 15 -1 17 -1 19 -1 21 22]
[ 7 6 5 -1 7 -1 9 10 -1 14 15 16 17 18 19 20 -1]
[-1 7 -1 11 -1 13 -1 11 12 13 -1 17 18 19 -1 21 22]
[ 9 8 9 10 11 12 13 12 -1 14 -1 18 19 -1 23 22 -1]
[-1 9 10 -1 12 13 -1 13 14 -1 24 -1 20 21 -1 23 -1]
[11 10 -1 24 -1 14 15 14 -1 24 23 22 21 -1 25 24 -1]
[12 -1 24 23 24 -1 16 -1 20 -1 24 -1 22 23 -1 25 26]
[13 14 -1 22 -1 18 17 18 19 20 -1 24 -1 24 25 26 -1]
[14 -1 22 21 20 19 -1 19 -1 21 22 23 24 -1 26 -1 38]
[-1 38 -1 22 -1 20 21 20 21 -1 23 -1 25 26 -1 38 37]
[38 37 38 -1 32 -1 22 -1 22 23 -1 27 26 27 28 -1 36]
[37 36 -1 32 31 30 -1 28 -1 24 25 26 -1 28 -1 36 35]
[-1 35 34 33 -1 29 28 27 26 25 -1 27 28 29 30 -1 34]
[37 36 -1 34 35 -1 29 28 27 -1 29 28 -1 30 31 32 33]
[38 -1 36 35 -1 31 30 -1 28 29 -1 29 30 -1 32 -1 34]]
> Is Connected: True
***************
Addendum
This is the result for a 20x20 Hitori on Colab:
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
integers: 207
booleans: 1
conflicts: 0
branches: 1
propagations: 0
integer_propagations: 420
restarts: 0
lp_iterations: 0
walltime: 0.0888325
usertime: 0.0888326
deterministic_time: 0.00802969
gap_integral: 0
solution_fingerprint: 0x88d601be3265190d
Input:
[ 7 10 16 11 2 9 2 3 14 14 17 19 4 4 13 17 13 15 5 11 ]
[ 1 18 13 19 9 8 2 1 7 5 15 7 3 4 9 3 11 11 16 16 ]
[ 9 19 6 11 7 16 14 2 3 20 11 5 12 11 15 4 18 3 14 10 ]
[ 5 6 9 11 12 8 17 17 4 8 18 8 15 14 12 19 13 18 2 1 ]
[ 20 4 11 9 7 17 13 5 1 15 19 15 7 10 12 11 2 8 18 11 ]
[ 18 12 4 17 15 20 11 7 14 1 11 9 18 13 2 6 3 18 19 17 ]
[ 11 13 3 16 18 2 20 15 7 3 2 13 8 19 2 1 17 11 2 9 ]
[ 4 8 10 3 15 14 1 4 9 13 12 2 13 18 19 2 6 18 20 6 ]
[ 10 16 5 20 1 15 9 9 18 19 6 10 6 11 7 13 12 3 17 4 ]
[ 5 14 8 1 13 4 12 7 11 16 9 11 17 19 3 5 16 19 1 12 ]
[ 4 7 7 15 11 6 19 12 11 2 19 4 16 10 18 20 9 11 8 5 ]
[ 3 17 5 6 5 11 15 19 8 18 1 12 3 15 10 10 3 13 9 13 ]
[ 10 5 9 18 20 4 10 18 16 10 7 12 6 12 1 2 8 7 3 19 ]
[ 14 3 18 2 16 5 12 13 1 9 16 11 20 11 10 8 17 7 4 17 ]
[ 1 20 19 4 17 8 16 19 19 6 17 8 16 7 4 11 1 12 13 18 ]
[ 11 12 19 2 8 18 9 14 19 17 11 19 1 5 20 5 10 8 6 14 ]
[ 19 19 17 13 4 12 3 12 15 10 14 1 18 15 8 2 4 5 16 7 ]
[ 16 11 7 18 4 3 14 7 7 15 1 3 10 9 6 17 19 2 12 2 ]
[ 7 15 12 14 3 19 4 10 17 17 20 16 1 18 14 9 1 2 7 8 ]
[ 13 3 16 8 10 2 10 4 10 16 3 18 19 8 17 12 12 6 2 20 ]
Mask (unshaded == 1):
[ 1 1 1 0 1 1 0 1 0 1 0 1 1 0 1 1 0 1 1 1 ]
[ 0 1 1 1 0 1 1 1 1 1 1 0 1 1 1 0 1 0 1 0 ]
[ 1 1 1 0 1 1 0 1 0 1 0 1 1 0 1 1 1 1 1 1 ]
[ 0 1 0 1 1 0 1 0 1 1 1 0 1 1 0 1 1 0 1 1 ]
[ 1 1 1 1 0 1 1 1 1 0 1 1 1 0 1 0 1 1 1 0 ]
[ 1 0 1 0 1 1 1 0 1 1 0 1 0 1 1 1 1 0 1 1 ]
[ 0 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 0 1 0 1 ]
[ 1 1 1 1 0 1 1 0 1 0 1 1 1 0 1 0 1 1 1 0 ]
[ 0 1 0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 ]
[ 1 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1 1 0 ]
[ 0 1 0 1 1 1 1 1 0 1 0 1 1 1 1 1 1 0 1 1 ]
[ 1 1 1 1 0 1 0 1 1 1 1 1 0 1 0 1 0 1 1 0 ]
[ 0 1 1 0 1 1 1 1 1 0 1 0 1 1 1 0 1 0 1 1 ]
[ 1 0 1 1 0 1 0 1 0 1 1 1 1 0 1 1 1 1 1 0 ]
[ 1 1 0 1 1 0 1 0 1 1 0 1 0 1 0 1 0 1 1 1 ]
[ 0 1 1 0 1 1 1 1 0 1 1 0 1 1 1 0 1 0 1 0 ]
[ 1 0 1 1 0 1 1 0 1 1 1 1 1 0 1 1 1 1 0 1 ]
[ 1 1 0 1 1 0 1 1 0 1 0 1 1 1 1 0 1 0 1 1 ]
[ 0 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 1 ]
[ 1 1 0 1 0 1 0 1 1 1 0 1 1 0 1 1 0 1 0 1 ]
Solution:
[ 7 10 16 _ 2 9 _ 3 _ 14 _ 19 4 _ 13 17 _ 15 5 11 ]
[ _ 18 13 19 _ 8 2 1 7 5 15 _ 3 4 9 _ 11 _ 16 _ ]
[ 9 19 6 _ 7 16 _ 2 _ 20 _ 5 12 _ 15 4 18 3 14 10 ]
[ _ 6 _ 11 12 _ 17 _ 4 8 18 _ 15 14 _ 19 13 _ 2 1 ]
[ 20 4 11 9 _ 17 13 5 1 _ 19 15 7 _ 12 _ 2 8 18 _ ]
[ 18 _ 4 _ 15 20 11 _ 14 1 _ 9 _ 13 2 6 3 _ 19 17 ]
[ _ 13 _ 16 18 _ 20 15 _ 3 2 _ 8 19 _ 1 _ 11 _ 9 ]
[ 4 8 10 3 _ 14 1 _ 9 _ 12 2 13 _ 19 _ 6 18 20 _ ]
[ _ 16 _ 20 1 15 _ 9 18 19 6 10 _ 11 7 13 12 _ 17 4 ]
[ 5 14 8 _ 13 _ 12 _ 11 _ 9 _ 17 _ 3 _ 16 19 1 _ ]
[ _ 7 _ 15 11 6 19 12 _ 2 _ 4 16 10 18 20 9 _ 8 5 ]
[ 3 17 5 6 _ 11 _ 19 8 18 1 12 _ 15 _ 10 _ 13 9 _ ]
[ _ 5 9 _ 20 4 10 18 16 _ 7 _ 6 12 1 _ 8 _ 3 19 ]
[ 14 _ 18 2 _ 5 _ 13 _ 9 16 11 20 _ 10 8 17 7 4 _ ]
[ 1 20 _ 4 17 _ 16 _ 19 6 _ 8 _ 7 _ 11 _ 12 13 18 ]
[ _ 12 19 _ 8 18 9 14 _ 17 11 _ 1 5 20 _ 10 _ 6 _ ]
[ 19 _ 17 13 _ 12 3 _ 15 10 14 1 18 _ 8 2 4 5 _ 7 ]
[ 16 11 _ 18 4 _ 14 7 _ 15 _ 3 10 9 6 _ 19 _ 12 2 ]
[ _ 15 12 14 3 19 4 10 17 _ 20 16 _ 18 _ 9 1 2 7 8 ]
[ 13 3 _ 8 _ 2 _ 4 10 16 _ 18 19 _ 17 12 _ 6 _ 20 ]
Distances:
[ 0 1 2 _ 14 13 _ 15 _ 17 _ 27 26 _ 28 29 _ 35 34 35 ]
[ _ 2 3 4 _ 12 13 14 15 16 17 _ 25 26 27 _ 31 _ 33 _ ]
[ 4 3 4 _ 10 11 _ 15 _ 17 _ 25 24 _ 28 29 30 31 32 33 ]
[ _ 4 _ 8 9 _ 23 _ 19 18 19 _ 23 24 _ 30 31 _ 33 34 ]
[ 6 5 6 7 _ 23 22 21 20 _ 20 21 22 _ 32 _ 32 33 34 _ ]
[ 7 _ 7 _ 25 24 23 _ 21 22 _ 22 _ 30 31 32 33 _ 35 36 ]
[ _ 31 _ 27 26 _ 24 25 _ 23 24 _ 28 29 _ 33 _ 49 _ 37 ]
[ 31 30 29 28 _ 26 25 _ 29 _ 25 26 27 _ 45 _ 47 48 49 _ ]
[ _ 31 _ 29 28 27 _ 29 28 27 26 27 _ 45 44 45 46 _ 48 49 ]
[ 33 32 33 _ 29 _ 33 _ 29 _ 27 _ 41 _ 43 _ 45 46 47 _ ]
[ _ 33 _ 31 30 31 32 33 _ 37 _ 39 40 41 42 43 44 _ 48 49 ]
[ 35 34 33 32 _ 32 _ 34 35 36 37 38 _ 42 _ 44 _ 50 49 _ ]
[ _ 35 34 _ 34 33 34 35 36 _ 38 _ 42 43 44 _ 48 _ 50 51 ]
[ 55 _ 35 36 _ 34 _ 36 _ 40 39 40 41 _ 45 46 47 48 49 _ ]
[ 54 53 _ 37 38 _ 42 _ 42 41 _ 41 _ 49 _ 47 _ 49 50 51 ]
[ _ 52 51 _ 39 40 41 42 _ 42 43 _ 47 48 49 _ 53 _ 51 _ ]
[ 52 _ 50 49 _ 41 42 _ 44 43 44 45 46 _ 50 51 52 53 _ 59 ]
[ 51 50 _ 48 47 _ 43 44 _ 44 _ 46 47 48 49 _ 53 _ 57 58 ]
[ _ 49 48 47 46 45 44 45 46 _ 48 47 _ 49 _ 55 54 55 56 57 ]
[ 51 50 _ 48 _ 46 _ 46 47 48 _ 48 49 _ 57 56 _ 56 _ 58 ]
Upvotes: 1
Reputation: 67
Not often that somebody cares about my code. So here is the code of my MiniZinc model with two data files. I also include a snippet of the output for the smaller data file run with Gecode.
Solvers and Running Times: I typically start with Gecode, because it is the default. It also gives useful output statistics. However, when Gecode tests my patience, I try Chuffed, COIN, HiGHS and OR Tools. For the 7x7 instance these five solvers all did the job in a few 100 ms. HiGHS solved the 15x15 problem in less than half a second, COIN took 6 seconds and Chuffed needed 10 seconds. I stopped Gecode after 1 minute as it had not given any output yet. OR Tools produced some output after a few seconds, but it made very slow progress and I stopped it after 2 minutes.
----------
sum(flowremain): 34 obj: 47
hitorimask hitorisoln flowremain
0 1 1 1 0 1 0 . 1 5 7 . 3 . 0 2 1 1 0 1 0
1 1 0 1 1 1 1 1 3 . 2 7 5 6 1 1 0 1 1 1 1
0 1 1 0 1 1 0 . 2 3 . 6 7 . 0 1 1 0 1 1 0
1 1 0 1 0 1 1 5 7 . 6 . 2 4 1 1 0 1 0 1 1
1 0 1 1 1 1 0 2 . 7 4 1 6 . 1 0 1 1 1 1 0
0 1 0 1 0 1 1 . 6 . 1 . 4 7 0 1 0 1 0 1 1
1 1 1 1 1 0 1 7 4 6 3 2 . 5 1 1 1 1 1 0 0
% time elapsed: 132msec
----------
sum(flowremain): 34 obj: 49
hitorimask hitorisoln flowremain
0 1 1 1 0 1 0 . 1 5 7 . 3 . 0 1 1 1 0 1 0
1 1 0 1 1 1 1 1 3 . 2 7 5 6 1 1 0 1 1 1 1
0 1 1 0 1 1 0 . 2 3 . 6 7 . 0 1 1 0 1 1 0
1 1 0 1 0 1 1 5 7 . 6 . 2 4 1 1 0 1 0 1 1
1 0 1 1 1 1 0 2 . 7 4 1 6 . 1 0 1 1 1 1 0
0 1 0 1 0 1 1 . 6 . 1 . 4 7 0 1 0 1 0 1 1
1 1 1 1 1 0 1 7 4 6 3 2 . 5 1 1 1 1 1 0 1
% time elapsed: 133msec
----------
==========
%%%mzn-stat: failures=1716
%%%mzn-stat: initTime=0.00951658
%%%mzn-stat: nodes=3504
%%%mzn-stat: peakDepth=69
%%%mzn-stat: propagations=382613
%%%mzn-stat: propagators=917
%%%mzn-stat: restarts=0
%%%mzn-stat: solutions=34
%%%mzn-stat: solveTime=0.0152159
%%%mzn-stat: variables=731
%%%mzn-stat-end
%%%mzn-stat: nSolutions=34
%%%mzn-stat-end
Finished in 159msec.
Data files
% hitori_7x7.dzn
nsize = 7 ;
hitoritabl = [|
5, 1, 5, 7, 7, 3, 3 |
1, 3, 3, 2, 7, 5, 6 |
2, 2, 3, 3, 6, 7, 6 |
5, 7, 7, 6, 6, 2, 4 |
2, 7, 7, 4, 1, 6, 4 |
1, 6, 6, 1, 4, 4, 7 |
7, 4, 6, 3, 2, 4, 5 |] ;
% hitori_15x15.dzn
nsize = 15 ;
hitoritabl = [|
15, 3, 1, 4, 15, 13, 1, 6, 15, 5, 7, 2, 1, 7, 11 |
8, 13, 6, 13, 4, 14, 13, 7, 11, 12, 10, 1, 15, 5, 1 |
6, 1, 15, 8, 6, 15, 4, 6, 13, 14, 15, 3, 5, 6, 9 |
11, 12, 3, 12, 12, 1, 12, 15, 14, 12, 5, 12, 8, 2, 4 |
6, 8, 3, 1, 13, 6, 5, 9, 3, 15, 3, 4, 10, 11, 9 |
10, 6, 11, 10, 7, 9, 5, 5, 3, 6, 2, 15, 10, 4, 8 |
12, 11, 7, 2, 12, 6, 15, 4, 7, 9, 13, 7, 3, 12, 10 |
3, 11, 4, 11, 5, 11, 9, 11, 6, 2, 11, 14, 2, 10, 11 |
6, 5, 6, 9, 6, 12, 6, 14, 6, 10, 15, 6, 11, 6, 3 |
4, 14, 2, 13, 15, 5, 11, 3, 7, 8, 1, 10, 5, 9, 3 |
11, 2, 14, 3, 8, 10, 11, 1, 11, 7, 3, 11, 13, 3, 6 |
2, 4, 12, 15, 1, 12, 6, 12, 8, 12, 11, 5, 12, 14, 12 |
5, 13, 12, 14, 13, 8, 7, 11, 13, 4, 9, 13, 2, 13, 15 |
14, 15, 5, 12, 9, 4, 12, 10, 2, 11, 1, 13, 1, 12, 7 |
13, 7, 13, 6, 4, 3, 2, 4, 10, 13, 14, 3, 4, 15, 4 |] ;
and finally my MiniZinc model
%% hitori.mzn
% Given square a grid of sidelength n
% The grid is filled with numbers 1..n, some of them repeated in rows and columns
% Goal: Cover some of the repeated numbers so that
% i) there are no duplicates in any row of column
% ii) no covered square is vertically of horizontally adjacent to another coverde square
% iii) the uncovered blocks are connected
% data instance defined in *.dzn
int: nsize ;
array[SDX,SDX] of SDX: hitoritabl ;
int : maxinputflow = nsize*nsize ;
var 1..2 : toprowstartcol ;
var int : sumflowremain ;
var int: inputflow = sum(hitorimask) ;
set of int: SDX = 1..nsize ;
set of int: SDXmin1 = 1..nsize-1 ; % index set over pairs of adjacent cells
set of int: SDX0 = 0..nsize ;
set of int: SDX2 = 2..nsize ;
% primary unknowns, 0 or 1 whether a square is covered or open
array[SDX,SDX] of var 0..1 : hitorimask ;
% modelling variable to display the solution, 0 is used for a covered/black cell
array[SDX,SDX] of var SDX0: hitorisoln ;
constraint forall( ir in SDX, jc in SDX) ( hitorisoln[ir,jc] = hitorimask[ir,jc] * hitoritabl[ir,jc] ) ;
% model constraint
% use primary variables to construct quantity on which nonrepeat constraint is imposed
include "alldifferent_except_0.mzn" ;
constraint forall ( ir in SDX) ( alldifferent_except_0 ( [ hitorisoln[ir,jc]| jc in SDX ] ) ) ;
constraint forall ( jc in SDX) ( alldifferent_except_0 ( [ hitorisoln[ir,jc]| ir in SDX ] ) ) ;
% impose condition on non-adjacent masked cells
% no black squares directly below another black square
constraint forall (ir in SDX2, jc in SDX) ( 0 < hitorimask[ir-1,jc] + hitorimask[ir,jc] ) ;
% no black square directly to the right of another black square
constraint forall (ir in SDX, jc in SDX2) ( 0 < hitorimask[ir,jc-1] + hitorimask[ir,jc] ) ;
% this constraint weeds out numerous sibgle isolated open cells
constraint forall (ir in SDX, jc in SDX)
( if hitorimask[ir,jc]=1 then
( if ir > 1 then hitorimask[ir-1,jc ] else 0 endif
+ if jc < nsize then hitorimask[ir ,jc+1] else 0 endif
+ if ir < nsize then hitorimask[ir+1,jc ] else 0 endif
+ if jc > 1 then hitorimask[ir ,jc-1] else 0 endif
> 0
)
endif
) ;
% the following implements the connectedness condition of the open cells
array[SDX,SDX] of var -maxinputflow..maxinputflow : xtop ; % inflows from top v
array[SDX,SDX] of var -maxinputflow..maxinputflow : xright ; % inflows from right <-
array[SDX,SDX] of var -maxinputflow..maxinputflow : xbottom ; % inflows from bottom ^
array[SDX,SDX] of var -maxinputflow..maxinputflow : xleft ; % inflows from left ->
array[SDX,SDX] of var 0..maxinputflow: flowremain ;
% initialize top, right, bottom and left rows and columns of xl, xr, xb and xl
constraint if hitorimask[1,1] = 0 % then necessarily hitorimask[1,2]=1
then toprowstartcol = 2 /\ xtop[1,2]=inputflow /\ xtop[1,1] = 0
else toprowstartcol = 1 /\ xtop[1,1]=inputflow
endif ;
constraint forall (jc in 3..nsize) (xtop[1,jc] = 0) ; % make zero the inflow from the top for the rest of the top row
constraint forall (ir in SDX) (xright[ir,nsize] = 0) ;
constraint forall (jc in SDX) (xbottom[nsize,jc] = 0) ;
constraint forall (ir in SDX) (xleft[ir,1] = 0) ;
% now make zero the inflow from all directions for all black cells
constraint forall (ir in SDX, jc in SDX)
( if hitorimask[ir, jc] = 0
then (xtop[ir,jc]=0) /\ (xright[ir,jc]=0) /\ (xbottom[ir,jc]=0) /\ (xleft[ir,jc]=0)
else true
endif
) ;
% flow balance across cell boundaries, first horizontally then vertically
constraint forall (ir in SDXmin1, jc in SDX) ( xbottom[ir,jc] = -xtop[ir+1,jc] ) ;
constraint forall (ir in SDX, jc in SDXmin1) ( xright[ir,jc] = -xleft[ir,jc+1] ) ;
% sum of flows into a cell
constraint forall (ir in SDX, jc in SDX)
(flowremain[ir,jc] = xtop[ir,jc] + xright[ir,jc] + xbottom[ir,jc] + xleft[ir,jc] ) ;
% objective function
var int: obj ;
constraint obj = sum (ir in SDX, jc in SDX) (hitorimask[ir,jc]=flowremain[ir,jc]) ;
solve maximize obj ;
output["\nsum(flowremain): \(sum(flowremain)) obj: \(obj) " ] ++ ["\n"] ;
output ["\n hitorimask hitorisoln flowremain \n" ] ;
output [
if jc <= nsize then "" ++ show_int(3, hitorimask[ir,jc] ) else "" endif
++ if jc = nsize then " " else "" endif
++ if nsize < jc /\ jc <= 2*nsize
then if fix(hitorisoln[ir,jc-nsize])>0
then show_int(3, hitorisoln[ir,jc-nsize] )
else " ."
endif
else ""
endif
++ if jc = 2*nsize then " " else "" endif
++ if 2*nsize < jc /\ jc <= 3*nsize
then show_int(3, flowremain[ir,jc-2*nsize] )
else ""
endif
++ if jc = 3*nsize then "\n" else "" endif
| ir in SDX, jc in 1..3*nsize ] ;
% the following displays first the four matrices xt, xr, xb and xl and then their sum flowremain
% set of int: SD5n = 1..5*nsize ; % index set along columns of printout
%output ["row xtop x <- xbottom x -> flowremain " ] ++ ["\n"] ;
%
%output [if jc=1 then show_int(3,ir) ++ " " else "" endif
% ++ if jc <= nsize then show_int(3,xtop[ir,jc]) else "" endif
%
% ++ if jc = nsize then " " else "" endif
% ++ if nsize < jc /\ jc <= 2*nsize then show_int(3, xright[ir,jc-nsize]) else "" endif
%
% ++ if jc = 2*nsize then " " else "" endif
% ++ if 2*nsize < jc /\ jc <= 3*nsize then show_int(3,xbottom[ir,jc-2*nsize]) else "" endif
%
% ++ if jc = 3*nsize then " " else "" endif
% ++ if 3*nsize < jc /\ jc <= 4*nsize then show_int(3,xleft[ir,jc-3*nsize]) else "" endif
%
% ++ if jc = 4*nsize then " " else "" endif
% ++ if 4*nsize < jc /\ jc <= 5*nsize then show_int(3,flowremain[ir,jc-4*nsize]) else "" endif
%
% ++ if jc = 5*nsize then "\n" else "" endif
% | ir in SDX, jc in SD5n ] ++ ["\n"] ;
Upvotes: 0
Reputation: 67
Your suggestion of flow into the network from the top-left corner works great!
I implemented it using four tables for the flows into a cell, from the top, the right, the bottom and from the cell to the left. These flows are constrained by zero flow through the outside boundaries, except for the single flow at the top left or second from left, by continuity across interior boundaries and zero flows through the boundaries of the black cells. I then “discovered” that the condition of a single unit of flow remaining in each uncovered cell and the zero flow condition for the covered cells correspond exactly to the mask I used to identify covered and uncovered cells. Then, by maximizing the number of cells where the remaining flow corresponds to whether the cells is covered or not, one finds the solution.
I tested my MiniZinc model on several instances, including a 15x15 hard one. Most of the solvers find the unique solution within half a second.
BTW, you must have quite some experience with network flow models to have realised that the connectedness condition could be modelled as network flow.
So, with great appreciation, thank you for your suggestions.
Upvotes: 0