@@ -0,0 +1,949 @@
1+{-# LANGUAGE TypeApplications #-}
2+{-# LANGUAGE LambdaCase #-}
3+{-# LANGUAGE StrictData #-}
4+{-# LANGUAGE TypeFamilyDependencies #-}
5+{-# LANGUAGE ConstraintKinds #-}
6+{-# LANGUAGE AllowAmbiguousTypes #-}
7+{-# LANGUAGE DataKinds #-}
8+{-# LANGUAGE FlexibleContexts #-}
9+{-# LANGUAGE FlexibleInstances #-}
10+{-# LANGUAGE GADTs #-}
11+{-# LANGUAGE InstanceSigs #-}
12+{-# LANGUAGE MultiParamTypeClasses #-}
13+{-# LANGUAGE NoStarIsType #-}
14+{-# LANGUAGE ScopedTypeVariables #-}
15+{-# LANGUAGE StandaloneDeriving #-}
16+{-# LANGUAGE TypeOperators #-}
17+{-# LANGUAGE UndecidableInstances #-}
19+-----------------------------------------------------------------------------
21+-- Module : LAoP.Matrix.Internal
22+-- Copyright : (c) Armando Santos 2019-2020
23+-- Maintainer : armandoifsantos@gmail.com
24+-- Stability : experimental
26+-- The LAoP discipline generalises relations and functions treating them as
27+-- Boolean matrices and in turn consider these as arrows.
29+-- __LAoP__ is a library for algebraic (inductive) construction and manipulation of matrices
30+-- in Haskell. See <https://github.com/bolt12/master-thesis my Msc Thesis> for the
31+-- motivation behind the library, the underlying theory, and implementation details.
33+-- This module offers many of the combinators mentioned in the work of
34+-- Macedo (2012) and Oliveira (2012).
36+-- This is an Internal module and it is no supposed to be imported.
38+-----------------------------------------------------------------------------
40+module LAoP.Matrix.Alternative
41+ ( -- | This definition makes use of the fact that 'Void' is
42+ -- isomorphic to 0 and '()' to 1 and captures matrix
43+ -- dimensions as stacks of 'Either's.
45+ -- There exists two type families that make it easier to write
46+ -- matrix dimensions: 'FromNat' and 'Count'. This approach
47+ -- leads to a very straightforward implementation
48+ -- of LAoP combinators.
50+ -- * Type safe matrix representation
59+ -- * Auxiliary type families
64+ -- * Matrix construction and conversion
86+ -- ** Matrix Transposition
89+ -- ** Scalar multiplication/division of matrices
93+ -- ** Selective operator
97+ -- ** McCarthy's Conditional
100+ -- ** Matrix "abiding"
107+ -- * Biproduct approach
122+ -- ** Applicative matrix combinators
124+ -- | Note that given the restrictions imposed it is not possible to
125+ -- implement the standard type classes present in standard Haskell.
126+ -- *** Matrix pairing projections
130+ -- *** Matrix pairing
133+ -- * Matrix composition and lifting
135+ -- ** Arrow matrix combinators
137+ -- | Note that given the restrictions imposed it is not possible to
138+ -- implement the standard type classes present in standard Haskell.
144+ -- * Matrix printing
164+import LAoP.Utils.Internal
166+import Data.Functor.Contravariant
173+import Data.Type.Equality
175+import Control.DeepSeq
176+import Prelude hiding (id, (.))
178+-- | LAoP (Linear Algebra of Programming) Inductive Matrix definition.
179+data Matrix e cols rows where
180+ Empty :: Matrix e Void Void
181+ One :: e -> Matrix e () ()
182+ Join :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
183+ Fork :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
185+deriving instance (Show e) => Show (Matrix e cols rows)
187+-- | Type family that computes the cardinality of a given type dimension.
189+-- It can also count the cardinality of custom types that implement the
190+-- 'Generic' instance.
191+type family Count (d :: Type) :: Nat where
192+ Count (Natural n m) = (m - n) + 1
193+ Count (List a) = (^) 2 (Count a)
194+ Count (Either a b) = (+) (Count a) (Count b)
195+ Count (a, b) = (*) (Count a) (Count b)
196+ Count (a -> b) = (^) (Count b) (Count a)
198+ Count (M1 _ _ f p) = Count (f p)
199+ Count (K1 _ _ _) = 1
202+ Count ((:*:) a b p) = Count (a p) * Count (b p)
203+ Count ((:+:) a b p) = Count (a p) + Count (b p)
204+ Count d = Count (Rep d R)
206+-- | Type family that computes of a given type dimension from a given natural
208+-- Thanks to Li-Yao Xia this type family is super fast.
209+type family FromNat (n :: Nat) :: Type where
212+ FromNat n = FromNat' (Mod n 2 == 0) (FromNat (Div n 2))
214+type family FromNat' (b :: Bool) (m :: Type) :: Type where
215+ FromNat' 'True m = Either m m
216+ FromNat' 'False m = Either () (Either m m)
218+-- | Type family that normalizes the representation of a given data
220+type family Normalize (d :: Type) :: Type where
221+ Normalize (Either a b) = Either (Normalize a) (Normalize b)
222+ Normalize d = FromNat (Count d)
224+-- | Constraint type synonyms to keep the type signatures less convoluted
225+type Countable a = KnownNat (Count a)
226+type CountableN a = KnownNat (Count (Normalize a))
227+type CountableDimensions a b = (Countable a, Countable b)
228+type Enumerable a = (Enum a, Bounded a)
229+type ConstructNorm a = (Enum a, Enum (Normalize a))
230+type ConstructableNorm a b = (ConstructNorm a, ConstructNorm b)
231+type ConstructN a = Construct (Normalize a)
232+type Constructable a b = (Construct a, Construct b)
233+type ConstructableN a b = (Construct (Normalize a), Construct (Normalize b))
236+-- | It isn't possible to implement the 'id' function so it's
237+-- implementation is 'undefined'. However 'comp' can be and this partial
238+-- class implementation exists just to make the code more readable.
240+-- Please use 'iden' instead.
241+instance (Num e) => Category (Matrix e) where
245+instance NFData e => NFData (Matrix e cols rows) where
248+ rnf (Join a b) = rnf a `seq` rnf b
249+ rnf (Fork a b) = rnf a `seq` rnf b
251+instance Eq e => Eq (Matrix e cols rows) where
252+ Empty == Empty = True
253+ (One a) == (One b) = a == b
254+ (Join a b) == (Join c d) = a == c && b == d
255+ (Fork a b) == (Fork c d) = a == c && b == d
256+ x@(Fork _ _) == y@(Join _ _) = x == abideJF y
257+ x@(Join _ _) == y@(Fork _ _) = abideJF x == y
259+instance Num e => Num (Matrix e cols rows) where
261+ a + b = zipWithM (+) a b
263+ a - b = zipWithM (-) a b
265+ a * b = zipWithM (*) a b
268+ abs (One a) = One (abs a)
269+ abs (Join a b) = Join (abs a) (abs b)
270+ abs (Fork a b) = Fork (abs a) (abs b)
272+ signum Empty = Empty
273+ signum (One a) = One (signum a)
274+ signum (Join a b) = Join (signum a) (signum b)
275+ signum (Fork a b) = Fork (signum a) (signum b)
277+instance Ord e => Ord (Matrix e cols rows) where
278+ Empty <= Empty = True
279+ (One a) <= (One b) = a <= b
280+ (Join a b) <= (Join c d) = (a <= c) && (b <= d)
281+ (Fork a b) <= (Fork c d) = (a <= c) && (b <= d)
282+ x@(Fork _ _) <= y@(Join _ _) = x <= abideJF y
283+ x@(Join _ _) <= y@(Fork _ _) = abideJF x <= y
287+-- | Empty matrix constructor
288+empty :: Matrix e Void Void
291+-- | Unit matrix constructor
292+one :: e -> Matrix e () ()
295+-- | Matrix 'Join' constructor
296+join :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
301+-- | Matrix 'Join' constructor
302+(|||) :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
305+-- | Matrix 'Fork' constructor
306+fork :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
311+-- | Matrix 'Fork' constructor
312+(===) :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
317+-- | Type class for defining the 'fromList' conversion function.
319+-- Given that it is not possible to branch on types at the term level type
320+-- classes are needed very much like an inductive definition but on types.
321+class FromLists e cols rows where
322+ -- | Build a matrix out of a list of list of elements. Throws a runtime
323+ -- error if the dimensions do not match.
324+ fromLists :: [[e]] -> Matrix e cols rows
326+instance FromLists e Void Void where
327+ fromLists [] = Empty
328+ fromLists _ = error "Wrong dimensions"
330+instance {-# OVERLAPPING #-} FromLists e () () where
331+ fromLists [[e]] = One e
332+ fromLists _ = error "Wrong dimensions"
334+instance {-# OVERLAPPING #-} (FromLists e cols ()) => FromLists e (Either () cols) () where
335+ fromLists [h : t] = Join (One h) (fromLists [t])
336+ fromLists _ = error "Wrong dimensions"
338+instance {-# OVERLAPPABLE #-} (FromLists e a (), FromLists e b (), Countable a) => FromLists e (Either a b) () where
340+ let rowsA = fromInteger (natVal (Proxy :: Proxy (Count a)))
341+ in Join (fromLists [take rowsA l]) (fromLists [drop rowsA l])
342+ fromLists _ = error "Wrong dimensions"
344+instance {-# OVERLAPPING #-} (FromLists e () rows) => FromLists e () (Either () rows) where
345+ fromLists ([h] : t) = Fork (One h) (fromLists t)
346+ fromLists _ = error "Wrong dimensions"
348+instance {-# OVERLAPPABLE #-} (FromLists e () a, FromLists e () b, Countable a) => FromLists e () (Either a b) where
349+ fromLists l@([_] : _) =
350+ let rowsA = fromInteger (natVal (Proxy :: Proxy (Count a)))
351+ in Fork (fromLists (take rowsA l)) (fromLists (drop rowsA l))
352+ fromLists _ = error "Wrong dimensions"
354+instance {-# OVERLAPPABLE #-} (FromLists e (Either a b) c, FromLists e (Either a b) d, Countable c) => FromLists e (Either a b) (Either c d) where
355+ fromLists l@(h : t) =
357+ rowsC = fromInteger (natVal (Proxy :: Proxy (Count c)))
358+ condition = all ((== lh) . length) t
359+ in if lh > 0 && condition
360+ then Fork (fromLists (take rowsC l)) (fromLists (drop rowsC l))
361+ else error "Not all rows have the same length"
363+-- | Constructs a column vector matrix
364+colL :: (FromLists e () rows) => [e] -> Matrix e () rows
365+colL = fromLists . map (: [])
367+-- | Constructs a row vector matrix
368+rowL :: (FromLists e cols ()) => [e] -> Matrix e cols ()
369+rowL = fromLists . (: [])
373+-- | Converts a matrix to a list of lists of elements.
374+toLists :: Matrix e cols rows -> [[e]]
376+toLists (One e) = [[e]]
377+toLists (Fork l r) = toLists l ++ toLists r
378+toLists (Join l r) = zipWith (++) (toLists l) (toLists r)
380+-- | Converts a matrix to a list of elements.
381+toList :: Matrix e cols rows -> [e]
382+toList = concat . toLists
386+-- | The zero matrix. A matrix wholly filled with zeros.
387+zeros :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols rows
388+zeros = matrixBuilder (const 0)
392+-- | The ones matrix. A matrix wholly filled with ones.
394+-- Also known as T (Top) matrix.
395+ones :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols rows
396+ones = matrixBuilder (const 1)
400+-- | The constant matrix constructor. A matrix wholly filled with a given
402+constant :: (Num e, Constructable cols rows, Enumerable cols) => e -> Matrix e cols rows
403+constant e = matrixBuilder (const e)
407+-- | The T (Top) row vector matrix.
408+bang :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols ()
411+-- | Point constant matrix
412+point :: (Num e, Eq a, Construct a) => a -> Matrix e () a
413+point = fromF . const
418+iden :: (Num e, Eq cols, Construct cols, Enumerable cols) => Matrix e cols cols
421+-- Matrix composition (MMM)
423+-- | Matrix composition. Equivalent to matrix-matrix multiplication.
425+-- This definition takes advantage of divide-and-conquer and fusion laws
427+comp :: (Num e) => Matrix e cr rows -> Matrix e cols cr -> Matrix e cols rows
428+comp Empty Empty = Empty
429+comp (One a) (One b) = One (a * b)
430+comp (Join a b) (Fork c d) = comp a c + comp b d -- Divide-and-conquer law
431+comp (Fork a b) c = Fork (comp a c) (comp b c) -- Fork fusion law
432+comp c (Join a b) = Join (comp c a) (comp c b) -- Join fusion law
433+{-# NOINLINE comp #-}
435+ "comp/iden1" forall m. comp m iden = m ;
436+ "comp/iden2" forall m. comp iden m = m
439+-- Scalar multiplication of matrices
442+-- | Scalar multiplication of matrices.
443+(.|) :: Num e => e -> Matrix e cols rows -> Matrix e cols rows
445+(.|) e (One a) = One (e * a)
446+(.|) e (Join a b) = Join (e .| a) (e .| b)
447+(.|) e (Fork a b) = Fork (e .| a) (e .| b)
449+-- Scalar division of matrices
452+-- | Scalar multiplication of matrices.
453+(./) :: Fractional e => Matrix e cols rows -> e -> Matrix e cols rows
455+(./) (One a) e = One (a / e)
456+(./) (Join a b) e = Join (a ./ e) (b ./ e)
457+(./) (Fork a b) e = Fork (a ./ e) (b ./ e)
461+-- | Biproduct first component projection
462+p1 :: (Num e, Eq m, Constructable m n, Enumerable m, Enumerable n) => Matrix e (Either m n) m
465+-- | Biproduct second component projection
466+p2 :: (Num e, Eq n, Constructable m n, Enumerable m, Enumerable n) => Matrix e (Either m n) n
471+-- | Biproduct first component injection
472+i1 :: (Num e, Eq m, Constructable m n, Enumerable m, Enumerable n) => Matrix e m (Either m n)
475+-- | Biproduct second component injection
476+i2 :: (Num e, Eq n, Constructable m n, Enumerable m, Enumerable n) => Matrix e n (Either m n)
481+-- | Obtain the number of rows.
483+-- NOTE: The 'KnownNat' constaint is needed in order to obtain the
484+-- dimensions in constant time.
486+-- TODO: A 'rows' function that does not need the 'KnownNat' constraint in
487+-- exchange for performance.
488+rows :: forall e cols rows. (Countable rows) => Matrix e cols rows -> Int
489+rows _ = fromInteger $ natVal (Proxy :: Proxy (Count rows))
491+-- | Obtain the number of columns.
493+-- NOTE: The 'KnownNat' constaint is needed in order to obtain the
494+-- dimensions in constant time.
496+-- TODO: A 'columns' function that does not need the 'KnownNat' constraint in
497+-- exchange for performance.
498+columns :: forall e cols rows. (Countable cols) => Matrix e cols rows -> Int
499+columns _ = fromInteger $ natVal (Proxy :: Proxy (Count cols))
501+-- Coproduct Bifunctor
505+-- | Matrix coproduct functor also known as matrix direct sum.
513+ ) => Matrix e n k -> Matrix e m j -> Matrix e (Either n m) (Either k j)
514+(-|-) a b = Join (i1 . a) (i2 . b)
516+-- Khatri Rao Product and projections
518+-- | Khatri Rao product first component projection matrix.
523+ ConstructableN (a, b) a,
524+ ConstructableNorm (a, b) a,
527+ ) => Matrix e (Normalize (a, b)) (Normalize a)
528+fstM = fromFN (fst :: (a, b) -> a)
530+-- | Khatri Rao product second component projection matrix.
535+ ConstructableN (a, b) b,
536+ ConstructableNorm (a, b) b,
539+ ) => Matrix e (Normalize (a, b)) (Normalize b)
540+sndM = fromFN (snd :: (a, b) -> b)
542+-- | Khatri Rao Matrix product also known as matrix pairing.
544+-- NOTE: That this is not a true categorical product, see for instance:
547+-- | fstM . kr a b == a
549+-- | sndM . kr a b == b
552+-- __Emphasis__ on the implication symbol.
558+ ConstructableN a (a, b),
560+ ConstructableNorm a b,
561+ ConstructNorm (a, b),
564+ ) => Matrix e cols (Normalize a) -> Matrix e cols (Normalize b) -> Matrix e cols (Normalize (a, b))
566+ let fstM' = fstM @e @a @b
567+ sndM' = sndM @e @a @b
568+ in (tr fstM' . a) * (tr sndM' . b)
570+-- Product Bifunctor (Kronecker)
574+type Kronecker e m n p q =
581+ ConstructableNorm m (m, n),
585+ ConstructableNorm p (p, q),
594+-- | Matrix product functor also known as kronecker product
598+ => Matrix e (Normalize m) (Normalize p)
599+ -> Matrix e (Normalize n) (Normalize q)
600+ -> Matrix e (Normalize (m, n)) (Normalize (p, q))
602+ let fstM' = fstM @e @m @n
603+ sndM' = sndM @e @m @n
604+ in kr @e @(Normalize (m, n)) @p @q (a . fstM') (b . sndM')
606+-- Matrix abide Join Fork
608+-- | Matrix "abiding" following the 'Join'-'Fork' exchange law.
613+-- 'Join' ('Fork' a c) ('Fork' b d) == 'Fork' ('Join' a b) ('Join' c d)
615+abideJF :: Matrix e cols rows -> Matrix e cols rows
616+abideJF (Join (Fork a c) (Fork b d)) = Fork (Join (abideJF a) (abideJF b)) (Join (abideJF c) (abideJF d)) -- Join-Fork abide law
617+abideJF Empty = Empty
618+abideJF (One e) = One e
619+abideJF (Join a b) = Join (abideJF a) (abideJF b)
620+abideJF (Fork a b) = Fork (abideJF a) (abideJF b)
622+-- Matrix abide Fork Join
624+-- | Matrix "abiding" followin the 'Fork'-'Join' abide law.
627+-- 'Fork' ('Join' a b) ('Join' c d) == 'Join' ('Fork' a c) ('Fork' b d)
629+abideFJ :: Matrix e cols rows -> Matrix e cols rows
630+abideFJ (Fork (Join a b) (Join c d)) = Join (Fork (abideFJ a) (abideFJ c)) (Fork (abideFJ b) (abideFJ d)) -- Fork-Join abide law
631+abideFJ Empty = Empty
632+abideFJ (One e) = One e
633+abideFJ (Join a b) = Join (abideFJ a) (abideFJ b)
634+abideFJ (Fork a b) = Fork (abideFJ a) (abideFJ b)
636+-- Matrix transposition
638+-- | Matrix transposition.
639+tr :: Matrix e cols rows -> Matrix e rows cols
642+tr (Join a b) = Fork (tr a) (tr b)
643+tr (Fork a b) = Join (tr a) (tr b)
645+-- Selective 'select' operator
647+-- | Selective functors 'select' operator equivalent inspired by the
648+-- ArrowMonad solution presented in the paper.
649+select :: (Num e) => Matrix e cols (Either a b) -> Matrix e a b -> Matrix e cols b
650+select (Fork a b) y = y . a + b -- Divide-and-conquer law
651+select (Join (Fork a c) (Fork b d)) y = join (y . a + c) (y . b + d) -- Pattern matching + DnC law
652+select m y = join y id . m
665+ => Matrix e cols (Either a b) -> Matrix e a c -> Matrix e b c -> Matrix e cols c
666+branch x l r = f x `select` g l `select` r
668+ f m = fork (tr i1) (i1 . tr i2) . m
673+prettyAux :: Show e => [[e]] -> [[e]] -> String
675+prettyAux [[e]] m = "│ " ++ fill (show e) ++ " │\n"
678+ widest = maximum $ fmap length v
679+ fill str = replicate (widest - length str - 2) ' ' ++ str
680+prettyAux [h] m = "│ " ++ fill (unwords $ map show h) ++ " │\n"
683+ widest = maximum $ fmap length v
684+ fill str = replicate (widest - length str - 2) ' ' ++ str
685+prettyAux (h : t) l = "│ " ++ fill (unwords $ map show h) ++ " │\n" ++
689+ widest = maximum $ fmap length v
690+ fill str = replicate (widest - length str - 2) ' ' ++ str
692+-- | Matrix pretty printer
693+pretty :: (CountableDimensions cols rows, Show e) => Matrix e cols rows -> String
695+ [ "┌ ", unwords (replicate (columns m) blank), " ┐\n"
697+ [ "│ " ++ unwords (fmap (\j -> fill $ show $ getElem i j m) [1..columns m]) ++ " │" | i <- [1..rows m] ]
698+ , "└ ", unwords (replicate (columns m) blank), " ┘"
701+ strings = map show (toList m)
702+ widest = maximum $ map length strings
703+ fill str = replicate (widest - length str) ' ' ++ str
706+ | i > rows m || j > columns m || i < 1 || j < 1 = Nothing
707+ | otherwise = Just $ unsafeGet i j m (toList m)
708+ unsafeGet i j m l = l !! encode (columns m) (i,j)
709+ encode m (i,j) = (i-1)*m + j - 1
713+ "getElem: Trying to get the "
715+ ++ " element from a "
716+ ++ show (rows m) ++ "x" ++ show (columns m)
721+-- | Matrix pretty printer
722+prettyPrint :: (CountableDimensions cols rows, Show e) => Matrix e cols rows -> IO ()
723+prettyPrint = putStrLn . pretty
725+-- | Zip two matrices with a given binary function
726+zipWithM :: (e -> f -> g) -> Matrix e cols rows -> Matrix f cols rows -> Matrix g cols rows
727+zipWithM _ Empty Empty = Empty
728+zipWithM f (One a) (One b) = One (f a b)
729+zipWithM f (Join a b) (Join c d) = Join (zipWithM f a c) (zipWithM f b d)
730+zipWithM f (Fork a b) (Fork c d) = Fork (zipWithM f a c) (zipWithM f b d)
731+zipWithM f x@(Fork _ _) y@(Join _ _) = zipWithM f x (abideJF y)
732+zipWithM f x@(Join _ _) y@(Fork _ _) = zipWithM f (abideJF x) y
734+-- Relational operators functions
736+type Boolean = Natural 0 1
737+type Relation a b = Matrix Boolean a b
739+-- | Helper conversion function
740+toBool :: (Num e, Eq e) => e -> Bool
745+-- | Helper conversion function
746+fromBool :: Bool -> Natural 0 1
747+fromBool True = nat 1
748+fromBool False = nat 0
750+-- | Relational negation
751+negateM :: Relation cols rows -> Relation cols rows
752+negateM Empty = Empty
753+negateM (One (Nat p)) = One (Nat (negate p))
754+negateM (Join a b) = Join (negateM a) (negateM b)
755+negateM (Fork a b) = Fork (negateM a) (negateM b)
757+-- | Relational addition
758+orM :: Relation cols rows -> Relation cols rows -> Relation cols rows
759+orM Empty Empty = Empty
760+orM (One a) (One b) = One (fromBool (toBool a || toBool b))
761+orM (Join a b) (Join c d) = Join (orM a c) (orM b d)
762+orM (Fork a b) (Fork c d) = Fork (orM a c) (orM b d)
763+orM x@(Fork _ _) y@(Join _ _) = orM x (abideJF y)
764+orM x@(Join _ _) y@(Fork _ _) = orM (abideJF x) y
766+-- | Relational multiplication
767+andM :: Relation cols rows -> Relation cols rows -> Relation cols rows
768+andM Empty Empty = Empty
769+andM (One a) (One b) = One (fromBool (toBool a && toBool b))
770+andM (Join a b) (Join c d) = Join (andM a c) (andM b d)
771+andM (Fork a b) (Fork c d) = Fork (andM a c) (andM b d)
772+andM x@(Fork _ _) y@(Join _ _) = andM x (abideJF y)
773+andM x@(Join _ _) y@(Fork _ _) = andM (abideJF x) y
775+-- | Relational subtraction
776+subM :: Relation cols rows -> Relation cols rows -> Relation cols rows
777+subM Empty Empty = Empty
778+subM (One a) (One b) = if a - b < nat 0 then One (nat 0) else One (a - b)
779+subM (Join a b) (Join c d) = Join (subM a c) (subM b d)
780+subM (Fork a b) (Fork c d) = Fork (subM a c) (subM b d)
781+subM x@(Fork _ _) y@(Join _ _) = subM x (abideJF y)
782+subM x@(Join _ _) y@(Fork _ _) = subM (abideJF x) y
784+-- | Matrix relational composition.
785+compRel :: Relation cr rows -> Relation cols cr -> Relation cols rows
786+compRel Empty Empty = Empty
787+compRel (One a) (One b) = One (fromBool (toBool a && toBool b))
788+compRel (Join a b) (Fork c d) = orM (compRel a c) (compRel b d) -- Divide-and-conquer law
789+compRel (Fork a b) c = Fork (compRel a c) (compRel b c) -- Fork fusion law
790+compRel c (Join a b) = Join (compRel c a) (compRel c b) -- Join fusion law
792+-- | Matrix relational right division
793+divR :: Relation b c -> Relation b a -> Relation a c
794+divR Empty Empty = Empty
795+divR (One a) (One b) = One (fromBool (not (toBool b) || toBool a)) -- b implies a
796+divR (Join a b) (Join c d) = andM (divR a c) (divR b d)
797+divR (Fork a b) c = Fork (divR a c) (divR b c)
798+divR c (Fork a b) = Join (divR c a) (divR c b)
800+-- | Matrix relational left division
801+divL :: Relation c b -> Relation a b -> Relation a c
802+divL x y = tr (divR (tr y) (tr x))
804+-- | Matrix relational symmetric division
805+divS :: Relation c a -> Relation b a -> Relation c b
806+divS s r = divL r s `intersection` divR (tr r) (tr s)
810+-- | Lifts functions to relations with dimensions matching @a@ and @b@
814+ ConstructableNorm a b,
818+ => (a -> b) -> Relation (Normalize a) (Normalize b)
821+-- | Lifts a relation function to a Boolean Matrix
823+ ( ConstructableN a b,
824+ ConstructableNorm a b,
827+ => (a -> b -> Bool) -> Relation (Normalize a) (Normalize b)
828+toRelN f = matrixBuilderN (\a -> if uncurry f a then 1 else 0)
830+----------------------------- Linear map semantics -----------------------------
831+newtype Vector e a = Vector { at :: a -> e }
833+instance Contravariant (Vector e) where
834+ contramap f (Vector g) = Vector (g . f)
836+instance Num e => Num (Vector e a) where
837+ fromInteger = Vector . const . fromInteger
843+ negate = liftV1 negate
844+ signum = error "No sensible definition"
846+liftV1 :: (e -> e) -> Vector e a -> Vector e a
847+liftV1 f x = Vector (f . at x)
849+liftV2 :: (e -> e -> e) -> Vector e a -> Vector e a -> Vector e a
850+liftV2 f x y = Vector (\a -> f (at x a) (at y a))
852+-- Semantics of Matrix e a b
853+type LinearMap e a b = Vector e a -> Vector e b
855+semantics :: Num e => Matrix e a b -> LinearMap e a b
856+semantics m = case m of
858+ One e -> const (Vector (const e))
859+ Join x y -> \v -> semantics x (Left >$< v) + semantics y (Right >$< v)
860+ Fork x y -> \v -> Vector $ either (at (semantics x v)) (at (semantics y v))
862+padLeft :: Num e => Vector e b -> Vector e (Either a b)
863+padLeft v = Vector $ \case Left _ -> 0
866+padRight :: Num e => Vector e a -> Vector e (Either a b)
867+padRight v = Vector $ \case Left a -> at v a
870+dot :: (Num e, Enumerable a) => Vector e a -> Vector e a -> e
871+dot x y = sum [ at x a * at y a | a <- enumerate ]
873+--------------------------------- Construction ---------------------------------
875+class Construct a where
876+ row :: Num e => Vector e a -> Matrix e a ()
878+ linearMap :: (Construct b, Num e) => LinearMap e a b -> Matrix e a b
880+instance Construct () where
881+ row v = one (at v ())
883+ linearMap m = col (m 1)
885+instance (Constructable a b) => Construct (Either a b) where
886+ row v = row (Left >$< v) ||| row (Right >$< v)
888+ linearMap m = linearMap (m . padRight) ||| linearMap (m . padLeft)
890+col :: (Construct a, Num e) => Vector e a -> Matrix e () a
893+-- | Matrix builder function. Constructs a matrix provided with
894+-- a construction function that operates with canonical types.
895+matrixBuilder :: (Num e, Constructable a b, Enumerable a) => ((a, b) -> e) -> Matrix e a b
896+matrixBuilder f = linearMap $ \v -> Vector $ \b -> dot v $ Vector $ \a -> f (a, b)
898+-- | Lifts functions to matrices.
899+fromF :: (Num e, Eq b, Constructable a b, Enumerable a) => (a -> b) -> Matrix e a b
900+fromF f = matrixBuilder (\(a, b) -> if f a == b then 1 else 0)
902+-------------------------------- Deconstruction --------------------------------
903+enumerate :: Enumerable a => [a]
904+enumerate = [minBound .. maxBound]
906+basis :: (Enumerable a, Eq a, Num e) => [Vector e a]
907+basis = [ Vector (bool 0 1 . (==a)) | a <- enumerate ]
909+toLists' :: (Enumerable a, Enumerable b, Eq a, Num e) => Matrix e a b -> [[e]]
910+toLists' m = transpose
911+ [ [ at r i | i <- enumerate ] | c <- basis, let r = semantics m c ]
913+dump :: (Enumerable a, Enumerable b, Eq a, Num e, Show e) => Matrix e a b -> IO ()
914+dump = mapM_ print . toLists'
916+------------------------------- Normalised types -------------------------------
917+class Profunctor p where
918+ dimap :: (a -> b) -> (c -> d) -> p b c -> p a d
920+instance Profunctor (->) where
921+ dimap f g h = g . h . f
923+toNorm :: ConstructNorm a => a -> Normalize a
924+toNorm = toEnum . fromEnum
926+fromNorm :: ConstructNorm a => Normalize a -> a
927+fromNorm = toEnum . fromEnum
929+rowN :: (ConstructN a, ConstructNorm a, Num e) => Vector e a -> Matrix e (Normalize a) ()
930+rowN = row . contramap fromNorm
932+linearMapN :: (ConstructableN a b, ConstructableNorm a b, Num e)
933+ => LinearMap e a b -> Matrix e (Normalize a) (Normalize b)
934+linearMapN = linearMap . dimap (contramap toNorm) (contramap fromNorm)
936+colN :: (ConstructN a, ConstructNorm a, Num e) => Vector e a -> Matrix e () (Normalize a)
939+-- | Matrix builder function. Constructs a matrix provided with
940+-- a construction function that operates with arbitrary types.
941+matrixBuilderN :: (Num e, ConstructableN a b, ConstructableNorm a b, Enumerable a)
942+ => ((a, b) -> e) -> Matrix e (Normalize a) (Normalize b)
943+matrixBuilderN f = linearMapN $ \v -> Vector $ \b -> dot v $ Vector $ \a -> f (a, b)
945+-- | Lifts functions to matrices with dimensions matching @a@ and @b@
947+fromFN :: (Num e, Eq b, ConstructableN a b, ConstructableNorm a b, Enumerable a)
948+ => (a -> b) -> Matrix e (Normalize a) (Normalize b)
949+fromFN f = matrixBuilderN (\(a, b) -> if f a == b then 1 else 0)