## Modelling Stochastically Independent Processes with F# Computation Expressions: Part 1

The idea for doing this is not new. There is an excellent series of posts closely tracing an article on applications of functional programming to probability.

A colleague of mine has recently called my attention to his own post of two years ago, where he describes a monad that models stochastically independent events in Clojure. I thought the latter implementation actually went to the heart of the whole idea of monads and that is why, once I started writing my own in F# from scratch, event though I naturally/brute-force-like, chose something described below and which corresponds exactly to the series of posts above, I eventually was compelled to do it my colleague’s way.

Here is the natural choice: a distribution (p.m.f) is just a sequence of events. Each event is simply a record/tuple/map of two values: the signifier of the event and its probability.

//underlying type type 'a Outcome = { Value: 'a Probability : BigRational } //Monadic type type 'a Distribution = 'a Outcome seq

This is great, and it works very well, as the posts show. My slight dissatisfaction with this approach is due to the fact that the competing one appears much more transparent, extremely easy to read, easy to use. In fact, anyone who knows Clojure syntax can use it right away to model independent processes simply by writing them out in Clojure, AND there is no need for any helper functions to extract the results. It just works (from the post, with some minor corrections):

(domonad probability-m [die1 (uniform [1 2 3 4 5 6]) die2 (uniform [1 2 3 4 5 6])] (+ die1 die2))

This code models an experiment of rolling two dice and returns a probability mass function that describes the experiment. All we had to do as users was describe the experiment in the most natural way possible: what is the probability of getting an exact value as the sum of both values of 2 rolled dice. Hence the expression: extract the value from the first pmf, add it to the value extracted from the second. Don’t even need to think about probabilities, as the magic happens behind the monadic scene.

The code above gives us the result: {2 1/36, 3 1/18, 4 1/12, 5 1/9, 6 5/36, 7 1/6, 8 5/36, 9 1/9, 10 1/12, 11 1/18, 12 1/36} – which is a map of the sum of dice values to the probability of getting them. In F#, I believe, the Algol-like syntax actually works to our advantage (I call the computation expression “independent” for obvious reasons):

independent { let! x = uniform [1..6] let! y = uniform [1..6] return x + y }

When executed in F# interactive, and using the FSharp.PowerPack:

val mp : Map<int,BigRational> = map [(2, 1/36N); (3, 1/18N); (4, 1/12N); (5, 1/9N); (6, 5/36N); (7, 1/6N); (8, 5/36N); (9, 1/9N); (10, 1/12N); ...]

We, amazingly enough, get the same answers. Next time: better examples and nuts & bolts of the code.

## Generating Permutations: Clojure or F#: Part 2

Marching on from the last post.

### Lazy Sequences

This is my favorite feature ever. If I want to generate just a few of 10! (nobody even knows how much that is) permutations, I could:

(take 10 (permute [1 2 3 4 5 6 7 8 9 10]))

provided, the function is defined (as described in the first post):

(defn permute [v] (when-let [[pos2 pos1] (findStartingPos v)] (let [nxt (sort-remainder (swapDigits v pos2 pos1) (inc pos1))] (cons nxt (lazy-seq (permute nxt))))))

Here I am not sure which language I like more. Clojure has easier syntax: everything fits nicely within the recursive function call. Returning nil terminates the loop, while in F# you need to know to return an option type where None terminates iteration. On the other hand, I like the fact that everything is neatly wrapped in the “unfold” function: seems more “natural” to me: fold/unfold – there is a certain symmetry here. Also, everything exists nicely in this LINQ-like world…

let permute (v : 'a array when 'a: comparison) = Seq.unfold (fun prev -> match findStartingPos prev with | None -> None | Some (cur, pos) -> Some(prev, sortRemainder (swapPositions prev cur pos) (pos + 1))) v

### Weak Typing

I really like Clojure weak typing. And I like the F# strong type system:

