## A monadic constraint solver in F#

Monads, a.k.a. computation expressions, permits the programmer to define a custom domain specific language (dsl) within the host language. When thinking about this I got the, in my opinion, great idea to build a dsl for constraint programming in F# using monads. However, as is common with great ideas, someone else has already though of it before. In this case I found out that David Overton was the man. He uses a specialized monad that allows constraint programming within Haskell. I definitely recommend reading his excellent blog post on this and also how he uses his library to solve Sudoku puzzles. Still, I thought it would be interesting to see how this worked out in F#.

Using the state transformer from an earlier post, with a few minor modifications, provided the core framework that I needed. David also uses a few monad control functions that made his code short and readable. Mostly due to the fact that they cannot be implemented generically, I think no F# equivalent exists in any library. However, non-generic versions are easy to implement – and that is what I did. Here is the entire monadic framework:

// ref: http://www.randomhacks.net/articles/2007/03/12/monads-in-15-minutes

type ListMonad() =

// Semantics

member o.Return(x) = [x]

member o.Bind(m, f) = List.concat( List.map f m )

member o.Zero() = []

// Helpers

member o.Guard(c) = if c then o.Return(()) else o.Zero()

// ref: http://matthewmanela.com/blog/parameterized-state-transformer-monad-in-f-2/

// ref: http://monads.haskell.cz/html/xformeranatomy.html

// ref: http://cs.hubfs.net/forums/thread/7472.aspx (combine, delay)

// ref: http://tomasp.net/articles/imperative-i-return.aspx (combine, delay)

type StateMonadT(listM : ListMonad) =

// Semantics

member o.Return(v) = fun st -> listM.Return(v, st)

member o.ReturnFrom(x) = x

member o.Bind(m, f) = fun st -> listM.Bind(m st, fun (v, st’) -> f v st’)

member o.Zero() = fun st -> listM.Return((), st)

member o.Combine(p1,dp2) = o.Bind(p1, fun _ -> dp2)

member b.Delay(f) = f()

// Helpers

member o.Update(f) = fun st -> listM.Return((), f st)

member o.Get() = fun st -> listM.Return(st, st)

member o.Set(st) = fun _ -> listM.Return((), st)

member o.Execute(m, s0) = m s0 |> List.map fst

member o.Lift(c) = fun st -> listM.Bind(c, fun x -> listM.Return(x,st))

member o.Guard(c) = o.Lift(listM.Guard(c))

let listM = ListMonad()

let stateMT = StateMonadT(listM)

let rec mapM m xs = stateMT {

match xs with

| x::rest ->

let! value = m x

let! values = mapM m rest

return value :: values

| _ -> return []

}

let rec mapM_ m xs = stateMT {

match xs with

| x::rest ->

do! m x

do! mapM_ m rest

| _ -> return ()

}

let rec replicateM n m = stateMT {

match n with

| 0 -> return []

| _ ->

let! x = m

let! rest = replicateM (n-1) m

return x::rest

}

let rec zipWithM_ f xs ys = stateMT {

match xs, ys with

| x::rxs, y::rys ->

let m = f x y

do! m

do! zipWithM_ f rxs rys

| _ -> return ()

}

Having this in place the constraint solver can be constructed. First of all the types needed:

type FDVar = { id : int }

and VarInfo = { values : Set<int>; delayedConstraints : FDState -> List<unit*FDState> }

and FDState = { varMap : Map<FDVar, VarInfo>; nextId : int }

I kept it a little less verbose then the original but the semantics is the same. The function for executing a constraint program was short:

let runFD fd = stateMT.Execute(fd, { varMap = Map.empty; nextId = 0 })

The next step was to add functions for declarations of variables:

let newVar domain =

let nextVar = stateMT {

let! s = stateMT.Get()

let vid = s.nextId

do! stateMT.Set({ s with nextId = vid + 1})

return { id = vid }

}

