psykotic

> import Data.List
> import Control.Monad

We start with some basic types and functions for vectors and matrices. Nothing there
is particularly interesting, so if you already know this material then you are encouraged
to skip to the section on zippers.

---

Vectors can be scaled by a number and added in pairs when they have the same length.

> vScale :: Num a => a -> Vector a -> Vector a
> vAdd   :: Num a => Vector a -> Vector a -> Vector a

Vectors also have a dot product:

> vDot :: Num a => Vector a -> Vector a -> a

A vector will be represented as a list of numbers.

> type Vector a = [a]

Addition and scaling is then easily defined entrywise.

> vScale a = map (a *)
> vAdd = zipWith (+)

The dot product is computed by multiplying vectors entrywise and summing.

> vDot v w = sum (zipWith (*) v w)

Subtraction can be defined in terms of addition and scaling.

> vSub v w = v `vAdd` ((-1) `vScale` w)

A core operation in algorithms of linear algebra like Gaussian elimination is to
eliminate the first entry in a vector by subtracting an appropriate multiple of
another vector. This is only well-defined when the other vector's leading entry
is non-zero, or there will be a division by zero.

> vEliminate (v:vs) (w:ws) = 0 : ws `vSub` ((w / v) `vScale` vs)

---

A matrix is an mxn rectangular array of numbers. As with vectors, matrices can be
scaled by a number and added in pairs when equal in length.

> mScale :: Num a => a -> Matrix a -> Matrix a
> mAdd :: Num a => Matrix a -> Matrix a -> Matrix a

An mxn matrix can also be transposed to yield an nxm matrix.

> mTranspose :: Num a => Matrix a -> Matrix a

Unlike vectors, matrices can be multiplied, which is generally non-commutative.
An mxp matrix can be multiplied by a pxn matrix to yield a mxn matrix.

> mMul :: Num a => Matrix a -> Matrix a -> Matrix a

For every n, there is a special kind of nxn matrix, analogous to 1 for numbers,
which acts like a left and right identity for multiplication.

> mIdentity :: Num a => Int -> Matrix a

It is natural to think of a matrix in terms of its individual rows or columns.

> mRows    :: Num a => Matrix a -> [Vector a]
> mColumns :: Num a => Matrix a -> [Vector a]

A matrix will be represented as a list of row vectors.

> type Matrix a = [Vector a]

Its representation as a list of rows is then trivial.

> mRows = id

Scaling and adding is a simple matter of scaling and adding the rows.

> mScale a = map (a `vScale`)
> mAdd = zipWith vAdd

To transpose a matrix, transpose the list of lists that is its representation.

> mTranspose = transpose

The columns of a matrix are the rows of the transposed matrix.

> mColumns = mRows . mTranspose

Multiplication is only a little trickier. The entries in the product consist of
every dot product of a row in the left operand with a column in the right operand.

> mMul m1 m2 = [[row `vDot` column | row <- mRows m1] | column <- mColumns m2]

Finally, a finite identity matrix is computed by taking a finite slice of an infinite
identity matrix.

> mIdentity n = (map (take n) . take n) infiniteIdentity
>     where infiniteIdentity = iterate (0:) (1 : repeat 0)

> mNumRows = length

> mNumColumns []     = 0
> mNumColumns (x:xs) = length x

TODO: Talk about the shape invariant. No empty rows, etc.

---

The representation of a matrix as a list of rows automatically induces a data type of
matrix zippers that will turn out be very useful. In the sequel, for brevity 'zipper' will
without further qualification refer to matrix zippers rather than zippers generally.

A matrix zipper is a  four-way block decomposition of a matrix that efficiently supports
incremental motion of the horizontal and vertical dividers in addition to easy access to
the entries straddling the dividers. In diagrammatic form it looks like this:

    / A | B \
    |---+---|
    \ C | D /

The intersection between the dividers, marked by + in the diagram, is called the focus.

Note that despite the symmetry of the diagram, the matrices A, B, C and D are not generally
all the same size. A and B have the same height, as do C and D. Likewise, A and C have the
same width, as do B and D.

We will call A, B, C and D the upper-left, upper-right, lower-left and lower-right blocks
of the zipper, respectively.

