@@ -0,0 +1,954 @@
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
167+import LAoP.Utils.Internal
169+import Data.Functor.Contravariant
176+import Data.Type.Equality
178+import Control.DeepSeq
179+import Prelude hiding (id, (.))
181+-- | LAoP (Linear Algebra of Programming) Inductive Matrix definition.
182+data Matrix e cols rows where
183+ Empty :: Matrix e Void Void
184+ One :: e -> Matrix e () ()
185+ Join :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
186+ Fork :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
188+deriving instance (Show e) => Show (Matrix e cols rows)
190+-- | Type family that computes the cardinality of a given type dimension.
192+-- It can also count the cardinality of custom types that implement the
193+-- 'Generic' instance.
194+type family Count (d :: Type) :: Nat where
195+ Count (Natural n m) = (m - n) + 1
196+ Count (List a) = (^) 2 (Count a)
197+ Count (Either a b) = (+) (Count a) (Count b)
198+ Count (a, b) = (*) (Count a) (Count b)
199+ Count (a -> b) = (^) (Count b) (Count a)
201+ Count (M1 _ _ f p) = Count (f p)
202+ Count (K1 _ _ _) = 1
205+ Count ((:*:) a b p) = Count (a p) * Count (b p)
206+ Count ((:+:) a b p) = Count (a p) + Count (b p)
207+ Count d = Count (Rep d R)
209+-- | Type family that computes of a given type dimension from a given natural
211+-- Thanks to Li-Yao Xia this type family is super fast.
212+type family FromNat (n :: Nat) :: Type where
215+ FromNat n = FromNat' (Mod n 2 == 0) (FromNat (Div n 2))
217+type family FromNat' (b :: Bool) (m :: Type) :: Type where
218+ FromNat' 'True m = Either m m
219+ FromNat' 'False m = Either () (Either m m)
221+-- | Type family that normalizes the representation of a given data
223+type family Normalize (d :: Type) :: Type where
224+ Normalize (Either a b) = Either (Normalize a) (Normalize b)
225+ Normalize d = FromNat (Count d)
227+-- | Constraint type synonyms to keep the type signatures less convoluted
228+type Countable a = KnownNat (Count a)
229+type CountableDimensions a b = (Countable a, Countable b)
230+type Liftable a b = (Eq b, Constructable a b, Enumerable a)
231+type Enumerable a = (Enum a, Bounded a)
232+type ConstructNorm a = (Enum a, Enum (Normalize a))
233+type ConstructableNorm a b = (ConstructNorm a, ConstructNorm b)
234+type ConstructN a = Construct (Normalize a)
235+type Constructable a b = (Construct a, Construct b)
236+type ConstructableN a b = (Construct (Normalize a), Construct (Normalize b))
238+-- | It isn't possible to implement the 'id' function so it's
239+-- implementation is 'undefined'. However 'comp' can be and this partial
240+-- class implementation exists just to make the code more readable.
242+-- Please use 'iden' instead.
243+instance (Num e) => Category (Matrix e) where
244+ type Object (Matrix e) a = Liftable a a
248+instance NFData e => NFData (Matrix e cols rows) where
251+ rnf (Join a b) = rnf a `seq` rnf b
252+ rnf (Fork a b) = rnf a `seq` rnf b
254+instance Eq e => Eq (Matrix e cols rows) where
255+ Empty == Empty = True
256+ (One a) == (One b) = a == b
257+ (Join a b) == (Join c d) = a == c && b == d
258+ (Fork a b) == (Fork c d) = a == c && b == d
259+ x@(Fork _ _) == y@(Join _ _) = x == abideJF y
260+ x@(Join _ _) == y@(Fork _ _) = abideJF x == y
262+instance Num e => Num (Matrix e cols rows) where
264+ a + b = zipWithM (+) a b
266+ a - b = zipWithM (-) a b
268+ a * b = zipWithM (*) a b
271+ abs (One a) = One (abs a)
272+ abs (Join a b) = Join (abs a) (abs b)
273+ abs (Fork a b) = Fork (abs a) (abs b)
275+ signum Empty = Empty
276+ signum (One a) = One (signum a)
277+ signum (Join a b) = Join (signum a) (signum b)
278+ signum (Fork a b) = Fork (signum a) (signum b)
280+instance Ord e => Ord (Matrix e cols rows) where
281+ Empty <= Empty = True
282+ (One a) <= (One b) = a <= b
283+ (Join a b) <= (Join c d) = (a <= c) && (b <= d)
284+ (Fork a b) <= (Fork c d) = (a <= c) && (b <= d)
285+ x@(Fork _ _) <= y@(Join _ _) = x <= abideJF y
286+ x@(Join _ _) <= y@(Fork _ _) = abideJF x <= y
290+-- | Empty matrix constructor
291+empty :: Matrix e Void Void
294+-- | Unit matrix constructor
295+one :: e -> Matrix e () ()
298+-- | Matrix 'Join' constructor
299+join :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
304+-- | Matrix 'Join' constructor
305+(|||) :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
308+-- | Matrix 'Fork' constructor
309+fork :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
314+-- | Matrix 'Fork' constructor
315+(===) :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
320+-- | Constructs a column vector matrix
321+colL :: (FromLists e () rows) => [e] -> Matrix e () rows
322+colL = fromLists . map (: [])
324+-- | Constructs a row vector matrix
325+rowL :: (FromLists e cols ()) => [e] -> Matrix e cols ()
326+rowL = fromLists . (: [])
330+-- | Converts a matrix to a list of lists of elements.
331+toLists :: Matrix e cols rows -> [[e]]
333+toLists (One e) = [[e]]
334+toLists (Fork l r) = toLists l ++ toLists r
335+toLists (Join l r) = zipWith (++) (toLists l) (toLists r)
337+-- | Converts a matrix to a list of elements.
338+toList :: Matrix e cols rows -> [e]
339+toList = concat . toLists
343+-- | The zero matrix. A matrix wholly filled with zeros.
344+zeros :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols rows
345+zeros = matrixBuilder (const 0)
349+-- | The ones matrix. A matrix wholly filled with ones.
351+-- Also known as T (Top) matrix.
352+ones :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols rows
353+ones = matrixBuilder (const 1)
357+-- | The constant matrix constructor. A matrix wholly filled with a given
359+constant :: (Num e, Constructable cols rows, Enumerable cols) => e -> Matrix e cols rows
360+constant e = matrixBuilder (const e)
364+-- | The T (Top) row vector matrix.
365+bang :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols ()
368+-- | Point constant matrix
369+point :: (Num e, Liftable () a) => a -> Matrix e () a
370+point = fromF . const
375+iden :: (Num e, Liftable cols cols) => Matrix e cols cols
378+-- Matrix composition (MMM)
380+-- | Matrix composition. Equivalent to matrix-matrix multiplication.
382+-- This definition takes advantage of divide-and-conquer and fusion laws
384+comp :: (Num e) => Matrix e cr rows -> Matrix e cols cr -> Matrix e cols rows
385+comp Empty Empty = Empty
386+comp (One a) (One b) = One (a * b)
387+comp (Join a b) (Fork c d) = comp a c + comp b d -- Divide-and-conquer law
388+comp (Fork a b) c = Fork (comp a c) (comp b c) -- Fork fusion law
389+comp c (Join a b) = Join (comp c a) (comp c b) -- Join fusion law
390+{-# NOINLINE comp #-}
392+ "comp/iden1" forall m. comp m iden = m ;
393+ "comp/iden2" forall m. comp iden m = m
396+-- Scalar multiplication of matrices
399+-- | Scalar multiplication of matrices.
400+(.|) :: Num e => e -> Matrix e cols rows -> Matrix e cols rows
402+(.|) e (One a) = One (e * a)
403+(.|) e (Join a b) = Join (e .| a) (e .| b)
404+(.|) e (Fork a b) = Fork (e .| a) (e .| b)
406+-- Scalar division of matrices
409+-- | Scalar multiplication of matrices.
410+(./) :: Fractional e => Matrix e cols rows -> e -> Matrix e cols rows
412+(./) (One a) e = One (a / e)
413+(./) (Join a b) e = Join (a ./ e) (b ./ e)
414+(./) (Fork a b) e = Fork (a ./ e) (b ./ e)
418+-- | Biproduct first component projection
419+p1 :: (Num e, Eq m, Constructable m n, Enumerable m, Enumerable n) => Matrix e (Either m n) m
422+-- | Biproduct second component projection
423+p2 :: (Num e, Eq n, Constructable m n, Enumerable m, Enumerable n) => Matrix e (Either m n) n
428+-- | Biproduct first component injection
429+i1 :: (Num e, Eq m, Constructable m n, Enumerable m, Enumerable n) => Matrix e m (Either m n)
432+-- | Biproduct second component injection
433+i2 :: (Num e, Eq n, Constructable m n, Enumerable m, Enumerable n) => Matrix e n (Either m n)
438+-- | Obtain the number of rows.
440+-- NOTE: The 'KnownNat' constaint is needed in order to obtain the
441+-- dimensions in constant time.
443+-- TODO: A 'rows' function that does not need the 'KnownNat' constraint in
444+-- exchange for performance.
445+rows :: forall e cols rows. (Countable rows) => Matrix e cols rows -> Int
446+rows _ = fromInteger $ natVal (Proxy :: Proxy (Count rows))
448+-- | Obtain the number of columns.
450+-- NOTE: The 'KnownNat' constaint is needed in order to obtain the
451+-- dimensions in constant time.
453+-- TODO: A 'columns' function that does not need the 'KnownNat' constraint in
454+-- exchange for performance.
455+columns :: forall e cols rows. (Countable cols) => Matrix e cols rows -> Int
456+columns _ = fromInteger $ natVal (Proxy :: Proxy (Count cols))
458+-- Coproduct Bifunctor
462+-- | Matrix coproduct functor also known as matrix direct sum.
470+ ) => Matrix e n k -> Matrix e m j -> Matrix e (Either n m) (Either k j)
471+(-|-) a b = Join (i1 . a) (i2 . b)
473+-- Khatri Rao Product and projections
475+-- | Khatri Rao product first component projection matrix.
480+ ConstructableN (a, b) a,
481+ ConstructableNorm (a, b) a,
484+ ) => Matrix e (Normalize (a, b)) (Normalize a)
485+fstM = fromFN (fst :: (a, b) -> a)
487+-- | Khatri Rao product second component projection matrix.
492+ ConstructableN (a, b) b,
493+ ConstructableNorm (a, b) b,
496+ ) => Matrix e (Normalize (a, b)) (Normalize b)
497+sndM = fromFN (snd :: (a, b) -> b)
499+-- | Khatri Rao Matrix product also known as matrix pairing.
501+-- NOTE: That this is not a true categorical product, see for instance:
504+-- | fstM . kr a b == a
506+-- | sndM . kr a b == b
509+-- __Emphasis__ on the implication symbol.
515+ ConstructableN a (a, b),
517+ ConstructableNorm a b,
518+ ConstructNorm (a, b),
521+ ) => Matrix e cols (Normalize a) -> Matrix e cols (Normalize b) -> Matrix e cols (Normalize (a, b))
523+ let fstM' = fstM @e @a @b
524+ sndM' = sndM @e @a @b
525+ in (tr fstM' . a) * (tr sndM' . b)
527+-- Product Bifunctor (Kronecker)
531+type Kronecker e m n p q =
538+ ConstructableNorm m (m, n),
542+ ConstructableNorm p (p, q),
551+-- | Matrix product functor also known as kronecker product
555+ => Matrix e (Normalize m) (Normalize p)
556+ -> Matrix e (Normalize n) (Normalize q)
557+ -> Matrix e (Normalize (m, n)) (Normalize (p, q))
559+ let fstM' = fstM @e @m @n
560+ sndM' = sndM @e @m @n
561+ in kr @e @(Normalize (m, n)) @p @q (a . fstM') (b . sndM')
563+-- Matrix abide Join Fork
565+-- | Matrix "abiding" following the 'Join'-'Fork' exchange law.
570+-- 'Join' ('Fork' a c) ('Fork' b d) == 'Fork' ('Join' a b) ('Join' c d)
572+abideJF :: Matrix e cols rows -> Matrix e cols rows
573+abideJF (Join (Fork a c) (Fork b d)) = Fork (Join (abideJF a) (abideJF b)) (Join (abideJF c) (abideJF d)) -- Join-Fork abide law
574+abideJF Empty = Empty
575+abideJF (One e) = One e
576+abideJF (Join a b) = Join (abideJF a) (abideJF b)
577+abideJF (Fork a b) = Fork (abideJF a) (abideJF b)
579+-- Matrix abide Fork Join
581+-- | Matrix "abiding" followin the 'Fork'-'Join' abide law.
584+-- 'Fork' ('Join' a b) ('Join' c d) == 'Join' ('Fork' a c) ('Fork' b d)
586+abideFJ :: Matrix e cols rows -> Matrix e cols rows
587+abideFJ (Fork (Join a b) (Join c d)) = Join (Fork (abideFJ a) (abideFJ c)) (Fork (abideFJ b) (abideFJ d)) -- Fork-Join abide law
588+abideFJ Empty = Empty
589+abideFJ (One e) = One e
590+abideFJ (Join a b) = Join (abideFJ a) (abideFJ b)
591+abideFJ (Fork a b) = Fork (abideFJ a) (abideFJ b)
593+-- Matrix transposition
595+-- | Matrix transposition.
596+tr :: Matrix e cols rows -> Matrix e rows cols
599+tr (Join a b) = Fork (tr a) (tr b)
600+tr (Fork a b) = Join (tr a) (tr b)
602+-- Selective 'select' operator
604+-- | Selective functors 'select' operator equivalent inspired by the
605+-- ArrowMonad solution presented in the paper.
606+select :: (Num e, Liftable b b) => Matrix e cols (Either a b) -> Matrix e a b -> Matrix e cols b
607+select (Fork a b) y = y . a + b -- Divide-and-conquer law
608+select (Join (Fork a c) (Fork b d)) y = join (y . a + c) (y . b + d) -- Pattern matching + DnC law
609+select m y = join y id . m
622+ => Matrix e cols (Either a b) -> Matrix e a c -> Matrix e b c -> Matrix e cols c
623+branch x l r = f x `select` g l `select` r
625+ f m = fork (tr i1) (i1 . tr i2) . m
630+prettyAux :: Show e => [[e]] -> [[e]] -> String
632+prettyAux [[e]] m = "│ " ++ fill (show e) ++ " │\n"
635+ widest = maximum $ fmap length v
636+ fill str = replicate (widest - length str - 2) ' ' ++ str
637+prettyAux [h] m = "│ " ++ fill (unwords $ map show h) ++ " │\n"
640+ widest = maximum $ fmap length v
641+ fill str = replicate (widest - length str - 2) ' ' ++ str
642+prettyAux (h : t) l = "│ " ++ fill (unwords $ map show h) ++ " │\n" ++
646+ widest = maximum $ fmap length v
647+ fill str = replicate (widest - length str - 2) ' ' ++ str
649+-- | Matrix pretty printer
650+pretty :: (CountableDimensions cols rows, Show e) => Matrix e cols rows -> String
652+ [ "┌ ", unwords (replicate (columns m) blank), " ┐\n"
654+ [ "│ " ++ unwords (fmap (\j -> fill $ show $ getElem i j m) [1..columns m]) ++ " │" | i <- [1..rows m] ]
655+ , "└ ", unwords (replicate (columns m) blank), " ┘"
658+ strings = map show (toList m)
659+ widest = maximum $ map length strings
660+ fill str = replicate (widest - length str) ' ' ++ str
663+ | i > rows m || j > columns m || i < 1 || j < 1 = Nothing
664+ | otherwise = Just $ unsafeGet i j m (toList m)
665+ unsafeGet i j m l = l !! encode (columns m) (i,j)
666+ encode m (i,j) = (i-1)*m + j - 1
670+ "getElem: Trying to get the "
672+ ++ " element from a "
673+ ++ show (rows m) ++ "x" ++ show (columns m)
678+-- | Matrix pretty printer
679+prettyPrint :: (CountableDimensions cols rows, Show e) => Matrix e cols rows -> IO ()
680+prettyPrint = putStrLn . pretty
682+-- | Zip two matrices with a given binary function
683+zipWithM :: (e -> f -> g) -> Matrix e cols rows -> Matrix f cols rows -> Matrix g cols rows
684+zipWithM _ Empty Empty = Empty
685+zipWithM f (One a) (One b) = One (f a b)
686+zipWithM f (Join a b) (Join c d) = Join (zipWithM f a c) (zipWithM f b d)
687+zipWithM f (Fork a b) (Fork c d) = Fork (zipWithM f a c) (zipWithM f b d)
688+zipWithM f x@(Fork _ _) y@(Join _ _) = zipWithM f x (abideJF y)
689+zipWithM f x@(Join _ _) y@(Fork _ _) = zipWithM f (abideJF x) y
691+-- Relational operators functions
693+type Boolean = Natural 0 1
694+type Relation a b = Matrix Boolean a b
696+-- | Helper conversion function
697+toBool :: (Num e, Eq e) => e -> Bool
702+-- | Helper conversion function
703+fromBool :: Bool -> Natural 0 1
704+fromBool True = nat 1
705+fromBool False = nat 0
707+-- | Relational negation
708+negateM :: Relation cols rows -> Relation cols rows
709+negateM Empty = Empty
710+negateM (One (Nat p)) = One (Nat (negate p))
711+negateM (Join a b) = Join (negateM a) (negateM b)
712+negateM (Fork a b) = Fork (negateM a) (negateM b)
714+-- | Relational addition
715+orM :: Relation cols rows -> Relation cols rows -> Relation cols rows
716+orM Empty Empty = Empty
717+orM (One a) (One b) = One (fromBool (toBool a || toBool b))
718+orM (Join a b) (Join c d) = Join (orM a c) (orM b d)
719+orM (Fork a b) (Fork c d) = Fork (orM a c) (orM b d)
720+orM x@(Fork _ _) y@(Join _ _) = orM x (abideJF y)
721+orM x@(Join _ _) y@(Fork _ _) = orM (abideJF x) y
723+-- | Relational multiplication
724+andM :: Relation cols rows -> Relation cols rows -> Relation cols rows
725+andM Empty Empty = Empty
726+andM (One a) (One b) = One (fromBool (toBool a && toBool b))
727+andM (Join a b) (Join c d) = Join (andM a c) (andM b d)
728+andM (Fork a b) (Fork c d) = Fork (andM a c) (andM b d)
729+andM x@(Fork _ _) y@(Join _ _) = andM x (abideJF y)
730+andM x@(Join _ _) y@(Fork _ _) = andM (abideJF x) y
732+-- | Relational subtraction
733+subM :: Relation cols rows -> Relation cols rows -> Relation cols rows
734+subM Empty Empty = Empty
735+subM (One a) (One b) = if a - b < nat 0 then One (nat 0) else One (a - b)
736+subM (Join a b) (Join c d) = Join (subM a c) (subM b d)
737+subM (Fork a b) (Fork c d) = Fork (subM a c) (subM b d)
738+subM x@(Fork _ _) y@(Join _ _) = subM x (abideJF y)
739+subM x@(Join _ _) y@(Fork _ _) = subM (abideJF x) y
741+-- | Matrix relational composition.
742+compRel :: Relation cr rows -> Relation cols cr -> Relation cols rows
743+compRel Empty Empty = Empty
744+compRel (One a) (One b) = One (fromBool (toBool a && toBool b))
745+compRel (Join a b) (Fork c d) = orM (compRel a c) (compRel b d) -- Divide-and-conquer law
746+compRel (Fork a b) c = Fork (compRel a c) (compRel b c) -- Fork fusion law
747+compRel c (Join a b) = Join (compRel c a) (compRel c b) -- Join fusion law
749+-- | Matrix relational right division
750+divR :: Relation b c -> Relation b a -> Relation a c
751+divR Empty Empty = Empty
752+divR (One a) (One b) = One (fromBool (not (toBool b) || toBool a)) -- b implies a
753+divR (Join a b) (Join c d) = andM (divR a c) (divR b d)
754+divR (Fork a b) c = Fork (divR a c) (divR b c)
755+divR c (Fork a b) = Join (divR c a) (divR c b)
757+-- | Matrix relational left division
758+divL :: Relation c b -> Relation a b -> Relation a c
759+divL x y = tr (divR (tr y) (tr x))
761+-- | Matrix relational symmetric division
762+divS :: Relation c a -> Relation b a -> Relation c b
763+divS s r = divL r s `intersection` divR (tr r) (tr s)
767+-- | Lifts functions to relations with dimensions matching @a@ and @b@
771+ ConstructableNorm a b,
775+ => (a -> b) -> Relation (Normalize a) (Normalize b)
778+-- | Lifts a relation function to a Boolean Matrix
780+ ( ConstructableN a b,
781+ ConstructableNorm a b,
784+ => (a -> b -> Bool) -> Relation (Normalize a) (Normalize b)
785+toRelN f = matrixBuilderN (\a -> if uncurry f a then 1 else 0)
787+----------------------------- Linear map semantics -----------------------------
788+newtype Vector e a = Vector { at :: a -> e }
790+instance Contravariant (Vector e) where
791+ contramap f (Vector g) = Vector (g . f)
793+instance Num e => Num (Vector e a) where
794+ fromInteger = Vector . const . fromInteger
800+ negate = liftV1 negate
801+ signum = error "No sensible definition"
803+liftV1 :: (e -> e) -> Vector e a -> Vector e a
804+liftV1 f x = Vector (f . at x)
806+liftV2 :: (e -> e -> e) -> Vector e a -> Vector e a -> Vector e a
807+liftV2 f x y = Vector (\a -> f (at x a) (at y a))
809+-- Semantics of Matrix e a b
810+type LinearMap e a b = Vector e a -> Vector e b
812+semantics :: Num e => Matrix e a b -> LinearMap e a b
813+semantics m = case m of
815+ One e -> const (Vector (const e))
816+ Join x y -> \v -> semantics x (Left >$< v) + semantics y (Right >$< v)
817+ Fork x y -> \v -> Vector $ either (at (semantics x v)) (at (semantics y v))
819+padLeft :: Num e => Vector e b -> Vector e (Either a b)
820+padLeft v = Vector $ \case Left _ -> 0
823+padRight :: Num e => Vector e a -> Vector e (Either a b)
824+padRight v = Vector $ \case Left a -> at v a
827+dot :: (Num e, Enumerable a) => Vector e a -> Vector e a -> e
828+dot x y = sum [ at x a * at y a | a <- enumerate ]
830+--------------------------------- Construction ---------------------------------
832+class Construct a where
833+ row :: Num e => Vector e a -> Matrix e a ()
835+ linearMap :: (Construct b, Num e) => LinearMap e a b -> Matrix e a b
837+instance Construct () where
838+ row v = one (at v ())
840+ linearMap m = col (m 1)
842+instance (Constructable a b) => Construct (Either a b) where
843+ row v = row (Left >$< v) ||| row (Right >$< v)
845+ linearMap m = linearMap (m . padRight) ||| linearMap (m . padLeft)
847+col :: (Construct a, Num e) => Vector e a -> Matrix e () a
850+-- | Matrix builder function. Constructs a matrix provided with
851+-- a construction function that operates with canonical types.
852+matrixBuilder :: (Num e, Constructable a b, Enumerable a) => ((a, b) -> e) -> Matrix e a b
853+matrixBuilder f = linearMap $ \v -> Vector $ \b -> dot v $ Vector $ \a -> f (a, b)
855+-- | Lifts functions to matrices.
856+fromF :: (Num e, Liftable a b) => (a -> b) -> Matrix e a b
857+fromF f = matrixBuilder (\(a, b) -> if f a == b then 1 else 0)
859+-------------------------------- Deconstruction --------------------------------
860+enumerate :: Enumerable a => [a]
861+enumerate = [minBound .. maxBound]
863+basis :: (Enumerable a, Eq a, Num e) => [Vector e a]
864+basis = [ Vector (bool 0 1 . (==a)) | a <- enumerate ]
866+toLists' :: (Enumerable a, Enumerable b, Eq a, Num e) => Matrix e a b -> [[e]]
867+toLists' m = transpose
868+ [ [ at r i | i <- enumerate ] | c <- basis, let r = semantics m c ]
870+dump :: (Enumerable a, Enumerable b, Eq a, Num e, Show e) => Matrix e a b -> IO ()
871+dump = mapM_ print . toLists'
873+------------------------------- Normalised types -------------------------------
874+class Profunctor p where
875+ dimap :: (a -> b) -> (c -> d) -> p b c -> p a d
877+instance Profunctor (->) where
878+ dimap f g h = g . h . f
880+toNorm :: ConstructNorm a => a -> Normalize a
881+toNorm = toEnum . fromEnum
883+fromNorm :: ConstructNorm a => Normalize a -> a
884+fromNorm = toEnum . fromEnum
886+rowN :: (ConstructN a, ConstructNorm a, Num e) => Vector e a -> Matrix e (Normalize a) ()
887+rowN = row . contramap fromNorm
889+linearMapN :: (ConstructableN a b, ConstructableNorm a b, Num e)
890+ => LinearMap e a b -> Matrix e (Normalize a) (Normalize b)
891+linearMapN = linearMap . dimap (contramap toNorm) (contramap fromNorm)
893+colN :: (ConstructN a, ConstructNorm a, Num e) => Vector e a -> Matrix e () (Normalize a)
896+-- | Matrix builder function. Constructs a matrix provided with
897+-- a construction function that operates with arbitrary types.
898+matrixBuilderN :: (Num e, ConstructableN a b, ConstructableNorm a b, Enumerable a)
899+ => ((a, b) -> e) -> Matrix e (Normalize a) (Normalize b)
900+matrixBuilderN f = linearMapN $ \v -> Vector $ \b -> dot v $ Vector $ \a -> f (a, b)
902+-- | Lifts functions to matrices with dimensions matching @a@ and @b@
904+fromFN :: (Num e, Eq b, ConstructableN a b, ConstructableNorm a b, Enumerable a)
905+ => (a -> b) -> Matrix e (Normalize a) (Normalize b)
906+fromFN f = matrixBuilderN (\(a, b) -> if f a == b then 1 else 0)
908+------------------------------- Unsafe fromLists matrix construction -------------------------------
910+-- | Type class for defining the 'fromList' conversion function.
912+-- Given that it is not possible to branch on types at the term level type
913+-- classes are needed very much like an inductive definition but on types.
914+class FromLists e cols rows where
915+ -- | Build a matrix out of a list of list of elements. Throws a runtime
916+ -- error if the dimensions do not match.
917+ fromLists :: [[e]] -> Matrix e cols rows
919+instance FromLists e Void Void where
920+ fromLists [] = Empty
921+ fromLists _ = error "Wrong dimensions"
923+instance {-# OVERLAPPING #-} FromLists e () () where
924+ fromLists [[e]] = One e
925+ fromLists _ = error "Wrong dimensions"
927+instance {-# OVERLAPPING #-} (FromLists e cols ()) => FromLists e (Either () cols) () where
928+ fromLists [h : t] = Join (One h) (fromLists [t])
929+ fromLists _ = error "Wrong dimensions"
931+instance {-# OVERLAPPABLE #-} (FromLists e a (), FromLists e b (), Countable a) => FromLists e (Either a b) () where
933+ let rowsA = fromInteger (natVal (Proxy :: Proxy (Count a)))
934+ in Join (fromLists [take rowsA l]) (fromLists [drop rowsA l])
935+ fromLists _ = error "Wrong dimensions"
937+instance {-# OVERLAPPING #-} (FromLists e () rows) => FromLists e () (Either () rows) where
938+ fromLists ([h] : t) = Fork (One h) (fromLists t)
939+ fromLists _ = error "Wrong dimensions"
941+instance {-# OVERLAPPABLE #-} (FromLists e () a, FromLists e () b, Countable a) => FromLists e () (Either a b) where
942+ fromLists l@([_] : _) =
943+ let rowsA = fromInteger (natVal (Proxy :: Proxy (Count a)))
944+ in Fork (fromLists (take rowsA l)) (fromLists (drop rowsA l))
945+ fromLists _ = error "Wrong dimensions"
947+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
948+ fromLists l@(h : t) =
950+ rowsC = fromInteger (natVal (Proxy :: Proxy (Count c)))
951+ condition = all ((== lh) . length) t
952+ in if lh > 0 && condition
953+ then Fork (fromLists (take rowsC l)) (fromLists (drop rowsC l))
954+ else error "Not all rows have the same length"