Costas Array

This is an example of solving the costas array problem with Numberjack.
The nxn Boolean matrix from the description below is represented as a set of n variables encoding a permuation, much like for the n-Queens problem.

Description:

(From Wikipedia:) A Costas array can be regarded geometrically as a set of n points lying on the squares of a nÃ—n checkerboard, such that each row or column contains only one point, and that all of the n(n-1)/2 displacement vectors between each pair of dots are distinct. This result in an ideal 'thumbtack' auto-ambiguity function, making the arrays useful in applications such as Sonar and Radar. A Costas array may be represented numerically as an nÃ—n array of numbers, where each entry is either 1, for a point, or 0, for the absence of a point. When interpreted as binary matrices, these arrays of numbers have the property that, since each row and column has the constraint that it only has one point on it, they are therefore also permutation matrices. Thus, the Costas arrays for any given n are a subset of the permutation matrices of order n. Costas arrays can be regarded as two-dimensional cousins of the one-dimensional Golomb ruler construction, and, as well as being of mathematical interest, have similar applications in experimental design and phased array radar engineering. All Costas arrays of size up to and including 26x26 are known. Costas arrays are known for infinitely many sizes because there exist two methods of generating them, both of which work near primes (of which there are infinitely many) and powers of primes. These are known as the Welch and Lempel-Golomb generation methods, and come from a branch of mathematics known as finite field theory. It is not known whether Costas arrays exist for all sizes. Currently, the smallest sizes for which no arrays are known are 32x32 and 33x33.

Code:

from Numberjack import *

def model_costas_array(N):
sequence = VarArray(N,1,N)
model = Model(
AllDiff(sequence),
[AllDiff([sequence[j] - sequence[j+i+1] for j in range(N-i-1)]) for i in range(N-2)]
)
return(sequence,model)

def solve_costas_array(param):
(sequence,model) = model_costas_array(param['N'])
solver = model.load(param['solver'])
if solver.solve():
print_costas_triangle(sequence)
print 'Nodes:', solver.getNodes(), ' Time:', solver.getTime()

def print_costas_triangle(seq):
N = len(seq)
print ''.join([str(int(var.get_value())).rjust(3) for var in seq])
for i in range(N-1):
print ''.join([str(int(seq[j].get_value())-int(seq[j+i+1].get_value())).rjust(3)
for j in range(N-i-1)])

solve_costas_array(input({'solver':'Mistral','N':10}))