Totally Spatially Balanced Latin Square

From CometPublic

Contents 1 The problem 2 The model 2.1 Avoiding differentiable invariants 2.1.1 A naive version 2.1.2 More elaborated version 3 Performances 4 Reference

[edit] The problem

A Latin Square of size n is a n × n matrix in which each number 1..n appears exactly once in each row and column. The Latin Square is said totally balanced if values are spread uniformly in this matrix. Expressing it more formally requires the following definitions :

The distance d_r(v,w) of a pair (v,w) in a row r is the absolute difference of the column indices in which v and w appear in row r.

The total distance d(v,w) of a pair (v,w) is the sum of the distances in every line :

.

Then, a Latin Square is totally spatially balanced if

Here is an example of a totally spatially balanced Latin square with size 5

[edit] The model

We use integer variables col[r,v] which represent the column index of value v inside row r.

In the example above, col[2,1]=5 because the element with value 1 in the 2nd row is the 5th one.

The constraints for Latin Square are then

forall(v in R)
    S.post(alldifferent(all(r in R) col[r,v]));

The totally spatially balanced objective is

It will be implemented as a differentiable invariant

And the final objective is to minimize this and the violations of Latin square constraints :

We use a Tabu search.

The complete code is :

include "LocalSolver";

LocalSolver m();

int n=5;
range R = 1..n;
// balance is the optimal distance between 2 values
int balance = n*(n+1)/3;

// The variables of the model are
// col[r,v], the column index of value v in row r
var{int} col[R,R](m,R);

// for each r, initialize col[r,*] as a random permutation of 1..n
forall(r in R) {
    RandomPermutation p(R);
    forall(v in R) col[r,v] := p.get();
}

// The Latin Square constraints
ConstraintSystem S(m);
forall(v in R)
    S.post(alldifferent(all(r in R) col[r,v]));

// The totally spatially balanced objective
ObjectiveSum OS(m);
forall(v in R, w in R: v < w)
    OS.post(((sum(r in R) (abs(col[r,v] - col[r,w])))-balance)^2);

// Total objective is a combination of both (arbitrary coefficients)
Objective O = n * S + OS;

m.close();


// Search part : tabu search
int tabu[R,R,R] = -1;
int tabuLength = 10;
int it = 0;
while (O.evaluation() > 0) {
    // Differentiable invariants provide a method getSwapDelta
    selectMin(r in R,v in R, w in R: v < w && tabu[r,v,w] <= it)
        (O.getSwapDelta(col[r,v],col[r,w])) {
            col[r,v] :=: col[r,w];
            tabu[r,v,w] = it + tabuLength;
        }
    it++;
}

// Print the solution found
forall ( r in R ) {
    int row[R];
    forall ( v in R )
        row[col[r,v]]=v;
    forall ( i in R )
        cout << " " << row[i];
    cout << endl;
}

[edit] Avoiding differentiable invariants

The use of differentiable invariants (ObjectiveSum OS in the code) could be avoided by using simple invariants and providing a custom implementation of getSwapDelta.

[edit] A naive version

A first and naive version consists in remplacing the two Objective instances O and OS by invariants:

// totally spatially balanced
var{int} OS(m) <- sum(v in R, w in R: v < w)( ((sum(r in R) (abs(col[r,v] - col[r,w])))-balance)^2 );
// final objective 
var{int} O(m) <- n*S.violations() + OS;

We implement the function getSwapDelta by simply swapping the values, getting the value of O, swapping back and getting the difference of values O

function int getSwapDeltaNaive(var{int} O, var{int} v, var{int} w) {
 v:=:w;
 int futureVal = O;
 v:=:w;
 return futureVal - O;
}

Note : A better way would be to use lookahead

And the O.getSwapDelta(col[r,v],col[r,w]) must be remplaced by getSwapDeltaNaive(O,col[r,v],col[r,w]) in the selectMin of the search

[edit] More elaborated version

A second (and more elaborated) method consists in defining differentiable invariants that keep the increase of the objective OS obtained when swapping some value in some row.

The objectives are not changed, except that OS is not anymore a differentiable invariant.

// The totally spatially balanced objective, as a simple invariant.
var{int} OS(m) <- sum(v in R, w in R: v < w)
      ( ((sum(r in R) (abs(col[r,v] - col[r,w])))-balance)^2 );

// Total objective, as a simple invariant too.
var{int} O(m) <- n*S.violations() + OS;

We introduce many invariants to compute the change in OS when performing a swap.

// d_r[r,v,w] is an invariant storing the distance between values v and w in
// row r (written d_r(v,w) in the referenced article)
var{int} d_r[R,R,R](m,R);
forall(r in R, v in R, w in R : v < w)
    d_r[r,v,w] <- abs(col[r,v] - col[r,w]);
// exploit symmetry of distance
forall(r in R, v in R, w in R : v > w)
    d_r[r,v,w] = d_r[r,w,v];

// d[v,w] is an invariant storing the total distance between values v and w
// (written d(v,w) in the referenced article)
var{int} d[R,R](m,n..n*(n-1));
forall(v in R, w in R : v < w)
    d[v,w] <- sum(r in R) d_r[r,v,w];
// exploit symmetry of distance
forall(v in R, w in R : v > w)
    d[v,w] = d[w,v];

// delta_d[r0,v0,w0,w] is the increase of d(v0,w) when swapping v0 and w0
// in row r0. It is equal to the increase in d_r0(v0,w) when swapping v0 and
// w0 in row r0 (since rows other than r0 are not affected by that swap) :
// delta_d[r0,v0,w0,w] = d_r0(w0,w) - d_r0(v0,w)
var{int} delta_d[R,R,R,R](m,R);
forall(r0 in R, v0 in R, w0 in R, w in R : v0<w0)
    delta_d[r0,v0,w0,w] <- d_r[r0,w0,w] - d_r[r0,v0,w];
// no need to define it for v0>w0

// delta_OS[r0,v0,w0,w] is the increase of the objective function OS when
// swapping v0 and w0 in row r0. The computation details are a bit long, but
// simple and will not be presented here.
var{int} delta_OS[R,R,R](m);
forall(r0 in R, v0 in R, w0 in R : v0<w0)
    delta_OS[r0,v0,w0] <- sum(w in R : w!=w0 && w!=v0)
            ( delta_d[r0,v0,w0,w] * (delta_d[r0,v0,w0,w] + 2*d[v0,w] - 2*balance) )
        + sum(v in R : v!=v0 && v!=w0)
            ( -delta_d[r0,v0,w0,v] * (-delta_d[r0,v0,w0,v] + 2*d[v,w0] - 2*balance) );

With these, the original call to O.getSwapDelta(col[r,v],col[r,w]) can be replaced with

delta_OS[r,v,w] + n * S.getSwapDelta(col[r,v],col[r,w])

[edit] Performances

The following graphs compare the performances of each method.

• Running time :

Image:RunningTimes.png

• Number of iterations :

Image:NbIterations.png

• Time spent at each iteration (dividing the total time by the number of iterations) :

Image:TimePerIteration.png

[edit] Reference

Van Hentenryck and Michel, Differentiable invariants, In Proceedings of CP2006, Springer.