Monads from Two Perspectives

Just random scribblings.  This is pretty standard stuff, but I worked this out for someone recently, so I’m typing it out.

Monads For Mathematicians

If you talk to a category theorist and ask for a definition of a monad, it looks something like this (taken from Wikipedia): a monad is a triple consisting of an endofunctor , and natural transformations and , satisfying the following laws:

• 
• 

Enough mumbo-jumbo.  It’s worth taking a second to think about what that means.  The two natural transformations just give us a single way to treat any of the infinite family of objects (for any integer ) as just .  We can flatten them, in essence.  If there are too few ’s, you apply .  If there are too many, you apply .  The monad laws above just guarantee that it doesn’t matter at what level (number of applictaions of the functor) you work; you’ll get the same result.

It’s definitely time for an example.  Suppose you have a fixed monoid .  Now, in the category , there is a functor , which just takes any set to the cartesian product of itself and the monoid – i.e., consists of values from the underlying set, plus a “side value” from the monoid.  We can also define natural transformations , and where multiplication is taken from the monoid structure.  In other words, if you don’t give me a side value (that’s ), then I’ll take the side value to be the identity from the monoid.  If you give me several side values in a given order (that’s ) I’ll multiply them in the monoid to decide what the side value is.

Verifying the monad laws, stated above, is fairly straightforward.  Really, the tricky part is just to understand the meaning of the pieces.  Generally speaking, for a natural transformation , think of as applying the transformation on the left hand side of a tuple.  Then think of as being the specialization of to values of type .  Here’s what happens:

• If , then .  Then .
• If , then .  Then .  By associativity in the monoid, this is the same as the above, satisfying the first law.
• If , then and .  Verifying the second monad law is left as an exercise (just write out what the function composition means, and apply ).

Monads for Haskell Programmers

In Haskell, of course, we prefer to think of monads differently, in terms of a type constructor and two functions “return” and “(>>=)”, satisfying the following laws:

return a >>= f           = f a
m >>= return             = m
(m >>= f) >>= g          = m >>= (\x -> f x >>= g)

The “return” function takes any value, and wraps it in the monad.  The second operation, often read as “bind”, is a little trickier.  It takes one value wrapped in the monad, and function from a plain value of the first type to a second value wrapped in the monad, and produces the second value wrapped in the monad.  Here are a few examples:

1. Containers can be seen as monads.  In this case, if I have a crate of thingamabobs, and each of those thingamabobs can be unpacked into a crate of whatsits, then I may as well (by unpacking the thingamabobs as needed) have a very large crate full of whatsits.
2. Famously, I/O actions are monads.  If I have an I/O action that produces a character string, and for any character string, I could get an I/O action that produces a number… then I can string them together into a single interaction producing a number.  (The resulting I/O action will get the character string, compute the second action, and perform it to get the number, all in one go.)

A lot of this corresponds exactly to the mathematical definition above.  The type constructor corresponds with an endofunctor (in the category of Haskell types), and the “return” function with the natural transformation .  That is, “return” takes a value of any possible type T, and turns it into a value of type M T (M being our monad type constructor).  The idea of a natural transformation is partially captured by polymorphic functions in Haskell.  (Partially because we’d also need to verify that the tranformation commutes with morphisms in the category; i.e., that “return (f x)” is the same as “fmap f (return x)” for all values of f and x.)

It turns out that is actually the standard library operation called “join”, which is defined in terms of the above, as follows:

join a = a >>= id

In other words, if I have a value of type (M (M T)), then I can turn it into a value of type (M T), by grabbing the encapsulated value of type (M T) and, well, doing nothing at all with it.

There’s actually one more piece we’re missing.  In mathematical language, a functor simultaneously maps objects to other objects in a category, and also morphisms to other morphisms.  So to recover a functor from the Haskell definition of a monad, we also need to decide how to map functions from type X to type Y into other functions from type M X into type M Y.  It turns out we can recover this from the two operations, return and bind, that are standard for monads.

fmap f x = x >>= (return . f)

In other words, if I have a function , then I can easily obtain a function as follows.  Use the bind operation to grab the encapsulated value, apply the function f to it, and then wrap it back up again with the return operation.

In the other direction, given an appropriate “fmap” and “join”, the bind operation is easy to implement as well:

a >>= f = join (fmap f a)

Basically, the idea here is that because M is a functor, a function can be turned into a function .  Then one can just apply the join () operation on the result to obtain the desired result.

The Equivalence Between the Two

It is a little tedious, but not really difficult, to show that these two structures are actually completely equivalent.  In other words, using the definitions and correspondences above, if you have a Haskell monad satisfying the monad laws, then you also have a mathematical monad for the category of Haskell types.  Similarly, if you have a mathematical monad for the category of Haskell types, then you also have a Haskell monad.  In fact, as many people probably noticed, my first example of a monad at the beginning of this article is none other than Haskell’s Writer monad.

If you haven’t done it before, it’s really worth doing.

View Patterns as Pattern Matching for Records

I’ve learned a lot in the last day about record systems for Haskell built as libraries.  I came up with what seemed like an obvious answer, and indeed it was essentially the same as two existing packages on Hackage — fclabels and data-accessor — the second of which is among the more popular downloads out there.  That’s always a good sign, that a significant number of people needed the same thing, and liked the same answer.  But there are some things, still, that are provided by the built-in package system, either as standard features or GHC extensions, that it’s unclear how to accomplish as a library.

The one that’s intriguing me most right now is pattern matching, because I think it almost got inadvertently solved a couple years ago.

Right now, I can write this for a record type in Haskell.

foo :: MyRecord -> Int
foo (MyRecord { field1 = Just k, field2 = 4, field3 = (x:_) }) = 4*k - x

That’s a pattern match that disassembles the record, by named fields, and matches each field using the full variety of patterns available in any other context.  Pretty powerful stuff.  Unfortunately, it’s built in to Haskell’s rather limited record system, and is among the casualties incurred when switching to a different system.  One needs to fall back to the underlying “native” record.  In the case of the fclabels and data-accessor packages, for example, this means using a whole separate set of names involving underscores, further polluting the namespace with symbols related only by similarity of name.  Not a purely practical disadvantage, perhaps, but one that I’m sure grates against the nerves of anyone looking for neat, elegant answers.

Fortunately, as I mentioned briefly in my last post, the answer already exists, and has already been added to GHC.  The answer here is view patterns.  In data-accessor package language, one can actually write the following to express part of the above.

foo ((^. field1) -> Just k) = k

(Admittedly, one does obtain a warning about pattern matching overlaps when doing this.  I’m unsure why.)

Unfortunately, this very nice approach doesn’t appear to currently extend as nicely to pattern matching on several fields at the same time.  A view pattern has one view, which matches an argument; a second view pattern matches a second argument, which isn’t what we want.  The best we can do is to define a view that selects several fields…

foo ((^. field1) &&& (^. field3) -> (Just k, (x:_))) = 4*k

Barely passable for two fields… but add a third field, and suddenly you have to know about the associativity of &&& (which, by the way, we took from the Control.Arrow package).  Besides, the distance and way that conditions on different fields are mixed together isn’t terribly appealing.  We’d really like to just be able to pattern match the same argument multiple times, requiring that all of them succeed.

