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:
BoundedInt
type with a Bounded
instance that can be set at runtime when we know the size of the array, orApplicative (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?
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With