let isOneOf x domain = stateMT {

return! stateMT.Update(fun s ->

let vm = s.varMap

let vi = { values = Set.ofList domain; delayedConstraints = stateMT.Return(()) }

{ s with varMap = Map.add x vi vm }

)

}

stateMT {

let! v = nextVar

do! isOneOf v domain

return v

}

let rec newVars n domain = replicateM n (newVar domain)

A few helper methods (following the same path as David):

let lookup x = stateMT {

let! s = stateMT.Get()

return s.varMap.[x].values

}

let update x i = stateMT {

let! s = stateMT.Get()

let vm = s.varMap

let vi = vm.[x]

do! stateMT.Set({ s with varMap = Map.add x { vi with values = i } vm })

do! vi.delayedConstraints

}

let addConstraint x constr = stateMT {

let! s = stateMT.Get()

let vm = s.varMap

let vi = vm.[x]

let delayedConstraints’ = stateMT {

do! constr

do! vi.delayedConstraints

}

do! stateMT.Set({ s with varMap = Map.add x { vi with delayedConstraints = delayedConstraints’ } vm })

}

let addBinaryConstraint f x y = stateMT {

let constr = f x y

do! constr

do! addConstraint x constr

do! addConstraint y constr

}

Now, using the helper methods, we start defining constraints:

let hasValue var value = stateMT {

let! vals = lookup var

do! stateMT.Guard(Set.contains value vals)

let i = Set.singleton value

if i <> vals then

do! update var i

}

let same =

let sameConstraint x y = stateMT {

let! xv = lookup x

let! yv = lookup y

let i = Set.intersect xv yv

do! stateMT.Guard(Set.isEmpty i |> not)

if i <> xv then

do! update x i

if i <> yv then

do! update y i

}

addBinaryConstraint sameConstraint

let different =

let differentConstraint x y = stateMT {

let! xv = lookup x

let! yv = lookup y

do! stateMT.Guard(Set.count xv > 1 || Set.count yv > 1 || xv <> yv)

if Set.count xv = 1 && Set.isProperSubset xv yv then

do! update y (Set.difference yv xv)

if Set.count yv = 1 && Set.isProperSubset yv xv then

do! update x (Set.difference xv yv)

}

addBinaryConstraint differentConstraint

let rec allDifferent xs =

stateMT {

match xs with

| x::rest ->

do! mapM_ (different x) rest

do! allDifferent rest

| _ -> return ()

}

I use the optimized different constraint that David suggests in his later post. The only thing that remains now is the labelling function for backtrack search.

let rec labelling vars =

let label var = stateMT {

let! vals = lookup var

let! value = stateMT.Lift(Set.toList vals)

do! hasValue var value

return value

}

mapM label vars

That’s all we need. A couple of samples may be in place:

let sample1 = stateMT {

let! a = newVar [1..4]

let! b = newVar [1..4]

do! different a b

return! labelling [a; b]

}

runFD sample1 |> printf “sample1: %A\n”

let sample2 = stateMT {

let! vars = newVars 3 [1..4]

do! allDifferent vars

return! labelling vars

}

runFD sample2 |> printf “sample2: %A\n”

let sample3 = stateMT {

let! a = newVar [1..9]

let! b = newVar [1..3]

let! c = newVar [1..2]

do! hasValue a 7

return! labelling [a;b;c]

}

runFD sample3 |> printf “sample3: %A\n”

let sample4 = stateMT {

let! a = newVar [1..9]

let! b = newVar [2..7]

do! same a b

return! labelling [a;b]

}

runFD sample4 |> printf “sample4: %A\n”

And of course also a Sudoku solver:

type Puzzle = int list

let test : Puzzle = [

0; 0; 0; 0; 8; 0; 0; 0; 0;

0; 0; 0; 1; 0; 6; 5; 0; 7;

4; 0; 2; 7; 0; 0; 0; 0; 0;

0; 8; 0; 3; 0; 0; 1; 0; 0;

0; 0; 3; 0; 0; 0; 8; 0; 0;

0; 0; 5; 0; 0; 9; 0; 7; 0;

0; 5; 0; 0; 0; 8; 0; 0; 6;

3; 0; 1; 2; 0; 4; 0; 0; 0;

0; 0; 6; 0; 1; 0; 0; 0; 0 ];

