start = time.time()
N = 9
cons = cons_big
#set up the solver and variables.
s = Tactic("qffd").solver()
# H = hooks
H = np.array(IntVector("h", N ** 2), dtype=object).reshape((N,N))
row_fix =IntVector("r", N)
col_fix =IntVector("c", N)
s += [And(e > 0 , e <= N) for (i,j), e in np.ndenumerate(H)] # each hook is an integer in 1...n
s += [And(e >= 0, e < N) for e in row_fix + col_fix] # each row_fix and col_fix varibale is an integer 1..n
s += Distinct(row_fix) # The row_fix and col_fix variables are distinct
s += Distinct(col_fix)
for n in range(1, N+1): # For each of the hooks 1...n
s += PbEq([(e == n, 1) for _ , e in np.ndenumerate(H)], 2 * n - 1)
if n != 1:
s+=PbEq([(If(e==int(n),Sum([H[k,l] == int(n) for k,l in neigh(i,j,N)]),0)==1,1)for (i,j),e in np.ndenumerate(H)],2)
for n in range(N):
s += [Implies(e == int(n+1),Or(i == row_fix[n],j == col_fix[n])) for (i,j),e in np.ndenumerate(H)]
# X = nums
X = np.array(IntVector("x",N**2),dtype=object).reshape((N,N))
hook_num = IntVector("n",N)
s += [And(e>0,e<N+1) for e in hook_num] # hook nums range from 1...N
s += [And(e>0,e<N+1) for (i,j),e in np.ndenumerate(X)] # Each non-zero entry ranges from 1...N
s += Distinct(hook_num) # The hook numbers must be distinct/ one of each
for n in range(N):
s += [Implies(e==n+1,X[i,j]==hook_num[n]) for (i,j),e in np.ndenumerate(H)] # X entries must match their hook number
# Y = filled
Y = np.array(IntVector("y",N**2),dtype=object).reshape((N,N)) # define Y va
s += [e == 0 if (i,j) in cons.keys() else Or(e==0,e==1) for (i,j),e in np.ndenumerate(Y)]
for n in range(1,N+1):
s+=PbEq([(And(e==n,Y[i,j]==1),1) for (i,j),e in np.ndenumerate(X)],n)
for (i,j), v in cons.items():
s += [sum(X[a,b] * Y[a,b] for a,b in neigh(i,j,N)) == v]
###
# no filled 2x2 Constraint
s += [Or(Y[i,j] ==0,Y[i+1,j] ==0,Y[i,j+1] ==0,Y[i+1,j+1] ==0) for j in range(N-1) for i in range(N-1)]
################### Set up connectivity constraints ################################################################
edge ={}
Z = np.array(IntVector('z', N*N),dtype=object).reshape(N,N)
# create edges and variable to define order in connectivity tree
for i in range(N):
for j in range(N):
for (k,l) in neigh(i,j,N):
edge[(i,j,k,l)] = Int("e%d%d%d%d" % (i,j,k,l))
s+=Or(edge[(i,j,k,l)] ==0,edge[(i,j,k,l)] ==1)
# no edges into or out of unfilled cells. Only one edge per cell
for i in range(N):
for j in range(N):
s += Implies(Y[i,j] ==0 , Sum([edge[(k,l,i,j)] for (k,l) in neigh(i,j,N)]) == 0)
s += Implies(Y[i,j] ==0 , Sum([edge[(i,j,k,l)] for (k,l) in neigh(i,j,N)]) == 0)
s += Implies(Y[i,j] > 0 , Sum([edge[(k,l,i,j)] for (k,l) in neigh(i,j,N)]) <= 1)
for (k,l) in neigh(i,j,N):
# Make edges one way to form a tree
s += (edge[(i,j,k,l)] + edge[(k,l,i,j)]) <= 1
# limit the order values
s += [And(e>= 0 ,e <=int(N*N)) for (i,j),e in np.ndenumerate(Z)]
# order is ascending along edges
for i in range(N):
for j in range(N):
s += [Implies(And(Y[i,j]!=0,edge[(k,l,i,j)]==1),Z[i,j] > Z[k,l]) for (k,l) in neigh(i,j,N)]
# only one cell with no feed in => root
s += PbEq([(And(Y[i,j] != 0,Sum([edge[(k,l,i,j)] for (k,l) in neigh(i,j,N) ]) == 0),1) for i in range(N) for j in range(N)],1)
################### End of connectivity constraints ################################################################
print("setup done in {:.2f} secs".format(time.time()-start))
stop = False
iteration = 0
while stop == False:
iteration += 1
if iteration % 1000 == 0:
print(f'On iteration {iteration}')
if s.check() == sat:
m = s.model()
evalu = np.vectorize(lambda x:m.evaluate(x).as_long())
result_x = evalu(X) # The filled numbers
result_y = evalu(Y) # The filled entries
if test_solution(result_x,result_y, cons):
grid_print(result_y,result_x,cons)
print("Took {:.4f} seconds".format(time.time()- start))
soln = result_x*result_y;
print("solution is {:,.0f}".format(areas(soln)))
print("Connectivity: ",np.max(label(soln>0,connectivity=1)) == 1)
stop=True
else:
s+=Or([X[i,j]*Y[i,j] != int(result_x[i,j]*result_y[i,j]) for (i,j),e in np.ndenumerate(X)])
