## 2010 April 4

### On powersets and folds

Filed under: algorithms, corecursion, haskell — lpsmith @ 10:03 pm

A while ago, somebody on reddit brought up a particularly cute and fairly well known definition of powerset:

```powerset :: [a] -> [[a]]
powerset = filterM (const [False,True])
```

Unfortunately, it is also fairly well known that this is a very space-inefficient way of producing the powerset, for example evaluating `length (powerset [1..32])` will allocate gigabytes of memory within seconds on modern computers, and is an effective denial-of-service attack against most default configurations.

Fortunately, there are solutions that avoid this space leak while maintaining the same style of definition. One solution is to use a revised definition of Control.Monad.filterM. Here is the standard definition:

```filterM          :: (Monad m) => (a -> m Bool) -> [a] -> m [a]
filterM _ []     =  return []
filterM p (x:xs) =  do
flg <- p x
ys  <- filterM p xs
return (if flg then x:ys else ys)
```

Here is the revised definition:

```filterM'          :: (Monad m) => (a -> m Bool) -> [a] -> m [a]
filterM' _ []     =  return []
filterM' p (x:xs) =  do
ys  <- filterM' p xs
flg <- p x
return (if flg then x:ys else ys)
```

The only difference is that we changed the order of two lines. This in turn changes the order in which the subsets are produced:

```powerset  "abc" == ["","c","b","bc","a","ac","ab","abc"]
powerset' "abc" == ["","a","b","ab","c","ac","bc","abc"]
```

These orderings are analogous to counting in binary, with the original definition having the least significant bit at the end of the list, whereas the revised versions have the least signficant bit at the beginning:

```pow      pow'

000 *    000 *
001      100 *
010      010
011      110
100 *    001
101      101
110      011
111      111
```

Both solutions have maximal sharing among the tails of the sublists, but the second solution arranges the sharing so that shared sublists appear adjacent to each other in the overall order, as illustrated by the asterisks. This allows `length (powerset' [1..32])` to be evaluated in a very modest amount of space.

Note that this space leak is inherent to any solution that produces the first ordering, and has maximal sharing! If the first ordering is needed, it is better to recompute the sublists than to share them. As Cale Gibbard pointed out, one way to eliminate the sharing is to use another monad, such as the Logic monad.

```import Control.Applicative ((<|>))

powerset2 :: [a] -> [[a]]
powerset2 = observeAll . filterM (const (return False <|> return True))
```

## Making powerset lazier

The definition of `powerset` can be generalized to produce all the finite sublists of a potentially infinite list. However, it is not obvious how to express this via tweaking `filterM` without breaking the abstraction, so let’s start with another efficient formulation based on concatMap.

```powerset' [] = [[]]
powerset' (x:xs) = concatMap (\ys -> [ys, x:ys]) (powerset' xs)
```

On an infinite list, this definition gets stuck in an non-productive loop because `powerset'` must traverse the entire list before it returns anything. Note that on finite cases, the first sublist returned is always `[]`. Formally, we can state this as `powerset' xs == [] : tail (powerset' xs)`. It is sensible to extend this equality to the infinite case, due to the analogy to counting in binary. By making this substitution, we produce a lazier definition that can handle infinite lists, from which we can calculate a definition that’s more aesthetically pleasing and slightly more efficient:

```powerset' []     = [[]]
powerset' (x:xs) = concatMap (\ys -> [ys, x:ys]) ([] : tail (powerset' xs))
-- apply definition of concatMap
= [] : [x] : concatMap (\ys -> [ys, x:ys]) (tail (powerset' xs))
```

We can clean this definition up by calculating definitions for `tail . powerset'`. We start by applying `tail` to both sides of the two cases:

```tail (powerset' [])     = tail [[]]
= []
tail (powerset' (x:xs)) = tail ([] : [x] : concatMap (\ys -> [ys, x:ys]) (tail (powerset' xs)))
= [x] : concatMap (\ys -> [ys, x:ys]) (tail (powerset' xs))
```

Next, we can choose a new name for `tail . powerset'`, such as `nonempties`, and substitute it through the definition. Finally, we can recover something equivalent to the original definition by appending the empty list to `nonempties`

```lazyPowerset xs = [] : nonempties xs
where nonempties [] = []
nonempties (x:xs) = [x] : concatMap (\ys -> [ys, x:ys]) (nonempties xs)
```

This is equivalent to Data.List.subsequences as added in GHC 6.10.1. I thank Bertram Felgenhauer for (unknowingly) providing some comments that were rather helpful in creating this presentation.

## Generalizing powerset to an enriched fold

By parameterizing the instances of `(:)` and `[]` that create the inner sublists, we can turn `lazyPowerset` into a kind of fold:

