How to create Bounded instance based on runtime values?

Question

I've been playing with O'Connor's matrix implementation based on *-semirings, allowing very neat solutions for graph algorithms:

import Data.Array
newtype Matrix i e = Matrix (Array (i,i) e)

matrix :: (Ix i, Bounded i) => ((i,i) -> e) -> Matrix i e
matrix f = Matrix . listArray (minBound, maxBound) . map f $ entireRange

However, I'd like to read in adjacency matrices of arbitrary sizes from files in the outside world, so having an enumerated type that the matrix is indexed on (like Matrix Node :: Matrix Node2 (Maybe Integer) from the same paper) doesn't really work for me.

My first thought was something like

toMatrix :: [[a]] -> Matrix Int a
toMatrix list = Matrix (listArray ((0,0),(l-1,l-1)) $ concat list)
  where l = length list

but of course this doesn't work either: trying to actually use this matrix blows up when various typeclass instances try to access index (minBound :: Int, minBound :: Int).

Parameterizing the matrix type with a size like

newtype Matrix i e = Matrix i (Array (i,i) e)

doesn't quite work either: although I can change the matrix function to build matrices this way, now I have trouble writing pure for the Applicative (Matrix i e) instance or one for the Semiring (Matrix i e) instance, as the correct one :: Matrix i e depends on the size of the matrices in context.

Conceptually, I can think of two ways out of this:

1. Define a new BoundedInt type with a Bounded instance that can be set at runtime when we know the size of the array, or
2. Find a way to declare instances of Applicative (Matrix i e) parameterized on the size of the matrix.

But I don't know how to implement either of these, and searches around the subject seem to turn up gnarly complicated things. This question also looks relevant, but I don't think it solves the problem (though it would let me use the Bounded i constructor on matrices of fixed Int size).

What's the simplest solution here? Is there one without having to learn how to use the singleton library/some kind of dependent typing?

Answer 1

Today's a bank holiday here in the UK, so I had time to finish that answer about statically sized matrices. I wouldn't necessarily recommend doing this in production - even leaving aside how silly the code is, this is a terrible representation of matrices if you want to do efficient linear algebra on real hardware - but it's kinda fun to mess around with.

From Hasochism:

-- Natural numbers and their singletons in explicit and implicit varieties
data Nat = Z | S Nat  -- page 2 of the paper
intToNat :: Int -> Maybe Nat  -- paraphrased from page 10
intToNat n
    | n < 0     = Nothing
    | n == 0    = Just Z
    | otherwise = S <$> intToNat (n-1)
data Natty n where  -- page 2
    Zy :: Natty Z
    Sy :: Natty n -> Natty (S n)
-- page 3
class NATTY n where
    natty :: Natty n
instance NATTY Z where
    natty = Zy
instance NATTY n => NATTY (S n) where
    natty = Sy natty

-- turn an explicit Natty into an implicit one
natter :: Natty n -> (NATTY n => r) -> r  -- page 4
natter Zy r = r
natter (Sy n) r = natter n r

-- vectors, matrices in row-major order
data Vec n a where  -- page 2
    V0 :: Vec Z a
    (:>) :: a -> Vec n a -> Vec (S n) a
newtype Mat w h a = Mat { unMat :: Vec h (Vec w a) }  -- page 4

-- vector addition, in the form of an Applicative instance
vcopies :: Natty n -> a -> Vec n a  -- page 4
vcopies Zy x = V0
vcopies (Sy n) x =  :> vcopies n x
vapp :: Vec n (a -> b) -> Vec n a -> Vec n b  -- page 4
vapp V0 V0 = V0
vapp (f :> fs) (x :> xs) = f x :> vapp fs xs
instance NATTY n => Applicative (Vec n) where  -- page 4
    pure = vcopies natty
    (<*>) = vapp

-- iterating vectors
instance Traversable (Vec n) where  -- page 4
    traverse f V0 = pure V0
    traverse f (x :> xs) = liftA2 (:>) (f x) (traverse f xs)
instance Foldable (Vec n) where  -- page 4
    foldMap = foldMapDefault
instance Functor (Vec n) where  -- page 4
    fmap = fmapDefault

