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