```powerfoldr :: (a -> b -> b) -> b -> [a] -> [b]
powerfoldr f b as = b : nonempties as
where nonempties [] = []
nonempties (a:as) = f a b : concatMap (\b' -> [b', f a b']) (nonempties as)
```

For example, we can calculate the sum of every finite combination of the powers of 2:

```powerfoldr (+) 0 (iterate (2*) 1)  == [0,1,2,3,4..]
```

An equivalent “higher-level” definition of `powerfoldr` in terms of `map`, `foldr`, and `lazyPowerset` is as follows:

```powerfoldr' f b = map (foldr f b) . lazyPowerset
```

This definition, however, recomputes the right fold over common sublists; the sharing provided by `lazyPowerset` is not preserved by this latter definition. The first definition has the same level of sharing as the original, reducing the number of applications of `f` by half. I found the first definition of `powerfoldr` to be quite helpful in solving Project Euler Problem 268; it certainly was not optimal for the purpose, but it was more than adequate.

## 2009 December 30

### Fun with the Lazy State Monad

Filed under: continuations, corecursion, haskell, monads — lpsmith @ 9:20 pm

The lazy state monad doesn’t work the way you think it works. Thinking of the lazy state monad in terms of imperative programming is a very useful first approximation, at least in terms of what is computed, but in terms how things are computed this intuition is beyond useless. I hope to dissuade you of the latter part of the intuition in this post, by demonstrating two of the more interesting things that can be done with the lazy state monad.

Albert Y.C. Lai recently shared a program that demonstrates the lazy state monad in a rather interesting way:

```pro :: State [Bool] ()
pro = do
pro
s <- get
put (True : s)
```

Instead of tail recursion, this employs head recursion. The question is, what does `pro` return when you run it, and why? I found it easy to guess the correct answer, but my reasoning was completely wrong. Of course, if viewed through a wholly imperative mindset, this leads to non-termination, but the lazy state monad extracts usable information from this definition.

In my recent Monad Reader article, I disentangled the bits of circular programming from Lloyd Allison’s queue, and for the last year I have been obsessed with pulling the same trick with Jones and Gibbons’ breadth-first labeling algorithm. At times, this obsession has been a bit of a disease. It’s also worth pointing out Chris Okasaki’s discussion of a particular instance of this algorithm, and the mention of this algorithm on the FP Lunch Weblog. Here is code for the special case of breadth-first numbering:

```data Tree a b = Leaf a | Branch b (Tree a b) (Tree a b)

lazyRenum :: Num n => Tree a b -> Tree n n
lazyRenum t = t'
where
(ns, t') = renumber (0:ns, t)

renumber (n:ns,  Leaf    _    ) = (n+1:ns  , Leaf    n      )
renumber (n:ns,  Branch  _ l r) = (n+1:ns'', Branch  n l' r')
where
(ns' , l')  = renumber (ns , l)
(ns'', r')  = renumber (ns', r)
```

I finally disentangled the corecursive bits from this example today. The circular programming occurs in the list argument, not the tree. Note that the flow of the list `ns` from one call of `renumber` to the next is like the state monad. From this observation, I wrote the following whimsical code:

```lazyRenum :: Num n => Tree a b -> Tree n n
lazyRenum = runFresh . renumber
where
renumber (Leaf    _    )
= fresh (\n -> do
return (Leaf n))
renumber (Branch  _ l r)
= fresh (\n -> do
l' <- renumber l
r' <- renumber r
return (Branch n l' r'))
```

Once I had this right, it was pretty easy to fill in the definitions for `fresh` and `runFresh`, by cribbing off Chris Okasaki’s simplification of Jones and Gibbons’ algorithm:

```type Fresh n a = State [n] a

runFresh :: Num n => Fresh n a -> a
runFresh m = a
where
(a, ns) = runState m (0:ns)

fresh :: Num n => (n -> Fresh n a) -> Fresh n a
fresh f = do
(n:ns) <- get
put ns
a <- f n
ns' <- get
put ((n+1) : ns')
return a
```

And finally, we have arrived at a way to disentangle Jones and Gibbons’ algorithm. This easily generalizes from breadth-first numbering to breadth-first labeling, and like the original, it is capable of renumbering the Stern-Brocot tree. The key insights here are the use of the lazy state monad, and getting the type of `fresh` correct. Everything else is relatively straightforward, once this groundwork is laid.

It’s interesting to perform a post-mortem analysis of why coming up with this proved to be so difficult for me. I’ll admit that I spent a few weeks trying to decipher the operational characteristics of Jones and Gibbons’ labeling algorithm, and while I think I have a reasonable grasp on it, I’m still not completely comfortable with it. However, this monadic abstraction seems perfectly obvious in hindsight.