transpose :: NATTY w => Mat w h a -> Mat h w a  -- page 4
transpose = Mat . sequenceA . unMat

I've taken the liberty of renaming the authors' Matrix type to Mat, rearranging its type arguments, and changing it from a GADT to a newtype. Forgive me for skipping over the explanation of the above - the paper does a better job than I could, and I want to get to the part where I answer your question.

---

Mat w h is an h-vector of w-vectors. It's the type-level composition of two Vec functors. Its Applicative instance, which implements matrix addition, reflects that structure,

instance (NATTY w, NATTY h) => Applicative (Mat w h) where
    pure = Mat . pure . pure
    Mat fss <*> Mat xss = Mat $ liftA2 (<*>) fss xss

as does its Traversable instance.

instance Traversable (Mat w h) where
    traverse f = fmap Mat . traverse (traverse f) . unMat
instance Foldable (Mat w h) where
    foldMap = foldMapDefault
instance Functor (Mat w h) where
    fmap = fmapDefault

We also need a bit of equipment to work with indexes of vectors. To identify a particular element in an n-vector, you have to give a number less than n.

data Fin n where
    FZ :: Fin (S n)
    FS :: Fin n -> Fin (S n)

The type Fin n has exactly n elements, so Fin is the family of finite sets. A value of type Fin n is structurally a natural number less then n (compare FS FZ with S Z), so FS FZ :: Fin (S (S Z)), or FS FZ :: Fin (S (S (S Z))), but FS FZ :: Fin (S Z) will fail to type check.

Here's a higher order function which builds a vector containing all the possible results of its argument.

tabulate :: Natty n -> (Fin n -> a) -> Vec n a
tabulate Zy f = V0
tabulate (Sy n) f = f FZ :> tabulate n (f . FS)

---

Now we can start working with semirings. Taking the dot product of two vectors involves multiplying their elements and then summing the result.

dot :: Semiring a => Vec n a -> Vec n a -> a
dot xs ys = foldr (<+>) zero $ vapp (fmap (<.>) xs) ys

Here's a vector which is zero everywhere except the specified index.

oneAt :: Semiring a => Natty n -> Fin n -> Vec n a
oneAt (Sy n) FZ = one :> vcopies n zero
oneAt (Sy n) (FS f) = zero :> oneAt n f

We'll use oneAt and tabulate to make an identity matrix.

type Square n = Mat n n

identity :: Semiring a => Natty n -> Square n a
identity n = Mat $ tabulate n (oneAt n)

ghci> identity (Sy (Sy Zy)) :: Square (S (S Z)) Int
Mat {unMat = (1 :> (0 :> V0)) :> ((0 :> (1 :> V0)) :> V0)}
-- ┌      ┐
-- │ 1, 0 │
-- │ 0, 1 │
-- └      ┘

And transpose comes in useful for matrix multiplication.

mul :: (NATTY w, Semiring a) => Mat r h a -> Mat w r a -> Mat w h a
mul m n =
    let mRows = unMat m
        nCols = unMat $ transpose n
    in Mat $ fmap (\r -> dot r <$> nCols) mRows

ghci> let m = Mat $ (1 :> 2 :> V0) :> (3 :> 4 :> V0) :> V0 :: Square (S (S Z)) Int in mul m m
Mat {unMat = (7 :> (10 :> V0)) :> ((15 :> (22 :> V0)) :> V0)}
-- ┌      ┐2   ┌        ┐
-- │ 1, 2 │  = │  7, 10 │
-- │ 3, 4 │    │ 15, 22 │
-- └      ┘    └        ┘

So that's the Semiring instance for square matrices sorted. Phew!

instance (NATTY n, Semiring a) => Semiring (Square n a) where
    zero = pure zero
    (<+>) = liftA2 (<+>)
    one = identity natty
    (<.>) = mul

The thing to notice about this implementation is that zero and one dynamically build matrices of a statically known size, typically based on contextual type information at the call site. They get the runtime representation of that size (a Natty) from the NATTY dictionary, which the elaborator builds based on the matrix's inferred type.