As the matrix representation is row-oriented, the induced zipper representation will
likewise be row-oriented, taking the form ((A, B), (C, D)).

> type Zipper a = ((Matrix a, Matrix a), (Matrix a, Matrix a))

One of the most important properties of zippers are their symmetries. The symmetry of a
plain old pair is swapping.

> swap (a, b) = (b, a)

A zipper has two levels of pairing and therefore two kinds of swapping. The first is vertical:

    / A | B \     / C | D \
    |---+---| ==> |---+---|
    \ C | D /     \ A | B /

> swapVertically = swap

The second is horizontal:

    / A | B \     / B | A \
    |---+---| ==> |---+---|
    \ C | D /     \ D | C /

> swapHorizontally ((a, b), (c, d)) = ((b, a), (d, c))

These swapping symmetries can be composed in two ways. The first is vertical then horizontal:

    / A | B \     / C | D \     / D | C \
    |---+---| ==> |---+---| ==> |---+---|
    \ C | D /     \ A | B /     \ B | A /

The second is the opposite, horizontal then vertical:

    / A | B \     / B | A \     / D | C \
    |---+---| ==> |---+---| ==> |---+---|
    \ C | D /     \ D | C /     \ B | A /

They give the same result!

> swapDiagonally = swapVertically . swapHorizontally

It's important to realize that though the diagram might suggest otherwise, 90-degree rotations are not
natural symmetries of zippers. In order to match up the invariants on width and height, we would need
to transpose the individual blocks:

    / A | B \     / C^T | A^T \
    |---+---| ==> |---+-------|
    \ C | D /     \ D^T | B^T /

This is a non-incremental (requiring O(n^2) time) transformation of the zipper's structure,
and therefore not considered a proper symmetry. In a sense this is the underlying reason for why a zipper
should not be considered exactly synonymous with a four-way block matrix.

Thus a zipper has a total of four symmetries. The first is the identity, the trivial symmetry, which
does nothing. The next two are horizontal and vertical reflection. The final one is a rotational
symmetry that is associated with the diagonal and can be generated by composing horizontal and vertical
reflection in either order.

Mathematicians know these symmetries as Klein's four-group. The underlying reason why the two orders
of composition for swapping are identical is that a 180-degree and a -180-degree rotation are the same.

It was mentioned earlier that one of the purposes of zippers is to allow incremental motion of the dividers
and hence the focus. We will decompose motion into single steps along the the horizontal and vertical.

Take the case of moving the vertical right by one step. With our representation of zippers this will entail
stripping off the column of B and adding it to A, and similarly for C and D. This is fundamentally similar to
moving the focus in a one-dimensional list zipper (A, B). There we would uncons a single entry from A and cons
that single entry onto B, assuming A is non-empty. Here our notion of unconsing and consing will have to be
more general because we are not dealing with single elements but with whole columns whose entries are scattered
throughout multiple rows.

For this task a useful function will be unconsing a list into its head and tail and returning them as a pair:

> uncons []     = Nothing
> uncons (x:xs) = Just (x, xs)

We can extend this function over a matrix by unconsing all leading row elements and returning the remainder
rows as a matrix. When we uncons we have to ensure we normalize blocks so we don't leave behind any empty rows,
which would violate the invariant for matrix representation.

> normalizeBlock ([]:m) = []
> normalizeBlock m      = m