This contrasts starkly to my own understanding of Lloyd Allison’s queue: I was completely comfortable with the operational aspects of the algorithm, but the monadic abstraction was quite difficult to come to grips with. So my difficulties with Jones and Gibbons’ algorithm was an over-emphasis on the operational aspects of the algorithm, and too much focus on the continuation monad as part of the solution. Basically, I was hopeful that the same methodologies in my Monad Reader article would be directly useful, so much so that I had a case of tunnel vision.

While it is not obvious how the continuation monad might be applicable to this problem, continuations still did play a role in the solution: look at the type of `fresh`, it’s the same `(a -> r) -> r` pattern that the continuation monad uses. It seems to me that there is something deeper here, although I don’t know what that would be.

These examples might make a stronger case for the lazy state monad than the example in my previous post. While `foldrByLevel` is relatively easy to adapt to the continuation passing state monad, I don’t know how to do the same with either of these two examples.

## 2009 June 22

Filed under: continuations, corecursion, haskell, monads, queues — lpsmith @ 8:41 pm

Haskell aficionados, take note! My library for corecursive queues has now been uploaded to Hackage. You can now cabal-install it.

I also have a substantially revised draft of the associated paper, Lloyd Allison’s Corecursive Queues, available. It has been re-organized so that it is hopefully easier to follow, it includes a careful performance comparison, and a tentative proof that `mapCont` cannot be expressed in terms of `callCC`, `(>>=)`, `return`.

The library includes a somewhat crude micro-benchmarking program in the `tests/` directory. Those who have read previous drafts, be warned that the few brief statements about performance were based on past notes, and I found some several issues with the testing methodology contained in the notes. Here the revised results:

Description Time (ms) -H500M Bytes allocated
GHC 6.10.3 mean σ mean σ per Branch
levelOrder’ 446 5 172 15 44.0
CorecQ 555 5 619 4 133.5
CorecQW _ 696 5 1128 6 213.6
CorecQW () 907 56 2235 11 213.6
Side Channel _ 959 3 1171 7 228.7
Side Channel () 1500 56 2171 7 276.4
STQ 1140 8 1087 14 371.2
TwoStack 1158 4 778 10 185.8
Okasaki 1553 7 1574 12 209.0
Data.Sequence 962 5 1308 5 348.1
GHC 6.8.3
levelOrder’ 461 2 173 15 44.1
CorecQ 458 4 267 13 67.5
CorecQW _ 526 5 713 5 141.2
CorecQW () 781 62 1775 62 141.3

These benchmarks come from performing breadth-first traversals repeatedly on the 34th fibonacci tree, on an Intel Core 2 Duo T9550. The first few data points were discarded, and the mean and standard deviation of the remaining times were computed. Note that `getCPUTime` was used to time each run, and this has a resolution of only 10 milliseconds.

If you would like to play with the queue transformer, which doesn’t appear in the library, or other bits of code exactly as they appear in the paper, you can download the source code here.

## 2009 April 7

### Polynomial multiplication

Filed under: algorithms, corecursion, haskell — lpsmith @ 6:27 am

Alright, it’s been almost a month since I’ve last posted, and I wanted to post weekly. Perhaps bi-weekly would be a better goal, as long as I really can stick to it.

So, a few months ago, while working on a combinatorial problem for Project Euler, I was re-implementing polynomial arithmetic, and made an interesting observation. I was using a sparse list representation, with each element of the list being `(M c k)`, representing the monomial `cxk`.

```> data Mono = M !Integer !Integer  deriving (Eq, Ord, Show)
```

A polynomial is represented by a list of monomials, which is assumed to be sorted in ascending order of the power `k` above, and every power is assumed to be unique in our representation. Thus we can easily define addition and subtraction using standard techniques for dealing with ordered lists:

```> mergeBy :: (Integer -> Integer -> Integer) -> [Mono] -> [Mono] -> [Mono]
> mergeBy f xs ys = loop xs ys
>   where
>     loop [] ys = ys
>     loop xs [] = xs
>     loop (M c0 k0:xs) (M c1 k1:ys)
>            = case compare k0 k1 of
>                LT -> M c0 k0 : loop xs (M c1 k1:ys)
>                EQ -> let c' = f c0 c1
>                       in if c' == 0
>                           then loop xs ys
>                           else M c' k0 : loop xs ys
>                GT -> M c1 k1 : loop (M c0 k0:xs) ys

> add = mergeBy (+)
> sub = mergeBy (-)
```

It’s also very straightforward to multiply two monomials, and to perform a “scalar” multiplication of a single monomial with a polynomial. Note that the degree of each term is changed by a constant amount, so we don’t have to worry about re-ordering the list, and that the first case of `smul` is a completely optional minor optimization.