This is a totally different approach than that of the reflection library (which I outlined in my other answer). reflection is about stuffing explicit runtime values into implicit instance dictionaries, whereas this style takes information that would otherwise only be known at runtime - the size of the matrix - and makes it static, using singletons to make the type information available in the world of values. Of course a real dependently typed language would dispense with the Natty ceremony: n would be a plain old value and we could use it directly, rather than having to go via a singleton hidden in an instance dictionary.

I'll leave the Kleene algebra stuff to you because I'm lazy and I want to get on to the question of synthesising type information based on runtime input.

---

How can we use these statically sized matrices when we don't know the size statically? You mentioned that your program asks the user how big their graph is (and thus how big the adjacency matrix used to represent the graph is). So the user types a number (a Nat-the-value, not Nat-the-type) and we're somehow expected to know statically what the user was going to type?

The trick is to write blocks of code which are agnostic to the value of the matrix's size. Then, no matter what the input was, as long as it was a natural number, we know that that block of code will work. We can force a function to be polymorphic using higher-rank types.

withNatty :: Nat -> (forall n. Natty n -> r) -> r
withNatty Z r = r Zy
withNatty (S n) r = withNatty n (r . Sy)

withNatty n r applies the function r to the singleton representation of the natural number n. r has the Natty n available at runtime, so it can recover static knowledge of n by pattern matching the Natty, but n can't leak to the outside of the block. (You can also use existential quantification, which is covered briefly in Hasochism, to wrap a Natty and pass it around. It amounts to the same thing.)

So, for example, supposing we want to print an identity matrix of a dynamically-determined size:

main = do
    Just size <- fmap intToNat readLn
    withNatty size (print . mkIdentity)
        where mkIdentity :: Natty n -> Square n Int
              mkIdentity n = natter n one

ghci> main
4
Mat {unMat = (1 :> (0 :> (0 :> (0 :> V0)))) :> ((0 :> (1 :> (0 :> (0 :> V0)))) :> ((0 :> (0 :> (1 :> (0 :> V0)))) :> ((0 :> (0 :> (0 :> (1 :> V0)))) :> V0)))}

The same technique applies if you want to, say, build a matrix from a list of lists. This time it's a bit trickier because you have to prove to GHC that all of the lists have the same length, by measuring them.

withVec :: [a] -> (forall n. NATTY n => Vec n a -> r) -> r
withVec [] r = r V0
withVec (x:xs) r = withVec xs (r . (x :>))


-- this operation can fail because the input lists may not all be the same length
withMat :: [[a]] -> (forall w h. (NATTY w, NATTY h) => Mat w h a -> r) -> Maybe r
withMat xss r = assertEqualLengths xss (\vs -> withVec vs (r . Mat))
    where assertEqualLengths :: [[a]] -> (forall n. NATTY n => [Vec n a] -> r) -> Maybe r
          assertEqualLengths [] r = Just (r noVecs)
          assertEqualLengths xss@(xs:_) r = withLen xs (\n -> natter n $ r <$> traverse (assertLength n) xss)

          noVecs :: [Vec Z a]
          noVecs = []

          assertLength :: Natty n -> [a] -> Maybe (Vec n a)
          assertLength Zy [] = Just V0
          assertLength (Sy n) (x:xs) = fmap (x :>) (assertLength n xs)
          assertLength _ _ = Nothing

          withLen :: [a] -> (forall n. Natty n -> r) -> r
          withLen [] r = r Zy
          withLen (x:xs) r = withLen xs (r . Sy)

ghci> withMat [[1,2], [3,4]] show
Just "Mat {unMat = (1 :> (2 :> V0)) :> ((3 :> (4 :> V0)) :> V0)}"
ghci> withMat [[1,2], [3]] show  -- a ragged input list
Nothing

And if you want to work with square matrices, you have to prove to GHC that the matrix's height is equal to its width.

withEqual :: Natty n -> Natty m -> (n ~ m => r) -> Maybe r
withEqual Zy Zy r = Just r
withEqual (Sy n) (Sy m) r = withEqual n m r
withEqual _ _ _ = Nothing

square :: Natty w -> Natty h -> Mat w h a -> Maybe (Square w a)
square = withEqual