> unconsMatrix [] = Nothing
> unconsMatrix m  = do (v, m') <- mapAndUnzipM uncons m
>                      return (v, normalizeBlock m')

The opposite of unconsMatrix is of course consMatrix, which conses a column vector onto a matrix. We add
an infinite supply of virtual empty rows onto the matrix to make sure no column entries are lost. Keeping
in mind what was said about block normalization for unconsing, this can be considered block denormalization.

> denormalizeBlock = (++ repeat [])

> consMatrix v m = zipWith (:) v (denormalizeBlock m)

To step (A, B) right, we first check to make sure B is not empty. If it is then we do nothing. Otherwise we
uncons from B and cons onto A.

> goRight' (a, b) = case unconsMatrix b of
>                   Just (v, b') -> (consMatrix v a, b')
>                   Nothing      -> (a, b)

This is now easily extended to zippers pairwise.

> goRight ((a, b), (c, d)) = (goRight' (a, b), goRight' (c, d))

What about goLeft, goDown and goUp? Do we have to take pains and code them separately? This is where the
symmetry group of the zipper comes into play. Let's represent a right move as a diagram.

    / A a | b B \     / A a b | B \
    |-----+-----| ==> |-------+---|
    \ C c | d D /     \ C c d | D /

What happens if we first swap horizontally and then go right?

    / A a | b B \     / B b | a A \     / B b a | A \
    |-----+-----| ==> |-----+-----| ==> |-------+---|
    \ C c | d D /     \ D d | c C /     \ D d c | C /

This is a horizontal mirror image of a left move! Swapping horizontally again undoes the mirroring.

    / A a | b B \     / B b | a A \     / B b a | A \     / A | a b B \
    |-----+-----| ==> |-----+-----| ==> |-------+---| ==> |---+-------|
    \ C c | d D /     \ D d | c C /     \ D d c | C /     \ C | c d D /

This trick of applying a symmetry, then a function, then an inverse symmetry is called conjugation.
In the most general case the symmetry and its inverse will be distinct. Fortunately, all our symmetries
will be involutionary, meaning they are their own inverses.

> conjugate g f = g . f . g

> flipHorizontally = conjugate swapHorizontally
> flipVertically   = conjugate swapVertically
> flipDiagonally   = conjugate swapDiagonally

By flipping goRight we therefore get goLeft for free.

> goLeft = flipHorizontally goRight

All the different motions of a zipper are generated by these single-step motions in cardinal directions.
That conjugating motion by symmetries gives another kind of motion (as opposed to some function that is a
non-motion) can be expressed mathematically by saying that the subgroup of motions is isotropic (or "normal")
within the zipper's group of automorphisms inside of which the zipper's symmetries reside.

That takes care of horizontal motions. What about the vertical motions? Can we get them from horizontal motions
by conjugation? The only thing we haven't tried yet is conjugation by vertical swapping. Unfortunately horizontal
motions and vertically swapping commute for much the same reason that horizontal and vertical motions commute.
Not too surprising when you think about it. The issue is that the subgroup of horizontal motions and the subgroup
of vertical motions are individually isotropic.

We are in fact stuck. Earlier we talked about how 90-degree rotations are not real zipper symmetries. If they
were, we could easily express goDown as the conjugation of goRight by a clockwise 90-degree rotation. Here the
gods of mathematics led us down the right path by excluding these as symmetries. If we had implemented goDown
in terms of 90-degree conjugation of goRight then it would have taken O(n^2) rather than O(n) time because of
the O(n^2) cost of matrix transposition.

Fortunately there is a simple direct and efficient O(n) implementation of goDown. Let's look at a diagram.

    / A | B \     / A | B \
    |---+---|     | c | b |
    | c | d | ==> |---+---|
    \ C | D /     \ C | D /

When we implemented goRight, we only had to uncons the first column, which took a tiny bit of work. Here we
only have to uncons the first row from C and cons it onto A. It's exactly like shifting the focus in a
one-dimensional list-based zipper.

> goDown' (a, b) = case uncons b of
>                   Just (v, b') -> (v : a, b')
>                   Nothing      -> (a, b)

As with goRight', it is easy to extend this to full zippers.

> goDown ((a, b), (c, d)) = ((a', b'), (c', d'))
>     where (a', c') = goDown' (a, c)
>           (b', d') = goDown' (b, d)

Now we get goUp for free by vertical conjugation.

> goUp = flipVertically goDown

That takes care of most of the groundwork we needed to lay down. Only a few more things remain.

We've already seen how important the vertical and horizontal symmetries are in understanding zippers.
Each symmetry axis is associated with an interesting kind of structure, something close to a bialgebra.

The usual kind of multiplication is something that takes two elements and produces a new element:

    (a, a) -> a

What happens if you turn the arrow around?

    a -> (a, a)

You get a comultiplication! Here we are of course dealing with zippers, so it takes this form:

    Zipper a -> (Zipper a, Zipper a)

What's a way to produce two zippers from a zipper? We cut it two! There are two axes and therefore
two comultiplications. For the horizontal axis, this is what it looks like:

                    / A | B \
    / A | B \       \---+---/
    |---+---|  ==>
    \ C | D /       /---+---\
                    \ C | D /

It looks simple and is even simpler to implement.

> splitHorizontally ((a, b), (c, d)) = (((a, b), ([], [])), (([], []), (c, d)))

Vertical splitting is no more difficult.

    / A | B \       / A |\  /| B \
    |---+---|  ==>  |---+|  |+---|
    \ C | D /       \ C |/  \| D /

> splitVertically ((a, b), (c, d)) = (((a, []), (c, [])), (([], b), ([], d)))

What's the natural kind of multiplications associated with these comultiplications? It
seems pretty obvious. If you read the diagrams in reverse and it looks like concatenation.
The only hitch is that this assumes the focuses (the + marks in the diagrams) match up.

If we want to define multiplication in this way for arbitrary zippers (at least of compatible
dimensions) then we will need to somehow normalize the focus before trying to concatenate.
It's clear there are only four possible distinguished focus positions for a zipper: the four corners.

Moving the focus to a given corner is most easily accomplished by moving diagonally until an
empty block in that direction is reached. To go to the upper-right corner, for example,
we would iterate goRight . goUp until we stop making progress. For this a simple utility
function will be useful.

> iterateUntil p f x = if p x then x else (iterateUntil p f) (f x)

It turns out we need only define one of the four functions. We pick goUpperLeft.

> goUpperLeft  = iterateUntil isCorner (goUp . goLeft)
>     where isCorner ((a, b), (c, d)) = all null [a, b, c]

Earlier we discovered to our delight that we could get goLeft for free from goRight by
symmetric conjugations. We were not quite so delighted to discover that we couldn't get the
vertical motions this way, so we were forced to manually define at least one of the vertical motions.
However, in the present case the direction of motion for goUpperLeft is diagonal, and a
diagonal vector is non-invariant under both vertical and horizontal reflections. These
reflections generate all four 90-degree rotations of the diagonal motion!

> goUpperRight = conjugate swapHorizontally goUpperLeft
> goLowerRight = conjugate swapVertically   goUpperRight
> goLowerLeft  = conjugate swapHorizontally goLowerRight

With these functions, we can now easily implement conversion between matrices and zippers.
Conversion from a matrix to a zipper happens by putting the matrix in D's place and filling
A, B and C with blanks.

     /+---\
     || M |
     \|   /

> toZipper m = (([], []), ([], m))

Converting a zipper to matrix conceptually erases the dividers between A, B, C and D.

    / A B \
    \ C D /

Equivalently, it moves the focus to the upper-left corner and returns the lower-right block.

    /+-----\
    || A B |
    \| C D /

> fromZipper z = let (([], []), ([], m)) = goUpperLeft z in m

That was a small but crucial detour. Now we can get back to splitting and concatenating.

To horizontally concatenate two zippers of the same height, we will have to move the
focus of the left zipper one of its two right corners, and conversely for the right zipper.
That leaves a choice of a top and bottom corner. Though a consistent choice of either
would yield a well-defined concatenation, the upper-left corner of a zipper is in a sense
the most canonical corner. That is to say, we are really working with pointed zippers, those
that have a distinguished origin. Thus we will always pick top over bottom and left over right
whenever there is a choice like that to be made.

With those preliminaries out of the way, concatenation is surprisingly easy to implement.

> concatHorizontally zLeft zRight = (([], []), (c, d))
>     where (([], []), (c, [])) = goUpperRight zLeft
>           (([], []), ([], d)) = goUpperLeft  zRight

> concatVertically zTop zBottom = (([], b), ([], d))
>     where (([], b), ([], [])) = goLowerLeft zTop
>           (([], []), ([], d)) = goUpperLeft zBottom

TODO: More about bialgebra properties that intertwine concatenation and splitting.

---

> iterateTimes 0 f = id
> iterateTimes n f = (iterateTimes (n-1) f) . f

TODO: Gauss, Laplace, Gauss-Seidel, etc.