TRANSPOSITION FOR FIXED-SIZED MATRICES

Problem-Solving Club
2nd December 2005


We generalize the approach shown in "Transposition for Three-by-Three
Matrices" to matrices of arbitrary size, but with that size enforced
at the type level. 

This script needs GHC extensions - specifically, GADTs. (So it
requires GHC 6.4 or later.)

> {-# OPTIONS -fglasgow-exts #-}


The types Z and S are uninhabited, but are used to index GADTs (by
size). We call types of the form S(S(Z)) "type-level naturals".

> data Z
> data S n


The type Vector is a polymorphic GADT, indexed by a type-level natural
(the length of the vector) and an element type. Note that the type
system enforces the property that all inhabitants of a type Vector n a
have the same length.

> data Vector n a where
>   VNil :: Vector Z a
>   VCons :: a -> Vector n a -> Vector (S n) a

For simplicity, we allow ourselves to flatten a vector to a list.

> toList :: Vector n a -> [a]
> toList VNil = []
> toList (VCons a x) = a : toList x

> instance Show a => Show (Vector n a) where
>   show = show . toList


Two simple shapely "wholemeal" operations on vectors that we will use
are map and fork on vectors. You can think of both of these as being
instances of "zip with apply", where in each case one of the two lists
is constant.

> vmap :: (a->b) -> Vector n a -> Vector n b
> vmap f VNil = VNil
> vmap f (VCons a x) = VCons (f a) (vmap f x)

> vfork :: Vector n (a->b) -> a -> Vector n b
> vfork VNil a = VNil
> vfork (VCons f fs) a = VCons (f a) (vfork fs a)


Other operations on vectors are less wholesome; they work on single
elements. Indexing is a typical example. For this to be safe, the
index selected has to be less than the vector length (we index from
zero); we specify this at the type level by a indexing parameter. We
therefore need some machinery for specifying and manipulating those
index parameters.

The type Nat is a GADT, indexed by a type-level natural. For size n,
the type Nat n has a single inhabitant: a term constructed from the
constructors Z and S of the same structure as the type index itself.

> data Nat n where
>   Z :: Nat Z
>   S :: Nat n -> Nat (S n)

For example, here is a nine-by-nine matrix, with the numbers from 1..81.

> type Nine = S(S(S(S(S(S(S(S(S Z))))))))
> nine :: Nat Nine
> nine = S(S(S(S(S(S(S(S(S Z))))))))

> viter :: Nat n -> (a->a) -> a -> Vector n a
> viter Z f a = VNil
> viter (S n) f a = VCons a (viter n f (f a))

> mat :: Vector Nine (Vector Nine Int)
> mat = vmap (viter nine (+1)) (viter nine (+9) 1)

Given a matrix of at least one row, we can compute its width:

> width :: Vector (S n) (Vector m a) -> Nat m
> width (VCons x xs) = len x

> len :: Vector n a -> Nat n
> len VNil = Z
> len (VCons a x) = S (len x)

The type LT is a GADT, indexed by a pair of type-level naturals. For
given type-level naturals m and n, the type LT m n has a single
inhabitant: a proof of the fact that m<n. This proof is a term
constructed from Base and Step, which in effect constructs a witness,
namely the difference n-m.

> data LT m n where
>   Base :: LT Z (S n)
>   Step :: LT m n -> LT (S m) (S n)

As well as the type Nat n, which has a single inhabitant n, we need a
type BNat n of bounded naturals, which has inhabitants naturals less
than n. These inhabitants are constructed from proofs of the
inequality.

> data BNat n where
>   Inject :: LT m n -> BNat n

Now we can write the list indexing operation:

> index :: BNat n -> Vector n a -> a
> index (Inject Base) (VCons a x) = a
> index (Inject (Step u)) (VCons a x) = index (Inject u) x


We need to convince the type checker of various properties of
naturals. These proofs take the form of simple coercions. Firstly, any
natural n is less than its successor:

 > ordered :: Nat n -> LT n (S n)
 > ordered Z = Base
 > ordered (S n) = Step (ordered n)

Secondly, successor increases: 

 > weaken :: LT m n -> LT m (S n)
 > weaken Base = Base
 > weaken (Step u) = Step (weaken u)

 Finally, any type of bounded naturals can be injected into the next
 larger such type:

 > include :: BNat n -> BNat (S n)
 > include (Inject l) = Inject (weaken l)


We will also need to reason about sums of type-level naturals.

> data Sum m n p where
>   SumZ :: Sum Z n n
>   SumS :: Sum m n p -> Sum (S m) n (S p)

For example (though we don't actually need this in what follows), sums
allow us to propagate vector lengths through append.

> vapp :: Sum m n p -> Vector m a -> Vector n a -> Vector p a
> vapp SumZ VNil y = y
> vapp (SumS u) (VCons a x) y = VCons a (vapp u x y)

We will need two properties of sum. One is a kind of associativity:
that if m+(1+n)=p then also (1+m)+n=p.

> sumAssoc :: Sum m (S n) p -> Sum (S m) n p
> sumAssoc SumZ = SumS SumZ
> sumAssoc (SumS s) = SumS (sumAssoc s)

The other need is to extract a proof that m<p from a witness that
m+(1+n)=p.

> sumLT :: Sum m (S n) p -> LT m p
> sumLT SumZ = Base
> sumLT (SumS s) = Step (sumLT s)


The final ingredient we will need for matrix transposition is to
compute the vector of predecessors of a type-level natural. We do this
using an accumulating parameter; preds' n m s produces the vector of n
elements counting up from m, making use of the proof s that all these
elements (and in particular the last) are less than an implicit bound
b.

> preds :: Nat n -> Vector n (BNat n)
> preds b = preds' b Z SumZ

> preds' :: Nat n -> Nat m -> Sum m n b -> Vector n (BNat b)
> preds' Z m s = VNil
> preds' (S n) m s 
>   = VCons (Inject (sumLT s)) (preds' n (S m) (sumAssoc s))

(The definition we constructed of this function during the meeting did
not work; this is thebest I have come up with. I'm rather unhappy with
it, because it is somewhat inefficient: both sumLT and sumAssoc take
time linear in the size of the proof s, which grows with each
step. But I suspect that those type-level manipulations are all
necessary, even though they are ultimately unproductive.)

Given predecessors, we can at last generalize the approach we took for
size three:

  tr3 = fork3 (Trip (fmap fst3) (fmap snd3) (fmap thd3))

The functions fst3, snd3, thd3 are the indexing operators for triples;
we generate a vector of such operators for indexing a vector.

> tr :: Vector (S n) (Vector (S m) a) -> Vector (S m) (Vector (S n) a)
> tr x = vfork (vmap (vmap . index) (preds (width x))) x

I have yet to check that the proof that tr is its own inverse goes
through smoothly, though, because I think I have a better plan.