Ironically, this is hardly a new problem.  It’s actually the same thing encountered in @-patterns.  One wants to pattern match twice there, too — once to bind a name to the entire value, and again to deconstruct it and bind names or match pieces.  Without view patterns, that was the only situation in which one might want to pattern match the same value in several ways; the only choice was to either open up a data constructor and match on it, or else don’t.  With view patterns, though, the possibility arises any number of ways.  (Data.FingerTree contains both the functions viewl and viewr, for instance; I’m not personally aware of an application that would like to pattern match on both at once, but it doesn’t seem terribly hard to imagine that it might occur.)

It makes sense to propose, in lieu of trying to invent a way to extend record pattern-matching syntax to a new kind of records, to simply provide a mechanism for simultaneously pattern-matching the same argument in several ways.  What we would obtain from this is to generalize the problems solved by @-patterns and record syntax at once, as well as potential future as-yet-unidentified problems.

In other words, the real feature we want is simultaneous patterns — multiple patterns that are matched side-by-side with a single value.  The semantics are that the pattern match succeeds only if all of them do, in which case all of the variables occurring in the pattern are appropriately bound.

In my last post, I proposed using a semicolon to separate the patterns.  Here’s the previous example in this form, using data-accessor syntax.

foo :: MyRecord -> Int
foo ((^.field1) -> Just k ; (^.field2) -> 4 ; (^.field3) -> (x:_)) = 4*k - x

With that, we’d have solved one of the outstanding issues in building a library replacement for Haskell’s distinguished record syntax.  Essentially, there’d just be no more need for a “special” record syntax for pattern matching.  We’d have done it not by adding a new distinguished record syntax, but rather by opening up the existing pattern matching features (including view patterns, of course) to accomodate a rather straight-forward extension.

If I don’t see an obvious reason not to pursue this, I may try to determine how difficult it would be to implement this as a GHC extension.

Playing With Records

Haskell has an existing record syntax that works for some purposes, and doesn’t work very well for some others.  For this blog entry today, I’ll try to build a better alternative as a library within the language.  It will have some advantages (mainly in flexibility), and some disadvantages (mainly in difficulty of type checking, and in having poorer syntactic niceties).  I’m also about 100% certain that this isn’t the first time someone has done these things; but perhaps it’s the first time someone has described them in this way, or perhaps you can discuss existing systems in the comments.

The Current System

Here’s an example of the existing record system:

data Employee = Employee { name :: String, salary :: Double }

That defines a new type for employees, and two named fields called name and salary, each of which has a different type.  I can define an employee in one of two ways:

chris = Employee "Chris" 0.01
bob   = Employee { name = "Bob", salary = 1000000 }

The first one uses positional values, and the second uses named values.  Both define exactly the same sort of value, and I can use the field names to access their properties

name chris {- The result is "Chris" -}
salary bob {- The result is 1000000 -}

The property names can also be used to set values.  For example, I can double my salary.

chris' = chris { salary = 0.02 }

The important thing to note here is that currently, the field names for a record are used in two ways: first, they operate as labels in the syntax that uses braces; and second, they work as functions to get the value of a field.

Building the Alternative

We’ll get started with a few language extensions:

