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
Applicative
; invent new function names that take your extra size parameters. (Whether you stuff those functions in a class or not depends on your needs.) But depending on what library operations are already available for things that supportApplicative
this may be painful... you'll know more than we will about whether that's a non-starter or not. What do you think?Applicative
,Semiring
, andStarSemiring
is certainly an option, though I think a moderately painful and inelegant one. I'll keep this as a backup option depending on the other answers.reflection
.) They're very different requirements!data Matrix :: (Nat,Nat) -> * -> * where Mat :: {unMat :: Vec h (Vec w a)} -> Matrix '(w,h) a
: What does theApplicative
instance look like, for example? If I want to define a type synonym liketype AdjMat = Matrix Int (Tropical Double)
from O'Connor, can I? What's the equivalenttoMatrix :: NATTY n => [[a]] -> Matrix (Pair n n) a
? Thanks!