Arc Consistency in J

NB. ----------------------------------------------------------------------------
NB. Arc consistency for constraint satisfaction problems (CSP).
NB. ----------------------------------------------------------------------------
NB. Michal Dobrogost 2012-11-03
NB.
NB. A CSP can model various satisfaction problems such as sudoku puzzles.
NB. A CSP consists of a set of variables X, each with an associated domain D of
NB. values. In sudoku, X = {1,2,...,81} and for each x in X, Dx = {1,2,...,9}.
NB. There are also a set of constraints C between pairs of variables.
NB. In sudoku, Cxy is simply the 'not equal' constraint and applies to all pairs
NB. of variables in rows, columns, and 3x3 squares.
NB.
NB. Arc consistency (AC) is a filtering algorithm that removes values, from the
NB. domain of a variable, that can not participate in any solution. In sudoku,
NB. if a variable x1 has domain {1,2,6} and is on the same row as x2 with domain
NB. {6} we can remove the value 6 from the domain of x1 because there can not
NB. be a solution where Dx1 = Dx2 = {6}.
NB.
NB. Filtering algorithms are usually combined with backtracking search which
NB. guesses a specific value for a variable with multiple available values.
NB. Guesses are interleaved with filtering (commonly refered to as maintaining
NB. arc consistency (MAC) search).
NB.
NB. Thanks to:
NB.   Raul M., Mike D., programming@jsoftware.com for many suggestions.

load 'format/printf'

NB. ----------------------------------------------------------------------------
NB. Generate a random CSP instance and filter it.

NB. Constants
n=. 4 NB. # variables
d=. 4 NB. domain size

NB.* 'D' is the domain of each variable, one row per variable.  We generate the
NB.*     domains so that each one is half full.  An empty domain would signify
NB.*     that the instance is not consistent to begin with so we avoid that.
] D=. (n%2) > (?]) n$d

NB. a single constraint (note: constraints may be asymmetric, for example x<y).
c1=. ? (d,d)$2
2 2 $ '' ; 'Dx' ; 'Dy' ; c1

NB.* 'C' is all of our constraints with the 0'th entry being the empty.
] C=. a: , c1 ; (|:c1)

NB.* 'A' an adjacency matrix (which constraint is between a pair of variables).
A=. ? (2$n)$2          NB. generate random matrix of [0,1]
A=. A *. (i.n) </ i.n  NB. make it upper diagonal, zeros on diagonal
] A=. A + |: 2*A       NB. make it symmetrix referencing transpose in C.
2 2 $ '' ; 'y' ; 'x' ; A

NB. ----------------------------------------------------------------------------
NB. Utility functions.

NB.* 'adj A' computes adjacency list from an adjacency matrix.
adj=: < @ I. @: *
NB.* 'a' is an adjacency list: which variables are adjacent to a given variable.
] a=. adj A

NB.* 'xs arcsX A' computes a list of all adjacent variables to those in xs.
NB.* 'ys arcsY A' computes a list of all adjacent variables to those in ys.
arcsX=: [  ((#~ #@>) ,. ;@]) {
arcsY=: [  (;@] ,. (#~ #@>)) {
] (i.n) arcsX a
] (i.>.n%2) arcsY a

NB. Mike D.'s compute arcs directly from A
NB. arcsA=: ((]#~ (e.~ {."1))   ($#: I.@,)) *

NB.* 'c revDom (Dx,Dy)' will return the domain of y supported by x.
getx=: 0{ ]
gety=: 1{ ]
revDom=: getx *. +./ @ (gety # [)

NB.* 'isValid D' decides if all domains are non-empty.
isValid=: <./ @: (>./"1)

NB.* 'allAssigned D' is true if all domains have exactly one value.
allAssigned=: 1&= @ */ @: (+/"1)

NB. ----------------------------------------------------------------------------
NB. The arc consistency algorithm itself.

NB.* '(A;a;<C) revise (ys;D)' filter domains of all variables (xs) adjacent to
NB.*                          ys. Returns (newXs;newD) where newXs are those
NB.*                          variables that are different across newD and D.
revise=: 4 : 0
    'A a C'=. x
    'ys D'=. y

    if. 0 < # ys
    do.
        arcs=. ys arcsY a
        ax=. 0{"1 arcs
        xs=. ~. ax
        if. 0 = # xs do.
            (0#0);D
        else.


            NB. revLookup [x,y]  <=>  Cxy revDom (Dx;Dy)
            revLookup=. (> @ {&C @ {::&A) revDom ({&D)

            NB. produce modified domains and the variables they correspond to.
            newD=. ax *.//. revLookup"1 arcs

            ((newD ([: >./"1 ~:) xs{D) # xs) ; (newD xs} D)
        end.
    else.
        (0$0) ; D
    end.
)

NB.* '(A;a;<C) ac (ys;D)' return D filtered into an arc consistent state given
NB.*                      that the the variables in ys have been changed.
NB.*                      Use (i.#D);D as right argument to filter from scratch.
ac=: 4 : 0
    > (1&{) ((x&revise)^:_) y
)

NB. ----------------------------------------------------------------------------
NB. Filter the random instance.

NB. Recall the random instance.
(A;a;<C)
D

NB. Filter the random instance.
(A;a;<C) ac ((i.#D);D)

NB. Are all the filtered domains non-empty?
isValid ((A;a;<C) ac ((i.#D);D))


NB. ----------------------------------------------------------------------------
NB. Search

search=: 4 : 0
    if. allAssigned >1{y do.
        > 1{y
    else.
        D=. x ac y

        if. -. isValid D do.
            D
        else.
            if. allAssigned D do.
                D
            else.
                nonEmpties=. 1< +/"1
                var=. {. I. nonEmpties D
                ret=. ($D) $ 0
                for_i. I. var{D do.
                    D2=. (1 i} ((1{$D) # 0)) var} D
                    D2=. x search ((1#var) ; D2)
                    if. allAssigned D2 do.
                        ret =. D2
                        break.
                    end.
                end.
                ret
            end.
        end.
    end.
)

(A;a;<C) search (i.#D);D


NB. ----------------------------------------------------------------------------
NB. Raul M.'s brute force sudoku instance. (doesn't work for me).

NB. C=: , |: 9#,:i.9 NB. identity within columns (C)
NB. R=: ,9#,:i.9 NB. identity within rows (R)
NB. B=: ,3#3#"1 i.3 3 NB. identity within boxes
NB. A=: ((2*>/~i.81)+</~i.81) * (+. |:) (C =/ R) +. (C =/B ) +. (R =/ B)
NB. D=. 81 9 $ 1
NB.
NB. (A;(adj A);<C) search (i.#D);D