NB. ----------------------------------------------------------------------------
NB. Arc consistency for constraint satisfaction problems (CSP).
NB. ----------------------------------------------------------------------------
NB. Michal Dobrogost 2012-11-21
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. ----------------------------------------------------------------------------
NB. Globals

NB. the number of digits in the sudoku (usually 9) has to be a perfect square.
N=: 4

NB. The ratio of set cells in the randomly generated sudoku.
R=: 0.3

NB. Set a specific RNG state for nice results... remove this for randomness.
(9!:43) 1
(9!:1) 598354

NB. ----------------------------------------------------------------------------
NB. Misc helpers

NB.* 'allZero X' returns 1 if all cells of X are 0.
allZero=: 1 - (>./ @ ,)

NB.* 'X xor Y' is the exclusive or of X and Y.  The result has the rank of the
NB.*           ranks of the arguments added together.
xor=: 2 | 1&+

NB.* 'N nOnes X' creates a vector of N 0's with the X'th 0-based entry set to 1.
NB.*             This is somewhat like the inverse of I. except it supports only
NB.*             a single index in X.
nOnes=: [ $!.0 (}. @ (,&1) @ #&0 @ ])

NB. ----------------------------------------------------------------------------
NB. Sudoku specific stuff (ie. generate a sudoku and convert it into a CSP).
NB. Note that the generated sudoku puzzles may be inconsistent (ie. have no
NB. solutions).

NB.* 'N genSudoku R' generates a suduko over 1..N, with R*N^2 set entries.
genSudoku=: 13 : '((2#x) $ ((y*x^2)> (?]) x^2)) * 1 + x ? (x#x)'
] S=: N genSudoku R

NB.* 'sudokuToD S' converts a sudoku puzzle S into an AC set of domains D.
sudokuToD=: (xor allZero)"1 @: (# (nOnes"0) ,)
D=: sudokuToD S
9 {. D NB. show the domains for the first row of the sudoku.


NB.* 'dToN D' converts each domain d of D into a visual representation.
NB.*          If d is assigned then the single number of the domain is shown.
NB.*          If d has multiple possible values or none, 0 is displayed.
dToN=: {. @ I. @ ((1 ~: +/) , ])"1

NB.* 'dToSudoku D' is the inverse of sudokuToD creating a sudoku representation
NB.*               from the given domains.
dToSudoku=: (%: @ #) ((2#[) $ [:dToN]) ]
dToSudoku D

NB.* 'sudokuToC S' generates the constraint tables to use (we need just the
NB.*               not-equal table in addition to the empty table at index 0).
sudokuToC=: (a:,<) @ (~:/)~ @ (i. @ #)
] C=: sudokuToC S


NB.* 'sudokuToArcs S' generates the arcs between variables who have a not-equal
NB.*                  constraint between them.  A row of the output will have
NB.*                  x y 1, where x is one variable, y another, and 1 is the
NB.*                  entry in C to use.
sudokuToArcs=: 3 : 0
    N=: #y  NB. # of digits in the sudoku
    B=: %:N NB. length of a side of a box

    NB.* 'n toRowCol x' works on a sudoku of n digits, converts variable x to
    NB.*                it's row and column representation.
    toRowCol=: (2#N) #: ]

    NB.* '(rx,cx) isAdj (ry,cy)' returns true if there is a not-equal constraint
    NB.*                         for two variables given in row-column notation.
    isAdj=: (-.@-:) *. (+./ @: =) +. ((B <.@%~ [) -: (B <.@%~ ]))

    NB. The adjacency matrix for the variables.
    X=: (isAdj"1)/~ toRowCol i.N^2

    NB. The arcs (ie. 1 entries in X).
    Arcs=: ,/ (i.#X) (,."0 1) I. X

    NB. Keep only one direction of each arc (eg. 0 1 not 1 0 )and append 1
    NB. indicating use the 1'th entry of C.
    (,.&1) (#~(</"1)) Arcs
)
15 {. A=: sudokuToArcs S

NB. ----------------------------------------------------------------------------
NB. Arc consistency stuff.

NB.* '(Dx,Dy) relevantCRows c' gets the rows from c specified by Dx and the
NB.*                           columns specified by Dy.
relevantCRows=: (#"1 2) ((,:|:)@:])

NB.* '(Dx,Dy) revDom c' compute (newDx,newDy) where the new versions are a subset
NB.*                    of the originals with only the mutually supported values.
revDom=: [ *. (1 0&{ @: (+./"2) @: relevantCRows)

NB.* '(A;C) revise 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 C'=. x
    'D'=. y
    getCIndex=: 2&{
    getVars=: 0 1&{

    ] RevDs=: ,/ (({&D @ getVars) revDom (> @ {&C @ getCIndex))"1 A
    ] RevAs=: , (getVars"1) A

    (RevAs *//. RevDs) (~.RevAs)} D
)

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
    ((x&revise)^:_) y
)

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. Run arc consistency on sudoku

S
dToSudoku (A;<C) revise D
] DSolved=: dToSudoku (A;<C) ac D

isValid DSolved
allAssigned DSolved