@@ -0,0 +1,956 @@
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 #-}
18+
19+-----------------------------------------------------------------------------
20+-- |
21+-- Module : LAoP.Matrix.Internal
22+-- Copyright : (c) Armando Santos 2019-2020
23+-- Maintainer : armandoifsantos@gmail.com
24+-- Stability : experimental
25+--
26+-- The LAoP discipline generalises relations and functions treating them as
27+-- Boolean matrices and in turn consider these as arrows.
28+--
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.
32+--
33+-- This module offers many of the combinators mentioned in the work of
34+-- Macedo (2012) and Oliveira (2012).
35+--
36+-- This is an Internal module and it is no supposed to be imported.
37+--
38+-----------------------------------------------------------------------------
39+
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.
44+ --
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.
49+
50+ -- * Type safe matrix representation
51+ Matrix (..),
52+
53+ -- * Primitives
54+ empty,
55+ one,
56+ join,
57+ fork,
58+
59+ -- * Auxiliary type families
60+ FromNat,
61+ Count,
62+ Normalize,
63+
64+ -- * Matrix construction and conversion
65+ FromLists,
66+ fromLists,
67+ toLists,
68+ toList,
69+ matrixBuilder,
70+ matrixBuilderN,
71+ rowL,
72+ rowN,
73+ colL,
74+ colN,
75+ zeros,
76+ ones,
77+ bang,
78+ point,
79+ constant,
80+
81+ -- * Misc
82+ -- ** Get dimensions
83+ columns,
84+ rows,
85+
86+ -- ** Matrix Transposition
87+ tr,
88+
89+ -- ** Scalar multiplication/division of matrices
90+ (.|),
91+ (./),
92+
93+ -- ** Selective operator
94+ select,
95+ branch,
96+
97+ -- ** McCarthy's Conditional
98+ -- cond,
99+
100+ -- ** Matrix "abiding"
101+ abideJF,
102+ abideFJ,
103+
104+ -- ** Zip matrices
105+ zipWithM,
106+
107+ -- * Biproduct approach
108+ -- ** Fork
109+ (===),
110+ -- *** Projections
111+ p1,
112+ p2,
113+ -- ** Join
114+ (|||),
115+ -- *** Injections
116+ i1,
117+ i2,
118+ -- ** Bifunctors
119+ (-|-),
120+ (><),
121+
122+ -- ** Applicative matrix combinators
123+
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
127+ fstM,
128+ sndM,
129+
130+ -- *** Matrix pairing
131+ kr,
132+
133+ -- * Matrix composition and lifting
134+
135+ -- ** Arrow matrix combinators
136+
137+ -- | Note that given the restrictions imposed it is not possible to
138+ -- implement the standard type classes present in standard Haskell.
139+ iden,
140+ comp,
141+ fromF,
142+ fromFN,
143+
144+ -- * Matrix printing
145+ pretty,
146+ prettyPrint,
147+
148+ -- * Other
149+ toBool,
150+ fromBool,
151+ compRel,
152+ divR,
153+ divL,
154+ divS,
155+ fromFRelN,
156+ toRelN,
157+ negateM,
158+ orM,
159+ andM,
160+ subM,
161+
162+ -- * Semantics
163+ Vector(..),
164+ LinearMap,
165+ Construct
166+ )
167+ where
168+
169+import LAoP.Utils.Internal
170+import Data.Bool
171+import Data.Functor.Contravariant
172+import Data.Kind
173+import Data.List
174+import Data.Maybe
175+import Data.Proxy
176+import Data.Void
177+import GHC.TypeLits
178+import Data.Type.Equality
179+import GHC.Generics
180+import Control.DeepSeq
181+import Prelude hiding (id, (.))
182+
183+-- | LAoP (Linear Algebra of Programming) Inductive Matrix definition.
184+data Matrix e cols rows where
185+ Empty :: Matrix e Void Void
186+ One :: e -> Matrix e () ()
187+ Join :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
188+ Fork :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
189+
190+deriving instance (Show e) => Show (Matrix e cols rows)
191+
192+-- | Type family that computes the cardinality of a given type dimension.
193+--
194+-- It can also count the cardinality of custom types that implement the
195+-- 'Generic' instance.
196+type family Count (d :: Type) :: Nat where
197+ Count (Natural n m) = (m - n) + 1
198+ Count (List a) = (^) 2 (Count a)
199+ Count (Either a b) = (+) (Count a) (Count b)
200+ Count (a, b) = (*) (Count a) (Count b)
201+ Count (a -> b) = (^) (Count b) (Count a)
202+ -- Generics
203+ Count (M1 _ _ f p) = Count (f p)
204+ Count (K1 _ _ _) = 1
205+ Count (V1 _) = 0
206+ Count (U1 _) = 1
207+ Count ((:*:) a b p) = Count (a p) * Count (b p)
208+ Count ((:+:) a b p) = Count (a p) + Count (b p)
209+ Count d = Count (Rep d R)
210+
211+-- | Type family that computes of a given type dimension from a given natural
212+--
213+-- Thanks to Li-Yao Xia this type family is super fast.
214+type family FromNat (n :: Nat) :: Type where
215+ FromNat 0 = Void
216+ FromNat 1 = ()
217+ FromNat n = FromNat' (Mod n 2 == 0) (FromNat (Div n 2))
218+
219+type family FromNat' (b :: Bool) (m :: Type) :: Type where
220+ FromNat' 'True m = Either m m
221+ FromNat' 'False m = Either () (Either m m)
222+
223+-- | Type family that normalizes the representation of a given data
224+-- structure
225+type family Normalize (d :: Type) :: Type where
226+ Normalize (Either a b) = Either (Normalize a) (Normalize b)
227+ Normalize d = FromNat (Count d)
228+
229+-- | Constraint type synonyms to keep the type signatures less convoluted
230+type Countable a = KnownNat (Count a)
231+type CountableDimensions a b = (Countable a, Countable b)
232+type Liftable a b = (Eq b, Constructable a b, Enumerable a)
233+type Enumerable a = (Enum a, Bounded a)
234+type ConstructNorm a = (Enum a, Enum (Normalize a))
235+type ConstructableNorm a b = (ConstructNorm a, ConstructNorm b)
236+type ConstructN a = Construct (Normalize a)
237+type Constructable a b = (Construct a, Construct b)
238+type ConstructableN a b = (Construct (Normalize a), Construct (Normalize b))
239+
240+-- | It isn't possible to implement the 'id' function so it's
241+-- implementation is 'undefined'. However 'comp' can be and this partial
242+-- class implementation exists just to make the code more readable.
243+--
244+-- Please use 'iden' instead.
245+instance (Num e) => Category (Matrix e) where
246+ type Object (Matrix e) a = Liftable a a
247+ id = iden
248+ (.) = comp
249+
250+instance NFData e => NFData (Matrix e cols rows) where
251+ rnf Empty = ()
252+ rnf (One e) = rnf e
253+ rnf (Join a b) = rnf a `seq` rnf b
254+ rnf (Fork a b) = rnf a `seq` rnf b
255+
256+instance Eq e => Eq (Matrix e cols rows) where
257+ Empty == Empty = True
258+ (One a) == (One b) = a == b
259+ (Join a b) == (Join c d) = a == c && b == d
260+ (Fork a b) == (Fork c d) = a == c && b == d
261+ x@(Fork _ _) == y@(Join _ _) = x == abideJF y
262+ x@(Join _ _) == y@(Fork _ _) = abideJF x == y
263+
264+instance Num e => Num (Matrix e cols rows) where
265+
266+ a + b = zipWithM (+) a b
267+
268+ a - b = zipWithM (-) a b
269+
270+ a * b = zipWithM (*) a b
271+
272+ abs Empty = Empty
273+ abs (One a) = One (abs a)
274+ abs (Join a b) = Join (abs a) (abs b)
275+ abs (Fork a b) = Fork (abs a) (abs b)
276+
277+ signum Empty = Empty
278+ signum (One a) = One (signum a)
279+ signum (Join a b) = Join (signum a) (signum b)
280+ signum (Fork a b) = Fork (signum a) (signum b)
281+
282+instance Ord e => Ord (Matrix e cols rows) where
283+ Empty <= Empty = True
284+ (One a) <= (One b) = a <= b
285+ (Join a b) <= (Join c d) = (a <= c) && (b <= d)
286+ (Fork a b) <= (Fork c d) = (a <= c) && (b <= d)
287+ x@(Fork _ _) <= y@(Join _ _) = x <= abideJF y
288+ x@(Join _ _) <= y@(Fork _ _) = abideJF x <= y
289+
290+-- Primitives
291+
292+-- | Empty matrix constructor
293+empty :: Matrix e Void Void
294+empty = Empty
295+
296+-- | Unit matrix constructor
297+one :: e -> Matrix e () ()
298+one = One
299+
300+-- | Matrix 'Join' constructor
301+join :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
302+join = Join
303+
304+infixl 3 |||
305+
306+-- | Matrix 'Join' constructor
307+(|||) :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
308+(|||) = Join
309+
310+-- | Matrix 'Fork' constructor
311+fork :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
312+fork = Fork
313+
314+infixl 2 ===
315+
316+-- | Matrix 'Fork' constructor
317+(===) :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
318+(===) = Fork
319+
320+-- Construction
321+
322+-- | Constructs a column vector matrix
323+colL :: (FromLists e () rows) => [e] -> Matrix e () rows
324+colL = fromLists . map (: [])
325+
326+-- | Constructs a row vector matrix
327+rowL :: (FromLists e cols ()) => [e] -> Matrix e cols ()
328+rowL = fromLists . (: [])
329+
330+-- Conversion
331+
332+-- | Converts a matrix to a list of lists of elements.
333+toLists :: Matrix e cols rows -> [[e]]
334+toLists Empty = []
335+toLists (One e) = [[e]]
336+toLists (Fork l r) = toLists l ++ toLists r
337+toLists (Join l r) = zipWith (++) (toLists l) (toLists r)
338+
339+-- | Converts a matrix to a list of elements.
340+toList :: Matrix e cols rows -> [e]
341+toList = concat . toLists
342+
343+-- Zeros Matrix
344+
345+-- | The zero matrix. A matrix wholly filled with zeros.
346+zeros :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols rows
347+zeros = matrixBuilder (const 0)
348+
349+-- Ones Matrix
350+
351+-- | The ones matrix. A matrix wholly filled with ones.
352+--
353+-- Also known as T (Top) matrix.
354+ones :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols rows
355+ones = matrixBuilder (const 1)
356+
357+-- Const Matrix
358+
359+-- | The constant matrix constructor. A matrix wholly filled with a given
360+-- value.
361+constant :: (Num e, Constructable cols rows, Enumerable cols) => e -> Matrix e cols rows
362+constant e = matrixBuilder (const e)
363+
364+-- Bang Matrix
365+
366+-- | The T (Top) row vector matrix.
367+bang :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols ()
368+bang = ones
369+
370+-- | Point constant matrix
371+point :: (Num e, Liftable () a) => a -> Matrix e () a
372+point = fromF . const
373+
374+-- iden Matrix
375+
376+-- | iden matrix.
377+iden :: (Num e, Liftable cols cols) => Matrix e cols cols
378+iden = fromF id
379+
380+-- Matrix composition (MMM)
381+
382+-- | Matrix composition. Equivalent to matrix-matrix multiplication.
383+--
384+-- This definition takes advantage of divide-and-conquer and fusion laws
385+-- from LAoP.
386+comp :: (Num e) => Matrix e cr rows -> Matrix e cols cr -> Matrix e cols rows
387+comp Empty Empty = Empty
388+comp (One a) (One b) = One (a * b)
389+comp (Join a b) (Fork c d) = comp a c + comp b d -- Divide-and-conquer law
390+comp (Fork a b) c = Fork (comp a c) (comp b c) -- Fork fusion law
391+comp c (Join a b) = Join (comp c a) (comp c b) -- Join fusion law
392+{-# NOINLINE comp #-}
393+{-# RULES
394+ "comp/iden1" forall m. comp m iden = m ;
395+ "comp/iden2" forall m. comp iden m = m
396+#-}
397+
398+-- Scalar multiplication of matrices
399+
400+infixl 7 .|
401+-- | Scalar multiplication of matrices.
402+(.|) :: Num e => e -> Matrix e cols rows -> Matrix e cols rows
403+(.|) _ Empty = Empty
404+(.|) e (One a) = One (e * a)
405+(.|) e (Join a b) = Join (e .| a) (e .| b)
406+(.|) e (Fork a b) = Fork (e .| a) (e .| b)
407+
408+-- Scalar division of matrices
409+
410+infixl 7 ./
411+-- | Scalar multiplication of matrices.
412+(./) :: Fractional e => Matrix e cols rows -> e -> Matrix e cols rows
413+(./) Empty _ = Empty
414+(./) (One a) e = One (a / e)
415+(./) (Join a b) e = Join (a ./ e) (b ./ e)
416+(./) (Fork a b) e = Fork (a ./ e) (b ./ e)
417+
418+-- Projections
419+
420+-- | Biproduct first component projection
421+p1 :: (Num e, Eq m, Constructable m n, Enumerable m, Enumerable n) => Matrix e (Either m n) m
422+p1 = join iden zeros
423+
424+-- | Biproduct second component projection
425+p2 :: (Num e, Eq n, Constructable m n, Enumerable m, Enumerable n) => Matrix e (Either m n) n
426+p2 = join zeros iden
427+
428+-- Injections
429+
430+-- | Biproduct first component injection
431+i1 :: (Num e, Eq m, Constructable m n, Enumerable m, Enumerable n) => Matrix e m (Either m n)
432+i1 = tr p1
433+
434+-- | Biproduct second component injection
435+i2 :: (Num e, Eq n, Constructable m n, Enumerable m, Enumerable n) => Matrix e n (Either m n)
436+i2 = tr p2
437+
438+-- Dimensions
439+
440+-- | Obtain the number of rows.
441+--
442+-- NOTE: The 'KnownNat' constaint is needed in order to obtain the
443+-- dimensions in constant time.
444+--
445+-- TODO: A 'rows' function that does not need the 'KnownNat' constraint in
446+-- exchange for performance.
447+rows :: forall e cols rows. (Countable rows) => Matrix e cols rows -> Int
448+rows _ = fromInteger $ natVal (Proxy :: Proxy (Count rows))
449+
450+-- | Obtain the number of columns.
451+--
452+-- NOTE: The 'KnownNat' constaint is needed in order to obtain the
453+-- dimensions in constant time.
454+--
455+-- TODO: A 'columns' function that does not need the 'KnownNat' constraint in
456+-- exchange for performance.
457+columns :: forall e cols rows. (Countable cols) => Matrix e cols rows -> Int
458+columns _ = fromInteger $ natVal (Proxy :: Proxy (Count cols))
459+
460+-- Coproduct Bifunctor
461+
462+infixl 5 -|-
463+
464+-- | Matrix coproduct functor also known as matrix direct sum.
465+(-|-) ::
466+ ( Num e,
467+ Eq k,
468+ Eq j,
469+ Constructable k j,
470+ Enumerable k,
471+ Enumerable j
472+ ) => Matrix e n k -> Matrix e m j -> Matrix e (Either n m) (Either k j)
473+(-|-) a b = Join (i1 . a) (i2 . b)
474+
475+-- Khatri Rao Product and projections
476+
477+-- | Khatri Rao product first component projection matrix.
478+fstM ::
479+ forall e a b .
480+ ( Num e,
481+ Eq a,
482+ ConstructableN (a, b) a,
483+ ConstructableNorm (a, b) a,
484+ Enumerable a,
485+ Enumerable b
486+ ) => Matrix e (Normalize (a, b)) (Normalize a)
487+fstM = fromFN (fst :: (a, b) -> a)
488+
489+-- | Khatri Rao product second component projection matrix.
490+sndM ::
491+ forall e a b .
492+ ( Num e,
493+ Eq b,
494+ ConstructableN (a, b) b,
495+ ConstructableNorm (a, b) b,
496+ Enumerable a,
497+ Enumerable b
498+ ) => Matrix e (Normalize (a, b)) (Normalize b)
499+sndM = fromFN (snd :: (a, b) -> b)
500+
501+-- | Khatri Rao Matrix product also known as matrix pairing.
502+--
503+-- NOTE: That this is not a true categorical product, see for instance:
504+--
505+-- @
506+-- | fstM . kr a b == a
507+-- kr a b ==> |
508+-- | sndM . kr a b == b
509+-- @
510+--
511+-- __Emphasis__ on the implication symbol.
512+kr ::
513+ forall e cols a b .
514+ ( Num e,
515+ Eq a,
516+ Eq b,
517+ ConstructableN a (a, b),
518+ ConstructN b,
519+ ConstructableNorm a b,
520+ ConstructNorm (a, b),
521+ Enumerable a,
522+ Enumerable b
523+ ) => Matrix e cols (Normalize a) -> Matrix e cols (Normalize b) -> Matrix e cols (Normalize (a, b))
524+kr a b =
525+ let fstM' = fstM @e @a @b
526+ sndM' = sndM @e @a @b
527+ in (tr fstM' . a) * (tr sndM' . b)
528+
529+-- Product Bifunctor (Kronecker)
530+
531+infixl 4 ><
532+
533+type Kronecker e m n p q =
534+ ( Num e,
535+ Eq m,
536+ Eq n,
537+ Eq p,
538+ Eq q,
539+ ConstructableN m n,
540+ ConstructableNorm m (m, n),
541+ ConstructNorm n,
542+ ConstructN (m, n),
543+ ConstructableN p q,
544+ ConstructableNorm p (p, q),
545+ ConstructNorm q,
546+ ConstructN (p, q),
547+ Enumerable m,
548+ Enumerable n,
549+ Enumerable p,
550+ Enumerable q
551+ )
552+
553+-- | Matrix product functor also known as kronecker product
554+(><) ::
555+ forall e m p n q .
556+ Kronecker e m n p q
557+ => Matrix e (Normalize m) (Normalize p)
558+ -> Matrix e (Normalize n) (Normalize q)
559+ -> Matrix e (Normalize (m, n)) (Normalize (p, q))
560+(><) a b =
561+ let fstM' = fstM @e @m @n
562+ sndM' = sndM @e @m @n
563+ in kr @e @(Normalize (m, n)) @p @q (a . fstM') (b . sndM')
564+
565+-- Matrix abide Join Fork
566+
567+-- | Matrix "abiding" following the 'Join'-'Fork' exchange law.
568+--
569+-- Law:
570+--
571+-- @
572+-- 'Join' ('Fork' a c) ('Fork' b d) == 'Fork' ('Join' a b) ('Join' c d)
573+-- @
574+abideJF :: Matrix e cols rows -> Matrix e cols rows
575+abideJF (Join (Fork a c) (Fork b d)) = Fork (Join (abideJF a) (abideJF b)) (Join (abideJF c) (abideJF d)) -- Join-Fork abide law
576+abideJF Empty = Empty
577+abideJF (One e) = One e
578+abideJF (Join a b) = Join (abideJF a) (abideJF b)
579+abideJF (Fork a b) = Fork (abideJF a) (abideJF b)
580+
581+-- Matrix abide Fork Join
582+
583+-- | Matrix "abiding" followin the 'Fork'-'Join' abide law.
584+--
585+-- @
586+-- 'Fork' ('Join' a b) ('Join' c d) == 'Join' ('Fork' a c) ('Fork' b d)
587+-- @
588+abideFJ :: Matrix e cols rows -> Matrix e cols rows
589+abideFJ (Fork (Join a b) (Join c d)) = Join (Fork (abideFJ a) (abideFJ c)) (Fork (abideFJ b) (abideFJ d)) -- Fork-Join abide law
590+abideFJ Empty = Empty
591+abideFJ (One e) = One e
592+abideFJ (Join a b) = Join (abideFJ a) (abideFJ b)
593+abideFJ (Fork a b) = Fork (abideFJ a) (abideFJ b)
594+
595+-- Matrix transposition
596+
597+-- | Matrix transposition.
598+tr :: Matrix e cols rows -> Matrix e rows cols
599+tr Empty = Empty
600+tr (One e) = One e
601+tr (Join a b) = Fork (tr a) (tr b)
602+tr (Fork a b) = Join (tr a) (tr b)
603+
604+-- Selective 'select' operator
605+
606+-- | Selective functors 'select' operator equivalent inspired by the
607+-- ArrowMonad solution presented in the paper.
608+select :: (Num e, Liftable b b) => Matrix e cols (Either a b) -> Matrix e a b -> Matrix e cols b
609+select (Fork a b) y = y . a + b -- Divide-and-conquer law
610+select (Join (Fork a c) (Fork b d)) y = join (y . a + c) (y . b + d) -- Pattern matching + DnC law
611+select m y = join y id . m
612+
613+branch ::
614+ ( Num e,
615+ Eq a,
616+ Eq b,
617+ Eq c,
618+ Constructable a b,
619+ Construct c,
620+ Enumerable a,
621+ Enumerable b,
622+ Enumerable c
623+ )
624+ => Matrix e cols (Either a b) -> Matrix e a c -> Matrix e b c -> Matrix e cols c
625+branch x l r = f x `select` g l `select` r
626+ where
627+ f m = fork (tr i1) (i1 . tr i2) . m
628+ g m = i2 . m
629+
630+-- Pretty print
631+
632+prettyAux :: Show e => [[e]] -> [[e]] -> String
633+prettyAux [] _ = ""
634+prettyAux [[e]] m = "│ " ++ fill (show e) ++ " │\n"
635+ where
636+ v = fmap show m
637+ widest = maximum $ fmap length v
638+ fill str = replicate (widest - length str - 2) ' ' ++ str
639+prettyAux [h] m = "│ " ++ fill (unwords $ map show h) ++ " │\n"
640+ where
641+ v = fmap show m
642+ widest = maximum $ fmap length v
643+ fill str = replicate (widest - length str - 2) ' ' ++ str
644+prettyAux (h : t) l = "│ " ++ fill (unwords $ map show h) ++ " │\n" ++
645+ prettyAux t l
646+ where
647+ v = fmap show l
648+ widest = maximum $ fmap length v
649+ fill str = replicate (widest - length str - 2) ' ' ++ str
650+
651+-- | Matrix pretty printer
652+pretty :: (CountableDimensions cols rows, Show e) => Matrix e cols rows -> String
653+pretty m = concat
654+ [ "┌ ", unwords (replicate (columns m) blank), " ┐\n"
655+ , unlines
656+ [ "│ " ++ unwords (fmap (\j -> fill $ show $ getElem i j m) [1..columns m]) ++ " │" | i <- [1..rows m] ]
657+ , "└ ", unwords (replicate (columns m) blank), " ┘"
658+ ]
659+ where
660+ strings = map show (toList m)
661+ widest = maximum $ map length strings
662+ fill str = replicate (widest - length str) ' ' ++ str
663+ blank = fill ""
664+ safeGet i j m
665+ | i > rows m || j > columns m || i < 1 || j < 1 = Nothing
666+ | otherwise = Just $ unsafeGet i j m (toList m)
667+ unsafeGet i j m l = l !! encode (columns m) (i,j)
668+ encode m (i,j) = (i-1)*m + j - 1
669+ getElem i j m =
670+ fromMaybe
671+ (error $
672+ "getElem: Trying to get the "
673+ ++ show (i, j)
674+ ++ " element from a "
675+ ++ show (rows m) ++ "x" ++ show (columns m)
676+ ++ " matrix."
677+ )
678+ (safeGet i j m)
679+
680+-- | Matrix pretty printer
681+prettyPrint :: (CountableDimensions cols rows, Show e) => Matrix e cols rows -> IO ()
682+prettyPrint = putStrLn . pretty
683+
684+-- | Zip two matrices with a given binary function
685+zipWithM :: (e -> f -> g) -> Matrix e cols rows -> Matrix f cols rows -> Matrix g cols rows
686+zipWithM _ Empty Empty = Empty
687+zipWithM f (One a) (One b) = One (f a b)
688+zipWithM f (Join a b) (Join c d) = Join (zipWithM f a c) (zipWithM f b d)
689+zipWithM f (Fork a b) (Fork c d) = Fork (zipWithM f a c) (zipWithM f b d)
690+zipWithM f x@(Fork _ _) y@(Join _ _) = zipWithM f x (abideJF y)
691+zipWithM f x@(Join _ _) y@(Fork _ _) = zipWithM f (abideJF x) y
692+
693+-- Relational operators functions
694+
695+type Boolean = Natural 0 1
696+type Relation a b = Matrix Boolean a b
697+
698+-- | Helper conversion function
699+toBool :: (Num e, Eq e) => e -> Bool
700+toBool n
701+ | n == 0 = False
702+ | n == 1 = True
703+
704+-- | Helper conversion function
705+fromBool :: Bool -> Natural 0 1
706+fromBool True = nat 1
707+fromBool False = nat 0
708+
709+-- | Relational negation
710+negateM :: Relation cols rows -> Relation cols rows
711+negateM Empty = Empty
712+negateM (One (Nat p)) = One (Nat (negate p))
713+negateM (Join a b) = Join (negateM a) (negateM b)
714+negateM (Fork a b) = Fork (negateM a) (negateM b)
715+
716+-- | Relational addition
717+orM :: Relation cols rows -> Relation cols rows -> Relation cols rows
718+orM Empty Empty = Empty
719+orM (One a) (One b) = One (fromBool (toBool a || toBool b))
720+orM (Join a b) (Join c d) = Join (orM a c) (orM b d)
721+orM (Fork a b) (Fork c d) = Fork (orM a c) (orM b d)
722+orM x@(Fork _ _) y@(Join _ _) = orM x (abideJF y)
723+orM x@(Join _ _) y@(Fork _ _) = orM (abideJF x) y
724+
725+-- | Relational multiplication
726+andM :: Relation cols rows -> Relation cols rows -> Relation cols rows
727+andM Empty Empty = Empty
728+andM (One a) (One b) = One (fromBool (toBool a && toBool b))
729+andM (Join a b) (Join c d) = Join (andM a c) (andM b d)
730+andM (Fork a b) (Fork c d) = Fork (andM a c) (andM b d)
731+andM x@(Fork _ _) y@(Join _ _) = andM x (abideJF y)
732+andM x@(Join _ _) y@(Fork _ _) = andM (abideJF x) y
733+
734+-- | Relational subtraction
735+subM :: Relation cols rows -> Relation cols rows -> Relation cols rows
736+subM Empty Empty = Empty
737+subM (One a) (One b) = if a - b < nat 0 then One (nat 0) else One (a - b)
738+subM (Join a b) (Join c d) = Join (subM a c) (subM b d)
739+subM (Fork a b) (Fork c d) = Fork (subM a c) (subM b d)
740+subM x@(Fork _ _) y@(Join _ _) = subM x (abideJF y)
741+subM x@(Join _ _) y@(Fork _ _) = subM (abideJF x) y
742+
743+-- | Matrix relational composition.
744+compRel :: Relation cr rows -> Relation cols cr -> Relation cols rows
745+compRel Empty Empty = Empty
746+compRel (One a) (One b) = One (fromBool (toBool a && toBool b))
747+compRel (Join a b) (Fork c d) = orM (compRel a c) (compRel b d) -- Divide-and-conquer law
748+compRel (Fork a b) c = Fork (compRel a c) (compRel b c) -- Fork fusion law
749+compRel c (Join a b) = Join (compRel c a) (compRel c b) -- Join fusion law
750+
751+-- | Matrix relational right division
752+divR :: Relation b c -> Relation b a -> Relation a c
753+divR Empty Empty = Empty
754+divR (One a) (One b) = One (fromBool (not (toBool b) || toBool a)) -- b implies a
755+divR (Join a b) (Join c d) = andM (divR a c) (divR b d)
756+divR (Fork a b) c = Fork (divR a c) (divR b c)
757+divR c (Fork a b) = Join (divR c a) (divR c b)
758+
759+-- | Matrix relational left division
760+divL :: Relation c b -> Relation a b -> Relation a c
761+divL x y = tr (divR (tr y) (tr x))
762+
763+-- | Matrix relational symmetric division
764+divS :: Relation c a -> Relation b a -> Relation c b
765+divS s r = divL r s `intersection` divR (tr r) (tr s)
766+ where
767+ intersection = andM
768+
769+-- | Lifts functions to relations with dimensions matching @a@ and @b@
770+-- cardinality's.
771+fromFRelN ::
772+ ( Eq b,
773+ ConstructableNorm a b,
774+ ConstructableN a b,
775+ Enumerable a
776+ )
777+ => (a -> b) -> Relation (Normalize a) (Normalize b)
778+fromFRelN = fromFN
779+
780+-- | Lifts a relation function to a Boolean Matrix
781+toRelN ::
782+ ( ConstructableN a b,
783+ ConstructableNorm a b,
784+ Enumerable a
785+ )
786+ => (a -> b -> Bool) -> Relation (Normalize a) (Normalize b)
787+toRelN f = matrixBuilderN (\a -> if uncurry f a then 1 else 0)
788+
789+----------------------------- Linear map semantics -----------------------------
790+newtype Vector e a = Vector { at :: a -> e }
791+
792+instance Contravariant (Vector e) where
793+ contramap f (Vector g) = Vector (g . f)
794+
795+instance Num e => Num (Vector e a) where
796+ fromInteger = Vector . const . fromInteger
797+
798+ (+) = liftV2 (+)
799+ (-) = liftV2 (-)
800+ (*) = liftV2 (*)
801+ abs = liftV1 abs
802+ negate = liftV1 negate
803+ signum = error "No sensible definition"
804+
805+liftV1 :: (e -> e) -> Vector e a -> Vector e a
806+liftV1 f x = Vector (f . at x)
807+
808+liftV2 :: (e -> e -> e) -> Vector e a -> Vector e a -> Vector e a
809+liftV2 f x y = Vector (\a -> f (at x a) (at y a))
810+
811+-- Semantics of Matrix e a b
812+type LinearMap e a b = Vector e a -> Vector e b
813+
814+semantics :: Num e => Matrix e a b -> LinearMap e a b
815+semantics m = case m of
816+ Empty -> id
817+ One e -> const (Vector (const e))
818+ Join x y -> \v -> semantics x (Left >$< v) + semantics y (Right >$< v)
819+ Fork x y -> \v -> Vector $ either (at (semantics x v)) (at (semantics y v))
820+
821+padLeft :: Num e => Vector e b -> Vector e (Either a b)
822+padLeft v = Vector $ \case Left _ -> 0
823+ Right b -> at v b
824+
825+padRight :: Num e => Vector e a -> Vector e (Either a b)
826+padRight v = Vector $ \case Left a -> at v a
827+ Right _ -> 0
828+
829+dot :: (Num e, Enumerable a) => Vector e a -> Vector e a -> e
830+dot x y = sum [ at x a * at y a | a <- enumerate ]
831+
832+--------------------------------- Construction ---------------------------------
833+
834+class Construct a where
835+ row :: Num e => Vector e a -> Matrix e a ()
836+
837+ linearMap :: (Construct b, Num e) => LinearMap e a b -> Matrix e a b
838+
839+instance Construct () where
840+ row v = one (at v ())
841+
842+ linearMap m = col (m 1)
843+
844+instance (Constructable a b) => Construct (Either a b) where
845+ row v = row (Left >$< v) ||| row (Right >$< v)
846+
847+ linearMap m = linearMap (m . padRight) ||| linearMap (m . padLeft)
848+
849+col :: (Construct a, Num e) => Vector e a -> Matrix e () a
850+col = tr . row
851+
852+-- | Matrix builder function. Constructs a matrix provided with
853+-- a construction function that operates with canonical types.
854+matrixBuilder :: (Num e, Constructable a b, Enumerable a) => ((a, b) -> e) -> Matrix e a b
855+matrixBuilder f = linearMap $ \v -> Vector $ \b -> dot v $ Vector $ \a -> f (a, b)
856+
857+-- | Lifts functions to matrices.
858+fromF :: (Num e, Liftable a b) => (a -> b) -> Matrix e a b
859+fromF f = matrixBuilder (\(a, b) -> if f a == b then 1 else 0)
860+
861+-------------------------------- Deconstruction --------------------------------
862+enumerate :: Enumerable a => [a]
863+enumerate = [minBound .. maxBound]
864+
865+basis :: (Enumerable a, Eq a, Num e) => [Vector e a]
866+basis = [ Vector (bool 0 1 . (==a)) | a <- enumerate ]
867+
868+toLists' :: (Enumerable a, Enumerable b, Eq a, Num e) => Matrix e a b -> [[e]]
869+toLists' m = transpose
870+ [ [ at r i | i <- enumerate ] | c <- basis, let r = semantics m c ]
871+
872+dump :: (Enumerable a, Enumerable b, Eq a, Num e, Show e) => Matrix e a b -> IO ()
873+dump = mapM_ print . toLists'
874+
875+------------------------------- Normalised types -------------------------------
876+class Profunctor p where
877+ dimap :: (a -> b) -> (c -> d) -> p b c -> p a d
878+
879+instance Profunctor (->) where
880+ dimap f g h = g . h . f
881+
882+toNorm :: ConstructNorm a => a -> Normalize a
883+toNorm = toEnum . fromEnum
884+
885+fromNorm :: ConstructNorm a => Normalize a -> a
886+fromNorm = toEnum . fromEnum
887+
888+rowN :: (ConstructN a, ConstructNorm a, Num e) => Vector e a -> Matrix e (Normalize a) ()
889+rowN = row . contramap fromNorm
890+
891+linearMapN :: (ConstructableN a b, ConstructableNorm a b, Num e)
892+ => LinearMap e a b -> Matrix e (Normalize a) (Normalize b)
893+linearMapN = linearMap . dimap (contramap toNorm) (contramap fromNorm)
894+
895+colN :: (ConstructN a, ConstructNorm a, Num e) => Vector e a -> Matrix e () (Normalize a)
896+colN = tr . rowN
897+
898+-- | Matrix builder function. Constructs a matrix provided with
899+-- a construction function that operates with arbitrary types.
900+matrixBuilderN :: (Num e, ConstructableN a b, ConstructableNorm a b, Enumerable a)
901+ => ((a, b) -> e) -> Matrix e (Normalize a) (Normalize b)
902+matrixBuilderN f = linearMapN $ \v -> Vector $ \b -> dot v $ Vector $ \a -> f (a, b)
903+
904+-- | Lifts functions to matrices with dimensions matching @a@ and @b@
905+-- cardinality's.
906+fromFN :: (Num e, Eq b, ConstructableN a b, ConstructableNorm a b, Enumerable a)
907+ => (a -> b) -> Matrix e (Normalize a) (Normalize b)
908+fromFN f = matrixBuilderN (\(a, b) -> if f a == b then 1 else 0)
909+
910+------------------------------- Unsafe fromLists matrix construction -------------------------------
911+
912+-- | Type class for defining the 'fromList' conversion function.
913+--
914+-- Given that it is not possible to branch on types at the term level type
915+-- classes are needed very much like an inductive definition but on types.
916+class FromLists e cols rows where
917+ -- | Build a matrix out of a list of list of elements. Throws a runtime
918+ -- error if the dimensions do not match.
919+ fromLists :: [[e]] -> Matrix e cols rows
920+
921+instance FromLists e Void Void where
922+ fromLists [] = Empty
923+ fromLists _ = error "Wrong dimensions"
924+
925+instance {-# OVERLAPPING #-} FromLists e () () where
926+ fromLists [[e]] = One e
927+ fromLists _ = error "Wrong dimensions"
928+
929+instance {-# OVERLAPPING #-} (FromLists e cols ()) => FromLists e (Either () cols) () where
930+ fromLists [h : t] = Join (One h) (fromLists [t])
931+ fromLists _ = error "Wrong dimensions"
932+
933+instance {-# OVERLAPPABLE #-} (FromLists e a (), FromLists e b (), Countable a) => FromLists e (Either a b) () where
934+ fromLists [l] =
935+ let rowsA = fromInteger (natVal (Proxy :: Proxy (Count a)))
936+ in Join (fromLists [take rowsA l]) (fromLists [drop rowsA l])
937+ fromLists _ = error "Wrong dimensions"
938+
939+instance {-# OVERLAPPING #-} (FromLists e () rows) => FromLists e () (Either () rows) where
940+ fromLists ([h] : t) = Fork (One h) (fromLists t)
941+ fromLists _ = error "Wrong dimensions"
942+
943+instance {-# OVERLAPPABLE #-} (FromLists e () a, FromLists e () b, Countable a) => FromLists e () (Either a b) where
944+ fromLists l@([_] : _) =
945+ let rowsA = fromInteger (natVal (Proxy :: Proxy (Count a)))
946+ in Fork (fromLists (take rowsA l)) (fromLists (drop rowsA l))
947+ fromLists _ = error "Wrong dimensions"
948+
949+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
950+ fromLists l@(h : t) =
951+ let lh = length h
952+ rowsC = fromInteger (natVal (Proxy :: Proxy (Count c)))
953+ condition = all ((== lh) . length) t
954+ in if lh > 0 && condition
955+ then Fork (fromLists (take rowsC l)) (fromLists (drop rowsC l))
956+ else error "Not all rows have the same length"