let sortRemainder (v : 'a array) pos = if v.Length - 1 = pos then v else [| yield! v.[0..pos - 1] yield! Array.sort v.[pos..v.Length - 1]; |]

F# type system requires that the first argument be qualified, but it is happy with this abbreviation, while the full qualification should be:

let sortRemainder (v : 'a array when 'a: comparison) pos =

Since we are sorting a subvector, the array has to be of a “comparable” type. Which is the condition of the applicability of the algorithm.

In Clojure it looks simpler, but it’s essentially the same:

(defn sort-remainder [v pos1] (if (= (dec (count v)) pos1) v (into (subvec v 0 pos1) (sort (subvec v pos1)))))

### Tail Recursion

One more cool feature of functional languages. I think it’s another tie once you use it, although the “loop” construct that demands it is very nice.

The following function returns a tuple (current, found) of two positions within the array: one of the element that is being “promoted” up (current), and the other – of the smaller element being pushed back. (So, current > found && v[current] < v[found]). Or nil/None if no such pair can be found. This is the key function of the algorithm:

(defn findStartingPos [v] (loop [cur (dec (count v)) acc [-1 -1]] (let [maxPos (second acc)] (if (or (< cur maxPos) (< cur 0)) (if (= maxPos -1) nil acc) (if-let [pos (findFirstLessThan v cur)] (recur (dec cur) (if (< maxPos pos) [cur pos] acc)) (recur (dec cur) acc))))))

And F#:

let findStartingPos v = let rec findStartingPosRec cur acc = let maxPos = snd acc if cur < 0 || cur < maxPos then if maxPos < 0 then None else Some acc else let pos = findFirstLessThan v cur match pos with | Some pos -> findStartingPosRec (cur - 1) (if maxPos < pos then (cur, pos) else acc) | None -> findStartingPosRec (cur - 1) acc findStartingPosRec (v.Length - 1) (-1, -1)

It’s nice that we have a “loop” keyword in Clojure to provide cleaner syntax and more discipline for defining tail-recursive functions, but I am not appalled with the way we do it in F# either.

(The above functions contain obvious optimizations: we stop scanning once we have a pair of “swappable” elements and we have moved beyond the “found” position. Also, we discard a valid pair if we already have a pair where “found” position is larger than the “found” position of the current iteration).

## Doing it in a Massively Parallel Way

Of course, I like everything parallel… So what about doing it on a GPU, using, say CUDA? It is definitely possible, although probably not very practical. Even if we only have an array of 10 distinct elements, the number of permutations is already ridiculously large (although who knows what we are going to be using them for)… In any event, this is solvable if we can get “random access” to permutations. Instead of unfolding them as a lazy sequence, generate them all at once in a massively parallel fashion.

This is possible because permutations are neatly tied to factoradic numbers, as this Wikipedia article explains. So, it is always possible to generate “permutation #10” to be guaranteed different from “permutation #5” for distinct, fully ordered sets. (Any sets where ordering relationship is not defined can still be easily permuted as long as its elements are stored in indexed data structures, such as arrays, by simply generating permutations of indices). Thus, taking CUDA “single data multiple threads” computation model it is easy to generate all (or however many) permutations in parallel. Naturally, if we are not just outputting the results but need to store them, the exponential nature of the problem memory growth, as well as the number of threads required, and the limited amount of GPU memory (a single computer RAM for that matter) will quickly become a problem. I guess the CUDA C++ version of this will have to wait until the next job interview…

## Generating Permutations: Clojure or F#: Part 1

## The Alogirthm

Recently, I have entered a brave (new?) world of Clojure and was looking for a small project to take it for a ride. I stopped on a popular/interesting enough little problem, that subsumed a certain interview question which I was once unfortunate enough to stumble through with half my brain tied behind my back (long story).

The algorithm to generate permutations in sequence is actually quite simple and intuitive, I believe it originates somewhere in India thousands of years ago. At least that was the impression I got from a Wikipedia article. Anyway, given a fully ordered set S, generate all permutations of it. The algorithm works for sets with repetition, but we can assume no repetitions for simplicity:

- (Since the set is fully ordered), sort it in ascending order, represent it as a sequence. Output it.
- Find the next sequence consisting of the elements of S and following the current one immediately based on the ordering relationship. Output it.
- Repeat 2, until no such sequence can be found.

The juice here is of course 2. To illustrate with an example. Suppose our sequence is just numbers “{}” here mean sequence, order matters: S ={1, 2, 3, 4}. The next one in lexicographical order is: {1, 2, 4, 3}. The one after it: {1, 3, 2, 4}. The one after: {1, 3, 4, 2}. Shall we do one more? {1, 4, 2, 3} Maybe just one more, so it becomes even clearer: {1, 4, 3, 2}.

So the key here, in case you have not caught on yet:

- Starting from the end of the sequence, going backwards, find the first smallest position in the sequence, if such exists, where the element in current position can be swapped with the element in the found position to make the sequence “larger” than the previous one. This is accomplished if the element at the current position is less than the element at the found position. Swap elements in the current and the found positions.
- Sort the sub-sequence starting at the found + 1 position, ascending. So, in the above example we get: {1, 3, 2, 4}

For instance, looking at {1, 2, 4, 3} we find: (2 < 3), (2 < 4), (1 < 2). But 2 < 3 is the very first pair that fits the bill (we are starting from the back), so we don't go any further and swap them. We get {1, 3, 4, 2}

Et c’est tout!

(How simple is that? I bet the brute force recursive “solution” feels very stupid right about now).

## Implementation

At this point it’s all but unnecessary to actually show the implementations. However, I just had fun with some of the Clojure concepts/structures, and wanted to compare them to what I am used to in F#. So this is the point. Meta-programming. The Clojure code is here and F# – here.

Observations:

### Immutability

The way the algorithm is structured, immutability is intuitive: we are producing a new sequence in each step, and old sequences should also be preserved and output. Should have been a no-brainer, since we are dealing with functional languages with immutable structures. And here was a slight jolt. Consider a function that swaps elements in the sequence (the algorithm calls for swapping, so…):

Clojure:

(defn swapDigits [v pos1 pos2] (assoc v pos2 (v pos1) pos1 (v pos2)))

(I was using strings of digits as input sets in my initial awkward steps through the language, hence the function name).

Now, I wanted something as short and as sweet as this for F#, so I cheated! .NET arrays are mutable and used I .NET arrays instead of F# lists (hope functional programming purists will commute a very harsh sentence I deserve for writing the following, and not condemn me to Java programming for life):

let swapPositions (v : 'a array) pos1 pos2 = let res = Array.copy v res.[pos1] <- v.[pos2] res.[pos2] <- v.[pos1] res

So yes. The code at line 2 is pure sin and I should have used F# lists which are immutable, and the code should have been:

let swapPositions (v : 'a list) pos1 pos2 = v |> List.permute (fun i -> if i = pos1 then pos2 elif i = pos2 then pos1 else i)

but frankly, this is not the (very) first thing that comes to mind if you haven’t been using the language for a long time, and not as clear as the Clojure solution, which was close enough to the surface for a noob like myself to scoop it up.

Alright. This has gone on long enough, more soon…