Imperative Haskell

Posted on 29 May 2017

Get Colorings has kindly translated this post to Russian.

This post covers essentially the same material as a 5-minute presentation I gave at RC, because giving that talk over and over again doesn’t scale and there are things I would like to cover that are difficult within that time limit.

I was working through Tim Roughgarden’s Algorithms 1 (which has now been replaced by two smaller courses) and attempting to do all the exercises in Haskell when I bumped up against an uncomfortable truth. Haskell’s ‘quicksort’:

= [] qsort [][] : xs) = lt ++ [x] ++ gt qsort (xxs)lt[x]gt where lt = qsort [e | e <- xs, e < x] ltqsort [exs, ex] = qsort [e | e <- xs, e >= x] gtqsort [exs, ex]

isn’t a true quicksort! Specifically, it doesn’t sort the elements in place, and the assignment I was working on involved counting the number of comparisons, so I couldn’t get away with my fake quicksort. With my tail between my legs, I gave up on my pure Haskell approach and implemented a solution in Python:

import sys sys 10000 ) sys.setrecursionlimit( def partition_first(array, l, r): partition_first(array, l, r): p = array[l] array[l] i = l + 1 for j in range (l + 1 , r): (l, r): if array[j] < p: array[j]p: = array[i], array[j] array[j], array[i]array[i], array[j] i += 1 -1 ] = array[i -1 ], array[l] array[l], array[iarray[i], array[l] return (i -1 ) (i def partition_last(array, l, r): partition_last(array, l, r): -1 ], array[l] = array[l], array[r -1 ] array[r], array[l]array[l], array[r return partition_first(array, l, r) partition_first(array, l, r) def partition_median(array, l, r): partition_median(array, l, r): = choose_median(array, l, r) p_idxchoose_median(array, l, r) = array[l], array[p_idx] array[p_idx], array[l]array[l], array[p_idx] return partition_first(array, l, r) partition_first(array, l, r) def choose_median(array, l, r): choose_median(array, l, r): = array[l] headarray[l] = array[r -1 ] lastarray[r = r - l length if length % 2 == 0 : length = l + (length // 2 ) - 1 mid_idx(length else : = l + (length // 2 ) mid_idx(length = array[mid_idx] midarray[mid_idx] = [(l, head), (mid_idx, mid), (r -1 , last)] options[(l, head), (mid_idx, mid), (r, last)] max (options, key = lambda v: v[ 1 ])) options.remove((options, keyv: v[])) min (options, key = lambda v: v[ 1 ])) options.remove((options, keyv: v[])) return options[ 0 ][ 0 ] options[][ def quicksort(array, start, end, partition): quicksort(array, start, end, partition): global comparisons comparisons if end <= start: return endstart: else : = partition(array, start, end) p_idxpartition(array, start, end) += (end - start -1 ) comparisons(endstart quicksort(array, start, p_idx, partition) + 1 , end, partition) quicksort(array, p_idx, end, partition) = 0 comparisons = contents.copy() inp1contents.copy() 0 , len (inp1), partition_first) quicksort(inp1,(inp1), partition_first) print (comparisons) (comparisons) = 0 comparisons = contents.copy() inp2contents.copy() 0 , len (inp2), partition_last) quicksort(inp2,(inp2), partition_last) print (comparisons) (comparisons) = 0 comparisons = contents.copy() inp3contents.copy() 0 , len (inp3), partition_median) quicksort(inp3,(inp3), partition_median) print (comparisons) (comparisons)

This implementation is not particularly Pythonic: note the recursion limit and the use of a global variable. I actually forgot to reset the variable to 0 between iterations, which was fun to track down. But it works!

So far, so good. This isn’t something we’d be able to do in Haskell, right? And even if we could, the equivalent implementation would be so different as to be unrecognisable. At least this is what I thought until I took a closer look at Control.Monad.ST and Data.STRef.

One of my biggest gripes with Haskell is the quality of the documentation. Control.Monad.ST is introduced as

This library provides support for strict state threads, as described in the PLDI ’94 paper by John Launchbury and Simon Peyton Jones Lazy Functional State Threads.

and Data.STRef is introduced as

Mutable references in the (strict) ST monad.

I don’t want to read a paper to figure out how to use these libraries, and in fact I don’t have to! In recognition of this, I humbly present alternative descriptions for Control.Monad.ST :

You asked for mutable state, here it is!

and Data.STRef :

Variables that you can actually vary!!!1!1!one!1eleventyone

Code utilising these libraries can look very familiar to people used to imperative languages, e.g. past me. Here’s the above quicksort rewritten in Haskell:

{-# LANGUAGE RankNTypes #-} import Control.Monad.ST import Data.STRef import Data.Vector (fromList, toList, freeze, thaw) (fromList, toList, freeze, thaw) import Control.Monad import Data.Vector.Mutable ( STVector , read, write, swap) , read, write, swap) import qualified Data.Vector as V ( Vector , length) , length) import Data.List (sortOn) (sortOn) import Prelude hiding (read) (read) = fromList contents vectorfromList contents = do partitionFirst array l r p <- read array l array l i <- newSTRef (l + 1 ) newSTRef (l + 1 .. (r - 1 )] $ \j -> do forM_ [l(r)]\j <- read array j arrayJarray j <- readSTRef i i'readSTRef i < p) $ do when (arrayJp) swap array i' j + 1 ) modifySTRef' i ( <- readSTRef i i'readSTRef i - 1 ) l swap array (i') l return (i' - 1 ) (i' = do partitionLast array l r - 1 ) l swap array (r) l partitionFirst array l r = do partitionMedian array l r p <- chooseMedian array l r chooseMedian array l r swap array p l partitionFirst array l r = do chooseMedian array l r h <- read array l array l t <- read array (r - 1 ) array (r let len = r - l len let mid = if (len `mod` 2 ) == 0 mid(len then l + (len `div` 2 ) - 1 (len else l + (len `div` 2 ) (len m <- read array mid array mid let options = sortOn snd [(l, h), (mid, m), (r - 1 , t)] optionssortOn[(l, h), (mid, m), (r, t)] return ( fst (options !! 1 )) (options)) = when (start < end) $ do quicksort array start end partition comparisonswhen (startend) i <- partition array start end partition array start end + (end - start - 1 )) modifySTRef' comparisons ((endstart)) quicksort array start i partition comparisons + 1 ) end partition comparisons quicksort array (i) end partition comparisons quicksort' :: Ord a => V.Vector a -> ( forall s a . ( Ord a) => STVector s a -> Int -> Int -> ST s Int ) -> Int s aa)s a = runST $ do quicksort' vector partitionrunST <- thaw vector arraythaw vector <- newSTRef 0 compsnewSTRef 0 (V.length vector) partition comps quicksort array(V.length vector) partition comps readSTRef comps quicksort' vector partitionFirst quicksort' vector partitionLast quicksort' vector partitionMedian

This is roughly the same length as the Python implementation, and even improves on it in some ways: no recursion limit fiddling and no global variables.

If we can write Haskell that resembles Python, and Python is executable pseudocode, can we cut out the middleman and translate pseudocode directly to Haskell? Let’s take a look at another problem.

I needed to calculate the size of the strongly connected components of a graph for another assignment, and I decided to use Tarjan’s Strongly Connected Components algorithm. The pseudocode for that (as taken from Wikipedia) is:

algorithm tarjan is input: graph G = (V, E) output: set of strongly connected components (sets of vertices) index := 0 S := empty array for each v in V do if (v.index is undefined) then strongconnect(v) end if end for function strongconnect(v) // Set the depth index for v to the smallest unused index v.index := index v.lowlink := index index := index + 1 S.push(v) v.onStack := true // Consider successors of v for each (v, w) in E do if (w.index is undefined) then // Successor w has not yet been visited; recurse on it strongconnect(w) v.lowlink := min(v.lowlink, w.lowlink) else if (w.onStack) then // Successor w is in stack S and hence in the current SCC // Note: The next line may look odd - but is correct. // It says w.index not w.lowlink; that is deliberate and from the original paper v.lowlink := min(v.lowlink, w.index) end if end for // If v is a root node, pop the stack and generate an SCC if (v.lowlink = v.index) then start a new strongly connected component repeat w := S.pop() w.onStack := false add w to current strongly connected component while (w != v) output the current strongly connected component end if end function

and here’s what that looks like in Haskell:

import qualified Data.Array as A import qualified Data.Graph as G import Control.Monad (forM_, when) (forM_, when) import Control.Monad.ST import Data.STRef import Data.Vector.Mutable ( STVector , read, replicate, write) , read, replicate, write) import Prelude hiding (read, replicate) (read, replicate) = runST $ do tarjan graphrunST index <- newSTRef 0 newSTRef <- newSTRef [] stacknewSTRef [] <- replicate size False stackSetsize <- replicate size Nothing indicessize <- replicate size Nothing lowlinkssize <- newSTRef [] outputnewSTRef [] $ \v -> do forM_ (G.vertices graph)\v <- read indices v vIndexindices v == Nothing ) $ when (vIndex index stack stackSet indices lowlinks output strongConnect v graphstack stackSet indices lowlinks output reverse <$> readSTRef output readSTRef output where size = snd (A.bounds graph) + 1 size(A.bounds graph) index stack stackSet indices lowlinks output = do strongConnect v graphstack stackSet indices lowlinks output i <- readSTRef index readSTRef Just i) write indices v (i) Just i) write lowlinks v (i) index ( + 1 ) modifySTRef' push v A.! v) $ \w -> read indices w >>= \found -> case found of forM_ (graphv)\windices w\foundfound Nothing -> do index stack stackSet indices lowlinks output strongConnect w graphstack stackSet indices lowlinks output =<< ( min <$> read lowlinks v <*> read lowlinks w) write lowlinks vlowlinks vlowlinks w) Just {} -> read stackSet w >>= \wOnStack -> when wOnStack $ {}stackSet w\wOnStackwhen wOnStack =<< ( min <$> read lowlinks v <*> read indices w) write lowlinks vlowlinks vindices w) <- read lowlinks v vLowLinklowlinks v <- read indices v vIndexindices v == vIndex) $ modifySTRef' output . ( : ) =<< addSCC v [] when (vLowLinkvIndex)modifySTRef' outputaddSCC v [] where = do addSCC v scc w <- pop pop let scc' = w : scc scc'scc if w == v then return scc' else addSCC v scc' scc'addSCC v scc' = do push e : ) modifySTRef' stack (e True write stackSet e = do pop e <- head <$> readSTRef stack readSTRef stack tail modifySTRef' stack False write stackSet e return e

Aside from explicitly declaring our variables and passing them around, I think this looks pretty close.

How do we square this with Haskell’s reputation for purity and referential transparency? That’s the subject of the paper mentioned above that you don’t have to read (but totally can if you want)! They figured out a way to provide a principled pure interface to mutable state by passing the references as arguments into each function that makes use of them and leveraging the type system to make sure any impurity is well contained. The correctness of this approach was very recently verified. If desired, we can replace any of the functions with purer and more idiomatic definitions without changing the output, and that satisfies the definition of referential transparency!

Why don’t we do this all the time, when Haskell is at least a serviceable imperative language? Because writing imperative programs is hard! They don’t compose as well, have less useful type signatures, and are harder to reason about. Getting away from those things is why we have Haskell to begin with! The real question should be: how can we avoid doing things this way as much as possible?

Before I discovered this part of Haskell, I had this perception of Haskell (and declarative programming more generally) as “imperative programming but less” from a practical perspective. I thought that although writing declarative code in Python was purely (heh) a matter of discipline, writing imperative code in Haskell required completely reconceptualising the algorithm. Thanks to ST , I now know that this not the case, which is a huge relief. If required, I can do a literal translation of the algorithm, and clean it up (or not) later. In fact Haskell is “imperative programming and more”, and that’s awesome!

Thanks to Peter Fraenkel, Julia Evans, and Michelle Steigerwalt for feedback.

If you’d rather try to make sense of the set of disconnected files that constitutes my slides for that presentation, you can do that instead, although I wouldn’t recommend it.