{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE FunctionalDependencies    #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE ExistentialQuantification #-}

To start, we’ll need some data type to represent a property of a record.  The type will depend on the types of the record, and the field.  A property consists of a way to retrieve the value, and a way to change the value.

data Property rec a = Property { set :: a -> rec -> rec,
                                 get :: rec -> a }

Yes, I’m aware of the irony of using the existing record syntax to define a replacement!

At least for backward compatibility purposes, I’d like to retain the ability to use the property name alone as the getter function.  I’ll present this in the order I thought when I wrote it, so hopefully the thought process can be duplicated.  I started with an example.

data Employee = Employee {- name   :: -} String
                         {- salary :: -} Double

First, I know that I want a type class, because that’s the only way that one name can stand for both a function and a named data type like Property.  That type class needs to encode the type of the record, and the type of the field, and also needs a type for the property name to map to.

name   :: Accessor Employee String prop => prop
salary :: Accessor Employee Double prop => prop

When I use one of those variables, I want it to possibly represent a value of the appropriate Property type, or possibly a value of the getter function type.  To define them, though, I want to just start with the property.  So the purpose of this type class is to convert from a Property to either itself, or else a getter function.  Then, as is often needed with multi-parameter type classes, I need some functional dependencies to help the type checker.

class Accessor rec a prop | prop -> rec, prop -> a where
    fromProperty :: Property rec a -> prop

instance Accessor rec a (Property rec a) where fromProperty = id
instance Accessor rec a (rec -> a)       where fromProperty = get

I can now define the variables I declared types for earlier.

name   = fromProperty $ Property setter getter
    where setter n' (Employee n s) = Employee n' s
          getter (Employee n s)    = n

salary = fromProperty $ Property setter getter
    where setter s' (Employee n s) = Employee n s'
          getter (Employee n s)    = s

Using the New Records

Using the new record types starts to look somewhat similar to the old ones.  I can either use “get”, or I can just use the field name as an accessor function.

name chris {- The result is "Chris" -}
salary bob {- The result is 1000000 -}
get name chris {- The result is also "Chris" -}
get salary bob {- The result is also 1000000 -}

Setting a value is a little different.

chris' = set salary 0.02 chris

Not too bad.  This also has the advantage, versus the current record syntax, that it can be partially applied; in other words, it’s a first class function.

rewardLoyalty = set salary 1000000

If I want to set multiple fields at once, though, it gets a little uglier.

bobette = set name "Bobette" (set salary 2000000 bob)

Hmm, definitely not pretty.  Fortunately, since properties are a first-class value, we can fix this now.  I’ll exploit the very convenient fact that := is a constructor in Haskell.  I’ll also need existential types to contain the variation in types between different properties of an object.

data NameValuePair rec = forall a. Property rec a := a

setAll :: [NameValuePair rec] -> rec -> rec
setAll pairs r = foldl setOne r pairs
    where setOne r (p := v) = set p v r

and now the earlier example looks much nicer.

bobette = setAll [ name := "Bobette", salary := 20000 ] bob

It’s a nice property that we were able to solve the problem within the language just by virtue of having a first-class representation of the record in the language.

Definition of a new record type is a tad more complex.  The code above to do so for Employee was very ugly, to say the least.  However, this is the one area where we’d simply expect some kind of language support to help out.  We could build such a thing using Template Haskell.  In fact, it would very nearly not cause any kind of incompatibility if the existing record syntax were just modified so that instead of producing the field names as functions, it produced them as polymorphic Accessor types as described above.

Limitations

This is certainly not perfect, as record systems go.  The Accessor type class requires all manner of type system extensions (namely, MPTCs, fundeps, and flexible instances and contexts), and even so, I can’t currently convince it to fully realize all of the types.  For example, take this exchange with GHCi:

*Main> :t name chris
name chris :: (Accessor Employee String (Employee -> t)) => t

Even though the only type that could possibly be used there is String, GHCi fails to recognize this, and gives us a long type involving the Accessor type class.  Hopefully, there’s a way around this within the type system, which I’ve merely missed.  Without caution, though, the code above could lead to lots of unintended polymorphism that could result in a serious performance hit.  In fact, there’s perhaps quite a strong argument to be made that if a new record system is being defined, perhaps one ought to eschew backward compatibility, or treat it as deprecated and eventually to be removed… and instead require that everyone write “get name” instead of merely “name”.  This eliminates most of the really weird type stuff that’s going on here, and thereby improves the error messages considerably.  On the other hand, it feels a but unusual to write verbs like that in a functional programming language.

Pattern matching is another very serious limitation.  Since the individual values in a record can’t be accessed except by evaluating functions, they can’t be referred to in a pattern matching definition of a function.  While I tend to never do this with the existing record syntax anyway, it remains a serious limitation.  It’s worth mentioning, though, that view patterns alleviate this difficulty when pattern matching on just one field value.  With multiple field values, view patterns can be used in conjunction with &&& from the Control.Arrow module to accomplish the same thing; it’s just ugly.  It would be nice to see a minimal extension of the view pattern feature to get around this by allowing multiple views to be applied side-by-side to the same parameter.  Perhaps something like (and no, I haven’t checked what this does to grammar conflicts):

displayName :: Employee -> String
displayName (name -> n ; salary -> s) | s > 100000 = map toUpper n
                                      | otherwise  = map toLower n

GHC also provides language extensions for treating record field names as unambiguous in places where the intended field can be determined by the types of values in pattern matching.  Indeed, there have been other proposals for type-directed resolution of ambiguous types; but I’ve always seen it mixed in with dot-record syntax.  This would need to be a language extension still.  Record punning and wildcards, similarly, don’t look to be definable within a library.  Whether they are worth recovering is a different discussion on which I know people have varying opinions.  Perhaps, though, punning at least could be viewed as a feature of view patterns rather than of records.

More Fun With the Property Type

There are two more interesting things we can do with the Property type defined earlier.

One of the features of record syntax in Haskell is that data types with multiple data constructors can still be records, and different constructors can provide fields by the same names.  It’s not hard to see that the same can happen here.

data Person = Customer {- name :: -} String {- , creditRating :: -} Double
            | Employee {- name :: -} String {- , salary :: -} Double

name = fromProperty $ Property setter getter
    where setter n' (Customer n c) = Customer n' c
          setter n' (Employee n s) = Employee n' s
          getter (Customer n c)    = n
          getter (Employee n s)    = n

creditRating = fromProperty $ Property getter setter
    where setter c' (Customer n c) = Customer n c'
          getter (Customer n c)    = c

salary = fromProperty $ Property getter setter
    where setter s' (Employee n s) = Employee n s'
          getter (Employee n s)    = s

This behaves the same way as the corresponding instance would in normal Haskell record syntax.  The treatment of properties as first-class language-definable values, though, opens up additional possibilities.  For example, take this definition of a complex number:

data MyComplex = MyComplex {-  real :: -} Double {- , imag :: -} Double

real = fromProperty $ Property setter getter
    where setter r' (MyComplex r i) = MyComplex r' i
          getter (MyComplex r i)    = r

imag = fromProperty $ Property setter getter
    where setter i' (MyComplex r i) = MyComplex r i'
          getter (MyComplex r i)    = i

That’s all the same so far… but, one can also manipulate the magnitude and argument of a complex number as properties.

magnitude = fromProperty $ Property setter getter
    where setter m c             = let a = get argument c
                                   in MyComplex (m * cos a) (m * sin a)
          getter (MyComplex r i) = sqrt (r^2 + i^2)

argument = fromProperty $ Property setter getter
    where setter a c             = let m = get magnitude c
                                   in MyComplex (m * cos a) (m * sin a)
          getter (MyComplex r i) = atan2 i r

We’ve now defined virtual properties of this record type.  (Yes, it seems silly to do it this way, but mainly just because it’s silly to define complex numbers as a record type to begin with.)

Thoughts on Hackage and the Haskell Platform

I just watched this video of the discussion on the future of Haskell from the Haskell Symposium.  I wanted to write a response on a series of points made about Hackage, the Haskell platform, and the need for a Debian-model “unstable” or a development level of releases.

When comparing Hackage and the Haskell Platform to the Debian model, it’s important to realize that Debian does two things, and there’s no compelling reason, in general, to intertwine them.

• The first is packaging.  Packaging includes putting everything together so that it’s all available from the same place, declares its dependencies on other packages, and is easily downloaded and installed onto the system.  Haskell’s solution to packaging, of course, is Cabal, Hackage, and cabal-install.  It’s an excellent system — easy to use, and quite versatile.  (Debian’s system, apt, is also excellent, by the way.)

• The second issue is integration testing and quality assurance.  This is the part that says we make sure that all the different packages work together nicely, don’t have incompatible version dependencies, that A doesn’t depend on some old behavior of B that recently got changed, and so on.  In the Haskell world, GHC used to do an ad hoc job of this for a few commonly used packages, in the form of extralibs… and now the Haskell Platform is attempting to solve these issues.

People talking about adopting a Debian-like model are, in general, talking about the quality assurance picture.  It’s important to keep in mind the role played by various pieces in this correspondence.  For instance, Hackage — from a quality assurance standpoint only — is actually playing the role of upstream.  Don’t look at Hackage as if it were analogous to Debian unstable.  It’s not!  Packages may get put onto Hackage with no intention of ever becoming a part of the Haskell Platform.  Uploading a package to Hackage does not incur any sort of continuing obligation to maintain that package for future GHC versions.  It just means that someone might see my description and it’s Cabal properties, and that should they want to install my package they can do it easily with cabal-install.  It’s rather important, to me, that Hackage remain that way.  Otherwise, we may sadly see developers hesitating to put their packages on Hackage for fear of the maintenance obligation.  That would be quite a disaster.  (As Duncan mentioned in the video linked above, a sort of non-centralized community moderation may be a good thing in Hackage; e.g., ranking packages by downloads or dependencies.  These would be for user information only, and would not negatively impact someone who uploads a package.  User ratings are a tougher question; they would need to be on a per-version basis to avoid penalizing developers for uploading early… and that may well make ratings useless in the long run, as any new release would essentially wipe them clean.)

Where the Debian model is more worthy of discussion, then, is not so much in Hackage, but rather in the Haskell Platform.  There it might actually make some sense, off in the future, to separate an unstable versus a stable version of the Haskell Platform.  (I doubt we’ll see so much activity in the near future that we have any need for an intermediate testing stage as Debian has now, especially since many of the basic packaging problems can be resolved in Hackage prior to accepting the packages into the Haskell Platform anyway.)  A beta release of the Haskell Platform could well be something worth doing in the long term.

Ultimately, though, I think the main thing is to remember that: (a) the Haskell Platform is far and away better than what we had before this point, and (b) a lot of smart people will be watching what happens, and the community will learn a LOT about quality assurance and blessing of collections of libraries in the next year, and that will help the community make better decisions going forward.  So I guess, in the end, I think things are going in the right way, and hope people continue to keep their eyes and ears open, and be patient so that we can learn from the process instead of going in assuming that we’ve got all the answers.

Type Classes With An Easier Example

I recently wrote this post.  It contained a lot of rather obtuse mathematics, just to introduce what was basically the example problem for the article.  That’s because it was a real honest-to-goodness problem I was solving and writing papers on for algebra journals… admittedly, not the best choice for a blog aimed partly at non-mathematicians.  Here, I do something similar, but with a much simpler and more general purpose problem.

Overview

We’ll be playing with converting matrices to an upper triangular form, essentially using Gaussian elimination.  As a reminder, here’s what that means (I’m simplifying a little):

Goal: I have a matrix, and I’d like it to be upper triangular.  (Upper triangular means that everything below the main diagonal is zero.)

Rules: I’ll allow myself to do these things to the matrix: swap any two rows or columns, or add any multiple of a row or column to another one.

Strategy: I’ll look for a non-zero element in the first column.  If there is one, I’ll swap rows to move it to the first row of the matrix.  Then I’ll add multiples of the first row to all of the other rows, until the entire rest of the first column is equal to zero.  Then (here’s where we get a little tricky) I can just do the same thing on the rest of the matrix, ignoring the first row and column.  In other words, I’ve reduced the problem to a smaller version of itself.  (Yep, it’s a recursive algorithm.)

The Trick Up Our Sleeve: Instead of writing our function to operate directly on some representation of matrices, we’ll make it work on a type class.  This will let us play all sorts of cool tricks.

Preliminary Stuff

I’d like to be able to declare instance for my matrices without jumping through newtype hoops, so I’ll start with a language extension.

{-# LANGUAGE TypeSynonymInstances #-}

Imports:

import Data.Maybe

Next, I need a few easy utility functions on lists.  There’s nothing terribly interesting here; just a swap function, and a function to apply a function to the nth element of  a list.

swap :: Int -> Int -> [a] -> [a]
swap i j xs | i == j    = xs
            | i > j     = swap j i xs
            | otherwise = swap' i j xs
    where swap'  0 j (x:xs) = let (b,xs') = swap'' x (j-1) xs in b : xs'
          swap'  i j (x:xs) = x : swap' (i-1) (j-1) xs
          swap'' a 0 (x:xs) = (x, a:xs)
          swap'' a j (x:xs) = let (b,xs') = swap'' a (j-1) xs in (b, x:xs')

modifynth :: Int -> (a -> a) -> [a] -> [a]
modifynth _ _ []     = []
modifynth 0 f (x:xs) = f x : xs
modifynth n f (x:xs) = x : modifynth (n-1) f xs

Finally, I need a type to represent matrices.  Since this is just toy code where I don’t need really high performance, a list of lists will do just fine.

type Matrix = [[Double]]

All done.  On to the interesting stuff.

Building a Type Class

I already mentioned that I don’t want to operate directly on the representation of a matrix as a list of lists.  Instead, I’ll declare a type class capturing all of the operations that I’d like to be able to perform.

class Eliminable a where
    (@@)     :: a -> (Int,Int) -> Double
    size     :: a -> Int
    swapRows :: Int -> Int -> a -> a
    swapCols :: Int -> Int -> a -> a
    addRow   :: Double -> Int -> Int -> a -> a
    addCol   :: Double -> Int -> Int -> a -> a

I’ve reserved the operator @@ to examine an entry of a matrix, size to give me its size, and then included functions to swap rows, swap columns, add a multiple of one row to another, and add a multiple of one column to another.  Now I just need an implementation for the concrete matrix type I declared earlier.

instance Eliminable Matrix where
    m @@ (i,j)     = m !! i !! j
    size m         = length m

    swapRows p q m = swap p q m
    swapCols p q m = map (swap p q) m
    addRow k p q m = modifynth q (zipWith comb (m!!p)) m
        where comb a b = k*a + b
    addCol k p q m = map (\row -> modifynth q (comb (row!!p)) row) m
        where comb a b = k*a + b

Done.

Programming With Our Type Class

Recall that the substantial portion of the elimination algorithm earlier was to zero out most of the first column of the matrix, leaving only the top element possibly non-zero.  We’re now in a position to implement this piece of the algorithm.  It’s not all that tricky.

zeroCol :: Eliminable a => a -> a
zeroCol m = let clearRow j m' = addRow (-(m' @@ (j,0) / m' @@ (0,0))) 0 j m'
                clearCol m'   = foldr clearRow m' [1 .. size m - 1]
            in  case listToMaybe [ i | i <- [0 .. size m - 1],
                                       m @@ (i,0) /= 0 ]
                of   Nothing -> m
                     Just 0  -> clearCol m
                     Just i  -> clearCol (swapRows 0 i m)

This function looks only at the first column of the matrix, and clears it out my moving a non-zero element to the top, and then adding the right multiple of that column to all those below it.  The important thing to notice is that this was implemented for any arbitrary instance of the type class I called “Eliminable”.  This will be incredibly useful in the next few steps.

Using the Type Class

The next part of the algorithm is to ignore the first row and column, and perform the same operation on the submatrix obtained by deleting them.  It’s actually a bit unclear how to implement this.  We have a few options:

1. Modify zeroCol above, to have it take a parameter representing the current column, and do everything relative to the current column.  This is pretty messy.  It actually might not be too messy in this case, but if the algorithm I were implementing were a little less trivial to begin with, It might definitely be quite messy.
2. Actually perform the elimination on a separate matrix, and then somehow graft the first row and column from this matrix onto that one.  Again, this could get pretty messy in general.
3. Change the representation of the submatrices.

I’ll choose the third.  Luckily, this isn’t too tough, since we have a type class.  I’ll just define a newtype, and a new instance, to encapsulate the idea of a matrix with the first row and column deleted.

newtype SubMatrix a = SubMatrix { unwrap :: a }

instance Eliminable a => Eliminable (SubMatrix a) where
    (SubMatrix m) @@ (i,j) = m @@ (i+1,j+1)
    size (SubMatrix m)     = size m - 1

    swapRows p q (SubMatrix m) = SubMatrix (swapRows (p+1) (q+1) m)
    swapCols p q (SubMatrix m) = SubMatrix (swapCols (p+1) (q+1) m)
    addRow k p q (SubMatrix m) = SubMatrix (addRow k (p+1) (q+1) m)
    addCol k p q (SubMatrix m) = SubMatrix (addCol k (p+1) (q+1) m)

Using this new instance, I can easily complete the elimination algorithm.

eliminate :: Eliminable a => a -> a
eliminate m | size m <= 1 = m
            | otherwise   = unwrap . eliminate . SubMatrix . zeroCol $ m

Yep, that’s all there is to it, and we have a working elimination algorithm.

Type Class Games

Suppose, now, that I want a lower triangular matrix.  It might initially seem that I’m out of luck; I need to write all this code again.  That turns out not to be the case, though.  If I just teach the existing code how to operate on the transpose of a matrix instead of the matrix I’ve given it, then all is well!  Here goes.

newtype Transposed a = Transposed { untranspose :: a }

instance Eliminable a => Eliminable (Transposed a) where
    (Transposed m) @@ (i,j) = m @@ (j,i)
    size (Transposed m)     = size m

    swapRows p q (Transposed m) = Transposed (swapCols p q m)
    swapCols p q (Transposed m) = Transposed (swapRows p q m)
    addRow k p q (Transposed m) = Transposed (addCol k p q m)
    addCol k p q (Transposed m) = Transposed (addRow k p q m)

To implement the lower triangular conversion, now, is simple.

lowerTriang :: Eliminable a => a -> a
lowerTriang = untranspose . eliminate . Transposed

Any number of changes to the operation we’re trying to perform can often be expressed by simply substituting a different representation for the type on which we’re performing the operation.  (Thinking about this fact can actually get pretty deep.)

Side Calculations

There’s a fairly common problem that many people run into when moving from an imperative language to a functional one.  This can apply to learning functional programming, converting existing imperative code, or even just translating the concepts in one’s mind when talking to someone who thinks imperatively.  The problem goes something like this: you have some code that performs some computation, and now you want to change the code to add some new concept to the existing computation.  Often, the new idea you’re trying to add could be performed trivially in an imperative language, by adding print statements to some function seven layers in, or by keeping track of some value in a global variable, or in some field of some object.  In the functional setting, these aren’t available to you.

The minimalist answer is simply to add all the plumbing code; new parameters and return values, etc. to every function in the entire call tree.  To say the least, this is unappealing!  To a new Haskell programmer, the obvious answer often looks like monads.  However, again, the entire call tree has to be rewritten in a monadic style, and besides, this is a tad like using a rocket launcher to rid the house of mice.

The solution I propose here is that many times, it’s sufficient to use a type class.  Here’s an example.

Problem: Calculate the determinant of a matrix efficiently.

Determinants can be calculated in a lot of different ways, but one of the most common uses elimination.  The interesting fact here is that once you’ve got a triangular matrix (lower or upper; doesn’t matter), then its determinant is just the product of its diagonal elements.  Furthermore, we know precisely what happens to the determinant when you swap rows or columns (it flips sign, but the magnitude stays the same), or when you add a multiple of one row or column to the other (it stays the same).  So a (very fast) way to calculate a determinant is to perform elimination, but also keep track, at each step, of what you’ve done to the determinant so far.

So now we need, not merely a matrix, but a pair consisting of a matrix and some side information – namely, which change we’ve made so far to the determinant.

data WithDeterminant a = WithDeterminant Double a

instance Eliminable a => Eliminable (WithDeterminant a) where
    (WithDeterminant _ m) @@ (i,j) = m @@ (i,j)
    size (WithDeterminant _ m)     = size m

    swapRows p q (WithDeterminant d m) = WithDeterminant (-d) (swapRows p q m)
    swapCols p q (WithDeterminant d m) = WithDeterminant (-d) (swapCols p q m)
    addRow k p q (WithDeterminant d m) = WithDeterminant d    (addRow k p q m)
    addCol k p q (WithDeterminant d m) = WithDeterminant d    (addCol k p q m)

As before, once we’ve defined the appropriate instance, the implementation is actually quite easy.

diags :: Matrix -> [Double]
diags []     = []
diags (r:rs) = head r : diags (map tail rs)

determinant :: Matrix -> Double
determinant m = let WithDeterminant d m' = eliminate (WithDeterminant 1 m)
                in  d * product (diags m')

The resulting determinant function actually performs quite well.  For example, calculation of the determinant of a 100 by 100 matrix is done in 1.61 seconds, fully interpreted in GHCi.  I didn’t bother compiling with optimization to see how well that does, nor replacing the inefficient list-of-lists representation of matrices with one based on contiguous memory or arrays.  (Edit: Compiled and optimized with GHC, but still using the list-of-list representation, the time is around a third of a second.)

(It’s worth pointing out that automatic differentiation is another very impressive example of this same technique, except that it uses the standard numeric type classes instead of a custom type class.)

Conclusion

The point of this article is that plain old type classes in Haskell can be used to make your purely functional code very flexible and versatile.  By defining a type class to capture the concept of a set of related operations, I was able to achieve:

1. Choice of data structures.  Had I wanted to use a contiguous array instead of a list of lists, I could have easily done so.
2. Easier programming.  For example, operating on a submatrix of the original matrix became much easier.
3. Flexible code.  I was able to get lower triangular matrices, too, without rewriting the code.
4. Better composability.  I easily reused by upper triangular matrix calculation to find determinants, even though additional calculations had to be traced through the original.
5. Separation of concerns.  When I started, I never even dreamed that I might need to trace determinant calculations through the process.  That got added later on, in it’s own separate bit of code.  If someone else wanted different plumbing… say, for logging, or precondition checking, or estimating the possible rounding error, all they’d need to do is define a new instance of the type class.

None of this is new, of course.  But type classes are definitely an underused language feature by many Haskell programmers.

On Inverses of Haskell Functions

Recall that given a function , a function is called a left inverse of f in case .  In other words, we want g to undo whatever f did to its parameter.  A function has a left inverse if and only if it is injective, a.k.a. one-to-one.

Question: How can we obtain left inverses of Haskell functions automatically?

In other words, I’d like to have a function

f :: a -> b

and obtain, without explicitly writing it, another function

g :: b -> a
g . f = (id :: a -> a)

I might first set out to accomplish the gold standard; to do this automatically, for an arbitrary domain.  Alas I wouldn’t get very far.  To succeed in that task would imply that the proposition is constructively valid.  But it isn’t even classically valid, so that’s hopeless.  Of course, if the domain of f is enumerable, I can in principle simply try all inputs to f in parallel until I find one that works… this is not terribly satisfying, however; it’s completely impractical for most purposes.

So what if I restrict the domain of f further? In particular, I will require that f operate polymorphically on values of some type class, and then restrict the type class to ensure that I can get inverses easily.  By requiring that f be polymorphic, I can ensure that only “good” (i.e., easily invertible) things happen in the implementation of f.  (By the way, we need rank-2 types for this.)

class Fooable a where
    foo :: Int -> a -> a

foo' n = foo (-n) -- an inverse for (foo n)

I now proceed to implement an automatic inverter for functions defined from and to Fooable types.

newtype FooInversion = FooInversion { unInversion :: Fooable a => a -> a }

instance Fooable FooInversion where
    foo n (FooInversion inv) = FooInversion (inv . foo' n)

invertFooable :: Fooable a => (forall b. Fooable b => b -> b) -> a -> a
invertFooable f = unInversion (f (FooInversion id))

Finished!  I did something a little tricky there, I suppose — I defined functions from Fooable things to Fooable things to actually be Fooable things themselves.  Not too unusual a trick to play in abstract algebra, really.

-- Declare a Fooable instance for testing
instance Fooable Int where foo = (+)

-- Define a QuickCheck property to ensure it's working
prop a b = invertFooable (foo b) a == a - b

Unfortunately, while this worked out, it’s fairly fragile.  It turns out that Fooable isn’t really a terribly useful type class.  The only fully polymorphic functions that can be defined from Fooable things to other Fooable things are actually finite iterations of one operation.  I can add other operations, sure, but the thing I’m missing is any kind of choice.  Since I have no way to inspect the value I’ve got as input to my function, I can’t make good use of if, case, pattern matching, or anything else of that form.  While it might be interesting to characterize precisely what type signatures can exist in the Fooable type class without making inversion impossible, I’m going in a different direction.

Suppose that Fooable contained functions to inspect properties of the value you’re working with.  Then we can’t play the trick above quite as cleanly as one might hope.

class Barrable a where
    look :: a -> Int
    bar  :: Int -> a -> a

bar' n = bar (-n)

Now if we tried to proceed as before, how should we define look in the inversion instance?  Really, we can only get an inverse for specific sequences of operations that we end up in, when computing in the forward direction.  To implement that, we need to know the input (for the forward function) in advance, and compute with it.

data BarInversion a = BarInversion { unBarInversion :: Barrable b => b -> b,
                                     barInvertible  :: a }

instance Barrable a => Barrable (BarInversion a) where
    look (BarInversion f x)  = look x
    bar n (BarInversion f x) = BarInversion (f . bar'') (bar n x)
        where bar'' y = let y' = bar' n y
                        in  if look y' == look x then y' else undefined

invertBarrable :: (Barrable a, Barrable b) => a -> (forall c. Barrable c => c -> c) -> b -> b
invertBarrable x f = unBarInversion (f (BarInversion id x))

This is similar to what was done above, except for two changes:

1. In addition to an inverse function, BarInversion carries around a current value at which the inverse is defined.
2. The inverse function checks to ensure that all the properties of the current value that are observable via the type class are the same, since the original function might have relied on those properties in deciding what to do.  If they don’t match, it gives up.

The intention is that our look function doesn’t really give away the house; i.e., it provides partial, not complete, information about the value.  So for our test case, we’ll only give partial information:

instance Barrable Int where
    look n = n `div` 5
    bar    = (+)

Now:

let f = invertBarrable 5 (bar 2)
f 7 == 5
f 8 == 6
f 9 == 7
f 10 == 8
f 11 == 9
f 12 == *** undefined ***

In other words, as long as the input is observably the same as the point at which we performed the forward calculation, the inverse is available.  If not, then the result is undefined.  We’ve got a partial inverse, at least.

It’s worth noting that even if the inverse were only defined at the one point we started with, the inverse function could be valuable.  Note that the type of invertBarrable only requires that x be of some Barrable type, and that the inverse operation on som Barrable type.  They need not actually be the same type!  This is quite useful if, for example, you want to trace some other computation through the calculation of the inverse.  I think I’ll write another blog post on the practice of tracing one calculation through another one using type classes… but, some other time.

Flow Equivalence Code in Haskell

I’ve talked about this code for a few blog entries… so in case anyone was interested in seeing the end result:

http://patch-tag.com/r/flowequiv/snapshot/current/content/pretty

There it is.  The main function of interest is canonicalForm, defined at the very end.  This is a computational version of Franks’ construction of canonical forms of flow equivalence classes.  Sorry about the GPL license; it was to avoid an argument, and it’s unlikely this would ever be commercially useful anyway.

Catching a Mathematical Error Using Haskell’s Type System

I found this worth sharing.  It’s the story of an error in a mathematics paper.  Perhaps not too terribly surprising, you might think.  But it really is surprising.  The paper in question is quite well-written and clear.  It’s central to its field, and has been read in detail by, I’m sure, hundreds if not thousands of people.  It’s 25 years old.  It doesn’t involve a lot of references to other results, but rather builds its entire argument right there in this one paper.  In other words, it’s exactly the sort of mathematics that you’d expect to be correct.  (The error, by the way, is not major and is easily overcome, and the result is still true.  So I’m not claiming that 25 years of research is invalid or anything like that!)

Here’s the story.  John Franks published a central result in the theory of subshifts of finite type, in 1983.  Here’s everything you need to know about that field:

Subshifts of finite type can be described by directed graphs (loops and parallel edges allowed.)  There are some interesting operations that we can perform on these graphs: in-splitting, out-splitting, in-amalgamation, out-amalgamation, expansion, and contraction.  These operations arise in a really cool way, but you don’t need to worry about that.  The main point, for the moment, is that although they come about quite naturally, they are not easy to work with directly.

Franks noticed that a certain matrix of integers seems to come up quite a bit in understanding the structure of these subshifts: that matrix is A – I, where A is the adjacency matrix of the graph, and I is the identity matrix of the right size.  (In other words, the (i,j) entry of A is the number of edges from vertex i to vertex j.)  He also noticed that he can perform row and column operations on that matrix — adding or subtracting a row or column to/from another — by doing those other fancy operations above in a certain way.  This was huge.  Now if I want to get from one graph G to another graph H, using those splits and amalgamations and expansions and such, it actually suffices to get from the matrix A – I to the matrix B – I (where A and B are the adjacency matrices of G and H) by using elementary row and column operations!

(There are a few technicalities, but skip this paragraph if you don’t care about them.  You first have to arrange for G and H to have a loop at every vertex, so that A – I and B – I are non-negative.  Then you have to make sure that your row and column operations don’t leave you with negative entries in your matrices.  These are handled correctly in the paper, but are not necessary for the point I’m making, so I won’t mention them again.)

Then Franks sets out to do this.

First, he shows that he can arrange so that the entire first column of the matrix A – I is equal to the greatest common divisor of all the entries.  Starting from this, he uses an elimination algorithm similar to what you often see in introductory linear algebra courses: you clear out some row and column except for the single cell where they coincide, and then ignore that row and column and recursively perform the same operation on the remaining smaller matrix.  Franks chooses, in particular, to clear out the first column and the second row.

But here’s the error.  Elementary row and column operations on this sub-matrix carry over to the larger matrix… but when arranging for the first column to be equal to the g.c.d. of all of the entries, Franks accidentally used some different operation — namely, conjugation by a permutation matrix.  A seemingly trivial operation, since it corresponds to a graph isomorphism.  In a sense, it doesn’t even change the graph at all, so it’s hardly even an operation.  However, when you delete the first column and second row, and then conjugate the resulting matrix by a permutation matrix, this does not correspond to any kind of good operation on the larger graph.  The proof of this central result in the theory of subshifts of finite type, therefore, contains an error.

How did I find this?  Actually, I found it because I was implementing Franks’ construction in Haskell, and Haskell has a type system that was remarkably helpful in reasoning about the result.  I started by writing a type class (called FlowEquiv) for representations of graphs on which I can perform those weird operations: splits, amalgamations, etc.  As I described in my last blog entry, I was able to capture a lot of higher level concepts by writing wrapper types that are instances of this type class.  For example, these operations can be performed on the transpose of a graph, and when they are, they correspond to other operations on the original graph.  I was able to write down that correspondence explicitly by having a wrapper type called Transposed implement the same type class.

Naturally, when I came to the point of needing to delete the first column and second row, I looked for a way to do the same thing.  I couldn’t do it.  I poked at this for a couple days, and reached the conclusion that, in fact, there is no obvious correspondence between these operations (splits and amalgamations and such) on a sub-matrix, and the corresponding operations on the larger matrix.  Clearly, a new abstraction was needed, so I defined a new type class, ReducibleMatrix, defining just those operations that carry over nicely to sub-matrices: namely, the row and column operations.  Then the huge insight of Franks, the ability to compose splits and amalgamations and expansions and the like to form elementary row and column operations, corresponds to this line of code:

instance FlowEquiv a => ReducibleMatrix a where

and Franks proof is realized by the definitions of the functions from ReducibleMatrix.  So then I set out to implement those operations of Franks that need to be performed on submatrices in terms of the ReducibleMatrix type class.  And that’s when I saw it.  I could not convert my implementation of that one piece (turning the first column of A – I into the gcd of the entries) on the ReducibleMatrix type class in the way that Franks described it.  Why?  Because of that conjugation by a permutation matrix.  In a sense, Haskell’s type system found the error in Franks’ proof.

Epilogue: it turns out there’s an easy modification of Franks’ proof that works around this issue.  That corresponded to implementing the g.c.d. piece in Haskell in a different way.  The correction is minor, and doesn’t really have an effect on any other results, so it will appear as a brief remark in an upcoming paper I’m writing on Leavitt path algebras.

Anyway, I thought that was worth sharing.  In this case, Haskell’s type system did catch not just a programming error, but an error in a mathematical result that had to the best of my knowledge gone undetected by likely hundreds of people for 25 years.

The Magic of Type Classes

Type classes are amazing.  Just wanted to point that out.  I don’t even mean magical stuff like multiparameter type classes or fancy newfangled nonsense kids these days are using.  (I will, however, use FlexibleInstances in this blog entry, just because it lets me avoid a bunch of type conversions.)  But really, I just mean ordinary old type classes.  Some of the least interesting uses of type classes in Haskell are, for example, the Num and Fractional and such classes and their instances in the standard library.  To be sure, they are sort of interesting — they take things that are part of the language definition in other languages (the “/” operator works on either single or double precision floats, for example), and describe them within the language itself.  That’s cool,but it’s nowhere near the limit of how cool type classes can be.

There’s a much better example of how cool type classes can be in a blog post I wrote some time ago (and the dozens of other people who have done the same thing, and probably blogged about it, too).  But I think this one, while a little more complicated, is even cooler.  Here’s the problem I’m working on this morning.  I’m going to oversimplify shamelessly, in hopes of getting to the point some time in the foreseeable future.  I hope no one reads this who knows the real mathematics involved, because I’ll get angry comments and corrections.

Some people in a branch of mathematics called symbolic dynamics are occasionally interested in certain square matrices of non-negative integers (if one feels like it, one could say they are interested in graphs, where the matrices are just the adjacency matrices of those graphs; it’s a distinction without a difference).  In particular, they are sometimes interested in whether it’s possible to change one of these matrices into a different one by making only certain kinds of changes, of which they are rather fond.  They’ve even given these changes cutesy nicknames: in-splitting, in-amalgamation, out-splitting, out-amalgamation, expansion, and contraction.  When matrices can be converted into each other using these operations, they are called “flow equivalent”.

These operations aren’t hard to implement.  Here are type signatures for functions to implement the operations.  If we represent matrices as lists of lists and allow the use of a few simple utility functions on lists, the actual implementations are all one or two lines long, but are not really relevant here.  The point is that that each of these operations takes some values, and a non-negative square integer matrix, and gives back another non-negative square integer matrix.

    type Vertex = Int

    outsplit   :: Vertex -> [Vertex] -> [Vertex] -> Matrix Int -> Matrix Int
    insplit    :: Vertex -> [Vertex] -> [Vertex] -> Matrix Int -> Matrix Int
    outamalg   :: Vertex -> Vertex -> Matrix Int -> Matrix Int
    inamalg    :: Vertex -> Vertex -> Matrix Int -> Matrix Int
    expand     :: Vertex -> Bool -> Matrix Int -> Matrix Int
    contract   :: Vertex -> Vertex -> Bool -> Matrix Int -> Matrix Int
    dropVertex :: Vertex -> Matrix Int -> Matrix Int

And those are called the “generators” of flow equivalence.  (I feel compelled, in the name of defending myself against the hordes of symolic dynamics researchers who will do doubt try me for crimes against their field, to put in here a disclaimer that deleting sources and sinks are not normally considered among the generators of flow equivalence… but for certain technical reasons, it is convenient to include them here.)

Now, here is where the problem arises.  My functions above will take a matrix, perform an in-amalgamation, for example, and tell me what matrix I’ve got back.  But it’s quite rare, actually, that I’m really interested in merely trying out an in-amalgamation to see what happens.  Actually, I’m far more interested in following someone’s construction.  For example, there are some constructions by some smart guys named John Franks and Danrun Huang that both do interesting things.  And here’s the interesting bit – by the very nature of the construction, I already know what the result is going to be.  What I really care about is what sequence of flow equivalence generators I used to get there.  And the functions I’ve written above simply can’t tell me that.

One thought might be to use the Writer monad.  I myself considered such a path, but it turns out that I would have been quite mistaken to do so.  Consider the following:

fiddleWith x = do
    a <- tryOneThing x
    b <- tryAnotherThing x
    if a `betterThan` b
        then return a
        else return b

See what’s going to happen here?  When I calculate the two candidates, I’ll produce the sequence of flow equivalence generators for each of them in turn, even though only one of them ends up coming back as the result of the computation.  This is quite wrong, and in fact will give me non-sensical results.  Since I anticipate eventually trying to optimize the constructions to try and get better flow equivalences, I certainly don’t want to create a situation where trying two different things will give nonsensical results!

What I really want is for the sequence of flow equivalence generators to be carried along with the non-negative integer matrix as I perform the construction.  So I want:

    outsplit   :: Vertex -> [Vertex] -> [Vertex] -> ([FlowOp], Matrix Int) -> ([FlowOp], Matrix Int)
    insplit    :: Vertex -> [Vertex] -> [Vertex] -> ([FlowOp], Matrix Int) -> ([FlowOp], Matrix Int)
    outamalg   :: Vertex -> Vertex -> ([FlowOp], Matrix Int) -> ([FlowOp], Matrix Int)
    inamalg    :: Vertex -> Vertex -> ([FlowOp], Matrix Int) -> ([FlowOp], Matrix Int)
    expand     :: Vertex -> Bool -> ([FlowOp], Matrix Int) -> ([FlowOp], Matrix Int)
    contract   :: Vertex -> Vertex -> Bool -> ([FlowOp], Matrix Int) -> ([FlowOp], Matrix Int)
    dropVertex :: Vertex -> ([FlowOp], Matrix Int) -> ([FlowOp], Matrix Int)

But I anticipate needing to change some other things, as well.  Eventually, I might want to represent my non-negative integer matrices differently — perhaps as a graph, with labelled edges and vertices.  (It turns out I will want this!)  I may decide that lists of lists are quite inefficient, and I need some data structure better suited to my task.  I may have other needs that I can’t even predict.  So, at this point, I decide that instead of rewriting all my types for this particular variation, I could try a type class, so what I actually write is:

class FlowEquiv a where
    matrix     :: a -> Matrix Int
    size       :: a -> Int

    outsplit   :: Vertex -> [Vertex] -> [Vertex] -> a -> a
    insplit    :: Vertex -> [Vertex] -> [Vertex] -> a -> a
    outamalg   :: Vertex -> Vertex -> a -> a
    inamalg    :: Vertex -> Vertex -> a -> a
    expand     :: Vertex -> Bool -> a -> a
    contract   :: Vertex -> Vertex -> Bool -> a -> a
    dropVertex :: Vertex -> a -> a

This is great, because it allows me to proceed with my implementation in little steps.  I’m actually going to build myself four (yes, four!) little instances of my type class, and they are going to look like the following.

instance FlowEquiv (Matrix Int) where
    matrix a = a
    ...

I am still refusing to give you the implementation of this one, not because it’s a secret, but because it would obscure the point.  If you care, feel free to ask me.  But then adding the logging that I need is simple enough.  First, I need a data type to represent these operations.

data FlowOp = OutSplit   Vertex [Vertex] [Vertex]
            | InSplit    Vertex [Vertex] [Vertex]
            | OutAmalg   Vertex Vertex
            | InAmalg    Vertex Vertex
            | Expand     Vertex Bool
            | Contract   Vertex Vertex Bool
            | DropVertex Vertex
    deriving Show

I’ll now define a new wrapper type to add such a list to some other FlowEquiv type.

newtype WithFlowOps a = WithFlowOps ([FlowOp], a) deriving Show

starting :: FlowEquiv a => a -> WithFlowOps a
starting a = WithFlowOps ([], a)

flowOps :: FlowEquiv a => WithFlowOps a -> [FlowOp]
flowOps (WithFlowOps (ops,_)) = ops

And finally, I’ll build myself an instance:

instance FlowEquiv a => FlowEquiv (WithFlowOps a) where
    matrix (WithFlowOps (_,a)) = matrix a
    size   (WithFlowOps (_,a)) = size a

    outsplit v r1 r2 (WithFlowOps (ops,a))
        = WithFlowOps (ops ++ [ OutSplit v r1 r2 ], outsplit v r1 r2 a)
    insplit v c1 c2 (WithFlowOps (ops,a))
        = WithFlowOps (ops ++ [ InSplit v c1 c2 ],  insplit v c1 c2 a)
    outamalg v w (WithFlowOps (ops,a))
        = WithFlowOps (ops ++ [ OutAmalg v w ],     outamalg v w a)
    inamalg v w (WithFlowOps (ops,a))
        = WithFlowOps (ops ++ [ InAmalg v w ],      inamalg v w a)
    expand v d (WithFlowOps (ops,a))
        = WithFlowOps (ops ++ [ Expand v d ],       expand v d a)
    contract v w d (WithFlowOps (ops,a))
        = WithFlowOps (ops ++ [ Contract v w d ],   contract v w d a)
    dropVertex v (WithFlowOps (ops,a))
        = WithFlowOps (ops ++ [ DropVertex v ],     dropVertex v a)

Now, so matter what my FlowEquiv type, I can simply wrap it in WithFlowOps (typically using the “starting” function I’ve defined that starts with an empty list of flow ops), and then perform my construction!  When I get back the result, I can ask it for its list of flow equivalence generators, and peruse them as I like.

But, as television infomercials are fond of saying, THAT’S NOT ALL!  It also turns out that there are some preconditions for performing many of the flow equivalence generating operations. I can only delete a vertex if it’s a source or a sink.  I can only contract an edge if it’s the only edge leaving its source, and the only edge entering its target.  There are more.  So far, I’ve neglected to check those, so that if I do things correctly, everything is fine; but if ‘be made a mistake, I may just get the wrong answer.  I sure don’t like that!  No problem, I’ll just make another instance:

newtype WithErrorChecks a = WithErrorChecks { unwrapErrorChecks :: a } deriving Show

instance FlowEquiv a => FlowEquiv (WithErrorChecks a) where
    matrix (WithErrorChecks a)           = matrix a
    size   (WithErrorChecks a)           = size a

    outsplit v r1 r2 (WithErrorChecks a)
      | outsplitIsOkay           = WithErrorChecks (outsplit v r1 r2 a)
      | otherwise                = error "Invalid outsplit"

    insplit v c1 c2 (WithErrorChecks a)
      | insplitIsOkay            = WithErrorChecks (insplit v c1 c2 a)
      | otherwise                = error "Invalid insplit"

    ...

You get the idea.  Well, that’s sort of cool, too.  I’m able to perform my error checking in one place, rather than having to build it into every possible representation I might try to implement this for.  But here’s something really cool!  It turns out there’s an operation, well known to those symbolic dynamics sorts, called a Franks row move.  It’s really just a sequence of four of those generators… and not hard to express, but it takes some thinking.  The hardest part, really, is keeping track of the indices of all those vertices.

primitiveRowMove :: FlowEquiv a => Vertex -> Vertex -> a -> a
primitiveRowMove v w a = let r1   = replacenth w 1 (replicate (size a) 0)
                             r2   = modifynth  w (subtract 1) (matrix a !! v)
                             a'   = outsplit v r1 r2 a
                             w'   = if w > v then w + 1 else w
                             c1   = replacenth v 1 (replicate (size a') 0)
                             c2   = modifynth  v (subtract 1) (map (!! w') (matrix a'))
                             a''  = insplit w' c1 c2 a'
                             v'   = if v > w then v + 1 else v
                             a''' = contract v' w' False a''
                         in  outamalg v (v+1) a'''

There it is (using some utility functions from elsewhere… don’t worry about the details of what’s going on, unless you like worrying).

Now there’s another operation, called a Franks column move.  This is a pattern that will show up again and again, though.  Every time I can do something to a matrix, I can also do the same thing to the transpose of that matrix!  That simply amounts to replacing in-splitting with out-splitting, in-amalgamations with out-amalgamations, and so on.  I don’t want to implement everything twice, so I’ll automate this process with yet another instance of FlowEquiv.

newtype Transposed a = Transposed { untranspose :: a } deriving Show

instance FlowEquiv a => FlowEquiv (Transposed a) where
    matrix (Transposed a) = transpose (matrix a)
    size   (Transposed a) = size a
    outsplit v r1 r2 (Transposed a) = Transposed (insplit  v r1 r2     a)
    insplit  v c1 c2 (Transposed a) = Transposed (outsplit v c1 c2     a)
    outamalg v w     (Transposed a) = Transposed (inamalg  v w         a)
    inamalg  v w     (Transposed a) = Transposed (outamalg v w         a)
    expand   v d     (Transposed a) = Transposed (expand   v (not d)   a)
    contract v w d   (Transposed a) = Transposed (contract w v (not d) a)
    dropVertex v     (Transposed a) = Transposed (dropVertex v         a)

Now implementing the Franks column move is a piece of cake!

primitiveColMove :: FlowEquiv a => Vertex -> Vertex -> a -> a
primitiveColMove v w = untranspose . primitiveRowMove w v . Transposed

And that is indisputably cool.  In addition to implementing this type class for the underlying type, Matrix Int, we’ve managed to craft three different secondary instances, each of which is solved a very different problem, and each of which saves a fair amount of work by removing some secondary task from the things we need to think about, instead separating it cleanly into a wrapper type.

Calculating Multiplicative Inverses in Modular Arithmetic

I’m sure plenty of people already know this, but I ran into it today, and it’s interesting… so here it is.  A common and very easy result in abstract algebra is that the ring (where is any positive integer) contains a multiplicative inverse for (with ), if and only if .  But what if you want to know what the multiplicative inverse is?  Well, for small values of , you could just try numbers until one of them works.  It turns out that also has a multiplicative inverse (namely, ), so it must also be relatively prime to … that helps when we’re doing the math mentally, but less so with a computer, since every possibility still requires some amount of work.

It turns out there’s a slightly faster way to find multiplicative inverses in , and it ban be done by doing the calculation of by the Euclidean algorithm, and tracing our own calculation through the steps.  For those who may not remember, the Euclidean algorithm for the greatest common divisor depends on the fact that the greatest common divisor of and is the same as the greatest common divisor of and the remainder when is divided by .  So, for example: .  In our case, we start with the assumption that , so ultimately, the process will end with .

That’s all the setup we need.  Here are the two cases for our calculation:

Case 1: .  Clearly, the muliplicative inverse of one is one, so we are done.

Case 2: .  Here, we reduce the problem to the next step in the chain of the Euclidean algorithm.  Restating the problem, we wish to find a such that , for any integer .

Using simple division, we know that we can write , where are integers and .  But, and this is important, we also know that (the next step in the Euclidean algorithm), so we can actually refine our inequality to .  Now, by solving a smaller instance of this same problem, we can find such that .  Another way of putting this is that divides .

Now we are done.  We simply choose .  Then .  By examining each term separately, it’s clear that this is evenly divisible by , so we have , and we’re done.

Here’s some code to compute this in Haskell:

    inverse :: Integral a => a -> a -> a
    inverse q 1 = 1
    inverse q p = (n * q + 1) `div` p
        where n = p - inverse p (q `mod` p)

And there you have it, multiplicative inverses in .