// ref: http://fsharpcode.blogspot.com/2010/03/splitevery-chunk.html

let rec chunk n xs =

let splitAt n (xs:#seq<‘a>) = (Seq.take n xs |> Seq.toList, Seq.skip n xs |> Seq.toList)

match xs with

| [] -> []

| _ ->

let ys, zs = splitAt n xs

ys :: chunk n zs

let rows xs = chunk 9 xs

let columns xs =

// ref: http://stackoverflow.com/questions/3016139/help-me-to-explain-the-f-matrix-transpose-function

let rec transpose = function

| (_::_)::_ as M -> List.map List.head M :: transpose (List.map List.tail M)

| _ -> []

xs |> rows |> transpose

let boxes xs =

let box (xa : array<‘a>) bi bj : list<‘a> =

[

for j in 0..2 do

let r = bj*3 + j

for i in 0..2 do

let c = bi*3 + i

yield xa.[r*9+c]

]

[

let xa = Array.ofList xs

for bj in 0..2 do

for bi in 0..2 do

yield box xa bi bj

]

let sudoku (puzzle : Puzzle) = stateMT {

let! vars = newVars 81 [1..9]

do! zipWithM_ (fun x n -> stateMT { if n > 0 then do! hasValue x n }) vars puzzle

do! mapM_ allDifferent (rows vars)

do! mapM_ allDifferent (columns vars)

do! mapM_ allDifferent (boxes vars)

return! labelling vars

}

let result = runFD (sudoku test)

for r in result do

r |> rows |> List.map (fun row -> printf “%A\n” row) |> ignore

printf “———–\n”

## Poor mans state monad transformer in F#

It seems like the type system of F# (and .NET) is currently not powerful enough to capture the semantics of (generic) monad transformers. The issue has been up for discussion at hubFS. One alternative is to give up generic and make explicit instantiation for the inner monads of interest. I think this approach is used by The F# Monad Library. Taking this approach for combining the State monad and the List monad has also been publish by Matthew Manela in his blog. I added methods for lift and execute and ended up having something like this:

//ref: http://www.randomhacks.net/articles/2007/03/12/monads-in-15-minutes

type ListMonad() =

// Semantics

member o.Return(x) = [x]

member o.Bind( (m:’a list), (f: ‘a -> ‘b list) ) = List.concat( List.map f m )

member o.Zero() = []

// Helpers

member o.Guard(c) = if c then o.Return(()) else o.Zero()

// ref: http://matthewmanela.com/blog/parameterized-state-transformer-monad-in-f-2/

// ref: http://monads.haskell.cz/html/xformeranatomy.html

type StateMonadT(listM : ListMonad) =

// Semantics

member o.Return(v) = fun st -> listM.Return(v, st)

member o.ReturnFrom(x) = x

member o.Bind(m, f) =

fun st -> listM.Bind(m st, fun (v, st’) -> f v st’)

member o.Zero() = fun st -> listM.Zero()

// Helpers

member o.Update(f) = fun st -> listM.Return(st, f st)

member o.Get() = fun st -> listM.Return(st, st)

member o.Set(st) = fun _ -> listM.Return((), st)

member o.Execute(m, s0) = m s0 |> List.map fst

member o.Lift(c) = fun st -> listM.Bind(c, fun x -> listM.Return(x,st))

let listM = ListMonad()

let stateMT = StateMonadT(listM)

I used the state monad transformer to write an n-queen program.

let rec checkHorizontal board q =

match board with

| x::xs -> q <> x && checkHorizontal xs q

| _ -> true

let checkDiagonals board q =

let rec checkDiagonals’ board q i =

match board with

| x::xs ->

let j = i + 1

q <> x + j && q <> x – j && checkDiagonals’ xs q j

| _ -> true

checkDiagonals’ (List.rev board) q 0

let validPositions board n = listM {

let! q = [1..n]

do! listM.Guard(checkHorizontal board q && checkDiagonals board q)

return q

}

#nowarn “40”

let rec addQueen = stateMT {

let! board, n = stateMT.Get()

if List.length board = n then

return board

else

let! position = stateMT.Lift(validPositions board n)

do! stateMT.Set(board @ [position], n)

return! addQueen

}

stateMT.Execute(addQueen, ([],4)) |> printf “solutions: %A\n”

Do not take this as the best way to solve the N-queen problem (in fact I had to discard a compiler complaint, no 40, about recursion being implicit or something) but only as a sample usage of the monad transformer.

## Constraint Programming

Wikipedia defines constraint programming as:

Constraint programming is a programming paradigm wherein relations between variables are stated in the form of constraints. Constraints differ from the common primitives of imperative programming languages in that they do not specify a step or sequence of steps to execute, but rather the properties of a solution to be found. This makes constraint programming a form of declarative programming.

The clpfd library for Prolog makes it possible to express constraint programs. This way of embedding constraints in a logic program is called constraint logic programming and clpfd is indeed a shorthand for “Constraint Logic Program over Finite Domains”. Finite domains refers to the fact that the variables are limited to a finite set of values.

Other libraries that support constraint programming exists for other languages. Microsoft provides one that is accessible from .net languages that is called Microsoft Solver Foundation. Robert Pickering has published a sample where he uses it to build a Sudoku solver. The host language is F# but of course C# or Visual Basic .NET would have worked out as well. Fasail has also done the same thing.

Writing a full featured, high performance constraint solver requires lots of work but building a basic yet functional one is in fact not an overwhelming task. As an exercise I wrote one consisting of about 140 lines of F# code. It provides constraint primitives for all_different, validate_pair and has_value which makes it possible to solve both Sudoku and N-queen problems.

Ignoring a few utility functions for row, column and block extraction the Sudoku solver boils down to:

type Puzzle = int list

let sudoku (puzzle : Puzzle) : Puzzle list =

let vars = newVars 81 (Set.ofList [1..9])

let givenConds =

[ for n, var in List.zip puzzle vars do

if n > 0 then

yield HasValue(var,n) :> FDCond ]

let gameConds = List.map (fun vars -> AllDifferent(vars) :> FDCond) (rows vars @ columns vars @ boxes vars)

labelling vars (givenConds @ gameConds)

The N-queen problem can be solved using:

let nqueen n =

let allvars = newVars n (Set.ofList [1..n])

let rows = AllDifferent(allvars) :> FDCond

let rec diagonals vars : FDCond list =

let rec diagonals’ v vrest i =

match vrest with

| vr::vrs ->

let j = i + 1

(ValidatePair(v, vr, fun v1 v2 -> v1 <> v2 + j) :> FDCond) ::

(ValidatePair(v, vr, fun v1 v2 -> v1 <> v2 – j) :> FDCond) ::

diagonals’ v vrs j

| [] -> []

match vars with

| v::vs -> diagonals’ v vs 0 @ diagonals vs

| [] -> []

labelling allvars (rows :: diagonals allvars)

The core of the solver consists of the following functions:

let rec solve (state : FDState) (conds : FDCond list) : FDState list =

let state’ = propagateConstraints state conds

if state’.Unique then

if evalConstraints state’ conds then

[ state’ ]

else

[]

else if state’.Invalid then

[]

else

[

let leastAmbiguousVar = state’.LeastAmbigousVar

for value in leastAmbiguousVar.domain do

let state” : FDState = state’.UpdateVar { leastAmbiguousVar with domain = Set.singleton value }

yield solve state” conds

] |> List.concat

let labelling (vars : FDVar list) (conds : FDCond list) : int list list =

[

let initialState = new FDState(vars)

for solution in solve initialState conds do

yield solution.Export

]

which means, algorithmically, that we need one mechanism for constraint propagation and one for (backtrack) search. The efficiency of the solver is related to how good it propagates constraints and this an area of active research that is beyond the scope of this post. Christian Bessier provides a lot of information on this subject in his paper, simply called constraint propagation.

## Prolog

One way to gain better skills in programming is to learn a new language. The pragmatic programmer recommends one new language every year and I guess that might be a good goal. Readers of Seven languages in seven weeks might learn even more, though it should be kept in mind that quantity and quality is not the same thing, I guess. The author expresses the idea of learning a new language (programming or spoken) using the following inspirational words:

A second language can help you encounter new worlds. You may even seek enlightenment, knowing every new language can shape the way you think.

The author covers Ruby, Io, Prolog, Scala, Erlang, Clojure and Haskell. One that has also been on my todo list is Prolog. Prolog belongs to the domain of logic programming and thus differs quite a bit from what I am used to, which of course contributes to my curiousness. I downloaded SWI-Prolog which was said to be mature and maybe more importantly, free. I decided to try to write a Sudoku solver, guessing that Prolog might be able to solve this with elegance, but also since it was also a problem I was being familiar with, having written solvers in both F#, C# and Powershell before. My task was prematurely aborted though when I searched for help about how to implement a matrix transpose method. The first hit when googling for “prolog tranpose” lead to a library called clpfd (Constraint Logic Programming over Finite Domains) which defined the function I was trying to express. However, as a sample usage of this function, the documentation also provided a – that’s right – a sudoku solver.

:- use_module(library(clpfd)). sudoku(Rows) :- length(Rows, 9), maplist(length_(9), Rows), append(Rows, Vs), Vs ins 1..9, maplist(all_distinct, Rows), transpose(Rows, Columns), maplist(all_distinct, Columns), Rows = [A,B,C,D,E,F,G,H,I], blocks(A, B, C), blocks(D, E, F), blocks(G, H, I). length_(L, Ls) :- length(Ls, L). blocks([], [], []). blocks([A,B,C|Bs1], [D,E,F|Bs2], [G,H,I|Bs3]) :- all_distinct([A,B,C,D,E,F,G,H,I]), blocks(Bs1, Bs2, Bs3).

Magnificent! A kind soul at stackoverflow had also published an explanation for prolog dummies (like myself). One important row was not explained though which had been pointed out but not corrected. The row not explained was

“append(Rows, Vs), Vs ins 1..9”

With a bit of trouble I was able to figure it out myself;

- append(Rows, Vs) is predefined and takes a list of lists and returns the concatenation. For example append([[1,2],[3],[4,5]],X) yields X=[1,2,3,4,5].
- Vs ins 1..9 means that every element in Vs should be constrained to values in the interval 1 to 9. Ins is defined I clpfd.

Nice, but I still wanted to write something myself. So… here is my n-queen solver instead:

:- use_module(library(clpfd)).

% usage: solution([X1,X2,X3,X4,X5,X6,X7,X8]).

solution(X) :-

length(X, Size),

X ins 1..Size,

all_distinct(X),

diagonal(X),

label(X).diagonal([]).

diagonal([X|XS]) :- diagonal(X, XS, 0), diagonal(XS).diagonal(_, [], _).

diagonal(X, [Y|YS], I) :-

J is I+1,

X #\= Y + J,

X #\= Y – J,

diagonal(X, YS, J).

I close this session by giving a reference to a prolog tutorial that I liked:

http://www.csupomona.edu/~jrfisher/www/prolog_tutorial/contents.html

(You might enjoy reading section 2.11 provides another n-queen solver and compare it to mine. The example uses a clever encoding to avoids checks for queens on the same row.)

## C++ AMP

Interesting video about the next(?) step in GPU programming from Microsoft.

Daniel Moth: Blazing-fast code using GPUs and more, with C++ AMP