@@ -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 #-}
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+ where
163+
164+import LAoP.Utils.Internal
165+import Data.Bool
166+import Data.Functor.Contravariant
167+import Data.Kind
168+import Data.List
169+import Data.Maybe
170+import Data.Proxy
171+import Data.Void
172+import GHC.TypeLits
173+import Data.Type.Equality
174+import GHC.Generics
175+import Control.DeepSeq
176+import Prelude hiding (id, (.))
177+
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)
184+
185+deriving instance (Show e) => Show (Matrix e cols rows)
186+
187+-- | Type family that computes the cardinality of a given type dimension.
188+--
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)
197+ -- Generics
198+ Count (M1 _ _ f p) = Count (f p)
199+ Count (K1 _ _ _) = 1
200+ Count (V1 _) = 0
201+ Count (U1 _) = 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)
205+
206+-- | Type family that computes of a given type dimension from a given natural
207+--
208+-- Thanks to Li-Yao Xia this type family is super fast.
209+type family FromNat (n :: Nat) :: Type where
210+ FromNat 0 = Void
211+ FromNat 1 = ()
212+ FromNat n = FromNat' (Mod n 2 == 0) (FromNat (Div n 2))
213+
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)
217+
218+-- | Type family that normalizes the representation of a given data
219+-- structure
220+type family Normalize (d :: Type) :: Type where
221+ Normalize (Either a b) = Either (Normalize a) (Normalize b)
222+ Normalize d = FromNat (Count d)
223+
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))
234+
235+
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.
239+--
240+-- Please use 'iden' instead.
241+instance (Num e) => Category (Matrix e) where
242+ id = undefined
243+ (.) = comp
244+
245+instance NFData e => NFData (Matrix e cols rows) where
246+ rnf Empty = ()
247+ rnf (One e) = rnf e
248+ rnf (Join a b) = rnf a `seq` rnf b
249+ rnf (Fork a b) = rnf a `seq` rnf b
250+
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
258+
259+instance Num e => Num (Matrix e cols rows) where
260+
261+ a + b = zipWithM (+) a b
262+
263+ a - b = zipWithM (-) a b
264+
265+ a * b = zipWithM (*) a b
266+
267+ abs Empty = Empty
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 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)
276+
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
284+
285+-- Primitives
286+
287+-- | Empty matrix constructor
288+empty :: Matrix e Void Void
289+empty = Empty
290+
291+-- | Unit matrix constructor
292+one :: e -> Matrix e () ()
293+one = One
294+
295+-- | Matrix 'Join' constructor
296+join :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
297+join = Join
298+
299+infixl 3 |||
300+
301+-- | Matrix 'Join' constructor
302+(|||) :: Matrix e a rows -> Matrix e b rows -> Matrix e (Either a b) rows
303+(|||) = Join
304+
305+-- | Matrix 'Fork' constructor
306+fork :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
307+fork = Fork
308+
309+infixl 2 ===
310+
311+-- | Matrix 'Fork' constructor
312+(===) :: Matrix e cols a -> Matrix e cols b -> Matrix e cols (Either a b)
313+(===) = Fork
314+
315+-- Construction
316+
317+-- | Type class for defining the 'fromList' conversion function.
318+--
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
325+
326+instance FromLists e Void Void where
327+ fromLists [] = Empty
328+ fromLists _ = error "Wrong dimensions"
329+
330+instance {-# OVERLAPPING #-} FromLists e () () where
331+ fromLists [[e]] = One e
332+ fromLists _ = error "Wrong dimensions"
333+
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"
337+
338+instance {-# OVERLAPPABLE #-} (FromLists e a (), FromLists e b (), Countable a) => FromLists e (Either a b) () where
339+ fromLists [l] =
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"
343+
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"
347+
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"
353+
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) =
356+ let lh = length h
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"
362+
363+-- | Constructs a column vector matrix
364+colL :: (FromLists e () rows) => [e] -> Matrix e () rows
365+colL = fromLists . map (: [])
366+
367+-- | Constructs a row vector matrix
368+rowL :: (FromLists e cols ()) => [e] -> Matrix e cols ()
369+rowL = fromLists . (: [])
370+
371+-- Conversion
372+
373+-- | Converts a matrix to a list of lists of elements.
374+toLists :: Matrix e cols rows -> [[e]]
375+toLists Empty = []
376+toLists (One e) = [[e]]
377+toLists (Fork l r) = toLists l ++ toLists r
378+toLists (Join l r) = zipWith (++) (toLists l) (toLists r)
379+
380+-- | Converts a matrix to a list of elements.
381+toList :: Matrix e cols rows -> [e]
382+toList = concat . toLists
383+
384+-- Zeros Matrix
385+
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)
389+
390+-- Ones Matrix
391+
392+-- | The ones matrix. A matrix wholly filled with ones.
393+--
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)
397+
398+-- Const Matrix
399+
400+-- | The constant matrix constructor. A matrix wholly filled with a given
401+-- value.
402+constant :: (Num e, Constructable cols rows, Enumerable cols) => e -> Matrix e cols rows
403+constant e = matrixBuilder (const e)
404+
405+-- Bang Matrix
406+
407+-- | The T (Top) row vector matrix.
408+bang :: (Num e, Constructable cols rows, Enumerable cols) => Matrix e cols ()
409+bang = ones
410+
411+-- | Point constant matrix
412+point :: (Num e, Eq a, Construct a) => a -> Matrix e () a
413+point = fromF . const
414+
415+-- iden Matrix
416+
417+-- | iden matrix.
418+iden :: (Num e, Eq cols, Construct cols, Enumerable cols) => Matrix e cols cols
419+iden = fromF id
420+
421+-- Matrix composition (MMM)
422+
423+-- | Matrix composition. Equivalent to matrix-matrix multiplication.
424+--
425+-- This definition takes advantage of divide-and-conquer and fusion laws
426+-- from LAoP.
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 #-}
434+{-# RULES
435+ "comp/iden1" forall m. comp m iden = m ;
436+ "comp/iden2" forall m. comp iden m = m
437+#-}
438+
439+-- Scalar multiplication of matrices
440+
441+infixl 7 .|
442+-- | Scalar multiplication of matrices.
443+(.|) :: Num e => e -> Matrix e cols rows -> Matrix e cols rows
444+(.|) _ Empty = Empty
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)
448+
449+-- Scalar division of matrices
450+
451+infixl 7 ./
452+-- | Scalar multiplication of matrices.
453+(./) :: Fractional e => Matrix e cols rows -> e -> Matrix e cols rows
454+(./) Empty _ = Empty
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)
458+
459+-- Projections
460+
461+-- | Biproduct first component projection
462+p1 :: (Num e, Eq m, Constructable m n, Enumerable m, Enumerable n) => Matrix e (Either m n) m
463+p1 = join iden zeros
464+
465+-- | Biproduct second component projection
466+p2 :: (Num e, Eq n, Constructable m n, Enumerable m, Enumerable n) => Matrix e (Either m n) n
467+p2 = join zeros iden
468+
469+-- Injections
470+
471+-- | Biproduct first component injection
472+i1 :: (Num e, Eq m, Constructable m n, Enumerable m, Enumerable n) => Matrix e m (Either m n)
473+i1 = tr p1
474+
475+-- | Biproduct second component injection
476+i2 :: (Num e, Eq n, Constructable m n, Enumerable m, Enumerable n) => Matrix e n (Either m n)
477+i2 = tr p2
478+
479+-- Dimensions
480+
481+-- | Obtain the number of rows.
482+--
483+-- NOTE: The 'KnownNat' constaint is needed in order to obtain the
484+-- dimensions in constant time.
485+--
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))
490+
491+-- | Obtain the number of columns.
492+--
493+-- NOTE: The 'KnownNat' constaint is needed in order to obtain the
494+-- dimensions in constant time.
495+--
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))
500+
501+-- Coproduct Bifunctor
502+
503+infixl 5 -|-
504+
505+-- | Matrix coproduct functor also known as matrix direct sum.
506+(-|-) ::
507+ ( Num e,
508+ Eq k,
509+ Eq j,
510+ Constructable k j,
511+ Enumerable k,
512+ Enumerable j
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)
515+
516+-- Khatri Rao Product and projections
517+
518+-- | Khatri Rao product first component projection matrix.
519+fstM ::
520+ forall e a b .
521+ ( Num e,
522+ Eq a,
523+ ConstructableN (a, b) a,
524+ ConstructableNorm (a, b) a,
525+ Enumerable a,
526+ Enumerable b
527+ ) => Matrix e (Normalize (a, b)) (Normalize a)
528+fstM = fromFN (fst :: (a, b) -> a)
529+
530+-- | Khatri Rao product second component projection matrix.
531+sndM ::
532+ forall e a b .
533+ ( Num e,
534+ Eq b,
535+ ConstructableN (a, b) b,
536+ ConstructableNorm (a, b) b,
537+ Enumerable a,
538+ Enumerable b
539+ ) => Matrix e (Normalize (a, b)) (Normalize b)
540+sndM = fromFN (snd :: (a, b) -> b)
541+
542+-- | Khatri Rao Matrix product also known as matrix pairing.
543+--
544+-- NOTE: That this is not a true categorical product, see for instance:
545+--
546+-- @
547+-- | fstM . kr a b == a
548+-- kr a b ==> |
549+-- | sndM . kr a b == b
550+-- @
551+--
552+-- __Emphasis__ on the implication symbol.
553+kr ::
554+ forall e cols a b .
555+ ( Num e,
556+ Eq a,
557+ Eq b,
558+ ConstructableN a (a, b),
559+ ConstructN b,
560+ ConstructableNorm a b,
561+ ConstructNorm (a, b),
562+ Enumerable a,
563+ Enumerable b
564+ ) => Matrix e cols (Normalize a) -> Matrix e cols (Normalize b) -> Matrix e cols (Normalize (a, b))
565+kr a b =
566+ let fstM' = fstM @e @a @b
567+ sndM' = sndM @e @a @b
568+ in (tr fstM' . a) * (tr sndM' . b)
569+
570+-- Product Bifunctor (Kronecker)
571+
572+infixl 4 ><
573+
574+type Kronecker e m n p q =
575+ ( Num e,
576+ Eq m,
577+ Eq n,
578+ Eq p,
579+ Eq q,
580+ ConstructableN m n,
581+ ConstructableNorm m (m, n),
582+ ConstructNorm n,
583+ ConstructN (m, n),
584+ ConstructableN p q,
585+ ConstructableNorm p (p, q),
586+ ConstructNorm q,
587+ ConstructN (p, q),
588+ Enumerable m,
589+ Enumerable n,
590+ Enumerable p,
591+ Enumerable q
592+ )
593+
594+-- | Matrix product functor also known as kronecker product
595+(><) ::
596+ forall e m p n q .
597+ Kronecker e m n p q
598+ => Matrix e (Normalize m) (Normalize p)
599+ -> Matrix e (Normalize n) (Normalize q)
600+ -> Matrix e (Normalize (m, n)) (Normalize (p, q))
601+(><) a b =
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')
605+
606+-- Matrix abide Join Fork
607+
608+-- | Matrix "abiding" following the 'Join'-'Fork' exchange law.
609+--
610+-- Law:
611+--
612+-- @
613+-- 'Join' ('Fork' a c) ('Fork' b d) == 'Fork' ('Join' a b) ('Join' c d)
614+-- @
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)
621+
622+-- Matrix abide Fork Join
623+
624+-- | Matrix "abiding" followin the 'Fork'-'Join' abide law.
625+--
626+-- @
627+-- 'Fork' ('Join' a b) ('Join' c d) == 'Join' ('Fork' a c) ('Fork' b d)
628+-- @
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)
635+
636+-- Matrix transposition
637+
638+-- | Matrix transposition.
639+tr :: Matrix e cols rows -> Matrix e rows cols
640+tr Empty = Empty
641+tr (One e) = One e
642+tr (Join a b) = Fork (tr a) (tr b)
643+tr (Fork a b) = Join (tr a) (tr b)
644+
645+-- Selective 'select' operator
646+
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
653+
654+branch ::
655+ ( Num e,
656+ Eq a,
657+ Eq b,
658+ Eq c,
659+ Constructable a b,
660+ Construct c,
661+ Enumerable a,
662+ Enumerable b,
663+ Enumerable c
664+ )
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
667+ where
668+ f m = fork (tr i1) (i1 . tr i2) . m
669+ g m = i2 . m
670+
671+-- Pretty print
672+
673+prettyAux :: Show e => [[e]] -> [[e]] -> String
674+prettyAux [] _ = ""
675+prettyAux [[e]] m = "│ " ++ fill (show e) ++ " │\n"
676+ where
677+ v = fmap show m
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"
681+ where
682+ v = fmap show m
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" ++
686+ prettyAux t l
687+ where
688+ v = fmap show l
689+ widest = maximum $ fmap length v
690+ fill str = replicate (widest - length str - 2) ' ' ++ str
691+
692+-- | Matrix pretty printer
693+pretty :: (CountableDimensions cols rows, Show e) => Matrix e cols rows -> String
694+pretty m = concat
695+ [ "┌ ", unwords (replicate (columns m) blank), " ┐\n"
696+ , unlines
697+ [ "│ " ++ unwords (fmap (\j -> fill $ show $ getElem i j m) [1..columns m]) ++ " │" | i <- [1..rows m] ]
698+ , "└ ", unwords (replicate (columns m) blank), " ┘"
699+ ]
700+ where
701+ strings = map show (toList m)
702+ widest = maximum $ map length strings
703+ fill str = replicate (widest - length str) ' ' ++ str
704+ blank = fill ""
705+ safeGet i j m
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
710+ getElem i j m =
711+ fromMaybe
712+ (error $
713+ "getElem: Trying to get the "
714+ ++ show (i, j)
715+ ++ " element from a "
716+ ++ show (rows m) ++ "x" ++ show (columns m)
717+ ++ " matrix."
718+ )
719+ (safeGet i j m)
720+
721+-- | Matrix pretty printer
722+prettyPrint :: (CountableDimensions cols rows, Show e) => Matrix e cols rows -> IO ()
723+prettyPrint = putStrLn . pretty
724+
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
733+
734+-- Relational operators functions
735+
736+type Boolean = Natural 0 1
737+type Relation a b = Matrix Boolean a b
738+
739+-- | Helper conversion function
740+toBool :: (Num e, Eq e) => e -> Bool
741+toBool n
742+ | n == 0 = False
743+ | n == 1 = True
744+
745+-- | Helper conversion function
746+fromBool :: Bool -> Natural 0 1
747+fromBool True = nat 1
748+fromBool False = nat 0
749+
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)
756+
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
765+
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
774+
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
783+
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
791+
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)
799+
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))
803+
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)
807+ where
808+ intersection = andM
809+
810+-- | Lifts functions to relations with dimensions matching @a@ and @b@
811+-- cardinality's.
812+fromFRelN ::
813+ ( Eq b,
814+ ConstructableNorm a b,
815+ ConstructableN a b,
816+ Enumerable a
817+ )
818+ => (a -> b) -> Relation (Normalize a) (Normalize b)
819+fromFRelN = fromFN
820+
821+-- | Lifts a relation function to a Boolean Matrix
822+toRelN ::
823+ ( ConstructableN a b,
824+ ConstructableNorm a b,
825+ Enumerable a
826+ )
827+ => (a -> b -> Bool) -> Relation (Normalize a) (Normalize b)
828+toRelN f = matrixBuilderN (\a -> if uncurry f a then 1 else 0)
829+
830+----------------------------- Linear map semantics -----------------------------
831+newtype Vector e a = Vector { at :: a -> e }
832+
833+instance Contravariant (Vector e) where
834+ contramap f (Vector g) = Vector (g . f)
835+
836+instance Num e => Num (Vector e a) where
837+ fromInteger = Vector . const . fromInteger
838+
839+ (+) = liftV2 (+)
840+ (-) = liftV2 (-)
841+ (*) = liftV2 (*)
842+ abs = liftV1 abs
843+ negate = liftV1 negate
844+ signum = error "No sensible definition"
845+
846+liftV1 :: (e -> e) -> Vector e a -> Vector e a
847+liftV1 f x = Vector (f . at x)
848+
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))
851+
852+-- Semantics of Matrix e a b
853+type LinearMap e a b = Vector e a -> Vector e b
854+
855+semantics :: Num e => Matrix e a b -> LinearMap e a b
856+semantics m = case m of
857+ Empty -> id
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))
861+
862+padLeft :: Num e => Vector e b -> Vector e (Either a b)
863+padLeft v = Vector $ \case Left _ -> 0
864+ Right b -> at v b
865+
866+padRight :: Num e => Vector e a -> Vector e (Either a b)
867+padRight v = Vector $ \case Left a -> at v a
868+ Right _ -> 0
869+
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 ]
872+
873+--------------------------------- Construction ---------------------------------
874+
875+class Construct a where
876+ row :: Num e => Vector e a -> Matrix e a ()
877+
878+ linearMap :: (Construct b, Num e) => LinearMap e a b -> Matrix e a b
879+
880+instance Construct () where
881+ row v = one (at v ())
882+
883+ linearMap m = col (m 1)
884+
885+instance (Constructable a b) => Construct (Either a b) where
886+ row v = row (Left >$< v) ||| row (Right >$< v)
887+
888+ linearMap m = linearMap (m . padRight) ||| linearMap (m . padLeft)
889+
890+col :: (Construct a, Num e) => Vector e a -> Matrix e () a
891+col = tr . row
892+
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)
897+
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)
901+
902+-------------------------------- Deconstruction --------------------------------
903+enumerate :: Enumerable a => [a]
904+enumerate = [minBound .. maxBound]
905+
906+basis :: (Enumerable a, Eq a, Num e) => [Vector e a]
907+basis = [ Vector (bool 0 1 . (==a)) | a <- enumerate ]
908+
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 ]
912+
913+dump :: (Enumerable a, Enumerable b, Eq a, Num e, Show e) => Matrix e a b -> IO ()
914+dump = mapM_ print . toLists'
915+
916+------------------------------- Normalised types -------------------------------
917+class Profunctor p where
918+ dimap :: (a -> b) -> (c -> d) -> p b c -> p a d
919+
920+instance Profunctor (->) where
921+ dimap f g h = g . h . f
922+
923+toNorm :: ConstructNorm a => a -> Normalize a
924+toNorm = toEnum . fromEnum
925+
926+fromNorm :: ConstructNorm a => Normalize a -> a
927+fromNorm = toEnum . fromEnum
928+
929+rowN :: (ConstructN a, ConstructNorm a, Num e) => Vector e a -> Matrix e (Normalize a) ()
930+rowN = row . contramap fromNorm
931+
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)
935+
936+colN :: (ConstructN a, ConstructNorm a, Num e) => Vector e a -> Matrix e () (Normalize a)
937+colN = tr . rowN
938+
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)
944+
945+-- | Lifts functions to matrices with dimensions matching @a@ and @b@
946+-- cardinality's.
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)