```> mmul (M c0 k0) (M c1 k1) = M (c0 * c1) (k0 + k1)

> smul (M 1 0) ys = ys
> smul x ys = map (mmul x) ys
```

Taking into account the possibility of polynomials with an infinite number of terms, it’s now fairly easy to implement a polynomial multiplication program in terms of `add` and `mmul`:

```> pmul _  []    = []
> pmul [] (_:_) = []
> pmul (x:xs) (y:ys) = mmul x y : add (smul x ys) (pmul xs (y:ys))
```

This is based around the algebraic identity `(x + xs) * (y + ys) == x * y + x * ys + xs * (y + ys)`. It works on infinite polynomials because `x` and `y` are the terms of least degree to be multiplied, and every other product will have a higher degree. Unfortunately, this code turns out to be hopelessly inefficient.

The second thing I tried was simply changing the order in which arguments are passed to the recursive call to `pmul`:

```> pmul2 _  []    = []
> pmul2 [] (_:_) = []
> pmul2 (x:xs) (y:ys) = mmul x y : add (smul x ys) (pmul2 (y:ys) xs)
```

Empirically, this is asymptotically faster than the first attempt. I do not know why. Any insight would be fun to read about!

## 2009 March 9

### Lloyd Allison’s Corecursive Queues

Filed under: continuations, corecursion, haskell, monads, queues — lpsmith @ 9:34 pm

I’m proud to announce that a draft of the release of my paper, “Lloyd Allison’s Corecursive Queues: Why Continuations Matter“, is now available. (Source code available here, with a hackage package soon to come available here) Wouter Swierstra has informed me that he will publish it in the Monad Reader. However, it will appear in the next issue after the upcoming one, due to an unexpectedly large number of submissions this time around. Here is the abstract:

In a purely functional setting, real-time queues are traditionally thought to be much harder to implement than either real-time stacks or amortized O(1) queues. In “Circular Programs and Self-Referential Structures,” Lloyd Allison uses corecursion to implement a queue by defining a lazy list in terms of itself. This provides a simple, efficient, and attractive implementation of real-time queues.

While Allison’s queues are general, in the sense it is straightforward to adapt his technique to a new algorithm, a significant problem has been the lack of a reusable library implementation. This paper solves this problem through the use of a monadic interface and continuations.

Because Allison’s queues are not fully persistent, they cannot be first class values. Rather, they are encoded in particular algorithms written in an extended continuation passing style. In direct style, this extension corresponds to `mapCont`, a control operator found in `Control.Monad.Cont`, part of the Monad Template Library for Haskell. This paper conjectures that `mapCont` cannot be expressed in terms of `callCC`, `return`, and `(>>=)`.

I intend to include a careful performance comparison before this becomes an official Monad Reader article. Allison’s queues come out very well; often better than two stack queues. I have conducted a careful performance comparison in the past, although with older versions of GHC, and older versions of my code. While I did take reasonably careful notes, things have changed. Haskell being what it is, figuring out why is often a challenge. In the meantime I am interested in feedback.

For fun, here is something I wrote right after I first came up with the basic idea behind the paper. It’s still the best error message I’ve gotten out of GHC. Kudos to whomever came up with that strategically placed bit of humor!

Thursday, August 25th, 2005, 5:22 am: Back to the Future

I’ve been up all night, but I now have a working fragment of computer code that is entirely too cute. It’s easily the cleverest bit I’ve written in years. I managed to implement… a queue.

Yes, a queue. One queue, not two. One purely functional queue, with one esoteric implementation! On my penultimate attempt, which was an utter failure except that it got me thinking in the right direction, I came by the most amusing error message I’ve seen to date out of GHC:

```leon@deleon:~/Programs/snippets \$  ghci -fglasgow-exts Queue.hs
___         ___ _
/ _ \ /\  /\/ __(_)
/ /_\// /_/ / /  | |      GHC Interactive, version 6.2.2, for Haskell 98.
/ /_\\/ __  / /___| |      http://www.haskell.org/ghc/
\____/\/ /_/\____/|_|      Type :? for help.

Compiling Queue            ( Queue.hs, interpreted )

Queue.hs:84:
My brain just exploded.
I can't handle pattern bindings for existentially-quantified constructors.

...

Failed, modules loaded: none.
Prelude>```

Yeah, mine did too. Even I don’t fully understand my own code yet.

It should be noted that I knew full well that the code I was trying wouldn’t work… but after hours of bewilderment, not even trying to load anything into GHCi, for amusement’s sake I simply had to try something.

Update: (March 23)
— Data.Sequence is not a real time queue: rather, they are amortized.
— Added citation to Chris Okasaki’s Purely Functional Data Structures
— Other minor changes